From e4020b792701f326e69ac58fb8245fe305b51239 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sun, 10 Jan 2021 22:44:15 +0100
Subject: [PATCH] Enabling mixed meshes in GridIn::read_vkt() amd read_msh()

---
 source/grid/grid_in.cc | 30 ++++++++++++------------------
 source/grid/tria.cc    | 26 ++++++++++++--------------
 2 files changed, 24 insertions(+), 32 deletions(-)

diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc
index 50d432be0f..2eec0e3195 100644
--- a/source/grid/grid_in.cc
+++ b/source/grid/grid_in.cc
@@ -529,17 +529,12 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
       Assert(subcelldata.check_consistency(dim), ExcInternalError());
 
 
-      // make sure that only either simplex or hypercube cells are available
-      //
       // TODO: the functions below (GridTools::delete_unused_vertices(),
       // GridTools::invert_all_cells_of_negative_grid(),
-      // GridReordering::reorder_cells(),
-      // Triangulation::create_triangulation_compatibility()) need to be
-      // revisited for simplex meshes
-      AssertThrow(dim == 1 || (is_tria_or_tet_mesh ^ is_quad_or_hex_mesh),
-                  ExcNotImplemented());
+      // GridReordering::reorder_cells()) need to be
+      // revisited for simplex/mixed meshes
 
-      if (dim == 1 || is_quad_or_hex_mesh)
+      if (dim == 1 || (is_quad_or_hex_mesh && !is_tria_or_tet_mesh))
         {
           GridTools::delete_unused_vertices(vertices, cells, subcelldata);
 
@@ -556,7 +551,10 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
         }
       else
         {
-          tria->create_triangulation(vertices, cells, subcelldata);
+          // simplex or mixed mesh
+          tria->create_triangulation_compatibility(vertices,
+                                                   cells,
+                                                   subcelldata);
         }
     }
   else
@@ -2150,17 +2148,12 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
   // check that we actually read some cells.
   AssertThrow(cells.size() > 0, ExcGmshNoCellInformation());
 
-  // make sure that only either simplex or hypercube cells are available
-  //
   // TODO: the functions below (GridTools::delete_unused_vertices(),
   // GridTools::invert_all_cells_of_negative_grid(),
-  // GridReordering::reorder_cells(),
-  // Triangulation::create_triangulation_compatibility()) need to be revisited
-  // for simplex meshes
-  AssertThrow(dim == 1 || (is_tria_or_tet_mesh ^ is_quad_or_hex_mesh),
-              ExcNotImplemented());
+  // GridReordering::reorder_cells()) need to be revisited
+  // for simplex/mixed meshes
 
-  if (dim == 1 || is_quad_or_hex_mesh)
+  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);
@@ -2173,7 +2166,8 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
     }
   else
     {
-      tria->create_triangulation(vertices, cells, subcelldata);
+      // simplex or mixed mesh
+      tria->create_triangulation_compatibility(vertices, cells, subcelldata);
     }
 
   // in 1d, we also have to attach boundary ids to vertices, which does not
diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index 964bb3b598..dc48caea96 100644
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -481,7 +481,8 @@ namespace
                              const SubCellData &)
   {
     for (auto &cell : cells)
-      std::swap(cell.vertices[2], cell.vertices[3]);
+      if (cell.vertices.size() == GeometryInfo<2>::vertices_per_cell)
+        std::swap(cell.vertices[2], cell.vertices[3]);
   }
 
 
@@ -490,21 +491,18 @@ namespace
   {
     unsigned int tmp[GeometryInfo<3>::vertices_per_cell];
     for (auto &cell : cells)
-      {
-        for (const unsigned int i : GeometryInfo<3>::vertex_indices())
-          tmp[i] = cell.vertices[i];
-        for (const unsigned int i : GeometryInfo<3>::vertex_indices())
-          cell.vertices[GeometryInfo<3>::ucd_to_deal[i]] = tmp[i];
-      }
+      if (cell.vertices.size() == GeometryInfo<3>::vertices_per_cell)
+        {
+          for (const unsigned int i : GeometryInfo<3>::vertex_indices())
+            tmp[i] = cell.vertices[i];
+          for (const unsigned int i : GeometryInfo<3>::vertex_indices())
+            cell.vertices[GeometryInfo<3>::ucd_to_deal[i]] = tmp[i];
+        }
 
     // now points in boundary quads
-    std::vector<CellData<2>>::iterator boundary_quad =
-      subcelldata.boundary_quads.begin();
-    std::vector<CellData<2>>::iterator end_quad =
-      subcelldata.boundary_quads.end();
-    for (unsigned int quad_no = 0; boundary_quad != end_quad;
-         ++boundary_quad, ++quad_no)
-      std::swap(boundary_quad->vertices[2], boundary_quad->vertices[3]);
+    for (auto &boundary_quad : subcelldata.boundary_quads)
+      if (boundary_quad.vertices.size() == GeometryInfo<2>::vertices_per_cell)
+        std::swap(boundary_quad.vertices[2], boundary_quad.vertices[3]);
   }
 
 
-- 
2.39.5