From 15f66452f4ec7b27de84693ccbc8940d7d1b69a1 Mon Sep 17 00:00:00 2001 From: David Wells Date: Fri, 20 May 2022 13:45:52 -0400 Subject: [PATCH] Remove std::abs() from the tet volume calculation. --- source/grid/grid_tools.cc | 26 ++--- source/grid/grid_tools_nontemplates.cc | 2 +- tests/simplex/grid_in_msh.cc | 5 + tests/simplex/grid_in_msh.output | 64 +++++++++++ .../simplex/grid_in_msh/tet_half_twisted.msh | 107 ++++++++++++++++++ 5 files changed, 184 insertions(+), 20 deletions(-) create mode 100644 tests/simplex/grid_in_msh/tet_half_twisted.msh diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 85ba1aa885..d15f5b4c11 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -871,16 +871,8 @@ namespace GridTools const std::vector> &all_vertices, std::vector> & 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 index 03e867d3cc..cf2a3be76a 100644 --- 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() == 6) // wedge { diff --git a/tests/simplex/grid_in_msh.cc b/tests/simplex/grid_in_msh.cc index 087581a009..75fc6bcf09 100644 --- 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 index b1fa345b50..a76b63f017 100644 --- 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 new file mode 100644 index 0000000000..720e270446 --- /dev/null +++ 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 -- 2.39.5