469 *
const std::array<double, 2>
polar =
472 *
constexpr const double alpha = 2. / 3.;
482 * <a name=
"step_75-Parameters"></a>
504 *
std::string type =
"cg_with_amg";
505 *
unsigned int maxiter = 10000;
508 *
unsigned int smoother_sweeps = 1;
509 *
unsigned int n_cycles = 1;
510 *
std::string smoother_type =
"ILU";
515 *
std::string type =
"chebyshev";
516 *
double smoothing_range = 20;
517 *
unsigned int degree = 5;
518 *
unsigned int eig_cg_n_iterations = 20;
525 *
PolynomialCoarseningSequenceType::decrease_by_one;
544 *
unsigned int n_cycles = 8;
569 * <a name=
"step_75-MatrixfreeLaplaceoperator"></a>
570 * <
h3>Matrix-free Laplace
operator</
h3>
590 *
template <
int dim,
typename number>
598 *
LaplaceOperator() =
default;
614 *
number el(
unsigned int,
unsigned int)
const;
616 *
void initialize_dof_vector(VectorType &
vec)
const;
618 *
void vmult(VectorType &dst,
const VectorType &src)
const;
620 *
void Tvmult(VectorType &dst,
const VectorType &src)
const;
631 *
const VectorType &src)
const;
637 *
const VectorType &src,
638 *
const std::pair<unsigned int, unsigned int> &
range)
const;
665 * right-
hand-side vector.
668 *
template <
int dim,
typename number>
669 *
LaplaceOperator<dim, number>::LaplaceOperator(
681 *
template <
int dim,
typename number>
682 *
void LaplaceOperator<dim, number>::reinit(
694 *
this->system_matrix.
clear();
702 *
this->constraints.copy_from(constraints);
714 *
matrix_free.reinit(mapping, dof_handler, constraints, quad,
data);
727 *
dof_handler.locally_owned_dofs(),
739 *
matrix_free.reinit(
742 *
matrix_free.initialize_dof_vector(b);
743 *
matrix_free.initialize_dof_vector(
x);
745 *
constraints.distribute(
x);
747 *
matrix_free.cell_loop(&LaplaceOperator::do_cell_integral_range,
752 *
constraints.set_zero(b);
771 *
template <
int dim,
typename number>
774 *
return matrix_free.get_dof_handler().n_dofs();
785 *
template <
int dim,
typename number>
786 *
number LaplaceOperator<dim, number>::el(
unsigned int,
unsigned int)
const
800 *
template <
int dim,
typename number>
802 *
LaplaceOperator<dim, number>::initialize_dof_vector(VectorType &
vec)
const
804 *
matrix_free.initialize_dof_vector(
vec);
820 *
this->matrix_free.cell_loop(
832 *
template <
int dim,
typename number>
833 *
void LaplaceOperator<dim, number>::Tvmult(VectorType &dst,
834 *
const VectorType &src)
const
836 *
this->vmult(dst, src);
850 *
template <
int dim,
typename number>
851 *
void LaplaceOperator<dim, number>::compute_inverse_diagonal(
852 *
VectorType &diagonal)
const
854 *
this->matrix_free.initialize_dof_vector(diagonal);
857 *
&LaplaceOperator::do_cell_integral_local,
861 *
i = (
std::
abs(i) > 1.0e-10) ? (1.0 / i) : 1.0;
881 *
LaplaceOperator<dim, number>::get_system_matrix()
const
883 *
if (system_matrix.m() == 0 && system_matrix.n() == 0)
885 *
const auto &dof_handler = this->matrix_free.get_dof_handler();
888 *
dof_handler.locally_owned_dofs(),
889 *
dof_handler.get_triangulation().get_mpi_communicator());
894 *
system_matrix.reinit(
dsp);
900 *
&LaplaceOperator::do_cell_integral_local,
904 *
return this->system_matrix;
916 *
template <
int dim,
typename number>
917 *
void LaplaceOperator<dim, number>::do_cell_integral_local(
922 *
for (
const unsigned int q :
integrator.quadrature_point_indices())
935 *
template <
int dim,
typename number>
936 *
void LaplaceOperator<dim, number>::do_cell_integral_global(
939 *
const VectorType &src)
const
943 *
for (
const unsigned int q :
integrator.quadrature_point_indices())
957 *
template <
int dim,
typename number>
958 *
void LaplaceOperator<dim, number>::do_cell_integral_range(
961 *
const VectorType &src,
962 *
const std::pair<unsigned int, unsigned int> &
range)
const
966 *
for (
unsigned cell =
range.first; cell <
range.second; ++cell)
979 * <a name=
"step_75-Solverandpreconditioner"></a>
985 * <a name=
"step_75-Conjugategradientsolverwithmultigridpreconditioner"></a>
1003 *
const VectorType &src,
1014 *
const unsigned int min_level =
mg_matrices.min_level();
1015 *
const unsigned int max_level =
mg_matrices.max_level();
1031 *
min_level, max_level);
1036 *
std::make_shared<SmootherPreconditionerType>();
1042 *
mg_data.smoother.eig_cg_n_iterations;
1056 *
mg_data.coarse_solver.abstol,
1057 *
mg_data.coarse_solver.reltol,
1062 *
std::unique_ptr<MGCoarseGridBase<VectorType>>
mg_coarse;
1068 *
amg_data.smoother_type =
mg_data.coarse_solver.smoother_type.c_str();
1101 * <a name=
"step_75-Hybridpolynomialgeometricglobalcoarseningmultigridpreconditioner"></a>
1112 *
template <
typename VectorType,
typename OperatorType,
int dim>
1116 *
const VectorType &src,
1142 *
std::vector<std::shared_ptr<const Triangulation<dim>>>
1144 *
if (
mg_data.transfer.perform_h_transfer)
1147 *
dof_handler.get_triangulation());
1150 *
&(dof_handler.get_triangulation()), [](
auto *) {});
1161 *
unsigned int max = 0;
1163 *
for (
auto &cell : dof_handler.active_cell_iterators())
1164 *
if (cell->is_locally_owned())
1166 *
std::
max(
max, dof_handler.get_fe(cell->active_fe_index()).degree);
1177 *
for (
unsigned int i = 0; i < dof_handler.get_fe_collection().
size(); ++i)
1179 *
const unsigned int degree = dof_handler.get_fe(i).degree;
1181 *
ExcMessage(
"FECollection does not contain unique degrees."));
1185 *
unsigned int minlevel = 0;
1188 *
dof_handlers.resize(minlevel, maxlevel);
1202 *
dof_handlers[
l].distribute_dofs(dof_handler.get_fe_collection());
1213 *
for (
unsigned int i = 0, l = maxlevel; i <
n_p_levels; ++i, --
l)
1215 *
dof_handlers[
l].reinit(dof_handler.get_triangulation());
1217 *
if (l == maxlevel)
1221 *
auto cell_other = dof_handler.begin_active();
1224 *
if (cell->is_locally_owned())
1225 *
cell->set_active_fe_index(
cell_other->active_fe_index());
1231 *
auto &dof_handler_fine = dof_handlers[
l + 1];
1234 *
auto cell_other = dof_handler_fine.begin_active();
1237 *
if (cell->is_locally_owned())
1246 *
ExcMessage(
"Next polynomial degree in sequence "
1247 *
"does not exist in FECollection."));
1255 *
dof_handlers[
l].distribute_dofs(dof_handler.get_fe_collection());
1267 *
constraints(minlevel, maxlevel);
1271 *
const auto &dof_handler = dof_handlers[
level];
1274 *
constraint.reinit(dof_handler.locally_owned_dofs(),
1287 *
operators[
level] = std::make_unique<OperatorType>(mapping_collection,
1302 *
dof_handlers[
level],
1303 *
constraints[
level + 1],
1304 *
constraints[
level]);
1331 * <a name=
"step_75-ThecodeLaplaceProblemcodeclasstemplate"></a>
1347 *
template <
int dim>
1385 *
std::unique_ptr<FESeries::Legendre<dim>>
legendre;
1386 *
std::unique_ptr<parallel::CellWeights<dim>>
cell_weights;
1400 * <a name=
"step_75-ThecodeLaplaceProblemcodeclassimplementation"></a>
1406 * <a name=
"step_75-Constructor"></a>
1417 *
template <
int dim>
1423 *
,
pcout(std::cout,
1430 *
Assert(prm.min_h_level <= prm.max_h_level,
1432 *
"Triangulation level limits have been incorrectly set up."));
1433 *
Assert(prm.min_p_degree <= prm.max_p_degree,
1434 *
ExcMessage(
"FECollection degrees have been incorrectly set up."));
1456 *
for (
unsigned int degree = 1; degree <= prm.max_p_degree; ++degree)
1458 *
fe_collection.push_back(
FE_Q<dim>(degree));
1476 *
const unsigned int min_fe_index = prm.min_p_degree - 1;
1477 *
fe_collection.set_hierarchy(
1480 *
const unsigned int fe_index) ->
unsigned int {
1481 *
return ((fe_index + 1) < fe_collection.
size()) ? fe_index + 1 :
1486 *
const unsigned int fe_index) ->
unsigned int {
1488 *
ExcMessage(
"Finite element is not part of hierarchy!"));
1498 *
legendre = std::make_unique<FESeries::Legendre<dim>>(
1532 *
cell_weights = std::make_unique<parallel::CellWeights<dim>>(
1535 *
{prm.weighting_factor, prm.weighting_exponent}));
1566 *
prm.max_p_level_difference,
1569 *
boost::signals2::at_front);
1577 * <a name=
"step_75-LaplaceProbleminitialize_grid"></a>
1589 * in
the positive z-direction.
1603 * remove one cell in
every cell
from the negative direction,
but remove one
1623 *
for (
unsigned int d = 0;
d < dim; ++
d)
1643 *
const unsigned int min_fe_index = prm.min_p_degree - 1;
1644 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1645 *
if (cell->is_locally_owned())
1648 *
dof_handler.distribute_dofs(fe_collection);
1658 * <a name=
"step_75-LaplaceProblemsetup_system"></a>
1670 *
template <
int dim>
1675 *
dof_handler.distribute_dofs(fe_collection);
1677 *
locally_owned_dofs = dof_handler.locally_owned_dofs();
1683 *
mpi_communicator);
1684 *
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1686 *
constraints.
clear();
1690 *
mapping_collection, dof_handler, 0,
Solution<dim>(), constraints);
1691 *
constraints.close();
1705 * <a name=
"step_75-LaplaceProblemprint_diagnostics"></a>
1737 *
std::min<unsigned int>(8,
1743 *
pcout <<
" Number of active cells: "
1745 *
<<
" by partition: ";
1755 *
pcout << std::endl;
1759 *
pcout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1761 *
<<
" by partition: ";
1765 *
dof_handler.n_locally_owned_dofs());
1771 *
pcout << std::endl;
1778 *
pcout <<
" Number of constraints: "
1783 *
<<
" by partition: ";
1789 *
pcout << std::endl;
1794 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1795 *
if (cell->is_locally_owned())
1800 *
pcout <<
" Frequencies of poly. degrees:";
1801 *
for (
unsigned int i = 0; i < fe_collection.
size(); ++i)
1804 *
pcout << std::endl;
1813 * <a name=
"step_75-LaplaceProblemsolve_system"></a>
1824 *
template <
int dim>
1833 *
prm.tolerance_factor *
system_rhs.l2_norm());
1840 *
mapping_collection,
1844 *
pcout <<
" Solved in " << solver_control.last_step() <<
" iterations."
1859 * <a name=
"step_75-LaplaceProblemcompute_indicators"></a>
1889 *
template <
int dim>
1897 *
face_quadrature_collection,
1921 * <a name=
"step_75-LaplaceProblemadapt_resolution"></a>
1932 *
template <
int dim>
1953 *
prm.refine_fraction,
1954 *
prm.coarsen_fraction);
1961 * polynomial degree.
1981 *
prm.p_refine_fraction,
1982 *
prm.p_coarsen_fraction);
2006 *
for (
const auto &cell :
2008 *
cell->clear_refine_flag();
2010 *
for (
const auto &cell :
2012 *
cell->clear_coarsen_flag();
2030 *
hp::Refinement::choose_p_over_h(dof_handler);
2075 *
if (cell->is_locally_owned())
2076 *
fe_degrees(cell->active_cell_index()) = cell->get_fe().degree;
2083 *
data_out.attach_dof_handler(dof_handler);
2085 *
data_out.add_data_vector(
fe_degrees,
"fe_degree");
2086 *
data_out.add_data_vector(
subdomain,
"subdomain");
2089 *
data_out.build_patches(mapping_collection);
2091 *
data_out.write_vtu_with_pvtu_record(
2092 *
"./",
"solution", cycle, mpi_communicator, 2, 1);
2100 * <a name=
"step_75-LaplaceProblemrun"></a>
2101 * <
h4>LaplaceProblem::run</
h4>
2114 *
template <
int dim>
2117 *
pcout <<
"Running with Trilinos on "
2119 *
<<
" MPI rank(s)..." << std::endl;
2122 *
pcout <<
"Calculating transformation matrices..." << std::endl;
2124 *
legendre->precalculate_all_transformation_matrices();
2127 *
for (
unsigned int cycle = 0; cycle < prm.n_cycles; ++cycle)
2129 *
pcout <<
"Cycle " << cycle <<
':' << std::endl;
2150 *
pcout << std::endl;
2160 * <a name=
"step_75-main"></a>
2174 *
using namespace dealii;
2175 *
using namespace Step75;
2183 *
catch (std::exception &exc)
2185 *
std::cerr << std::endl
2187 *
<<
"----------------------------------------------------"
2189 *
std::cerr <<
"Exception on processing: " << std::endl
2190 *
<< exc.what() << std::endl
2191 *
<<
"Aborting!" << std::endl
2192 *
<<
"----------------------------------------------------"
2199 *
std::cerr << std::endl
2201 *
<<
"----------------------------------------------------"
2203 *
std::cerr <<
"Unknown exception!" << std::endl
2204 *
<<
"Aborting!" << std::endl
2205 *
<<
"----------------------------------------------------"
2226 Number
of constraints: 542
2232+---------------------------------------------+------------+------------+
2236+---------------------------------+-----------+------------+------------+
2239| initialize grid | 1 | 0.011s | 6.4% |
2240| output results | 1 | 0.0343s | 20% |
2241| setup
system | 1 | 0.00839s | 4.9% |
2242| solve
system | 1 | 0.0896s | 52% |
2243+---------------------------------+-----------+------------+------------+
2251 Number
of constraints: 1202
2257+---------------------------------------------+------------+------------+
2261+---------------------------------+-----------+------------+------------+
2264| output results | 1 | 0.0243s | 15% |
2265| setup
system | 1 | 0.00662s | 4% |
2266| solve
system | 1 | 0.121s | 74% |
2267+---------------------------------+-----------+------------+------------+
2278 Number
of constraints: 14357
2280 Frequencies of poly. degrees: 2:1126 3:1289 4:2725 5:465 6:5
2284+---------------------------------------------+------------+------------+
2288+---------------------------------+-----------+------------+------------+
2291| output results | 1 | 0.0434s | 2.4% |
2292| setup
system | 1 | 0.0165s | 0.9% |
2293| solve
system | 1 | 1.74s | 95% |
2294+---------------------------------+-----------+------------+------------+
2318<
div class=
"twocolumn" style=
"width: 80%; text-align: center;">
2320 <
img src=
"https://www.dealii.org/images/steps/developer/step-75.subdomains-07.svg"
2321 alt=
"Partitioning after seven refinements.">
2324 <
img src=
"https://www.dealii.org/images/steps/developer/step-75.fedegrees-07.svg"
2325 alt=
"Local approximation degrees after seven refinements.">
2331<a name=
"step-75-extensions"></a>
2332<a name=
"step_75-Possibilitiesforextensions"></a><
h3>Possibilities
for extensions</
h3>
2365hp::Refinement::limit_p_level_difference() function.
At this stage, all
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)
static constexpr unsigned int rank
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
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
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#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::vector< index_type > data
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.
constexpr types::blas_int one
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)
constexpr unsigned int invalid_unsigned_int
constexpr types::material_id invalid_material_id
constexpr types::subdomain_id invalid_subdomain_id
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
boost::signals2::signal< void()> post_p4est_refinement