Error when loading Tetgen tet10 meshes

[Edward-Eth]

I have encountered a bug that I believe is occurring in libmesh (when used with MOOSE). When I load a tetgen tet10 mesh, the node orders are not read correctly leading to element self intersection:

I have made a tetgen file containing a single tet10 element, drawing the nodes out on paper shows that they're in the right places, and the nodes visible in paraview are in the correct place, it just seems that the order expected by libmesh and the order supplied by tetgen differ.
The ele file as generated by tetgen:

ele file contents

6  10  0
    1       3     8     1     4      9    10    11    12    13    14
# Generated by tetgen -p ModelGeneration/cube.poly -o2 

The node file as generated by tetgen:

node file contents

27  3  0  0
   1    0  0  1
   3    1  1  1
   4    0  1  1
   8    0  1  0
   9    0  0.5  1
  10    0.5  1  1
  11    0.5  1  0.5
  12    0  0.5  0.5
  13    0  1  0.5
  14    0.5  0.5  1
# Generated by tetgen -p ModelGeneration/cube.poly -o2 

Not sure if this is any help, but this is the order the nodes are written in the tetgen file.

I have attempted a few "logical" feeling orders to rewrite these, but given there are several million permutations of 10 values I'm not holding out much hope for this scattershot approach...

[jwpeterson]

In tetgen_io.C, it looks like we just assume there is an identity mapping between the Tetgen numbering and our Tet10 numbering, but apparently that's not correct...

libmesh/src/mesh/tetgen_io.C

Lines 209 to 213 in db05ce9

	// Vector that assigns element nodes to their correct position.
	// TetGen is normally 0-based
	// (right now this is strictly not necessary since it is the identity map,
	// but in the future TetGen could change their numbering scheme.)
	static const unsigned int assign_elm_nodes[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9};

Our Tet10 documentation shows the element numbering that we use:

libmesh/include/geom/cell_tet10.h

Lines 30 to 63 in db05ce9

	* The \p Tet10 is an element in 3D composed of 10 nodes. It is
	* numbered like this:
	*
	* \verbatim
	* 3
	* TET10: o
	* /|\
	* / | \
	* 7 / | \9
	* o | o zeta
	* / |8 \ ^
	* / o \ |
	* / 6 | \ |
	* 0 o.....o.|.......o 2 o---> eta
	* \ | / \
	* \ | / \
	* \ | / xi (out of page)
	* 4 o | o 5
	* \ | /
	* \ | /
	* \|/
	* o
	* 1
	* \endverbatim
	*
	* (xi, eta, zeta): { 0 <= xi <= 1
	* { 0 <= eta <= 1
	* { 0 <= zeta <= 1
	* { xi + eta + zeta <= 1
	*
	* \author Benjamin S. Kirk
	* \date 2002
	* \brief A 3D tetrahedral element with 10 nodes.
	*/

and to me it looks like the main difference is that we number the mid-edge nodes between the {0,1,2} face of the Tet first whereas Tetgen numbers the mid-edge nodes in a different way that isn't immediately obvious to me. It should be possible to directly set up the correct assign_elm_nodes array given these two numberings.

[Edward-Eth]

Figured out the mapping:
If a tetgen element is described in the ele file using ten nodes:
0 1 2 3 4 5 6 7 8 9
then reordering them in the file to:
0 1 2 3 6 7 9 5 8 4
gives me the correct element. I'm not sure how that would translate to your assign_elm_nodes system, (maybe literally just like that?) but I hope that helps! In the mean time I'm just going to write a script to reorder each line of my tetgen ele file into that format and hopefully that should work.

[jwpeterson]

We just use assign_elm_nodes to map like this:

// Assign node to element
elem->set_node(assign_elm_nodes[j]) =
  mesh.node_ptr(/*libmesh global node id*/);

So based on (a 0-based version of) your screenshot above I would say

assign_elm_nodes = [0, 1, 2, 3, 9, 7, 4, 5, 8, 6]

Should do the trick. For example, this would map (0-based) Tetgen node 4 to libmesh node 9. To check this, observe that (0-based) Tetgen node 4 and libmesh node 9 are situated between (0-based) vertices {2,3} in both diagrams.

Would it be possible for you to provide a non-reordered Tetgen input file that we could use as the basis of a unit test for this change in libmesh?

[Edward-Eth]

I redrew my diagram in the same orientation as the documentation diagram in cell_tet10.h and with zero indexing.

Using this, and testing with a script to rewrite the order I have confirmed that the mapping is definitely:
[0, 1, 2, 3, 6, 7, 9, 5, 8, 4]
As for a unit test, you need the .ele file:

1  10  0
    1       1     2     3     4      5    6    7    8    9    10
# Generated by tetgen -p ModelGeneration/pyramid.poly -o2

and the .node file:

10  3  0  0
  1    1  1  1
  2    0  1  0
  3    0  0  1
  4    0  1  1
  5    0  0.5  1
  6    0.5  1  1
  7    0.5  1  0.5
  8    0  0.5  0.5
  9    0  1  0.5
  10    0.5  0.5  1
# Generated by tetgen -p ModelGeneration/pyramid.poly -o2 

this should create a pyramid with corners at (1,1,1) (0,0,1) (0,1,0) and (0,1,1) and a total volume of 0.166... Does that work as something you can implement? If the ordering is off then the volume would be.

[jwpeterson]

Using this, and testing with a script to rewrite the order I have confirmed that the mapping is definitely:
[0, 1, 2, 3, 6, 7, 9, 5, 8, 4]

Yes, this is the "inverse" to the assign_elm_nodes map that I posted in my previous comment, in other words,

inverse_map[assign_elm_nodes[j]] = j

for each j, so we are in agreement.

Does that work as something you can implement?

Yes, I think so. I will hopefully make a PR some time this weekend that fixes the bug and adds a unit test. Thanks!

[jwpeterson]

#3338 was merged including the test case attached here.