]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Map boundary to bulk dof iterators.
authorLuca Heltai <luca.heltai@sissa.it>
Wed, 12 Apr 2023 08:15:38 +0000 (11:15 +0300)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 4 Jul 2023 07:48:28 +0000 (09:48 +0200)
doc/news/changes/minor/20230412LucaHeltai [new file with mode: 0644]
include/deal.II/dofs/dof_tools.h
source/dofs/dof_tools.cc
source/dofs/dof_tools.inst.in
tests/dofs/map_boundary_to_bulk_dof_iterators_01.cc [new file with mode: 0644]
tests/dofs/map_boundary_to_bulk_dof_iterators_01.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20230412LucaHeltai b/doc/news/changes/minor/20230412LucaHeltai
new file mode 100644 (file)
index 0000000..5633177
--- /dev/null
@@ -0,0 +1,7 @@
+New: Added a function DoFTools::map_boundary_to_bulk_dof_iterators() that 
+generates a mapping of codimension-1 active DoFHandler cell iterators to
+codimension-0 cells and face indices, to couple DoFHandler objects of different
+co-dimensions, initialized on grids generated with
+GridTools::extract_boundary_mesh()
+<br>
+(Luca Heltai, 2023/04/12)
index fb0397d8db9a80272fcd5f36e261c7c3b343f298..22d954c0e64c098fd9beec5c79df817de9d4a4f9 100644 (file)
@@ -92,120 +92,133 @@ namespace DoFTools
 #endif
 
 /**
- * This is a collection of functions operating on, and manipulating the
- * numbers of degrees of freedom. The documentation of the member functions
- * will provide more information, but for functions that exist in multiple
- * versions, there are sections in this global documentation stating some
- * commonalities.
+ * This is a collection of functions operating on, and manipulating the numbers
+ * of degrees of freedom. The documentation of the member functions will provide
+ * more information, but for functions that exist in multiple versions, there
+ * are sections in this global documentation stating some commonalities.
  *
  * <h3>Setting up sparsity patterns</h3>
  *
- * When assembling system matrices, the entries are usually of the form
- * $a_{ij} = a(\phi_i, \phi_j)$, where $a$ is a bilinear functional, often an
- * integral. When using sparse matrices, we therefore only need to reserve
- * space for those $a_{ij}$ only, which are nonzero, which is the same as to
- * say that the basis functions $\phi_i$ and $\phi_j$ have a nonempty
- * intersection of their support. Since the support of basis functions is
- * bound only on cells on which they are located or to which they are
- * adjacent, to determine the sparsity pattern it is sufficient to loop over
- * all cells and connect all basis functions on each cell with all other basis
- * functions on that cell.  There may be finite elements for which not all
- * basis functions on a cell connect with each other, but no use of this case
- * is made since no examples where this occurs are known to the author.
+ * When assembling system matrices, the entries are usually of the form $a_{ij}
+ * = a(\phi_i, \phi_j)$, where $a$ is a bilinear functional, often an integral.
+ * When using sparse matrices, we therefore only need to reserve space for those
+ * $a_{ij}$ only, which are nonzero, which is the same as to say that the basis
+ * functions $\phi_i$ and $\phi_j$ have a nonempty intersection of their
+ * support. Since the support of basis functions is bound only on cells on which
+ * they are located or to which they are adjacent, to determine the sparsity
+ * pattern it is sufficient to loop over all cells and connect all basis
+ * functions on each cell with all other basis functions on that cell.  There
+ * may be finite elements for which not all basis functions on a cell connect
+ * with each other, but no use of this case is made since no examples where this
+ * occurs are known to the author.
  *
  *
  * <h3>DoF numberings on boundaries</h3>
  *
- * When projecting the traces of functions to the boundary or parts thereof,
- * one needs to build matrices and vectors that act only on those degrees of
- * freedom that are located on the boundary, rather than on all degrees of
- * freedom. One could do that by simply building matrices in which the entries
- * for all interior DoFs are zero, but such matrices are always very rank
- * deficient and not very practical to work with.
+ * When projecting the traces of functions to the boundary or parts thereof, one
+ * needs to build matrices and vectors that act only on those degrees of freedom
+ * that are located on the boundary, rather than on all degrees of freedom. One
+ * could do that by simply building matrices in which the entries for all
+ * interior DoFs are zero, but such matrices are always very rank deficient and
+ * not very practical to work with.
  *
- * What is needed instead in this case is a numbering of the boundary degrees
- * of freedom, i.e. we should enumerate all the degrees of freedom that are
- * sitting on the boundary, and exclude all other (interior) degrees of
- * freedom. The map_dof_to_boundary_indices() function does exactly this: it
- * provides a vector with as many entries as there are degrees of freedom on
- * the whole domain, with each entry being the number in the numbering of the
- * boundary or numbers::invalid_dof_index if the dof is not on the
- * boundary.
+ * What is needed instead in this case is a numbering of the boundary degrees of
+ * freedom, i.e. we should enumerate all the degrees of freedom that are sitting
+ * on the boundary, and exclude all other (interior) degrees of freedom. The
+ * map_dof_to_boundary_indices() function does exactly this: it provides a
+ * vector with as many entries as there are degrees of freedom on the whole
+ * domain, with each entry being the number in the numbering of the boundary or
+ * numbers::invalid_dof_index if the dof is not on the boundary.
  *
  * With this vector, one can get, for any given degree of freedom, a unique
  * number among those DoFs that sit on the boundary; or, if your DoF was
- * interior to the domain, the result would be numbers::invalid_dof_index.
- * We need this mapping, for example, to build the @ref GlossMassMatrix "mass matrix" on the boundary
- * (for this, see make_boundary_sparsity_pattern() function, the corresponding
- * section below, as well as the MatrixCreator namespace documentation).
+ * interior to the domain, the result would be numbers::invalid_dof_index. We
+ * need this mapping, for example, to build the @ref GlossMassMatrix "mass
+ * matrix" on the boundary (for this, see make_boundary_sparsity_pattern()
+ * function, the corresponding section below, as well as the MatrixCreator
+ * namespace documentation).
  *
  * Actually, there are two map_dof_to_boundary_indices() functions, one
- * producing a numbering for all boundary degrees of freedom and one producing
- * a numbering for only parts of the boundary, namely those parts for which
- * the boundary indicator is listed in a set of indicators given to the
- * function. The latter case is needed if, for example, we would only want to
- * project the boundary values for the Dirichlet part of the boundary. You
- * then give the function a list of boundary indicators referring to Dirichlet
- * parts on which the projection is to be performed. The parts of the boundary
- * on which you want to project need not be contiguous; however, it is not
- * guaranteed that the indices of each of the boundary parts are continuous,
- * i.e. the indices of degrees of freedom on different parts may be
- * intermixed.
+ * producing a numbering for all boundary degrees of freedom and one producing a
+ * numbering for only parts of the boundary, namely those parts for which the
+ * boundary indicator is listed in a set of indicators given to the function.
+ * The latter case is needed if, for example, we would only want to project the
+ * boundary values for the Dirichlet part of the boundary. You then give the
+ * function a list of boundary indicators referring to Dirichlet parts on which
+ * the projection is to be performed. The parts of the boundary on which you
+ * want to project need not be contiguous; however, it is not guaranteed that
+ * the indices of each of the boundary parts are continuous, i.e. the indices of
+ * degrees of freedom on different parts may be intermixed.
  *
  * Degrees of freedom on the boundary but not on one of the specified boundary
- * parts are given the index numbers::invalid_dof_index, as if they were in
- * the interior. If no boundary indicator was given or if no face of a cell
- * has a boundary indicator contained in the given list, the vector of new
- * indices consists solely of numbers::invalid_dof_index.
+ * parts are given the index numbers::invalid_dof_index, as if they were in the
+ * interior. If no boundary indicator was given or if no face of a cell has a
+ * boundary indicator contained in the given list, the vector of new indices
+ * consists solely of numbers::invalid_dof_index.
  *
  * (As a side note, for corner cases: The question what a degree of freedom on
- * the boundary is, is not so easy.  It should really be a degree of freedom
- * of which the respective basis function has nonzero values on the boundary.
- * At least for Lagrange elements this definition is equal to the statement
- * that the off-point, or what deal.II calls support_point, of the shape
- * function, i.e. the point where the function assumes its nominal value (for
- * Lagrange elements this is the point where it has the function value 1), is
- * located on the boundary. We do not check this directly, the criterion is
- * rather defined through the information the finite element class gives: the
- * FiniteElement class defines the numbers of basis functions per vertex, per
- * line, and so on and the basis functions are numbered after this
- * information; a basis function is to be considered to be on the face of a
- * cell (and thus on the boundary if the cell is at the boundary) according to
- * it belonging to a vertex, line, etc but not to the interior of the cell.
- * The finite element uses the same cell-wise numbering so that we can say
- * that if a degree of freedom was numbered as one of the dofs on lines, we
- * assume that it is located on the line. Where the off-point actually is, is
- * a secret of the finite element (well, you can ask it, but we don't do it
- * here) and not relevant in this context.)
+ * the boundary is, is not so easy.  It should really be a degree of freedom of
+ * which the respective basis function has nonzero values on the boundary. At
+ * least for Lagrange elements this definition is equal to the statement that
+ * the off-point, or what deal.II calls support_point, of the shape function,
+ * i.e. the point where the function assumes its nominal value (for Lagrange
+ * elements this is the point where it has the function value 1), is located on
+ * the boundary. We do not check this directly, the criterion is rather defined
+ * through the information the finite element class gives: the FiniteElement
+ * class defines the numbers of basis functions per vertex, per line, and so on
+ * and the basis functions are numbered after this information; a basis function
+ * is to be considered to be on the face of a cell (and thus on the boundary if
+ * the cell is at the boundary) according to it belonging to a vertex, line, etc
+ * but not to the interior of the cell. The finite element uses the same
+ * cell-wise numbering so that we can say that if a degree of freedom was
+ * numbered as one of the dofs on lines, we assume that it is located on the
+ * line. Where the off-point actually is, is a secret of the finite element
+ * (well, you can ask it, but we don't do it here) and not relevant in this
+ * context.)
  *
  *
  * <h3>Setting up sparsity patterns for boundary matrices</h3>
  *
- * In some cases, one wants to only work with DoFs that sit on the boundary.
- * One application is, for example, if rather than interpolating non-
- * homogeneous boundary values, one would like to project them. For this, we
- * need two things: a way to identify nodes that are located on (parts of) the
- * boundary, and a way to build matrices out of only degrees of freedom that
- * are on the boundary (i.e. much smaller matrices, in which we do not even
- * build the large zero block that stems from the fact that most degrees of
- * freedom have no support on the boundary of the domain). The first of these
- * tasks is done by the map_dof_to_boundary_indices() function (described
- * above).
+ * In some cases, one wants to only work with DoFs that sit on the boundary. One
+ * application is, for example, if rather than interpolating non- homogeneous
+ * boundary values, one would like to project them. For this, we need two
+ * things: a way to identify nodes that are located on (parts of) the boundary,
+ * and a way to build matrices out of only degrees of freedom that are on the
+ * boundary (i.e. much smaller matrices, in which we do not even build the large
+ * zero block that stems from the fact that most degrees of freedom have no
+ * support on the boundary of the domain). The first of these tasks is done by
+ * the map_dof_to_boundary_indices() function (described above).
  *
  * The second part requires us first to build a sparsity pattern for the
  * couplings between boundary nodes, and then to actually build the components
- * of this matrix. While actually computing the entries of these small
- * boundary matrices is discussed in the MatrixCreator namespace, the creation
- * of the sparsity pattern is done by the create_boundary_sparsity_pattern()
- * function. For its work, it needs to have a numbering of all those degrees
- * of freedom that are on those parts of the boundary that we are interested
- * in. You can get this from the map_dof_to_boundary_indices() function. It
- * then builds the sparsity pattern corresponding to integrals like
- * $\int_\Gamma \varphi_{b2d(i)} \varphi_{b2d(j)} dx$, where $i$ and $j$ are
- * indices into the matrix, and $b2d(i)$ is the global DoF number of a degree
- * of freedom sitting on a boundary (i.e., $b2d$ is the inverse of the mapping
- * returned by map_dof_to_boundary_indices() function).
+ * of this matrix. While actually computing the entries of these small boundary
+ * matrices is discussed in the MatrixCreator namespace, the creation of the
+ * sparsity pattern is done by the create_boundary_sparsity_pattern() function.
+ * For its work, it needs to have a numbering of all those degrees of freedom
+ * that are on those parts of the boundary that we are interested in. You can
+ * get this from the map_dof_to_boundary_indices() function. It then builds the
+ * sparsity pattern corresponding to integrals like $\int_\Gamma
+ * \varphi_{b2d(i)} \varphi_{b2d(j)} dx$, where $i$ and $j$ are indices into the
+ * matrix, and $b2d(i)$ is the global DoF number of a degree of freedom sitting
+ * on a boundary (i.e., $b2d$ is the inverse of the mapping returned by
+ * map_dof_to_boundary_indices() function).
  *
+ * <h3>DoF coupling between surface triangulations and bulk triangulations</h3>
+ *
+ * When working with Triangulation and DoFHandler objects of different
+ * co-dimension, such as a `Triangulation<2,3>`, describing (part of) the
+ * boundary of a `Triangulation<3>`, and their corresponding DoFHandler objects,
+ * one often needs to build a one-to-one matching between the degrees of freedom
+ * that live on the surface Triangulation and those that live on the boundary of
+ * the bulk Triangulation. The GridGenerator::extract_boundary_mesh() function
+ * returns a mapping of surface cell iterators to face iterators, that can be
+ * used by the function map_boundary_to_bulk_dof_iterators() to construct a map
+ * between cell iterators of the surface DoFHandler, and the corresponding pair
+ * of cell iterator and face index of the bulk DoFHandler. Such map can be used
+ * to initialize FEValues and FEFaceValues for the corresponding DoFHandler
+ * objects. Notice that one must still ensure that the ordering of the
+ * quadrature points coincide in the two objects, in order to build a coupling
+ * matrix between the two sytesm.
  *
  * @ingroup dofs
  */
@@ -1411,6 +1424,105 @@ namespace DoFTools
                          std::vector<std::vector<bool>> & constant_modes);
   /** @} */
 
+  /**
+   * @name Coupling between DoFHandler objects on different dimensions
+   * @{
+   */
+
+  /**
+   * This function generates a mapping of codimension-1 active DoFHandler cell
+   * iterators to codimension-0 cells and face indices, for DoFHandler objects
+   * built on top of the boundary of a given `Triangulation<spacedim>`.
+   *
+   * If you need to couple a PDE defined on the surface of an existing
+   * Triangulation (as in step-38) with the PDE defined on the bulk (say, for
+   * example, a hyper ball and its boundary), the information that is returned
+   * by the function GridGenerator::extract_boundary_mesh() is not enough, since
+   * you need to build FEValues objects on active cell iterators of the
+   * DoFHandler defined on the surface mesh, and FEFaceValues objects on face
+   * iterators of the DoFHandler defined on the bulk mesh. This second step
+   * requires knowledge of the bulk cell iterator and of the corresponding face
+   * index.
+   *
+   * This function examines the map `c1_to_c0` returned by
+   * GridGenerator::extract_boundary_mesh() (when used with two Triangulation
+   * objects as input), and associates to each active cell iterator of the
+   * `c1_dh` DoFHandler that is also contained in the map `c1_to_c0`, a pair of
+   * corresponding active bulk cell iterator of `c0_dh`, and face index
+   * corresponding to the cell iterator on the surface mesh.
+   *
+   * An example usage of this function is the following:
+   *
+   * @code
+   * Triangulation<dim>          triangulation;
+   * Triangulation<dim - 1, dim> surface_triangulation;
+   *
+   * FE_Q<dim>       fe(1);
+   * DoFHandler<dim> dof_handler(triangulation);
+   *
+   * FE_Q<dim - 1, dim>       surface_fe(1);
+   * DoFHandler<dim - 1, dim> surface_dof_handler(surface_triangulation);
+   *
+   * GridGenerator::half_hyper_ball(triangulation);
+   * triangulation.refine_global(4);
+   *
+   * surface_triangulation.set_manifold(0, SphericalManifold<dim - 1, dim>());
+   * const auto surface_to_bulk_map =
+   *   GridGenerator::extract_boundary_mesh(triangulation,
+   *                                        surface_triangulation,
+   *                                        {0});
+   *
+   * dof_handler.distribute_dofs(fe);
+   * surface_dof_handler.distribute_dofs(surface_fe);
+   *
+   * // Extract the mapping between surface and bulk degrees of freedom:
+   * const auto surface_to_bulk_dof_iterator_map =
+   *   DoFTools::map_boundary_to_bulk_dof_iterators(surface_to_bulk_map,
+   *                                                dof_handler,
+   *                                                surface_dof_handler);
+   *
+   * // Loop over the map, and print some information:
+   * for (const auto &p : surface_to_bulk_dof_iterator_map)
+   *   {
+   *     const auto &surface_cell = p.first;
+   *     const auto &bulk_cell    = p.second.first;
+   *     const auto &bulk_face    = p.second.second;
+   *     deallog << "Surface cell " << surface_cell << " coincides with face "
+   *             << bulk_face << " of bulk cell " << bulk_cell << std::endl;
+   *   }
+   * @endcode
+   *
+   * \tparam dim The dimension of the codimension-0 mesh.
+   *
+   * \tparam spacedim The dimension of the underlying space.
+   *
+   * \param[in] c1_to_c0 A map from codimension-1 triangulation cell iterators
+   * to codimension-0 face iterators, as generated by the
+   * GridGenerators::extract_boundary_mesh() function.
+   *
+   * \param[in] c0_dh The DoFHandler object of the codimension-0 mesh.
+   *
+   * \param[in] c1_dh The DoFHandler object of the codimension-1 mesh.
+   *
+   * \return A std::map object that maps codimension-1 active DoFHandler cell
+   * iterators to a pair consisting of the corresponding codimension-0 cell
+   * iterator and face index.
+   */
+  template <int dim, int spacedim>
+  std::map<typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
+           std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
+                     unsigned int>>
+  map_boundary_to_bulk_dof_iterators(
+    const std::map<typename Triangulation<dim - 1, spacedim>::cell_iterator,
+                   typename Triangulation<dim, spacedim>::face_iterator>
+      &                                  c1_to_c0,
+    const DoFHandler<dim, spacedim> &    c0_dh,
+    const DoFHandler<dim - 1, spacedim> &c1_dh);
+
+  /**
+   * @}
+   */
+
   /**
    * @name Parallelization and domain decomposition
    * @{
index bd2f44dae1b286776a423d7147a62c4d2048199d..cf2d1b5280ccb0a540048f647dbea66317095d41 100644 (file)
@@ -1353,6 +1353,56 @@ namespace DoFTools
 
 
 
+  template <int dim, int spacedim>
+  std::map<typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
+           std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
+                     unsigned int>>
+  map_boundary_to_bulk_dof_iterators(
+    const std::map<typename Triangulation<dim - 1, spacedim>::cell_iterator,
+                   typename Triangulation<dim, spacedim>::face_iterator>
+      &                                  c1_to_c0,
+    const DoFHandler<dim, spacedim> &    c0_dh,
+    const DoFHandler<dim - 1, spacedim> &c1_dh)
+  {
+    // This is the returned object: a map of codimension-1 active dof cell
+    // iterators to codimension-0 cells and face indices
+    std::map<typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
+             std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
+                       unsigned int>>
+      c1_to_c0_cells_and_faces;
+
+    // Shortcut if there are no faces to check
+    if (c1_to_c0.empty())
+      return c1_to_c0_cells_and_faces;
+
+    // This is the partial inverse of the map passed as input, for dh
+    std::map<typename Triangulation<dim, spacedim>::face_iterator,
+             typename DoFHandler<dim - 1, spacedim>::active_cell_iterator>
+      c0_to_c1;
+
+    // map volume mesh face -> codimension 1 dof cell
+    for (auto c1_cell : c1_dh.active_cell_iterators())
+      if (c1_to_c0.find(c1_cell) != c1_to_c0.end())
+        c0_to_c1[c1_to_c0.at(c1_cell)] = c1_cell;
+
+    // generate a mapping that maps codimension-1 cells
+    // to codimension-0 cells and faces
+    for (auto cell :
+         c0_dh.active_cell_iterators()) // disp_dof.active_cell_iterators())
+      for (const auto f : cell->face_indices())
+        if (cell->face(f)->at_boundary() &&
+            c0_to_c1.find(cell->face(f)) != c0_to_c1.end())
+          {
+            const auto &c1_cell               = c0_to_c1[cell->face(f)];
+            c1_to_c0_cells_and_faces[c1_cell] = {cell, f};
+          }
+    // Check the dimensions: make sure all active cells we had have been mapped.
+    AssertDimension(c1_to_c0_cells_and_faces.size(), c0_to_c1.size());
+    return c1_to_c0_cells_and_faces;
+  }
+
+
+
   template <int dim, int spacedim>
   void
   get_active_fe_indices(const DoFHandler<dim, spacedim> &dof_handler,
index 9b9efa853a6cdb8371ab7cc9772563110c9f2cc6..456ea317c8ffa2dbb081a2b11e735cf1f9d09ec7 100644 (file)
@@ -484,3 +484,28 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
     \}
 #endif
   }
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
+  {
+#if deal_II_dimension >= 2 && (deal_II_dimension <= deal_II_space_dimension)
+    namespace DoFTools
+    \{
+      template std::map<
+        typename DoFHandler<deal_II_dimension - 1,
+                            deal_II_space_dimension>::active_cell_iterator,
+        std::pair<
+          typename DoFHandler<deal_II_dimension,
+                              deal_II_space_dimension>::active_cell_iterator,
+          unsigned int>>
+      map_boundary_to_bulk_dof_iterators<deal_II_dimension,
+                                         deal_II_space_dimension>(
+        const std::map<
+          typename Triangulation<deal_II_dimension - 1,
+                                 deal_II_space_dimension>::cell_iterator,
+          typename Triangulation<deal_II_dimension,
+                                 deal_II_space_dimension>::face_iterator> &,
+        const DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+        const DoFHandler<deal_II_dimension - 1, deal_II_space_dimension> &);
+    \}
+#endif
+  }
\ No newline at end of file
diff --git a/tests/dofs/map_boundary_to_bulk_dof_iterators_01.cc b/tests/dofs/map_boundary_to_bulk_dof_iterators_01.cc
new file mode 100644 (file)
index 0000000..36e8ecb
--- /dev/null
@@ -0,0 +1,90 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2023 by the deal.II authors
+//
+//    This file is subject to LGPL and may not be distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------
+
+// Test map_boundary_to_bulk_dof_iterators
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test()
+{
+  Triangulation<dim>          triangulation;
+  Triangulation<dim - 1, dim> surface_triangulation;
+
+  FE_Q<dim>       fe(1);
+  DoFHandler<dim> dof_handler(triangulation);
+
+  FE_Q<dim - 1, dim>       surface_fe(1);
+  DoFHandler<dim - 1, dim> surface_dof_handler(surface_triangulation);
+
+  GridGenerator::half_hyper_ball(triangulation);
+  triangulation.refine_global(4 - dim);
+
+  surface_triangulation.set_manifold(0, SphericalManifold<dim - 1, dim>());
+  const auto surface_to_bulk_map =
+    GridGenerator::extract_boundary_mesh(triangulation,
+                                         surface_triangulation,
+                                         {0});
+
+  deallog << "Bulk mesh active cells:" << triangulation.n_active_cells()
+          << std::endl
+          << "Surface mesh active cells:"
+          << surface_triangulation.n_active_cells() << std::endl;
+
+  dof_handler.distribute_dofs(fe);
+  surface_dof_handler.distribute_dofs(surface_fe);
+
+  // Log degrees of freedom:
+  deallog << "Bulk mesh degrees of freedom:" << dof_handler.n_dofs()
+          << std::endl
+          << "Surface mesh degrees of freedom:" << surface_dof_handler.n_dofs()
+          << std::endl;
+
+  // Extract the mapping between surface and bulk degrees of freedom:
+  const auto surface_to_bulk_dof_iterator_map =
+    DoFTools::map_boundary_to_bulk_dof_iterators(surface_to_bulk_map,
+                                                 dof_handler,
+                                                 surface_dof_handler);
+
+  // Loop over the map, and print some information:
+  for (const auto &p : surface_to_bulk_dof_iterator_map)
+    {
+      const auto &surface_cell = p.first;
+      const auto &bulk_cell    = p.second.first;
+      const auto &bulk_face    = p.second.second;
+      deallog << "Surface cell " << surface_cell << " coincides with face "
+              << bulk_face << " of bulk cell " << bulk_cell << std::endl;
+    }
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  test<2>();
+  test<3>();
+}
diff --git a/tests/dofs/map_boundary_to_bulk_dof_iterators_01.output b/tests/dofs/map_boundary_to_bulk_dof_iterators_01.output
new file mode 100644 (file)
index 0000000..cf4ceb0
--- /dev/null
@@ -0,0 +1,41 @@
+
+DEAL::Bulk mesh active cells:64
+DEAL::Surface mesh active cells:12
+DEAL::Bulk mesh degrees of freedom:77
+DEAL::Surface mesh degrees of freedom:13
+DEAL::Surface cell 2.0 coincides with face 2 of bulk cell 2.5
+DEAL::Surface cell 2.1 coincides with face 2 of bulk cell 2.4
+DEAL::Surface cell 2.2 coincides with face 2 of bulk cell 2.1
+DEAL::Surface cell 2.3 coincides with face 2 of bulk cell 2.0
+DEAL::Surface cell 2.4 coincides with face 2 of bulk cell 2.37
+DEAL::Surface cell 2.5 coincides with face 2 of bulk cell 2.36
+DEAL::Surface cell 2.6 coincides with face 2 of bulk cell 2.33
+DEAL::Surface cell 2.7 coincides with face 2 of bulk cell 2.32
+DEAL::Surface cell 2.8 coincides with face 0 of bulk cell 2.48
+DEAL::Surface cell 2.9 coincides with face 0 of bulk cell 2.50
+DEAL::Surface cell 2.10 coincides with face 0 of bulk cell 2.56
+DEAL::Surface cell 2.11 coincides with face 0 of bulk cell 2.58
+DEAL::Bulk mesh active cells:48
+DEAL::Surface mesh active cells:20
+DEAL::Bulk mesh degrees of freedom:77
+DEAL::Surface mesh degrees of freedom:25
+DEAL::Surface cell 1.0 coincides with face 4 of bulk cell 1.0
+DEAL::Surface cell 1.1 coincides with face 4 of bulk cell 1.2
+DEAL::Surface cell 1.2 coincides with face 4 of bulk cell 1.1
+DEAL::Surface cell 1.3 coincides with face 4 of bulk cell 1.3
+DEAL::Surface cell 1.4 coincides with face 0 of bulk cell 1.8
+DEAL::Surface cell 1.5 coincides with face 0 of bulk cell 1.12
+DEAL::Surface cell 1.6 coincides with face 0 of bulk cell 1.10
+DEAL::Surface cell 1.7 coincides with face 0 of bulk cell 1.14
+DEAL::Surface cell 1.8 coincides with face 4 of bulk cell 1.24
+DEAL::Surface cell 1.9 coincides with face 4 of bulk cell 1.26
+DEAL::Surface cell 1.10 coincides with face 4 of bulk cell 1.25
+DEAL::Surface cell 1.11 coincides with face 4 of bulk cell 1.27
+DEAL::Surface cell 1.12 coincides with face 0 of bulk cell 1.32
+DEAL::Surface cell 1.13 coincides with face 0 of bulk cell 1.36
+DEAL::Surface cell 1.14 coincides with face 0 of bulk cell 1.34
+DEAL::Surface cell 1.15 coincides with face 0 of bulk cell 1.38
+DEAL::Surface cell 1.16 coincides with face 0 of bulk cell 1.40
+DEAL::Surface cell 1.17 coincides with face 0 of bulk cell 1.44
+DEAL::Surface cell 1.18 coincides with face 0 of bulk cell 1.42
+DEAL::Surface cell 1.19 coincides with face 0 of bulk cell 1.46

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.