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.fccx_2.20/src/gauss.fccx_2.20/src/shape4tet.f.
I re-implemented this error in Golang:
nonpositive jacobian determinant
My first function is this:
import (
v3 "github.com/deadsy/sdfx/vec/v3"
)
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)
}
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 [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
}
Input
When the 4 nodes of a tetrahedron are like this:
import (
"github.com/deadsy/sdfx/vec/v3"
)
coords := [4]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
Investigating bad 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:
You found a bug 
