941 *
double product = 1;
942 *
for (
unsigned int d = 0;
d < dim; ++
d)
943 * product *= (p[d] + 1);
952 * <a name=
"step_27-Implementationofthemainclass"></a>
953 * <h3>Implementation of the main
class</h3>
958 * <a name=
"step_27-LaplaceProblemLaplaceProblemconstructor"></a>
959 * <h4>LaplaceProblem::LaplaceProblem constructor</h4>
963 * The constructor of
this class is fairly straightforward. It associates
965 * maximal polynomial degree to 7 (in 1
d and 2
d) or 5 (in 3d and higher). We
966 *
do so because
using higher order polynomial degrees becomes prohibitively
967 * expensive, especially in higher space dimensions.
971 * Following
this, we fill the collections of finite element, and cell and
972 * face quadrature objects. We start with quadratic elements, and each
973 * quadrature formula is chosen so that it is appropriate
for the
matching
978 * LaplaceProblem<dim>::LaplaceProblem()
980 * , max_degree(dim <= 2 ? 7 : 5)
982 * for (unsigned
int degree = 2; degree <= max_degree; ++degree)
985 * quadrature_collection.push_back(
QGauss<dim>(degree + 1));
994 * <a name=
"step_27-LaplaceProblemLaplaceProblemdestructor"></a>
995 * <h4>LaplaceProblem::~LaplaceProblem destructor</h4>
999 * The destructor is unchanged from what we already did in @ref step_6
"step-6":
1002 *
template <
int dim>
1003 * LaplaceProblem<dim>::~LaplaceProblem()
1005 * dof_handler.clear();
1012 * <a name=
"step_27-LaplaceProblemsetup_system"></a>
1013 * <h4>LaplaceProblem::setup_system</h4>
1017 * This function is again a verbatim
copy of what we already did in
1018 * @ref step_6
"step-6". Despite function calls with exactly the same names and arguments,
1019 * the algorithms used internally are different in some aspect since the
1020 * dof_handler variable here is in @f$hp@f$-mode.
1023 *
template <
int dim>
1024 *
void LaplaceProblem<dim>::setup_system()
1026 * dof_handler.distribute_dofs(fe_collection);
1028 * solution.reinit(dof_handler.n_dofs());
1029 * system_rhs.reinit(dof_handler.n_dofs());
1031 * constraints.clear();
1037 * constraints.close();
1041 * sparsity_pattern.copy_from(dsp);
1043 * system_matrix.reinit(sparsity_pattern);
1051 * <a name=
"step_27-LaplaceProblemassemble_system"></a>
1052 * <h4>LaplaceProblem::assemble_system</h4>
1056 * This is the function that assembles the global
matrix and right hand side
1057 * vector from the local contributions of each cell. Its main working is as
1058 * has been described in many of the tutorial programs before. The
1059 * significant deviations are the ones necessary for <i>
hp</i> finite
1060 * element methods. In particular, that we need to use a collection of
1062 * have to eliminate constrained degrees of freedom already when copying
1063 * local contributions into global objects. Both of these are explained in
1064 * detail in the introduction of
this program.
1068 * One other slight complication is the fact that because we use different
1069 * polynomial degrees on different cells, the matrices and vectors holding
1070 * local contributions
do not have the same size on all cells. At the
1071 * beginning of the
loop over all cells, we therefore each time have to
1072 * resize them to the correct size (given by <code>dofs_per_cell</code>).
1073 * Because these classes are implemented in such a way that reducing the size
1074 * of a
matrix or vector does not release the currently allocated memory
1075 * (unless the
new size is zero), the process of resizing at the beginning of
1076 * the loop will only require re-allocation of memory during the
first few
1077 * iterations. Once we have found in a cell with the maximal finite element
1078 * degree, no more re-allocations will happen because all subsequent
1079 * <code>reinit</code> calls will only set the size to something that fits the
1080 * currently allocated memory. This is important since allocating memory is
1081 * expensive, and doing so every time we visit a
new cell would take
1082 * significant compute time.
1085 * template <int dim>
1086 *
void LaplaceProblem<dim>::assemble_system()
1089 * quadrature_collection,
1094 * RightHandSide<dim> rhs_function;
1099 * std::vector<types::global_dof_index> local_dof_indices;
1101 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1103 * const unsigned
int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1105 *
cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1108 * cell_rhs.reinit(dofs_per_cell);
1111 * hp_fe_values.reinit(cell);
1115 * std::vector<double> rhs_values(fe_values.n_quadrature_points);
1116 * rhs_function.value_list(fe_values.get_quadrature_points(), rhs_values);
1118 *
for (
unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
1120 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1122 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1124 * (fe_values.shape_grad(i, q_point) *
1125 * fe_values.shape_grad(j, q_point) *
1126 * fe_values.JxW(q_point));
1128 * cell_rhs(i) += (fe_values.shape_value(i, q_point) *
1129 * rhs_values[q_point] *
1130 * fe_values.JxW(q_point));
1133 * local_dof_indices.resize(dofs_per_cell);
1134 * cell->get_dof_indices(local_dof_indices);
1136 * constraints.distribute_local_to_global(
1137 * cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1146 * <a name=
"step_27-LaplaceProblemsolve"></a>
1147 * <h4>LaplaceProblem::solve</h4>
1151 * The function solving the linear system is entirely unchanged from
1152 * previous examples. We simply
try to
reduce the
initial residual (which
1153 * equals the @f$l_2@f$ norm of the right hand side) by a certain factor:
1156 *
template <
int dim>
1157 *
void LaplaceProblem<dim>::solve()
1160 * 1e-12 * system_rhs.l2_norm());
1164 * preconditioner.
initialize(system_matrix, 1.2);
1166 * cg.solve(system_matrix, solution, system_rhs, preconditioner);
1168 * constraints.distribute(solution);
1176 * <a name=
"step_27-LaplaceProblempostprocess"></a>
1177 * <h4>LaplaceProblem::postprocess</h4>
1181 * After solving the linear system, we will want to postprocess the
1182 * solution. Here, all we
do is to estimate the error, estimate the local
1183 * smoothness of the solution as described in the introduction, then write
1184 * graphical output, and
finally refine the mesh in both @f$h@f$ and @f$p@f$
1185 * according to the indicators computed before. We
do all
this in the same
1186 * function because we want the estimated error and smoothness indicators
1187 * not only
for refinement, but also include them in the graphical output.
1190 *
template <
int dim>
1191 *
void LaplaceProblem<dim>::postprocess(
const unsigned int cycle)
1195 * Let us start with computing estimated error and smoothness indicators,
1196 * which each are one number
for each active cell of our
1204 * face_quadrature_collection,
1207 * estimated_error_per_cell);
1211 * Estimating the smoothness is performed with the method of decaying
1212 * expansion coefficients as outlined in the introduction. We will
first
1213 * need to create an
object capable of transforming the finite element
1214 * solution on every single cell into a sequence of Fourier series
1217 * estimating smoothness. The actual determination of the decay of Fourier
1218 * coefficients on every individual cell then happens in the last function.
1227 * smoothness_indicators);
1231 * Next we want to generate graphical output. In addition to the two
1232 * estimated quantities derived above, we would also like to output the
1233 * polynomial degree of the finite elements used on each of the elements
1238 * The way to
do that
requires that we
loop over all cells and poll the
1239 * active finite element
index of them
using
1240 * <code>cell-@>active_fe_index()</code>. We then use the result of
this
1241 * operation and query the finite element collection
for the finite
1242 * element with that
index, and
finally determine the polynomial degree of
1243 * that element. The result we put into a vector with one element per
1244 * cell. The
DataOut class requires this to be a vector of
1245 * <code>
float</code> or <code>
double</code>, even though our
values are
1246 * all integers, so that is what we use:
1251 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1252 * fe_degrees(cell->active_cell_index()) =
1253 * fe_collection[cell->active_fe_index()].degree;
1257 * With now all data vectors available -- solution, estimated errors and
1258 * smoothness indicators, and finite element degrees --, we create a
1259 *
DataOut object for graphical output and attach all data:
1265 * data_out.add_data_vector(solution,
"solution");
1266 * data_out.add_data_vector(estimated_error_per_cell,
"error");
1267 * data_out.add_data_vector(smoothness_indicators,
"smoothness");
1268 * data_out.add_data_vector(fe_degrees,
"fe_degree");
1269 * data_out.build_patches();
1273 * The
final step in generating output is to determine a file name, open
1274 * the file, and write the data into it (here, we use VTK format):
1277 * const
std::string filename =
1279 * std::ofstream output(filename);
1280 * data_out.write_vtk(output);
1285 * After
this, we would like to actually
refine the mesh, in both @f$h@f$ and
1286 * @f$p@f$. The way we are going to
do this is as follows:
first, we use the
1287 * estimated error to flag those cells
for refinement that have the
1288 * largest error. This is what we have
always done:
1293 * estimated_error_per_cell,
1299 * Next we would like to figure out which of the cells that have been
1300 * flagged
for refinement should actually have @f$p@f$ increased instead of
1301 * @f$h@f$ decreased. The strategy we choose here is that we look at the
1302 * smoothness indicators of those cells that are flagged
for refinement,
1303 * and increase @f$p@f$
for those with a smoothness larger than a certain
1304 * relative threshold. In other words,
for every cell
for which (i) the
1305 * refinement flag is set, (ii) the smoothness indicator is larger than
1306 * the threshold, and (iii) we still have a finite element with a
1307 * polynomial degree higher than the current one in the finite element
1308 * collection, we will assign a future FE
index that corresponds to a
1309 * polynomial with degree one higher than it currently is. The following
1310 * function is capable of doing exactly this. Absent any better
1311 * strategies, we will set the threshold via interpolation between the
1312 * minimal and maximal smoothness indicators on cells flagged
for
1313 * refinement. Since the corner singularities are strongly localized, we
1314 * will favor @f$p@f$- over @f$h@f$-refinement quantitatively. We achieve
this
1315 * with a low threshold by setting a small interpolation factor of 0.2. In
1316 * the same way, we deal with cells that are going to be coarsened and
1317 * decrease their polynomial degree when their smoothness indicator is
1318 * below the corresponding threshold determined on cells to be coarsened.
1322 * dof_handler, smoothness_indicators, 0.2, 0.2);
1326 * The above function only determines whether the polynomial degree will
1327 * change via future FE indices, but does not manipulate the
1328 * @f$h@f$-refinement flags. So
for cells that are flagged
for both refinement
1329 * categories, we prefer @f$p@f$- over @f$h@f$-refinement. The following function
1330 * call ensures that only one of @f$p@f$- or @f$h@f$-refinement is imposed, and
1338 * For grid adaptive refinement, we ensure a 2:1 mesh balance by limiting
1339 * the difference of refinement levels of neighboring cells to one by
1341 * like to achieve something similar
for the p-levels of neighboring
1342 * cells: levels of future finite elements are not allowed to differ by
1343 * more than a specified difference. With its
default parameters, a call
1345 * difference is limited to one. This will not necessarily decrease the
1346 * number of hanging nodes in the domain, but makes sure that high order
1347 * polynomials are not constrained to much lower polynomials on faces,
1348 * e.g., fifth order to
second order polynomials.
1352 *
hp::Refinement::limit_p_level_difference(dof_handler);
1356 * At the end of this procedure, we then refine the mesh. During this
1357 * process, children of cells undergoing bisection inherit their mother
1358 * cell's finite element index. Further, future finite element indices
1359 * will turn into active ones, so that the new finite elements will be
1360 * assigned to cells after the next call of
DoFHandler::distribute_dofs().
1371 * <a name="step_27-LaplaceProblemcreate_coarse_grid"></a>
1372 * <h4>LaplaceProblem::create_coarse_grid</h4>
1376 * The following function is used when creating the initial grid. The grid we
1377 * would like to create is actually similar to the one from @ref step_14 "step-14", i.e., the
1378 * square domain with the square hole in the middle. It can be generated by
1379 * exactly the same function. However, since its implementation is only a
1380 * specialization of the 2d case, we will present a different way of creating
1386 * already holds our desired domain @f$[-1,1]^d@f$, subdivided into @f$4^d@f$ cells.
1387 * We then remove those cells in the
center of the domain by testing the
1388 * coordinate values of the
vertices on each cell. In the end, we refine the
1389 * so created grid globally as usual.
1392 * template <
int dim>
1393 *
void LaplaceProblem<dim>::create_coarse_grid()
1398 * std::set<typename Triangulation<dim>::active_cell_iterator> cells_to_remove;
1399 *
for (
const auto &cell : cube.active_cell_iterators())
1400 * for (unsigned
int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
1401 *
if (cell->vertex(v).square() < .1)
1402 * cells_to_remove.insert(cell);
1416 * <a name=
"step_27-LaplaceProblemrun"></a>
1417 * <h4>LaplaceProblem::run</h4>
1421 * This function implements the logic of the program, as did the respective
1422 * function in most of the previous programs already, see
for example @ref step_6
"step-6".
1426 * Basically, it contains the adaptive
loop: in the
first iteration create a
1427 * coarse grid, and then set up the linear system,
assemble it, solve, and
1428 * postprocess the solution including mesh refinement. Then start over
1429 * again. In the meantime, also output some information
for those staring at
1430 * the screen trying to figure out what the program does:
1433 *
template <
int dim>
1434 *
void LaplaceProblem<dim>::run()
1436 *
for (
unsigned int cycle = 0; cycle < 6; ++cycle)
1438 * std::cout <<
"Cycle " << cycle <<
':' << std::endl;
1441 * create_coarse_grid();
1445 * std::cout <<
" Number of active cells : "
1447 * <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1449 * <<
" Number of constraints : "
1450 * << constraints.n_constraints() << std::endl;
1452 * assemble_system();
1454 * postprocess(cycle);
1463 * <a name=
"step_27-Themainfunction"></a>
1464 * <h3>The main function</h3>
1468 * The main function is again verbatim what we had before: wrap creating and
1469 * running an
object of the main
class into a <code>try</code> block and catch
1470 * whatever exceptions are thrown, thereby producing meaningful output
if
1471 * anything should go wrong:
1478 *
using namespace Step27;
1480 * LaplaceProblem<2> laplace_problem;
1481 * laplace_problem.run();
1483 *
catch (std::exception &exc)
1485 * std::cerr << std::endl
1487 * <<
"----------------------------------------------------"
1489 * std::cerr <<
"Exception on processing: " << std::endl
1490 * << exc.what() << std::endl
1491 * <<
"Aborting!" << std::endl
1492 * <<
"----------------------------------------------------"
1499 * std::cerr << std::endl
1501 * <<
"----------------------------------------------------"
1503 * std::cerr <<
"Unknown exception!" << std::endl
1504 * <<
"Aborting!" << std::endl
1505 * <<
"----------------------------------------------------"
1513<a name=
"step_27-Results"></a><h1>Results</h1>
1516In
this section, we discuss a few results produced from running the
1517current tutorial program. More results, in particular the extension to
15183
d calculations and determining how much compute time the individual
1519components of the program take, are given in the @ref hp_paper
"hp-paper".
1521When
run,
this is what the program produces:
1525[ 66%] Built target @ref step_27
"step-27"
1526[100%] Run @ref step_27
"step-27" with Release configuration
1528 Number of active cells : 768
1529 Number of degrees of freedom: 3264
1530 Number of constraints : 384
1532 Number of active cells : 807
1533 Number of degrees of freedom: 4764
1534 Number of constraints : 756
1536 Number of active cells : 927
1537 Number of degrees of freedom: 8226
1538 Number of constraints : 1856
1540 Number of active cells : 978
1541 Number of degrees of freedom: 12146
1542 Number of constraints : 2944
1544 Number of active cells : 1104
1545 Number of degrees of freedom: 16892
1546 Number of constraints : 3998
1548 Number of active cells : 1149
1549 Number of degrees of freedom: 22078
1550 Number of constraints : 5230
1553The
first thing we learn from
this is that the number of constrained degrees
1554of freedom is on the order of 20-25% of the total number of degrees of
1555freedom, at least on the later grids when we have elements of relatively
1556high order (in 3d, the fraction of constrained degrees of freedom can be up
1557to 30%). This is, in fact, on the same order of magnitude as
for
1558non-@f$hp@f$-discretizations. For example, in the last step of the @ref step_6
"step-6"
1559program, we have 18353 degrees of freedom, 4432 of which are
1560constrained. The difference is that in the latter program, each constrained
1561hanging node is constrained against only the two adjacent degrees of
1562freedom, whereas in the @f$hp@f$-
case, constrained nodes are constrained against
1563many more degrees of freedom. Note also that the current program also
1564includes nodes subject to Dirichlet boundary conditions in the list of
1565constraints. In cycle 0, all the constraints are actually because of
1568Of maybe more interest is to look at the graphical output. First, here is the
1569solution of the problem:
1571<img src=
"https://www.dealii.org/images/steps/developer/step-27-solution.png"
1572 alt=
"Elevation plot of the solution, showing the lack of regularity near
1573 the interior (reentrant) corners."
1574 width=
"200" height=
"200">
1576Secondly, let us look at the sequence of meshes generated:
1578<div
class=
"threecolumn" style=
"width: 80%">
1580 <img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-00.svg"
1581 alt=
"Triangulation containing reentrant corners without adaptive refinement."
1582 width=
"200" height=
"200">
1585 <img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-01.svg"
1586 alt=
"Triangulation containing reentrant corners with one level of
1587 refinement. New cells are placed near the corners."
1588 width=
"200" height=
"200">
1591 <img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-02.svg"
1592 alt=
"Triangulation containing reentrant corners with two levels of
1593 refinement. New cells are placed near the corners."
1594 width=
"200" height=
"200">
1597 <img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-03.svg"
1598 alt=
"Triangulation containing reentrant corners with three levels of
1599 refinement. New cells are placed near the corners."
1600 width=
"200" height=
"200">
1603 <img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-04.svg"
1604 alt=
"Triangulation containing reentrant corners with four levels of
1605 refinement. New cells are placed near the corners."
1606 width=
"200" height=
"200">
1609 <img src=
"https://www.dealii.org/images/steps/developer/step-27.mesh-05.svg"
1610 alt=
"Triangulation containing reentrant corners with five levels of
1611 refinement. New cells are placed near the corners."
1612 width=
"200" height=
"200">
1616It is clearly visible how the mesh is refined near the corner singularities,
1617as one would expect it. More interestingly, we should be curious to see the
1618distribution of finite element polynomial degrees to these mesh cells, where
1619the lightest color corresponds to degree two and the darkest one corresponds
1622<div
class=
"threecolumn" style=
"width: 80%">
1624 <img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-00.svg"
1625 alt=
"Initial grid where all cells contain just biquadratic functions."
1626 width=
"200" height=
"200">
1629 <img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-01.svg"
1630 alt=
"Depiction of local approximation degrees after one refinement."
1631 width=
"200" height=
"200">
1634 <img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-02.svg"
1635 alt=
"Depiction of local approximation degrees after two refinements."
1636 width=
"200" height=
"200">
1639 <img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-03.svg"
1640 alt=
"Depiction of local approximation degrees after three refinements."
1641 width=
"200" height=
"200">
1644 <img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-04.svg"
1645 alt=
"Depiction of local approximation degrees after four refinements."
1646 width=
"200" height=
"200">
1649 <img src=
"https://www.dealii.org/images/steps/developer/step-27.fe_degree-05.svg"
1650 alt=
"Depiction of local approximation degrees after five refinements."
1651 width=
"200" height=
"200">
1655While
this is certainly not a perfect arrangement, it does make some sense: we
1656use low order elements close to boundaries and corners where regularity is
1657low. On the other hand, higher order elements are used where (i) the error was
1658at one
point fairly large, i.e., mainly in the
general area around the corner
1659singularities and in the top right corner where the solution is large, and
1660(ii) where the solution is smooth, i.e., far away from the boundary.
1662This arrangement of polynomial degrees of course follows from our smoothness
1663estimator. Here is the estimated smoothness of the solution, with darker colors
1664indicating least smoothness and lighter indicating the smoothest areas:
1666<div
class=
"threecolumn" style=
"width: 80%">
1668 <img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-00.svg"
1669 alt=
"Estimated regularity per cell on the initial grid."
1670 width=
"200" height=
"200">
1673 <img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-01.svg"
1674 alt=
"Depiction of the estimated regularity per cell after one refinement."
1675 width=
"200" height=
"200">
1678 <img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-02.svg"
1679 alt=
"Depiction of the estimated regularity per cell after two refinements."
1680 width=
"200" height=
"200">
1683 <img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-03.svg"
1684 alt=
"Depiction of the estimated regularity per cell after three refinements."
1685 width=
"200" height=
"200">
1688 <img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-04.svg"
1689 alt=
"Depiction of the estimated regularity per cell after four refinements."
1690 width=
"200" height=
"200">
1693 <img src=
"https://www.dealii.org/images/steps/developer/step-27.smoothness-05.svg"
1694 alt=
"Depiction of the estimated regularity per cell after five refinements."
1695 width=
"200" height=
"200">
1699The primary conclusion one can draw from
this is that the loss of regularity at
1700the
internal corners is a highly localized phenomenon; it only seems to impact
1701the cells adjacent to the corner itself, so when we
refine the mesh the black
1702coloring is no longer visible. Besides the corners,
this sequence of plots
1703implies that the smoothness estimates are somewhat
independent of the mesh
1704refinement, particularly when we are far away from boundaries.
1705It is also obvious that the smoothness estimates are
independent of the actual
1706size of the solution (see the picture of the solution above), as it should be.
1707A
point of larger concern, however, is that one realizes on closer inspection
1708that the estimator we have overestimates the smoothness of the solution on
1709cells with hanging nodes. This in turn leads to higher polynomial degrees in
1710these areas, skewing the allocation of finite elements onto cells.
1712We have no good explanation
for this effect at the moment. One theory is that
1713the numerical solution on cells with hanging nodes is, of course, constrained
1714and therefore not entirely free to explore the function space to get close to
1715the exact solution. This lack of degrees of freedom may manifest itself by
1716yielding numerical solutions on these cells with suppressed oscillation,
1717meaning a higher degree of smoothness. The estimator picks
this signal up and
1718the estimated smoothness overestimates the actual
value. However, a definite
1719answer to what is going on currently eludes the authors of
this program.
1721The bigger question is, of course, how to avoid
this problem. Possibilities
1722include estimating the smoothness not on single cells, but cell assemblies or
1723patches surrounding each cell. It may also be possible to find simple
1724correction factors
for each cell depending on the number of constrained
1725degrees of freedom it has. In either
case, there are ample opportunities
for
1726further research on finding good @f$hp@f$-refinement criteria. On the other hand,
1727the main
point of the current program was to demonstrate
using the
1728@f$hp@f$-technology in deal.II, which is unaffected by our use of a possible
1729sub-optimal refinement criterion.
1733<a name=
"step-27-extensions"></a>
1734<a name=
"step_27-Possibilitiesforextensions"></a><h3>Possibilities
for extensions</h3>
1737<a name=
"step_27-Differenthpdecisionstrategies"></a><h4>Different
hp-decision strategies</h4>
1740This tutorial demonstrates only one particular strategy to decide between @f$h@f$- and
1741@f$p@f$-adaptation. In fact, there are many more ways to automatically decide on the
1742adaptation type, of which a few are already implemented in deal.II:
1744 <li><i>Fourier coefficient decay:</i> This is the strategy currently
1745 implemented in
this tutorial. For more information on
this strategy, see
1748 <li><i>Legendre coefficient decay:</i> This strategy is quite similar
1749 to the current one, but uses Legendre series expansion rather than the
1750 Fourier one: instead of sinusoids as basis
functions,
this strategy uses
1751 Legendre polynomials. Of course, since we
approximate the solution
using a
1752 finite-dimensional polynomial on each cell, the expansion of the solution in
1753 Legendre polynomials is also finite and, consequently, when we talk about the
1754 "decay" of
this expansion, we can only consider the finitely many
nonzero
1755 coefficients of
this expansion, rather than think about it in asymptotic terms.
1756 But,
if we have enough of these coefficients, we can certainly think of the
1757 decay of these coefficients as characteristic of the decay of the coefficients
1758 of the exact solution (which is, in general, not polynomial and so will have an
1759 infinite Legendre expansion), and considering the coefficients we have should
1760 reveal something about the properties of the exact solution.
1762 The transition from the Fourier strategy to the Legendre one is quite simple:
1763 You just need to change the series expansion
class and the corresponding
1764 smoothness estimation function to be part of the proper namespaces
1766 in @ref step_75 "step-75". For the theoretical background of this strategy, consult the
1768 as @cite mavriplis1994hp , @cite eibner2007hp and @cite davydov2017hp .</li>
1770 <li><i>Refinement history:</i> The last strategy is quite different
1771 from the other two. In theory, we know how the error will converge
1772 after changing the discretization of the function space. With
1773 @f$h@f$-refinement the solution converges algebraically as already pointed
1774 out in @ref step_7
"step-7". If the solution is sufficiently smooth, though, we
1775 expect that the solution will converge exponentially with increasing
1776 polynomial degree of the finite element. We can compare a proper
1777 prediction of the error with the actual error in the following step to
1778 see
if our choice of adaptation type was justified.
1780 The transition to
this strategy is a bit more complicated. For
this, we need
1781 an initialization step with pure @f$h@f$- or @f$p@f$-refinement and we need to
1782 transfer the predicted errors over adapted meshes. The extensive
1784 only the theoretical details of this approach, but also presents a blueprint
1785 on how to implement this strategy in your code. For more information, see
1786 @cite melenk2001hp .
1788 Note that with this particular function you cannot predict the error for
1789 the next time step in time-dependent problems. Therefore, this strategy
1790 cannot be applied to this type of problem without further ado. Alternatively,
1791 the following approach could be used, which works for all the other
1792 strategies as well: start each time step with a coarse mesh, keep refining
1793 until happy with the result, and only then move on to the next time step.</li>
1796Try implementing one of these strategies into this tutorial and observe the
1797subtle changes to the results. You will notice that all strategies are
1798capable of identifying the singularities near the reentrant corners and
1799will perform @f$h@f$-refinement in these regions, while preferring @f$p@f$-refinement
1800in the bulk domain. A detailed comparison of these strategies is presented
1801in @cite fehling2020 .
1804<a name="step_27-Parallelhpadaptivefiniteelements"></a><h4>Parallel
hp-adaptive finite elements</h4>
1807All functionality presented in this tutorial already works for both
1808sequential and
parallel applications. It is possible without too much
1811it, we recommend reading @ref step_18 "step-18" for the former and @ref step_40 "step-40" for the
1812latter case
first for further background information on the topic, and
1813then come back to this tutorial to try out your newly acquired skills.
1815We go one step further in @ref step_75 "step-75": Here, we combine
hp-adaptive and
1820<a name="step_27-PlainProg"></a>
1821<h1> The plain program</h1>
1822@include "step-27.cc"
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
const FEValues< dim, spacedim > & get_present_fe_values() const
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 MatrixType &A, const AdditionalData ¶meters=AdditionalData())
virtual bool prepare_coarsening_and_refinement()
void push_back(const FiniteElement< dim, spacedim > &new_fe)
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.
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.
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