496 *
return (-8. * p(0) * p(1));
502 *
const unsigned int )
const
510 * hessian[2][2] =
sin(
PI * p(0)) *
cos(
PI * p(1)) *
exp(p(2));
527 * normal /= p.
norm();
529 *
return (-
trace(hessian) + 2 * (gradient * normal) +
530 * (hessian * normal) * normal);
537 * <a name=
"ImplementationofthecodeLaplaceBeltramiProblemcodeclass"></a>
538 * <h3>Implementation of the <code>LaplaceBeltramiProblem</code>
class</h3>
542 * The rest of the program is actually quite unspectacular
if you know
543 * @ref step_4
"step-4". Our
first step is to define the constructor, setting the
544 * polynomial degree of the finite element and mapping, and associating the
548 *
template <
int spacedim>
549 * LaplaceBeltramiProblem<spacedim>::LaplaceBeltramiProblem(
550 *
const unsigned degree)
560 * <a name=
"LaplaceBeltramiProblemmake_grid_and_dofs"></a>
561 * <h4>LaplaceBeltramiProblem::make_grid_and_dofs</h4>
565 * The next step is to create the mesh, distribute degrees of freedom, and
566 *
set up the various variables that describe the linear system. All of
567 * these steps are standard with the exception of how to create a mesh that
568 * describes a surface. We could generate a mesh
for the domain we are
569 * interested in, generate a
triangulation using a mesh generator, and read
570 * it in
using the
GridIn class. Or, as we
do here, we generate the mesh
575 * In particular, what we
're going to do is this (enclosed between the set
576 * of braces below): we generate a <code>spacedim</code> dimensional mesh
577 * for the half disk (in 2d) or half ball (in 3d), using the
578 * GridGenerator::half_hyper_ball function. This function sets the boundary
579 * indicators of all faces on the outside of the boundary to zero for the
580 * ones located on the perimeter of the disk/ball, and one on the straight
581 * part that splits the full disk/ball into two halves. The next step is the
582 * main point: The GridGenerator::extract_boundary_mesh function creates a
583 * mesh that consists of those cells that are the faces of the previous mesh,
584 * i.e. it describes the <i>surface</i> cells of the original (volume)
585 * mesh. However, we do not want all faces: only those on the perimeter of
586 * the disk or ball which carry boundary indicator zero; we can select these
587 * cells using a set of boundary indicators that we pass to
588 * GridGenerator::extract_boundary_mesh.
592 * There is one point that needs to be mentioned. In order to refine a
593 * surface mesh appropriately if the manifold is curved (similarly to
594 * refining the faces of cells that are adjacent to a curved boundary), the
595 * triangulation has to have an object attached to it that describes where
596 * new vertices should be located. If you don't attach such a boundary
597 * object, they will be located halfway between existing
vertices;
this is
598 * appropriate
if you have a domain with straight boundaries (
e.g. a
599 * polygon) but not when, as here, the manifold has curvature. So
for things
600 * to work properly, we need to attach a manifold
object to our (surface)
601 *
triangulation, in much the same way as we
've already done in 1d for the
602 * boundary. We create such an object and attach it to the triangulation.
606 * The final step in creating the mesh is to refine it a number of
607 * times. The rest of the function is the same as in previous tutorial
611 * template <int spacedim>
612 * void LaplaceBeltramiProblem<spacedim>::make_grid_and_dofs()
615 * Triangulation<spacedim> volume_mesh;
616 * GridGenerator::half_hyper_ball(volume_mesh);
618 * std::set<types::boundary_id> boundary_ids;
619 * boundary_ids.insert(0);
621 * GridGenerator::extract_boundary_mesh(volume_mesh,
625 * triangulation.set_all_manifold_ids(0);
626 * triangulation.set_manifold(0, SphericalManifold<dim, spacedim>());
628 * triangulation.refine_global(4);
630 * std::cout << "Surface mesh has " << triangulation.n_active_cells()
631 * << " cells." << std::endl;
633 * dof_handler.distribute_dofs(fe);
635 * std::cout << "Surface mesh has " << dof_handler.n_dofs()
636 * << " degrees of freedom." << std::endl;
638 * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
639 * DoFTools::make_sparsity_pattern(dof_handler, dsp);
640 * sparsity_pattern.copy_from(dsp);
642 * system_matrix.reinit(sparsity_pattern);
644 * solution.reinit(dof_handler.n_dofs());
645 * system_rhs.reinit(dof_handler.n_dofs());
652 * <a name="LaplaceBeltramiProblemassemble_system"></a>
653 * <h4>LaplaceBeltramiProblem::assemble_system</h4>
657 * The following is the central function of this program, assembling the
658 * matrix that corresponds to the surface Laplacian (Laplace-Beltrami
659 * operator). Maybe surprisingly, it actually looks exactly the same as for
660 * the regular Laplace operator discussed in, for example, @ref step_4 "step-4". The key
661 * is that the FEValues::shape_grad() function does the magic: It returns
662 * the surface gradient @f$\nabla_K \phi_i(x_q)@f$ of the @f$i@f$th shape function
663 * at the @f$q@f$th quadrature point. The rest then does not need any changes
667 * template <int spacedim>
668 * void LaplaceBeltramiProblem<spacedim>::assemble_system()
673 * const QGauss<dim> quadrature_formula(2 * fe.degree);
674 * FEValues<dim, spacedim> fe_values(mapping,
676 * quadrature_formula,
677 * update_values | update_gradients |
678 * update_quadrature_points |
679 * update_JxW_values);
681 * const unsigned int dofs_per_cell = fe.dofs_per_cell;
682 * const unsigned int n_q_points = quadrature_formula.size();
684 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
685 * Vector<double> cell_rhs(dofs_per_cell);
687 * std::vector<double> rhs_values(n_q_points);
688 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
690 * RightHandSide<spacedim> rhs;
692 * for (const auto &cell : dof_handler.active_cell_iterators())
697 * fe_values.reinit(cell);
699 * rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
701 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
702 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
703 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
704 * cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
705 * fe_values.shape_grad(j, q_point) *
706 * fe_values.JxW(q_point);
708 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
709 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
710 * cell_rhs(i) += fe_values.shape_value(i, q_point) *
711 * rhs_values[q_point] * fe_values.JxW(q_point);
713 * cell->get_dof_indices(local_dof_indices);
714 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
716 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
717 * system_matrix.add(local_dof_indices[i],
718 * local_dof_indices[j],
719 * cell_matrix(i, j));
721 * system_rhs(local_dof_indices[i]) += cell_rhs(i);
725 * std::map<types::global_dof_index, double> boundary_values;
726 * VectorTools::interpolate_boundary_values(
727 * mapping, dof_handler, 0, Solution<spacedim>(), boundary_values);
729 * MatrixTools::apply_boundary_values(
730 * boundary_values, system_matrix, solution, system_rhs, false);
738 * <a name="LaplaceBeltramiProblemsolve"></a>
739 * <h4>LaplaceBeltramiProblem::solve</h4>
743 * The next function is the one that solves the linear system. Here, too, no
744 * changes are necessary:
747 * template <int spacedim>
748 * void LaplaceBeltramiProblem<spacedim>::solve()
750 * SolverControl solver_control(solution.size(), 1e-7 * system_rhs.l2_norm());
751 * SolverCG<Vector<double>> cg(solver_control);
753 * PreconditionSSOR<SparseMatrix<double>> preconditioner;
754 * preconditioner.initialize(system_matrix, 1.2);
756 * cg.solve(system_matrix, solution, system_rhs, preconditioner);
764 * <a name="LaplaceBeltramiProblemoutput_result"></a>
765 * <h4>LaplaceBeltramiProblem::output_result</h4>
769 * This is the function that generates graphical output from the
770 * solution. Most of it is boilerplate code, but there are two points worth
775 * - The DataOut::add_data_vector() function can take two kinds of vectors:
776 * Either vectors that have one value per degree of freedom defined by the
777 * DoFHandler object previously attached via DataOut::attach_dof_handler();
778 * and vectors that have one value for each cell of the triangulation, for
779 * example to output estimated errors for each cell. Typically, the
780 * DataOut class knows to tell these two kinds of vectors apart: there are
781 * almost always more degrees of freedom than cells, so we can
782 * differentiate by the two kinds looking at the length of a vector. We
783 * could do the same here, but only because we got lucky: we use a half
784 * sphere. If we had used the whole sphere as domain and @f$Q_1@f$ elements,
785 * we would have the same number of cells as vertices and consequently the
786 * two kinds of vectors would have the same number of elements. To avoid
787 * the resulting confusion, we have to tell the DataOut::add_data_vector()
788 * function which kind of vector we have: DoF data. This is what the third
789 * argument to the function does.
790 * - The DataOut::build_patches() function can generate output that subdivides
791 * each cell so that visualization programs can resolve curved manifolds
792 * or higher polynomial degree shape functions better. We here subdivide
793 * each element in each coordinate direction as many times as the
794 * polynomial degree of the finite element in use.
797 * template <int spacedim>
798 * void LaplaceBeltramiProblem<spacedim>::output_results() const
800 * DataOut<dim, DoFHandler<dim, spacedim>> data_out;
801 * data_out.attach_dof_handler(dof_handler);
802 * data_out.add_data_vector(
805 * DataOut<dim, DoFHandler<dim, spacedim>>::type_dof_data);
806 * data_out.build_patches(mapping, mapping.get_degree());
808 * const std::string filename =
809 * "solution-" + std::to_string(spacedim) + "d.vtk";
810 * std::ofstream output(filename);
811 * data_out.write_vtk(output);
819 * <a name="LaplaceBeltramiProblemcompute_error"></a>
820 * <h4>LaplaceBeltramiProblem::compute_error</h4>
824 * This is the last piece of functionality: we want to compute the error in
825 * the numerical solution. It is a verbatim copy of the code previously
826 * shown and discussed in @ref step_7 "step-7". As mentioned in the introduction, the
827 * <code>Solution</code> class provides the (tangential) gradient of the
828 * solution. To avoid evaluating the error only a superconvergence points,
829 * we choose a quadrature rule of sufficiently high order.
832 * template <int spacedim>
833 * void LaplaceBeltramiProblem<spacedim>::compute_error() const
835 * Vector<float> difference_per_cell(triangulation.n_active_cells());
836 * VectorTools::integrate_difference(mapping,
839 * Solution<spacedim>(),
840 * difference_per_cell,
841 * QGauss<dim>(2 * fe.degree + 1),
842 * VectorTools::H1_norm);
844 * double h1_error = VectorTools::compute_global_error(triangulation,
845 * difference_per_cell,
846 * VectorTools::H1_norm);
847 * std::cout << "H1 error = " << h1_error << std::endl;
855 * <a name="LaplaceBeltramiProblemrun"></a>
856 * <h4>LaplaceBeltramiProblem::run</h4>
860 * The last function provides the top-level logic. Its contents are
864 * template <int spacedim>
865 * void LaplaceBeltramiProblem<spacedim>::run()
867 * make_grid_and_dofs();
873 * } // namespace Step38
879 * <a name="Themainfunction"></a>
880 * <h3>The main() function</h3>
884 * The remainder of the program is taken up by the <code>main()</code>
885 * function. It follows exactly the general layout first introduced in @ref step_6 "step-6"
886 * and used in all following tutorial programs:
893 * using namespace Step38;
895 * LaplaceBeltramiProblem<3> laplace_beltrami;
896 * laplace_beltrami.run();
898 * catch (std::exception &exc)
900 * std::cerr << std::endl
902 * << "----------------------------------------------------"
904 * std::cerr << "Exception on processing: " << std::endl
905 * << exc.what() << std::endl
906 * << "Aborting!" << std::endl
907 * << "----------------------------------------------------"
913 * std::cerr << std::endl
915 * << "----------------------------------------------------"
917 * std::cerr << "Unknown exception!" << std::endl
918 * << "Aborting!" << std::endl
919 * << "----------------------------------------------------"
927 <a name="Results"></a><h1>Results</h1>
930 When you run the program, the following output should be printed on screen:
933 Surface mesh has 1280 cells.
934 Surface mesh has 5185 degrees of freedom.
939 By playing around with the number of global refinements in the
940 <code>LaplaceBeltrami::make_grid_and_dofs</code> function you increase or decrease mesh
941 refinement. For example, doing one more refinement and only running the 3d surface
942 problem yields the following
946 Surface mesh has 5120 cells.
947 Surface mesh has 20609 degrees of freedom.
948 H1 error = 0.00543481
951 This is what we expect: make the mesh size smaller by a factor of two and the
952 error goes down by a factor of four (remember that we use bi-quadratic
953 elements). The full sequence of errors from one to five refinements looks like
954 this, neatly following the theoretically predicted pattern:
963 Finally, the program produces graphical output that we can visualize. Here is
964 a plot of the results:
966 <img src="https://www.dealii.org/images/steps/developer/step-38.solution-3d.png" alt="">
968 The program also works for 1d curves in 2d, not just 2d surfaces in 3d. You
969 can test this by changing the template argument in <code>main()</code> like
972 LaplaceBeltramiProblem<2> laplace_beltrami;
974 The domain is a curve in 2d, and we can visualize the solution by using the
975 third dimension (and color) to denote the value of the function @f$u(x)@f$. This
976 then looks like so (the white curve is the domain, the colored curve is the
977 solution extruded into the third dimension, clearly showing the change in sign
978 as the curve moves from one quadrant of the domain into the adjacent one):
980 <img src="https://www.dealii.org/images/steps/developer/step-38.solution-2d.png" alt="">
983 <a name="extensions"></a>
984 <a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
987 Computing on surfaces only becomes interesting if the surface is more
988 interesting than just a half sphere. To achieve this, deal.II can read
989 meshes that describe surfaces through the usual GridIn class. Or, in case you
990 have an analytic description, a simple mesh can sometimes be stretched and
991 bent into a shape we are interested in.
993 Let us consider a relatively simple example: we take the half sphere we used
994 before, we stretch it by a factor of 10 in the z-direction, and then we jumble
995 the x- and y-coordinates a bit. Let's show the computational domain and the
996 solution
first before we go into details of the implementation below:
998 <img src=
"https://www.dealii.org/images/steps/developer/step-38.warp-1.png" alt=
"">
1000 <img src=
"https://www.dealii.org/images/steps/developer/step-38.warp-2.png" alt=
"">
1003 function. It needs a way to
transform each individual mesh
point to a
1004 different position. Let us here use the following, rather simple function
1005 (remember: stretch in
one direction, jumble in the other two):
1008 template <
int spacedim>
1009 Point<spacedim> warp(const
Point<spacedim> &p)
1012 q[spacedim-1] *= 10;
1023 If we followed the <code>LaplaceBeltrami::make_grid_and_dofs</code>
function, we would
1024 extract the half spherical surface mesh as before, warp it into the shape we
1025 want, and
refine as often as necessary. This is not quite as simple as we
'd
1026 like here, though: refining requires that we have an appropriate manifold
1027 object attached to the triangulation that describes where new vertices of the
1028 mesh should be located upon refinement. I'm sure it
's possible to describe
1029 this manifold in a not-too-complicated way by simply undoing the
1030 transformation above (yielding the spherical surface again), finding the
1031 location of a new point on the sphere, and then re-warping the result. But I'm
1032 a lazy person, and since doing
this is not really the
point here, let
's just
1033 make our lives a bit easier: we'll
extract the half sphere,
refine it as
1034 often as necessary, get rid of the
object that describes the manifold since we
1035 now no longer need it, and then
finally warp the mesh. With the
function
1036 above,
this would look as follows:
1039 template <
int spacedim>
1040 void LaplaceBeltrami<spacedim>::make_grid_and_dofs()
1048 std::set<types::boundary_id> boundary_ids;
1049 boundary_ids.insert(0);
1054 std::ofstream x(
"x"), y(
"y");
1059 std::cout <<
"Surface mesh has " <<
triangulation.n_active_cells()
1066 Note that the only essential addition is the line marked with
1067 asterisks. It is worth pointing out
one other thing here, though: because we
1068 detach the manifold description from the surface mesh, whenever we use a
1069 mapping
object in the rest of the program, it has no curves boundary
1070 description to go on any more. Rather, it will have to use the implicit,
1071 FlatManifold class that is used on all parts of the domain not
1072 explicitly assigned a different manifold object. Consequently, whether we use
1074 using a bilinear approximation.
1076 All these drawbacks aside, the resulting pictures are still pretty. The only
1077 other differences to what's in @ref step_38 "step-38" is that we changed the right hand side
1078 to @f$f(\mathbf x)=
\sin x_3@f$ and the boundary values (through the
1079 <code>Solution</code>
class) to @f$u(\mathbf x)|_{\partial\Omega}=
\cos x_3@f$. Of
1080 course, we now no longer know the exact solution, so the computation of the
1085 <a name=
"PlainProg"></a>
1086 <h1> The plain program</h1>
1087 @include
"step-38.cc"