]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix dim = 1 mesh loading with gmsh. 15307/head
authorDavid Wells <drwells@email.unc.edu>
Mon, 5 Jun 2023 18:27:23 +0000 (14:27 -0400)
committerDavid Wells <drwells@email.unc.edu>
Mon, 5 Jun 2023 18:33:18 +0000 (14:33 -0400)
source/grid/grid_in.cc
tests/grid/grid_in_gmsh_04.cc [new file with mode: 0644]
tests/grid/grid_in_gmsh_04.msh [new file with mode: 0644]
tests/grid/grid_in_gmsh_04.output [new file with mode: 0644]

index 03a3ae82083b92a546bef7f02dcd03823ecd3f6a..24e456b551a280dfc631654a62d77afa1911062b 100644 (file)
@@ -64,35 +64,41 @@ namespace
    * boundary indicators on vertices after the triangulation has already been
    * created.
    *
-   * TODO: Fix this properly via SubcellData
+   * This function ignores all non-boundary vertices. Boundary ids can only be
+   * assigned to internal vertices - however, the grid fixup routines (e.g.,
+   * combining duplicated vertices) may switch vertices from boundary to
+   * non-boundary. Get around this in two ways:
+   *
+   * 1. Ignore any attempt to set a boundary id on a non-boundary vertex.
+   * 2. Work around possible vertex renumberings by using Point<spacedim>
+   *   instead of the corresponding vertex number (which may be out of date)
+   *
+   * TODO: fix this properly via SubCellData
    */
   template <int spacedim>
   void
   assign_1d_boundary_ids(
-    const std::map<unsigned int, types::boundary_id> &boundary_ids,
-    Triangulation<1, spacedim> &                      triangulation)
+    const std::vector<std::pair<Point<spacedim>, types::boundary_id>>
+      &                         boundary_ids,
+    Triangulation<1, spacedim> &triangulation)
   {
-    if (boundary_ids.size() > 0)
-      for (const auto &cell : triangulation.active_cell_iterators())
-        for (unsigned int f : GeometryInfo<1>::face_indices())
-          if (boundary_ids.find(cell->vertex_index(f)) != boundary_ids.end())
-            {
-              AssertThrow(
-                cell->at_boundary(f),
-                ExcMessage(
-                  "You are trying to prescribe boundary ids on the face "
-                  "of a 1d cell (i.e., on a vertex), but this face is not actually at "
-                  "the boundary of the mesh. This is not allowed."));
-              cell->face(f)->set_boundary_id(
-                boundary_ids.find(cell->vertex_index(f))->second);
-            }
+    for (auto &cell : triangulation.active_cell_iterators())
+      for (unsigned int face_no : ReferenceCells::Line.face_indices())
+        if (cell->face(face_no)->at_boundary())
+          for (const auto &pair : boundary_ids)
+            if (cell->face(face_no)->vertex(0) == pair.first)
+              {
+                cell->face(face_no)->set_boundary_id(pair.second);
+                break;
+              }
   }
 
 
   template <int dim, int spacedim>
   void
-  assign_1d_boundary_ids(const std::map<unsigned int, types::boundary_id> &,
-                         Triangulation<dim, spacedim> &)
+  assign_1d_boundary_ids(
+    const std::vector<std::pair<Point<spacedim>, types::boundary_id>> &,
+    Triangulation<dim, spacedim> &)
   {
     // we shouldn't get here since boundary ids are not assigned to
     // vertices except in 1d
@@ -2319,6 +2325,12 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
   SubCellData                                subcelldata;
   std::map<unsigned int, types::boundary_id> boundary_ids_1d;
 
+  // Track the number of times each vertex is used in 1D. This determines
+  // whether or not we can assign a boundary id to a vertex. This is necessary
+  // because sometimes gmsh saves internal vertices in the $ELEM list in codim
+  // 1 or codim 2.
+  std::map<unsigned int, unsigned int> vertex_counts;
+
   {
     static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
       {0, 1, 5, 4, 2, 3, 7, 6}};
@@ -2485,22 +2497,18 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
 
                 // allocate and read indices
                 cells.emplace_back();
-                cells.back().vertices.resize(vertices_per_cell);
+                CellData<dim> &cell = cells.back();
+                cell.vertices.resize(vertices_per_cell);
                 for (unsigned int i = 0; i < vertices_per_cell; ++i)
                   {
                     // hypercube cells need to be reordered
                     if (vertices_per_cell ==
                         GeometryInfo<dim>::vertices_per_cell)
-                      {
-                        in >> cells.back()
-                                .vertices[dim == 3 ?
+                      in >> cell.vertices[dim == 3 ?
                                             local_vertex_numbering[i] :
                                             GeometryInfo<dim>::ucd_to_deal[i]];
-                      }
                     else
-                      {
-                        in >> cells.back().vertices[i];
-                      }
+                      in >> cell.vertices[i];
                   }
 
                 // to make sure that the cast won't fail
@@ -2514,21 +2522,21 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
                 // numbers::invalid_material_id-1
                 AssertIndexRange(material_id, numbers::invalid_material_id);
 
-                cells.back().material_id = material_id;
+                cell.material_id = material_id;
 
                 // transform from gmsh to consecutive numbering
                 for (unsigned int i = 0; i < vertices_per_cell; ++i)
                   {
-                    AssertThrow(
-                      vertex_indices.find(cells.back().vertices[i]) !=
-                        vertex_indices.end(),
-                      ExcInvalidVertexIndexGmsh(cell_per_entity,
-                                                elm_number,
-                                                cells.back().vertices[i]));
-
-                    // vertex with this index exists
-                    cells.back().vertices[i] =
-                      vertex_indices[cells.back().vertices[i]];
+                    AssertThrow(vertex_indices.find(cell.vertices[i]) !=
+                                  vertex_indices.end(),
+                                ExcInvalidVertexIndexGmsh(cell_per_entity,
+                                                          elm_number,
+                                                          cell.vertices[i]));
+
+                    const auto vertex = vertex_indices[cell.vertices[i]];
+                    if (dim == 1)
+                      vertex_counts[vertex] += 1u;
+                    cell.vertices[i] = vertex;
                   }
               }
             else if ((cell_type == 1) &&
@@ -2662,13 +2670,23 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
               ExcGmshNoCellInformation(subcelldata.boundary_lines.size(),
                                        subcelldata.boundary_quads.size()));
 
+  // apply_grid_fixup_functions() may invalidate the vertex indices (since it
+  // will delete duplicated or unused vertices). Get around this by storing
+  // Points directly in that case so that the comparisons are valid.
+  std::vector<std::pair<Point<spacedim>, types::boundary_id>> boundary_id_pairs;
+  if (dim == 1)
+    for (const auto &pair : vertex_counts)
+      if (pair.second == 1u)
+        boundary_id_pairs.emplace_back(vertices[pair.first],
+                                       boundary_ids_1d[pair.first]);
+
   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
-  // currently work through the call above
+  // currently work through the call above.
   if (dim == 1)
-    assign_1d_boundary_ids(boundary_ids_1d, *tria);
+    assign_1d_boundary_ids(boundary_id_pairs, *tria);
 }
 
 
@@ -2699,6 +2717,12 @@ GridIn<dim, spacedim>::read_msh(const std::string &fname)
   SubCellData                                subcelldata;
   std::map<unsigned int, types::boundary_id> boundary_ids_1d;
 
+  // Track the number of times each vertex is used in 1D. This determines
+  // whether or not we can assign a boundary id to a vertex. This is necessary
+  // because sometimes gmsh saves internal vertices in the $ELEM list in codim
+  // 1 or codim 2.
+  std::map<unsigned int, unsigned int> vertex_counts;
+
   gmsh::initialize();
   gmsh::option::setNumber("General.Verbosity", 0);
   gmsh::open(fname);
@@ -2832,8 +2856,12 @@ GridIn<dim, spacedim>::read_msh(const std::string &fname)
                   cells.emplace_back(n_vertices);
                   auto &cell = cells.back();
                   for (unsigned int v = 0; v < n_vertices; ++v)
-                    cell.vertices[v] =
-                      nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
+                    {
+                      cell.vertices[v] =
+                        nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
+                      if (dim == 1)
+                        vertex_counts[cell.vertices[v]] += 1u;
+                    }
                   cell.manifold_id = manifold_id;
                   cell.material_id = boundary_id;
                 }
@@ -2870,13 +2898,23 @@ GridIn<dim, spacedim>::read_msh(const std::string &fname)
         }
     }
 
+  // apply_grid_fixup_functions() may invalidate the vertex indices (since it
+  // will delete duplicated or unused vertices). Get around this by storing
+  // Points directly in that case so that the comparisons are valid.
+  std::vector<std::pair<Point<spacedim>, types::boundary_id>> boundary_id_pairs;
+  if (dim == 1)
+    for (const auto &pair : vertex_counts)
+      if (pair.second == 1u)
+        boundary_id_pairs.emplace_back(vertices[pair.first],
+                                       boundary_ids_1d[pair.first]);
+
   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
-  // currently work through the call above
+  // currently work through the call above.
   if (dim == 1)
-    assign_1d_boundary_ids(boundary_ids_1d, *tria);
+    assign_1d_boundary_ids(boundary_id_pairs, *tria);
 
   gmsh::clear();
   gmsh::finalize();
diff --git a/tests/grid/grid_in_gmsh_04.cc b/tests/grid/grid_in_gmsh_04.cc
new file mode 100644 (file)
index 0000000..f8eac02
--- /dev/null
@@ -0,0 +1,75 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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 that we can read dim = 1, codim = 2 meshes with boundary ids correctly.
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+#include <map>
+
+#include "../tests.h"
+
+template <int dim, int spacedim = dim>
+void
+print_mesh_info(const Triangulation<dim, spacedim> &triangulation,
+                const std::string &                 filename)
+{
+  deallog << "Mesh info:" << std::endl
+          << " dimension: " << dim << std::endl
+          << " no. of cells: " << triangulation.n_active_cells() << std::endl;
+
+  {
+    std::map<types::boundary_id, unsigned int> boundary_count;
+    for (const auto &cell : triangulation.active_cell_iterators())
+      for (unsigned int face_no : cell->face_indices())
+        if (cell->face(face_no)->at_boundary())
+          boundary_count[cell->face(face_no)->boundary_id()]++;
+
+    deallog << " boundary indicators: ";
+    for (const std::pair<const types::boundary_id, unsigned int> &pair :
+         boundary_count)
+      {
+        deallog << pair.first << "(" << pair.second << " times) ";
+      }
+    deallog << std::endl;
+  }
+
+  // Finally, produce a graphical representation of the mesh to an output
+  // file:
+  std::ofstream out(filename);
+  GridOut       grid_out;
+  grid_out.write_vtu(triangulation, out);
+  deallog << " written to " << filename << std::endl << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+  std::ifstream in(SOURCE_DIR "/grid_in_gmsh_04.msh");
+
+  Triangulation<1, 2> tria;
+  GridIn<1, 2>        gi;
+  gi.attach_triangulation(tria);
+  gi.read_msh(in);
+
+  print_mesh_info(tria, "grid_in_gmsh_04.vtu");
+}
diff --git a/tests/grid/grid_in_gmsh_04.msh b/tests/grid/grid_in_gmsh_04.msh
new file mode 100644 (file)
index 0000000..8582326
--- /dev/null
@@ -0,0 +1,113 @@
+$NOD
+49
+1 -6 -4 0
+2 -6 -4 0
+3 0 0 0
+4 6 -4 0
+5 -4 -4 0
+6 -6 0 0
+7 6 0 0
+8 4 -4 0
+9 4 -2 0
+10 -4 -2 0
+11 -4 -2.5 0
+12 -4 -3 0
+13 -4 -3.5 0
+14 -4.5 -4 0
+15 -5 -4 0
+16 -5.5 -4 0
+17 -6 -3 0
+18 -6 -2 0
+19 -6 -1 0
+20 -5.25 0 0
+21 -4.5 0 0
+22 -3.75 0 0
+23 -3 0 0
+24 -2.25 0 0
+25 -1.5 0 0
+26 -0.75 0 0
+27 0.75 0 0
+28 1.5 0 0
+29 2.25 0 0
+30 3 0 0
+31 3.75 0 0
+32 4.5 0 0
+33 5.25 0 0
+34 6 -1 0
+35 6 -2 0
+36 6 -3 0
+37 5.5 -4 0
+38 5 -4 0
+39 4.5 -4 0
+40 4 -3.5 0
+41 4 -3 0
+42 4 -2.5 0
+43 3 -2 0
+44 2 -2 0
+45 1 -2 0
+46 0 -2 0
+47 -1 -2 0
+48 -2 -2 0
+49 -3 -2 0
+$ENDNOD
+$ELM
+58
+1 15 0 1 1 1
+2 15 0 2 1 2
+3 15 0 3 1 3
+4 15 0 4 1 4
+5 15 0 6 1 5
+6 15 0 7 1 6
+7 15 0 8 1 7
+8 15 0 9 1 8
+9 15 0 10 1 9
+10 15 0 11 1 10
+11 1 0 1 2 10 11
+12 1 0 1 2 11 12
+13 1 0 1 2 12 13
+14 1 0 1 2 13 5
+15 1 0 2 2 5 14
+16 1 0 2 2 14 15
+17 1 0 2 2 15 16
+18 1 0 2 2 16 1
+19 1 0 3 2 1 17
+20 1 0 3 2 17 18
+21 1 0 3 2 18 19
+22 1 0 3 2 19 6
+23 1 0 4 2 6 20
+24 1 0 4 2 20 21
+25 1 0 4 2 21 22
+26 1 0 4 2 22 23
+27 1 0 4 2 23 24
+28 1 0 4 2 24 25
+29 1 0 4 2 25 26
+30 1 0 4 2 26 3
+31 1 0 5 2 3 27
+32 1 0 5 2 27 28
+33 1 0 5 2 28 29
+34 1 0 5 2 29 30
+35 1 0 5 2 30 31
+36 1 0 5 2 31 32
+37 1 0 5 2 32 33
+38 1 0 5 2 33 7
+39 1 0 6 2 7 34
+40 1 0 6 2 34 35
+41 1 0 6 2 35 36
+42 1 0 6 2 36 4
+43 1 0 7 2 4 37
+44 1 0 7 2 37 38
+45 1 0 7 2 38 39
+46 1 0 7 2 39 8
+47 1 0 8 2 8 40
+48 1 0 8 2 40 41
+49 1 0 8 2 41 42
+50 1 0 8 2 42 9
+51 1 0 9 2 9 43
+52 1 0 9 2 43 44
+53 1 0 9 2 44 45
+54 1 0 9 2 45 46
+55 1 0 9 2 46 47
+56 1 0 9 2 47 48
+57 1 0 9 2 48 49
+58 1 0 9 2 49 10
+$ENDELM
diff --git a/tests/grid/grid_in_gmsh_04.output b/tests/grid/grid_in_gmsh_04.output
new file mode 100644 (file)
index 0000000..5b71810
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::Mesh info:
+DEAL:: dimension: 1
+DEAL:: no. of cells: 48
+DEAL:: boundary indicators: 
+DEAL:: written to grid_in_gmsh_04.vtu
+DEAL::

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.