From: David Wells <drwells@email.unc.edu>
Date: Wed, 25 May 2022 14:37:57 +0000 (-0400)
Subject: Add pyramids and wedges to invert_cells_with_negative_measure().
X-Git-Tag: v9.5.0-rc1~1001^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e63e78a41b10308c17c07a593ec085a73914e176;p=dealii.git

Add pyramids and wedges to invert_cells_with_negative_measure().
---

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 101671a0ae..c9da816fc3 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -872,8 +872,8 @@ namespace GridTools
     const std::vector<Point<spacedim>> &all_vertices,
     std::vector<CellData<dim>> &        cells)
   {
-    // This function is presently only implemented for hypercube and simplex
-    // volumetric (codimension 0) elements.
+    // This function is presently only implemented for volumetric (codimension
+    // 0) elements.
 
     if (dim == 1)
       return 0;
@@ -887,10 +887,10 @@ namespace GridTools
         const ArrayView<const unsigned int> vertices(cell.vertices);
         if (GridTools::cell_measure(all_vertices, vertices) < 0)
           {
-            const unsigned int n_vertices = vertices.size();
+            const auto reference_cell =
+              ReferenceCell::n_vertices_to_type(dim, vertices.size());
 
-            if (ReferenceCell::n_vertices_to_type(dim, n_vertices)
-                  .is_hyper_cube())
+            if (reference_cell.is_hyper_cube())
               {
                 ++n_negative_cells;
 
@@ -908,16 +908,27 @@ namespace GridTools
                     std::swap(cell.vertices[5], cell.vertices[7]);
                   }
               }
-
-            else if (ReferenceCell::n_vertices_to_type(dim, n_vertices)
-                       .is_simplex())
+            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.
-                std::swap(cell.vertices[n_vertices - 2],
-                          cell.vertices[n_vertices - 1]);
+                std::swap(cell.vertices[cell.vertices.size() - 2],
+                          cell.vertices[cell.vertices.size() - 1]);
+              }
+            else if (reference_cell == ReferenceCells::Wedge)
+              {
+                // swap the two triangular faces
+                std::swap(cell.vertices[0], cell.vertices[3]);
+                std::swap(cell.vertices[1], cell.vertices[4]);
+                std::swap(cell.vertices[2], cell.vertices[5]);
+              }
+            else if (reference_cell == ReferenceCells::Pyramid)
+              {
+                // Try swapping two vertices in the base - perhaps things were
+                // read in the UCD (counter-clockwise) order instead of lexical
+                std::swap(cell.vertices[2], cell.vertices[3]);
               }
             else
               {