From: Nils Margenberg <margenbn@hsu-hh.de>
Date: Thu, 9 Feb 2023 08:09:30 +0000 (+0100)
Subject: Introduce local vertex numberings in grid_in and grid_out
X-Git-Tag: v9.5.0-rc1~488^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9afc08af31b9ee0f42506c410e7cbf6abb5af9b6;p=dealii.git

Introduce local vertex numberings in grid_in and grid_out
---

diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc
index 4a8f7ba8f5..6bd7b24f9e 100644
--- a/source/grid/grid_in.cc
+++ b/source/grid/grid_in.cc
@@ -1327,6 +1327,8 @@ GridIn<dim, spacedim>::read_dbmesh(std::istream &in)
   getline(in, line);
   AssertThrow(line == "Quadrilaterals", ExcInvalidDBMESHInput(line));
 
+  static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
+    {0, 1, 5, 4, 2, 3, 7, 6}};
   std::vector<CellData<dim>> cells;
   SubCellData                subcelldata;
   unsigned int               n_cells;
@@ -1338,7 +1340,9 @@ GridIn<dim, spacedim>::read_dbmesh(std::istream &in)
       cells.emplace_back();
       for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
         {
-          in >> cells.back().vertices[GeometryInfo<dim>::ucd_to_deal[i]];
+          in >>
+            cells.back().vertices[dim == 3 ? local_vertex_numbering[i] :
+                                             GeometryInfo<dim>::ucd_to_deal[i]];
 
           AssertThrow((cells.back().vertices[i] >= 1) &&
                         (static_cast<unsigned int>(cells.back().vertices[i]) <=
@@ -2306,6 +2310,8 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
   std::map<unsigned int, types::boundary_id> boundary_ids_1d;
 
   {
+    static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
+      {0, 1, 5, 4, 2, 3, 7, 6}};
     unsigned int global_cell = 0;
     for (int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
       {
@@ -2477,7 +2483,9 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
                         GeometryInfo<dim>::vertices_per_cell)
                       {
                         in >> cells.back()
-                                .vertices[GeometryInfo<dim>::ucd_to_deal[i]];
+                                .vertices[dim == 3 ?
+                                            local_vertex_numbering[i] :
+                                            GeometryInfo<dim>::ucd_to_deal[i]];
                       }
                     else
                       {
@@ -3389,6 +3397,8 @@ GridIn<dim, spacedim>::read_assimp(const std::string &filename,
   unsigned int v_offset = 0;
   unsigned int c_offset = 0;
 
+  static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
+    {0, 1, 5, 4, 2, 3, 7, 6}};
   // The index of the mesh will be used as a material index.
   for (unsigned int m = start_mesh; m < end_mesh; ++m)
     {
@@ -3433,7 +3443,8 @@ GridIn<dim, spacedim>::read_assimp(const std::string &filename,
               for (const unsigned int f : GeometryInfo<dim>::vertex_indices())
                 {
                   cells[valid_cell]
-                    .vertices[GeometryInfo<dim>::ucd_to_deal[f]] =
+                    .vertices[dim == 3 ? local_vertex_numbering[f] :
+                                         GeometryInfo<dim>::ucd_to_deal[f]] =
                     mFaces[i].mIndices[f] + v_offset;
                 }
               cells[valid_cell].material_id = m;
diff --git a/source/grid/grid_out.cc b/source/grid/grid_out.cc
index 2f72079b9a..7dcfcb5658 100644
--- a/source/grid/grid_out.cc
+++ b/source/grid/grid_out.cc
@@ -1069,6 +1069,9 @@ GridOut::write_msh(const Triangulation<dim, spacedim> &tria,
             (msh_flags.write_lines ? n_boundary_lines(tria) : 0))
       << '\n';
 
+  static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
+    {0, 1, 5, 4, 2, 3, 7, 6}};
+
   // write cells. Enumerate cells
   // consecutively, starting with 1
   for (const auto &cell : tria.active_cell_iterators())
@@ -1083,7 +1086,9 @@ GridOut::write_msh(const Triangulation<dim, spacedim> &tria,
       for (const unsigned int vertex : cell->vertex_indices())
         {
           if (cell->reference_cell() == ReferenceCells::get_hypercube<dim>())
-            out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex]) +
+            out << cell->vertex_index(
+                     dim == 3 ? local_vertex_numbering[vertex] :
+                                GeometryInfo<dim>::ucd_to_deal[vertex]) +
                      1
                 << ' ';
           else if (cell->reference_cell() == ReferenceCells::get_simplex<dim>())
@@ -4237,6 +4242,8 @@ namespace internal
             ReferenceCells::get_hypercube<dim>(), quadrature);
         }
 
+      static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
+        {0, 1, 5, 4, 2, 3, 7, 6}};
       for (const auto &cell : tria.active_cell_iterators())
         {
           if (gnuplot_flags.write_cell_numbers)
@@ -4253,8 +4260,11 @@ namespace internal
               // points (+ the initial point again) in a row and lifting the
               // drawing pencil at the end
               for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
-                out << cell->vertex(GeometryInfo<dim>::ucd_to_deal[i]) << ' '
-                    << cell->level() << ' ' << cell->material_id() << '\n';
+                out << cell->vertex(dim == 3 ?
+                                      local_vertex_numbering[i] :
+                                      GeometryInfo<dim>::ucd_to_deal[i])
+                    << ' ' << cell->level() << ' ' << cell->material_id()
+                    << '\n';
               out << cell->vertex(0) << ' ' << cell->level() << ' '
                   << cell->material_id() << '\n'
                   << '\n' // double new line for gnuplot 3d plots