FE - Finite Elements

examples/finite_elements/README.md

@@ -0,0 +1,76 @@
+# Output
+
+## Nodes and elements
+
+Output is generated as `inp` files for ABAQUS or CalculiX. Nodes and elements are saved in separate files, which are then included inside a single, final `inp` file. This final `inp` file can be opened and viewed using [FreeCAD](https://en.wikipedia.org/wiki/FreeCAD).
+
+## Loading
+
+The final `inp` file can apply a distributed *gravity* load to all elements. Additional point loads can be defined by a JSON file.
+
+## Boundary
+
+In the final `inp` file, the boundary is defined by the point restraints inside a JSON file.
+
+## Modification
+
+Currently, a convenient way for making adjustments is provided. Simply modify the load, restraint, and specification JSON files. Examining the unit tests is helpful to figure out how.
+
+# CCX and CGX
+
+The results can be consumed by ABAQUS or CalculiX.
+
+CalculiX and CalculiX GraphiX binaries are available for different platforms, like Linux distributions.
+
+## openSUSE
+
+openSUSE has [CCX](https://software.opensuse.org/package/ccx) package and also [CGX](https://software.opensuse.org/package/cgx) one.
+
+To install CCX and CGX on openSUSE Leap 15.5 you can run as root:
+
+```bash
+zypper addrepo https://download.opensuse.org/repositories/science/15.5/science.repo
+zypper refresh
+zypper install ccx
+zypper install cgx
+```
+
+# Visualize `inp` file
+
+To visualize the `inp` file by CalculiX GraphiX:
+
+```bash
+cgx -c hex8.inp
+```
+
+# Analyze `inp` file
+
+To run the `inp` files by FEA engines like CalculiX:
+
+```bash
+ccx -i hex8
+```
+
+The above `-i` flag expects a `hex8.inp` file.
+
+The above command creates `frd` files containing the results. They can be viewed by CalculiX GraphiX:
+
+```bash
+cgx hex8.frd
+```
+
+The boundary conditions and loads used in the calculation will be available together with the results if you run:
+
+```bash
+cgx hex8.frd hex8.inp
+```
+
+## Math solver
+
+The default CCX math solver is `SPOOLES` which is slow. Apparently `PARDISO` is faster and `PaStiX` is the fastest. But it's needed to build the CCX with `PARDISO` or `PaStiX` math libraries.
+
+### PARDISO
+
+#### Linux executable with the Intel Pardiso Solver
+
+You can download [here](https://www.dropbox.com/s/x8axi53l9dk9w4g/ccx_2.19_MT?dl=1) an executable with the Intel Pardiso solver for x86_64 Linux systems. The executable has all the libraries statically linked into it. So it should run by itself without any dependency. Thanks to [these guys](https://www.feacluster.com/calculix.php).

examples/finite_elements/main.go

@@ -0,0 +1,190 @@
+//-----------------------------------------------------------------------------
+/*
+
+Finite elements from triangle mesh.
+The result `inp` file is consumable by ABAQUS or CalculiX.
+
+*/
+//-----------------------------------------------------------------------------
+
+package main
+
+import (
+	"encoding/json"
+	"log"
+	"os"
+
+	"github.com/deadsy/sdfx/obj"
+	"github.com/deadsy/sdfx/render"
+	"github.com/deadsy/sdfx/sdf/finiteelements/mesh"
+	v3 "github.com/deadsy/sdfx/vec/v3"
+)
+
+type Specs struct {
+	PathResult             string // Result file, consumable by ABAQUS or CalculiX.
+	PathReport             string // Report some details after finite elements are generated.
+	PathStl                string // Input STL file.
+	PathLoadPoints         string // File containing point loads.
+	PathRestraintPoints    string // File containing point restraints.
+	MassDensity            float64
+	YoungModulus           float64
+	PoissonRatio           float64
+	GravityDirectionX      float64
+	GravityDirectionY      float64
+	GravityDirectionZ      float64
+	GravityMagnitude       float64
+	GravityIsNeeded        bool
+	Resolution             int  // Number of voxels on the longest axis of 3D model AABB.
+	NonlinearConsidered    bool // If true, nonlinear finite elements are generated.
+	ExactSurfaceConsidered bool // If true, surface is approximated by tetrahedral finite elements.
+}
+
+type Restraint struct {
+	LocX     float64
+	LocY     float64
+	LocZ     float64
+	IsFixedX bool
+	IsFixedY bool
+	IsFixedZ bool
+}
+
+type Load struct {
+	LocX float64
+	LocY float64
+	LocZ float64
+	MagX float64
+	MagY float64
+	MagZ float64
+}
+
+type Component struct {
+	VoxelCount int
+}
+
+type Report struct {
+	VoxelsX        int
+	VoxelsY        int
+	VoxelsZ        int
+	ComponentCount int
+	Components     []Component
+}
+
+// Render STL to SDF3 to finite elements.
+// Write finite elements to an `inp` file.
+// Written file can be used by ABAQUS or CalculiX.
+func main() {
+	if len(os.Args) != 2 {
+		log.Fatalf("usage: wrong argument count")
+	}
+
+	pthSpecs := os.Args[1]
+
+	jsonData, err := os.ReadFile(pthSpecs)
+	if err != nil {
+		log.Fatalf(err.Error())
+	}
+
+	var specs Specs
+	err = json.Unmarshal(jsonData, &specs)
+	if err != nil {
+		log.Fatalf(err.Error())
+	}
+
+	jsonData, err = os.ReadFile(specs.PathLoadPoints)
+	if err != nil {
+		log.Fatalf(err.Error())
+	}
+
+	var loads []Load
+	err = json.Unmarshal(jsonData, &loads)
+	if err != nil {
+		log.Fatalf(err.Error())
+	}
+
+	jsonData, err = os.ReadFile(specs.PathRestraintPoints)
+	if err != nil {
+		log.Fatalf(err.Error())
+	}
+
+	var restraints []Restraint
+	err = json.Unmarshal(jsonData, &restraints)
+	if err != nil {
+		log.Fatalf(err.Error())
+	}
+
+	// create the SDF from the STL mesh
+	inSdf, err := obj.ImportSTL(specs.PathStl, 20, 3, 5)
+	if err != nil {
+		log.Fatalf(err.Error())
+	}
+
+	var order render.Order
+	if specs.NonlinearConsidered {
+		order = render.Quadratic
+	} else {
+		order = render.Linear
+	}
+
+	var shape render.Shape
+	if specs.ExactSurfaceConsidered {
+		shape = render.HexAndTet
+	} else {
+		shape = render.Hexahedral
+	}
+
+	// Create a mesh of finite elements.
+	m, voxelsX, voxelsY, voxelsZ := mesh.NewFem(inSdf, render.NewMarchingCubesFeUniform(specs.Resolution, order, shape))
+
+	components := m.Components()
+	report := Report{
+		VoxelsX:        voxelsX,
+		VoxelsY:        voxelsY,
+		VoxelsZ:        voxelsZ,
+		ComponentCount: len(components),
+		Components:     make([]Component, len(components)),
+	}
+	for i, component := range components {
+		report.Components[i] = Component{VoxelCount: component.VoxelCount()}
+	}
+
+	jsonData, err = json.MarshalIndent(report, "", "    ")
+	if err != nil {
+		log.Fatalf(err.Error())
+	}
+	err = os.WriteFile(specs.PathReport, jsonData, 0644)
+	if err != nil {
+		log.Fatalf(err.Error())
+	}
+
+	// Generate finite elements for all layers of mesh.
+	err = m.WriteInp(
+		specs.PathResult,
+		float32(specs.MassDensity), float32(specs.YoungModulus), float32(specs.PoissonRatio),
+		restraintsConvert(restraints),
+		loadsConvert(loads),
+		v3.Vec{X: specs.GravityDirectionX, Y: specs.GravityDirectionY, Z: specs.GravityDirectionZ}, specs.GravityMagnitude,
+		specs.GravityIsNeeded,
+	)
+
+	if err != nil {
+		log.Fatalf(err.Error())
+	}
+}
+
+func restraintsConvert(rs []Restraint) []*mesh.Restraint {
+	restraints := make([]*mesh.Restraint, len(rs))
+	for i, r := range rs {
+		restraint := mesh.NewRestraint(v3.Vec{X: r.LocX, Y: r.LocY, Z: r.LocZ}, r.IsFixedX, r.IsFixedY, r.IsFixedZ)
+		restraints[i] = restraint
+	}
+	return restraints
+}
+
+func loadsConvert(ls []Load) []*mesh.Load {
+	loads := make([]*mesh.Load, len(ls))
+	for i, l := range ls {
+		load := mesh.NewLoad(v3.Vec{X: l.LocX, Y: l.LocY, Z: l.LocZ}, v3.Vec{X: l.MagX, Y: l.MagY, Z: l.MagZ})
+		loads[i] = load
+	}
+	return loads
+}