Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

FE - Finite Elements #68

Open
wants to merge 798 commits into
base: master
Choose a base branch
from
Open

FE - Finite Elements #68

wants to merge 798 commits into from

Conversation

Megidd
Copy link
Contributor

@Megidd Megidd commented Feb 21, 2023

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

Screenshot_20230221_110954

FEA volume mesh

Screenshot_20230221_111047

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:

Marching Cubes

  • 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 🙄

@Megidd Megidd marked this pull request as draft February 21, 2023 08:18
@deadsy
Copy link
Owner

deadsy commented Feb 21, 2023

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
Copy link
Contributor Author

Megidd commented Feb 21, 2023

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 👍

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 ✔️

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 🙄

@deadsy
Copy link
Owner

deadsy commented Feb 22, 2023

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
Copy link
Contributor Author

Megidd commented Feb 22, 2023

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
Copy link
Contributor Author

Megidd commented Feb 22, 2023

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.

Screenshot_20230222_140608

10-node

From the same document:

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

Screenshot_20230222_140623

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
Copy link
Owner

deadsy commented Feb 22, 2023

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
Copy link
Contributor Author

Megidd commented Feb 23, 2023

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:

Screenshot_20230223_125445

Screenshot_20230223_125519

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

@Megidd
Copy link
Contributor Author

Megidd commented Feb 28, 2023

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
Copy link
Owner

deadsy commented Feb 28, 2023

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
Copy link
Contributor Author

Megidd commented Mar 2, 2023

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
Copy link
Owner

deadsy commented Mar 9, 2023

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
Copy link
Contributor Author

Megidd commented Mar 12, 2023

Test

4-node tetrahedra

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

8-node hexahedra

Screenshot_20230312_173136

Screenshot_20230312_173707

20-node hexahedra

Screenshot_20230312_155101

Screenshot_20230312_155454

@Megidd Megidd marked this pull request as ready for review March 12, 2023 14:11
@Megidd
Copy link
Contributor Author

Megidd commented Mar 12, 2023

Possible optimization

  • 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
Copy link
Contributor Author

Megidd commented Mar 12, 2023

Hex20 ✔️

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
Copy link
Contributor Author

Megidd commented Mar 14, 2023

AMR

Among the optimizations here: #68 (comment)

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 Megidd marked this pull request as draft March 30, 2023 08:35
@Megidd
Copy link
Contributor Author

Megidd commented May 13, 2023

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

@Megidd
Copy link
Contributor Author

Megidd commented Jun 15, 2023

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.

image

image

@Megidd
Copy link
Contributor Author

Megidd commented Jun 15, 2023

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
Copy link
Contributor Author

Megidd commented Jun 15, 2023

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
Copy link
Contributor Author

Megidd commented Jun 19, 2023

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.

2023-06-19 05_43_39-Window
2023-06-19 06_01_27-Window
2023-06-19 06_02_15-Window

Artifacat: only hex

Artifact when only hex elements are generated.

2023-06-19 06_10_46-Window
2023-06-19 06_11_31-Window

Artifact: only tet

Artifact when only tet elements are generated.

2023-06-19 06_12_01-Window
2023-06-19 06_12_48-Window

Reproduce

To reproduce the artifact, commit b7973b8 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
Copy link
Owner

deadsy commented Jun 21, 2023

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

@deadsy
Copy link
Owner

deadsy commented Jun 21, 2023

Is this code good for a commit to master?

@Megidd
Copy link
Contributor Author

Megidd commented Jun 22, 2023

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 Megidd marked this pull request as ready for review June 22, 2023 05:55
@Megidd
Copy link
Contributor Author

Megidd commented Jul 8, 2023

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

2023-07-08 11_35_43-FreeCAD 0 20 2
2023-07-08 11_35_07-FreeCAD 0 20 2
2023-07-08 11_36_11-FreeCAD 0 20 2

Not expected

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

2023-07-08 11_41_31-FreeCAD 0 20 2

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants