FE - Finite Elements

[Megidd]

Objective

Render an SDF3 to finite elements in the shape of tetrahedra.

Philosohpy / why

There are some commercial software to convert STL surface mesh to FEA volume mesh. Like WELSIM. But as far as I tested, they are unreliable, limited, slow, ...

Problem

I have been only able to do STL mesh to FEA mesh conversion for WELSIM's own sample STL files like the following hinge.stl sample. For any other STL file, the conversion is throwing errors. Looks like even the commercial apps only work for perfect STL meshes without any fault. Not sure, maybe?

STL surface mesh

FEA volume mesh

Methodology

This looks like a feasible approach:

STL triangles -> SDF3 -> MT (marching tetrahedra) -> FE (finite elements) -> ABAQUS or CalculiX input file

Marching cubes algorithm and marching tetrahedra are closely related.

It's needed to develop a code like below. But for marching tetrahedra rather than marching cubes:

https://github.com/deadsy/sdfx/blob/5bd4805907429608b176ed6ddabe7063e48e87f3/render/march3.go#L4

• The above code extracts a triangle isosurface with all the 0 values.
• We need a code to extract tetrahedra elements with all the non-positive <=0 values.
  • Non-positive means 3D space on and inside the surface geometry.

Output API

Feeding the finite elements in the shape of tetrahedra - or any other shape - into the FEA engines could be done in various ways. But looks like the ABAQUS way makes sense. Popular opensource FEA engines like CalculiX are using it too.

...naming conventions and input style formats for CalculiX are based on those used by ABAQUS...

Question

Does the whole thing make sense? It looks like a solution to me :roll_eyes:

[deadsy]

So- if I understand this, the idea is to generate a set of tetrahedra who collectively make up the volume of the model. If you want to use it tor FEA I assume the tetraheda should be small and roughly the same size.

There's a terminology issue that needs to be addressed:

Marching cubes generates surface triangles. Marching cubes could also generate tetrahedra throughout the volume of the model (I think).

That is: the "cubes" is a reference to the sampling of the field - not to the result. Likewise I think "marching tetrahedra" is a reference to the sampling technique.

So: sample the field with "marching X" and generate tetrahedra instead of surface triangles. That is: within the volume of the model (all samples <=0) you just generate N tetrahedra that fill the sample volume (cube or tetrahedron).

Works for 2D as well. ie: generate 2D triangle instead of lines.

[Megidd]

So- if I understand this, the idea is to generate a set of tetrahedra who collectively make up the volume of the model.

@deadsy Yes, that's right :+1:

If you want to use it tor FEA I assume the tetraheda should be small and roughly the same size.

I believe it depends on the desired accuracy and the complexity of the shape. To achieve the same accuracy throughout a 3D model, simple regions may need fewer finite elements. Complex regions may need more. As the initial implementation, the same size for all finite elements might be a good start :heavy_check_mark:

There's a terminology issue that needs to be addressed...

Right. Great point. I'm still not sure which sampling is more suitable. It is mentioned that:

Marching tetrahedra ... clarifies a minor ambiguity problem of the marching cubes algorithm with some cube configurations.

Maybe the above sentence means marching tetrahedra is the preferred sampling. I'm not sure :roll_eyes:

[deadsy]

Marching tetrahedra ... clarifies a minor ambiguity

You need to take a punt at saddle points and you may get it wrong. In practice it doesn't seem to be too much of an issue. At least I haven't had to deal with it.

See case 10.... https://en.wikipedia.org/wiki/Marching_squares

If you want to keep things simple you could just adapt marching cubes/squares to generate tetrahedra/triangles.

I haven't looked at your code yet...

[Megidd]

To be able to easily debug the code, I guess I'm going to implement the output API first, then work on the sampling algorithm. To develop the output API, some constant hard-coded tetrahedra vertices may be used.

[Megidd]

4-node tetrahedral element

CalculiX solver documentation mentions:

Four-node tetrahedral element (C3D4) ... This element is included for completeness, however, it is not suited for structural calculations unless a lot of them are used (the element is too stiff). Please use the 10-node tetrahedral element instead.

10-node

From the same document:

Ten-node tetrahedral element (C3D10) ... The element behaves very well and is a good general purpose element...

Pick 4-node

Let's keep it simple. Let's consider 4-node tetrahedral element (C3D4) for the initial implementation. Later, it can be expanded to 10-node.

[deadsy]

TIL... about FEA nodes. I still don't understand why adding edge nodes to the TET that appear to be linear combinations of the corner nodes makes them better.

[Megidd]

TIL... about FEA nodes. I still don't understand why adding edge nodes to the TET that appear to be linear combinations of the corner nodes makes them better.

CalculiX solver documentation mentions:

As far as I know, stress and strain are distributed linear for one and quadratic for the other :slightly_smiling_face:

[Megidd]

Trying to figure out how to adapt the marching cubes algorithm: https://math.stackexchange.com/q/4648068/197913

UPDATE:

About the tools to do so: https://softwarerecs.stackexchange.com/q/86486/12269

[deadsy]

Conceptually it seems simple enough. For each of those 256 cases work, break up the defined volume into tetrahedra. But yes- you don't want to do it by hand. fwiw- doing the 2d form of the same thing can lead to insights.

[Megidd]

MC tables

Tables of marching cubes algorithm are better understood now.

Edge table

To determine which edges are involved with triangle creation:

// 8 vertices -> 256 possible inside/outside combinations
// A 1 bit in the value indicates an edge with a line end point.
// 12 edges -> 12 bit values, note the fwd/rev symmetry
var mcEdgeTable = [256]int{ /* ... */ }

// Example: mcEdgeTable[3] is:
//
// 0x030a
//
// It means:
// 0b 0000 0011 0000 1010
// It means edges are:
// 9, 8, 3, 1

// Min possible:
// 0b 0000 0000 0000 0000 i.e. 0x0000
// Max possible:
// 0b 0000 1111 1111 1111 i.e. 0x0fff

Pair table

To define how 8 cube nodes are related to 12 cube edges.

// These are the vertex pairs for the edges
var mcPairTable = [12][2]int{
    {0, 1}, // Edge 0
    {1, 2}, // Edge 1
    {2, 3}, // Edge 2
    {3, 0}, // Edge 3
    {4, 5}, // Edge 4
    {5, 6}, // Edge 5
    {6, 7}, // Edge 6
    {7, 4}, // Edge 7
    {0, 4}, // Edge 8
    {1, 5}, // Edge 9
    {2, 6}, // Edge 10
    {3, 7}, // Edge 11
}

Triangle table

To determine the edge order for each triangle being created. Each triangle is created by connecting 3 edges.

// specify the edges used to create the triangle(s)
var mcTriangleTable = [256][]int{ /* ... */ }

// Example: mcTriangleTable[3] is:
//
// {1, 8, 3, 9, 8, 1},
//
// It means:
// 1 of 2 triangles connecting edge 1 to edge 8 to edge 3
// 1 of 2 triangles connecting edge 9 to edge 8 to edge 1

Notes

Repeated edge

An edge can be repeated multiple times by one of the 256 cases on mcTriangleTable.

However, each edge can have only a single point on it. An edge cannot have multiple points on it. Because each edge can only have a single point with a 0 value.

[deadsy]

btw- the rest of the code has been written to follow the strictures of golint. I realise this is deprecated, but it is a standard.

[Megidd]

Test

4-node tetrahedra

Implementation of Tet4 FE is not finished. Function mcToTet4 has a related TODO.

8-node hexahedra

20-node hexahedra

[Megidd]

Possible optimization

• [x] Use faster math solvers like PaStiX or PARDISO rather than slow default SPOOLES.
  1. It requires building CCX with those math dependencies.
  2. You can download here an executable with the Intel Pardiso solver for x86_64 Linux systems. The executable has all the libraries statically linked to it. So it should run by itself without any dependency. Thanks to these guys.
  3. Actually I realized that the official executables are built with PaStiX which is the fastest solver available at the moment.
• [ ] Adaptive mesh refinement or AMR.
• [ ] ...

[Megidd]

Hex20 :heavy_check_mark:

CalculiX solver documentation mentions about 20-node hexahedral element with reduced integration:

The C3D20R element ... The element behaves very well and is an excellent general purpose element (if you are setting off for a long journey and you are allowed to take only one element type with you, that’s the one to take).

This PR implements this finite element. Hexahedral elements would define a pixelated/pixeled surface. But maybe that's good enough.

[Megidd]

AMR

Among the optimizations here: https://github.com/deadsy/sdfx/pull/68#issuecomment-1465210911

Adaptive mesh refinement AMR may have the highest priority. It would be helpful regardless of the FEA engine.

Octree space subdivision

File render/march3x.go implements an octree space subdivision to convert an SDF3 to a triangle mesh by the marching cubes algorithm. Can it be leveraged to implement a kind of AMR?

[Megidd]

Those 256 cases from 0 to 255 are being worked out here: https://github.com/Megidd/tetrahedron-table :wrench:

[Megidd]

Tetrahedron

The tetrahedron element is implemented using this repository: https://github.com/Megidd/tetrahedron-table

Now both hexahedron and tetrahedron can be used to represent a 3D mode.

Test

An STL file of WELSIM software i.e. hinge.stl is tested. It is modeled with:

• Combination of hex8 + tet4 finite elements.
• Combination of hex20 + tet10 finite elements.

[Megidd]

Good STL, good FEA

The last comment tested a good STL file of WELSIM software: hinge.stl

This STL had no faults, so the FEA had no problem too.

Faulty STL, faulty FEA

When input STL model has faults, the output FEA model can be difficult for math solvers. Throwing an error like this:

*ERROR in e_c3d: nonpositive jacobian determinant in element

Or this error:

matrix found to be singular

Repair

A repair stage may be required. I mean we may need a repair stage after finite elements are generated. The repair main objectives might be:

• To connect the disconnected or improperly connected elements.
• To replace the elements with non-positive Jacobian determinant
• ...

[Megidd]

Voxel

We are currently storing the finite elements by their containing voxel. Every voxel is corresponding to a cube of the marching cubes algorithm.

Voxel & repair

A repair logic requires knowing the neighborhood of elements. Our voxel data structure may help with repair logic.

[Megidd]

Artifact

Inside the 3D model some artifacts are observed. The artifact starts at some layers and it grows layer by layer. I ran out if ideas. Is the artifact a result of the marching cubes algorithm's famous ambiguity? 🙄

Artifact: hex + tet

Artifact when a combination of hex an tet elements are generated.

Artifacat: only hex

Artifact when only hex elements are generated.

Artifact: only tet

Artifact when only tet elements are generated.

Reproduce

To reproduce the artifact, commit b7973b8bbda55e5671ceef8c631a4efcb7796b31 sets the layer where the artifact begins. Then, just like this commit, layer can be incremented to see how artifact grows.

FreeCAD

The FreeCAD software can be used to open the ouput *.inp files.

[deadsy]

The marching cubes issue shows up when you have saddle points. These look like plain surfaces. Therefore: A bug?

[deadsy]

Is this code good for a commit to master?

[Megidd]

Is this code good for a commit to master?

If you ask me, I'd say let's do a code review and merge, if it's fine. Of course, renaming, formatting, cleaning up, and more could be done with golint and other tools.

The marching cubes issue shows up when you have saddle points. These look like plain surfaces. Therefore: A bug?

Right, it looks like a bug 🐞 I try to pinpoint the bug by stepping through the code with the debugger. Maybe I can fix it by submitting other PRs.

[Megidd]

4 FEA benchmarks

I have picked four benchmark 3D models from here to verify the FEA results. The square, circular, and I-shaped 3D beams are OK, but the pipe one is a bit strange.

As expected

Not expected

Why are some elements missing from the pipe? 🙄 Is it due to the same bug that was mentioned by the last comment? 🐞

[Megidd]

Previously, it was mentioned that our pipe benchmark is not as expected: https://github.com/deadsy/sdfx/pull/68#issuecomment-1626907610

Further investigation indicates that it's due to resolution sensitivity.

Resolution sensitivity

When resolution is bad

Our pipe benchmark is bad with this resolution:

voxel counts of marching algorithm are: (51 x 8 x 8)

When resolution is fine

It is fixed by just incrementing the resolution:

voxel counts of marching algorithm are: (52 x 9 x 9)

[Megidd]

However, playing around with resolution value couldn't get rid of the artifacts here yet: https://github.com/deadsy/sdfx/pull/68#issuecomment-1597156034