step-38.h

Go to the documentation of this file.

) const
 *     {
 *       return (-8. * p(0) * p(1));
 *     }
 *   
 *   
 *     template <>
 *     double RightHandSide<3>::value(const Point<3> &p,
 *                                    const unsigned int /*component*/) const
 *     {
 *       using numbers::PI;
 *   
 *       Tensor<2, 3> hessian;
 *   
 *       hessian[0][0] = -PI * PI * sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2));
 *       hessian[1][1] = -PI * PI * sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2));
 *       hessian[2][2] = sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2));
 *   
 *       hessian[0][1] = -PI * PI * cos(PI * p(0)) * sin(PI * p(1)) * exp(p(2));
 *       hessian[1][0] = -PI * PI * cos(PI * p(0)) * sin(PI * p(1)) * exp(p(2));
 *   
 *       hessian[0][2] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2));
 *       hessian[2][0] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2));
 *   
 *       hessian[1][2] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2));
 *       hessian[2][1] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2));
 *   
 *       Tensor<1, 3> gradient;
 *       gradient[0] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2));
 *       gradient[1] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2));
 *       gradient[2] = sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2));
 *   
 *       Point<3> normal = p;
 *       normal /= p.norm();
 *   
 *       return (-trace(hessian) + 2 * (gradient * normal) +
 *               (hessian * normal) * normal);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="ImplementationofthecodeLaplaceBeltramiProblemcodeclass"></a> 
 * <h3>Implementation of the <code>LaplaceBeltramiProblem</code> class</h3>
 * 
 
 * 
 * The rest of the program is actually quite unspectacular if you know
 * @ref step_4 "step-4". Our first step is to define the constructor, setting the
 * polynomial degree of the finite element and mapping, and associating the
 * DoF handler to the triangulation:
 * 
 * @code
 *     template <int spacedim>
 *     LaplaceBeltramiProblem<spacedim>::LaplaceBeltramiProblem(
 *       const unsigned degree)
 *       : fe(degree)
 *       , dof_handler(triangulation)
 *       , mapping(degree)
 *     {}
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaplaceBeltramiProblemmake_grid_and_dofs"></a> 
 * <h4>LaplaceBeltramiProblem::make_grid_and_dofs</h4>
 * 
 
 * 
 * The next step is to create the mesh, distribute degrees of freedom, and
 * set up the various variables that describe the linear system. All of
 * these steps are standard with the exception of how to create a mesh that
 * describes a surface. We could generate a mesh for the domain we are
 * interested in, generate a triangulation using a mesh generator, and read
 * it in using the GridIn class. Or, as we do here, we generate the mesh
 * using the facilities in the GridGenerator namespace.
 *   
 
 * 
 * In particular, what we're going to do is this (enclosed between the set
 * of braces below): we generate a <code>spacedim</code> dimensional mesh
 * for the half disk (in 2d) or half ball (in 3d), using the
 * GridGenerator::half_hyper_ball function. This function sets the boundary
 * indicators of all faces on the outside of the boundary to zero for the
 * ones located on the perimeter of the disk/ball, and one on the straight
 * part that splits the full disk/ball into two halves. The next step is the
 * main point: The GridGenerator::extract_boundary_mesh function creates a
 * mesh that consists of those cells that are the faces of the previous mesh,
 * i.e. it describes the <i>surface</i> cells of the original (volume)
 * mesh. However, we do not want all faces: only those on the perimeter of
 * the disk or ball which carry boundary indicator zero; we can select these
 * cells using a set of boundary indicators that we pass to
 * GridGenerator::extract_boundary_mesh.
 *   
 
 * 
 * There is one point that needs to be mentioned. In order to refine a
 * surface mesh appropriately if the manifold is curved (similarly to
 * refining the faces of cells that are adjacent to a curved boundary), the
 * triangulation has to have an object attached to it that describes where
 * new vertices should be located. If you don't attach such a boundary
 * object, they will be located halfway between existing vertices; this is
 * appropriate if you have a domain with straight boundaries (e.g. a
 * polygon) but not when, as here, the manifold has curvature. So for things
 * to work properly, we need to attach a manifold object to our (surface)
 * triangulation, in much the same way as we've already done in 1d for the
 * boundary. We create such an object and attach it to the triangulation.
 *   
 
 * 
 * The final step in creating the mesh is to refine it a number of
 * times. The rest of the function is the same as in previous tutorial
 * programs.
 * 
 * @code
 *     template <int spacedim>
 *     void LaplaceBeltramiProblem<spacedim>::make_grid_and_dofs()
 *     {
 *       {
 *         Triangulation<spacedim> volume_mesh;
 *         GridGenerator::half_hyper_ball(volume_mesh);
 *   
 *         const std::set<types::boundary_id> boundary_ids = {0};
 *   
 *         GridGenerator::extract_boundary_mesh(volume_mesh,
 *                                              triangulation,
 *                                              boundary_ids);
 *       }
 *       triangulation.set_all_manifold_ids(0);
 *       triangulation.set_manifold(0, SphericalManifold<dim, spacedim>());
 *   
 *       triangulation.refine_global(4);
 *   
 *       std::cout << "Surface mesh has " << triangulation.n_active_cells()
 *                 << " cells." << std::endl;
 *   
 *       dof_handler.distribute_dofs(fe);
 *   
 *       std::cout << "Surface mesh has " << dof_handler.n_dofs()
 *                 << " degrees of freedom." << std::endl;
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *       sparsity_pattern.copy_from(dsp);
 *   
 *       system_matrix.reinit(sparsity_pattern);
 *   
 *       solution.reinit(dof_handler.n_dofs());
 *       system_rhs.reinit(dof_handler.n_dofs());
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaplaceBeltramiProblemassemble_system"></a> 
 * <h4>LaplaceBeltramiProblem::assemble_system</h4>
 * 
 
 * 
 * The following is the central function of this program, assembling the
 * matrix that corresponds to the surface Laplacian (Laplace-Beltrami
 * operator). Maybe surprisingly, it actually looks exactly the same as for
 * the regular Laplace operator discussed in, for example, @ref step_4 "step-4". The key
 * is that the FEValues::shape_grad() function does the magic: It returns
 * the surface gradient @f$\nabla_K \phi_i(x_q)@f$ of the @f$i@f$th shape function
 * at the @f$q@f$th quadrature point. The rest then does not need any changes
 * either:
 * 
 * @code
 *     template <int spacedim>
 *     void LaplaceBeltramiProblem<spacedim>::assemble_system()
 *     {
 *       system_matrix = 0;
 *       system_rhs    = 0;
 *   
 *       const QGauss<dim>       quadrature_formula(2 * fe.degree);
 *       FEValues<dim, spacedim> fe_values(mapping,
 *                                         fe,
 *                                         quadrature_formula,
 *                                         update_values | update_gradients |
 *                                           update_quadrature_points |
 *                                           update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *   
 *       FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *       Vector<double>     cell_rhs(dofs_per_cell);
 *   
 *       std::vector<double>                  rhs_values(n_q_points);
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       RightHandSide<spacedim> rhs;
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           cell_matrix = 0;
 *           cell_rhs    = 0;
 *   
 *           fe_values.reinit(cell);
 *   
 *           rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
 *   
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *               for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *                 cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
 *                                      fe_values.shape_grad(j, q_point) *
 *                                      fe_values.JxW(q_point);
 *   
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *               cell_rhs(i) += fe_values.shape_value(i, q_point) *
 *                              rhs_values[q_point] * fe_values.JxW(q_point);
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             {
 *               for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                 system_matrix.add(local_dof_indices[i],
 *                                   local_dof_indices[j],
 *                                   cell_matrix(i, j));
 *   
 *               system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *             }
 *         }
 *   
 *       std::map<types::global_dof_index, double> boundary_values;
 *       VectorTools::interpolate_boundary_values(
 *         mapping, dof_handler, 0, Solution<spacedim>(), boundary_values);
 *   
 *       MatrixTools::apply_boundary_values(
 *         boundary_values, system_matrix, solution, system_rhs, false);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaplaceBeltramiProblemsolve"></a> 
 * <h4>LaplaceBeltramiProblem::solve</h4>
 * 
 
 * 
 * The next function is the one that solves the linear system. Here, too, no
 * changes are necessary:
 * 
 * @code
 *     template <int spacedim>
 *     void LaplaceBeltramiProblem<spacedim>::solve()
 *     {
 *       SolverControl solver_control(solution.size(), 1e-7 * system_rhs.l2_norm());
 *       SolverCG<Vector<double>> cg(solver_control);
 *   
 *       PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *       preconditioner.initialize(system_matrix, 1.2);
 *   
 *       cg.solve(system_matrix, solution, system_rhs, preconditioner);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaplaceBeltramiProblemoutput_result"></a> 
 * <h4>LaplaceBeltramiProblem::output_result</h4>
 * 
 
 * 
 * This is the function that generates graphical output from the
 * solution. Most of it is boilerplate code, but there are two points worth
 * pointing out:
 *   
 
 * 
 * - The DataOut::add_data_vector() function can take two kinds of vectors:
 * Either vectors that have one value per degree of freedom defined by the
 * DoFHandler object previously attached via DataOut::attach_dof_handler();
 * and vectors that have one value for each cell of the triangulation, for
 * example to output estimated errors for each cell. Typically, the
 * DataOut class knows to tell these two kinds of vectors apart: there are
 * almost always more degrees of freedom than cells, so we can
 * differentiate by the two kinds looking at the length of a vector. We
 * could do the same here, but only because we got lucky: we use a half
 * sphere. If we had used the whole sphere as domain and @f$Q_1@f$ elements,
 * we would have the same number of cells as vertices and consequently the
 * two kinds of vectors would have the same number of elements. To avoid
 * the resulting confusion, we have to tell the DataOut::add_data_vector()
 * function which kind of vector we have: DoF data. This is what the third
 * argument to the function does.
 * - The DataOut::build_patches() function can generate output that subdivides
 * each cell so that visualization programs can resolve curved manifolds
 * or higher polynomial degree shape functions better. We here subdivide
 * each element in each coordinate direction as many times as the
 * polynomial degree of the finite element in use.
 * 
 * @code
 *     template <int spacedim>
 *     void LaplaceBeltramiProblem<spacedim>::output_results() const
 *     {
 *       DataOut<dim, spacedim> data_out;
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution,
 *                                "solution",
 *                                DataOut<dim, spacedim>::type_dof_data);
 *       data_out.build_patches(mapping, mapping.get_degree());
 *   
 *       const std::string filename =
 *         "solution-" + std::to_string(spacedim) + "d.vtk";
 *       std::ofstream output(filename);
 *       data_out.write_vtk(output);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaplaceBeltramiProblemcompute_error"></a> 
 * <h4>LaplaceBeltramiProblem::compute_error</h4>
 * 
 
 * 
 * This is the last piece of functionality: we want to compute the error in
 * the numerical solution. It is a verbatim copy of the code previously
 * shown and discussed in @ref step_7 "step-7". As mentioned in the introduction, the
 * <code>Solution</code> class provides the (tangential) gradient of the
 * solution. To avoid evaluating the error only a superconvergence points,
 * we choose a quadrature rule of sufficiently high order.
 * 
 * @code
 *     template <int spacedim>
 *     void LaplaceBeltramiProblem<spacedim>::compute_error() const
 *     {
 *       Vector<float> difference_per_cell(triangulation.n_active_cells());
 *       VectorTools::integrate_difference(mapping,
 *                                         dof_handler,
 *                                         solution,
 *                                         Solution<spacedim>(),
 *                                         difference_per_cell,
 *                                         QGauss<dim>(2 * fe.degree + 1),
 *                                         VectorTools::H1_norm);
 *   
 *       double h1_error = VectorTools::compute_global_error(triangulation,
 *                                                           difference_per_cell,
 *                                                           VectorTools::H1_norm);
 *       std::cout << "H1 error = " << h1_error << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaplaceBeltramiProblemrun"></a> 
 * <h4>LaplaceBeltramiProblem::run</h4>
 * 
 
 * 
 * The last function provides the top-level logic. Its contents are
 * self-explanatory:
 * 
 * @code
 *     template <int spacedim>
 *     void LaplaceBeltramiProblem<spacedim>::run()
 *     {
 *       make_grid_and_dofs();
 *       assemble_system();
 *       solve();
 *       output_results();
 *       compute_error();
 *     }
 *   } // namespace Step38
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Themainfunction"></a> 
 * <h3>The main() function</h3>
 * 
 
 * 
 * The remainder of the program is taken up by the <code>main()</code>
 * function. It follows exactly the general layout first introduced in @ref step_6 "step-6"
 * and used in all following tutorial programs:
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         using namespace Step38;
 *   
 *         LaplaceBeltramiProblem<3> laplace_beltrami;
 *         laplace_beltrami.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *   
 *     return 0;
 *   }
 * @endcode
<a name="Results"></a><h1>Results</h1>
 
 
When you run the program, the following output should be printed on screen:
 
@verbatim
Surface mesh has 1280 cells.
Surface mesh has 5185 degrees of freedom.
H1 error = 0.0217136
@endverbatim
 
 
By playing around with the number of global refinements in the
<code>LaplaceBeltrami::make_grid_and_dofs</code> function you increase or decrease mesh
refinement. For example, doing one more refinement and only running the 3d surface
problem yields the following
output:
 
@verbatim
Surface mesh has 5120 cells.
Surface mesh has 20609 degrees of freedom.
H1 error = 0.00543481
@endverbatim
 
This is what we expect: make the mesh size smaller by a factor of two and the
error goes down by a factor of four (remember that we use bi-quadratic
elements). The full sequence of errors from one to five refinements looks like
this, neatly following the theoretically predicted pattern:
@verbatim
0.339438
0.0864385
0.0217136
0.00543481
0.00135913
@endverbatim
 
Finally, the program produces graphical output that we can visualize. Here is
a plot of the results:
 
<img src="https://www.dealii.org/images/steps/developer/step-38.solution-3d.png" alt="">
 
The program also works for 1d curves in 2d, not just 2d surfaces in 3d. You
can test this by changing the template argument in <code>main()</code> like
so:
@code
      LaplaceBeltramiProblem<2> laplace_beltrami;
@endcode
The domain is a curve in 2d, and we can visualize the solution by using the
third dimension (and color) to denote the value of the function @f$u(x)@f$. This
then looks like so (the white curve is the domain, the colored curve is the
solution extruded into the third dimension, clearly showing the change in sign
as the curve moves from one quadrant of the domain into the adjacent one):
 
<img src="https://www.dealii.org/images/steps/developer/step-38.solution-2d.png" alt="">
 
 
<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
Computing on surfaces only becomes interesting if the surface is more
interesting than just a half sphere. To achieve this, deal.II can read
meshes that describe surfaces through the usual GridIn class. Or, in case you
have an analytic description, a simple mesh can sometimes be stretched and
bent into a shape we are interested in.
 
Let us consider a relatively simple example: we take the half sphere we used
before, we stretch it by a factor of 10 in the z-direction, and then we jumble
the x- and y-coordinates a bit. Let's show the computational domain and the
solution first before we go into details of the implementation below:
 
<img src="https://www.dealii.org/images/steps/developer/step-38.warp-1.png" alt="">
 
<img src="https://www.dealii.org/images/steps/developer/step-38.warp-2.png" alt="">
 
The way to produce such a mesh is by using the GridTools::transform()
function. It needs a way to transform each individual mesh point to a
different position. Let us here use the following, rather simple function
(remember: stretch in one direction, jumble in the other two):
 
@code
template <int spacedim>
Point<spacedim> warp(const Point<spacedim> &p)
{
  Point<spacedim> q = p;
  q[spacedim-1] *= 10;
 
  if (spacedim >= 2)
    q[0] += 2*std::sin(q[spacedim-1]);
  if (spacedim >= 3)
    q[1] += 2*std::cos(q[spacedim-1]);
 
  return q;
}
@endcode
 
If we followed the <code>LaplaceBeltrami::make_grid_and_dofs</code> function, we would
extract the half spherical surface mesh as before, warp it into the shape we
want, and refine as often as necessary. This is not quite as simple as we'd
like here, though: refining requires that we have an appropriate manifold
object attached to the triangulation that describes where new vertices of the
mesh should be located upon refinement. I'm sure it's possible to describe
this manifold in a not-too-complicated way by simply undoing the
transformation above (yielding the spherical surface again), finding the
location of a new point on the sphere, and then re-warping the result. But I'm
a lazy person, and since doing this is not really the point here, let's just
make our lives a bit easier: we'll extract the half sphere, refine it as
often as necessary, get rid of the object that describes the manifold since we
now no longer need it, and then finally warp the mesh. With the function
above, this would look as follows:
 
@code
template <int spacedim>
void LaplaceBeltrami<spacedim>::make_grid_and_dofs()
{
  {
    Triangulation<spacedim> volume_mesh;
    GridGenerator::half_hyper_ball(volume_mesh);
 
    volume_mesh.refine_global(4);
 
    const std::set<types::boundary_id> boundary_ids = {0};
 
    GridGenerator::extract_boundary_mesh(volume_mesh, triangulation,
                                         boundary_ids);
    GridTools::transform(&warp<spacedim>, triangulation);       /* ** */
    std::ofstream x("x"), y("y");
    GridOut().write_gnuplot(volume_mesh, x);
    GridOut().write_gnuplot(triangulation, y);
  }
 
  std::cout << "Surface mesh has " << triangulation.n_active_cells()
            << " cells."
            << std::endl;
  ...
}
@endcode
 
Note that the only essential addition is the line marked with
asterisks. It is worth pointing out one other thing here, though: because we
detach the manifold description from the surface mesh, whenever we use a
mapping object in the rest of the program, it has no curves boundary
description to go on any more. Rather, it will have to use the implicit,
FlatManifold class that is used on all parts of the domain not
explicitly assigned a different manifold object. Consequently, whether we use
MappingQ(2), MappingQ(15) or MappingQ1, each cell of our mesh will be mapped
using a bilinear approximation.
 
All these drawbacks aside, the resulting pictures are still pretty. The only
other differences to what's in @ref step_38 "step-38" is that we changed the right hand side
to @f$f(\mathbf x)=\sin x_3@f$ and the boundary values (through the
<code>Solution</code> class) to @f$u(\mathbf x)|_{\partial\Omega}=\cos x_3@f$. Of
course, we now no longer know the exact solution, so the computation of the
error at the end of <code>LaplaceBeltrami::run</code> will yield a meaningless
number.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-38.cc"
*/