Node orderings after refinement bug

[ackirby88]

Our team is having issues with node orderings after refinement of tets, prisms, and pyramids (hexes are fine). Below, I pasted some code that can be placed in the t8_step2_uniform_forest.cxx tutorial (paste after the forest construction and refinement) for demonstration:

  // paste at ~line 167 of t8_step2_uniform_forest.cxx
  int negvol_counts[T8_ECLASS_COUNT] = {0};
  double vertex_coords[3];
  double coords[3];
  int icorner;

  t8_locidx_t num_local_trees = t8_forest_get_num_local_trees(forest);
  for (t8_locidx_t itree = 0,elem_index = 0; itree < num_local_trees; ++itree) {
    t8_eclass_t tree_class = t8_forest_get_tree_class(forest,itree);
    t8_locidx_t num_elements_in_tree = t8_forest_get_tree_num_elements(forest,itree);
    t8_eclass_scheme_c *eclass_scheme = t8_forest_get_eclass_scheme(forest,tree_class);
    t8_gloidx_t gtreeid = t8_forest_global_tree_id(forest,itree);
    t8_cmesh_t cmesh = t8_forest_get_cmesh(forest);

    /* loop over all local elements in the tree. */
    for (t8_locidx_t ielement = 0; ielement < num_elements_in_tree; ++ielement, ++elem_index) {
      t8_element_t *element = t8_forest_get_element_in_tree(forest,itree,ielement);
      const t8_element_shape_t element_shape = eclass_scheme->t8_element_shape(element);
      const int num_corners = t8_eclass_num_vertices[element_shape];

      /* insert vertices */
      double tree_vertices[T8_ECLASS_MAX_CORNERS * 3];
      for (icorner = 0; icorner < num_corners; icorner++) {
        eclass_scheme->t8_element_vertex_reference_coords(element, icorner, vertex_coords);
        t8_geometry_evaluate(cmesh,gtreeid,vertex_coords,1,coords);
        memcpy(&tree_vertices[3*icorner],coords,3*sizeof(double));
      }
      if(t8_cmesh_tree_vertices_negative_volume(element_shape,tree_vertices,num_corners)){
        negvol_counts[element_shape]++;
      }
    }
  }

  printf("NEGATIVE VOLUME COUNTS: tets %d, prism %d, pyramid %d, hex %d\n",
          negvol_counts[T8_ECLASS_TET],
          negvol_counts[T8_ECLASS_PRISM],
          negvol_counts[T8_ECLASS_PYRAMID],
          negvol_counts[T8_ECLASS_HEX]);

Using the default case with prisms, refining three levels, I get the following result:
NEGATIVE VOLUME COUNTS: tets 0, prism 448, pyramid 0, hex 0.
I get similar results using just one level of refinement and for tet and pyramid element types on different cmeshes.

My question is:
(a) is there bug in the assembly of the node orderings here in my example, or
(b) is there a bug in element node ordering after refinement?

[ackirby88]

Also, you can test with the other element types by adjusting the cmesh construction:

static t8_cmesh_t
t8_step2_build_prismcube_coarse_mesh (sc_MPI_Comm comm, t8_eclass_t eclass)
{
  t8_cmesh_t cmesh;

  /* Build a coarse mesh of 2 prism trees that form a cube. */
  cmesh = t8_cmesh_new_hypercube (eclass, comm, 0, 0, 0);
  t8_global_productionf (" [step2] Constructed coarse mesh with 2 prism trees.\n");

  return cmesh;
}

Then, e.g., cmesh = t8_step2_build_prismcube_coarse_mesh (comm, T8_ECLASS_HEX);

[ackirby88]

Testing this case for each element type, I found the following incorrect children node orderings on refinement:
tetrahedra elements: children elements #3 and #5
prism elements: children elements #2 and #6
pyramid elements: children elements #1, #6, and #8
hexahedral elements: no bad children elements
Note: children ID numbers here are base 0.

For the t8_cmesh_new_hypercube cmesh construction with multiple trees, the incorrect children are consistent in each tree refined, e.g., tree 1's bad children for tetrahedral refinement are global element id #11 (local element #3) and global element id #13 (local element #5).

Example: Following the above example with Tets, children 3 and 5 have the following node coordinates

ELEM 3 is NEGATIVE
         X                   Y                   Z
Node 0   0.500000000000000   0.000000000000000   0.000000000000000
Node 1   0.500000000000000   0.000000000000000   0.500000000000000
Node 2   1.000000000000000   0.000000000000000   0.500000000000000
Node 3   1.000000000000000   0.500000000000000   0.500000000000000

ELEM 5 is NEGATIVE
         X                   Y                   Z
Node 0   0.500000000000000   0.000000000000000   0.500000000000000
Node 1   1.000000000000000   0.000000000000000   0.500000000000000
Node 2   1.000000000000000   0.500000000000000   0.500000000000000
Node 3   1.000000000000000   0.500000000000000   1.000000000000000

Drawing these points out on paper show these nodes follow a right-hand-rule, which is incorrect relative to t8code's TET node ordering (https://github.com/DLR-AMR/t8code/wiki/Build-Cmesh).

[Davknapp]

Hey @ackirby88 , thank you for submitting this issue! Let me first tackle the volume issue. This test is meant to be used to check the volume of a coarse mesh element. If I understand your code correctly, you compute the coordinates of the nodes of a refined element in a forest and pass them to the function 't8_cmesh_tree_vertices_negative_volume'.

This function does not compute the actual volume of the element, but checks that the order of the nodes will give a positive volume. For hexes and pyramids this is done by checking if the fourth node is above the plane given by nodes 0-3. For Tets and prisms the check is similar, but we check for the third node and the plane given by nodes 0-2. (check the comment here).

This assumptions holds for tree / coarse mesh elements. After refining the forest given by the coarse mesh, elements of different types arise. For prisms for example we have two types, defined by their triangular base. The type of the root of the tree (the coarse mesh element) is 0. All children in the corner have type 0, too. The prism-child in the middle has type 1. For all children of type 1 the function 't8_cmesh_tree_vertices_negative_volume' will fail, because its computation is based on the vectors between node_0-node_1 and node_0-node_2 and does not follow the node-ordering you described correctly.

base of a         base of a
type 0 prism     type 1 prism
    2              1--2
  / |              | /
 /  |              |/
0 --1              0

This order is given by the triangles, their node order has been defined here: https://arxiv.org/pdf/1509.04627.pdf, see section 2.1, defintion 1 (2a) This also the reason why you detect 448 prisms with negative volume on level 3 (per prism, there are 224 prisms of type 1).

All the child-ids of the different element classes given by you have in common that they refer to children in the "middle" of their parent, for which the node-ordering does not follow the rule. We will discuss, if these children should also follow the right-hand-rule and will keep you up to date!

If you want to compute the volume of an element, you can use t8_forest_element_volume in the meantime

[ackirby88]

@Davknapp This is really helpful!

Originally, I was reassembling the element-node connectivities by following the VTK plotting function, which I assumed the element vertex iteration orderings were consistent. But now that is not the case, I should be able to look at the t8_forest_element_volume function to reconstruct them correctly.

My guess is more users will get tripped up by this in the future. On our side, we are offloading all of the geometry related information to the application solver to be solver-consistent, following the p4est modus operandi. I would recommend either following the right-hand-rule for children or providing a function to the user to get a right-hand-rule consistent iteration of the vertices.

Thanks, again, for the very insightful and helpful information!

[ackirby88]

@Davknapp For tets, (following https://arxiv.org/pdf/1509.04627.pdf, section 2.1, definition 1 (2b)), child #3 - T3 is supposed to be similar to the parent (same node orientation). Why is it then that we get a negative volume using the cmesh function if the parent element corresponds to a cmesh tree (which follows the t8code node ordering correctly)?

[Davknapp]

@ackirby88 First of all, thank you for your tipps on the user experience. It helps a lot to improve our code!

On the subject: I guess that the child ids you mention here are given in morton-order. We do not use this order in definition 1 (2b). This numbers are used to help to define the subsimplices. Even though the corner-elements are similar to the parent, they are not traversered first. Similar to a triangle, some of the middle tetrahedra have a lower morton-index than the corner elements (see Fig 5 in the same paper, page 10). For all tetrahedra the children are traversered similar to the picture:

The third tetrahedra is here:

and is a tetrahedron of type 5. You can check table 2 in the same paper to see the relation between local-ids, and the type of their parents. If I am correct the node order of the third child (with respect to the SFC/morton order) is not similar to its parent and we therefore still get a negative volume using the cmesh function. Please feel free to ask further questions if something has been left unclear or if you observe further/other bugs