Merge pull request #13765 from drwells/signed-tet-volume

Remove std::abs() from the tet volume calculation.
dealii · May 22, 2022 · cb36d2c · cb36d2c
2 parents 30f04b7 + 15f6645
commit cb36d2c
Show file tree

Hide file tree

Showing 5 changed files with 184 additions and 20 deletions.
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
@@ -871,16 +871,8 @@ namespace GridTools
     const std::vector<Point<spacedim>> &all_vertices,
     std::vector<CellData<dim>> &        cells)
   {
-    // This function only works for quads and hexes,
-    // and triangles. From the examples we tested,
-    // tetrahedron meshes loaded from Gmsh don't have
-    // negative measures, so we skip them for now.
-
-    // Since we use std::abs() in GridTools::cell_measure() to
-    // compute the measures of tetrahdra, the measures
-    // are always positive. But if the orientation is wrong,
-    // we may get negative determinants of the Jacobians in
-    // MappingFE, which will triger the assert there.
+    // This function is presently only implemented for hypercube and simplex
+    // volumetric (codimension 0) elements.
 
     if (dim == 1)
       return 0;
@@ -919,15 +911,11 @@ namespace GridTools
                        .is_simplex())
               {
                 ++n_negative_cells;
-                if (dim == 2)
-                  {
-                    // Triangular mesh, swap any two vertices
-                    std::swap(cell.vertices[1], cell.vertices[2]);
-                  }
-                else
-                  {
-                    AssertThrow(false, ExcNotImplemented());
-                  }
+                // By basic rules for computing determinants we can just swap
+                // two vertices to fix a negative volume. Arbitrarily pick the
+                // last two.
+                std::swap(cell.vertices[n_vertices - 2],
+                          cell.vertices[n_vertices - 1]);
               }
             else
               {

diff --git a/source/grid/grid_tools_nontemplates.cc b/source/grid/grid_tools_nontemplates.cc
@@ -128,7 +128,7 @@ namespace GridTools
         const auto &c = all_vertices[vertex_indices[2]];
         const auto &d = all_vertices[vertex_indices[3]];
 
-        return (1.0 / 6.0) * std::abs((a - d) * cross_product_3d(b - d, c - d));
+        return (1.0 / 6.0) * (d - a) * cross_product_3d(b - a, c - a);
       }
     else if (vertex_indices.size() == 5) // pyramid
       {

diff --git a/tests/simplex/grid_in_msh.cc b/tests/simplex/grid_in_msh.cc
@@ -70,4 +70,9 @@ main()
   deallog.push("tetrahedral_elements_dim3_spacedim3: ");
   check_file<3, 3>(std::string(SOURCE_DIR "/grid_in_msh/tet.msh"));
   deallog.pop();
+
+  // Make sure things work when half of the tets are twisted too:
+  deallog.push("tetrahedral_elements_dim3_spacedim3_twisted: ");
+  check_file<3, 3>(std::string(SOURCE_DIR "/grid_in_msh/tet_half_twisted.msh"));
+  deallog.pop();
 }
diff --git a/tests/simplex/grid_in_msh.output b/tests/simplex/grid_in_msh.output
@@ -133,3 +133,67 @@ LOOKUP_TABLE default
 
 
 DEAL:tetrahedral_elements_dim3_spacedim3: ::OK!
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 14 double
+0.00000 0.00000 1.00000
+0.00000 0.00000 0.00000
+0.00000 1.00000 1.00000
+0.00000 1.00000 0.00000
+1.00000 0.00000 1.00000
+1.00000 0.00000 0.00000
+1.00000 1.00000 1.00000
+1.00000 1.00000 0.00000
+0.00000 0.500000 0.500000
+1.00000 0.500000 0.500000
+0.500000 0.00000 0.500000
+0.500000 1.00000 0.500000
+0.500000 0.500000 0.00000
+0.500000 0.500000 1.00000
+
+CELLS 24 120
+4 13 10 11 9
+4 9 10 11 12
+4 11 8 10 13
+4 8 11 10 12
+4 4 13 6 9
+4 3 8 11 2
+4 1 8 0 10
+4 0 8 2 13
+4 4 10 0 13
+4 9 4 10 5
+4 11 2 13 6
+4 9 6 7 11
+4 1 3 8 12
+4 12 3 11 7
+4 7 12 5 9
+4 12 1 5 10
+4 13 4 10 9
+4 13 11 8 2
+4 13 6 9 11
+4 8 10 13 0
+4 10 1 8 12
+4 12 8 11 3
+4 9 7 12 11
+4 12 10 5 9
+
+CELL_TYPES 24
+10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 
+
+
+
+CELL_DATA 24
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 
+
+
+DEAL:tetrahedral_elements_dim3_spacedim3_twisted: ::OK!
diff --git a/tests/simplex/grid_in_msh/tet_half_twisted.msh b/tests/simplex/grid_in_msh/tet_half_twisted.msh
@@ -0,0 +1,107 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Entities
+8 12 6 1
+1 0 0 1 0 
+2 0 0 0 0 
+3 0 1 1 0 
+4 0 1 0 0 
+5 1 0 1 0 
+6 1 0 0 0 
+7 1 1 1 0 
+8 1 1 0 0 
+1 -1e-07 -1e-07 -9.999999994736442e-08 1e-07 1e-07 1.0000001 0 2 2 -1 
+2 -1e-07 -9.999999994736442e-08 0.9999999000000001 1e-07 1.0000001 1.0000001 0 2 1 -3 
+3 -1e-07 0.9999999000000001 -9.999999994736442e-08 1e-07 1.0000001 1.0000001 0 2 4 -3 
+4 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 1.0000001 1e-07 0 2 2 -4 
+5 0.9999999000000001 -1e-07 -9.999999994736442e-08 1.0000001 1e-07 1.0000001 0 2 6 -5 
+6 0.9999999000000001 -9.999999994736442e-08 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 2 5 -7 
+7 0.9999999000000001 0.9999999000000001 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 2 8 -7 
+8 0.9999999000000001 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 2 6 -8 
+9 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 2 -6 
+10 -9.999999994736442e-08 -1e-07 0.9999999000000001 1.0000001 1e-07 1.0000001 0 2 1 -5 
+11 -9.999999994736442e-08 0.9999999000000001 -1e-07 1.0000001 1.0000001 1e-07 0 2 4 -8 
+12 -9.999999994736442e-08 0.9999999000000001 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 2 3 -7 
+1 -1e-07 -9.999999994736442e-08 -9.999999994736442e-08 1e-07 1.0000001 1.0000001 0 4 1 2 -3 -4 
+2 0.9999999000000001 -9.999999994736442e-08 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 4 5 6 -7 -8 
+3 -9.999999994736442e-08 -1e-07 -9.999999994736442e-08 1.0000001 1e-07 1.0000001 0 4 9 5 -10 -1 
+4 -9.999999994736442e-08 0.9999999000000001 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 4 11 7 -12 -3 
+5 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 4 4 11 -8 -9 
+6 -9.999999994736442e-08 -9.999999994736442e-08 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 4 2 12 -6 -10 
+1 -9.999999994736442e-08 -9.999999994736442e-08 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 1 1 6 1 2 3 4 5 6 
+$EndEntities
+$Nodes
+15 14 1 14
+0 1 0 1
+1
+0 0 1
+0 2 0 1
+2
+0 0 0
+0 3 0 1
+3
+0 1 1
+0 4 0 1
+4
+0 1 0
+0 5 0 1
+5
+1 0 1
+0 6 0 1
+6
+1 0 0
+0 7 0 1
+7
+1 1 1
+0 8 0 1
+8
+1 1 0
+2 1 0 1
+9
+0 0.5 0.5
+2 2 0 1
+10
+1 0.5 0.5
+2 3 0 1
+11
+0.5 0 0.5
+2 4 0 1
+12
+0.5 1 0.5
+2 5 0 1
+13
+0.5 0.5 0
+2 6 0 1
+14
+0.5 0.5 1
+3 1 0 0
+$EndNodes
+$Elements
+1 24 1 24
+3 1 4 24
+1 14 11 12 10 
+2 10 11 12 13 
+3 12 9 11 14 
+4 9 12 11 13 
+5 5 14 7 10 
+6 4 9 12 3 
+7 2 9 1 11 
+8 1 9 3 14 
+9 5 11 1 14 
+10 10 5 11 6 
+11 12 3 14 7 
+12 10 7 12 8 
+13 2 4 13 9
+14 13 4 8 12 
+15 8 13 10 6 
+16 13 2 11 6 
+17 14 5 10 11 
+18 14 12 3 9 
+19 14 7 12 10 
+20 9 11 1 14 
+21 11 2 13 9
+22 13 9 4 12
+23 10 8 12 13 
+24 13 11 10 6 
+$EndElements