# Detect bad finite elements: 4-node tetrahedral

**URL:** <https://calculix.discourse.group/t/detect-bad-finite-elements-4-node-tetrahedral/1700>\
**Category:** Uncategorized\
**Created:** [June 10, 2023, 1:54pm UTC](https://calculix.discourse.group/t/detect-bad-finite-elements-4-node-tetrahedral/1700 "2023-06-10T13:54:09Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Megidd](https://yyz2.discourse-cdn.com/free1/user_avatar/calculix.discourse.group/megidd/32/1875_2.png) [@Megidd](https://calculix.discourse.group/u/Megidd)\
**Post date:** [June 10, 2023, 1:54pm UTC](https://calculix.discourse.group/t/detect-bad-finite-elements-4-node-tetrahedral/1700/1 "2023-06-10T13:54:09Z")

</div>

# 4-node tetrahedron

This is a following question to this one: [Filter bad finite elements](https://calculix.discourse.group/t/filter-bad-finite-elements/1689)

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

# Re-implementing FORTRAN in Golang

I studied some files of the [source code](http://www.dhondt.de/ccx_2.20.src.tar.bz2) 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:

```go

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:

```go

// 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:

```go
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:

 ![image](https://global.discourse-cdn.com/free1/uploads/calculix/original/2X/f/f9d8437c1d21bcc1876200e0786162dfed2be42b.png)

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:

> **[Go Playground - The Go Programming Language](https://go.dev/play/p/ZANLpu2lMUW)**

---

<div class="post-metadata">

**Author:** ![fgr](https://avatars.discourse-cdn.com/v4/letter/f/f4b2a3/32.png) [@fgr](https://calculix.discourse.group/u/fgr)\
**Post date:** [June 10, 2023, 5:10pm UTC](https://calculix.discourse.group/t/detect-bad-finite-elements-4-node-tetrahedral/1700/2 "2023-06-10T17:10:46Z")

</div>

@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.

```auto
*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
*Nset, Nset=LoadNode
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
*Cload, op=New
*Dload, op=New
*Cload
LoadNode, 3, -1
*Node file
RF, U
*El file
S, E
*End step

```

---

<div class="post-metadata">

**Author:** ![Megidd](https://yyz2.discourse-cdn.com/free1/user_avatar/calculix.discourse.group/megidd/32/1875_2.png) [@Megidd](https://calculix.discourse.group/u/Megidd)\
**Post date:** [June 11, 2023, 2:36am UTC](https://calculix.discourse.group/t/detect-bad-finite-elements-4-node-tetrahedral/1700/3 "2023-06-11T02:36:49Z")

</div>

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

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

> **[teapot-tet4.inp.zip](https://drive.google.com/file/d/1vG54mH0aFONJILctwt-jE--sGhjj-LBY/view?usp=sharing)**
>
> Google Drive file.

---

<div class="post-metadata">

**Author:** ![Megidd](https://yyz2.discourse-cdn.com/free1/user_avatar/calculix.discourse.group/megidd/32/1875_2.png) [@Megidd](https://calculix.discourse.group/u/Megidd)\
**Post date:** [June 11, 2023, 5:34am UTC](https://calculix.discourse.group/t/detect-bad-finite-elements-4-node-tetrahedral/1700/4 "2023-06-11T05:34:04Z")

</div>

# Single element analysis

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

```inp
*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
*DLOAD
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:

 ![Screenshot_20230611_084906](https://global.discourse-cdn.com/free1/uploads/calculix/original/2X/b/b84ced9bc70e49885b745a51fca7fb9f13fc6969.png)

# 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?

---

<div class="post-metadata">

**Author:** ![fgr](https://avatars.discourse-cdn.com/v4/letter/f/f4b2a3/32.png) [@fgr](https://calculix.discourse.group/u/fgr)\
**Post date:** [June 11, 2023, 6:29am UTC](https://calculix.discourse.group/t/detect-bad-finite-elements-4-node-tetrahedral/1700/5 "2023-06-11T06:29:27Z")

</div>

@Megidd  
Running your dataset, I got this

```auto
*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

```auto
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

---

<div class="post-metadata">

**Author:** ![Megidd](https://yyz2.discourse-cdn.com/free1/user_avatar/calculix.discourse.group/megidd/32/1875_2.png) [@Megidd](https://calculix.discourse.group/u/Megidd)\
**Post date:** [June 11, 2023, 7:02am UTC](https://calculix.discourse.group/t/detect-bad-finite-elements-4-node-tetrahedral/1700/6 "2023-06-11T07:02:45Z")

</div>

Thanks 🙂 You found a bug 🐞

---

<div class="post-metadata">

**Author:** ![Megidd](https://yyz2.discourse-cdn.com/free1/user_avatar/calculix.discourse.group/megidd/32/1875_2.png) [@Megidd](https://calculix.discourse.group/u/Megidd)\
**Post date:** [June 11, 2023, 9:08am UTC](https://calculix.discourse.group/t/detect-bad-finite-elements-4-node-tetrahedral/1700/7 "2023-06-11T09:08:43Z")

</div>

That bug is fixed now ✅