757 * <a name=
"step_14-TheLaplacesolverclasses"></a>
771 * <a name=
"step_14-TheLaplacesolverbaseclass"></a>
772 * <
h4>
The Laplace solver base
class</
h4>
790 *
virtual ~Base() =
default;
793 *
virtual void postprocess(
795 *
virtual void refine_grid() = 0;
796 *
virtual unsigned int n_dofs()
const = 0;
818 *
void Base<dim>::set_refinement_cycle(
const unsigned int cycle)
827 * <a name=
"step_14-TheLaplaceSolverclass"></a>
837 *
class Solver :
public virtual Base<dim>
849 *
virtual void postprocess(
852 *
virtual unsigned int n_dofs()
const override;
898 *
std::vector<types::global_dof_index> local_dof_indices;
924 *
, quadrature(&quadrature)
925 *
, face_quadrature(&face_quadrature)
934 *
dof_handler.
clear();
941 *
dof_handler.distribute_dofs(*fe);
942 *
solution.reinit(dof_handler.n_dofs());
954 *
postprocessor(dof_handler, solution);
961 *
return dof_handler.n_dofs();
1012 *
template <
int dim>
1020 *
template <
int dim>
1023 *
: fe_values(scratch_data.fe_values.get_fe(),
1024 *
scratch_data.fe_values.get_quadrature(),
1029 *
template <
int dim>
1035 *
const unsigned int dofs_per_cell = fe->n_dofs_per_cell();
1036 *
const unsigned int n_q_points = quadrature->size();
1038 *
copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1040 *
copy_data.local_dof_indices.resize(dofs_per_cell);
1042 *
scratch_data.fe_values.reinit(cell);
1045 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1046 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1047 *
copy_data.cell_matrix(i,
j) +=
1048 *
(scratch_data.fe_values.shape_grad(i,
q_point) *
1049 *
scratch_data.fe_values.shape_grad(
j,
q_point) *
1050 *
scratch_data.fe_values.JxW(
q_point));
1052 *
cell->get_dof_indices(copy_data.local_dof_indices);
1057 *
template <
int dim>
1061 *
for (
unsigned int i = 0; i < copy_data.local_dof_indices.size(); ++i)
1062 *
for (
unsigned int j = 0;
j < copy_data.local_dof_indices.size(); ++
j)
1064 *
copy_data.local_dof_indices[
j],
1065 *
copy_data.cell_matrix(i,
j));
1082 * class here, but rather use the one created task object
1083 * directly to wait for this particular task's
exit.
The
1140 *
sparsity_pattern.copy_from(
dsp);
1142 *
matrix.reinit(sparsity_pattern);
1143 *
rhs.reinit(dof_handler.n_dofs());
1148 *
template <
int dim>
1155 *
preconditioner.initialize(matrix, 1.2);
1157 *
cg.solve(matrix, solution,
rhs, preconditioner);
1167 * <a name=
"step_14-ThePrimalSolverclass"></a>
1181 *
template <
int dim>
1200 *
template <
int dim>
1218 *
template <
int dim>
1222 *
data_out.attach_dof_handler(this->dof_handler);
1223 *
data_out.add_data_vector(this->solution,
"solution");
1224 *
data_out.build_patches();
1233 *
template <
int dim>
1237 *
*this->quadrature,
1241 *
const unsigned int dofs_per_cell = this->fe->n_dofs_per_cell();
1242 *
const unsigned int n_q_points = this->quadrature->size();
1245 *
std::vector<double>
rhs_values(n_q_points);
1246 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1248 *
for (
const auto &cell :
this->dof_handler.active_cell_iterators())
1252 *
fe_values.reinit(cell);
1254 *
rhs_function->value_list(fe_values.get_quadrature_points(),
1258 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1263 *
cell->get_dof_indices(local_dof_indices);
1264 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1273 * <a name=
"step_14-TheRefinementGlobalandRefinementKellyclasses"></a>
1293 *
virtual void refine_grid()
override;
1298 *
template <
int dim>
1317 *
template <
int dim>
1325 *
template <
int dim>
1336 *
virtual void refine_grid()
override;
1341 *
template <
int dim>
1360 *
template <
int dim>
1366 *
this->dof_handler,
1383 * <a name=
"step_14-TheRefinementWeightedKellyclass"></a>
1390 * indicator by some function. We include this class since the goal of
1391 * this example program is to demonstrate automatic refinement criteria
1392 * even for complex output quantities such as point values or stresses. If
1393 * we did not solve a dual problem and compute the weights thereof, we
1394 * would probably be tempted to give a hand-crafted weighting to the
1395 * indicators to account for the fact that we are going to evaluate these
1396 * quantities. This class accepts such a weighting function as argument to
1400 * template <int dim>
1401 * class RefinementWeightedKelly : public PrimalSolver<dim>
1404 * RefinementWeightedKelly(Triangulation<dim> &coarse_grid,
1405 * const FiniteElement<dim> &fe,
1406 * const Quadrature<dim> &quadrature,
1407 * const Quadrature<dim - 1> &face_quadrature,
1408 * const Function<dim> &rhs_function,
1409 * const Function<dim> &boundary_values,
1410 * const Function<dim> &weighting_function);
1412 * virtual void refine_grid() override;
1415 * const ObserverPointer<const Function<dim>> weighting_function;
1420 * template <int dim>
1421 * RefinementWeightedKelly<dim>::RefinementWeightedKelly(
1422 * Triangulation<dim> &coarse_grid,
1423 * const FiniteElement<dim> &fe,
1424 * const Quadrature<dim> &quadrature,
1425 * const Quadrature<dim - 1> &face_quadrature,
1426 * const Function<dim> &rhs_function,
1427 * const Function<dim> &boundary_values,
1428 * const Function<dim> &weighting_function)
1429 * : Base<dim>(coarse_grid)
1430 * , PrimalSolver<dim>(coarse_grid,
1436 * , weighting_function(&weighting_function)
1443 * Now, here comes the main function, including the weighting:
1446 * template <int dim>
1447 * void RefinementWeightedKelly<dim>::refine_grid()
1451 * First compute some residual based error indicators for all cells by a
1452 * method already implemented in the library. What exactly we compute
1453 * here is described in more detail in the documentation of that class.
1456 * Vector<float> estimated_error_per_cell(
1457 * this->triangulation->n_active_cells());
1458 * std::map<types::boundary_id, const Function<dim> *> dummy_function_map;
1459 * KellyErrorEstimator<dim>::estimate(this->dof_handler,
1460 * *this->face_quadrature,
1461 * dummy_function_map,
1463 * estimated_error_per_cell);
1467 * Next weigh each entry in the vector of indicators by the value of the
1468 * function given to the constructor, evaluated at the cell center. We
1469 * need to write the result into the vector entry that corresponds to the
1470 * current cell, which we can obtain by asking the cell what its index
1471 * among all active cells is using CellAccessor::active_cell_index(). (In
1472 * reality, this index is zero for the first cell we handle in the loop,
1473 * one for the second cell, etc., and we could as well just keep track of
1474 * this index using an integer counter; but using
1475 * CellAccessor::active_cell_index() makes this more explicit.)
1478 * for (const auto &cell : this->dof_handler.active_cell_iterators())
1479 * estimated_error_per_cell(cell->active_cell_index()) *=
1480 * weighting_function->value(cell->center());
1482 * GridRefinement::refine_and_coarsen_fixed_number(*this->triangulation,
1483 * estimated_error_per_cell,
1486 * this->triangulation->execute_coarsening_and_refinement();
1489 * } // namespace LaplaceSolver
1495 * <a name="step_14-Equationdata"></a>
1496 * <h3>Equation data</h3>
1500 * In this example program, we work with the same data sets as in the
1501 * previous one, but as it may so happen that someone wants to run the
1502 * program with different boundary values and right hand side functions, or
1503 * on a different grid, we show a simple technique to do exactly that. For
1504 * more clarity, we furthermore pack everything that has to do with equation
1505 * data into a namespace of its own.
1509 * The underlying assumption is that this is a research program, and that
1510 * there we often have a number of test cases that consist of a domain, a
1511 * right hand side, boundary values, possibly a specified coefficient, and a
1512 * number of other parameters. They often vary all at the same time when
1513 * shifting from one example to another. To make handling such sets of
1514 * problem description parameters simple is the goal of the following.
1518 * Basically, the idea is this: let us have a structure for each set of
1519 * data, in which we pack everything that describes a test case: here, these
1520 * are two subclasses, one called <code>BoundaryValues</code> for the
1521 * boundary values of the exact solution, and one called
1522 * <code>RightHandSide</code>, and then a way to generate the coarse
1523 * grid. Since the solution of the previous example program looked like
1524 * curved ridges, we use this name here for the enclosing class. Note that
1525 * the names of the two inner classes have to be the same for all enclosing
1526 * test case classes, and also that we have attached the dimension template
1527 * argument to the enclosing class rather than to the inner ones, to make
1528 * further processing simpler. (From a language viewpoint, a namespace
1529 * would be better to encapsulate these inner classes, rather than a
1530 * structure. However, namespaces cannot be given as template arguments, so
1531 * we use a structure to allow a second object to select from within its
1532 * given argument. The enclosing structure, of course, has no member
1533 * variables apart from the classes it declares, and a static function to
1534 * generate the coarse mesh; it will in general never be instantiated.)
1538 * The idea is then the following (this is the right time to also take a
1539 * brief look at the code below): we can generate objects for boundary
1540 * values and right hand side by simply giving the name of the outer class
1541 * as a template argument to a class which we call here
1542 * <code>Data::SetUp</code>, and it then creates objects for the inner
1543 * classes. In this case, to get all that characterizes the curved ridge
1544 * solution, we would simply generate an instance of
1545 * <code>Data::SetUp@<Data::CurvedRidge@></code>, and everything we need to
1546 * know about the solution would be static member variables and functions of
1551 * This approach might seem like overkill in this case, but will become very
1552 * handy once a certain set up is not only characterized by Dirichlet
1553 * boundary values and a right hand side function, but in addition by
1554 * material properties, Neumann values, different boundary descriptors,
1555 * etc. In that case, the <code>SetUp</code> class might consist of a dozen
1556 * or more objects, and each descriptor class (like the
1557 * <code>CurvedRidges</code> class below) would have to provide them. Then,
1558 * you will be happy to be able to change from one set of data to another by
1559 * only changing the template argument to the <code>SetUp</code> class at
1560 * one place, rather than at many.
1564 * With this framework for different test cases, we are almost finished, but
1565 * one thing remains: by now we can select statically, by changing one
1566 * template argument, which data set to choose. In order to be able to do
1567 * that dynamically, i.e. at run time, we need a base class. This we provide
1568 * in the obvious way, see below, with virtual abstract functions. It forces
1569 * us to introduce a second template parameter <code>dim</code> which we
1570 * need for the base class (which could be avoided using some template
1571 * magic, but we omit that), but that's all.
1576 * classes, only a structure like the <code>CurvedRidges</code> one is
1585 * <a name="step_14-TheSetUpBaseandSetUpclasses"></a>
1586 * <h4>The SetUpBase and SetUp classes</h4>
1590 * Based on the above description, the <code>SetUpBase</code> class then
1591 * looks as follows. To allow using the <code>ObserverPointer</code> class
1592 * with this class, we derived from the
1593 * <code>EnableObserverPointer</code> class.
1596 * template <int dim>
1597 * struct SetUpBase : public EnableObserverPointer
1599 * virtual const Function<dim> &get_boundary_values() const = 0;
1601 * virtual const Function<dim> &get_right_hand_side() const = 0;
1604 * create_coarse_grid(Triangulation<dim> &coarse_grid) const = 0;
1610 * And now for the derived class that takes the template argument as
1615 * Here we pack the data elements into private variables, and allow access
1616 * to them through the methods of the base class.
1619 * template <class Traits, int dim>
1620 * struct SetUp : public SetUpBase<dim>
1622 * virtual const Function<dim> &get_boundary_values() const override;
1624 * virtual const Function<dim> &get_right_hand_side() const override;
1628 * create_coarse_grid(Triangulation<dim> &coarse_grid) const override;
1631 * static const typename Traits::BoundaryValues boundary_values;
1632 * static const typename Traits::RightHandSide right_hand_side;
1637 * We have to provide definitions for the static member variables of the
1641 * template <class Traits, int dim>
1642 * const typename Traits::BoundaryValues SetUp<Traits, dim>::boundary_values;
1643 * template <class Traits, int dim>
1644 * const typename Traits::RightHandSide SetUp<Traits, dim>::right_hand_side;
1648 * And definitions of the member functions:
1651 * template <class Traits, int dim>
1652 * const Function<dim> &SetUp<Traits, dim>::get_boundary_values() const
1654 * return boundary_values;
1658 * template <class Traits, int dim>
1659 * const Function<dim> &SetUp<Traits, dim>::get_right_hand_side() const
1661 * return right_hand_side;
1665 * template <class Traits, int dim>
1666 * void SetUp<Traits, dim>::create_coarse_grid(
1667 * Triangulation<dim> &coarse_grid) const
1669 * Traits::create_coarse_grid(coarse_grid);
1676 * <a name="step_14-TheCurvedRidgesclass"></a>
1677 * <h4>The CurvedRidges class</h4>
1681 * The class that is used to describe the boundary values and right hand
1682 * side of the <code>curved ridge</code> problem already used in the
1683 * @ref step_13 "step-13" example program is then like so:
1686 * template <int dim>
1687 * struct CurvedRidges
1689 * class BoundaryValues : public Function<dim>
1692 * virtual double value(const Point<dim> &p,
1693 * const unsigned int component) const;
1697 * class RightHandSide : public Function<dim>
1700 * virtual double value(const Point<dim> &p,
1701 * const unsigned int component) const;
1704 * static void create_coarse_grid(Triangulation<dim> &coarse_grid);
1708 * template <int dim>
1709 * double CurvedRidges<dim>::BoundaryValues::value(
1710 * const Point<dim> &p,
1711 * const unsigned int /*component*/) const
1714 * for (unsigned int i = 1; i < dim; ++i)
1715 * q += std::sin(10 * p(i) + 5 * p(0) * p(0));
1716 * const double exponential = std::exp(q);
1717 * return exponential;
1722 * template <int dim>
1723 * double CurvedRidges<dim>::RightHandSide::value(
1724 * const Point<dim> &p,
1725 * const unsigned int /*component*/) const
1728 * for (unsigned int i = 1; i < dim; ++i)
1729 * q += std::sin(10 * p(i) + 5 * p(0) * p(0));
1730 * const double u = std::exp(q);
1731 * double t1 = 1, t2 = 0, t3 = 0;
1732 * for (unsigned int i = 1; i < dim; ++i)
1734 * t1 += std::cos(10 * p(i) + 5 * p(0) * p(0)) * 10 * p(0);
1735 * t2 += 10 * std::cos(10 * p(i) + 5 * p(0) * p(0)) -
1736 * 100 * std::sin(10 * p(i) + 5 * p(0) * p(0)) * p(0) * p(0);
1737 * t3 += 100 * std::cos(10 * p(i) + 5 * p(0) * p(0)) *
1738 * std::cos(10 * p(i) + 5 * p(0) * p(0)) -
1739 * 100 * std::sin(10 * p(i) + 5 * p(0) * p(0));
1743 * return -u * (t1 + t2 + t3);
1747 * template <int dim>
1748 * void CurvedRidges<dim>::create_coarse_grid(Triangulation<dim> &coarse_grid)
1750 * GridGenerator::hyper_cube(coarse_grid, -1, 1);
1751 * coarse_grid.refine_global(2);
1758 * <a name="step_14-TheExercise_2_3class"></a>
1759 * <h4>The Exercise_2_3 class</h4>
1763 * This example program was written while giving practical courses for a
1764 * lecture on adaptive finite element methods and duality based error
1765 * estimates. For these courses, we had one exercise, which required to
1766 * solve the Laplace equation with constant right hand side on a square
1767 * domain with a square hole in the center, and zero boundary
1768 * values. Since the implementation of the properties of this problem is
1769 * so particularly simple here, lets do it. As the number of the exercise
1770 * was 2.3, we take the liberty to retain this name for the class as well.
1773 * template <int dim>
1774 * struct Exercise_2_3
1778 * We need a class to denote the boundary values of the problem. In this
1779 * case, this is simple: it's
the zero function,
so don't even declare a
1780 * class, just an alias:
1783 * using BoundaryValues = Functions::ZeroFunction<dim>;
1787 * Second, a class that denotes the right hand side. Since they are
1788 * constant, just subclass the corresponding class of the library and be
1792 * class RightHandSide : public Functions::ConstantFunction<dim>
1796 * : Functions::ConstantFunction<dim>(1.)
1802 * Finally a function to generate the coarse grid. This is somewhat more
1803 * complicated here, see immediately below.
1806 * static void create_coarse_grid(Triangulation<dim> &coarse_grid);
1812 * As stated above, the grid for this example is the square [-1,1]^2 with
1813 * the square [-1/2,1/2]^2 as hole in it. We create the coarse grid as 4
1814 * times 4 cells with the middle four ones missing. To understand how
1815 * exactly the mesh is going to look, it may be simplest to just look
1816 * at the "Results" section of this tutorial program first. In general,
1827 * linker that this function is not implemented for 3d, and needs to be
1832 * For the creation of this geometry, the library has no predefined
1833 * method. In this case, the geometry is still simple enough to do the
1834 * creation by hand, rather than using a mesh generator.
1838 * void Exercise_2_3<2>::create_coarse_grid(Triangulation<2> &coarse_grid)
1842 * We first define the space dimension, to allow those parts of the
1843 * function that are actually dimension independent to use this
1844 * variable. That makes it simpler if you later take this as a starting
1845 * point to implement a 3d version of this mesh. The next step is then
1846 * to have a list of vertices. Here, they are 24 (5 times 5, with the
1847 * middle one omitted). It is probably best to draw a sketch here.
1850 * const unsigned int dim = 2;
1852 * const std::vector<Point<2>> vertices = {
1853 * {-1.0, -1.0}, {-0.5, -1.0}, {+0.0, -1.0}, {+0.5, -1.0}, {+1.0, -1.0},
1854 * {-1.0, -0.5}, {-0.5, -0.5}, {+0.0, -0.5}, {+0.5, -0.5}, {+1.0, -0.5},
1855 * {-1.0, +0.0}, {-0.5, +0.0}, {+0.5, +0.0}, {+1.0, +0.0},
1856 * {-1.0, +0.5}, {-0.5, +0.5}, {+0.0, +0.5}, {+0.5, +0.5}, {+1.0, +0.5},
1857 * {-1.0, +1.0}, {-0.5, +1.0}, {+0.0, +1.0}, {+0.5, +1.0}, {+1.0, +1.0}};
1861 * Next, we have to define the cells and the vertices they contain.
1864 * const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
1865 * cell_vertices = {{{0, 1, 5, 6}},
1871 * {{10, 11, 14, 15}},
1872 * {{12, 13, 17, 18}},
1873 * {{14, 15, 19, 20}},
1874 * {{15, 16, 20, 21}},
1875 * {{16, 17, 21, 22}},
1876 * {{17, 18, 22, 23}}};
1878 * const unsigned int n_cells = cell_vertices.size();
1882 * Again, we generate a C++ vector type from this, but this time by
1883 * looping over the cells (yes, this is boring). Additionally, we set
1884 * the material indicator to zero for all the cells:
1887 * std::vector<CellData<dim>> cells(n_cells, CellData<dim>());
1888 * for (unsigned int i = 0; i < n_cells; ++i)
1890 * for (unsigned int j = 0; j < cell_vertices[i].size(); ++j)
1891 * cells[i].vertices[j] = cell_vertices[i][j];
1892 * cells[i].material_id = 0;
1897 * Finally pass all this information to the library to generate a
1898 * triangulation. The right call for this is
1899 * Triangulation::create_triangulation(), but that function wants
1900 * its input in a format in which cells are consistently oriented
1901 * in some way. It turns out that the mesh we describe with the
1902 * `vertices` and `cells` objects above already is consistently
1903 * oriented, but if you modify the code in some way it may not
1904 * be any more, and so it is good practice to call a function
1905 * that ensures it is -- GridTools::consistently_order_cells()
1910 * The last parameter of the call to Triangulation::create_triangulation()
1911 * below describes what we want to do about boundary and manifold
1912 * indicators on boundary faces. Here, we don't
want to do anything
1919 *
coarse_grid.create_triangulation(vertices, cells,
SubCellData());
1927 *
coarse_grid.refine_global(1);
2016 * <a name=
"step_14-TheDualFunctionalBaseclass"></a>
2027 *
template <
int dim>
2039 * <a name=
"step_14-ThedualfunctionalPointValueEvaluationclass"></a>
2046 * assume to be a vertex. Apart from the constructor that takes and stores
2047 * the evaluation point, this class consists only of the function that
2048 * implements assembling the right hand side.
2051 * template <int dim>
2052 * class PointValueEvaluation : public DualFunctionalBase<dim>
2055 * PointValueEvaluation(const Point<dim> &evaluation_point);
2057 * virtual void assemble_rhs(const DoFHandler<dim> &dof_handler,
2058 * Vector<double> &rhs) const override;
2061 * ExcEvaluationPointNotFound,
2063 * << "The evaluation point " << arg1
2064 * << " was not found among the vertices of the present grid.");
2067 * const Point<dim> evaluation_point;
2071 * template <int dim>
2072 * PointValueEvaluation<dim>::PointValueEvaluation(
2073 * const Point<dim> &evaluation_point)
2074 * : evaluation_point(evaluation_point)
2080 * As for doing the main purpose of the class, assembling the right hand
2081 * side, let us first consider what is necessary: The right hand side of
2082 * the dual problem is a vector of values J(phi_i), where J is the error
2083 * functional, and phi_i is the i-th shape function. Here, J is the
2084 * evaluation at the point x0, i.e. J(phi_i)=phi_i(x0).
2088 * Now, we have assumed that the evaluation point is a vertex. Thus, for
2089 * the usual finite elements we might be using in this program, we can
2090 * take for granted that at such a point exactly one shape function is
2091 * nonzero, and in particular has the value one. Thus, we set the right
2092 * hand side vector to all-zeros, then seek for the shape function
2093 * associated with that point and set the corresponding value of the right
2094 * hand side vector to one:
2097 * template <int dim>
2099 * PointValueEvaluation<dim>::assemble_rhs(const DoFHandler<dim> &dof_handler,
2100 * Vector<double> &rhs) const
2104 * So, first set everything to zeros...
2107 * rhs.reinit(dof_handler.n_dofs());
2111 * ...then loop over cells and find the evaluation point among the
2112 * vertices (or very close to a vertex, which may happen due to floating
2116 * for (const auto &cell : dof_handler.active_cell_iterators())
2117 * for (const auto vertex : cell->vertex_indices())
2118 * if (cell->vertex(vertex).distance(evaluation_point) <
2119 * cell->diameter() * 1e-8)
2123 * Ok, found, so set corresponding entry, and leave function
2124 * since we are finished:
2127 * rhs(cell->vertex_dof_index(vertex, 0)) = 1;
2133 * Finally, a sanity check: if we somehow got here, then we must have
2134 * missed the evaluation point, so raise an exception unconditionally:
2137 * AssertThrow(false, ExcEvaluationPointNotFound(evaluation_point));
2144 * <a name="step_14-ThedualfunctionalPointXDerivativeEvaluationclass"></a>
2145 * <h4>The dual functional PointXDerivativeEvaluation class</h4>
2149 * As second application, we again consider the evaluation of the
2150 * x-derivative of the solution at one point. Again, the declaration of
2151 * the class, and the implementation of its constructor is not too
2155 * template <int dim>
2156 * class PointXDerivativeEvaluation : public DualFunctionalBase<dim>
2159 * PointXDerivativeEvaluation(const Point<dim> &evaluation_point);
2161 * virtual void assemble_rhs(const DoFHandler<dim> &dof_handler,
2162 * Vector<double> &rhs) const;
2165 * ExcEvaluationPointNotFound,
2167 * << "The evaluation point " << arg1
2168 * << " was not found among the vertices of the present grid.");
2171 * const Point<dim> evaluation_point;
2175 * template <int dim>
2176 * PointXDerivativeEvaluation<dim>::PointXDerivativeEvaluation(
2177 * const Point<dim> &evaluation_point)
2178 * : evaluation_point(evaluation_point)
2184 * What is interesting is the implementation of this functional: here,
2185 * J(phi_i)=d/dx phi_i(x0).
2189 * We could, as in the implementation of the respective evaluation object
2190 * take the average of the gradients of each shape function phi_i at this
2191 * evaluation point. However, we take a slightly different approach: we
2192 * simply take the average over all cells that surround this point. The
2193 * question which cells <code>surrounds</code> the evaluation point is
2194 * made dependent on the mesh width by including those cells for which the
2196 *
the cell
's diameter.
2200 * Taking the average of the gradient over the area/volume of these cells
2201 * leads to a dual solution which is very close to the one which would
2202 * result from the point evaluation of the gradient. It is simple to
2203 * justify theoretically that this does not change the method
2207 * template <int dim>
2208 * void PointXDerivativeEvaluation<dim>::assemble_rhs(
2209 * const DoFHandler<dim> &dof_handler,
2210 * Vector<double> &rhs) const
2214 * Again, first set all entries to zero:
2217 * rhs.reinit(dof_handler.n_dofs());
2221 * Initialize a <code>FEValues</code> object with a quadrature formula,
2222 * have abbreviations for the number of quadrature points and shape
2226 * const QGauss<dim> quadrature(dof_handler.get_fe().degree + 1);
2227 * FEValues<dim> fe_values(dof_handler.get_fe(),
2229 * update_gradients | update_quadrature_points |
2230 * update_JxW_values);
2231 * const unsigned int n_q_points = fe_values.n_quadrature_points;
2232 * const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
2236 * ...and have two objects that are used to store the global indices of
2237 * the degrees of freedom on a cell, and the values of the gradients of
2238 * the shape functions at the quadrature points:
2241 * Vector<double> cell_rhs(dofs_per_cell);
2242 * std::vector<unsigned int> local_dof_indices(dofs_per_cell);
2246 * Finally have a variable in which we will sum up the area/volume of
2247 * the cells over which we integrate, by integrating the unit functions
2251 * double total_volume = 0;
2255 * Then start the loop over all cells, and select those cells which are
2256 * close enough to the evaluation point:
2259 * for (const auto &cell : dof_handler.active_cell_iterators())
2260 * if (cell->center().distance(evaluation_point) <= cell->diameter())
2264 * If we have found such a cell, then initialize the
2265 * <code>FEValues</code> object and integrate the x-component of
2266 * the gradient of each shape function, as well as the unit
2267 * function for the total area/volume.
2270 * fe_values.reinit(cell);
2273 * for (unsigned int q = 0; q < n_q_points; ++q)
2275 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2277 * fe_values.shape_grad(i, q)[0] // (d/dx phi_i(x_q))
2278 * * fe_values.JxW(q); // * dx
2279 * total_volume += fe_values.JxW(q);
2284 * If we have the local contributions, distribute them to the
2288 * cell->get_dof_indices(local_dof_indices);
2289 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2290 * rhs(local_dof_indices[i]) += cell_rhs(i);
2295 * After we have looped over all cells, check whether we have found any
2296 * at all, by making sure that their volume is non-zero. If not, then
2297 * the results will be botched, as the right hand side should then still
2298 * be zero, so throw an exception:
2301 * AssertThrow(total_volume > 0,
2302 * ExcEvaluationPointNotFound(evaluation_point));
2306 * Finally, we have by now only integrated the gradients of the shape
2307 * functions, not taking their mean value. We fix this by dividing by
2308 * the measure of the volume over which we have integrated:
2311 * rhs /= total_volume;
2315 * } // namespace DualFunctional
2321 * <a name="step_14-ExtendingtheLaplaceSolvernamespace"></a>
2322 * <h3>Extending the LaplaceSolver namespace</h3>
2325 * namespace LaplaceSolver
2330 * <a name="step_14-TheDualSolverclass"></a>
2331 * <h4>The DualSolver class</h4>
2335 * In the same way as the <code>PrimalSolver</code> class above, we now
2336 * implement a <code>DualSolver</code>. It has all the same features, the
2337 * only difference is that it does not take a function object denoting a
2338 * right hand side object, but now takes a <code>DualFunctionalBase</code>
2339 * object that will assemble the right hand side vector of the dual
2340 * problem. The rest of the class is rather trivial.
2344 * Since both primal and dual solver will use the same triangulation, but
2345 * different discretizations, it now becomes clear why we have made the
2346 * <code>Base</code> class a virtual one: since the final class will be
2347 * derived from both <code>PrimalSolver</code> as well as
2348 * <code>DualSolver</code>, it would have two <code>Base</code> instances,
2349 * would we not have marked the inheritance as virtual. Since in many
2350 * applications the base class would store much more information than just
2351 * the triangulation which needs to be shared between primal and dual
2352 * solvers, we do not usually want to use two such base classes.
2355 * template <int dim>
2356 * class DualSolver : public Solver<dim>
2360 * Triangulation<dim> &triangulation,
2361 * const FiniteElement<dim> &fe,
2362 * const Quadrature<dim> &quadrature,
2363 * const Quadrature<dim - 1> &face_quadrature,
2364 * const DualFunctional::DualFunctionalBase<dim> &dual_functional);
2367 * const ObserverPointer<const DualFunctional::DualFunctionalBase<dim>>
2369 * virtual void assemble_rhs(Vector<double> &rhs) const override;
2371 * static const Functions::ZeroFunction<dim> boundary_values;
2374 * template <int dim>
2375 * const Functions::ZeroFunction<dim> DualSolver<dim>::boundary_values;
2377 * template <int dim>
2378 * DualSolver<dim>::DualSolver(
2379 * Triangulation<dim> &triangulation,
2380 * const FiniteElement<dim> &fe,
2381 * const Quadrature<dim> &quadrature,
2382 * const Quadrature<dim - 1> &face_quadrature,
2383 * const DualFunctional::DualFunctionalBase<dim> &dual_functional)
2384 * : Base<dim>(triangulation)
2385 * , Solver<dim>(triangulation,
2390 * , dual_functional(&dual_functional)
2395 * template <int dim>
2396 * void DualSolver<dim>::assemble_rhs(Vector<double> &rhs) const
2398 * dual_functional->assemble_rhs(this->dof_handler, rhs);
2405 * <a name="step_14-TheWeightedResidualclass"></a>
2406 * <h4>The WeightedResidual class</h4>
2410 * Here finally comes the main class of this program, the one that
2411 * implements the dual weighted residual error estimator. It joins the
2412 * primal and dual solver classes to use them for the computation of
2413 * primal and dual solutions, and implements the error representation
2414 * formula for use as error estimate and mesh refinement.
2418 * The first few of the functions of this class are mostly overriders of
2419 * the respective functions of the base class:
2422 * template <int dim>
2423 * class WeightedResidual : public PrimalSolver<dim>, public DualSolver<dim>
2427 * Triangulation<dim> &coarse_grid,
2428 * const FiniteElement<dim> &primal_fe,
2429 * const FiniteElement<dim> &dual_fe,
2430 * const Quadrature<dim> &quadrature,
2431 * const Quadrature<dim - 1> &face_quadrature,
2432 * const Function<dim> &rhs_function,
2433 * const Function<dim> &boundary_values,
2434 * const DualFunctional::DualFunctionalBase<dim> &dual_functional);
2436 * virtual void solve_problem() override;
2438 * virtual void postprocess(
2439 * const Evaluation::EvaluationBase<dim> &postprocessor) const override;
2441 * virtual unsigned int n_dofs() const override;
2443 * virtual void refine_grid() override;
2445 * virtual void output_solution() const override;
2450 * In the private section, we have two functions that are used to call
2451 * the <code>solve_problem</code> functions of the primal and dual base
2452 * classes. These two functions will be called in parallel by the
2453 * <code>solve_problem</code> function of this class.
2456 * void solve_primal_problem();
2457 * void solve_dual_problem();
2460 * Then declare abbreviations for active cell iterators, to avoid that
2461 * we have to write this lengthy name over and over again:
2467 * using active_cell_iterator =
2468 * typename DoFHandler<dim>::active_cell_iterator;
2472 * Next, declare a data type that we will us to store the contribution
2473 * of faces to the error estimator. The idea is that we can compute the
2474 * face terms from each of the two cells to this face, as they are the
2475 * same when viewed from both sides. What we will do is to compute them
2476 * only once, based on some rules explained below which of the two
2477 * adjacent cells will be in charge to do so. We then store the
2478 * contribution of each face in a map mapping faces to their values, and
2479 * only collect the contributions for each cell by looping over the
2480 * cells a second time and grabbing the values from the map.
2484 * The data type of this map is declared here:
2487 * using FaceIntegrals =
2488 * typename std::map<typename DoFHandler<dim>::face_iterator, double>;
2492 * In the computation of the error estimates on cells and faces, we need
2493 * a number of helper objects, such as <code>FEValues</code> and
2494 * <code>FEFaceValues</code> functions, but also temporary objects
2495 * storing the values and gradients of primal and dual solutions, for
2496 * example. These fields are needed in the three functions that do the
2497 * integration on cells, and regular and irregular faces, respectively.
2501 * There are three reasonable ways to provide these fields: first, as
2502 * local variables in the function that needs them; second, as member
2503 * variables of this class; third, as arguments passed to that function.
2507 * These three alternatives all have drawbacks: the third that their
2508 * number is not negligible and would make calling these functions a
2509 * lengthy enterprise. The second has the drawback that it disallows
2510 * parallelization, since the threads that will compute the error
2511 * estimate have to have their own copies of these variables each, so
2512 * member variables of the enclosing class will not work. The first
2513 * approach, although straightforward, has a subtle but important
2514 * drawback: we will call these functions over and over again, many
2515 * thousands of times maybe; it now turns out that allocating
2516 * vectors and other objects that need memory from the heap is an
2517 * expensive business in terms of run-time, since memory allocation is
2518 * expensive when several threads are involved. It is thus
2519 * significantly better to allocate the memory only once, and recycle
2520 * the objects as often as possible.
2524 * What to do? Our answer is to use a variant of the third strategy.
2525 * In fact, this is exactly what the WorkStream concept is supposed to
2526 * do (we have already introduced it above, but see also @ref threads).
2527 * To avoid that we have to give these functions a dozen or so
2528 * arguments, we pack all these variables into two structures, one which
2529 * is used for the computations on cells, the other doing them on the
2530 * faces. Both are then joined into the WeightedResidualScratchData class
2531 * that will serve as the "scratch data" class of the WorkStream concept:
2536 * FEValues<dim> fe_values;
2537 * const ObserverPointer<const Function<dim>> right_hand_side;
2539 * std::vector<double> cell_residual;
2540 * std::vector<double> rhs_values;
2541 * std::vector<double> dual_weights;
2542 * std::vector<double> cell_laplacians;
2543 * CellData(const FiniteElement<dim> &fe,
2544 * const Quadrature<dim> &quadrature,
2545 * const Function<dim> &right_hand_side);
2546 * CellData(const CellData &cell_data);
2551 * FEFaceValues<dim> fe_face_values_cell;
2552 * FEFaceValues<dim> fe_face_values_neighbor;
2553 * FESubfaceValues<dim> fe_subface_values_cell;
2555 * std::vector<double> jump_residual;
2556 * std::vector<double> dual_weights;
2557 * typename std::vector<Tensor<1, dim>> cell_grads;
2558 * typename std::vector<Tensor<1, dim>> neighbor_grads;
2559 * FaceData(const FiniteElement<dim> &fe,
2560 * const Quadrature<dim - 1> &face_quadrature);
2561 * FaceData(const FaceData &face_data);
2564 * struct WeightedResidualScratchData
2566 * WeightedResidualScratchData(
2567 * const FiniteElement<dim> &primal_fe,
2568 * const Quadrature<dim> &primal_quadrature,
2569 * const Quadrature<dim - 1> &primal_face_quadrature,
2570 * const Function<dim> &rhs_function,
2571 * const Vector<double> &primal_solution,
2572 * const Vector<double> &dual_weights);
2574 * WeightedResidualScratchData(
2575 * const WeightedResidualScratchData &scratch_data);
2577 * CellData cell_data;
2578 * FaceData face_data;
2579 * Vector<double> primal_solution;
2580 * Vector<double> dual_weights;
2586 * WorkStream::run generally wants both a scratch object and a copy
2587 * object. Here, for reasons similar to what we had in @ref step_9 "step-9" when
2588 * discussing the computation of an approximation of the gradient, we
2663 *
, cell_residual(quadrature.
size())
2671 *
template <
int dim>
2673 *
: fe_values(cell_data.fe_values.get_fe(),
2674 *
cell_data.fe_values.get_quadrature(),
2686 *
template <
int dim>
2694 *
, fe_face_values_neighbor(fe,
2710 *
template <
int dim>
2713 *
face_data.fe_face_values_cell.get_quadrature(),
2716 *
, fe_face_values_neighbor(
2717 *
face_data.fe_face_values_neighbor.get_fe(),
2718 *
face_data.fe_face_values_neighbor.get_quadrature(),
2722 *
face_data.fe_subface_values_cell.get_fe(),
2723 *
face_data.fe_subface_values_cell.get_quadrature(),
2733 *
template <
int dim>
2748 *
template <
int dim>
2752 *
: cell_data(scratch_data.cell_data)
2753 *
, face_data(scratch_data.face_data)
2760 *
template <
int dim>
2793 *
template <
int dim>
2805 *
template <
int dim>
2811 *
template <
int dim>
2818 *
template <
int dim>
2826 *
template <
int dim>
2841 *
template <
int dim>
2894 *
template <
int dim>
2921 *
data_out.build_patches();
2932 * <a name=
"step_14-Estimatingerrors"></a>
2938 * <a name=
"step_14-Errorestimationdriverfunctions"></a>
2944 * function that drives all this, i.e. calls those functions that actually
2945 * do the work, and finally collects the results.
2948 * template <int dim>
2950 * WeightedResidual<dim>::estimate_error(Vector<float> &error_indicators) const
2954 * The first task in computing the error is to set up vectors that
2955 * denote the primal solution, and the weights (z-z_h)=(z-I_hz), both in
2956 * the finite element space for which we have computed the dual
2957 * solution. For this, we have to interpolate the primal solution to the
2958 * dual finite element space, and to subtract the interpolation of the
2959 * computed dual solution to the primal finite element
2960 * space. Fortunately, the library provides functions for the
2961 * interpolation into larger or smaller finite element spaces, so this
2962 * is mostly obvious.
3039 *
for (
const auto &cell :
3040 *
DualSolver<dim>::dof_handler.active_cell_iterators())
3044 *
auto worker = [
this,
3054 *
std::function<void(const WeightedResidualCopyData &)>();
3063 *
DualSolver<dim>::dof_handler.begin_active(),
3085 *
unsigned int present_cell = 0;
3087 *
DualSolver<dim>::dof_handler.active_cell_iterators())
3089 *
for (
const auto &face : cell->face_iterators())
3097 *
std::cout <<
" Estimated error: "
3108 * <a name=
"step_14-Estimatingonasinglecell"></a>
3118 *
template <
int dim>
3120 *
const active_cell_iterator &cell,
3143 *
scratch_data.primal_solution,
3144 *
scratch_data.dual_weights,
3145 *
scratch_data.cell_data,
3155 *
for (
const auto face_no : cell->face_indices())
3165 *
if (cell->face(
face_no)->at_boundary())
3184 * work
on this face:
3187 *
if ((cell->neighbor(
face_no)->has_children() ==
false) &&
3188 *
(cell->neighbor(
face_no)->level() == cell->level()) &&
3189 *
(cell->neighbor(
face_no)->index() < cell->index()))
3201 *
if (cell->at_boundary(
face_no) ==
false)
3202 *
if (cell->neighbor(
face_no)->level() < cell->level())
3210 *
other side
's cell is neither coarser not finer than this
3211 * cell, then call one function, and if the cell on the other
3212 * side is further refined, then use another function. Note that
3213 * the case that the cell on the other side is coarser cannot
3214 * happen since we have decided above that we handle this case
3215 * when we pass over that other cell.
3218 * if (cell->face(face_no)->has_children() == false)
3219 * integrate_over_regular_face(cell,
3221 * scratch_data.primal_solution,
3222 * scratch_data.dual_weights,
3223 * scratch_data.face_data,
3226 * integrate_over_irregular_face(cell,
3228 * scratch_data.primal_solution,
3229 * scratch_data.dual_weights,
3230 * scratch_data.face_data,
3239 * <a name="step_14-Computingcelltermerrorcontributions"></a>
3240 * <h4>Computing cell term error contributions</h4>
3244 * As for the actual computation of the error contributions, first turn to
3248 * template <int dim>
3249 * void WeightedResidual<dim>::integrate_over_cell(
3250 * const active_cell_iterator &cell,
3251 * const Vector<double> &primal_solution,
3252 * const Vector<double> &dual_weights,
3253 * CellData &cell_data,
3254 * Vector<float> &error_indicators) const
3258 * The tasks to be done are what appears natural from looking at the
3259 * error estimation formula: first get the right hand side and Laplacian
3260 * of the numerical solution at the quadrature points for the cell
3264 * cell_data.fe_values.reinit(cell);
3265 * cell_data.right_hand_side->value_list(
3266 * cell_data.fe_values.get_quadrature_points(), cell_data.rhs_values);
3267 * cell_data.fe_values.get_function_laplacians(primal_solution,
3268 * cell_data.cell_laplacians);
3272 * ...then get the dual weights...
3275 * cell_data.fe_values.get_function_values(dual_weights,
3276 * cell_data.dual_weights);
3280 * ...and finally build the sum over all quadrature points and store it
3281 * with the present cell:
3285 * for (unsigned int p = 0; p < cell_data.fe_values.n_quadrature_points; ++p)
3286 * sum += ((cell_data.rhs_values[p] + cell_data.cell_laplacians[p]) *
3287 * cell_data.dual_weights[p] * cell_data.fe_values.JxW(p));
3288 * error_indicators(cell->active_cell_index()) += sum;
3295 * <a name="step_14-Computingedgetermerrorcontributions1"></a>
3296 * <h4>Computing edge term error contributions -- 1</h4>
3300 * On the other hand, computation of the edge terms for the error estimate
3301 * is not so simple. First, we have to distinguish between faces with and
3302 * without hanging nodes. Because it is the simple case, we first consider
3303 * the case without hanging nodes on a face (let's call
this the `
regular'
3307 * template <int dim>
3308 * void WeightedResidual<dim>::integrate_over_regular_face(
3309 * const active_cell_iterator &cell,
3310 * const unsigned int face_no,
3311 * const Vector<double> &primal_solution,
3312 * const Vector<double> &dual_weights,
3313 * FaceData &face_data,
3314 * FaceIntegrals &face_integrals) const
3316 * const unsigned int n_q_points =
3317 * face_data.fe_face_values_cell.n_quadrature_points;
3321 * The first step is to get the values of the gradients at the
3322 * quadrature points of the finite element field on the present
3323 * cell. For this, initialize the <code>FEFaceValues</code> object
3324 * corresponding to this side of the face, and extract the gradients
3325 * using that object.
3328 * face_data.fe_face_values_cell.reinit(cell, face_no);
3329 * face_data.fe_face_values_cell.get_function_gradients(
3330 * primal_solution, face_data.cell_grads);
3334 * The second step is then to extract the gradients of the finite
3335 * element solution at the quadrature points on the other side of the
3336 * face, i.e. from the neighboring cell.
3340 * For this, do a sanity check before: make sure that the neighbor
3341 * actually exists (yes, we should not have come here if the neighbor
3342 * did not exist, but in complicated software there are bugs, so better
3343 * check this), and if this is not the case throw an error.
3346 * Assert(cell->neighbor(face_no).state() == IteratorState::valid,
3347 * ExcInternalError());
3350 * If we have that, then we need to find out with which face of the
3351 * neighboring cell we have to work, i.e. the <code>how-many'th</
code>
the
3358 *
cell->neighbor_of_neighbor(
face_no);
3366 *
const active_cell_iterator neighbor = cell->neighbor(
face_no);
3368 *
face_data.fe_face_values_neighbor.get_function_gradients(
3378 *
for (
unsigned int p = 0; p < n_q_points; ++p)
3379 *
face_data.jump_residual[p] =
3380 *
((face_data.cell_grads[p] - face_data.neighbor_grads[p]) *
3381 *
face_data.fe_face_values_cell.normal_vector(p));
3388 *
face_data.fe_face_values_cell.get_function_values(
dual_weights,
3389 *
face_data.dual_weights);
3394 * weights,
and quadrature weights,
to get
the result for this face:
3398 *
for (
unsigned int p = 0; p < n_q_points; ++p)
3400 *
(face_data.jump_residual[p] * face_data.dual_weights[p] *
3401 *
face_data.fe_face_values_cell.JxW(p));
3432 * <a name=
"step_14-Computingedgetermerrorcontributions2"></a>
3441 *
template <
int dim>
3443 *
const active_cell_iterator &cell,
3457 *
const unsigned int n_q_points =
3458 *
face_data.fe_face_values_cell.n_quadrature_points;
3476 *
cell->neighbor_of_neighbor(
face_no);
3481 * face
for all
the sub-faces now:
3514 *
face_data.fe_subface_values_cell.get_function_gradients(
3523 *
face_data.fe_face_values_neighbor.get_function_gradients(
3533 *
for (
unsigned int p = 0; p < n_q_points; ++p)
3534 *
face_data.jump_residual[p] =
3535 *
((face_data.neighbor_grads[p] - face_data.cell_grads[p]) *
3536 *
face_data.fe_face_values_neighbor.normal_vector(p));
3543 *
face_data.fe_face_values_neighbor.get_function_values(
3553 *
for (
unsigned int p = 0; p < n_q_points; ++p)
3555 *
(face_data.jump_residual[p] * face_data.dual_weights[p] *
3556 *
face_data.fe_face_values_neighbor.JxW(p));
3596 * <a name=
"step_14-Asimulationframework"></a>
3620 *
template <
int dim>
3692 *
std::unique_ptr<const DualFunctional::DualFunctionalBase<dim>>
3750 *
template <
int dim>
3765 *
template <
int dim>
3784 *
const QGauss<dim> quadrature(descriptor.dual_fe_degree + 1);
3785 *
const QGauss<dim - 1> face_quadrature(descriptor.dual_fe_degree + 1);
3793 *
std::unique_ptr<LaplaceSolver::Base<dim>> solver;
3794 *
switch (descriptor.refinement_criterion)
3798 *
solver = std::make_unique<LaplaceSolver::WeightedResidual<dim>>(
3804 *
descriptor.data->get_right_hand_side(),
3805 *
descriptor.data->get_boundary_values(),
3806 *
*descriptor.dual_functional);
3812 *
solver = std::make_unique<LaplaceSolver::RefinementGlobal<dim>>(
3817 *
descriptor.data->get_right_hand_side(),
3818 *
descriptor.data->get_boundary_values());
3824 *
solver = std::make_unique<LaplaceSolver::RefinementKelly<dim>>(
3829 *
descriptor.data->get_right_hand_side(),
3830 *
descriptor.data->get_boundary_values());
3837 *
std::make_unique<LaplaceSolver::RefinementWeightedKelly<dim>>(
3842 *
descriptor.data->get_right_hand_side(),
3843 *
descriptor.data->get_boundary_values(),
3844 *
*descriptor.kelly_weight);
3865 *
for (
unsigned int step = 0;
true; ++step)
3867 *
std::cout <<
"Refinement cycle: " << step << std::endl;
3869 *
solver->set_refinement_cycle(step);
3870 *
solver->solve_problem();
3871 *
solver->output_solution();
3873 *
std::cout <<
" Number of degrees of freedom: " << solver->n_dofs()
3883 *
if (solver->n_dofs() < descriptor.max_degrees_of_freedom)
3884 *
solver->refine_grid();
3894 *
std::cout << std::endl;
3904 * <a name=
"step_14-Themainfunction"></a>
3919 *
using namespace Step14;
3927 *
const unsigned int dim = 2;
3935 *
descriptor.refinement_criterion =
3951 *
descriptor.primal_fe_degree = 1;
3952 *
descriptor.dual_fe_degree = 2;
3963 *
std::make_unique<Data::SetUp<Data::Exercise_2_3<dim>, dim>>();
3985 *
descriptor.dual_functional =
3986 *
std::make_unique<DualFunctional::PointValueEvaluation<dim>>(
4001 *
descriptor.max_degrees_of_freedom = 20000;
4017 *
catch (std::exception &exc)
4019 *
std::cerr << std::endl
4021 *
<<
"----------------------------------------------------"
4023 *
std::cerr <<
"Exception on processing: " << std::endl
4024 *
<< exc.what() << std::endl
4025 *
<<
"Aborting!" << std::endl
4026 *
<<
"----------------------------------------------------"
4032 *
std::cerr << std::endl
4034 *
<<
"----------------------------------------------------"
4036 *
std::cerr <<
"Unknown exception!" << std::endl
4037 *
<<
"Aborting!" << std::endl
4038 *
<<
"----------------------------------------------------"
4097First let's look what the program actually computed. On the seventh
4098grid, primal and dual numerical solutions look like this (using a
4099color scheme intended to evoke the snow-capped mountains of
4100Colorado that the original author of this program now calls
4102<table align="center">
4105 <img src="https://www.dealii.org/images/steps/developer/step-14.point-value.solution-7.9.2.png" alt="">
4108 <img src="https://www.dealii.org/images/steps/developer/step-14.point-value.solution-7-dual.9.2.png" alt="">
4112Apparently, the region at the bottom left is so unimportant for the
4113point value evaluation at the top right that the grid is left entirely
4114unrefined there, even though the solution has singularities at the inner
4115corner of that cell! Due
4116to the symmetry in right hand side and domain, the solution should
4117actually look like at the top right in all four corners, but the mesh
4118refinement criterion involving the dual solution chose to refine them
4119differently -- because we said that we really only care about a single
4120function value somewhere at the top right.
4124Here are some of the meshes that are produced in refinement cycles 0,
41252, 4 (top row), and 5, 7, and 8 (bottom row):
4127<table width="80%" align="center">
4129 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-0.9.2.png" alt="" width="100%"></td>
4130 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-2.9.2.png" alt="" width="100%"></td>
4131 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-4.9.2.png" alt="" width="100%"></td>
4134 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-5.9.2.png" alt="" width="100%"></td>
4135 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-7.9.2.png" alt="" width="100%"></td>
4136 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.grid-8.9.2.png" alt="" width="100%"></td>
4140Note the subtle interplay between resolving the corner singularities,
4141and resolving around the point of evaluation. It will be rather
4142difficult to generate such a mesh by hand, as this would involve to
4143judge quantitatively how much which of the four corner singularities
4144should be resolved, and to set the weight compared to the vicinity of
4145the evaluation point.
4149The program prints the point value and the estimated error in this
4150quantity. From extrapolating it, we can guess that the exact value is
4151somewhere close to 0.0334473, plus or minus 0.0000001 (note that we get
4152almost 6 valid digits from only 22,000 (primal) degrees of
4153freedom. This number cannot be obtained from the value of the
4154functional alone, but I have used the assumption that the error
4155estimator is mostly exact, and extrapolated the computed value plus
4156the estimated error, to get an approximation of the true
4157value. Computing with more degrees of freedom shows that this
4158assumption is indeed valid.
4162From the computed results, we can generate two graphs: one that shows
4163the convergence of the error @f$J(u)-J(u_h)@f$ (taking the
4164extrapolated value as correct) in the point value, and the value that
4165we get by adding up computed value @f$J(u_h)@f$ and estimated
4166error eta (if the error estimator @f$eta@f$ were exact, then the value
4167@f$J(u_h)+\eta@f$ would equal the exact point value, and the error
4168in this quantity would always be zero; however, since the error
4169estimator is only a - good - approximation to the true error, we can
4170by this only reduce the size of the error). In this graph, we also
4171indicate the complexity @f${\cal O}(1/N)@f$ to show that mesh refinement
4172acts optimal in this case. The second chart compares
4173true and estimated error, and shows that the two are actually very
4174close to each other, even for such a complicated quantity as the point
4178<table width="80%" align="center">
4180 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.error.png" alt="" width="100%"></td>
4181 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-value.error-estimation.png" alt="" width="100%"></td>
4186<a name="step_14-Comparingrefinementcriteria"></a><h3>Comparing refinement criteria</h3>
4190Since we have accepted quite some effort when using the mesh
4191refinement driven by the dual weighted error estimator (for solving
4192the dual problem, and for evaluating the error representation), it is
4193worth while asking whether that effort was successful. To this end, we
4194first compare the achieved error levels for different mesh refinement
4195criteria. To generate this data, simply change the value of the mesh
4196refinement criterion variable in the main program. The results are
4197thus (for the weight in the Kelly indicator, we have chosen the
4198function @f$1/(r^2+0.1^2)@f$, where @f$r@f$
4199is the distance to the evaluation point; it can be shown that this is
4200the optimal weight if we neglect the effects of boundaries):
4202<img src="https://www.dealii.org/images/steps/developer/step-14.point-value.error-comparison.png" alt="">
4206Checking these numbers, we see that for global refinement, the error
4207is proportional to @f$O(1/(sqrt(N) log(N)))@f$, and for the dual
4208estimator @f$O(1/N)@f$. Generally speaking, we see that the dual
4209weighted error estimator is better than the other refinement
4210indicators, at least when compared with those that have a similarly
4211regular behavior. The Kelly indicator produces smaller errors, but
4212jumps about the picture rather irregularly, with the error also
4213changing signs sometimes. Therefore, its behavior does not allow to
4214extrapolate the results to larger values of N. Furthermore, if we
4215trust the error estimates of the dual weighted error estimator, the
4216results can be improved by adding the estimated error to the computed
4217values. In terms of reliability, the weighted estimator is thus better
4218than the Kelly indicator, although the latter sometimes produces
4223<a name="step_14-Evaluationofpointstresses"></a><h3>Evaluation of point stresses</h3>
4227Besides evaluating the values of the solution at a certain point, the
4228program also offers the possibility to evaluate the x-derivatives at a
4229certain point, and also to tailor mesh refinement for this. To let the
4230program compute these quantities, simply replace the two occurrences of
4231<code>PointValueEvaluation</code> in the main function by
4232<code>PointXDerivativeEvaluation</code>, and let the program run:
4235 Number of degrees of freedom: 72
4236 Point x-derivative: -0.0719397
4237 Estimated error: -0.0126173
4239 Number of degrees of freedom: 61
4240 Point x-derivative: -0.0707956
4241 Estimated error: -0.00774316
4243 Number of degrees of freedom: 131
4244 Point x-derivative: -0.0568671
4245 Estimated error: -0.00313426
4247 Number of degrees of freedom: 247
4248 Point x-derivative: -0.053033
4249 Estimated error: -0.00136114
4251 Number of degrees of freedom: 532
4252 Point x-derivative: -0.0526429
4253 Estimated error: -0.000558868
4255 Number of degrees of freedom: 1267
4256 Point x-derivative: -0.0526955
4257 Estimated error: -0.000220116
4259 Number of degrees of freedom: 2864
4260 Point x-derivative: -0.0527495
4261 Estimated error: -9.46731e-05
4263 Number of degrees of freedom: 6409
4264 Point x-derivative: -0.052785
4265 Estimated error: -4.21543e-05
4267 Number of degrees of freedom: 14183
4268 Point x-derivative: -0.0528028
4269 Estimated error: -2.04241e-05
4271 Number of degrees of freedom: 29902
4272 Point x-derivative: -0.052814
4277The solution looks roughly the same as before (the exact solution of
4278course <em>is</em> the same, only the grid changed a little), but the
4279dual solution is now different. A close-up around the point of
4280evaluation shows this:
4281<table align="center">
4284 <img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.solution-7-dual.png" alt="">
4287 <img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.solution-7-dual-close-up.png" alt="">
4290This time, the grids in refinement cycles 0, 5, 6, 7, 8, and 9 look
4293<table align="center" width="80%">
4295 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-0.9.2.png" alt="" width="100%"></td>
4296 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-5.9.2.png" alt="" width="100%"></td>
4297 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-6.9.2.png" alt="" width="100%"></td>
4300 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-7.9.2.png" alt="" width="100%"></td>
4301 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-8.9.2.png" alt="" width="100%"></td>
4302 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.grid-9.9.2.png" alt="" width="100%"></td>
4306Note the asymmetry of the grids compared with those we obtained for
4307the point evaluation. This is due to the fact that the domain and the primal
4308solution may be symmetric about the diagonal, but the @f$x@f$-derivative is
4309not, and the latter enters the refinement criterion.
4313Then, it is interesting to compare actually computed values of the
4314quantity of interest (i.e. the x-derivative of the solution at one
4315point) with a reference value of -0.0528223... plus or minus
43160.0000005. We get this reference value by computing on finer grid after
4317some more mesh refinements, with approximately 130,000 cells.
4318Recall that if the error is @f$O(1/N)@f$ in the optimal case, then
4319taking a mesh with ten times more cells gives us one additional digit
4324In the left part of the following chart, you again see the convergence
4325of the error towards this extrapolated value, while on the right you
4326see a comparison of true and estimated error:
4328<table width="80%" align="center">
4330 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.error.png" alt="" width="100%"></td>
4331 <td><img src="https://www.dealii.org/images/steps/developer/step-14.point-derivative.error-estimation.png" alt="" width="100%"></td>
4335After an initial phase where the true error changes its sign, the
4336estimated error matches it quite well, again. Also note the dramatic
4337improvement in the error when using the estimated error to correct the
4338computed value of @f$J(u_h)@f$.
4342<a name="step_14-step13revisited"></a><h3>step-13 revisited</h3>
4346If instead of the <code>Exercise_2_3</code> data set, we choose
4347<code>CurvedRidges</code> in the main function, and choose @f$(0.5,0.5)@f$
4348as the evaluation point, then we can redo the
4349computations of the previous example program, to compare whether the
4350results obtained with the help of the dual weighted error estimator
4351are better than those we had previously.
4355First, the meshes after 9 adaptive refinement cycles obtained with
4356the point evaluation and derivative evaluation refinement
4357criteria, respectively, look like this:
4359<table width="80%" align="center">
4361 <td><img src="https://www.dealii.org/images/steps/developer/step-14.step-13.point-value.png" alt="" width="100%"></td>
4362 <td><img src="https://www.dealii.org/images/steps/developer/step-14.step-13.point-derivative.png" alt="" width="100%"></td>
4366The features of the solution can still be seen in the mesh, but since the
4367solution is smooth, the singularities of the dual solution entirely
4368dominate the mesh refinement criterion, and lead to strongly
4369concentrated meshes. The solution after the seventh refinement step looks
4372<table width="40%" align="center">
4374 <td><img src="https://www.dealii.org/images/steps/developer/step-14.step-13.solution-7.9.2.png" alt="" width="100%"></td>
4378Obviously, the solution is worse at some places, but the mesh
4379refinement process should have taken care that these places are not
4380important for computing the point value.
4385The next point is to compare the new (duality based) mesh refinement
4386criterion with the old ones. These are the results:
4388<img src="https://www.dealii.org/images/steps/developer/step-14.step-13.error-comparison.png" alt="">
4392The results are, well, somewhat mixed. First, the Kelly indicator
4393disqualifies itself by its unsteady behavior, changing the sign of the
4394error several times, and with increasing errors under mesh
4395refinement. The dual weighted error estimator has a monotone decrease
4396in the error, and is better than the weighted Kelly and global
4397refinement, but the margin is not as large as expected. This is, here,
4398due to the fact the global refinement can exploit the regular
4399structure of the meshes around the point of evaluation, which leads to
4400a better order of convergence for the point error. However, if we had
4401a mesh that is not locally rectangular, for example because we had to
4402approximate curved boundaries, or if the coefficients were not
4403constant, then this advantage of globally refinement meshes would
4404vanish, while the good performance of the duality based estimator
4410<a name="step_14-Conclusionsandoutlook"></a><h3>Conclusions and outlook</h3>
4414The results here are not too clearly indicating the superiority of the
4415dual weighted error estimation approach for mesh refinement over other
4416mesh refinement criteria, such as the Kelly indicator. This is due to
4417the relative simplicity of the shown applications. If you are not
4418convinced yet that this approach is indeed superior, you are invited
4419to browse through the literature indicated in the introduction, where
4420plenty of examples are provided where the dual weighted approach can
4421reduce the necessary numerical work by orders of magnitude, making
4422this the only way to compute certain quantities to reasonable
4427Besides the objections you may raise against its use as a mesh
4428refinement criterion, consider that accurate knowledge of the error in
4429the quantity one might want to compute is of great use, since we can
4430stop computations when we are satisfied with the accuracy. Using more
4431traditional approaches, it is very difficult to get accurate estimates
4432for arbitrary quantities, except for, maybe, the error in the energy
4433norm, and we will then have no guarantee that the result we computed
4434satisfies any requirements on its accuracy. Also, as was shown for the
4435evaluation of point values and derivatives, the error estimate can be
4436used to extrapolate the results, yielding much higher accuracy in the
4437quantity we want to know.
4441Leaving these mathematical considerations, we tried to write the
4442program in a modular way, such that implementing another test case, or
4443another evaluation and dual functional is simple. You are encouraged
4444to take the program as a basis for your own experiments, and to play a
4448<a name="step_14-PlainProg"></a>
4449<h1> The plain program</h1>
4450@include "step-14.cc"
void write_svg(const Triangulation< 2, 2 > &tria, std::ostream &out) 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)
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
static constexpr unsigned int dimension
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
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_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void consistently_order_cells(std::vector< CellData< dim > > &cells)
Task< RT > new_task(const std::function< RT()> &function)
std::vector< index_type > data
Expression fabs(const Expression &x)
Expression sign(const Expression &x)
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())
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
void refine_and_coarsen_fixed_fraction(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double top_fraction, const double bottom_fraction, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max(), const VectorTools::NormType norm_type=VectorTools::L1_norm)
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
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.)
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
constexpr unsigned int invalid_unsigned_int
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation