diff --git a/examples/finite_elements/README.md b/examples/finite_elements/README.md new file mode 100644 index 000000000..80d682a0f --- /dev/null +++ b/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). diff --git a/examples/finite_elements/main.go b/examples/finite_elements/main.go new file mode 100644 index 000000000..2db3d4c2b --- /dev/null +++ b/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 +} diff --git a/examples/finite_elements/main_test.go b/examples/finite_elements/main_test.go new file mode 100644 index 000000000..50f2505f1 --- /dev/null +++ b/examples/finite_elements/main_test.go @@ -0,0 +1,322 @@ +//----------------------------------------------------------------------------- + +/* + +Finite elements from triangle mesh. +The result `inp` file is consumable by ABAQUS or CalculiX. + +*/ + +//----------------------------------------------------------------------------- + +package main + +import ( + "encoding/json" + "math" + "os" + "path/filepath" + "testing" +) + +// Benchmark reference: +// https://github.com/calculix/CalculiX-Examples/tree/master/NonLinear/Sections +func Test_main(t *testing.T) { + tests := []struct { + skip bool + name string + pathSpecs string // File to be created by test. + specs Specs + loads []Load // If load is zero, gravity would be the dominant force. + restraints []Restraint + }{ + { + skip: false, + name: "teapot", + pathSpecs: filepath.Join(os.TempDir(), "specs.json"), + specs: Specs{ + PathResult: filepath.Join(os.TempDir(), "result.inp"), + PathReport: filepath.Join(os.TempDir(), "report.json"), + PathStl: filepath.Join("..", "..", "files", "teapot.stl"), // Valid STL, Unit: mm + PathLoadPoints: filepath.Join(os.TempDir(), "load-points.json"), + PathRestraintPoints: filepath.Join(os.TempDir(), "restraint-points.json"), + MassDensity: 1130 * math.Pow(10, -12), // (N*s2/mm4) // Assumed: 1.13 g/cm3 + YoungModulus: 1.6 * 1000, // MPa (N/mm2) + PoissonRatio: 0.3, // Unitless. + GravityDirectionX: 0, + GravityDirectionY: 0, + GravityDirectionZ: -1, + GravityMagnitude: 9810, // mm/s2 + GravityIsNeeded: false, + Resolution: 60, + NonlinearConsidered: false, + ExactSurfaceConsidered: true, + }, + loads: []Load{ + {LocX: -7.7018147506213062, LocY: -0.4793329364029888, LocZ: 5.4655784011739659, MagX: -70.381474830032147, MagY: -174.42493975029208, MagZ: 59.390099907428898}, + {LocX: -0.011008272390835461, LocY: -0.7768798803556729, LocZ: 8.0940818810755175, MagX: -7.1696819796276845, MagY: -157.24707657594607, MagZ: -122.86811489950169}, + {LocX: 7.7771501865767299, LocY: -0.44676917365822177, LocZ: 6.1957182567021745, MagX: 8.5251191596139506, MagY: -198.29032531361477, MagZ: 18.196364257032524}, + }, + restraints: []Restraint{ + {LocX: 2.6121906631017695, LocY: 0.20348199936959829, LocZ: 0.050483960817894413, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: -1.3968227044257533, LocY: -2.035934011608322, LocZ: 0.04909315835598238, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: -1.8197506822193277, LocY: 2.2580011513606717, LocZ: 0.064527793304306025, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + }, + }, + { + skip: true, + name: "benchmarkSquare", + pathSpecs: filepath.Join(os.TempDir(), "bms-specs.json"), + specs: Specs{ + PathResult: filepath.Join(os.TempDir(), "bms-result.inp"), + PathReport: filepath.Join(os.TempDir(), "bms-report.json"), + PathStl: filepath.Join("..", "..", "files", "benchmark-square.stl"), + PathLoadPoints: filepath.Join(os.TempDir(), "bms-load-points.json"), + PathRestraintPoints: filepath.Join(os.TempDir(), "bms-restraint-points.json"), + MassDensity: 7.85e-9, + YoungModulus: 210000, + PoissonRatio: 0.3, + GravityDirectionX: 0, + GravityDirectionY: 0, + GravityDirectionZ: -1, + GravityMagnitude: 9810, + GravityIsNeeded: true, + Resolution: 50, + NonlinearConsidered: false, + ExactSurfaceConsidered: true, + }, + loads: []Load{ + { + LocX: 0, + LocY: 0, + LocZ: 0, + MagX: 0, + MagY: 0, + MagZ: 0, + }, + }, + restraints: func() []Restraint { + restraints := make([]Restraint, 0) + + gap := 1.0 + var y float64 + for y <= 17.32 { + restraint := Restraint{ + LocX: 0, + LocY: y, + LocZ: 0, + IsFixedX: true, + IsFixedY: true, + IsFixedZ: true, + } + restraints = append(restraints, restraint) + y += gap + } + + y = 0 + for y <= 17.32 { + restraint := Restraint{ + LocX: 200, + LocY: y, + LocZ: 0, + IsFixedX: false, + IsFixedY: true, + IsFixedZ: true, + } + restraints = append(restraints, restraint) + y += gap + } + return restraints + }(), + }, + { + skip: true, + name: "benchmarkCircle", + pathSpecs: filepath.Join(os.TempDir(), "bmc-specs.json"), + specs: Specs{ + PathResult: filepath.Join(os.TempDir(), "bmc-result.inp"), + PathReport: filepath.Join(os.TempDir(), "bmc-report.json"), + PathStl: filepath.Join("..", "..", "files", "benchmark-circle.stl"), + PathLoadPoints: filepath.Join(os.TempDir(), "bmc-load-points.json"), + PathRestraintPoints: filepath.Join(os.TempDir(), "bmc-restraint-points.json"), + MassDensity: 7.85e-9, + YoungModulus: 210000, + PoissonRatio: 0.3, + GravityDirectionX: 0, + GravityDirectionY: 0, + GravityDirectionZ: -1, + GravityMagnitude: 9810, + GravityIsNeeded: true, + Resolution: 50, + NonlinearConsidered: false, + ExactSurfaceConsidered: true, + }, + loads: []Load{ + { + LocX: 0, + LocY: 0, + LocZ: 0, + MagX: 0, + MagY: 0, + MagZ: 0, + }, + }, + restraints: []Restraint{ + {LocX: 0, LocY: 0, LocZ: 0, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: 0, LocY: -2.0313, LocZ: 0.213498, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: 0, LocY: -3.97382, LocZ: 0.844661, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: 0, LocY: 2.0313, LocZ: 0.213498, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: 0, LocY: 3.97382, LocZ: 0.844661, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: 200, LocY: 0, LocZ: 0, IsFixedX: false, IsFixedY: true, IsFixedZ: true}, + {LocX: 200, LocY: -2.0313, LocZ: 0.213498, IsFixedX: false, IsFixedY: true, IsFixedZ: true}, + {LocX: 200, LocY: -3.97382, LocZ: 0.844661, IsFixedX: false, IsFixedY: true, IsFixedZ: true}, + {LocX: 200, LocY: 2.0313, LocZ: 0.213498, IsFixedX: false, IsFixedY: true, IsFixedZ: true}, + {LocX: 200, LocY: 3.97382, LocZ: 0.844661, IsFixedX: false, IsFixedY: true, IsFixedZ: true}, + }, + }, + { + skip: true, + name: "benchmarkPipe", + pathSpecs: filepath.Join(os.TempDir(), "bmp-specs.json"), + specs: Specs{ + PathResult: filepath.Join(os.TempDir(), "bmp-result.inp"), + PathReport: filepath.Join(os.TempDir(), "bmp-report.json"), + PathStl: filepath.Join("..", "..", "files", "benchmark-pipe.stl"), + PathLoadPoints: filepath.Join(os.TempDir(), "bmp-load-points.json"), + PathRestraintPoints: filepath.Join(os.TempDir(), "bmp-restraint-points.json"), + MassDensity: 7.85e-9, + YoungModulus: 210000, + PoissonRatio: 0.3, + GravityDirectionX: 0, + GravityDirectionY: 0, + GravityDirectionZ: -1, + GravityMagnitude: 9810, + GravityIsNeeded: true, + Resolution: 50, + NonlinearConsidered: false, + ExactSurfaceConsidered: true, + }, + loads: []Load{ + { + LocX: 0, + LocY: 0, + LocZ: 0, + MagX: 0, + MagY: 0, + MagZ: 0, + }, + }, + restraints: []Restraint{ + {LocX: 0, LocY: 0, LocZ: 0, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: 0, LocY: -2.0313, LocZ: 0.213498, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: 0, LocY: -3.97382, LocZ: 0.844661, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: 0, LocY: 2.0313, LocZ: 0.213498, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: 0, LocY: 3.97382, LocZ: 0.844661, IsFixedX: true, IsFixedY: true, IsFixedZ: true}, + {LocX: 200, LocY: 0, LocZ: 0, IsFixedX: false, IsFixedY: true, IsFixedZ: true}, + {LocX: 200, LocY: -2.0313, LocZ: 0.213498, IsFixedX: false, IsFixedY: true, IsFixedZ: true}, + {LocX: 200, LocY: -3.97382, LocZ: 0.844661, IsFixedX: false, IsFixedY: true, IsFixedZ: true}, + {LocX: 200, LocY: 2.0313, LocZ: 0.213498, IsFixedX: false, IsFixedY: true, IsFixedZ: true}, + {LocX: 200, LocY: 3.97382, LocZ: 0.844661, IsFixedX: false, IsFixedY: true, IsFixedZ: true}, + }, + }, + { + skip: true, + name: "benchmarkI", + pathSpecs: filepath.Join(os.TempDir(), "bmi-specs.json"), + specs: Specs{ + PathResult: filepath.Join(os.TempDir(), "bmi-result.inp"), + PathReport: filepath.Join(os.TempDir(), "bmi-report.json"), + PathStl: filepath.Join("..", "..", "files", "benchmark-I.stl"), + PathLoadPoints: filepath.Join(os.TempDir(), "bmi-load-points.json"), + PathRestraintPoints: filepath.Join(os.TempDir(), "bmi-restraint-points.json"), + MassDensity: 7.85e-9, + YoungModulus: 210000, + PoissonRatio: 0.3, + GravityDirectionX: 0, + GravityDirectionY: 0, + GravityDirectionZ: -1, + GravityMagnitude: 9810, + GravityIsNeeded: true, + Resolution: 50, + NonlinearConsidered: false, + ExactSurfaceConsidered: true, + }, + loads: []Load{ + { + LocX: 0, + LocY: 0, + LocZ: 0, + MagX: 0, + MagY: 0, + MagZ: 0, + }, + }, + restraints: func() []Restraint { + restraints := make([]Restraint, 0) + + gap := 1.0 + var y float64 + for y <= 25 { + restraints = append(restraints, Restraint{LocX: 0, LocY: y, LocZ: 0, IsFixedX: true, IsFixedY: true, IsFixedZ: true}) + y += gap + } + + y = 0 + for y <= 25 { + restraints = append(restraints, Restraint{LocX: 200, LocY: y, LocZ: 0, IsFixedX: false, IsFixedY: true, IsFixedZ: true}) + y += gap + } + + return restraints + }(), + }, + } + + for _, tt := range tests { + if tt.skip { + continue + } + t.Run(tt.name, func(t *testing.T) { + jsonData, err := json.MarshalIndent(tt.specs, "", " ") + if err != nil { + t.Error(err) + return + } + err = os.WriteFile(tt.pathSpecs, jsonData, 0644) + if err != nil { + t.Error(err) + return + } + + jsonData, err = json.MarshalIndent(tt.loads, "", " ") + if err != nil { + t.Error(err) + return + } + + err = os.WriteFile(tt.specs.PathLoadPoints, jsonData, 0644) + if err != nil { + t.Error(err) + return + } + + jsonData, err = json.MarshalIndent(tt.restraints, "", " ") + if err != nil { + t.Error(err) + return + } + + err = os.WriteFile(tt.specs.PathRestraintPoints, jsonData, 0644) + if err != nil { + t.Error(err) + return + } + + os.Args = []string{ + "executable-name-dummy", + tt.pathSpecs, + } + main() + }) + } +} diff --git a/examples/finite_elements_print3d/README.md b/examples/finite_elements_print3d/README.md new file mode 100644 index 000000000..1180e5d23 --- /dev/null +++ b/examples/finite_elements_print3d/README.md @@ -0,0 +1,5 @@ +# Finite elements for 3D print analysis + +## Layer-by-layer along Z axis + +The result is generated layer-by-layer along the Z axis. Just like the 3D print process itself. The result is in the form of `inp` files for ABAQUS or CalculiX diff --git a/examples/finite_elements_print3d/main.go b/examples/finite_elements_print3d/main.go new file mode 100644 index 000000000..972a341ae --- /dev/null +++ b/examples/finite_elements_print3d/main.go @@ -0,0 +1,174 @@ +//----------------------------------------------------------------------------- +/* + +Finite elements from triangle mesh. +The result `inp` file is consumable by ABAQUS or CalculiX. + +*/ +//----------------------------------------------------------------------------- + +package main + +import ( + "encoding/json" + "fmt" + "log" + "os" + "strings" + + "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 { + PathResultWithPlaceholder string // Result file, consumable by ABAQUS or CalculiX. Must include "#" character as placeholder for layer number. + PathReport string // Report some details after finite elements are generated. + PathStl string // Input STL file. + MassDensity float64 + YoungModulus float64 + PoissonRatio float64 + GravityDirectionX float64 + GravityDirectionY float64 + GravityDirectionZ float64 + GravityMagnitude float64 + 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()) + } + + // 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 layer-by-layer. + // Applicable to 3D print analysis that is done layer-by-layer. + + if voxelsZ < 3 { + log.Fatalf("not enough voxel layers along the Z axis %d < %d", voxelsZ, 3) + } + + // Skip 0 and start from 1, to be sensible. + for z := 1; z <= voxelsZ; z++ { + err := m.WriteInpLayers( + strings.Replace(specs.PathResultWithPlaceholder, "#", fmt.Sprintf("%d", z), 1), + 0, z, // Note that the start layer is included, the end layer is excluded. + float32(specs.MassDensity), float32(specs.YoungModulus), float32(specs.PoissonRatio), + restraintsPrintFloor(m), + []*mesh.Load{}, // Load is empty since only gravity is assumed responsible for 3D print collapse. + v3.Vec{X: specs.GravityDirectionX, Y: specs.GravityDirectionY, Z: specs.GravityDirectionZ}, specs.GravityMagnitude, + true, + ) + if err != nil { + log.Fatalf(err.Error()) + } + fmt.Printf("Finite elements are generated from layer 0 to layer %v out of %v total.\n", z, voxelsZ) + } + + if err != nil { + log.Fatalf(err.Error()) + } +} + +// For 3D print analysis, all the voxels at the first layer along Z axis are considered as restraint. +// Since, the 3D print floor is at the first Z level. +func restraintsPrintFloor(m *mesh.Fem) []*mesh.Restraint { + voxels := m.VoxelsOn1stLayerZ() + restraints := make([]*mesh.Restraint, len(voxels)) + for i, voxel := range voxels { + restraints[i] = mesh.NewRestraintByVoxel(voxel, true, true, true) + } + return restraints +} diff --git a/examples/finite_elements_print3d/main_test.go b/examples/finite_elements_print3d/main_test.go new file mode 100644 index 000000000..2f6a4d0d0 --- /dev/null +++ b/examples/finite_elements_print3d/main_test.go @@ -0,0 +1,70 @@ +//----------------------------------------------------------------------------- + +/* + +Finite elements from triangle mesh. +The result `inp` file is consumable by ABAQUS or CalculiX. + +*/ + +//----------------------------------------------------------------------------- + +package main + +import ( + "encoding/json" + "os" + "path/filepath" + "testing" +) + +func Test_main(t *testing.T) { + tests := []struct { + skip bool + name string + pathSpecs string // File to be created by test. + specs Specs + }{ + { + skip: false, + name: "teapot", + pathSpecs: filepath.Join(os.TempDir(), "teapot-specs.json"), + specs: Specs{ + PathResultWithPlaceholder: filepath.Join(os.TempDir(), "teapot-result-layer0-to-layer#.inp"), + PathReport: filepath.Join(os.TempDir(), "teapot-report.json"), + PathStl: filepath.Join("..", "..", "files", "teapot.stl"), + MassDensity: 7.85e-9, + YoungModulus: 210000, + PoissonRatio: 0.3, + GravityDirectionX: 0, + GravityDirectionY: 0, + GravityDirectionZ: +1, // SLA 3D print is usually done upside down. + GravityMagnitude: 9810, // mm unit. + Resolution: 50, + NonlinearConsidered: false, + ExactSurfaceConsidered: true, + }, + }, + } + + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + jsonData, err := json.MarshalIndent(tt.specs, "", " ") + if err != nil { + t.Error(err) + return + } + err = os.WriteFile(tt.pathSpecs, jsonData, 0644) + if err != nil { + t.Error(err) + return + } + + os.Args = []string{ + "executable-name-dummy", + tt.pathSpecs, + } + main() + }) + } +} diff --git a/examples/stl_to_sdf_evaluation/main.go b/examples/stl_to_sdf_evaluation/main.go new file mode 100644 index 000000000..87e0ad825 --- /dev/null +++ b/examples/stl_to_sdf_evaluation/main.go @@ -0,0 +1,26 @@ +package main + +import ( + "fmt" + "log" + + "github.com/deadsy/sdfx/obj" + v3 "github.com/deadsy/sdfx/vec/v3" +) + +func main() { + // create the SDF from the STL file mesh. + inSdf, err := obj.ImportSTL("../../files/teapot.stl", 20, 3, 5) + if err != nil { + log.Fatalf("error: %s", err) + } + + // This point is definitely inside the teapot model, + // so SDF value should be negative. + value := inSdf.Evaluate(v3.Vec{X: -0.8164382918936324, Y: 2.542909114087213, Z: 5.006102143191411}) + if value >= 0 { + fmt.Println("not expected") + } else { + fmt.Println("as expected") + } +} diff --git a/files/benchmark-I.scad b/files/benchmark-I.scad new file mode 100644 index 000000000..a202d6782 --- /dev/null +++ b/files/benchmark-I.scad @@ -0,0 +1,16 @@ +t = 4; // Thickness +length = 200; // Length of the beam +width = 25; // Width of the beam +height = 25 + t + t; // Overall height of the beam + +flange_height = t; // Height of the top and bottom flanges +web_height = height - 2 * flange_height; // Height of the vertical web + +// Top flange +translate([ 0, 0, height - flange_height ]) cube([ length, width, flange_height ]); + +// Bottom flange +cube([ length, width, flange_height ]); + +// Vertical web +translate([ 0, (width - t) / 2, flange_height ]) cube([ length, t, web_height ]); \ No newline at end of file diff --git a/files/benchmark-I.stl b/files/benchmark-I.stl new file mode 100644 index 000000000..73169f864 --- /dev/null +++ b/files/benchmark-I.stl @@ -0,0 +1,310 @@ +solid OpenSCAD_Model + facet normal 1 0 0 + outer loop + vertex 200 25 33 + vertex 200 14.5 29 + vertex 200 25 29 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 14.5 29 + vertex 200 10.5 29 + vertex 200 14.5 4 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 25 33 + vertex 200 10.5 29 + vertex 200 14.5 29 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 0 33 + vertex 200 10.5 29 + vertex 200 25 33 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 10.5 29 + vertex 200 0 33 + vertex 200 0 29 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 14.5 4 + vertex 200 25 0 + vertex 200 25 4 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 10.5 4 + vertex 200 14.5 4 + vertex 200 10.5 29 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 10.5 4 + vertex 200 25 0 + vertex 200 14.5 4 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 0 0 + vertex 200 10.5 4 + vertex 200 0 4 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 10.5 4 + vertex 200 0 0 + vertex 200 25 0 + endloop + endfacet + facet normal -0 0 1 + outer loop + vertex 0 25 33 + vertex 200 0 33 + vertex 200 25 33 + endloop + endfacet + facet normal 0 0 1 + outer loop + vertex 200 0 33 + vertex 0 25 33 + vertex 0 0 33 + endloop + endfacet + facet normal 0 0 -1 + outer loop + vertex 0 14.5 29 + vertex 200 25 29 + vertex 200 14.5 29 + endloop + endfacet + facet normal -0 0 -1 + outer loop + vertex 200 25 29 + vertex 0 14.5 29 + vertex 0 25 29 + endloop + endfacet + facet normal 0 0 -1 + outer loop + vertex 0 0 29 + vertex 200 10.5 29 + vertex 200 0 29 + endloop + endfacet + facet normal -0 0 -1 + outer loop + vertex 200 10.5 29 + vertex 0 0 29 + vertex 0 10.5 29 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 25 0 + vertex 0 14.5 4 + vertex 0 25 4 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 14.5 4 + vertex 0 10.5 4 + vertex 0 14.5 29 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 25 0 + vertex 0 10.5 4 + vertex 0 14.5 4 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 0 + vertex 0 10.5 4 + vertex 0 25 0 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 10.5 4 + vertex 0 0 0 + vertex 0 0 4 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 14.5 29 + vertex 0 25 33 + vertex 0 25 29 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 10.5 29 + vertex 0 14.5 29 + vertex 0 10.5 4 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 10.5 29 + vertex 0 25 33 + vertex 0 14.5 29 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 33 + vertex 0 10.5 29 + vertex 0 0 29 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 10.5 29 + vertex 0 0 33 + vertex 0 25 33 + endloop + endfacet + facet normal 0 1 -0 + outer loop + vertex 200 25 29 + vertex 0 25 33 + vertex 200 25 33 + endloop + endfacet + facet normal 0 1 0 + outer loop + vertex 0 25 33 + vertex 200 25 29 + vertex 0 25 29 + endloop + endfacet + facet normal 0 -1 0 + outer loop + vertex 0 0 29 + vertex 200 0 33 + vertex 0 0 33 + endloop + endfacet + facet normal 0 -1 -0 + outer loop + vertex 200 0 33 + vertex 0 0 29 + vertex 200 0 29 + endloop + endfacet + facet normal -0 0 1 + outer loop + vertex 0 25 4 + vertex 200 14.5 4 + vertex 200 25 4 + endloop + endfacet + facet normal 0 0 1 + outer loop + vertex 200 14.5 4 + vertex 0 25 4 + vertex 0 14.5 4 + endloop + endfacet + facet normal -0 0 1 + outer loop + vertex 0 10.5 4 + vertex 200 0 4 + vertex 200 10.5 4 + endloop + endfacet + facet normal 0 0 1 + outer loop + vertex 200 0 4 + vertex 0 10.5 4 + vertex 0 0 4 + endloop + endfacet + facet normal 0 0 -1 + outer loop + vertex 0 0 0 + vertex 200 25 0 + vertex 200 0 0 + endloop + endfacet + facet normal -0 0 -1 + outer loop + vertex 200 25 0 + vertex 0 0 0 + vertex 0 25 0 + endloop + endfacet + facet normal 0 1 -0 + outer loop + vertex 200 25 0 + vertex 0 25 4 + vertex 200 25 4 + endloop + endfacet + facet normal 0 1 0 + outer loop + vertex 0 25 4 + vertex 200 25 0 + vertex 0 25 0 + endloop + endfacet + facet normal 0 -1 0 + outer loop + vertex 0 0 0 + vertex 200 0 4 + vertex 0 0 4 + endloop + endfacet + facet normal 0 -1 -0 + outer loop + vertex 200 0 4 + vertex 0 0 0 + vertex 200 0 0 + endloop + endfacet + facet normal 0 1 -0 + outer loop + vertex 200 14.5 4 + vertex 0 14.5 29 + vertex 200 14.5 29 + endloop + endfacet + facet normal 0 1 0 + outer loop + vertex 0 14.5 29 + vertex 200 14.5 4 + vertex 0 14.5 4 + endloop + endfacet + facet normal 0 -1 0 + outer loop + vertex 0 10.5 4 + vertex 200 10.5 29 + vertex 0 10.5 29 + endloop + endfacet + facet normal 0 -1 -0 + outer loop + vertex 200 10.5 29 + vertex 0 10.5 4 + vertex 200 10.5 4 + endloop + endfacet +endsolid OpenSCAD_Model diff --git a/files/benchmark-circle.scad b/files/benchmark-circle.scad new file mode 100644 index 000000000..b33360fc3 --- /dev/null +++ b/files/benchmark-circle.scad @@ -0,0 +1,4 @@ +length = 200; +diameter = 19.54; + +translate([ 0, 0, diameter / 2 ]) rotate([ 0, 90, 0 ]) cylinder(h = length, d = diameter); \ No newline at end of file diff --git a/files/benchmark-circle.stl b/files/benchmark-circle.stl new file mode 100644 index 000000000..f5be69c2e --- /dev/null +++ b/files/benchmark-circle.stl @@ -0,0 +1,814 @@ +solid OpenSCAD_Model + facet normal 0 0.104528 -0.994522 + outer loop + vertex 0 0 0 + vertex 200 2.0313 0.213498 + vertex 200 0 0 + endloop + endfacet + facet normal 0 0.104528 -0.994522 + outer loop + vertex 200 2.0313 0.213498 + vertex 0 0 0 + vertex 0 2.0313 0.213498 + endloop + endfacet + facet normal 0 0.309017 -0.951057 + outer loop + vertex 0 2.0313 0.213498 + vertex 200 3.97382 0.844661 + vertex 200 2.0313 0.213498 + endloop + endfacet + facet normal 0 0.309017 -0.951057 + outer loop + vertex 200 3.97382 0.844661 + vertex 0 2.0313 0.213498 + vertex 0 3.97382 0.844661 + endloop + endfacet + facet normal 0 0.5 -0.866026 + outer loop + vertex 0 3.97382 0.844661 + vertex 200 5.74266 1.8659 + vertex 200 3.97382 0.844661 + endloop + endfacet + facet normal 0 0.5 -0.866026 + outer loop + vertex 200 5.74266 1.8659 + vertex 0 3.97382 0.844661 + vertex 0 5.74266 1.8659 + endloop + endfacet + facet normal 0 0.669131 -0.743144 + outer loop + vertex 0 5.74266 1.8659 + vertex 200 7.26052 3.23259 + vertex 200 5.74266 1.8659 + endloop + endfacet + facet normal 0 0.669131 -0.743144 + outer loop + vertex 200 7.26052 3.23259 + vertex 0 5.74266 1.8659 + vertex 0 7.26052 3.23259 + endloop + endfacet + facet normal 0 0.809016 -0.587786 + outer loop + vertex 200 7.26052 3.23259 + vertex 0 8.46107 4.885 + vertex 200 8.46107 4.885 + endloop + endfacet + facet normal 0 0.809016 -0.587786 + outer loop + vertex 0 8.46107 4.885 + vertex 200 7.26052 3.23259 + vertex 0 7.26052 3.23259 + endloop + endfacet + facet normal 0 0.913546 -0.406736 + outer loop + vertex 200 8.46107 4.885 + vertex 0 9.29182 6.7509 + vertex 200 9.29182 6.7509 + endloop + endfacet + facet normal 0 0.913546 -0.406736 + outer loop + vertex 0 9.29182 6.7509 + vertex 200 8.46107 4.885 + vertex 0 8.46107 4.885 + endloop + endfacet + facet normal 0 0.978147 -0.207913 + outer loop + vertex 200 9.29182 6.7509 + vertex 0 9.71648 8.74876 + vertex 200 9.71648 8.74876 + endloop + endfacet + facet normal 0 0.978147 -0.207913 + outer loop + vertex 0 9.71648 8.74876 + vertex 200 9.29182 6.7509 + vertex 0 9.29182 6.7509 + endloop + endfacet + facet normal 0 1 -0 + outer loop + vertex 200 9.71648 8.74876 + vertex 0 9.71648 10.7912 + vertex 200 9.71648 10.7912 + endloop + endfacet + facet normal 0 1 0 + outer loop + vertex 0 9.71648 10.7912 + vertex 200 9.71648 8.74876 + vertex 0 9.71648 8.74876 + endloop + endfacet + facet normal 0 0.978148 0.207909 + outer loop + vertex 200 9.71648 10.7912 + vertex 0 9.29182 12.7891 + vertex 200 9.29182 12.7891 + endloop + endfacet + facet normal 0 0.978148 0.207909 + outer loop + vertex 0 9.29182 12.7891 + vertex 200 9.71648 10.7912 + vertex 0 9.71648 10.7912 + endloop + endfacet + facet normal 0 0.913546 0.406736 + outer loop + vertex 200 9.29182 12.7891 + vertex 0 8.46107 14.655 + vertex 200 8.46107 14.655 + endloop + endfacet + facet normal 0 0.913546 0.406736 + outer loop + vertex 0 8.46107 14.655 + vertex 200 9.29182 12.7891 + vertex 0 9.29182 12.7891 + endloop + endfacet + facet normal 0 0.809014 0.587789 + outer loop + vertex 200 8.46107 14.655 + vertex 0 7.26052 16.3074 + vertex 200 7.26052 16.3074 + endloop + endfacet + facet normal 0 0.809014 0.587789 + outer loop + vertex 0 7.26052 16.3074 + vertex 200 8.46107 14.655 + vertex 0 8.46107 14.655 + endloop + endfacet + facet normal -0 0.669134 0.743142 + outer loop + vertex 0 7.26052 16.3074 + vertex 200 5.74266 17.6741 + vertex 200 7.26052 16.3074 + endloop + endfacet + facet normal 0 0.669134 0.743142 + outer loop + vertex 200 5.74266 17.6741 + vertex 0 7.26052 16.3074 + vertex 0 5.74266 17.6741 + endloop + endfacet + facet normal -0 0.499985 0.866034 + outer loop + vertex 0 5.74266 17.6741 + vertex 200 3.97382 18.6953 + vertex 200 5.74266 17.6741 + endloop + endfacet + facet normal 0 0.499985 0.866034 + outer loop + vertex 200 3.97382 18.6953 + vertex 0 5.74266 17.6741 + vertex 0 3.97382 18.6953 + endloop + endfacet + facet normal -0 0.309033 0.951051 + outer loop + vertex 0 3.97382 18.6953 + vertex 200 2.0313 19.3265 + vertex 200 3.97382 18.6953 + endloop + endfacet + facet normal 0 0.309033 0.951051 + outer loop + vertex 200 2.0313 19.3265 + vertex 0 3.97382 18.6953 + vertex 0 2.0313 19.3265 + endloop + endfacet + facet normal -0 0.104529 0.994522 + outer loop + vertex 0 2.0313 19.3265 + vertex 200 0 19.54 + vertex 200 2.0313 19.3265 + endloop + endfacet + facet normal 0 0.104529 0.994522 + outer loop + vertex 200 0 19.54 + vertex 0 2.0313 19.3265 + vertex 0 0 19.54 + endloop + endfacet + facet normal 0 -0.104529 0.994522 + outer loop + vertex 0 0 19.54 + vertex 200 -2.0313 19.3265 + vertex 200 0 19.54 + endloop + endfacet + facet normal 0 -0.104529 0.994522 + outer loop + vertex 200 -2.0313 19.3265 + vertex 0 0 19.54 + vertex 0 -2.0313 19.3265 + endloop + endfacet + facet normal 0 -0.309033 0.951051 + outer loop + vertex 0 -2.0313 19.3265 + vertex 200 -3.97382 18.6953 + vertex 200 -2.0313 19.3265 + endloop + endfacet + facet normal 0 -0.309033 0.951051 + outer loop + vertex 200 -3.97382 18.6953 + vertex 0 -2.0313 19.3265 + vertex 0 -3.97382 18.6953 + endloop + endfacet + facet normal 0 -0.499985 0.866034 + outer loop + vertex 0 -3.97382 18.6953 + vertex 200 -5.74266 17.6741 + vertex 200 -3.97382 18.6953 + endloop + endfacet + facet normal 0 -0.499985 0.866034 + outer loop + vertex 200 -5.74266 17.6741 + vertex 0 -3.97382 18.6953 + vertex 0 -5.74266 17.6741 + endloop + endfacet + facet normal 0 -0.669134 0.743142 + outer loop + vertex 0 -5.74266 17.6741 + vertex 200 -7.26052 16.3074 + vertex 200 -5.74266 17.6741 + endloop + endfacet + facet normal 0 -0.669134 0.743142 + outer loop + vertex 200 -7.26052 16.3074 + vertex 0 -5.74266 17.6741 + vertex 0 -7.26052 16.3074 + endloop + endfacet + facet normal 0 -0.809014 0.587789 + outer loop + vertex 0 -8.46107 14.655 + vertex 200 -7.26052 16.3074 + vertex 0 -7.26052 16.3074 + endloop + endfacet + facet normal 0 -0.809014 0.587789 + outer loop + vertex 200 -7.26052 16.3074 + vertex 0 -8.46107 14.655 + vertex 200 -8.46107 14.655 + endloop + endfacet + facet normal 0 -0.913546 0.406736 + outer loop + vertex 0 -9.29182 12.7891 + vertex 200 -8.46107 14.655 + vertex 0 -8.46107 14.655 + endloop + endfacet + facet normal 0 -0.913546 0.406736 + outer loop + vertex 200 -8.46107 14.655 + vertex 0 -9.29182 12.7891 + vertex 200 -9.29182 12.7891 + endloop + endfacet + facet normal 0 -0.978148 0.207909 + outer loop + vertex 0 -9.71648 10.7912 + vertex 200 -9.29182 12.7891 + vertex 0 -9.29182 12.7891 + endloop + endfacet + facet normal 0 -0.978148 0.207909 + outer loop + vertex 200 -9.29182 12.7891 + vertex 0 -9.71648 10.7912 + vertex 200 -9.71648 10.7912 + endloop + endfacet + facet normal 0 -1 0 + outer loop + vertex 0 -9.71648 8.74876 + vertex 200 -9.71648 10.7912 + vertex 0 -9.71648 10.7912 + endloop + endfacet + facet normal 0 -1 -0 + outer loop + vertex 200 -9.71648 10.7912 + vertex 0 -9.71648 8.74876 + vertex 200 -9.71648 8.74876 + endloop + endfacet + facet normal 0 -0.978147 -0.207913 + outer loop + vertex 0 -9.29182 6.7509 + vertex 200 -9.71648 8.74876 + vertex 0 -9.71648 8.74876 + endloop + endfacet + facet normal 0 -0.978147 -0.207913 + outer loop + vertex 200 -9.71648 8.74876 + vertex 0 -9.29182 6.7509 + vertex 200 -9.29182 6.7509 + endloop + endfacet + facet normal 0 -0.913546 -0.406736 + outer loop + vertex 0 -8.46107 4.885 + vertex 200 -9.29182 6.7509 + vertex 0 -9.29182 6.7509 + endloop + endfacet + facet normal 0 -0.913546 -0.406736 + outer loop + vertex 200 -9.29182 6.7509 + vertex 0 -8.46107 4.885 + vertex 200 -8.46107 4.885 + endloop + endfacet + facet normal 0 -0.809016 -0.587786 + outer loop + vertex 0 -7.26052 3.23259 + vertex 200 -8.46107 4.885 + vertex 0 -8.46107 4.885 + endloop + endfacet + facet normal 0 -0.809016 -0.587786 + outer loop + vertex 200 -8.46107 4.885 + vertex 0 -7.26052 3.23259 + vertex 200 -7.26052 3.23259 + endloop + endfacet + facet normal 0 -0.669131 -0.743144 + outer loop + vertex 0 -7.26052 3.23259 + vertex 200 -5.74266 1.8659 + vertex 200 -7.26052 3.23259 + endloop + endfacet + facet normal -0 -0.669131 -0.743144 + outer loop + vertex 200 -5.74266 1.8659 + vertex 0 -7.26052 3.23259 + vertex 0 -5.74266 1.8659 + endloop + endfacet + facet normal 0 -0.5 -0.866026 + outer loop + vertex 0 -5.74266 1.8659 + vertex 200 -3.97382 0.844661 + vertex 200 -5.74266 1.8659 + endloop + endfacet + facet normal -0 -0.5 -0.866026 + outer loop + vertex 200 -3.97382 0.844661 + vertex 0 -5.74266 1.8659 + vertex 0 -3.97382 0.844661 + endloop + endfacet + facet normal 0 -0.309017 -0.951057 + outer loop + vertex 0 -3.97382 0.844661 + vertex 200 -2.0313 0.213498 + vertex 200 -3.97382 0.844661 + endloop + endfacet + facet normal -0 -0.309017 -0.951057 + outer loop + vertex 200 -2.0313 0.213498 + vertex 0 -3.97382 0.844661 + vertex 0 -2.0313 0.213498 + endloop + endfacet + facet normal 0 -0.104528 -0.994522 + outer loop + vertex 0 -2.0313 0.213498 + vertex 200 0 0 + vertex 200 -2.0313 0.213498 + endloop + endfacet + facet normal -0 -0.104528 -0.994522 + outer loop + vertex 200 0 0 + vertex 0 -2.0313 0.213498 + vertex 0 0 0 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 9.29182 6.7509 + vertex 0 9.71648 10.7912 + vertex 0 9.71648 8.74876 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 9.29182 6.7509 + vertex 0 9.29182 12.7891 + vertex 0 9.71648 10.7912 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 8.46107 4.885 + vertex 0 9.29182 12.7891 + vertex 0 9.29182 6.7509 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 8.46107 4.885 + vertex 0 8.46107 14.655 + vertex 0 9.29182 12.7891 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 7.26052 3.23259 + vertex 0 8.46107 14.655 + vertex 0 8.46107 4.885 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 7.26052 3.23259 + vertex 0 7.26052 16.3074 + vertex 0 8.46107 14.655 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 5.74266 1.8659 + vertex 0 7.26052 16.3074 + vertex 0 7.26052 3.23259 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 5.74266 1.8659 + vertex 0 5.74266 17.6741 + vertex 0 7.26052 16.3074 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 3.97382 0.844661 + vertex 0 5.74266 17.6741 + vertex 0 5.74266 1.8659 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 3.97382 0.844661 + vertex 0 3.97382 18.6953 + vertex 0 5.74266 17.6741 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 2.0313 0.213498 + vertex 0 3.97382 18.6953 + vertex 0 3.97382 0.844661 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 2.0313 0.213498 + vertex 0 2.0313 19.3265 + vertex 0 3.97382 18.6953 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 0 + vertex 0 2.0313 19.3265 + vertex 0 2.0313 0.213498 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 0 + vertex 0 0 19.54 + vertex 0 2.0313 19.3265 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -2.0313 0.213498 + vertex 0 0 19.54 + vertex 0 0 0 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -2.0313 0.213498 + vertex 0 -2.0313 19.3265 + vertex 0 0 19.54 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -3.97382 0.844661 + vertex 0 -2.0313 19.3265 + vertex 0 -2.0313 0.213498 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -3.97382 0.844661 + vertex 0 -3.97382 18.6953 + vertex 0 -2.0313 19.3265 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -5.74266 1.8659 + vertex 0 -3.97382 18.6953 + vertex 0 -3.97382 0.844661 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -5.74266 1.8659 + vertex 0 -5.74266 17.6741 + vertex 0 -3.97382 18.6953 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -7.26052 3.23259 + vertex 0 -5.74266 17.6741 + vertex 0 -5.74266 1.8659 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -7.26052 3.23259 + vertex 0 -7.26052 16.3074 + vertex 0 -5.74266 17.6741 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -8.46107 4.885 + vertex 0 -7.26052 16.3074 + vertex 0 -7.26052 3.23259 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -8.46107 4.885 + vertex 0 -8.46107 14.655 + vertex 0 -7.26052 16.3074 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -9.29182 6.7509 + vertex 0 -8.46107 14.655 + vertex 0 -8.46107 4.885 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -9.29182 6.7509 + vertex 0 -9.29182 12.7891 + vertex 0 -8.46107 14.655 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -9.71648 8.74876 + vertex 0 -9.29182 12.7891 + vertex 0 -9.29182 6.7509 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -9.29182 12.7891 + vertex 0 -9.71648 8.74876 + vertex 0 -9.71648 10.7912 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 9.29182 12.7891 + vertex 200 9.71648 8.74876 + vertex 200 9.71648 10.7912 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 9.29182 12.7891 + vertex 200 9.29182 6.7509 + vertex 200 9.71648 8.74876 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 8.46107 14.655 + vertex 200 9.29182 6.7509 + vertex 200 9.29182 12.7891 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 8.46107 14.655 + vertex 200 8.46107 4.885 + vertex 200 9.29182 6.7509 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 7.26052 16.3074 + vertex 200 8.46107 4.885 + vertex 200 8.46107 14.655 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 7.26052 16.3074 + vertex 200 7.26052 3.23259 + vertex 200 8.46107 4.885 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 5.74266 17.6741 + vertex 200 7.26052 3.23259 + vertex 200 7.26052 16.3074 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 5.74266 17.6741 + vertex 200 5.74266 1.8659 + vertex 200 7.26052 3.23259 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 3.97382 18.6953 + vertex 200 5.74266 1.8659 + vertex 200 5.74266 17.6741 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 3.97382 18.6953 + vertex 200 3.97382 0.844661 + vertex 200 5.74266 1.8659 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 2.0313 19.3265 + vertex 200 3.97382 0.844661 + vertex 200 3.97382 18.6953 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 2.0313 19.3265 + vertex 200 2.0313 0.213498 + vertex 200 3.97382 0.844661 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 0 19.54 + vertex 200 2.0313 0.213498 + vertex 200 2.0313 19.3265 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 0 19.54 + vertex 200 0 0 + vertex 200 2.0313 0.213498 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -2.0313 19.3265 + vertex 200 0 0 + vertex 200 0 19.54 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -2.0313 19.3265 + vertex 200 -2.0313 0.213498 + vertex 200 0 0 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -3.97382 18.6953 + vertex 200 -2.0313 0.213498 + vertex 200 -2.0313 19.3265 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -3.97382 18.6953 + vertex 200 -3.97382 0.844661 + vertex 200 -2.0313 0.213498 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -5.74266 17.6741 + vertex 200 -3.97382 0.844661 + vertex 200 -3.97382 18.6953 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -5.74266 17.6741 + vertex 200 -5.74266 1.8659 + vertex 200 -3.97382 0.844661 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -7.26052 16.3074 + vertex 200 -5.74266 1.8659 + vertex 200 -5.74266 17.6741 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -7.26052 16.3074 + vertex 200 -7.26052 3.23259 + vertex 200 -5.74266 1.8659 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -8.46107 14.655 + vertex 200 -7.26052 3.23259 + vertex 200 -7.26052 16.3074 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -8.46107 14.655 + vertex 200 -8.46107 4.885 + vertex 200 -7.26052 3.23259 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -9.29182 12.7891 + vertex 200 -8.46107 4.885 + vertex 200 -8.46107 14.655 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -9.29182 12.7891 + vertex 200 -9.29182 6.7509 + vertex 200 -8.46107 4.885 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -9.71648 10.7912 + vertex 200 -9.29182 6.7509 + vertex 200 -9.29182 12.7891 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -9.29182 6.7509 + vertex 200 -9.71648 10.7912 + vertex 200 -9.71648 8.74876 + endloop + endfacet +endsolid OpenSCAD_Model diff --git a/files/benchmark-pipe.scad b/files/benchmark-pipe.scad new file mode 100644 index 000000000..c109a5d9b --- /dev/null +++ b/files/benchmark-pipe.scad @@ -0,0 +1,12 @@ +l = 200; +d = 27.87; +t = 4; + +translate([ 0, 0, d / 2 ]) rotate([ 0, 90, 0 ]) difference() +{ + // Outer cylinder + cylinder(d = d, h = l); + + // Inner cylinder + cylinder(d = d - 2 * t, h = l); +} \ No newline at end of file diff --git a/files/benchmark-pipe.stl b/files/benchmark-pipe.stl new file mode 100644 index 000000000..a884a9324 --- /dev/null +++ b/files/benchmark-pipe.stl @@ -0,0 +1,1682 @@ +solid OpenSCAD_Model + facet normal -0 0.309008 0.951059 + outer loop + vertex 0 5.66787 26.6653 + vertex 200 2.89725 27.5655 + vertex 200 5.66787 26.6653 + endloop + endfacet + facet normal 0 0.309008 0.951059 + outer loop + vertex 200 2.89725 27.5655 + vertex 0 5.66787 26.6653 + vertex 0 2.89725 27.5655 + endloop + endfacet + facet normal 0 0.309016 -0.951057 + outer loop + vertex 0 2.89725 0.304514 + vertex 200 5.66787 1.20474 + vertex 200 2.89725 0.304514 + endloop + endfacet + facet normal 0 0.309016 -0.951057 + outer loop + vertex 200 5.66787 1.20474 + vertex 0 2.89725 0.304514 + vertex 0 5.66787 1.20474 + endloop + endfacet + facet normal 0 0.104529 -0.994522 + outer loop + vertex 0 0 5.34058e-07 + vertex 200 2.89725 0.304514 + vertex 200 0 5.34058e-07 + endloop + endfacet + facet normal 0 0.104529 -0.994522 + outer loop + vertex 200 2.89725 0.304514 + vertex 0 0 5.34058e-07 + vertex 0 2.89725 0.304514 + endloop + endfacet + facet normal 0 0.5 -0.866025 + outer loop + vertex 0 5.66787 1.20474 + vertex 200 8.19079 2.66135 + vertex 200 5.66787 1.20474 + endloop + endfacet + facet normal 0 0.5 -0.866025 + outer loop + vertex 200 8.19079 2.66135 + vertex 0 5.66787 1.20474 + vertex 0 8.19079 2.66135 + endloop + endfacet + facet normal 0 -0.669136 -0.74314 + outer loop + vertex 0 -10.3557 4.61067 + vertex 200 -8.19079 2.66135 + vertex 200 -10.3557 4.61067 + endloop + endfacet + facet normal -0 -0.669136 -0.74314 + outer loop + vertex 200 -8.19079 2.66135 + vertex 0 -10.3557 4.61067 + vertex 0 -8.19079 2.66135 + endloop + endfacet + facet normal 0 0.913547 -0.406734 + outer loop + vertex 200 12.0681 6.9675 + vertex 0 13.253 9.62885 + vertex 200 13.253 9.62885 + endloop + endfacet + facet normal 0 0.913547 -0.406734 + outer loop + vertex 0 13.253 9.62885 + vertex 200 12.0681 6.9675 + vertex 0 12.0681 6.9675 + endloop + endfacet + facet normal 0 -0.104529 -0.994522 + outer loop + vertex 0 -2.89725 0.304514 + vertex 200 0 5.34058e-07 + vertex 200 -2.89725 0.304514 + endloop + endfacet + facet normal -0 -0.104529 -0.994522 + outer loop + vertex 200 0 5.34058e-07 + vertex 0 -2.89725 0.304514 + vertex 0 0 5.34058e-07 + endloop + endfacet + facet normal 0 -0.5 -0.866025 + outer loop + vertex 0 -8.19079 2.66135 + vertex 200 -5.66787 1.20474 + vertex 200 -8.19079 2.66135 + endloop + endfacet + facet normal -0 -0.5 -0.866025 + outer loop + vertex 200 -5.66787 1.20474 + vertex 0 -8.19079 2.66135 + vertex 0 -5.66787 1.20474 + endloop + endfacet + facet normal -0 0.499998 0.866027 + outer loop + vertex 0 8.19079 25.2087 + vertex 200 5.66787 26.6653 + vertex 200 8.19079 25.2087 + endloop + endfacet + facet normal 0 0.499998 0.866027 + outer loop + vertex 200 5.66787 26.6653 + vertex 0 8.19079 25.2087 + vertex 0 5.66787 26.6653 + endloop + endfacet + facet normal 0 0.669136 -0.74314 + outer loop + vertex 0 8.19079 2.66135 + vertex 200 10.3557 4.61067 + vertex 200 8.19079 2.66135 + endloop + endfacet + facet normal 0 0.669136 -0.74314 + outer loop + vertex 200 10.3557 4.61067 + vertex 0 8.19079 2.66135 + vertex 0 10.3557 4.61067 + endloop + endfacet + facet normal 0 -0.499998 0.866027 + outer loop + vertex 0 -5.66787 26.6653 + vertex 200 -8.19079 25.2087 + vertex 200 -5.66787 26.6653 + endloop + endfacet + facet normal 0 -0.499998 0.866027 + outer loop + vertex 200 -8.19079 25.2087 + vertex 0 -5.66787 26.6653 + vertex 0 -8.19079 25.2087 + endloop + endfacet + facet normal 0 -0.913547 -0.406734 + outer loop + vertex 0 -12.0681 6.9675 + vertex 200 -13.253 9.62885 + vertex 0 -13.253 9.62885 + endloop + endfacet + facet normal 0 -0.913547 -0.406734 + outer loop + vertex 200 -13.253 9.62885 + vertex 0 -12.0681 6.9675 + vertex 200 -12.0681 6.9675 + endloop + endfacet + facet normal 0 -0.978148 0.207911 + outer loop + vertex 0 -13.8587 15.3916 + vertex 200 -13.253 18.2412 + vertex 0 -13.253 18.2412 + endloop + endfacet + facet normal 0 -0.978148 0.207911 + outer loop + vertex 200 -13.253 18.2412 + vertex 0 -13.8587 15.3916 + vertex 200 -13.8587 15.3916 + endloop + endfacet + facet normal -0 0.104524 0.994522 + outer loop + vertex 0 2.89725 27.5655 + vertex 200 0 27.87 + vertex 200 2.89725 27.5655 + endloop + endfacet + facet normal 0 0.104524 0.994522 + outer loop + vertex 200 0 27.87 + vertex 0 2.89725 27.5655 + vertex 0 0 27.87 + endloop + endfacet + facet normal 0 -0.309016 -0.951057 + outer loop + vertex 0 -5.66787 1.20474 + vertex 200 -2.89725 0.304514 + vertex 200 -5.66787 1.20474 + endloop + endfacet + facet normal -0 -0.309016 -0.951057 + outer loop + vertex 200 -2.89725 0.304514 + vertex 0 -5.66787 1.20474 + vertex 0 -2.89725 0.304514 + endloop + endfacet + facet normal 0 1 -0 + outer loop + vertex 200 13.8587 12.4784 + vertex 0 13.8587 15.3916 + vertex 200 13.8587 15.3916 + endloop + endfacet + facet normal 0 1 0 + outer loop + vertex 0 13.8587 15.3916 + vertex 200 13.8587 12.4784 + vertex 0 13.8587 12.4784 + endloop + endfacet + facet normal 0 0.809003 0.587804 + outer loop + vertex 200 12.0681 20.9025 + vertex 0 10.3557 23.2593 + vertex 200 10.3557 23.2593 + endloop + endfacet + facet normal 0 0.809003 0.587804 + outer loop + vertex 0 10.3557 23.2593 + vertex 200 12.0681 20.9025 + vertex 0 12.0681 20.9025 + endloop + endfacet + facet normal 0 -0.913544 0.40674 + outer loop + vertex 0 -13.253 18.2412 + vertex 200 -12.0681 20.9025 + vertex 0 -12.0681 20.9025 + endloop + endfacet + facet normal 0 -0.913544 0.40674 + outer loop + vertex 200 -12.0681 20.9025 + vertex 0 -13.253 18.2412 + vertex 200 -13.253 18.2412 + endloop + endfacet + facet normal 0 -0.809007 -0.587799 + outer loop + vertex 0 -10.3557 4.61067 + vertex 200 -12.0681 6.9675 + vertex 0 -12.0681 6.9675 + endloop + endfacet + facet normal 0 -0.809007 -0.587799 + outer loop + vertex 200 -12.0681 6.9675 + vertex 0 -10.3557 4.61067 + vertex 200 -10.3557 4.61067 + endloop + endfacet + facet normal -0 0.669151 0.743127 + outer loop + vertex 0 10.3557 23.2593 + vertex 200 8.19079 25.2087 + vertex 200 10.3557 23.2593 + endloop + endfacet + facet normal 0 0.669151 0.743127 + outer loop + vertex 200 8.19079 25.2087 + vertex 0 10.3557 23.2593 + vertex 0 8.19079 25.2087 + endloop + endfacet + facet normal 0 0.913544 0.40674 + outer loop + vertex 200 13.253 18.2412 + vertex 0 12.0681 20.9025 + vertex 200 12.0681 20.9025 + endloop + endfacet + facet normal 0 0.913544 0.40674 + outer loop + vertex 0 12.0681 20.9025 + vertex 200 13.253 18.2412 + vertex 0 13.253 18.2412 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 9.88057 14.9735 + vertex 200 13.8587 15.3916 + vertex 200 13.253 18.2412 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 13.8587 15.3916 + vertex 200 9.88057 14.9735 + vertex 200 13.8587 12.4784 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 9.44875 17.0051 + vertex 200 13.253 18.2412 + vertex 200 12.0681 20.9025 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 9.88057 12.8965 + vertex 200 13.8587 12.4784 + vertex 200 9.88057 14.9735 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 8.60396 18.9025 + vertex 200 12.0681 20.9025 + vertex 200 10.3557 23.2593 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 13.8587 12.4784 + vertex 200 9.88057 12.8965 + vertex 200 13.253 9.62885 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 9.44875 10.8649 + vertex 200 13.253 9.62885 + vertex 200 9.88057 12.8965 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 13.253 18.2412 + vertex 200 9.44875 17.0051 + vertex 200 9.88057 14.9735 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 12.0681 20.9025 + vertex 200 8.60396 18.9025 + vertex 200 9.44875 17.0051 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 7.38314 20.5828 + vertex 200 10.3557 23.2593 + vertex 200 8.19079 25.2087 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 10.3557 23.2593 + vertex 200 7.38314 20.5828 + vertex 200 8.60396 18.9025 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 8.19079 25.2087 + vertex 200 5.83965 21.9726 + vertex 200 7.38314 20.5828 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 5.66787 26.6653 + vertex 200 5.83965 21.9726 + vertex 200 8.19079 25.2087 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 5.66787 26.6653 + vertex 200 4.04093 23.0111 + vertex 200 5.83965 21.9726 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 2.89725 27.5655 + vertex 200 4.04093 23.0111 + vertex 200 5.66787 26.6653 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 2.89725 27.5655 + vertex 200 2.0656 23.6529 + vertex 200 4.04093 23.0111 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 0 27.87 + vertex 200 2.0656 23.6529 + vertex 200 2.89725 27.5655 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 0 27.87 + vertex 200 0 23.87 + vertex 200 2.0656 23.6529 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 0 27.87 + vertex 200 -2.0656 23.6529 + vertex 200 0 23.87 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -2.89725 27.5655 + vertex 200 -2.0656 23.6529 + vertex 200 0 27.87 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -2.89725 27.5655 + vertex 200 -4.04093 23.0111 + vertex 200 -2.0656 23.6529 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -5.66787 26.6653 + vertex 200 -4.04093 23.0111 + vertex 200 -2.89725 27.5655 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -5.66787 26.6653 + vertex 200 -5.83965 21.9726 + vertex 200 -4.04093 23.0111 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -8.19079 25.2087 + vertex 200 -5.83965 21.9726 + vertex 200 -5.66787 26.6653 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -5.83965 21.9726 + vertex 200 -8.19079 25.2087 + vertex 200 -7.38314 20.5828 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -10.3557 23.2593 + vertex 200 -7.38314 20.5828 + vertex 200 -8.19079 25.2087 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -7.38314 20.5828 + vertex 200 -10.3557 23.2593 + vertex 200 -8.60396 18.9025 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -12.0681 20.9025 + vertex 200 -8.60396 18.9025 + vertex 200 -10.3557 23.2593 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -13.253 18.2412 + vertex 200 -9.44875 17.0051 + vertex 200 -12.0681 20.9025 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -8.60396 18.9025 + vertex 200 -12.0681 20.9025 + vertex 200 -9.44875 17.0051 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 13.253 9.62885 + vertex 200 9.44875 10.8649 + vertex 200 12.0681 6.9675 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 8.60396 8.9675 + vertex 200 12.0681 6.9675 + vertex 200 9.44875 10.8649 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 12.0681 6.9675 + vertex 200 8.60396 8.9675 + vertex 200 10.3557 4.61067 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 7.38314 7.28719 + vertex 200 10.3557 4.61067 + vertex 200 8.60396 8.9675 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 10.3557 4.61067 + vertex 200 7.38314 7.28719 + vertex 200 8.19079 2.66135 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 5.83965 5.89742 + vertex 200 8.19079 2.66135 + vertex 200 7.38314 7.28719 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 5.83965 5.89742 + vertex 200 5.66787 1.20474 + vertex 200 8.19079 2.66135 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 4.04093 4.85893 + vertex 200 5.66787 1.20474 + vertex 200 5.83965 5.89742 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 4.04093 4.85893 + vertex 200 2.89725 0.304514 + vertex 200 5.66787 1.20474 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 2.0656 4.2171 + vertex 200 2.89725 0.304514 + vertex 200 4.04093 4.85893 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 0 4 + vertex 200 2.89725 0.304514 + vertex 200 2.0656 4.2171 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 0 4 + vertex 200 0 5.34058e-07 + vertex 200 2.89725 0.304514 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -2.0656 4.2171 + vertex 200 0 5.34058e-07 + vertex 200 0 4 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -2.0656 4.2171 + vertex 200 -2.89725 0.304514 + vertex 200 0 5.34058e-07 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -4.04093 4.85893 + vertex 200 -2.89725 0.304514 + vertex 200 -2.0656 4.2171 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -4.04093 4.85893 + vertex 200 -5.66787 1.20474 + vertex 200 -2.89725 0.304514 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -5.83965 5.89742 + vertex 200 -5.66787 1.20474 + vertex 200 -4.04093 4.85893 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -8.19079 2.66135 + vertex 200 -5.83965 5.89742 + vertex 200 -7.38314 7.28719 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -5.83965 5.89742 + vertex 200 -8.19079 2.66135 + vertex 200 -5.66787 1.20474 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -10.3557 4.61067 + vertex 200 -7.38314 7.28719 + vertex 200 -8.60396 8.9675 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -12.0681 6.9675 + vertex 200 -8.60396 8.9675 + vertex 200 -9.44875 10.8649 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -13.253 9.62885 + vertex 200 -9.44875 10.8649 + vertex 200 -9.88057 12.8965 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -7.38314 7.28719 + vertex 200 -10.3557 4.61067 + vertex 200 -8.19079 2.66135 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -9.44875 17.0051 + vertex 200 -13.253 18.2412 + vertex 200 -9.88057 14.9735 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 -13.8587 15.3916 + vertex 200 -9.88057 14.9735 + vertex 200 -13.253 18.2412 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -8.60396 8.9675 + vertex 200 -12.0681 6.9675 + vertex 200 -10.3557 4.61067 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -9.88057 14.9735 + vertex 200 -13.8587 15.3916 + vertex 200 -9.88057 12.8965 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -9.44875 10.8649 + vertex 200 -13.253 9.62885 + vertex 200 -12.0681 6.9675 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -13.8587 12.4784 + vertex 200 -9.88057 12.8965 + vertex 200 -13.8587 15.3916 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 -9.88057 12.8965 + vertex 200 -13.8587 12.4784 + vertex 200 -13.253 9.62885 + endloop + endfacet + facet normal 0 0.978147 -0.207915 + outer loop + vertex 200 13.253 9.62885 + vertex 0 13.8587 12.4784 + vertex 200 13.8587 12.4784 + endloop + endfacet + facet normal 0 0.978147 -0.207915 + outer loop + vertex 0 13.8587 12.4784 + vertex 200 13.253 9.62885 + vertex 0 13.253 9.62885 + endloop + endfacet + facet normal 0 0.809007 -0.587799 + outer loop + vertex 200 10.3557 4.61067 + vertex 0 12.0681 6.9675 + vertex 200 12.0681 6.9675 + endloop + endfacet + facet normal 0 0.809007 -0.587799 + outer loop + vertex 0 12.0681 6.9675 + vertex 200 10.3557 4.61067 + vertex 0 10.3557 4.61067 + endloop + endfacet + facet normal 0 0.978148 0.207911 + outer loop + vertex 200 13.8587 15.3916 + vertex 0 13.253 18.2412 + vertex 200 13.253 18.2412 + endloop + endfacet + facet normal 0 0.978148 0.207911 + outer loop + vertex 0 13.253 18.2412 + vertex 200 13.8587 15.3916 + vertex 0 13.8587 15.3916 + endloop + endfacet + facet normal 0 -0.309008 0.951059 + outer loop + vertex 0 -2.89725 27.5655 + vertex 200 -5.66787 26.6653 + vertex 200 -2.89725 27.5655 + endloop + endfacet + facet normal 0 -0.309008 0.951059 + outer loop + vertex 200 -5.66787 26.6653 + vertex 0 -2.89725 27.5655 + vertex 0 -5.66787 26.6653 + endloop + endfacet + facet normal 0 -0.104524 0.994522 + outer loop + vertex 0 0 27.87 + vertex 200 -2.89725 27.5655 + vertex 200 0 27.87 + endloop + endfacet + facet normal 0 -0.104524 0.994522 + outer loop + vertex 200 -2.89725 27.5655 + vertex 0 0 27.87 + vertex 0 -2.89725 27.5655 + endloop + endfacet + facet normal 0 -0.978147 -0.207915 + outer loop + vertex 0 -13.253 9.62885 + vertex 200 -13.8587 12.4784 + vertex 0 -13.8587 12.4784 + endloop + endfacet + facet normal 0 -0.978147 -0.207915 + outer loop + vertex 200 -13.8587 12.4784 + vertex 0 -13.253 9.62885 + vertex 200 -13.253 9.62885 + endloop + endfacet + facet normal 0 -1 0 + outer loop + vertex 0 -13.8587 12.4784 + vertex 200 -13.8587 15.3916 + vertex 0 -13.8587 15.3916 + endloop + endfacet + facet normal 0 -1 -0 + outer loop + vertex 200 -13.8587 15.3916 + vertex 0 -13.8587 12.4784 + vertex 200 -13.8587 12.4784 + endloop + endfacet + facet normal 0 -0.669151 0.743127 + outer loop + vertex 0 -8.19079 25.2087 + vertex 200 -10.3557 23.2593 + vertex 200 -8.19079 25.2087 + endloop + endfacet + facet normal 0 -0.669151 0.743127 + outer loop + vertex 200 -10.3557 23.2593 + vertex 0 -8.19079 25.2087 + vertex 0 -10.3557 23.2593 + endloop + endfacet + facet normal 0 -0.809003 0.587804 + outer loop + vertex 0 -12.0681 20.9025 + vertex 200 -10.3557 23.2593 + vertex 0 -10.3557 23.2593 + endloop + endfacet + facet normal 0 -0.809003 0.587804 + outer loop + vertex 200 -10.3557 23.2593 + vertex 0 -12.0681 20.9025 + vertex 200 -12.0681 20.9025 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 9.88057 12.8965 + vertex 0 13.8587 12.4784 + vertex 0 13.253 9.62885 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 13.8587 12.4784 + vertex 0 9.88057 12.8965 + vertex 0 13.8587 15.3916 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 9.44875 10.8649 + vertex 0 13.253 9.62885 + vertex 0 12.0681 6.9675 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 9.88057 14.9735 + vertex 0 13.8587 15.3916 + vertex 0 9.88057 12.8965 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 8.60396 8.9675 + vertex 0 12.0681 6.9675 + vertex 0 10.3557 4.61067 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 13.8587 15.3916 + vertex 0 9.88057 14.9735 + vertex 0 13.253 18.2412 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 9.44875 17.0051 + vertex 0 13.253 18.2412 + vertex 0 9.88057 14.9735 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 13.253 9.62885 + vertex 0 9.44875 10.8649 + vertex 0 9.88057 12.8965 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 12.0681 6.9675 + vertex 0 8.60396 8.9675 + vertex 0 9.44875 10.8649 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 7.38314 7.28719 + vertex 0 10.3557 4.61067 + vertex 0 8.19079 2.66135 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 10.3557 4.61067 + vertex 0 7.38314 7.28719 + vertex 0 8.60396 8.9675 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 8.19079 2.66135 + vertex 0 5.83965 5.89742 + vertex 0 7.38314 7.28719 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 5.66787 1.20474 + vertex 0 5.83965 5.89742 + vertex 0 8.19079 2.66135 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 5.66787 1.20474 + vertex 0 4.04093 4.85893 + vertex 0 5.83965 5.89742 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 2.89725 0.304514 + vertex 0 4.04093 4.85893 + vertex 0 5.66787 1.20474 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 2.89725 0.304514 + vertex 0 2.0656 4.2171 + vertex 0 4.04093 4.85893 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 5.34058e-07 + vertex 0 2.0656 4.2171 + vertex 0 2.89725 0.304514 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 5.34058e-07 + vertex 0 0 4 + vertex 0 2.0656 4.2171 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 5.34058e-07 + vertex 0 -2.0656 4.2171 + vertex 0 0 4 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -2.89725 0.304514 + vertex 0 -2.0656 4.2171 + vertex 0 0 5.34058e-07 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -2.89725 0.304514 + vertex 0 -4.04093 4.85893 + vertex 0 -2.0656 4.2171 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -5.66787 1.20474 + vertex 0 -4.04093 4.85893 + vertex 0 -2.89725 0.304514 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -5.66787 1.20474 + vertex 0 -5.83965 5.89742 + vertex 0 -4.04093 4.85893 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -8.19079 2.66135 + vertex 0 -5.83965 5.89742 + vertex 0 -5.66787 1.20474 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 -5.83965 5.89742 + vertex 0 -8.19079 2.66135 + vertex 0 -7.38314 7.28719 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -10.3557 4.61067 + vertex 0 -7.38314 7.28719 + vertex 0 -8.19079 2.66135 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 -7.38314 7.28719 + vertex 0 -10.3557 4.61067 + vertex 0 -8.60396 8.9675 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -12.0681 6.9675 + vertex 0 -8.60396 8.9675 + vertex 0 -10.3557 4.61067 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -13.253 9.62885 + vertex 0 -9.44875 10.8649 + vertex 0 -12.0681 6.9675 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 -8.60396 8.9675 + vertex 0 -12.0681 6.9675 + vertex 0 -9.44875 10.8649 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 13.253 18.2412 + vertex 0 9.44875 17.0051 + vertex 0 12.0681 20.9025 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 8.60396 18.9025 + vertex 0 12.0681 20.9025 + vertex 0 9.44875 17.0051 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 12.0681 20.9025 + vertex 0 8.60396 18.9025 + vertex 0 10.3557 23.2593 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 7.38314 20.5828 + vertex 0 10.3557 23.2593 + vertex 0 8.60396 18.9025 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 10.3557 23.2593 + vertex 0 7.38314 20.5828 + vertex 0 8.19079 25.2087 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 5.83965 21.9726 + vertex 0 8.19079 25.2087 + vertex 0 7.38314 20.5828 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 5.83965 21.9726 + vertex 0 5.66787 26.6653 + vertex 0 8.19079 25.2087 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 4.04093 23.0111 + vertex 0 5.66787 26.6653 + vertex 0 5.83965 21.9726 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 4.04093 23.0111 + vertex 0 2.89725 27.5655 + vertex 0 5.66787 26.6653 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 2.0656 23.6529 + vertex 0 2.89725 27.5655 + vertex 0 4.04093 23.0111 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 23.87 + vertex 0 2.89725 27.5655 + vertex 0 2.0656 23.6529 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 23.87 + vertex 0 0 27.87 + vertex 0 2.89725 27.5655 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -2.0656 23.6529 + vertex 0 0 27.87 + vertex 0 0 23.87 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -2.0656 23.6529 + vertex 0 -2.89725 27.5655 + vertex 0 0 27.87 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -4.04093 23.0111 + vertex 0 -2.89725 27.5655 + vertex 0 -2.0656 23.6529 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -4.04093 23.0111 + vertex 0 -5.66787 26.6653 + vertex 0 -2.89725 27.5655 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -5.83965 21.9726 + vertex 0 -5.66787 26.6653 + vertex 0 -4.04093 23.0111 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -8.19079 25.2087 + vertex 0 -5.83965 21.9726 + vertex 0 -7.38314 20.5828 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -5.83965 21.9726 + vertex 0 -8.19079 25.2087 + vertex 0 -5.66787 26.6653 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -10.3557 23.2593 + vertex 0 -7.38314 20.5828 + vertex 0 -8.60396 18.9025 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -12.0681 20.9025 + vertex 0 -8.60396 18.9025 + vertex 0 -9.44875 17.0051 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -13.253 18.2412 + vertex 0 -9.44875 17.0051 + vertex 0 -9.88057 14.9735 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -7.38314 20.5828 + vertex 0 -10.3557 23.2593 + vertex 0 -8.19079 25.2087 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 -9.44875 10.8649 + vertex 0 -13.253 9.62885 + vertex 0 -9.88057 12.8965 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -13.8587 12.4784 + vertex 0 -9.88057 12.8965 + vertex 0 -13.253 9.62885 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -8.60396 18.9025 + vertex 0 -12.0681 20.9025 + vertex 0 -10.3557 23.2593 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 -9.88057 12.8965 + vertex 0 -13.8587 12.4784 + vertex 0 -9.88057 14.9735 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -9.44875 17.0051 + vertex 0 -13.253 18.2412 + vertex 0 -12.0681 20.9025 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -13.8587 15.3916 + vertex 0 -9.88057 14.9735 + vertex 0 -13.8587 12.4784 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 -9.88057 14.9735 + vertex 0 -13.8587 15.3916 + vertex 0 -13.253 18.2412 + endloop + endfacet + facet normal 0 0.66914 -0.743136 + outer loop + vertex 0 -7.38314 20.5828 + vertex 200 -5.83965 21.9726 + vertex 200 -7.38314 20.5828 + endloop + endfacet + facet normal 0 0.66914 -0.743136 + outer loop + vertex 200 -5.83965 21.9726 + vertex 0 -7.38314 20.5828 + vertex 0 -5.83965 21.9726 + endloop + endfacet + facet normal 0 -0.30902 0.951056 + outer loop + vertex 0 4.04093 4.85893 + vertex 200 2.0656 4.2171 + vertex 200 4.04093 4.85893 + endloop + endfacet + facet normal 0 -0.30902 0.951056 + outer loop + vertex 200 2.0656 4.2171 + vertex 0 4.04093 4.85893 + vertex 0 2.0656 4.2171 + endloop + endfacet + facet normal 0 -0.104527 0.994522 + outer loop + vertex 0 2.0656 4.2171 + vertex 200 0 4 + vertex 200 2.0656 4.2171 + endloop + endfacet + facet normal 0 -0.104527 0.994522 + outer loop + vertex 200 0 4 + vertex 0 2.0656 4.2171 + vertex 0 0 4 + endloop + endfacet + facet normal 0 -0.499999 0.866026 + outer loop + vertex 0 5.83965 5.89742 + vertex 200 4.04093 4.85893 + vertex 200 5.83965 5.89742 + endloop + endfacet + facet normal 0 -0.499999 0.866026 + outer loop + vertex 200 4.04093 4.85893 + vertex 0 5.83965 5.89742 + vertex 0 4.04093 4.85893 + endloop + endfacet + facet normal -0 0.30902 0.951056 + outer loop + vertex 0 -2.0656 4.2171 + vertex 200 -4.04093 4.85893 + vertex 200 -2.0656 4.2171 + endloop + endfacet + facet normal 0 0.30902 0.951056 + outer loop + vertex 200 -4.04093 4.85893 + vertex 0 -2.0656 4.2171 + vertex 0 -4.04093 4.85893 + endloop + endfacet + facet normal 0 0.809016 0.587786 + outer loop + vertex 200 -7.38314 7.28719 + vertex 0 -8.60396 8.9675 + vertex 200 -8.60396 8.9675 + endloop + endfacet + facet normal 0 0.809016 0.587786 + outer loop + vertex 0 -8.60396 8.9675 + vertex 200 -7.38314 7.28719 + vertex 0 -7.38314 7.28719 + endloop + endfacet + facet normal -0 0.669132 0.743144 + outer loop + vertex 0 -5.83965 5.89742 + vertex 200 -7.38314 7.28719 + vertex 200 -5.83965 5.89742 + endloop + endfacet + facet normal 0 0.669132 0.743144 + outer loop + vertex 200 -7.38314 7.28719 + vertex 0 -5.83965 5.89742 + vertex 0 -7.38314 7.28719 + endloop + endfacet + facet normal -0 0.104527 0.994522 + outer loop + vertex 0 0 4 + vertex 200 -2.0656 4.2171 + vertex 200 0 4 + endloop + endfacet + facet normal 0 0.104527 0.994522 + outer loop + vertex 200 -2.0656 4.2171 + vertex 0 0 4 + vertex 0 -2.0656 4.2171 + endloop + endfacet + facet normal 0 0.913543 -0.406742 + outer loop + vertex 200 -9.44875 17.0051 + vertex 0 -8.60396 18.9025 + vertex 200 -8.60396 18.9025 + endloop + endfacet + facet normal 0 0.913543 -0.406742 + outer loop + vertex 0 -8.60396 18.9025 + vertex 200 -9.44875 17.0051 + vertex 0 -9.44875 17.0051 + endloop + endfacet + facet normal 0 -0.913543 0.406742 + outer loop + vertex 0 8.60396 8.9675 + vertex 200 9.44875 10.8649 + vertex 0 9.44875 10.8649 + endloop + endfacet + facet normal 0 -0.913543 0.406742 + outer loop + vertex 200 9.44875 10.8649 + vertex 0 8.60396 8.9675 + vertex 200 8.60396 8.9675 + endloop + endfacet + facet normal 0 -0.669132 0.743144 + outer loop + vertex 0 7.38314 7.28719 + vertex 200 5.83965 5.89742 + vertex 200 7.38314 7.28719 + endloop + endfacet + facet normal 0 -0.669132 0.743144 + outer loop + vertex 200 5.83965 5.89742 + vertex 0 7.38314 7.28719 + vertex 0 5.83965 5.89742 + endloop + endfacet + facet normal 0 0.809015 -0.587789 + outer loop + vertex 200 -8.60396 18.9025 + vertex 0 -7.38314 20.5828 + vertex 200 -7.38314 20.5828 + endloop + endfacet + facet normal 0 0.809015 -0.587789 + outer loop + vertex 0 -7.38314 20.5828 + vertex 200 -8.60396 18.9025 + vertex 0 -8.60396 18.9025 + endloop + endfacet + facet normal 0 -0.66914 -0.743136 + outer loop + vertex 0 5.83965 21.9726 + vertex 200 7.38314 20.5828 + vertex 200 5.83965 21.9726 + endloop + endfacet + facet normal -0 -0.66914 -0.743136 + outer loop + vertex 200 7.38314 20.5828 + vertex 0 5.83965 21.9726 + vertex 0 7.38314 20.5828 + endloop + endfacet + facet normal 0 -1 0 + outer loop + vertex 0 9.88057 12.8965 + vertex 200 9.88057 14.9735 + vertex 0 9.88057 14.9735 + endloop + endfacet + facet normal 0 -1 -0 + outer loop + vertex 200 9.88057 14.9735 + vertex 0 9.88057 12.8965 + vertex 200 9.88057 12.8965 + endloop + endfacet + facet normal 0 0.913543 0.406742 + outer loop + vertex 200 -8.60396 8.9675 + vertex 0 -9.44875 10.8649 + vertex 200 -9.44875 10.8649 + endloop + endfacet + facet normal 0 0.913543 0.406742 + outer loop + vertex 0 -9.44875 10.8649 + vertex 200 -8.60396 8.9675 + vertex 0 -8.60396 8.9675 + endloop + endfacet + facet normal 0 0.978149 0.207907 + outer loop + vertex 200 -9.44875 10.8649 + vertex 0 -9.88057 12.8965 + vertex 200 -9.88057 12.8965 + endloop + endfacet + facet normal 0 0.978149 0.207907 + outer loop + vertex 0 -9.88057 12.8965 + vertex 200 -9.44875 10.8649 + vertex 0 -9.44875 10.8649 + endloop + endfacet + facet normal 0 -0.309007 -0.95106 + outer loop + vertex 0 2.0656 23.6529 + vertex 200 4.04093 23.0111 + vertex 200 2.0656 23.6529 + endloop + endfacet + facet normal -0 -0.309007 -0.95106 + outer loop + vertex 200 4.04093 23.0111 + vertex 0 2.0656 23.6529 + vertex 0 4.04093 23.0111 + endloop + endfacet + facet normal 0 0.309007 -0.95106 + outer loop + vertex 0 -4.04093 23.0111 + vertex 200 -2.0656 23.6529 + vertex 200 -4.04093 23.0111 + endloop + endfacet + facet normal 0 0.309007 -0.95106 + outer loop + vertex 200 -2.0656 23.6529 + vertex 0 -4.04093 23.0111 + vertex 0 -2.0656 23.6529 + endloop + endfacet + facet normal 0 -0.913543 -0.406742 + outer loop + vertex 0 9.44875 17.0051 + vertex 200 8.60396 18.9025 + vertex 0 8.60396 18.9025 + endloop + endfacet + facet normal 0 -0.913543 -0.406742 + outer loop + vertex 200 8.60396 18.9025 + vertex 0 9.44875 17.0051 + vertex 200 9.44875 17.0051 + endloop + endfacet + facet normal 0 0.978149 -0.207907 + outer loop + vertex 200 -9.88057 14.9735 + vertex 0 -9.44875 17.0051 + vertex 200 -9.44875 17.0051 + endloop + endfacet + facet normal 0 0.978149 -0.207907 + outer loop + vertex 0 -9.44875 17.0051 + vertex 200 -9.88057 14.9735 + vertex 0 -9.88057 14.9735 + endloop + endfacet + facet normal 0 -0.978149 0.207907 + outer loop + vertex 0 9.44875 10.8649 + vertex 200 9.88057 12.8965 + vertex 0 9.88057 12.8965 + endloop + endfacet + facet normal 0 -0.978149 0.207907 + outer loop + vertex 200 9.88057 12.8965 + vertex 0 9.44875 10.8649 + vertex 200 9.44875 10.8649 + endloop + endfacet + facet normal 0 0.104527 -0.994522 + outer loop + vertex 0 -2.0656 23.6529 + vertex 200 0 23.87 + vertex 200 -2.0656 23.6529 + endloop + endfacet + facet normal 0 0.104527 -0.994522 + outer loop + vertex 200 0 23.87 + vertex 0 -2.0656 23.6529 + vertex 0 0 23.87 + endloop + endfacet + facet normal 0 -0.104527 -0.994522 + outer loop + vertex 0 0 23.87 + vertex 200 2.0656 23.6529 + vertex 200 0 23.87 + endloop + endfacet + facet normal -0 -0.104527 -0.994522 + outer loop + vertex 200 2.0656 23.6529 + vertex 0 0 23.87 + vertex 0 2.0656 23.6529 + endloop + endfacet + facet normal 0 0.500003 -0.866024 + outer loop + vertex 0 -5.83965 21.9726 + vertex 200 -4.04093 23.0111 + vertex 200 -5.83965 21.9726 + endloop + endfacet + facet normal 0 0.500003 -0.866024 + outer loop + vertex 200 -4.04093 23.0111 + vertex 0 -5.83965 21.9726 + vertex 0 -4.04093 23.0111 + endloop + endfacet + facet normal 0 -0.500003 -0.866024 + outer loop + vertex 0 4.04093 23.0111 + vertex 200 5.83965 21.9726 + vertex 200 4.04093 23.0111 + endloop + endfacet + facet normal -0 -0.500003 -0.866024 + outer loop + vertex 200 5.83965 21.9726 + vertex 0 4.04093 23.0111 + vertex 0 5.83965 21.9726 + endloop + endfacet + facet normal -0 0.499999 0.866026 + outer loop + vertex 0 -4.04093 4.85893 + vertex 200 -5.83965 5.89742 + vertex 200 -4.04093 4.85893 + endloop + endfacet + facet normal 0 0.499999 0.866026 + outer loop + vertex 200 -5.83965 5.89742 + vertex 0 -4.04093 4.85893 + vertex 0 -5.83965 5.89742 + endloop + endfacet + facet normal 0 -0.809016 0.587786 + outer loop + vertex 0 7.38314 7.28719 + vertex 200 8.60396 8.9675 + vertex 0 8.60396 8.9675 + endloop + endfacet + facet normal 0 -0.809016 0.587786 + outer loop + vertex 200 8.60396 8.9675 + vertex 0 7.38314 7.28719 + vertex 200 7.38314 7.28719 + endloop + endfacet + facet normal 0 -0.809015 -0.587789 + outer loop + vertex 0 8.60396 18.9025 + vertex 200 7.38314 20.5828 + vertex 0 7.38314 20.5828 + endloop + endfacet + facet normal 0 -0.809015 -0.587789 + outer loop + vertex 200 7.38314 20.5828 + vertex 0 8.60396 18.9025 + vertex 200 8.60396 18.9025 + endloop + endfacet + facet normal 0 -0.978149 -0.207907 + outer loop + vertex 0 9.88057 14.9735 + vertex 200 9.44875 17.0051 + vertex 0 9.44875 17.0051 + endloop + endfacet + facet normal 0 -0.978149 -0.207907 + outer loop + vertex 200 9.44875 17.0051 + vertex 0 9.88057 14.9735 + vertex 200 9.88057 14.9735 + endloop + endfacet + facet normal 0 1 -0 + outer loop + vertex 200 -9.88057 12.8965 + vertex 0 -9.88057 14.9735 + vertex 200 -9.88057 14.9735 + endloop + endfacet + facet normal 0 1 0 + outer loop + vertex 0 -9.88057 14.9735 + vertex 200 -9.88057 12.8965 + vertex 0 -9.88057 12.8965 + endloop + endfacet +endsolid OpenSCAD_Model diff --git a/files/benchmark-square.scad b/files/benchmark-square.scad new file mode 100644 index 000000000..d03d154a1 --- /dev/null +++ b/files/benchmark-square.scad @@ -0,0 +1,6 @@ +length = 200; +width = 17.32; +height = 17.32; + +// translate([0, -width/2, 0]) +cube([ length, width, height ]); \ No newline at end of file diff --git a/files/benchmark-square.stl b/files/benchmark-square.stl new file mode 100644 index 000000000..379e3838e --- /dev/null +++ b/files/benchmark-square.stl @@ -0,0 +1,86 @@ +solid OpenSCAD_Model + facet normal -0 0 1 + outer loop + vertex 0 17.32 17.32 + vertex 200 0 17.32 + vertex 200 17.32 17.32 + endloop + endfacet + facet normal 0 0 1 + outer loop + vertex 200 0 17.32 + vertex 0 17.32 17.32 + vertex 0 0 17.32 + endloop + endfacet + facet normal 0 0 -1 + outer loop + vertex 0 0 0 + vertex 200 17.32 0 + vertex 200 0 0 + endloop + endfacet + facet normal -0 0 -1 + outer loop + vertex 200 17.32 0 + vertex 0 0 0 + vertex 0 17.32 0 + endloop + endfacet + facet normal 0 -1 0 + outer loop + vertex 0 0 0 + vertex 200 0 17.32 + vertex 0 0 17.32 + endloop + endfacet + facet normal 0 -1 -0 + outer loop + vertex 200 0 17.32 + vertex 0 0 0 + vertex 200 0 0 + endloop + endfacet + facet normal 1 -0 0 + outer loop + vertex 200 0 17.32 + vertex 200 17.32 0 + vertex 200 17.32 17.32 + endloop + endfacet + facet normal 1 0 0 + outer loop + vertex 200 17.32 0 + vertex 200 0 17.32 + vertex 200 0 0 + endloop + endfacet + facet normal 0 1 -0 + outer loop + vertex 200 17.32 0 + vertex 0 17.32 17.32 + vertex 200 17.32 17.32 + endloop + endfacet + facet normal 0 1 0 + outer loop + vertex 0 17.32 17.32 + vertex 200 17.32 0 + vertex 0 17.32 0 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 0 + vertex 0 17.32 17.32 + vertex 0 17.32 0 + endloop + endfacet + facet normal -1 -0 0 + outer loop + vertex 0 17.32 17.32 + vertex 0 0 0 + vertex 0 0 17.32 + endloop + endfacet +endsolid OpenSCAD_Model diff --git a/go.mod b/go.mod index bb45d5454..dc43f2161 100644 --- a/go.mod +++ b/go.mod @@ -3,6 +3,7 @@ module github.com/deadsy/sdfx go 1.13 require ( + github.com/Megidd/tetrahedron-table v0.0.0-20230901074008-499ef3fcef87 github.com/ajstarks/svgo v0.0.0-20211024235047-1546f124cd8b github.com/dhconnelly/rtreego v1.1.0 github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0 diff --git a/go.sum b/go.sum index 444f00cb0..9a2c2d840 100644 --- a/go.sum +++ b/go.sum @@ -3,6 +3,8 @@ gioui.org v0.0.0-20210308172011-57750fc8a0a6/go.mod h1:RSH6KIUZ0p2xy5zHDxgAM4zum git.sr.ht/~sbinet/gg v0.3.1/go.mod h1:KGYtlADtqsqANL9ueOFkWymvzUvLMQllU5Ixo+8v3pc= github.com/BurntSushi/toml v0.3.1/go.mod h1:xHWCNGjB5oqiDr8zfno3MHue2Ht5sIBksp03qcyfWMU= github.com/BurntSushi/xgb v0.0.0-20160522181843-27f122750802/go.mod h1:IVnqGOEym/WlBOVXweHU+Q+/VP0lqqI8lqeDx9IjBqo= +github.com/Megidd/tetrahedron-table v0.0.0-20230901074008-499ef3fcef87 h1:Q46IyMIux8zFiWM1cP4mXhPldL+z0jLGeKy/2tQTlXQ= +github.com/Megidd/tetrahedron-table v0.0.0-20230901074008-499ef3fcef87/go.mod h1:TlCbcNZKX8/kXirrkWTy7AUgEJoryHM9ahRHLHZtsBc= github.com/ajstarks/deck v0.0.0-20200831202436-30c9fc6549a9/go.mod h1:JynElWSGnm/4RlzPXRlREEwqTHAN3T56Bv2ITsFT3gY= github.com/ajstarks/deck/generate v0.0.0-20210309230005-c3f852c02e19/go.mod h1:T13YZdzov6OU0A1+RfKZiZN9ca6VeKdBdyDV+BY97Tk= github.com/ajstarks/svgo v0.0.0-20180226025133-644b8db467af/go.mod h1:K08gAheRH3/J6wwsYMMT4xOr94bZjxIelGM0+d/wbFw= diff --git a/obj/stl.go b/obj/stl.go index 0dcc36e8a..1c8689303 100644 --- a/obj/stl.go +++ b/obj/stl.go @@ -39,6 +39,10 @@ func (t *triMeshSdf) Evaluate(p v3.Vec) float64 { closestTriangle := math.MaxFloat64 // Quickly skip checking most triangles by only checking the N closest neighbours (AABB based) neighbors := t.rtree.NearestNeighbors(t.numNeighbors, v3ToPoint(p)) + + // To check if all the distances have the same sign. + dists := make([]float64, 0, t.numNeighbors) + for _, neighbor := range neighbors { triangle := neighbor.(*sdf.Triangle3) testPointToTriangle := p.Sub(triangle[0]) @@ -50,10 +54,67 @@ func (t *triMeshSdf) Evaluate(p v3.Vec) float64 { closestTriangle = distToTri signedDistanceResult = signedDistanceToTriPlane } + dists = append(dists, signedDistanceToTriPlane) } + + // Does the approach of this paper make sense: + // https://www2.imm.dtu.dk/pubdb/edoc/imm1289.pdf + // TODO: If so, try to implement it in the future. + + if !sameSign(dists) { + // Sometimes the sign of the final result is not consistent. + signedDistanceResult = signConsistency(dists, signedDistanceResult) + } + return signedDistanceResult } +func sameSign(values []float64) bool { + + positive, negative := false, false + + for _, v := range values { + if v > 0 { + positive = true + } else if v < 0 { + negative = true + } + + // If we've seen both positive and negative, return early + if positive && negative { + return false + } + } + + // All values must have been the same sign + return true +} + +func signConsistency(values []float64, value float64) float64 { + positive := 0 + negative := 0 + + for _, v := range values { + if v > 0 { + positive++ + } else if v < 0 { + negative++ + } + } + + if positive > negative { + if value < 0 { + return -value + } + } else if negative > positive { + if value > 0 { + return -value + } + } + + return value +} + func (t *triMeshSdf) BoundingBox() sdf.Box3 { return t.bb } diff --git a/render/march3.go b/render/march3.go index 892151e25..602c2c9fb 100644 --- a/render/march3.go +++ b/render/march3.go @@ -357,261 +357,517 @@ var mcEdgeTable = [256]int{ // specify the edges used to create the triangle(s) var mcTriangleTable = [256][]int{ + // 0b00000000 case 0 {}, + // 0b00000001 case 1 {0, 8, 3}, + // 0b00000010 case 2 {0, 1, 9}, + // 0b00000011 case 3 {1, 8, 3, 9, 8, 1}, + // 0b00000100 case 4 {1, 2, 10}, + // 0b00000101 case 5 {0, 8, 3, 1, 2, 10}, + // 0b00000110 case 6 {9, 2, 10, 0, 2, 9}, + // 0b00000111 case 7 {2, 8, 3, 2, 10, 8, 10, 9, 8}, + // 0b00001000 case 8 {3, 11, 2}, + // 0b00001001 case 9 {0, 11, 2, 8, 11, 0}, + // 0b00001010 case 10 {1, 9, 0, 2, 3, 11}, + // 0b00001011 case 11 {1, 11, 2, 1, 9, 11, 9, 8, 11}, + // 0b00001100 case 12 {3, 10, 1, 11, 10, 3}, + // 0b00001101 case 13 {0, 10, 1, 0, 8, 10, 8, 11, 10}, + // 0b00001110 case 14 {3, 9, 0, 3, 11, 9, 11, 10, 9}, + // 0b00001111 case 15 {9, 8, 10, 10, 8, 11}, + // 0b00010000 case 16 {4, 7, 8}, + // 0b00010001 case 17 {4, 3, 0, 7, 3, 4}, + // 0b00010010 case 18 {0, 1, 9, 8, 4, 7}, + // 0b00010011 case 19 {4, 1, 9, 4, 7, 1, 7, 3, 1}, + // 0b00010100 case 20 {1, 2, 10, 8, 4, 7}, + // 0b00010101 case 21 {3, 4, 7, 3, 0, 4, 1, 2, 10}, + // 0b00010110 case 22 {9, 2, 10, 9, 0, 2, 8, 4, 7}, + // 0b00010111 case 23 {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4}, + // 0b00011000 case 24 {8, 4, 7, 3, 11, 2}, + // 0b00011001 case 25 {11, 4, 7, 11, 2, 4, 2, 0, 4}, + // 0b00011010 case 26 {9, 0, 1, 8, 4, 7, 2, 3, 11}, + // 0b00011011 case 27 {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1}, + // 0b00011100 case 28 {3, 10, 1, 3, 11, 10, 7, 8, 4}, + // 0b00011101 case 29 {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4}, + // 0b00011110 case 30 {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3}, + // 0b00011111 case 31 {4, 7, 11, 4, 11, 9, 9, 11, 10}, + // 0b00100000 case 32 {9, 5, 4}, + // 0b00100001 case 33 {9, 5, 4, 0, 8, 3}, + // 0b00100010 case 34 {0, 5, 4, 1, 5, 0}, + // 0b00100011 case 35 {8, 5, 4, 8, 3, 5, 3, 1, 5}, + // 0b00100100 case 36 {1, 2, 10, 9, 5, 4}, + // 0b00100101 case 37 {3, 0, 8, 1, 2, 10, 4, 9, 5}, + // 0b00100110 case 38 {5, 2, 10, 5, 4, 2, 4, 0, 2}, + // 0b00100111 case 39 {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8}, + // 0b00101000 case 40 {9, 5, 4, 2, 3, 11}, + // 0b00101001 case 41 {0, 11, 2, 0, 8, 11, 4, 9, 5}, + // 0b00101010 case 42 {0, 5, 4, 0, 1, 5, 2, 3, 11}, + // 0b00101011 case 43 {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5}, + // 0b00101100 case 44 {10, 3, 11, 10, 1, 3, 9, 5, 4}, + // 0b00101101 case 45 {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10}, + // 0b00101110 case 46 {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3}, + // 0b00101111 case 47 {5, 4, 8, 5, 8, 10, 10, 8, 11}, + // 0b00110000 case 48 {9, 7, 8, 5, 7, 9}, + // 0b00110001 case 49 {9, 3, 0, 9, 5, 3, 5, 7, 3}, + // 0b00110010 case 50 {0, 7, 8, 0, 1, 7, 1, 5, 7}, + // 0b00110011 case 51 {1, 5, 3, 3, 5, 7}, + // 0b00110100 case 52 {9, 7, 8, 9, 5, 7, 10, 1, 2}, + // 0b00110101 case 53 {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3}, + // 0b00110110 case 54 {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2}, + // 0b00110111 case 55 {2, 10, 5, 2, 5, 3, 3, 5, 7}, + // 0b00111000 case 56 {7, 9, 5, 7, 8, 9, 3, 11, 2}, + // 0b00111001 case 57 {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11}, + // 0b00111010 case 58 {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7}, + // 0b00111011 case 59 {11, 2, 1, 11, 1, 7, 7, 1, 5}, + // 0b00111100 case 60 {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11}, + // 0b00111101 case 61 {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0}, + // 0b00111110 case 62 {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0}, + // 0b00111111 case 63 {11, 10, 5, 7, 11, 5}, + // 0b01000000 case 64 {10, 6, 5}, + // 0b01000001 case 65 {0, 8, 3, 5, 10, 6}, + // 0b01000010 case 66 {9, 0, 1, 5, 10, 6}, + // 0b01000011 case 67 {1, 8, 3, 1, 9, 8, 5, 10, 6}, + // 0b01000100 case 68 {1, 6, 5, 2, 6, 1}, + // 0b01000101 case 69 {1, 6, 5, 1, 2, 6, 3, 0, 8}, + // 0b01000110 case 70 {9, 6, 5, 9, 0, 6, 0, 2, 6}, + // 0b01000111 case 71 {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8}, + // 0b01001000 case 72 {2, 3, 11, 10, 6, 5}, + // 0b01001001 case 73 {11, 0, 8, 11, 2, 0, 10, 6, 5}, + // 0b01001010 case 74 {0, 1, 9, 2, 3, 11, 5, 10, 6}, + // 0b01001011 case 75 {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11}, + // 0b01001100 case 76 {6, 3, 11, 6, 5, 3, 5, 1, 3}, + // 0b01001101 case 77 {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6}, + // 0b01001110 case 78 {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9}, + // 0b01001111 case 79 {6, 5, 9, 6, 9, 11, 11, 9, 8}, + // 0b01010000 case 80 {5, 10, 6, 4, 7, 8}, + // 0b01010001 case 81 {4, 3, 0, 4, 7, 3, 6, 5, 10}, + // 0b01010010 case 82 {1, 9, 0, 5, 10, 6, 8, 4, 7}, + // 0b01010011 case 83 {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4}, + // 0b01010100 case 84 {6, 1, 2, 6, 5, 1, 4, 7, 8}, + // 0b01010101 case 85 {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7}, + // 0b01010110 case 86 {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6}, + // 0b01010111 case 87 {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9}, + // 0b01011000 case 88 {3, 11, 2, 7, 8, 4, 10, 6, 5}, + // 0b01011001 case 89 {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11}, + // 0b01011010 case 90 {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6}, + // 0b01011011 case 91 {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6}, + // 0b01011100 case 92 {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6}, + // 0b01011101 case 93 {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11}, + // 0b01011110 case 94 {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7}, + // 0b01011111 case 95 {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9}, + // 0b01100000 case 96 {10, 4, 9, 6, 4, 10}, + // 0b01100001 case 97 {4, 10, 6, 4, 9, 10, 0, 8, 3}, + // 0b01100010 case 98 {10, 0, 1, 10, 6, 0, 6, 4, 0}, + // 0b01100011 case 99 {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10}, + // 0b01100100 case 100 {1, 4, 9, 1, 2, 4, 2, 6, 4}, + // 0b01100101 case 101 {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4}, + // 0b01100110 case 102 {0, 2, 4, 4, 2, 6}, + // 0b01100111 case 103 {8, 3, 2, 8, 2, 4, 4, 2, 6}, + // 0b01101000 case 104 {10, 4, 9, 10, 6, 4, 11, 2, 3}, + // 0b01101001 case 105 {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6}, + // 0b01101010 case 106 {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10}, + // 0b01101011 case 107 {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1}, + // 0b01101100 case 108 {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3}, + // 0b01101101 case 109 {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1}, + // 0b01101110 case 110 {3, 11, 6, 3, 6, 0, 0, 6, 4}, + // 0b01101111 case 111 {6, 4, 8, 11, 6, 8}, + // 0b01110000 case 112 {7, 10, 6, 7, 8, 10, 8, 9, 10}, + // 0b01110001 case 113 {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10}, + // 0b01110010 case 114 {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0}, + // 0b01110011 case 115 {10, 6, 7, 10, 7, 1, 1, 7, 3}, + // 0b01110100 case 116 {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7}, + // 0b01110101 case 117 {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9}, + // 0b01110110 case 118 {7, 8, 0, 7, 0, 6, 6, 0, 2}, + // 0b01110111 case 119 {7, 3, 2, 6, 7, 2}, + // 0b01111000 case 120 {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7}, + // 0b01111001 case 121 {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7}, + // 0b01111010 case 122 {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11}, + // 0b01111011 case 123 {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1}, + // 0b01111100 case 124 {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6}, + // 0b01111101 case 125 {0, 9, 1, 11, 6, 7}, + // 0b01111110 case 126 {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0}, + // 0b01111111 case 127 {7, 11, 6}, + // 0b10000000 case 128 {7, 6, 11}, + // 0b10000001 case 129 {3, 0, 8, 11, 7, 6}, + // 0b10000010 case 130 {0, 1, 9, 11, 7, 6}, + // 0b10000011 case 131 {8, 1, 9, 8, 3, 1, 11, 7, 6}, + // 0b10000100 case 132 {10, 1, 2, 6, 11, 7}, + // 0b10000101 case 133 {1, 2, 10, 3, 0, 8, 6, 11, 7}, + // 0b10000110 case 134 {2, 9, 0, 2, 10, 9, 6, 11, 7}, + // 0b10000111 case 135 {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8}, + // 0b10001000 case 136 {7, 2, 3, 6, 2, 7}, + // 0b10001001 case 137 {7, 0, 8, 7, 6, 0, 6, 2, 0}, + // 0b10001010 case 138 {2, 7, 6, 2, 3, 7, 0, 1, 9}, + // 0b10001011 case 139 {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6}, + // 0b10001100 case 140 {10, 7, 6, 10, 1, 7, 1, 3, 7}, + // 0b10001101 case 141 {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8}, + // 0b10001110 case 142 {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7}, + // 0b10001111 case 143 {7, 6, 10, 7, 10, 8, 8, 10, 9}, + // 0b10010000 case 144 {6, 8, 4, 11, 8, 6}, + // 0b10010001 case 145 {3, 6, 11, 3, 0, 6, 0, 4, 6}, + // 0b10010010 case 146 {8, 6, 11, 8, 4, 6, 9, 0, 1}, + // 0b10010011 case 147 {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6}, + // 0b10010100 case 148 {6, 8, 4, 6, 11, 8, 2, 10, 1}, + // 0b10010101 case 149 {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6}, + // 0b10010110 case 150 {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9}, + // 0b10010111 case 151 {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3}, + // 0b10011000 case 152 {8, 2, 3, 8, 4, 2, 4, 6, 2}, + // 0b10011001 case 153 {0, 4, 2, 4, 6, 2}, + // 0b10011010 case 154 {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8}, + // 0b10011011 case 155 {1, 9, 4, 1, 4, 2, 2, 4, 6}, + // 0b10011100 case 156 {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1}, + // 0b10011101 case 157 {10, 1, 0, 10, 0, 6, 6, 0, 4}, + // 0b10011110 case 158 {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3}, + // 0b10011111 case 159 {10, 9, 4, 6, 10, 4}, + // 0b10100000 case 160 {4, 9, 5, 7, 6, 11}, + // 0b10100001 case 161 {0, 8, 3, 4, 9, 5, 11, 7, 6}, + // 0b10100010 case 162 {5, 0, 1, 5, 4, 0, 7, 6, 11}, + // 0b10100011 case 163 {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5}, + // 0b10100100 case 164 {9, 5, 4, 10, 1, 2, 7, 6, 11}, + // 0b10100101 case 165 {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5}, + // 0b10100110 case 166 {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2}, + // 0b10100111 case 167 {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6}, + // 0b10101000 case 168 {7, 2, 3, 7, 6, 2, 5, 4, 9}, + // 0b10101001 case 169 {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7}, + // 0b10101010 case 170 {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0}, + // 0b10101011 case 171 {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8}, + // 0b10101100 case 172 {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7}, + // 0b10101101 case 173 {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4}, + // 0b10101110 case 174 {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10}, + // 0b10101111 case 175 {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10}, + // 0b10110000 case 176 {6, 9, 5, 6, 11, 9, 11, 8, 9}, + // 0b10110001 case 177 {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5}, + // 0b10110010 case 178 {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11}, + // 0b10110011 case 179 {6, 11, 3, 6, 3, 5, 5, 3, 1}, + // 0b10110100 case 180 {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6}, + // 0b10110101 case 181 {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10}, + // 0b10110110 case 182 {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5}, + // 0b10110111 case 183 {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3}, + // 0b10111000 case 184 {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2}, + // 0b10111001 case 185 {9, 5, 6, 9, 6, 0, 0, 6, 2}, + // 0b10111010 case 186 {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8}, + // 0b10111011 case 187 {1, 5, 6, 2, 1, 6}, + // 0b10111100 case 188 {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6}, + // 0b10111101 case 189 {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0}, + // 0b10111110 case 190 {0, 3, 8, 5, 6, 10}, + // 0b10111111 case 191 {10, 5, 6}, + // 0b11000000 case 192 {11, 5, 10, 7, 5, 11}, + // 0b11000001 case 193 {11, 5, 10, 11, 7, 5, 8, 3, 0}, + // 0b11000010 case 194 {5, 11, 7, 5, 10, 11, 1, 9, 0}, + // 0b11000011 case 195 {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1}, + // 0b11000100 case 196 {11, 1, 2, 11, 7, 1, 7, 5, 1}, + // 0b11000101 case 197 {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11}, + // 0b11000110 case 198 {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7}, + // 0b11000111 case 199 {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2}, + // 0b11001000 case 200 {2, 5, 10, 2, 3, 5, 3, 7, 5}, + // 0b11001001 case 201 {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5}, + // 0b11001010 case 202 {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2}, + // 0b11001011 case 203 {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2}, + // 0b11001100 case 204 {1, 3, 5, 3, 7, 5}, + // 0b11001101 case 205 {0, 8, 7, 0, 7, 1, 1, 7, 5}, + // 0b11001110 case 206 {9, 0, 3, 9, 3, 5, 5, 3, 7}, + // 0b11001111 case 207 {9, 8, 7, 5, 9, 7}, + // 0b11010000 case 208 {5, 8, 4, 5, 10, 8, 10, 11, 8}, + // 0b11010001 case 209 {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0}, + // 0b11010010 case 210 {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5}, + // 0b11010011 case 211 {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4}, + // 0b11010100 case 212 {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8}, + // 0b11010101 case 213 {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11}, + // 0b11010110 case 214 {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5}, + // 0b11010111 case 215 {9, 4, 5, 2, 11, 3}, + // 0b11011000 case 216 {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4}, + // 0b11011001 case 217 {5, 10, 2, 5, 2, 4, 4, 2, 0}, + // 0b11011010 case 218 {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9}, + // 0b11011011 case 219 {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2}, + // 0b11011100 case 220 {8, 4, 5, 8, 5, 3, 3, 5, 1}, + // 0b11011101 case 221 {0, 4, 5, 1, 0, 5}, + // 0b11011110 case 222 {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5}, + // 0b11011111 case 223 {9, 4, 5}, + // 0b11100000 case 224 {4, 11, 7, 4, 9, 11, 9, 10, 11}, + // 0b11100001 case 225 {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11}, + // 0b11100010 case 226 {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11}, + // 0b11100011 case 227 {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4}, + // 0b11100100 case 228 {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2}, + // 0b11100101 case 229 {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3}, + // 0b11100110 case 230 {11, 7, 4, 11, 4, 2, 2, 4, 0}, + // 0b11100111 case 231 {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4}, + // 0b11101000 case 232 {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9}, + // 0b11101001 case 233 {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7}, + // 0b11101010 case 234 {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10}, + // 0b11101011 case 235 {1, 10, 2, 8, 7, 4}, + // 0b11101100 case 236 {4, 9, 1, 4, 1, 7, 7, 1, 3}, + // 0b11101101 case 237 {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1}, + // 0b11101110 case 238 {4, 0, 3, 7, 4, 3}, + // 0b11101111 case 239 {4, 8, 7}, + // 0b11110000 case 240 {9, 10, 8, 10, 11, 8}, + // 0b11110001 case 241 {3, 0, 9, 3, 9, 11, 11, 9, 10}, + // 0b11110010 case 242 {0, 1, 10, 0, 10, 8, 8, 10, 11}, + // 0b11110011 case 243 {3, 1, 10, 11, 3, 10}, + // 0b11110100 case 244 {1, 2, 11, 1, 11, 9, 9, 11, 8}, + // 0b11110101 case 245 {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9}, + // 0b11110110 case 246 {0, 2, 11, 8, 0, 11}, + // 0b11110111 case 247 {3, 2, 11}, + // 0b11111000 case 248 {2, 3, 8, 2, 8, 10, 10, 8, 9}, + // 0b11111001 case 249 {9, 10, 2, 0, 9, 2}, + // 0b11111010 case 250 {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8}, + // 0b11111011 case 251 {1, 10, 2}, + // 0b11111100 case 252 {1, 3, 8, 9, 1, 8}, + // 0b11111101 case 253 {0, 9, 1}, + // 0b11111110 case 254 {0, 3, 8}, + // 0b11111111 case 255 {}, } diff --git a/render/marchfe.go b/render/marchfe.go new file mode 100644 index 000000000..889ed61a3 --- /dev/null +++ b/render/marchfe.go @@ -0,0 +1,124 @@ +package render + +import ( + "fmt" + + "github.com/deadsy/sdfx/sdf" + "github.com/deadsy/sdfx/vec/conv" + v3 "github.com/deadsy/sdfx/vec/v3" + "github.com/deadsy/sdfx/vec/v3i" +) + +//----------------------------------------------------------------------------- + +// A finite element can be linear or non-linear. +type Order int + +const ( + Linear Order = iota + 1 // 4-node tetrahedron and 8-node hexahedron + Quadratic // 10-node tetrahedron and 20-node hexahedron +) + +//----------------------------------------------------------------------------- + +// Two shapes of finite element can be generated: tetrahedral and hexahedral. +type Shape int + +const ( + Hexahedral Shape = iota + 1 + Tetrahedral + HexAndTet +) + +//----------------------------------------------------------------------------- + +// MarchingCubesFeUniform renders using marching cubes with uniform space sampling. +type MarchingCubesFeUniform struct { + meshCells int // number of cells on the longest axis of bounding box. e.g 200 + order Order // Linear or quadratic. + shape Shape // Hexahedral, tetrahedral, or both. +} + +// NewMarchingCubesFeUniform returns a RenderHex8 object. +func NewMarchingCubesFeUniform(meshCells int, order Order, shape Shape) *MarchingCubesFeUniform { + return &MarchingCubesFeUniform{ + meshCells: meshCells, + order: order, + shape: shape, + } +} + +// Info returns a string describing the rendered volume. +func (r *MarchingCubesFeUniform) Info(s sdf.SDF3) string { + bb0 := s.BoundingBox() + bb0Size := bb0.Size() + meshInc := bb0Size.MaxComponent() / float64(r.meshCells) + bb1Size := bb0Size.DivScalar(meshInc) + bb1Size = bb1Size.Ceil().AddScalar(1) + cells := conv.V3ToV3i(bb1Size) + return fmt.Sprintf("%dx%dx%d", cells.X, cells.Y, cells.Z) +} + +// Render produces a finite elements mesh over the bounding volume of an sdf3. +// Order and shape of finite elements are selectable. +func (r *MarchingCubesFeUniform) RenderFe(s sdf.SDF3, output sdf.FeWriter) { + // work out the region we will sample + bb0 := s.BoundingBox() + bb0Size := bb0.Size() + meshInc := bb0Size.MaxComponent() / float64(r.meshCells) + bb1Size := bb0Size.DivScalar(meshInc) + bb1Size = bb1Size.Ceil().AddScalar(1) + bb1Size = bb1Size.MulScalar(meshInc) + bb := sdf.NewBox3(bb0.Center(), bb1Size) + marchingCubesFe(s, bb, meshInc, r.order, r.shape, output) +} + +//----------------------------------------------------------------------------- + +// To get the voxel count, dimension, and min/max corner which are consistent with loops of marching algorithm. +// This func loops are exactly like `marchingCubesFe` loops. We have to be consistant. +func (r *MarchingCubesFeUniform) Voxels(s sdf.SDF3) (v3i.Vec, v3.Vec, []v3.Vec, []v3.Vec) { + // work out the region we will sample + bb0 := s.BoundingBox() + bb0Size := bb0.Size() + meshInc := bb0Size.MaxComponent() / float64(r.meshCells) + bb1Size := bb0Size.DivScalar(meshInc) + bb1Size = bb1Size.Ceil().AddScalar(1) + bb1Size = bb1Size.MulScalar(meshInc) + bb := sdf.NewBox3(bb0.Center(), bb1Size) + + size := bb.Size() + base := bb.Min + steps := conv.V3ToV3i(size.DivScalar(meshInc).Ceil()) + inc := size.Div(conv.V3iToV3(steps)) + + nx, ny, nz := steps.X, steps.Y, steps.Z + dx, dy, dz := inc.X, inc.Y, inc.Z + + mins := make([]v3.Vec, 0, nz*nx*ny) + maxs := make([]v3.Vec, 0, nz*nx*ny) + + var p v3.Vec + p.X = base.X + for x := 0; x < nx; x++ { + p.Y = base.Y + for y := 0; y < ny; y++ { + p.Z = base.Z + for z := 0; z < nz; z++ { + x0, y0, z0 := p.X, p.Y, p.Z + x1, y1, z1 := x0+dx, y0+dy, z0+dz + + mins = append(mins, v3.Vec{X: x0, Y: y0, Z: z0}) + maxs = append(maxs, v3.Vec{X: x1, Y: y1, Z: z1}) + + p.Z += dz + } + p.Y += dy + } + p.X += dx + } + + return v3i.Vec{X: nx, Y: ny, Z: nz}, v3.Vec{X: dx, Y: dy, Z: dz}, mins, maxs +} + +//----------------------------------------------------------------------------- diff --git a/render/marchfehelper.go b/render/marchfehelper.go new file mode 100644 index 000000000..17150e8a6 --- /dev/null +++ b/render/marchfehelper.go @@ -0,0 +1,330 @@ +package render + +import ( + "math" + + "github.com/deadsy/sdfx/sdf" + v3 "github.com/deadsy/sdfx/vec/v3" +) + +//----------------------------------------------------------------------------- + +// Specify the point to create the tetrahedra. +// Point can be on edges or corners. +// Index from 0 to 11 means an edge. +// Index from 12 to 19 means a corner. +// Corners were originally indexed from 0 to 7 but they are shifted by 12. +// So, corners are from 0+12 to 7+12 i.e. from 12 to 19. +func point(edges [12]v3.Vec, corners [8]v3.Vec, index int) v3.Vec { + if 0 <= index && index < 12 { + return edges[index] + } else if index < 20 { + return corners[index-12] + } else { + // Should never reach here. + return v3.Vec{} + } +} + +// ----------------------------------------------------------------------------- +// Tries to avoid distorted tetrahedron with +// non-positive Jacobian or with flat shape or with degenerate faces. +// Points that are close to corners, are handled differently. +func mcInterpolateFE(p1, p2 v3.Vec, v1, v2, x float64) v3.Vec { + // Try 80/20 rule. + closeToV1 := math.Abs(x-v1)/math.Abs(v2-v1) < 0.2 + closeToV2 := math.Abs(x-v1)/math.Abs(v2-v1) > 0.8 + + var t float64 + + if closeToV1 { + // Pick a point away from corner + t = 0.2 + } else if closeToV2 { + // Pick a point away from corner + t = 0.8 + } else { + // linear interpolation + t = (x - v1) / (v2 - v1) + } + + return v3.Vec{ + X: p1.X + t*(p2.X-p1.X), + Y: p1.Y + t*(p2.Y-p1.Y), + Z: p1.Z + t*(p2.Z-p1.Z), + } +} + +//----------------------------------------------------------------------------- + +// Is a tetrahedron almost flat? +// To avoid mathematical problem while FEA is run. +// MATHEMATICA script is available here: +// https://math.stackexchange.com/a/4709610/197913 +func almostFlat(a, b, c, d v3.Vec) (bool, float64) { + ab := b.Sub(a) + ac := c.Sub(a) + ad := d.Sub(a) + + // Note that the `Norm` function of MATHEMATICA is equivalent to our `Length()` function. + nab := ab.Length() + ncd := d.Sub(c).Length() + nbd := d.Sub(b).Length() + nbc := c.Sub(b).Length() + nac := ac.Length() + nad := ad.Length() + + // Check for 0 edge lengths + if nab == 0 || ncd == 0 || + nbd == 0 || nbc == 0 || + nac == 0 || nad == 0 { + return true, 0 + } + + volume := 1.0 / 6.0 * math.Abs(ab.Cross(ac).Dot(ad)) + denom := (nab + ncd) * (nac + nbd) * (nad + nbc) + + // Tolerance derived from here would be `480.0`: + // https://math.stackexchange.com/a/4709610/197913 + // A different value is calibrated according to observations. + // TODO: Could be further calibrated. + tolerance := 1000.0 + + rho := tolerance * volume / denom + + return rho < 1, volume +} + +//----------------------------------------------------------------------------- + +// Reference: +// CCX source code: +// ccx_2.20/src/shape4tet.f +func isBadGaussTet4(coords [4]v3.Vec, xi, et, ze float64) (bool, float64) { + // Coordinates of the nodes. + var xl [3][4]float64 + + for i := 0; i < 4; i++ { + xl[0][i] = coords[i].X + xl[1][i] = coords[i].Y + xl[2][i] = coords[i].Z + } + + // Shape functions. + var shp [4][4]float64 + + shp[3][0] = 1.0 - xi - et - ze + shp[3][1] = xi + shp[3][2] = et + shp[3][3] = ze + + // local derivatives of the shape functions: xi-derivative + + shp[0][0] = -1.0 + shp[0][1] = 1.0 + shp[0][2] = 0.0 + shp[0][3] = 0.0 + + // local derivatives of the shape functions: eta-derivative + + shp[1][0] = -1.0 + shp[1][1] = 0.0 + shp[1][2] = 1.0 + shp[1][3] = 0.0 + + // local derivatives of the shape functions: zeta-derivative + + shp[2][0] = -1.0 + shp[2][1] = 0.0 + shp[2][2] = 0.0 + shp[2][3] = 1.0 + + // computation of the local derivative of the global coordinates (xs) + xs := [3][3]float64{} + for i := 0; i < 3; i++ { + for j := 0; j < 3; j++ { + xs[i][j] = 0.0 + for k := 0; k < 4; k++ { + xs[i][j] += xl[i][k] * shp[j][k] + } + } + } + + // computation of the jacobian determinant + xsj := xs[0][0]*(xs[1][1]*xs[2][2]-xs[1][2]*xs[2][1]) - + xs[0][1]*(xs[1][0]*xs[2][2]-xs[1][2]*xs[2][0]) + + xs[0][2]*(xs[1][0]*xs[2][1]-xs[1][1]*xs[2][0]) + + // According to CCX source code to detect nonpositive jacobian determinant in element + // Fortran threshold for non-positive Jacobian determinant is 1e-20. + return xsj < 1e-20, xsj +} + +// Reference: +// CCX source code: +// ccx_2.20/src/shape10tet.f +func isBadGaussTet10(coords [10]v3.Vec, xi, et, ze float64) (bool, float64) { + // Coordinates of the nodes. + var xl [3][10]float64 + + for i := 0; i < 10; i++ { + xl[0][i] = coords[i].X + xl[1][i] = coords[i].Y + xl[2][i] = coords[i].Z + } + + // Shape functions. + var shp [4][10]float64 + + // Shape functions + a := 1.0 - xi - et - ze + shp[3][0] = (2.0*a - 1.0) * a + shp[3][1] = xi * (2.0*xi - 1.0) + shp[3][2] = et * (2.0*et - 1.0) + shp[3][3] = ze * (2.0*ze - 1.0) + shp[3][4] = 4.0 * xi * a + shp[3][5] = 4.0 * xi * et + shp[3][6] = 4.0 * et * a + shp[3][7] = 4.0 * ze * a + shp[3][8] = 4.0 * xi * ze + shp[3][9] = 4.0 * et * ze + + // Local derivatives of the shape functions: xi-derivative + shp[0][0] = 1.0 - 4.0*a + shp[0][1] = 4.0*xi - 1.0 + shp[0][2] = 0.0 + shp[0][3] = 0.0 + shp[0][4] = 4.0 * (a - xi) + shp[0][5] = 4.0 * et + shp[0][6] = -4.0 * et + shp[0][7] = -4.0 * ze + shp[0][8] = 4.0 * ze + shp[0][9] = 0.0 + + // Local derivatives of the shape functions: eta-derivative + shp[1][0] = 1.0 - 4.0*a + shp[1][1] = 0.0 + shp[1][2] = 4.0*et - 1.0 + shp[1][3] = 0.0 + shp[1][4] = -4.0 * xi + shp[1][5] = 4.0 * xi + shp[1][6] = 4.0 * (a - et) + shp[1][7] = -4.0 * ze + shp[1][8] = 0.0 + shp[1][9] = 4.0 * ze + + // Local derivatives of the shape functions: zeta-derivative + shp[2][0] = 1.0 - 4.0*a + shp[2][1] = 0.0 + shp[2][2] = 0.0 + shp[2][3] = 4.0*ze - 1.0 + shp[2][4] = -4.0 * xi + shp[2][5] = 0.0 + shp[2][6] = -4.0 * et + shp[2][7] = 4.0 * (a - ze) + shp[2][8] = 4.0 * xi + shp[2][9] = 4.0 * et + + // Computation of the local derivative of the global coordinates (xs) + var xs [3][3]float64 + for i := 0; i < 3; i++ { + for j := 0; j < 3; j++ { + xs[i][j] = 0.0 + for k := 0; k < 10; k++ { + xs[i][j] = xs[i][j] + xl[i][k]*shp[j][k] + } + } + } + + // computation of the jacobian determinant + xsj := xs[0][0]*(xs[1][1]*xs[2][2]-xs[1][2]*xs[2][1]) - + xs[0][1]*(xs[1][0]*xs[2][2]-xs[1][2]*xs[2][0]) + + xs[0][2]*(xs[1][0]*xs[2][1]-xs[1][1]*xs[2][0]) + + // According to CCX source code to detect nonpositive jacobian determinant in element + // Fortran threshold for non-positive Jacobian determinant is 1e-20. + return xsj < 1e-20, xsj +} + +//----------------------------------------------------------------------------- + +// Exactly follow CCX source code leading to this error: +// *ERROR in e_c3d: nonpositive jacobian determinant in element +func isBadTet4(coords [4]v3.Vec) (bool, float64) { + + // xi, et, and ze are the coordinates of the Gauss point + // in the integration scheme for the 4-node tetrahedral element. + // For this element type, there is typically only 1 Gauss point used, + // which is located at the centroid of the tetrahedron. + // The coordinates of this Gauss point are (xi, et, ze) = (1/4, 1/4, 1/4). + // Reference: + // ccx_2.20/src/gauss.f + var xi float64 = 0.25 + var et float64 = 0.25 + var ze float64 = 0.25 + + return isBadGaussTet4(coords, xi, et, ze) +} + +// Exactly follow CCX source code leading to this error: +// *ERROR in e_c3d: nonpositive jacobian determinant in element +func isBadTet10(coords [10]v3.Vec) (bool, float64) { + // Gause points are according to CCX source code. + // Reference: + // ccx_2.20/src/gauss.f + var gaussPoints [4]v3.Vec + gaussPoints[0] = v3.Vec{0.138196601125011, 0.138196601125011, 0.138196601125011} + gaussPoints[1] = v3.Vec{0.585410196624968, 0.138196601125011, 0.138196601125011} + gaussPoints[2] = v3.Vec{0.138196601125011, 0.585410196624968, 0.138196601125011} + gaussPoints[3] = v3.Vec{0.138196601125011, 0.138196601125011, 0.585410196624968} + + var bad bool + var jacobianDeterminant float64 + + for i := 0; i < 4; i++ { + bad, jacobianDeterminant = isBadGaussTet10(coords, gaussPoints[i].X, gaussPoints[i].Y, gaussPoints[i].Z) + if bad { + return true, jacobianDeterminant + } + } + + return false, jacobianDeterminant +} + +//----------------------------------------------------------------------------- + +// If triangles are degenerate, then tetrahedra will be bad. +// To filter bad tetrahedra. +func degenerateTriangles(a, b, c, d v3.Vec) bool { + // 4 triangles are possible. + // Each triangle is a tetrahedron side. + t := sdf.Triangle3{} + t[0] = a + t[1] = b + t[2] = c + // Use the epsilon value of `vertexbuffer.go` + if t.Degenerate(0.0001) { + return true + } + t[0] = a + t[1] = b + t[2] = d + // Use the epsilon value of `vertexbuffer.go` + if t.Degenerate(0.0001) { + return true + } + t[0] = a + t[1] = c + t[2] = d + // Use the epsilon value of `vertexbuffer.go` + if t.Degenerate(0.0001) { + return true + } + t[0] = b + t[1] = c + t[2] = d + // Use the epsilon value of `vertexbuffer.go` + return t.Degenerate(0.0001) +} + +//----------------------------------------------------------------------------- diff --git a/render/marchfelogic.go b/render/marchfelogic.go new file mode 100644 index 000000000..dfee25d08 --- /dev/null +++ b/render/marchfelogic.go @@ -0,0 +1,342 @@ +package render + +import ( + "github.com/Megidd/tetrahedron-table/src/gotable" + "github.com/deadsy/sdfx/sdf" + "github.com/deadsy/sdfx/vec/conv" + v3 "github.com/deadsy/sdfx/vec/v3" +) + +func marchingCubesFe(s sdf.SDF3, box sdf.Box3, step float64, order Order, shape Shape, output sdf.FeWriter) { + size := box.Size() + base := box.Min + steps := conv.V3ToV3i(size.DivScalar(step).Ceil()) + inc := size.Div(conv.V3iToV3(steps)) + + // start the evaluation routines + evalRoutines() + + // create the SDF layer cache + l := newLayerYZ(base, inc, steps) + // evaluate the SDF for x = 0 + l.Evaluate(s, 0) + + nx, ny, nz := steps.X, steps.Y, steps.Z + dx, dy, dz := inc.X, inc.Y, inc.Z + + var p v3.Vec + p.X = base.X + for x := 0; x < nx; x++ { + // read the x + 1 layer + l.Evaluate(s, x+1) + // process all cubes in the x and x + 1 layers + p.Y = base.Y + for y := 0; y < ny; y++ { + p.Z = base.Z + for z := 0; z < nz; z++ { + x0, y0, z0 := p.X, p.Y, p.Z + x1, y1, z1 := x0+dx, y0+dy, z0+dz + corners := [8]v3.Vec{ + {x0, y0, z0}, + {x1, y0, z0}, + {x1, y1, z0}, + {x0, y1, z0}, + {x0, y0, z1}, + {x1, y0, z1}, + {x1, y1, z1}, + {x0, y1, z1}} + values := [8]float64{ + l.Get(0, y, z), + l.Get(1, y, z), + l.Get(1, y+1, z), + l.Get(0, y+1, z), + l.Get(0, y, z+1), + l.Get(1, y, z+1), + l.Get(1, y+1, z+1), + l.Get(0, y+1, z+1)} + output.Write(mcToFE(corners, values, x, y, z, order, shape)) + p.Z += dz + } + p.Y += dy + } + p.X += dx + } +} + +//----------------------------------------------------------------------------- + +func mcToFE(corners [8]v3.Vec, values [8]float64, x, y, z int, order Order, shape Shape) []*sdf.Fe { + var fes []*sdf.Fe + switch order { + case Linear: + { + switch shape { + case Hexahedral: + { + fes = append(fes, mcToHex8(corners, values, 0, x, y, z)...) + } + case Tetrahedral: + { + fes = append(fes, mcToTet4(corners, values, 0, x, y, z)...) + } + case HexAndTet: + { + // If all cube corners are inside surface mesh, a single hexahedral element is generated. + // If all cube corners are outside surface mesh, no element is generated. + // If cube is colliding with surface mesh, one or more tetrahedral elements are generated. + tmp := mcToHex8(corners, values, 0, x, y, z) + if len(tmp) < 1 { + tmp = mcToTet4(corners, values, 0, x, y, z) + } + fes = append(fes, tmp...) + } + } + } + case Quadratic: + { + switch shape { + case Hexahedral: + { + fes = append(fes, mcToHex20(corners, values, 0, x, y, z)...) + } + case Tetrahedral: + { + fes = append(fes, mcToTet10(corners, values, 0, x, y, z)...) + } + case HexAndTet: + { + // If all cube corners are inside surface mesh, a single hexahedral element is generated. + // If all cube corners are outside surface mesh, no element is generated. + // If cube is colliding with surface mesh, one or more tetrahedral elements are generated. + tmp := mcToHex20(corners, values, 0, x, y, z) + if len(tmp) < 1 { + tmp = mcToTet10(corners, values, 0, x, y, z) + } + fes = append(fes, tmp...) + } + } + } + } + + return fes +} + +//----------------------------------------------------------------------------- + +func mcToHex8(p [8]v3.Vec, v [8]float64, x float64, layerX, layerY, layerZ int) []*sdf.Fe { + result := make([]*sdf.Fe, 0) + + anyPositive := false + for i := 0; i < 8; i++ { + if v[i] > 0 { + anyPositive = true + break + } + } + + // Create a finite element if all 8 values are non-positive. + // Finite element is inside the 3D model if all values are non-positive. + + if !anyPositive { + fe := sdf.Fe{ + V: make([]v3.Vec, 8), + X: layerX, + Y: layerY, + Z: layerZ, + } + + // Refer to CalculiX solver documentation: + // http://www.dhondt.de/ccx_2.20.pdf + + fe.V[7] = p[7] + fe.V[6] = p[6] + fe.V[5] = p[5] + fe.V[4] = p[4] + fe.V[3] = p[3] + fe.V[2] = p[2] + fe.V[1] = p[1] + fe.V[0] = p[0] + result = append(result, &fe) + } + + return result +} + +//----------------------------------------------------------------------------- + +func mcToHex20(p [8]v3.Vec, v [8]float64, x float64, layerX, layerY, layerZ int) []*sdf.Fe { + result := make([]*sdf.Fe, 0) + + anyPositive := false + for i := 0; i < 8; i++ { + if v[i] > 0 { + anyPositive = true + break + } + } + + // Create a finite element if all 8 values are non-positive. + // Finite element is inside the 3D model if all values are non-positive. + + if !anyPositive { + fe := sdf.Fe{ + V: make([]v3.Vec, 20), + X: layerX, + Y: layerY, + Z: layerZ, + } + + // Refer to CalculiX solver documentation: + // http://www.dhondt.de/ccx_2.20.pdf + + // Points on cube corners: + fe.V[7] = p[7] + fe.V[6] = p[6] + fe.V[5] = p[5] + fe.V[4] = p[4] + fe.V[3] = p[3] + fe.V[2] = p[2] + fe.V[1] = p[1] + fe.V[0] = p[0] + + // Points on cube edges: + fe.V[8] = p[0].Add(p[1]).MulScalar(0.5) + fe.V[9] = p[1].Add(p[2]).MulScalar(0.5) + fe.V[10] = p[2].Add(p[3]).MulScalar(0.5) + fe.V[11] = p[3].Add(p[0]).MulScalar(0.5) + + fe.V[12] = p[4].Add(p[5]).MulScalar(0.5) + fe.V[13] = p[5].Add(p[6]).MulScalar(0.5) + fe.V[14] = p[6].Add(p[7]).MulScalar(0.5) + fe.V[15] = p[7].Add(p[4]).MulScalar(0.5) + + fe.V[16] = p[0].Add(p[4]).MulScalar(0.5) + fe.V[17] = p[1].Add(p[5]).MulScalar(0.5) + fe.V[18] = p[2].Add(p[6]).MulScalar(0.5) + fe.V[19] = p[3].Add(p[7]).MulScalar(0.5) + + result = append(result, &fe) + } + + return result +} + +//----------------------------------------------------------------------------- + +func mcToTet4(p [8]v3.Vec, v [8]float64, x float64, layerX, layerY, layerZ int) []*sdf.Fe { + // which of the 0..255 patterns do we have? + index := 0 + for i := 0; i < 8; i++ { + if v[i] < x { + index |= 1 << uint(i) + } + } + + // work out the interpolated points on the edges + var points [12]v3.Vec + for i := 0; i < 12; i++ { + bit := 1 << uint(i) + if mcEdgeTable[index]&bit != 0 { + a := mcPairTable[i][0] + b := mcPairTable[i][1] + points[i] = mcInterpolateFE(p[a], p[b], v[a], v[b], x) + } + } + + // Create the tetrahedra. + table := gotable.TetrahedronTable[index] + count := len(table) / 4 + result := make([]*sdf.Fe, 0, count) + for i := 0; i < count; i++ { + t := sdf.Fe{ + V: make([]v3.Vec, 4), + X: layerX, + Y: layerY, + Z: layerZ, + } + + t.V[0] = point(points, p, table[i*4+0]) + t.V[1] = point(points, p, table[i*4+1]) + t.V[2] = point(points, p, table[i*4+2]) + t.V[3] = point(points, p, table[i*4+3]) + degenerated := degenerateTriangles(t.V[0], t.V[1], t.V[2], t.V[3]) + flat, _ := almostFlat(t.V[0], t.V[1], t.V[2], t.V[3]) + + // In the case of marching cubes algorithm to generate triangle, it's avoiding zero-area triangles by `!t.Degenerate(0)` check. + // In our case of marching cubes algorithm to generate tetrahedron, we can do a check too: + bad, _ := isBadTet4([4]v3.Vec{t.V[0], t.V[1], t.V[2], t.V[3]}) + if !degenerated && !bad && !flat { + result = append(result, &t) + } else { + // CCX solver may throw error for this element. So, skip it. + // *ERROR in e_c3d: nonpositive jacobian determinant in element + } + } + + return result +} + +//----------------------------------------------------------------------------- + +func mcToTet10(p [8]v3.Vec, v [8]float64, x float64, layerX, layerY, layerZ int) []*sdf.Fe { + // which of the 0..255 patterns do we have? + index := 0 + for i := 0; i < 8; i++ { + if v[i] < x { + index |= 1 << uint(i) + } + } + + // work out the interpolated points on the edges + var points [12]v3.Vec + for i := 0; i < 12; i++ { + bit := 1 << uint(i) + if mcEdgeTable[index]&bit != 0 { + a := mcPairTable[i][0] + b := mcPairTable[i][1] + points[i] = mcInterpolateFE(p[a], p[b], v[a], v[b], x) + } + } + + // Create the tetrahedra. + table := gotable.TetrahedronTable[index] + count := len(table) / 4 + result := make([]*sdf.Fe, 0, count) + for i := 0; i < count; i++ { + t := sdf.Fe{ + V: make([]v3.Vec, 10), + X: layerX, + Y: layerY, + Z: layerZ, + } + + // Points on tetrahedron corners. + t.V[0] = point(points, p, table[i*4+0]) + t.V[1] = point(points, p, table[i*4+1]) + t.V[2] = point(points, p, table[i*4+2]) + t.V[3] = point(points, p, table[i*4+3]) + degenerated := degenerateTriangles(t.V[0], t.V[1], t.V[2], t.V[3]) + flat, _ := almostFlat(t.V[0], t.V[1], t.V[2], t.V[3]) + // Points on tetrahedron edges. + // Followoing CalculiX node numbering. + t.V[4] = t.V[0].Add(t.V[1]).MulScalar(0.5) + t.V[5] = t.V[1].Add(t.V[2]).MulScalar(0.5) + t.V[6] = t.V[0].Add(t.V[2]).MulScalar(0.5) + t.V[7] = t.V[0].Add(t.V[3]).MulScalar(0.5) + t.V[8] = t.V[1].Add(t.V[3]).MulScalar(0.5) + t.V[9] = t.V[2].Add(t.V[3]).MulScalar(0.5) + // In the case of marching cubes algorithm to generate triangle, it's avoiding zero-area triangles by `!t.Degenerate(0)` check. + // In our case of marching cubes algorithm to generate tetrahedron, we can do a check too: + bad, _ := isBadTet10([10]v3.Vec{t.V[0], t.V[1], t.V[2], t.V[3], t.V[4], t.V[5], t.V[6], t.V[7], t.V[8], t.V[9]}) + if !degenerated && !bad && !flat { + result = append(result, &t) + } else { + // CCX solver may throw error for this element. So, skip it. + // *ERROR in e_c3d: nonpositive jacobian determinant in element + } + } + + return result +} + +//----------------------------------------------------------------------------- diff --git a/render/render.go b/render/render.go index 1a5501ace..beebf41c9 100644 --- a/render/render.go +++ b/render/render.go @@ -13,6 +13,8 @@ import ( "sync" "github.com/deadsy/sdfx/sdf" + v3 "github.com/deadsy/sdfx/vec/v3" + "github.com/deadsy/sdfx/vec/v3i" ) //----------------------------------------------------------------------------- @@ -29,6 +31,13 @@ type Render2 interface { Info(s sdf.SDF2) string } +// RenderFe renders a finite element mesh over the bounding volume of an sdf3. +type RenderFe interface { + RenderFe(sdf3 sdf.SDF3, output sdf.FeWriter) + Info(sdf3 sdf.SDF3) string + Voxels(sdf3 sdf.SDF3) (v3i.Vec, v3.Vec, []v3.Vec, []v3.Vec) +} + //----------------------------------------------------------------------------- // ToTriangles renders an SDF3 to a triangle mesh. @@ -52,6 +61,36 @@ func ToTriangles( //----------------------------------------------------------------------------- +// ToFem renders an SDF3 to finite elements. +func ToFem( + s sdf.SDF3, // sdf3 to render + r RenderFe, // rendering method +) []sdf.Fe { + fmt.Printf("rendering %s\n", r.Info(s)) + + voxelCount, _, _, _ := r.Voxels(s) + fmt.Printf("voxel counts of marching algorithm are: (%v x %v x %v)\n", voxelCount.X, voxelCount.Y, voxelCount.Z) + + // Will be filled by the rendering. + fes := make([]sdf.Fe, 0) + + var wg sync.WaitGroup + + // To write the finite elements. + output := sdf.WriteFes(&wg, &fes) + + // run the renderer + r.RenderFe(s, sdf.NewFeBuffer(output)) + // stop the writer reading on the channel + close(output) + // wait for the file write to complete + wg.Wait() + + return fes +} + +//----------------------------------------------------------------------------- + // ToSTL renders an SDF3 to an STL file. func ToSTL( s sdf.SDF3, // sdf3 to render diff --git a/sdf/fe.go b/sdf/fe.go new file mode 100644 index 000000000..59b04bdad --- /dev/null +++ b/sdf/fe.go @@ -0,0 +1,103 @@ +package sdf + +import ( + "sync" + + v3 "github.com/deadsy/sdfx/vec/v3" +) + +// Fe is a finite element. +type Fe struct { + // Coordinates of nodes or vertices. + V []v3.Vec + // Coordinates of the voxel to which the element belongs. + X int + Y int + Z int +} + +//----------------------------------------------------------------------------- + +// WriteFes writes a stream of finite elements to a slice. +func WriteFes(wg *sync.WaitGroup, elements *[]Fe) chan<- []*Fe { + // External code writes to this channel. + // This goroutine reads the channel and stores finite elements. + c := make(chan []*Fe) + + wg.Add(1) + go func() { + defer wg.Done() + // read finite elements from the channel and handle them + for fes := range c { + for _, fe := range fes { + *elements = append(*elements, *fe) + } + } + }() + + return c +} + +//----------------------------------------------------------------------------- +// Finite element Buffering + +// We write finite elements to a channel to decouple the rendering routines from the +// routine that writes file output. We have a lot of finite elements and channels +// are not very fast, so it's best to bundle many finite elements into a single channel +// write. The renderer doesn't naturally do that, so we buffer finite elements before +// writing them to the channel. + +// FeWriter is the interface of a finite element writer/closer object. +type FeWriter interface { + Write(in []*Fe) error + Close() error +} + +// size the buffer to avoid re-allocations when appending. +const feBufferSize = 256 + +// marching cubes produces 0 or 1 finite element type of hex. +// marching cubes produces 0 to less than 20 finite element of type tet. +// According to https://github.com/Megidd/tetrahedron-table +// +// TODO: can value be further calibrated? +const feBufferMargin = 20 + +// FeBuffer buffers finite elements before writing them to a channel. +type FeBuffer struct { + buf []*Fe // finite element buffer + out chan<- []*Fe // output channel + lock sync.Mutex // lock the the buffer during access +} + +// NewFeBuffer returns a FeBuffer. +func NewFeBuffer(out chan<- []*Fe) FeWriter { + return &FeBuffer{ + buf: make([]*Fe, 0, feBufferSize+feBufferMargin), + out: out, + } +} + +func (a *FeBuffer) Write(in []*Fe) error { + a.lock.Lock() + a.buf = append(a.buf, in...) + if len(a.buf) >= tBufferSize { + a.out <- a.buf + a.buf = make([]*Fe, 0, feBufferSize+feBufferMargin) + } + a.lock.Unlock() + return nil +} + +// Close flushes out any remaining finite elements in the buffer. +func (a *FeBuffer) Close() error { + a.lock.Lock() + if len(a.buf) != 0 { + a.out <- a.buf + a.buf = nil + } + a.lock.Unlock() + return nil +} + +//----------------------------------------------------------------------------- diff --git a/sdf/finiteelements/buffer/buffer.go b/sdf/finiteelements/buffer/buffer.go new file mode 100644 index 000000000..5a11ed477 --- /dev/null +++ b/sdf/finiteelements/buffer/buffer.go @@ -0,0 +1,2 @@ +// Package buffer implements index and vertex buffers for finite element mesh. +package buffer diff --git a/sdf/finiteelements/buffer/clean.go b/sdf/finiteelements/buffer/clean.go new file mode 100644 index 000000000..65d6eeb74 --- /dev/null +++ b/sdf/finiteelements/buffer/clean.go @@ -0,0 +1,32 @@ +package buffer + +import "fmt" + +// Rather than connecting disconnected components, we can keep the largest one and delete the rest. +func (vg *VoxelGrid) CleanDisconnections(components []*Component) { + // Find largest component. Consider volume criterion. + maxComponentIndex := -1 + maxVoxelCount := 0 + for i, component := range components { + if component.VoxelCount() > maxVoxelCount { + maxVoxelCount = component.VoxelCount() + maxComponentIndex = i + } + } + + if maxComponentIndex != -1 { + fmt.Printf("Component %v has the largest voxel count: %v\n", maxComponentIndex, maxVoxelCount) + } else { + fmt.Printf("No components found.") + return + } + + // Remove elements inside the voxels of smaller components. + for i, component := range components { + if i != maxComponentIndex { + for v := range component.Voxels { + vg.DelAll(v[0], v[1], v[2]) + } + } + } +} diff --git a/sdf/finiteelements/buffer/component.go b/sdf/finiteelements/buffer/component.go new file mode 100644 index 000000000..580e1b7c5 --- /dev/null +++ b/sdf/finiteelements/buffer/component.go @@ -0,0 +1,127 @@ +package buffer + +// Component represents a set of connected elements. +type Component struct { + Voxels map[[3]int]struct{} // A map avoids repeated voxels. +} + +func (c *Component) VoxelCount() int { + return len(c.Voxels) +} + +// Count separate components consisting of disconnected finite elements. +// They cause FEA solver to throw error. +func (vg *VoxelGrid) Components() []*Component { + visited := make(map[*Element]bool) + components := make([]*Component, 0) + process := func(x, y, z int, els []*Element) { + for _, el := range els { + if !visited[el] { + component := vg.bfs(visited, el, [3]int{x, y, z}) + components = append(components, component) + } + } + } + vg.Iterate(process) + return components +} + +// Algorithm: breadth-first search (bfs). +func (vg *VoxelGrid) bfs(visited map[*Element]bool, start *Element, startV [3]int) *Component { + queue := []*Element{start} + quVox := [][3]int{startV} // To store the voxel of each element. + visited[start] = true + + component := &Component{ + Voxels: make(map[[3]int]struct{}, 0), + } + component.Voxels[startV] = struct{}{} + + for len(queue) > 0 { + e := queue[0] + v := quVox[0] + queue = queue[1:] + quVox = quVox[1:] + + neighbors, neighVoxs := vg.neighbors(e, v) + + for i := range neighbors { + n := neighbors[i] + nv := neighVoxs[i] + if !visited[n] { + visited[n] = true + queue = append(queue, n) + quVox = append(quVox, nv) + + component.Voxels[nv] = struct{}{} + } + } + } + + return component +} + +// It returns a list of neighbors. +func (vg *VoxelGrid) neighbors(e *Element, v [3]int) ([]*Element, [][3]int) { + var neighbors []*Element + var neighVoxs [][3]int + + for i := -1; i <= 1; i++ { + for j := -1; j <= 1; j++ { + for k := -1; k <= 1; k++ { + x := v[0] + i + y := v[1] + j + z := v[2] + k + + if !vg.isValid(x, y, z) { + continue + } + + for _, el := range vg.Get(x, y, z) { + if i == 0 && j == 0 && k == 0 { + // The same voxel: skip the same element. + if el == e { + continue + } + } + if sharesNode(e, el) { + neighbors = append(neighbors, el) + neighVoxs = append(neighVoxs, [3]int{x, y, z}) + } + } + } + } + } + + return neighbors, neighVoxs +} + +func sharesNode(e1, e2 *Element) bool { + // The node count doesn't need to be equal for the two elements. + // Since, the two elements could be of different types. + + for _, n1 := range e1.Nodes { + if contains(e2.Nodes, n1) { + return true + } + } + return false +} + +func contains(arr []uint32, i uint32) bool { + for _, n := range arr { + if n == i { + return true + } + } + return false +} + +//----------------------------------------------------------------------------- + +// Is voxel index inside a valid range? +func (vg *VoxelGrid) isValid(x, y, z int) bool { + return x >= 0 && y >= 0 && z >= 0 && x < vg.Len.X && y < vg.Len.Y && z < vg.Len.Z +} + +//----------------------------------------------------------------------------- diff --git a/sdf/finiteelements/buffer/indexbuffer.go b/sdf/finiteelements/buffer/indexbuffer.go new file mode 100644 index 000000000..d38e712ea --- /dev/null +++ b/sdf/finiteelements/buffer/indexbuffer.go @@ -0,0 +1,41 @@ +package buffer + +import ( + v3 "github.com/deadsy/sdfx/vec/v3" + "github.com/deadsy/sdfx/vec/v3i" +) + +// Index buffer for a mesh of finite elements. +type IB struct { + Grid *VoxelGrid +} + +func NewIB(voxelLen v3i.Vec, voxelDim v3.Vec, mins, maxs []v3.Vec) *IB { + ib := IB{ + Grid: NewVoxelGrid(voxelLen, voxelDim, mins, maxs), + } + + return &ib +} + +// Add a finite element to buffer. +// Voxel coordinate and nodes are input. +// The node numbering should follow the convention of CalculiX. +// http://www.dhondt.de/ccx_2.20.pdf +func (ib *IB) AddFE(x, y, z int, nodes []uint32) { + ib.Grid.Append(x, y, z, NewElement(nodes)) +} + +// To delete all elements off of a voxel. +func (ib *IB) DelAll(x, y, z int) { + ib.Grid.DelAll(x, y, z) +} + +// To iterate over all voxels and get elements inside each voxel and do stuff with them. +func (ib *IB) Iterate(f func(int, int, int, []*Element)) { + ib.Grid.Iterate(f) +} + +func (ib *IB) Size() (int, int, int) { + return ib.Grid.Size() +} diff --git a/sdf/finiteelements/buffer/vertexbuffer.go b/sdf/finiteelements/buffer/vertexbuffer.go new file mode 100644 index 000000000..84c36bfad --- /dev/null +++ b/sdf/finiteelements/buffer/vertexbuffer.go @@ -0,0 +1,65 @@ +package buffer + +import ( + "runtime" + + v3 "github.com/deadsy/sdfx/vec/v3" +) + +// Vertex buffer for finite elements. +// Vertex buffer avoids repeating vertices by a hash table. +// The same vertex buffer is used for 4-node tetrahedra, 8-node hexahedra, and others. +type VB struct { + // To store index of vertices. Repeated vertices would have the same index. + hashTable map[[3]float64]uint32 + // To store coordinates of vertices. + V []v3.Vec +} + +func NewVB() *VB { + b := VB{ + hashTable: map[[3]float64]uint32{}, + V: []v3.Vec{}, + } + + return &b +} + +// Add vertex to buffer and get vertex ID. +// If vertex is already available on the buffer, its ID is just returned. +// So, all vertices will be unique. Not repeated. +func (b *VB) Id(v v3.Vec) uint32 { + // Do NOT remove any small details. Small details matter. + // Removing small details would cause this bug: + // https://calculix.discourse.group/t/detect-bad-finite-elements-4-node-tetrahedral/1700/5?u=megidd + key := [3]float64{v.X, v.Y, v.Z} + if vID, ok := b.hashTable[key]; ok { + // Vertex already exists. It's repeated. + return vID + } + + // Vertex is new, so append it. + b.V = append(b.V, v) + + // Store index of the appended vertex. + b.hashTable[key] = uint32(b.VertexCount() - 1) + + // Return index of the appended vertex. + return uint32(b.VertexCount() - 1) +} + +func (b *VB) VertexCount() int { + return len(b.V) +} + +// To be called after adding all vertices to the vertex buffer. +// Call if you are sure that no new vertex will be added to the vertex buffer. +func (b *VB) DestroyHashTable() { + // Clear memory. + b.hashTable = nil + runtime.GC() +} + +func (b *VB) Vertex(i uint32) v3.Vec { + return b.V[i] +} diff --git a/sdf/finiteelements/buffer/voxel.go b/sdf/finiteelements/buffer/voxel.go new file mode 100644 index 000000000..5fd1a0c33 --- /dev/null +++ b/sdf/finiteelements/buffer/voxel.go @@ -0,0 +1,222 @@ +package buffer + +import ( + v3 "github.com/deadsy/sdfx/vec/v3" + "github.com/deadsy/sdfx/vec/v3i" +) + +type Element struct { + Nodes []uint32 // Node indices +} + +// Declare the enum using iota and const +type ElementType int + +const ( + C3D4 ElementType = iota + 1 + C3D10 + C3D8 + C3D20R + Unknown +) + +func (e *Element) Type() ElementType { + if len(e.Nodes) == 4 { + return C3D4 + } else if len(e.Nodes) == 10 { + return C3D10 + } else if len(e.Nodes) == 8 { + return C3D8 + } else if len(e.Nodes) == 20 { + return C3D20R + } + return Unknown +} + +func NewElement(nodes []uint32) *Element { + e := Element{ + Nodes: nodes, + } + return &e +} + +type Voxel struct { + data []*Element // Each voxel stores multiple elements. + Min v3.Vec // Min corner of voxel. + Max v3.Vec // Max corner of voxel. +} + +func NewVoxel(min, max v3.Vec) *Voxel { + return &Voxel{ + data: make([]*Element, 0), + Min: min, + Max: max, + } +} + +// Acts like a three-dimensional nested slice using +// a one-dimensional slice under the hood. +// To increase performance. +type VoxelGrid struct { + Voxels []*Voxel // + Len v3i.Vec // Voxel count in 3 directions. + Dim v3.Vec // Voxel dimension in 3 directions. +} + +func NewVoxelGrid(len v3i.Vec, dim v3.Vec, mins, maxs []v3.Vec) *VoxelGrid { + vg := &VoxelGrid{ + Voxels: make([]*Voxel, len.X*len.Y*len.Z), + Len: len, + Dim: dim, + } + + // Assign the min corner and max corner of each voxel. + for i := range vg.Voxels { + vg.Voxels[i] = NewVoxel(mins[i], maxs[i]) + } + + return vg +} + +// This func must be consistent with `(r *MarchingCubesFEUniform) Voxels` func. +// This func must be consistent with `marchingCubesFE` func too. +func (vg *VoxelGrid) index1Dto3D(i int) (int, int, int) { + z := i % vg.Len.Z + y := (i / vg.Len.Z) % vg.Len.Y + x := i / (vg.Len.Z * vg.Len.Y) + return x, y, z +} + +// This func must be consistent with `(r *MarchingCubesFEUniform) Voxels` func. +// This func must be consistent with `marchingCubesFE` func too. +func (vg *VoxelGrid) index3Dto1D(x, y, z int) int { + return x*vg.Len.Y*vg.Len.Z + y*vg.Len.Z + z +} + +func (vg *VoxelGrid) Size() (int, int, int) { + return vg.Len.X, vg.Len.Y, vg.Len.Z +} + +// To get all the elements inside a voxel. +func (vg *VoxelGrid) Get(x, y, z int) []*Element { + return vg.Voxels[vg.index3Dto1D(x, y, z)].data +} + +// To set all the elements inside a voxel at once. +func (vg *VoxelGrid) Set(x, y, z int, value []*Element) { + vg.Voxels[vg.index3Dto1D(x, y, z)].data = value +} + +// To append a single element to the elements inside a voxel. +func (vg *VoxelGrid) Append(x, y, z int, value *Element) { + vg.Voxels[vg.index3Dto1D(x, y, z)].data = append(vg.Voxels[vg.index3Dto1D(x, y, z)].data, value) +} + +// To delete all elements off of a voxel. +func (vg *VoxelGrid) DelAll(x, y, z int) { + vg.Voxels[vg.index3Dto1D(x, y, z)].data = nil + vg.Voxels[vg.index3Dto1D(x, y, z)].data = make([]*Element, 0) +} + +// Compute the bounding box of all the input points. +// Return all the voxels that are intersecting with that bounding box. +func (vg *VoxelGrid) VoxelsIntersecting(points []v3.Vec) ([]v3i.Vec, v3.Vec, v3.Vec) { + if len(points) == 0 { + return nil, v3.Vec{}, v3.Vec{} + } + + // compute the bounding box of all the input points + min, max := points[0], points[0] + for _, point := range points { + min = min.Min(point) + max = max.Max(point) + } + + var intersectingVoxels []v3i.Vec + + // iterate over all the voxels + for i, voxel := range vg.Voxels { + // check if the voxel intersects with the bounding box + if !DoesIntersect(voxel.Min, voxel.Max, min, max) { + continue + } + + // convert the 1D index to a 3D index + x, y, z := vg.index1Dto3D(i) + + intersectingVoxels = append(intersectingVoxels, v3i.Vec{X: x, Y: y, Z: z}) + } + + return intersectingVoxels, min, max +} + +// Does two voxels or b-boxes intersect with each other? +func DoesIntersect(aMin, aMax, bMin, bMax v3.Vec) bool { + if aMin.X > bMax.X || aMin.Y > bMax.Y || aMin.Z > bMax.Z { + return false + } + if bMin.X > aMax.X || bMin.Y > aMax.Y || bMin.Z > aMax.Z { + return false + } + return true +} + +// Return all the voxels that are intersecting with a single point. +func (vg *VoxelGrid) VoxelsIntersectingWithPoint(point v3.Vec) []v3i.Vec { + var intersectingVoxels []v3i.Vec + + // iterate over all the voxels + for i, voxel := range vg.Voxels { + // check if the voxel intersects with the bounding box + if !DoesIntersectPointWithBox(point, voxel.Min, voxel.Max) { + continue + } + + // convert the 1D index to a 3D index + x, y, z := vg.index1Dto3D(i) + + intersectingVoxels = append(intersectingVoxels, v3i.Vec{X: x, Y: y, Z: z}) + } + + return intersectingVoxels +} + +// Does a point and a bounding box intersect with each other? +func DoesIntersectPointWithBox(point, boxMin, boxMax v3.Vec) bool { + if point.X < boxMin.X || point.Y < boxMin.Y || point.Z < boxMin.Z { + return false + } + if point.X > boxMax.X || point.Y > boxMax.Y || point.Z > boxMax.Z { + return false + } + return true +} + +// To iterate over all voxels and get elements inside each voxel and do stuff with them. +func (vg *VoxelGrid) Iterate(f func(int, int, int, []*Element)) { + for z := 0; z < vg.Len.Z; z++ { + for y := 0; y < vg.Len.Y; y++ { + for x := 0; x < vg.Len.X; x++ { + value := vg.Get(x, y, z) + f(x, y, z, value) + } + } + } +} + +//----------------------------------------------------------------------------- + +func (vg *VoxelGrid) VoxelsOn1stLayerZ() []v3i.Vec { + voxels := make([]v3i.Vec, 0) + z := 0 // 1st layer along Z axis. + for y := 0; y < vg.Len.Y; y++ { + for x := 0; x < vg.Len.X; x++ { + elements := vg.Get(x, y, z) + if len(elements) < 1 { + continue + } + voxels = append(voxels, v3i.Vec{X: x, Y: y, Z: z}) + } + } + return voxels +} diff --git a/sdf/finiteelements/mesh/fem.go b/sdf/finiteelements/mesh/fem.go new file mode 100644 index 000000000..1a9e2f64f --- /dev/null +++ b/sdf/finiteelements/mesh/fem.go @@ -0,0 +1,253 @@ +package mesh + +import ( + "log" + "math" + + "github.com/deadsy/sdfx/render" + "github.com/deadsy/sdfx/sdf" + "github.com/deadsy/sdfx/sdf/finiteelements/buffer" + v3 "github.com/deadsy/sdfx/vec/v3" + "github.com/deadsy/sdfx/vec/v3i" +) + +// Fem is a mesh of finite elements. +// A sophisticated data structure for mesh is required. +// The repeated nodes would be removed. +// The element connectivity would be created with unique nodes. +type Fem struct { + // Index buffer. + IBuff *buffer.IB + // Vertex buffer. + VBuff *buffer.VB +} + +// NewFem returns a new mesh and number of its voxel layers along X, Y, Z axis. +func NewFem(s sdf.SDF3, r render.RenderFe) (*Fem, int, int, int) { + fes := render.ToFem(s, r) + + voxelLen, voxelDim, mins, maxs := r.Voxels(s) + + m := newFem(voxelLen, voxelDim, mins, maxs) + + // Fill out the mesh with finite elements. + for _, fe := range fes { + m.addFE(fe.X, fe.Y, fe.Z, fe.V) + } + + defer m.VBuff.DestroyHashTable() + + return m, voxelLen.X, voxelLen.Y, voxelLen.Z +} + +func newFem(voxelLen v3i.Vec, voxelDim v3.Vec, mins, maxs []v3.Vec) *Fem { + return &Fem{ + IBuff: buffer.NewIB(voxelLen, voxelDim, mins, maxs), + VBuff: buffer.NewVB(), + } +} + +func (m *Fem) Size() (int, int, int) { + return m.IBuff.Size() +} + +// Add a finite element. +// Voxel coordinate and nodes are input. +// The node numbering should follow the convention of CalculiX. +// http://www.dhondt.de/ccx_2.20.pdf +func (m *Fem) addFE(x, y, z int, nodes []v3.Vec) { + indices := make([]uint32, len(nodes)) + for n := 0; n < len(nodes); n++ { + indices[n] = m.addVertex(nodes[n]) + } + m.IBuff.AddFE(x, y, z, indices) +} + +func (m *Fem) addVertex(vert v3.Vec) uint32 { + return m.VBuff.Id(vert) +} + +func (m *Fem) vertexCount() int { + return m.VBuff.VertexCount() +} + +func (m *Fem) vertex(i uint32) v3.Vec { + return m.VBuff.Vertex(i) +} + +// To iterate over all voxels and get elements inside each voxel and do stuff with them. +func (m *Fem) iterate(f func(int, int, int, []*buffer.Element)) { + m.IBuff.Iterate(f) +} + +// The closest vertex/node is identified. +// Also, the containing voxel is figured out. +// +// This logic has to be here, since we need access to any node vertex. +func (m *Fem) Locate(location v3.Vec) (v3.Vec, v3i.Vec) { + // Calculating voxel indices. + idxX := int(math.Floor((location.X - m.IBuff.Grid.Voxels[0].Min.X) / (m.IBuff.Grid.Dim.X))) + idxY := int(math.Floor((location.Y - m.IBuff.Grid.Voxels[0].Min.Y) / (m.IBuff.Grid.Dim.Y))) + idxZ := int(math.Floor((location.Z - m.IBuff.Grid.Voxels[0].Min.Z) / (m.IBuff.Grid.Dim.Z))) + + // Ensure indices are within bounds + if idxX >= m.IBuff.Grid.Len.X { + idxX = m.IBuff.Grid.Len.X - 1 + } + if idxY >= m.IBuff.Grid.Len.Y { + idxY = m.IBuff.Grid.Len.Y - 1 + } + if idxZ >= m.IBuff.Grid.Len.Z { + idxZ = m.IBuff.Grid.Len.Z - 1 + } + // Sometimes a location on the STL geometry, is outside the voxels. + // Since the marching cubes algorithm might not be consistent with the exact STL geometry boundaries. + if idxX <= 0 { + idxX = 0 + } + if idxY <= 0 { + idxY = 0 + } + if idxZ <= 0 { + idxZ = 0 + } + + // Get elements in the voxel + elements := m.IBuff.Grid.Get(idxX, idxY, idxZ) + + // Find the closest node + var closestNode v3.Vec + minDistance := math.Inf(1) + + for _, element := range elements { + for _, node := range element.Nodes { + // A function that gives you the position of a node. + nodePos := m.vertex(node) + + distance := location.Sub(nodePos).Length() + if distance < minDistance { + minDistance = distance + closestNode = nodePos + } + } + } + + return closestNode, v3i.Vec{X: idxX, Y: idxY, Z: idxZ} +} + +// Compute the bounding box of all the input points. +// Return all the voxels that are intersecting with that bounding box. +func (m *Fem) VoxelsIntersecting(points []v3.Vec) ([]v3i.Vec, v3.Vec, v3.Vec) { + return m.IBuff.Grid.VoxelsIntersecting(points) +} + +//----------------------------------------------------------------------------- + +func (m *Fem) VoxelsOn1stLayerZ() []v3i.Vec { + return m.IBuff.Grid.VoxelsOn1stLayerZ() +} + +//----------------------------------------------------------------------------- + +// Count separate components consisting of disconnected finite elements. +// They cause FEA solver to throw error. +func (m *Fem) Components() []*buffer.Component { + return m.IBuff.Grid.Components() +} + +func (m *Fem) CleanDisconnections(components []*buffer.Component) { + m.IBuff.Grid.CleanDisconnections(components) +} + +//----------------------------------------------------------------------------- + +// WriteInp writes mesh to ABAQUS or CalculiX `inp` file. +func (m *Fem) WriteInp( + path string, + massDensity float32, + youngModulus float32, + poissonRatio float32, + restraints []*Restraint, + loads []*Load, + gravityDirection v3.Vec, + gravityMagnitude float64, + gravityIsNeeded bool, +) error { + _, _, layersZ := m.IBuff.Size() + return m.WriteInpLayers(path, 0, layersZ, massDensity, youngModulus, poissonRatio, restraints, loads, gravityDirection, gravityMagnitude, gravityIsNeeded) +} + +// WriteInpLayers writes specific layers of mesh to ABAQUS or CalculiX `inp` file. +// Result would include start layer. +// Result would exclude end layer. +func (m *Fem) WriteInpLayers( + path string, + layerStart, layerEnd int, + massDensity float32, + youngModulus float32, + poissonRatio float32, + restraints []*Restraint, + loads []*Load, + gravityDirection v3.Vec, + gravityMagnitude float64, + gravityIsNeeded bool, +) error { + restraints = restraintSetup(m, restraints) + loads = loadSetup(m, loads) + inp := NewInp(m, path, layerStart, layerEnd, massDensity, youngModulus, poissonRatio, restraints, loads, gravityDirection, gravityMagnitude, gravityIsNeeded) + return inp.Write() +} + +//----------------------------------------------------------------------------- + +func restraintSetup(m *Fem, restraints []*Restraint) []*Restraint { + // Figure out voxel for each. + for _, r := range restraints { + // Set voxel, if not already set. + // If voxel is already set, it means the caller has decided about voxel. + // It means all the nodes inside the voxel will be restraint. + if r.voxel.X == -1 && r.voxel.Y == -1 && r.voxel.Z == -1 { + voxels := m.IBuff.Grid.VoxelsIntersectingWithPoint(r.Location) + _, closestVoxel := m.Locate(r.Location) + + if len(voxels) < 1 { + log.Printf("no voxel is intersecting with the point restraint: %f, %f, %f\n", r.Location.X, r.Location.Y, r.Location.Z) + log.Println("...no worries, the closest voxel will be assumed.") + r.voxel = closestVoxel + } else { + if voxels[0].X != closestVoxel.X && voxels[0].Y != closestVoxel.Y && voxels[0].Z != closestVoxel.Z { + log.Println("point restraint is outside voxel boundaries: m.VoxelsIntersecting() != m.Locate()") + log.Println("...no worries, the closest voxel will be assumed.") + r.voxel = closestVoxel + } else { + r.voxel = voxels[0] + } + } + } + } + return restraints +} + +func loadSetup(m *Fem, loads []*Load) []*Load { + // Figure out voxel for each. + for _, l := range loads { + voxels, _, _ := m.VoxelsIntersecting([]v3.Vec{l.Location}) + closestVertex, closestVoxel := m.Locate(l.Location) + + if len(voxels) < 1 { + log.Printf("no voxel is intersecting with the point load: %f, %f, %f\n", l.Location.X, l.Location.Y, l.Location.Z) + log.Println("...no worries, the closest voxel will be assumed.") + } else { + if voxels[0].X != closestVoxel.X && voxels[0].Y != closestVoxel.Y && voxels[0].Z != closestVoxel.Z { + log.Println("point load is outside voxel boundaries: m.VoxelsIntersecting() != m.Locate()") + log.Println("...no worries, the closest voxel will be assumed.") + } + } + + l.voxel = closestVoxel + l.nodeREF = closestVertex + } + return loads +} + +//----------------------------------------------------------------------------- diff --git a/sdf/finiteelements/mesh/inp.go b/sdf/finiteelements/mesh/inp.go new file mode 100644 index 000000000..b5829157a --- /dev/null +++ b/sdf/finiteelements/mesh/inp.go @@ -0,0 +1,706 @@ +package mesh + +import ( + "fmt" + "os" + "time" + + "github.com/deadsy/sdfx/sdf/finiteelements/buffer" + v3 "github.com/deadsy/sdfx/vec/v3" +) + +// Inp writes different types of finite elements as ABAQUS or CalculiX `inp` file. +type Inp struct { + // Finite elements mesh. + Mesh *Fem + // Output `inp` file path. + Path string + // For writing nodes to a separate file. + PathNodes string + // For writing elements to a separate file. + PathElsC3D4 string + // For writing elements to a separate file. + PathElsC3D10 string + // For writing elements to a separate file. + PathElsC3D8 string + // For writing elements to a separate file. + PathElsC3D20R string + // For writing boundary conditions to a separate file. + PathBou string + // For writing loads to a separate file. + PathLoad string + // For writing gravity to a separate file. + PathGravity string + GravityIsNeeded bool + // Output `inp` file would include start layer. + LayerStart int + // Output `inp` file would exclude end layer. + LayerEnd int + // To write only required nodes to `inp` file. + TempVBuff *buffer.VB + // Mechanical properties of 3D print resin. + MassDensity float32 + YoungModulus float32 + PoissonRatio float32 + // Just a counter to keep track of written elements + eleID uint32 + // Just a counter to keep track of written nodes + nextNode uint32 + // Single point constraint: one or more degrees of freedom are fixed for a given node. Passed in by the caller. + Restraints []*Restraint + // Point loads are applied to the nodes of the mesh. Passed in by the caller. + Loads []*Load + // Assigns gravity loading in this direction to all elements. + GravityDirection v3.Vec + // Assigns gravity loading with magnitude to all elements. + GravityMagnitude float64 +} + +// NewInp sets up a new writer. +func NewInp( + m *Fem, + path string, + layerStart, layerEnd int, + massDensity float32, youngModulus float32, poissonRatio float32, + restraints []*Restraint, + loads []*Load, + gravityDirection v3.Vec, + gravityMagnitude float64, + gravityIsNeeded bool, +) *Inp { + inp := &Inp{ + Mesh: m, + Path: path, + PathNodes: path + ".nodes", + PathElsC3D4: path + ".elements_C3D4", + PathElsC3D10: path + ".elements_C3D10", + PathElsC3D8: path + ".elements_C3D8", + PathElsC3D20R: path + ".elements_C3D20R", + PathBou: path + ".boundary", + PathLoad: path + ".load", + PathGravity: path + ".gravity", + GravityIsNeeded: gravityIsNeeded, + LayerStart: layerStart, + LayerEnd: layerEnd, + TempVBuff: buffer.NewVB(), + MassDensity: massDensity, + YoungModulus: youngModulus, + PoissonRatio: poissonRatio, + Restraints: restraints, + Loads: loads, + GravityDirection: gravityDirection, + GravityMagnitude: gravityMagnitude, + } + + return inp +} + +// Write starts writing to `inp` file. +func (inp *Inp) Write() error { + f, err := os.Create(inp.Path) + if err != nil { + return err + } + defer f.Close() + + err = inp.writeHeader(f) + if err != nil { + return err + } + + // Write nodes. + + // Include a separate file to avoid cluttering the `inp` file. + _, err = f.WriteString(fmt.Sprintf("*INCLUDE,INPUT=%s\n", inp.PathNodes)) + if err != nil { + return err + } + + // Temp buffer is just to avoid writing repeated nodes into the `inpt` file. + defer inp.TempVBuff.DestroyHashTable() + + err = inp.writeNodes() + if err != nil { + return err + } + + // Write elements. + + // Include a separate file to avoid cluttering the `inp` file. + _, err = f.WriteString(fmt.Sprintf("*INCLUDE,INPUT=%s\n", inp.PathElsC3D4)) + if err != nil { + return err + } + // Include a separate file to avoid cluttering the `inp` file. + _, err = f.WriteString(fmt.Sprintf("*INCLUDE,INPUT=%s\n", inp.PathElsC3D10)) + if err != nil { + return err + } + // Include a separate file to avoid cluttering the `inp` file. + _, err = f.WriteString(fmt.Sprintf("*INCLUDE,INPUT=%s\n", inp.PathElsC3D8)) + if err != nil { + return err + } + // Include a separate file to avoid cluttering the `inp` file. + _, err = f.WriteString(fmt.Sprintf("*INCLUDE,INPUT=%s\n", inp.PathElsC3D20R)) + if err != nil { + return err + } + + err = inp.writeElements() + if err != nil { + return err + } + + // Fix the degrees of freedom one through three for all nodes on specific layers. + + // Include a separate file to avoid cluttering the `inp` file. + _, err = f.WriteString(fmt.Sprintf("*INCLUDE,INPUT=%s\n", inp.PathBou)) + if err != nil { + return err + } + + err = inp.writeBoundary() + if err != nil { + return err + } + + return inp.writeFooter(f) +} + +func (inp *Inp) writeHeader(f *os.File) error { + _, _, layersZ := inp.Mesh.Size() + if 0 <= inp.LayerStart && inp.LayerStart < inp.LayerEnd && inp.LayerEnd <= layersZ { + // Good. + } else { + return fmt.Errorf("start or end layer is beyond range") + } + + _, err := f.WriteString("**\n** Structure: finite elements of a 3D model.\n** Generated by: https://github.com/deadsy/sdfx\n**\n") + if err != nil { + return err + } + + _, err = f.WriteString("*HEADING\nModel: 3D model Date: " + time.Now().UTC().Format("2006-Jan-02 MST") + "\n") + if err != nil { + return err + } + + return nil +} + +func (inp *Inp) writeNodes() error { + // Write to a separate file to avoid cluttering the `inp` file. + f, err := os.Create(inp.PathNodes) + if err != nil { + return err + } + defer f.Close() + + _, err = f.WriteString("*NODE\n") + if err != nil { + return err + } + + var process func(int, int, int, []*buffer.Element) + + inp.nextNode = 1 // ID starts from one not zero. + + process = func(x, y, z int, els []*buffer.Element) { + if z >= inp.LayerStart && z < inp.LayerEnd { + // Good. + } else { + return + } + + for _, el := range els { + vertices := make([]v3.Vec, len(el.Nodes)) + ids := make([]uint32, len(el.Nodes)) + for n := 0; n < len(el.Nodes); n++ { + vertices[n] = inp.Mesh.vertex(el.Nodes[n]) + ids[n] = inp.TempVBuff.Id(vertices[n]) + } + + // Write the node IDs. + for n := 0; n < len(el.Nodes); n++ { + // Only write node if it's not already written to file. + if ids[n]+1 == inp.nextNode { + // ID starts from one not zero. + _, err = f.WriteString(fmt.Sprintf("%d,%f,%f,%f\n", ids[n]+1, float32(vertices[n].X), float32(vertices[n].Y), float32(vertices[n].Z))) + if err != nil { + panic("Couldn't write node to file: " + err.Error()) + } + inp.nextNode++ + } + + } + } + } + + inp.Mesh.iterate(process) + + return nil +} + +func (inp *Inp) writeElements() error { + // Write to a separate file to avoid cluttering the `inp` file. + fC3D4, err := os.Create(inp.PathElsC3D4) + if err != nil { + return err + } + defer fC3D4.Close() + + // Write to a separate file to avoid cluttering the `inp` file. + fC3D10, err := os.Create(inp.PathElsC3D10) + if err != nil { + return err + } + defer fC3D10.Close() + + // Write to a separate file to avoid cluttering the `inp` file. + fC3D8, err := os.Create(inp.PathElsC3D8) + if err != nil { + return err + } + defer fC3D8.Close() + + // Write to a separate file to avoid cluttering the `inp` file. + fC3D20R, err := os.Create(inp.PathElsC3D20R) + if err != nil { + return err + } + defer fC3D20R.Close() + + _, err = fC3D4.WriteString(fmt.Sprintf("*ELEMENT, TYPE=%s, ELSET=eC3D4\n", "C3D4")) + if err != nil { + return err + } + + _, err = fC3D10.WriteString(fmt.Sprintf("*ELEMENT, TYPE=%s, ELSET=eC3D10\n", "C3D10")) + if err != nil { + return err + } + + _, err = fC3D8.WriteString(fmt.Sprintf("*ELEMENT, TYPE=%s, ELSET=eC3D8\n", "C3D8")) + if err != nil { + return err + } + + _, err = fC3D20R.WriteString(fmt.Sprintf("*ELEMENT, TYPE=%s, ELSET=eC3D20R\n", "C3D20R")) + if err != nil { + return err + } + + // Define a function variable with the signature + var process func(int, int, int, []*buffer.Element) + // Assign a function literal to the variable + process = func(x, y, z int, els []*buffer.Element) { + if z >= inp.LayerStart && z < inp.LayerEnd { + // Good. + } else { + return + } + for _, el := range els { + ids := make([]uint32, len(el.Nodes)) + for n := 0; n < len(el.Nodes); n++ { + vertex := inp.Mesh.vertex(el.Nodes[n]) + ids[n] = inp.TempVBuff.Id(vertex) + } + + // ID starts from one not zero. + + switch el.Type() { + case buffer.C3D4: + { + _, err = fC3D4.WriteString(fmt.Sprintf("%d,%d,%d,%d,%d\n", inp.eleID+1, ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1)) + } + case buffer.C3D10: + { + _, err = fC3D10.WriteString(fmt.Sprintf("%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d\n", inp.eleID+1, ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, ids[4]+1, ids[5]+1, ids[6]+1, ids[7]+1, ids[8]+1, ids[9]+1)) + } + case buffer.C3D8: + { + _, err = fC3D8.WriteString(fmt.Sprintf("%d,%d,%d,%d,%d,%d,%d,%d,%d\n", inp.eleID+1, ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, ids[4]+1, ids[5]+1, ids[6]+1, ids[7]+1)) + } + case buffer.C3D20R: + { + // There should not be more than 16 entries in a line; + // That's why there is new line in the middle. + // Refer to CalculiX solver documentation: + // http://www.dhondt.de/ccx_2.20.pdf + _, err = fC3D20R.WriteString(fmt.Sprintf("%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,\n%d,%d,%d,%d,%d\n", inp.eleID+1, ids[0]+1, ids[1]+1, ids[2]+1, ids[3]+1, ids[4]+1, ids[5]+1, ids[6]+1, ids[7]+1, ids[8]+1, ids[9]+1, ids[10]+1, ids[11]+1, ids[12]+1, ids[13]+1, ids[14]+1, ids[15]+1, ids[16]+1, ids[17]+1, ids[18]+1, ids[19]+1)) + } + case buffer.Unknown: + { + fmt.Println("Element has unknown type :(") + } + } + + if err != nil { + panic("Couldn't write finite element to file: " + err.Error()) + } + + inp.eleID++ + } + } + + inp.Mesh.iterate(process) + + return nil +} + +func (inp *Inp) writeBoundary() error { + // Write to a separate file to avoid cluttering the `inp` file. + f, err := os.Create(inp.PathBou) + if err != nil { + return err + } + defer f.Close() + + for i, r := range inp.Restraints { + isFixedX, isFixedY, isFixedZ := r.IsFixedX, r.IsFixedY, r.IsFixedZ + if !isFixedX && !isFixedY && !isFixedZ { + return fmt.Errorf("restraint has no fixed degree of freedom") + } + + nodeSet := make([]uint32, 0) + + elements := inp.Mesh.IBuff.Grid.Get(r.voxel.X, r.voxel.Y, r.voxel.Z) + for _, element := range elements { + for _, node := range element.Nodes { + // Node ID should be consistant with the temp vertex buffer. + // Node ID is different on these two: (1) original vertex buffer, (2) temp vertex buffer. + vertex := inp.Mesh.vertex(node) + id := inp.TempVBuff.Id(vertex) + nodeSet = append(nodeSet, id) + } + } + + // Write node set for this restraint. + _, err = f.WriteString(fmt.Sprintf("*NSET,NSET=restraint%d\n", i+1)) + if err != nil { + return err + } + + for j, id := range nodeSet { + if j == len(nodeSet)-1 { + // The last one is written differently. + _, err = f.WriteString(fmt.Sprintf("%d\n", id+1)) + if err != nil { + return err + } + } else if j != 0 && j%15 == 0 { + // According to CCX manual: maximum 16 entries per line. + _, err = f.WriteString(fmt.Sprintf("%d,\n", id+1)) + if err != nil { + return err + } + } else { + _, err = f.WriteString(fmt.Sprintf("%d,", id+1)) + if err != nil { + return err + } + } + } + + // Put the boundary constraints on the reference node. + _, err = f.WriteString("*BOUNDARY\n") + if err != nil { + return err + } + + // To be written: + // + // 1) Node number/ID or node set label. + // 2) First degree of freedom constrained. + // 3) Last degree of freedom constrained. This field may be left blank if only + // one degree of freedom is constrained. + // + // Note: written node ID would start from one not zero. + + if isFixedX && isFixedY && isFixedZ { + _, err = f.WriteString(fmt.Sprintf("restraint%d,1\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + _, err = f.WriteString(fmt.Sprintf("restraint%d,2\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + _, err = f.WriteString(fmt.Sprintf("restraint%d,3\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + } else if isFixedX && isFixedY && !isFixedZ { + _, err = f.WriteString(fmt.Sprintf("restraint%d,1\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + _, err = f.WriteString(fmt.Sprintf("restraint%d,2\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + } else if !isFixedX && isFixedY && isFixedZ { + _, err = f.WriteString(fmt.Sprintf("restraint%d,2\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + _, err = f.WriteString(fmt.Sprintf("restraint%d,3\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + } else if isFixedX && !isFixedY && isFixedZ { + _, err = f.WriteString(fmt.Sprintf("restraint%d,1\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + _, err = f.WriteString(fmt.Sprintf("restraint%d,3\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + } else if isFixedX && !isFixedY && !isFixedZ { + _, err = f.WriteString(fmt.Sprintf("restraint%d,1\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + } else if !isFixedX && isFixedY && !isFixedZ { + _, err = f.WriteString(fmt.Sprintf("restraint%d,2\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + } else if !isFixedX && !isFixedY && isFixedZ { + _, err = f.WriteString(fmt.Sprintf("restraint%d,3\n", i+1)) + if err != nil { + panic("Couldn't write boundary to file: " + err.Error()) + } + } + } + + return nil +} + +func (inp *Inp) writeLoad() error { + // Write to a separate file to avoid cluttering the `inp` file. + f, err := os.Create(inp.PathLoad) + if err != nil { + return err + } + defer f.Close() + + _, err = f.WriteString("*CLOAD\n") + if err != nil { + return err + } + + // The closest node to any restraint is already computed. + for _, l := range inp.Loads { + // Node ID should be consistant with the temp vertex buffer. + // Node ID is different on these two: (1) original vertex buffer, (2) temp vertex buffer. + id := inp.TempVBuff.Id(l.nodeREF) + + // To be written: + // + // 1) Node ID. + // 2) Degree of freedom. + // 3) Magnitude of the load. + // + // Note: written node ID would start from one not zero. + + _, err = f.WriteString(fmt.Sprintf("%d,1,%f\n", id+1, l.Magnitude.X)) + if err != nil { + panic("Couldn't write load to file: " + err.Error()) + } + _, err = f.WriteString(fmt.Sprintf("%d,2,%f\n", id+1, l.Magnitude.Y)) + if err != nil { + panic("Couldn't write load to file: " + err.Error()) + } + _, err = f.WriteString(fmt.Sprintf("%d,3,%f\n", id+1, l.Magnitude.Z)) + if err != nil { + panic("Couldn't write load to file: " + err.Error()) + } + } + + return nil +} + +func (inp *Inp) writeGravity() error { + // Write to a separate file to avoid cluttering the `inp` file. + f, err := os.Create(inp.PathGravity) + if err != nil { + return err + } + defer f.Close() + + // Write distributed loads. + + _, err = f.WriteString("*DLOAD\n") + if err != nil { + return err + } + + // Assign gravity loading in any direction with any magnitude to all elements. + // + // 9810 could be gravity magnitude in mm/sec^2 units. + // + // SLA 3D printing is done upside-down. 3D model is hanging from the print floor. + // That's why gravity could be in "positive" z-direction. + // Here ”gravity” really stands for any acceleration vector. + // + // Refer to CalculiX solver documentation: + // http://www.dhondt.de/ccx_2.20.pdf + _, err = f.WriteString( + fmt.Sprintf( + "eC3D4,GRAV,%v,%v,%v,%v\n", + inp.GravityMagnitude, + inp.GravityDirection.X, inp.GravityDirection.Y, inp.GravityDirection.Z, + ), + ) + if err != nil { + return err + } + _, err = f.WriteString( + fmt.Sprintf( + "eC3D10,GRAV,%v,%v,%v,%v\n", + inp.GravityMagnitude, + inp.GravityDirection.X, inp.GravityDirection.Y, inp.GravityDirection.Z, + ), + ) + if err != nil { + return err + } + _, err = f.WriteString( + fmt.Sprintf( + "eC3D8,GRAV,%v,%v,%v,%v\n", + inp.GravityMagnitude, + inp.GravityDirection.X, inp.GravityDirection.Y, inp.GravityDirection.Z, + ), + ) + if err != nil { + return err + } + _, err = f.WriteString( + fmt.Sprintf( + "eC3D20R,GRAV,%v,%v,%v,%v\n", + inp.GravityMagnitude, + inp.GravityDirection.X, inp.GravityDirection.Y, inp.GravityDirection.Z, + ), + ) + + return err +} + +func (inp *Inp) writeFooter(f *os.File) error { + + // Define material. + // Units of measurement are mm,N,s,K. + // Refer to: + // https://engineering.stackexchange.com/q/54454/15178 + // Refer to: + // Units chapter of CalculiX solver documentation: + // http://www.dhondt.de/ccx_2.20.pdf + + _, err := f.WriteString("*MATERIAL, name=resin\n") + if err != nil { + return err + } + + _, err = f.WriteString(fmt.Sprintf("*ELASTIC,TYPE=ISO\n%e,%e,0\n", inp.YoungModulus, inp.PoissonRatio)) + if err != nil { + return err + } + + _, err = f.WriteString(fmt.Sprintf("*DENSITY\n%e\n", inp.MassDensity)) + if err != nil { + return err + } + + // Assign material to all elements + _, err = f.WriteString("*SOLID SECTION,MATERIAL=resin,ELSET=eC3D4\n") + if err != nil { + return err + } + + // Assign material to all elements + _, err = f.WriteString("*SOLID SECTION,MATERIAL=resin,ELSET=eC3D10\n") + if err != nil { + return err + } + + // Assign material to all elements + _, err = f.WriteString("*SOLID SECTION,MATERIAL=resin,ELSET=eC3D8\n") + if err != nil { + return err + } + + // Assign material to all elements + _, err = f.WriteString("*SOLID SECTION,MATERIAL=resin,ELSET=eC3D20R\n") + if err != nil { + return err + } + + // Write analysis + + _, err = f.WriteString("*STEP\n*STATIC\n") + if err != nil { + return err + } + + // Write point loads. + + // Include a separate file to avoid cluttering the `inp` file. + _, err = f.WriteString(fmt.Sprintf("*INCLUDE,INPUT=%s\n", inp.PathLoad)) + if err != nil { + return err + } + + err = inp.writeLoad() + if err != nil { + return err + } + + if inp.GravityIsNeeded { + // Include a separate file to avoid cluttering the `inp` file. + _, err = f.WriteString(fmt.Sprintf("*INCLUDE,INPUT=%s\n", inp.PathGravity)) + if err != nil { + return err + } + + err = inp.writeGravity() + if err != nil { + return err + } + } + + // Pick element results. + + _, err = f.WriteString("*EL FILE\n") + if err != nil { + return err + } + + _, err = f.WriteString("S\n") + if err != nil { + return err + } + + // Pick node results. + + _, err = f.WriteString("*NODE FILE\n") + if err != nil { + return err + } + + _, err = f.WriteString("U\n") + if err != nil { + return err + } + + // Conclude. + + _, err = f.WriteString("*END STEP\n") + if err != nil { + return err + } + + return nil +} diff --git a/sdf/finiteelements/mesh/mesh.go b/sdf/finiteelements/mesh/mesh.go new file mode 100644 index 000000000..9115fd5a1 --- /dev/null +++ b/sdf/finiteelements/mesh/mesh.go @@ -0,0 +1,5 @@ +// Package mesh provides convenient types & functions for meshes consisting of finite elements. +// Like 4-node & 10-node tetrahedra, 8-node & 20-node hexahedra. +package mesh + +//----------------------------------------------------------------------------- diff --git a/sdf/finiteelements/mesh/restraint_load.go b/sdf/finiteelements/mesh/restraint_load.go new file mode 100644 index 000000000..57bb63677 --- /dev/null +++ b/sdf/finiteelements/mesh/restraint_load.go @@ -0,0 +1,58 @@ +package mesh + +import ( + v3 "github.com/deadsy/sdfx/vec/v3" + "github.com/deadsy/sdfx/vec/v3i" +) + +// A voxel is detected that intersects with the restraint location. +// ... or that is closest to the restraint location. +// Boundary is created for all the nodes inside that voxel. +// This way, the stress concentration at the restraint may be alleviated by distributing it. +type Restraint struct { + Location v3.Vec // Exact coordinates inside restraint. + IsFixedX bool // Is X degree of freedom fixed? + IsFixedY bool // Is Y degree of freedom fixed? + IsFixedZ bool // Is Z degree of freedom fixed? + voxel v3i.Vec // Intersecting voxel: to be computed by logic. +} + +// A voxel is detected that intersects with the point load location. +// ... or that is closest the point load location. +// A vertex inside that voxel is selected, i.e. the vertex closest to load location. +// The vertex is assigned the load. +// This way, point load is applied to a single node/vertex of the mesh. +type Load struct { + Location v3.Vec // Exact coordinates inside restraint. + Magnitude v3.Vec // X, Y, Z magnitude. + voxel v3i.Vec // Intersecting voxel: to be computed by logic. + nodeREF v3.Vec // Eventual vertex/node to which the load is applied. To be computed. +} + +func NewRestraint(location v3.Vec, isFixedX, isFixedY, isFixedZ bool) *Restraint { + return &Restraint{ + Location: location, + IsFixedX: isFixedX, + IsFixedY: isFixedY, + IsFixedZ: isFixedZ, + voxel: v3i.Vec{X: -1, Y: -1, Z: -1}, // We depend on -1 value to see if voxel is valid. + } +} + +// All the nodes inside the input voxel will be considered as restraint. +func NewRestraintByVoxel(voxel v3i.Vec, isFixedX, isFixedY, isFixedZ bool) *Restraint { + return &Restraint{ + IsFixedX: isFixedX, + IsFixedY: isFixedY, + IsFixedZ: isFixedZ, + voxel: voxel, + } +} + +func NewLoad(location v3.Vec, magnitude v3.Vec) *Load { + return &Load{ + Location: location, + Magnitude: magnitude, + voxel: v3i.Vec{X: -1, Y: -1, Z: -1}, // We depend on -1 value to see if voxel is valid. + } +}