 |
deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
|
Loading...
Searching...
No Matches
|
deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
|
This program was contributed by Marco Feder <[email protected]>.
It comes without any warranty or support by its authors or the authors of deal.II.
This program is part of the deal.II code gallery and consists of the following files (click to inspect):
 |  |  |  |
 |  |  |  |
 |  |  |
This program solves a Poisson problem on an agglomerated polytopal mesh using a symmetric interior penalty discontinuous Galerkin (SIPDG) method. Agglomerates are constructed by an R-tree based spatial indexing strategy, following the approach proposed in [1]. In addition, a graph-based METIS partitioner is also provided in this program for comparison.
As in the tutorial programs, type
cmake -DDEAL_II_DIR=/path/to/deal.II .
on the command line to configure the program. After that, you can compile with make and run with either make run or with
./agglomeration_poisson
on the command line.
Running the program produces two kinds of output:
.vtu).The program prints a short summary including:
FE degree);Size of tria);N subdomains);#DoFs, L2 error, and H1 error.The program writes the following .vtu files for each run:
grid_input_mesh.vtu, containing the input mesh from gmsh;grid_fine_mesh_refined.vtu, containing the globally refined fine mesh;grid_<partitioner>_<n_subdomains>.vtu, containing the agglomerated mesh partition information (cell-wise agglomeration labels);interpolated_solution_<partitioner>_<n_subdomains>.vtu, containing the numerical solution interpolated to the fine grid together with agglomerate labels for visualization.These files can be visualized in ParaView to inspect both the agglomeration structure and the computed solution.
We consider the Poisson problem in a bounded, simply connected domain \(\Omega \subset \mathbb{R}^d\), \(d = 2,3\). The strong formulation reads
\begin{align*} -\Delta u &= f && \text{in } \Omega, \\ u &= u_D && \text{on } \partial\Omega, \end{align*}
where the right-hand side satisfies \(f \in L^2(\Omega)\) and the prescribed Dirichlet data satisfy \(u_D \in H^{1/2}(\partial\Omega)\).
The corresponding weak formulation is: find \(u \in H^1(\Omega)\) with \(u = u_D\) on \(\partial\Omega\) such that
\begin{align*} \int_{\Omega} \nabla u \cdot \nabla v \,\mathrm d\mathbf{x} =\int_{\Omega} f\, v \,\mathrm d\mathbf{x} \qquad \text{for all } v \in H_0^1(\Omega). \end{align*}
We discretize the weak formulation by a symmetric interior penalty discontinuous Galerkin (SIPDG) method on the agglomerated polytopal mesh \(T_h\), whose elements \(K \in T_h\) are mutually disjoint open polygons (for \(d=2\)) or polyhedra (for \(d=3\)). The mesh skeleton is defined by
\begin{align*} \Gamma := \bigcup_{K \in T_h} \partial K. \end{align*}
The mesh skeleton \(\Gamma\) is decomposed into \((d-1)\)–dimensional simplices \(F\) representing the mesh faces, shared by at most two elements. These are distinct from elemental interfaces, which are defined as the simply connected components of the intersection between the boundary of an element and either a neighboring element or \(\partial \Omega\). As such, an interface between two elements may consist of more than one face, separated by hanging nodes/edges shared by those two elements only. We denote by \(\Gamma_{\mathrm{int}}\) the union of all interior faces, and by \(\Gamma_{\mathrm D} := \Gamma \cap \partial\Omega\) the union of Dirichlet boundary faces.
In practice, the polytopic mesh is obtained by agglomeration, so that each element \(K \in T_h\) is the union of a collection of leaf cells. For each agglomerated element \(K\), we associate an axis-aligned bounding box \(B_K\). On \(B_K\) we define the standard polynomial space \(Q_p(B_K)\) spanned by tensor-product Lagrange polynomials of degree \(p\) in each coordinate direction. Since \(K \subset B_K\), the basis on \(K\) is defined by restricting each basis function to \(K\). In the implementation, this corresponds to using the deal.II finite element FE_DGQ on the bounding box and taking its restriction to the agglomerated element. The global discrete space \(V_h\) is then obtained by assembling these local spaces in a discontinuous manner over all \(K \in T_h\). For \(u_h, v_h \in V_h\) we use the broken gradient \(\nabla_h\) and the standard jump and average operators \([\![\cdot]\!]\) and \(\{\!\!\{\cdot\}\!\!\}\) on faces.
The DG formulation reads: find \(u_h \in V_h\) such that
\begin{align*} B(u_h,v_h) = l(v_h) \qquad \forall\, v_h \in V_h, \end{align*}
with
\begin{align*} B(u_h,v_h) &= \int_{\Omega} \nabla_h u_h \cdot \nabla_h v_h \,\mathrm d\mathbf{x} \\ &\quad - \int_{\Gamma} \Bigl( \{\!\!\{\nabla u_h\}\!\!\} \cdot [\![v_h]\!] + \{\!\!\{\nabla v_h\}\!\!\} \cdot [\![u_h]\!] \Bigr)\,\mathrm d s \\ &\quad + \int_{\Gamma} \sigma \,[\![u_h]\!] \cdot [\![v_h]\!] \,\mathrm d s, \end{align*}
and
\begin{align*} l(v_h) = \int_\Omega f\, v_h \,\mathrm d\mathbf{x} + \int_{\Gamma_{\mathrm D}} u_D \bigl(\sigma v_h - \nabla v_h \cdot \mathbf n\bigr)\,\mathrm d s. \end{align*}
The penalty parameter is chosen as
\begin{align*} \sigma(\mathbf x) = C_\sigma \begin{cases} \dfrac{(p+1)(p+d)}{ h_{B_K} }, & \text{if } \mathbf x \in \partial K \cap \partial\Omega, \\[0.5em] \dfrac{(p+1)(p+d)}{\min\{h_{B_K}^+,h_{B_K}^-\}}, & \text{if } \mathbf x \in \Gamma_{\mathrm{int}}, \end{cases} \end{align*}
where \(h_{B_K}^\pm\) denote the diameters of the bounding boxes associated with the two elements sharing the interior face. In this program, we set \(C_\sigma = 10\).
This scheme is well posed and admits optimal-order a priori error estimates. More precisely, assuming that \(u|_K \in H^{s+1}(K)\) for all \(K \in T_h\) and some \(1 \le s \le p\), there exists a constant \(C > 0\), independent of \(h\), such that
\begin{align*} \|u - u_h\|_{L^2(\Omega)} \le C\, h^{s+1} \, |u|_{H^{s+1}(\Omega)}, \end{align*}
and
\begin{align*} \|\nabla(u - u_h)\|_{L^2(\Omega)} \le C\, h^{s} \, |u|_{H^{s+1}(\Omega)}. \end{align*}
We refer to [2] for details of the analysis.
Agglomeration is a natural mechanism for constructing polytopic meshes. This program supports two strategies for generating agglomerates, corresponding to the choices metis and rtree. In this example, we mainly focus on the rtree strategy, which is the method developed and introduced in the work [1].
In the rtree option, axis-aligned bounding boxes of all fine cells are inserted into a spatial R-tree. Agglomerates are obtained by grouping the cells whose bounding boxes belong to the same node at a user-selected extraction level of the tree. This purely geometric strategy does not require external graph partitioners and is typically fast and scalable. The number and shape of the agglomerates are determined by the R-tree structure and the chosen level.
This approach is particularly suitable for multilevel methods, where a nested agglomerated hierarchy is desirable.
At the data-structure level, we distinguish leaf nodes and internal nodes:
As a result, each internal node represents a spatial grouping of the objects below it.
The R-tree is used here as a geometry-aware structure for organizing cell bounding boxes into hierarchical groups. Our construction is guided by the classical R*-tree criteria of Beckmann et al. [3], namely:
These criteria improve the spatial quality of the hierarchy and typically lead to better grouping and query behavior.
To make the R-tree construction more intuitive, we first illustrate the relation between geometric objects, their minimum bounding rectangles (MBRs), and the corresponding R-tree hierarchy. The left image shows geometric objects together with their enclosing MBRs, while the right image shows the associated tree structure (leaf and internal nodes). This visual example helps explain how the hierarchical grouping is later used to extract agglomerates.
The construction of an R-tree spatial index on an arbitrary fine grid provides a natural agglomeration strategy with the following features:
These properties make the R-tree approach an attractive alternative to graph-based agglomeration methods.
Given a collection of fine-level cells, or of cut-cell bounding boxes, the agglomeration procedure is as follows:
This yields a nested hierarchy with natural parent-child relations across levels and, by repeating the extraction at different levels, produces a sequence of nested agglomerated meshes that can be used in multilevel solvers and preconditioners.
For the agglomeration workflow considered here, the R-tree-based extraction has the following practical features:
Boost.Geometry R-tree.In the present setting, this makes the R-tree approach particularly convenient for constructing nested agglomerated meshes in multilevel finite element and DG settings.
For illustration, the following images show R-tree-based agglomeration on a structured fine mesh:
Figures (3)-(5) show, respectively, the original fine mesh, the blocks induced by the R-tree on the cell bounding boxes, and the corresponding tree structure.
In the metis option, the adjacency graph of the fine mesh is constructed with one vertex per cell and edges between face-neighboring cells. This graph is then partitioned by the multilevel graph partitioner METIS into a prescribed number of parts, and each part defines one agglomerate; see [4] for details.
We consider the Poisson problem on the unit square \(\Omega = (0,1)^2\) with the manufactured exact solution
\begin{align*} u(x,y) = \sin(\pi x)\sin(\pi y). \end{align*}
The corresponding right-hand side is
\begin{align*} f(x,y) = 2\pi^2 \sin(\pi x)\sin(\pi y). \end{align*}
This manufactured solution allows us to compute the global \(L^2\)- and \(H^1\)-seminorm errors of the discrete solution in order to assess the quality of the numerical approximation.
In this example, we start from an unstructured initial mesh and then perform five global refinement steps. The resulting meshes are shown below.
Agglomerates are then constructed from the fine level mesh by METIS and by the R-tree strategy, leading to different polytopal meshes. The following images compare the resulting agglomerates corresponding to two different agglomeration levels (91 and 364 agglomerates).
These plots illustrate how the two strategies distribute and shape the agglomerates starting from the same fine level mesh. In particular, the R-tree approach produces geometry-driven groupings induced by the spatial hierarchy, while METIS produces graph-based partitions of the cell adjacency graph.
To assess the discretization accuracy, we next compare the error behavior under mesh refinement or increasing the polynomial degree. The plots below are obtained by collecting the program outputs over multiple runs and post-processing the reported error data.
(12) h-convergence for Q_p elements (p=1,2,3) with METIS and R-tree agglomeration*
The figure reports the \(L^2\)-norm error with respect to the manufactured solution. Optimal convergence rates are observed for all polynomial degrees and for both agglomeration strategies. In addition, the curves associated with the R-tree approach are marginally lower than or comparable to those obtained with METIS-based partitioning.
We next compare p-convergence for the two agglomeration strategies in both the \(L^2\)-norm and the \(H^1\)-seminorm.
(13) Convergence under p-refinement for METIS and R-tree agglomeration (p = 1,2,3,4,5)*
In addition to accuracy, the cost of constructing the agglomerated polytopal meshes is also relevant in practice. The following timing plot compares the wall-clock time required by the R-tree and METIS strategies. The timing values are collected from the program outputs and summarized in post-processing.
(14) Wall-clock time (seconds) for building polytopal grids with R-tree and METIS*
Finally, the program also outputs VTU files for visualization in ParaView. Figures (15) and (16) show the solution computed on 91 agglomerates and then interpolated back onto the fine mesh, for the R-tree and METIS strategies, respectively, visualized in ParaView using Warp By Scalar together with the Surface representation.
[1] Marco Feder, Andrea Cangiani and Luca Heltai (2025), R3MG: R-tree based agglomeration of polytopal grids with applications to multilevel methods. DOI: 10.1016/j.jcp.2025.113773
[2] Daniele Antonio Di Pietro and Alexandre Ern (2012), Mathematical Aspects of Discontinuous Galerkin Methods. DOI: 10.1007/978-3-642-22980-0
[3] Norbert Beckmann, Hans-Peter Kriegel, Ralf Schneider and Bernhard Seeger (1990), The R*-tree: an efficient and robust access method for points and rectangles. DOI: 10.1145/93597.98741
[4] George Karypis and Vipin Kumar (1998), A fast and high quality multilevel scheme for partitioning irregular graphs. DOI: 10.1137/S1064827595287997
The following headers provide the deal.II functionality needed in this example. Most of them are standard components for mesh handling, finite element mappings, linear algebra, and graphical output. In addition, we include the agglomeration-specific headers that define the data structures and utilities used to construct and manage polytopal agglomerates.
deal.II base utilities.
Finite element mappings.
Grid generation, mesh input/output, and mesh-related utilities.
Linear algebra objects and sparse direct solvers.
Output of finite element data for visualization.
Agglomeration-specific headers used in this example.
C++ standard library headers.
We use the struct ConvergenceInfo to store the number of degrees of freedom together with the corresponding L2 and H1 errors, and print a simple convergence table to the console.
We will compare the performance of three different partitioning strategies: using METIS, using an R-tree based agglomeration, or not performing any partitioning at all.
We then implement the manufactured right-hand side f(x, y) = 2 π² sin(π x) sin(π y), which corresponds to the exact solution u(x, y) = sin(π x) sin(π y).
Exact solution is set as u(x,y) = sin(pi x) sin(pi y). It is used to impose Dirichlet boundary conditions and to evaluate the L2 and H1-seminorm errors. Its gradient is also provided for the computation of the H1 error.
The Poisson<dim> class encapsulates the solution of the model Poisson problem
\[ -\Delta u = f \quad \text{in } \Omega, \qquad u = u_D \quad \text{on } \partial\Omega. \]
It sets up a fine triangulation, constructs agglomerated polytopal cells according to the chosen partitioning strategy, assembles the symmetric interior penalty DG discretization on the agglomerated mesh, solves the resulting linear system, and finally postprocesses the numerical solution by writing visualization output and computing global error norms.
The constructor initializes the Poisson<dim> solver with the selected partitioning strategy, agglomeration parameters, polynomial degree, and the manufactured exact solution and right-hand side.
Initialize manufactured solution.
Build the fine triangulation from a Gmsh mesh, apply a global refinement, initialize the cache and agglomeration handler, and define agglomerates according to the selected partitioning strategy.
To finalize the agglomeration. In the no-partition case, each fine cell is declared as its own agglomerate. The function then distributes the degrees of freedom on the agglomerated mesh, builds the corresponding sparsity pattern, and writes a VTU file visualizing the agglomeration and the partitioning of the fine grid.
Assemble the global SIPG matrix and right-hand side on the agglomerated mesh.
It initializes the system matrix and right-hand side, sets up FEValues objects on polytopal cells and interfaces, and then adds the volume, boundary, and interior face contributions of the symmetric interior penalty formulation.
Next, we define the four dofsxdofs matrices needed to assemble jumps and averages.
Distribute the local contributions to the global system.
Solve the linear system by means of a sparse direct solver.
Write VTU output and compute the global \(L^2\) and \(H^1\)-seminorm errors of the agglomerated DG approximation.
Mark fine cells belonging to the same agglomerate.
Return the number of degrees of freedom on the agglomerated mesh.
Return the pair consisting of the \(L^2\) error and the \(H^1\)-seminorm error of the numerical solution.
Run the full workflow: mesh generation, agglomeration setup, assembly, solution, and postprocessing.
Driver code.
Forward declarations
The following path is needed when the present function is called from neighbor_of_neighbor()
Use the id of the master cell to uniquely identify the neighboring agglomerate
Get master_id from the neighboring ghost polytope. This uniquely identifies the neighboring polytope among all processors.
Use the id of the master cell to uniquely identify the neighboring agglomerate
First, make sure it's not a boundary face.
if it is locally owned, retrieve the number of faces
The neighboring polytope is not locally owned. We need to get the number of its faces from the neighboring rank.
First, retrieve the CellId of the neighboring polytope.
Then, get the neighboring rank
From the neighboring rank, use the CellId of the neighboring polytope to get the number of its faces.
Loop over all faces of neighboring agglomerate
Check if same CellId
Face is at boundary
---------------------------— inline functions ----------------------—
neigh_cell is ghosted
neigh_cell is ghosted, use the CellId of that agglomerate
Forward the call to the master cell
Get the bounding box associated with the master cell
Standard deal.II way to get the measure of a cell.
Increment the present index and update the polytope
Make sure not to query the CellId if it's past the last
Decrement the present index and update the polytope
Forward the call to the master cell using the right DoFHandler.
Forward declarations
clear all the members
disconnect the signal
TODO: move it to private interface
First disconnect existing connections
////////////////////////////////////////////////////
n_faces
CellId (including slaves)
send to neighborign rank the information that
CellIds from neighboring rank
Exchange neighboring bounding boxes
Exchange DoF indices with ghosted polytopes
Exchange qpoints
Exchange jxws
Exchange normals
Exchange values
////////////////////////////////////////////////////
The FiniteElement space we have on each cell. Currently supported types are FE_DGQ and FE_DGP elements.
Associate the master cell to the slaves.
Map the master cell index with the polytope index
Dummy FiniteElement objects needed only to generate quadratures
Support for hp::FECollection
Stores quadrature rules; these QCollections should have the same size as hp_fe_collection
containing only dummy_fe Note: The above two variables provide an hp::FECollection interface but actually contain only one element each.
Analogous to no_values and no_face_values, but used when different cells employ different FEs or quadratures
---------------------------— inline functions ----------------------—
if different subdomain, then by construction they will not be together if (cell->subdomain_id() != other_cell->subdomain_id()) return false; else
---------------------------— inline functions ----------------------—
node
Done with node number 'node_counter' on level target_level.
I am on a child (internal) node on a deeper level.
Keep visiting until you go to the leafs.
looping through entries of node
done with node on level l > target_level (not just "target_level+1). @code n_nodes_per_level[level_backup]++; const types::global_cell_index node_idx = n_nodes_per_level[level_backup] - 1; // so to start from 0 parent_node_to_children_nodes[{n_nodes_per_level[level_backup - 1], level_backup - 1}] .push_back(node_idx); level = level_backup; } } template <typename Value, typename Options, typename Translator, typename Box, typename Allocators> void Rtree_visitor<Value, Options, Translator, Box, Allocators>::operator()( const Rtree_visitor::Leaf &leaf) { using elements_type = typename boost::geometry::index::detail::rtree::elements_type< Leaf>::type; // pairs of bounding box and pointer to child node const elements_type &elements = boost::geometry::index::detail::rtree::elements(leaf); if (level == target_level) { @endcode If I want to extract from leaf node, i.e. the target_level is the last one where leafs are grouped together. @code const auto offset = agglomerates.size(); agglomerates.resize(offset + 1); for (const auto &it : elements) agglomerates[node_counter].push_back(it.second); ++node_counter; n_nodes_per_level[target_level]++; } else { for (const auto &it : elements) agglomerates[node_counter].push_back(it.second); if (level == target_level + 1) { const unsigned int node_idx = n_nodes_per_level[level]; parent_node_to_children_nodes[{n_nodes_per_level[level - 1], level - 1}] .push_back(node_idx); n_nodes_per_level[level]++; } } } } // namespace internal /** * Helper class which handles agglomeration based on the R-tree data * structure. Notice that the R-tree type is assumed to be an R-star-tree. */ template <int dim, typename RtreeType> class CellsAgglomerator { public: template <int, int> friend class ::AgglomerationHandler; /** * Constructor. It takes a given rtree and an integer representing the * index of the level to be extracted. */ CellsAgglomerator(const RtreeType &rtree, const unsigned int extraction_level); /** * Extract agglomerates based on the current tree and the extraction level. * This function returns a reference to */ const std::vector< std::vector<typename Triangulation<dim>::active_cell_iterator>> & extract_agglomerates(); /** * Get total number of levels. */ inline unsigned int get_n_levels() const; /** * Return the number of nodes present in level @p level. */ inline types::global_cell_index get_n_nodes_per_level(const unsigned int level) const; /** * This function returns a map which associates to each node on level * @p extraction_level a list of children. */ inline const std::map< std::pair<types::global_cell_index, types::global_cell_index>, std::vector<types::global_cell_index>> & get_hierarchy() const; private: /** * Raw pointer to the actual R-tree. */ RtreeType *rtree; /** * Extraction level. */ const unsigned int extraction_level; /** * Store agglomerates obtained after recursive extraction on nodes of * level @p extraction_level. */ std::vector<std::vector<typename Triangulation<dim>::active_cell_iterator>> agglomerates_on_level; /** * Vector storing the number of nodes (and, ultimately, agglomerates) for * each level. */ std::vector<types::global_cell_index> n_nodes_per_level; /** * Map which maps a node parent @n on level @p l to a vector of integers * which stores the index of children. */ std::map<std::pair<types::global_cell_index, types::global_cell_index>, std::vector<types::global_cell_index>> parent_node_to_children_nodes; }; template <int dim, typename RtreeType> CellsAgglomerator<dim, RtreeType>::CellsAgglomerator( const RtreeType &tree, const unsigned int extraction_level_) : extraction_level(extraction_level_) { rtree = const_cast<RtreeType *>(&tree); Assert(n_levels(*rtree), ExcMessage("At least two levels are needed.")); } template <int dim, typename RtreeType> const std::vector< std::vector<typename Triangulation<dim>::active_cell_iterator>> & CellsAgglomerator<dim, RtreeType>::extract_agglomerates() { AssertThrow(extraction_level <= n_levels(*rtree), ExcInternalError("You are trying to extract level " + std::to_string(extraction_level) + " of the tree, but it only has a total of " + std::to_string(n_levels(*rtree)) + " levels.")); using RtreeView = boost::geometry::index::detail::rtree::utilities::view<RtreeType>; RtreeView rtv(*rtree); n_nodes_per_level.resize(rtv.depth() + 1); // store how many nodes we have for each level. if (rtv.depth() == 0) { @endcode The below algorithm does not work for <tt>rtv.depth()==0</tt>, which might happen if the number entries in the tree is too small. @code agglomerates_on_level.resize(1); agglomerates_on_level[0].resize(1); } else { const unsigned int target_level = std::min<unsigned int>(extraction_level, rtv.depth()); internal::Rtree_visitor<typename RtreeView::value_type, typename RtreeView::options_type, typename RtreeView::translator_type, typename RtreeView::box_type, typename RtreeView::allocators_type> extractor_visitor(rtv.translator(), target_level, agglomerates_on_level, n_nodes_per_level, parent_node_to_children_nodes); rtv.apply_visitor(extractor_visitor); } return agglomerates_on_level; } @endcode ---------------------------— inline functions ----------------------— @code template <int dim, typename RtreeType> inline unsigned int CellsAgglomerator<dim, RtreeType>::get_n_levels() const { return n_levels(*rtree); } template <int dim, typename RtreeType> inline types::global_cell_index CellsAgglomerator<dim, RtreeType>::get_n_nodes_per_level( const unsigned int level) const { return n_nodes_per_level[level]; } template <int dim, typename RtreeType> inline const std::map< std::pair<types::global_cell_index, types::global_cell_index>, std::vector<types::global_cell_index>> & CellsAgglomerator<dim, RtreeType>::get_hierarchy() const { Assert(parent_node_to_children_nodes.size(), ExcMessage( "The hierarchy has not been computed. Did you forget to call" " extract_agglomerates() first?")); return parent_node_to_children_nodes; } } // namespace dealii #endif @endcode <a name="ann-include/mapping_box.h"></a> <h1>Annotated version of include/mapping_box.h</h1> @code /* ----------------------------------------------------------------------------- * * SPDX-License-Identifier: LGPL-2.1-or-later * Copyright (C) 2025 by Marco Feder, Pasquale Claudio Africa, Xinping Gui, * Andrea Cangiani * * This file is part of the deal.II code gallery. * * ----------------------------------------------------------------------------- */ #ifndef dealii_mapping_box_h #define dealii_mapping_box_h #include <deal.II/base/config.h> #include <deal.II/base/bounding_box.h> #include <deal.II/base/qprojector.h> #include <deal.II/fe/mapping.h> #include <cmath> DEAL_II_NAMESPACE_OPEN /** * @addtogroup mapping * @{ */ /** * A class providing a mapping from the reference cell to cells that are * axiparallel, i.e., that have the shape of rectangles (in 2d) or * boxes (in 3d) with edges parallel to the coordinate directions. The * class therefore provides functionality that is equivalent to what, * for example, MappingQ would provide for such cells. However, knowledge * of the shape of cells allows this class to be substantially more * efficient. * * Specifically, the mapping is meant for cells for which the mapping from * the reference to the real cell is a scaling along the coordinate * directions: The transformation from reference coordinates \hat {\mathbf * x} to real coordinates \mathbf x on each cell is of the form * @f{align*}{ * {\mathbf x}(\hat {\mathbf x}) * = * \begin{pmatrix} * h_x & 0 \\ * 0 & h_y * \end{pmatrix} * \hat{\mathbf x} * + {\mathbf v}_0 * @f} * in 2d, and * @f{align*}{ * {\mathbf x}(\hat {\mathbf x}) * = * \begin{pmatrix} * h_x & 0 & 0 \\ * 0 & h_y & 0 \\ * 0 & 0 & h_z * \end{pmatrix} * \hat{\mathbf x} * + {\mathbf v}_0 * @f} * in 3d, where {\mathbf v}_0 is the bottom left vertex and h_x,h_y,h_z * are the extents of the cell along the axes. * * The class is intended for efficiency, and it does not do a whole lot of * error checking. If you apply this mapping to a cell that does not conform * to the requirements above, you will get strange results. */ template <int dim, int spacedim = dim> class MappingBox : public Mapping<dim, spacedim> { public: MappingBox(const std::vector<BoundingBox<dim>> &local_boxes, const std::map<types::global_cell_index, types::global_cell_index> &polytope_translator); @endcode for documentation, see the Mapping base class @code virtual std::unique_ptr<Mapping<dim, spacedim>> clone() const override; /** * Return @p true because MappingBox preserves vertex * locations. */ virtual bool preserves_vertex_locations() const override; virtual bool is_compatible_with( #if DEAL_II_VERSION_GTE(9, 8, 0) const ReferenceCell<dim> &reference_cell #else const ReferenceCell &reference_cell #endif ) const override; /** * @name Mapping points between reference and real cells * @{ */ @endcode for documentation, see the Mapping base class @code virtual Point<spacedim> transform_unit_to_real_cell( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const Point<dim> &p) const override; @endcode for documentation, see the Mapping base class @code virtual Point<dim> transform_real_to_unit_cell( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const Point<spacedim> &p) const override; @endcode for documentation, see the Mapping base class @code virtual void transform_points_real_to_unit_cell( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const ArrayView<const Point<spacedim>> &real_points, const ArrayView<Point<dim>> &unit_points) const override; /** * @} */ /** * @name Functions to transform tensors from reference to real coordinates * @{ */ @endcode for documentation, see the Mapping base class @code virtual void transform(const ArrayView<const Tensor<1, dim>> &input, const MappingKind kind, const typename Mapping<dim, spacedim>::InternalDataBase &internal, const ArrayView<Tensor<1, spacedim>> &output) const override; @endcode for documentation, see the Mapping base class @code virtual void transform(const ArrayView<const DerivativeForm<1, dim, spacedim>> &input, const MappingKind kind, const typename Mapping<dim, spacedim>::InternalDataBase &internal, const ArrayView<Tensor<2, spacedim>> &output) const override; @endcode for documentation, see the Mapping base class @code virtual void transform(const ArrayView<const Tensor<2, dim>> &input, const MappingKind kind, const typename Mapping<dim, spacedim>::InternalDataBase &internal, const ArrayView<Tensor<2, spacedim>> &output) const override; @endcode for documentation, see the Mapping base class @code virtual void transform(const ArrayView<const DerivativeForm<2, dim, spacedim>> &input, const MappingKind kind, const typename Mapping<dim, spacedim>::InternalDataBase &internal, const ArrayView<Tensor<3, spacedim>> &output) const override; @endcode for documentation, see the Mapping base class @code virtual void transform(const ArrayView<const Tensor<3, dim>> &input, const MappingKind kind, const typename Mapping<dim, spacedim>::InternalDataBase &internal, const ArrayView<Tensor<3, spacedim>> &output) const override; /** * @} */ /** * @name Interface with FEValues * @{ */ /** * Storage for internal data of the mapping. See Mapping::InternalDataBase * for an extensive description. * * This includes data that is computed once when the object is created (in * get_data()) as well as data the class wants to store from between the * call to fill_fe_values(), fill_fe_face_values(), or * fill_fe_subface_values() until possible later calls from the finite * element to functions such as transform(). The latter class of member * variables are marked as 'mutable'. */ class InternalData : public Mapping<dim, spacedim>::InternalDataBase { public: /** * Default constructor. */ InternalData() = default; /** * Constructor that initializes the object with a quadrature. */ InternalData(const Quadrature<dim> &quadrature); @endcode Documentation see Mapping::InternalDataBase. @code virtual void reinit(const UpdateFlags update_flags, const Quadrature<dim> &quadrature) override; /** * Return an estimate (in bytes) for the memory consumption of this object. */ virtual std::size_t memory_consumption() const override; /** * Extents of the last cell we have seen in the coordinate directions, * i.e., <i>h<sub>x</sub></i>, <i>h<sub>y</sub></i>, <i>h<sub>z</sub></i>. */ mutable Tensor<1, dim> cell_extents; /** * Traslation term in F(\hat{x})=J\hat{x} + c. */ mutable Tensor<1, dim> traslation; /** * Reciprocal of the extents of the last cell we have seen in the * coordinate directions, i.e., <i>h<sub>x</sub></i>, * <i>h<sub>y</sub></i>, <i>h<sub>z</sub></i>. */ mutable Tensor<1, dim> inverse_cell_extents; /** * The volume element */ mutable double volume_element; /** * Location of quadrature points of faces or subfaces in 3d with all * possible orientations. Can be accessed with the correct offset provided * via QProjector::DataSetDescriptor. Not needed/used for cells. */ std::vector<Point<dim>> quadrature_points; }; private: @endcode documentation can be found in Mapping::requires_update_flags() @code virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override; @endcode documentation can be found in Mapping::get_data() @code virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> get_data(const UpdateFlags, const Quadrature<dim> &quadrature) const override; using Mapping<dim, spacedim>::get_face_data; @endcode documentation can be found in Mapping::get_subface_data() @code virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> get_subface_data(const UpdateFlags flags, const Quadrature<dim - 1> &quadrature) const override; @endcode documentation can be found in Mapping::fill_fe_values() @code virtual CellSimilarity::Similarity fill_fe_values( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature<dim> &quadrature, const typename Mapping<dim, spacedim>::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const override; using Mapping<dim, spacedim>::fill_fe_face_values; @endcode documentation can be found in Mapping::fill_fe_subface_values() @code virtual void fill_fe_subface_values( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature<dim - 1> &quadrature, const typename Mapping<dim, spacedim>::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const override; @endcode documentation can be found in Mapping::fill_fe_immersed_surface_values() @code virtual void fill_fe_immersed_surface_values( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature, const typename Mapping<dim, spacedim>::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const override; /** * @} */ /** * Update the cell_extents field of the incoming InternalData object with the * size of the incoming cell. */ void update_cell_extents( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const InternalData &data) const; /** * Compute the quadrature points if the UpdateFlags of the incoming * InternalData object say that they should be updated. * * Called from fill_fe_values. */ void maybe_update_cell_quadrature_points( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const InternalData &data, const ArrayView<const Point<dim>> &unit_quadrature_points, std::vector<Point<dim>> &quadrature_points) const; /** * Compute the normal vectors if the UpdateFlags of the incoming InternalData * object say that they should be updated. */ void maybe_update_normal_vectors( const unsigned int face_no, const InternalData &data, std::vector<Tensor<1, dim>> &normal_vectors) const; /** * Since the Jacobian is constant for this mapping all derivatives of the * Jacobian are identically zero. Fill these quantities with zeros if the * corresponding update flags say that they should be updated. */ void maybe_update_jacobian_derivatives( const InternalData &data, const CellSimilarity::Similarity cell_similarity, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const; /** * Compute the volume elements if the UpdateFlags of the incoming * InternalData object say that they should be updated. */ void maybe_update_volume_elements(const InternalData &data) const; /** * Compute the Jacobians if the UpdateFlags of the incoming * InternalData object say that they should be updated. */ void maybe_update_jacobians( const InternalData &data, const CellSimilarity::Similarity cell_similarity, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const; /** * Compute the inverse Jacobians if the UpdateFlags of the incoming * InternalData object say that they should be updated. */ void maybe_update_inverse_jacobians( const InternalData &data, const CellSimilarity::Similarity cell_similarity, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const; /** * Vector of (local) bounding boxes */ std::vector<BoundingBox<dim>> boxes; /** * Map from global cell index to bounding box index */ std::map<types::global_cell_index, types::global_cell_index> polytope_translator; }; /** @} */ DEAL_II_NAMESPACE_CLOSE #endif @endcode <a name="ann-include/poly_utils.h"></a> <h1>Annotated version of include/poly_utils.h</h1> @code /* ----------------------------------------------------------------------------- * * SPDX-License-Identifier: LGPL-2.1-or-later * Copyright (C) 2025 by Marco Feder, Pasquale Claudio Africa, Xinping Gui, * Andrea Cangiani * * This file is part of the deal.II code gallery. * * ----------------------------------------------------------------------------- */ #ifndef poly_utils_h #define poly_utils_h #include <deal.II/base/config.h> #include <deal.II/base/conditional_ostream.h> #include <deal.II/base/point.h> #include <deal.II/base/quadrature.h> #include <deal.II/base/std_cxx20/iota_view.h> #include <deal.II/boost_adaptors/bounding_box.h> #include <deal.II/boost_adaptors/point.h> #include <deal.II/boost_adaptors/segment.h> #include <deal.II/distributed/tria.h> #include <deal.II/dofs/dof_handler.h> #include <deal.II/fe/fe_dgq.h> #include <deal.II/fe/fe_values.h> #include <deal.II/grid/grid_tools.h> #include <deal.II/lac/dynamic_sparsity_pattern.h> #include <deal.II/lac/sparse_matrix.h> #include <deal.II/lac/sparsity_pattern.h> #include <deal.II/lac/sparsity_tools.h> #include <deal.II/lac/trilinos_sparse_matrix.h> #include <deal.II/numerics/vector_tools_common.h> #include <boost/geometry/algorithms/distance.hpp> #include <boost/geometry/index/detail/rtree/utilities/print.hpp> #include <boost/geometry/index/rtree.hpp> #include <boost/geometry/strategies/strategies.hpp> #include <memory> namespace ::PolyUtils::internal { /** * Helper function to compute the position of index @p index in vector @p v. */ inline types::global_cell_index get_index(const std::vector<types::global_cell_index> &v, const types::global_cell_index index) { return std::distance(v.begin(), std::find(v.begin(), v.end(), index)); } /** * Compute the connectivity graph for locally owned regions of a distributed * triangulation. */ template <int dim, int spacedim> void get_face_connectivity_of_cells( const parallel::fullydistributed::Triangulation<dim, spacedim> &triangulation, DynamicSparsityPattern &cell_connectivity, const std::vector<types::global_cell_index> locally_owned_cells) { cell_connectivity.reinit(triangulation.n_locally_owned_active_cells(), triangulation.n_locally_owned_active_cells()); @endcode loop over all cells and their neighbors to build the sparsity pattern. note that it's a bit hard to enter all the connections when a neighbor has children since we would need to find out which of its children is adjacent to the current cell. this problem can be omitted if we only do something if the neighbor has no children – in that case it is either on the same or a coarser level than we are. in return, we have to add entries in both directions for both cells @code for (const auto &cell : triangulation.active_cell_iterators()) { if (cell->is_locally_owned()) { const unsigned int index = cell->active_cell_index(); cell_connectivity.add(get_index(locally_owned_cells, index), get_index(locally_owned_cells, index)); for (auto f : cell->face_indices()) if ((cell->at_boundary(f) == false) && (cell->neighbor(f)->has_children() == false) && cell->neighbor(f)->is_locally_owned()) { const unsigned int other_index = cell->neighbor(f)->active_cell_index(); cell_connectivity.add(get_index(locally_owned_cells, index), get_index(locally_owned_cells, other_index)); cell_connectivity.add(get_index(locally_owned_cells, other_index), get_index(locally_owned_cells, index)); } } } } } // namespace ::PolyUtils::internal namespace ::PolyUtils { template <typename Value, typename Options, typename Translator, typename Box, typename Allocators> struct Rtree_visitor : public boost::geometry::index::detail::rtree::visitor< Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag, true>::type { inline Rtree_visitor( const Translator &translator, unsigned int target_level, std::vector<std::vector<typename Triangulation< boost::geometry::dimension<Box>::value>::active_cell_iterator>> &boxes, std::vector<std::vector<unsigned int>> &csr); /** * An alias that identifies an InternalNode of the tree. */ using InternalNode = typename boost::geometry::index::detail::rtree::internal_node< Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type; /** * An alias that identifies a Leaf of the tree. */ using Leaf = typename boost::geometry::index::detail::rtree::leaf< Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type; /** * Implements the visitor interface for InternalNode objects. If the node * belongs to the level next to @p target_level, then fill the bounding box * vector for that node. */ inline void operator()(const InternalNode &node); /** * Implements the visitor interface for Leaf objects. */ inline void operator()(const Leaf &); /** * Translator interface, required by the boost implementation of the rtree. */ const Translator &translator; /** * Store the level we are currently visiting. */ size_t level; /** * Index used to keep track of the number of different visited nodes during * recursion/ */ size_t node_counter; size_t next_level_leafs_processed; /** * The level where children are living. * Before: "we want to extract from the RTree object." */ const size_t target_level; /** * A reference to the input vector of vector of BoundingBox objects. This * vector v has the following property: v[i] = vector with all * of the BoundingBox bounded by the i-th node of the Rtree. */ std::vector<std::vector<typename Triangulation< boost::geometry::dimension<Box>::value>::active_cell_iterator>> &agglomerates; std::vector<std::vector<unsigned int>> &row_ptr; }; template <typename Value, typename Options, typename Translator, typename Box, typename Allocators> Rtree_visitor<Value, Options, Translator, Box, Allocators>::Rtree_visitor( const Translator &translator, const unsigned int target_level, std::vector<std::vector<typename Triangulation< boost::geometry::dimension<Box>::value>::active_cell_iterator>> &bb_in_boxes, std::vector<std::vector<unsigned int>> &csr) : translator(translator) , level(0) , node_counter(0) , next_level_leafs_processed(0) , target_level(target_level) , agglomerates(bb_in_boxes) , row_ptr(csr) {} template <typename Value, typename Options, typename Translator, typename Box, typename Allocators> void Rtree_visitor<Value, Options, Translator, Box, Allocators>::operator()( const Rtree_visitor::InternalNode &node) { using elements_type = typename boost::geometry::index::detail::rtree::elements_type< InternalNode>::type; // pairs of bounding box and pointer to child @endcode node @code const elements_type &elements = boost::geometry::index::detail::rtree::elements(node); if (level < target_level) { size_t level_backup = level; ++level; for (typename elements_type::const_iterator it = elements.begin(); it != elements.end(); ++it) { boost::geometry::index::detail::rtree::apply_visitor(*this, *it->second); } level = level_backup; } else if (level == target_level) { @endcode const unsigned int n_children = elements.size(); @code const auto offset = agglomerates.size(); agglomerates.resize(offset + 1); row_ptr.resize(row_ptr.size() + 1); next_level_leafs_processed = 0; row_ptr.back().push_back( next_level_leafs_processed); // convention: row_ptr[0]=0 size_t level_backup = level; ++level; for (const auto &child : elements) { boost::geometry::index::detail::rtree::apply_visitor(*this, *child.second); } @endcode Done with node number 'node_counter' @code ++node_counter; // visited all children of an internal node level = level_backup; } else if (level > target_level) { @endcode Keep visiting until you go to the leafs. @code size_t level_backup = level; ++level; for (const auto &child : elements) { boost::geometry::index::detail::rtree::apply_visitor(*this, *child.second); } level = level_backup; row_ptr[node_counter].push_back(next_level_leafs_processed); } } template <typename Value, typename Options, typename Translator, typename Box, typename Allocators> void Rtree_visitor<Value, Options, Translator, Box, Allocators>::operator()( const Rtree_visitor::Leaf &leaf) { using elements_type = typename boost::geometry::index::detail::rtree::elements_type< Leaf>::type; // pairs of bounding box and pointer to child node const elements_type &elements = boost::geometry::index::detail::rtree::elements(leaf); for (const auto &it : elements) { agglomerates[node_counter].push_back(it.second); } next_level_leafs_processed += elements.size(); } template <typename T> inline constexpr T constexpr_pow(T num, unsigned int pow) { return (pow >= sizeof(unsigned int) * 8) ? 0 : pow == 0 ? 1 : num * constexpr_pow(num, pow - 1); } namespace internal { /** * Same as the public free function with the same name, but storing * explicitly the interpolation matrix and performing interpolation through * matrix-vector product. */ template <int dim, int spacedim, typename VectorType> void interpolate_to_fine_grid( const AgglomerationHandler<dim, spacedim> &agglomeration_handler, VectorType &dst, const VectorType &src) { Assert((dim == spacedim), ExcNotImplemented()); Assert( dst.size() == 0, ExcMessage( "The destination vector must the empt upon calling this function.")); using NumberType = typename VectorType::value_type; constexpr bool is_trilinos_vector = std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector>; using MatrixType = std::conditional_t<is_trilinos_vector, TrilinosWrappers::SparseMatrix, SparseMatrix<NumberType>>; MatrixType interpolation_matrix; [[maybe_unused]] typename std::conditional_t<!is_trilinos_vector, SparsityPattern, void *> sp; @endcode Get some info from the handler @code const DoFHandler<dim, spacedim> &agglo_dh = agglomeration_handler.agglo_dh; DoFHandler<dim, spacedim> *output_dh = const_cast<DoFHandler<dim, spacedim> *>( &agglomeration_handler.output_dh); const FiniteElement<dim, spacedim> &fe = agglomeration_handler.get_fe(); const Mapping<dim> &mapping = agglomeration_handler.get_mapping(); const Triangulation<dim, spacedim> &tria = agglomeration_handler.get_triangulation(); const auto &bboxes = agglomeration_handler.get_local_bboxes(); std::unique_ptr<FiniteElement<dim>> output_fe; if (tria.all_reference_cells_are_hyper_cube()) output_fe = std::make_unique<FE_DGQ<dim>>(fe.degree); else if (tria.all_reference_cells_are_simplex()) output_fe = std::make_unique<FE_SimplexDGP<dim>>(fe.degree); else AssertThrow(false, ExcNotImplemented()); @endcode Setup an auxiliary DoFHandler for output purposes @code output_dh->reinit(tria); output_dh->distribute_dofs(*output_fe); const IndexSet &locally_owned_dofs = output_dh->locally_owned_dofs(); const IndexSet locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(*output_dh); const IndexSet &locally_owned_dofs_agglo = agglo_dh.locally_owned_dofs(); DynamicSparsityPattern dsp(output_dh->n_dofs(), agglo_dh.n_dofs(), locally_relevant_dofs); std::vector<types::global_dof_index> agglo_dof_indices(fe.dofs_per_cell); std::vector<types::global_dof_index> standard_dof_indices( fe.dofs_per_cell); std::vector<types::global_dof_index> output_dof_indices( output_fe->dofs_per_cell); Quadrature<dim> quad(output_fe->get_unit_support_points()); FEValues<dim, spacedim> output_fe_values(mapping, *output_fe, quad, update_quadrature_points); for (const auto &cell : agglo_dh.active_cell_iterators()) if (cell->is_locally_owned()) { if (agglomeration_handler.is_master_cell(cell)) { auto slaves = agglomeration_handler.get_slaves_of_idx( cell->active_cell_index()); slaves.emplace_back(cell); cell->get_dof_indices(agglo_dof_indices); for (const auto &slave : slaves) { @endcode addd master-slave relationship @code const auto slave_output = slave->as_dof_handler_iterator(*output_dh); slave_output->get_dof_indices(output_dof_indices); for (const auto row : output_dof_indices) dsp.add_entries(row, agglo_dof_indices.begin(), agglo_dof_indices.end()); } } } const auto assemble_interpolation_matrix = [&]() { FullMatrix<NumberType> local_matrix(fe.dofs_per_cell, fe.dofs_per_cell); std::vector<Point<dim>> reference_q_points(fe.dofs_per_cell); @endcode Dummy AffineConstraints, only needed for loc2glb @code AffineConstraints<NumberType> c; c.close(); for (const auto &cell : agglo_dh.active_cell_iterators()) if (cell->is_locally_owned()) { if (agglomeration_handler.is_master_cell(cell)) { auto slaves = agglomeration_handler.get_slaves_of_idx( cell->active_cell_index()); slaves.emplace_back(cell); cell->get_dof_indices(agglo_dof_indices); const types::global_cell_index polytope_index = agglomeration_handler.cell_to_polytope_index(cell); @endcode Get the box of this agglomerate. @code const BoundingBox<dim> &box = bboxes[polytope_index]; for (const auto &slave : slaves) { @endcode add master-slave relationship @code const auto slave_output = slave->as_dof_handler_iterator(*output_dh); slave_output->get_dof_indices(output_dof_indices); output_fe_values.reinit(slave_output); local_matrix = 0.; const auto &q_points = output_fe_values.get_quadrature_points(); for (const auto i : output_fe_values.dof_indices()) { const auto &p = box.real_to_unit(q_points[i]); for (const auto j : output_fe_values.dof_indices()) { local_matrix(i, j) = fe.shape_value(j, p); } } c.distribute_local_to_global(local_matrix, output_dof_indices, agglo_dof_indices, interpolation_matrix); } } } }; if constexpr (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix>) { const MPI_Comm &communicator = tria.get_mpi_communicator(); SparsityTools::distribute_sparsity_pattern(dsp, locally_owned_dofs, communicator, locally_relevant_dofs); interpolation_matrix.reinit(locally_owned_dofs, locally_owned_dofs_agglo, dsp, communicator); dst.reinit(locally_owned_dofs); assemble_interpolation_matrix(); } else if constexpr (std::is_same_v<MatrixType, SparseMatrix<NumberType>>) { sp.copy_from(dsp); interpolation_matrix.reinit(sp); dst.reinit(output_dh->n_dofs()); assemble_interpolation_matrix(); } else { @endcode PETSc, LA::d::v options not implemented. @code (void)agglomeration_handler; (void)dst; (void)src; AssertThrow(false, ExcNotImplemented()); } @endcode If tria is distributed @code if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>( &tria) != nullptr) interpolation_matrix.compress(VectorOperation::add); @endcode Finally, perform the interpolation. @code interpolation_matrix.vmult(dst, src); } } // namespace internal /** * Given a vector @p src, typically the solution stemming after the * agglomerate problem has been solved, this function interpolates @p src * onto the finer grid and stores the result in vector @p dst. The last * argument @p on_the_fly does not build any interpolation matrix and allows * computing the entries in @p dst in a matrix-free fashion. * * @note Supported parallel types are TrilinosWrappers::SparseMatrix and * TrilinosWrappers::MPI::Vector. */ template <int dim, int spacedim, typename VectorType> void interpolate_to_fine_grid( const AgglomerationHandler<dim, spacedim> &agglomeration_handler, VectorType &dst, const VectorType &src, const bool on_the_fly = true) { Assert((dim == spacedim), ExcNotImplemented()); Assert( dst.size() == 0, ExcMessage( "The destination vector must the empt upon calling this function.")); using NumberType = typename VectorType::value_type; static constexpr bool is_trilinos_vector = std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector>; static constexpr bool is_supported_vector = std::is_same_v<VectorType, Vector<NumberType>> || is_trilinos_vector; static_assert(is_supported_vector); @endcode First, check for an easy return @code if (on_the_fly == false) { return internal::interpolate_to_fine_grid(agglomeration_handler, dst, src); } else { @endcode otherwise, do not create any matrix @code if (!agglomeration_handler.used_fe_collection()) { @endcode Original version: handle case without hp::FECollection @code const Triangulation<dim, spacedim> &tria = agglomeration_handler.get_triangulation(); const Mapping<dim> &mapping = agglomeration_handler.get_mapping(); const FiniteElement<dim, spacedim> &original_fe = agglomeration_handler.get_fe(); @endcode We use DGQ (on tensor-product meshes) or DGP (on simplex meshes) nodal elements of the same degree as the ones in the agglomeration handler to interpolate the solution onto the finer grid. @code std::unique_ptr<FiniteElement<dim>> output_fe; if (tria.all_reference_cells_are_hyper_cube()) output_fe = std::make_unique<FE_DGQ<dim>>(original_fe.degree); else if (tria.all_reference_cells_are_simplex()) output_fe = std::make_unique<FE_SimplexDGP<dim>>(original_fe.degree); else AssertThrow(false, ExcNotImplemented()); DoFHandler<dim> &output_dh = const_cast<DoFHandler<dim> &>(agglomeration_handler.output_dh); output_dh.reinit(tria); output_dh.distribute_dofs(*output_fe); if constexpr (std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector>) { const IndexSet &locally_owned_dofs = output_dh.locally_owned_dofs(); dst.reinit(locally_owned_dofs); } else if constexpr (std::is_same_v<VectorType, Vector<NumberType>>) { dst.reinit(output_dh.n_dofs()); } else { @endcode PETSc, LA::d::v options not implemented. @code (void)agglomeration_handler; (void)dst; (void)src; AssertThrow(false, ExcNotImplemented()); } const unsigned int dofs_per_cell = agglomeration_handler.n_dofs_per_cell(); const unsigned int output_dofs_per_cell = output_fe->n_dofs_per_cell(); Quadrature<dim> quad(output_fe->get_unit_support_points()); FEValues<dim> output_fe_values(mapping, *output_fe, quad, update_quadrature_points); std::vector<types::global_dof_index> local_dof_indices( dofs_per_cell); std::vector<types::global_dof_index> local_dof_indices_output( output_dofs_per_cell); const auto &bboxes = agglomeration_handler.get_local_bboxes(); for (const auto &polytope : agglomeration_handler.polytope_iterators()) { if (polytope->is_locally_owned()) { polytope->get_dof_indices(local_dof_indices); const BoundingBox<dim> &box = bboxes[polytope->index()]; const auto &deal_cells = polytope->get_agglomerate(); // fine deal.II cells for (const auto &cell : deal_cells) { const auto slave_output = cell->as_dof_handler_iterator( agglomeration_handler.output_dh); slave_output->get_dof_indices(local_dof_indices_output); output_fe_values.reinit(slave_output); const auto &qpoints = output_fe_values.get_quadrature_points(); for (unsigned int j = 0; j < output_dofs_per_cell; ++j) { const auto &ref_qpoint = box.real_to_unit(qpoints[j]); for (unsigned int i = 0; i < dofs_per_cell; ++i) dst(local_dof_indices_output[j]) += src(local_dof_indices[i]) * original_fe.shape_value(i, ref_qpoint); } } } } } else { @endcode Handle the hp::FECollection case @code const Triangulation<dim, spacedim> &tria = agglomeration_handler.get_triangulation(); const Mapping<dim> &mapping = agglomeration_handler.get_mapping(); const hp::FECollection<dim, spacedim> &original_fe_collection = agglomeration_handler.get_fe_collection(); @endcode We use DGQ (on tensor-product meshes) or DGP (on simplex meshes) nodal elements of the same degree as the ones in the agglomeration handler to interpolate the solution onto the finer grid. @code hp::FECollection<dim, spacedim> output_fe_collection; Assert(original_fe_collection[0].n_components() >= 1, ExcMessage("Invalid FE: must have at least one component.")); if (original_fe_collection[0].n_components() == 1) { @endcode Scalar case @code for (unsigned int i = 0; i < original_fe_collection.size(); ++i) { std::unique_ptr<FiniteElement<dim>> output_fe; if (tria.all_reference_cells_are_hyper_cube()) output_fe = std::make_unique<FE_DGQ<dim>>( original_fe_collection[i].degree); else if (tria.all_reference_cells_are_simplex()) output_fe = std::make_unique<FE_SimplexDGP<dim>>( original_fe_collection[i].degree); else AssertThrow(false, ExcNotImplemented()); output_fe_collection.push_back(*output_fe); } } else if (original_fe_collection[0].n_components() > 1) { @endcode System case @code for (unsigned int i = 0; i < original_fe_collection.size(); ++i) { std::vector<const FiniteElement<dim, spacedim> *> base_elements; std::vector<unsigned int> multiplicities; for (unsigned int b = 0; b < original_fe_collection[i].n_base_elements(); ++b) { if (dynamic_cast<const FE_Nothing<dim> *>( &original_fe_collection[i].base_element(b))) base_elements.push_back( new FE_Nothing<dim, spacedim>()); else { if (tria.all_reference_cells_are_hyper_cube()) base_elements.push_back(new FE_DGQ<dim, spacedim>( original_fe_collection[i] .base_element(b) .degree)); else if (tria.all_reference_cells_are_simplex()) base_elements.push_back( new FE_SimplexDGP<dim, spacedim>( original_fe_collection[i] .base_element(b) .degree)); else AssertThrow(false, ExcNotImplemented()); } multiplicities.push_back( original_fe_collection[i].element_multiplicity(b)); } FESystem<dim, spacedim> output_fe_system(base_elements, multiplicities); for (const auto *ptr : base_elements) delete ptr; output_fe_collection.push_back(output_fe_system); } } DoFHandler<dim> &output_dh = const_cast<DoFHandler<dim> &>(agglomeration_handler.output_dh); output_dh.reinit(tria); for (const auto &polytope : agglomeration_handler.polytope_iterators()) { if (polytope->is_locally_owned()) { const auto &deal_cells = polytope->get_agglomerate(); // fine deal.II cells const unsigned int active_fe_idx = polytope->active_fe_index(); for (const auto &cell : deal_cells) { const typename DoFHandler<dim>::active_cell_iterator slave_cell_dh_iterator = cell->as_dof_handler_iterator(output_dh); slave_cell_dh_iterator->set_active_fe_index( active_fe_idx); } } } output_dh.distribute_dofs(output_fe_collection); if constexpr (std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector>) { const IndexSet &locally_owned_dofs = output_dh.locally_owned_dofs(); dst.reinit(locally_owned_dofs); } else if constexpr (std::is_same_v<VectorType, Vector<NumberType>>) { dst.reinit(output_dh.n_dofs()); } else { @endcode PETSc, LA::d::v options not implemented. @code (void)agglomeration_handler; (void)dst; (void)src; AssertThrow(false, ExcNotImplemented()); } const auto &bboxes = agglomeration_handler.get_local_bboxes(); for (const auto &polytope : agglomeration_handler.polytope_iterators()) { if (polytope->is_locally_owned()) { const unsigned int active_fe_idx = polytope->active_fe_index(); const unsigned int dofs_per_cell = polytope->get_fe().dofs_per_cell; const unsigned int output_dofs_per_cell = output_fe_collection[active_fe_idx].n_dofs_per_cell(); Quadrature<dim> quad(output_fe_collection[active_fe_idx] .get_unit_support_points()); FEValues<dim> output_fe_values( mapping, output_fe_collection[active_fe_idx], quad, update_quadrature_points); std::vector<types::global_dof_index> local_dof_indices( dofs_per_cell); std::vector<types::global_dof_index> local_dof_indices_output(output_dofs_per_cell); polytope->get_dof_indices(local_dof_indices); const BoundingBox<dim> &box = bboxes[polytope->index()]; const auto &deal_cells = polytope->get_agglomerate(); // fine deal.II cells for (const auto &cell : deal_cells) { const auto slave_output = cell->as_dof_handler_iterator( agglomeration_handler.output_dh); slave_output->get_dof_indices(local_dof_indices_output); output_fe_values.reinit(slave_output); const auto &qpoints = output_fe_values.get_quadrature_points(); for (unsigned int j = 0; j < output_dofs_per_cell; ++j) { const unsigned int component_idx_of_this_dof = slave_output->get_fe() .system_to_component_index(j) .first; const auto &ref_qpoint = box.real_to_unit(qpoints[j]); for (unsigned int i = 0; i < dofs_per_cell; ++i) dst(local_dof_indices_output[j]) += src(local_dof_indices[i]) * original_fe_collection[active_fe_idx] .shape_value_component( i, ref_qpoint, component_idx_of_this_dof); } } } } } } } /** * Similar to VectorTools::compute_global_error(), but customized for * polytopic elements. Aside from the solution vector and a reference * function, this function takes in addition a vector @p norms with types * VectorTools::NormType to be computed and later stored in the last * argument @p global_errors. * In case of a parallel vector, the local errors are collected over each * processor and later a classical reduction operation is performed. */ template <int dim, typename Number, typename VectorType> void compute_global_error(const AgglomerationHandler<dim> &agglomeration_handler, const VectorType &solution, const Function<dim, Number> &exact_solution, const std::vector<VectorTools::NormType> &norms, std::vector<double> &global_errors) { Assert(solution.size() > 0, ExcNotImplemented( "Solution vector must be non-empty upon calling this function.")); Assert(std::any_of(norms.cbegin(), norms.cend(), [](VectorTools::NormType norm_type) { return (norm_type == VectorTools::NormType::H1_seminorm || norm_type == VectorTools::NormType::L2_norm); }), ExcMessage("Norm type not supported")); global_errors.resize(norms.size()); std::fill(global_errors.begin(), global_errors.end(), 0.); @endcode Vector storing errors local to the current processor. @code std::vector<double> local_errors(norms.size()); std::fill(local_errors.begin(), local_errors.end(), 0.); @endcode Get some info from the handler @code const unsigned int dofs_per_cell = agglomeration_handler.n_dofs_per_cell(); const bool compute_semi_H1 = std::any_of(norms.cbegin(), norms.cend(), [](VectorTools::NormType norm_type) { return norm_type == VectorTools::NormType::H1_seminorm; }); std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell); for (const auto &polytope : agglomeration_handler.polytope_iterators()) { if (polytope->is_locally_owned()) { const auto &agglo_values = agglomeration_handler.reinit(polytope); polytope->get_dof_indices(local_dof_indices); const auto &q_points = agglo_values.get_quadrature_points(); const unsigned int n_qpoints = q_points.size(); std::vector<double> analyical_sol_at_qpoints(n_qpoints); exact_solution.value_list(q_points, analyical_sol_at_qpoints); std::vector<Tensor<1, dim>> grad_analyical_sol_at_qpoints( n_qpoints); if (compute_semi_H1) exact_solution.gradient_list(q_points, grad_analyical_sol_at_qpoints); for (unsigned int q_index : agglo_values.quadrature_point_indices()) { double solution_at_qpoint = 0.; Tensor<1, dim> grad_solution_at_qpoint; for (unsigned int i = 0; i < dofs_per_cell; ++i) { solution_at_qpoint += solution(local_dof_indices[i]) * agglo_values.shape_value(i, q_index); if (compute_semi_H1) grad_solution_at_qpoint += solution(local_dof_indices[i]) * agglo_values.shape_grad(i, q_index); } @endcode L2 @code local_errors[0] += std::pow((analyical_sol_at_qpoints[q_index] - solution_at_qpoint), 2) * agglo_values.JxW(q_index); @endcode H1 seminorm @code if (compute_semi_H1) for (unsigned int d = 0; d < dim; ++d) local_errors[1] += std::pow((grad_analyical_sol_at_qpoints[q_index][d] - grad_solution_at_qpoint[d]), 2) * agglo_values.JxW(q_index); } } } @endcode Perform reduction and take sqrt of each error @code global_errors[0] = Utilities::MPI::reduce<double>( local_errors[0], agglomeration_handler.get_triangulation().get_mpi_communicator(), [](const double a, const double b) { return a + b; }); global_errors[0] = std::sqrt(global_errors[0]); if (compute_semi_H1) { global_errors[1] = Utilities::MPI::reduce<double>( local_errors[1], agglomeration_handler.get_triangulation().get_mpi_communicator(), [](const double a, const double b) { return a + b; }); global_errors[1] = std::sqrt(global_errors[1]); } } /** * Utility function that builds the multilevel hierarchy from the tree level * @p starting_level. This function fills the vector of * @p AgglomerationHandlers objects by distributing degrees of freedom on * each level of the hierarchy. It returns the total number of levels in the * hierarchy. */ template <int dim> unsigned int construct_agglomerated_levels( const Triangulation<dim> &tria, std::vector<std::unique_ptr<AgglomerationHandler<dim>>> &agglomeration_handlers, const FE_DGQ<dim> &fe_dg, const Mapping<dim> &mapping, const unsigned int starting_tree_level) { const auto parallel_tria = dynamic_cast<const parallel::TriangulationBase<dim> *>(&tria); GridTools::Cache<dim> cached_tria(tria); Assert(parallel_tria->n_active_cells() > 0, ExcInternalError()); const MPI_Comm comm = parallel_tria->get_mpi_communicator(); ConditionalOStream pcout(std::cout, (Utilities::MPI::this_mpi_process(comm) == 0)); @endcode Start building R-tree @code namespace bgi = boost::geometry::index; static constexpr unsigned int max_elem_per_node = constexpr_pow(2, dim); // 2^dim std::vector<std::pair<BoundingBox<dim>, typename Triangulation<dim>::active_cell_iterator>> boxes(parallel_tria->n_locally_owned_active_cells()); unsigned int i = 0; for (const auto &cell : parallel_tria->active_cell_iterators()) if (cell->is_locally_owned()) boxes[i++] = std::make_pair(mapping.get_bounding_box(cell), cell); auto tree = pack_rtree<bgi::rstar<max_elem_per_node>>(boxes); Assert(n_levels(tree) >= 2, ExcMessage("At least two levels are needed.")); pcout << "Total number of available levels: " << n_levels(tree) << std::endl; pcout << "Starting level: " << starting_tree_level << std::endl; const unsigned int total_tree_levels = n_levels(tree) - starting_tree_level + 1; @endcode Resize the agglomeration handlers to the right size @code agglomeration_handlers.resize(total_tree_levels); @endcode Loop through the available levels and set AgglomerationHandlers up. @code for (unsigned int extraction_level = starting_tree_level; extraction_level <= n_levels(tree); ++extraction_level) { agglomeration_handlers[extraction_level - starting_tree_level] = std::make_unique<AgglomerationHandler<dim>>(cached_tria); CellsAgglomerator<dim, decltype(tree)> agglomerator{tree, extraction_level}; const auto agglomerates = agglomerator.extract_agglomerates(); agglomeration_handlers[extraction_level - starting_tree_level] ->connect_hierarchy(agglomerator); @endcode Flag elements for agglomeration @code unsigned int agglo_index = 0; for (unsigned int i = 0; i < agglomerates.size(); ++i) { const auto &agglo = agglomerates[i]; // i-th agglomerate for (const auto &el : agglo) { el->set_material_id(agglo_index); } ++agglo_index; } const unsigned int n_local_agglomerates = agglo_index; unsigned int total_agglomerates = Utilities::MPI::sum(n_local_agglomerates, comm); pcout << "Total agglomerates per (tree) level: " << extraction_level << ": " << total_agglomerates << std::endl; @endcode Now, perform agglomeration within each locally owned partition @code std::vector< std::vector<typename Triangulation<dim>::active_cell_iterator>> cells_per_subdomain(n_local_agglomerates); for (const auto &cell : parallel_tria->active_cell_iterators()) if (cell->is_locally_owned()) cells_per_subdomain[cell->material_id()].push_back(cell); @endcode For every subdomain, agglomerate elements together @code for (std::size_t i = 0; i < cells_per_subdomain.size(); ++i) agglomeration_handlers[extraction_level - starting_tree_level] ->define_agglomerate(cells_per_subdomain[i]); agglomeration_handlers[extraction_level - starting_tree_level] ->initialize_fe_values(QGauss<dim>(fe_dg.degree + 1), update_values | update_gradients | update_JxW_values | update_quadrature_points, QGauss<dim - 1>(fe_dg.degree + 1), update_JxW_values); agglomeration_handlers[extraction_level - starting_tree_level] ->distribute_agglomerated_dofs(fe_dg); } return total_tree_levels; } /** * Utility to compute jump terms when the interface is locally owned, i.e. * both elements are locally owned. */ template <int dim> void assemble_local_jumps_and_averages(FullMatrix<double> &M11, FullMatrix<double> &M12, FullMatrix<double> &M21, FullMatrix<double> &M22, const FEValuesBase<dim> &fe_faces0, const FEValuesBase<dim> &fe_faces1, const double penalty_constant, const double h_f) { const std::vector<Tensor<1, dim>> &normals = fe_faces0.get_normal_vectors(); const unsigned int dofs_per_cell = M11.m(); // size of local matrices equals the #DoFs for (unsigned int q_index : fe_faces0.quadrature_point_indices()) { const Tensor<1, dim> &normal = normals[q_index]; for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { M11(i, j) += (-0.5 * fe_faces0.shape_grad(i, q_index) * normal * fe_faces0.shape_value(j, q_index) - 0.5 * fe_faces0.shape_grad(j, q_index) * normal * fe_faces0.shape_value(i, q_index) + (penalty_constant / h_f) * fe_faces0.shape_value(i, q_index) * fe_faces0.shape_value(j, q_index)) * fe_faces0.JxW(q_index); M12(i, j) += (0.5 * fe_faces0.shape_grad(i, q_index) * normal * fe_faces1.shape_value(j, q_index) - 0.5 * fe_faces1.shape_grad(j, q_index) * normal * fe_faces0.shape_value(i, q_index) - (penalty_constant / h_f) * fe_faces0.shape_value(i, q_index) * fe_faces1.shape_value(j, q_index)) * fe_faces1.JxW(q_index); M21(i, j) += (-0.5 * fe_faces1.shape_grad(i, q_index) * normal * fe_faces0.shape_value(j, q_index) + 0.5 * fe_faces0.shape_grad(j, q_index) * normal * fe_faces1.shape_value(i, q_index) - (penalty_constant / h_f) * fe_faces1.shape_value(i, q_index) * fe_faces0.shape_value(j, q_index)) * fe_faces1.JxW(q_index); M22(i, j) += (0.5 * fe_faces1.shape_grad(i, q_index) * normal * fe_faces1.shape_value(j, q_index) + 0.5 * fe_faces1.shape_grad(j, q_index) * normal * fe_faces1.shape_value(i, q_index) + (penalty_constant / h_f) * fe_faces1.shape_value(i, q_index) * fe_faces1.shape_value(j, q_index)) * fe_faces1.JxW(q_index); } } } } /** * Same as above, but for a ghosted neighbor. */ template <int dim> void assemble_local_jumps_and_averages_ghost( FullMatrix<double> &M11, FullMatrix<double> &M12, FullMatrix<double> &M21, FullMatrix<double> &M22, const FEValuesBase<dim> &fe_faces0, const std::vector<std::vector<double>> &recv_values, const std::vector<std::vector<Tensor<1, dim>>> &recv_gradients, const std::vector<double> &recv_jxws, const double penalty_constant, const double h_f) { Assert( (recv_values.size() > 0 && recv_gradients.size() && recv_jxws.size()), ExcMessage("Not possible to assemble jumps and averages at a ghosted " "interface.")); const unsigned int dofs_per_cell = M11.m(); const std::vector<Tensor<1, dim>> &normals = fe_faces0.get_normal_vectors(); for (unsigned int q_index : fe_faces0.quadrature_point_indices()) { const Tensor<1, dim> &normal = normals[q_index]; for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { M11(i, j) += (-0.5 * fe_faces0.shape_grad(i, q_index) * normal * fe_faces0.shape_value(j, q_index) - 0.5 * fe_faces0.shape_grad(j, q_index) * normal * fe_faces0.shape_value(i, q_index) + (penalty_constant / h_f) * fe_faces0.shape_value(i, q_index) * fe_faces0.shape_value(j, q_index)) * fe_faces0.JxW(q_index); M12(i, j) += (0.5 * fe_faces0.shape_grad(i, q_index) * normal * recv_values[j][q_index] - 0.5 * recv_gradients[j][q_index] * normal * fe_faces0.shape_value(i, q_index) - (penalty_constant / h_f) * fe_faces0.shape_value(i, q_index) * recv_values[j][q_index]) * recv_jxws[q_index]; M21(i, j) += (-0.5 * recv_gradients[i][q_index] * normal * fe_faces0.shape_value(j, q_index) + 0.5 * fe_faces0.shape_grad(j, q_index) * normal * recv_values[i][q_index] - (penalty_constant / h_f) * recv_values[i][q_index] * fe_faces0.shape_value(j, q_index)) * recv_jxws[q_index]; M22(i, j) += (0.5 * recv_gradients[i][q_index] * normal * recv_values[j][q_index] + 0.5 * recv_gradients[j][q_index] * normal * recv_values[i][q_index] + (penalty_constant / h_f) * recv_values[i][q_index] * recv_values[j][q_index]) * recv_jxws[q_index]; } } } } /** * Utility function to assemble the SIPDG Laplace matrix. * @note Supported matrix types are Trilinos types and native SparseMatrix * objects provided by deal.II. */ template <int dim, typename MatrixType> void assemble_dg_matrix(MatrixType &system_matrix, const FiniteElement<dim> &fe_dg, const AgglomerationHandler<dim> &ah) { static_assert( (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix> || std::is_same_v<MatrixType, SparseMatrix<typename MatrixType::value_type>>)); Assert((dynamic_cast<const FE_DGQ<dim> *>(&fe_dg) || dynamic_cast<const FE_DGP<dim> *>(&fe_dg) || dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_dg)), ExcMessage("FE type not supported.")); AffineConstraints constraints; constraints.close(); const double penalty_constant = 10 * (fe_dg.degree + dim) * (fe_dg.degree + 1); TrilinosWrappers::SparsityPattern dsp; const_cast<AgglomerationHandler<dim> &>(ah) .create_agglomeration_sparsity_pattern(dsp); system_matrix.reinit(dsp); const unsigned int dofs_per_cell = fe_dg.n_dofs_per_cell(); FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell); FullMatrix<double> M11(dofs_per_cell, dofs_per_cell); FullMatrix<double> M12(dofs_per_cell, dofs_per_cell); FullMatrix<double> M21(dofs_per_cell, dofs_per_cell); FullMatrix<double> M22(dofs_per_cell, dofs_per_cell); std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell); std::vector<types::global_dof_index> local_dof_indices_neighbor( dofs_per_cell); for (const auto &polytope : ah.polytope_iterators()) { if (polytope->is_locally_owned()) { cell_matrix = 0.; const auto &agglo_values = ah.reinit(polytope); for (unsigned int q_index : agglo_values.quadrature_point_indices()) { for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { cell_matrix(i, j) += agglo_values.shape_grad(i, q_index) * agglo_values.shape_grad(j, q_index) * agglo_values.JxW(q_index); } } } @endcode get volumetric DoFs @code polytope->get_dof_indices(local_dof_indices); @endcode Assemble face terms @code unsigned int n_faces = polytope->n_faces(); const double h_f = polytope->diameter(); for (unsigned int f = 0; f < n_faces; ++f) { if (polytope->at_boundary(f)) { @endcode Get normal vectors seen from each agglomeration. @code const auto &fe_face = ah.reinit(polytope, f); const auto &normals = fe_face.get_normal_vectors(); for (unsigned int q_index : fe_face.quadrature_point_indices()) { const Tensor<1, dim> &normal = normals[q_index]; for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { cell_matrix(i, j) += (-fe_face.shape_value(i, q_index) * fe_face.shape_grad(j, q_index) * normal - fe_face.shape_grad(i, q_index) * normal * fe_face.shape_value(j, q_index) + (penalty_constant / h_f) * fe_face.shape_value(i, q_index) * fe_face.shape_value(j, q_index)) * fe_face.JxW(q_index); } } } } else { const auto &neigh_polytope = polytope->neighbor(f); if (polytope->id() < neigh_polytope->id()) { unsigned int nofn = polytope->neighbor_of_agglomerated_neighbor(f); Assert(neigh_polytope->neighbor(nofn)->id() == polytope->id(), ExcMessage("Mismatch.")); const auto &fe_faces = ah.reinit_interface( polytope, neigh_polytope, f, nofn); const auto &fe_faces0 = fe_faces.first; if (neigh_polytope->is_locally_owned()) { @endcode use both fevalues @code const auto &fe_faces1 = fe_faces.second; M11 = 0.; M12 = 0.; M21 = 0.; M22 = 0.; assemble_local_jumps_and_averages(M11, M12, M21, M22, fe_faces0, fe_faces1, penalty_constant, h_f); @endcode distribute DoFs accordingly fluxes @code neigh_polytope->get_dof_indices( local_dof_indices_neighbor); constraints.distribute_local_to_global( M11, local_dof_indices, system_matrix); constraints.distribute_local_to_global( M12, local_dof_indices, local_dof_indices_neighbor, system_matrix); constraints.distribute_local_to_global( M21, local_dof_indices_neighbor, local_dof_indices, system_matrix); constraints.distribute_local_to_global( M22, local_dof_indices_neighbor, system_matrix); } else { @endcode neigh polytope is ghosted, so retrieve necessary metadata. @code types::subdomain_id neigh_rank = neigh_polytope->subdomain_id(); const auto &recv_jxws = ah.recv_jxws.at(neigh_rank) .at({neigh_polytope->id(), nofn}); const auto &recv_values = ah.recv_values.at(neigh_rank) .at({neigh_polytope->id(), nofn}); const auto &recv_gradients = ah.recv_gradients.at(neigh_rank) .at({neigh_polytope->id(), nofn}); M11 = 0.; M12 = 0.; M21 = 0.; M22 = 0.; @endcode there's no FEFaceValues on the other side (it's ghosted), so we just pass the actual data we have recevied from the neighboring ghosted polytope @code assemble_local_jumps_and_averages_ghost( M11, M12, M21, M22, fe_faces0, recv_values, recv_gradients, recv_jxws, penalty_constant, h_f); @endcode distribute DoFs accordingly fluxes @code neigh_polytope->get_dof_indices( local_dof_indices_neighbor); constraints.distribute_local_to_global( M11, local_dof_indices, system_matrix); constraints.distribute_local_to_global( M12, local_dof_indices, local_dof_indices_neighbor, system_matrix); constraints.distribute_local_to_global( M21, local_dof_indices_neighbor, local_dof_indices, system_matrix); constraints.distribute_local_to_global( M22, local_dof_indices_neighbor, system_matrix); } // ghosted polytope case } // only once } // internal face } // face loop constraints.distribute_local_to_global(cell_matrix, local_dof_indices, system_matrix); } // locally owned polytopes } system_matrix.compress(VectorOperation::add); } /** * Compute SIPDG matrix as well as rhs vector. * @note Hardcoded for f=1 and simplex elements. * TODO: Pass Function object for boundary conditions and forcing term. */ template <int dim, typename MatrixType, typename VectorType> void assemble_dg_matrix_on_standard_mesh(MatrixType &system_matrix, VectorType &system_rhs, const Mapping<dim> &mapping, const FiniteElement<dim> &fe_dg, const DoFHandler<dim> &dof_handler) { static_assert( (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix> || std::is_same_v<MatrixType, SparseMatrix<typename MatrixType::value_type>>)); Assert((dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_dg) != nullptr), ExcNotImplemented( "Implemented only for simplex meshes for the time being.")); Assert(dof_handler.get_triangulation().all_reference_cells_are_simplex(), ExcNotImplemented()); const double penalty_constant = .5 * fe_dg.degree * (fe_dg.degree + 1); AffineConstraints<typename MatrixType::value_type> constraints; constraints.close(); const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs(); const IndexSet locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler); DynamicSparsityPattern dsp(locally_relevant_dofs); DoFTools::make_flux_sparsity_pattern(dof_handler, dsp); SparsityTools::distribute_sparsity_pattern(dsp, dof_handler.locally_owned_dofs(), dof_handler.get_communicator(), locally_relevant_dofs); system_matrix.reinit(locally_owned_dofs, locally_owned_dofs, dsp, dof_handler.get_communicator()); system_rhs.reinit(locally_owned_dofs, dof_handler.get_communicator()); const unsigned int quadrature_degree = fe_dg.degree + 1; FEFaceValues<dim> fe_faces0(mapping, fe_dg, QGaussSimplex<dim - 1>(quadrature_degree), update_values | update_JxW_values | update_gradients | update_quadrature_points | update_normal_vectors); FEValues<dim> fe_values(mapping, fe_dg, QGaussSimplex<dim>(quadrature_degree), update_values | update_JxW_values | update_gradients | update_quadrature_points); FEFaceValues<dim> fe_faces1(mapping, fe_dg, QGaussSimplex<dim - 1>(quadrature_degree), update_values | update_JxW_values | update_gradients | update_quadrature_points | update_normal_vectors); const unsigned int dofs_per_cell = fe_dg.n_dofs_per_cell(); FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell); Vector<double> cell_rhs(dofs_per_cell); FullMatrix<double> M11(dofs_per_cell, dofs_per_cell); FullMatrix<double> M12(dofs_per_cell, dofs_per_cell); FullMatrix<double> M21(dofs_per_cell, dofs_per_cell); FullMatrix<double> M22(dofs_per_cell, dofs_per_cell); std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell); @endcode Loop over standard deal.II cells @code for (const auto &cell : dof_handler.active_cell_iterators()) { if (cell->is_locally_owned()) { cell_matrix = 0.; cell_rhs = 0.; fe_values.reinit(cell); @endcode const auto &q_points = fe_values.get_quadrature_points(); const unsigned int n_qpoints = q_points.size(); @code for (unsigned int q_index : fe_values.quadrature_point_indices()) { for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { cell_matrix(i, j) += fe_values.shape_grad(i, q_index) * fe_values.shape_grad(j, q_index) * fe_values.JxW(q_index); } cell_rhs(i) += fe_values.shape_value(i, q_index) * 1. * fe_values.JxW(q_index); // TODO: pass functional } } @endcode distribute volumetric DoFs @code cell->get_dof_indices(local_dof_indices); double hf = 0.; for (const auto f : cell->face_indices()) { const double extent1 = cell->measure() / cell->face(f)->measure(); if (cell->face(f)->at_boundary()) { hf = (1. / extent1 + 1. / extent1); fe_faces0.reinit(cell, f); const auto &normals = fe_faces0.get_normal_vectors(); for (unsigned int q_index : fe_faces0.quadrature_point_indices()) { for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { cell_matrix(i, j) += (-fe_faces0.shape_value(i, q_index) * fe_faces0.shape_grad(j, q_index) * normals[q_index] - fe_faces0.shape_grad(i, q_index) * normals[q_index] * fe_faces0.shape_value(j, q_index) + (penalty_constant * hf) * fe_faces0.shape_value(i, q_index) * fe_faces0.shape_value(j, q_index)) * fe_faces0.JxW(q_index); } cell_rhs(i) += 0.; // TODO: add bdary conditions functional } } } else { const auto &neigh_cell = cell->neighbor(f); if (cell->global_active_cell_index() < neigh_cell->global_active_cell_index()) { const double extent2 = neigh_cell->measure() / neigh_cell->face(cell->neighbor_of_neighbor(f)) ->measure(); hf = (1. / extent1 + 1. / extent2); fe_faces0.reinit(cell, f); fe_faces1.reinit(neigh_cell, cell->neighbor_of_neighbor(f)); std::vector<types::global_dof_index> local_dof_indices_neighbor(dofs_per_cell); M11 = 0.; M12 = 0.; M21 = 0.; M22 = 0.; const auto &normals = fe_faces0.get_normal_vectors(); @endcode M11 @code for (unsigned int q_index : fe_faces0.quadrature_point_indices()) { for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { M11(i, j) += (-0.5 * fe_faces0.shape_grad(i, q_index) * normals[q_index] * fe_faces0.shape_value(j, q_index) - 0.5 * fe_faces0.shape_grad(j, q_index) * normals[q_index] * fe_faces0.shape_value(i, q_index) + (penalty_constant * hf) * fe_faces0.shape_value(i, q_index) * fe_faces0.shape_value(j, q_index)) * fe_faces0.JxW(q_index); M12(i, j) += (0.5 * fe_faces0.shape_grad(i, q_index) * normals[q_index] * fe_faces1.shape_value(j, q_index) - 0.5 * fe_faces1.shape_grad(j, q_index) * normals[q_index] * fe_faces0.shape_value(i, q_index) - (penalty_constant * hf) * fe_faces0.shape_value(i, q_index) * fe_faces1.shape_value(j, q_index)) * fe_faces1.JxW(q_index); @endcode A10 @code M21(i, j) += (-0.5 * fe_faces1.shape_grad(i, q_index) * normals[q_index] * fe_faces0.shape_value(j, q_index) + 0.5 * fe_faces0.shape_grad(j, q_index) * normals[q_index] * fe_faces1.shape_value(i, q_index) - (penalty_constant * hf) * fe_faces1.shape_value(i, q_index) * fe_faces0.shape_value(j, q_index)) * fe_faces1.JxW(q_index); @endcode A11 @code M22(i, j) += (0.5 * fe_faces1.shape_grad(i, q_index) * normals[q_index] * fe_faces1.shape_value(j, q_index) + 0.5 * fe_faces1.shape_grad(j, q_index) * normals[q_index] * fe_faces1.shape_value(i, q_index) + (penalty_constant * hf) * fe_faces1.shape_value(i, q_index) * fe_faces1.shape_value(j, q_index)) * fe_faces1.JxW(q_index); } } } @endcode distribute DoFs accordingly @code neigh_cell->get_dof_indices(local_dof_indices_neighbor); constraints.distribute_local_to_global( M11, local_dof_indices, system_matrix); constraints.distribute_local_to_global( M12, local_dof_indices, local_dof_indices_neighbor, system_matrix); constraints.distribute_local_to_global( M21, local_dof_indices_neighbor, local_dof_indices, system_matrix); constraints.distribute_local_to_global( M22, local_dof_indices_neighbor, system_matrix); } // check idx neighbors } // over faces } constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); } } system_matrix.compress(VectorOperation::add); system_rhs.compress(VectorOperation::add); } } // namespace ::PolyUtils #endif @endcode <a name="ann-source/agglomeration_handler.cc"></a> <h1>Annotated version of source/agglomeration_handler.cc</h1> @code /* ----------------------------------------------------------------------------- * * SPDX-License-Identifier: LGPL-2.1-or-later * Copyright (C) 2025 by Marco Feder, Pasquale Claudio Africa, Xinping Gui, * Andrea Cangiani * * This file is part of the deal.II code gallery. * * ----------------------------------------------------------------------------- */ #include <deal.II/base/quadrature_lib.h> #include <deal.II/lac/sparsity_tools.h> #include <agglomeration_handler.h> template <int dim, int spacedim> AgglomerationHandler<dim, spacedim>::AgglomerationHandler( const GridTools::Cache<dim, spacedim> &cache_tria) : cached_tria(std::make_unique<GridTools::Cache<dim, spacedim>>( cache_tria.get_triangulation(), cache_tria.get_mapping())) , communicator(cache_tria.get_triangulation().get_mpi_communicator()) { Assert(dim == spacedim, ExcNotImplemented("Not available with codim > 0")); Assert(dim == 2 || dim == 3, ExcImpossibleInDim(1)); Assert((dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>( &cached_tria->get_triangulation()) == nullptr), ExcNotImplemented()); Assert(cached_tria->get_triangulation().n_active_cells() > 0, ExcMessage( "The triangulation must not be empty upon calling this function.")); n_agglomerations = 0; hybrid_mesh = false; initialize_agglomeration_data(cached_tria); } template <int dim, int spacedim> typename AgglomerationHandler<dim, spacedim>::agglomeration_iterator AgglomerationHandler<dim, spacedim>::define_agglomerate( const AgglomerationContainer &cells) { Assert(cells.size() > 0, ExcMessage("No cells to be agglomerated.")); if (cells.size() == 1) hybrid_mesh = true; // mesh is made also by classical cells @endcode First index drives the selection of the master cell. After that, store the master cell. @code const types::global_cell_index global_master_idx = cells[0]->global_active_cell_index(); const types::global_cell_index master_idx = cells[0]->active_cell_index(); master_cells_container.push_back(cells[0]); master_slave_relationships[global_master_idx] = -1; const typename DoFHandler<dim>::active_cell_iterator cell_dh = cells[0]->as_dof_handler_iterator(agglo_dh); cell_dh->set_active_fe_index(CellAgglomerationType::master); @endcode Store slave cells and save the relationship with the parent @code std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator> slaves; slaves.reserve(cells.size() - 1); @endcode exclude first cell since it's the master cell @code for (auto it = ++cells.begin(); it != cells.end(); ++it) { slaves.push_back(*it); master_slave_relationships[(*it)->global_active_cell_index()] = global_master_idx; // mark each slave master_slave_relationships_iterators[(*it)->active_cell_index()] = cells[0]; const typename DoFHandler<dim>::active_cell_iterator cell = (*it)->as_dof_handler_iterator(agglo_dh); cell->set_active_fe_index(CellAgglomerationType::slave); // slave cell @endcode If we have a p::d::T, check that all cells are in the same subdomain. If serial, just check that the subdomain_id is invalid. @code Assert(((*it)->subdomain_id() == tria->locally_owned_subdomain() || tria->locally_owned_subdomain() == numbers::invalid_subdomain_id), ExcInternalError()); } master_slave_relationships_iterators[master_idx] = cells[0]; // set iterator to master cell @endcode Store the slaves of each master @code master2slaves[master_idx] = slaves; @endcode Save to which polygon this agglomerate correspond @code master2polygon[master_idx] = n_agglomerations; ++n_agglomerations; // an agglomeration has been performed, record it create_bounding_box(cells); // fill the vector of bboxes @endcode Finally, return a polygonal iterator to the polytope just constructed. @code return {cells[0], this}; } template <int dim, int spacedim> typename AgglomerationHandler<dim, spacedim>::agglomeration_iterator AgglomerationHandler<dim, spacedim>::define_agglomerate( const AgglomerationContainer &cells, const unsigned int fecollection_size) { Assert(cells.size() > 0, ExcMessage("No cells to be agglomerated.")); if (cells.size() == 1) hybrid_mesh = true; // mesh is made also by classical cells @endcode First index drives the selection of the master cell. After that, store the master cell. @code const types::global_cell_index global_master_idx = cells[0]->global_active_cell_index(); const types::global_cell_index master_idx = cells[0]->active_cell_index(); master_cells_container.push_back(cells[0]); master_slave_relationships[global_master_idx] = -1; const typename DoFHandler<dim>::active_cell_iterator cell_dh = cells[0]->as_dof_handler_iterator(agglo_dh); cell_dh->set_active_fe_index(CellAgglomerationType::master); @endcode Store slave cells and save the relationship with the parent @code std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator> slaves; slaves.reserve(cells.size() - 1); @endcode exclude first cell since it's the master cell @code for (auto it = ++cells.begin(); it != cells.end(); ++it) { slaves.push_back(*it); master_slave_relationships[(*it)->global_active_cell_index()] = global_master_idx; // mark each slave master_slave_relationships_iterators[(*it)->active_cell_index()] = cells[0]; const typename DoFHandler<dim>::active_cell_iterator cell = (*it)->as_dof_handler_iterator(agglo_dh); cell->set_active_fe_index( fecollection_size); // slave cell (the last index) @endcode If we have a p::d::T, check that all cells are in the same subdomain. If serial, just check that the subdomain_id is invalid. @code Assert(((*it)->subdomain_id() == tria->locally_owned_subdomain() || tria->locally_owned_subdomain() == numbers::invalid_subdomain_id), ExcInternalError()); } master_slave_relationships_iterators[master_idx] = cells[0]; // set iterator to master cell @endcode Store the slaves of each master @code master2slaves[master_idx] = slaves; @endcode Save to which polygon this agglomerate correspond @code master2polygon[master_idx] = n_agglomerations; ++n_agglomerations; // an agglomeration has been performed, record it create_bounding_box(cells); // fill the vector of bboxes @endcode Finally, return a polygonal iterator to the polytope just constructed. @code return {cells[0], this}; } template <int dim, int spacedim> void AgglomerationHandler<dim, spacedim>::initialize_fe_values( const Quadrature<dim> &cell_quadrature, const UpdateFlags &flags, const Quadrature<dim - 1> &face_quadrature, const UpdateFlags &face_flags) { agglomeration_quad = cell_quadrature; agglomeration_flags = flags; agglomeration_face_quad = face_quadrature; agglomeration_face_flags = face_flags | internal_agglomeration_face_flags; no_values = std::make_unique<FEValues<dim>>(*mapping, dummy_fe, agglomeration_quad, update_quadrature_points | update_JxW_values); // only for quadrature no_face_values = std::make_unique<FEFaceValues<dim>>( *mapping, dummy_fe, agglomeration_face_quad, update_quadrature_points | update_JxW_values | update_normal_vectors); // only for quadrature } template <int dim, int spacedim> void AgglomerationHandler<dim, spacedim>::initialize_fe_values( const hp::QCollection<dim> &cell_qcollection, const UpdateFlags &flags, const hp::QCollection<dim - 1> &face_qcollection, const UpdateFlags &face_flags) { agglomeration_quad_collection = cell_qcollection; agglomeration_flags = flags; agglomeration_face_quad_collection = face_qcollection; agglomeration_face_flags = face_flags | internal_agglomeration_face_flags; mapping_collection = hp::MappingCollection<dim>(*mapping); dummy_fe_collection = hp::FECollection<dim, spacedim>(dummy_fe); hp_no_values = std::make_unique<hp::FEValues<dim>>( mapping_collection, dummy_fe_collection, agglomeration_quad_collection, update_quadrature_points | update_JxW_values); // only for quadrature hp_no_face_values = std::make_unique<hp::FEFaceValues<dim>>( mapping_collection, dummy_fe_collection, agglomeration_face_quad_collection, update_quadrature_points | update_JxW_values | update_normal_vectors); // only for quadrature } template <int dim, int spacedim> unsigned int AgglomerationHandler<dim, spacedim>::n_agglomerated_faces_per_cell( const typename Triangulation<dim, spacedim>::active_cell_iterator &cell) const { unsigned int n_neighbors = 0; for (const auto &f : cell->face_indices()) { const auto &neighboring_cell = cell->neighbor(f); if ((cell->face(f)->at_boundary()) || (neighboring_cell->is_active() && !are_cells_agglomerated(cell, neighboring_cell))) { ++n_neighbors; } } return n_neighbors; } template <int dim, int spacedim> void AgglomerationHandler<dim, spacedim>::initialize_agglomeration_data( const std::unique_ptr<GridTools::Cache<dim, spacedim>> &cache_tria) { tria = &(cache_tria->get_triangulation()); mapping = &(cache_tria->get_mapping()); agglo_dh.reinit(*tria); if (const auto parallel_tria = dynamic_cast< const ::parallel::TriangulationBase<dim, spacedim> *>(&*tria)) { const std::weak_ptr<const Utilities::MPI::Partitioner> cells_partitioner = parallel_tria->global_active_cell_index_partitioner(); master_slave_relationships.reinit( cells_partitioner.lock()->locally_owned_range(), communicator); } else { master_slave_relationships.reinit(tria->n_active_cells(), MPI_COMM_SELF); } polytope_cache.clear(); bboxes.clear(); @endcode First, update the pointer @code cached_tria = std::make_unique<GridTools::Cache<dim, spacedim>>( cache_tria->get_triangulation(), cache_tria->get_mapping()); connect_to_tria_signals(); n_agglomerations = 0; } template <int dim, int spacedim> void AgglomerationHandler<dim, spacedim>::distribute_agglomerated_dofs( const FiniteElement<dim> &fe_space) { if (dynamic_cast<const FE_DGQ<dim> *>(&fe_space)) fe = std::make_unique<FE_DGQ<dim>>(fe_space.degree); else if (dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_space)) fe = std::make_unique<FE_SimplexDGP<dim>>(fe_space.degree); else AssertThrow( false, ExcNotImplemented( "Currently, this interface supports only DGQ and DGP bases.")); box_mapping = std::make_unique<MappingBox<dim>>( bboxes, master2polygon); // construct bounding box mapping if (hybrid_mesh) { @endcode the mesh is composed by standard and agglomerate cells. initialize classes needed for standard cells in order to treat that finite element space as defined on a standard shape and not on the BoundingBox. @code standard_scratch = std::make_unique<ScratchData>(*mapping, *fe, QGauss<dim>(2 * fe_space.degree + 2), internal_agglomeration_flags); } fe_collection.push_back(*fe); // master fe_collection.push_back( FE_Nothing<dim, spacedim>(fe->reference_cell())); // slave initialize_hp_structure(); @endcode in case the tria is distributed, communicate ghost information with neighboring ranks @code const bool needs_ghost_info = dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&*tria) != nullptr; if (needs_ghost_info) setup_ghost_polytopes(); setup_connectivity_of_agglomeration(); if (needs_ghost_info) exchange_interface_values(); } template <int dim, int spacedim> void AgglomerationHandler<dim, spacedim>::distribute_agglomerated_dofs( const hp::FECollection<dim, spacedim> &fe_collection_in) { is_hp_collection = true; hp_fe_collection = std::make_unique<hp::FECollection<dim, spacedim>>( fe_collection_in); // copy the input collection box_mapping = std::make_unique<MappingBox<dim>>( bboxes, master2polygon); // construct bounding box mapping if (hybrid_mesh) { AssertThrow(false, ExcNotImplemented( "Hybrid mesh is not implemented for hp::FECollection.")); } for (unsigned int i = 0; i < fe_collection_in.size(); ++i) { if (dynamic_cast<const FESystem<dim> *>(&fe_collection_in[i])) { @endcode System case @code for (unsigned int b = 0; b < fe_collection_in[i].n_base_elements(); ++b) { if (!(dynamic_cast<const FE_DGQ<dim> *>( &fe_collection_in[i].base_element(b)) || dynamic_cast<const FE_SimplexDGP<dim> *>( &fe_collection_in[i].base_element(b)) || dynamic_cast<const FE_Nothing<dim> *>( &fe_collection_in[i].base_element(b)))) AssertThrow( false, ExcNotImplemented( "Currently, this interface supports only DGQ and DGP bases.")); } } else { @endcode Scalar case @code if (!(dynamic_cast<const FE_DGQ<dim> *>(&fe_collection_in[i]) || dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_collection_in[i]))) AssertThrow( false, ExcNotImplemented( "Currently, this interface supports only DGQ and DGP bases.")); } fe_collection.push_back(fe_collection_in[i]); } Assert(fe_collection[0].n_components() >= 1, ExcMessage("Invalid FE: must have at least one component.")); if (fe_collection[0].n_components() == 1) { fe_collection.push_back(FE_Nothing<dim, spacedim>()); } else if (fe_collection[0].n_components() > 1) { std::vector<const FiniteElement<dim, spacedim> *> base_elements; std::vector<unsigned int> multiplicities; for (unsigned int b = 0; b < fe_collection[0].n_base_elements(); ++b) { base_elements.push_back(new FE_Nothing<dim, spacedim>()); multiplicities.push_back(fe_collection[0].element_multiplicity(b)); } FESystem<dim, spacedim> fe_system_nothing(base_elements, multiplicities); for (const auto *ptr : base_elements) delete ptr; fe_collection.push_back(fe_system_nothing); } initialize_hp_structure(); @endcode in case the tria is distributed, communicate ghost information with neighboring ranks @code const bool needs_ghost_info = dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&*tria) != nullptr; if (needs_ghost_info) setup_ghost_polytopes(); setup_connectivity_of_agglomeration(); if (needs_ghost_info) exchange_interface_values(); } template <int dim, int spacedim> void AgglomerationHandler<dim, spacedim>::create_bounding_box( const AgglomerationContainer &polytope) { Assert(n_agglomerations > 0, ExcMessage("No agglomeration has been performed.")); Assert(dim > 1, ExcNotImplemented()); std::vector<Point<spacedim>> pts; // store all the vertices for (const auto &cell : polytope) for (const auto i : cell->vertex_indices()) pts.push_back(cell->vertex(i)); bboxes.emplace_back(pts); } template <int dim, int spacedim> void AgglomerationHandler<dim, spacedim>::setup_connectivity_of_agglomeration() { Assert(master_cells_container.size() > 0, ExcMessage("No agglomeration has been performed.")); Assert( agglo_dh.n_dofs() > 0, ExcMessage( "The DoFHandler associated to the agglomeration has not been initialized." "It's likely that you forgot to distribute the DoFs. You may want" "to check if a call to `initialize_hp_structure()` has been done.")); number_of_agglomerated_faces.resize(master2polygon.size(), 0); for (const auto &cell : master_cells_container) { internal::AgglomerationHandlerImplementation<dim, spacedim>:: setup_master_neighbor_connectivity(cell, *this); } if (Utilities::MPI::job_supports_mpi()) { @endcode communicate the number of faces @code recv_n_faces = Utilities::MPI::some_to_some(communicator, local_n_faces); @endcode send information about boundaries and neighboring polytopes id @code recv_bdary_info = Utilities::MPI::some_to_some(communicator, local_bdary_info); recv_ghosted_master_id = Utilities::MPI::some_to_some(communicator, local_ghosted_master_id); } } template <int dim, int spacedim> void AgglomerationHandler<dim, spacedim>::exchange_interface_values() { const unsigned int dofs_per_cell = fe->dofs_per_cell; for (const auto &polytope : polytope_iterators()) { if (polytope->is_locally_owned()) { const unsigned int n_faces = polytope->n_faces(); for (unsigned int f = 0; f < n_faces; ++f) { if (!polytope->at_boundary(f)) { const auto &neigh_polytope = polytope->neighbor(f); if (!neigh_polytope->is_locally_owned()) { @endcode Neighboring polytope is ghosted. Compute shape functions at the interface @code const auto ¤t_fe = reinit(polytope, f); std::vector<Point<spacedim>> qpoints_to_send = current_fe.get_quadrature_points(); const std::vector<double> &jxws_to_send = current_fe.get_JxW_values(); const std::vector<Tensor<1, spacedim>> &normals_to_send = current_fe.get_normal_vectors(); const types::subdomain_id neigh_rank = neigh_polytope->subdomain_id(); std::pair<CellId, unsigned int> cell_and_face{ polytope->id(), f}; @endcode Prepare data to send @code local_qpoints[neigh_rank].emplace(cell_and_face, qpoints_to_send); local_jxws[neigh_rank].emplace(cell_and_face, jxws_to_send); local_normals[neigh_rank].emplace(cell_and_face, normals_to_send); const unsigned int n_qpoints = qpoints_to_send.size(); @endcode TODO: check <tt>agglomeration_flags</tt> before computing values and gradients. @code std::vector<std::vector<double>> values_per_qpoints( dofs_per_cell); std::vector<std::vector<Tensor<1, spacedim>>> gradients_per_qpoints(dofs_per_cell); for (unsigned int i = 0; i < dofs_per_cell; ++i) { values_per_qpoints[i].resize(n_qpoints); gradients_per_qpoints[i].resize(n_qpoints); for (unsigned int q = 0; q < n_qpoints; ++q) { values_per_qpoints[i][q] = current_fe.shape_value(i, q); gradients_per_qpoints[i][q] = current_fe.shape_grad(i, q); } } local_values[neigh_rank].emplace(cell_and_face, values_per_qpoints); local_gradients[neigh_rank].emplace( cell_and_face, gradients_per_qpoints); } } } } } @endcode Finally, exchange with neighboring ranks @code recv_qpoints = Utilities::MPI::some_to_some(communicator, local_qpoints); recv_jxws = Utilities::MPI::some_to_some(communicator, local_jxws); recv_normals = Utilities::MPI::some_to_some(communicator, local_normals); recv_values = Utilities::MPI::some_to_some(communicator, local_values); recv_gradients = Utilities::MPI::some_to_some(communicator, local_gradients); } template <int dim, int spacedim> Quadrature<dim> AgglomerationHandler<dim, spacedim>::agglomerated_quadrature( const typename AgglomerationHandler<dim, spacedim>::AgglomerationContainer &cells, const typename Triangulation<dim, spacedim>::active_cell_iterator &master_cell) const { Assert(is_master_cell(master_cell), ExcMessage("This must be a master cell.")); std::vector<Point<dim>> vec_pts; std::vector<double> vec_JxWs; if (!is_hp_collection) { @endcode Original version: handle case without hp::FECollection @code for (const auto &dummy_cell : cells) { no_values->reinit(dummy_cell); auto q_points = no_values->get_quadrature_points(); // real qpoints const auto &JxWs = no_values->get_JxW_values(); std::transform(q_points.begin(), q_points.end(), std::back_inserter(vec_pts), [&](const Point<spacedim> &p) { return p; }); std::transform(JxWs.begin(), JxWs.end(), std::back_inserter(vec_JxWs), [&](const double w) { return w; }); } } else { @endcode Handle the hp::FECollection case @code const auto &master_cell_as_dh_iterator = master_cell->as_dof_handler_iterator(agglo_dh); for (const auto &dummy_cell : cells) { @endcode The following verbose call is necessary to handle cases where different slave cells on different polytopes use different quadrature rules. If the hp::QCollection contains multiple elements, calling hp_no_values->reinit(dummy_cell) won't work because it cannot infer the correct quadrature rule. By explicitly passing the active FE index as q_index, and setting mapping_index and fe_index to 0, we ensure that the dummy cell uses the same quadrature rule as its corresponding master cell. This assumes a one-to-one correspondence between hp::QCollection and hp::FECollection, which is the convention in deal.II. However, this implementation does not support cases where hp::QCollection and hp::FECollection have different sizes. TODO: Refactor the architecture to better handle numerical integration for hp::QCollection. @code hp_no_values->reinit(dummy_cell, master_cell_as_dh_iterator->active_fe_index(), 0, 0); auto q_points = hp_no_values->get_present_fe_values() .get_quadrature_points(); // real qpoints const auto &JxWs = hp_no_values->get_present_fe_values().get_JxW_values(); std::transform(q_points.begin(), q_points.end(), std::back_inserter(vec_pts), [&](const Point<spacedim> &p) { return p; }); std::transform(JxWs.begin(), JxWs.end(), std::back_inserter(vec_JxWs), [&](const double w) { return w; }); } } @endcode Map back each point in real space by using the map associated to the bounding box. @code std::vector<Point<dim>> unit_points(vec_pts.size()); const auto &bbox = bboxes[master2polygon.at(master_cell->active_cell_index())]; unit_points.reserve(vec_pts.size()); for (unsigned int i = 0; i < vec_pts.size(); i++) unit_points[i] = bbox.real_to_unit(vec_pts[i]); return Quadrature<dim>(unit_points, vec_JxWs); } template <int dim, int spacedim> void AgglomerationHandler<dim, spacedim>::initialize_hp_structure() { Assert(agglo_dh.get_triangulation().n_cells() > 0, ExcMessage( "Triangulation must not be empty upon calling this function.")); Assert(n_agglomerations > 0, ExcMessage("No agglomeration has been performed.")); agglo_dh.distribute_dofs(fe_collection); @endcode euler_mapping = std::make_unique< MappingFEField<dim, spacedim, LinearAlgebra::distributed::Vector<double>>>( euler_dh, euler_vector); @code } template <int dim, int spacedim> const FEValues<dim, spacedim> & AgglomerationHandler<dim, spacedim>::reinit( const AgglomerationIterator<dim, spacedim> &polytope) const { @endcode Assert(euler_mapping, ExcMessage("The mapping describing the physical element stemming from " "agglomeration has not been set up.")); @code const auto &deal_cell = polytope->as_dof_handler_iterator(agglo_dh); @endcode First check if the polytope is made just by a single cell. If so, use classical FEValues if (polytope->n_background_cells() == 1) return standard_scratch->reinit(deal_cell); @code const auto &agglo_cells = polytope->get_agglomerate(); Quadrature<dim> agglo_quad = agglomerated_quadrature(agglo_cells, deal_cell); if (!is_hp_collection) { @endcode Original version: handle case without hp::FECollection @code agglomerated_scratch = std::make_unique<ScratchData>(*box_mapping, fe_collection[0], agglo_quad, agglomeration_flags); } else { @endcode Handle the hp::FECollection case @code agglomerated_scratch = std::make_unique<ScratchData>(*box_mapping, polytope->get_fe(), agglo_quad, agglomeration_flags); } return agglomerated_scratch->reinit(deal_cell); } template <int dim, int spacedim> const FEValuesBase<dim, spacedim> & AgglomerationHandler<dim, spacedim>::reinit_master( const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell, const unsigned int face_index, std::unique_ptr<NonMatching::FEImmersedSurfaceValues<spacedim>> &agglo_isv_ptr) const { return internal::AgglomerationHandlerImplementation<dim, spacedim>:: reinit_master(cell, face_index, agglo_isv_ptr, *this); } template <int dim, int spacedim> const FEValuesBase<dim, spacedim> & AgglomerationHandler<dim, spacedim>::reinit( const AgglomerationIterator<dim, spacedim> &polytope, const unsigned int face_index) const { @endcode Assert(euler_mapping, ExcMessage("The mapping describing the physical element stemming from " "agglomeration has not been set up.")); @code const auto &deal_cell = polytope->as_dof_handler_iterator(agglo_dh); Assert(is_master_cell(deal_cell), ExcMessage("This should be true.")); return internal::AgglomerationHandlerImplementation<dim, spacedim>:: reinit_master(deal_cell, face_index, agglomerated_isv_bdary, *this); } template <int dim, int spacedim> std::pair<const FEValuesBase<dim, spacedim> &, const FEValuesBase<dim, spacedim> &> AgglomerationHandler<dim, spacedim>::reinit_interface( const AgglomerationIterator<dim, spacedim> &polytope_in, const AgglomerationIterator<dim, spacedim> &neigh_polytope, const unsigned int local_in, const unsigned int local_neigh) const { @endcode If current and neighboring polytopes are both locally owned, then compute the jump in the classical way without needing information about ghosted entities. @code if (polytope_in->is_locally_owned() && neigh_polytope->is_locally_owned()) { const auto &cell_in = polytope_in->as_dof_handler_iterator(agglo_dh); const auto &neigh_cell = neigh_polytope->as_dof_handler_iterator(agglo_dh); const auto &fe_in = internal::AgglomerationHandlerImplementation<dim, spacedim>:: reinit_master(cell_in, local_in, agglomerated_isv, *this); const auto &fe_out = internal::AgglomerationHandlerImplementation<dim, spacedim>:: reinit_master(neigh_cell, local_neigh, agglomerated_isv_neigh, *this); std::pair<const FEValuesBase<dim, spacedim> &, const FEValuesBase<dim, spacedim> &> my_p(fe_in, fe_out); return my_p; } else { Assert((polytope_in->is_locally_owned() && !neigh_polytope->is_locally_owned()), ExcInternalError()); const auto &cell = polytope_in->as_dof_handler_iterator(agglo_dh); const auto &bbox = bboxes[master2polygon.at(cell->active_cell_index())]; @endcode const double bbox_measure = bbox.volume(); @code const unsigned int neigh_rank = neigh_polytope->subdomain_id(); const CellId &neigh_id = neigh_polytope->id(); @endcode Retrieve qpoints,JxWs, normals sent previously from the neighboring rank. @code std::vector<Point<spacedim>> &real_qpoints = recv_qpoints.at(neigh_rank).at({neigh_id, local_neigh}); const auto &JxWs = recv_jxws.at(neigh_rank).at({neigh_id, local_neigh}); std::vector<Tensor<1, spacedim>> &normals = recv_normals.at(neigh_rank).at({neigh_id, local_neigh}); @endcode Apply the necessary scalings due to the bbox. @code std::vector<Point<spacedim>> final_unit_q_points; std::transform(real_qpoints.begin(), real_qpoints.end(), std::back_inserter(final_unit_q_points), [&](const Point<spacedim> &p) { return bbox.real_to_unit(p); }); @endcode std::vector<double> scale_factors(final_unit_q_points.size()); std::vector<double> scaled_weights(final_unit_q_points.size()); std::vector<Tensor<1, dim>> scaled_normals(final_unit_q_points.size()); Since we received normal vectors from a neighbor, we have to swap the // sign of the vector in order to have outward normals. for (unsigned int q = 0; q < final_unit_q_points.size(); ++q) { for (unsigned int direction = 0; direction < spacedim; ++direction) scaled_normals[q][direction] = normals[q][direction] * (bbox.side_length(direction)); scaled_normals[q] *= -1; scaled_weights[q] = (JxWs[q] * scaled_normals[q].norm()) / bbox_measure; scaled_normals[q] /= scaled_normals[q].norm(); } @code for (unsigned int q = 0; q < final_unit_q_points.size(); ++q) normals[q] *= -1; NonMatching::ImmersedSurfaceQuadrature<dim, spacedim> surface_quad( final_unit_q_points, JxWs, normals); agglomerated_isv = std::make_unique<NonMatching::FEImmersedSurfaceValues<spacedim>>( *box_mapping, *fe, surface_quad, agglomeration_face_flags); agglomerated_isv->reinit(cell); std::pair<const FEValuesBase<dim, spacedim> &, const FEValuesBase<dim, spacedim> &> my_p(*agglomerated_isv, *agglomerated_isv); return my_p; } } template <int dim, int spacedim> template <typename SparsityPatternType, typename Number> void AgglomerationHandler<dim, spacedim>::create_agglomeration_sparsity_pattern( SparsityPatternType &dsp, const AffineConstraints<Number> &constraints, const bool keep_constrained_dofs, const types::subdomain_id subdomain_id) { Assert(n_agglomerations > 0, ExcMessage("The agglomeration has not been set up correctly.")); Assert(dsp.empty(), ExcMessage( "The Sparsity pattern must be empty upon calling this function.")); const IndexSet &locally_owned_dofs = agglo_dh.locally_owned_dofs(); const IndexSet locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(agglo_dh); if constexpr (std::is_same_v<SparsityPatternType, DynamicSparsityPattern>) dsp.reinit(locally_owned_dofs.size(), locally_owned_dofs.size(), locally_relevant_dofs); else if constexpr (std::is_same_v<SparsityPatternType, TrilinosWrappers::SparsityPattern>) dsp.reinit(locally_owned_dofs, communicator); else AssertThrow(false, ExcNotImplemented()); @endcode Create the sparsity pattern corresponding only to volumetric terms. The fluxes needed by DG methods will be filled later. @code DoFTools::make_sparsity_pattern( agglo_dh, dsp, constraints, keep_constrained_dofs, subdomain_id); if (!is_hp_collection) { @endcode Original version: handle case without hp::FECollection @code const unsigned int dofs_per_cell = agglo_dh.get_fe(0).n_dofs_per_cell(); std::vector<types::global_dof_index> current_dof_indices(dofs_per_cell); std::vector<types::global_dof_index> neighbor_dof_indices(dofs_per_cell); @endcode Loop over all locally owned polytopes, find the neighbor (also ghosted) and add fluxes to the sparsity pattern. @code for (const auto &polytope : polytope_iterators()) { if (polytope->is_locally_owned()) { const unsigned int n_current_faces = polytope->n_faces(); polytope->get_dof_indices(current_dof_indices); for (unsigned int f = 0; f < n_current_faces; ++f) { const auto &neigh_polytope = polytope->neighbor(f); if (neigh_polytope.state() == IteratorState::valid) { neigh_polytope->get_dof_indices(neighbor_dof_indices); constraints.add_entries_local_to_global( current_dof_indices, neighbor_dof_indices, dsp, keep_constrained_dofs, {}); } } } } } else { @endcode Handle the hp::FECollection case Loop over all locally owned polytopes, find the neighbor (also ghosted) and add fluxes to the sparsity pattern. @code for (const auto &polytope : polytope_iterators()) { if (polytope->is_locally_owned()) { const unsigned int current_dofs_per_cell = polytope->get_fe().dofs_per_cell; std::vector<types::global_dof_index> current_dof_indices( current_dofs_per_cell); const unsigned int n_current_faces = polytope->n_faces(); polytope->get_dof_indices(current_dof_indices); for (unsigned int f = 0; f < n_current_faces; ++f) { const auto &neigh_polytope = polytope->neighbor(f); if (neigh_polytope.state() == IteratorState::valid) { const unsigned int neighbor_dofs_per_cell = neigh_polytope->get_fe().dofs_per_cell; std::vector<types::global_dof_index> neighbor_dof_indices( neighbor_dofs_per_cell); neigh_polytope->get_dof_indices(neighbor_dof_indices); constraints.add_entries_local_to_global( current_dof_indices, neighbor_dof_indices, dsp, keep_constrained_dofs, {}); } } } } } if constexpr (std::is_same_v<SparsityPatternType, TrilinosWrappers::SparsityPattern>) dsp.compress(); } template <int dim, int spacedim> void AgglomerationHandler<dim, spacedim>::setup_ghost_polytopes() { [[maybe_unused]] const auto parallel_triangulation = dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&*tria); Assert(parallel_triangulation != nullptr, ExcInternalError()); const unsigned int n_dofs_per_cell = fe->dofs_per_cell; std::vector<types::global_dof_index> global_dof_indices(n_dofs_per_cell); for (const auto &polytope : polytope_iterators()) if (polytope->is_locally_owned()) { const CellId &master_cell_id = polytope->id(); const auto polytope_dh = polytope->as_dof_handler_iterator(agglo_dh); polytope_dh->get_dof_indices(global_dof_indices); const auto &agglomerate = polytope->get_agglomerate(); for (const auto &cell : agglomerate) { @endcode interior, locally owned, cell @code for (const auto &f : cell->face_indices()) { if (!cell->at_boundary(f)) { const auto &neighbor = cell->neighbor(f); if (neighbor->is_ghost()) { @endcode key of the map: the rank to which send the data @code const types::subdomain_id neigh_rank = neighbor->subdomain_id(); @endcode inform the "standard" neighbor about the neighboring id and its master cell @code local_cell_ids_neigh_cell[neigh_rank].emplace( cell->id(), master_cell_id); @endcode inform the neighboring rank that this master cell (hence polytope) has the following DoF indices @code local_ghost_dofs[neigh_rank].emplace( master_cell_id, global_dof_indices); @endcode ...same for bounding boxes @code const auto &bbox = bboxes[polytope->index()]; local_ghosted_bbox[neigh_rank].emplace(master_cell_id, bbox); } } } } } recv_cell_ids_neigh_cell = Utilities::MPI::some_to_some(communicator, local_cell_ids_neigh_cell); @endcode Exchange with neighboring ranks the neighboring bounding boxes @code recv_ghosted_bbox = Utilities::MPI::some_to_some(communicator, local_ghosted_bbox); @endcode Exchange with neighboring ranks the neighboring ghosted DoFs @code recv_ghost_dofs = Utilities::MPI::some_to_some(communicator, local_ghost_dofs); } namespace dealii { namespace internal { template <int dim, int spacedim> class AgglomerationHandlerImplementation { public: static const FEValuesBase<dim, spacedim> & reinit_master( const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell, const unsigned int face_index, std::unique_ptr<NonMatching::FEImmersedSurfaceValues<spacedim>> &agglo_isv_ptr, const AgglomerationHandler<dim, spacedim> &handler) { Assert(handler.is_master_cell(cell), ExcMessage("This cell must be a master one.")); AgglomerationIterator<dim, spacedim> it{cell, &handler}; const auto &neigh_polytope = it->neighbor(face_index); const CellId polytope_in_id = cell->id(); @endcode Retrieve the bounding box of the agglomeration @code const auto &bbox = handler.bboxes[handler.master2polygon.at(cell->active_cell_index())]; CellId polytope_out_id; if (neigh_polytope.state() == IteratorState::valid) polytope_out_id = neigh_polytope->id(); else polytope_out_id = polytope_in_id; // on the boundary. Same id const auto &common_face = handler.polytope_cache.interface.at( {polytope_in_id, polytope_out_id}); std::vector<Point<spacedim>> final_unit_q_points; std::vector<double> final_weights; std::vector<Tensor<1, dim>> final_normals; if (!handler.is_hp_collection) { @endcode Original version: handle case without hp::FECollection @code const unsigned int expected_qpoints = common_face.size() * handler.agglomeration_face_quad.size(); final_unit_q_points.reserve(expected_qpoints); final_weights.reserve(expected_qpoints); final_normals.reserve(expected_qpoints); for (const auto &[deal_cell, local_face_idx] : common_face) { handler.no_face_values->reinit(deal_cell, local_face_idx); const auto &q_points = handler.no_face_values->get_quadrature_points(); const auto &JxWs = handler.no_face_values->get_JxW_values(); const auto &normals = handler.no_face_values->get_normal_vectors(); const unsigned int n_qpoints_agglo = q_points.size(); for (unsigned int q = 0; q < n_qpoints_agglo; ++q) { final_unit_q_points.push_back( bbox.real_to_unit(q_points[q])); final_weights.push_back(JxWs[q]); final_normals.push_back(normals[q]); } } } else { @endcode Handle the hp::FECollection case @code unsigned int higher_order_quad_index = cell->active_fe_index(); if (neigh_polytope.state() == IteratorState::valid) if (handler .agglomeration_face_quad_collection[cell->active_fe_index()] .size() < handler .agglomeration_face_quad_collection[neigh_polytope ->active_fe_index()] .size()) higher_order_quad_index = neigh_polytope->active_fe_index(); const unsigned int expected_qpoints = common_face.size() * handler .agglomeration_face_quad_collection[higher_order_quad_index] .size(); final_unit_q_points.reserve(expected_qpoints); final_weights.reserve(expected_qpoints); final_normals.reserve(expected_qpoints); for (const auto &[deal_cell, local_face_idx] : common_face) { handler.hp_no_face_values->reinit( deal_cell, local_face_idx, higher_order_quad_index, 0, 0); const auto &q_points = handler.hp_no_face_values->get_present_fe_values() .get_quadrature_points(); const auto &JxWs = handler.hp_no_face_values->get_present_fe_values() .get_JxW_values(); const auto &normals = handler.hp_no_face_values->get_present_fe_values() .get_normal_vectors(); const unsigned int n_qpoints_agglo = q_points.size(); for (unsigned int q = 0; q < n_qpoints_agglo; ++q) { final_unit_q_points.push_back( bbox.real_to_unit(q_points[q])); final_weights.push_back(JxWs[q]); final_normals.push_back(normals[q]); } } } NonMatching::ImmersedSurfaceQuadrature<dim, spacedim> surface_quad( final_unit_q_points, final_weights, final_normals); if (!handler.is_hp_collection) { agglo_isv_ptr = std::make_unique<NonMatching::FEImmersedSurfaceValues<spacedim>>( *(handler.box_mapping), *(handler.fe), surface_quad, handler.agglomeration_face_flags); } else { agglo_isv_ptr = std::make_unique<NonMatching::FEImmersedSurfaceValues<spacedim>>( *(handler.box_mapping), cell->get_fe(), surface_quad, handler.agglomeration_face_flags); } agglo_isv_ptr->reinit(cell); return *agglo_isv_ptr; } /** * Given an agglomeration described by the master cell `master_cell`, * this function: * - enumerates the faces of the agglomeration * - stores who is the neighbor, the local face indices from outside and * inside*/ static void setup_master_neighbor_connectivity( const typename Triangulation<dim, spacedim>::active_cell_iterator &master_cell, const AgglomerationHandler<dim, spacedim> &handler) { Assert( handler.master_slave_relationships[master_cell ->global_active_cell_index()] == -1, ExcMessage("The present cell with index " + std::to_string(master_cell->global_active_cell_index()) + "is not a master one.")); const auto &agglomeration = handler.get_agglomerate(master_cell); const types::global_cell_index current_polytope_index = handler.master2polygon.at(master_cell->active_cell_index()); CellId current_polytope_id = master_cell->id(); std::set<types::global_cell_index> visited_polygonal_neighbors; std::map<unsigned int, CellId> face_to_neigh_id; std::map<unsigned int, bool> is_face_at_boundary; @endcode same as above, but with CellId @code std::set<CellId> visited_polygonal_neighbors_id; unsigned int ghost_counter = 0; for (const auto &cell : agglomeration) { const types::global_cell_index cell_index = cell->active_cell_index(); const CellId cell_id = cell->id(); for (const auto f : cell->face_indices()) { const auto &neighboring_cell = cell->neighbor(f); const bool valid_neighbor = neighboring_cell.state() == IteratorState::valid; if (valid_neighbor) { if (neighboring_cell->is_locally_owned() && !handler.are_cells_agglomerated(cell, neighboring_cell)) { @endcode - cell is not on the boundary, - it's not agglomerated with the neighbor. If so, it's a neighbor of the present agglomeration std::cout << " (from rank) " << Utilities::MPI::this_mpi_process( handler.communicator) << std::endl; std::cout << "neighbor locally owned? " << std::boolalpha << neighboring_cell->is_locally_owned() << std::endl; if (neighboring_cell->is_ghost()) handler.ghosted_indices.push_back( neighboring_cell->active_cell_index()); a new face of the agglomeration has been discovered. @code handler.polygon_boundary[master_cell].push_back( cell->face(f)); @endcode global index of neighboring deal.II cell @code const types::global_cell_index neighboring_cell_index = neighboring_cell->active_cell_index(); @endcode master cell for the neighboring polytope @code const auto &master_of_neighbor = handler.master_slave_relationships_iterators.at( neighboring_cell_index); const auto nof = cell->neighbor_of_neighbor(f); if (handler.is_slave_cell(neighboring_cell)) { @endcode index of the neighboring polytope @code const types::global_cell_index neighbor_polytope_index = handler.master2polygon.at( master_of_neighbor->active_cell_index()); CellId neighbor_polytope_id = master_of_neighbor->id(); if (visited_polygonal_neighbors.find( neighbor_polytope_index) == std::end(visited_polygonal_neighbors)) { @endcode found a neighbor @code const unsigned int n_face = handler.number_of_agglomerated_faces [current_polytope_index]; handler.polytope_cache.cell_face_at_boundary[{ current_polytope_index, n_face}] = { false, master_of_neighbor}; is_face_at_boundary[n_face] = true; ++handler.number_of_agglomerated_faces [current_polytope_index]; visited_polygonal_neighbors.insert( neighbor_polytope_index); } if (handler.polytope_cache.visited_cell_and_faces .find({cell_index, f}) == std::end(handler.polytope_cache .visited_cell_and_faces)) { handler.polytope_cache .interface[{current_polytope_id, neighbor_polytope_id}] .emplace_back(cell, f); handler.polytope_cache.visited_cell_and_faces .insert({cell_index, f}); } if (handler.polytope_cache.visited_cell_and_faces .find({neighboring_cell_index, nof}) == std::end(handler.polytope_cache .visited_cell_and_faces)) { handler.polytope_cache .interface[{neighbor_polytope_id, current_polytope_id}] .emplace_back(neighboring_cell, nof); handler.polytope_cache.visited_cell_and_faces .insert({neighboring_cell_index, nof}); } } else { @endcode neighboring cell is a master save the pair of neighboring cells @code const types::global_cell_index neighbor_polytope_index = handler.master2polygon.at( neighboring_cell_index); CellId neighbor_polytope_id = neighboring_cell->id(); if (visited_polygonal_neighbors.find( neighbor_polytope_index) == std::end(visited_polygonal_neighbors)) { @endcode found a neighbor @code const unsigned int n_face = handler.number_of_agglomerated_faces [current_polytope_index]; handler.polytope_cache.cell_face_at_boundary[{ current_polytope_index, n_face}] = { false, neighboring_cell}; is_face_at_boundary[n_face] = true; ++handler.number_of_agglomerated_faces [current_polytope_index]; visited_polygonal_neighbors.insert( neighbor_polytope_index); } if (handler.polytope_cache.visited_cell_and_faces .find({cell_index, f}) == std::end(handler.polytope_cache .visited_cell_and_faces)) { handler.polytope_cache .interface[{current_polytope_id, neighbor_polytope_id}] .emplace_back(cell, f); handler.polytope_cache.visited_cell_and_faces .insert({cell_index, f}); } if (handler.polytope_cache.visited_cell_and_faces .find({neighboring_cell_index, nof}) == std::end(handler.polytope_cache .visited_cell_and_faces)) { handler.polytope_cache .interface[{neighbor_polytope_id, current_polytope_id}] .emplace_back(neighboring_cell, nof); handler.polytope_cache.visited_cell_and_faces .insert({neighboring_cell_index, nof}); } } } else if (neighboring_cell->is_ghost()) { const auto nof = cell->neighbor_of_neighbor(f); @endcode from neighboring rank,receive the association between standard cell ids and neighboring polytope. This tells to the current rank that the neighboring cell has the following CellId as master cell. @code const auto &check_neigh_poly_ids = handler.recv_cell_ids_neigh_cell.at( neighboring_cell->subdomain_id()); const CellId neighboring_cell_id = neighboring_cell->id(); const CellId &check_neigh_polytope_id = check_neigh_poly_ids.at(neighboring_cell_id); @endcode const auto master_index = master_indices[ghost_counter]; @code if (visited_polygonal_neighbors_id.find( check_neigh_polytope_id) == std::end(visited_polygonal_neighbors_id)) { handler.polytope_cache.cell_face_at_boundary[{ current_polytope_index, handler.number_of_agglomerated_faces [current_polytope_index]}] = {false, neighboring_cell}; @endcode record the cell id of the neighboring polytope @code handler.polytope_cache.ghosted_master_id[{ current_polytope_id, handler.number_of_agglomerated_faces [current_polytope_index]}] = check_neigh_polytope_id; const unsigned int n_face = handler.number_of_agglomerated_faces [current_polytope_index]; face_to_neigh_id[n_face] = check_neigh_polytope_id; is_face_at_boundary[n_face] = false; @endcode increment number of faces @code ++handler.number_of_agglomerated_faces [current_polytope_index]; visited_polygonal_neighbors_id.insert( check_neigh_polytope_id); @endcode ghosted polytope has been found, increment ghost counter @code ++ghost_counter; } if (handler.polytope_cache.visited_cell_and_faces_id .find({cell_id, f}) == std::end( handler.polytope_cache.visited_cell_and_faces_id)) { handler.polytope_cache .interface[{current_polytope_id, check_neigh_polytope_id}] .emplace_back(cell, f); @endcode std::cout << "ADDED (" << cell->active_cell_index() << ") BETWEEN " << current_polytope_id << " e " << check_neigh_polytope_id << std::endl; @code handler.polytope_cache.visited_cell_and_faces_id .insert({cell_id, f}); } if (handler.polytope_cache.visited_cell_and_faces_id .find({neighboring_cell_id, nof}) == std::end( handler.polytope_cache.visited_cell_and_faces_id)) { handler.polytope_cache .interface[{check_neigh_polytope_id, current_polytope_id}] .emplace_back(neighboring_cell, nof); handler.polytope_cache.visited_cell_and_faces_id .insert({neighboring_cell_id, nof}); } } } else if (cell->face(f)->at_boundary()) { @endcode Boundary face of a boundary cell. Note that the neighboring cell must be invalid. @code handler.polygon_boundary[master_cell].push_back( cell->face(f)); if (visited_polygonal_neighbors.find( std::numeric_limits<unsigned int>::max()) == std::end(visited_polygonal_neighbors)) { @endcode boundary face. Notice that <tt>neighboring_cell</tt> is invalid here. @code handler.polytope_cache.cell_face_at_boundary[{ current_polytope_index, handler.number_of_agglomerated_faces [current_polytope_index]}] = {true, neighboring_cell}; const unsigned int n_face = handler.number_of_agglomerated_faces [current_polytope_index]; is_face_at_boundary[n_face] = true; ++handler.number_of_agglomerated_faces [current_polytope_index]; visited_polygonal_neighbors.insert( std::numeric_limits<unsigned int>::max()); } if (handler.polytope_cache.visited_cell_and_faces.find( {cell_index, f}) == std::end(handler.polytope_cache.visited_cell_and_faces)) { handler.polytope_cache .interface[{current_polytope_id, current_polytope_id}] .emplace_back(cell, f); handler.polytope_cache.visited_cell_and_faces.insert( {cell_index, f}); } } } // loop over faces } // loop over all cells of agglomerate if (ghost_counter > 0) { const auto parallel_triangulation = dynamic_cast< const ::parallel::TriangulationBase<dim, spacedim> *>( &(*handler.tria)); const unsigned int n_faces_current_poly = handler.number_of_agglomerated_faces[current_polytope_index]; @endcode Communicate to neighboring ranks that current_polytope_id has a number of faces equal to n_faces_current_poly faces: current_polytope_id -> n_faces_current_poly @code for (const unsigned int neigh_rank : parallel_triangulation->ghost_owners()) { handler.local_n_faces[neigh_rank].emplace(current_polytope_id, n_faces_current_poly); handler.local_bdary_info[neigh_rank].emplace( current_polytope_id, is_face_at_boundary); handler.local_ghosted_master_id[neigh_rank].emplace( current_polytope_id, face_to_neigh_id); } } } }; } // namespace internal } // namespace dealii template class AgglomerationHandler<1>; template void AgglomerationHandler<1>::create_agglomeration_sparsity_pattern( DynamicSparsityPattern &sparsity_pattern, const AffineConstraints<double> &constraints, const bool keep_constrained_dofs, const types::subdomain_id subdomain_id); template void AgglomerationHandler<1>::create_agglomeration_sparsity_pattern( TrilinosWrappers::SparsityPattern &sparsity_pattern, const AffineConstraints<double> &constraints, const bool keep_constrained_dofs, const types::subdomain_id subdomain_id); template class AgglomerationHandler<2>; template void AgglomerationHandler<2>::create_agglomeration_sparsity_pattern( DynamicSparsityPattern &sparsity_pattern, const AffineConstraints<double> &constraints, const bool keep_constrained_dofs, const types::subdomain_id subdomain_id); template void AgglomerationHandler<2>::create_agglomeration_sparsity_pattern( TrilinosWrappers::SparsityPattern &sparsity_pattern, const AffineConstraints<double> &constraints, const bool keep_constrained_dofs, const types::subdomain_id subdomain_id); template class AgglomerationHandler<3>; template void AgglomerationHandler<3>::create_agglomeration_sparsity_pattern( DynamicSparsityPattern &sparsity_pattern, const AffineConstraints<double> &constraints, const bool keep_constrained_dofs, const types::subdomain_id subdomain_id); template void AgglomerationHandler<3>::create_agglomeration_sparsity_pattern( TrilinosWrappers::SparsityPattern &sparsity_pattern, const AffineConstraints<double> &constraints, const bool keep_constrained_dofs, const types::subdomain_id subdomain_id); @endcode <a name="ann-source/mapping_box.cc"></a> <h1>Annotated version of source/mapping_box.cc</h1> @code /* ----------------------------------------------------------------------------- * * SPDX-License-Identifier: LGPL-2.1-or-later * Copyright (C) 2025 by Marco Feder, Pasquale Claudio Africa, Xinping Gui, * Andrea Cangiani * * This file is part of the deal.II code gallery. * * ----------------------------------------------------------------------------- */ #include <deal.II/base/array_view.h> #include <deal.II/base/memory_consumption.h> #include <deal.II/base/qprojector.h> #include <deal.II/base/quadrature.h> #include <deal.II/base/signaling_nan.h> #include <deal.II/base/tensor.h> #include <deal.II/dofs/dof_accessor.h> #include <deal.II/fe/fe_values.h> #include <deal.II/grid/tria.h> #include <deal.II/grid/tria_iterator.h> #include <deal.II/lac/full_matrix.h> #include <mapping_box.h> #include <algorithm> DEAL_II_NAMESPACE_OPEN DeclExceptionMsg( ExcCellNotAssociatedWithBox, "You are using MappingBox, but the incoming element is not associated with a" "Bounding Box Cartesian."); /** * Return whether the incoming element has a BoundingBox associated to it. * Simplicial and quad-hex meshes are supported. */ template <typename CellType> bool has_box(const CellType &cell, const std::map<types::global_cell_index, types::global_cell_index> &translator) { Assert((cell->reference_cell().is_hyper_cube() || cell->reference_cell().is_simplex()), ExcNotImplemented()); Assert((translator.find(cell->active_cell_index()) != translator.cend()), ExcCellNotAssociatedWithBox()); return true; } template <int dim, int spacedim> MappingBox<dim, spacedim>::MappingBox( const std::vector<BoundingBox<dim>> &input_boxes, const std::map<types::global_cell_index, types::global_cell_index> &global_to_polytope) { Assert(input_boxes.size() > 0, ExcMessage("Invalid number of bounding boxes.")); @endcode copy boxes and map @code boxes.resize(input_boxes.size()); for (unsigned int i = 0; i < input_boxes.size(); ++i) boxes[i] = input_boxes[i]; polytope_translator = global_to_polytope; } template <int dim, int spacedim> MappingBox<dim, spacedim>::InternalData::InternalData(const Quadrature<dim> &q) : cell_extents(numbers::signaling_nan<Tensor<1, dim>>()) , traslation(numbers::signaling_nan<Tensor<1, dim>>()) , inverse_cell_extents(numbers::signaling_nan<Tensor<1, dim>>()) , volume_element(numbers::signaling_nan<double>()) , quadrature_points(q.get_points()) {} template <int dim, int spacedim> void MappingBox<dim, spacedim>::InternalData::reinit(const UpdateFlags update_flags, const Quadrature<dim> &) { @endcode store the flags in the internal data object so we can access them in fill_fe_*_values(). use the transitive hull of the required flags @code this->update_each = update_flags; } template <int dim, int spacedim> std::size_t MappingBox<dim, spacedim>::InternalData::memory_consumption() const { return (Mapping<dim, spacedim>::InternalDataBase::memory_consumption() + MemoryConsumption::memory_consumption(cell_extents) + MemoryConsumption::memory_consumption(traslation) + MemoryConsumption::memory_consumption(inverse_cell_extents) + MemoryConsumption::memory_consumption(volume_element)); } template <int dim, int spacedim> bool MappingBox<dim, spacedim>::preserves_vertex_locations() const { return true; } template <int dim, int spacedim> bool MappingBox<dim, spacedim>::is_compatible_with( #if DEAL_II_VERSION_GTE(9, 8, 0) const ReferenceCell<dim> &reference_cell #else const ReferenceCell &reference_cell #endif ) const { Assert(dim == reference_cell.get_dimension(), ExcMessage("The dimension of your mapping (" + Utilities::to_string(dim) + ") and the reference cell cell_type (" + Utilities::to_string(reference_cell.get_dimension()) + " ) do not agree.")); return reference_cell.is_hyper_cube() || reference_cell.is_simplex(); } template <int dim, int spacedim> UpdateFlags MappingBox<dim, spacedim>::requires_update_flags(const UpdateFlags in) const { @endcode this mapping is pretty simple in that it can basically compute every piece of information wanted by FEValues without requiring computing any other quantities. boundary forms are one exception since they can be computed from the normal vectors without much further ado @code UpdateFlags out = in; if (out & update_boundary_forms) out |= update_normal_vectors; return out; } template <int dim, int spacedim> std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> MappingBox<dim, spacedim>::get_data(const UpdateFlags update_flags, const Quadrature<dim> &q) const { std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr = std::make_unique<InternalData>(); data_ptr->reinit(requires_update_flags(update_flags), q); return data_ptr; } template <int dim, int spacedim> std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> MappingBox<dim, spacedim>::get_subface_data( const UpdateFlags update_flags, const Quadrature<dim - 1> &quadrature) const { (void)update_flags; (void)quadrature; DEAL_II_NOT_IMPLEMENTED(); return {}; } template <int dim, int spacedim> void MappingBox<dim, spacedim>::update_cell_extents( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const InternalData &data) const { @endcode Compute start point and sizes along axes. The vertices to be looked at are 1, 2, 4 compared to the base vertex 0. @code if (cell_similarity != CellSimilarity::translation) { const BoundingBox<dim> ¤t_box = boxes[polytope_translator.at(cell->active_cell_index())]; const std::pair<Point<dim>, Point<dim>> &bdary_points = current_box.get_boundary_points(); for (unsigned int d = 0; d < dim; ++d) { const double cell_extent_d = current_box.side_length(d); data.cell_extents[d] = cell_extent_d; data.traslation[d] = .5 * (bdary_points.first[d] + bdary_points.second[d]); // midpoint of each interval Assert(cell_extent_d != 0., ExcMessage("Cell does not appear to be Cartesian!")); data.inverse_cell_extents[d] = 1. / cell_extent_d; } } } namespace { template <int dim> void transform_quadrature_points( const BoundingBox<dim> &box, const ArrayView<const Point<dim>> &unit_quadrature_points, std::vector<Point<dim>> &quadrature_points) { for (unsigned int i = 0; i < quadrature_points.size(); ++i) quadrature_points[i] = box.unit_to_real(unit_quadrature_points[i]); } } // namespace template <int dim, int spacedim> void MappingBox<dim, spacedim>::maybe_update_cell_quadrature_points( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const InternalData &data, const ArrayView<const Point<dim>> &unit_quadrature_points, std::vector<Point<dim>> &quadrature_points) const { if (data.update_each & update_quadrature_points) transform_quadrature_points( boxes[polytope_translator.at(cell->active_cell_index())], unit_quadrature_points, quadrature_points); } template <int dim, int spacedim> void MappingBox<dim, spacedim>::maybe_update_normal_vectors( const unsigned int face_no, const InternalData &data, std::vector<Tensor<1, dim>> &normal_vectors) const { @endcode compute normal vectors. All normals on a face have the same value. @code if (data.update_each & update_normal_vectors) { Assert(face_no < GeometryInfo<dim>::faces_per_cell, ExcInternalError()); std::fill(normal_vectors.begin(), normal_vectors.end(), GeometryInfo<dim>::unit_normal_vector[face_no]); } } template <int dim, int spacedim> void MappingBox<dim, spacedim>::maybe_update_jacobian_derivatives( const InternalData &data, const CellSimilarity::Similarity cell_similarity, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const { if (cell_similarity != CellSimilarity::translation) { if (data.update_each & update_jacobian_grads) for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i) output_data.jacobian_grads[i] = DerivativeForm<2, dim, spacedim>(); if (data.update_each & update_jacobian_pushed_forward_grads) for (unsigned int i = 0; i < output_data.jacobian_pushed_forward_grads.size(); ++i) output_data.jacobian_pushed_forward_grads[i] = Tensor<3, spacedim>(); if (data.update_each & update_jacobian_2nd_derivatives) for (unsigned int i = 0; i < output_data.jacobian_2nd_derivatives.size(); ++i) output_data.jacobian_2nd_derivatives[i] = DerivativeForm<3, dim, spacedim>(); if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives) for (unsigned int i = 0; i < output_data.jacobian_pushed_forward_2nd_derivatives.size(); ++i) output_data.jacobian_pushed_forward_2nd_derivatives[i] = Tensor<4, spacedim>(); if (data.update_each & update_jacobian_3rd_derivatives) for (unsigned int i = 0; i < output_data.jacobian_3rd_derivatives.size(); ++i) output_data.jacobian_3rd_derivatives[i] = DerivativeForm<4, dim, spacedim>(); if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives) for (unsigned int i = 0; i < output_data.jacobian_pushed_forward_3rd_derivatives.size(); ++i) output_data.jacobian_pushed_forward_3rd_derivatives[i] = Tensor<5, spacedim>(); } } template <int dim, int spacedim> void MappingBox<dim, spacedim>::maybe_update_volume_elements( const InternalData &data) const { if (data.update_each & update_volume_elements) { double volume = data.cell_extents[0]; for (unsigned int d = 1; d < dim; ++d) volume *= data.cell_extents[d]; data.volume_element = volume; } } template <int dim, int spacedim> void MappingBox<dim, spacedim>::maybe_update_jacobians( const InternalData &data, const CellSimilarity::Similarity cell_similarity, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const { @endcode "compute" Jacobian at the quadrature points, which are all the same @code if (data.update_each & update_jacobians) if (cell_similarity != CellSimilarity::translation) for (unsigned int i = 0; i < output_data.jacobians.size(); ++i) { output_data.jacobians[i] = DerivativeForm<1, dim, spacedim>(); for (unsigned int j = 0; j < dim; ++j) output_data.jacobians[i][j][j] = data.cell_extents[j]; } } template <int dim, int spacedim> void MappingBox<dim, spacedim>::maybe_update_inverse_jacobians( const InternalData &data, const CellSimilarity::Similarity cell_similarity, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const { @endcode "compute" inverse Jacobian at the quadrature points, which are all the same @code if (data.update_each & update_inverse_jacobians) if (cell_similarity != CellSimilarity::translation) for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i) { output_data.inverse_jacobians[i] = Tensor<2, dim>(); for (unsigned int j = 0; j < dim; ++j) output_data.inverse_jacobians[i][j][j] = data.inverse_cell_extents[j]; } } template <int dim, int spacedim> CellSimilarity::Similarity MappingBox<dim, spacedim>::fill_fe_values( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature<dim> &quadrature, const typename Mapping<dim, spacedim>::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const { Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); @endcode convert data object to internal data for this class. fails with an exception if that is not possible @code Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr, ExcInternalError()); const InternalData &data = static_cast<const InternalData &>(internal_data); update_cell_extents(cell, cell_similarity, data); maybe_update_cell_quadrature_points(cell, data, quadrature.get_points(), output_data.quadrature_points); @endcode compute Jacobian determinant. all values are equal and are the product of the local lengths in each coordinate direction @code if (data.update_each & (update_JxW_values | update_volume_elements)) if (cell_similarity != CellSimilarity::translation) { double J = data.cell_extents[0]; for (unsigned int d = 1; d < dim; ++d) J *= data.cell_extents[d]; data.volume_element = J; if (data.update_each & update_JxW_values) for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i) output_data.JxW_values[i] = quadrature.weight(i); } maybe_update_jacobians(data, cell_similarity, output_data); maybe_update_jacobian_derivatives(data, cell_similarity, output_data); maybe_update_inverse_jacobians(data, cell_similarity, output_data); return cell_similarity; } template <int dim, int spacedim> void MappingBox<dim, spacedim>::fill_fe_subface_values( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature<dim - 1> &quadrature, const typename Mapping<dim, spacedim>::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const { (void)cell; (void)face_no; (void)subface_no; (void)quadrature; (void)internal_data; (void)output_data; DEAL_II_NOT_IMPLEMENTED(); } template <int dim, int spacedim> void MappingBox<dim, spacedim>::fill_fe_immersed_surface_values( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature, const typename Mapping<dim, spacedim>::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &output_data) const { AssertDimension(dim, spacedim); Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); @endcode Convert data object to internal data for this class. Fails with an exception if that is not possible. @code Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr, ExcInternalError()); const InternalData &data = static_cast<const InternalData &>(internal_data); update_cell_extents(cell, CellSimilarity::none, data); maybe_update_cell_quadrature_points(cell, data, quadrature.get_points(), output_data.quadrature_points); if (data.update_each & update_normal_vectors) for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i) output_data.normal_vectors[i] = quadrature.normal_vector(i); if (data.update_each & update_JxW_values) for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i) output_data.JxW_values[i] = quadrature.weight(i); maybe_update_volume_elements(data); maybe_update_jacobians(data, CellSimilarity::none, output_data); maybe_update_jacobian_derivatives(data, CellSimilarity::none, output_data); maybe_update_inverse_jacobians(data, CellSimilarity::none, output_data); } template <int dim, int spacedim> void MappingBox<dim, spacedim>::transform( const ArrayView<const Tensor<1, dim>> &input, const MappingKind mapping_kind, const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data, const ArrayView<Tensor<1, spacedim>> &output) const { AssertDimension(input.size(), output.size()); Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr, ExcInternalError()); const InternalData &data = static_cast<const InternalData &>(mapping_data); switch (mapping_kind) { case mapping_covariant: { Assert(data.update_each & update_covariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_covariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d = 0; d < dim; ++d) output[i][d] = input[i][d] * data.inverse_cell_extents[d]; return; } case mapping_contravariant: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d = 0; d < dim; ++d) output[i][d] = input[i][d] * data.cell_extents[d]; return; } case mapping_piola: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); Assert(data.update_each & update_volume_elements, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_volume_elements")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d = 0; d < dim; ++d) output[i][d] = input[i][d] * data.cell_extents[d] / data.volume_element; return; } default: DEAL_II_NOT_IMPLEMENTED(); } } template <int dim, int spacedim> void MappingBox<dim, spacedim>::transform( const ArrayView<const DerivativeForm<1, dim, spacedim>> &input, const MappingKind mapping_kind, const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data, const ArrayView<Tensor<2, spacedim>> &output) const { AssertDimension(input.size(), output.size()); Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr, ExcInternalError()); const InternalData &data = static_cast<const InternalData &>(mapping_data); switch (mapping_kind) { case mapping_covariant: { Assert(data.update_each & update_covariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_covariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.inverse_cell_extents[d2]; return; } case mapping_contravariant: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2]; return; } case mapping_covariant_gradient: { Assert(data.update_each & update_covariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_covariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.inverse_cell_extents[d2] * data.inverse_cell_extents[d1]; return; } case mapping_contravariant_gradient: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * data.inverse_cell_extents[d1]; return; } case mapping_piola: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); Assert(data.update_each & update_volume_elements, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_volume_elements")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] / data.volume_element; return; } case mapping_piola_gradient: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); Assert(data.update_each & update_volume_elements, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_volume_elements")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * data.inverse_cell_extents[d1] / data.volume_element; return; } default: DEAL_II_NOT_IMPLEMENTED(); } } template <int dim, int spacedim> void MappingBox<dim, spacedim>::transform( const ArrayView<const Tensor<2, dim>> &input, const MappingKind mapping_kind, const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data, const ArrayView<Tensor<2, spacedim>> &output) const { AssertDimension(input.size(), output.size()); Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr, ExcInternalError()); const InternalData &data = static_cast<const InternalData &>(mapping_data); switch (mapping_kind) { case mapping_covariant: { Assert(data.update_each & update_covariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_covariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.inverse_cell_extents[d2]; return; } case mapping_contravariant: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2]; return; } case mapping_covariant_gradient: { Assert(data.update_each & update_covariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_covariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.inverse_cell_extents[d2] * data.inverse_cell_extents[d1]; return; } case mapping_contravariant_gradient: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * data.inverse_cell_extents[d1]; return; } case mapping_piola: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); Assert(data.update_each & update_volume_elements, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_volume_elements")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] / data.volume_element; return; } case mapping_piola_gradient: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); Assert(data.update_each & update_volume_elements, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_volume_elements")); for (unsigned int i = 0; i < output.size(); ++i) for (unsigned int d1 = 0; d1 < dim; ++d1) for (unsigned int d2 = 0; d2 < dim; ++d2) output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * data.inverse_cell_extents[d1] / data.volume_element; return; } default: DEAL_II_NOT_IMPLEMENTED(); } } template <int dim, int spacedim> void MappingBox<dim, spacedim>::transform( const ArrayView<const DerivativeForm<2, dim, spacedim>> &input, const MappingKind mapping_kind, const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data, const ArrayView<Tensor<3, spacedim>> &output) const { AssertDimension(input.size(), output.size()); Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr, ExcInternalError()); const InternalData &data = static_cast<const InternalData &>(mapping_data); switch (mapping_kind) { case mapping_covariant_gradient: { Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_covariant_transformation")); for (unsigned int q = 0; q < output.size(); ++q) for (unsigned int i = 0; i < spacedim; ++i) for (unsigned int j = 0; j < spacedim; ++j) for (unsigned int k = 0; k < spacedim; ++k) { output[q][i][j][k] = input[q][i][j][k] * data.inverse_cell_extents[j] * data.inverse_cell_extents[k]; } return; } default: DEAL_II_NOT_IMPLEMENTED(); } } template <int dim, int spacedim> void MappingBox<dim, spacedim>::transform( const ArrayView<const Tensor<3, dim>> &input, const MappingKind mapping_kind, const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data, const ArrayView<Tensor<3, spacedim>> &output) const { AssertDimension(input.size(), output.size()); Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr, ExcInternalError()); const InternalData &data = static_cast<const InternalData &>(mapping_data); switch (mapping_kind) { case mapping_contravariant_hessian: { Assert(data.update_each & update_covariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_covariant_transformation")); Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); for (unsigned int q = 0; q < output.size(); ++q) for (unsigned int i = 0; i < spacedim; ++i) for (unsigned int j = 0; j < spacedim; ++j) for (unsigned int k = 0; k < spacedim; ++k) { output[q][i][j][k] = input[q][i][j][k] * data.cell_extents[i] * data.inverse_cell_extents[j] * data.inverse_cell_extents[k]; } return; } case mapping_covariant_hessian: { Assert(data.update_each & update_covariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_covariant_transformation")); for (unsigned int q = 0; q < output.size(); ++q) for (unsigned int i = 0; i < spacedim; ++i) for (unsigned int j = 0; j < spacedim; ++j) for (unsigned int k = 0; k < spacedim; ++k) { output[q][i][j][k] = input[q][i][j][k] * (data.inverse_cell_extents[i] * data.inverse_cell_extents[j]) * data.inverse_cell_extents[k]; } return; } case mapping_piola_hessian: { Assert(data.update_each & update_covariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_covariant_transformation")); Assert(data.update_each & update_contravariant_transformation, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_contravariant_transformation")); Assert(data.update_each & update_volume_elements, typename FEValuesBase<dim>::ExcAccessToUninitializedField( "update_volume_elements")); for (unsigned int q = 0; q < output.size(); ++q) for (unsigned int i = 0; i < spacedim; ++i) for (unsigned int j = 0; j < spacedim; ++j) for (unsigned int k = 0; k < spacedim; ++k) { output[q][i][j][k] = input[q][i][j][k] * (data.cell_extents[i] / data.volume_element * data.inverse_cell_extents[j]) * data.inverse_cell_extents[k]; } return; } default: DEAL_II_NOT_IMPLEMENTED(); } } template <int dim, int spacedim> Point<spacedim> MappingBox<dim, spacedim>::transform_unit_to_real_cell( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const Point<dim> &p) const { Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); Assert(dim == spacedim, ExcNotImplemented()); return boxes[polytope_translator.at(cell->active_cell_index())].unit_to_real( p); } template <int dim, int spacedim> Point<dim> MappingBox<dim, spacedim>::transform_real_to_unit_cell( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const Point<spacedim> &p) const { Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); Assert(dim == spacedim, ExcNotImplemented()); return boxes[polytope_translator.at(cell->active_cell_index())].real_to_unit( p); } template <int dim, int spacedim> void MappingBox<dim, spacedim>::transform_points_real_to_unit_cell( const typename Triangulation<dim, spacedim>::cell_iterator &cell, const ArrayView<const Point<spacedim>> &real_points, const ArrayView<Point<dim>> &unit_points) const { Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); AssertDimension(real_points.size(), unit_points.size()); if (dim != spacedim) DEAL_II_NOT_IMPLEMENTED(); for (unsigned int i = 0; i < real_points.size(); ++i) unit_points[i] = boxes[polytope_translator.at(cell->active_cell_index())].real_to_unit( real_points[i]); } template <int dim, int spacedim> std::unique_ptr<Mapping<dim, spacedim>> MappingBox<dim, spacedim>::clone() const { return std::make_unique<MappingBox<dim, spacedim>>(*this); }
explicit instantiations
1.9.8