]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove std::abs() from the tet volume calculation. 13765/head
authorDavid Wells <drwells@email.unc.edu>
Fri, 20 May 2022 17:45:52 +0000 (13:45 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sat, 21 May 2022 21:25:58 +0000 (17:25 -0400)
source/grid/grid_tools.cc
source/grid/grid_tools_nontemplates.cc
tests/simplex/grid_in_msh.cc
tests/simplex/grid_in_msh.output
tests/simplex/grid_in_msh/tet_half_twisted.msh [new file with mode: 0644]

index 85ba1aa885af184efa846bf510b8c6a3b757d799..d15f5b4c11831a4d451e4f3da38b2aaa6d5c7cfb 100644 (file)
@@ -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
               {
index 03e867d3cc5bfc401110d709d2c47d71ade4aca8..cf2a3be76af457ee2c833f6a8459309b0b8b59a2 100644 (file)
@@ -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() == 6) // wedge
       {
index 087581a00984eb80a81589fdf3544a1ef49f24f2..75fc6bcf097dfceb2b6d3d6737476fbe4ddf235b 100644 (file)
@@ -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();
 }
index b1fa345b509c0d34589e645ff8fb62cec33f64ab..a76b63f017e5fd1dbf66a36909528424baece836 100644 (file)
@@ -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
new file mode 100644 (file)
index 0000000..720e270
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.