diff --git a/agglomeration_poisson/CMakeLists.txt b/agglomeration_poisson/CMakeLists.txt new file mode 100644 index 00000000..fca2079c --- /dev/null +++ b/agglomeration_poisson/CMakeLists.txt @@ -0,0 +1,58 @@ +# Set the name of the project and target: +SET(TARGET "agglomeration_poisson") + +# Declare all source files the target consists of. +FILE(GLOB TARGET_SRC CONFIGURE_DEPENDS ${TARGET}.cc ${CMAKE_CURRENT_SOURCE_DIR}/source/*.cc) + +CMAKE_MINIMUM_REQUIRED(VERSION 3.13.4) + +FIND_PACKAGE(deal.II 9.6.0 + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() + +TARGET_INCLUDE_DIRECTORIES(${TARGET} PRIVATE include) +TARGET_COMPILE_DEFINITIONS(${TARGET} PRIVATE MESH_DIR="${CMAKE_SOURCE_DIR}/mesh") + + + +# ------------------------------------------------------------------- +# Coverage target: "make coverage" +# ------------------------------------------------------------------- +target_compile_options(${TARGET} PRIVATE + -O0 + -g + --coverage + -fprofile-arcs + -fprofile-update=atomic + -ftest-coverage +) +target_link_options(${TARGET} PRIVATE "--coverage") + +find_program(LCOV_EXEC lcov) +find_program(GENHTML_EXEC genhtml) + +if(LCOV_EXEC AND GENHTML_EXEC) + add_custom_target(coverage + COMMAND ${CMAKE_COMMAND} -E echo "Generating coverage report..." + COMMAND ${CMAKE_COMMAND} -E rm -rf coverage_report coverage.info + COMMAND ${LCOV_EXEC} --capture --directory . --output-file coverage.info + COMMAND ${LCOV_EXEC} --extract coverage.info "*${CMAKE_SOURCE_DIR}*" --output-file coverage.info + COMMAND ${GENHTML_EXEC} --function-coverage --output-directory coverage_report coverage.info + WORKING_DIRECTORY ${CMAKE_BINARY_DIR} + COMMENT "Running LCOV + genhtml" + VERBATIM + ) +else() + message(WARNING "lcov/genhtml not found - coverage target will not be available.") +endif() diff --git a/agglomeration_poisson/README.md b/agglomeration_poisson/README.md new file mode 100644 index 00000000..d01b5c73 --- /dev/null +++ b/agglomeration_poisson/README.md @@ -0,0 +1,135 @@ +# A Discontinuous Galerkin solver for the Poisson problem on general polytopal meshes generated through mesh agglomeration + +## Running the code: + +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 using + +`./DG_advection_reaction` + +on the command line. + +## The problem: +This program solves the problem, for $\Omega \in \mathbb{R^2}$ + +@f[ +\begin{cases} b \cdot \nabla u + c u = f \qquad \text{in } \Omega \\ +\qquad \qquad u=g \qquad \text{on } \partial_{-}\Omega \end{cases} +@f] + +where $g \in L^2(\partial_{-}\Omega)$ and $\partial_{-}\Omega=\{ x \in +\partial \Omega: b(x)\cdot n(x) <0\}$ is the inflow part of the +boundary, with $b=(b_1,b_2) \in \mathbb{R^2}$. As we know from +classical DG theory, we need to ensure that +@f[ +c(x) - \frac{1}{2}\nabla \cdot b \geq \gamma_0 >0 +@f] +for some positive $\gamma_0$ so that we have coercivity in $L^2$ at the continuous level. Discrete coercivity is achieved by using a stronger norm which takes care of jumps, see Di Pietro and Ern [2] for details. + + +## The weak formulation: + +As trial space we choose $V_h = \{ v_h \in L^2(\Omega): v_h \in P^1(\mathbb{T_h})\} \notin H^1(\Omega)$. If we integrate by parts and sum over all cells + +@f[ +\sum_{T \in \mathbb{T}_h} \Bigl( (-u,\beta \cdot \nabla v_h) _T + (c +u,v_h)_T + \bigl<(b \cdot n) u ,v_h \bigr>_{\partial T} \Bigr) = +(f,v_h)_{\Omega} +@f] + +and use the so-called DG magic formula and exploit the property $[bu]_{\mathbb{F}^i} = 0$ where $\mathbb{F}^i$ are set of internal faces we obtain the (unstable!) formulation: + +Find $u_h \in V_h$: + +@f[ + a_h(u_h,v_h) + b_h(u_h,v_h)=l(v_h) \qquad \forall v_h \in V_h +@f] +where +@f[ +a_h(u,v_h)=\sum_{T \in \mathbb{T}_h} \Bigl( (-u,b \cdot \nabla v_h) _T + (c u,v_h)_T \Bigr) +@f] + +@f[ +b_h(u,v_h)= \sum_{F \not \in \partial_{-}\Omega} \bigl< \{ bu\}, [v_h]\bigr>_F +@f] + +@f[ + l(v_h)= (f,v_h)_{\Omega} - \sum_{F \in \partial_{-}\Omega} \bigl< (b \cdot n) g,v_h \bigr>_F +@f] + +It's well known this formulation is coercive only in $L^2$, hence the formulation is unstable as we don't "see" the derivatives. To stabilize this, we can use a jump-penalty term, i.e. our $b_h$ is replaced by: + +@f[ +b_h^s(u_h,v_h)=b_h(u_h,v_h)+ \sum_{F \in \mathbb{F}^i} \bigl< c_F +[u_h],[v_h] \bigr> +@f] + +where $c_F>0$ is a function on each edge such that $c_F \geq \theta |b \cdot n|$ for some positive $\theta$. In this program, $\theta=\frac{1}{2}$ and $c_F = \frac{1}{2} |b \cdot n|$, which corresponds to an upwind formulation. Notice that consistency is trivially achieved, as $[u]_{\mathbb{F}^i} =0$. This formulation is stable in the energy norm + +@f[ + |||\cdot ||| = \Bigl(||\cdot||_{0,\Omega}^2 + \sum_{F \in + \mathbb{F}}||c_F^{\frac{1}{2}}[\cdot] ||_{0,F}^2 + \Bigr)^{\frac{1}{2}} +@f] + +(well defined on $H^1(\Omega) + V_h$) and moreover we have the a-priori bound: + +@f[ +|||u-u_h||| \leq C h^{k+\frac{1}{2}}||u||_{k+1,\Omega} +@f] + +valid for $u \in H^{k+1}(\Omega)$. + +See Brezzi-Marini-Süli [3] for more details. + +## A-posteriori error estimator: + +The estimator is the one proposed by Georgoulis, Edward Hall and Charalambos Makridakis in [3]. This approach is quite different with respect to other works in the field, as the authors are trying to develop an estimator for the original hyperbolic problem, rather than taking the hyperbolic regime as the vanishing diffusivity limit. + +The reliability is: + +@f[ +|||u-u_h|||^2 \leq C || \sqrt{b \cdot n}[u_h]||_{\Gamma^{-}}^2 + C +\sum_{T \in \mathbb{T}_h}\Bigl( ||\beta (g-u_h^+)||_{\partial_{-}T +\cap \partial_{-} \Omega}^2 +||f-c u_h - \Pi(f- cu_h)||_T^2 \Bigr) +@f] + +where: + +- $\Pi$ is the (local) $L^2$ orthogonal projection onto $V_h$ + +- $\Gamma$ is the skeleton of the mesh + +- $c$ is constant + +- $\beta = |b \cdot n|$ + +- $u_h^+$ is the interior trace from the current cell $T$ of a the finite element function $u_h$. + +## Test case: + +The following test case has been taken from [3]. Consider: +- $c=1$ +- $b=(1,1)$ +- $f$ to be such that the exact solution is $u(x,y)=\tanh(100(x+y-\frac{1}{2}))$ +This solution has an internal layer along the line $y=\frac{1}{2} -x$, hence we would like to see that part of the domain to be much more refined than the rest. + +The next image is the 3D view of the numerical solution: + +![Screenshot](./doc/images/warp_by_scalar_solution_layer.png) + +More interestingly, we see that the estimator has been able to capture the layer. Here a bulk-chasing criterion is used, with bottom fraction ´0.5´ and no coarsening. This mesh is obtained after 12 refinement cycles. +![Screenshot](./doc/images/refined_mesh_internal_layer.png) + +If we look at the decrease of the energy norm of the error in the globally refined case and in the adaptively case, with respect to the DoFs, we obtain: + +![Screenshot](./doc/images/adaptive_vs_global_refinement.png) + +## References +* [1] Emmanuil H. Georgoulis, Edward Hall and Charalambos Makridakis (2013), Error Control for Discontinuous Galerkin Methods for First Order Hyperbolic Problems. DOI: [10.1007/978-3-319-01818-8_8 +](https://link.springer.com/chapter/10.1007%2F978-3-319-01818-8_8) +* [2] Di Pietro, Daniele Antonio and Ern, Alexandre (2012), Mathematical Aspects of Discontinuous Galerkin Methods. ISBN: [978-3-642-22980-0](https://www.springer.com/gp/book/9783642229794) +* [3] Franco Brezzi, Luisa Donatella Marini and Endre Süli (2004) Discontinuous Galerkin Methods for First-Order Hyperbolic Problems. DOI: [10.1142/S0218202504003866](https://doi.org/10.1142/S0218202504003866) diff --git a/agglomeration_poisson/agglomeration_poisson.cc b/agglomeration_poisson/agglomeration_poisson.cc new file mode 100644 index 00000000..8fb029f5 --- /dev/null +++ b/agglomeration_poisson/agglomeration_poisson.cc @@ -0,0 +1,676 @@ +/* ----------------------------------------------------------------------------- + * + * 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 + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include + +struct ConvergenceInfo +{ + ConvergenceInfo() = default; + void + add(const std::pair> + &dofs_and_errs) + { + vec_data.push_back(dofs_and_errs); + } + + void + print() + { + Assert(vec_data.size() > 0, ExcInternalError()); + std::cout << std::left << "#DoFs, L2 error, H1 error" << std::endl; + + for (const auto &dof_and_errs : vec_data) + std::cout << std::scientific << dof_and_errs.first << ", " + << dof_and_errs.second.first << ", " + << dof_and_errs.second.second << std::endl; + } + + std::vector>> + vec_data; +}; + +enum class PartitionerType +{ + metis, + rtree, + no_partition +}; + +template +class RightHandSide : public Function +{ +public: + RightHandSide() + : Function() + {} + + virtual void + value_list(const std::vector> &points, + std::vector &values, + const unsigned int /*component*/) const override + { + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = 2 * numbers::PI * numbers::PI * + std::sin(numbers::PI * points[i][0]) * + std::sin(numbers::PI * points[i][1]); + } +}; + +template +class ExactSolution : public Function +{ +public: + ExactSolution() + : Function() + { + Assert(dim == 2, ExcNotImplemented()); + } + + virtual double + value(const Point &p, + const unsigned int /* component */ = 0) const override + { + return std::sin(numbers::PI * p[0]) * std::sin(numbers::PI * p[1]); + } + + + virtual void + value_list(const std::vector> &points, + std::vector &values, + const unsigned int /*component*/) const override + { + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = this->value(points[i]); + } + + virtual Tensor<1, dim> + gradient(const Point &p, + const unsigned int /* component */ = 0) const override + { + Tensor<1, dim> return_value; + return_value[0] = + numbers::PI * std::cos(numbers::PI * p[0]) * std::sin(numbers::PI * p[1]); + return_value[1] = + numbers::PI * std::cos(numbers::PI * p[1]) * std::sin(numbers::PI * p[0]); + return return_value; + } +}; + + +template +class Poisson +{ +private: + void + make_grid(); + void + setup_agglomeration(); + void + assemble_system(); + void + solve(); + void + output_results(); + + Triangulation tria; + MappingQ1 mapping; + FE_DGQ dg_fe; + std::unique_ptr> ah; + AffineConstraints constraints; + SparsityPattern sparsity; + DynamicSparsityPattern dsp; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + std::unique_ptr> cached_tria; + std::unique_ptr> rhs_function; + std::unique_ptr> analytical_solution; + +public: + Poisson(const PartitionerType &partitioner_type = PartitionerType::rtree, + const unsigned int = 0, + const unsigned int = 0, + const unsigned int fe_degree = 1); + void + run(); + + types::global_dof_index + get_n_dofs() const; + + std::pair + get_error() const; + + PartitionerType partitioner_type; + unsigned int extraction_level; + unsigned int n_subdomains; + double penalty_constant = 60.; // 10*(p+1)(p+d) for p = 1 and d = 2 => 60 + double l2_err; + double semih1_err; +}; + +template +Poisson::Poisson(const PartitionerType &partitioner_type, + const unsigned int extraction_level, + const unsigned int n_subdomains, + const unsigned int fe_degree) + : mapping() + , dg_fe(fe_degree) + , partitioner_type(partitioner_type) + , extraction_level(extraction_level) + , n_subdomains(n_subdomains) + , penalty_constant(10. * (fe_degree + 1) * (fe_degree + dim)) +{ + // Initialize manufactured solution. + analytical_solution = std::make_unique>(); + rhs_function = std::make_unique>(); + constraints.close(); +} + +template +void +Poisson::make_grid() +{ + GridIn grid_in; + grid_in.attach_triangulation(tria); + std::ifstream gmsh_file(std::string(MESH_DIR) + + "/unit_square_quad_unstructured.msh"); + grid_in.read_msh(gmsh_file); + tria.refine_global(2); + + std::cout << "Size of tria: " << tria.n_active_cells() << std::endl; + cached_tria = std::make_unique>(tria, mapping); + ah = std::make_unique>(*cached_tria); + + if (partitioner_type == PartitionerType::metis) + { + // Partition the triangulation with graph partitioner. + auto start = std::chrono::system_clock::now(); + GridTools::partition_triangulation(n_subdomains, + tria, + SparsityTools::Partitioner::metis); + + std::vector< + std::vector::active_cell_iterator>> + cells_per_subdomain(n_subdomains); + for (const auto &cell : tria.active_cell_iterators()) + cells_per_subdomain[cell->subdomain_id()].push_back(cell); + + // For every subdomain, agglomerate elements together + for (std::size_t i = 0; i < n_subdomains; ++i) + ah->define_agglomerate(cells_per_subdomain[i]); + + std::chrono::duration wctduration = + (std::chrono::system_clock::now() - start); + std::cout << "METIS built in " << wctduration.count() + << " seconds [wall clock]" << std::endl; + } + else if (partitioner_type == PartitionerType::rtree) + { + // Partition with Rtree + + namespace bgi = boost::geometry::index; + static constexpr unsigned int max_elem_per_node = + PolyUtils::constexpr_pow(2, dim); + std::vector, + typename Triangulation::active_cell_iterator>> + boxes(tria.n_active_cells()); + unsigned int i = 0; + for (const auto &cell : tria.active_cell_iterators()) + boxes[i++] = std::make_pair(mapping.get_bounding_box(cell), cell); + + auto start = std::chrono::system_clock::now(); + auto tree = pack_rtree>(boxes); + + CellsAgglomerator agglomerator{tree, + extraction_level}; + const auto vec_agglomerates = agglomerator.extract_agglomerates(); + + // Flag elements for agglomeration + for (const auto &agglo : vec_agglomerates) + ah->define_agglomerate(agglo); + + std::chrono::duration wctduration = + (std::chrono::system_clock::now() - start); + std::cout << "R-tree agglomerates built in " << wctduration.count() + << " seconds [wall clock]" << std::endl; + } + else if (partitioner_type == PartitionerType::no_partition) + { + } + else + { + Assert(false, ExcMessage("Wrong partitioning.")); + } + n_subdomains = ah->n_agglomerates(); + std::cout << "N subdomains = " << n_subdomains << std::endl; +} + +template +void +Poisson::setup_agglomeration() +{ + if (partitioner_type == PartitionerType::no_partition) + { + // No partitioning means that each cell is a master cell + for (const auto &cell : tria.active_cell_iterators()) + ah->define_agglomerate({cell}); + } + + ah->distribute_agglomerated_dofs(dg_fe); + ah->create_agglomeration_sparsity_pattern(dsp); + sparsity.copy_from(dsp); + + { + std::string partitioner; + if (partitioner_type == PartitionerType::metis) + partitioner = "metis"; + else if (partitioner_type == PartitionerType::rtree) + partitioner = "rtree"; + else + partitioner = "no_partitioning"; + + const std::string filename = + "grid_" + partitioner + "_" + std::to_string(n_subdomains) + ".vtu"; + std::ofstream output(filename); + + DataOut data_out; + data_out.attach_dof_handler(ah->agglo_dh); + + Vector agglomerated(tria.n_active_cells()); + Vector agglo_idx(tria.n_active_cells()); + for (const auto &cell : tria.active_cell_iterators()) + { + agglomerated[cell->active_cell_index()] = + ah->get_relationships()[cell->active_cell_index()]; + agglo_idx[cell->active_cell_index()] = cell->subdomain_id(); + } + data_out.add_data_vector(agglomerated, + "agglo_relationships", + DataOut::type_cell_data); + data_out.add_data_vector(agglo_idx, + "agglomerated_idx", + DataOut::type_cell_data); + data_out.build_patches(mapping); + data_out.write_vtu(output); + } +} + +template +void +Poisson::assemble_system() +{ + system_matrix.reinit(sparsity); + solution.reinit(ah->n_dofs()); + system_rhs.reinit(ah->n_dofs()); + + const unsigned int quadrature_degree = dg_fe.get_degree() + 1; + const unsigned int face_quadrature_degree = dg_fe.get_degree() + 1; + + ah->initialize_fe_values(QGauss(quadrature_degree), + update_gradients | update_JxW_values | + update_quadrature_points | update_JxW_values | + update_values, + QGauss(face_quadrature_degree)); + + const unsigned int dofs_per_cell = ah->n_dofs_per_cell(); + std::cout << "DoFs per cell: " << dofs_per_cell << std::endl; + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + // Next, we define the four dofsxdofs matrices needed to assemble jumps and + // averages. + FullMatrix M11(dofs_per_cell, dofs_per_cell); + FullMatrix M12(dofs_per_cell, dofs_per_cell); + FullMatrix M21(dofs_per_cell, dofs_per_cell); + FullMatrix M22(dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + for (const auto &polytope : ah->polytope_iterators()) + { + cell_matrix = 0; + cell_rhs = 0; + const auto &agglo_values = ah->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 rhs(n_qpoints); + rhs_function->value_list(q_points, rhs); + + 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); + } + cell_rhs(i) += agglo_values.shape_value(i, q_index) * + rhs[q_index] * agglo_values.JxW(q_index); + } + } + + // Face terms + const unsigned int n_faces = polytope->n_faces(); + AssertThrow(n_faces > 0, + ExcMessage( + "Invalid element: at least 4 faces are required.")); + + auto polygon_boundary_vertices = polytope->polytope_boundary(); + for (unsigned int f = 0; f < n_faces; ++f) + { + if (polytope->at_boundary(f)) + { + // std::cout << "at boundary!" << std::endl; + const auto &fe_face = ah->reinit(polytope, f); + + const unsigned int dofs_per_cell = fe_face.dofs_per_cell; + + const auto &face_q_points = fe_face.get_quadrature_points(); + std::vector analytical_solution_values( + face_q_points.size()); + analytical_solution->value_list(face_q_points, + analytical_solution_values, + 1); + + // Get normal vectors seen from each agglomeration. + const auto &normals = fe_face.get_normal_vectors(); + + const double penalty = + penalty_constant / std::fabs(polytope->diameter()); + + for (unsigned int q_index : fe_face.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_face.shape_value(i, q_index) * + fe_face.shape_grad(j, q_index) * + normals[q_index] - + fe_face.shape_grad(i, q_index) * normals[q_index] * + fe_face.shape_value(j, q_index) + + (penalty)*fe_face.shape_value(i, q_index) * + fe_face.shape_value(j, q_index)) * + fe_face.JxW(q_index); + } + cell_rhs(i) += + (penalty * analytical_solution_values[q_index] * + fe_face.shape_value(i, q_index) - + fe_face.shape_grad(i, q_index) * normals[q_index] * + analytical_solution_values[q_index]) * + fe_face.JxW(q_index); + } + } + } + else + { + const auto &neigh_polytope = polytope->neighbor(f); + + // This is necessary to loop over internal faces only once. + if (polytope->index() < neigh_polytope->index()) + { + unsigned int nofn = + polytope->neighbor_of_agglomerated_neighbor(f); + + const auto &fe_faces = + ah->reinit_interface(polytope, neigh_polytope, f, nofn); + const auto &fe_faces0 = fe_faces.first; + const auto &fe_faces1 = fe_faces.second; + + std::vector + local_dof_indices_neighbor(dofs_per_cell); + + M11 = 0.; + M12 = 0.; + M21 = 0.; + M22 = 0.; + + const auto &normals = fe_faces0.get_normal_vectors(); + + const double penalty = + penalty_constant / std::fabs(polytope->diameter()); + + // M11 + 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)*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)*fe_faces0.shape_value(i, q_index) * + fe_faces1.shape_value(j, q_index)) * + fe_faces1.JxW(q_index); + + // A10 + 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)*fe_faces1.shape_value(i, q_index) * + fe_faces0.shape_value(j, q_index)) * + fe_faces1.JxW(q_index); + + // A11 + 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)*fe_faces1.shape_value(i, q_index) * + fe_faces1.shape_value(j, q_index)) * + fe_faces1.JxW(q_index); + } + } + } + + 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); + } // Loop only once trough internal faces + } + } // Loop over faces of current cell + + // distribute DoFs + constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); + } // Loop over cells +} + +template +void +Poisson::solve() +{ + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + A_direct.vmult(solution, system_rhs); +} + +template +void +Poisson::output_results() +{ + { + std::string partitioner; + if (partitioner_type == PartitionerType::metis) + partitioner = "metis"; + else if (partitioner_type == PartitionerType::rtree) + partitioner = "rtree"; + else + partitioner = "no_partitioning"; + + const std::string filename = "interpolated_solution_" + partitioner + "_" + + std::to_string(n_subdomains) + ".vtu"; + std::ofstream output(filename); + + DataOut data_out; + Vector interpolated_solution; + PolyUtils::interpolate_to_fine_grid(*ah, + interpolated_solution, + solution, + true /*on_the_fly*/); + data_out.attach_dof_handler(ah->output_dh); + data_out.add_data_vector(interpolated_solution, + "u", + DataOut::type_dof_data); + + Vector agglo_idx(tria.n_active_cells()); + + // Mark fine cells belonging to the same agglomerate. + for (const auto &polytope : ah->polytope_iterators()) + { + const types::global_cell_index polytope_index = polytope->index(); + const auto &patch_of_cells = polytope->get_agglomerate(); // fine cells + // Flag them + for (const auto &cell : patch_of_cells) + agglo_idx[cell->active_cell_index()] = polytope_index; + } + + data_out.add_data_vector(agglo_idx, + "agglo_idx", + DataOut::type_cell_data); + + data_out.build_patches(mapping); + data_out.write_vtu(output); + + // Compute L2 and semiH1 norm of error + std::vector errors; + PolyUtils::compute_global_error(*ah, + solution, + *analytical_solution, + {VectorTools::L2_norm, + VectorTools::H1_seminorm}, + errors); + l2_err = errors[0]; + semih1_err = errors[1]; + } +} + +template +inline types::global_dof_index +Poisson::get_n_dofs() const +{ + return ah->n_dofs(); +} + +template +inline std::pair +Poisson::get_error() const +{ + return std::make_pair(l2_err, semih1_err); +} + +template +void +Poisson::run() +{ + make_grid(); + setup_agglomeration(); + auto start = std::chrono::high_resolution_clock::now(); + assemble_system(); + auto stop = std::chrono::high_resolution_clock::now(); + auto duration = + std::chrono::duration_cast(stop - start); + + std::cout << "Time taken by assemble_system(): " << duration.count() + << " seconds" << std::endl; + solve(); + output_results(); +} + +int +main() +{ + ConvergenceInfo convergence_info; + + for (unsigned int fe_degree : {1}) //, 2, 3}) + { + std::cout << "Running with FE degree: " << fe_degree << std::endl; + Poisson<2> poisson_problem{PartitionerType::rtree, + 4 /* extraction_level */, + 256 /* n_subdomains */, + fe_degree}; + poisson_problem.run(); + convergence_info.add( + std::make_pair>( + poisson_problem.get_n_dofs(), poisson_problem.get_error())); + std::cout << std::endl; + } + + std::cout << "Convergence table:" << std::endl; + convergence_info.print(); + std::cout << std::endl; + + return 0; +} diff --git a/agglomeration_poisson/doc/author b/agglomeration_poisson/doc/author new file mode 100644 index 00000000..74d91324 --- /dev/null +++ b/agglomeration_poisson/doc/author @@ -0,0 +1,4 @@ +Marco Feder +Pasquale Claudio Africa +Xinping Gui +Andrea Cangiani diff --git a/agglomeration_poisson/doc/builds-on b/agglomeration_poisson/doc/builds-on new file mode 100644 index 00000000..154f8bf8 --- /dev/null +++ b/agglomeration_poisson/doc/builds-on @@ -0,0 +1 @@ +step-12 step-74 diff --git a/agglomeration_poisson/doc/entry-name b/agglomeration_poisson/doc/entry-name new file mode 100644 index 00000000..0cdacd6e --- /dev/null +++ b/agglomeration_poisson/doc/entry-name @@ -0,0 +1 @@ +An agglomeration-based solver for the Poisson problem diff --git a/agglomeration_poisson/doc/images/adaptive_vs_global_refinement.png b/agglomeration_poisson/doc/images/adaptive_vs_global_refinement.png new file mode 100644 index 00000000..70d30393 Binary files /dev/null and b/agglomeration_poisson/doc/images/adaptive_vs_global_refinement.png differ diff --git a/agglomeration_poisson/doc/images/refined_mesh_internal_layer.png b/agglomeration_poisson/doc/images/refined_mesh_internal_layer.png new file mode 100644 index 00000000..bd38f0ec Binary files /dev/null and b/agglomeration_poisson/doc/images/refined_mesh_internal_layer.png differ diff --git a/agglomeration_poisson/doc/images/warp_by_scalar_solution_layer.png b/agglomeration_poisson/doc/images/warp_by_scalar_solution_layer.png new file mode 100644 index 00000000..f0c30f3e Binary files /dev/null and b/agglomeration_poisson/doc/images/warp_by_scalar_solution_layer.png differ diff --git a/agglomeration_poisson/doc/tooltip b/agglomeration_poisson/doc/tooltip new file mode 100644 index 00000000..8f0328bc --- /dev/null +++ b/agglomeration_poisson/doc/tooltip @@ -0,0 +1 @@ +A Discontinuous Galerkin solver for the Poisson problem on general polytopal meshes generated through mesh agglomeration diff --git a/agglomeration_poisson/include/agglomeration_accessor.h b/agglomeration_poisson/include/agglomeration_accessor.h new file mode 100644 index 00000000..8eb803c8 --- /dev/null +++ b/agglomeration_poisson/include/agglomeration_accessor.h @@ -0,0 +1,794 @@ +/* ----------------------------------------------------------------------------- + * + * 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 agglomeration_accessor_h +#define agglomeration_accessor_h + +#include + +#include +#include + +#include + +#include + +using namespace dealii; + + +// Forward declarations +#ifndef DOXYGEN +template +class AgglomerationHandler; +template +class AgglomerationIterator; +#endif + + +/** + * Accessor class used by AgglomerationIterator to access agglomeration data. + */ +template +class AgglomerationAccessor +{ +public: + /** + * Type for storing the polygons in an agglomerate. + */ + using AgglomerationContainer = + std::vector::active_cell_iterator>; + + + /** + * Get the DoFs indices associated to the agglomerate. + */ + void + get_dof_indices(std::vector &) const; + + /** + * Return, for a cell, the number of faces. In case the cell is a standard + * cell, then the number of faces is the classical one. If it's a master cell, + * then it returns the number of faces of the agglomeration identified by the + * master cell itself. + */ + unsigned int + n_faces() const; + + /** + * Return the number of deal.II faces that is building a polygon face. + */ + unsigned int + n_agglomerated_faces() const; + + /** + * Return the agglomerate which shares face f. + */ + const AgglomerationIterator + neighbor(const unsigned int f) const; + + /** + * Return the present index (seen from the neighboring agglomerate) of the + * present face f. + */ + unsigned int + neighbor_of_agglomerated_neighbor(const unsigned int f) const; + + /** + * + * This function generalizes the behaviour of cell->face(f)->at_boundary() + * in the case where f is an index out of the range [0,..., n_faces). + * In practice, if you call this function with a standard deal.II cell, you + * have precisely the same result as calling cell->face(f)->at_boundary(). + * Otherwise, if the cell is a master one, you have a boolean returning true + * is that face for the agglomeration is on the boundary or not. + */ + bool + at_boundary(const unsigned int f) const; + + /** + * Return a vector of face iterators describing the boundary of agglomerate. + */ + const std::vector::active_face_iterator> & + polytope_boundary() const; + + /** + * + * Return the volume of a polytope. + */ + double + volume() const; + + /** + * Return the diameter of the present polytopal element. + */ + double + diameter() const; + + /** + * Returns the deal.II cells that build the agglomerate. + */ + AgglomerationContainer + get_agglomerate() const; + + /** + * Return the BoundingBox which bounds the present polytope. In case the + * present polytope is not locally owned, it returns the BoundingBox of that + * ghosted polytope. + */ + const BoundingBox & + get_bounding_box() const; + + /** + * Return the index of the present polytope. + */ + types::global_cell_index + index() const; + + /** + * Returns an active cell iterator for the dof_handler, matching the polytope + * referenced by the input iterator. The type of the returned object is a + * DoFHandler::active_cell_iterator which can be used to initialize + * FiniteElement data. + */ + typename DoFHandler::active_cell_iterator + as_dof_handler_iterator(const DoFHandler &dof_handler) const; + + /** + * Returns the number of classical deal.II cells that are building the present + * polygon. + */ + unsigned int + n_background_cells() const; + + /* Returns true if this polygon is owned by the current processor. On a serial + * Triangulation this returs always true, but may yield false for a + * parallel::distributed::Triangulation. + */ + bool + is_locally_owned() const; + + /** + * The polytopal analogue of CellAccessor::id(). It provides a way to uniquely + * identify cells in a parallel Triangulation such as a + * parallel::distributed::Triangulation. + */ + CellId + id() const; + + /** + * The polytopal analogue of CellAccessor::subdomain_id(). In case of a serial + * Triangulation, it returns the numbers::invalid_subdomain_id. + */ + types::subdomain_id + subdomain_id() const; + + /** + * Returns a vector of indices identifying the children polytopes. + */ + inline const std::vector & + children() const; + + /** + * Returns the FiniteElement object used by the current polytope. + * This function should only be called after the corresponding agglomeration + * handler has invoked distribute_agglomerated_dofs(). + */ + const FiniteElement & + get_fe() const; + + /** + * Sets the active finite element index. + * This function should be called when using hp::FECollection to specify + * which finite element in the collection is assigned to the current polytope. + */ + void + set_active_fe_index(const types::fe_index index) const; + + /** + * Returns the index of the active finite element. + * When using hp::FECollection, this function retrieves the index + * of the finite element assigned to the current polytope. + */ + types::fe_index + active_fe_index() const; + +private: + /** + * Private default constructor. This is not supposed to be used and hence will + * throw. + */ + AgglomerationAccessor(); + + /** + * Private constructor for an agglomerate. This is meant to be invoked by + * the AgglomerationIterator class. It takes as input the master cell of the + * agglomerate and a pointer to the handler. + */ + AgglomerationAccessor( + const typename Triangulation::active_cell_iterator + &master_cell, + const AgglomerationHandler *ah); + + /** + * Same as above, but needed when the argument @p cells is a ghost cell. + */ + AgglomerationAccessor( + const typename Triangulation::active_cell_iterator &cell, + const CellId &cell_id, + const AgglomerationHandler *ah); + + /** + * Default destructor. + */ + ~AgglomerationAccessor() = default; + + + /** + * The unique deal.II cell associated to the present polytope. + */ + typename Triangulation::active_cell_iterator master_cell; + + /** + * The index of the present polytope. + */ + types::global_cell_index present_index; + + /** + * The index of the present polytope. + */ + CellId present_id; + + /** + * The rank owning of the present polytope. + */ + types::subdomain_id present_subdomain_id; + + /** + * A pointer to the Handler. + */ + AgglomerationHandler *handler; + + /** + * Comparison operator for Accessor. Two accessors are equal if they refer to + * the same polytopal element. + */ + bool + operator==(const AgglomerationAccessor &other) const; + + /** + * Compare for inequality. + */ + bool + operator!=(const AgglomerationAccessor &other) const; + + /** + * Move to the next cell in the polytopal mesh. + */ + void + next(); + + /** + * Move to the previous cell in the polytopal mesh. + */ + void + prev(); + + /** + * Returns the slaves of the present agglomeration. + */ + const AgglomerationContainer & + get_slaves() const; + + unsigned int + n_agglomerated_faces_per_cell( + const typename Triangulation::active_cell_iterator &cell) + const; + + template + friend class AgglomerationIterator; +}; + + + +template +unsigned int +AgglomerationAccessor::n_agglomerated_faces_per_cell( + const typename Triangulation::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() && + !handler->are_cells_agglomerated(cell, neighboring_cell))) + { + ++n_neighbors; + } + } + return n_neighbors; +} + + + +template +unsigned int +AgglomerationAccessor::n_faces() const +{ + Assert(!handler->is_slave_cell(master_cell), + ExcMessage("You cannot pass a slave cell.")); + return handler->number_of_agglomerated_faces[present_index]; +} + + + +template +const AgglomerationIterator +AgglomerationAccessor::neighbor(const unsigned int f) const +{ + if (!at_boundary(f)) + { + if (master_cell->is_ghost()) + { + // The following path is needed when the present function is called + // from neighbor_of_neighbor() + + const unsigned int sender_rank = master_cell->subdomain_id(); + + const CellId &master_id_ghosted_neighbor = + handler->recv_ghosted_master_id.at(sender_rank) + .at(present_id) + .at(f); + + // Use the id of the master cell to uniquely identify the neighboring + // agglomerate + + return {master_cell, + master_id_ghosted_neighbor, + handler}; // dummy master? + } + + const types::global_cell_index polytope_index = + handler->master2polygon.at(master_cell->active_cell_index()); + + const auto &neigh = + handler->polytope_cache.cell_face_at_boundary.at({polytope_index, f}) + .second; + + + if (neigh->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator cell_dh( + *neigh, &(handler->agglo_dh)); + return {cell_dh, handler}; + } + else + { + // Get master_id from the neighboring ghost polytope. This uniquely + // identifies the neighboring polytope among all processors. + const CellId &master_id_neighbor = + handler->polytope_cache.ghosted_master_id.at({present_id, f}); + + // Use the id of the master cell to uniquely identify the neighboring + // agglomerate + return {neigh, master_id_neighbor, handler}; + } + } + else + { + return {}; + } +} + + + +template +unsigned int +AgglomerationAccessor::neighbor_of_agglomerated_neighbor( + const unsigned int f) const +{ + // First, make sure it's not a boundary face. + if (!at_boundary(f)) + { + const auto &neigh_polytope = + neighbor(f); // returns the neighboring master and id + + AssertThrow(neigh_polytope.state() == IteratorState::valid, + ExcInternalError()); + + unsigned int n_faces_agglomerated_neighbor; + + // if it is locally owned, retrieve the number of faces + if (neigh_polytope->is_locally_owned()) + { + n_faces_agglomerated_neighbor = neigh_polytope->n_faces(); + } + else + { + // 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. + const CellId &master_id_neighbor = neigh_polytope->id(); + + // Then, get the neighboring rank + const unsigned int sender_rank = neigh_polytope->subdomain_id(); + + // From the neighboring rank, use the CellId of the neighboring + // polytope to get the number of its faces. + n_faces_agglomerated_neighbor = + handler->recv_n_faces.at(sender_rank).at(master_id_neighbor); + } + + + // Loop over all faces of neighboring agglomerate + for (unsigned int f_out = 0; f_out < n_faces_agglomerated_neighbor; + ++f_out) + { + // Check if same CellId + if (neigh_polytope->neighbor(f_out).state() == IteratorState::valid) + if (neigh_polytope->neighbor(f_out)->id() == present_id) + return f_out; + } + return numbers::invalid_unsigned_int; + } + else + { + // Face is at boundary + return numbers::invalid_unsigned_int; + } +} + +// ------------------------------ inline functions ------------------------- + +template +inline AgglomerationAccessor::AgglomerationAccessor() +{} + + + +template +inline AgglomerationAccessor::AgglomerationAccessor( + const typename Triangulation::active_cell_iterator &cell, + const AgglomerationHandler *ah) +{ + handler = const_cast *>(ah); + if (&(*handler->master_cells_container.end()) == std::addressof(cell)) + { + present_index = handler->master_cells_container.size(); + master_cell = *handler->master_cells_container.end(); + present_id = CellId(); // invalid id (TODO) + present_subdomain_id = numbers::invalid_subdomain_id; + } + else + { + present_index = handler->master2polygon.at(cell->active_cell_index()); + master_cell = cell; + present_id = master_cell->id(); + present_subdomain_id = master_cell->subdomain_id(); + } +} + + + +template +inline AgglomerationAccessor::AgglomerationAccessor( + const typename Triangulation::active_cell_iterator &neigh_cell, + const CellId &master_cell_id, + const AgglomerationHandler *ah) +{ + Assert(neigh_cell->is_ghost(), ExcInternalError()); + // neigh_cell is ghosted + + handler = const_cast *>(ah); + master_cell = neigh_cell; + present_index = numbers::invalid_unsigned_int; + // neigh_cell is ghosted, use the CellId of that agglomerate + present_id = master_cell_id; + present_subdomain_id = master_cell->subdomain_id(); +} + + + +template +inline void +AgglomerationAccessor::get_dof_indices( + std::vector &dof_indices) const +{ + Assert(dof_indices.size() > 0, + ExcMessage( + "The vector of DoFs indices must be already properly resized.")); + if (is_locally_owned()) + { + // Forward the call to the master cell + typename DoFHandler::cell_iterator master_cell_dh( + *master_cell, &(handler->agglo_dh)); + master_cell_dh->get_dof_indices(dof_indices); + } + else + { + const std::vector &recv_dof_indices = + handler->recv_ghost_dofs.at(present_subdomain_id).at(present_id); + + std::copy(recv_dof_indices.cbegin(), + recv_dof_indices.cend(), + dof_indices.begin()); + } +} + + + +template +inline typename AgglomerationAccessor::AgglomerationContainer +AgglomerationAccessor::get_agglomerate() const +{ + auto agglomeration = get_slaves(); + agglomeration.push_back(master_cell); + return agglomeration; +} + + + +template +inline const std::vector::active_face_iterator> & +AgglomerationAccessor::polytope_boundary() const +{ + return handler->polygon_boundary[master_cell]; +} + + + +template +inline double +AgglomerationAccessor::diameter() const +{ + Assert(!handler->is_slave_cell(master_cell), + ExcMessage("The present function cannot be called for slave cells.")); + + if (handler->is_master_cell(master_cell)) + { + // Get the bounding box associated with the master cell + const auto &bdary_pts = + handler->bboxes[present_index].get_boundary_points(); + return (bdary_pts.second - bdary_pts.first).norm(); + } + else + { + // Standard deal.II way to get the measure of a cell. + return master_cell->diameter(); + } +} + + + +template +inline const BoundingBox & +AgglomerationAccessor::get_bounding_box() const +{ + if (is_locally_owned()) + return handler->bboxes[present_index]; + else + return handler->recv_ghosted_bbox.at(present_subdomain_id).at(present_id); +} + + + +template +inline double +AgglomerationAccessor::volume() const +{ + Assert(!handler->is_slave_cell(master_cell), + ExcMessage("The present function cannot be called for slave cells.")); + + if (handler->is_master_cell(master_cell)) + { + return handler->bboxes[present_index].volume(); + } + else + { + return master_cell->measure(); + } +} + + + +template +inline void +AgglomerationAccessor::next() +{ + // Increment the present index and update the polytope + ++present_index; + + // Make sure not to query the CellId if it's past the last + if (present_index < handler->master_cells_container.size()) + { + master_cell = handler->master_cells_container[present_index]; + present_id = master_cell->id(); + present_subdomain_id = master_cell->subdomain_id(); + } +} + + + +template +inline void +AgglomerationAccessor::prev() +{ + // Decrement the present index and update the polytope + --present_index; + master_cell = handler->master_cells_container[present_index]; + present_id = master_cell->id(); +} + + +template +inline bool +AgglomerationAccessor::operator==( + const AgglomerationAccessor &other) const +{ + return present_index == other.present_index; +} + +template +inline bool +AgglomerationAccessor::operator!=( + const AgglomerationAccessor &other) const +{ + return !(*this == other); +} + + + +template +inline types::global_cell_index +AgglomerationAccessor::index() const +{ + return present_index; +} + + + +template +typename DoFHandler::active_cell_iterator +AgglomerationAccessor::as_dof_handler_iterator( + const DoFHandler &dof_handler) const +{ + // Forward the call to the master cell using the right DoFHandler. + return master_cell->as_dof_handler_iterator(dof_handler); +} + + + +template +inline const typename AgglomerationAccessor::AgglomerationContainer & +AgglomerationAccessor::get_slaves() const +{ + return handler->master2slaves.at(master_cell->active_cell_index()); +} + + + +template +inline unsigned int +AgglomerationAccessor::n_background_cells() const +{ + AssertThrow(get_agglomerate().size() > 0, ExcMessage("Empty agglomeration.")); + return get_agglomerate().size(); +} + + + +template +unsigned int +AgglomerationAccessor::n_agglomerated_faces() const +{ + const auto &agglomeration = get_agglomerate(); + unsigned int n_neighbors = 0; + for (const auto &cell : agglomeration) + n_neighbors += n_agglomerated_faces_per_cell(cell); + return n_neighbors; +} + + + +template +inline bool +AgglomerationAccessor::at_boundary(const unsigned int f) const +{ + if (master_cell->is_ghost()) + { + const unsigned int sender_rank = master_cell->subdomain_id(); + return handler->recv_bdary_info.at(sender_rank).at(present_id).at(f); + } + else + { + Assert(!handler->is_slave_cell(master_cell), + ExcMessage( + "This function should not be called for a slave cell.")); + + + typename DoFHandler::active_cell_iterator cell_dh( + *master_cell, &(handler->agglo_dh)); + return handler->at_boundary(cell_dh, f); + } +} + + + +template +inline bool +AgglomerationAccessor::is_locally_owned() const +{ + return master_cell->is_locally_owned(); +} + + + +template +inline CellId +AgglomerationAccessor::id() const +{ + return present_id; +} + + + +template +inline types::subdomain_id +AgglomerationAccessor::subdomain_id() const +{ + return present_subdomain_id; +} + +template +inline const std::vector & +AgglomerationAccessor::children() const +{ + Assert(!handler->parent_child_info.empty(), ExcInternalError()); + return handler->parent_child_info.at( + {present_index, handler->present_extraction_level}); +} + +template +inline const FiniteElement & +AgglomerationAccessor::get_fe() const +{ + typename DoFHandler::active_cell_iterator + master_cell_as_dof_handler_iterator = + master_cell->as_dof_handler_iterator(handler->agglo_dh); + return master_cell_as_dof_handler_iterator->get_fe(); +} + +template +inline void +AgglomerationAccessor::set_active_fe_index( + const types::fe_index index) const +{ + Assert(!handler->is_slave_cell(master_cell), + ExcMessage("The present function cannot be called for slave cells.")); + typename DoFHandler::active_cell_iterator + master_cell_as_dof_handler_iterator = + master_cell->as_dof_handler_iterator(handler->agglo_dh); + master_cell_as_dof_handler_iterator->set_active_fe_index(index); +} + +template +inline types::fe_index +AgglomerationAccessor::active_fe_index() const +{ + typename DoFHandler::active_cell_iterator + master_cell_as_dof_handler_iterator = + master_cell->as_dof_handler_iterator(handler->agglo_dh); + return master_cell_as_dof_handler_iterator->active_fe_index(); +} + +#endif diff --git a/agglomeration_poisson/include/agglomeration_handler.h b/agglomeration_poisson/include/agglomeration_handler.h new file mode 100644 index 00000000..26df6303 --- /dev/null +++ b/agglomeration_poisson/include/agglomeration_handler.h @@ -0,0 +1,1250 @@ +/* ----------------------------------------------------------------------------- + * + * 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 agglomeration_handler_h +#define agglomeration_handler_h + +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include + +#include +#include + +using namespace dealii; + +// Forward declarations +template +class AgglomerationHandler; + +namespace dealii +{ + namespace internal + { + /** + * Helper class to reinit finite element spaces on polytopal cells. + */ + template + class AgglomerationHandlerImplementation; + } // namespace internal +} // namespace dealii + + + +/** + * Helper class for the storage of connectivity information of the polytopal + * grid. + */ +namespace dealii +{ + namespace internal + { + template + class PolytopeCache + { + public: + /** + * Default constructor. + */ + PolytopeCache() = default; + + /** + * Destructor. It simply calls clear() for all of its members. + */ + ~PolytopeCache() = default; + + void + clear() + { + // clear all the members + cell_face_at_boundary.clear(); + interface.clear(); + visited_cell_and_faces.clear(); + } + + /** + * Standard std::set for recording the standard cells and faces (in the + * deal.II lingo) that have been already visited. The first argument of + * the pair identifies the global index of a deal.II cell, while the + * second its local face number. + * + */ + mutable std::set> + visited_cell_and_faces; + + + mutable std::set> + visited_cell_and_faces_id; + + + + /** + * Map that associate the pair of (polytopal index, polytopal face) to + * (b,cell). The latter pair indicates whether or not the present face is + * on boundary. If it's on the boundary, then b is true and cell is an + * invalid cell iterator. Otherwise, b is false and cell points to the + * neighboring polytopal cell. + * + */ + mutable std::map< + std::pair, + std::pair::active_cell_iterator>> + cell_face_at_boundary; + + /** + * Map that associate the **local** pair of (polytope id, polytopal face) + * to the master id of the neighboring **ghosted** cell. + */ + mutable std::map, CellId> + ghosted_master_id; + + /** + * Standard std::map that associated to a pair of neighboring polytopic + * cells (current_polytope, neighboring_polytope) a sequence of + * ({deal_cell,deal_face_index}) which is meant to describe their + * interface. + * Indeed, the pair is identified by the two polytopic global indices, + * while the interface is described by a std::vector of deal.II cells and + * faces. + * + */ + mutable std::map< + std::pair, + std::vector< + std::pair::active_cell_iterator, + unsigned int>>> + interface; + }; + } // namespace internal +} // namespace dealii + + +/** + * + */ +template +class AgglomerationHandler : public Subscriptor +{ +public: + using agglomeration_iterator = AgglomerationIterator; + + using AgglomerationContainer = + typename AgglomerationIterator::AgglomerationContainer; + + + enum CellAgglomerationType + { + master = 0, + slave = 1 + }; + + + + explicit AgglomerationHandler( + const GridTools::Cache &cached_tria); + + AgglomerationHandler() = default; + + ~AgglomerationHandler() + { + // disconnect the signal + tria_listener.disconnect(); + } + + /** + * Iterator to the first polytope. + */ + agglomeration_iterator + begin() const; + + /** + * Iterator to the first polytope. + */ + agglomeration_iterator + begin(); + + /** + * Iterator to one past the last polygonal element. + */ + agglomeration_iterator + end() const; + + /** + * Iterator to one past the last polygonal element. + */ + agglomeration_iterator + end(); + + /** + * Iterator to the last polygonal element. + */ + agglomeration_iterator + last(); + + /** + * Returns an IteratorRange that makes up all the polygonal elements in the + * mesh. + */ + IteratorRange + polytope_iterators() const; + + template + friend class AgglomerationIterator; + + template + friend class AgglomerationAccessor; + + /** + * Distribute degrees of freedom on a grid where some cells have been + * agglomerated. + */ + void + distribute_agglomerated_dofs(const FiniteElement &fe_space); + + /** + * Overload for hp::FECollection. + */ + void + distribute_agglomerated_dofs( + const hp::FECollection &fe_collection_in); + + /** + * + * Set the degree of the quadrature formula to be used and the proper flags + * for the FEValues object on the agglomerated cell. + */ + void + initialize_fe_values( + const Quadrature &cell_quadrature = QGauss(1), + const UpdateFlags &flags = UpdateFlags::update_default, + const Quadrature &face_quadrature = QGauss(1), + const UpdateFlags &face_flags = UpdateFlags::update_default); + + /** + * Overload for hp::FECollection. + */ + void + initialize_fe_values( + const hp::QCollection &cell_qcollection = + hp::QCollection(QGauss(1)), + const UpdateFlags &flags = UpdateFlags::update_default, + const hp::QCollection &face_qcollection = + hp::QCollection(QGauss(1)), + const UpdateFlags &face_flags = UpdateFlags::update_default); + + /** + * Given a Triangulation with some agglomerated cells, create the sparsity + * pattern corresponding to a Discontinuous Galerkin discretization where the + * agglomerated cells are seen as one **unique** cell, with only the DoFs + * associated to the master cell of the agglomeration. + */ + template + void + create_agglomeration_sparsity_pattern( + SparsityPatternType &sparsity_pattern, + const AffineConstraints &constraints = AffineConstraints(), + const bool keep_constrained_dofs = true, + const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id); + + /** + * Store internally that the given cells are agglomerated. The convenction we + * take is the following: + * -1: cell is a master cell + * + * @note cells are assumed to be adjacent one to each other, and no check + * about this is done. + */ + agglomeration_iterator + define_agglomerate(const AgglomerationContainer &cells); + + /** + * Overload for hp::FECollection. + * + * The parameter @p fecollection_size provides the number of finite elements + * in the collection, allowing Polydeal to insert an empty element for + * slave cells internally. + * + * When @p fecollection_size equals 1, this function behaves identically to + * define_agglomerate(const AgglomerationContainer &cells). + */ + agglomeration_iterator + define_agglomerate(const AgglomerationContainer &cells, + const unsigned int fecollection_size); + + + inline const Triangulation & + get_triangulation() const; + + inline const FiniteElement & + get_fe() const; + + inline const Mapping & + get_mapping() const; + + inline const MappingBox & + get_agglomeration_mapping() const; + + inline const std::vector> & + get_local_bboxes() const; + + /** + * Return the mesh size of the polytopal mesh. It simply takes the maximum + * diameter over all the polytopes. + */ + double + get_mesh_size() const; + + inline types::global_cell_index + cell_to_polytope_index( + const typename Triangulation::active_cell_iterator &cell) + const; + + inline decltype(auto) + get_interface() const; + + /** + * Helper function to determine whether or not a cell is a master or a slave + */ + template + inline bool + is_master_cell(const CellIterator &cell) const; + + /** + * Find (if any) the cells that have the given master index. Note that `idx` + * is as it can be equal to -1 (meaning that the cell is a master one). + */ + inline const std::vector< + typename Triangulation::active_cell_iterator> & + get_slaves_of_idx(types::global_cell_index idx) const; + + + inline const LinearAlgebra::distributed::Vector & + get_relationships() const; + + /** + * + * @param master_cell + * @return std::vector< + * typename Triangulation::active_cell_iterator> + */ + inline std::vector< + typename Triangulation::active_cell_iterator> + get_agglomerate( + const typename Triangulation::active_cell_iterator + &master_cell) const; + + /** + * + * Return a constant reference to the DoFHandler underlying the + * agglomeration. It knows which cell have been agglomerated, and which FE + * spaces are present on each cell of the triangulation. + */ + inline const DoFHandler & + get_dof_handler() const; + + /** + * Returns the number of agglomerate cells in the grid. + */ + unsigned int + n_agglomerates() const; + + /** + * Return the number of agglomerated faces for a generic deal.II cell. + */ + unsigned int + n_agglomerated_faces_per_cell( + const typename Triangulation::active_cell_iterator &cell) + const; + + /** + * Construct a finite element space on the agglomeration. + */ + const FEValues & + reinit(const AgglomerationIterator &polytope) const; + + /** + * For a given polytope and face index, initialize shape functions, normals + * and quadratures rules to integrate there. + */ + const FEValuesBase & + reinit(const AgglomerationIterator &polytope, + const unsigned int face_index) const; + + /** + * + * Return a pair of FEValuesBase object reinited from the two sides of the + * agglomeration. + */ + std::pair &, + const FEValuesBase &> + reinit_interface(const AgglomerationIterator &polytope_in, + const AgglomerationIterator &neigh_polytope, + const unsigned int local_in, + const unsigned int local_outside) const; + + /** + * Return the agglomerated quadrature for the given agglomeration. This + * amounts to loop over all cells in an agglomeration and collecting together + * all the rules. + */ + Quadrature + agglomerated_quadrature( + const AgglomerationContainer &cells, + const typename Triangulation::active_cell_iterator + &master_cell) const; + + + /** + * + * This function generalizes the behaviour of cell->face(f)->at_boundary() + * in the case where f is an index out of the range [0,..., n_faces). + * In practice, if you call this function with a standard deal.II cell, you + * have precisely the same result as calling cell->face(f)->at_boundary(). + * Otherwise, if the cell is a master one, you have a boolean returning true + * is that face for the agglomeration is on the boundary or not. + */ + inline bool + at_boundary( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int f) const; + + inline unsigned int + n_dofs_per_cell() const noexcept; + + inline types::global_dof_index + n_dofs() const noexcept; + + + + /** + * Return the collection of vertices describing the boundary of the polytope + * associated to the master cell `cell`. The return type is meant to describe + * a sequence of edges (in 2D) or faces (in 3D). + */ + inline const std::vector::active_face_iterator> & + polytope_boundary( + const typename Triangulation::active_cell_iterator &cell); + + + /** + * DoFHandler for the agglomerated space + */ + DoFHandler agglo_dh; + + /** + * DoFHandler for the finest space: classical deal.II space + */ + DoFHandler output_dh; + + std::unique_ptr> box_mapping; + + /** + * This function stores the information needed to identify which polytopes are + * ghosted w.r.t the local partition. The issue this function addresses is due + * to the fact that the layer of ghost cells is made by just one layer of + * deal.II cells. Therefore, the neighboring polytopes will always be made by + * some ghost cells and **artificial** ones. This implies that we need to + * communicate the missing information from the neighboring rank. + */ + void + setup_ghost_polytopes(); + + void + exchange_interface_values(); + + // TODO: move it to private interface + mutable std::map< + types::subdomain_id, + std::map, std::vector>>> + recv_qpoints; + + mutable std::map< + types::subdomain_id, + std::map, std::vector>> + recv_jxws; + + mutable std::map< + types::subdomain_id, + std::map, std::vector>>> + recv_normals; + + mutable std::map< + types::subdomain_id, + std::map, std::vector>>> + recv_values; + + mutable std::map, + std::vector>>>> + recv_gradients; + + /** + * Given the index of a polytopic element, return a DoFHandler iterator + * for which DoFs associated to that polytope can be queried. + */ + inline const typename DoFHandler::active_cell_iterator + polytope_to_dh_iterator(const types::global_cell_index polytope_index) const; + + /** + * + */ + template + void + connect_hierarchy(const CellsAgglomerator &agglomerator); + + /** + * Return the finite element collection passed to + * distribute_agglomerated_dofs(). + */ + inline const hp::FECollection & + get_fe_collection() const; + + /** + * Return whether a hp::FECollection is being used. + */ + inline bool + used_fe_collection() const; + +private: + /** + * Initialize connectivity informations + */ + void + initialize_agglomeration_data( + const std::unique_ptr> &cache_tria); + + void + update_agglomerate( + AgglomerationContainer &polytope, + const typename Triangulation::active_cell_iterator + &master_cell); + + /** + * Reinitialize the agglomeration data. + */ + void + connect_to_tria_signals() + { + // First disconnect existing connections + tria_listener.disconnect(); + tria_listener = tria->signals.any_change.connect( + [&]() { this->initialize_agglomeration_data(this->cached_tria); }); + } + + /** + * Helper function to determine whether or not a cell is a slave cell. + * Instead of returning a boolean, it gives the index of the master cell. If + * it's a master cell, then the it returns -1, by construction. + */ + + inline typename Triangulation::active_cell_iterator & + is_slave_cell_of( + const typename Triangulation::active_cell_iterator &cell); + + /** + * Construct bounding boxes for an agglomeration described by a sequence of + * cells. This fills also the euler vector + */ + void + create_bounding_box(const AgglomerationContainer &polytope); + + + inline types::global_cell_index + get_master_idx_of_cell( + const typename Triangulation::active_cell_iterator &cell) + const; + + /** + * Returns true if the two given cells are agglomerated together. + */ + inline bool + are_cells_agglomerated( + const typename Triangulation::active_cell_iterator &cell, + const typename Triangulation::active_cell_iterator + &other_cell) const; + + /** + * Assign a finite element index on each cell of a triangulation, depending + * if it is a master cell, a slave cell, or a standard deal.II cell. A user + * doesn't need to know the internals of this, the only thing that is + * relevant is that after the call to the present function, DoFs are + * distributed in a different way if a cell is a master, slave, or standard + * cell. + */ + void + initialize_hp_structure(); + + + /** + * Helper function to call reinit on a master cell. + */ + const FEValuesBase & + reinit_master( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_number, + std::unique_ptr> + &agglo_isv_ptr) const; + + + /** + * Helper function to determine whether or not a cell is a slave cell, without + * any information about his parents. + */ + template + inline bool + is_slave_cell(const CellIterator &cell) const; + + + /** + * Initialize all the necessary connectivity information for an + * agglomeration. + */ + void + setup_connectivity_of_agglomeration(); + + + /** + * Record the number of agglomerations on the grid. + */ + unsigned int n_agglomerations; + + + /** + * Vector of indices such that v[cell->active_cell_index()] returns + * { -1 if `cell` is a master cell + * { `cell_master->active_cell_index()`, i.e. the index of the master cell if + * `cell` is a slave cell. + */ + LinearAlgebra::distributed::Vector master_slave_relationships; + + /** + * Same as the one above, but storing cell iterators rather than indices. + * + */ + std::map::active_cell_iterator> + master_slave_relationships_iterators; + + using ScratchData = MeshWorker::ScratchData; + + mutable std::vector number_of_agglomerated_faces; + + /** + * Associate a master cell (hence, a given polytope) to its boundary faces. + * The boundary is described through a vector of face iterators. + * + */ + mutable std::map< + const typename Triangulation::active_cell_iterator, + std::vector::active_face_iterator>> + polygon_boundary; + + + /** + * Vector of `BoundingBoxes` s.t. `bboxes[idx]` equals BBOx associated to the + * agglomeration with master cell indexed by ìdx`. Othwerwise default BBox is + * empty + * + */ + std::vector> bboxes; + + //////////////////////////////////////////////////////// + + + // n_faces + mutable std::map> + local_n_faces; + + mutable std::map> + recv_n_faces; + + + // CellId (including slaves) + mutable std::map> + local_cell_ids_neigh_cell; + + mutable std::map> + recv_cell_ids_neigh_cell; + + + // send to neighborign rank the information that + // - current polytope id + // - face f + // has the following neighboring id. + mutable std::map>> + local_ghosted_master_id; + + mutable std::map>> + recv_ghosted_master_id; + + // CellIds from neighboring rank + mutable std::map>> + local_bdary_info; + + mutable std::map>> + recv_bdary_info; + + // Exchange neighboring bounding boxes + mutable std::map>> + local_ghosted_bbox; + + mutable std::map>> + recv_ghosted_bbox; + + // Exchange DoF indices with ghosted polytopes + mutable std::map>> + local_ghost_dofs; + + mutable std::map>> + recv_ghost_dofs; + + // Exchange qpoints + mutable std::map< + types::subdomain_id, + std::map, std::vector>>> + local_qpoints; + + // Exchange jxws + mutable std::map< + types::subdomain_id, + std::map, std::vector>> + local_jxws; + + // Exchange normals + mutable std::map< + types::subdomain_id, + std::map, std::vector>>> + local_normals; + + // Exchange values + mutable std::map< + types::subdomain_id, + std::map, std::vector>>> + local_values; + + mutable std::map, + std::vector>>>> + local_gradients; + + + + //////////////////////////////////////////////////////// + + const Triangulation *tria; + + const Mapping *mapping; + + std::unique_ptr> cached_tria; + + const MPI_Comm communicator; + + // The FiniteElement space we have on each cell. Currently supported types are + // FE_DGQ and FE_DGP elements. + std::unique_ptr> fe; + + hp::FECollection fe_collection; + + /** + * Eulerian vector describing the new cells obtained by the bounding boxes + */ + LinearAlgebra::distributed::Vector euler_vector; + + + /** + * Use this in reinit(cell) for (non-agglomerated, standard) cells, + * and return the result of scratch.reinit(cell) for cells + */ + mutable std::unique_ptr standard_scratch; + + /** + * Fill this up in reinit(cell), for agglomerated cells, using the custom + * quadrature, and return the result of + * scratch.reinit(cell); + */ + mutable std::unique_ptr agglomerated_scratch; + + + mutable std::unique_ptr> + agglomerated_isv; + + mutable std::unique_ptr> + agglomerated_isv_neigh; + + mutable std::unique_ptr> + agglomerated_isv_bdary; + + boost::signals2::connection tria_listener; + + UpdateFlags agglomeration_flags; + + const UpdateFlags internal_agglomeration_flags = + update_values | update_gradients | update_JxW_values | + update_quadrature_points; + + UpdateFlags agglomeration_face_flags; + + const UpdateFlags internal_agglomeration_face_flags = + update_quadrature_points | update_normal_vectors | update_values | + update_gradients | update_JxW_values | update_inverse_jacobians; + + Quadrature agglomeration_quad; + + Quadrature agglomeration_face_quad; + + // Associate the master cell to the slaves. + std::unordered_map< + types::global_cell_index, + std::vector::active_cell_iterator>> + master2slaves; + + // Map the master cell index with the polytope index + std::map master2polygon; + + + std::vector::active_cell_iterator> + master_disconnected; + + // Dummy FiniteElement objects needed only to generate quadratures + + /** + * Dummy FE_Nothing + */ + FE_Nothing dummy_fe; + + /** + * Dummy FEValues, needed for cell quadratures. + */ + std::unique_ptr> no_values; + + /** + * Dummy FEFaceValues, needed for face quadratures. + */ + std::unique_ptr> no_face_values; + + /** + * A contiguous container for all of the master cells. + */ + std::vector::active_cell_iterator> + master_cells_container; + + friend class internal::AgglomerationHandlerImplementation; + + internal::PolytopeCache polytope_cache; + + /** + * Bool that keeps track whether the mesh is composed also by standard deal.II + * cells as (trivial) polytopes. + */ + bool hybrid_mesh; + + std::map, + std::vector> + parent_child_info; + + unsigned int present_extraction_level; + + // Support for hp::FECollection + bool is_hp_collection = false; // Indicates whether hp::FECollection is used + std::unique_ptr> + hp_fe_collection; // External input FECollection + + // Stores quadrature rules; these QCollections should have the same size as + // hp_fe_collection + hp::QCollection agglomeration_quad_collection; + hp::QCollection agglomeration_face_quad_collection; + + hp::MappingCollection + mapping_collection; // Contains only one mapping object + hp::FECollection + dummy_fe_collection; // Similar to dummy_fe, but as an FECollection + // 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 + std::unique_ptr> hp_no_values; + std::unique_ptr> hp_no_face_values; +}; + + + +// ------------------------------ inline functions ------------------------- +template +inline const FiniteElement & +AgglomerationHandler::get_fe() const +{ + return *fe; +} + + + +template +inline const Mapping & +AgglomerationHandler::get_mapping() const +{ + return *mapping; +} + + + +template +inline const MappingBox & +AgglomerationHandler::get_agglomeration_mapping() const +{ + return *box_mapping; +} + + + +template +inline const Triangulation & +AgglomerationHandler::get_triangulation() const +{ + return *tria; +} + + +template +inline const std::vector> & +AgglomerationHandler::get_local_bboxes() const +{ + return bboxes; +} + + + +template +inline types::global_cell_index +AgglomerationHandler::cell_to_polytope_index( + const typename Triangulation::active_cell_iterator &cell) const +{ + return master2polygon.at(cell->active_cell_index()); +} + + + +template +inline decltype(auto) +AgglomerationHandler::get_interface() const +{ + return polytope_cache.interface; +} + + + +template +inline const LinearAlgebra::distributed::Vector & +AgglomerationHandler::get_relationships() const +{ + return master_slave_relationships; +} + + + +template +inline std::vector::active_cell_iterator> +AgglomerationHandler::get_agglomerate( + const typename Triangulation::active_cell_iterator + &master_cell) const +{ + Assert(is_master_cell(master_cell), ExcInternalError()); + auto agglomeration = get_slaves_of_idx(master_cell->active_cell_index()); + agglomeration.push_back(master_cell); + return agglomeration; +} + + + +template +inline const DoFHandler & +AgglomerationHandler::get_dof_handler() const +{ + return agglo_dh; +} + + + +template +inline const std::vector< + typename Triangulation::active_cell_iterator> & +AgglomerationHandler::get_slaves_of_idx( + types::global_cell_index idx) const +{ + return master2slaves.at(idx); +} + + + +template +template +inline bool +AgglomerationHandler::is_master_cell( + const CellIterator &cell) const +{ + return master_slave_relationships[cell->global_active_cell_index()] == -1; +} + + + +/** + * Helper function to determine whether or not a cell is a slave cell, without + * any information about his parents. + */ +template +template +inline bool +AgglomerationHandler::is_slave_cell( + const CellIterator &cell) const +{ + return master_slave_relationships[cell->global_active_cell_index()] >= 0; +} + + + +template +inline bool +AgglomerationHandler::at_boundary( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_index) const +{ + Assert(!is_slave_cell(cell), + ExcMessage("This function should not be called for a slave cell.")); + + return polytope_cache.cell_face_at_boundary + .at({master2polygon.at(cell->active_cell_index()), face_index}) + .first; +} + + +template +inline unsigned int +AgglomerationHandler::n_dofs_per_cell() const noexcept +{ + return fe->n_dofs_per_cell(); +} + + + +template +inline types::global_dof_index +AgglomerationHandler::n_dofs() const noexcept +{ + return agglo_dh.n_dofs(); +} + + + +template +inline const std::vector::active_face_iterator> & +AgglomerationHandler::polytope_boundary( + const typename Triangulation::active_cell_iterator &cell) +{ + return polygon_boundary[cell]; +} + + + +template +inline typename Triangulation::active_cell_iterator & +AgglomerationHandler::is_slave_cell_of( + const typename Triangulation::active_cell_iterator &cell) +{ + return master_slave_relationships_iterators.at(cell->active_cell_index()); +} + + + +template +inline types::global_cell_index +AgglomerationHandler::get_master_idx_of_cell( + const typename Triangulation::active_cell_iterator &cell) const +{ + auto idx = master_slave_relationships[cell->global_active_cell_index()]; + if (idx == -1) + return cell->global_active_cell_index(); + else + return static_cast(idx); +} + + + +template +inline bool +AgglomerationHandler::are_cells_agglomerated( + const typename Triangulation::active_cell_iterator &cell, + const typename Triangulation::active_cell_iterator &other_cell) + const +{ + // if different subdomain, then **by construction** they will not be together + // if (cell->subdomain_id() != other_cell->subdomain_id()) + // return false; + // else + return (get_master_idx_of_cell(cell) == get_master_idx_of_cell(other_cell)); +} + + + +template +inline unsigned int +AgglomerationHandler::n_agglomerates() const +{ + return n_agglomerations; +} + + + +template +inline const typename DoFHandler::active_cell_iterator +AgglomerationHandler::polytope_to_dh_iterator( + const types::global_cell_index polytope_index) const +{ + return master_cells_container[polytope_index]->as_dof_handler_iterator( + agglo_dh); +} + + + +template +AgglomerationIterator +AgglomerationHandler::begin() const +{ + Assert(n_agglomerations > 0, + ExcMessage("No agglomeration has been performed.")); + return {*master_cells_container.begin(), this}; +} + + + +template +AgglomerationIterator +AgglomerationHandler::begin() +{ + Assert(n_agglomerations > 0, + ExcMessage("No agglomeration has been performed.")); + return {*master_cells_container.begin(), this}; +} + + + +template +AgglomerationIterator +AgglomerationHandler::end() const +{ + Assert(n_agglomerations > 0, + ExcMessage("No agglomeration has been performed.")); + return {*master_cells_container.end(), this}; +} + + + +template +AgglomerationIterator +AgglomerationHandler::end() +{ + Assert(n_agglomerations > 0, + ExcMessage("No agglomeration has been performed.")); + return {*master_cells_container.end(), this}; +} + + + +template +AgglomerationIterator +AgglomerationHandler::last() +{ + Assert(n_agglomerations > 0, + ExcMessage("No agglomeration has been performed.")); + return {master_cells_container.back(), this}; +} + + + +template +IteratorRange< + typename AgglomerationHandler::agglomeration_iterator> +AgglomerationHandler::polytope_iterators() const +{ + return IteratorRange< + typename AgglomerationHandler::agglomeration_iterator>( + begin(), end()); +} + +template +template +void +AgglomerationHandler::connect_hierarchy( + const CellsAgglomerator &agglomerator) +{ + parent_child_info = agglomerator.parent_node_to_children_nodes; + present_extraction_level = agglomerator.extraction_level; +} + +template +inline const hp::FECollection & +AgglomerationHandler::get_fe_collection() const +{ + return *hp_fe_collection; +} + +template +inline bool +AgglomerationHandler::used_fe_collection() const +{ + return is_hp_collection; +} + + +#endif diff --git a/agglomeration_poisson/include/agglomeration_iterator.h b/agglomeration_poisson/include/agglomeration_iterator.h new file mode 100644 index 00000000..878e74f8 --- /dev/null +++ b/agglomeration_poisson/include/agglomeration_iterator.h @@ -0,0 +1,304 @@ +/* ----------------------------------------------------------------------------- + * + * 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 agglomeration_iterator_h +#define agglomeration_iterator_h + + +#include + + +/** + * A class that is used to iterate over polygons. Together with the + * AgglomerationAccessor class this is used to hide the internal implementation + * of the particle class and the particle container. + */ +template +class AgglomerationIterator +{ +public: + using AgglomerationContainer = + typename AgglomerationAccessor::AgglomerationContainer; + + /** + * Empty constructor. This constructore creates an iterator pointing to an + * invalid object. + */ + AgglomerationIterator(); + + /** + * Constructor of the iterator. Takes a reference to the master cell encoding + * the actual polytope. + */ + AgglomerationIterator( + const typename Triangulation::active_cell_iterator &cell, + const AgglomerationHandler *handler); + + /** + * Same as above, needed for ghosted elements. + */ + AgglomerationIterator( + const typename Triangulation::active_cell_iterator + &master_cell, + const CellId &cell_id, + const AgglomerationHandler *handler); + + /** + * Dereferencing operator, returns a reference to an accessor. Usage is thus + * like (*i).get_dof_indices (); + */ + const AgglomerationAccessor & + operator*() const; + + /** + * Dereferencing operator, non-@p const version. + */ + AgglomerationAccessor & + operator*(); + + /** + * Dereferencing operator, returns a pointer of the particle pointed to. + * Usage is thus like i->get_dof_indices (); + * + * There is a @p const and a non-@p const version. + */ + const AgglomerationAccessor * + operator->() const; + + /** + * Dereferencing operator, non-@p const version. + */ + AgglomerationAccessor * + operator->(); + + /** + * Compare for equality. + */ + bool + operator==(const AgglomerationIterator &) const; + + /** + * Compare for inequality. + */ + bool + operator!=(const AgglomerationIterator &) const; + + /** + * Prefix ++ operator: ++iterator. This operator advances + * the iterator to the next element and returns a reference to + * *this. + */ + AgglomerationIterator & + operator++(); + + /** + * Postfix ++ operator: iterator++. This operator advances + * the iterator to the next element, but returns an iterator to the element + * previously pointed to. + */ + AgglomerationIterator + operator++(int); + + /** + * Prefix \-- operator: \--iterator. This operator moves + * the iterator to the previous element and returns a reference to + * *this. + */ + AgglomerationIterator & + operator--(); + + /** + * Postfix \-- operator: iterator\--. This operator moves + * the iterator to the previous element, but returns an iterator to the + * element previously pointed to. + */ + AgglomerationIterator + operator--(int); + + /** + * Return the state of the present iterator. + */ + IteratorState::IteratorStates + state() const; + + /** + * Return the master cell associated to the present polytope. + */ + const typename Triangulation::active_cell_iterator & + master_cell() const; + + /** + * Mark the class as bidirectional iterator and declare some alias which + * are standard for iterators and are used by algorithms to enquire about + * the specifics of the iterators they work on. + */ + using iterator_category = std::bidirectional_iterator_tag; + using value_type = AgglomerationAccessor; + using difference_type = std::ptrdiff_t; + using pointer = AgglomerationAccessor *; + using reference = AgglomerationAccessor &; + +private: + /** + * The accessor to the actual polytope. + */ + AgglomerationAccessor accessor; +}; + + + +// ------------------------------ inline functions ------------------------- + +template +inline AgglomerationIterator::AgglomerationIterator() + : accessor() +{} + + + +template +inline AgglomerationIterator::AgglomerationIterator( + const typename Triangulation::active_cell_iterator + &master_cell, + const AgglomerationHandler *handler) + : accessor(master_cell, handler) +{} + +template +inline AgglomerationIterator::AgglomerationIterator( + const typename Triangulation::active_cell_iterator + &master_cell, + const CellId &cell_id, + const AgglomerationHandler *handler) + : accessor(master_cell, cell_id, handler) +{} + + + +template +inline AgglomerationAccessor & +AgglomerationIterator::operator*() +{ + return accessor; +} + + + +template +inline AgglomerationAccessor * +AgglomerationIterator::operator->() +{ + return &(this->operator*()); +} + + + +template +inline const AgglomerationAccessor & +AgglomerationIterator::operator*() const +{ + return accessor; +} + + + +template +inline const AgglomerationAccessor * +AgglomerationIterator::operator->() const +{ + return &(this->operator*()); +} + + + +template +inline bool +AgglomerationIterator::operator!=( + const AgglomerationIterator &other) const +{ + return accessor != other.accessor; +} + + + +template +inline bool +AgglomerationIterator::operator==( + const AgglomerationIterator &other) const +{ + return accessor == other.accessor; +} + + + +template +inline AgglomerationIterator & +AgglomerationIterator::operator++() +{ + accessor.next(); + return *this; +} + + + +template +inline AgglomerationIterator +AgglomerationIterator::operator++(int) +{ + AgglomerationIterator tmp(*this); + operator++(); + + return tmp; +} + + + +template +inline AgglomerationIterator & +AgglomerationIterator::operator--() +{ + accessor.prev(); + return *this; +} + + + +template +inline AgglomerationIterator +AgglomerationIterator::operator--(int) +{ + AgglomerationIterator tmp(*this); + operator--(); + + return tmp; +} + + + +template +inline IteratorState::IteratorStates +AgglomerationIterator::state() const +{ + return accessor.master_cell.state(); +} + + + +template +inline const typename Triangulation::active_cell_iterator & +AgglomerationIterator::master_cell() const +{ + return accessor.master_cell; +} + + + +#endif diff --git a/agglomeration_poisson/include/agglomerator.h b/agglomeration_poisson/include/agglomerator.h new file mode 100644 index 00000000..79a4d5f4 --- /dev/null +++ b/agglomeration_poisson/include/agglomerator.h @@ -0,0 +1,471 @@ +/* ----------------------------------------------------------------------------- + * + * 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 agglomerator_h +#define agglomerator_h + + +#include + +#include + +#include +#include +#include + +template +class AgglomerationHandler; + +namespace dealii +{ + namespace internal + { + template + 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, + const unsigned int target_level, + std::vector::value>::active_cell_iterator>> + &agglomerates_, + std::vector &n_nodes_per_level, + std::map, + std::vector> &parent_to_children); + + /** + * 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; + + /** + * 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 the mesh + * iterators composing the i-th agglomerate. + */ + std::vector::value>::active_cell_iterator>> + &agglomerates; + + /** + * Store the total number of nodes on each level. + */ + std::vector &n_nodes_per_level; + + /** + * Map that associates to a given node on level l its children, identified + * by their integer index. + */ + std::map, + std::vector> + &parent_node_to_children_nodes; + }; + + + + template + Rtree_visitor::Rtree_visitor( + const Translator &translator, + const unsigned int target_level, + std::vector::value>::active_cell_iterator>> + &agglomerates_, + std::vector &n_nodes_per_level_, + std::map, + std::vector> &parent_to_children) + : translator(translator) + , level(0) + , node_counter(0) + , target_level(target_level) + , agglomerates(agglomerates_) + , n_nodes_per_level(n_nodes_per_level_) + , parent_node_to_children_nodes(parent_to_children) + {} + + + + template + void + Rtree_visitor::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 + // node + 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) + { + const auto offset = agglomerates.size(); + agglomerates.resize(offset + 1); + size_t level_backup = level; + + ++level; + for (const auto &entry : elements) + { + boost::geometry::index::detail::rtree::apply_visitor( + *this, *entry.second); + } + // Done with node number 'node_counter' on level target_level. + + ++node_counter; // visited all children of an internal node + n_nodes_per_level[target_level]++; + + level = level_backup; + } + else if (level > target_level) + { + // I am on a child (internal) node on a deeper level. + + // Keep visiting until you go to the leafs. + size_t level_backup = level; + + ++level; + + // looping through entries of node + for (const auto &entry : elements) + { + boost::geometry::index::detail::rtree::apply_visitor( + *this, *entry.second); + } + // done with node on level l > target_level (not just + // "target_level+1). + 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 + void + Rtree_visitor::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) + { + // If I want to extract from leaf node, i.e. the target_level is the + // last one where leafs are grouped together. + 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 + class CellsAgglomerator + { + public: + template + 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::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, + std::vector> & + 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::active_cell_iterator>> + agglomerates_on_level; + + /** + * Vector storing the number of nodes (and, ultimately, agglomerates) for + * each level. + */ + std::vector 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::vector> + parent_node_to_children_nodes; + }; + + + + template + CellsAgglomerator::CellsAgglomerator( + const RtreeType &tree, + const unsigned int extraction_level_) + : extraction_level(extraction_level_) + { + rtree = const_cast(&tree); + Assert(n_levels(*rtree), ExcMessage("At least two levels are needed.")); + } + + + + template + const std::vector< + std::vector::active_cell_iterator>> & + CellsAgglomerator::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; + RtreeView rtv(*rtree); + + n_nodes_per_level.resize(rtv.depth() + + 1); // store how many nodes we have for each level. + + if (rtv.depth() == 0) + { + // The below algorithm does not work for `rtv.depth()==0`, which might + // happen if the number entries in the tree is too small. + agglomerates_on_level.resize(1); + agglomerates_on_level[0].resize(1); + } + else + { + const unsigned int target_level = + std::min(extraction_level, rtv.depth()); + + internal::Rtree_visitor + 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; + } + + + + // ------------------------------ inline functions ------------------------- + + + template + inline unsigned int + CellsAgglomerator::get_n_levels() const + { + return n_levels(*rtree); + } + + + + template + inline types::global_cell_index + CellsAgglomerator::get_n_nodes_per_level( + const unsigned int level) const + { + return n_nodes_per_level[level]; + } + + + + template + inline const std::map< + std::pair, + std::vector> & + CellsAgglomerator::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 diff --git a/agglomeration_poisson/include/mapping_box.h b/agglomeration_poisson/include/mapping_box.h new file mode 100644 index 00000000..6183e5c1 --- /dev/null +++ b/agglomeration_poisson/include/mapping_box.h @@ -0,0 +1,382 @@ +/* ----------------------------------------------------------------------------- + * + * 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 + +#include +#include + +#include + +#include + + +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 +class MappingBox : public Mapping +{ +public: + MappingBox(const std::vector> &local_boxes, + const std::map + &polytope_translator); + // for documentation, see the Mapping base class + virtual std::unique_ptr> + clone() const override; + + /** + * Return @p true because MappingBox preserves vertex + * locations. + */ + virtual bool + preserves_vertex_locations() const override; + + virtual bool + is_compatible_with(const ReferenceCell &reference_cell) const override; + + /** + * @name Mapping points between reference and real cells + * @{ + */ + + // for documentation, see the Mapping base class + virtual Point + transform_unit_to_real_cell( + const typename Triangulation::cell_iterator &cell, + const Point &p) const override; + + // for documentation, see the Mapping base class + virtual Point + transform_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const Point &p) const override; + + // for documentation, see the Mapping base class + virtual void + transform_points_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const ArrayView> &real_points, + const ArrayView> &unit_points) const override; + + /** + * @} + */ + + /** + * @name Functions to transform tensors from reference to real coordinates + * @{ + */ + + // for documentation, see the Mapping base class + virtual void + transform(const ArrayView> &input, + const MappingKind kind, + const typename Mapping::InternalDataBase &internal, + const ArrayView> &output) const override; + + // for documentation, see the Mapping base class + virtual void + transform(const ArrayView> &input, + const MappingKind kind, + const typename Mapping::InternalDataBase &internal, + const ArrayView> &output) const override; + + // for documentation, see the Mapping base class + virtual void + transform(const ArrayView> &input, + const MappingKind kind, + const typename Mapping::InternalDataBase &internal, + const ArrayView> &output) const override; + + // for documentation, see the Mapping base class + virtual void + transform(const ArrayView> &input, + const MappingKind kind, + const typename Mapping::InternalDataBase &internal, + const ArrayView> &output) const override; + + // for documentation, see the Mapping base class + virtual void + transform(const ArrayView> &input, + const MappingKind kind, + const typename Mapping::InternalDataBase &internal, + const ArrayView> &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::InternalDataBase + { + public: + /** + * Default constructor. + */ + InternalData() = default; + + /** + * Constructor that initializes the object with a quadrature. + */ + InternalData(const Quadrature &quadrature); + + // Documentation see Mapping::InternalDataBase. + virtual void + reinit(const UpdateFlags update_flags, + const Quadrature &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., hx, hy, hz. + */ + 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., hx, + * hy, hz. + */ + 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> quadrature_points; + }; + +private: + // documentation can be found in Mapping::requires_update_flags() + virtual UpdateFlags + requires_update_flags(const UpdateFlags update_flags) const override; + + // documentation can be found in Mapping::get_data() + virtual std::unique_ptr::InternalDataBase> + get_data(const UpdateFlags, const Quadrature &quadrature) const override; + + using Mapping::get_face_data; + + // documentation can be found in Mapping::get_subface_data() + virtual std::unique_ptr::InternalDataBase> + get_subface_data(const UpdateFlags flags, + const Quadrature &quadrature) const override; + + // documentation can be found in Mapping::fill_fe_values() + virtual CellSimilarity::Similarity + fill_fe_values( + const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const override; + + using Mapping::fill_fe_face_values; + + // documentation can be found in Mapping::fill_fe_subface_values() + virtual void + fill_fe_subface_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const override; + + // documentation can be found in Mapping::fill_fe_immersed_surface_values() + virtual void + fill_fe_immersed_surface_values( + const typename Triangulation::cell_iterator &cell, + const NonMatching::ImmersedSurfaceQuadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &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::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::cell_iterator &cell, + const InternalData &data, + const ArrayView> &unit_quadrature_points, + std::vector> &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> &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 + &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 + &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 + &output_data) const; + + /** + * Vector of (local) bounding boxes + */ + std::vector> boxes; + + /** + * Map from global cell index to bounding box index + */ + std::map + polytope_translator; +}; + +/** @} */ + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/agglomeration_poisson/include/poly_utils.h b/agglomeration_poisson/include/poly_utils.h new file mode 100644 index 00000000..f34fe23e --- /dev/null +++ b/agglomeration_poisson/include/poly_utils.h @@ -0,0 +1,1667 @@ +/* ----------------------------------------------------------------------------- + * + * 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 + +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include + +namespace dealii::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 &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 + void + get_face_connectivity_of_cells( + const parallel::fullydistributed::Triangulation + &triangulation, + DynamicSparsityPattern &cell_connectivity, + const std::vector locally_owned_cells) + { + cell_connectivity.reinit(triangulation.n_locally_owned_active_cells(), + triangulation.n_locally_owned_active_cells()); + + // 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 + 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 dealii::PolyUtils::internal + +namespace dealii::PolyUtils +{ + template + 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::value>::active_cell_iterator>> &boxes, + std::vector> &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::value>::active_cell_iterator>> + &agglomerates; + + std::vector> &row_ptr; + }; + + template + Rtree_visitor::Rtree_visitor( + const Translator &translator, + const unsigned int target_level, + std::vector::value>::active_cell_iterator>> + &bb_in_boxes, + std::vector> &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 + void + Rtree_visitor::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 + // node + 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) + { + // const unsigned int n_children = elements.size(); + 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); + } + // Done with node number 'node_counter' + + ++node_counter; // visited all children of an internal node + + level = level_backup; + } + else if (level > target_level) + { + // Keep visiting until you go to the leafs. + 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 + void + Rtree_visitor::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 + 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 + void + interpolate_to_fine_grid( + const AgglomerationHandler &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; + using MatrixType = std::conditional_t>; + + MatrixType interpolation_matrix; + + [[maybe_unused]] + typename std::conditional_t + sp; + + // Get some info from the handler + const DoFHandler &agglo_dh = + agglomeration_handler.agglo_dh; + + DoFHandler *output_dh = + const_cast *>( + &agglomeration_handler.output_dh); + const FiniteElement &fe = agglomeration_handler.get_fe(); + const Mapping &mapping = agglomeration_handler.get_mapping(); + const Triangulation &tria = + agglomeration_handler.get_triangulation(); + const auto &bboxes = agglomeration_handler.get_local_bboxes(); + + std::unique_ptr> output_fe; + if (tria.all_reference_cells_are_hyper_cube()) + output_fe = std::make_unique>(fe.degree); + else if (tria.all_reference_cells_are_simplex()) + output_fe = std::make_unique>(fe.degree); + else + AssertThrow(false, ExcNotImplemented()); + + // Setup an auxiliary DoFHandler for output purposes + 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 agglo_dof_indices(fe.dofs_per_cell); + std::vector standard_dof_indices( + fe.dofs_per_cell); + std::vector output_dof_indices( + output_fe->dofs_per_cell); + + Quadrature quad(output_fe->get_unit_support_points()); + FEValues 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) + { + // addd master-slave relationship + 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 local_matrix(fe.dofs_per_cell, fe.dofs_per_cell); + std::vector> reference_q_points(fe.dofs_per_cell); + + // Dummy AffineConstraints, only needed for loc2glb + AffineConstraints 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); + + // Get the box of this agglomerate. + const BoundingBox &box = bboxes[polytope_index]; + + for (const auto &slave : slaves) + { + // add master-slave relationship + 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) + { + const MPI_Comm &communicator = tria.get_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>) + { + sp.copy_from(dsp); + interpolation_matrix.reinit(sp); + dst.reinit(output_dh->n_dofs()); + assemble_interpolation_matrix(); + } + else + { + // PETSc, LA::d::v options not implemented. + (void)agglomeration_handler; + (void)dst; + (void)src; + AssertThrow(false, ExcNotImplemented()); + } + + // If tria is distributed + if (dynamic_cast *>( + &tria) != nullptr) + interpolation_matrix.compress(VectorOperation::add); + + // Finally, perform the interpolation. + 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 + void + interpolate_to_fine_grid( + const AgglomerationHandler &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; + + static constexpr bool is_supported_vector = + std::is_same_v> || is_trilinos_vector; + static_assert(is_supported_vector); + + // First, check for an easy return + if (on_the_fly == false) + { + return internal::interpolate_to_fine_grid(agglomeration_handler, + dst, + src); + } + else + { + // otherwise, do not create any matrix + if (!agglomeration_handler.used_fe_collection()) + { + // Original version: handle case without hp::FECollection + const Triangulation &tria = + agglomeration_handler.get_triangulation(); + const Mapping &mapping = agglomeration_handler.get_mapping(); + const FiniteElement &original_fe = + agglomeration_handler.get_fe(); + + // 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. + std::unique_ptr> output_fe; + if (tria.all_reference_cells_are_hyper_cube()) + output_fe = std::make_unique>(original_fe.degree); + else if (tria.all_reference_cells_are_simplex()) + output_fe = + std::make_unique>(original_fe.degree); + else + AssertThrow(false, ExcNotImplemented()); + + DoFHandler &output_dh = + const_cast &>(agglomeration_handler.output_dh); + output_dh.reinit(tria); + output_dh.distribute_dofs(*output_fe); + + if constexpr (std::is_same_v) + { + const IndexSet &locally_owned_dofs = + output_dh.locally_owned_dofs(); + dst.reinit(locally_owned_dofs); + } + else if constexpr (std::is_same_v>) + { + dst.reinit(output_dh.n_dofs()); + } + else + { + // PETSc, LA::d::v options not implemented. + (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 quad(output_fe->get_unit_support_points()); + FEValues output_fe_values(mapping, + *output_fe, + quad, + update_quadrature_points); + + std::vector local_dof_indices( + dofs_per_cell); + std::vector 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 &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 + { + // Handle the hp::FECollection case + const Triangulation &tria = + agglomeration_handler.get_triangulation(); + const Mapping &mapping = agglomeration_handler.get_mapping(); + const hp::FECollection &original_fe_collection = + agglomeration_handler.get_fe_collection(); + + // 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. + hp::FECollection 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) + { + // Scalar case + for (unsigned int i = 0; i < original_fe_collection.size(); ++i) + { + std::unique_ptr> output_fe; + if (tria.all_reference_cells_are_hyper_cube()) + output_fe = std::make_unique>( + original_fe_collection[i].degree); + else if (tria.all_reference_cells_are_simplex()) + output_fe = std::make_unique>( + 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) + { + // System case + for (unsigned int i = 0; i < original_fe_collection.size(); ++i) + { + std::vector *> + base_elements; + std::vector multiplicities; + for (unsigned int b = 0; + b < original_fe_collection[i].n_base_elements(); + ++b) + { + if (dynamic_cast *>( + &original_fe_collection[i].base_element(b))) + base_elements.push_back( + new FE_Nothing()); + else + { + if (tria.all_reference_cells_are_hyper_cube()) + base_elements.push_back(new FE_DGQ( + original_fe_collection[i] + .base_element(b) + .degree)); + else if (tria.all_reference_cells_are_simplex()) + base_elements.push_back( + new FE_SimplexDGP( + original_fe_collection[i] + .base_element(b) + .degree)); + else + AssertThrow(false, ExcNotImplemented()); + } + multiplicities.push_back( + original_fe_collection[i].element_multiplicity(b)); + } + + FESystem output_fe_system(base_elements, + multiplicities); + for (const auto *ptr : base_elements) + delete ptr; + output_fe_collection.push_back(output_fe_system); + } + } + + DoFHandler &output_dh = + const_cast &>(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::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) + { + const IndexSet &locally_owned_dofs = + output_dh.locally_owned_dofs(); + dst.reinit(locally_owned_dofs); + } + else if constexpr (std::is_same_v>) + { + dst.reinit(output_dh.n_dofs()); + } + else + { + // PETSc, LA::d::v options not implemented. + (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 quad(output_fe_collection[active_fe_idx] + .get_unit_support_points()); + FEValues output_fe_values( + mapping, + output_fe_collection[active_fe_idx], + quad, + update_quadrature_points); + std::vector local_dof_indices( + dofs_per_cell); + std::vector + local_dof_indices_output(output_dofs_per_cell); + + polytope->get_dof_indices(local_dof_indices); + const BoundingBox &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 + void + compute_global_error(const AgglomerationHandler &agglomeration_handler, + const VectorType &solution, + const Function &exact_solution, + const std::vector &norms, + std::vector &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.); + + // Vector storing errors local to the current processor. + std::vector local_errors(norms.size()); + std::fill(local_errors.begin(), local_errors.end(), 0.); + + // Get some info from the handler + 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 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 analyical_sol_at_qpoints(n_qpoints); + exact_solution.value_list(q_points, analyical_sol_at_qpoints); + std::vector> 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); + } + // L2 + local_errors[0] += std::pow((analyical_sol_at_qpoints[q_index] - + solution_at_qpoint), + 2) * + agglo_values.JxW(q_index); + + // H1 seminorm + 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); + } + } + } + + // Perform reduction and take sqrt of each error + global_errors[0] = Utilities::MPI::reduce( + local_errors[0], + agglomeration_handler.get_triangulation().get_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( + local_errors[1], + agglomeration_handler.get_triangulation().get_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 + unsigned int + construct_agglomerated_levels( + const Triangulation &tria, + std::vector>> + &agglomeration_handlers, + const FE_DGQ &fe_dg, + const Mapping &mapping, + const unsigned int starting_tree_level) + { + const auto parallel_tria = + dynamic_cast *>(&tria); + + GridTools::Cache cached_tria(tria); + Assert(parallel_tria->n_active_cells() > 0, ExcInternalError()); + + const MPI_Comm comm = parallel_tria->get_communicator(); + ConditionalOStream pcout(std::cout, + (Utilities::MPI::this_mpi_process(comm) == 0)); + + // Start building R-tree + namespace bgi = boost::geometry::index; + static constexpr unsigned int max_elem_per_node = + constexpr_pow(2, dim); // 2^dim + std::vector, + typename Triangulation::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>(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; + + // Resize the agglomeration handlers to the right size + + agglomeration_handlers.resize(total_tree_levels); + // Loop through the available levels and set AgglomerationHandlers up. + 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>(cached_tria); + CellsAgglomerator agglomerator{tree, + extraction_level}; + const auto agglomerates = agglomerator.extract_agglomerates(); + agglomeration_handlers[extraction_level - starting_tree_level] + ->connect_hierarchy(agglomerator); + + // Flag elements for agglomeration + 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; + + // Now, perform agglomeration within each locally owned partition + std::vector< + std::vector::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); + + // For every subdomain, agglomerate elements together + 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(fe_dg.degree + 1), + update_values | update_gradients | + update_JxW_values | update_quadrature_points, + QGauss(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 + void + assemble_local_jumps_and_averages(FullMatrix &M11, + FullMatrix &M12, + FullMatrix &M21, + FullMatrix &M22, + const FEValuesBase &fe_faces0, + const FEValuesBase &fe_faces1, + const double penalty_constant, + const double h_f) + { + const std::vector> &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 + void + assemble_local_jumps_and_averages_ghost( + FullMatrix &M11, + FullMatrix &M12, + FullMatrix &M21, + FullMatrix &M22, + const FEValuesBase &fe_faces0, + const std::vector> &recv_values, + const std::vector>> &recv_gradients, + const std::vector &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> &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 + void + assemble_dg_matrix(MatrixType &system_matrix, + const FiniteElement &fe_dg, + const AgglomerationHandler &ah) + { + static_assert( + (std::is_same_v || + std::is_same_v>)); + + Assert((dynamic_cast *>(&fe_dg) || + dynamic_cast *>(&fe_dg) || + dynamic_cast *>(&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 &>(ah) + .create_agglomeration_sparsity_pattern(dsp); + system_matrix.reinit(dsp); + const unsigned int dofs_per_cell = fe_dg.n_dofs_per_cell(); + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix M11(dofs_per_cell, dofs_per_cell); + FullMatrix M12(dofs_per_cell, dofs_per_cell); + FullMatrix M21(dofs_per_cell, dofs_per_cell); + FullMatrix M22(dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + std::vector 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); + } + } + } + // get volumetric DoFs + polytope->get_dof_indices(local_dof_indices); + // Assemble face terms + 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)) + { + // Get normal vectors seen from each agglomeration. + 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()) + { + // use both fevalues + 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); + // distribute DoFs accordingly + // fluxes + 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 + { + // neigh polytope is ghosted, so retrieve necessary + // metadata. + 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.; + // 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 + assemble_local_jumps_and_averages_ghost( + M11, + M12, + M21, + M22, + fe_faces0, + recv_values, + recv_gradients, + recv_jxws, + penalty_constant, + h_f); + // distribute DoFs accordingly + // fluxes + 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 + void + assemble_dg_matrix_on_standard_mesh(MatrixType &system_matrix, + VectorType &system_rhs, + const Mapping &mapping, + const FiniteElement &fe_dg, + const DoFHandler &dof_handler) + { + static_assert( + (std::is_same_v || + std::is_same_v>)); + + Assert((dynamic_cast *>(&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 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 fe_faces0(mapping, + fe_dg, + QGaussSimplex(quadrature_degree), + update_values | update_JxW_values | + update_gradients | update_quadrature_points | + update_normal_vectors); + + FEValues fe_values(mapping, + fe_dg, + QGaussSimplex(quadrature_degree), + update_values | update_JxW_values | + update_gradients | update_quadrature_points); + + FEFaceValues fe_faces1(mapping, + fe_dg, + QGaussSimplex(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 cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + FullMatrix M11(dofs_per_cell, dofs_per_cell); + FullMatrix M12(dofs_per_cell, dofs_per_cell); + FullMatrix M21(dofs_per_cell, dofs_per_cell); + FullMatrix M22(dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + // Loop over standard deal.II cells + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + cell_matrix = 0.; + cell_rhs = 0.; + + fe_values.reinit(cell); + + // const auto &q_points = fe_values.get_quadrature_points(); + // const unsigned int n_qpoints = q_points.size(); + + 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 + } + } + + // distribute volumetric DoFs + 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 + local_dof_indices_neighbor(dofs_per_cell); + + M11 = 0.; + M12 = 0.; + M21 = 0.; + M22 = 0.; + + const auto &normals = fe_faces0.get_normal_vectors(); + // M11 + 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); + + // A10 + 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); + + // A11 + 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); + } + } + } + + // distribute DoFs accordingly + + 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 dealii::PolyUtils + +#endif diff --git a/agglomeration_poisson/mesh/unit_square_quad_unstructured.geo b/agglomeration_poisson/mesh/unit_square_quad_unstructured.geo new file mode 100644 index 00000000..54c2eab5 --- /dev/null +++ b/agglomeration_poisson/mesh/unit_square_quad_unstructured.geo @@ -0,0 +1,27 @@ +// ----------------------------------------------------------------------------- +// +// Gmsh GEO tutorial 11 +// +// Unstructured quadrangular meshes +// +// ----------------------------------------------------------------------------- + +// We have seen in tutorials `t3.geo' and `t6.geo' that extruded and transfinite +// meshes can be "recombined" into quads, prisms or hexahedra by using the +// "Recombine" keyword. Unstructured meshes can be recombined in the same +// way. Let's define a simple geometry with an analytical mesh size field: + +Point(1) = {0., 0., 0}; Point(2) = {1., 0., 0}; +Point(3) = {1., 1., 0}; Point(4) = {0., 1., 0}; + +Line(1) = {1, 2}; Line(2) = {2, 3}; +Line(3) = {3, 4}; Line(4) = {4, 1}; + +Curve Loop(4) = {1, 2, 3, 4}; Plane Surface(100) = {4}; + +// Field[1] = MathEval; +// Field[1].F = "0.01*(1.0+30.*(y-x*x)*(y-x*x) + (1-x)*(1-x))"; +Background Field = 1; + +// To generate quadrangles instead of triangles, we can simply add +Recombine Surface{100}; diff --git a/agglomeration_poisson/mesh/unit_square_quad_unstructured.msh b/agglomeration_poisson/mesh/unit_square_quad_unstructured.msh new file mode 100644 index 00000000..08ec6b77 --- /dev/null +++ b/agglomeration_poisson/mesh/unit_square_quad_unstructured.msh @@ -0,0 +1,382 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +4 4 1 0 +1 0 0 0 0 +2 1 0 0 0 +3 1 1 0 0 +4 0 1 0 0 +1 0 0 0 1 0 0 0 2 1 -2 +2 1 0 0 1 1 0 0 2 2 -3 +3 0 1 0 1 1 0 0 2 3 -4 +4 0 0 0 0 1 0 0 2 4 -1 +100 0 0 0 1 1 0 0 4 1 2 3 4 +$EndEntities +$Nodes +9 108 1 108 +0 1 0 1 +1 +0 0 0 +0 2 0 1 +2 +1 0 0 +0 3 0 1 +3 +1 1 0 +0 4 0 1 +4 +0 1 0 +1 1 0 7 +5 +6 +7 +8 +9 +10 +11 +0.1249999999997731 0 0 +0.2499999999994107 0 0 +0.374999999999046 0 0 +0.4999999999986922 0 0 +0.6249999999990107 0 0 +0.7499999999993404 0 0 +0.8749999999996702 0 0 +1 2 0 7 +12 +13 +14 +15 +16 +17 +18 +1 0.1249999999997731 0 +1 0.2499999999994107 0 +1 0.374999999999046 0 +1 0.4999999999986922 0 +1 0.6249999999990107 0 +1 0.7499999999993404 0 +1 0.8749999999996702 0 +1 3 0 7 +19 +20 +21 +22 +23 +24 +25 +0.8749999999995015 1 0 +0.7500000000003476 1 0 +0.6250000000012155 1 0 +0.5000000000020615 1 0 +0.3750000000015629 1 0 +0.2500000000010419 1 0 +0.1250000000005209 1 0 +1 4 0 7 +26 +27 +28 +29 +30 +31 +32 +0 0.8749999999995015 0 +0 0.7500000000003476 0 +0 0.6250000000012155 0 +0 0.5000000000020615 0 +0 0.3750000000015629 0 +0 0.2500000000010419 0 +0 0.1250000000005209 0 +2 100 0 76 +33 +34 +35 +36 +37 +38 +39 +40 +41 +42 +43 +44 +45 +46 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +58 +59 +60 +61 +62 +63 +64 +65 +66 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +81 +82 +83 +84 +85 +86 +87 +88 +89 +90 +91 +92 +93 +94 +95 +96 +97 +98 +99 +100 +101 +102 +103 +104 +105 +106 +107 +108 +0.503501210952674 0.5305839901209262 0 +0.6942407700434752 0.6798621195953681 0 +0.3206341550738262 0.6968598453070916 0 +0.6769393310403389 0.3069458642867186 0 +0.3060189717879025 0.3229547154352804 0 +0.4867599217093032 0.7950271064161099 0 +0.2040175824919035 0.4835528617084321 0 +0.7926509560572483 0.5126839902182272 0 +0.5128646815587427 0.2073338332504069 0 +0.8392508566777073 0.8139831352234977 0 +0.1860844363657626 0.8392712019596287 0 +0.8140435920095974 0.1607975324360289 0 +0.1607966467592726 0.1862067607410803 0 +0.5508382332704997 0.6244740025302096 0 +0.3601191571983056 0.5684523279534492 0 +0.5969176378959661 0.4622951541503035 0 +0.4455186052276263 0.3799205198266247 0 +0.6222813900350821 0.8487323410419246 0 +0.1513741490112091 0.6225906059959003 0 +0.84784383128208 0.378128153976362 0 +0.3779481325802844 0.1523890841137806 0 +0.3719857633831656 0.8893223089876053 0 +0.889128196678698 0.6280386633289134 0 +0.1107096319053588 0.3719118844951486 0 +0.6280967071932586 0.1110055038827389 0 +0.4308509710180405 0.6710747048951622 0 +0.334061493425355 0.4344021003471804 0 +0.6594246712171535 0.5715285257732502 0 +0.5602136872407 0.3385736243784459 0 +0.4462815025779339 0.8913638390052673 0 +0.1085594630173113 0.4458179219241999 0 +0.8908681638441919 0.5536706834588853 0 +0.5537152226704879 0.1092900472702152 0 +0.731826769150998 0.8988216172351036 0 +0.1011935139974098 0.7318503656403255 0 +0.8988167872544847 0.2682261532693396 0 +0.2682236232881479 0.1012750181307538 0 +0.8730092650990028 0.7240647505683416 0 +0.7241400641280884 0.1270837420319705 0 +0.2760172229683994 0.8730643835132147 0 +0.1270275798675799 0.2761499260468754 0 +0.5919509445920679 0.7345204078151798 0 +0.2640773889525423 0.5956657197441078 0 +0.7277737674301493 0.4104521545407531 0 +0.4069046912558426 0.2688610821336825 0 +0.5438374419139871 0.8791137588528191 0 +0.1206797703020509 0.5436966939992822 0 +0.8778584871418256 0.4563929386001843 0 +0.4560772247802358 0.122015046692696 0 +0.367529277632855 0.7784363054940432 0 +0.7886783313038435 0.6398938192552009 0 +0.2105357569132224 0.3605795881872284 0 +0.6389334988909361 0.2114819951281795 0 +0.9028848636438731 0.902884865794588 0 +0.09711513422304896 0.9028848636462924 0 +0.9028848658136241 0.0971151363591552 0 +0.09711513635072225 0.09711513424948652 0 +0.7370185046385402 0.7941855903152235 0 +0.2062631004804109 0.7370141125252966 0 +0.7934931167346435 0.2633360516249019 0 +0.2630124346220526 0.2066494067601173 0 +0.801524995972853 0.7148893548287539 0 +0.2865416927366418 0.8000908126573411 0 +0.7144340528033952 0.1987431376655361 0 +0.1984233033966108 0.2857651340123553 0 +0.4847726232734851 0.5689026898049103 0 +0.5409513541683244 0.5361388899809862 0 +0.4833125112932423 0.4597674230337673 0 +0.7995128737200753 0.904499878503802 0 +0.09554631046526975 0.7995200971207616 0 +0.904470753876745 0.2005266226161684 0 +0.2005056778287199 0.09561066399559177 0 +0.9064543058193955 0.8264398800501755 0 +0.82645870073321 0.09356716827651582 0 +0.1735867980507265 0.906463920764619 0 +0.09355871740558798 0.1736313945069894 0 +$EndNodes +$Elements +9 127 1 127 +0 1 15 1 +1 1 +0 2 15 1 +2 2 +0 3 15 1 +3 3 +0 4 15 1 +4 4 +1 1 1 8 +5 1 5 +6 5 6 +7 6 7 +8 7 8 +9 8 9 +10 9 10 +11 10 11 +12 11 2 +1 2 1 8 +13 2 12 +14 12 13 +15 13 14 +16 14 15 +17 15 16 +18 16 17 +19 17 18 +20 18 3 +1 3 1 8 +21 3 19 +22 19 20 +23 20 21 +24 21 22 +25 22 23 +26 23 24 +27 24 25 +28 25 4 +1 4 1 8 +29 4 26 +30 26 27 +31 27 28 +32 28 29 +33 29 30 +34 30 31 +35 31 32 +36 32 1 +2 100 3 91 +37 56 30 31 73 +38 54 23 24 72 +39 57 9 10 71 +40 55 16 17 70 +41 58 46 74 38 +42 59 47 75 39 +43 60 48 76 40 +44 61 49 77 41 +45 60 34 74 46 +46 58 35 75 47 +47 61 36 76 48 +48 59 37 77 49 +49 62 54 82 38 +50 64 55 83 40 +51 63 56 84 39 +52 65 57 85 41 +53 78 38 74 50 +54 79 39 75 51 +55 80 40 76 52 +56 81 41 77 53 +57 50 74 34 90 +58 51 75 35 91 +59 52 76 36 92 +60 53 77 37 93 +61 54 62 22 23 +62 56 63 29 30 +63 55 64 15 16 +64 57 65 8 9 +65 12 88 11 2 +66 5 89 32 1 +67 19 86 18 3 +68 26 87 25 4 +69 35 58 38 82 +70 34 60 40 83 +71 37 59 39 84 +72 36 61 41 85 +73 50 21 22 78 +74 51 28 29 79 +75 52 14 15 80 +76 53 7 8 81 +77 47 100 33 98 +78 48 99 33 100 +79 47 98 46 58 +80 48 60 46 99 +81 49 100 47 59 +82 49 61 48 100 +83 38 78 22 62 +84 39 79 29 63 +85 40 80 15 64 +86 41 81 8 65 +87 21 50 90 66 +88 28 51 91 67 +89 14 52 92 68 +90 7 53 93 69 +91 70 42 90 94 +92 71 44 92 96 +93 72 43 91 95 +94 73 45 93 97 +95 101 66 90 42 +96 102 67 91 43 +97 103 68 92 44 +98 104 69 93 45 +99 21 66 101 20 +100 28 67 102 27 +101 14 68 103 13 +102 7 69 104 6 +103 42 70 17 105 +104 43 72 24 107 +105 44 71 10 106 +106 45 73 31 108 +107 83 55 70 94 +108 85 57 71 96 +109 82 54 72 95 +110 84 56 73 97 +111 101 42 105 86 +112 102 43 107 87 +113 103 44 106 88 +114 104 45 108 89 +115 83 94 90 34 +116 82 95 91 35 +117 85 96 92 36 +118 84 97 93 37 +119 20 101 86 19 +120 27 102 87 26 +121 13 103 88 12 +122 6 104 89 5 +123 17 18 86 105 +124 24 25 87 107 +125 10 11 88 106 +126 31 32 89 108 +127 33 99 46 98 +$EndElements diff --git a/agglomeration_poisson/source/agglomeration_handler.cc b/agglomeration_poisson/source/agglomeration_handler.cc new file mode 100644 index 00000000..96e3ef92 --- /dev/null +++ b/agglomeration_poisson/source/agglomeration_handler.cc @@ -0,0 +1,1649 @@ +/* ----------------------------------------------------------------------------- + * + * 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 + +#include + +template +AgglomerationHandler::AgglomerationHandler( + const GridTools::Cache &cache_tria) + : cached_tria(std::make_unique>( + cache_tria.get_triangulation(), + cache_tria.get_mapping())) + , communicator(cache_tria.get_triangulation().get_communicator()) +{ + Assert(dim == spacedim, ExcNotImplemented("Not available with codim > 0")); + Assert(dim == 2 || dim == 3, ExcImpossibleInDim(1)); + Assert((dynamic_cast *>( + &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 +typename AgglomerationHandler::agglomeration_iterator +AgglomerationHandler::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 + + // First index drives the selection of the master cell. After that, store the + // master cell. + 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::active_cell_iterator cell_dh = + cells[0]->as_dof_handler_iterator(agglo_dh); + cell_dh->set_active_fe_index(CellAgglomerationType::master); + + // Store slave cells and save the relationship with the parent + std::vector::active_cell_iterator> + slaves; + slaves.reserve(cells.size() - 1); + // exclude first cell since it's the master cell + 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::active_cell_iterator cell = + (*it)->as_dof_handler_iterator(agglo_dh); + cell->set_active_fe_index(CellAgglomerationType::slave); // slave cell + + // 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. + 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 + + // Store the slaves of each master + master2slaves[master_idx] = slaves; + // Save to which polygon this agglomerate correspond + master2polygon[master_idx] = n_agglomerations; + + ++n_agglomerations; // an agglomeration has been performed, record it + + create_bounding_box(cells); // fill the vector of bboxes + + // Finally, return a polygonal iterator to the polytope just constructed. + return {cells[0], this}; +} + +template +typename AgglomerationHandler::agglomeration_iterator +AgglomerationHandler::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 + + // First index drives the selection of the master cell. After that, store the + // master cell. + 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::active_cell_iterator cell_dh = + cells[0]->as_dof_handler_iterator(agglo_dh); + cell_dh->set_active_fe_index(CellAgglomerationType::master); + + // Store slave cells and save the relationship with the parent + std::vector::active_cell_iterator> + slaves; + slaves.reserve(cells.size() - 1); + // exclude first cell since it's the master cell + 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::active_cell_iterator cell = + (*it)->as_dof_handler_iterator(agglo_dh); + cell->set_active_fe_index( + fecollection_size); // slave cell (the last index) + + // 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. + 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 + + // Store the slaves of each master + master2slaves[master_idx] = slaves; + // Save to which polygon this agglomerate correspond + master2polygon[master_idx] = n_agglomerations; + + ++n_agglomerations; // an agglomeration has been performed, record it + + create_bounding_box(cells); // fill the vector of bboxes + + // Finally, return a polygonal iterator to the polytope just constructed. + return {cells[0], this}; +} + + +template +void +AgglomerationHandler::initialize_fe_values( + const Quadrature &cell_quadrature, + const UpdateFlags &flags, + const Quadrature &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>(*mapping, + dummy_fe, + agglomeration_quad, + update_quadrature_points | + update_JxW_values); // only for quadrature + no_face_values = std::make_unique>( + *mapping, + dummy_fe, + agglomeration_face_quad, + update_quadrature_points | update_JxW_values | + update_normal_vectors); // only for quadrature +} + +template +void +AgglomerationHandler::initialize_fe_values( + const hp::QCollection &cell_qcollection, + const UpdateFlags &flags, + const hp::QCollection &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(*mapping); + dummy_fe_collection = hp::FECollection(dummy_fe); + hp_no_values = std::make_unique>( + mapping_collection, + dummy_fe_collection, + agglomeration_quad_collection, + update_quadrature_points | update_JxW_values); // only for quadrature + + hp_no_face_values = std::make_unique>( + mapping_collection, + dummy_fe_collection, + agglomeration_face_quad_collection, + update_quadrature_points | update_JxW_values | + update_normal_vectors); // only for quadrature +} + + + +template +unsigned int +AgglomerationHandler::n_agglomerated_faces_per_cell( + const typename Triangulation::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 +void +AgglomerationHandler::initialize_agglomeration_data( + const std::unique_ptr> &cache_tria) +{ + tria = &(cache_tria->get_triangulation()); + mapping = &(cache_tria->get_mapping()); + + agglo_dh.reinit(*tria); + + if (const auto parallel_tria = dynamic_cast< + const dealii::parallel::TriangulationBase *>(&*tria)) + { + const std::weak_ptr 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(); + + // First, update the pointer + cached_tria = std::make_unique>( + cache_tria->get_triangulation(), cache_tria->get_mapping()); + + connect_to_tria_signals(); + n_agglomerations = 0; +} + + + +template +void +AgglomerationHandler::distribute_agglomerated_dofs( + const FiniteElement &fe_space) +{ + if (dynamic_cast *>(&fe_space)) + fe = std::make_unique>(fe_space.degree); + else if (dynamic_cast *>(&fe_space)) + fe = std::make_unique>(fe_space.degree); + else + AssertThrow( + false, + ExcNotImplemented( + "Currently, this interface supports only DGQ and DGP bases.")); + + box_mapping = std::make_unique>( + bboxes, + master2polygon); // construct bounding box mapping + + if (hybrid_mesh) + { + // 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. + standard_scratch = + std::make_unique(*mapping, + *fe, + QGauss(2 * fe_space.degree + 2), + internal_agglomeration_flags); + } + + + fe_collection.push_back(*fe); // master + fe_collection.push_back( + FE_Nothing(fe->reference_cell())); // slave + + initialize_hp_structure(); + + // in case the tria is distributed, communicate ghost information with + // neighboring ranks + const bool needs_ghost_info = + dynamic_cast *>(&*tria) != + nullptr; + if (needs_ghost_info) + setup_ghost_polytopes(); + + setup_connectivity_of_agglomeration(); + + if (needs_ghost_info) + exchange_interface_values(); +} + +template +void +AgglomerationHandler::distribute_agglomerated_dofs( + const hp::FECollection &fe_collection_in) +{ + is_hp_collection = true; + + hp_fe_collection = std::make_unique>( + fe_collection_in); // copy the input collection + + box_mapping = std::make_unique>( + 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 *>(&fe_collection_in[i])) + { + // System case + for (unsigned int b = 0; b < fe_collection_in[i].n_base_elements(); + ++b) + { + if (!(dynamic_cast *>( + &fe_collection_in[i].base_element(b)) || + dynamic_cast *>( + &fe_collection_in[i].base_element(b)) || + dynamic_cast *>( + &fe_collection_in[i].base_element(b)))) + AssertThrow( + false, + ExcNotImplemented( + "Currently, this interface supports only DGQ and DGP bases.")); + } + } + else + { + // Scalar case + if (!(dynamic_cast *>(&fe_collection_in[i]) || + dynamic_cast *>(&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()); + } + else if (fe_collection[0].n_components() > 1) + { + std::vector *> base_elements; + std::vector multiplicities; + for (unsigned int b = 0; b < fe_collection[0].n_base_elements(); ++b) + { + base_elements.push_back(new FE_Nothing()); + multiplicities.push_back(fe_collection[0].element_multiplicity(b)); + } + FESystem fe_system_nothing(base_elements, multiplicities); + for (const auto *ptr : base_elements) + delete ptr; + fe_collection.push_back(fe_system_nothing); + } + + initialize_hp_structure(); + + // in case the tria is distributed, communicate ghost information with + // neighboring ranks + const bool needs_ghost_info = + dynamic_cast *>(&*tria) != + nullptr; + if (needs_ghost_info) + setup_ghost_polytopes(); + + setup_connectivity_of_agglomeration(); + + if (needs_ghost_info) + exchange_interface_values(); +} + +template +void +AgglomerationHandler::create_bounding_box( + const AgglomerationContainer &polytope) +{ + Assert(n_agglomerations > 0, + ExcMessage("No agglomeration has been performed.")); + Assert(dim > 1, ExcNotImplemented()); + + std::vector> 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 +void +AgglomerationHandler::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:: + setup_master_neighbor_connectivity(cell, *this); + } + + if (Utilities::MPI::job_supports_mpi()) + { + // communicate the number of faces + recv_n_faces = Utilities::MPI::some_to_some(communicator, local_n_faces); + + // send information about boundaries and neighboring polytopes id + 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 +void +AgglomerationHandler::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()) + { + // Neighboring polytope is ghosted. + + // Compute shape functions at the interface + const auto ¤t_fe = reinit(polytope, f); + + std::vector> qpoints_to_send = + current_fe.get_quadrature_points(); + + const std::vector &jxws_to_send = + current_fe.get_JxW_values(); + + const std::vector> &normals_to_send = + current_fe.get_normal_vectors(); + + + const types::subdomain_id neigh_rank = + neigh_polytope->subdomain_id(); + + std::pair cell_and_face{ + polytope->id(), f}; + // Prepare data to send + 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(); + + // TODO: check `agglomeration_flags` before computing + // values and gradients. + std::vector> values_per_qpoints( + dofs_per_cell); + + std::vector>> + 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); + } + } + } + } + } + + // Finally, exchange with neighboring ranks + 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 +Quadrature +AgglomerationHandler::agglomerated_quadrature( + const typename AgglomerationHandler::AgglomerationContainer + &cells, + const typename Triangulation::active_cell_iterator + &master_cell) const +{ + Assert(is_master_cell(master_cell), + ExcMessage("This must be a master cell.")); + + std::vector> vec_pts; + std::vector vec_JxWs; + + if (!is_hp_collection) + { + // Original version: handle case without hp::FECollection + 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 &p) { return p; }); + std::transform(JxWs.begin(), + JxWs.end(), + std::back_inserter(vec_JxWs), + [&](const double w) { return w; }); + } + } + else + { + // Handle the hp::FECollection case + const auto &master_cell_as_dh_iterator = + master_cell->as_dof_handler_iterator(agglo_dh); + for (const auto &dummy_cell : cells) + { + // 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. + 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 &p) { return p; }); + std::transform(JxWs.begin(), + JxWs.end(), + std::back_inserter(vec_JxWs), + [&](const double w) { return w; }); + } + } + + // Map back each point in real space by using the map associated to the + // bounding box. + std::vector> 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(unit_points, vec_JxWs); +} + + + +template +void +AgglomerationHandler::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); + // euler_mapping = std::make_unique< + // MappingFEField>>( euler_dh, euler_vector); +} + + + +template +const FEValues & +AgglomerationHandler::reinit( + const AgglomerationIterator &polytope) const +{ + // Assert(euler_mapping, + // ExcMessage("The mapping describing the physical element stemming + // from " + // "agglomeration has not been set up.")); + + const auto &deal_cell = polytope->as_dof_handler_iterator(agglo_dh); + + // 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); + + const auto &agglo_cells = polytope->get_agglomerate(); + + Quadrature agglo_quad = agglomerated_quadrature(agglo_cells, deal_cell); + + if (!is_hp_collection) + { + // Original version: handle case without hp::FECollection + agglomerated_scratch = std::make_unique(*box_mapping, + fe_collection[0], + agglo_quad, + agglomeration_flags); + } + else + { + // Handle the hp::FECollection case + agglomerated_scratch = std::make_unique(*box_mapping, + polytope->get_fe(), + agglo_quad, + agglomeration_flags); + } + return agglomerated_scratch->reinit(deal_cell); +} + + + +template +const FEValuesBase & +AgglomerationHandler::reinit_master( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_index, + std::unique_ptr> + &agglo_isv_ptr) const +{ + return internal::AgglomerationHandlerImplementation:: + reinit_master(cell, face_index, agglo_isv_ptr, *this); +} + + + +template +const FEValuesBase & +AgglomerationHandler::reinit( + const AgglomerationIterator &polytope, + const unsigned int face_index) const +{ + // Assert(euler_mapping, + // ExcMessage("The mapping describing the physical element stemming + // from " + // "agglomeration has not been set up.")); + + 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:: + reinit_master(deal_cell, face_index, agglomerated_isv_bdary, *this); +} + + + +template +std::pair &, + const FEValuesBase &> +AgglomerationHandler::reinit_interface( + const AgglomerationIterator &polytope_in, + const AgglomerationIterator &neigh_polytope, + const unsigned int local_in, + const unsigned int local_neigh) const +{ + // If current and neighboring polytopes are both locally owned, then compute + // the jump in the classical way without needing information about ghosted + // entities. + 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:: + reinit_master(cell_in, local_in, agglomerated_isv, *this); + const auto &fe_out = + internal::AgglomerationHandlerImplementation:: + reinit_master(neigh_cell, local_neigh, agglomerated_isv_neigh, *this); + std::pair &, + const FEValuesBase &> + 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())]; + // const double bbox_measure = bbox.volume(); + + const unsigned int neigh_rank = neigh_polytope->subdomain_id(); + const CellId &neigh_id = neigh_polytope->id(); + + // Retrieve qpoints,JxWs, normals sent previously from the neighboring + // rank. + std::vector> &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> &normals = + recv_normals.at(neigh_rank).at({neigh_id, local_neigh}); + + // Apply the necessary scalings due to the bbox. + std::vector> final_unit_q_points; + std::transform(real_qpoints.begin(), + real_qpoints.end(), + std::back_inserter(final_unit_q_points), + [&](const Point &p) { + return bbox.real_to_unit(p); + }); + + // std::vector scale_factors(final_unit_q_points.size()); + // std::vector scaled_weights(final_unit_q_points.size()); + // std::vector> 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(); + // } + for (unsigned int q = 0; q < final_unit_q_points.size(); ++q) + normals[q] *= -1; + + + NonMatching::ImmersedSurfaceQuadrature surface_quad( + final_unit_q_points, JxWs, normals); + + agglomerated_isv = + std::make_unique>( + *box_mapping, *fe, surface_quad, agglomeration_face_flags); + + + agglomerated_isv->reinit(cell); + + std::pair &, + const FEValuesBase &> + my_p(*agglomerated_isv, *agglomerated_isv); + + return my_p; + } +} + + + +template +template +void +AgglomerationHandler::create_agglomeration_sparsity_pattern( + SparsityPatternType &dsp, + const AffineConstraints &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) + dsp.reinit(locally_owned_dofs.size(), + locally_owned_dofs.size(), + locally_relevant_dofs); + else if constexpr (std::is_same_v) + dsp.reinit(locally_owned_dofs, communicator); + else + AssertThrow(false, ExcNotImplemented()); + + // Create the sparsity pattern corresponding only to volumetric terms. The + // fluxes needed by DG methods will be filled later. + DoFTools::make_sparsity_pattern( + agglo_dh, dsp, constraints, keep_constrained_dofs, subdomain_id); + + + if (!is_hp_collection) + { + // Original version: handle case without hp::FECollection + const unsigned int dofs_per_cell = agglo_dh.get_fe(0).n_dofs_per_cell(); + std::vector current_dof_indices(dofs_per_cell); + std::vector neighbor_dof_indices(dofs_per_cell); + + // Loop over all locally owned polytopes, find the neighbor (also ghosted) + // and add fluxes to the sparsity pattern. + 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 + { + // Handle the hp::FECollection case + + // Loop over all locally owned polytopes, find the neighbor (also ghosted) + // and add fluxes to the sparsity pattern. + 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 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 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) + dsp.compress(); +} + + + +template +void +AgglomerationHandler::setup_ghost_polytopes() +{ + [[maybe_unused]] const auto parallel_triangulation = + dynamic_cast *>(&*tria); + Assert(parallel_triangulation != nullptr, ExcInternalError()); + + const unsigned int n_dofs_per_cell = fe->dofs_per_cell; + std::vector 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) + { + // interior, locally owned, cell + for (const auto &f : cell->face_indices()) + { + if (!cell->at_boundary(f)) + { + const auto &neighbor = cell->neighbor(f); + if (neighbor->is_ghost()) + { + // key of the map: the rank to which send the data + const types::subdomain_id neigh_rank = + neighbor->subdomain_id(); + + // inform the "standard" neighbor about the neighboring + // id and its master cell + local_cell_ids_neigh_cell[neigh_rank].emplace( + cell->id(), master_cell_id); + + // inform the neighboring rank that this master cell + // (hence polytope) has the following DoF indices + local_ghost_dofs[neigh_rank].emplace( + master_cell_id, global_dof_indices); + + // ...same for bounding boxes + 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); + + // Exchange with neighboring ranks the neighboring bounding boxes + recv_ghosted_bbox = + Utilities::MPI::some_to_some(communicator, local_ghosted_bbox); + + // Exchange with neighboring ranks the neighboring ghosted DoFs + recv_ghost_dofs = + Utilities::MPI::some_to_some(communicator, local_ghost_dofs); +} + + + +namespace dealii +{ + namespace internal + { + template + class AgglomerationHandlerImplementation + { + public: + static const FEValuesBase & + reinit_master( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face_index, + std::unique_ptr> + &agglo_isv_ptr, + const AgglomerationHandler &handler) + { + Assert(handler.is_master_cell(cell), + ExcMessage("This cell must be a master one.")); + + AgglomerationIterator it{cell, &handler}; + const auto &neigh_polytope = it->neighbor(face_index); + + const CellId polytope_in_id = cell->id(); + + // Retrieve the bounding box of the agglomeration + 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> final_unit_q_points; + std::vector final_weights; + std::vector> final_normals; + + if (!handler.is_hp_collection) + { + // Original version: handle case without hp::FECollection + 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 + { + // Handle the hp::FECollection case + 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 surface_quad( + final_unit_q_points, final_weights, final_normals); + + if (!handler.is_hp_collection) + { + agglo_isv_ptr = + std::make_unique>( + *(handler.box_mapping), + *(handler.fe), + surface_quad, + handler.agglomeration_face_flags); + } + else + { + agglo_isv_ptr = + std::make_unique>( + *(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::active_cell_iterator + &master_cell, + const AgglomerationHandler &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 visited_polygonal_neighbors; + + std::map face_to_neigh_id; + + std::map is_face_at_boundary; + + // same as above, but with CellId + std::set 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)) + { + // - 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. + handler.polygon_boundary[master_cell].push_back( + cell->face(f)); + + // global index of neighboring deal.II cell + const types::global_cell_index neighboring_cell_index = + neighboring_cell->active_cell_index(); + + // master cell for the neighboring polytope + 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)) + { + // index of the neighboring polytope + 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)) + { + // found a neighbor + + 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 + { + // neighboring cell is a master + + // save the pair of neighboring cells + 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)) + { + // found a neighbor + 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); + + // 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. + 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); + + // const auto master_index = + // master_indices[ghost_counter]; + + 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}; + + + // record the cell id of the neighboring polytope + 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; + + + // increment number of faces + ++handler.number_of_agglomerated_faces + [current_polytope_index]; + + visited_polygonal_neighbors_id.insert( + check_neigh_polytope_id); + + // ghosted polytope has been found, increment + // ghost counter + ++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); + + // std::cout << "ADDED (" + // << cell->active_cell_index() << ") + // BETWEEN " + // << current_polytope_id << " e " + // << check_neigh_polytope_id << + // std::endl; + + 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()) + { + // Boundary face of a boundary cell. + // Note that the neighboring cell must be invalid. + + handler.polygon_boundary[master_cell].push_back( + cell->face(f)); + + if (visited_polygonal_neighbors.find( + std::numeric_limits::max()) == + std::end(visited_polygonal_neighbors)) + { + // boundary face. Notice that `neighboring_cell` is + // invalid here. + 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::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 dealii::parallel::TriangulationBase *>( + &(*handler.tria)); + + const unsigned int n_faces_current_poly = + handler.number_of_agglomerated_faces[current_polytope_index]; + + // 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 + 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 &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 &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 &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 &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 &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 &constraints, + const bool keep_constrained_dofs, + const types::subdomain_id subdomain_id); diff --git a/agglomeration_poisson/source/mapping_box.cc b/agglomeration_poisson/source/mapping_box.cc new file mode 100644 index 00000000..7e65a2eb --- /dev/null +++ b/agglomeration_poisson/source/mapping_box.cc @@ -0,0 +1,983 @@ +/* ----------------------------------------------------------------------------- + * + * 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 +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include + +#include + +#include + +#include + +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 +bool +has_box(const CellType &cell, + const std::map + &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 +MappingBox::MappingBox( + const std::vector> &input_boxes, + const std::map + &global_to_polytope) +{ + Assert(input_boxes.size() > 0, + ExcMessage("Invalid number of bounding boxes.")); + + // copy boxes and map + 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 +MappingBox::InternalData::InternalData(const Quadrature &q) + : cell_extents(numbers::signaling_nan>()) + , traslation(numbers::signaling_nan>()) + , inverse_cell_extents(numbers::signaling_nan>()) + , volume_element(numbers::signaling_nan()) + , quadrature_points(q.get_points()) +{} + + + +template +void +MappingBox::InternalData::reinit(const UpdateFlags update_flags, + const Quadrature &) +{ + // 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 + this->update_each = update_flags; +} + + + +template +std::size_t +MappingBox::InternalData::memory_consumption() const +{ + return (Mapping::InternalDataBase::memory_consumption() + + MemoryConsumption::memory_consumption(cell_extents) + + MemoryConsumption::memory_consumption(traslation) + + MemoryConsumption::memory_consumption(inverse_cell_extents) + + MemoryConsumption::memory_consumption(volume_element)); +} + + + +template +bool +MappingBox::preserves_vertex_locations() const +{ + return true; +} + + + +template +bool +MappingBox::is_compatible_with( + const ReferenceCell &reference_cell) 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 +UpdateFlags +MappingBox::requires_update_flags(const UpdateFlags in) const +{ + // 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 + UpdateFlags out = in; + if (out & update_boundary_forms) + out |= update_normal_vectors; + + return out; +} + + + +template +std::unique_ptr::InternalDataBase> +MappingBox::get_data(const UpdateFlags update_flags, + const Quadrature &q) const +{ + std::unique_ptr::InternalDataBase> data_ptr = + std::make_unique(); + data_ptr->reinit(requires_update_flags(update_flags), q); + + return data_ptr; +} + + + +template +std::unique_ptr::InternalDataBase> +MappingBox::get_subface_data( + const UpdateFlags update_flags, + const Quadrature &quadrature) const +{ + (void)update_flags; + (void)quadrature; + DEAL_II_NOT_IMPLEMENTED(); + return {}; +} + + + +template +void +MappingBox::update_cell_extents( + const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const InternalData &data) const +{ + // Compute start point and sizes along axes. The vertices to be looked at + // are 1, 2, 4 compared to the base vertex 0. + if (cell_similarity != CellSimilarity::translation) + { + const BoundingBox ¤t_box = + boxes[polytope_translator.at(cell->active_cell_index())]; + const std::pair, Point> &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 + void + transform_quadrature_points( + const BoundingBox &box, + const ArrayView> &unit_quadrature_points, + std::vector> &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 +void +MappingBox::maybe_update_cell_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const InternalData &data, + const ArrayView> &unit_quadrature_points, + std::vector> &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 +void +MappingBox::maybe_update_normal_vectors( + const unsigned int face_no, + const InternalData &data, + std::vector> &normal_vectors) const +{ + // compute normal vectors. All normals on a face have the same value. + if (data.update_each & update_normal_vectors) + { + Assert(face_no < GeometryInfo::faces_per_cell, ExcInternalError()); + std::fill(normal_vectors.begin(), + normal_vectors.end(), + GeometryInfo::unit_normal_vector[face_no]); + } +} + + + +template +void +MappingBox::maybe_update_jacobian_derivatives( + const InternalData &data, + const CellSimilarity::Similarity cell_similarity, + internal::FEValuesImplementation::MappingRelatedData + &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 +void +MappingBox::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 +void +MappingBox::maybe_update_jacobians( + const InternalData &data, + const CellSimilarity::Similarity cell_similarity, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + // "compute" Jacobian at the quadrature points, which are all the + // same + 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 +void +MappingBox::maybe_update_inverse_jacobians( + const InternalData &data, + const CellSimilarity::Similarity cell_similarity, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + // "compute" inverse Jacobian at the quadrature points, which are + // all the same + 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 +CellSimilarity::Similarity +MappingBox::fill_fe_values( + const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); + + // convert data object to internal data for this class. fails with + // an exception if that is not possible + Assert(dynamic_cast(&internal_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(internal_data); + + + update_cell_extents(cell, cell_similarity, data); + + maybe_update_cell_quadrature_points(cell, + data, + quadrature.get_points(), + output_data.quadrature_points); + + // compute Jacobian determinant. all values are equal and are the + // product of the local lengths in each coordinate direction + 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 +void +MappingBox::fill_fe_subface_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + (void)cell; + (void)face_no; + (void)subface_no; + (void)quadrature; + (void)internal_data; + (void)output_data; + DEAL_II_NOT_IMPLEMENTED(); +} + + + +template +void +MappingBox::fill_fe_immersed_surface_values( + const typename Triangulation::cell_iterator &cell, + const NonMatching::ImmersedSurfaceQuadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + AssertDimension(dim, spacedim); + Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); + + // Convert data object to internal data for this class. Fails with an + // exception if that is not possible. + Assert(dynamic_cast(&internal_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(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 +void +MappingBox::transform( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) const +{ + AssertDimension(input.size(), output.size()); + Assert(dynamic_cast(&mapping_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_kind) + { + case mapping_covariant: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::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::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::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::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 +void +MappingBox::transform( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) const +{ + AssertDimension(input.size(), output.size()); + Assert(dynamic_cast(&mapping_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_kind) + { + case mapping_covariant: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::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::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::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::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::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::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::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::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 +void +MappingBox::transform( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) const +{ + AssertDimension(input.size(), output.size()); + Assert(dynamic_cast(&mapping_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_kind) + { + case mapping_covariant: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::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::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::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::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::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::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::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::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 +void +MappingBox::transform( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) const +{ + AssertDimension(input.size(), output.size()); + Assert(dynamic_cast(&mapping_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_kind) + { + case mapping_covariant_gradient: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::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 +void +MappingBox::transform( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) const +{ + AssertDimension(input.size(), output.size()); + Assert(dynamic_cast(&mapping_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_kind) + { + case mapping_contravariant_hessian: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::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::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::ExcAccessToUninitializedField( + "update_covariant_transformation")); + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::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 +Point +MappingBox::transform_unit_to_real_cell( + const typename Triangulation::cell_iterator &cell, + const Point &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 +Point +MappingBox::transform_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const Point &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 +void +MappingBox::transform_points_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const ArrayView> &real_points, + const ArrayView> &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 +std::unique_ptr> +MappingBox::clone() const +{ + return std::make_unique>(*this); +} + + +//--------------------------------------------------------------------------- +// explicit instantiations +template class MappingBox<1>; +template class MappingBox<2>; +template class MappingBox<3>; + + +DEAL_II_NAMESPACE_CLOSE