]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix bugs in extract_boundary_mesh.
authorSebastian Stark <stark.sebastian@gmx.de>
Mon, 26 Nov 2018 08:27:18 +0000 (10:27 +0200)
committerSebastian Stark <stark.sebastian@gmx.de>
Tue, 8 Jan 2019 19:07:55 +0000 (21:07 +0200)
In 3d, extract_boundary_mesh swaps vertices to get consistent boundary cell orientations. However, this swapping was so far not taken into account during the assembly of the map between boundary cells and corresponding faces of the volume cells. Fixing the latter required substantial changes (a quick and dirty fix would have been possible in principle, but certainly not clean) to extract_boundary_mesh because the (local) face indices of the volume cells must be known when the map is assembled. Also, the algorithm to do the refinement of the boundary mesh has been modified such that map lookups are no longer necessary. Moreover, vertex swapping has been introduced also for the 2d case in order to make the 2d behavior consistent with the 3d behavior (in previous versions of the function, the normal vector of the boundary mesh pointed always into the domain in 3d, whereas it pointed either into or out of the domain in 2d depending on the situation). Further, the normals have been adjusted to generally face outward.  Fixes #7467.

include/deal.II/grid/grid_generator.h
source/grid/grid_generator.cc

index 661d9cc5771beef9adb55875a0d2fee16a74de9e..354037d6f20ac9f909e472193388754e988852cd 100644 (file)
@@ -1579,13 +1579,12 @@ namespace GridGenerator
    * @return A map that for each cell of the surface mesh (key) returns an
    * iterator to the corresponding face of a cell of the volume mesh (value).
    * The keys include both active and non-active cells of the surface mesh.
-   * For dim=2 (i.e., where volume cells are quadrilaterals and surface cells
-   * are lines), the order of vertices of surface cells and the corresponding
-   * volume faces match. For dim=3 (i.e., where volume cells are hexahedra and
-   * surface cells are quadrilaterals), the order of vertices may not match in
-   * order to ensure that each surface cell has a right-handed coordinate
-   * system when viewed from one of the two sides of the surface connecting
-   * the cells of the surface mesh.
+   * The order of vertices of surface cells and the corresponding
+   * volume faces may not match in order to ensure that each surface cell is
+   * associated with an outward facing normal.
+   * As a consequence, if you want to match quantities on the faces of the
+   * domain cells and on the cells of the surface mesh, you may have to
+   * translate between vertex locations or quadrature points.
    *
    * @note The algorithm outlined above assumes that all faces on higher
    * refinement levels always have exactly the same boundary indicator as
index 44e8a50ed0b19a4ed42d739c65d4b412ac8691bf..97e567900d4c84c9400a79212bb515e8498ee9ec 100644 (file)
@@ -5358,27 +5358,64 @@ namespace GridGenerator
     //    preserve the order of cells passed in using the CellData argument;
     //    also, that it will not reorder the vertices.
 
-    std::map<typename MeshType<dim - 1, spacedim>::cell_iterator,
-             typename MeshType<dim, spacedim>::face_iterator>
-      surface_to_volume_mapping;
-
-    const unsigned int boundary_dim = dim - 1; // dimension of the boundary mesh
-
-    // First create surface mesh and mapping
-    // from only level(0) cells of volume_mesh
-    std::vector<typename MeshType<dim, spacedim>::face_iterator>
-      mapping; // temporary map for level==0
-
-
+    // dimension of the boundary mesh
+    const unsigned int boundary_dim = dim - 1;
+
+    // temporary map for level==0
+    // iterator to face is stored along with face number
+    // (this is required by the algorithm to adjust the normals of the
+    // cells of the boundary mesh)
+    std::vector<
+      std::pair<typename MeshType<dim, spacedim>::face_iterator, unsigned int>>
+      temporary_mapping_level0;
+
+    // vector indicating whether a vertex of the volume mesh has
+    // already been visited (necessary to avoid duplicate vertices in
+    // boundary mesh)
     std::vector<bool> touched(volume_mesh.get_triangulation().n_vertices(),
                               false);
+
+    // data structures required for creation of boundary mesh
     std::vector<CellData<boundary_dim>> cells;
     SubCellData                         subcell_data;
     std::vector<Point<spacedim>>        vertices;
 
-    std::map<unsigned int, unsigned int>
-      map_vert_index; // volume vertex indices to surf ones
+    // volume vertex indices to surf ones
+    std::map<unsigned int, unsigned int> map_vert_index;
+
+    // define swapping of vertices to get proper normal orientation of boundary
+    // mesh;
+    // the entry (i,j) of swap_matrix stores the index of the vertex of
+    // the boundary cell corresponding to the j-th vertex on the i-th face
+    // of the underlying volume cell
+    // if e.g. face 3 of a volume cell is considered and vertices 1 and 2 of the
+    // corresponding boundary cell are swapped to get
+    // proper normal orientation, swap_matrix[3]=( 0, 2, 1, 3 )
+    Table<2, unsigned int> swap_matrix(
+      GeometryInfo<spacedim>::faces_per_cell,
+      GeometryInfo<dim - 1>::vertices_per_cell);
+    for (unsigned int i1 = 0; i1 < GeometryInfo<spacedim>::faces_per_cell; i1++)
+      {
+        for (unsigned int i2 = 0; i2 < GeometryInfo<dim - 1>::vertices_per_cell;
+             i2++)
+          swap_matrix[i1][i2] = i2;
+      }
+    // vertex swapping such that normals on the surface mesh point out of the
+    // underlying volume
+    if (dim == 3)
+      {
+        std::swap(swap_matrix[0][1], swap_matrix[0][2]);
+        std::swap(swap_matrix[2][1], swap_matrix[2][2]);
+        std::swap(swap_matrix[4][1], swap_matrix[4][2]);
+      }
+    else if (dim == 2)
+      {
+        std::swap(swap_matrix[1][0], swap_matrix[1][1]);
+        std::swap(swap_matrix[2][0], swap_matrix[2][1]);
+      }
 
+    // Create boundary mesh and mapping
+    // from only level(0) cells of volume_mesh
     for (typename MeshType<dim, spacedim>::cell_iterator cell =
            volume_mesh.begin(0);
          cell != volume_mesh.end(0);
@@ -5407,26 +5444,14 @@ namespace GridGenerator
                       touched[v_index]        = true;
                     }
 
-                  c_data.vertices[j] = map_vert_index[v_index];
-                  c_data.material_id =
-                    static_cast<types::material_id>(face->boundary_id());
-                  c_data.manifold_id = face->manifold_id();
+                  c_data.vertices[swap_matrix[i][j]] = map_vert_index[v_index];
                 }
+              c_data.material_id =
+                static_cast<types::material_id>(face->boundary_id());
+              c_data.manifold_id = face->manifold_id();
 
-              // if we start from a 3d mesh, then we have copied the
-              // vertex information in the same order in which they
-              // appear in the face; however, this means that we
-              // impart a coordinate system that is right-handed when
-              // looked at *from the outside* of the cell if the
-              // current face has index 0, 2, 4 within a 3d cell, but
-              // right-handed when looked at *from the inside* for the
-              // other faces. we fix this by flipping opposite
-              // vertices if we are on a face 1, 3, 5
-              if (dim == 3)
-                if (i % 2 == 1)
-                  std::swap(c_data.vertices[1], c_data.vertices[2]);
 
-              // in 3d, we also need to make sure we copy the manifold
+              // in 3d, we need to make sure we copy the manifold
               // indicators from the edges of the volume mesh to the
               // edges of the surface mesh
               //
@@ -5459,7 +5484,8 @@ namespace GridGenerator
                             break;
                           }
                       if (edge_found == true)
-                        continue; // try next edge of current face
+                        // try next edge of current face
+                        continue;
                     }
 
                     CellData<1> edge;
@@ -5473,9 +5499,8 @@ namespace GridGenerator
                     subcell_data.boundary_lines.push_back(edge);
                   }
 
-
               cells.push_back(c_data);
-              mapping.push_back(face);
+              temporary_mapping_level0.push_back(std::make_pair(face, i));
             }
         }
 
@@ -5485,48 +5510,119 @@ namespace GridGenerator
       surface_mesh.get_triangulation())
       .create_triangulation(vertices, cells, subcell_data);
 
-    // Make the actual mapping
-    for (typename MeshType<dim - 1, spacedim>::active_cell_iterator cell =
-           surface_mesh.begin(0);
-         cell != surface_mesh.end(0);
-         ++cell)
-      surface_to_volume_mapping[cell] = mapping.at(cell->index());
+    // in 2d: set default boundary ids for "boundary vertices"
+    if (dim == 2)
+      {
+        for (const auto &cell : surface_mesh.active_cell_iterators())
+          for (unsigned int vertex = 0; vertex < 2; vertex++)
+            if (cell->face(vertex)->at_boundary())
+              cell->face(vertex)->set_boundary_id(0);
+      }
 
+    // Make mapping for level 0
+
+    // temporary map between cells on the boundary and corresponding faces of
+    // domain mesh (each face is characterized by an iterator to the face and
+    // the face number within the underlying cell)
+    std::vector<std::pair<
+      const typename MeshType<dim - 1, spacedim>::cell_iterator,
+      std::pair<typename MeshType<dim, spacedim>::face_iterator, unsigned int>>>
+      temporary_map_boundary_cell_face;
+    for (const auto &cell : surface_mesh.active_cell_iterators())
+      temporary_map_boundary_cell_face.push_back(
+        std::make_pair(cell, temporary_mapping_level0.at(cell->index())));
+
+
+    // refine the boundary mesh according to the refinement of the underlying
+    // volume mesh,
+    // algorithm:
+    //   (1) check which cells on refinement level i need to be refined
+    //   (2) do refinement (yields cells on level i+1)
+    //   (3) repeat for the next level (i+1->i) until refinement is completed
+
+    // stores the index into temporary_map_boundary_cell_face at which
+    // presently deepest refinement level of boundary mesh begins
+    unsigned int index_cells_deepest_level = 0;
     do
       {
         bool changed = false;
 
-        for (typename MeshType<dim - 1, spacedim>::active_cell_iterator cell =
-               surface_mesh.begin_active();
-             cell != surface_mesh.end();
-             ++cell)
-          if (surface_to_volume_mapping[cell]->has_children() == true)
-            {
-              cell->set_refine_flag();
-              changed = true;
-            }
+        // vector storing cells which have been marked for
+        // refinement
+        std::vector<unsigned int> cells_refined;
+
+        // loop over cells of presently deepest level of boundary triangulation
+        for (unsigned int cell_n = index_cells_deepest_level;
+             cell_n < temporary_map_boundary_cell_face.size();
+             cell_n++)
+          {
+            // mark boundary cell for refinement if underlying volume face has
+            // children
+            if (temporary_map_boundary_cell_face[cell_n]
+                  .second.first->has_children())
+              {
+                // algorithm only works for
+                // isotropic refinement!
+                Assert(temporary_map_boundary_cell_face[cell_n]
+                           .second.first->refinement_case() ==
+                         RefinementCase<dim - 1>::isotropic_refinement,
+                       ExcNotImplemented());
+                temporary_map_boundary_cell_face[cell_n]
+                  .first->set_refine_flag();
+                cells_refined.push_back(cell_n);
+                changed = true;
+              }
+          }
 
+        // if cells have been marked for refinement (i.e., presently deepest
+        // level is not the deepest level of the volume mesh)
         if (changed)
           {
+            // do actual refinement
             const_cast<Triangulation<dim - 1, spacedim> &>(
               surface_mesh.get_triangulation())
               .execute_coarsening_and_refinement();
 
-            for (typename MeshType<dim - 1, spacedim>::cell_iterator
-                   surface_cell = surface_mesh.begin();
-                 surface_cell != surface_mesh.end();
-                 ++surface_cell)
-              for (unsigned int c = 0; c < surface_cell->n_children(); c++)
-                if (surface_to_volume_mapping.find(surface_cell->child(c)) ==
-                    surface_to_volume_mapping.end())
-                  surface_to_volume_mapping[surface_cell->child(c)] =
-                    surface_to_volume_mapping[surface_cell]->child(c);
+            // add new level of cells to temporary_map_boundary_cell_face
+            index_cells_deepest_level = temporary_map_boundary_cell_face.size();
+            for (const auto &refined_cell_n : cells_refined)
+              {
+                const typename MeshType<dim - 1, spacedim>::cell_iterator
+                  refined_cell =
+                    temporary_map_boundary_cell_face[refined_cell_n].first;
+                const typename MeshType<dim,
+                                        spacedim>::face_iterator refined_face =
+                  temporary_map_boundary_cell_face[refined_cell_n].second.first;
+                const unsigned int refined_face_number =
+                  temporary_map_boundary_cell_face[refined_cell_n]
+                    .second.second;
+                for (unsigned int child_n = 0;
+                     child_n < refined_cell->n_children();
+                     ++child_n)
+                  // at this point, the swapping of vertices done earlier must
+                  // be taken into account to get the right association between
+                  // volume faces and boundary cells!
+                  temporary_map_boundary_cell_face.push_back(
+                    std::make_pair(refined_cell->child(
+                                     swap_matrix[refined_face_number][child_n]),
+                                   std::make_pair(refined_face->child(child_n),
+                                                  refined_face_number)));
+              }
           }
+        // we are at the deepest level of refinement of the volume mesh
         else
           break;
       }
     while (true);
 
+    // generate the final mapping from the temporary mapping
+    std::map<typename MeshType<dim - 1, spacedim>::cell_iterator,
+             typename MeshType<dim, spacedim>::face_iterator>
+      surface_to_volume_mapping;
+    for (unsigned int i = 0; i < temporary_map_boundary_cell_face.size(); i++)
+      surface_to_volume_mapping[temporary_map_boundary_cell_face[i].first] =
+        temporary_map_boundary_cell_face[i].second.first;
+
     return surface_to_volume_mapping;
   }
 

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.