]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Refactor the way we apply grid fixup functions. 13820/head
authorDavid Wells <drwells@email.unc.edu>
Wed, 25 May 2022 14:43:17 +0000 (10:43 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sat, 27 Aug 2022 17:15:23 +0000 (13:15 -0400)
This exposes a few issues with the ExodusII test grids.

doc/news/changes/incompatibilities/20220325DavidWells [new file with mode: 0644]
source/grid/grid_in.cc
source/grid/grid_tools.cc
tests/grid/grid_in_exodusii.with_trilinos_with_seacas=on.output
tests/simplex/negative_measure_01.cc [new file with mode: 0644]
tests/simplex/negative_measure_01.output [new file with mode: 0644]

diff --git a/doc/news/changes/incompatibilities/20220325DavidWells b/doc/news/changes/incompatibilities/20220325DavidWells
new file mode 100644 (file)
index 0000000..acb0802
--- /dev/null
@@ -0,0 +1,4 @@
+Changed: All GridIn functions now remove unused vertices and will attempt to fix
+pyramids and wedges with negative volumes.
+<br>
+(David Wells, 2022/05/25)
index b9a8edae923fbc4fc125487b9e42f774b418666e..882d1b0bdae38abe9e9204e9fde00b7e4027bbb1 100644 (file)
@@ -97,6 +97,35 @@ namespace
     // vertices except in 1d
     Assert(dim != 1, ExcInternalError());
   }
+
+  /**
+   * Apply each of the grid fixup routines in the correct sequence.
+   */
+  template <int dim, int spacedim>
+  void
+  apply_grid_fixup_functions(std::vector<Point<spacedim>> &vertices,
+                             std::vector<CellData<dim>> &  cells,
+                             SubCellData &                 subcelldata)
+  {
+    // check that no forbidden arrays are used
+    Assert(subcelldata.check_consistency(dim), ExcInternalError());
+    const auto n_hypercube_vertices =
+      ReferenceCells::get_hypercube<dim>().n_vertices();
+    bool is_only_hypercube = true;
+    for (const CellData<dim> &cell : cells)
+      if (cell.vertices.size() != n_hypercube_vertices)
+        {
+          is_only_hypercube = false;
+          break;
+        }
+
+    GridTools::delete_unused_vertices(vertices, cells, subcelldata);
+    if (dim == spacedim)
+      GridTools::invert_cells_with_negative_measure(vertices, cells);
+
+    if (is_only_hypercube)
+      GridTools::consistently_order_cells(cells);
+  }
 } // namespace
 
 template <int dim, int spacedim>
@@ -195,9 +224,6 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
   unsigned int n_geometric_objects = 0;
   unsigned int n_ints;
 
-  bool is_quad_or_hex_mesh = false;
-  bool is_tria_or_tet_mesh = false;
-
   if (keyword == "CELLS")
     {
       // jump to the `CELL_TYPES` section and read in cell types
@@ -235,11 +261,6 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
               // VTK_TETRA is 10, VTK_HEXAHEDRON is 12
               if (cell_types[count] == 10 || cell_types[count] == 12)
                 {
-                  if (cell_types[count] == 10)
-                    is_tria_or_tet_mesh = true;
-                  if (cell_types[count] == 12)
-                    is_quad_or_hex_mesh = true;
-
                   // we assume that the file contains first all cells,
                   // and only then any faces or lines
                   AssertThrow(subcelldata.boundary_quads.size() == 0 &&
@@ -267,11 +288,6 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
               // VTK_TRIANGLE is 5, VTK_QUAD is 9
               else if (cell_types[count] == 5 || cell_types[count] == 9)
                 {
-                  if (cell_types[count] == 5)
-                    is_tria_or_tet_mesh = true;
-                  if (cell_types[count] == 9)
-                    is_quad_or_hex_mesh = true;
-
                   // we assume that the file contains first all cells,
                   // then all faces, and finally all lines
                   AssertThrow(subcelldata.boundary_lines.size() == 0,
@@ -319,11 +335,6 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
                   AssertThrow(subcelldata.boundary_lines.size() == 0,
                               ExcNotImplemented());
 
-                  if (cell_types[count] == 5)
-                    is_tria_or_tet_mesh = true;
-                  if (cell_types[count] == 9)
-                    is_quad_or_hex_mesh = true;
-
                   cells.emplace_back(n_vertices);
 
                   for (unsigned int j = 0; j < n_vertices;
@@ -556,29 +567,8 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
               }
           }
 
-      Assert(subcelldata.check_consistency(dim), ExcInternalError());
-
-
-      // TODO: the functions below (GridTools::delete_unused_vertices(),
-      // GridTools::invert_all_negative_measure_cells(),
-      // GridTools::consistently_order_cells()) need to be
-      // revisited for simplex/mixed meshes
-
-      if (dim == 1 || (is_quad_or_hex_mesh && !is_tria_or_tet_mesh))
-        {
-          GridTools::delete_unused_vertices(vertices, cells, subcelldata);
-
-          if (dim == spacedim)
-            GridTools::invert_all_negative_measure_cells(vertices, cells);
-
-          GridTools::consistently_order_cells(cells);
-          tria->create_triangulation(vertices, cells, subcelldata);
-        }
-      else
-        {
-          // simplex or mixed mesh
-          tria->create_triangulation(vertices, cells, subcelldata);
-        }
+      apply_grid_fixup_functions(vertices, cells, subcelldata);
+      tria->create_triangulation(vertices, cells, subcelldata);
     }
   else
     AssertThrow(false,
@@ -879,15 +869,7 @@ GridIn<dim, spacedim>::read_unv(std::istream &in)
         }
     }
 
-  Assert(subcelldata.check_consistency(dim), ExcInternalError());
-
-  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
-
-  if (dim == spacedim)
-    GridTools::invert_all_negative_measure_cells(vertices, cells);
-
-  GridTools::consistently_order_cells(cells);
-
+  apply_grid_fixup_functions(vertices, cells, subcelldata);
   tria->create_triangulation(vertices, cells, subcelldata);
 }
 
@@ -1108,18 +1090,9 @@ GridIn<dim, spacedim>::read_ucd(std::istream &in,
         AssertThrow(false, ExcUnknownIdentifier(cell_type));
     }
 
-
-  // check that no forbidden arrays are used
-  Assert(subcelldata.check_consistency(dim), ExcInternalError());
-
   AssertThrow(in.fail() == false, ExcIO());
 
-  // do some clean-up on vertices...
-  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
-  // ... and cells
-  if (dim == spacedim)
-    GridTools::invert_all_negative_measure_cells(vertices, cells);
-  GridTools::consistently_order_cells(cells);
+  apply_grid_fixup_functions(vertices, cells, subcelldata);
   tria->create_triangulation(vertices, cells, subcelldata);
 }
 
@@ -1361,18 +1334,9 @@ GridIn<dim, spacedim>::read_dbmesh(std::istream &in)
     ;
   // ok, so we are not at the end of
   // the file, that's it, mostly
-
-
-  // check that no forbidden arrays are used
-  Assert(subcelldata.check_consistency(dim), ExcInternalError());
-
   AssertThrow(in.fail() == false, ExcIO());
 
-  // do some clean-up on vertices...
-  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
-  // ...and cells
-  GridTools::invert_all_negative_measure_cells(vertices, cells);
-  GridTools::consistently_order_cells(cells);
+  apply_grid_fixup_functions(vertices, cells, subcelldata);
   tria->create_triangulation(vertices, cells, subcelldata);
 }
 
@@ -1436,11 +1400,7 @@ GridIn<dim, spacedim>::read_xda(std::istream &in)
     }
   AssertThrow(in.fail() == false, ExcIO());
 
-  // do some clean-up on vertices...
-  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
-  // ... and cells
-  GridTools::invert_all_negative_measure_cells(vertices, cells);
-  GridTools::consistently_order_cells(cells);
+  apply_grid_fixup_functions(vertices, cells, subcelldata);
   tria->create_triangulation(vertices, cells, subcelldata);
 }
 
@@ -2313,8 +2273,6 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
   std::vector<CellData<dim>>                 cells;
   SubCellData                                subcelldata;
   std::map<unsigned int, types::boundary_id> boundary_ids_1d;
-  bool                                       is_quad_or_hex_mesh = false;
-  bool                                       is_tria_or_tet_mesh = false;
 
   {
     unsigned int global_cell = 0;
@@ -2465,25 +2423,13 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
                 if (cell_type == 1) // line
                   vertices_per_cell = 2;
                 else if (cell_type == 2) // tri
-                  {
-                    vertices_per_cell   = 3;
-                    is_tria_or_tet_mesh = true;
-                  }
+                  vertices_per_cell = 3;
                 else if (cell_type == 3) // quad
-                  {
-                    vertices_per_cell   = 4;
-                    is_quad_or_hex_mesh = true;
-                  }
+                  vertices_per_cell = 4;
                 else if (cell_type == 4) // tet
-                  {
-                    vertices_per_cell   = 4;
-                    is_tria_or_tet_mesh = true;
-                  }
+                  vertices_per_cell = 4;
                 else if (cell_type == 5) // hex
-                  {
-                    vertices_per_cell   = 8;
-                    is_quad_or_hex_mesh = true;
-                  }
+                  vertices_per_cell = 8;
 
                 AssertThrow(nod_num == vertices_per_cell,
                             ExcMessage(
@@ -2582,15 +2528,9 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
                 unsigned int vertices_per_cell = 0;
                 // check cell type
                 if (cell_type == 2) // tri
-                  {
-                    vertices_per_cell   = 3;
-                    is_tria_or_tet_mesh = true;
-                  }
+                  vertices_per_cell = 3;
                 else if (cell_type == 3) // quad
-                  {
-                    vertices_per_cell   = 4;
-                    is_quad_or_hex_mesh = true;
-                  }
+                  vertices_per_cell = 4;
 
                 subcelldata.boundary_quads.emplace_back();
 
@@ -2666,10 +2606,6 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
   static const std::string end_elements_marker[] = {"$ENDELM", "$EndElements"};
   AssertThrow(line == end_elements_marker[gmsh_file_format == 10 ? 0 : 1],
               ExcInvalidGMSHInput(line));
-
-  // check that no forbidden arrays are used
-  Assert(subcelldata.check_consistency(dim), ExcInternalError());
-
   AssertThrow(in.fail() == false, ExcIO());
 
   // check that we actually read some cells.
@@ -2677,25 +2613,7 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
               ExcGmshNoCellInformation(subcelldata.boundary_lines.size(),
                                        subcelldata.boundary_quads.size()));
 
-  // TODO: the functions below (GridTools::delete_unused_vertices(),
-  // GridTools::invert_all_negative_measure_cells(),
-  // GridTools::consistently_order_cells()) need to be revisited
-  // for simplex/mixed meshes
-
-  if (dim == 1 || (is_quad_or_hex_mesh && !is_tria_or_tet_mesh))
-    {
-      // do some clean-up on vertices...
-      GridTools::delete_unused_vertices(vertices, cells, subcelldata);
-      // ... and cells
-      if (dim == spacedim)
-        GridTools::invert_cells_with_negative_measure(vertices, cells);
-      GridTools::consistently_order_cells(cells);
-    }
-  else if (is_tria_or_tet_mesh)
-    {
-      if (dim == spacedim)
-        GridTools::invert_cells_with_negative_measure(vertices, cells);
-    }
+  apply_grid_fixup_functions(vertices, cells, subcelldata);
   tria->create_triangulation(vertices, cells, subcelldata);
 
   // in 1d, we also have to attach boundary ids to vertices, which does not
@@ -2903,8 +2821,7 @@ GridIn<dim, spacedim>::read_msh(const std::string &fname)
         }
     }
 
-  Assert(subcelldata.check_consistency(dim), ExcInternalError());
-
+  apply_grid_fixup_functions(vertices, cells, subcelldata);
   tria->create_triangulation(vertices, cells, subcelldata);
 
   // in 1d, we also have to attach boundary ids to vertices, which does not
@@ -3375,19 +3292,10 @@ GridIn<2>::read_tecplot(std::istream &in)
           for (const unsigned int j : GeometryInfo<dim>::vertex_indices())
             in >> cells[i].vertices[GeometryInfo<dim>::ucd_to_deal[j]];
         }
-      // do some clean-up on vertices
-      GridTools::delete_unused_vertices(vertices, cells, subcelldata);
     }
-
-  // check that no forbidden arrays are
-  // used. as we do not read in any
-  // subcelldata, nothing should happen here.
-  Assert(subcelldata.check_consistency(dim), ExcInternalError());
   AssertThrow(in.fail() == false, ExcIO());
 
-  // do some cleanup on cells
-  GridTools::invert_all_negative_measure_cells(vertices, cells);
-  GridTools::consistently_order_cells(cells);
+  apply_grid_fixup_functions(vertices, cells, subcelldata);
   tria->create_triangulation(vertices, cells, subcelldata);
 }
 
@@ -3539,10 +3447,7 @@ GridIn<dim, spacedim>::read_assimp(const std::string &filename,
         }
     }
 
-  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
-  if (dim == spacedim)
-    GridTools::invert_all_negative_measure_cells(vertices, cells);
-  GridTools::consistently_order_cells(cells);
+  apply_grid_fixup_functions(vertices, cells, subcelldata);
   tria->create_triangulation(vertices, cells, subcelldata);
 
 #else
@@ -3842,8 +3747,6 @@ GridIn<dim, spacedim>::read_exodusii(
   ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data());
   AssertThrowExodusII(ierr);
 
-  bool is_only_quad_or_hex = true;
-
   std::vector<CellData<dim>> cells;
   // Elements are grouped together by same reference cell type in element
   // blocks. There may be multiple blocks for a single reference cell type,
@@ -3871,9 +3774,6 @@ GridIn<dim, spacedim>::read_exodusii(
       const ReferenceCell type =
         exodusii_name_to_type(string_temp.data(), n_nodes_per_element);
 
-      if (type.is_simplex())
-        is_only_quad_or_hex = false;
-
       // The number of nodes per element may be larger than what we want to
       // read - for example, if the Exodus file contains a QUAD9 element, we
       // only want to read the first four values and ignore the rest.
@@ -3908,16 +3808,7 @@ GridIn<dim, spacedim>::read_exodusii(
   ierr = ex_close(ex_id);
   AssertThrowExodusII(ierr);
 
-  if (is_only_quad_or_hex)
-    {
-      // do some clean-up on vertices...
-      GridTools::delete_unused_vertices(vertices, cells, pair.first);
-      // ... and cells
-      if (dim == spacedim)
-        GridTools::invert_all_negative_measure_cells(vertices, cells);
-      GridTools::consistently_order_cells(cells);
-    }
-
+  apply_grid_fixup_functions(vertices, cells, pair.first);
   tria->create_triangulation(vertices, cells, pair.first);
   ExodusIIData out;
   out.id_to_sideset_ids = std::move(pair.second);
index c9da816fc3a0e0e13b2c7ad47e079caf21b954cb..f6f642ab1fbed2b4de64616c72596c1cf89b6a08 100644 (file)
@@ -885,15 +885,16 @@ namespace GridTools
     for (auto &cell : cells)
       {
         const ArrayView<const unsigned int> vertices(cell.vertices);
-        if (GridTools::cell_measure(all_vertices, vertices) < 0)
+        // Some pathologically twisted cells can have exactly zero measure but
+        // we can still fix them
+        if (GridTools::cell_measure(all_vertices, vertices) <= 0)
           {
+            ++n_negative_cells;
             const auto reference_cell =
               ReferenceCell::n_vertices_to_type(dim, vertices.size());
 
             if (reference_cell.is_hyper_cube())
               {
-                ++n_negative_cells;
-
                 if (dim == 2)
                   {
                     // flip the cell across the y = x line in 2D
@@ -910,7 +911,6 @@ namespace GridTools
               }
             else if (reference_cell.is_simplex())
               {
-                ++n_negative_cells;
                 // By basic rules for computing determinants we can just swap
                 // two vertices to fix a negative volume. Arbitrarily pick the
                 // last two.
index 14d6585b069ba5b8c4c2752bd0fb81498c53a4d4..b1b43d02d4b5101089ae9ee34e38b5acb171b07b 100644 (file)
@@ -123,18 +123,18 @@ DEAL::face center = 1.50000 3.00000 1.50000   face boundary id = 0
 DEAL::face center = 1.33333 0.333333 3.66667   face boundary id = 0
 DEAL::face center = 0.333333 1.33333 3.66667   face boundary id = 0
 DEAL::face center = 2.33333 0.666667 4.33333   face boundary id = 0
-DEAL::face center = 5.50000 0.500000 4.00000   face boundary id = 0
-DEAL::face center = 5.50000 1.50000 3.00000   face boundary id = 0
 DEAL::face center = 1.33333 2.33333 3.66667   face boundary id = 0
 DEAL::face center = 2.33333 1.66667 4.33333   face boundary id = 0
-DEAL::face center = 5.50000 2.00000 4.00000   face boundary id = 0
+DEAL::face center = 5.50000 0.500000 4.00000   face boundary id = 0
+DEAL::face center = 5.50000 1.50000 3.00000   face boundary id = 0
 DEAL::face center = 8.00000 1.33333 3.66667   face boundary id = 0
+DEAL::face center = 5.50000 2.00000 4.00000   face boundary id = 0
 DEAL::Number of vertices: 13
 DEAL::Number of cells: 4
 DEAL::cell 0 type = Hex volume = 27.0000
 DEAL::cell 1 type = Pyramid volume = 6.00000
 DEAL::cell 2 type = Tet volume = 2.00000
-DEAL::cell 3 type = Wedge volume = -15.0000
+DEAL::cell 3 type = Wedge volume = 15.0000
 # vtk DataFile Version 3.0
 Triangulation generated with deal.II
 ASCII
@@ -158,7 +158,7 @@ CELLS 4 27
 8 0 1 2 3 4 5 6 7
 5 4 5 6 7 8
 4 5 9 6 8
-6 6 5 9 11 10 12
+6 11 10 12 6 5 9
 
 CELL_TYPES 4
 12 14 10 13 
@@ -189,23 +189,23 @@ DEAL::face center = 1.50000 3.00000 1.50000   face boundary id = 0
 DEAL::face center = 1.33333 0.333333 3.66667   face boundary id = 0
 DEAL::face center = 0.333333 1.33333 3.66667   face boundary id = 0
 DEAL::face center = 2.33333 0.666667 4.33333   face boundary id = 0
-DEAL::face center = 5.50000 0.500000 4.00000   face boundary id = 0
-DEAL::face center = 5.50000 1.50000 3.00000   face boundary id = 0
 DEAL::face center = 1.33333 2.33333 3.66667   face boundary id = 0
 DEAL::face center = 2.33333 1.66667 4.33333   face boundary id = 0
-DEAL::face center = 5.50000 2.00000 4.00000   face boundary id = 0
+DEAL::face center = 5.50000 0.500000 4.00000   face boundary id = 0
+DEAL::face center = 5.50000 1.50000 3.00000   face boundary id = 0
 DEAL::face center = 8.00000 1.33333 3.66667   face boundary id = 0
-DEAL::Number of vertices: 29
+DEAL::face center = 5.50000 2.00000 4.00000   face boundary id = 0
+DEAL::Number of vertices: 13
 DEAL::Number of cells: 4
 DEAL::cell 0 type = Hex volume = 27.0000
 DEAL::cell 1 type = Pyramid volume = 6.00000
 DEAL::cell 2 type = Tet volume = 2.00000
-DEAL::cell 3 type = Wedge volume = -15.0000
+DEAL::cell 3 type = Wedge volume = 15.0000
 # vtk DataFile Version 3.0
 Triangulation generated with deal.II
 ASCII
 DATASET UNSTRUCTURED_GRID
-POINTS 29 double
+POINTS 13 double
 0.0 0.0 0.0
 3.0 0.0 0.0
 3.0 3.0 0.0
@@ -219,28 +219,12 @@ POINTS 29 double
 8.0 0.0 3.0
 8.0 3.0 3.0
 8.0 1.0 5.0
-3.0 1.5 0.0
-1.5 3.0 0.0
-0.0 1.5 0.0
-1.5 0.0 0.0
-3.0 0.0 1.5
-3.0 3.0 1.5
-0.0 3.0 1.5
-0.0 0.0 1.5
-3.0 1.5 3.0
-1.5 3.0 3.0
-0.0 1.5 3.0
-1.5 0.0 3.0
-0.50 0.50 4.0
-2.0 0.50 4.0
-2.0 2.0 4.0
-0.50 2.0 4.0
 
 CELLS 4 27
 8 0 1 2 3 4 5 6 7
 5 4 5 6 7 8
 4 5 9 6 8
-6 6 5 9 11 10 12
+6 11 10 12 6 5 9
 
 CELL_TYPES 4
 12 14 10 13 
@@ -277,18 +261,18 @@ DEAL::face center = 1.50000 3.00000 1.50000   face boundary id = 4
 DEAL::face center = 1.33333 0.333333 3.66667   face boundary id = 1
 DEAL::face center = 0.333333 1.33333 3.66667   face boundary id = 5
 DEAL::face center = 2.33333 0.666667 4.33333   face boundary id = 1
-DEAL::face center = 5.50000 0.500000 4.00000   face boundary id = 1
-DEAL::face center = 5.50000 1.50000 3.00000   face boundary id = 4
 DEAL::face center = 1.33333 2.33333 3.66667   face boundary id = 4
 DEAL::face center = 2.33333 1.66667 4.33333   face boundary id = 3
-DEAL::face center = 5.50000 2.00000 4.00000   face boundary id = 3
+DEAL::face center = 5.50000 0.500000 4.00000   face boundary id = 1
+DEAL::face center = 5.50000 1.50000 3.00000   face boundary id = 4
 DEAL::face center = 8.00000 1.33333 3.66667   face boundary id = 6
+DEAL::face center = 5.50000 2.00000 4.00000   face boundary id = 3
 DEAL::Number of vertices: 13
 DEAL::Number of cells: 4
 DEAL::cell 0 type = Hex volume = 27.0000
 DEAL::cell 1 type = Pyramid volume = 6.00000
 DEAL::cell 2 type = Tet volume = 2.00000
-DEAL::cell 3 type = Wedge volume = -15.0000
+DEAL::cell 3 type = Wedge volume = 15.0000
 # vtk DataFile Version 3.0
 Triangulation generated with deal.II
 ASCII
@@ -312,7 +296,7 @@ CELLS 18 91
 8 0 1 2 3 4 5 6 7
 5 4 5 6 7 8
 4 5 9 6 8
-6 6 5 9 11 10 12
+6 11 10 12 6 5 9
 4 0 4 5 1
 4 0 1 2 3
 4 0 3 7 4
@@ -321,23 +305,23 @@ CELLS 18 91
 3 5 4 8
 3 4 7 8
 3 9 5 8
-4 5 9 12 10
-4 6 5 10 11
 3 7 6 8
 3 6 9 8
-4 9 6 11 12
-3 11 10 12
+4 10 12 9 5
+4 11 10 5 6
+3 10 11 12
+4 12 11 6 9
 
 CELL_TYPES 18
 12 14 10 13 
-9 9 9 9 9 5 5 5 9 9 5 5 9 5 
+9 9 9 9 9 5 5 5 5 5 9 9 5 9 
 
 
 CELL_DATA 18
 SCALARS MaterialID int 1
 LOOKUP_TABLE default
 1 10 100 1000 
-1 2 5 3 4 1 5 1 1 4 4 3 3 6 
+1 2 5 3 4 1 5 1 4 3 1 4 6 3 
 
 
 SCALARS ManifoldID int 1
@@ -363,18 +347,18 @@ DEAL::face center = 1.50000 3.00000 1.50000   face manifold id = 4
 DEAL::face center = 1.33333 0.333333 3.66667   face manifold id = 1
 DEAL::face center = 0.333333 1.33333 3.66667   face manifold id = 5
 DEAL::face center = 2.33333 0.666667 4.33333   face manifold id = 1
-DEAL::face center = 5.50000 0.500000 4.00000   face manifold id = 1
-DEAL::face center = 5.50000 1.50000 3.00000   face manifold id = 4
 DEAL::face center = 1.33333 2.33333 3.66667   face manifold id = 4
 DEAL::face center = 2.33333 1.66667 4.33333   face manifold id = 3
-DEAL::face center = 5.50000 2.00000 4.00000   face manifold id = 3
+DEAL::face center = 5.50000 0.500000 4.00000   face manifold id = 1
+DEAL::face center = 5.50000 1.50000 3.00000   face manifold id = 4
 DEAL::face center = 8.00000 1.33333 3.66667   face manifold id = 6
+DEAL::face center = 5.50000 2.00000 4.00000   face manifold id = 3
 DEAL::Number of vertices: 13
 DEAL::Number of cells: 4
 DEAL::cell 0 type = Hex volume = 27.0000
 DEAL::cell 1 type = Pyramid volume = 6.00000
 DEAL::cell 2 type = Tet volume = 2.00000
-DEAL::cell 3 type = Wedge volume = -15.0000
+DEAL::cell 3 type = Wedge volume = 15.0000
 # vtk DataFile Version 3.0
 Triangulation generated with deal.II
 ASCII
@@ -398,7 +382,7 @@ CELLS 18 91
 8 0 1 2 3 4 5 6 7
 5 4 5 6 7 8
 4 5 9 6 8
-6 6 5 9 11 10 12
+6 11 10 12 6 5 9
 4 0 4 5 1
 4 0 1 2 3
 4 0 3 7 4
@@ -407,16 +391,16 @@ CELLS 18 91
 3 5 4 8
 3 4 7 8
 3 9 5 8
-4 5 9 12 10
-4 6 5 10 11
 3 7 6 8
 3 6 9 8
-4 9 6 11 12
-3 11 10 12
+4 10 12 9 5
+4 11 10 5 6
+3 10 11 12
+4 12 11 6 9
 
 CELL_TYPES 18
 12 14 10 13 
-9 9 9 9 9 5 5 5 9 9 5 5 9 5 
+9 9 9 9 9 5 5 5 5 5 9 9 5 9 
 
 
 CELL_DATA 18
@@ -429,6 +413,6 @@ LOOKUP_TABLE default
 SCALARS ManifoldID int 1
 LOOKUP_TABLE default
 -1 -1 -1 -1 
-1 2 5 3 4 1 5 1 1 4 4 3 3 6 
+1 2 5 3 4 1 5 1 4 3 1 4 6 3 
 
 DEAL::OK
diff --git a/tests/simplex/negative_measure_01.cc b/tests/simplex/negative_measure_01.cc
new file mode 100644 (file)
index 0000000..f3c4097
--- /dev/null
@@ -0,0 +1,162 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2002 - 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test GridTools::invert_cells_with_negative_measure() with a mixed mesh.
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <string>
+
+#include "../tests.h"
+
+void
+test()
+{
+  deal_II_exceptions::disable_abort_on_exception();
+
+  // test a grid with a pyramid that cannot be fixed
+  const static int dim = 3;
+  try
+    {
+      std::vector<Point<dim>> vertices(10);
+      // hypercube
+      vertices[0] = Point<dim>(0, 0, 0);
+      vertices[1] = Point<dim>(1, 0, 0);
+      vertices[2] = Point<dim>(0, 1, 0);
+      vertices[3] = Point<dim>(1, 1, 0);
+      vertices[4] = Point<dim>(0, 0, 1);
+      vertices[5] = Point<dim>(1, 0, 1);
+      vertices[6] = Point<dim>(0, 1, 1);
+      vertices[7] = Point<dim>(1, 1, 1);
+      // pyramid on top
+      vertices[8] = Point<dim>(0.75, 0.5, 1.5);
+      // pyramid on right - a mess of an element
+      vertices[9] = Point<dim>(1.0, 1.25, 0.5);
+
+
+      std::vector<CellData<dim>> cells(3);
+      cells[0].vertices = {0, 1, 2, 3, 4, 5, 6, 7};
+      // this one is twisted
+      cells[1].vertices = {4, 5, 7, 6, 8};
+      // this one is twisted in an unfixable way
+      cells[2].vertices = {1, 3, 5, 7, 9};
+
+      SubCellData       subcelldata;
+      const std::size_t n_cells_inverted =
+        GridTools::invert_cells_with_negative_measure(vertices, cells);
+
+      deallog << "We inverted " << n_cells_inverted << " cell(s)." << std::endl;
+    }
+  catch (const ExcGridHasInvalidCell &exception)
+    {
+      deallog << "Successfully failed to set up a mesh with a pyramid "
+              << "pointing into another element" << std::endl;
+    }
+
+  // test with a fixable pyramid
+  {
+    std::vector<Point<dim>> vertices(10);
+    // hypercube
+    vertices[0] = Point<dim>(0, 0, 0);
+    vertices[1] = Point<dim>(1, 0, 0);
+    vertices[2] = Point<dim>(0, 1, 0);
+    vertices[3] = Point<dim>(1, 1, 0);
+    vertices[4] = Point<dim>(0, 0, 1);
+    vertices[5] = Point<dim>(1, 0, 1);
+    vertices[6] = Point<dim>(0, 1, 1);
+    vertices[7] = Point<dim>(1, 1, 1);
+    // pyramid on top
+    vertices[8] = Point<dim>(0.75, 0.5, 1.5);
+    // pyramid on right
+    vertices[9] = Point<dim>(1.25, 0.75, 0.5);
+
+
+    std::vector<CellData<dim>> cells(3);
+    cells[0].vertices = {0, 1, 2, 3, 4, 5, 6, 7};
+    // this one is twisted
+    cells[1].vertices = {4, 5, 7, 6, 8};
+    // this one is not twisted
+    cells[2].vertices = {5, 1, 7, 3, 9};
+
+    SubCellData       subcelldata;
+    const std::size_t n_cells_inverted =
+      GridTools::invert_cells_with_negative_measure(vertices, cells);
+
+    deallog << "We inverted " << n_cells_inverted << " cell(s)." << std::endl;
+
+    Triangulation<dim> tria;
+    tria.create_triangulation(vertices, cells, subcelldata);
+
+    std::ostream &logfile = deallog.get_file_stream();
+    logfile << "---------------------------------------------" << std::endl
+            << std::endl
+            << std::endl;
+    GridOut grid_out;
+    grid_out.write_vtk(tria, logfile);
+  }
+
+  // test with a twisted and untwisted wedge
+  {
+    std::vector<Point<dim>> vertices(12);
+    // hypercube
+    vertices[0] = Point<dim>(0, 0, 0);
+    vertices[1] = Point<dim>(1, 0, 0);
+    vertices[2] = Point<dim>(0, 1, 0);
+    vertices[3] = Point<dim>(1, 1, 0);
+    vertices[4] = Point<dim>(0, 0, 1);
+    vertices[5] = Point<dim>(1, 0, 1);
+    vertices[6] = Point<dim>(0, 1, 1);
+    vertices[7] = Point<dim>(1, 1, 1);
+    // wedge on top
+    vertices[8] = Point<dim>(0.0, 0.5, 1.5);
+    vertices[9] = Point<dim>(1.0, 0.6, 1.6);
+    // wedge on right
+    vertices[10] = Point<dim>(1.0, 0.5, 0.0);
+    vertices[11] = Point<dim>(1.0, 0.6, 0.9);
+
+
+    std::vector<CellData<dim>> cells(3);
+    cells[0].vertices = {0, 1, 2, 3, 4, 5, 6, 7};
+    // this one is twisted
+    cells[1].vertices = {5, 7, 9, 4, 6, 8};
+    // this one is not twisted
+    cells[2].vertices = {1, 10, 3, 5, 11, 7};
+
+    SubCellData       subcelldata;
+    const std::size_t n_cells_inverted =
+      GridTools::invert_cells_with_negative_measure(vertices, cells);
+
+    deallog << "We inverted " << n_cells_inverted << " cell(s)." << std::endl;
+
+    Triangulation<dim> tria;
+    tria.create_triangulation(vertices, cells, subcelldata);
+
+    std::ostream &logfile = deallog.get_file_stream();
+    logfile << "---------------------------------------------" << std::endl
+            << std::endl
+            << std::endl;
+    GridOut grid_out;
+    grid_out.write_vtk(tria, logfile);
+  }
+}
+
+int
+main()
+{
+  initlog();
+  test();
+}
diff --git a/tests/simplex/negative_measure_01.output b/tests/simplex/negative_measure_01.output
new file mode 100644 (file)
index 0000000..20e2135
--- /dev/null
@@ -0,0 +1,88 @@
+
+DEAL::Successfully failed to set up a mesh with a pyramid pointing into another element
+DEAL::We inverted 1 cell(s).
+---------------------------------------------
+
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 10 double
+0.00000 0.00000 0.00000
+1.00000 0.00000 0.00000
+0.00000 1.00000 0.00000
+1.00000 1.00000 0.00000
+0.00000 0.00000 1.00000
+1.00000 0.00000 1.00000
+0.00000 1.00000 1.00000
+1.00000 1.00000 1.00000
+0.750000 0.500000 1.50000
+1.25000 0.750000 0.500000
+
+CELLS 3 21
+8 0 1 3 2 4 5 7 6
+5 4 5 7 6 8
+5 5 1 3 7 9
+
+CELL_TYPES 3
+12 14 14 
+
+
+
+CELL_DATA 3
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 
+
+
+DEAL::We inverted 1 cell(s).
+---------------------------------------------
+
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 12 double
+0.00000 0.00000 0.00000
+1.00000 0.00000 0.00000
+0.00000 1.00000 0.00000
+1.00000 1.00000 0.00000
+0.00000 0.00000 1.00000
+1.00000 0.00000 1.00000
+0.00000 1.00000 1.00000
+1.00000 1.00000 1.00000
+0.00000 0.500000 1.50000
+1.00000 0.600000 1.60000
+1.00000 0.500000 0.00000
+1.00000 0.600000 0.900000
+
+CELLS 3 23
+8 0 1 3 2 4 5 7 6
+6 4 6 8 5 7 9
+6 1 10 3 5 11 7
+
+CELL_TYPES 3
+12 13 13 
+
+
+
+CELL_DATA 3
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 
+
+

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.