# 4-node tetrahedron

This is a following question to this one: Filter bad finite elements

I focus only on 4-node tetrahedral element for now.

# Re-implementing FORTRAN in Golang

I studied some files of the source code like these:

• `ccx_2.20/src/e_c3d.f`
• `ccx_2.20/src/gauss.f`
• `ccx_2.20/src/shape4tet.f`.

I re-implemented this error in Golang:

nonpositive jacobian determinant

My first function is this:

``````
import (
)

func isBadTet4(coords 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

}
``````

And this is my second function that is called by the first function:

``````
// Reference:
// CCX source code:
// ccx_2.20/src/shape4tet.f
func isBadGaussTet4(coords v3.Vec, xi, et, ze float64) (bool, float64) {
// Coordinates of the nodes.
var xl float64

for i := 0; i < 4; i++ {
xl[i] = coords[i].X
xl[i] = coords[i].Y
xl[i] = coords[i].Z
}

// Shape functions.
var shp float64

shp = 1.0 - xi - et - ze
shp = xi
shp = et
shp = ze

// local derivatives of the shape functions: xi-derivative

shp = -1.0
shp = 1.0
shp = 0.0
shp = 0.0

// local derivatives of the shape functions: eta-derivative

shp = -1.0
shp = 0.0
shp = 1.0
shp = 0.0

// local derivatives of the shape functions: zeta-derivative

shp = -1.0
shp = 0.0
shp = 0.0
shp = 1.0

// computation of the local derivative of the global coordinates (xs)
xs := 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*(xs*xs-xs*xs) -
xs*(xs*xs-xs*xs) +
xs*(xs*xs-xs*xs)

// 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
}
``````

# Input

When the 4 nodes of a tetrahedron are like this:

``````import (
)

coords := v3.Vec{
{-3.2205676734447444, -0.01609302163123899, 6.644957738369699},
{-3.4127413321494826, -0.01609302163123899, 6.644957738369699},
{-3.2205676734447444, -0.235796231031416,  6.7280195367106765},
{-3.2205676734447444, -0.01609302163123899, 6.732770196393929},
}
``````

My Golang code doesn’t detect it as a bad element. But, CCX solver throws this error for it:

*ERROR in e_c3d: nonpositive jacobian
determinant in element

Investigations show that the bad tetrahedron has actually a fine shape. I visualized the bad element with OpenSCAD, it is just fine:

But the bad element is very small with a volume of `0.000617924113573828`.

# Question

Despite the fact that I have exactly followed the CCX source code to filter bad elements, I’m seeing CCX solver throwing errors of nonpositive jacobian determinant. I’m just looking for some hints and pointers to continue my investigations. Thanks.

# Test code

I have shared the Go code by this playground so that you can run and test it:

@Megidd
Just for your info, when I use your coordinates in the following dataset and run it in CCX_2.20, I don’t get any error at all, so I believe something must be wrong in your setup.

``````*Node
1,-3.2205676734447444, -0.01609302163123899, 6.644957738369699
2,-3.4127413321494826, -0.01609302163123899, 6.644957738369699
3,-3.2205676734447444, -0.235796231031416,  6.7280195367106765
4,-3.2205676734447444, -0.01609302163123899, 6.732770196393929
*Element, Type=C3D4
1, 1, 2, 3, 4
*Nset, Nset=FixedNodes
1, 2, 3
4
*Elset, Elset=Eall
1
*Material, Name=S185
*Density
7.8E-09
*Elastic
210000, 0.28
*Expansion, Zero=20
1.1E-05
*Conductivity
14
*Specific heat
440000000
*Solid section, Elset=Eall, Material=S185
*Step
*Static
*Boundary, op=New
*Boundary
FixedNodes, 1, 6, 0
*Node file
RF, U
*El file
S, E
*End step
``````
2 Likes

Thanks. I’m going to double check my setup.

Just for reference, my `inp` setup is shared here:

# Single element analysis

Inspired by your help, I modified my huge `inp` file. So that it only contains a single element:

``````*HEADING
Model: 3D model Date: 2023-Jun-10 UTC
*NODE
1, -3.2205676734447444, -0.01609302163123899, 6.644957738369699
2, -3.4127413321494826, -0.01609302163123899, 6.644957738369699
3, -3.2205676734447444, -0.235796231031416,  6.7280195367106765
4, -3.2205676734447444, -0.01609302163123899, 6.732770196393929
*ELEMENT, TYPE=C3D4, ELSET=Eall
1, 1, 2, 3, 4
*BOUNDARY
1, 1, 3
2, 1, 3
3, 1, 3
*MATERIAL, name=resin
*ELASTIC,TYPE=ISO
9.000000e+02,3.000000e-01,0
*DENSITY
1.250000e-09
*SOLID SECTION,MATERIAL=resin,ELSET=Eall
*ELSET, ELSET=ECheck, GENERATE
1, 330696
*STEP
*STATIC
Eall,GRAV,9810.,0.,0.,1.
*EL FILE
S
*NODE FILE
U
*END STEP
``````

# No error

The element analysis is just fine and it throws no error at all:

# Implication

### Element connectivity

One possibility is that my huge `inp` file contains tetrahedral elements that their connectivity to neighbor elements is problematic. I mean, the element itself might be fine, but for example it’s not properly connected to neighbor elements or boundary. Maybe # Question

Is it a feasible implication? Is there any convenient tool that lets me figure out and visualize such problems?

@Megidd
Running your dataset, I got this

``````*ERROR in e_c3d: nonpositive jacobian
determinant in element       80983
*ERROR in e_c3d: nonpositive jacobian
determinant in element       81013
*ERROR in e_c3d: nonpositive jacobian
determinant in element       39028
*ERROR in e_c3d: nonpositive jacobian
determinant in element       39261
*ERROR in e_c3d: nonpositive jacobian
determinant in element      145748
``````

Isolation the topology for these elements I got this

``````39028,8791,8791,8791,8791
39261,8832,8832,8832,8832
80983,16581,16579,16599,16581
81013,16593,16612,16595,16593
145748,26102,26103,26102,28375

``````

Four different coordinates are needed to make a proper tetrahedral element, so you have some errors in your topology routine

2 Likes

Thanks You found a bug 1 Like

That bug is fixed now 1 Like