943 *
for (
unsigned int d = 0;
d < dim; ++
d)
953 * <a name=
"step_27-Implementationofthemainclass"></a>
959 * <a name=
"step_27-LaplaceProblemLaplaceProblemconstructor"></a>
981 *
, max_degree(dim <= 2 ? 7 : 5)
983 *
for (unsigned
int degree = 2; degree <= max_degree; ++degree)
985 *
fe_collection.push_back(
FE_Q<dim>(degree));
995 * <a name=
"step_27-LaplaceProblemLaplaceProblemdestructor"></a>
1003 *
template <
int dim>
1006 *
dof_handler.
clear();
1013 * <a name=
"step_27-LaplaceProblemsetup_system"></a>
1024 *
template <
int dim>
1027 *
dof_handler.distribute_dofs(fe_collection);
1029 *
solution.reinit(dof_handler.n_dofs());
1032 *
constraints.
clear();
1038 *
constraints.close();
1042 *
sparsity_pattern.copy_from(
dsp);
1044 *
system_matrix.reinit(sparsity_pattern);
1052 * <a name=
"step_27-LaplaceProblemassemble_system"></a>
1100 *
std::vector<types::global_dof_index> local_dof_indices;
1102 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1104 *
const unsigned
int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1106 *
cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1112 *
hp_fe_values.reinit(cell);
1114 *
const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
1116 *
std::vector<double>
rhs_values(fe_values.n_quadrature_points);
1119 *
for (
unsigned int q_point = 0;
q_point < fe_values.n_quadrature_points;
1121 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1123 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1125 *
(fe_values.shape_grad(i,
q_point) *
1126 *
fe_values.shape_grad(
j,
q_point) *
1134 *
local_dof_indices.resize(dofs_per_cell);
1135 *
cell->get_dof_indices(local_dof_indices);
1137 *
constraints.distribute_local_to_global(
1147 * <a name=
"step_27-LaplaceProblemsolve"></a>
1148 * <
h4>LaplaceProblem::solve</
h4>
1157 *
template <
int dim>
1165 *
preconditioner.initialize(system_matrix, 1.2);
1167 *
cg.solve(system_matrix, solution,
system_rhs, preconditioner);
1169 *
constraints.distribute(solution);
1177 * <a name=
"step_27-LaplaceProblempostprocess"></a>
1178 * <
h4>LaplaceProblem::postprocess</
h4>
1191 *
template <
int dim>
1205 *
face_quadrature_collection,
1252 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1254 *
fe_collection[cell->active_fe_index()].degree;
1265 *
data_out.attach_dof_handler(dof_handler);
1266 *
data_out.add_data_vector(solution,
"solution");
1269 *
data_out.add_data_vector(
fe_degrees,
"fe_degree");
1270 *
data_out.build_patches();
1281 *
data_out.write_vtk(output);
1353 *
hp::Refinement::limit_p_level_difference(dof_handler);
1399 *
std::set<typename Triangulation<dim>::active_cell_iterator>
cells_to_remove;
1400 *
for (
const auto &cell :
cube.active_cell_iterators())
1402 *
if (cell->vertex(v).square() < .1)
1417 * <a name=
"step_27-LaplaceProblemrun"></a>
1418 * <
h4>LaplaceProblem::run</
h4>
1434 *
template <
int dim>
1437 *
for (
unsigned int cycle = 0; cycle < 6; ++cycle)
1439 *
std::cout <<
"Cycle " << cycle <<
':' << std::endl;
1446 *
std::cout <<
" Number of active cells : "
1448 *
<<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1450 *
<<
" Number of constraints : "
1451 *
<< constraints.n_constraints() << std::endl;
1455 *
postprocess(cycle);
1464 * <a name=
"step_27-Themainfunction"></a>
1479 *
using namespace Step27;
1484 *
catch (std::exception &exc)
1486 *
std::cerr << std::endl
1488 *
<<
"----------------------------------------------------"
1490 *
std::cerr <<
"Exception on processing: " << std::endl
1491 *
<< exc.what() << std::endl
1492 *
<<
"Aborting!" << std::endl
1493 *
<<
"----------------------------------------------------"
1500 *
std::cerr << std::endl
1502 *
<<
"----------------------------------------------------"
1504 *
std::cerr <<
"Unknown exception!" << std::endl
1505 *
<<
"Aborting!" << std::endl
1506 *
<<
"----------------------------------------------------"
1531 Number
of constraints : 384
1535 Number
of constraints : 756
1539 Number
of constraints : 1856
1543 Number
of constraints : 2944
1547 Number
of constraints : 3998
1551 Number
of constraints : 5230
1572<
img src=
"https://www.dealii.org/images/steps/developer/step-27-solution.png"
1573 alt=
"Elevation plot of the solution, showing the lack of regularity near
1574 the interior (reentrant) corners."
1575 width=
"200" height=
"200">
1579<
div class=
"threecolumn" style=
"width: 80%">
1581 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-00.svg"
1582 alt=
"Triangulation containing reentrant corners without adaptive refinement."
1583 width=
"200" height=
"200">
1586 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-01.svg"
1587 alt=
"Triangulation containing reentrant corners with one level of
1588 refinement. New cells are placed near the corners."
1589 width=
"200" height=
"200">
1592 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-02.svg"
1593 alt=
"Triangulation containing reentrant corners with two levels of
1594 refinement. New cells are placed near the corners."
1595 width=
"200" height=
"200">
1598 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-03.svg"
1599 alt=
"Triangulation containing reentrant corners with three levels of
1600 refinement. New cells are placed near the corners."
1601 width=
"200" height=
"200">
1604 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-04.svg"
1605 alt=
"Triangulation containing reentrant corners with four levels of
1606 refinement. New cells are placed near the corners."
1607 width=
"200" height=
"200">
1610 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-05.svg"
1611 alt=
"Triangulation containing reentrant corners with five levels of
1612 refinement. New cells are placed near the corners."
1613 width=
"200" height=
"200">
1623<
div class=
"threecolumn" style=
"width: 80%">
1625 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-00.svg"
1626 alt=
"Initial grid where all cells contain just biquadratic functions."
1627 width=
"200" height=
"200">
1630 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-01.svg"
1631 alt=
"Depiction of local approximation degrees after one refinement."
1632 width=
"200" height=
"200">
1635 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-02.svg"
1636 alt=
"Depiction of local approximation degrees after two refinements."
1637 width=
"200" height=
"200">
1640 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-03.svg"
1641 alt=
"Depiction of local approximation degrees after three refinements."
1642 width=
"200" height=
"200">
1645 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-04.svg"
1646 alt=
"Depiction of local approximation degrees after four refinements."
1647 width=
"200" height=
"200">
1650 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-05.svg"
1651 alt=
"Depiction of local approximation degrees after five refinements."
1652 width=
"200" height=
"200">
1667<
div class=
"threecolumn" style=
"width: 80%">
1669 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-00.svg"
1670 alt=
"Estimated regularity per cell on the initial grid."
1671 width=
"200" height=
"200">
1674 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-01.svg"
1675 alt=
"Depiction of the estimated regularity per cell after one refinement."
1676 width=
"200" height=
"200">
1679 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-02.svg"
1680 alt=
"Depiction of the estimated regularity per cell after two refinements."
1681 width=
"200" height=
"200">
1684 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-03.svg"
1685 alt=
"Depiction of the estimated regularity per cell after three refinements."
1686 width=
"200" height=
"200">
1689 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-04.svg"
1690 alt=
"Depiction of the estimated regularity per cell after four refinements."
1691 width=
"200" height=
"200">
1694 <
img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-05.svg"
1695 alt=
"Depiction of the estimated regularity per cell after five refinements."
1696 width=
"200" height=
"200">
1734<a name=
"step-27-extensions"></a>
1735<a name=
"step_27-Possibilitiesforextensions"></a><
h3>Possibilities
for extensions</
h3>
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)
numbers::NumberTraits< Number >::real_type norm() const
static constexpr unsigned int dimension
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
virtual bool prepare_coarsening_and_refinement()
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_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
void create_triangulation_with_removed_cells(const Triangulation< dim, spacedim > &input_triangulation, const std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_remove, Triangulation< dim, spacedim > &result)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
@ general
No special properties.
constexpr types::blas_int zero
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.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
FESeries::Fourier< dim, spacedim > default_fe_series(const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)
void coefficient_decay(FESeries::Fourier< dim, spacedim > &fe_fourier, 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 reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
std::string int_to_string(const unsigned int 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)
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 choose_p_over_h(const DoFHandler< dim, spacedim > &dof_handler)
void p_adaptivity_from_relative_threshold(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 >())
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)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation