-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2019 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-// This tutorial program was contributed by Martin Kronbichler
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2019 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+
+ * This tutorial program was contributed by Martin Kronbichler
+ */
+
+// @sect3{Include files}
+
+// The include files for this tutorial are essentially the same as in step-6.
#include <deal.II/base/timer.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/mapping_q_cache.h>
#include <deal.II/fe/mapping_q_generic.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/vector_tools.h>
+#include <fstream>
+
+// The only new include file is the one for the MappingQCache class.
+#include <deal.II/fe/mapping_q_cache.h>
namespace step65
{
using namespace dealii;
+
+ // @sect3{Analytical solution and coefficient}
+
+ // In this tutorial program, we want to solve the Poisson equation with a
+ // coefficient that jumps along a sphere of radius 0.5 and constant right
+ // hand side of value $f(\mathbf{x}) = -3$. Due to the jump in the
+ // coefficient, the analytical solution must have a kink where the
+ // coefficient switches from one value to the other. To keep things simple,
+ // we select an analytical solution that is quadratic in all components,
+ // i.e., $u(x,y,z) = x^2 + y^2 + z^2$ in the ball of radius 0.5 and
+ // $u(x,y,z) = 0.1(x^2 + y^2 + z^2) + 0.25-0.025$ in the outer part of the
+ // domain. This analytical solution is compatible with the right hand side
+ // in case the coefficient is 0.5 in the inner ball and 5 outside.
template <int dim>
class ExactSolution : public Function<dim>
{
+ // @sect3{The PoissonProblem class}
+ //
+ // The implementation of the Poisson problem is very similar to what we used
+ // in the step-5 tutorial program. The two main differences are that we pass
+ // a mapping object to the various steps in the program in order to switch
+ // between two mapping representations as explained in the introduction, and
+ // the time object that will be used for measuring the run times in the
+ // various cases.
template <int dim>
class PoissonProblem
{
void assemble_system(const Mapping<dim> &mapping);
void solve();
void postprocess(const Mapping<dim> &mapping);
- void output_results(const Mapping<dim> &mapping);
Triangulation<dim> triangulation;
FE_Q<dim> fe;
+ // In the constructor, we set up the timer object to record wall times but
+ // be quiet during the normal execution. We will query it for timing details
+ // in the PoissonProblem::run() function. Furthermore, we select a
+ // relatively high polynomial degree of three.
template <int dim>
PoissonProblem<dim>::PoissonProblem()
: fe(3)
+ // @sect3{Grid creation and initialization of the manifolds}
+ //
+ // This function presents the typical usage of
+ // TransfiniteInterpolationManifold. The first step is to create the desired
+ // grid, which can be done by composition of two grids from
+ // GridGenerator. The inner ball mesh is simple enough: We run
+ // GridGenerator::hyper_cube centered in the origin with radius 0.5 (third
+ // function argument). The second mesh is more interesting and constructed
+ // as follows: We want to have a mesh that is spherical in the interior but
+ // flat on the outer surface. Furthermore, the mesh topology of the inner
+ // ball should be compatible with the outer grid in the sense that their
+ // vertices coincide to allow the two grid to be merged. The grid coming out
+ // of GridGenerator::hyper_shell fulfills the requirements on the inner side
+ // in case it is created with $2d$ coarse cells (6 coarse cells in 3D which
+ // we are going to use) &endash; this is the same number of cells as there
+ // are boundary faces for the ball. For the outer surface, we use the fact
+ // that the 6 faces on the surface of the shell without a manifold attached
+ // would degenerate to the surface of a cube. What we are still missing is
+ // the radius of the outer shell boundary. Since we desire a cube of extent
+ // $[-1, 1]$ and the 6-cell shell puts its 8 outer vertices at the 8
+ // opposing diagonals, we must translate the points $(\pm 1, \pm 1, \pm 1)$
+ // into a radius: Clearly, the radius must be $\sqrt{d}$ in $d$ dimensions,
+ // i.e., $\sqrt{3}$ for the three-dimensional case we want to consider.
+ //
+ // Thus, we have a plan: After creating the inner triangulation for the ball
+ // and the one for the outer shell, we merge those two grids but remove all
+ // manifolds from the resulting triangulation to ensure that we have full
+ // control over manifolds. In particular, we want additional points added on
+ // the boundary during refinement to follow a flat manifold description. To
+ // start the process of adding more appropriate manifold ids, we assign the
+ // manifold id 0 to all mesh entities (cells, faces, lines), which will
+ // later be associated with the TransfiniteInterpolationManifold. Then, we
+ // must identify the faces and lines that are along the sphere of radius 0.5
+ // and mark them with a different manifold id to then assign a
+ // SphericalManifold to those. We will choose the manifold id of 1. Since we
+ // have thrown away all manifolds that pre-existed after calling
+ // GridGenerator::hyper_ball(), we manually go through the cells of the mesh
+ // and all their faces. We have found a face on the sphere if all four
+ // vertices have a radius of 0.5, or, as we write in the program, have
+ // $r^2=0.25$. Note that we call `cell->face(f)->set_all_manifold_ids(1)` to
+ // set the manifold id both on the faces and the surrounding
+ // lines. Furthermore, we want to distinguish the cells inside the ball and
+ // outside the ball by a material id for visualization according to the
+ // picture in the introduction.
template <int dim>
void PoissonProblem<dim>::create_grid()
{
Triangulation<dim> tria_outer, tria_inner;
+ GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.5);
GridGenerator::hyper_shell(
tria_outer, Point<dim>(), 0.5, std::sqrt(dim), 2 * dim);
- GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.5);
GridGenerator::merge_triangulations(tria_inner, tria_outer, triangulation);
triangulation.reset_all_manifolds();
triangulation.set_all_manifold_ids(0);
for (unsigned int v = 0;
v < GeometryInfo<dim - 1>::vertices_per_cell;
++v)
- if (std::abs(cell->face(f)->vertex(v).norm_square() - 0.25) >
- 1e-12)
- face_at_sphere_boundary = false;
+ {
+ if (std::abs(cell->face(f)->vertex(v).norm_square() - 0.25) >
+ 1e-12)
+ face_at_sphere_boundary = false;
+ }
if (face_at_sphere_boundary)
cell->face(f)->set_all_manifold_ids(1);
}
else
cell->set_material_id(0);
}
+
+ // With all cells, faces and lines marked appropriately, we can attach the
+ // Manifold objects to those numbers. The entities with manifold id 1 will
+ // get a spherical manifold, whereas the other entities, which have the
+ // manifold id 0, will be assigned the
+ // TransfiniteInterpolationManifold. As mentioned in the introduction, we
+ // must explicitly initialize the manifold to the mesh by a call to
+ // TransfiniteInterpolationManifold::initialize() in order to pick up the
+ // coarse mesh cells and the manifolds attached to the boundaries of those
+ // cells. Note that the manifolds are allowed to go out of scope, because
+ // the Triangulation object internally calls Manifold::clone() to have a
+ // valid object around.
+ //
+ // With all manifolds attached, we will finally go about and refine the
+ // mesh a few times to create a sufficiently large test case.
triangulation.set_manifold(1, SphericalManifold<dim>());
TransfiniteInterpolationManifold<dim> transfinite_manifold;
transfinite_manifold.initialize(triangulation);
+ // @sect3{Setup of data structures}
+ //
+ // This function is well-known from other tutorials in that it enumerates
+ // the degrees of freedom, creates a constraint object and sets up a sparse
+ // matrix for the linear system. The only thing worth mentioning is the fact
+ // that we pass the mapping to the
+ // VectorTools::interpolate_boundary_values() function to ensure that our
+ // boundary values are evaluated on the high-order mesh used for
+ // assembly. In the present example, it does not really matter because the
+ // outer surfaces are flat, but for curved outer cells this is mandatory.
template <int dim>
void PoissonProblem<dim>::setup_system(const Mapping<dim> &mapping)
{
}
-
+ // @sect3{Assembly of the system matrix and right hand side}
+ //
+ // This function is also well-known from the previous tutorial programs. One
+ // thing to note is that we set the number of quadrature points to the
+ // polynomial degree plus two, not the degree plus one as in most other
+ // tutorials. This is because we expect some extra accuracy as the mapping
+ // also involves a degree one more than the polynomials for the
+ // solution.
+ //
+ // The only somewhat unusual code in the assembly is the way we compute the
+ // cell matrix. Rather than using three nested loop over the quadrature
+ // point index, the row and column of the matrix, we first collect the
+ // derivatives of the shape function, multiplied by the square root of the
+ // product of the coefficient and the integration factor `JxW` in a separate
+ // matrix `partial_matrix`. To compute the cell matrix, we then execute
+ // `cell_matrix = partial_matrix * transpose(partial_matrix)` in the line
+ // `partial_matrix.mTmult(cell_matrix, partial_matrix);`. To understand why
+ // this works, we realize that the matrix-matrix multiplication performs a
+ // summation over the columns of `partial_matrix`. If we denote the
+ // coefficient by $a(\mathbf{x}_q)$, the entries in the temporary matrix are
+ // $\sqrt{\text{det}(J) w_q a(x)} \frac{\partial \varphi_i(\boldsymbol
+ // \xi_q)}{\partial x_k}$. If we take the product of the <i>i</i>th row with
+ // the <i>j</i>th column of that matrix, we compute a nested sum involving
+ // $\sum_q \sum_{k=1}^d \sqrt{\text{det}(J) w_q a(x)} \frac{\partial
+ // \varphi_i(\boldsymbol \xi_q)}{\partial x_k} \sqrt{\text{det}(J) w_q a(x)}
+ // \frac{\partial \varphi_j(\boldsymbol \xi_q)}{\partial x_k} = \sum_q
+ // \sum_{k=1}^d\text{det}(J) w_q a(x)\frac{\partial \varphi_i(\boldsymbol
+ // \xi_q)}{\partial x_k} \frac{\partial \varphi_j(\boldsymbol
+ // \xi_q)}{\partial x_k}$, which is exactly the terms needed for the
+ // Laplacian.
+ //
+ // The reason for choosing this somewhat unusual scheme is due to the heavy
+ // work involved in computing the cell matrix for a relatively high
+ // polynomial degree in 3D. As we want to highlight the cost of the mapping
+ // in this tutorial program, we better do the assembly in an optimized way
+ // in order to not chase bottlenecks that have been solved by the community
+ // already. Matrix-matrix multiplication is one of the best optimized
+ // kernels in the HPC context, and the FullMatrix::mTmult() function will
+ // call into those optimized BLAS functions. If the user has provided a good
+ // BLAS library when configuring deal.II (like OpenBLAS or Intel's MKL), the
+ // computation of the cell matrix will execute close to the processor's peak
+ // arithmetic performance. As a side note, we mention that despite an
+ // optimized matrix-matrix multiplication, the current strategy is
+ // sub-optimal in terms of complexity as the work to be done is proportional
+ // to $(p+1)^9$ operations for degree $p$ (this also applies to the usual
+ // evaluation with FEValues). One could compute the cell matrix with
+ // $\mathcal O((p+1)^7)$ operations by utilizing the tensor product
+ // structure of the shape functions, as is done by the matrix-free framework
+ // in deal.II. We refer to step-37 and the documentation of the
+ // tensor-product-aware evaluators FEEvaluation for details on how an even
+ // more efficient cell matrix computation could be realized.
template <int dim>
void PoissonProblem<dim>::assemble_system(const Mapping<dim> &mapping)
{
TimerOutput::Scope scope(timer, "Assemble linear system");
- QGauss<dim> quadrature_formula(fe.degree + 1);
+ QGauss<dim> quadrature_formula(fe.degree + 2);
FEValues<dim> fe_values(mapping,
fe,
quadrature_formula,
Vector<double> cell_rhs(dofs_per_cell);
FullMatrix<double> partial_matrix(dofs_per_cell, dim * n_q_points);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell_rhs = 0.;
+ // @sect3{Solution of the linear system}
+ //
+ // For solving the linear system, we pick a simple Jacobi-preconditioned
+ // variant, similar to the settings in the early tutorials.
template <int dim>
void PoissonProblem<dim>::solve()
{
+ // @sect3{Output of the solution and computation of errors}
+ //
+ // In this function, we do various post-processing steps with the
+ // solution, all of which involve the mapping in one way or the other.
+ //
+ // The first operation we do is to write the solution as well as the
+ // material ids to a vtu file. This is similar to what was done in many
+ // other tutorial programs. The new ingredient presented in this tutorial
+ // program is the way we ensure that a high-order representation is written
+ // to the file. We need to consider two particular topics. Firstly, we tell
+ // the DataOut object via the DataOutBase::VtkFlags that we intend to
+ // interpret the subdivisions of the elements as a high-order Lagrange
+ // polynomial. Recent visualization programs, like ParaView of version 5.5
+ // or newer, can then render a high-order solution (one typically needs to
+ // adjust a parameter called "nonlinear subdivision level"). Secondly, we
+ // need to make sure that the mapping is passed to the
+ // DataOut::build_patches() method. Finally, the DataOut class only prints
+ // curved boundary cells by default, so we need to ensure that even inner
+ // cells are printed in a curved representation via the mapping.
template <int dim>
void PoissonProblem<dim>::postprocess(const Mapping<dim> &mapping)
{
{
TimerOutput::Scope scope(timer, "Write output");
- Timer time;
- DataOut<dim> data_out;
- data_out.attach_dof_handler(dof_handler);
- data_out.add_data_vector(solution, "solution");
+
+ DataOut<dim> data_out;
DataOutBase::VtkFlags flags;
flags.write_higher_order_cells = true;
data_out.set_flags(flags);
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "solution");
+
Vector<double> material_ids(triangulation.n_active_cells());
for (const auto &cell : triangulation.active_cell_iterators())
material_ids[cell->active_cell_index()] = cell->material_id();
data_out.write_vtu(file);
}
+ // The next operation in the postprocessing function is to compute the L2
+ // and H1 errors against the analytical solution. As the analytical
+ // solution is a quadratic polynomial, we expect a very accurate result at
+ // this point. If we were solving on a simple mesh with affine element
+ // shapes, we would expect the numerical result to coincide with the
+ // analyical solution up to roundoff accuracy. However, since we are used
+ // deformed cells following a sphere, which are only tracked by
+ // polynomials of degree 4 (one more than the degree for the finite
+ // elements), we will make an error around 1e-7. We could get more
+ // accuracy by increasing the polynomial degree or refining the mesh.
{
TimerOutput::Scope scope(timer, "Compute error norms");
Vector<double> norm_per_cell_p(triangulation.n_active_cells());
<< norm_per_cell_p.l2_norm() << std::endl;
}
+ // The final post-processing operation we do here is to compute an error
+ // estimation with the KellyErrorEstimator. We use the exact same settings
+ // as in the step-6 tutorial program, except for the fact that we also
+ // hand in the mapping to ensure that errors are evaluated along the
+ // curved element, consistent with the remainder of the program. However,
+ // we do not really use the result here to drive a mesh adaptation step
+ // (that would refine the mesh around the material interface along the
+ // sphere), as the focus here is on the cost of this operation.
{
TimerOutput::Scope scope(timer, "Compute error estimator");
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell);
+ std::cout << " Max cell-wise error estimate: "
+ << estimated_error_per_cell.linfty_norm() << std::endl;
}
}
- template <int dim>
- void PoissonProblem<dim>::output_results(const Mapping<dim> &mapping)
- {
- TimerOutput::Scope scope(timer, "Write output");
- Timer time;
- DataOut<dim> data_out;
- data_out.attach_dof_handler(dof_handler);
- data_out.add_data_vector(solution, "solution");
-
- DataOutBase::VtkFlags flags;
- flags.write_higher_order_cells = true;
- data_out.set_flags(flags);
-
- Vector<double> material_ids(triangulation.n_active_cells());
- for (const auto &cell : triangulation.active_cell_iterators())
- material_ids[cell->active_cell_index()] = cell->material_id();
- data_out.add_data_vector(material_ids, "material_ids");
-
- data_out.build_patches(mapping,
- fe.degree,
- DataOut<dim>::curved_inner_cells);
-
- std::ofstream file(
- ("solution-" +
- std::to_string(triangulation.n_global_levels() - 10 + 2 * dim) + ".vtu")
- .c_str());
- data_out.write_vtu(file);
- }
-
-
-
+ // @sect3{The PoissonProblem::run() function}
+ //
+ // Finally, we define the `run()` function that controls how we want to
+ // execute this program (which is called by the main() function in the usual
+ // way). We start by calling the `create_grid()` function that sets up our
+ // geometry with the appropriate manifolds. We then run two instances of a
+ // solver chain, starting from the setup of the equations, the assembly of
+ // the linear system, its solution with a simple iterative solver, and the
+ // postprocessing discussed above. The two instances differ in the way they
+ // use the mapping. The first uses a conventional MappingQGeneric mapping
+ // object which we initialize to a degree one more than we use for the
+ // finite element &endash; after all, we expect the geometry representation
+ // to be the bottleneck as the analytic solution is only a quadratic
+ // polynomial. (In reality, things are interlinked to quite some extent
+ // because the evaluation of the polynomials in real coordinates involves
+ // the mapping of a higher-degree polynomials, which represent some smooth
+ // rational functions. As a consequence, higher-degree polynomials still pay
+ // off, so it does not make sense to increase the degree of the mapping
+ // further.) Once the first pass is completed, we let the timer print a
+ // summary of the compute times of the individual stages.
template <int dim>
void PoissonProblem<dim>::run()
{
create_grid();
+
{
std::cout << std::endl
<< "====== Running with the basic MappingQGeneric class ====== "
assemble_system(mapping);
solve();
postprocess(mapping);
- output_results(mapping);
timer.print_summary();
+ timer.reset();
}
- timer.reset();
-
+ // For the second instance, we instead set up the MappingQCache class. Its
+ // use is very simple: After constructing it (with the degree, given that
+ // we want it to show the correct degree functionality in other contexts),
+ // we fill the cache via the MappingQCache::initialize() function. At this
+ // stage, we specify which mapping we want to use (obviously, the same
+ // MappingQGeneric as previously in order to repeat the same computations)
+ // for the cache, and then run through the same functions again, now
+ // handing in the modified mapping. In the end, we again print the
+ // accumulated wall times since the reset to see how the times compare to
+ // the original setting.
{
std::cout
<< "====== Running with the optimized MappingQCache class ====== "
assemble_system(mapping);
solve();
postprocess(mapping);
- output_results(mapping);
timer.print_summary();
}
}
-
} // namespace step65
+
int main()
{
step65::PoissonProblem<3> test_program;
projects / dealii.git / commitdiff