488 *
return (-8. * p[0] * p[1]);
493 *
double RightHandSide<3>::value(
const Point<3> &p,
494 *
const unsigned int )
const
519 * normal /= p.
norm();
521 *
return (-
trace(hessian) + 2 * (gradient * normal) +
522 * (hessian * normal) * normal);
529 * <a name=
"step_38-ImplementationofthecodeLaplaceBeltramiProblemcodeclass"></a>
530 * <h3>Implementation of the <code>LaplaceBeltramiProblem</code>
class</h3>
534 * The rest of the program is actually quite unspectacular
if you know
535 * @ref step_4
"step-4". Our
first step is to define the constructor, setting the
536 * polynomial degree of the finite element and mapping, and associating the
540 *
template <
int spacedim>
541 * LaplaceBeltramiProblem<spacedim>::LaplaceBeltramiProblem(
542 *
const unsigned degree)
552 * <a name=
"step_38-LaplaceBeltramiProblemmake_grid_and_dofs"></a>
553 * <h4>LaplaceBeltramiProblem::make_grid_and_dofs</h4>
557 * The next step is to create the mesh, distribute degrees of freedom, and
558 * set up the various variables that describe the linear system. All of
559 * these steps are standard with the exception of how to create a mesh that
560 * describes a surface. We could generate a mesh
for the domain we are
561 * interested in, generate a
triangulation using a mesh generator, and read
562 * it in
using the
GridIn class. Or, as we
do here, we generate the mesh
567 * In particular, what we
're going to do is this (enclosed between the set
568 * of braces below): we generate a <code>spacedim</code> dimensional mesh
569 * for the half disk (in 2d) or half ball (in 3d), using the
570 * GridGenerator::half_hyper_ball function. This function sets the boundary
571 * indicators of all faces on the outside of the boundary to zero for the
572 * ones located on the perimeter of the disk/ball, and one on the straight
573 * part that splits the full disk/ball into two halves. The next step is the
574 * main point: The GridGenerator::extract_boundary_mesh function creates a
575 * mesh that consists of those cells that are the faces of the previous mesh,
576 * i.e. it describes the <i>surface</i> cells of the original (volume)
577 * mesh. However, we do not want all faces: only those on the perimeter of
578 * the disk or ball which carry boundary indicator zero; we can select these
579 * cells using a set of boundary indicators that we pass to
580 * GridGenerator::extract_boundary_mesh.
584 * There is one point that needs to be mentioned. In order to refine a
585 * surface mesh appropriately if the manifold is curved (similarly to
586 * refining the faces of cells that are adjacent to a curved boundary), the
587 * triangulation has to have an object attached to it that describes where
588 * new vertices should be located. If you don't attach such a boundary
589 * object, they will be located halfway between existing
vertices;
this is
590 * appropriate
if you have a domain with straight boundaries (
e.g. a
591 * polygon) but not when, as here, the manifold has curvature. So
for things
592 * to work properly, we need to attach a manifold
object to our (surface)
593 *
triangulation, in much the same way as we
've already done in 1d for the
594 * boundary. We create such an object and attach it to the triangulation.
598 * The final step in creating the mesh is to refine it a number of
599 * times. The rest of the function is the same as in previous tutorial
603 * template <int spacedim>
604 * void LaplaceBeltramiProblem<spacedim>::make_grid_and_dofs()
607 * Triangulation<spacedim> volume_mesh;
608 * GridGenerator::half_hyper_ball(volume_mesh);
610 * const std::set<types::boundary_id> boundary_ids = {0};
612 * GridGenerator::extract_boundary_mesh(volume_mesh,
616 * triangulation.set_all_manifold_ids(0);
617 * triangulation.set_manifold(0, SphericalManifold<dim, spacedim>());
619 * triangulation.refine_global(4);
621 * std::cout << "Surface mesh has " << triangulation.n_active_cells()
622 * << " cells." << std::endl;
624 * dof_handler.distribute_dofs(fe);
626 * std::cout << "Surface mesh has " << dof_handler.n_dofs()
627 * << " degrees of freedom." << std::endl;
629 * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
630 * DoFTools::make_sparsity_pattern(dof_handler, dsp);
631 * sparsity_pattern.copy_from(dsp);
633 * system_matrix.reinit(sparsity_pattern);
635 * solution.reinit(dof_handler.n_dofs());
636 * system_rhs.reinit(dof_handler.n_dofs());
643 * <a name="step_38-LaplaceBeltramiProblemassemble_system"></a>
644 * <h4>LaplaceBeltramiProblem::assemble_system</h4>
648 * The following is the central function of this program, assembling the
649 * matrix that corresponds to the surface Laplacian (Laplace-Beltrami
650 * operator). Maybe surprisingly, it actually looks exactly the same as for
651 * the regular Laplace operator discussed in, for example, @ref step_4 "step-4". The key
652 * is that the FEValues::shape_grad() function does the magic: It returns
653 * the surface gradient @f$\nabla_K \phi_i(x_q)@f$ of the @f$i@f$th shape function
654 * at the @f$q@f$th quadrature point. The rest then does not need any changes
658 * template <int spacedim>
659 * void LaplaceBeltramiProblem<spacedim>::assemble_system()
664 * const QGauss<dim> quadrature_formula(2 * fe.degree);
665 * FEValues<dim, spacedim> fe_values(mapping,
667 * quadrature_formula,
668 * update_values | update_gradients |
669 * update_quadrature_points |
670 * update_JxW_values);
672 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
673 * const unsigned int n_q_points = quadrature_formula.size();
675 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
676 * Vector<double> cell_rhs(dofs_per_cell);
678 * std::vector<double> rhs_values(n_q_points);
679 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
681 * RightHandSide<spacedim> rhs;
683 * for (const auto &cell : dof_handler.active_cell_iterators())
688 * fe_values.reinit(cell);
690 * rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
692 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
693 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
694 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
695 * cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
696 * fe_values.shape_grad(j, q_point) *
697 * fe_values.JxW(q_point);
699 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
700 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
701 * cell_rhs(i) += fe_values.shape_value(i, q_point) *
702 * rhs_values[q_point] * fe_values.JxW(q_point);
704 * cell->get_dof_indices(local_dof_indices);
705 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
707 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
708 * system_matrix.add(local_dof_indices[i],
709 * local_dof_indices[j],
710 * cell_matrix(i, j));
712 * system_rhs(local_dof_indices[i]) += cell_rhs(i);
716 * std::map<types::global_dof_index, double> boundary_values;
717 * VectorTools::interpolate_boundary_values(
718 * mapping, dof_handler, 0, Solution<spacedim>(), boundary_values);
720 * MatrixTools::apply_boundary_values(
721 * boundary_values, system_matrix, solution, system_rhs, false);
729 * <a name="step_38-LaplaceBeltramiProblemsolve"></a>
730 * <h4>LaplaceBeltramiProblem::solve</h4>
734 * The next function is the one that solves the linear system. Here, too, no
735 * changes are necessary:
738 * template <int spacedim>
739 * void LaplaceBeltramiProblem<spacedim>::solve()
741 * SolverControl solver_control(solution.size(), 1e-7 * system_rhs.l2_norm());
742 * SolverCG<Vector<double>> cg(solver_control);
744 * PreconditionSSOR<SparseMatrix<double>> preconditioner;
745 * preconditioner.initialize(system_matrix, 1.2);
747 * cg.solve(system_matrix, solution, system_rhs, preconditioner);
755 * <a name="step_38-LaplaceBeltramiProblemoutput_result"></a>
756 * <h4>LaplaceBeltramiProblem::output_result</h4>
760 * This is the function that generates graphical output from the
761 * solution. Most of it is boilerplate code, but there are two points worth
766 * - The DataOut::add_data_vector() function can take two kinds of vectors:
767 * Either vectors that have one value per degree of freedom defined by the
768 * DoFHandler object previously attached via DataOut::attach_dof_handler();
769 * and vectors that have one value for each cell of the triangulation, for
770 * example to output estimated errors for each cell. Typically, the
771 * DataOut class knows to tell these two kinds of vectors apart: there are
772 * almost always more degrees of freedom than cells, so we can
773 * differentiate by the two kinds looking at the length of a vector. We
774 * could do the same here, but only because we got lucky: we use a half
775 * sphere. If we had used the whole sphere as domain and @f$Q_1@f$ elements,
776 * we would have the same number of cells as vertices and consequently the
777 * two kinds of vectors would have the same number of elements. To avoid
778 * the resulting confusion, we have to tell the DataOut::add_data_vector()
779 * function which kind of vector we have: DoF data. This is what the third
780 * argument to the function does.
781 * - The DataOut::build_patches() function can generate output that subdivides
782 * each cell so that visualization programs can resolve curved manifolds
783 * or higher polynomial degree shape functions better. We here subdivide
784 * each element in each coordinate direction as many times as the
785 * polynomial degree of the finite element in use.
788 * template <int spacedim>
789 * void LaplaceBeltramiProblem<spacedim>::output_results() const
791 * DataOut<dim, spacedim> data_out;
792 * data_out.attach_dof_handler(dof_handler);
793 * data_out.add_data_vector(solution,
795 * DataOut<dim, spacedim>::type_dof_data);
796 * data_out.build_patches(mapping, mapping.get_degree());
798 * const std::string filename =
799 * "solution-" + std::to_string(spacedim) + "d.vtk";
800 * std::ofstream output(filename);
801 * data_out.write_vtk(output);
809 * <a name="step_38-LaplaceBeltramiProblemcompute_error"></a>
810 * <h4>LaplaceBeltramiProblem::compute_error</h4>
814 * This is the last piece of functionality: we want to compute the error in
815 * the numerical solution. It is a verbatim copy of the code previously
816 * shown and discussed in @ref step_7 "step-7". As mentioned in the introduction, the
817 * <code>Solution</code> class provides the (tangential) gradient of the
818 * solution. To avoid evaluating the error only a superconvergence points,
819 * we choose a quadrature rule of sufficiently high order.
822 * template <int spacedim>
823 * void LaplaceBeltramiProblem<spacedim>::compute_error() const
825 * Vector<float> difference_per_cell(triangulation.n_active_cells());
826 * VectorTools::integrate_difference(mapping,
829 * Solution<spacedim>(),
830 * difference_per_cell,
831 * QGauss<dim>(2 * fe.degree + 1),
832 * VectorTools::H1_norm);
834 * double h1_error = VectorTools::compute_global_error(triangulation,
835 * difference_per_cell,
836 * VectorTools::H1_norm);
837 * std::cout << "H1 error = " << h1_error << std::endl;
845 * <a name="step_38-LaplaceBeltramiProblemrun"></a>
846 * <h4>LaplaceBeltramiProblem::run</h4>
850 * The last function provides the top-level logic. Its contents are
854 * template <int spacedim>
855 * void LaplaceBeltramiProblem<spacedim>::run()
857 * make_grid_and_dofs();
863 * } // namespace Step38
869 * <a name="step_38-Themainfunction"></a>
870 * <h3>The main() function</h3>
874 * The remainder of the program is taken up by the <code>main()</code>
875 * function. It follows exactly the general layout first introduced in @ref step_6 "step-6"
876 * and used in all following tutorial programs:
883 * using namespace Step38;
885 * LaplaceBeltramiProblem<3> laplace_beltrami;
886 * laplace_beltrami.run();
888 * catch (std::exception &exc)
890 * std::cerr << std::endl
892 * << "----------------------------------------------------"
894 * std::cerr << "Exception on processing: " << std::endl
895 * << exc.what() << std::endl
896 * << "Aborting!" << std::endl
897 * << "----------------------------------------------------"
903 * std::cerr << std::endl
905 * << "----------------------------------------------------"
907 * std::cerr << "Unknown exception!" << std::endl
908 * << "Aborting!" << std::endl
909 * << "----------------------------------------------------"
917<a name="step_38-Results"></a><h1>Results</h1>
920When you run the program, the following output should be printed on screen:
923Surface mesh has 1280 cells.
924Surface mesh has 5185 degrees of freedom.
929By playing around with the number of global refinements in the
930<code>LaplaceBeltrami::make_grid_and_dofs</code> function you increase or decrease mesh
931refinement. For example, doing one more refinement and only running the 3d surface
932problem yields the following
936Surface mesh has 5120 cells.
937Surface mesh has 20609 degrees of freedom.
941This is what we expect: make the mesh size smaller by a factor of two and the
942error goes down by a factor of four (remember that we use bi-quadratic
943elements). The full sequence of errors from one to five refinements looks like
944this, neatly following the theoretically predicted pattern:
953Finally, the program produces graphical output that we can visualize. Here is
954a plot of the results:
956<img src="https://www.dealii.org/images/steps/developer/step-38.solution-3d.png" alt="">
958The program also works for 1d curves in 2d, not just 2d surfaces in 3d. You
959can test this by changing the template argument in <code>main()</code> like
962 LaplaceBeltramiProblem<2> laplace_beltrami;
964The domain is a curve in 2d, and we can visualize the solution by using the
965third dimension (and color) to denote the value of the function @f$u(x)@f$. This
966then looks like so (the white curve is the domain, the colored curve is the
967solution extruded into the third dimension, clearly showing the change in sign
968as the curve moves from one quadrant of the domain into the adjacent one):
970<img src="https://www.dealii.org/images/steps/developer/step-38.solution-2d.png" alt="">
973<a name="step-38-extensions"></a>
974<a name="step_38-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
977Computing on surfaces only becomes interesting if the surface is more
978interesting than just a half sphere. To achieve this, deal.II can read
979meshes that describe surfaces through the usual GridIn class. Or, in case you
980have an analytic description, a simple mesh can sometimes be stretched and
981bent into a shape we are interested in.
983Let us consider a relatively simple example: we take the half sphere we used
984before, we stretch it by a factor of 10 in the z-direction, and then we jumble
985the x- and y-coordinates a bit. Let's show the computational domain and the
986solution
first before we go into details of the implementation below:
988<img src=
"https://www.dealii.org/images/steps/developer/step-38.warp-1.png" alt=
"">
990<img src=
"https://www.dealii.org/images/steps/developer/step-38.warp-2.png" alt=
"">
993function. It needs a way to transform each individual mesh point to a
994different position. Let us here use the following, rather simple function
995(remember: stretch in one direction, jumble in the other two):
998template <
int spacedim>
1002 q[spacedim-1] *= 10;
1013If we followed the <code>LaplaceBeltrami::make_grid_and_dofs</code> function, we would
1014extract the half spherical surface mesh as before, warp it into the shape we
1015want, and
refine as often as necessary. This is not quite as simple as we
'd
1016like here, though: refining requires that we have an appropriate manifold
1017object attached to the triangulation that describes where new vertices of the
1018mesh should be located upon refinement. I'm sure it
's possible to describe
1019this manifold in a not-too-complicated way by simply undoing the
1020transformation above (yielding the spherical surface again), finding the
1021location of a new point on the sphere, and then re-warping the result. But I'm
1022a lazy person, and since doing
this is not really the
point here, let
's just
1023make our lives a bit easier: we'll
extract the half sphere,
refine it as
1024often as necessary, get rid of the
object that describes the manifold since we
1025now no longer need it, and then
finally warp the mesh. With the function
1026above,
this would look as follows:
1029template <
int spacedim>
1030void LaplaceBeltrami<spacedim>::make_grid_and_dofs()
1038 const std::set<types::boundary_id> boundary_ids = {0};
1043 std::ofstream x(
"x"), y(
"y");
1048 std::cout <<
"Surface mesh has " <<
triangulation.n_active_cells()
1055Note that the only essential addition is the line marked with
1056asterisks. It is worth pointing out one other thing here, though: because we
1057detach the manifold description from the surface mesh, whenever we use a
1058mapping
object in the rest of the program, it has no curves boundary
1059description to go on any more. Rather, it will have to use the implicit,
1060FlatManifold class that is used on all parts of the domain not
1061explicitly assigned a different manifold object. Consequently, whether we use
1063using a bilinear approximation.
1065All these drawbacks aside, the resulting pictures are still pretty. The only
1066other differences to what's in @ref step_38
"step-38" is that we changed the right hand side
1067to @f$f(\mathbf x)=\
sin x_3@f$ and the boundary
values (through the
1068<code>Solution</code>
class) to @f$u(\mathbf x)|_{\partial\Omega}=\
cos x_3@f$. Of
1069course, we now no longer know the exact solution, so the computation of the
1070error at the
end of <code>LaplaceBeltrami::run</code> will yield a meaningless
1074<a name=
"step_38-PlainProg"></a>
1075<h1> The plain program</h1>
1076@include
"step-38.cc"
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
numbers::NumberTraits< Number >::real_type norm() const
void refine_global(const unsigned int times=1)
return_type extract_boundary_mesh(const MeshType< dim, spacedim > &volume_mesh, MeshType< dim - 1, spacedim > &surface_mesh, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())
void half_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
static constexpr double PI
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)