From 0963f163d7c2c09bdd2366c8d7a85d8e4601e0b0 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 12 Apr 2023 11:15:38 +0300 Subject: [PATCH] Map boundary to bulk dof iterators. --- doc/news/changes/minor/20230412LucaHeltai | 7 + include/deal.II/dofs/dof_tools.h | 290 ++++++++++++------ source/dofs/dof_tools.cc | 50 +++ source/dofs/dof_tools.inst.in | 25 ++ .../map_boundary_to_bulk_dof_iterators_01.cc | 90 ++++++ ...p_boundary_to_bulk_dof_iterators_01.output | 41 +++ 6 files changed, 414 insertions(+), 89 deletions(-) create mode 100644 doc/news/changes/minor/20230412LucaHeltai create mode 100644 tests/dofs/map_boundary_to_bulk_dof_iterators_01.cc create mode 100644 tests/dofs/map_boundary_to_bulk_dof_iterators_01.output diff --git a/doc/news/changes/minor/20230412LucaHeltai b/doc/news/changes/minor/20230412LucaHeltai new file mode 100644 index 0000000000..56331779d9 --- /dev/null +++ b/doc/news/changes/minor/20230412LucaHeltai @@ -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() +
+(Luca Heltai, 2023/04/12) diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index fb0397d8db..22d954c0e6 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -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. * *

Setting up sparsity patterns

* - * 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. * * *

DoF numberings on boundaries

* - * 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.) * * *

Setting up sparsity patterns for boundary matrices

* - * 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). * + *

DoF coupling between surface triangulations and bulk triangulations

+ * + * 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> & 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`. + * + * 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 triangulation; + * Triangulation surface_triangulation; + * + * FE_Q fe(1); + * DoFHandler dof_handler(triangulation); + * + * FE_Q surface_fe(1); + * DoFHandler surface_dof_handler(surface_triangulation); + * + * GridGenerator::half_hyper_ball(triangulation); + * triangulation.refine_global(4); + * + * surface_triangulation.set_manifold(0, SphericalManifold()); + * 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 + std::map::active_cell_iterator, + std::pair::active_cell_iterator, + unsigned int>> + map_boundary_to_bulk_dof_iterators( + const std::map::cell_iterator, + typename Triangulation::face_iterator> + & c1_to_c0, + const DoFHandler & c0_dh, + const DoFHandler &c1_dh); + + /** + * @} + */ + /** * @name Parallelization and domain decomposition * @{ diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc index bd2f44dae1..cf2d1b5280 100644 --- a/source/dofs/dof_tools.cc +++ b/source/dofs/dof_tools.cc @@ -1353,6 +1353,56 @@ namespace DoFTools + template + std::map::active_cell_iterator, + std::pair::active_cell_iterator, + unsigned int>> + map_boundary_to_bulk_dof_iterators( + const std::map::cell_iterator, + typename Triangulation::face_iterator> + & c1_to_c0, + const DoFHandler & c0_dh, + const DoFHandler &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::active_cell_iterator, + std::pair::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::face_iterator, + typename DoFHandler::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 void get_active_fe_indices(const DoFHandler &dof_handler, diff --git a/source/dofs/dof_tools.inst.in b/source/dofs/dof_tools.inst.in index 9b9efa853a..456ea317c8 100644 --- a/source/dofs/dof_tools.inst.in +++ b/source/dofs/dof_tools.inst.in @@ -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::active_cell_iterator, + std::pair< + typename DoFHandler::active_cell_iterator, + unsigned int>> + map_boundary_to_bulk_dof_iterators( + const std::map< + typename Triangulation::cell_iterator, + typename Triangulation::face_iterator> &, + const DoFHandler &, + const DoFHandler &); + \} +#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 index 0000000000..36e8ecb8e9 --- /dev/null +++ b/tests/dofs/map_boundary_to_bulk_dof_iterators_01.cc @@ -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 +#include +#include + +#include + +#include +#include +#include + +#include "../tests.h" + + + +template +void +test() +{ + Triangulation triangulation; + Triangulation surface_triangulation; + + FE_Q fe(1); + DoFHandler dof_handler(triangulation); + + FE_Q surface_fe(1); + DoFHandler surface_dof_handler(surface_triangulation); + + GridGenerator::half_hyper_ball(triangulation); + triangulation.refine_global(4 - dim); + + surface_triangulation.set_manifold(0, SphericalManifold()); + 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 index 0000000000..cf4ceb0bcf --- /dev/null +++ b/tests/dofs/map_boundary_to_bulk_dof_iterators_01.output @@ -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 -- 2.39.5