]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Allow copying manifold ids in GridGenerator::merge_triangulations
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Sat, 7 Jul 2018 23:30:54 +0000 (01:30 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Wed, 11 Jul 2018 08:37:10 +0000 (10:37 +0200)
doc/news/changes/minor/20180708DanielArndt [new file with mode: 0644]
include/deal.II/grid/grid_generator.h
source/grid/grid_generator.cc
source/grid/grid_generator.inst.in
source/grid/grid_tools.cc
source/grid/tria.cc

diff --git a/doc/news/changes/minor/20180708DanielArndt b/doc/news/changes/minor/20180708DanielArndt
new file mode 100644 (file)
index 0000000..af38198
--- /dev/null
@@ -0,0 +1,3 @@
+New: GridGenerator_merge_triangulations gained an option to copy manifold ids.
+<br>
+(Daniel Arndt, 2018/07/08)
index e9f7e9f214feff0dcef90e3766ca86ab8ce21d58..b801236da0f465df02f1ae7fbc11ffdb1a297310 100644 (file)
@@ -984,11 +984,11 @@ namespace GridGenerator
    * refined cells.
    *
    * @note The function copies the material ids of the cells of the two input
-   * triangulations into the output triangulation but it currently makes no
-   * attempt to do the same for boundary ids. In other words, if the two
-   * coarse meshes have anything but the default boundary indicators, then you
-   * will currently have to set boundary indicators again by hand in the
-   * output triangulation.
+   * triangulations into the output triangulation. If @p copy_manifold_ids is
+   * activated, manifold ids will be copied. Currently, boundary_ids are never
+   * copied. In other words, if the two coarse meshes have anything but the
+   * default boundary indicators, then you will currently have to set boundary
+   * indicators again by hand in the output triangulation.
    *
    * @note Unlike most GridGenerator functions, this function does not attach
    * any manifolds to @p result, nor does it set any manifold ids.
@@ -1002,7 +1002,8 @@ namespace GridGenerator
   merge_triangulations(const Triangulation<dim, spacedim> &triangulation_1,
                        const Triangulation<dim, spacedim> &triangulation_2,
                        Triangulation<dim, spacedim> &      result,
-                       const double duplicated_vertex_tolerance = 1.0e-12);
+                       const double duplicated_vertex_tolerance = 1.0e-12,
+                       const bool   copy_manifold_ids           = false);
 
   /**
    * Given the two triangulations specified as the first two arguments, create
index 94f28d20fa2f07fcf163f98a9bcdca4fe1915a6b..3a55b2cf04894bc9381993ac7e31bb19e80c86d1 100644 (file)
@@ -3793,7 +3793,8 @@ namespace GridGenerator
   merge_triangulations(const Triangulation<dim, spacedim> &triangulation_1,
                        const Triangulation<dim, spacedim> &triangulation_2,
                        Triangulation<dim, spacedim> &      result,
-                       const double duplicated_vertex_tolerance)
+                       const double duplicated_vertex_tolerance,
+                       const bool   copy_manifold_ids)
   {
     Assert(triangulation_1.n_levels() == 1,
            ExcMessage("The input triangulations must be coarse meshes."));
@@ -3809,36 +3810,121 @@ namespace GridGenerator
     // now form the union of the set of cells
     std::vector<CellData<dim>> cells;
     cells.reserve(triangulation_1.n_cells() + triangulation_2.n_cells());
-    for (typename Triangulation<dim, spacedim>::cell_iterator cell =
-           triangulation_1.begin();
-         cell != triangulation_1.end();
-         ++cell)
+    for (const auto &cell : triangulation_1.cell_iterators())
       {
         CellData<dim> this_cell;
         for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
           this_cell.vertices[v] = cell->vertex_index(v);
         this_cell.material_id = cell->material_id();
-        cells.push_back(this_cell);
+        this_cell.manifold_id = cell->manifold_id();
+        cells.push_back(std::move(this_cell));
       }
 
     // now do the same for the other other mesh. note that we have to
     // translate the vertex indices
-    for (typename Triangulation<dim, spacedim>::cell_iterator cell =
-           triangulation_2.begin();
-         cell != triangulation_2.end();
-         ++cell)
+    for (const auto &cell : triangulation_2.cell_iterators())
       {
         CellData<dim> this_cell;
         for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
           this_cell.vertices[v] =
             cell->vertex_index(v) + triangulation_1.n_vertices();
         this_cell.material_id = cell->material_id();
-        cells.push_back(this_cell);
+        this_cell.manifold_id = cell->manifold_id();
+        cells.push_back(std::move(this_cell));
+      }
+
+    // Now for the SubCellData
+    SubCellData subcell_data;
+    if (copy_manifold_ids)
+      {
+        if (dim == 2)
+          {
+            subcell_data.boundary_lines.reserve(triangulation_1.n_lines() +
+                                                triangulation_2.n_lines());
+
+            auto copy_line_manifold_ids =
+              [&](const Triangulation<dim, spacedim> &tria,
+                  const unsigned int                  offset) {
+                for (typename Triangulation<dim, spacedim>::line_iterator line =
+                       tria.begin_face();
+                     line != tria.end_face();
+                     ++line)
+                  if (line->manifold_id() != numbers::flat_manifold_id)
+                    {
+                      CellData<1> boundary_line;
+                      boundary_line.vertices[0] =
+                        line->vertex_index(0) + offset;
+                      boundary_line.vertices[1] =
+                        line->vertex_index(1) + offset;
+                      boundary_line.manifold_id = line->manifold_id();
+                      subcell_data.boundary_lines.push_back(
+                        std::move(boundary_line));
+                    }
+              };
+
+            copy_line_manifold_ids(triangulation_1, 0);
+            copy_line_manifold_ids(triangulation_2,
+                                   triangulation_1.n_vertices());
+          }
+
+        if (dim == 3)
+          {
+            subcell_data.boundary_quads.reserve(triangulation_1.n_quads() +
+                                                triangulation_2.n_quads());
+            // we can't do better here than to loop over all the lines bounding
+            // a face. For regular meshes an (interior) line in 3D is part of
+            // four cells. So this should be an appropriate guess.
+            subcell_data.boundary_lines.reserve(triangulation_1.n_cells() * 4 +
+                                                triangulation_2.n_cells() * 4);
+
+            auto copy_face_and_line_manifold_ids =
+              [&](const Triangulation<dim, spacedim> &tria,
+                  const unsigned int                  offset) {
+                for (typename Triangulation<dim, spacedim>::face_iterator face =
+                       tria.begin_face();
+                     face != tria.end_face();
+                     ++face)
+                  if (face->manifold_id() != numbers::flat_manifold_id)
+                    {
+                      CellData<2> boundary_quad;
+                      boundary_quad.vertices[0] =
+                        face->vertex_index(0) + offset;
+                      boundary_quad.vertices[1] =
+                        face->vertex_index(1) + offset;
+                      boundary_quad.vertices[2] =
+                        face->vertex_index(2) + offset;
+                      boundary_quad.vertices[3] =
+                        face->vertex_index(3) + offset;
+
+                      boundary_quad.manifold_id = face->manifold_id();
+                      subcell_data.boundary_quads.push_back(
+                        std::move(boundary_quad));
+                    }
+                for (const auto &cell : tria.cell_iterators())
+                  for (unsigned int l = 0; l < 12; ++l)
+                    if (cell->line(l)->manifold_id() !=
+                        numbers::flat_manifold_id)
+                      {
+                        CellData<1> boundary_line;
+                        boundary_line.vertices[0] =
+                          cell->line(l)->vertex_index(0) + offset;
+                        boundary_line.vertices[1] =
+                          cell->line(l)->vertex_index(1) + offset;
+                        boundary_line.manifold_id =
+                          cell->line(l)->manifold_id();
+                        subcell_data.boundary_lines.push_back(
+                          std::move(boundary_line));
+                      }
+              };
+
+            copy_face_and_line_manifold_ids(triangulation_1, 0);
+            copy_face_and_line_manifold_ids(triangulation_2,
+                                            triangulation_1.n_vertices());
+          }
       }
 
     // throw out duplicated vertices from the two meshes, reorder vertices as
     // necessary and create the triangulation
-    SubCellData               subcell_data;
     std::vector<unsigned int> considered_vertices;
     GridTools::delete_duplicated_vertices(vertices,
                                           cells,
index f581285c88720901ca0387b584d4d4e577343bb3..b64c8f60a4a8d21f6b045b2fd96db50489da2ebe 100644 (file)
@@ -76,7 +76,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
         const Triangulation<deal_II_dimension, deal_II_space_dimension>
           &triangulation_2,
         Triangulation<deal_II_dimension, deal_II_space_dimension> &result,
-        const double duplicated_vertex_tolerance);
+        const double duplicated_vertex_tolerance,
+        const bool   copy_manifold_ids);
 
       template void
       create_union_triangulation(
index a62d82612ce7d8917feca6a198f177231d9c5dd6..188aafd91e3e02d3daef7db1c5fcc5aaf8f584df 100644 (file)
@@ -586,16 +586,19 @@ namespace GridTools
           }
       }
 
-    // now we got a renumbering list. simply
-    // renumber all vertices (non-duplicate
-    // vertices get renumbered to themselves, so
-    // nothing bad happens). after that, the
-    // duplicate vertices will be unused, so call
-    // delete_unused_vertices() to do that part
-    // of the job.
-    for (unsigned int c = 0; c < cells.size(); ++c)
-      for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
-        cells[c].vertices[v] = new_vertex_numbers[cells[c].vertices[v]];
+    // now we got a renumbering list. simply renumber all vertices
+    // (non-duplicate vertices get renumbered to themselves, so nothing bad
+    // happens). after that, the duplicate vertices will be unused, so call
+    // delete_unused_vertices() to do that part of the job.
+    for (auto &cell : cells)
+      for (auto &vertex_index : cell.vertices)
+        vertex_index = new_vertex_numbers[vertex_index];
+    for (auto &quad : subcelldata.boundary_quads)
+      for (auto &vertex_index : quad.vertices)
+        vertex_index = new_vertex_numbers[vertex_index];
+    for (auto &line : subcelldata.boundary_lines)
+      for (auto &vertex_index : line.vertices)
+        vertex_index = new_vertex_numbers[vertex_index];
 
     delete_unused_vertices(vertices, cells, subcelldata);
   }
index c9b2e2bfd0c58675dcf8784e09fd1a6fcfbc3377..90a1bf268f5c6f67f9f494f7a4e89f349233dcad 100644 (file)
@@ -3106,28 +3106,25 @@ namespace internal
                               ExcLineInexistant(line_vertices.first,
                                                 line_vertices.second));
               }
-            // Assert that only exterior
-            // lines are given a boundary
-            // indicator
-            AssertThrow(
-              line->at_boundary(),
-              ExcInteriorLineCantBeBoundary(line->vertex_index(0),
-                                            line->vertex_index(1),
-                                            boundary_line->boundary_id));
-
-            // and make sure that we don't
-            // attempt to reset the
-            // boundary indicator to a
-            // different than the
-            // previously set value
-            if (line->boundary_id() != 0)
-              AssertThrow(line->boundary_id() == boundary_line->boundary_id,
-                          ExcMessage(
-                            "Duplicate boundary lines are only allowed "
-                            "if they carry the same boundary indicator."));
+            // Only exterior lines can be given a boundary indicator
+            if (line->at_boundary())
+              {
+                // make sure that we don't attempt to reset the boundary
+                // indicator to a different than the previously set value
+                if (line->boundary_id() != 0)
+                  AssertThrow(line->boundary_id() == boundary_line->boundary_id,
+                              ExcMessage(
+                                "Duplicate boundary lines are only allowed "
+                                "if they carry the same boundary indicator."));
 
-            line->set_boundary_id_internal(boundary_line->boundary_id);
+                line->set_boundary_id_internal(boundary_line->boundary_id);
+              }
             // Set manifold id if given
+            if (line->manifold_id() != numbers::flat_manifold_id)
+              AssertThrow(line->manifold_id() == boundary_line->manifold_id,
+                          ExcMessage(
+                            "Duplicate boundary lines are only allowed "
+                            "if they carry the same manifold indicator."));
             line->set_manifold_id(boundary_line->manifold_id);
           }
 
@@ -3262,26 +3259,25 @@ namespace internal
 
             // check whether this face is
             // really an exterior one
-            AssertThrow(
-              quad->at_boundary(),
-              ExcInteriorQuadCantBeBoundary(quad->vertex_index(0),
-                                            quad->vertex_index(1),
-                                            quad->vertex_index(2),
-                                            quad->vertex_index(3),
-                                            boundary_quad->boundary_id));
-
-            // and make sure that we don't
-            // attempt to reset the
-            // boundary indicator to a
-            // different than the
-            // previously set value
-            if (quad->boundary_id() != 0)
-              AssertThrow(quad->boundary_id() == boundary_quad->boundary_id,
+            if (quad->at_boundary())
+              {
+                // and make sure that we don't attempt to reset the boundary
+                // indicator to a different than the previously set value
+                if (quad->boundary_id() != 0)
+                  AssertThrow(quad->boundary_id() == boundary_quad->boundary_id,
+                              ExcMessage(
+                                "Duplicate boundary quads are only allowed "
+                                "if they carry the same boundary indicator."));
+
+                quad->set_boundary_id_internal(boundary_quad->boundary_id);
+              }
+            // Set manifold id if given
+            if (quad->manifold_id() != numbers::flat_manifold_id)
+              AssertThrow(quad->manifold_id() == boundary_quad->manifold_id,
                           ExcMessage(
                             "Duplicate boundary quads are only allowed "
-                            "if they carry the same boundary indicator."));
+                            "if they carry the same manifold indicator."));
 
-            quad->set_boundary_id_internal(boundary_quad->boundary_id);
             quad->set_manifold_id(boundary_quad->manifold_id);
           }
 

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.