435 * The following headers provide the deal.II functionality needed in
this
436 * example. Most of them are standard components
for mesh handling, finite
437 * element mappings, linear algebra, and graphical output. In addition, we
438 * include the agglomeration-specific headers that define the
data structures
439 * and utilities used to construct and manage polytopal agglomerates.
443 * deal.II base utilities.
446 *
#include <deal.II/base/exceptions.h>
450 * Finite element mappings.
453 *
#include <deal.II/fe/mapping_fe.h>
457 * Grid generation, mesh input/output, and mesh-related utilities.
460 *
#include <deal.II/grid/grid_generator.h>
461 *
#include <deal.II/grid/grid_in.h>
462 *
#include <deal.II/grid/grid_out.h>
463 *
#include <deal.II/grid/grid_tools.h>
467 * Linear algebra objects and sparse direct solvers.
470 *
#include <deal.II/lac/precondition.h>
471 *
#include <deal.II/lac/solver_cg.h>
472 *
#include <deal.II/lac/sparse_direct.h>
473 *
#include <deal.II/lac/sparse_matrix.h>
477 * Output of finite element
data for visualization.
480 *
#include <deal.II/numerics/data_out.h>
484 * Agglomeration-specific headers used in
this example.
487 *
#include <agglomeration_handler.h>
488 *
#include <poly_utils.h>
492 *
C++ standard library headers.
495 *
#include <algorithm>
501 * We use the
struct ConvergenceInfo to store the number of degrees of freedom together
502 * with the corresponding
L2 and H1 errors, and print a simple
503 * convergence table to the console.
506 *
struct ConvergenceInfo
508 *
ConvergenceInfo() =
default;
513 *
vec_data.push_back(dofs_and_errs);
519 *
Assert(vec_data.size() > 0, ExcInternalError());
520 *
std::cout << std::left <<
"#DoFs, L2 error, H1 error" << std::endl;
522 *
for (
const auto &dof_and_errs : vec_data)
523 *
std::cout <<
std::scientific << dof_and_errs.
first <<
", "
528 *
std::vector<std::pair<types::global_dof_index, std::pair<double, double>>>
534 * We will compare the performance of three different partitioning strategies:
535 *
using METIS,
using an R-tree based agglomeration, or not performing any
536 * partitioning at all.
539 *
enum class PartitionerType
548 * We then implement the manufactured right-hand side
549 * f(x, y) = 2 π²
sin(π x)
sin(π y),
550 * which corresponds to the exact solution
551 * u(x, y) =
sin(π x)
sin(π y).
555 *
class RightHandSide :
public Function<dim>
564 *
std::vector<double> &values,
565 *
const unsigned int )
const override
567 *
for (
unsigned int i = 0; i <
values.size(); ++i)
577 * Exact solution is
set as u(x,y) =
sin(pi x)
sin(pi y).
578 * It is used to impose Dirichlet boundary conditions and to evaluate
579 * the
L2 and H1-seminorm errors. Its
gradient is also provided
for
580 * the computation of the H1 error.
584 *
class ExactSolution :
public Function<dim>
595 *
const unsigned int = 0) const override
602 *
std::vector<double> &values,
603 *
const unsigned int )
const override
605 *
for (
unsigned int i = 0; i <
values.size(); ++i)
606 *
values[i] = this->
value(points[i]);
611 *
const unsigned int = 0) const override
618 *
return return_value;
624 * The Poisson<dim>
class encapsulates the solution of the model Poisson
626 * @f[ -\Delta u = f \quad \text{in } \Omega, \qquad u = u_D \quad \text{on } \partial\Omega. @f]
627 * It sets up a fine triangulation, constructs agglomerated polytopal
628 * cells according to the chosen partitioning strategy, assembles the
629 *
symmetric interior penalty DG discretization on the agglomerated mesh,
630 * solves the resulting linear system, and
finally postprocesses the
631 * numerical solution by writing visualization output and computing
632 * global error norms.
642 *
setup_agglomeration();
653 *
std::unique_ptr<AgglomerationHandler<dim>> ah;
660 *
std::unique_ptr<GridTools::Cache<dim>> cached_tria;
661 *
std::unique_ptr<const Function<dim>> rhs_function;
662 *
std::unique_ptr<const Function<dim>> analytical_solution;
665 *
Poisson(
const PartitionerType &partitioner_type = PartitionerType::rtree,
666 *
const unsigned int = 0,
667 *
const unsigned int = 0,
668 *
const unsigned int fe_degree = 1);
673 *
get_n_dofs()
const;
675 *
std::pair<double, double>
678 *
PartitionerType partitioner_type;
679 *
unsigned int extraction_level;
680 *
unsigned int n_subdomains;
681 *
double penalty_constant = 60.;
689 * The constructor initializes the Poisson<dim> solver with the selected partitioning
690 * strategy, agglomeration parameters, polynomial degree, and the
691 * manufactured exact solution and right-hand side.
695 *
Poisson<dim>::Poisson(
const PartitionerType &partitioner_type,
696 *
const unsigned int extraction_level,
697 *
const unsigned int n_subdomains,
698 *
const unsigned int fe_degree)
701 *
, partitioner_type(partitioner_type)
702 *
, extraction_level(extraction_level)
703 *
, n_subdomains(n_subdomains)
704 *
, penalty_constant(10. * (fe_degree + 1) * (fe_degree + dim))
708 * Initialize manufactured solution.
711 *
analytical_solution = std::make_unique<ExactSolution<dim>>();
712 *
rhs_function = std::make_unique<const RightHandSide<dim>>();
713 *
constraints.close();
720 * Build the fine triangulation from a
Gmsh mesh,
apply a global
721 * refinement, initialize the cache and agglomeration handler, and
722 * define agglomerates according to the selected partitioning strategy.
730 *
Poisson<dim>::make_grid()
734 *
std::ifstream gmsh_file(std::string(MESH_DIR) +
735 *
"/unit_square_quad_unstructured.msh");
736 *
grid_in.read_msh(gmsh_file);
740 *
std::ofstream out(
"grid_input_mesh.vtu");
741 *
grid_out.write_vtu(tria, out);
744 *
tria.refine_global(5);
748 *
std::ofstream out(
"grid_fine_mesh_refined.vtu");
749 *
grid_out.write_vtu(tria, out);
753 *
std::cout <<
"Size of tria: " << tria.n_active_cells() << std::endl;
754 *
cached_tria = std::make_unique<GridTools::Cache<dim>>(tria, mapping);
755 *
ah = std::make_unique<AgglomerationHandler<dim>>(*cached_tria);
757 *
if (partitioner_type == PartitionerType::metis)
759 *
auto start = std::chrono::system_clock::now();
765 *
std::vector<typename Triangulation<dim>::active_cell_iterator>>
766 *
cells_per_subdomain(n_subdomains);
767 *
for (
const auto &cell : tria.active_cell_iterators())
770 *
for (std::size_t i = 0; i < n_subdomains; ++i)
771 *
ah->define_agglomerate(cells_per_subdomain[i]);
773 *
std::chrono::duration<double> wctduration =
774 *
(std::chrono::system_clock::now() - start);
775 *
std::cout <<
"METIS built in " << wctduration.count()
776 *
<<
" seconds [wall clock]" << std::endl;
778 *
else if (partitioner_type == PartitionerType::rtree)
782 *
static constexpr unsigned int max_elem_per_node =
783 *
PolyUtils::constexpr_pow(2, dim);
784 *
std::vector<std::pair<BoundingBox<dim>,
786 *
boxes(tria.n_active_cells());
787 *
unsigned int i = 0;
788 *
for (
const auto &cell : tria.active_cell_iterators())
789 *
boxes[i++] =
std::make_pair(mapping.get_bounding_box(cell), cell);
791 *
auto start = std::chrono::system_clock::now();
792 *
auto tree = pack_rtree<bgi::rstar<max_elem_per_node>>(boxes);
794 *
CellsAgglomerator<dim,
decltype(tree)> agglomerator{tree,
796 *
const auto vec_agglomerates = agglomerator.extract_agglomerates();
798 *
for (
const auto &agglo : vec_agglomerates)
799 *
ah->define_agglomerate(agglo);
801 *
std::chrono::duration<double> wctduration =
802 *
(std::chrono::system_clock::now() - start);
803 *
std::cout <<
"R-tree agglomerates built in " << wctduration.count()
804 *
<<
" seconds [wall clock]" << std::endl;
806 *
else if (partitioner_type == PartitionerType::no_partition)
811 *
Assert(
false, ExcMessage(
"Wrong partitioning."));
813 *
n_subdomains = ah->n_agglomerates();
814 *
std::cout <<
"N subdomains = " << n_subdomains << std::endl;
820 * To finalize the agglomeration. In the no-
partition case, each fine cell is declared as its own
821 * agglomerate. The function then distributes the degrees of freedom
822 * on the agglomerated mesh, builds the corresponding sparsity pattern,
823 * and writes a VTU file visualizing the agglomeration and the
824 * partitioning of the fine grid.
829 *
Poisson<dim>::setup_agglomeration()
831 *
if (partitioner_type == PartitionerType::no_partition)
833 *
for (
const auto &cell : tria.active_cell_iterators())
834 *
ah->define_agglomerate({cell});
837 *
ah->distribute_agglomerated_dofs(dg_fe);
838 *
ah->create_agglomeration_sparsity_pattern(dsp);
839 *
sparsity.copy_from(dsp);
842 *
std::string partitioner;
843 *
if (partitioner_type == PartitionerType::metis)
844 *
partitioner =
"metis";
845 *
else if (partitioner_type == PartitionerType::rtree)
846 *
partitioner =
"rtree";
848 *
partitioner =
"no_partitioning";
851 *
const std::string filename =
852 *
"grid_" + partitioner +
"_" + std::to_string(n_subdomains) +
".vtu";
853 *
std::ofstream output(filename);
859 *
const auto &rel = ah->get_relationships();
862 *
for (
const auto &cell : tria.active_cell_iterators())
864 *
const unsigned
int i = cell->active_cell_index();
865 *
agglo_relationships[i] = rel[i];
870 *
for (
const auto &polytope : ah->polytope_iterators())
872 *
const float id = static_cast<float>(polytope->
index());
873 *
const auto &patch_of_cells = polytope->get_agglomerate();
874 *
for (
const auto &cell : patch_of_cells)
875 *
agglo_idx[cell->active_cell_index()] = id;
878 *
data_out.add_data_vector(agglo_relationships,
879 *
"agglo_relationships",
881 *
data_out.add_data_vector(agglo_idx,
885 *
data_out.build_patches(mapping);
886 *
data_out.write_vtu(output);
894 * Assemble the global SIPG
matrix and right-hand side on the
899 * It initializes the system
matrix and right-hand side, sets up
FEValues
900 * objects on polytopal cells and interfaces, and then adds the
volume,
901 * boundary, and interior face contributions of the
symmetric interior
902 * penalty formulation.
907 *
Poisson<dim>::assemble_system()
909 *
system_matrix.reinit(sparsity);
910 *
solution.reinit(ah->n_dofs());
911 *
system_rhs.reinit(ah->n_dofs());
913 *
const unsigned int quadrature_degree = dg_fe.get_degree() + 1;
914 *
const unsigned int face_quadrature_degree = dg_fe.get_degree() + 1;
916 *
ah->initialize_fe_values(
QGauss<dim>(quadrature_degree),
922 *
const unsigned int dofs_per_cell = ah->n_dofs_per_cell();
923 *
std::cout <<
"DoFs per cell: " << dofs_per_cell << std::endl;
930 * Next, we define the four dofsxdofs matrices needed to
assemble jumps and
939 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
941 *
for (
const auto &polytope : ah->polytope_iterators())
945 *
const auto &agglo_values = ah->reinit(polytope);
946 *
polytope->get_dof_indices(local_dof_indices);
948 *
const auto &q_points = agglo_values.get_quadrature_points();
949 *
const unsigned int n_qpoints = q_points.size();
950 *
std::vector<double> rhs(n_qpoints);
951 *
rhs_function->value_list(q_points, rhs);
953 *
for (
unsigned int q_index : agglo_values.quadrature_point_indices())
955 *
for (unsigned
int i = 0; i < dofs_per_cell; ++i)
957 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
959 *
cell_matrix(i, j) += agglo_values.shape_grad(i, q_index) *
960 *
agglo_values.shape_grad(j, q_index) *
961 *
agglo_values.JxW(q_index);
963 *
cell_rhs(i) += agglo_values.shape_value(i, q_index) *
964 *
rhs[q_index] * agglo_values.JxW(q_index);
969 *
const unsigned int n_faces = polytope->n_faces();
972 *
"Invalid element: at least 4 faces are required."));
974 *
auto polygon_boundary_vertices = polytope->polytope_boundary();
975 *
for (
unsigned int f = 0; f < n_faces; ++f)
977 *
if (polytope->at_boundary(f))
979 *
const auto &fe_face = ah->reinit(polytope, f);
981 *
const unsigned int dofs_per_cell = fe_face.dofs_per_cell;
983 *
const auto &face_q_points = fe_face.get_quadrature_points();
984 *
std::vector<double> analytical_solution_values(
985 *
face_q_points.size());
986 *
analytical_solution->value_list(face_q_points,
987 *
analytical_solution_values,
990 *
const auto &normals = fe_face.get_normal_vectors();
992 *
const double penalty =
993 *
penalty_constant / std::fabs(polytope->diameter());
995 *
for (
unsigned int q_index : fe_face.quadrature_point_indices())
997 *
for (unsigned
int i = 0; i < dofs_per_cell; ++i)
999 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1002 *
(-fe_face.shape_value(i, q_index) *
1003 *
fe_face.shape_grad(j, q_index) *
1004 *
normals[q_index] -
1005 *
fe_face.shape_grad(i, q_index) * normals[q_index] *
1006 *
fe_face.shape_value(j, q_index) +
1007 *
(penalty)*fe_face.shape_value(i, q_index) *
1008 *
fe_face.shape_value(j, q_index)) *
1009 *
fe_face.JxW(q_index);
1012 *
(penalty * analytical_solution_values[q_index] *
1013 *
fe_face.shape_value(i, q_index) -
1014 *
fe_face.shape_grad(i, q_index) * normals[q_index] *
1015 *
analytical_solution_values[q_index]) *
1016 *
fe_face.JxW(q_index);
1022 *
const auto &neigh_polytope = polytope->neighbor(f);
1024 *
if (polytope->index() < neigh_polytope->index())
1026 *
unsigned int nofn =
1027 *
polytope->neighbor_of_agglomerated_neighbor(f);
1029 *
const auto &fe_faces =
1030 *
ah->reinit_interface(polytope, neigh_polytope, f, nofn);
1031 *
const auto &fe_faces0 = fe_faces.first;
1032 *
const auto &fe_faces1 = fe_faces.second;
1034 *
std::vector<types::global_dof_index>
1035 *
local_dof_indices_neighbor(dofs_per_cell);
1042 *
const auto &normals = fe_faces0.get_normal_vectors();
1044 *
const double penalty =
1045 *
penalty_constant /
std::min(polytope->diameter(), neigh_polytope->diameter());
1048 *
for (
unsigned int q_index :
1049 *
fe_faces0.quadrature_point_indices())
1051 *
for (unsigned
int i = 0; i < dofs_per_cell; ++i)
1053 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1056 *
(-0.5 * fe_faces0.shape_grad(i, q_index) *
1057 *
normals[q_index] *
1058 *
fe_faces0.shape_value(j, q_index) -
1059 *
0.5 * fe_faces0.shape_grad(j, q_index) *
1060 *
normals[q_index] *
1061 *
fe_faces0.shape_value(i, q_index) +
1062 *
(penalty)*fe_faces0.shape_value(i, q_index) *
1063 *
fe_faces0.shape_value(j, q_index)) *
1064 *
fe_faces0.JxW(q_index);
1067 *
(0.5 * fe_faces0.shape_grad(i, q_index) *
1068 *
normals[q_index] *
1069 *
fe_faces1.shape_value(j, q_index) -
1070 *
0.5 * fe_faces1.shape_grad(j, q_index) *
1071 *
normals[q_index] *
1072 *
fe_faces0.shape_value(i, q_index) -
1073 *
(penalty)*fe_faces0.shape_value(i, q_index) *
1074 *
fe_faces1.shape_value(j, q_index)) *
1075 *
fe_faces1.JxW(q_index);
1079 *
(-0.5 * fe_faces1.shape_grad(i, q_index) *
1080 *
normals[q_index] *
1081 *
fe_faces0.shape_value(j, q_index) +
1082 *
0.5 * fe_faces0.shape_grad(j, q_index) *
1083 *
normals[q_index] *
1084 *
fe_faces1.shape_value(i, q_index) -
1085 *
(penalty)*fe_faces1.shape_value(i, q_index) *
1086 *
fe_faces0.shape_value(j, q_index)) *
1087 *
fe_faces1.JxW(q_index);
1091 *
(0.5 * fe_faces1.shape_grad(i, q_index) *
1092 *
normals[q_index] *
1093 *
fe_faces1.shape_value(j, q_index) +
1094 *
0.5 * fe_faces1.shape_grad(j, q_index) *
1095 *
normals[q_index] *
1096 *
fe_faces1.shape_value(i, q_index) +
1097 *
(penalty)*fe_faces1.shape_value(i, q_index) *
1098 *
fe_faces1.shape_value(j, q_index)) *
1099 *
fe_faces1.JxW(q_index);
1104 *
neigh_polytope->get_dof_indices(local_dof_indices_neighbor);
1106 *
constraints.distribute_local_to_global(M11,
1107 *
local_dof_indices,
1109 *
constraints.distribute_local_to_global(
1111 *
local_dof_indices,
1112 *
local_dof_indices_neighbor,
1114 *
constraints.distribute_local_to_global(
1116 *
local_dof_indices_neighbor,
1117 *
local_dof_indices,
1119 *
constraints.distribute_local_to_global(
1120 *
M22, local_dof_indices_neighbor, system_matrix);
1127 * Distribute the local contributions to the global system.
1130 *
constraints.distribute_local_to_global(
1131 *
cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1137 * Solve the linear system by means of a sparse direct solver.
1140 *
template <
int dim>
1142 *
Poisson<dim>::solve()
1146 *
A_direct.vmult(solution, system_rhs);
1152 * Write VTU output and compute the global @f$L^2@f$ and @f$H^1@f$-seminorm
1153 * errors of the agglomerated DG approximation.
1156 *
template <
int dim>
1158 *
Poisson<dim>::output_results()
1161 *
std::string partitioner;
1162 *
if (partitioner_type == PartitionerType::metis)
1163 *
partitioner =
"metis";
1164 *
else if (partitioner_type == PartitionerType::rtree)
1165 *
partitioner =
"rtree";
1167 *
partitioner =
"no_partitioning";
1169 *
const std::string filename =
"interpolated_solution_" + partitioner +
"_" +
1170 *
std::to_string(n_subdomains) +
".vtu";
1171 *
std::ofstream output(filename);
1175 *
PolyUtils::interpolate_to_fine_grid(*ah,
1176 *
interpolated_solution,
1179 *
data_out.attach_dof_handler(ah->output_dh);
1180 *
data_out.add_data_vector(interpolated_solution,
1188 * Mark fine cells belonging to the same agglomerate.
1191 *
for (
const auto &polytope : ah->polytope_iterators())
1194 *
const auto &patch_of_cells = polytope->get_agglomerate();
1195 *
for (
const auto &cell : patch_of_cells)
1196 *
agglo_idx[cell->active_cell_index()] = polytope_index;
1199 *
data_out.add_data_vector(agglo_idx,
1203 *
data_out.build_patches(mapping);
1204 *
data_out.write_vtu(output);
1206 *
std::vector<double> errors;
1207 *
PolyUtils::compute_global_error(*ah,
1209 *
*analytical_solution,
1213 *
l2_err = errors[0];
1214 *
semih1_err = errors[1];
1221 * Return the number of degrees of freedom on the agglomerated mesh.
1224 *
template <
int dim>
1226 *
Poisson<dim>::get_n_dofs() const
1228 *
return ah->n_dofs();
1234 * Return the pair consisting of the @f$L^2@f$ error and the @f$H^1@f$-seminorm error of the numerical solution.
1237 *
template <
int dim>
1238 *
inline std::pair<double, double>
1239 *
Poisson<dim>::get_error() const
1241 *
return std::make_pair(l2_err, semih1_err);
1247 * Run the full workflow: mesh generation, agglomeration setup,
1248 * assembly, solution, and postprocessing.
1251 *
template <
int dim>
1253 *
Poisson<dim>::run()
1256 *
setup_agglomeration();
1257 *
auto start = std::chrono::high_resolution_clock::now();
1258 *
assemble_system();
1259 *
auto stop = std::chrono::high_resolution_clock::now();
1261 *
std::chrono::duration_cast<std::chrono::seconds>(stop - start);
1263 *
std::cout <<
"Time taken by assemble_system(): " << duration.count()
1264 *
<<
" seconds" << std::endl;
1277 *
ConvergenceInfo convergence_info;
1279 *
for (
unsigned int fe_degree : {1})
1281 *
std::cout <<
"Running with FE degree: " << fe_degree << std::endl;
1282 *
Poisson<2> poisson_problem{PartitionerType::rtree,
1286 *
poisson_problem.run();
1287 *
convergence_info.add(
1289 *
poisson_problem.get_n_dofs(), poisson_problem.get_error()));
1290 *
std::cout << std::endl;
1293 *
std::cout <<
"Convergence table:" << std::endl;
1294 *
convergence_info.print();
1295 *
std::cout << std::endl;
1302<a name=
"ann-include/agglomeration_accessor.h"></a>
1303<h1>Annotated version of include/agglomeration_accessor.h</h1>
1320 *
#ifndef agglomeration_accessor_h
1321 *
#define agglomeration_accessor_h
1323 *
#include <deal.II/base/config.h>
1325 *
#include <deal.II/base/bounding_box.h>
1326 *
#include <deal.II/base/iterator_range.h>
1328 *
#include <deal.II/grid/filtered_iterator.h>
1332 *
using namespace dealii;
1337 * Forward declarations
1341 *
template <
int,
int>
1342 *
class AgglomerationHandler;
1343 *
template <
int,
int>
1344 *
class AgglomerationIterator;
1351 *
template <
int dim,
int spacedim = dim>
1352 *
class AgglomerationAccessor
1358 *
using AgglomerationContainer =
1359 *
std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>;
1366 *
get_dof_indices(std::vector<types::global_dof_index> &)
const;
1381 *
n_agglomerated_faces()
const;
1386 *
const AgglomerationIterator<dim, spacedim>
1387 *
neighbor(
const unsigned int f)
const;
1394 *
neighbor_of_agglomerated_neighbor(
const unsigned int f)
const;
1406 *
at_boundary(
const unsigned int f)
const;
1411 *
const std::vector<typename Triangulation<dim>::active_face_iterator> &
1412 *
polytope_boundary()
const;
1430 *
AgglomerationContainer
1431 *
get_agglomerate()
const;
1439 *
get_bounding_box()
const;
1461 *
n_background_cells()
const;
1468 *
is_locally_owned()
const;
1488 *
inline const std::vector<types::global_cell_index> &
1513 *
active_fe_index()
const;
1520 *
AgglomerationAccessor();
1527 *
AgglomerationAccessor(
1530 *
const AgglomerationHandler<dim, spacedim> *ah);
1535 *
AgglomerationAccessor(
1538 *
const AgglomerationHandler<dim, spacedim> *ah);
1543 *
~AgglomerationAccessor() =
default;
1569 *
AgglomerationHandler<dim, spacedim> *handler;
1576 *
operator==(
const AgglomerationAccessor<dim, spacedim> &other)
const;
1582 *
operator!=(
const AgglomerationAccessor<dim, spacedim> &other)
const;
1599 *
const AgglomerationContainer &
1600 *
get_slaves()
const;
1603 *
n_agglomerated_faces_per_cell(
1607 *
template <
int,
int>
1608 *
friend class AgglomerationIterator;
1613 *
template <
int dim,
int spacedim>
1615 *
AgglomerationAccessor<dim, spacedim>::n_agglomerated_faces_per_cell(
1618 *
unsigned int n_neighbors = 0;
1619 *
for (
const auto &f : cell->face_indices())
1621 *
const auto &neighboring_cell = cell->neighbor(f);
1622 *
if ((cell->face(f)->at_boundary()) ||
1623 *
(neighboring_cell->is_active() &&
1624 *
!handler->are_cells_agglomerated(cell, neighboring_cell)))
1629 *
return n_neighbors;
1634 *
template <
int dim,
int spacedim>
1636 *
AgglomerationAccessor<dim, spacedim>::n_faces() const
1638 *
Assert(!handler->is_slave_cell(master_cell),
1639 *
ExcMessage(
"You cannot pass a slave cell."));
1640 *
return handler->number_of_agglomerated_faces[present_index];
1645 *
template <
int dim,
int spacedim>
1646 *
const AgglomerationIterator<dim, spacedim>
1647 *
AgglomerationAccessor<dim, spacedim>::neighbor(
const unsigned int f)
const
1649 *
if (!at_boundary(f))
1651 *
if (master_cell->is_ghost())
1655 * The following path is needed when the present function is called
1656 * from neighbor_of_neighbor()
1662 *
const unsigned int sender_rank = master_cell->subdomain_id();
1664 *
const CellId &master_id_ghosted_neighbor =
1665 *
handler->recv_ghosted_master_id.at(sender_rank)
1671 * Use the
id of the master cell to uniquely identify the neighboring
1678 *
return {master_cell,
1679 *
master_id_ghosted_neighbor,
1684 *
handler->master2polygon.at(master_cell->active_cell_index());
1686 *
const auto &neigh =
1687 *
handler->polytope_cache.cell_face_at_boundary.at({polytope_index, f})
1691 *
if (neigh->is_locally_owned())
1694 *
*neigh, &(handler->agglo_dh));
1695 *
return {cell_dh, handler};
1701 * Get master_id from the neighboring ghost polytope. This uniquely
1702 * identifies the neighboring polytope among all processors.
1705 *
const CellId &master_id_neighbor =
1706 *
handler->polytope_cache.ghosted_master_id.at({present_id, f});
1710 * Use the
id of the master cell to uniquely identify the neighboring
1714 *
return {neigh, master_id_neighbor, handler};
1725 *
template <
int dim,
int spacedim>
1727 *
AgglomerationAccessor<dim, spacedim>::neighbor_of_agglomerated_neighbor(
1728 *
const unsigned int f)
const
1732 * First, make sure it
's not a boundary face.
1735 * if (!at_boundary(f))
1737 * const auto &neigh_polytope =
1738 * neighbor(f); // returns the neighboring master and id
1740 * AssertThrow(neigh_polytope.state() == IteratorState::valid,
1741 * ExcInternalError());
1743 * unsigned int n_faces_agglomerated_neighbor;
1747 * if it is locally owned, retrieve the number of faces
1750 * if (neigh_polytope->is_locally_owned())
1752 * n_faces_agglomerated_neighbor = neigh_polytope->n_faces();
1758 * The neighboring polytope is not locally owned. We need to get the
1759 * number of its faces from the neighboring rank.
1763 * First, retrieve the CellId of the neighboring polytope.
1766 * const CellId &master_id_neighbor = neigh_polytope->id();
1770 * Then, get the neighboring rank
1773 * const unsigned int sender_rank = neigh_polytope->subdomain_id();
1777 * From the neighboring rank, use the CellId of the neighboring
1778 * polytope to get the number of its faces.
1781 * n_faces_agglomerated_neighbor =
1782 * handler->recv_n_faces.at(sender_rank).at(master_id_neighbor);
1788 * Loop over all faces of neighboring agglomerate
1791 * for (unsigned int f_out = 0; f_out < n_faces_agglomerated_neighbor;
1796 * Check if same CellId
1799 * if (neigh_polytope->neighbor(f_out).state() == IteratorState::valid)
1800 * if (neigh_polytope->neighbor(f_out)->id() == present_id)
1803 * return numbers::invalid_unsigned_int;
1809 * Face is at boundary
1812 * return numbers::invalid_unsigned_int;
1818 * ------------------------------ inline functions -------------------------
1824 * template <int dim, int spacedim>
1825 * inline AgglomerationAccessor<dim, spacedim>::AgglomerationAccessor()
1830 * template <int dim, int spacedim>
1831 * inline AgglomerationAccessor<dim, spacedim>::AgglomerationAccessor(
1832 * const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
1833 * const AgglomerationHandler<dim, spacedim> *ah)
1835 * handler = const_cast<AgglomerationHandler<dim, spacedim> *>(ah);
1836 * if (&(*handler->master_cells_container.end()) == std::addressof(cell))
1838 * present_index = handler->master_cells_container.size();
1839 * master_cell = *handler->master_cells_container.end();
1840 * present_id = CellId(); // invalid id (TODO)
1841 * present_subdomain_id = numbers::invalid_subdomain_id;
1845 * present_index = handler->master2polygon.at(cell->active_cell_index());
1846 * master_cell = cell;
1847 * present_id = master_cell->id();
1848 * present_subdomain_id = master_cell->subdomain_id();
1854 * template <int dim, int spacedim>
1855 * inline AgglomerationAccessor<dim, spacedim>::AgglomerationAccessor(
1856 * const typename Triangulation<dim, spacedim>::active_cell_iterator &neigh_cell,
1857 * const CellId &master_cell_id,
1858 * const AgglomerationHandler<dim, spacedim> *ah)
1860 * Assert(neigh_cell->is_ghost(), ExcInternalError());
1863 * neigh_cell is ghosted
1869 * handler = const_cast<AgglomerationHandler<dim, spacedim> *>(ah);
1870 * master_cell = neigh_cell;
1871 * present_index = numbers::invalid_unsigned_int;
1874 * neigh_cell is ghosted, use the CellId of that agglomerate
1877 * present_id = master_cell_id;
1878 * present_subdomain_id = master_cell->subdomain_id();
1883 * template <int dim, int spacedim>
1885 * AgglomerationAccessor<dim, spacedim>::get_dof_indices(
1886 * std::vector<types::global_dof_index> &dof_indices) const
1888 * Assert(dof_indices.size() > 0,
1890 * "The vector of DoFs indices must be already properly resized."));
1891 * if (is_locally_owned())
1895 * Forward the call to the master cell
1898 * typename DoFHandler<dim, spacedim>::cell_iterator master_cell_dh(
1899 * *master_cell, &(handler->agglo_dh));
1900 * master_cell_dh->get_dof_indices(dof_indices);
1904 * const std::vector<types::global_dof_index> &recv_dof_indices =
1905 * handler->recv_ghost_dofs.at(present_subdomain_id).at(present_id);
1907 * std::copy(recv_dof_indices.cbegin(),
1908 * recv_dof_indices.cend(),
1909 * dof_indices.begin());
1915 * template <int dim, int spacedim>
1916 * inline typename AgglomerationAccessor<dim, spacedim>::AgglomerationContainer
1917 * AgglomerationAccessor<dim, spacedim>::get_agglomerate() const
1919 * auto agglomeration = get_slaves();
1920 * agglomeration.push_back(master_cell);
1921 * return agglomeration;
1926 * template <int dim, int spacedim>
1927 * inline const std::vector<typename Triangulation<dim>::active_face_iterator> &
1928 * AgglomerationAccessor<dim, spacedim>::polytope_boundary() const
1930 * return handler->polygon_boundary[master_cell];
1935 * template <int dim, int spacedim>
1937 * AgglomerationAccessor<dim, spacedim>::diameter() const
1939 * Assert(!handler->is_slave_cell(master_cell),
1940 * ExcMessage("The present function cannot be called for slave cells."));
1942 * if (handler->is_master_cell(master_cell))
1946 * Get the bounding box associated with the master cell
1949 * const auto &bdary_pts =
1950 * handler->bboxes[present_index].get_boundary_points();
1951 * return (bdary_pts.second - bdary_pts.first).norm();
1957 * Standard deal.II way to get the measure of a cell.
1960 * return master_cell->diameter();
1966 * template <int dim, int spacedim>
1967 * inline const BoundingBox<dim> &
1968 * AgglomerationAccessor<dim, spacedim>::get_bounding_box() const
1970 * if (is_locally_owned())
1971 * return handler->bboxes[present_index];
1973 * return handler->recv_ghosted_bbox.at(present_subdomain_id).at(present_id);
1978 * template <int dim, int spacedim>
1980 * AgglomerationAccessor<dim, spacedim>::volume() const
1982 * Assert(!handler->is_slave_cell(master_cell),
1983 * ExcMessage("The present function cannot be called for slave cells."));
1985 * if (handler->is_master_cell(master_cell))
1987 * return handler->bboxes[present_index].volume();
1991 * return master_cell->measure();
1997 * template <int dim, int spacedim>
1999 * AgglomerationAccessor<dim, spacedim>::next()
2003 * Increment the present index and update the polytope
2010 * Make sure not to query the CellId if it's past the last
2013 *
if (present_index < handler->master_cells_container.size())
2015 *
master_cell = handler->master_cells_container[present_index];
2016 *
present_id = master_cell->id();
2017 *
present_subdomain_id = master_cell->subdomain_id();
2023 *
template <
int dim,
int spacedim>
2025 *
AgglomerationAccessor<dim, spacedim>::prev()
2029 * Decrement the present
index and update the polytope
2033 *
master_cell = handler->master_cells_container[present_index];
2034 *
present_id = master_cell->id();
2038 *
template <
int dim,
int spacedim>
2040 *
AgglomerationAccessor<dim, spacedim>::operator==(
2041 *
const AgglomerationAccessor<dim, spacedim> &other)
const
2043 *
return present_index == other.present_index;
2046 *
template <
int dim,
int spacedim>
2048 *
AgglomerationAccessor<dim, spacedim>::operator!=(
2049 *
const AgglomerationAccessor<dim, spacedim> &other)
const
2051 *
return !(*
this == other);
2056 *
template <
int dim,
int spacedim>
2058 *
AgglomerationAccessor<dim, spacedim>::index() const
2060 *
return present_index;
2065 *
template <
int dim,
int spacedim>
2067 *
AgglomerationAccessor<dim, spacedim>::as_dof_handler_iterator(
2072 * Forward the call to the master cell
using the right
DoFHandler.
2075 *
return master_cell->as_dof_handler_iterator(dof_handler);
2080 *
template <
int dim,
int spacedim>
2081 *
inline const typename AgglomerationAccessor<dim,
2082 *
spacedim>::AgglomerationContainer &
2083 *
AgglomerationAccessor<dim, spacedim>::get_slaves() const
2085 *
return handler->master2slaves.at(master_cell->active_cell_index());
2090 *
template <
int dim,
int spacedim>
2091 *
inline unsigned int
2092 *
AgglomerationAccessor<dim, spacedim>::n_background_cells() const
2094 *
AssertThrow(get_agglomerate().
size() > 0, ExcMessage(
"Empty agglomeration."));
2095 *
return get_agglomerate().size();
2100 *
template <
int dim,
int spacedim>
2102 *
AgglomerationAccessor<dim, spacedim>::n_agglomerated_faces() const
2104 *
const auto &agglomeration = get_agglomerate();
2105 *
unsigned int n_neighbors = 0;
2106 *
for (
const auto &cell : agglomeration)
2107 *
n_neighbors += n_agglomerated_faces_per_cell(cell);
2108 *
return n_neighbors;
2113 *
template <
int dim,
int spacedim>
2115 *
AgglomerationAccessor<dim, spacedim>::at_boundary(
const unsigned int f)
const
2117 *
if (master_cell->is_ghost())
2119 *
const unsigned int sender_rank = master_cell->subdomain_id();
2120 *
return handler->recv_bdary_info.at(sender_rank).at(present_id).at(f);
2124 *
Assert(!handler->is_slave_cell(master_cell),
2126 *
"This function should not be called for a slave cell."));
2130 *
*master_cell, &(handler->agglo_dh));
2131 *
return handler->at_boundary(cell_dh, f);
2137 *
template <
int dim,
int spacedim>
2139 *
AgglomerationAccessor<dim, spacedim>::is_locally_owned() const
2141 *
return master_cell->is_locally_owned();
2146 *
template <
int dim,
int spacedim>
2148 *
AgglomerationAccessor<dim, spacedim>::id() const
2150 *
return present_id;
2155 *
template <
int dim,
int spacedim>
2157 *
AgglomerationAccessor<dim, spacedim>::subdomain_id() const
2159 *
return present_subdomain_id;
2162 *
template <
int dim,
int spacedim>
2163 *
inline const std::vector<types::global_cell_index> &
2164 *
AgglomerationAccessor<dim, spacedim>::children() const
2166 *
Assert(!handler->parent_child_info.empty(), ExcInternalError());
2167 *
return handler->parent_child_info.at(
2168 *
{present_index, handler->present_extraction_level});
2171 *
template <
int dim,
int spacedim>
2173 *
AgglomerationAccessor<dim, spacedim>::get_fe() const
2176 *
master_cell_as_dof_handler_iterator =
2177 *
master_cell->as_dof_handler_iterator(handler->agglo_dh);
2178 *
return master_cell_as_dof_handler_iterator->get_fe();
2181 *
template <
int dim,
int spacedim>
2183 *
AgglomerationAccessor<dim, spacedim>::set_active_fe_index(
2186 *
Assert(!handler->is_slave_cell(master_cell),
2187 *
ExcMessage(
"The present function cannot be called for slave cells."));
2189 *
master_cell_as_dof_handler_iterator =
2190 *
master_cell->as_dof_handler_iterator(handler->agglo_dh);
2191 *
master_cell_as_dof_handler_iterator->set_active_fe_index(index);
2194 *
template <
int dim,
int spacedim>
2196 *
AgglomerationAccessor<dim, spacedim>::active_fe_index() const
2199 *
master_cell_as_dof_handler_iterator =
2200 *
master_cell->as_dof_handler_iterator(handler->agglo_dh);
2201 *
return master_cell_as_dof_handler_iterator->active_fe_index();
2208<a name=
"ann-include/agglomeration_handler.h"></a>
2209<h1>Annotated version of include/agglomeration_handler.h</h1>
2226 *
#ifndef agglomeration_handler_h
2227 *
#define agglomeration_handler_h
2229 *
#include <deal.II/base/mpi.h>
2230 *
#include <deal.II/base/quadrature.h>
2231 *
#include <deal.II/base/enable_observer_pointer.h>
2233 *
#include <deal.II/distributed/shared_tria.h>
2234 *
#include <deal.II/distributed/tria.h>
2236 *
#include <deal.II/dofs/dof_handler.h>
2237 *
#include <deal.II/dofs/dof_tools.h>
2239 *
#include <deal.II/fe/fe_dgp.h>
2240 *
#include <deal.II/fe/fe_dgq.h>
2241 *
#include <deal.II/fe/fe_nothing.h>
2242 *
#include <deal.II/fe/fe_simplex_p.h>
2243 *
#include <deal.II/fe/fe_system.h>
2244 *
#include <deal.II/fe/fe_values.h>
2245 *
#include <deal.II/fe/mapping_fe_field.h>
2246 *
#include <deal.II/fe/mapping_q.h>
2248 *
#include <deal.II/grid/grid_tools_cache.h>
2249 *
#include <deal.II/grid/tria.h>
2251 *
#include <deal.II/
hp/fe_collection.h>
2253 *
#include <deal.II/lac/dynamic_sparsity_pattern.h>
2254 *
#include <deal.II/lac/la_parallel_vector.h>
2255 *
#include <deal.II/lac/sparse_matrix.h>
2256 *
#include <deal.II/lac/sparsity_pattern.h>
2257 *
#include <deal.II/lac/trilinos_sparse_matrix.h>
2258 *
#include <deal.II/lac/vector.h>
2260 *
#include <deal.II/meshworker/scratch_data.h>
2262 *
#include <deal.II/non_matching/fe_immersed_values.h>
2263 *
#include <deal.II/non_matching/immersed_surface_quadrature.h>
2265 *
#include <agglomeration_iterator.h>
2266 *
#include <agglomerator.h>
2267 *
#include <mapping_box.h>
2269 *
#include <fstream>
2272 *
using namespace dealii;
2276 * Forward declarations
2279 *
template <
int dim,
int spacedim>
2280 *
class AgglomerationHandler;
2289 *
template <
int,
int>
2290 *
class AgglomerationHandlerImplementation;
2304 *
template <
int dim,
int spacedim>
2305 *
class PolytopeCache
2311 *
PolytopeCache() =
default;
2316 *
~PolytopeCache() =
default;
2323 * clear all the members
2326 *
cell_face_at_boundary.clear();
2327 *
interface.clear();
2328 *
visited_cell_and_faces.clear();
2338 *
mutable std::set<std::pair<types::global_cell_index, unsigned int>>
2339 *
visited_cell_and_faces;
2342 *
mutable std::set<std::pair<CellId, unsigned int>>
2343 *
visited_cell_and_faces_id;
2356 *
std::pair<types::global_cell_index, unsigned int>,
2359 *
cell_face_at_boundary;
2365 *
mutable std::map<std::pair<CellId, unsigned int>,
CellId>
2366 *
ghosted_master_id;
2379 *
std::pair<CellId, CellId>,
2381 *
std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
2392 *
template <
int dim,
int spacedim = dim>
2396 *
using agglomeration_iterator = AgglomerationIterator<dim, spacedim>;
2398 *
using AgglomerationContainer =
2399 *
typename AgglomerationIterator<dim, spacedim>::AgglomerationContainer;
2402 *
enum CellAgglomerationType
2410 *
explicit AgglomerationHandler(
2413 *
AgglomerationHandler() =
default;
2415 *
~AgglomerationHandler()
2419 * disconnect the signal
2422 *
tria_listener.disconnect();
2428 *
agglomeration_iterator
2434 *
agglomeration_iterator
2440 *
agglomeration_iterator
2446 *
agglomeration_iterator
2452 *
agglomeration_iterator
2460 *
polytope_iterators()
const;
2462 *
template <
int,
int>
2463 *
friend class AgglomerationIterator;
2465 *
template <
int,
int>
2466 *
friend class AgglomerationAccessor;
2479 *
distribute_agglomerated_dofs(
2488 *
initialize_fe_values(
2498 *
initialize_fe_values(
2512 *
template <
typename SparsityPatternType,
typename Number =
double>
2514 *
create_agglomeration_sparsity_pattern(
2515 *
SparsityPatternType &sparsity_pattern,
2517 *
const bool keep_constrained_dofs =
true,
2528 *
agglomeration_iterator
2529 *
define_agglomerate(
const AgglomerationContainer &cells);
2541 *
agglomeration_iterator
2542 *
define_agglomerate(
const AgglomerationContainer &cells,
2543 *
const unsigned int fecollection_size);
2547 *
get_triangulation()
const;
2553 *
get_mapping()
const;
2555 *
inline const MappingBox<dim> &
2556 *
get_agglomeration_mapping()
const;
2558 *
inline const std::vector<BoundingBox<dim>> &
2559 *
get_local_bboxes()
const;
2566 *
get_mesh_size()
const;
2569 *
cell_to_polytope_index(
2573 *
inline decltype(
auto)
2574 *
get_interface()
const;
2579 *
template <
typename CellIterator>
2581 *
is_master_cell(
const CellIterator &cell)
const;
2587 *
inline const std::vector<
2593 *
get_relationships()
const;
2601 *
inline std::vector<
2605 *
&master_cell)
const;
2614 *
get_dof_handler()
const;
2620 *
n_agglomerates()
const;
2626 *
n_agglomerated_faces_per_cell(
2634 *
reinit(
const AgglomerationIterator<dim, spacedim> &polytope)
const;
2641 *
reinit(
const AgglomerationIterator<dim, spacedim> &polytope,
2642 *
const unsigned int face_index)
const;
2649 *
std::pair<const FEValuesBase<dim, spacedim> &,
2651 *
reinit_interface(
const AgglomerationIterator<dim, spacedim> &polytope_in,
2652 *
const AgglomerationIterator<dim, spacedim> &neigh_polytope,
2653 *
const unsigned int local_in,
2654 *
const unsigned int local_outside)
const;
2662 *
agglomerated_quadrature(
2663 *
const AgglomerationContainer &cells,
2665 *
&master_cell)
const;
2680 *
const unsigned int f)
const;
2682 *
inline unsigned int
2683 *
n_dofs_per_cell() const noexcept;
2685 *
inline
types::global_dof_index
2686 *
n_dofs() const noexcept;
2695 *
inline const
std::vector<typename
Triangulation<dim>::active_face_iterator> &
2696 *
polytope_boundary(
2697 *
const typename
Triangulation<dim>::active_cell_iterator &cell);
2710 *
std::unique_ptr<MappingBox<dim>> box_mapping;
2721 *
setup_ghost_polytopes();
2724 *
exchange_interface_values();
2728 * TODO: move it to private interface
2732 *
types::subdomain_id,
2737 *
types::subdomain_id,
2742 *
types::subdomain_id,
2747 *
types::subdomain_id,
2751 *
mutable
std::map<
types::subdomain_id,
2760 *
inline const typename
DoFHandler<dim, spacedim>::active_cell_iterator
2761 *
polytope_to_dh_iterator(const
types::global_cell_index polytope_index) const;
2766 *
template <typename RtreeType>
2768 *
connect_hierarchy(const CellsAgglomerator<dim, RtreeType> &agglomerator);
2774 *
inline const
hp::FECollection<dim, spacedim> &
2775 *
get_fe_collection() const;
2781 *
used_fe_collection() const;
2788 *
initialize_agglomeration_data(
2789 *
const
std::unique_ptr<
GridTools::Cache<dim, spacedim>> &cache_tria);
2792 *
update_agglomerate(
2793 *
AgglomerationContainer &polytope,
2794 *
const typename
Triangulation<dim, spacedim>::active_cell_iterator
2801 *
connect_to_tria_signals()
2805 * First disconnect existing connections
2808 *
tria_listener.disconnect();
2809 *
tria_listener = tria->signals.any_change.connect(
2810 *
[&]() { this->initialize_agglomeration_data(this->cached_tria); });
2828 *
create_bounding_box(
const AgglomerationContainer &polytope);
2832 *
get_master_idx_of_cell(
2840 *
are_cells_agglomerated(
2843 *
&other_cell)
const;
2854 *
initialize_hp_structure();
2863 *
const unsigned int face_number,
2865 *
&agglo_isv_ptr)
const;
2872 *
template <
typename CellIterator>
2874 *
is_slave_cell(
const CellIterator &cell)
const;
2882 *
setup_connectivity_of_agglomeration();
2888 *
unsigned int n_agglomerations;
2905 *
master_slave_relationships_iterators;
2909 *
mutable std::vector<types::global_cell_index> number_of_agglomerated_faces;
2918 *
std::vector<typename Triangulation<dim>::active_face_iterator>>
2928 *
std::vector<BoundingBox<spacedim>> bboxes;
2942 *
mutable std::map<types::subdomain_id, std::map<CellId, unsigned int>>
2945 *
mutable std::map<types::subdomain_id, std::map<CellId, unsigned int>>
2951 *
CellId (including slaves)
2954 *
mutable std::map<types::subdomain_id, std::map<CellId, CellId>>
2955 *
local_cell_ids_neigh_cell;
2957 *
mutable std::map<types::subdomain_id, std::map<CellId, CellId>>
2958 *
recv_cell_ids_neigh_cell;
2963 * send to neighborign rank the information that
2964 * - current polytope
id
2966 * has the following neighboring
id.
2970 *
std::map<CellId, std::map<unsigned int, CellId>>>
2971 *
local_ghosted_master_id;
2974 *
std::map<CellId, std::map<unsigned int, CellId>>>
2975 *
recv_ghosted_master_id;
2979 * CellIds from neighboring rank
2983 *
std::map<CellId, std::map<unsigned int, bool>>>
2987 *
std::map<CellId, std::map<unsigned int, bool>>>
2992 * Exchange neighboring bounding boxes
2995 *
mutable std::map<types::subdomain_id, std::map<CellId, BoundingBox<dim>>>
2996 *
local_ghosted_bbox;
2998 *
mutable std::map<types::subdomain_id, std::map<CellId, BoundingBox<dim>>>
2999 *
recv_ghosted_bbox;
3003 * Exchange DoF indices with ghosted polytopes
3007 *
std::map<CellId, std::vector<types::global_dof_index>>>
3011 *
std::map<CellId, std::vector<types::global_dof_index>>>
3021 *
std::map<std::pair<CellId, unsigned int>, std::vector<Point<spacedim>>>>
3031 *
std::map<std::pair<CellId, unsigned int>, std::vector<double>>>
3041 *
std::map<std::pair<CellId, unsigned int>, std::vector<Tensor<1, spacedim>>>>
3051 *
std::map<std::pair<CellId, unsigned int>, std::vector<std::vector<double>>>>
3055 *
std::map<std::pair<CellId, unsigned int>,
3056 *
std::vector<std::vector<Tensor<1, spacedim>>>>>
3073 *
std::unique_ptr<GridTools::Cache<dim, spacedim>> cached_tria;
3083 *
std::unique_ptr<FiniteElement<dim>> fe;
3097 *
mutable std::unique_ptr<ScratchData> standard_scratch;
3104 *
mutable std::unique_ptr<ScratchData> agglomerated_scratch;
3107 *
mutable std::unique_ptr<NonMatching::FEImmersedSurfaceValues<spacedim>>
3110 *
mutable std::unique_ptr<NonMatching::FEImmersedSurfaceValues<spacedim>>
3111 *
agglomerated_isv_neigh;
3113 *
mutable std::unique_ptr<NonMatching::FEImmersedSurfaceValues<spacedim>>
3114 *
agglomerated_isv_bdary;
3116 *
boost::signals2::connection tria_listener;
3120 *
const UpdateFlags internal_agglomeration_flags =
3126 *
const UpdateFlags internal_agglomeration_face_flags =
3132 *
Quadrature<dim - 1> agglomeration_face_quad;
3136 * Associate the master cell to the slaves.
3139 *
std::unordered_map<
3141 *
std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3146 * Map the master cell
index with the polytope
index
3149 *
std::map<types::global_cell_index, types::global_cell_index> master2polygon;
3152 *
std::vector<typename Triangulation<dim>::active_cell_iterator>
3153 *
master_disconnected;
3157 * Dummy
FiniteElement objects needed only to generate quadratures
3171 *
std::unique_ptr<FEValues<dim, spacedim>> no_values;
3176 *
std::unique_ptr<FEFaceValues<dim, spacedim>> no_face_values;
3181 *
std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
3182 *
master_cells_container;
3184 *
friend class internal::AgglomerationHandlerImplementation<dim, spacedim>;
3186 *
internal::PolytopeCache<dim, spacedim> polytope_cache;
3194 *
std::map<std::pair<types::global_cell_index, types::global_cell_index>,
3195 *
std::vector<types::global_cell_index>>
3196 *
parent_child_info;
3198 *
unsigned int present_extraction_level;
3205 *
bool is_hp_collection =
false;
3206 *
std::unique_ptr<hp::FECollection<dim, spacedim>>
3211 * Stores quadrature rules; these QCollections should have the same
size as
3219 *
mapping_collection;
3221 *
dummy_fe_collection;
3224 * containing only dummy_fe
3226 * actually contain only
one element each.
3230 * Analogous to no_values and no_face_values, but used when different cells
3231 * employ different FEs or quadratures
3234 *
std::unique_ptr<hp::FEValues<dim, spacedim>> hp_no_values;
3235 *
std::unique_ptr<hp::FEFaceValues<dim, spacedim>> hp_no_face_values;
3242 * ------------------------------
inline functions -------------------------
3245 *
template <
int dim,
int spacedim>
3247 *
AgglomerationHandler<dim, spacedim>::get_fe() const
3254 *
template <
int dim,
int spacedim>
3256 *
AgglomerationHandler<dim, spacedim>::get_mapping() const
3263 *
template <
int dim,
int spacedim>
3264 *
inline const MappingBox<dim> &
3265 *
AgglomerationHandler<dim, spacedim>::get_agglomeration_mapping() const
3267 *
return *box_mapping;
3272 *
template <
int dim,
int spacedim>
3274 *
AgglomerationHandler<dim, spacedim>::get_triangulation() const
3280 *
template <
int dim,
int spacedim>
3281 *
inline const std::vector<BoundingBox<dim>> &
3282 *
AgglomerationHandler<dim, spacedim>::get_local_bboxes() const
3289 *
template <
int dim,
int spacedim>
3291 *
AgglomerationHandler<dim, spacedim>::cell_to_polytope_index(
3294 *
return master2polygon.at(cell->active_cell_index());
3299 *
template <
int dim,
int spacedim>
3300 *
inline decltype(
auto)
3301 *
AgglomerationHandler<dim, spacedim>::get_interface() const
3303 *
return polytope_cache.interface;
3308 *
template <
int dim,
int spacedim>
3310 *
AgglomerationHandler<dim, spacedim>::get_relationships() const
3312 *
return master_slave_relationships;
3317 *
template <
int dim,
int spacedim>
3318 *
inline std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
3319 *
AgglomerationHandler<dim, spacedim>::get_agglomerate(
3321 *
&master_cell)
const
3323 *
Assert(is_master_cell(master_cell), ExcInternalError());
3324 *
auto agglomeration = get_slaves_of_idx(master_cell->active_cell_index());
3325 *
agglomeration.push_back(master_cell);
3326 *
return agglomeration;
3331 *
template <
int dim,
int spacedim>
3333 *
AgglomerationHandler<dim, spacedim>::get_dof_handler() const
3340 *
template <
int dim,
int spacedim>
3341 *
inline const std::vector<
3343 *
AgglomerationHandler<dim, spacedim>::get_slaves_of_idx(
3346 *
return master2slaves.at(idx);
3351 *
template <
int dim,
int spacedim>
3352 *
template <
typename CellIterator>
3354 *
AgglomerationHandler<dim, spacedim>::is_master_cell(
3355 *
const CellIterator &cell)
const
3357 *
return master_slave_relationships[cell->global_active_cell_index()] == -1;
3366 *
template <
int dim,
int spacedim>
3367 *
template <
typename CellIterator>
3369 *
AgglomerationHandler<dim, spacedim>::is_slave_cell(
3370 *
const CellIterator &cell)
const
3372 *
return master_slave_relationships[cell->global_active_cell_index()] >= 0;
3377 *
template <
int dim,
int spacedim>
3379 *
AgglomerationHandler<dim, spacedim>::at_boundary(
3381 *
const unsigned int face_index)
const
3383 *
Assert(!is_slave_cell(cell),
3384 *
ExcMessage(
"This function should not be called for a slave cell."));
3386 *
return polytope_cache.cell_face_at_boundary
3387 *
.at({master2polygon.at(cell->active_cell_index()), face_index})
3392 *
template <
int dim,
int spacedim>
3393 *
inline unsigned int
3394 *
AgglomerationHandler<dim, spacedim>::n_dofs_per_cell() const noexcept
3396 *
return fe->n_dofs_per_cell();
3401 *
template <
int dim,
int spacedim>
3403 *
AgglomerationHandler<dim, spacedim>::n_dofs() const noexcept
3405 *
return agglo_dh.n_dofs();
3410 *
template <
int dim,
int spacedim>
3411 *
inline const std::vector<typename Triangulation<dim>::active_face_iterator> &
3412 *
AgglomerationHandler<dim, spacedim>::polytope_boundary(
3415 *
return polygon_boundary[cell];
3420 *
template <
int dim,
int spacedim>
3422 *
AgglomerationHandler<dim, spacedim>::is_slave_cell_of(
3425 *
return master_slave_relationships_iterators.at(cell->active_cell_index());
3430 *
template <
int dim,
int spacedim>
3432 *
AgglomerationHandler<dim, spacedim>::get_master_idx_of_cell(
3435 *
auto idx = master_slave_relationships[cell->global_active_cell_index()];
3437 *
return cell->global_active_cell_index();
3444 *
template <
int dim,
int spacedim>
3446 *
AgglomerationHandler<dim, spacedim>::are_cells_agglomerated(
3453 *
if different subdomain, then **by construction** they will not be together
3454 *
if (cell->subdomain_id() != other_cell->subdomain_id())
3459 *
return (get_master_idx_of_cell(cell) == get_master_idx_of_cell(other_cell));
3464 *
template <
int dim,
int spacedim>
3465 *
inline unsigned int
3466 *
AgglomerationHandler<dim, spacedim>::n_agglomerates() const
3468 *
return n_agglomerations;
3473 *
template <
int dim,
int spacedim>
3475 *
AgglomerationHandler<dim, spacedim>::polytope_to_dh_iterator(
3478 *
return master_cells_container[polytope_index]->as_dof_handler_iterator(
3484 *
template <
int dim,
int spacedim>
3485 *
AgglomerationIterator<dim, spacedim>
3486 *
AgglomerationHandler<dim, spacedim>::begin() const
3488 *
Assert(n_agglomerations > 0,
3489 *
ExcMessage(
"No agglomeration has been performed."));
3490 *
return {*master_cells_container.begin(),
this};
3495 *
template <
int dim,
int spacedim>
3496 *
AgglomerationIterator<dim, spacedim>
3497 *
AgglomerationHandler<dim, spacedim>::begin()
3499 *
Assert(n_agglomerations > 0,
3500 *
ExcMessage(
"No agglomeration has been performed."));
3501 *
return {*master_cells_container.begin(),
this};
3506 *
template <
int dim,
int spacedim>
3507 *
AgglomerationIterator<dim, spacedim>
3508 *
AgglomerationHandler<dim, spacedim>::end() const
3510 *
Assert(n_agglomerations > 0,
3511 *
ExcMessage(
"No agglomeration has been performed."));
3512 *
return {*master_cells_container.end(),
this};
3517 *
template <
int dim,
int spacedim>
3518 *
AgglomerationIterator<dim, spacedim>
3519 *
AgglomerationHandler<dim, spacedim>::end()
3521 *
Assert(n_agglomerations > 0,
3522 *
ExcMessage(
"No agglomeration has been performed."));
3523 *
return {*master_cells_container.end(),
this};
3528 *
template <
int dim,
int spacedim>
3529 *
AgglomerationIterator<dim, spacedim>
3530 *
AgglomerationHandler<dim, spacedim>::last()
3532 *
Assert(n_agglomerations > 0,
3533 *
ExcMessage(
"No agglomeration has been performed."));
3534 *
return {master_cells_container.back(),
this};
3539 *
template <
int dim,
int spacedim>
3541 *
typename AgglomerationHandler<dim, spacedim>::agglomeration_iterator>
3542 *
AgglomerationHandler<dim, spacedim>::polytope_iterators() const
3545 *
typename AgglomerationHandler<dim, spacedim>::agglomeration_iterator>(
3549 *
template <
int dim,
int spacedim>
3550 *
template <
typename RtreeType>
3552 *
AgglomerationHandler<dim, spacedim>::connect_hierarchy(
3553 *
const CellsAgglomerator<dim, RtreeType> &agglomerator)
3555 *
parent_child_info = agglomerator.parent_node_to_children_nodes;
3556 *
present_extraction_level = agglomerator.extraction_level;
3559 *
template <
int dim,
int spacedim>
3561 *
AgglomerationHandler<dim, spacedim>::get_fe_collection() const
3563 *
return *hp_fe_collection;
3566 *
template <
int dim,
int spacedim>
3568 *
AgglomerationHandler<dim, spacedim>::used_fe_collection() const
3570 *
return is_hp_collection;
3578<a name=
"ann-include/agglomeration_iterator.h"></a>
3579<h1>Annotated version of include/agglomeration_iterator.h</h1>
3596 *
#ifndef agglomeration_iterator_h
3597 *
#define agglomeration_iterator_h
3600 *
#include <agglomeration_accessor.h>
3608 *
template <
int dim,
int spacedim = dim>
3609 *
class AgglomerationIterator
3612 *
using AgglomerationContainer =
3613 *
typename AgglomerationAccessor<dim, spacedim>::AgglomerationContainer;
3619 *
AgglomerationIterator();
3625 *
AgglomerationIterator(
3627 *
const AgglomerationHandler<dim, spacedim> *handler);
3632 *
AgglomerationIterator(
3636 *
const AgglomerationHandler<dim, spacedim> *handler);
3642 *
const AgglomerationAccessor<dim, spacedim> &
3648 *
AgglomerationAccessor<dim, spacedim> &
3657 *
const AgglomerationAccessor<dim, spacedim> *
3658 *
operator->()
const;
3663 *
AgglomerationAccessor<dim, spacedim> *
3670 *
operator==(
const AgglomerationIterator<dim, spacedim> &)
const;
3676 *
operator!=(
const AgglomerationIterator<dim, spacedim> &)
const;
3683 *
AgglomerationIterator &
3691 *
AgglomerationIterator
3699 *
AgglomerationIterator &
3707 *
AgglomerationIterator
3720 *
master_cell()
const;
3727 *
using iterator_category = std::bidirectional_iterator_tag;
3728 *
using value_type = AgglomerationAccessor<dim, spacedim>;
3730 *
using pointer = AgglomerationAccessor<dim, spacedim> *;
3731 *
using reference = AgglomerationAccessor<dim, spacedim> &;
3737 *
AgglomerationAccessor<dim, spacedim> accessor;
3744 * ------------------------------
inline functions -------------------------
3750 *
template <
int dim,
int spacedim>
3751 *
inline AgglomerationIterator<dim, spacedim>::AgglomerationIterator()
3757 *
template <
int dim,
int spacedim>
3758 *
inline AgglomerationIterator<dim, spacedim>::AgglomerationIterator(
3761 *
const AgglomerationHandler<dim, spacedim> *handler)
3762 *
: accessor(master_cell, handler)
3765 *
template <
int dim,
int spacedim>
3766 *
inline AgglomerationIterator<dim, spacedim>::AgglomerationIterator(
3770 *
const AgglomerationHandler<dim, spacedim> *handler)
3771 *
: accessor(master_cell, cell_id, handler)
3776 *
template <
int dim,
int spacedim>
3777 *
inline AgglomerationAccessor<dim, spacedim> &
3778 *
AgglomerationIterator<dim, spacedim>::operator*()
3785 *
template <
int dim,
int spacedim>
3786 *
inline AgglomerationAccessor<dim, spacedim> *
3787 *
AgglomerationIterator<dim, spacedim>::operator->()
3794 *
template <
int dim,
int spacedim>
3795 *
inline const AgglomerationAccessor<dim, spacedim> &
3796 *
AgglomerationIterator<dim, spacedim>::operator*() const
3803 *
template <
int dim,
int spacedim>
3804 *
inline const AgglomerationAccessor<dim, spacedim> *
3805 *
AgglomerationIterator<dim, spacedim>::operator->() const
3812 *
template <
int dim,
int spacedim>
3814 *
AgglomerationIterator<dim, spacedim>::operator!=(
3815 *
const AgglomerationIterator<dim, spacedim> &other)
const
3817 *
return accessor != other.accessor;
3822 *
template <
int dim,
int spacedim>
3824 *
AgglomerationIterator<dim, spacedim>::operator==(
3825 *
const AgglomerationIterator<dim, spacedim> &other)
const
3827 *
return accessor == other.accessor;
3832 *
template <
int dim,
int spacedim>
3833 *
inline AgglomerationIterator<dim, spacedim> &
3834 *
AgglomerationIterator<dim, spacedim>::operator++()
3842 *
template <
int dim,
int spacedim>
3843 *
inline AgglomerationIterator<dim, spacedim>
3844 *
AgglomerationIterator<dim, spacedim>::operator++(
int)
3846 *
AgglomerationIterator tmp(*
this);
3854 *
template <
int dim,
int spacedim>
3855 *
inline AgglomerationIterator<dim, spacedim> &
3856 *
AgglomerationIterator<dim, spacedim>::operator--()
3864 *
template <
int dim,
int spacedim>
3865 *
inline AgglomerationIterator<dim, spacedim>
3866 *
AgglomerationIterator<dim, spacedim>::operator--(
int)
3868 *
AgglomerationIterator tmp(*
this);
3876 *
template <
int dim,
int spacedim>
3878 *
AgglomerationIterator<dim, spacedim>::state() const
3880 *
return accessor.master_cell.state();
3885 *
template <
int dim,
int spacedim>
3887 *
AgglomerationIterator<dim, spacedim>::master_cell() const
3889 *
return accessor.master_cell;
3898<a name=
"ann-include/agglomerator.h"></a>
3899<h1>Annotated version of include/agglomerator.h</h1>
3916 *
#ifndef agglomerator_h
3917 *
#define agglomerator_h
3920 *
#include <deal.II/base/config.h>
3922 *
#include <deal.II/base/bounding_box.h>
3924 *
#include <
boost/geometry/algorithms/distance.hpp>
3926 *
#include <
boost/geometry/strategies/strategies.hpp>
3928 *
template <
int dim,
int spacedim>
3929 *
class AgglomerationHandler;
3935 *
template <
typename Value,
3937 *
typename Translator,
3939 *
typename Allocators>
3940 *
struct Rtree_visitor
3941 *
:
public boost::geometry::index::detail::rtree::visitor<
3943 *
typename Options::parameters_type,
3946 *
typename Options::node_tag,
3949 *
inline Rtree_visitor(
3950 *
const Translator &translator,
3951 *
const unsigned int target_level,
3953 *
boost::geometry::dimension<Box>::value>::active_cell_iterator>>
3955 *
std::vector<types::global_cell_index> &n_nodes_per_level,
3956 *
std::map<std::pair<types::global_cell_index, types::global_cell_index>,
3957 *
std::vector<types::global_cell_index>> &parent_to_children);
3962 *
using InternalNode =
3963 *
typename boost::geometry::index::detail::rtree::internal_node<
3965 *
typename Options::parameters_type,
3968 *
typename Options::node_tag>::type;
3973 *
using Leaf =
typename boost::geometry::index::detail::rtree::leaf<
3975 *
typename Options::parameters_type,
3978 *
typename Options::node_tag>::type;
3998 *
const Translator &translator;
4009 *
size_t node_counter;
4015 *
const size_t target_level;
4023 *
boost::geometry::dimension<Box>::value>::active_cell_iterator>>
4029 *
std::vector<types::global_cell_index> &n_nodes_per_level;
4035 *
std::map<std::pair<types::global_cell_index, types::global_cell_index>,
4036 *
std::vector<types::global_cell_index>>
4037 *
&parent_node_to_children_nodes;
4042 *
template <
typename Value,
4044 *
typename Translator,
4046 *
typename Allocators>
4047 *
Rtree_visitor<Value, Options, Translator, Box, Allocators>::Rtree_visitor(
4048 *
const Translator &translator,
4049 *
const unsigned int target_level,
4051 *
boost::geometry::dimension<Box>::value>::active_cell_iterator>>
4053 *
std::vector<types::global_cell_index> &n_nodes_per_level_,
4054 *
std::map<std::pair<types::global_cell_index, types::global_cell_index>,
4055 *
std::vector<types::global_cell_index>> &parent_to_children)
4056 *
: translator(translator)
4059 *
, target_level(target_level)
4060 *
, agglomerates(agglomerates_)
4061 *
, n_nodes_per_level(n_nodes_per_level_)
4062 *
, parent_node_to_children_nodes(parent_to_children)
4067 *
template <
typename Value,
4069 *
typename Translator,
4071 *
typename Allocators>
4073 *
Rtree_visitor<Value, Options, Translator, Box, Allocators>::operator()(
4074 *
const Rtree_visitor::InternalNode &node)
4076 *
using elements_type =
4077 *
typename boost::geometry::index::detail::rtree::elements_type<
4078 *
InternalNode>::type;
4084 *
const elements_type &elements =
4085 *
boost::geometry::index::detail::rtree::elements(node);
4087 *
if (
level < target_level)
4089 *
size_t level_backup =
level;
4092 *
for (
typename elements_type::const_iterator it = elements.begin();
4093 *
it != elements.end();
4096 *
boost::geometry::index::detail::rtree::apply_visitor(*
this,
4100 *
level = level_backup;
4102 *
else if (
level == target_level)
4104 *
const auto offset = agglomerates.size();
4105 *
agglomerates.resize(offset + 1);
4106 *
size_t level_backup =
level;
4109 *
for (
const auto &entry : elements)
4111 *
boost::geometry::
index::detail::rtree::apply_visitor(
4116 * Done with node number
'node_counter' on
level target_level.
4123 *
n_nodes_per_level[target_level]++;
4125 *
level = level_backup;
4127 *
else if (
level > target_level)
4135 * Keep visiting until you go to the leafs.
4138 *
size_t level_backup =
level;
4144 * looping through entries of node
4147 *
for (
const auto &entry : elements)
4149 *
boost::geometry::
index::detail::rtree::apply_visitor(
4154 * done with node on
level l > target_level (not just
4158 * n_nodes_per_level[level_backup]++;
4159 * const types::global_cell_index node_idx =
4160 * n_nodes_per_level[level_backup] - 1; // so to start from 0
4162 * parent_node_to_children_nodes[{n_nodes_per_level[level_backup - 1],
4163 * level_backup - 1}]
4164 * .push_back(node_idx);
4166 * level = level_backup;
4172 * template <typename Value,
4174 * typename Translator,
4176 * typename Allocators>
4178 * Rtree_visitor<Value, Options, Translator, Box, Allocators>::operator()(
4179 * const Rtree_visitor::Leaf &leaf)
4181 * using elements_type =
4182 * typename boost::geometry::index::detail::rtree::elements_type<
4183 * Leaf>::type; // pairs of bounding box and pointer to child node
4184 * const elements_type &elements =
4185 * boost::geometry::index::detail::rtree::elements(leaf);
4187 * if (level == target_level)
4191 * If I want to extract from leaf node, i.e. the target_level is the
4192 * last one where leafs are grouped together.
4195 * const auto offset = agglomerates.size();
4196 * agglomerates.resize(offset + 1);
4198 * for (const auto &it : elements)
4199 * agglomerates[node_counter].push_back(it.second);
4202 * n_nodes_per_level[target_level]++;
4206 * for (const auto &it : elements)
4207 * agglomerates[node_counter].push_back(it.second);
4210 * if (level == target_level + 1)
4212 * const unsigned int node_idx = n_nodes_per_level[level];
4214 * parent_node_to_children_nodes[{n_nodes_per_level[level - 1],
4216 * .push_back(node_idx);
4217 * n_nodes_per_level[level]++;
4221 * } // namespace internal
4226 * * Helper class which handles agglomeration based on the R-tree data
4227 * * structure. Notice that the R-tree type is assumed to be an R-star-tree.
4229 * template <int dim, typename RtreeType>
4230 * class CellsAgglomerator
4233 * template <int, int>
4234 * friend class ::AgglomerationHandler;
4237 * * Constructor. It takes a given rtree and an integer representing the
4238 * * index of the level to be extracted.
4240 * CellsAgglomerator(const RtreeType &rtree,
4241 * const unsigned int extraction_level);
4244 * * Extract agglomerates based on the current tree and the extraction level.
4245 * * This function returns a reference to
4247 * const std::vector<
4248 * std::vector<typename Triangulation<dim>::active_cell_iterator>> &
4249 * extract_agglomerates();
4252 * * Get total number of levels.
4254 * inline unsigned int
4255 * get_n_levels() const;
4258 * * Return the number of nodes present in level @p level.
4260 * inline types::global_cell_index
4261 * get_n_nodes_per_level(const unsigned int level) const;
4264 * * This function returns a map which associates to each node on level
4265 * * @p extraction_level a list of children.
4267 * inline const std::map<
4268 * std::pair<types::global_cell_index, types::global_cell_index>,
4269 * std::vector<types::global_cell_index>> &
4270 * get_hierarchy() const;
4274 * * Raw pointer to the actual R-tree.
4279 * * Extraction level.
4281 * const unsigned int extraction_level;
4284 * * Store agglomerates obtained after recursive extraction on nodes of
4285 * * level @p extraction_level.
4287 * std::vector<std::vector<typename Triangulation<dim>::active_cell_iterator>>
4288 * agglomerates_on_level;
4291 * * Vector storing the number of nodes (and, ultimately, agglomerates) for
4294 * std::vector<types::global_cell_index> n_nodes_per_level;
4297 * * Map which maps a node parent @n on level @p l to a vector of integers
4298 * * which stores the index of children.
4300 * std::map<std::pair<types::global_cell_index, types::global_cell_index>,
4301 * std::vector<types::global_cell_index>>
4302 * parent_node_to_children_nodes;
4307 * template <int dim, typename RtreeType>
4308 * CellsAgglomerator<dim, RtreeType>::CellsAgglomerator(
4309 * const RtreeType &tree,
4310 * const unsigned int extraction_level_)
4311 * : extraction_level(extraction_level_)
4313 * rtree = const_cast<RtreeType *>(&tree);
4314 * Assert(n_levels(*rtree), ExcMessage("At least two levels are needed.
"));
4319 * template <int dim, typename RtreeType>
4320 * const std::vector<
4321 * std::vector<typename Triangulation<dim>::active_cell_iterator>> &
4322 * CellsAgglomerator<dim, RtreeType>::extract_agglomerates()
4324 * AssertThrow(extraction_level <= n_levels(*rtree),
4325 * ExcInternalError("You are trying to extract
level " +
4326 * std::to_string(extraction_level) +
4327 * " of the tree, but it only has a total of
" +
4328 * std::to_string(n_levels(*rtree)) +
4331 * boost::geometry::index::detail::rtree::utilities::view<RtreeType>;
4332 * RtreeView rtv(*rtree);
4334 * n_nodes_per_level.resize(rtv.depth() +
4335 * 1); // store how many nodes we have for each level.
4337 * if (rtv.depth() == 0)
4341 * The below algorithm does not work for `rtv.depth()==0`, which might
4342 * happen if the number entries in the tree is too small.
4345 * agglomerates_on_level.resize(1);
4346 * agglomerates_on_level[0].resize(1);
4350 * const unsigned int target_level =
4351 * std::min<unsigned int>(extraction_level, rtv.depth());
4353 * internal::Rtree_visitor<typename RtreeView::value_type,
4354 * typename RtreeView::options_type,
4355 * typename RtreeView::translator_type,
4356 * typename RtreeView::box_type,
4357 * typename RtreeView::allocators_type>
4358 * extractor_visitor(rtv.translator(),
4360 * agglomerates_on_level,
4361 * n_nodes_per_level,
4362 * parent_node_to_children_nodes);
4365 * rtv.apply_visitor(extractor_visitor);
4367 * return agglomerates_on_level;
4374 * ------------------------------ inline functions -------------------------
4383 * template <int dim, typename RtreeType>
4384 * inline unsigned int
4385 * CellsAgglomerator<dim, RtreeType>::get_n_levels() const
4387 * return n_levels(*rtree);
4392 * template <int dim, typename RtreeType>
4393 * inline types::global_cell_index
4394 * CellsAgglomerator<dim, RtreeType>::get_n_nodes_per_level(
4395 * const unsigned int level) const
4397 * return n_nodes_per_level[level];
4402 * template <int dim, typename RtreeType>
4403 * inline const std::map<
4404 * std::pair<types::global_cell_index, types::global_cell_index>,
4405 * std::vector<types::global_cell_index>> &
4406 * CellsAgglomerator<dim, RtreeType>::get_hierarchy() const
4408 * Assert(parent_node_to_children_nodes.size(),
4410 * "The hierarchy has not been computed. Did you forget to call
"
4411 * " extract_agglomerates()
first?
"));
4412 * return parent_node_to_children_nodes;
4414 * } // namespace dealii
4419<a name="ann-include/mapping_box.h
"></a>
4420<h1>Annotated version of include/mapping_box.h</h1>
4426 * /* -----------------------------------------------------------------------------
4428 * * SPDX-License-Identifier: LGPL-2.1-or-later
4429 * * Copyright (C) 2025 by Marco Feder, Pasquale Claudio Africa, Xinping Gui,
4432 * * This file is part of the deal.II code gallery.
4434 * * -----------------------------------------------------------------------------
4437 * #ifndef dealii_mapping_box_h
4438 * #define dealii_mapping_box_h
4441 * #include <deal.II/base/config.h>
4443 * #include <deal.II/base/bounding_box.h>
4444 * #include <deal.II/base/qprojector.h>
4446 * #include <deal.II/fe/mapping.h>
4451 * DEAL_II_NAMESPACE_OPEN
4454 * * @addtogroup mapping
4459 * * A class providing a mapping from the reference cell to cells that are
4460 * * axiparallel, i.e., that have the shape of rectangles (in 2d) or
4461 * * boxes (in 3d) with edges parallel to the coordinate directions. The
4462 * * class therefore provides functionality that is equivalent to what,
4463 * * for example, MappingQ would provide for such cells. However, knowledge
4464 * * of the shape of cells allows this class to be substantially more
4467 * * Specifically, the mapping is meant for cells for which the mapping from
4468 * * the reference to the real cell is a scaling along the coordinate
4469 * * directions: The transformation from reference coordinates \hat {\mathbf
4470 * * x} to real coordinates \mathbf x on each cell is of the form
4472 * * {\mathbf x}(\hat {\mathbf x})
4483 * * {\mathbf x}(\hat {\mathbf x})
4493 * * in 3d, where {\mathbf v}_0 is the bottom left vertex and h_x,h_y,h_z
4494 * * are the extents of the cell along the axes.
4496 * * The class is intended for efficiency, and it does not do a whole lot of
4497 * * error checking. If you apply this mapping to a cell that does not conform
4498 * * to the requirements above, you will get strange results.
4500 * template <int dim, int spacedim = dim>
4501 * class MappingBox : public Mapping<dim, spacedim>
4504 * MappingBox(const std::vector<BoundingBox<dim>> &local_boxes,
4505 * const std::map<types::global_cell_index, types::global_cell_index>
4506 * &polytope_translator);
4509 * for documentation, see the Mapping base class
4512 * virtual std::unique_ptr<Mapping<dim, spacedim>>
4513 * clone() const override;
4516 * * Return @p true because MappingBox preserves vertex
4520 * preserves_vertex_locations() const override;
4523 * is_compatible_with(
4524 * #if DEAL_II_VERSION_GTE(9, 8, 0)
4525 * const ReferenceCell<dim> &reference_cell
4527 * const ReferenceCell &reference_cell
4532 * * @name Mapping points between reference and real cells
4538 * for documentation, see the Mapping base class
4541 * virtual Point<spacedim>
4542 * transform_unit_to_real_cell(
4543 * const typename Triangulation<dim, spacedim>::cell_iterator &cell,
4544 * const Point<dim> &p) const override;
4548 * for documentation, see the Mapping base class
4551 * virtual Point<dim>
4552 * transform_real_to_unit_cell(
4553 * const typename Triangulation<dim, spacedim>::cell_iterator &cell,
4554 * const Point<spacedim> &p) const override;
4558 * for documentation, see the Mapping base class
4562 * transform_points_real_to_unit_cell(
4563 * const typename Triangulation<dim, spacedim>::cell_iterator &cell,
4564 * const ArrayView<const Point<spacedim>> &real_points,
4565 * const ArrayView<Point<dim>> &unit_points) const override;
4572 * * @name Functions to transform tensors from reference to real coordinates
4578 * for documentation, see the Mapping base class
4582 * transform(const ArrayView<const Tensor<1, dim>> &input,
4583 * const MappingKind kind,
4584 * const typename Mapping<dim, spacedim>::InternalDataBase &internal,
4585 * const ArrayView<Tensor<1, spacedim>> &output) const override;
4589 * for documentation, see the Mapping base class
4593 * transform(const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
4594 * const MappingKind kind,
4595 * const typename Mapping<dim, spacedim>::InternalDataBase &internal,
4596 * const ArrayView<Tensor<2, spacedim>> &output) const override;
4600 * for documentation, see the Mapping base class
4604 * transform(const ArrayView<const Tensor<2, dim>> &input,
4605 * const MappingKind kind,
4606 * const typename Mapping<dim, spacedim>::InternalDataBase &internal,
4607 * const ArrayView<Tensor<2, spacedim>> &output) const override;
4611 * for documentation, see the Mapping base class
4615 * transform(const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
4616 * const MappingKind kind,
4617 * const typename Mapping<dim, spacedim>::InternalDataBase &internal,
4618 * const ArrayView<Tensor<3, spacedim>> &output) const override;
4622 * for documentation, see the Mapping base class
4626 * transform(const ArrayView<const Tensor<3, dim>> &input,
4627 * const MappingKind kind,
4628 * const typename Mapping<dim, spacedim>::InternalDataBase &internal,
4629 * const ArrayView<Tensor<3, spacedim>> &output) const override;
4636 * * @name Interface with FEValues
4641 * * Storage for internal data of the mapping. See Mapping::InternalDataBase
4642 * * for an extensive description.
4644 * * This includes data that is computed once when the object is created (in
4645 * * get_data()) as well as data the class wants to store from between the
4646 * * call to fill_fe_values(), fill_fe_face_values(), or
4647 * * fill_fe_subface_values() until possible later calls from the finite
4648 * * element to functions such as transform(). The latter class of member
4649 * * variables are marked as 'mutable'.
4651 * class InternalData : public Mapping<dim, spacedim>::InternalDataBase
4655 * * Default constructor.
4657 * InternalData() = default;
4660 * * Constructor that initializes the object with a quadrature.
4662 * InternalData(const Quadrature<dim> &quadrature);
4666 * Documentation see Mapping::InternalDataBase.
4670 * reinit(const UpdateFlags update_flags,
4671 * const Quadrature<dim> &quadrature) override;
4674 * * Return an estimate (in bytes) for the memory consumption of this object.
4676 * virtual std::size_t
4677 * memory_consumption() const override;
4680 * * Extents of the last cell we have seen in the coordinate directions,
4681 * * i.e., <i>h<sub>x</sub></i>, <i>h<sub>y</sub></i>, <i>h<sub>z</sub></i>.
4683 * mutable Tensor<1, dim> cell_extents;
4686 * * Traslation term in F(\hat{x})=J\hat{x} + c.
4688 * mutable Tensor<1, dim> traslation;
4691 * * Reciprocal of the extents of the last cell we have seen in the
4692 * * coordinate directions, i.e., <i>h<sub>x</sub></i>,
4693 * * <i>h<sub>y</sub></i>, <i>h<sub>z</sub></i>.
4695 * mutable Tensor<1, dim> inverse_cell_extents;
4698 * * The volume element
4700 * mutable double volume_element;
4703 * * Location of quadrature points of faces or subfaces in 3d with all
4704 * * possible orientations. Can be accessed with the correct offset provided
4705 * * via QProjector::DataSetDescriptor. Not needed/used for cells.
4707 * std::vector<Point<dim>> quadrature_points;
4713 * documentation can be found in Mapping::requires_update_flags()
4716 * virtual UpdateFlags
4717 * requires_update_flags(const UpdateFlags update_flags) const override;
4721 * documentation can be found in Mapping::get_data()
4724 * virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
4725 * get_data(const UpdateFlags, const Quadrature<dim> &quadrature) const override;
4727 * using Mapping<dim, spacedim>::get_face_data;
4731 * documentation can be found in Mapping::get_subface_data()
4734 * virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
4735 * get_subface_data(const UpdateFlags flags,
4736 * const Quadrature<dim - 1> &quadrature) const override;
4740 * documentation can be found in Mapping::fill_fe_values()
4743 * virtual CellSimilarity::Similarity
4745 * const typename Triangulation<dim, spacedim>::cell_iterator &cell,
4746 * const CellSimilarity::Similarity cell_similarity,
4747 * const Quadrature<dim> &quadrature,
4748 * const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
4749 * internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
4750 * &output_data) const override;
4752 * using Mapping<dim, spacedim>::fill_fe_face_values;
4756 * documentation can be found in Mapping::fill_fe_subface_values()
4760 * fill_fe_subface_values(
4761 * const typename Triangulation<dim, spacedim>::cell_iterator &cell,
4762 * const unsigned int face_no,
4763 * const unsigned int subface_no,
4764 * const Quadrature<dim - 1> &quadrature,
4765 * const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
4766 * internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
4767 * &output_data) const override;
4771 * documentation can be found in Mapping::fill_fe_immersed_surface_values()
4775 * fill_fe_immersed_surface_values(
4776 * const typename Triangulation<dim, spacedim>::cell_iterator &cell,
4777 * const NonMatching::ImmersedSurfaceQuadrature<dim> &quadrature,
4778 * const typename Mapping<dim, spacedim>::InternalDataBase &internal_data,
4779 * internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
4780 * &output_data) const override;
4787 * * Update the cell_extents field of the incoming InternalData object with the
4788 * * size of the incoming cell.
4791 * update_cell_extents(
4792 * const typename Triangulation<dim, spacedim>::cell_iterator &cell,
4793 * const CellSimilarity::Similarity cell_similarity,
4794 * const InternalData &data) const;
4797 * * Compute the quadrature points if the UpdateFlags of the incoming
4798 * * InternalData object say that they should be updated.
4800 * * Called from fill_fe_values.
4803 * maybe_update_cell_quadrature_points(
4804 * const typename Triangulation<dim, spacedim>::cell_iterator &cell,
4805 * const InternalData &data,
4806 * const ArrayView<const Point<dim>> &unit_quadrature_points,
4807 * std::vector<Point<dim>> &quadrature_points) const;
4810 * * Compute the normal vectors if the UpdateFlags of the incoming InternalData
4811 * * object say that they should be updated.
4814 * maybe_update_normal_vectors(
4815 * const unsigned int face_no,
4816 * const InternalData &data,
4817 * std::vector<Tensor<1, dim>> &normal_vectors) const;
4820 * * Since the Jacobian is constant for this mapping all derivatives of the
4821 * * Jacobian are identically zero. Fill these quantities with zeros if the
4822 * * corresponding update flags say that they should be updated.
4825 * maybe_update_jacobian_derivatives(
4826 * const InternalData &data,
4827 * const CellSimilarity::Similarity cell_similarity,
4828 * internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
4829 * &output_data) const;
4833 * * Compute the volume elements if the UpdateFlags of the incoming
4834 * * InternalData object say that they should be updated.
4837 * maybe_update_volume_elements(const InternalData &data) const;
4840 * * Compute the Jacobians if the UpdateFlags of the incoming
4841 * * InternalData object say that they should be updated.
4844 * maybe_update_jacobians(
4845 * const InternalData &data,
4846 * const CellSimilarity::Similarity cell_similarity,
4847 * internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
4848 * &output_data) const;
4851 * * Compute the inverse Jacobians if the UpdateFlags of the incoming
4852 * * InternalData object say that they should be updated.
4855 * maybe_update_inverse_jacobians(
4856 * const InternalData &data,
4857 * const CellSimilarity::Similarity cell_similarity,
4858 * internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
4859 * &output_data) const;
4862 * * Vector of (local) bounding boxes
4864 * std::vector<BoundingBox<dim>> boxes;
4867 * * Map from global cell index to bounding box index
4869 * std::map<types::global_cell_index, types::global_cell_index>
4870 * polytope_translator;
4875 * DEAL_II_NAMESPACE_CLOSE
4881<a name="ann-include/poly_utils.h
"></a>
4882<h1>Annotated version of include/poly_utils.h</h1>
4888 * /* -----------------------------------------------------------------------------
4890 * * SPDX-License-Identifier: LGPL-2.1-or-later
4891 * * Copyright (C) 2025 by Marco Feder, Pasquale Claudio Africa, Xinping Gui,
4894 * * This file is part of the deal.II code gallery.
4896 * * -----------------------------------------------------------------------------
4899 * #ifndef poly_utils_h
4900 * #define poly_utils_h
4902 * #include <deal.II/base/config.h>
4904 * #include <deal.II/base/conditional_ostream.h>
4905 * #include <deal.II/base/point.h>
4906 * #include <deal.II/base/quadrature.h>
4907 * #include <deal.II/base/std_cxx20/iota_view.h>
4909 * #include <deal.II/boost_adaptors/bounding_box.h>
4910 * #include <deal.II/boost_adaptors/point.h>
4911 * #include <deal.II/boost_adaptors/segment.h>
4913 * #include <deal.II/distributed/tria.h>
4915 * #include <deal.II/dofs/dof_handler.h>
4917 * #include <deal.II/fe/fe_dgq.h>
4918 * #include <deal.II/fe/fe_values.h>
4920 * #include <deal.II/grid/grid_tools.h>
4922 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
4923 * #include <deal.II/lac/sparse_matrix.h>
4924 * #include <deal.II/lac/sparsity_pattern.h>
4925 * #include <deal.II/lac/sparsity_tools.h>
4926 * #include <deal.II/lac/trilinos_sparse_matrix.h>
4928 * #include <deal.II/numerics/vector_tools_common.h>
4930 * #include <boost/geometry/algorithms/distance.hpp>
4931 * #include <boost/geometry/index/detail/rtree/utilities/print.hpp>
4932 * #include <boost/geometry/index/rtree.hpp>
4933 * #include <boost/geometry/strategies/strategies.hpp>
4937 * namespace ::PolyUtils::internal
4940 * * Helper function to compute the position of index @p index in vector @p v.
4942 * inline types::global_cell_index
4943 * get_index(const std::vector<types::global_cell_index> &v,
4944 * const types::global_cell_index index)
4946 * return std::distance(v.begin(), std::find(v.begin(), v.end(), index));
4950 * * Compute the connectivity graph for locally owned regions of a distributed
4953 * template <int dim, int spacedim>
4955 * get_face_connectivity_of_cells(
4956 * const parallel::fullydistributed::Triangulation<dim, spacedim>
4958 * DynamicSparsityPattern &cell_connectivity,
4959 * const std::vector<types::global_cell_index> locally_owned_cells)
4961 * cell_connectivity.reinit(triangulation.n_locally_owned_active_cells(),
4962 * triangulation.n_locally_owned_active_cells());
4966 * loop over all cells and their neighbors to build the sparsity
4967 * pattern. note that it's a bit hard to enter all the connections when
4968 * a neighbor has children since we would need to find out which of its
4969 * children is adjacent to the current cell. this problem can be omitted
4970 * if we only do something if the neighbor has no children -- in that
4971 * case it is either on the same or a coarser level than we are. in
4972 * return, we have to add entries in both directions for both cells
4975 * for (const auto &cell : triangulation.active_cell_iterators())
4977 * if (cell->is_locally_owned())
4979 * const unsigned int index = cell->active_cell_index();
4980 * cell_connectivity.add(get_index(locally_owned_cells, index),
4981 * get_index(locally_owned_cells, index));
4982 * for (auto f : cell->face_indices())
4983 * if ((cell->at_boundary(f) == false) &&
4984 * (cell->neighbor(f)->has_children() == false) &&
4985 * cell->neighbor(f)->is_locally_owned())
4987 * const unsigned int other_index =
4988 * cell->neighbor(f)->active_cell_index();
4990 * cell_connectivity.add(get_index(locally_owned_cells, index),
4991 * get_index(locally_owned_cells,
4993 * cell_connectivity.add(get_index(locally_owned_cells,
4995 * get_index(locally_owned_cells, index));
5000 * } // namespace ::PolyUtils::internal
5002 * namespace ::PolyUtils
5004 * template <typename Value,
5006 * typename Translator,
5008 * typename Allocators>
5009 * struct Rtree_visitor : public boost::geometry::index::detail::rtree::visitor<
5011 * typename Options::parameters_type,
5014 * typename Options::node_tag,
5017 * inline Rtree_visitor(
5018 * const Translator &translator,
5019 * unsigned int target_level,
5020 * std::vector<std::vector<typename Triangulation<
5021 * boost::geometry::dimension<Box>::value>::active_cell_iterator>> &boxes,
5022 * std::vector<std::vector<unsigned int>> &csr);
5025 * * An alias that identifies an InternalNode of the tree.
5027 * using InternalNode =
5028 * typename boost::geometry::index::detail::rtree::internal_node<
5030 * typename Options::parameters_type,
5033 * typename Options::node_tag>::type;
5036 * * An alias that identifies a Leaf of the tree.
5038 * using Leaf = typename boost::geometry::index::detail::rtree::leaf<
5040 * typename Options::parameters_type,
5043 * typename Options::node_tag>::type;
5046 * * Implements the visitor interface for InternalNode objects. If the node
5047 * * belongs to the level next to @p target_level, then fill the bounding box
5048 * * vector for that node.
5051 * operator()(const InternalNode &node);
5054 * * Implements the visitor interface for Leaf objects.
5057 * operator()(const Leaf &);
5060 * * Translator interface, required by the boost implementation of the rtree.
5062 * const Translator &translator;
5065 * * Store the level we are currently visiting.
5070 * * Index used to keep track of the number of different visited nodes during
5073 * size_t node_counter;
5075 * size_t next_level_leafs_processed;
5077 * * The level where children are living.
5078 * * Before: "we want to extract from the
RTree object.
"
5080 * const size_t target_level;
5083 * * A reference to the input vector of vector of BoundingBox objects. This
5084 * * vector v has the following property: v[i] = vector with all
5085 * * of the BoundingBox bounded by the i-th node of the Rtree.
5087 * std::vector<std::vector<typename Triangulation<
5088 * boost::geometry::dimension<Box>::value>::active_cell_iterator>>
5091 * std::vector<std::vector<unsigned int>> &row_ptr;
5094 * template <typename Value,
5096 * typename Translator,
5098 * typename Allocators>
5099 * Rtree_visitor<Value, Options, Translator, Box, Allocators>::Rtree_visitor(
5100 * const Translator &translator,
5101 * const unsigned int target_level,
5102 * std::vector<std::vector<typename Triangulation<
5103 * boost::geometry::dimension<Box>::value>::active_cell_iterator>>
5105 * std::vector<std::vector<unsigned int>> &csr)
5106 * : translator(translator)
5109 * , next_level_leafs_processed(0)
5110 * , target_level(target_level)
5111 * , agglomerates(bb_in_boxes)
5115 * template <typename Value,
5117 * typename Translator,
5119 * typename Allocators>
5121 * Rtree_visitor<Value, Options, Translator, Box, Allocators>::operator()(
5122 * const Rtree_visitor::InternalNode &node)
5124 * using elements_type =
5125 * typename boost::geometry::index::detail::rtree::elements_type<
5126 * InternalNode>::type; // pairs of bounding box and pointer to child
5132 * const elements_type &elements =
5133 * boost::geometry::index::detail::rtree::elements(node);
5135 * if (level < target_level)
5137 * size_t level_backup = level;
5140 * for (typename elements_type::const_iterator it = elements.begin();
5141 * it != elements.end();
5144 * boost::geometry::index::detail::rtree::apply_visitor(*this,
5148 * level = level_backup;
5150 * else if (level == target_level)
5154 * const unsigned int n_children = elements.size();
5157 * const auto offset = agglomerates.size();
5158 * agglomerates.resize(offset + 1);
5159 * row_ptr.resize(row_ptr.size() + 1);
5160 * next_level_leafs_processed = 0;
5161 * row_ptr.back().push_back(
5162 * next_level_leafs_processed); // convention: row_ptr[0]=0
5163 * size_t level_backup = level;
5166 * for (const auto &child : elements)
5168 * boost::geometry::index::detail::rtree::apply_visitor(*this,
5173 * Done with node number 'node_counter'
5179 * ++node_counter; // visited all children of an internal node
5181 * level = level_backup;
5183 * else if (level > target_level)
5187 * Keep visiting until you go to the leafs.
5190 * size_t level_backup = level;
5194 * for (const auto &child : elements)
5196 * boost::geometry::index::detail::rtree::apply_visitor(*this,
5199 * level = level_backup;
5200 * row_ptr[node_counter].push_back(next_level_leafs_processed);
5204 * template <typename Value,
5206 * typename Translator,
5208 * typename Allocators>
5210 * Rtree_visitor<Value, Options, Translator, Box, Allocators>::operator()(
5211 * const Rtree_visitor::Leaf &leaf)
5213 * using elements_type =
5214 * typename boost::geometry::index::detail::rtree::elements_type<
5215 * Leaf>::type; // pairs of bounding box and pointer to child node
5216 * const elements_type &elements =
5217 * boost::geometry::index::detail::rtree::elements(leaf);
5219 * for (const auto &it : elements)
5221 * agglomerates[node_counter].push_back(it.second);
5223 * next_level_leafs_processed += elements.size();
5226 * template <typename T>
5227 * inline constexpr T
5228 * constexpr_pow(T num, unsigned int pow)
5230 * return (pow >= sizeof(unsigned int) * 8) ? 0 :
5232 * num * constexpr_pow(num, pow - 1);
5235 * namespace internal
5238 * * Same as the public free function with the same name, but storing
5239 * * explicitly the interpolation matrix and performing interpolation through
5240 * * matrix-vector product.
5242 * template <int dim, int spacedim, typename VectorType>
5244 * interpolate_to_fine_grid(
5245 * const AgglomerationHandler<dim, spacedim> &agglomeration_handler,
5247 * const VectorType &src)
5249 * Assert((dim == spacedim), ExcNotImplemented());
5253 * "The destination vector must the empt upon calling
this function.
"));
5255 * using NumberType = typename VectorType::value_type;
5256 * constexpr bool is_trilinos_vector =
5257 * std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector>;
5258 * using MatrixType = std::conditional_t<is_trilinos_vector,
5259 * TrilinosWrappers::SparseMatrix,
5260 * SparseMatrix<NumberType>>;
5262 * MatrixType interpolation_matrix;
5265 * typename std::conditional_t<!is_trilinos_vector, SparsityPattern, void *>
5270 * Get some info from the handler
5273 * const DoFHandler<dim, spacedim> &agglo_dh =
5274 * agglomeration_handler.agglo_dh;
5276 * DoFHandler<dim, spacedim> *output_dh =
5277 * const_cast<DoFHandler<dim, spacedim> *>(
5278 * &agglomeration_handler.output_dh);
5279 * const FiniteElement<dim, spacedim> &fe = agglomeration_handler.get_fe();
5280 * const Mapping<dim> &mapping = agglomeration_handler.get_mapping();
5281 * const Triangulation<dim, spacedim> &tria =
5282 * agglomeration_handler.get_triangulation();
5283 * const auto &bboxes = agglomeration_handler.get_local_bboxes();
5285 * std::unique_ptr<FiniteElement<dim>> output_fe;
5286 * if (tria.all_reference_cells_are_hyper_cube())
5287 * output_fe = std::make_unique<FE_DGQ<dim>>(fe.degree);
5288 * else if (tria.all_reference_cells_are_simplex())
5289 * output_fe = std::make_unique<FE_SimplexDGP<dim>>(fe.degree);
5291 * AssertThrow(false, ExcNotImplemented());
5295 * Setup an auxiliary DoFHandler for output purposes
5298 * output_dh->reinit(tria);
5299 * output_dh->distribute_dofs(*output_fe);
5301 * const IndexSet &locally_owned_dofs = output_dh->locally_owned_dofs();
5302 * const IndexSet locally_relevant_dofs =
5303 * DoFTools::extract_locally_relevant_dofs(*output_dh);
5305 * const IndexSet &locally_owned_dofs_agglo = agglo_dh.locally_owned_dofs();
5307 * DynamicSparsityPattern dsp(output_dh->n_dofs(),
5308 * agglo_dh.n_dofs(),
5309 * locally_relevant_dofs);
5311 * std::vector<types::global_dof_index> agglo_dof_indices(fe.dofs_per_cell);
5312 * std::vector<types::global_dof_index> standard_dof_indices(
5313 * fe.dofs_per_cell);
5314 * std::vector<types::global_dof_index> output_dof_indices(
5315 * output_fe->dofs_per_cell);
5317 * Quadrature<dim> quad(output_fe->get_unit_support_points());
5318 * FEValues<dim, spacedim> output_fe_values(mapping,
5321 * update_quadrature_points);
5323 * for (const auto &cell : agglo_dh.active_cell_iterators())
5324 * if (cell->is_locally_owned())
5326 * if (agglomeration_handler.is_master_cell(cell))
5328 * auto slaves = agglomeration_handler.get_slaves_of_idx(
5329 * cell->active_cell_index());
5330 * slaves.emplace_back(cell);
5332 * cell->get_dof_indices(agglo_dof_indices);
5334 * for (const auto &slave : slaves)
5338 * addd master-slave relationship
5341 * const auto slave_output =
5342 * slave->as_dof_handler_iterator(*output_dh);
5343 * slave_output->get_dof_indices(output_dof_indices);
5344 * for (const auto row : output_dof_indices)
5345 * dsp.add_entries(row,
5346 * agglo_dof_indices.begin(),
5347 * agglo_dof_indices.end());
5352 * const auto assemble_interpolation_matrix = [&]() {
5353 * FullMatrix<NumberType> local_matrix(fe.dofs_per_cell, fe.dofs_per_cell);
5354 * std::vector<Point<dim>> reference_q_points(fe.dofs_per_cell);
5358 * Dummy AffineConstraints, only needed for loc2glb
5361 * AffineConstraints<NumberType> c;
5364 * for (const auto &cell : agglo_dh.active_cell_iterators())
5365 * if (cell->is_locally_owned())
5367 * if (agglomeration_handler.is_master_cell(cell))
5369 * auto slaves = agglomeration_handler.get_slaves_of_idx(
5370 * cell->active_cell_index());
5371 * slaves.emplace_back(cell);
5373 * cell->get_dof_indices(agglo_dof_indices);
5375 * const types::global_cell_index polytope_index =
5376 * agglomeration_handler.cell_to_polytope_index(cell);
5380 * Get the box of this agglomerate.
5383 * const BoundingBox<dim> &box = bboxes[polytope_index];
5385 * for (const auto &slave : slaves)
5389 * add master-slave relationship
5392 * const auto slave_output =
5393 * slave->as_dof_handler_iterator(*output_dh);
5395 * slave_output->get_dof_indices(output_dof_indices);
5396 * output_fe_values.reinit(slave_output);
5398 * local_matrix = 0.;
5400 * const auto &q_points =
5401 * output_fe_values.get_quadrature_points();
5402 * for (const auto i : output_fe_values.dof_indices())
5404 * const auto &p = box.real_to_unit(q_points[i]);
5405 * for (const auto j : output_fe_values.dof_indices())
5407 * local_matrix(i, j) = fe.shape_value(j, p);
5410 * c.distribute_local_to_global(local_matrix,
5411 * output_dof_indices,
5412 * agglo_dof_indices,
5413 * interpolation_matrix);
5419 * if constexpr (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix>)
5421 * const MPI_Comm &communicator = tria.get_mpi_communicator();
5422 * SparsityTools::distribute_sparsity_pattern(dsp,
5423 * locally_owned_dofs,
5425 * locally_relevant_dofs);
5427 * interpolation_matrix.reinit(locally_owned_dofs,
5428 * locally_owned_dofs_agglo,
5431 * dst.reinit(locally_owned_dofs);
5432 * assemble_interpolation_matrix();
5434 * else if constexpr (std::is_same_v<MatrixType, SparseMatrix<NumberType>>)
5436 * sp.copy_from(dsp);
5437 * interpolation_matrix.reinit(sp);
5438 * dst.reinit(output_dh->n_dofs());
5439 * assemble_interpolation_matrix();
5445 * PETSc, LA::d::v options not implemented.
5448 * (void)agglomeration_handler;
5451 * AssertThrow(false, ExcNotImplemented());
5456 * If tria is distributed
5459 * if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
5460 * &tria) != nullptr)
5461 * interpolation_matrix.compress(VectorOperation::add);
5465 * Finally, perform the interpolation.
5468 * interpolation_matrix.vmult(dst, src);
5470 * } // namespace internal
5473 * * Given a vector @p src, typically the solution stemming after the
5474 * * agglomerate problem has been solved, this function interpolates @p src
5475 * * onto the finer grid and stores the result in vector @p dst. The last
5476 * * argument @p on_the_fly does not build any interpolation matrix and allows
5477 * * computing the entries in @p dst in a matrix-free fashion.
5479 * * @note Supported parallel types are TrilinosWrappers::SparseMatrix and
5480 * * TrilinosWrappers::MPI::Vector.
5482 * template <int dim, int spacedim, typename VectorType>
5484 * interpolate_to_fine_grid(
5485 * const AgglomerationHandler<dim, spacedim> &agglomeration_handler,
5487 * const VectorType &src,
5488 * const bool on_the_fly = true)
5490 * Assert((dim == spacedim), ExcNotImplemented());
5494 * "The destination vector must the empt upon calling
this function.
"));
5496 * using NumberType = typename VectorType::value_type;
5497 * static constexpr bool is_trilinos_vector =
5498 * std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector>;
5500 * static constexpr bool is_supported_vector =
5501 * std::is_same_v<VectorType, Vector<NumberType>> || is_trilinos_vector;
5502 * static_assert(is_supported_vector);
5506 * First, check for an easy return
5509 * if (on_the_fly == false)
5511 * return internal::interpolate_to_fine_grid(agglomeration_handler,
5519 * otherwise, do not create any matrix
5522 * if (!agglomeration_handler.used_fe_collection())
5526 * Original version: handle case without hp::FECollection
5529 * const Triangulation<dim, spacedim> &tria =
5530 * agglomeration_handler.get_triangulation();
5531 * const Mapping<dim> &mapping = agglomeration_handler.get_mapping();
5532 * const FiniteElement<dim, spacedim> &original_fe =
5533 * agglomeration_handler.get_fe();
5537 * We use DGQ (on tensor-product meshes) or DGP (on simplex meshes)
5538 * nodal elements of the same degree as the ones in the
5539 * agglomeration handler to interpolate the solution onto the finer
5543 * std::unique_ptr<FiniteElement<dim>> output_fe;
5544 * if (tria.all_reference_cells_are_hyper_cube())
5545 * output_fe = std::make_unique<FE_DGQ<dim>>(original_fe.degree);
5546 * else if (tria.all_reference_cells_are_simplex())
5548 * std::make_unique<FE_SimplexDGP<dim>>(original_fe.degree);
5550 * AssertThrow(false, ExcNotImplemented());
5552 * DoFHandler<dim> &output_dh =
5553 * const_cast<DoFHandler<dim> &>(agglomeration_handler.output_dh);
5554 * output_dh.reinit(tria);
5555 * output_dh.distribute_dofs(*output_fe);
5557 * if constexpr (std::is_same_v<VectorType,
5558 * TrilinosWrappers::MPI::Vector>)
5560 * const IndexSet &locally_owned_dofs =
5561 * output_dh.locally_owned_dofs();
5562 * dst.reinit(locally_owned_dofs);
5564 * else if constexpr (std::is_same_v<VectorType, Vector<NumberType>>)
5566 * dst.reinit(output_dh.n_dofs());
5572 * PETSc, LA::d::v options not implemented.
5575 * (void)agglomeration_handler;
5578 * AssertThrow(false, ExcNotImplemented());
5581 * const unsigned int dofs_per_cell =
5582 * agglomeration_handler.n_dofs_per_cell();
5583 * const unsigned int output_dofs_per_cell =
5584 * output_fe->n_dofs_per_cell();
5585 * Quadrature<dim> quad(output_fe->get_unit_support_points());
5586 * FEValues<dim> output_fe_values(mapping,
5589 * update_quadrature_points);
5591 * std::vector<types::global_dof_index> local_dof_indices(
5593 * std::vector<types::global_dof_index> local_dof_indices_output(
5594 * output_dofs_per_cell);
5596 * const auto &bboxes = agglomeration_handler.get_local_bboxes();
5597 * for (const auto &polytope :
5598 * agglomeration_handler.polytope_iterators())
5600 * if (polytope->is_locally_owned())
5602 * polytope->get_dof_indices(local_dof_indices);
5603 * const BoundingBox<dim> &box = bboxes[polytope->index()];
5605 * const auto &deal_cells =
5606 * polytope->get_agglomerate(); // fine deal.II cells
5607 * for (const auto &cell : deal_cells)
5609 * const auto slave_output = cell->as_dof_handler_iterator(
5610 * agglomeration_handler.output_dh);
5611 * slave_output->get_dof_indices(local_dof_indices_output);
5612 * output_fe_values.reinit(slave_output);
5614 * const auto &qpoints =
5615 * output_fe_values.get_quadrature_points();
5617 * for (unsigned int j = 0; j < output_dofs_per_cell; ++j)
5619 * const auto &ref_qpoint =
5620 * box.real_to_unit(qpoints[j]);
5621 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
5622 * dst(local_dof_indices_output[j]) +=
5623 * src(local_dof_indices[i]) *
5624 * original_fe.shape_value(i, ref_qpoint);
5634 * Handle the hp::FECollection case
5637 * const Triangulation<dim, spacedim> &tria =
5638 * agglomeration_handler.get_triangulation();
5639 * const Mapping<dim> &mapping = agglomeration_handler.get_mapping();
5640 * const hp::FECollection<dim, spacedim> &original_fe_collection =
5641 * agglomeration_handler.get_fe_collection();
5645 * We use DGQ (on tensor-product meshes) or DGP (on simplex meshes)
5646 * nodal elements of the same degree as the ones in the
5647 * agglomeration handler to interpolate the solution onto the finer
5651 * hp::FECollection<dim, spacedim> output_fe_collection;
5653 * Assert(original_fe_collection[0].n_components() >= 1,
5654 * ExcMessage("Invalid FE: must have at least
one component.
"));
5655 * if (original_fe_collection[0].n_components() == 1)
5662 * for (unsigned int i = 0; i < original_fe_collection.size(); ++i)
5664 * std::unique_ptr<FiniteElement<dim>> output_fe;
5665 * if (tria.all_reference_cells_are_hyper_cube())
5666 * output_fe = std::make_unique<FE_DGQ<dim>>(
5667 * original_fe_collection[i].degree);
5668 * else if (tria.all_reference_cells_are_simplex())
5669 * output_fe = std::make_unique<FE_SimplexDGP<dim>>(
5670 * original_fe_collection[i].degree);
5672 * AssertThrow(false, ExcNotImplemented());
5673 * output_fe_collection.push_back(*output_fe);
5676 * else if (original_fe_collection[0].n_components() > 1)
5683 * for (unsigned int i = 0; i < original_fe_collection.size(); ++i)
5685 * std::vector<const FiniteElement<dim, spacedim> *>
5687 * std::vector<unsigned int> multiplicities;
5688 * for (unsigned int b = 0;
5689 * b < original_fe_collection[i].n_base_elements();
5692 * if (dynamic_cast<const FE_Nothing<dim> *>(
5693 * &original_fe_collection[i].base_element(b)))
5694 * base_elements.push_back(
5695 * new FE_Nothing<dim, spacedim>());
5698 * if (tria.all_reference_cells_are_hyper_cube())
5699 * base_elements.push_back(new FE_DGQ<dim, spacedim>(
5700 * original_fe_collection[i]
5703 * else if (tria.all_reference_cells_are_simplex())
5704 * base_elements.push_back(
5705 * new FE_SimplexDGP<dim, spacedim>(
5706 * original_fe_collection[i]
5710 * AssertThrow(false, ExcNotImplemented());
5712 * multiplicities.push_back(
5713 * original_fe_collection[i].element_multiplicity(b));
5716 * FESystem<dim, spacedim> output_fe_system(base_elements,
5718 * for (const auto *ptr : base_elements)
5720 * output_fe_collection.push_back(output_fe_system);
5724 * DoFHandler<dim> &output_dh =
5725 * const_cast<DoFHandler<dim> &>(agglomeration_handler.output_dh);
5726 * output_dh.reinit(tria);
5727 * for (const auto &polytope :
5728 * agglomeration_handler.polytope_iterators())
5730 * if (polytope->is_locally_owned())
5732 * const auto &deal_cells =
5733 * polytope->get_agglomerate(); // fine deal.II cells
5734 * const unsigned int active_fe_idx =
5735 * polytope->active_fe_index();
5737 * for (const auto &cell : deal_cells)
5739 * const typename DoFHandler<dim>::active_cell_iterator
5740 * slave_cell_dh_iterator =
5741 * cell->as_dof_handler_iterator(output_dh);
5742 * slave_cell_dh_iterator->set_active_fe_index(
5747 * output_dh.distribute_dofs(output_fe_collection);
5749 * if constexpr (std::is_same_v<VectorType,
5750 * TrilinosWrappers::MPI::Vector>)
5752 * const IndexSet &locally_owned_dofs =
5753 * output_dh.locally_owned_dofs();
5754 * dst.reinit(locally_owned_dofs);
5756 * else if constexpr (std::is_same_v<VectorType, Vector<NumberType>>)
5758 * dst.reinit(output_dh.n_dofs());
5764 * PETSc, LA::d::v options not implemented.
5767 * (void)agglomeration_handler;
5770 * AssertThrow(false, ExcNotImplemented());
5773 * const auto &bboxes = agglomeration_handler.get_local_bboxes();
5774 * for (const auto &polytope :
5775 * agglomeration_handler.polytope_iterators())
5777 * if (polytope->is_locally_owned())
5779 * const unsigned int active_fe_idx =
5780 * polytope->active_fe_index();
5781 * const unsigned int dofs_per_cell =
5782 * polytope->get_fe().dofs_per_cell;
5783 * const unsigned int output_dofs_per_cell =
5784 * output_fe_collection[active_fe_idx].n_dofs_per_cell();
5785 * Quadrature<dim> quad(output_fe_collection[active_fe_idx]
5786 * .get_unit_support_points());
5787 * FEValues<dim> output_fe_values(
5789 * output_fe_collection[active_fe_idx],
5791 * update_quadrature_points);
5792 * std::vector<types::global_dof_index> local_dof_indices(
5794 * std::vector<types::global_dof_index>
5795 * local_dof_indices_output(output_dofs_per_cell);
5797 * polytope->get_dof_indices(local_dof_indices);
5798 * const BoundingBox<dim> &box = bboxes[polytope->index()];
5800 * const auto &deal_cells =
5801 * polytope->get_agglomerate(); // fine deal.II cells
5802 * for (const auto &cell : deal_cells)
5804 * const auto slave_output = cell->as_dof_handler_iterator(
5805 * agglomeration_handler.output_dh);
5806 * slave_output->get_dof_indices(local_dof_indices_output);
5807 * output_fe_values.reinit(slave_output);
5809 * const auto &qpoints =
5810 * output_fe_values.get_quadrature_points();
5812 * for (unsigned int j = 0; j < output_dofs_per_cell; ++j)
5814 * const unsigned int component_idx_of_this_dof =
5815 * slave_output->get_fe()
5816 * .system_to_component_index(j)
5818 * const auto &ref_qpoint =
5819 * box.real_to_unit(qpoints[j]);
5820 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
5821 * dst(local_dof_indices_output[j]) +=
5822 * src(local_dof_indices[i]) *
5823 * original_fe_collection[active_fe_idx]
5824 * .shape_value_component(
5825 * i, ref_qpoint, component_idx_of_this_dof);
5835 * * Similar to VectorTools::compute_global_error(), but customized for
5836 * * polytopic elements. Aside from the solution vector and a reference
5837 * * function, this function takes in addition a vector @p norms with types
5838 * * VectorTools::NormType to be computed and later stored in the last
5839 * * argument @p global_errors.
5840 * * In case of a parallel vector, the local errors are collected over each
5841 * * processor and later a classical reduction operation is performed.
5843 * template <int dim, typename Number, typename VectorType>
5845 * compute_global_error(const AgglomerationHandler<dim> &agglomeration_handler,
5846 * const VectorType &solution,
5847 * const Function<dim, Number> &exact_solution,
5848 * const std::vector<VectorTools::NormType> &norms,
5849 * std::vector<double> &global_errors)
5851 * Assert(solution.size() > 0,
5852 * ExcNotImplemented(
5853 * "Solution vector must be non-empty upon calling this function.
"));
5854 * Assert(std::any_of(norms.cbegin(),
5856 * [](VectorTools::NormType norm_type) {
5857 * return (norm_type ==
5858 * VectorTools::NormType::H1_seminorm ||
5859 * norm_type == VectorTools::NormType::L2_norm);
5861 * ExcMessage("Norm type not supported
"));
5862 * global_errors.resize(norms.size());
5863 * std::fill(global_errors.begin(), global_errors.end(), 0.);
5867 * Vector storing errors local to the current processor.
5870 * std::vector<double> local_errors(norms.size());
5871 * std::fill(local_errors.begin(), local_errors.end(), 0.);
5875 * Get some info from the handler
5878 * const unsigned int dofs_per_cell = agglomeration_handler.n_dofs_per_cell();
5880 * const bool compute_semi_H1 =
5881 * std::any_of(norms.cbegin(),
5883 * [](VectorTools::NormType norm_type) {
5884 * return norm_type == VectorTools::NormType::H1_seminorm;
5887 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
5888 * for (const auto &polytope : agglomeration_handler.polytope_iterators())
5890 * if (polytope->is_locally_owned())
5892 * const auto &agglo_values = agglomeration_handler.reinit(polytope);
5893 * polytope->get_dof_indices(local_dof_indices);
5895 * const auto &q_points = agglo_values.get_quadrature_points();
5896 * const unsigned int n_qpoints = q_points.size();
5897 * std::vector<double> analyical_sol_at_qpoints(n_qpoints);
5898 * exact_solution.value_list(q_points, analyical_sol_at_qpoints);
5899 * std::vector<Tensor<1, dim>> grad_analyical_sol_at_qpoints(
5902 * if (compute_semi_H1)
5903 * exact_solution.gradient_list(q_points,
5904 * grad_analyical_sol_at_qpoints);
5906 * for (unsigned int q_index : agglo_values.quadrature_point_indices())
5908 * double solution_at_qpoint = 0.;
5909 * Tensor<1, dim> grad_solution_at_qpoint;
5910 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
5912 * solution_at_qpoint += solution(local_dof_indices[i]) *
5913 * agglo_values.shape_value(i, q_index);
5915 * if (compute_semi_H1)
5916 * grad_solution_at_qpoint +=
5917 * solution(local_dof_indices[i]) *
5918 * agglo_values.shape_grad(i, q_index);
5925 * local_errors[0] += std::pow((analyical_sol_at_qpoints[q_index] -
5926 * solution_at_qpoint),
5928 * agglo_values.JxW(q_index);
5935 * if (compute_semi_H1)
5936 * for (unsigned int d = 0; d < dim; ++d)
5937 * local_errors[1] +=
5938 * std::pow((grad_analyical_sol_at_qpoints[q_index][d] -
5939 * grad_solution_at_qpoint[d]),
5941 * agglo_values.JxW(q_index);
5948 * Perform reduction and take sqrt of each error
5951 * global_errors[0] = Utilities::MPI::reduce<double>(
5953 * agglomeration_handler.get_triangulation().get_mpi_communicator(),
5954 * [](const double a, const double b) { return a + b; });
5956 * global_errors[0] = std::sqrt(global_errors[0]);
5958 * if (compute_semi_H1)
5960 * global_errors[1] = Utilities::MPI::reduce<double>(
5962 * agglomeration_handler.get_triangulation().get_mpi_communicator(),
5963 * [](const double a, const double b) { return a + b; });
5964 * global_errors[1] = std::sqrt(global_errors[1]);
5969 * * Utility function that builds the multilevel hierarchy from the tree level
5970 * * @p starting_level. This function fills the vector of
5971 * * @p AgglomerationHandlers objects by distributing degrees of freedom on
5972 * * each level of the hierarchy. It returns the total number of levels in the
5975 * template <int dim>
5977 * construct_agglomerated_levels(
5978 * const Triangulation<dim> &tria,
5979 * std::vector<std::unique_ptr<AgglomerationHandler<dim>>>
5980 * &agglomeration_handlers,
5981 * const FE_DGQ<dim> &fe_dg,
5982 * const Mapping<dim> &mapping,
5983 * const unsigned int starting_tree_level)
5985 * const auto parallel_tria =
5986 * dynamic_cast<const parallel::TriangulationBase<dim> *>(&tria);
5988 * GridTools::Cache<dim> cached_tria(tria);
5989 * Assert(parallel_tria->n_active_cells() > 0, ExcInternalError());
5991 * const MPI_Comm comm = parallel_tria->get_mpi_communicator();
5992 * ConditionalOStream pcout(std::cout,
5993 * (Utilities::MPI::this_mpi_process(comm) == 0));
5997 * Start building R-tree
6000 * namespace bgi = boost::geometry::index;
6001 * static constexpr unsigned int max_elem_per_node =
6002 * constexpr_pow(2, dim); // 2^dim
6003 * std::vector<std::pair<BoundingBox<dim>,
6004 * typename Triangulation<dim>::active_cell_iterator>>
6005 * boxes(parallel_tria->n_locally_owned_active_cells());
6006 * unsigned int i = 0;
6007 * for (const auto &cell : parallel_tria->active_cell_iterators())
6008 * if (cell->is_locally_owned())
6009 * boxes[i++] = std::make_pair(mapping.get_bounding_box(cell), cell);
6011 * auto tree = pack_rtree<bgi::rstar<max_elem_per_node>>(boxes);
6012 * Assert(n_levels(tree) >= 2, ExcMessage("At least two levels are needed.
"));
6013 * pcout << "Total number of available levels:
" << n_levels(tree)
6016 * pcout << "Starting
level:
" << starting_tree_level << std::endl;
6017 * const unsigned int total_tree_levels =
6018 * n_levels(tree) - starting_tree_level + 1;
6022 * Resize the agglomeration handlers to the right size
6028 * agglomeration_handlers.resize(total_tree_levels);
6031 * Loop through the available levels and set AgglomerationHandlers up.
6034 * for (unsigned int extraction_level = starting_tree_level;
6035 * extraction_level <= n_levels(tree);
6036 * ++extraction_level)
6038 * agglomeration_handlers[extraction_level - starting_tree_level] =
6039 * std::make_unique<AgglomerationHandler<dim>>(cached_tria);
6040 * CellsAgglomerator<dim, decltype(tree)> agglomerator{tree,
6041 * extraction_level};
6042 * const auto agglomerates = agglomerator.extract_agglomerates();
6043 * agglomeration_handlers[extraction_level - starting_tree_level]
6044 * ->connect_hierarchy(agglomerator);
6048 * Flag elements for agglomeration
6051 * unsigned int agglo_index = 0;
6052 * for (unsigned int i = 0; i < agglomerates.size(); ++i)
6054 * const auto &agglo = agglomerates[i]; // i-th agglomerate
6055 * for (const auto &el : agglo)
6057 * el->set_material_id(agglo_index);
6062 * const unsigned int n_local_agglomerates = agglo_index;
6063 * unsigned int total_agglomerates =
6064 * Utilities::MPI::sum(n_local_agglomerates, comm);
6065 * pcout << "Total agglomerates per (tree)
level:
" << extraction_level
6066 * << ":
" << total_agglomerates << std::endl;
6070 * Now, perform agglomeration within each locally owned partition
6074 * std::vector<typename Triangulation<dim>::active_cell_iterator>>
6075 * cells_per_subdomain(n_local_agglomerates);
6076 * for (const auto &cell : parallel_tria->active_cell_iterators())
6077 * if (cell->is_locally_owned())
6078 * cells_per_subdomain[cell->material_id()].push_back(cell);
6082 * For every subdomain, agglomerate elements together
6085 * for (std::size_t i = 0; i < cells_per_subdomain.size(); ++i)
6086 * agglomeration_handlers[extraction_level - starting_tree_level]
6087 * ->define_agglomerate(cells_per_subdomain[i]);
6089 * agglomeration_handlers[extraction_level - starting_tree_level]
6090 * ->initialize_fe_values(QGauss<dim>(fe_dg.degree + 1),
6091 * update_values | update_gradients |
6092 * update_JxW_values | update_quadrature_points,
6093 * QGauss<dim - 1>(fe_dg.degree + 1),
6094 * update_JxW_values);
6095 * agglomeration_handlers[extraction_level - starting_tree_level]
6096 * ->distribute_agglomerated_dofs(fe_dg);
6099 * return total_tree_levels;
6103 * * Utility to compute jump terms when the interface is locally owned, i.e.
6104 * * both elements are locally owned.
6106 * template <int dim>
6108 * assemble_local_jumps_and_averages(FullMatrix<double> &M11,
6109 * FullMatrix<double> &M12,
6110 * FullMatrix<double> &M21,
6111 * FullMatrix<double> &M22,
6112 * const FEValuesBase<dim> &fe_faces0,
6113 * const FEValuesBase<dim> &fe_faces1,
6114 * const double penalty_constant,
6117 * const std::vector<Tensor<1, dim>> &normals = fe_faces0.get_normal_vectors();
6118 * const unsigned int dofs_per_cell =
6119 * M11.m(); // size of local matrices equals the #DoFs
6120 * for (unsigned int q_index : fe_faces0.quadrature_point_indices())
6122 * const Tensor<1, dim> &normal = normals[q_index];
6123 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
6125 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
6127 * M11(i, j) += (-0.5 * fe_faces0.shape_grad(i, q_index) * normal *
6128 * fe_faces0.shape_value(j, q_index) -
6129 * 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
6130 * fe_faces0.shape_value(i, q_index) +
6131 * (penalty_constant / h_f) *
6132 * fe_faces0.shape_value(i, q_index) *
6133 * fe_faces0.shape_value(j, q_index)) *
6134 * fe_faces0.JxW(q_index);
6135 * M12(i, j) += (0.5 * fe_faces0.shape_grad(i, q_index) * normal *
6136 * fe_faces1.shape_value(j, q_index) -
6137 * 0.5 * fe_faces1.shape_grad(j, q_index) * normal *
6138 * fe_faces0.shape_value(i, q_index) -
6139 * (penalty_constant / h_f) *
6140 * fe_faces0.shape_value(i, q_index) *
6141 * fe_faces1.shape_value(j, q_index)) *
6142 * fe_faces1.JxW(q_index);
6143 * M21(i, j) += (-0.5 * fe_faces1.shape_grad(i, q_index) * normal *
6144 * fe_faces0.shape_value(j, q_index) +
6145 * 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
6146 * fe_faces1.shape_value(i, q_index) -
6147 * (penalty_constant / h_f) *
6148 * fe_faces1.shape_value(i, q_index) *
6149 * fe_faces0.shape_value(j, q_index)) *
6150 * fe_faces1.JxW(q_index);
6151 * M22(i, j) += (0.5 * fe_faces1.shape_grad(i, q_index) * normal *
6152 * fe_faces1.shape_value(j, q_index) +
6153 * 0.5 * fe_faces1.shape_grad(j, q_index) * normal *
6154 * fe_faces1.shape_value(i, q_index) +
6155 * (penalty_constant / h_f) *
6156 * fe_faces1.shape_value(i, q_index) *
6157 * fe_faces1.shape_value(j, q_index)) *
6158 * fe_faces1.JxW(q_index);
6164 * * Same as above, but for a ghosted neighbor.
6166 * template <int dim>
6168 * assemble_local_jumps_and_averages_ghost(
6169 * FullMatrix<double> &M11,
6170 * FullMatrix<double> &M12,
6171 * FullMatrix<double> &M21,
6172 * FullMatrix<double> &M22,
6173 * const FEValuesBase<dim> &fe_faces0,
6174 * const std::vector<std::vector<double>> &recv_values,
6175 * const std::vector<std::vector<Tensor<1, dim>>> &recv_gradients,
6176 * const std::vector<double> &recv_jxws,
6177 * const double penalty_constant,
6181 * (recv_values.size() > 0 && recv_gradients.size() && recv_jxws.size()),
6182 * ExcMessage("Not possible to
assemble jumps and averages at a ghosted
"
6184 * const unsigned int dofs_per_cell = M11.m();
6185 * const std::vector<Tensor<1, dim>> &normals = fe_faces0.get_normal_vectors();
6186 * for (unsigned int q_index : fe_faces0.quadrature_point_indices())
6188 * const Tensor<1, dim> &normal = normals[q_index];
6189 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
6191 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
6193 * M11(i, j) += (-0.5 * fe_faces0.shape_grad(i, q_index) * normal *
6194 * fe_faces0.shape_value(j, q_index) -
6195 * 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
6196 * fe_faces0.shape_value(i, q_index) +
6197 * (penalty_constant / h_f) *
6198 * fe_faces0.shape_value(i, q_index) *
6199 * fe_faces0.shape_value(j, q_index)) *
6200 * fe_faces0.JxW(q_index);
6201 * M12(i, j) += (0.5 * fe_faces0.shape_grad(i, q_index) * normal *
6202 * recv_values[j][q_index] -
6203 * 0.5 * recv_gradients[j][q_index] * normal *
6204 * fe_faces0.shape_value(i, q_index) -
6205 * (penalty_constant / h_f) *
6206 * fe_faces0.shape_value(i, q_index) *
6207 * recv_values[j][q_index]) *
6208 * recv_jxws[q_index];
6210 * (-0.5 * recv_gradients[i][q_index] * normal *
6211 * fe_faces0.shape_value(j, q_index) +
6212 * 0.5 * fe_faces0.shape_grad(j, q_index) * normal *
6213 * recv_values[i][q_index] -
6214 * (penalty_constant / h_f) * recv_values[i][q_index] *
6215 * fe_faces0.shape_value(j, q_index)) *
6216 * recv_jxws[q_index];
6218 * (0.5 * recv_gradients[i][q_index] * normal *
6219 * recv_values[j][q_index] +
6220 * 0.5 * recv_gradients[j][q_index] * normal *
6221 * recv_values[i][q_index] +
6222 * (penalty_constant / h_f) * recv_values[i][q_index] *
6223 * recv_values[j][q_index]) *
6224 * recv_jxws[q_index];
6231 * * Utility function to assemble the SIPDG Laplace matrix.
6232 * * @note Supported matrix types are Trilinos types and native SparseMatrix
6233 * * objects provided by deal.II.
6235 * template <int dim, typename MatrixType>
6237 * assemble_dg_matrix(MatrixType &system_matrix,
6238 * const FiniteElement<dim> &fe_dg,
6239 * const AgglomerationHandler<dim> &ah)
6242 * (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix> ||
6243 * std::is_same_v<MatrixType,
6244 * SparseMatrix<typename MatrixType::value_type>>));
6246 * Assert((dynamic_cast<const FE_DGQ<dim> *>(&fe_dg) ||
6247 * dynamic_cast<const FE_DGP<dim> *>(&fe_dg) ||
6248 * dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_dg)),
6249 * ExcMessage("FE type not supported.
"));
6251 * AffineConstraints constraints;
6252 * constraints.close();
6253 * const double penalty_constant =
6254 * 10 * (fe_dg.degree + dim) * (fe_dg.degree + 1);
6255 * TrilinosWrappers::SparsityPattern dsp;
6256 * const_cast<AgglomerationHandler<dim> &>(ah)
6257 * .create_agglomeration_sparsity_pattern(dsp);
6258 * system_matrix.reinit(dsp);
6259 * const unsigned int dofs_per_cell = fe_dg.n_dofs_per_cell();
6260 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
6261 * FullMatrix<double> M11(dofs_per_cell, dofs_per_cell);
6262 * FullMatrix<double> M12(dofs_per_cell, dofs_per_cell);
6263 * FullMatrix<double> M21(dofs_per_cell, dofs_per_cell);
6264 * FullMatrix<double> M22(dofs_per_cell, dofs_per_cell);
6265 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
6266 * std::vector<types::global_dof_index> local_dof_indices_neighbor(
6269 * for (const auto &polytope : ah.polytope_iterators())
6271 * if (polytope->is_locally_owned())
6274 * const auto &agglo_values = ah.reinit(polytope);
6275 * for (unsigned int q_index : agglo_values.quadrature_point_indices())
6277 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
6279 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
6281 * cell_matrix(i, j) +=
6282 * agglo_values.shape_grad(i, q_index) *
6283 * agglo_values.shape_grad(j, q_index) *
6284 * agglo_values.JxW(q_index);
6290 * get volumetric DoFs
6293 * polytope->get_dof_indices(local_dof_indices);
6296 * Assemble face terms
6299 * unsigned int n_faces = polytope->n_faces();
6300 * const double h_f = polytope->diameter();
6301 * for (unsigned int f = 0; f < n_faces; ++f)
6303 * if (polytope->at_boundary(f))
6307 * Get normal vectors seen from each agglomeration.
6310 * const auto &fe_face = ah.reinit(polytope, f);
6311 * const auto &normals = fe_face.get_normal_vectors();
6312 * for (unsigned int q_index :
6313 * fe_face.quadrature_point_indices())
6315 * const Tensor<1, dim> &normal = normals[q_index];
6316 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
6318 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
6320 * cell_matrix(i, j) +=
6321 * (-fe_face.shape_value(i, q_index) *
6322 * fe_face.shape_grad(j, q_index) * normal -
6323 * fe_face.shape_grad(i, q_index) * normal *
6324 * fe_face.shape_value(j, q_index) +
6325 * (penalty_constant / h_f) *
6326 * fe_face.shape_value(i, q_index) *
6327 * fe_face.shape_value(j, q_index)) *
6328 * fe_face.JxW(q_index);
6335 * const auto &neigh_polytope = polytope->neighbor(f);
6336 * if (polytope->id() < neigh_polytope->id())
6338 * unsigned int nofn =
6339 * polytope->neighbor_of_agglomerated_neighbor(f);
6340 * Assert(neigh_polytope->neighbor(nofn)->id() ==
6342 * ExcMessage("Mismatch.
"));
6343 * const auto &fe_faces = ah.reinit_interface(
6344 * polytope, neigh_polytope, f, nofn);
6345 * const auto &fe_faces0 = fe_faces.first;
6346 * if (neigh_polytope->is_locally_owned())
6353 * const auto &fe_faces1 = fe_faces.second;
6358 * assemble_local_jumps_and_averages(M11,
6368 * distribute DoFs accordingly
6372 * neigh_polytope->get_dof_indices(
6373 * local_dof_indices_neighbor);
6374 * constraints.distribute_local_to_global(
6375 * M11, local_dof_indices, system_matrix);
6376 * constraints.distribute_local_to_global(
6378 * local_dof_indices,
6379 * local_dof_indices_neighbor,
6381 * constraints.distribute_local_to_global(
6383 * local_dof_indices_neighbor,
6384 * local_dof_indices,
6386 * constraints.distribute_local_to_global(
6387 * M22, local_dof_indices_neighbor, system_matrix);
6393 * neigh polytope is ghosted, so retrieve necessary
6397 * types::subdomain_id neigh_rank =
6398 * neigh_polytope->subdomain_id();
6399 * const auto &recv_jxws =
6400 * ah.recv_jxws.at(neigh_rank)
6401 * .at({neigh_polytope->id(), nofn});
6402 * const auto &recv_values =
6403 * ah.recv_values.at(neigh_rank)
6404 * .at({neigh_polytope->id(), nofn});
6405 * const auto &recv_gradients =
6406 * ah.recv_gradients.at(neigh_rank)
6407 * .at({neigh_polytope->id(), nofn});
6414 * there's no FEFaceValues on the other side (it's
6415 * ghosted), so we just pass the actual data we have
6416 * recevied from the neighboring ghosted polytope
6419 * assemble_local_jumps_and_averages_ghost(
6432 * distribute DoFs accordingly
6436 * neigh_polytope->get_dof_indices(
6437 * local_dof_indices_neighbor);
6438 * constraints.distribute_local_to_global(
6439 * M11, local_dof_indices, system_matrix);
6440 * constraints.distribute_local_to_global(
6442 * local_dof_indices,
6443 * local_dof_indices_neighbor,
6445 * constraints.distribute_local_to_global(
6447 * local_dof_indices_neighbor,
6448 * local_dof_indices,
6450 * constraints.distribute_local_to_global(
6451 * M22, local_dof_indices_neighbor, system_matrix);
6452 * } // ghosted polytope case
6454 * } // internal face
6456 * constraints.distribute_local_to_global(cell_matrix,
6457 * local_dof_indices,
6459 * } // locally owned polytopes
6461 * system_matrix.compress(VectorOperation::add);
6465 * * Compute SIPDG matrix as well as rhs vector.
6466 * * @note Hardcoded for f=1 and simplex elements.
6467 * * TODO: Pass Function object for boundary conditions and forcing term.
6469 * template <int dim, typename MatrixType, typename VectorType>
6471 * assemble_dg_matrix_on_standard_mesh(MatrixType &system_matrix,
6472 * VectorType &system_rhs,
6473 * const Mapping<dim> &mapping,
6474 * const FiniteElement<dim> &fe_dg,
6475 * const DoFHandler<dim> &dof_handler)
6478 * (std::is_same_v<MatrixType, TrilinosWrappers::SparseMatrix> ||
6479 * std::is_same_v<MatrixType,
6480 * SparseMatrix<typename MatrixType::value_type>>));
6482 * Assert((dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_dg) != nullptr),
6483 * ExcNotImplemented(
6484 * "Implemented only
for simplex meshes
for the time being.
"));
6486 * Assert(dof_handler.get_triangulation().all_reference_cells_are_simplex(),
6487 * ExcNotImplemented());
6489 * const double penalty_constant = .5 * fe_dg.degree * (fe_dg.degree + 1);
6490 * AffineConstraints<typename MatrixType::value_type> constraints;
6491 * constraints.close();
6493 * const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
6494 * const IndexSet locally_relevant_dofs =
6495 * DoFTools::extract_locally_relevant_dofs(dof_handler);
6497 * DynamicSparsityPattern dsp(locally_relevant_dofs);
6498 * DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
6499 * SparsityTools::distribute_sparsity_pattern(dsp,
6500 * dof_handler.locally_owned_dofs(),
6501 * dof_handler.get_communicator(),
6502 * locally_relevant_dofs);
6504 * system_matrix.reinit(locally_owned_dofs,
6505 * locally_owned_dofs,
6507 * dof_handler.get_communicator());
6509 * system_rhs.reinit(locally_owned_dofs, dof_handler.get_communicator());
6511 * const unsigned int quadrature_degree = fe_dg.degree + 1;
6512 * FEFaceValues<dim> fe_faces0(mapping,
6514 * QGaussSimplex<dim - 1>(quadrature_degree),
6515 * update_values | update_JxW_values |
6516 * update_gradients | update_quadrature_points |
6517 * update_normal_vectors);
6519 * FEValues<dim> fe_values(mapping,
6521 * QGaussSimplex<dim>(quadrature_degree),
6522 * update_values | update_JxW_values |
6523 * update_gradients | update_quadrature_points);
6525 * FEFaceValues<dim> fe_faces1(mapping,
6527 * QGaussSimplex<dim - 1>(quadrature_degree),
6528 * update_values | update_JxW_values |
6529 * update_gradients | update_quadrature_points |
6530 * update_normal_vectors);
6531 * const unsigned int dofs_per_cell = fe_dg.n_dofs_per_cell();
6533 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
6534 * Vector<double> cell_rhs(dofs_per_cell);
6536 * FullMatrix<double> M11(dofs_per_cell, dofs_per_cell);
6537 * FullMatrix<double> M12(dofs_per_cell, dofs_per_cell);
6538 * FullMatrix<double> M21(dofs_per_cell, dofs_per_cell);
6539 * FullMatrix<double> M22(dofs_per_cell, dofs_per_cell);
6541 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
6545 * Loop over standard deal.II cells
6548 * for (const auto &cell : dof_handler.active_cell_iterators())
6550 * if (cell->is_locally_owned())
6555 * fe_values.reinit(cell);
6559 * const auto &q_points = fe_values.get_quadrature_points();
6560 * const unsigned int n_qpoints = q_points.size();
6566 * for (unsigned int q_index : fe_values.quadrature_point_indices())
6568 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
6570 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
6572 * cell_matrix(i, j) += fe_values.shape_grad(i, q_index) *
6573 * fe_values.shape_grad(j, q_index) *
6574 * fe_values.JxW(q_index);
6577 * fe_values.shape_value(i, q_index) * 1. *
6578 * fe_values.JxW(q_index); // TODO: pass functional
6584 * distribute volumetric DoFs
6587 * cell->get_dof_indices(local_dof_indices);
6589 * for (const auto f : cell->face_indices())
6591 * const double extent1 =
6592 * cell->measure() / cell->face(f)->measure();
6594 * if (cell->face(f)->at_boundary())
6596 * hf = (1. / extent1 + 1. / extent1);
6597 * fe_faces0.reinit(cell, f);
6599 * const auto &normals = fe_faces0.get_normal_vectors();
6600 * for (unsigned int q_index :
6601 * fe_faces0.quadrature_point_indices())
6603 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
6605 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
6607 * cell_matrix(i, j) +=
6608 * (-fe_faces0.shape_value(i, q_index) *
6609 * fe_faces0.shape_grad(j, q_index) *
6610 * normals[q_index] -
6611 * fe_faces0.shape_grad(i, q_index) *
6612 * normals[q_index] *
6613 * fe_faces0.shape_value(j, q_index) +
6614 * (penalty_constant * hf) *
6615 * fe_faces0.shape_value(i, q_index) *
6616 * fe_faces0.shape_value(j, q_index)) *
6617 * fe_faces0.JxW(q_index);
6620 * 0.; // TODO: add bdary conditions functional
6626 * const auto &neigh_cell = cell->neighbor(f);
6627 * if (cell->global_active_cell_index() <
6628 * neigh_cell->global_active_cell_index())
6630 * const double extent2 =
6631 * neigh_cell->measure() /
6632 * neigh_cell->face(cell->neighbor_of_neighbor(f))
6634 * hf = (1. / extent1 + 1. / extent2);
6635 * fe_faces0.reinit(cell, f);
6636 * fe_faces1.reinit(neigh_cell,
6637 * cell->neighbor_of_neighbor(f));
6639 * std::vector<types::global_dof_index>
6640 * local_dof_indices_neighbor(dofs_per_cell);
6647 * const auto &normals = fe_faces0.get_normal_vectors();
6653 * for (unsigned int q_index :
6654 * fe_faces0.quadrature_point_indices())
6656 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
6658 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
6661 * (-0.5 * fe_faces0.shape_grad(i, q_index) *
6662 * normals[q_index] *
6663 * fe_faces0.shape_value(j, q_index) -
6664 * 0.5 * fe_faces0.shape_grad(j, q_index) *
6665 * normals[q_index] *
6666 * fe_faces0.shape_value(i, q_index) +
6667 * (penalty_constant * hf) *
6668 * fe_faces0.shape_value(i, q_index) *
6669 * fe_faces0.shape_value(j, q_index)) *
6670 * fe_faces0.JxW(q_index);
6673 * (0.5 * fe_faces0.shape_grad(i, q_index) *
6674 * normals[q_index] *
6675 * fe_faces1.shape_value(j, q_index) -
6676 * 0.5 * fe_faces1.shape_grad(j, q_index) *
6677 * normals[q_index] *
6678 * fe_faces0.shape_value(i, q_index) -
6679 * (penalty_constant * hf) *
6680 * fe_faces0.shape_value(i, q_index) *
6681 * fe_faces1.shape_value(j, q_index)) *
6682 * fe_faces1.JxW(q_index);
6690 * (-0.5 * fe_faces1.shape_grad(i, q_index) *
6691 * normals[q_index] *
6692 * fe_faces0.shape_value(j, q_index) +
6693 * 0.5 * fe_faces0.shape_grad(j, q_index) *
6694 * normals[q_index] *
6695 * fe_faces1.shape_value(i, q_index) -
6696 * (penalty_constant * hf) *
6697 * fe_faces1.shape_value(i, q_index) *
6698 * fe_faces0.shape_value(j, q_index)) *
6699 * fe_faces1.JxW(q_index);
6707 * (0.5 * fe_faces1.shape_grad(i, q_index) *
6708 * normals[q_index] *
6709 * fe_faces1.shape_value(j, q_index) +
6710 * 0.5 * fe_faces1.shape_grad(j, q_index) *
6711 * normals[q_index] *
6712 * fe_faces1.shape_value(i, q_index) +
6713 * (penalty_constant * hf) *
6714 * fe_faces1.shape_value(i, q_index) *
6715 * fe_faces1.shape_value(j, q_index)) *
6716 * fe_faces1.JxW(q_index);
6723 * distribute DoFs accordingly
6729 * neigh_cell->get_dof_indices(local_dof_indices_neighbor);
6731 * constraints.distribute_local_to_global(
6732 * M11, local_dof_indices, system_matrix);
6733 * constraints.distribute_local_to_global(
6735 * local_dof_indices,
6736 * local_dof_indices_neighbor,
6738 * constraints.distribute_local_to_global(
6740 * local_dof_indices_neighbor,
6741 * local_dof_indices,
6743 * constraints.distribute_local_to_global(
6744 * M22, local_dof_indices_neighbor, system_matrix);
6746 * } // check idx neighbors
6749 * constraints.distribute_local_to_global(cell_matrix,
6751 * local_dof_indices,
6756 * system_matrix.compress(VectorOperation::add);
6757 * system_rhs.compress(VectorOperation::add);
6760 * } // namespace ::PolyUtils
6766<a name="ann-source/agglomeration_handler.cc
"></a>
6767<h1>Annotated version of source/agglomeration_handler.cc</h1>
6773 * /* -----------------------------------------------------------------------------
6775 * * SPDX-License-Identifier: LGPL-2.1-or-later
6776 * * Copyright (C) 2025 by Marco Feder, Pasquale Claudio Africa, Xinping Gui,
6779 * * This file is part of the deal.II code gallery.
6781 * * -----------------------------------------------------------------------------
6784 * #include <deal.II/base/quadrature_lib.h>
6785 * #include <deal.II/lac/sparsity_tools.h>
6787 * #include <agglomeration_handler.h>
6789 * template <int dim, int spacedim>
6790 * AgglomerationHandler<dim, spacedim>::AgglomerationHandler(
6791 * const GridTools::Cache<dim, spacedim> &cache_tria)
6792 * : cached_tria(std::make_unique<GridTools::Cache<dim, spacedim>>(
6793 * cache_tria.get_triangulation(),
6794 * cache_tria.get_mapping()))
6795 * , communicator(cache_tria.get_triangulation().get_mpi_communicator())
6797 * Assert(dim == spacedim, ExcNotImplemented("Not available with codim > 0
"));
6798 * Assert(dim == 2 || dim == 3, ExcImpossibleInDim(1));
6799 * Assert((dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
6800 * &cached_tria->get_triangulation()) == nullptr),
6801 * ExcNotImplemented());
6802 * Assert(cached_tria->get_triangulation().n_active_cells() > 0,
6804 * "The triangulation must not be empty upon calling this function.
"));
6806 * n_agglomerations = 0;
6807 * hybrid_mesh = false;
6808 * initialize_agglomeration_data(cached_tria);
6813 * template <int dim, int spacedim>
6814 * typename AgglomerationHandler<dim, spacedim>::agglomeration_iterator
6815 * AgglomerationHandler<dim, spacedim>::define_agglomerate(
6816 * const AgglomerationContainer &cells)
6818 * Assert(cells.size() > 0, ExcMessage("No cells to be agglomerated.
"));
6820 * if (cells.size() == 1)
6821 * hybrid_mesh = true; // mesh is made also by classical cells
6825 * First index drives the selection of the master cell. After that, store the
6829 * const types::global_cell_index global_master_idx =
6830 * cells[0]->global_active_cell_index();
6831 * const types::global_cell_index master_idx = cells[0]->active_cell_index();
6832 * master_cells_container.push_back(cells[0]);
6833 * master_slave_relationships[global_master_idx] = -1;
6835 * const typename DoFHandler<dim>::active_cell_iterator cell_dh =
6836 * cells[0]->as_dof_handler_iterator(agglo_dh);
6837 * cell_dh->set_active_fe_index(CellAgglomerationType::master);
6841 * Store slave cells and save the relationship with the parent
6844 * std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
6846 * slaves.reserve(cells.size() - 1);
6849 * exclude first cell since it's the master cell
6852 * for (auto it = ++cells.begin(); it != cells.end(); ++it)
6854 * slaves.push_back(*it);
6855 * master_slave_relationships[(*it)->global_active_cell_index()] =
6856 * global_master_idx; // mark each slave
6857 * master_slave_relationships_iterators[(*it)->active_cell_index()] =
6860 * const typename DoFHandler<dim>::active_cell_iterator cell =
6861 * (*it)->as_dof_handler_iterator(agglo_dh);
6862 * cell->set_active_fe_index(CellAgglomerationType::slave); // slave cell
6866 * If we have a p::d::T, check that all cells are in the same subdomain.
6867 * If serial, just check that the subdomain_id is invalid.
6870 * Assert(((*it)->subdomain_id() == tria->locally_owned_subdomain() ||
6871 * tria->locally_owned_subdomain() == numbers::invalid_subdomain_id),
6872 * ExcInternalError());
6875 * master_slave_relationships_iterators[master_idx] =
6876 * cells[0]; // set iterator to master cell
6880 * Store the slaves of each master
6883 * master2slaves[master_idx] = slaves;
6886 * Save to which polygon this agglomerate correspond
6889 * master2polygon[master_idx] = n_agglomerations;
6891 * ++n_agglomerations; // an agglomeration has been performed, record it
6893 * create_bounding_box(cells); // fill the vector of bboxes
6897 * Finally, return a polygonal iterator to the polytope just constructed.
6900 * return {cells[0], this};
6903 * template <int dim, int spacedim>
6904 * typename AgglomerationHandler<dim, spacedim>::agglomeration_iterator
6905 * AgglomerationHandler<dim, spacedim>::define_agglomerate(
6906 * const AgglomerationContainer &cells,
6907 * const unsigned int fecollection_size)
6909 * Assert(cells.size() > 0, ExcMessage("No cells to be agglomerated.
"));
6911 * if (cells.size() == 1)
6912 * hybrid_mesh = true; // mesh is made also by classical cells
6916 * First index drives the selection of the master cell. After that, store the
6920 * const types::global_cell_index global_master_idx =
6921 * cells[0]->global_active_cell_index();
6922 * const types::global_cell_index master_idx = cells[0]->active_cell_index();
6923 * master_cells_container.push_back(cells[0]);
6924 * master_slave_relationships[global_master_idx] = -1;
6926 * const typename DoFHandler<dim>::active_cell_iterator cell_dh =
6927 * cells[0]->as_dof_handler_iterator(agglo_dh);
6928 * cell_dh->set_active_fe_index(CellAgglomerationType::master);
6932 * Store slave cells and save the relationship with the parent
6935 * std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
6937 * slaves.reserve(cells.size() - 1);
6940 * exclude first cell since it's the master cell
6943 * for (auto it = ++cells.begin(); it != cells.end(); ++it)
6945 * slaves.push_back(*it);
6946 * master_slave_relationships[(*it)->global_active_cell_index()] =
6947 * global_master_idx; // mark each slave
6948 * master_slave_relationships_iterators[(*it)->active_cell_index()] =
6951 * const typename DoFHandler<dim>::active_cell_iterator cell =
6952 * (*it)->as_dof_handler_iterator(agglo_dh);
6953 * cell->set_active_fe_index(
6954 * fecollection_size); // slave cell (the last index)
6958 * If we have a p::d::T, check that all cells are in the same subdomain.
6959 * If serial, just check that the subdomain_id is invalid.
6962 * Assert(((*it)->subdomain_id() == tria->locally_owned_subdomain() ||
6963 * tria->locally_owned_subdomain() == numbers::invalid_subdomain_id),
6964 * ExcInternalError());
6967 * master_slave_relationships_iterators[master_idx] =
6968 * cells[0]; // set iterator to master cell
6972 * Store the slaves of each master
6975 * master2slaves[master_idx] = slaves;
6978 * Save to which polygon this agglomerate correspond
6981 * master2polygon[master_idx] = n_agglomerations;
6983 * ++n_agglomerations; // an agglomeration has been performed, record it
6985 * create_bounding_box(cells); // fill the vector of bboxes
6989 * Finally, return a polygonal iterator to the polytope just constructed.
6992 * return {cells[0], this};
6996 * template <int dim, int spacedim>
6998 * AgglomerationHandler<dim, spacedim>::initialize_fe_values(
6999 * const Quadrature<dim> &cell_quadrature,
7000 * const UpdateFlags &flags,
7001 * const Quadrature<dim - 1> &face_quadrature,
7002 * const UpdateFlags &face_flags)
7004 * agglomeration_quad = cell_quadrature;
7005 * agglomeration_flags = flags;
7006 * agglomeration_face_quad = face_quadrature;
7007 * agglomeration_face_flags = face_flags | internal_agglomeration_face_flags;
7011 * std::make_unique<FEValues<dim>>(*mapping,
7013 * agglomeration_quad,
7014 * update_quadrature_points |
7015 * update_JxW_values); // only for quadrature
7016 * no_face_values = std::make_unique<FEFaceValues<dim>>(
7019 * agglomeration_face_quad,
7020 * update_quadrature_points | update_JxW_values |
7021 * update_normal_vectors); // only for quadrature
7024 * template <int dim, int spacedim>
7026 * AgglomerationHandler<dim, spacedim>::initialize_fe_values(
7027 * const hp::QCollection<dim> &cell_qcollection,
7028 * const UpdateFlags &flags,
7029 * const hp::QCollection<dim - 1> &face_qcollection,
7030 * const UpdateFlags &face_flags)
7032 * agglomeration_quad_collection = cell_qcollection;
7033 * agglomeration_flags = flags;
7034 * agglomeration_face_quad_collection = face_qcollection;
7035 * agglomeration_face_flags = face_flags | internal_agglomeration_face_flags;
7037 * mapping_collection = hp::MappingCollection<dim>(*mapping);
7038 * dummy_fe_collection = hp::FECollection<dim, spacedim>(dummy_fe);
7039 * hp_no_values = std::make_unique<hp::FEValues<dim>>(
7040 * mapping_collection,
7041 * dummy_fe_collection,
7042 * agglomeration_quad_collection,
7043 * update_quadrature_points | update_JxW_values); // only for quadrature
7045 * hp_no_face_values = std::make_unique<hp::FEFaceValues<dim>>(
7046 * mapping_collection,
7047 * dummy_fe_collection,
7048 * agglomeration_face_quad_collection,
7049 * update_quadrature_points | update_JxW_values |
7050 * update_normal_vectors); // only for quadrature
7055 * template <int dim, int spacedim>
7057 * AgglomerationHandler<dim, spacedim>::n_agglomerated_faces_per_cell(
7058 * const typename Triangulation<dim, spacedim>::active_cell_iterator &cell) const
7060 * unsigned int n_neighbors = 0;
7061 * for (const auto &f : cell->face_indices())
7063 * const auto &neighboring_cell = cell->neighbor(f);
7064 * if ((cell->face(f)->at_boundary()) ||
7065 * (neighboring_cell->is_active() &&
7066 * !are_cells_agglomerated(cell, neighboring_cell)))
7071 * return n_neighbors;
7076 * template <int dim, int spacedim>
7078 * AgglomerationHandler<dim, spacedim>::initialize_agglomeration_data(
7079 * const std::unique_ptr<GridTools::Cache<dim, spacedim>> &cache_tria)
7081 * tria = &(cache_tria->get_triangulation());
7082 * mapping = &(cache_tria->get_mapping());
7084 * agglo_dh.reinit(*tria);
7086 * if (const auto parallel_tria = dynamic_cast<
7087 * const ::parallel::TriangulationBase<dim, spacedim> *>(&*tria))
7089 * const std::weak_ptr<const Utilities::MPI::Partitioner> cells_partitioner =
7090 * parallel_tria->global_active_cell_index_partitioner();
7091 * master_slave_relationships.reinit(
7092 * cells_partitioner.lock()->locally_owned_range(), communicator);
7096 * master_slave_relationships.reinit(tria->n_active_cells(), MPI_COMM_SELF);
7099 * polytope_cache.clear();
7104 * First, update the pointer
7107 * cached_tria = std::make_unique<GridTools::Cache<dim, spacedim>>(
7108 * cache_tria->get_triangulation(), cache_tria->get_mapping());
7110 * connect_to_tria_signals();
7111 * n_agglomerations = 0;
7116 * template <int dim, int spacedim>
7118 * AgglomerationHandler<dim, spacedim>::distribute_agglomerated_dofs(
7119 * const FiniteElement<dim> &fe_space)
7121 * if (dynamic_cast<const FE_DGQ<dim> *>(&fe_space))
7122 * fe = std::make_unique<FE_DGQ<dim>>(fe_space.degree);
7123 * else if (dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_space))
7124 * fe = std::make_unique<FE_SimplexDGP<dim>>(fe_space.degree);
7128 * ExcNotImplemented(
7129 * "Currently, this interface supports only DGQ and DGP bases.
"));
7131 * box_mapping = std::make_unique<MappingBox<dim>>(
7133 * master2polygon); // construct bounding box mapping
7139 * the mesh is composed by standard and agglomerate cells. initialize
7140 * classes needed for standard cells in order to treat that finite
7141 * element space as defined on a standard shape and not on the
7145 * standard_scratch =
7146 * std::make_unique<ScratchData>(*mapping,
7148 * QGauss<dim>(2 * fe_space.degree + 2),
7149 * internal_agglomeration_flags);
7153 * fe_collection.push_back(*fe); // master
7154 * fe_collection.push_back(
7155 * FE_Nothing<dim, spacedim>(fe->reference_cell())); // slave
7157 * initialize_hp_structure();
7161 * in case the tria is distributed, communicate ghost information with
7165 * const bool needs_ghost_info =
7166 * dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&*tria) !=
7168 * if (needs_ghost_info)
7169 * setup_ghost_polytopes();
7171 * setup_connectivity_of_agglomeration();
7173 * if (needs_ghost_info)
7174 * exchange_interface_values();
7177 * template <int dim, int spacedim>
7179 * AgglomerationHandler<dim, spacedim>::distribute_agglomerated_dofs(
7180 * const hp::FECollection<dim, spacedim> &fe_collection_in)
7182 * is_hp_collection = true;
7184 * hp_fe_collection = std::make_unique<hp::FECollection<dim, spacedim>>(
7185 * fe_collection_in); // copy the input collection
7187 * box_mapping = std::make_unique<MappingBox<dim>>(
7189 * master2polygon); // construct bounding box mapping
7194 * AssertThrow(false,
7195 * ExcNotImplemented(
7196 * "Hybrid mesh is not implemented
for hp::FECollection.
"));
7199 * for (unsigned int i = 0; i < fe_collection_in.size(); ++i)
7201 * if (dynamic_cast<const FESystem<dim> *>(&fe_collection_in[i]))
7208 * for (unsigned int b = 0; b < fe_collection_in[i].n_base_elements();
7211 * if (!(dynamic_cast<const FE_DGQ<dim> *>(
7212 * &fe_collection_in[i].base_element(b)) ||
7213 * dynamic_cast<const FE_SimplexDGP<dim> *>(
7214 * &fe_collection_in[i].base_element(b)) ||
7215 * dynamic_cast<const FE_Nothing<dim> *>(
7216 * &fe_collection_in[i].base_element(b))))
7219 * ExcNotImplemented(
7220 * "Currently, this interface supports only DGQ and DGP bases.
"));
7230 * if (!(dynamic_cast<const FE_DGQ<dim> *>(&fe_collection_in[i]) ||
7231 * dynamic_cast<const FE_SimplexDGP<dim> *>(&fe_collection_in[i])))
7234 * ExcNotImplemented(
7235 * "Currently, this interface supports only DGQ and DGP bases.
"));
7237 * fe_collection.push_back(fe_collection_in[i]);
7240 * Assert(fe_collection[0].n_components() >= 1,
7242 * if (fe_collection[0].n_components() == 1)
7244 * fe_collection.push_back(FE_Nothing<dim, spacedim>());
7246 * else if (fe_collection[0].n_components() > 1)
7248 * std::vector<const FiniteElement<dim, spacedim> *> base_elements;
7249 * std::vector<unsigned int> multiplicities;
7250 * for (unsigned int b = 0; b < fe_collection[0].n_base_elements(); ++b)
7252 * base_elements.push_back(new FE_Nothing<dim, spacedim>());
7253 * multiplicities.push_back(fe_collection[0].element_multiplicity(b));
7255 * FESystem<dim, spacedim> fe_system_nothing(base_elements, multiplicities);
7256 * for (const auto *ptr : base_elements)
7258 * fe_collection.push_back(fe_system_nothing);
7261 * initialize_hp_structure();
7265 * in case the tria is distributed, communicate ghost information with
7269 * const bool needs_ghost_info =
7270 * dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&*tria) !=
7272 * if (needs_ghost_info)
7273 * setup_ghost_polytopes();
7275 * setup_connectivity_of_agglomeration();
7277 * if (needs_ghost_info)
7278 * exchange_interface_values();
7281 * template <int dim, int spacedim>
7283 * AgglomerationHandler<dim, spacedim>::create_bounding_box(
7284 * const AgglomerationContainer &polytope)
7286 * Assert(n_agglomerations > 0,
7287 * ExcMessage("No agglomeration has been performed.
"));
7288 * Assert(dim > 1, ExcNotImplemented());
7290 * std::vector<Point<spacedim>> pts; // store all the vertices
7291 * for (const auto &cell : polytope)
7292 * for (const auto i : cell->vertex_indices())
7293 * pts.push_back(cell->vertex(i));
7295 * bboxes.emplace_back(pts);
7300 * template <int dim, int spacedim>
7302 * AgglomerationHandler<dim, spacedim>::setup_connectivity_of_agglomeration()
7304 * Assert(master_cells_container.size() > 0,
7305 * ExcMessage("No agglomeration has been performed.
"));
7307 * agglo_dh.n_dofs() > 0,
7309 * "The
DoFHandler associated to the agglomeration has not been initialized.
"
7310 * "It
's likely that you forgot to distribute the DoFs. You may want"
7311 * "to check if a call to `initialize_hp_structure()` has been done."));
7313 * number_of_agglomerated_faces.resize(master2polygon.size(), 0);
7314 * for (const auto &cell : master_cells_container)
7316 * internal::AgglomerationHandlerImplementation<dim, spacedim>::
7317 * setup_master_neighbor_connectivity(cell, *this);
7320 * if (Utilities::MPI::job_supports_mpi())
7324 * communicate the number of faces
7327 * recv_n_faces = Utilities::MPI::some_to_some(communicator, local_n_faces);
7331 * send information about boundaries and neighboring polytopes id
7335 * Utilities::MPI::some_to_some(communicator, local_bdary_info);
7337 * recv_ghosted_master_id =
7338 * Utilities::MPI::some_to_some(communicator, local_ghosted_master_id);
7344 * template <int dim, int spacedim>
7346 * AgglomerationHandler<dim, spacedim>::exchange_interface_values()
7348 * const unsigned int dofs_per_cell = fe->dofs_per_cell;
7349 * for (const auto &polytope : polytope_iterators())
7351 * if (polytope->is_locally_owned())
7353 * const unsigned int n_faces = polytope->n_faces();
7354 * for (unsigned int f = 0; f < n_faces; ++f)
7356 * if (!polytope->at_boundary(f))
7358 * const auto &neigh_polytope = polytope->neighbor(f);
7359 * if (!neigh_polytope->is_locally_owned())
7363 * Neighboring polytope is ghosted.
7367 * Compute shape functions at the interface
7370 * const auto ¤t_fe = reinit(polytope, f);
7372 * std::vector<Point<spacedim>> qpoints_to_send =
7373 * current_fe.get_quadrature_points();
7375 * const std::vector<double> &jxws_to_send =
7376 * current_fe.get_JxW_values();
7378 * const std::vector<Tensor<1, spacedim>> &normals_to_send =
7379 * current_fe.get_normal_vectors();
7382 * const types::subdomain_id neigh_rank =
7383 * neigh_polytope->subdomain_id();
7385 * std::pair<CellId, unsigned int> cell_and_face{
7386 * polytope->id(), f};
7389 * Prepare data to send
7392 * local_qpoints[neigh_rank].emplace(cell_and_face,
7395 * local_jxws[neigh_rank].emplace(cell_and_face,
7398 * local_normals[neigh_rank].emplace(cell_and_face,
7402 * const unsigned int n_qpoints = qpoints_to_send.size();
7406 * TODO: check `agglomeration_flags` before computing
7407 * values and gradients.
7410 * std::vector<std::vector<double>> values_per_qpoints(
7413 * std::vector<std::vector<Tensor<1, spacedim>>>
7414 * gradients_per_qpoints(dofs_per_cell);
7416 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
7418 * values_per_qpoints[i].resize(n_qpoints);
7419 * gradients_per_qpoints[i].resize(n_qpoints);
7420 * for (unsigned int q = 0; q < n_qpoints; ++q)
7422 * values_per_qpoints[i][q] =
7423 * current_fe.shape_value(i, q);
7424 * gradients_per_qpoints[i][q] =
7425 * current_fe.shape_grad(i, q);
7429 * local_values[neigh_rank].emplace(cell_and_face,
7430 * values_per_qpoints);
7431 * local_gradients[neigh_rank].emplace(
7432 * cell_and_face, gradients_per_qpoints);
7441 * Finally, exchange with neighboring ranks
7444 * recv_qpoints = Utilities::MPI::some_to_some(communicator, local_qpoints);
7445 * recv_jxws = Utilities::MPI::some_to_some(communicator, local_jxws);
7446 * recv_normals = Utilities::MPI::some_to_some(communicator, local_normals);
7447 * recv_values = Utilities::MPI::some_to_some(communicator, local_values);
7448 * recv_gradients = Utilities::MPI::some_to_some(communicator, local_gradients);
7453 * template <int dim, int spacedim>
7455 * AgglomerationHandler<dim, spacedim>::agglomerated_quadrature(
7456 * const typename AgglomerationHandler<dim, spacedim>::AgglomerationContainer
7458 * const typename Triangulation<dim, spacedim>::active_cell_iterator
7459 * &master_cell) const
7461 * Assert(is_master_cell(master_cell),
7462 * ExcMessage("This must be a master cell."));
7464 * std::vector<Point<dim>> vec_pts;
7465 * std::vector<double> vec_JxWs;
7467 * if (!is_hp_collection)
7471 * Original version: handle case without hp::FECollection
7474 * for (const auto &dummy_cell : cells)
7476 * no_values->reinit(dummy_cell);
7477 * auto q_points = no_values->get_quadrature_points(); // real qpoints
7478 * const auto &JxWs = no_values->get_JxW_values();
7480 * std::transform(q_points.begin(),
7482 * std::back_inserter(vec_pts),
7483 * [&](const Point<spacedim> &p) { return p; });
7484 * std::transform(JxWs.begin(),
7486 * std::back_inserter(vec_JxWs),
7487 * [&](const double w) { return w; });
7494 * Handle the hp::FECollection case
7497 * const auto &master_cell_as_dh_iterator =
7498 * master_cell->as_dof_handler_iterator(agglo_dh);
7499 * for (const auto &dummy_cell : cells)
7503 * The following verbose call is necessary to handle cases where
7504 * different slave cells on different polytopes use different
7505 * quadrature rules. If the hp::QCollection contains multiple
7506 * elements, calling hp_no_values->reinit(dummy_cell) won't work
7507 * because it cannot infer the correct quadrature rule. By explicitly
7508 * passing the active FE
index as q_index, and setting mapping_index
7509 * and
fe_index to 0, we ensure that the dummy cell uses the same
7510 * quadrature rule as its corresponding master cell. This assumes a
7511 *
one-to-
one correspondence between
hp::QCollection and
7512 *
hp::FECollection, which is the convention in deal.II. However, this
7513 * implementation does not support cases where
hp::QCollection and
7514 *
hp::FECollection have different sizes.
7515 * TODO: Refactor the architecture to better handle numerical
7516 * integration
for hp::QCollection.
7519 *
hp_no_values->
reinit(dummy_cell,
7520 *
master_cell_as_dh_iterator->active_fe_index(),
7523 *
auto q_points = hp_no_values->get_present_fe_values()
7524 *
.get_quadrature_points();
7525 *
const auto &JxWs =
7526 *
hp_no_values->get_present_fe_values().get_JxW_values();
7528 *
std::transform(q_points.begin(),
7530 *
std::back_inserter(vec_pts),
7532 *
std::transform(JxWs.begin(),
7534 *
std::back_inserter(vec_JxWs),
7535 *
[&](
const double w) { return w; });
7541 * Map back each
point in real space by
using the map associated to the
7545 *
std::vector<Point<dim>> unit_points(vec_pts.size());
7546 *
const auto &bbox =
7547 *
bboxes[master2polygon.at(master_cell->active_cell_index())];
7548 *
unit_points.reserve(vec_pts.size());
7550 *
for (
unsigned int i = 0; i < vec_pts.size(); i++)
7551 *
unit_points[i] = bbox.real_to_unit(vec_pts[i]);
7558 *
template <
int dim,
int spacedim>
7560 *
AgglomerationHandler<dim, spacedim>::initialize_hp_structure()
7562 *
Assert(agglo_dh.get_triangulation().n_cells() > 0,
7564 *
"Triangulation must not be empty upon calling this function."));
7565 *
Assert(n_agglomerations > 0,
7566 *
ExcMessage(
"No agglomeration has been performed."));
7568 *
agglo_dh.distribute_dofs(fe_collection);
7571 * euler_mapping = std::make_unique<
7580 *
template <
int dim,
int spacedim>
7582 *
AgglomerationHandler<dim, spacedim>::reinit(
7583 *
const AgglomerationIterator<dim, spacedim> &polytope)
const
7588 * ExcMessage(
"The mapping describing the physical element stemming
7590 *
"agglomeration has not been set up."));
7596 *
const auto &deal_cell = polytope->as_dof_handler_iterator(agglo_dh);
7600 * First
check if the polytope is made just by a single cell. If so, use
7602 *
if (polytope->n_background_cells() == 1)
7603 *
return standard_scratch->
reinit(deal_cell);
7609 *
const auto &agglo_cells = polytope->get_agglomerate();
7611 *
Quadrature<dim> agglo_quad = agglomerated_quadrature(agglo_cells, deal_cell);
7613 *
if (!is_hp_collection)
7620 *
agglomerated_scratch = std::make_unique<ScratchData>(*box_mapping,
7623 *
agglomeration_flags);
7632 *
agglomerated_scratch = std::make_unique<ScratchData>(*box_mapping,
7633 *
polytope->get_fe(),
7635 *
agglomeration_flags);
7637 *
return agglomerated_scratch->reinit(deal_cell);
7642 *
template <
int dim,
int spacedim>
7644 *
AgglomerationHandler<dim, spacedim>::reinit_master(
7646 *
const unsigned int face_index,
7648 *
&agglo_isv_ptr)
const
7650 *
return internal::AgglomerationHandlerImplementation<dim, spacedim>::
7651 *
reinit_master(cell, face_index, agglo_isv_ptr, *
this);
7656 *
template <
int dim,
int spacedim>
7658 *
AgglomerationHandler<dim, spacedim>::reinit(
7659 *
const AgglomerationIterator<dim, spacedim> &polytope,
7660 *
const unsigned int face_index)
const
7665 * ExcMessage(
"The mapping describing the physical element stemming
7667 *
"agglomeration has not been set up."));
7673 *
const auto &deal_cell = polytope->as_dof_handler_iterator(agglo_dh);
7674 *
Assert(is_master_cell(deal_cell), ExcMessage(
"This should be true."));
7676 *
return internal::AgglomerationHandlerImplementation<dim, spacedim>::
7677 *
reinit_master(deal_cell, face_index, agglomerated_isv_bdary, *
this);
7682 *
template <
int dim,
int spacedim>
7683 *
std::pair<const FEValuesBase<dim, spacedim> &,
7685 *
AgglomerationHandler<dim, spacedim>::reinit_interface(
7686 *
const AgglomerationIterator<dim, spacedim> &polytope_in,
7687 *
const AgglomerationIterator<dim, spacedim> &neigh_polytope,
7688 *
const unsigned int local_in,
7689 *
const unsigned int local_neigh)
const
7693 * If current and neighboring polytopes are both locally owned, then compute
7694 * the jump in the classical way without needing information about ghosted
7698 *
if (polytope_in->is_locally_owned() && neigh_polytope->is_locally_owned())
7700 *
const auto &cell_in = polytope_in->as_dof_handler_iterator(agglo_dh);
7701 *
const auto &neigh_cell =
7702 *
neigh_polytope->as_dof_handler_iterator(agglo_dh);
7704 *
const auto &fe_in =
7705 *
internal::AgglomerationHandlerImplementation<dim, spacedim>::
7706 *
reinit_master(cell_in, local_in, agglomerated_isv, *
this);
7707 *
const auto &fe_out =
7708 *
internal::AgglomerationHandlerImplementation<dim, spacedim>::
7709 *
reinit_master(neigh_cell, local_neigh, agglomerated_isv_neigh, *
this);
7710 *
std::pair<const FEValuesBase<dim, spacedim> &,
7712 *
my_p(fe_in, fe_out);
7718 *
Assert((polytope_in->is_locally_owned() &&
7719 *
!neigh_polytope->is_locally_owned()),
7720 *
ExcInternalError());
7722 *
const auto &cell = polytope_in->as_dof_handler_iterator(agglo_dh);
7723 *
const auto &bbox = bboxes[master2polygon.at(cell->active_cell_index())];
7726 *
const double bbox_measure = bbox.volume();
7732 *
const unsigned int neigh_rank = neigh_polytope->subdomain_id();
7733 *
const CellId &neigh_id = neigh_polytope->id();
7737 * Retrieve qpoints,JxWs, normals sent previously from the neighboring
7741 *
std::vector<Point<spacedim>> &real_qpoints =
7742 *
recv_qpoints.at(neigh_rank).at({neigh_id, local_neigh});
7744 *
const auto &JxWs = recv_jxws.at(neigh_rank).at({neigh_id, local_neigh});
7746 *
std::vector<Tensor<1, spacedim>> &normals =
7747 *
recv_normals.at(neigh_rank).at({neigh_id, local_neigh});
7751 * Apply the necessary scalings due to the bbox.
7754 *
std::vector<Point<spacedim>> final_unit_q_points;
7755 *
std::transform(real_qpoints.begin(),
7756 *
real_qpoints.end(),
7757 *
std::back_inserter(final_unit_q_points),
7759 * return bbox.real_to_unit(p);
7764 * std::vector<double> scale_factors(final_unit_q_points.size());
7765 * std::vector<double> scaled_weights(final_unit_q_points.size());
7766 * std::vector<Tensor<1, dim>> scaled_normals(final_unit_q_points.size());
7770 * Since we received normal vectors from a neighbor, we have to
swap
7773 *
for (
unsigned int q = 0; q < final_unit_q_points.size(); ++q)
7775 *
for (
unsigned int direction = 0; direction < spacedim; ++direction)
7776 * scaled_normals[q][direction] =
7777 * normals[q][direction] * (bbox.side_length(direction));
7781 * scaled_normals[q] *= -1;
7785 * scaled_weights[q] =
7786 * (JxWs[q] * scaled_normals[q].norm()) / bbox_measure;
7787 * scaled_normals[q] /= scaled_normals[q].norm();
7791 *
for (
unsigned int q = 0; q < final_unit_q_points.size(); ++q)
7796 *
final_unit_q_points, JxWs, normals);
7798 *
agglomerated_isv =
7799 *
std::make_unique<NonMatching::FEImmersedSurfaceValues<spacedim>>(
7800 *
*box_mapping, *fe, surface_quad, agglomeration_face_flags);
7803 *
agglomerated_isv->reinit(cell);
7805 *
std::pair<const FEValuesBase<dim, spacedim> &,
7807 *
my_p(*agglomerated_isv, *agglomerated_isv);
7815 *
template <
int dim,
int spacedim>
7816 *
template <
typename SparsityPatternType,
typename Number>
7818 *
AgglomerationHandler<dim, spacedim>::create_agglomeration_sparsity_pattern(
7819 *
SparsityPatternType &dsp,
7821 *
const bool keep_constrained_dofs,
7824 *
Assert(n_agglomerations > 0,
7825 *
ExcMessage(
"The agglomeration has not been set up correctly."));
7828 *
"The Sparsity pattern must be empty upon calling this function."));
7830 *
const IndexSet &locally_owned_dofs = agglo_dh.locally_owned_dofs();
7831 *
const IndexSet locally_relevant_dofs =
7834 *
if constexpr (std::is_same_v<SparsityPatternType, DynamicSparsityPattern>)
7835 *
dsp.reinit(locally_owned_dofs.size(),
7836 *
locally_owned_dofs.size(),
7837 *
locally_relevant_dofs);
7838 *
else if constexpr (std::is_same_v<SparsityPatternType,
7840 *
dsp.reinit(locally_owned_dofs, communicator);
7846 * Create the sparsity pattern corresponding only to volumetric terms. The
7847 * fluxes needed by DG methods will be filled later.
7851 *
agglo_dh, dsp, constraints, keep_constrained_dofs, subdomain_id);
7854 *
if (!is_hp_collection)
7861 *
const unsigned int dofs_per_cell = agglo_dh.get_fe(0).n_dofs_per_cell();
7862 *
std::vector<types::global_dof_index> current_dof_indices(dofs_per_cell);
7863 *
std::vector<types::global_dof_index> neighbor_dof_indices(dofs_per_cell);
7867 * Loop over all locally owned polytopes, find the neighbor (also ghosted)
7868 * and add fluxes to the sparsity pattern.
7871 *
for (
const auto &polytope : polytope_iterators())
7873 *
if (polytope->is_locally_owned())
7875 *
const unsigned
int n_current_faces = polytope->n_faces();
7876 *
polytope->get_dof_indices(current_dof_indices);
7877 *
for (
unsigned int f = 0; f < n_current_faces; ++f)
7879 *
const auto &neigh_polytope = polytope->neighbor(f);
7882 *
neigh_polytope->get_dof_indices(neighbor_dof_indices);
7883 *
constraints.add_entries_local_to_global(
7884 *
current_dof_indices,
7885 *
neighbor_dof_indices,
7887 *
keep_constrained_dofs,
7902 * Loop over all locally owned polytopes, find the neighbor (also ghosted)
7903 * and add fluxes to the sparsity pattern.
7906 *
for (
const auto &polytope : polytope_iterators())
7908 *
if (polytope->is_locally_owned())
7910 *
const unsigned
int current_dofs_per_cell =
7911 *
polytope->get_fe().dofs_per_cell;
7912 *
std::vector<types::global_dof_index> current_dof_indices(
7913 *
current_dofs_per_cell);
7915 *
const unsigned int n_current_faces = polytope->n_faces();
7916 *
polytope->get_dof_indices(current_dof_indices);
7917 *
for (
unsigned int f = 0; f < n_current_faces; ++f)
7919 *
const auto &neigh_polytope = polytope->neighbor(f);
7922 *
const unsigned int neighbor_dofs_per_cell =
7923 *
neigh_polytope->get_fe().dofs_per_cell;
7924 *
std::vector<types::global_dof_index> neighbor_dof_indices(
7925 *
neighbor_dofs_per_cell);
7927 *
neigh_polytope->get_dof_indices(neighbor_dof_indices);
7928 *
constraints.add_entries_local_to_global(
7929 *
current_dof_indices,
7930 *
neighbor_dof_indices,
7932 *
keep_constrained_dofs,
7942 *
if constexpr (std::is_same_v<SparsityPatternType,
7949 *
template <
int dim,
int spacedim>
7951 *
AgglomerationHandler<dim, spacedim>::setup_ghost_polytopes()
7953 *
[[maybe_unused]]
const auto parallel_triangulation =
7955 *
Assert(parallel_triangulation !=
nullptr, ExcInternalError());
7957 *
const unsigned int n_dofs_per_cell = fe->dofs_per_cell;
7958 *
std::vector<types::global_dof_index> global_dof_indices(n_dofs_per_cell);
7959 *
for (
const auto &polytope : polytope_iterators())
7960 *
if (polytope->is_locally_owned())
7964 *
const auto polytope_dh = polytope->as_dof_handler_iterator(agglo_dh);
7965 *
polytope_dh->get_dof_indices(global_dof_indices);
7968 *
const auto &agglomerate = polytope->get_agglomerate();
7970 *
for (
const auto &cell : agglomerate)
7974 * interior, locally owned, cell
7977 *
for (
const auto &f : cell->face_indices())
7979 *
if (!cell->at_boundary(f))
7981 *
const auto &neighbor = cell->neighbor(f);
7982 *
if (neighbor->is_ghost())
7986 * key of the map: the rank to which send the
data
7990 *
neighbor->subdomain_id();
7994 * inform the
"standard" neighbor about the neighboring
7995 *
id and its master cell
7998 *
local_cell_ids_neigh_cell[neigh_rank].emplace(
7999 *
cell->id(), master_cell_id);
8003 * inform the neighboring rank that
this master cell
8004 * (hence polytope) has the following DoF indices
8007 *
local_ghost_dofs[neigh_rank].emplace(
8008 *
master_cell_id, global_dof_indices);
8012 * ...same
for bounding boxes
8015 *
const auto &bbox = bboxes[polytope->index()];
8016 *
local_ghosted_bbox[neigh_rank].emplace(master_cell_id,
8024 *
recv_cell_ids_neigh_cell =
8029 * Exchange with neighboring ranks the neighboring bounding boxes
8032 *
recv_ghosted_bbox =
8037 * Exchange with neighboring ranks the neighboring ghosted DoFs
8050 *
template <
int dim,
int spacedim>
8051 *
class AgglomerationHandlerImplementation
8057 *
const unsigned int face_index,
8060 *
const AgglomerationHandler<dim, spacedim> &handler)
8062 *
Assert(handler.is_master_cell(cell),
8063 *
ExcMessage(
"This cell must be a master one."));
8065 *
AgglomerationIterator<dim, spacedim> it{cell, &handler};
8066 *
const auto &neigh_polytope = it->neighbor(face_index);
8068 *
const CellId polytope_in_id = cell->id();
8072 * Retrieve the bounding box of the agglomeration
8075 *
const auto &bbox =
8076 *
handler.bboxes[handler.master2polygon.at(cell->active_cell_index())];
8078 *
CellId polytope_out_id;
8080 *
polytope_out_id = neigh_polytope->id();
8082 *
polytope_out_id = polytope_in_id;
8084 *
const auto &common_face = handler.polytope_cache.interface.at(
8085 *
{polytope_in_id, polytope_out_id});
8087 *
std::vector<Point<spacedim>> final_unit_q_points;
8088 *
std::vector<double> final_weights;
8089 *
std::vector<Tensor<1, dim>> final_normals;
8091 *
if (!handler.is_hp_collection)
8098 *
const unsigned int expected_qpoints =
8099 *
common_face.
size() * handler.agglomeration_face_quad.size();
8100 *
final_unit_q_points.reserve(expected_qpoints);
8101 *
final_weights.reserve(expected_qpoints);
8102 *
final_normals.reserve(expected_qpoints);
8105 *
for (
const auto &[deal_cell, local_face_idx] : common_face)
8107 *
handler.no_face_values->
reinit(deal_cell, local_face_idx);
8109 *
const auto &q_points =
8110 *
handler.no_face_values->get_quadrature_points();
8111 *
const auto &JxWs = handler.no_face_values->get_JxW_values();
8112 *
const auto &normals =
8113 *
handler.no_face_values->get_normal_vectors();
8115 *
const unsigned int n_qpoints_agglo = q_points.size();
8117 *
for (
unsigned int q = 0; q < n_qpoints_agglo; ++q)
8119 *
final_unit_q_points.push_back(
8120 *
bbox.real_to_unit(q_points[q]));
8121 *
final_weights.push_back(JxWs[q]);
8122 *
final_normals.push_back(normals[q]);
8133 *
unsigned int higher_order_quad_index = cell->active_fe_index();
8136 *
.agglomeration_face_quad_collection[cell->active_fe_index()]
8139 *
.agglomeration_face_quad_collection[neigh_polytope
8140 *
->active_fe_index()]
8142 *
higher_order_quad_index = neigh_polytope->active_fe_index();
8144 *
const unsigned int expected_qpoints =
8145 *
common_face.
size() *
8147 *
.agglomeration_face_quad_collection[higher_order_quad_index]
8149 *
final_unit_q_points.reserve(expected_qpoints);
8150 *
final_weights.reserve(expected_qpoints);
8151 *
final_normals.reserve(expected_qpoints);
8153 *
for (
const auto &[deal_cell, local_face_idx] : common_face)
8155 *
handler.hp_no_face_values->
reinit(
8156 *
deal_cell, local_face_idx, higher_order_quad_index, 0, 0);
8158 *
const auto &q_points =
8159 *
handler.hp_no_face_values->get_present_fe_values()
8160 *
.get_quadrature_points();
8161 *
const auto &JxWs =
8162 *
handler.hp_no_face_values->get_present_fe_values()
8163 *
.get_JxW_values();
8164 *
const auto &normals =
8165 *
handler.hp_no_face_values->get_present_fe_values()
8166 *
.get_normal_vectors();
8168 *
const unsigned int n_qpoints_agglo = q_points.size();
8170 *
for (
unsigned int q = 0; q < n_qpoints_agglo; ++q)
8172 *
final_unit_q_points.push_back(
8173 *
bbox.real_to_unit(q_points[q]));
8174 *
final_weights.push_back(JxWs[q]);
8175 *
final_normals.push_back(normals[q]);
8182 *
final_unit_q_points, final_weights, final_normals);
8184 *
if (!handler.is_hp_collection)
8187 *
std::make_unique<NonMatching::FEImmersedSurfaceValues<spacedim>>(
8188 *
*(handler.box_mapping),
8191 *
handler.agglomeration_face_flags);
8196 *
std::make_unique<NonMatching::FEImmersedSurfaceValues<spacedim>>(
8197 *
*(handler.box_mapping),
8200 *
handler.agglomeration_face_flags);
8203 *
agglo_isv_ptr->reinit(cell);
8205 *
return *agglo_isv_ptr;
8217 *
setup_master_neighbor_connectivity(
8220 *
const AgglomerationHandler<dim, spacedim> &handler)
8223 *
handler.master_slave_relationships[master_cell
8224 *
->global_active_cell_index()] ==
8226 *
ExcMessage(
"The present cell with index " +
8227 *
std::to_string(master_cell->global_active_cell_index()) +
8228 *
"is not a master one."));
8230 *
const auto &agglomeration = handler.get_agglomerate(master_cell);
8232 *
handler.master2polygon.at(master_cell->active_cell_index());
8234 *
CellId current_polytope_id = master_cell->id();
8237 *
std::set<types::global_cell_index> visited_polygonal_neighbors;
8239 *
std::map<unsigned int, CellId> face_to_neigh_id;
8241 *
std::map<unsigned int, bool> is_face_at_boundary;
8245 * same as above, but with
CellId
8248 *
std::set<CellId> visited_polygonal_neighbors_id;
8249 *
unsigned int ghost_counter = 0;
8251 *
for (
const auto &cell : agglomeration)
8254 *
cell->active_cell_index();
8256 *
const CellId cell_id = cell->id();
8258 *
for (
const auto f : cell->face_indices())
8260 *
const auto &neighboring_cell = cell->neighbor(f);
8262 *
const bool valid_neighbor =
8265 *
if (valid_neighbor)
8267 *
if (neighboring_cell->is_locally_owned() &&
8268 *
!handler.are_cells_agglomerated(cell, neighboring_cell))
8272 * - cell is not on the boundary,
8273 * - it
's not agglomerated with the neighbor. If so,
8274 * it's a neighbor of the present agglomeration
8275 * std::cout <<
" (from rank) "
8277 * handler.communicator)
8283 * <<
"neighbor locally owned? " << std::boolalpha
8284 * << neighboring_cell->is_locally_owned() <<
8286 *
if (neighboring_cell->is_ghost())
8287 * handler.ghosted_indices.push_back(
8288 * neighboring_cell->active_cell_index());
8292 * a
new face of the agglomeration has been
8296 *
handler.polygon_boundary[master_cell].push_back(
8301 * global
index of neighboring deal.II cell
8305 *
neighboring_cell->active_cell_index();
8309 * master cell
for the neighboring polytope
8312 *
const auto &master_of_neighbor =
8313 *
handler.master_slave_relationships_iterators.at(
8314 *
neighboring_cell_index);
8316 *
const auto nof = cell->neighbor_of_neighbor(f);
8318 *
if (handler.is_slave_cell(neighboring_cell))
8322 *
index of the neighboring polytope
8326 *
neighbor_polytope_index =
8327 *
handler.master2polygon.at(
8328 *
master_of_neighbor->active_cell_index());
8330 *
CellId neighbor_polytope_id =
8331 *
master_of_neighbor->id();
8333 *
if (visited_polygonal_neighbors.find(
8334 *
neighbor_polytope_index) ==
8335 *
std::end(visited_polygonal_neighbors))
8345 *
const unsigned int n_face =
8346 *
handler.number_of_agglomerated_faces
8347 *
[current_polytope_index];
8349 *
handler.polytope_cache.cell_face_at_boundary[{
8350 *
current_polytope_index, n_face}] = {
8351 *
false, master_of_neighbor};
8353 *
is_face_at_boundary[n_face] =
true;
8355 *
++handler.number_of_agglomerated_faces
8356 *
[current_polytope_index];
8358 *
visited_polygonal_neighbors.insert(
8359 *
neighbor_polytope_index);
8363 *
if (handler.polytope_cache.visited_cell_and_faces
8364 *
.find({cell_index, f}) ==
8365 *
std::end(handler.polytope_cache
8366 *
.visited_cell_and_faces))
8368 *
handler.polytope_cache
8369 *
.interface[{current_polytope_id,
8370 *
neighbor_polytope_id}]
8371 *
.emplace_back(cell, f);
8373 *
handler.polytope_cache.visited_cell_and_faces
8378 *
if (handler.polytope_cache.visited_cell_and_faces
8379 *
.find({neighboring_cell_index, nof}) ==
8380 *
std::end(handler.polytope_cache
8381 *
.visited_cell_and_faces))
8383 *
handler.polytope_cache
8384 *
.interface[{neighbor_polytope_id,
8385 *
current_polytope_id}]
8386 *
.emplace_back(neighboring_cell, nof);
8388 *
handler.polytope_cache.visited_cell_and_faces
8389 *
.insert({neighboring_cell_index, nof});
8396 * neighboring cell is a master
8400 *
save the pair of neighboring cells
8404 *
neighbor_polytope_index =
8405 *
handler.master2polygon.at(
8406 *
neighboring_cell_index);
8408 *
CellId neighbor_polytope_id =
8409 *
neighboring_cell->id();
8411 *
if (visited_polygonal_neighbors.find(
8412 *
neighbor_polytope_index) ==
8413 *
std::end(visited_polygonal_neighbors))
8420 *
const unsigned int n_face =
8421 *
handler.number_of_agglomerated_faces
8422 *
[current_polytope_index];
8425 *
handler.polytope_cache.cell_face_at_boundary[{
8426 *
current_polytope_index, n_face}] = {
8427 *
false, neighboring_cell};
8429 *
is_face_at_boundary[n_face] =
true;
8431 *
++handler.number_of_agglomerated_faces
8432 *
[current_polytope_index];
8434 *
visited_polygonal_neighbors.insert(
8435 *
neighbor_polytope_index);
8440 *
if (handler.polytope_cache.visited_cell_and_faces
8441 *
.find({cell_index, f}) ==
8442 *
std::end(handler.polytope_cache
8443 *
.visited_cell_and_faces))
8445 *
handler.polytope_cache
8446 *
.interface[{current_polytope_id,
8447 *
neighbor_polytope_id}]
8448 *
.emplace_back(cell, f);
8450 *
handler.polytope_cache.visited_cell_and_faces
8454 *
if (handler.polytope_cache.visited_cell_and_faces
8455 *
.find({neighboring_cell_index, nof}) ==
8456 *
std::end(handler.polytope_cache
8457 *
.visited_cell_and_faces))
8459 *
handler.polytope_cache
8460 *
.interface[{neighbor_polytope_id,
8461 *
current_polytope_id}]
8462 *
.emplace_back(neighboring_cell, nof);
8464 *
handler.polytope_cache.visited_cell_and_faces
8465 *
.insert({neighboring_cell_index, nof});
8469 *
else if (neighboring_cell->is_ghost())
8471 *
const auto nof = cell->neighbor_of_neighbor(f);
8475 * from neighboring rank,receive the association
8476 * between standard cell ids and neighboring polytope.
8477 * This tells to the current rank that the
8478 * neighboring cell has the following
CellId as master
8482 *
const auto &check_neigh_poly_ids =
8483 *
handler.recv_cell_ids_neigh_cell.at(
8484 *
neighboring_cell->subdomain_id());
8486 *
const CellId neighboring_cell_id =
8487 *
neighboring_cell->id();
8489 *
const CellId &check_neigh_polytope_id =
8490 *
check_neigh_poly_ids.at(neighboring_cell_id);
8494 *
const auto master_index =
8495 * master_indices[ghost_counter];
8501 *
if (visited_polygonal_neighbors_id.find(
8502 *
check_neigh_polytope_id) ==
8503 *
std::end(visited_polygonal_neighbors_id))
8505 *
handler.polytope_cache.cell_face_at_boundary[{
8506 *
current_polytope_index,
8507 *
handler.number_of_agglomerated_faces
8508 *
[current_polytope_index]}] = {
false,
8509 *
neighboring_cell};
8514 * record the cell
id of the neighboring polytope
8517 *
handler.polytope_cache.ghosted_master_id[{
8518 *
current_polytope_id,
8519 *
handler.number_of_agglomerated_faces
8520 *
[current_polytope_index]}] =
8521 *
check_neigh_polytope_id;
8524 *
const unsigned int n_face =
8525 *
handler.number_of_agglomerated_faces
8526 *
[current_polytope_index];
8528 *
face_to_neigh_id[n_face] = check_neigh_polytope_id;
8530 *
is_face_at_boundary[n_face] =
false;
8535 * increment number of faces
8538 *
++handler.number_of_agglomerated_faces
8539 *
[current_polytope_index];
8541 *
visited_polygonal_neighbors_id.insert(
8542 *
check_neigh_polytope_id);
8546 * ghosted polytope has been found, increment
8555 *
if (handler.polytope_cache.visited_cell_and_faces_id
8556 *
.find({cell_id, f}) ==
8558 *
handler.polytope_cache.visited_cell_and_faces_id))
8560 *
handler.polytope_cache
8561 *
.interface[{current_polytope_id,
8562 *
check_neigh_polytope_id}]
8563 *
.emplace_back(cell, f);
8567 * std::cout <<
"ADDED ("
8568 * << cell->active_cell_index() <<
")
8570 * << current_polytope_id <<
" e "
8571 * << check_neigh_polytope_id <<
8578 *
handler.polytope_cache.visited_cell_and_faces_id
8579 *
.insert({cell_id, f});
8583 *
if (handler.polytope_cache.visited_cell_and_faces_id
8584 *
.find({neighboring_cell_id, nof}) ==
8586 *
handler.polytope_cache.visited_cell_and_faces_id))
8588 *
handler.polytope_cache
8589 *
.interface[{check_neigh_polytope_id,
8590 *
current_polytope_id}]
8591 *
.emplace_back(neighboring_cell, nof);
8593 *
handler.polytope_cache.visited_cell_and_faces_id
8594 *
.insert({neighboring_cell_id, nof});
8598 *
else if (cell->face(f)->at_boundary())
8602 * Boundary face of a boundary cell.
8603 * Note that the neighboring cell must be
invalid.
8609 *
handler.polygon_boundary[master_cell].push_back(
8612 *
if (visited_polygonal_neighbors.find(
8613 *
std::numeric_limits<unsigned int>::max()) ==
8614 *
std::end(visited_polygonal_neighbors))
8618 * boundary face. Notice that `neighboring_cell` is
8622 *
handler.polytope_cache.cell_face_at_boundary[{
8623 *
current_polytope_index,
8624 *
handler.number_of_agglomerated_faces
8625 *
[current_polytope_index]}] = {
true,
8626 *
neighboring_cell};
8628 *
const unsigned int n_face =
8629 *
handler.number_of_agglomerated_faces
8630 *
[current_polytope_index];
8632 *
is_face_at_boundary[n_face] =
true;
8634 *
++handler.number_of_agglomerated_faces
8635 *
[current_polytope_index];
8637 *
visited_polygonal_neighbors.insert(
8638 *
std::numeric_limits<unsigned int>::max());
8643 *
if (handler.polytope_cache.visited_cell_and_faces.find(
8644 *
{cell_index, f}) ==
8645 *
std::end(handler.polytope_cache.visited_cell_and_faces))
8647 *
handler.polytope_cache
8648 *
.interface[{current_polytope_id, current_polytope_id}]
8649 *
.emplace_back(cell, f);
8651 *
handler.polytope_cache.visited_cell_and_faces.insert(
8660 *
if (ghost_counter > 0)
8662 *
const auto parallel_triangulation =
dynamic_cast<
8663 *
const ::parallel::TriangulationBase<dim, spacedim> *
>(
8664 *
&(*handler.tria));
8666 *
const unsigned int n_faces_current_poly =
8667 *
handler.number_of_agglomerated_faces[current_polytope_index];
8671 * Communicate to neighboring ranks that current_polytope_id has
8672 * a number of faces
equal to n_faces_current_poly faces:
8673 * current_polytope_id -> n_faces_current_poly
8676 *
for (
const unsigned int neigh_rank :
8677 *
parallel_triangulation->ghost_owners())
8679 *
handler.local_n_faces[neigh_rank].emplace(current_polytope_id,
8680 *
n_faces_current_poly);
8682 *
handler.local_bdary_info[neigh_rank].emplace(
8683 *
current_polytope_id, is_face_at_boundary);
8685 *
handler.local_ghosted_master_id[neigh_rank].emplace(
8686 *
current_polytope_id, face_to_neigh_id);
8699 *
template class AgglomerationHandler<1>;
8701 *
AgglomerationHandler<1>::create_agglomeration_sparsity_pattern(
8704 *
const bool keep_constrained_dofs,
8708 *
AgglomerationHandler<1>::create_agglomeration_sparsity_pattern(
8711 *
const bool keep_constrained_dofs,
8714 *
template class AgglomerationHandler<2>;
8716 *
AgglomerationHandler<2>::create_agglomeration_sparsity_pattern(
8719 *
const bool keep_constrained_dofs,
8723 *
AgglomerationHandler<2>::create_agglomeration_sparsity_pattern(
8726 *
const bool keep_constrained_dofs,
8729 *
template class AgglomerationHandler<3>;
8731 *
AgglomerationHandler<3>::create_agglomeration_sparsity_pattern(
8734 *
const bool keep_constrained_dofs,
8738 *
AgglomerationHandler<3>::create_agglomeration_sparsity_pattern(
8741 *
const bool keep_constrained_dofs,
8746<a name=
"ann-source/mapping_box.cc"></a>
8747<h1>Annotated version of source/mapping_box.cc</h1>
8764 *
#include <deal.II/base/array_view.h>
8766 *
#include <deal.II/base/qprojector.h>
8767 *
#include <deal.II/base/quadrature.h>
8769 *
#include <deal.II/base/tensor.h>
8771 *
#include <deal.II/dofs/dof_accessor.h>
8773 *
#include <deal.II/fe/fe_values.h>
8775 *
#include <deal.II/grid/tria.h>
8776 *
#include <deal.II/grid/tria_iterator.h>
8778 *
#include <deal.II/lac/full_matrix.h>
8780 *
#include <mapping_box.h>
8782 *
#include <algorithm>
8787 *
ExcCellNotAssociatedWithBox,
8788 *
"You are using MappingBox, but the incoming element is not associated with a"
8789 *
"Bounding Box Cartesian.");
8797 *
template <
typename CellType>
8799 *
has_box(
const CellType &cell,
8800 *
const std::map<types::global_cell_index, types::global_cell_index>
8803 *
Assert((cell->reference_cell().is_hyper_cube() ||
8804 *
cell->reference_cell().is_simplex()),
8805 *
ExcNotImplemented());
8806 *
Assert((translator.find(cell->active_cell_index()) != translator.cend()),
8807 *
ExcCellNotAssociatedWithBox());
8814 *
template <
int dim,
int spacedim>
8815 *
MappingBox<dim, spacedim>::MappingBox(
8817 *
const std::map<types::global_cell_index, types::global_cell_index>
8818 *
&global_to_polytope)
8820 *
Assert(input_boxes.size() > 0,
8821 *
ExcMessage(
"Invalid number of bounding boxes."));
8825 *
copy boxes and map
8828 *
boxes.resize(input_boxes.size());
8829 *
for (
unsigned int i = 0; i < input_boxes.size(); ++i)
8830 *
boxes[i] = input_boxes[i];
8831 *
polytope_translator = global_to_polytope;
8836 *
template <
int dim,
int spacedim>
8837 *
MappingBox<dim, spacedim>::InternalData::InternalData(
const Quadrature<dim> &q)
8841 *
, volume_element(numbers::signaling_nan<double>())
8847 *
template <
int dim,
int spacedim>
8849 *
MappingBox<dim, spacedim>::InternalData::reinit(
const UpdateFlags update_flags,
8854 * store the flags in the
internal data object so we can access them
8855 * in fill_fe_*_values(). use the transitive hull of the required
8859 *
this->update_each = update_flags;
8864 *
template <
int dim,
int spacedim>
8866 *
MappingBox<dim, spacedim>::InternalData::memory_consumption() const
8877 *
template <
int dim,
int spacedim>
8879 *
MappingBox<dim, spacedim>::preserves_vertex_locations() const
8886 *
template <
int dim,
int spacedim>
8888 *
MappingBox<dim, spacedim>::is_compatible_with(
8897 *
ExcMessage(
"The dimension of your mapping (" +
8899 *
") and the reference cell cell_type (" +
8901 *
" ) do not agree."));
8908 *
template <
int dim,
int spacedim>
8910 *
MappingBox<dim, spacedim>::requires_update_flags(
const UpdateFlags in)
const
8914 *
this mapping is pretty simple in that it can basically compute
8915 * every piece of information wanted by
FEValues without requiring
8916 * computing any other quantities. boundary forms are
one exception
8917 * since they can be computed from the normal vectors without much
8930 *
template <
int dim,
int spacedim>
8931 *
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
8932 *
MappingBox<dim, spacedim>::get_data(
const UpdateFlags update_flags,
8935 *
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
8936 *
std::make_unique<InternalData>();
8937 *
data_ptr->reinit(requires_update_flags(update_flags), q);
8944 *
template <
int dim,
int spacedim>
8945 *
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
8946 *
MappingBox<dim, spacedim>::get_subface_data(
8950 *
(
void)update_flags;
8958 *
template <
int dim,
int spacedim>
8960 *
MappingBox<dim, spacedim>::update_cell_extents(
8963 *
const InternalData &
data)
const
8967 * Compute start
point and sizes along axes. The vertices to be looked at
8968 * are 1, 2, 4 compared to the base vertex 0.
8974 *
boxes[polytope_translator.at(cell->active_cell_index())];
8975 *
const std::pair<Point<dim>,
Point<dim>> &bdary_points =
8976 *
current_box.get_boundary_points();
8978 *
for (
unsigned int d = 0;
d < dim; ++
d)
8980 *
const double cell_extent_d = current_box.side_length(d);
8981 *
data.cell_extents[
d] = cell_extent_d;
8983 *
data.traslation[
d] =
8984 *
.5 * (bdary_points.first[
d] +
8985 *
bdary_points.second[
d]);
8987 *
Assert(cell_extent_d != 0.,
8988 *
ExcMessage(
"Cell does not appear to be Cartesian!"));
8989 *
data.inverse_cell_extents[
d] = 1. / cell_extent_d;
8998 *
template <
int dim>
9000 *
transform_quadrature_points(
9003 *
std::vector<
Point<dim>> &quadrature_points)
9006 *
quadrature_points[i] = box.unit_to_real(unit_quadrature_points[i]);
9012 *
template <
int dim,
int spacedim>
9014 *
MappingBox<dim, spacedim>::maybe_update_cell_quadrature_points(
9016 *
const InternalData &
data,
9018 *
std::vector<
Point<dim>> &quadrature_points)
const
9021 *
transform_quadrature_points(
9022 *
boxes[polytope_translator.at(cell->active_cell_index())],
9023 *
unit_quadrature_points,
9024 *
quadrature_points);
9029 *
template <
int dim,
int spacedim>
9031 *
MappingBox<dim, spacedim>::maybe_update_normal_vectors(
9032 *
const unsigned int face_no,
9033 *
const InternalData &
data,
9038 * compute normal vectors. All normals on a face have the same
value.
9044 *
std::fill(normal_vectors.begin(),
9045 *
normal_vectors.end(),
9052 *
template <
int dim,
int spacedim>
9054 *
MappingBox<dim, spacedim>::maybe_update_jacobian_derivatives(
9055 *
const InternalData &
data,
9058 *
&output_data)
const
9063 *
for (
unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i)
9067 *
for (
unsigned int i = 0;
9068 *
i < output_data.jacobian_pushed_forward_grads.size();
9073 *
for (
unsigned int i = 0;
9074 *
i < output_data.jacobian_2nd_derivatives.size();
9076 *
output_data.jacobian_2nd_derivatives[i] =
9080 *
for (
unsigned int i = 0;
9081 *
i < output_data.jacobian_pushed_forward_2nd_derivatives.size();
9083 *
output_data.jacobian_pushed_forward_2nd_derivatives[i] =
9087 *
for (
unsigned int i = 0;
9088 *
i < output_data.jacobian_3rd_derivatives.size();
9090 *
output_data.jacobian_3rd_derivatives[i] =
9094 *
for (
unsigned int i = 0;
9095 *
i < output_data.jacobian_pushed_forward_3rd_derivatives.size();
9097 *
output_data.jacobian_pushed_forward_3rd_derivatives[i] =
9104 *
template <
int dim,
int spacedim>
9106 *
MappingBox<dim, spacedim>::maybe_update_volume_elements(
9107 *
const InternalData &
data)
const
9112 *
for (
unsigned int d = 1;
d < dim; ++
d)
9113 *
volume *=
data.cell_extents[d];
9120 *
template <
int dim,
int spacedim>
9122 *
MappingBox<dim, spacedim>::maybe_update_jacobians(
9123 *
const InternalData &
data,
9126 *
&output_data)
const
9130 *
"compute" Jacobian at the quadrature points, which are all the
9136 *
for (
unsigned int i = 0; i < output_data.jacobians.size(); ++i)
9139 *
for (
unsigned int j = 0; j < dim; ++j)
9140 *
output_data.jacobians[i][j][j] =
data.cell_extents[j];
9146 *
template <
int dim,
int spacedim>
9148 *
MappingBox<dim, spacedim>::maybe_update_inverse_jacobians(
9149 *
const InternalData &
data,
9152 *
&output_data)
const
9156 *
"compute" inverse Jacobian at the quadrature points, which are
9162 *
for (
unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i)
9165 *
for (
unsigned int j = 0; j < dim; ++j)
9166 *
output_data.inverse_jacobians[i][j][j] =
9167 *
data.inverse_cell_extents[j];
9173 *
template <
int dim,
int spacedim>
9175 *
MappingBox<dim, spacedim>::fill_fe_values(
9181 *
&output_data)
const
9183 *
Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox());
9188 * an exception
if that is not possible
9191 *
Assert(
dynamic_cast<const InternalData *
>(&internal_data) !=
nullptr,
9192 *
ExcInternalError());
9193 *
const InternalData &
data =
static_cast<const InternalData &
>(internal_data);
9196 *
update_cell_extents(cell, cell_similarity,
data);
9198 *
maybe_update_cell_quadrature_points(cell,
9200 *
quadrature.get_points(),
9201 *
output_data.quadrature_points);
9206 * product of the local lengths in each coordinate direction
9212 *
double J =
data.cell_extents[0];
9213 *
for (
unsigned int d = 1;
d < dim; ++
d)
9214 *
J *=
data.cell_extents[d];
9215 *
data.volume_element =
J;
9217 *
for (
unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
9218 *
output_data.JxW_values[i] = quadrature.weight(i);
9222 *
maybe_update_jacobians(
data, cell_similarity, output_data);
9223 *
maybe_update_jacobian_derivatives(
data, cell_similarity, output_data);
9224 *
maybe_update_inverse_jacobians(
data, cell_similarity, output_data);
9226 *
return cell_similarity;
9231 *
template <
int dim,
int spacedim>
9233 *
MappingBox<dim, spacedim>::fill_fe_subface_values(
9235 *
const unsigned int face_no,
9236 *
const unsigned int subface_no,
9240 *
&output_data)
const
9246 *
(
void)internal_data;
9247 *
(
void)output_data;
9253 *
template <
int dim,
int spacedim>
9255 *
MappingBox<dim, spacedim>::fill_fe_immersed_surface_values(
9260 *
&output_data)
const
9263 *
Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox());
9268 * exception
if that is not possible.
9271 *
Assert(
dynamic_cast<const InternalData *
>(&internal_data) !=
nullptr,
9272 *
ExcInternalError());
9273 *
const InternalData &
data =
static_cast<const InternalData &
>(internal_data);
9278 *
maybe_update_cell_quadrature_points(cell,
9280 *
quadrature.get_points(),
9281 *
output_data.quadrature_points);
9284 *
for (
unsigned int i = 0; i < output_data.normal_vectors.size(); ++i)
9285 *
output_data.normal_vectors[i] = quadrature.normal_vector(i);
9288 *
for (
unsigned int i = 0; i < output_data.JxW_values.size(); ++i)
9289 *
output_data.JxW_values[i] = quadrature.weight(i);
9291 *
maybe_update_volume_elements(
data);
9299 *
template <
int dim,
int spacedim>
9301 *
MappingBox<dim, spacedim>::transform(
9308 *
Assert(
dynamic_cast<const InternalData *
>(&mapping_data) !=
nullptr,
9309 *
ExcInternalError());
9310 *
const InternalData &
data =
static_cast<const InternalData &
>(mapping_data);
9312 *
switch (mapping_kind)
9318 *
"update_covariant_transformation"));
9320 *
for (
unsigned int i = 0; i < output.size(); ++i)
9321 *
for (
unsigned int d = 0;
d < dim; ++
d)
9322 *
output[i][d] = input[i][d] *
data.inverse_cell_extents[d];
9330 *
"update_contravariant_transformation"));
9332 *
for (
unsigned int i = 0; i < output.size(); ++i)
9333 *
for (
unsigned int d = 0;
d < dim; ++
d)
9334 *
output[i][d] = input[i][d] *
data.cell_extents[d];
9341 *
"update_contravariant_transformation"));
9344 *
"update_volume_elements"));
9346 *
for (
unsigned int i = 0; i < output.size(); ++i)
9347 *
for (
unsigned int d = 0;
d < dim; ++
d)
9349 *
input[i][d] *
data.cell_extents[d] /
data.volume_element;
9359 *
template <
int dim,
int spacedim>
9361 *
MappingBox<dim, spacedim>::transform(
9368 *
Assert(
dynamic_cast<const InternalData *
>(&mapping_data) !=
nullptr,
9369 *
ExcInternalError());
9370 *
const InternalData &
data =
static_cast<const InternalData &
>(mapping_data);
9372 *
switch (mapping_kind)
9378 *
"update_covariant_transformation"));
9380 *
for (
unsigned int i = 0; i < output.size(); ++i)
9381 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9382 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9383 *
output[i][d1][d2] =
9384 *
input[i][d1][d2] *
data.inverse_cell_extents[d2];
9392 *
"update_contravariant_transformation"));
9394 *
for (
unsigned int i = 0; i < output.size(); ++i)
9395 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9396 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9397 *
output[i][d1][d2] = input[i][d1][d2] *
data.cell_extents[d2];
9405 *
"update_covariant_transformation"));
9407 *
for (
unsigned int i = 0; i < output.size(); ++i)
9408 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9409 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9410 *
output[i][d1][d2] = input[i][d1][d2] *
9411 *
data.inverse_cell_extents[d2] *
9412 *
data.inverse_cell_extents[d1];
9420 *
"update_contravariant_transformation"));
9422 *
for (
unsigned int i = 0; i < output.size(); ++i)
9423 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9424 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9425 *
output[i][d1][d2] = input[i][d1][d2] *
data.cell_extents[d2] *
9426 *
data.inverse_cell_extents[d1];
9434 *
"update_contravariant_transformation"));
9437 *
"update_volume_elements"));
9439 *
for (
unsigned int i = 0; i < output.size(); ++i)
9440 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9441 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9442 *
output[i][d1][d2] = input[i][d1][d2] *
data.cell_extents[d2] /
9443 *
data.volume_element;
9451 *
"update_contravariant_transformation"));
9454 *
"update_volume_elements"));
9456 *
for (
unsigned int i = 0; i < output.size(); ++i)
9457 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9458 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9459 *
output[i][d1][d2] = input[i][d1][d2] *
data.cell_extents[d2] *
9460 *
data.inverse_cell_extents[d1] /
9461 *
data.volume_element;
9472 *
template <
int dim,
int spacedim>
9474 *
MappingBox<dim, spacedim>::transform(
9481 *
Assert(
dynamic_cast<const InternalData *
>(&mapping_data) !=
nullptr,
9482 *
ExcInternalError());
9483 *
const InternalData &
data =
static_cast<const InternalData &
>(mapping_data);
9485 *
switch (mapping_kind)
9491 *
"update_covariant_transformation"));
9493 *
for (
unsigned int i = 0; i < output.size(); ++i)
9494 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9495 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9496 *
output[i][d1][d2] =
9497 *
input[i][d1][d2] *
data.inverse_cell_extents[d2];
9505 *
"update_contravariant_transformation"));
9507 *
for (
unsigned int i = 0; i < output.size(); ++i)
9508 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9509 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9510 *
output[i][d1][d2] = input[i][d1][d2] *
data.cell_extents[d2];
9518 *
"update_covariant_transformation"));
9520 *
for (
unsigned int i = 0; i < output.size(); ++i)
9521 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9522 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9523 *
output[i][d1][d2] = input[i][d1][d2] *
9524 *
data.inverse_cell_extents[d2] *
9525 *
data.inverse_cell_extents[d1];
9533 *
"update_contravariant_transformation"));
9535 *
for (
unsigned int i = 0; i < output.size(); ++i)
9536 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9537 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9538 *
output[i][d1][d2] = input[i][d1][d2] *
data.cell_extents[d2] *
9539 *
data.inverse_cell_extents[d1];
9547 *
"update_contravariant_transformation"));
9550 *
"update_volume_elements"));
9552 *
for (
unsigned int i = 0; i < output.size(); ++i)
9553 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9554 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9555 *
output[i][d1][d2] = input[i][d1][d2] *
data.cell_extents[d2] /
9556 *
data.volume_element;
9564 *
"update_contravariant_transformation"));
9567 *
"update_volume_elements"));
9569 *
for (
unsigned int i = 0; i < output.size(); ++i)
9570 *
for (
unsigned int d1 = 0; d1 < dim; ++d1)
9571 *
for (
unsigned int d2 = 0; d2 < dim; ++d2)
9572 *
output[i][d1][d2] = input[i][d1][d2] *
data.cell_extents[d2] *
9573 *
data.inverse_cell_extents[d1] /
9574 *
data.volume_element;
9585 *
template <
int dim,
int spacedim>
9587 *
MappingBox<dim, spacedim>::transform(
9594 *
Assert(
dynamic_cast<const InternalData *
>(&mapping_data) !=
nullptr,
9595 *
ExcInternalError());
9596 *
const InternalData &
data =
static_cast<const InternalData &
>(mapping_data);
9598 *
switch (mapping_kind)
9604 *
"update_covariant_transformation"));
9606 *
for (
unsigned int q = 0; q < output.size(); ++q)
9607 *
for (
unsigned int i = 0; i < spacedim; ++i)
9608 *
for (
unsigned int j = 0; j < spacedim; ++j)
9609 *
for (
unsigned int k = 0; k < spacedim; ++k)
9611 *
output[q][i][j][k] = input[q][i][j][k] *
9612 *
data.inverse_cell_extents[j] *
9613 *
data.inverse_cell_extents[k];
9624 *
template <
int dim,
int spacedim>
9626 *
MappingBox<dim, spacedim>::transform(
9633 *
Assert(
dynamic_cast<const InternalData *
>(&mapping_data) !=
nullptr,
9634 *
ExcInternalError());
9635 *
const InternalData &
data =
static_cast<const InternalData &
>(mapping_data);
9637 *
switch (mapping_kind)
9643 *
"update_covariant_transformation"));
9646 *
"update_contravariant_transformation"));
9648 *
for (
unsigned int q = 0; q < output.size(); ++q)
9649 *
for (
unsigned int i = 0; i < spacedim; ++i)
9650 *
for (
unsigned int j = 0; j < spacedim; ++j)
9651 *
for (
unsigned int k = 0; k < spacedim; ++k)
9653 *
output[q][i][j][k] = input[q][i][j][k] *
9654 *
data.cell_extents[i] *
9655 *
data.inverse_cell_extents[j] *
9656 *
data.inverse_cell_extents[k];
9665 *
"update_covariant_transformation"));
9667 *
for (
unsigned int q = 0; q < output.size(); ++q)
9668 *
for (
unsigned int i = 0; i < spacedim; ++i)
9669 *
for (
unsigned int j = 0; j < spacedim; ++j)
9670 *
for (
unsigned int k = 0; k < spacedim; ++k)
9672 *
output[q][i][j][k] = input[q][i][j][k] *
9673 *
(
data.inverse_cell_extents[i] *
9674 *
data.inverse_cell_extents[j]) *
9675 *
data.inverse_cell_extents[k];
9685 *
"update_covariant_transformation"));
9688 *
"update_contravariant_transformation"));
9691 *
"update_volume_elements"));
9693 *
for (
unsigned int q = 0; q < output.size(); ++q)
9694 *
for (
unsigned int i = 0; i < spacedim; ++i)
9695 *
for (
unsigned int j = 0; j < spacedim; ++j)
9696 *
for (
unsigned int k = 0; k < spacedim; ++k)
9698 *
output[q][i][j][k] =
9699 *
input[q][i][j][k] *
9700 *
(
data.cell_extents[i] /
data.volume_element *
9701 *
data.inverse_cell_extents[j]) *
9702 *
data.inverse_cell_extents[k];
9715 *
template <
int dim,
int spacedim>
9717 *
MappingBox<dim, spacedim>::transform_unit_to_real_cell(
9721 *
Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox());
9722 *
Assert(dim == spacedim, ExcNotImplemented());
9724 *
return boxes[polytope_translator.at(cell->active_cell_index())].unit_to_real(
9730 *
template <
int dim,
int spacedim>
9732 *
MappingBox<dim, spacedim>::transform_real_to_unit_cell(
9736 *
Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox());
9737 *
Assert(dim == spacedim, ExcNotImplemented());
9739 *
return boxes[polytope_translator.at(cell->active_cell_index())].real_to_unit(
9745 *
template <
int dim,
int spacedim>
9747 *
MappingBox<dim, spacedim>::transform_points_real_to_unit_cell(
9752 *
Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox());
9755 *
if (dim != spacedim)
9757 *
for (
unsigned int i = 0; i < real_points.size(); ++i)
9759 *
boxes[polytope_translator.at(cell->active_cell_index())].real_to_unit(
9765 *
template <
int dim,
int spacedim>
9766 *
std::unique_ptr<Mapping<dim, spacedim>>
9767 *
MappingBox<dim, spacedim>::clone() const
9769 *
return std::make_unique<MappingBox<dim, spacedim>>(*this);
9775 * ---------------------------------------------------------------------------
9776 *
explicit instantiations
9779 *
template class MappingBox<1>;
9780 *
template class MappingBox<2>;
9781 *
template class MappingBox<3>;
* * for(const auto &cell :triangulation.active_cell_iterators())
* * int main(int argc, char **argv)
* x_component_mask set(0, true)
* * reference operator*() const
* * * struct InterferenceTaperTransform *
std::ptrdiff_t difference_type
* * Point< dim > operator()(const Point< dim > &p) const *
* * iterator & operator++()
***mech_lbc_system increment_interpolation_handlers push_back(scale_z_handler)
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void attach_triangulation(const Triangulation< dim, spacedim > &)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
virtual Tensor< 1, dim, RangeNumberType > gradient(const Point< dim > &p, const unsigned int component=0) const
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
void attach_triangulation(Triangulation< dim, spacedim > &tria)
Abstract base class for mapping classes.
void initialize(const SparsityPattern &sparsity_pattern)
unsigned int size() const
#define DEAL_II_VERSION_GTE(major, minor, subminor)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define DeclExceptionMsg(Exception, defaulttext)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_jacobian_pushed_forward_2nd_derivatives
@ update_volume_elements
Determinant of the Jacobian.
@ update_contravariant_transformation
Contravariant transformation.
@ update_jacobian_pushed_forward_grads
@ update_jacobian_3rd_derivatives
@ update_values
Shape function values.
@ update_jacobian_grads
Gradient of volume element.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_covariant_transformation
Covariant transformation.
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
@ update_default
No update.
@ update_jacobian_pushed_forward_3rd_derivatives
@ update_boundary_forms
Outer normal vector, not normalized.
@ update_jacobian_2nd_derivatives
@ mapping_covariant_gradient
@ mapping_contravariant_hessian
@ mapping_covariant_hessian
@ mapping_contravariant_gradient
std::vector< index_type > data
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell< dim > &reference_cell)
void simplex(Triangulation< dim, dim > &tria, const std::vector< Point< dim > > &vertices)
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
* * * ScaleZFunction< dim, Number, components >::ScaleZFunction * component(component)
* * if(update_pressure &update_flags) * compute_pressure(constitutive_request
* * * * std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters const
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, ShapeDataMemorySpace > shape_data, const ViewTypeIn in, ViewTypeOut out)
constexpr ReferenceCell< dim > Invalid
std::map< unsigned int, T > some_to_some(const MPI_Comm comm, const std::map< unsigned int, T > &objects_to_send)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void save(Archive &ar, const ::std_cxx26::inplace_vector< T, N > &vec, const unsigned int)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr types::subdomain_id invalid_subdomain_id
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
unsigned int subdomain_id
unsigned short int fe_index
unsigned int global_cell_index
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
SynchronousIterators< Iterators > & operator--(SynchronousIterators< Iterators > &a)
void prev(std::tuple< I1, I2 > &t)
1.9.8