468 *
const std::array<double, 2> polar =
471 *
constexpr const double alpha = 2. / 3.;
481 * <a name=
"step_75-Parameters"></a>
482 * <h3>Parameters</h3>
486 * For
this tutorial, we will use a simplified set of parameters. It is also
488 *
short we decided on
using simple structs. The actual intention of all these
489 * parameters will be described in the upcoming classes at their respective
490 * location where they are used.
494 * The following parameter set controls the coarse-grid solver, the smoothers,
495 * and the inter-grid transfer scheme of the multigrid mechanism.
496 * We populate it with
default parameters.
499 *
struct MultigridParameters
503 * std::string type =
"cg_with_amg";
504 *
unsigned int maxiter = 10000;
505 *
double abstol = 1
e-20;
506 *
double reltol = 1
e-4;
507 *
unsigned int smoother_sweeps = 1;
508 *
unsigned int n_cycles = 1;
509 * std::string smoother_type =
"ILU";
514 * std::string type =
"chebyshev";
515 *
double smoothing_range = 20;
516 *
unsigned int degree = 5;
517 *
unsigned int eig_cg_n_iterations = 20;
524 * PolynomialCoarseningSequenceType::decrease_by_one;
525 *
bool perform_h_transfer =
true;
533 * This is the
general parameter
struct for the problem class. You will find
534 *
this struct divided into several categories, including
general runtime
536 * parameters
for cell weighting. It also contains an instance of the above
537 *
struct for multigrid parameters which will be passed to the multigrid
543 *
unsigned int n_cycles = 8;
544 *
double tolerance_factor = 1
e-12;
546 * MultigridParameters mg_data;
548 *
unsigned int min_h_level = 5;
549 *
unsigned int max_h_level = 12;
550 *
unsigned int min_p_degree = 2;
551 *
unsigned int max_p_degree = 6;
552 *
unsigned int max_p_level_difference = 1;
554 *
double refine_fraction = 0.3;
555 *
double coarsen_fraction = 0.03;
556 *
double p_refine_fraction = 0.9;
557 *
double p_coarsen_fraction = 0.9;
559 *
double weighting_factor = 1.;
560 *
double weighting_exponent = 1.;
568 * <a name=
"step_75-MatrixfreeLaplaceoperator"></a>
569 * <h3>Matrix-free Laplace
operator</h3>
573 * This is a
matrix-free implementation of the Laplace
operator that will
574 * basically take over the part of the `assemble_system()` function from other
575 * tutorials. The meaning of all member
functions will be explained at their
580 * We will use the
FEEvaluation class to evaluate the solution vector
581 * at the quadrature points and to perform the integration. In contrast to
582 * other tutorials, the
template arguments `degree` is set to @f$-1@f$ and
583 * `number of quadrature in 1
d` to @f$0@f$. In
this case,
FEEvaluation selects
584 * dynamically the correct polynomial degree and number of quadrature
585 * points. Here, we introduce an alias to
FEEvaluation with the correct
586 *
template parameters so that we
do not have to worry about them later on.
589 *
template <
int dim,
typename number>
595 *
using FECellIntegrator =
FEEvaluation<dim, -1, 0, 1, number>;
597 * LaplaceOperator() =
default;
603 * VectorType &system_rhs);
609 * VectorType &system_rhs);
613 * number el(
unsigned int,
unsigned int)
const;
615 *
void initialize_dof_vector(VectorType &vec)
const;
617 *
void vmult(VectorType &dst,
const VectorType &src)
const;
619 *
void Tvmult(VectorType &dst,
const VectorType &src)
const;
623 *
void compute_inverse_diagonal(VectorType &diagonal)
const;
626 *
void do_cell_integral_local(FECellIntegrator &integrator)
const;
628 *
void do_cell_integral_global(FECellIntegrator &integrator,
630 *
const VectorType &src)
const;
633 *
void do_cell_integral_range(
636 *
const VectorType &src,
637 *
const std::pair<unsigned int, unsigned int> &range)
const;
643 * To solve the equation system on the coarsest
level with an AMG
644 * preconditioner, we need an actual system
matrix on the coarsest
level.
645 * For
this purpose, we provide a mechanism that optionally computes a
646 *
matrix from the
matrix-free formulation,
for which we introduce a
648 * empty. Once `get_system_matrix()` is called,
this matrix is filled (lazy
649 * allocation). Since
this is a `
const` function, we need the
"mutable"
650 * keyword here. We also need a the constraints
object to build the
matrix.
661 * The following section contains
functions to initialize and reinitialize
663 *
MatrixFree instance. For sake of simplicity, we also compute the system
664 * right-hand-side vector.
667 *
template <
int dim,
typename number>
668 * LaplaceOperator<dim, number>::LaplaceOperator(
673 * VectorType &system_rhs)
675 * this->
reinit(mapping, dof_handler, quad, constraints, system_rhs);
680 *
template <
int dim,
typename number>
681 *
void LaplaceOperator<dim, number>::reinit(
686 * VectorType &system_rhs)
690 * Clear
internal data structures (in the
case that the
operator is reused).
693 * this->system_matrix.clear();
697 * Copy the constraints, since they might be needed
for computation of the
701 * this->constraints.copy_from(constraints);
705 * Set up
MatrixFree. At the quadrature points, we only need to evaluate
713 * matrix_free.reinit(mapping, dof_handler, constraints, quad, data);
717 * Compute the right-hand side vector. For
this purpose, we set up a
second
719 * the constraints due to Dirichlet-boundary conditions. This modified
720 *
operator is applied to a vector with only the Dirichlet
values set. The
721 * result is the
negative right-hand-side vector.
726 * dof_handler.locally_owned_dofs(),
730 * constraints_without_dbc);
731 * constraints_without_dbc.close();
735 * this->initialize_dof_vector(system_rhs);
739 * mapping, dof_handler, constraints_without_dbc, quad, data);
741 * matrix_free.initialize_dof_vector(b);
742 * matrix_free.initialize_dof_vector(x);
744 * constraints.distribute(x);
746 * matrix_free.cell_loop(&LaplaceOperator::do_cell_integral_range,
751 * constraints.set_zero(b);
761 * The following
functions are implicitly needed by the multigrid algorithm,
762 * including the smoothers.
767 * degrees of freedom.
770 *
template <
int dim,
typename number>
773 *
return matrix_free.get_dof_handler().n_dofs();
780 * Access a particular element in the
matrix. This function is neither
781 * needed nor implemented, however, is required to compile the program.
784 *
template <
int dim,
typename number>
785 * number LaplaceOperator<dim, number>::el(
unsigned int,
unsigned int)
const
795 * Initialize the given vector. We simply delegate the task to the
799 *
template <
int dim,
typename number>
801 * LaplaceOperator<dim, number>::initialize_dof_vector(VectorType &vec)
const
810 * Perform an
operator evaluation by looping with the help of
MatrixFree
811 * over all cells and evaluating the effect of the cell integrals (see also:
812 * `do_cell_integral_local()` and `do_cell_integral_global()`).
815 * template <
int dim, typename number>
816 * void LaplaceOperator<dim, number>::vmult(
VectorType &dst,
819 * this->matrix_free.cell_loop(
820 * &LaplaceOperator::do_cell_integral_range, this, dst, src, true);
827 * Perform the transposed
operator evaluation. Since we are considering
828 *
symmetric "matrices",
this function can simply delegate it task to vmult().
831 *
template <
int dim,
typename number>
832 *
void LaplaceOperator<dim, number>::Tvmult(VectorType &dst,
833 *
const VectorType &src)
const
835 * this->vmult(dst, src);
842 * Since we
do not have a system
matrix, we cannot
loop over the the
844 * performing a sequence of
operator evaluations to unit basis vectors.
846 *
namespace is used. The inversion is performed manually afterwards.
849 *
template <
int dim,
typename number>
850 *
void LaplaceOperator<dim, number>::compute_inverse_diagonal(
851 * VectorType &diagonal)
const
853 * this->matrix_free.initialize_dof_vector(diagonal);
856 * &LaplaceOperator::do_cell_integral_local,
860 * i = (
std::
abs(i) > 1.0e-10) ? (1.0 / i) : 1.0;
867 * In the
matrix-free context, no system
matrix is set up during
868 * initialization of
this class. As a consequence, it has to be computed
869 * here
if it should be requested. Since the
matrix is only computed in
870 *
this tutorial
for linear elements (on the coarse grid),
this is
872 * The
matrix entries are obtained via sequence of
operator evaluations.
874 * is used. The
matrix will only be computed
if it has not been set up yet
878 * template <int dim, typename number>
880 * LaplaceOperator<dim, number>::get_system_matrix() const
882 *
if (system_matrix.m() == 0 && system_matrix.n() == 0)
884 *
const auto &dof_handler = this->matrix_free.get_dof_handler();
887 * dof_handler.locally_owned_dofs(),
888 * dof_handler.get_triangulation().get_mpi_communicator());
893 * system_matrix.reinit(dsp);
899 * &LaplaceOperator::do_cell_integral_local,
903 *
return this->system_matrix;
910 * Perform cell integral on a cell batch without gathering and scattering
915 *
template <
int dim,
typename number>
916 *
void LaplaceOperator<dim, number>::do_cell_integral_local(
917 * FECellIntegrator &integrator)
const
921 *
for (
const unsigned int q : integrator.quadrature_point_indices())
922 * integrator.submit_gradient(integrator.get_gradient(q), q);
931 * Same as above but with access to the global vectors.
934 *
template <
int dim,
typename number>
935 *
void LaplaceOperator<dim, number>::do_cell_integral_global(
936 * FECellIntegrator &integrator,
938 *
const VectorType &src)
const
942 *
for (
const unsigned int q : integrator.quadrature_point_indices())
943 * integrator.submit_gradient(integrator.get_gradient(q), q);
952 * This function loops over all cell batches within a cell-batch range and
953 * calls the above function.
956 *
template <
int dim,
typename number>
957 *
void LaplaceOperator<dim, number>::do_cell_integral_range(
960 *
const VectorType &src,
961 *
const std::pair<unsigned int, unsigned int> &range)
const
963 * FECellIntegrator integrator(matrix_free, range);
965 *
for (
unsigned cell = range.first; cell < range.second; ++cell)
967 * integrator.reinit(cell);
969 * do_cell_integral_global(integrator, dst, src);
978 * <a name=
"step_75-Solverandpreconditioner"></a>
979 * <h3>Solver and preconditioner</h3>
984 * <a name=
"step_75-Conjugategradientsolverwithmultigridpreconditioner"></a>
985 * <h4>Conjugate-
gradient solver with multigrid preconditioner</h4>
989 * This function solves the equation system with a sequence of provided
990 * multigrid objects. It is meant to be treated as
general as possible, hence
991 * the multitude of
template parameters.
996 *
typename SystemMatrixType,
997 *
typename LevelMatrixType,
998 *
typename MGTransferType>
1002 *
const VectorType &src,
1003 *
const MultigridParameters &mg_data,
1005 *
const SystemMatrixType &fine_matrix,
1006 *
const MGLevelObject<std::unique_ptr<LevelMatrixType>> &mg_matrices,
1007 *
const MGTransferType &mg_transfer)
1009 *
AssertThrow(mg_data.coarse_solver.type ==
"cg_with_amg",
1010 * ExcNotImplemented());
1011 *
AssertThrow(mg_data.smoother.type ==
"chebyshev", ExcNotImplemented());
1013 *
const unsigned int min_level = mg_matrices.min_level();
1014 *
const unsigned int max_level = mg_matrices.max_level();
1019 * SmootherPreconditionerType>;
1024 * We initialize
level operators and Chebyshev smoothers here.
1030 * min_level, max_level);
1034 * smoother_data[
level].preconditioner =
1035 * std::make_shared<SmootherPreconditionerType>();
1036 * mg_matrices[
level]->compute_inverse_diagonal(
1037 * smoother_data[
level].preconditioner->get_vector());
1038 * smoother_data[
level].smoothing_range = mg_data.smoother.smoothing_range;
1039 * smoother_data[
level].degree = mg_data.smoother.degree;
1040 * smoother_data[
level].eig_cg_n_iterations =
1041 * mg_data.smoother.eig_cg_n_iterations;
1046 * mg_smoother.
initialize(mg_matrices, smoother_data);
1050 * Next, we initialize the coarse-grid solver. We use conjugate-
gradient
1051 * method with AMG as preconditioner.
1054 *
ReductionControl coarse_grid_solver_control(mg_data.coarse_solver.maxiter,
1055 * mg_data.coarse_solver.abstol,
1056 * mg_data.coarse_solver.reltol,
1061 * std::unique_ptr<MGCoarseGridBase<VectorType>> mg_coarse;
1066 * amg_data.n_cycles = mg_data.coarse_solver.n_cycles;
1067 * amg_data.smoother_type = mg_data.coarse_solver.smoother_type.c_str();
1069 * precondition_amg.initialize(mg_matrices[min_level]->get_system_matrix(),
1076 *
decltype(precondition_amg)>>(
1077 * coarse_grid_solver, *mg_matrices[min_level], precondition_amg);
1081 * Finally, we create the
Multigrid object, convert it to a preconditioner,
1082 * and use it
inside of a conjugate-
gradient solver to solve the linear
1083 * system of equations.
1087 * mg_matrix, *mg_coarse, mg_transfer, mg_smoother, mg_smoother);
1089 * PreconditionerType preconditioner(dof,
mg, mg_transfer);
1092 * .
solve(fine_matrix, dst, src, preconditioner);
1100 * <a name=
"step_75-Hybridpolynomialgeometricglobalcoarseningmultigridpreconditioner"></a>
1101 * <h4>Hybrid polynomial/geometric-global-coarsening multigrid preconditioner</h4>
1105 * The above function deals with the actual solution
for a given sequence of
1106 * multigrid objects. This
functions creates the actual multigrid levels, in
1107 * particular the operators, and the transfer
operator as a
1111 *
template <
typename VectorType,
typename OperatorType,
int dim>
1113 *
const OperatorType &system_matrix,
1115 *
const VectorType &src,
1116 *
const MultigridParameters &mg_data,
1124 * as well as, create transfer operators. To be able to
1125 * set up the operators, we need a set of
DoFHandler that we create
1126 * via global coarsening of p or h. For latter, we need also a sequence
1132 * In
case no h-transfer is requested, we provide an empty deleter
for the
1134 * an external field and its destructor is called somewhere
else.
1141 * std::vector<std::shared_ptr<const Triangulation<dim>>>
1142 * coarse_grid_triangulations;
1143 *
if (mg_data.transfer.perform_h_transfer)
1144 * coarse_grid_triangulations =
1146 * dof_handler.get_triangulation());
1148 * coarse_grid_triangulations.emplace_back(
1149 * &(dof_handler.get_triangulation()), [](
auto *) {});
1153 * Determine the total number of levels
for the multigrid operation and
1154 * allocate sufficient memory
for all levels.
1157 *
const unsigned int n_h_levels = coarse_grid_triangulations.size() - 1;
1159 *
const auto get_max_active_fe_degree = [&](
const auto &dof_handler) {
1160 *
unsigned int max = 0;
1162 *
for (
auto &cell : dof_handler.active_cell_iterators())
1163 * if (cell->is_locally_owned())
1165 *
std::
max(
max, dof_handler.get_fe(cell->active_fe_index()).degree);
1170 *
const unsigned int n_p_levels =
1172 * get_max_active_fe_degree(dof_handler), mg_data.transfer.p_sequence)
1175 * std::map<unsigned int, unsigned int> fe_index_for_degree;
1176 *
for (
unsigned int i = 0; i < dof_handler.get_fe_collection().size(); ++i)
1178 *
const unsigned int degree = dof_handler.get_fe(i).degree;
1179 *
Assert(fe_index_for_degree.find(degree) == fe_index_for_degree.end(),
1180 * ExcMessage(
"FECollection does not contain unique degrees."));
1181 * fe_index_for_degree[degree] = i;
1184 *
unsigned int minlevel = 0;
1185 *
unsigned int maxlevel = n_h_levels + n_p_levels - 1;
1187 * dof_handlers.resize(minlevel, maxlevel);
1188 * operators.resize(minlevel, maxlevel);
1189 * transfers.resize(minlevel, maxlevel);
1193 * Loop from the minimum (coarsest) to the maximum (finest)
level and set up
1194 *
DoFHandler accordingly. We start with the h-levels, where we distribute
1195 * on increasingly finer meshes linear elements.
1198 *
for (
unsigned int l = 0;
l < n_h_levels; ++
l)
1200 * dof_handlers[
l].
reinit(*coarse_grid_triangulations[l]);
1201 * dof_handlers[
l].distribute_dofs(dof_handler.get_fe_collection());
1206 * After we reached the finest mesh, we will adjust the polynomial degrees
1207 * on each
level. We reverse iterate over our data structure and start at
1208 * the finest mesh that contains all information about the active FE
1209 * indices. We then lower the polynomial degree of each cell
level by
level.
1212 *
for (
unsigned int i = 0, l = maxlevel; i < n_p_levels; ++i, --
l)
1214 * dof_handlers[
l].reinit(dof_handler.get_triangulation());
1216 *
if (l == maxlevel)
1218 *
auto &dof_handler_mg = dof_handlers[
l];
1220 *
auto cell_other = dof_handler.begin_active();
1221 *
for (
auto &cell : dof_handler_mg.active_cell_iterators())
1223 * if (cell->is_locally_owned())
1224 * cell->set_active_fe_index(cell_other->active_fe_index());
1230 *
auto &dof_handler_fine = dof_handlers[
l + 1];
1231 *
auto &dof_handler_coarse = dof_handlers[
l + 0];
1233 *
auto cell_other = dof_handler_fine.begin_active();
1234 *
for (
auto &cell : dof_handler_coarse.active_cell_iterators())
1236 * if (cell->is_locally_owned())
1238 * const unsigned
int next_degree =
1241 * cell_other->get_fe().degree,
1242 * mg_data.transfer.p_sequence);
1243 *
Assert(fe_index_for_degree.find(next_degree) !=
1244 * fe_index_for_degree.end(),
1245 * ExcMessage(
"Next polynomial degree in sequence "
1246 *
"does not exist in FECollection."));
1248 * cell->set_active_fe_index(fe_index_for_degree[next_degree]);
1254 * dof_handlers[
l].distribute_dofs(dof_handler.get_fe_collection());
1259 * Next, we will create all data structures additionally needed on each
1260 * multigrid
level. This involves determining constraints with homogeneous
1261 * Dirichlet boundary conditions, and building the
operator just like on the
1266 * constraints(minlevel, maxlevel);
1270 *
const auto &dof_handler = dof_handlers[
level];
1271 *
auto &constraint = constraints[
level];
1273 * constraint.reinit(dof_handler.locally_owned_dofs(),
1282 * constraint.close();
1286 * operators[
level] = std::make_unique<OperatorType>(mapping_collection,
1288 * quadrature_collection,
1295 * Set up intergrid operators and collect transfer operators within a single
1296 *
operator as needed by the
Multigrid solver
class.
1301 * dof_handlers[
level],
1302 * constraints[
level + 1],
1303 * constraints[
level]);
1306 * transfers, [&](
const auto l,
auto &vec) {
1307 * operators[
l]->initialize_dof_vector(vec);
1312 * Finally, proceed to solve the problem with multigrid.
1315 * mg_solve(solver_control,
1330 * <a name=
"step_75-ThecodeLaplaceProblemcodeclasstemplate"></a>
1331 * <h3>The <code>LaplaceProblem</code>
class template</h3>
1335 * Now we will
finally declare the main
class of this program, which solves
1336 * the Laplace equation on subsequently refined function spaces. Its structure
1337 * will look familiar as it is similar to the main classes of @ref step_27
"step-27" and
1338 * @ref step_40
"step-40". There are basically just two additions:
1340 * replaced by an
object of the LaplaceOperator
class for the
MatrixFree
1343 * balancing, has been added.
1346 *
template <
int dim>
1347 *
class LaplaceProblem
1350 * LaplaceProblem(
const Parameters ¶meters);
1355 *
void initialize_grid();
1356 *
void setup_system();
1357 *
void print_diagnostics();
1358 *
void solve_system();
1359 *
void compute_indicators();
1360 *
void adapt_resolution();
1361 *
void output_results(
const unsigned int cycle);
1365 *
const Parameters prm;
1380 * LaplaceOperator<dim, double> laplace_operator;
1384 * std::unique_ptr<FESeries::Legendre<dim>>
legendre;
1385 * std::unique_ptr<parallel::CellWeights<dim>> cell_weights;
1399 * <a name=
"step_75-ThecodeLaplaceProblemcodeclassimplementation"></a>
1400 * <h3>The <code>LaplaceProblem</code>
class implementation</h3>
1405 * <a name=
"step_75-Constructor"></a>
1406 * <h4>Constructor</h4>
1410 * The constructor starts with an initializer list that looks similar to the
1411 * one of @ref step_40
"step-40". We again prepare the
ConditionalOStream object to allow
1412 * only the
first process to output anything over the console, and initialize
1413 * the computing timer properly.
1416 *
template <
int dim>
1417 * LaplaceProblem<dim>::LaplaceProblem(
const Parameters ¶meters)
1418 * : mpi_communicator(MPI_COMM_WORLD)
1422 * , pcout(std::cout,
1424 * , computing_timer(mpi_communicator,
1429 *
Assert(prm.min_h_level <= prm.max_h_level,
1431 *
"Triangulation level limits have been incorrectly set up."));
1432 *
Assert(prm.min_p_degree <= prm.max_p_degree,
1433 * ExcMessage(
"FECollection degrees have been incorrectly set up."));
1437 * We need to prepare the data structures
for the
hp-functionality in the
1438 * actual body of the constructor, and create corresponding objects
for
1439 * every degree in the specified range from the parameter
struct. As we are
1440 * only dealing with non-distorted rectangular cells, a linear mapping
1441 *
object is sufficient in
this context.
1445 * In the Parameters
struct, we provide ranges
for levels on which the
1446 * function space is operating with a reasonable resolution. The multigrid
1447 * algorithm
requires linear elements on the coarsest possible
level. So we
1448 * start with the lowest polynomial degree and fill the collection with
1449 * consecutively higher degrees until the user-specified maximum is
1455 *
for (
unsigned int degree = 1; degree <= prm.max_p_degree; ++degree)
1457 * fe_collection.push_back(
FE_Q<dim>(degree));
1458 * quadrature_collection.push_back(
QGauss<dim>(degree + 1));
1464 * As our FECollection contains more finite elements than we want to use
for
1465 * the finite element approximation of our solution, we would like to limit
1466 * the range on which active FE indices can operate on. For
this, the
1467 * FECollection
class allows to register a hierarchy that determines the
1468 * succeeding and preceding finite element in
case of of p-refinement and
1470 * consult
this hierarchy to determine future FE indices. We will
register
1471 * such a hierarchy that only works on finite elements with polynomial
1472 * degrees in the proposed range <code>[min_p_degree, max_p_degree]</code>.
1475 *
const unsigned int min_fe_index = prm.min_p_degree - 1;
1476 * fe_collection.set_hierarchy(
1479 *
const unsigned int fe_index) ->
unsigned int {
1480 *
return ((fe_index + 1) < fe_collection.
size()) ? fe_index + 1 :
1485 *
const unsigned int fe_index) ->
unsigned int {
1486 *
Assert(fe_index >= min_fe_index,
1487 * ExcMessage(
"Finite element is not part of hierarchy!"));
1494 *
for smoothness estimation.
1497 *
legendre = std::make_unique<FESeries::Legendre<dim>>(
1502 * The next part is going to be tricky. During execution of refinement, a
1503 * few
hp-algorithms need to interfere with the actual refinement process on
1506 * the actual refinement process and trigger all connected
functions. We
1507 * require
this functionality
for load balancing and to limit the polynomial
1508 * degrees of neighboring cells.
1512 * For the former, we would like to assign a weight to every cell that is
1513 * proportional to the number of degrees of freedom of its future finite
1515 * easily attach individual weights at the right place during the refinement
1516 * process, i.e., after all
refine and
coarsen flags have been set correctly
1517 *
for hp-adaptation and right before repartitioning
for load balancing is
1518 * about to happen.
Functions can be registered that will attach weights in
1519 * the form that @f$a (n_\text{dofs})^b@f$ with a provided pair of parameters
1520 * @f$(a,b)@f$. We
register such a function in the following.
1524 * For load balancing, efficient solvers like the one we use should
scale
1525 * linearly with the number of degrees of freedom owned. We set the
1526 * parameters
for cell weighting correspondingly: A weighting factor of @f$1@f$
1527 * and an exponent of @f$1@f$ (see the definitions of the `weighting_factor` and
1528 * `weighting_exponent` above).
1531 * cell_weights = std::make_unique<parallel::CellWeights<dim>>(
1534 * {prm.weighting_factor, prm.weighting_exponent}));
1538 * In h-adaptive applications, we ensure a 2:1 mesh balance by limiting the
1539 * difference of refinement levels of neighboring cells to one. With the
1540 *
second call in the following code snippet, we will ensure the same
for
1541 * p-levels on neighboring cells: levels of future finite elements are not
1542 * allowed to differ by more than a specified difference. The function
1544 * be connected to a very specific signal in the
parallel context. The issue
1545 * is that we need to know how the mesh will be actually refined to set
1546 * future FE indices accordingly. As we ask the p4est oracle to perform
1547 * refinement, we need to ensure that the
Triangulation has been updated
1548 * with the adaptation flags of the oracle
first. An instantiation of
1550 * that
for the duration of its life. Thus, we will create an
object of
this
1551 *
class right before limiting the p-
level difference, and connect the
1552 * corresponding
lambda function to the signal
1554 * after the oracle got refined, but before the
Triangulation is refined.
1555 * Furthermore, we specify that
this function will be connected to the front
1556 * of the signal, to ensure that the modification is performed before any
1557 * other function connected to the same signal.
1561 * [&, min_fe_index]() {
1565 * prm.max_p_level_difference,
1568 * boost::signals2::at_front);
1576 * <a name=
"step_75-LaplaceProbleminitialize_grid"></a>
1577 * <h4>LaplaceProblem::initialize_grid</h4>
1582 * as demonstrated in @ref step_50
"step-50". However in the 2
d case, that particular
1583 * function removes the
first quadrant,
while we need the fourth quadrant
1584 * removed in our scenario. Thus, we will use a different function
1586 * the mesh. Furthermore, we formulate that function in a way that it also
1587 * generates a 3d mesh: the 2d L-shaped domain will basically elongated by 1
1588 * in the positive z-direction.
1593 * The parameters that we need to provide are
Point objects for the lower left
1594 * and top right corners, as well as the number of repetitions that the base
1595 * mesh will have in each direction. We provide them for the
first two
1596 * dimensions and treat the higher third dimension separately.
1600 * To create a L-shaped domain, we need to remove the excess cells. For this,
1601 * we specify the <code>cells_to_remove</code> accordingly. We would like to
1602 * remove one cell in every cell from the negative direction, but remove one
1603 * from the positive x-direction.
1607 * On the coarse grid, we set the initial active FE indices and distribute the
1608 * degrees of freedom once. We do that in order to assign the
hp::FECollection
1609 * to the
DoFHandler, so that all cells know how many DoFs they are going to
1610 * have. This step is mandatory for the weighted load balancing algorithm,
1611 * which will be called implicitly in
1615 * template <
int dim>
1616 *
void LaplaceProblem<dim>::initialize_grid()
1620 * std::vector<unsigned int> repetitions(dim);
1622 *
for (
unsigned int d = 0;
d < dim; ++
d)
1625 * repetitions[
d] = 2;
1626 * bottom_left[
d] = -1.;
1627 * top_right[
d] = 1.;
1631 * repetitions[
d] = 1;
1632 * bottom_left[
d] = 0.;
1633 * top_right[
d] = 1.;
1636 * std::vector<int> cells_to_remove(dim, 1);
1637 * cells_to_remove[0] = -1;
1640 *
triangulation, repetitions, bottom_left, top_right, cells_to_remove);
1642 *
const unsigned int min_fe_index = prm.min_p_degree - 1;
1643 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1644 * if (cell->is_locally_owned())
1645 * cell->set_active_fe_index(min_fe_index);
1647 * dof_handler.distribute_dofs(fe_collection);
1657 * <a name=
"step_75-LaplaceProblemsetup_system"></a>
1658 * <h4>LaplaceProblem::setup_system</h4>
1662 * This function looks exactly the same to the one of @ref step_40
"step-40", but you will
1663 * notice the absence of the system
matrix as well as the scaffold that
1664 * surrounds it. Instead, we will initialize the
MatrixFree formulation of the
1665 * <code>laplace_operator</code> here. For boundary conditions, we will use
1666 * the Solution
class introduced earlier in this tutorial.
1669 *
template <
int dim>
1670 *
void LaplaceProblem<dim>::setup_system()
1674 * dof_handler.distribute_dofs(fe_collection);
1676 * locally_owned_dofs = dof_handler.locally_owned_dofs();
1677 * locally_relevant_dofs =
1680 * locally_relevant_solution.reinit(locally_owned_dofs,
1681 * locally_relevant_dofs,
1682 * mpi_communicator);
1683 * system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1685 * constraints.clear();
1686 * constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
1689 * mapping_collection, dof_handler, 0, Solution<dim>(), constraints);
1690 * constraints.close();
1692 * laplace_operator.reinit(mapping_collection,
1694 * quadrature_collection,
1704 * <a name=
"step_75-LaplaceProblemprint_diagnostics"></a>
1705 * <h4>LaplaceProblem::print_diagnostics</h4>
1709 * This is a function that prints additional diagnostics about the equation
1710 * system and its partitioning. In addition to the usual global number of
1711 * active cells and degrees of freedom, we also output their local
1712 * equivalents. For a regulated output, we will communicate the local
1714 * which will then output all information. Output of local quantities is
1715 * limited to the
first 8 processes to avoid cluttering the terminal.
1719 * Furthermore, we would like to print the frequencies of the polynomial
1720 * degrees in the numerical discretization. Since
this information is only
1721 * stored locally, we will count the finite elements on locally owned cells
1725 *
template <
int dim>
1726 *
void LaplaceProblem<dim>::print_diagnostics()
1728 *
const unsigned int first_n_processes =
1729 * std::min<unsigned int>(8,
1731 *
const bool output_cropped =
1735 * pcout <<
" Number of active cells: "
1737 * <<
" by partition: ";
1739 * std::vector<unsigned int> n_active_cells_per_subdomain =
1742 *
for (
unsigned int i = 0; i < first_n_processes; ++i)
1743 * pcout <<
' ' << n_active_cells_per_subdomain[i];
1744 *
if (output_cropped)
1746 * pcout << std::endl;
1750 * pcout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1752 * <<
" by partition: ";
1754 * std::vector<types::global_dof_index> n_dofs_per_subdomain =
1756 * dof_handler.n_locally_owned_dofs());
1757 *
for (
unsigned int i = 0; i < first_n_processes; ++i)
1758 * pcout <<
' ' << n_dofs_per_subdomain[i];
1759 *
if (output_cropped)
1761 * pcout << std::endl;
1765 * std::vector<types::global_dof_index> n_constraints_per_subdomain =
1768 * pcout <<
" Number of constraints: "
1769 * << std::accumulate(n_constraints_per_subdomain.begin(),
1770 * n_constraints_per_subdomain.end(),
1773 * <<
" by partition: ";
1774 *
for (
unsigned int i = 0; i < first_n_processes; ++i)
1775 * pcout <<
' ' << n_constraints_per_subdomain[i];
1776 *
if (output_cropped)
1778 * pcout << std::endl;
1782 * std::vector<unsigned int> n_fe_indices(fe_collection.
size(), 0);
1783 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1784 * if (cell->is_locally_owned())
1785 * n_fe_indices[cell->active_fe_index()]++;
1789 * pcout <<
" Frequencies of poly. degrees:";
1790 *
for (
unsigned int i = 0; i < fe_collection.
size(); ++i)
1791 *
if (n_fe_indices[i] > 0)
1792 * pcout <<
' ' << fe_collection[i].degree <<
':' << n_fe_indices[i];
1793 * pcout << std::endl;
1802 * <a name=
"step_75-LaplaceProblemsolve_system"></a>
1803 * <h4>LaplaceProblem::solve_system</h4>
1807 * The scaffold around the solution is similar to the one of @ref step_40
"step-40". We
1808 * prepare a vector that matches the requirements of
MatrixFree and collect
1809 * the locally-relevant degrees of freedoms we solved the equation system. The
1810 * solution happens with the function introduced earlier.
1813 *
template <
int dim>
1814 *
void LaplaceProblem<dim>::solve_system()
1819 * laplace_operator.initialize_dof_vector(completely_distributed_solution);
1822 * prm.tolerance_factor * system_rhs.l2_norm());
1824 * solve_with_gmg(solver_control,
1826 * completely_distributed_solution,
1829 * mapping_collection,
1831 * quadrature_collection);
1833 * pcout <<
" Solved in " << solver_control.last_step() <<
" iterations."
1836 * constraints.distribute(completely_distributed_solution);
1838 * locally_relevant_solution.copy_locally_owned_data_from(
1839 * completely_distributed_solution);
1840 * locally_relevant_solution.update_ghost_values();
1848 * <a name=
"step_75-LaplaceProblemcompute_indicators"></a>
1849 * <h4>LaplaceProblem::compute_indicators</h4>
1853 * This function contains only a part of the typical <code>refine_grid</code>
1854 * function from other tutorials and is
new in that sense. Here, we will only
1855 * calculate all indicators
for adaptation with actually refining the grid. We
1856 *
do this for the purpose of writing all indicators to the file system, so we
1857 * store them
for later.
1861 * Since we are dealing the an elliptic problem, we will make use of the
1863 * scaling factor of the underlying face integrals to be dependent on the
1864 * actual polynomial degree of the neighboring elements is favorable in
1865 *
hp-adaptive applications @cite davydov2017hp. We can
do this by specifying
1866 * the very last parameter from the additional ones you notices. The others
1867 * are actually just the defaults.
1871 * For the purpose of
hp-adaptation, we will calculate smoothness estimates
1872 * with the strategy presented in the tutorial introduction and use the
1874 * we set the minimal polynomial degree to 2 as it seems that the smoothness
1875 * estimation algorithms have trouble with linear elements.
1878 *
template <
int dim>
1879 *
void LaplaceProblem<dim>::compute_indicators()
1883 * estimated_error_per_cell.grow_or_shrink(
triangulation.n_active_cells());
1886 * face_quadrature_collection,
1888 * locally_relevant_solution,
1889 * estimated_error_per_cell,
1898 * hp_decision_indicators.grow_or_shrink(
triangulation.n_active_cells());
1901 * locally_relevant_solution,
1902 * hp_decision_indicators);
1910 * <a name=
"step_75-LaplaceProblemadapt_resolution"></a>
1911 * <h4>LaplaceProblem::adapt_resolution</h4>
1915 * With the previously calculated indicators, we will
finally flag all cells
1916 *
for adaptation and also execute refinement in
this function. As in previous
1917 * tutorials, we will use the
"fixed number" strategy, but now
for
1921 *
template <
int dim>
1922 *
void LaplaceProblem<dim>::adapt_resolution()
1928 * First, we will set
refine and
coarsen flags based on the error estimates
1929 * on each cell. There is
nothing new here.
1934 * elaborated in the other deal.II tutorials:
using the fixed number
1935 * strategy, we will flag 30% of all cells
for refinement and 3%
for
1936 * coarsening, as provided in the Parameters
struct.
1941 * estimated_error_per_cell,
1942 * prm.refine_fraction,
1943 * prm.coarsen_fraction);
1947 * Next, we will make all adjustments
for hp-adaptation. We want to
refine
1948 * and
coarsen those cells flagged in the previous step, but need to decide
1949 *
if we would like to
do it by adjusting the grid resolution or the
1950 * polynomial degree.
1954 * The next function call sets future FE indices according to the previously
1955 * calculated smoothness indicators as p-adaptation indicators. These
1956 * indices will only be set on those cells that have
refine or
coarsen flags
1961 * For the p-adaptation fractions, we will take an educated guess. Since we
1962 * only expect a single singularity in our scenario, i.e., in the origin of
1963 * the domain, and a smooth solution anywhere
else, we would like to
1964 * strongly prefer to use p-adaptation over h-adaptation. This reflects in
1965 * our choice of a fraction of 90%
for both p-refinement and p-coarsening.
1969 * hp_decision_indicators,
1970 * prm.p_refine_fraction,
1971 * prm.p_coarsen_fraction);
1975 * After setting all indicators, we will remove those that exceed the
1976 * specified limits of the provided
level ranges in the Parameters
struct.
1977 * This limitation naturally arises
for p-adaptation as the number of
1978 * supplied finite elements is limited. In addition, we registered a custom
1979 * hierarchy
for p-adaptation in the constructor. Now, we need to
do this
1980 * manually in the h-adaptive context like in @ref step_31
"step-31".
1984 * We will iterate over all cells on the designated
min and
max levels and
1985 * remove the corresponding flags. As an alternative, we could also flag
1986 * these cells
for p-adaptation by setting future FE indices accordingly
1992 * ExcInternalError());
1995 *
for (
const auto &cell :
1996 *
triangulation.active_cell_iterators_on_level(prm.max_h_level))
1997 * cell->clear_refine_flag();
1999 *
for (
const auto &cell :
2000 *
triangulation.active_cell_iterators_on_level(prm.min_h_level))
2001 * cell->clear_coarsen_flag();
2005 * At
this stage, we have both the future FE indices and the classic
refine
2006 * and
coarsen flags set. The latter will be interpreted by
2008 * our previous modification ensures that the resulting
Triangulation stays
2009 * within the specified
level range.
2013 * Now, we would like to only impose one type of adaptation on cells, which
2014 * is what the next function will sort out for us. In
short, on cells which
2015 * have both
types of indicators assigned, we will favor the p-adaptation
2016 * one and remove the h-adaptation one.
2019 *
hp::Refinement::choose_p_over_h(dof_handler);
2023 * In the end, we are left to execute coarsening and refinement. Here, not
2024 * only the grid will be updated, but also all previous future FE indices
2025 * will become active.
2029 * Remember that we have attached functions to
triangulation signals in the
2030 * constructor, will be triggered in this function call. So there is even
2031 * more happening: weighted repartitioning will be performed to ensure load
2032 * balancing, as well as we will limit the difference of p-levels between
2033 * neighboring cells.
2044 * <a name="step_75-LaplaceProblemoutput_results"></a>
2045 * <h4>LaplaceProblem::output_results</h4>
2049 * Writing results to the file system in
parallel applications works exactly
2050 * like in @ref step_40 "step-40". In addition to the data containers that we prepared
2051 * throughout the tutorial, we would also like to write out the polynomial
2052 * degree of each finite element on the grid as well as the subdomain each
2053 * cell belongs to. We prepare necessary containers for this in the scope of
2057 * template <
int dim>
2058 *
void LaplaceProblem<dim>::output_results(const
unsigned int cycle)
2064 * if (cell->is_locally_owned())
2065 * fe_degrees(cell->active_cell_index()) = cell->get_fe().degree;
2068 *
for (
auto &subd : subdomain)
2073 * data_out.add_data_vector(locally_relevant_solution,
"solution");
2074 * data_out.add_data_vector(fe_degrees,
"fe_degree");
2075 * data_out.add_data_vector(subdomain,
"subdomain");
2076 * data_out.add_data_vector(estimated_error_per_cell,
"error");
2077 * data_out.add_data_vector(hp_decision_indicators,
"hp_indicator");
2078 * data_out.build_patches(mapping_collection);
2080 * data_out.write_vtu_with_pvtu_record(
2081 *
"./",
"solution", cycle, mpi_communicator, 2, 1);
2089 * <a name=
"step_75-LaplaceProblemrun"></a>
2090 * <h4>LaplaceProblem::run</h4>
2094 * The actual
run function again looks very familiar to @ref step_40
"step-40". The only
2095 * addition is the bracketed section that precedes the actual cycle
loop.
2096 * Here, we will pre-calculate the Legendre transformation matrices. In
2097 *
general, these will be calculated on the fly via lazy allocation whenever a
2098 * certain
matrix is needed. For timing purposes however, we would like to
2099 * calculate them all at once before the actual time measurement begins. We
2100 * will thus designate their calculation to their own scope.
2103 *
template <
int dim>
2104 *
void LaplaceProblem<dim>::run()
2106 * pcout <<
"Running with Trilinos on "
2108 * <<
" MPI rank(s)..." << std::endl;
2111 * pcout <<
"Calculating transformation matrices..." << std::endl;
2113 *
legendre->precalculate_all_transformation_matrices();
2116 *
for (
unsigned int cycle = 0; cycle < prm.n_cycles; ++cycle)
2118 * pcout <<
"Cycle " << cycle <<
':' << std::endl;
2121 * initialize_grid();
2123 * adapt_resolution();
2127 * print_diagnostics();
2131 * compute_indicators();
2134 * output_results(cycle);
2136 * computing_timer.print_summary();
2137 * computing_timer.reset();
2139 * pcout << std::endl;
2149 * <a name=
"step_75-main"></a>
2154 * The
final function is the <code>main</code> function that will ultimately
2155 * create and
run a LaplaceOperator instantiation. Its structure is similar to
2156 * most other tutorial programs.
2159 *
int main(
int argc,
char *argv[])
2163 *
using namespace dealii;
2164 *
using namespace Step75;
2169 * LaplaceProblem<2> laplace_problem(prm);
2170 * laplace_problem.run();
2172 *
catch (std::exception &exc)
2174 * std::cerr << std::endl
2176 * <<
"----------------------------------------------------"
2178 * std::cerr <<
"Exception on processing: " << std::endl
2179 * << exc.what() << std::endl
2180 * <<
"Aborting!" << std::endl
2181 * <<
"----------------------------------------------------"
2188 * std::cerr << std::endl
2190 * <<
"----------------------------------------------------"
2192 * std::cerr <<
"Unknown exception!" << std::endl
2193 * <<
"Aborting!" << std::endl
2194 * <<
"----------------------------------------------------"
2202<a name=
"step_75-Results"></a><h1>Results</h1>
2205When you
run the program with the given parameters on four processes in
2206release mode, your terminal output should look like
this:
2208Running with Trilinos on 4
MPI rank(s)...
2209Calculating transformation matrices...
2211 Number of active cells: 3072
2213 Number of degrees of freedom: 12545
2215 Number of constraints: 542
2217 Frequencies of poly. degrees: 2:3072
2218 Solved in 7 iterations.
2221+---------------------------------------------+------------+------------+
2222| Total wallclock time elapsed since start | 0.172s | |
2224| Section | no. calls | wall time | % of total |
2225+---------------------------------+-----------+------------+------------+
2226| calculate transformation | 1 | 0.0194s | 11% |
2227| compute indicators | 1 | 0.00676s | 3.9% |
2228| initialize grid | 1 | 0.011s | 6.4% |
2229| output results | 1 | 0.0343s | 20% |
2230| setup system | 1 | 0.00839s | 4.9% |
2231| solve system | 1 | 0.0896s | 52% |
2232+---------------------------------+-----------+------------+------------+
2236 Number of active cells: 3351
2238 Number of degrees of freedom: 18228
2240 Number of constraints: 1202
2242 Frequencies of poly. degrees: 2:2522 3:829
2243 Solved in 7 iterations.
2246+---------------------------------------------+------------+------------+
2247| Total wallclock time elapsed since start | 0.165s | |
2249| Section | no. calls | wall time | % of total |
2250+---------------------------------+-----------+------------+------------+
2251| adapt resolution | 1 | 0.00473s | 2.9% |
2252| compute indicators | 1 | 0.00764s | 4.6% |
2253| output results | 1 | 0.0243s | 15% |
2254| setup system | 1 | 0.00662s | 4% |
2255| solve system | 1 | 0.121s | 74% |
2256+---------------------------------+-----------+------------+------------+
2263 Number of active cells: 5610
2265 Number of degrees of freedom: 82047
2267 Number of constraints: 14357
2269 Frequencies of poly. degrees: 2:1126 3:1289 4:2725 5:465 6:5
2270 Solved in 7 iterations.
2273+---------------------------------------------+------------+------------+
2274| Total wallclock time elapsed since start | 1.83s | |
2276| Section | no. calls | wall time | % of total |
2277+---------------------------------+-----------+------------+------------+
2278| adapt resolution | 1 | 0.00834s | 0.46% |
2279| compute indicators | 1 | 0.0178s | 0.97% |
2280| output results | 1 | 0.0434s | 2.4% |
2281| setup system | 1 | 0.0165s | 0.9% |
2282| solve system | 1 | 1.74s | 95% |
2283+---------------------------------+-----------+------------+------------+
2286When running the code with more processes, you will notice slight
2287differences in the number of active cells and degrees of freedom. This
2288is due to the fact that solver and preconditioner depend on the
2289partitioning of the problem, which might yield to slight differences of
2290the solution in the last digits and ultimately yields to different
2293Furthermore, the number of iterations
for the solver stays about the
2294same in all cycles despite
hp-adaptation, indicating the robustness of
2295the proposed algorithms and promising good scalability
for even larger
2296problem sizes and on more processes.
2298Let us have a look at the graphical output of the program. After all
2299refinement cycles in the given parameter configuration, the actual
2300discretized function space looks like the following with its
2301partitioning on twelve processes on the left and the polynomial degrees
2302of finite elements on the right. In the left picture, each color
2303represents a unique subdomain. In the right picture, the lightest color
2304corresponds to the polynomial degree two and the darkest one corresponds
2307<div
class=
"twocolumn" style=
"width: 80%; text-align: center;">
2309 <img src=
"https://www.dealii.org/images/steps/developer/step-75.subdomains-07.svg"
2310 alt=
"Partitioning after seven refinements.">
2313 <img src=
"https://www.dealii.org/images/steps/developer/step-75.fedegrees-07.svg"
2314 alt=
"Local approximation degrees after seven refinements.">
2320<a name=
"step-75-extensions"></a>
2321<a name=
"step_75-Possibilitiesforextensions"></a><h3>Possibilities
for extensions</h3>
2324This tutorial shows only one particular way how to use
parallel
2325hp-adaptive finite element methods. In the following paragraphs, you
2326will get to know which alternatives are possible. Most of these
2327extensions are already part of https:
2328which provides you with implementation examples that you can play
2332<a name=
"step_75-Differenthpdecisionstrategies"></a><h4>Different
hp-decision strategies</h4>
2335The deal.II library offers multiple strategies to decide which type of
2336adaptation to impose on cells: either adjust the grid resolution or
2337change the polynomial degree. We only presented the <i>Legendre
2338coefficient decay</i> strategy in
this tutorial,
while @ref step_27
"step-27"
2339demonstrated the <i>Fourier</i> equivalent of the same idea.
2341See the
"possibilities for extensions" section of @ref step_27
"step-27" for an
2342overview over these strategies, or the corresponding documentation
2343for a detailed description.
2345There, another strategy is mentioned that has not been shown in any
2346tutorial so far: the strategy based on <i>refinement history</i>. The
2347usage of
this method
for parallel distributed applications is more
2348tricky than the others, so we will highlight the challenges that come
2349along with it. We need information about the
final state of refinement
2350flags, and we need to transfer the solution across refined meshes. For
2352function to the
Triangulation::Signals::post_p4est_refinement signal in
2353a way that it will be called <i>after</i> the
2354hp::Refinement::limit_p_level_difference() function. At this stage, all
2355refinement flags and future FE indices are terminally set and a reliable
2356prediction of the error is possible. The predicted error then needs to
2357be transferred across refined meshes with the aid of
2360Try implementing one of these strategies into this tutorial and observe
2361the subtle changes to the results. You will notice that all strategies
2362are capable of identifying the singularities near the reentrant corners
2363and will perform @f$h@f$-refinement in these regions, while preferring
2364@f$p@f$-refinement in the bulk domain. A detailed comparison of these
2365strategies is presented in @cite fehling2020 .
2368<a name="step_75-Solvewithmatrixbasedmethods"></a><h4>Solve with matrix-based methods</h4>
2371This tutorial focuses solely on matrix-free strategies. All
hp-adaptive
2372algorithms however also work with matrix-based approaches in the
2375To create a system matrix, you can either use the
2376LaplaceOperator::get_system_matrix() function, or use an
2377<code>assemble_system()</code> function similar to the one of @ref step_27 "step-27".
2378You can then pass the system matrix to the solver as usual.
2380You can time the results of both matrix-based and matrix-free
2381implementations, quantify the speed-up, and convince yourself which
2385<a name="step_75-Multigridvariants"></a><h4>
Multigrid variants</h4>
2388For sake of simplicity, we have restricted ourselves to a single type of
2389coarse-grid solver (CG with AMG), smoother (Chebyshev smoother with
2390point Jacobi preconditioner), and geometric-coarsening scheme (global
2391coarsening) within the multigrid algorithm. Feel free to try out
2392alternatives and investigate their performance and robustness.
2395<a name="step_75-PlainProg"></a>
2396<h1> The plain program</h1>
2397@include "step-75.cc"
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void reinit(const Triangulation< dim, spacedim > &tria)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename PreconditionerType::AdditionalData &additional_data=typename PreconditionerType::AdditionalData())
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
virtual types::subdomain_id locally_owned_subdomain() const
void coarsen_global(const unsigned int times=1)
virtual void execute_coarsening_and_refinement()
unsigned int size() const
static WeightingFunction ndofs_weighting(const std::pair< float, float > &coefficients)
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define Assert(cond, exc)
#define AssertThrow(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
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_gradients
Shape function gradients.
#define DEAL_II_NOT_IMPLEMENTED()
std::array< double, dim > to_spherical(const Point< dim > &point)
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
void subdivided_hyper_L(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &bottom_left, const Point< dim > &top_right, const std::vector< int > &n_cells_to_remove)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
@ general
No special properties.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
SymmetricTensor< 2, dim, Number > e(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 > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
FESeries::Legendre< dim, spacedim > default_fe_series(const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)
void coefficient_decay(FESeries::Legendre< dim, spacedim > &fe_legendre, const DoFHandler< dim, spacedim > &dof_handler, const VectorType &solution, Vector< float > &smoothness_indicators, const VectorTools::NormType regression_strategy=VectorTools::Linfty_norm, const double smallest_abs_coefficient=1e-10, const bool only_flagged_cells=false)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::vector< T > gather(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
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)
bool limit_p_level_difference(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int max_difference=1, const unsigned int contains_fe_index=0)
void predict_error(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &error_indicators, Vector< Number > &predicted_errors, const double gamma_p=std::sqrt(0.4), const double gamma_h=2., const double gamma_n=1.)
void p_adaptivity_fixed_number(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const double p_refine_fraction=0.5, const double p_coarsen_fraction=0.5, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_coarsen=std::less_equal< Number >())
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const types::material_id invalid_material_id
const types::subdomain_id invalid_subdomain_id
static const unsigned int invalid_unsigned_int
void refine_and_coarsen_fixed_number(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const types::global_cell_index max_n_cells=std::numeric_limits< types::global_cell_index >::max())
double legendre(unsigned int l, double x)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
UpdateFlags mapping_update_flags
boost::signals2::signal< void()> post_p4est_refinement
unsigned int smoother_sweeps