495 *
return (-8. * p(0) * p(1));
500 *
double RightHandSide<3>::value(
const Point<3> &p,
501 *
const unsigned int )
const
526 * normal /= p.
norm();
528 *
return (-
trace(hessian) + 2 * (gradient * normal) +
529 * (hessian * normal) * normal);
536 * <a name=
"ImplementationofthecodeLaplaceBeltramiProblemcodeclass"></a>
537 * <h3>Implementation of the <code>LaplaceBeltramiProblem</code>
class</h3>
541 * The rest of the program is actually quite unspectacular
if you know
542 * @ref step_4
"step-4". Our
first step is to define the constructor, setting the
543 * polynomial degree of the finite element and mapping, and associating the
547 *
template <
int spacedim>
548 * LaplaceBeltramiProblem<spacedim>::LaplaceBeltramiProblem(
549 *
const unsigned degree)
559 * <a name=
"LaplaceBeltramiProblemmake_grid_and_dofs"></a>
560 * <h4>LaplaceBeltramiProblem::make_grid_and_dofs</h4>
564 * The next step is to create the mesh, distribute degrees of freedom, and
565 *
set up the various variables that describe the linear system. All of
566 * these steps are standard with the exception of how to create a mesh that
567 * describes a surface. We could generate a mesh
for the domain we are
568 * interested in, generate a
triangulation using a mesh generator, and read
569 * it in
using the
GridIn class. Or, as we
do here, we generate the mesh
574 * In particular, what we
're going to do is this (enclosed between the set
575 * of braces below): we generate a <code>spacedim</code> dimensional mesh
576 * for the half disk (in 2d) or half ball (in 3d), using the
577 * GridGenerator::half_hyper_ball function. This function sets the boundary
578 * indicators of all faces on the outside of the boundary to zero for the
579 * ones located on the perimeter of the disk/ball, and one on the straight
580 * part that splits the full disk/ball into two halves. The next step is the
581 * main point: The GridGenerator::extract_boundary_mesh function creates a
582 * mesh that consists of those cells that are the faces of the previous mesh,
583 * i.e. it describes the <i>surface</i> cells of the original (volume)
584 * mesh. However, we do not want all faces: only those on the perimeter of
585 * the disk or ball which carry boundary indicator zero; we can select these
586 * cells using a set of boundary indicators that we pass to
587 * GridGenerator::extract_boundary_mesh.
591 * There is one point that needs to be mentioned. In order to refine a
592 * surface mesh appropriately if the manifold is curved (similarly to
593 * refining the faces of cells that are adjacent to a curved boundary), the
594 * triangulation has to have an object attached to it that describes where
595 * new vertices should be located. If you don't attach such a boundary
596 * object, they will be located halfway between existing
vertices;
this is
597 * appropriate
if you have a domain with straight boundaries (
e.g. a
598 * polygon) but not when, as here, the manifold has curvature. So
for things
599 * to work properly, we need to attach a manifold
object to our (surface)
600 *
triangulation, in much the same way as we
've already done in 1d for the
601 * boundary. We create such an object and attach it to the triangulation.
605 * The final step in creating the mesh is to refine it a number of
606 * times. The rest of the function is the same as in previous tutorial
610 * template <int spacedim>
611 * void LaplaceBeltramiProblem<spacedim>::make_grid_and_dofs()
614 * Triangulation<spacedim> volume_mesh;
615 * GridGenerator::half_hyper_ball(volume_mesh);
617 * std::set<types::boundary_id> boundary_ids;
618 * boundary_ids.insert(0);
620 * GridGenerator::extract_boundary_mesh(volume_mesh,
624 * triangulation.set_all_manifold_ids(0);
625 * triangulation.set_manifold(0, SphericalManifold<dim, spacedim>());
627 * triangulation.refine_global(4);
629 * std::cout << "Surface mesh has " << triangulation.n_active_cells()
630 * << " cells." << std::endl;
632 * dof_handler.distribute_dofs(fe);
634 * std::cout << "Surface mesh has " << dof_handler.n_dofs()
635 * << " degrees of freedom." << std::endl;
637 * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
638 * DoFTools::make_sparsity_pattern(dof_handler, dsp);
639 * sparsity_pattern.copy_from(dsp);
641 * system_matrix.reinit(sparsity_pattern);
643 * solution.reinit(dof_handler.n_dofs());
644 * system_rhs.reinit(dof_handler.n_dofs());
651 * <a name="LaplaceBeltramiProblemassemble_system"></a>
652 * <h4>LaplaceBeltramiProblem::assemble_system</h4>
656 * The following is the central function of this program, assembling the
657 * matrix that corresponds to the surface Laplacian (Laplace-Beltrami
658 * operator). Maybe surprisingly, it actually looks exactly the same as for
659 * the regular Laplace operator discussed in, for example, @ref step_4 "step-4". The key
660 * is that the FEValues::shape_grad() function does the magic: It returns
661 * the surface gradient @f$\nabla_K \phi_i(x_q)@f$ of the @f$i@f$th shape function
662 * at the @f$q@f$th quadrature point. The rest then does not need any changes
666 * template <int spacedim>
667 * void LaplaceBeltramiProblem<spacedim>::assemble_system()
672 * const QGauss<dim> quadrature_formula(2 * fe.degree);
673 * FEValues<dim, spacedim> fe_values(mapping,
675 * quadrature_formula,
676 * update_values | update_gradients |
677 * update_quadrature_points |
678 * update_JxW_values);
680 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
681 * const unsigned int n_q_points = quadrature_formula.size();
683 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
684 * Vector<double> cell_rhs(dofs_per_cell);
686 * std::vector<double> rhs_values(n_q_points);
687 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
689 * RightHandSide<spacedim> rhs;
691 * for (const auto &cell : dof_handler.active_cell_iterators())
696 * fe_values.reinit(cell);
698 * rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
700 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
701 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
702 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
703 * cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
704 * fe_values.shape_grad(j, q_point) *
705 * fe_values.JxW(q_point);
707 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
708 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
709 * cell_rhs(i) += fe_values.shape_value(i, q_point) *
710 * rhs_values[q_point] * fe_values.JxW(q_point);
712 * cell->get_dof_indices(local_dof_indices);
713 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
715 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
716 * system_matrix.add(local_dof_indices[i],
717 * local_dof_indices[j],
718 * cell_matrix(i, j));
720 * system_rhs(local_dof_indices[i]) += cell_rhs(i);
724 * std::map<types::global_dof_index, double> boundary_values;
725 * VectorTools::interpolate_boundary_values(
726 * mapping, dof_handler, 0, Solution<spacedim>(), boundary_values);
728 * MatrixTools::apply_boundary_values(
729 * boundary_values, system_matrix, solution, system_rhs, false);
737 * <a name="LaplaceBeltramiProblemsolve"></a>
738 * <h4>LaplaceBeltramiProblem::solve</h4>
742 * The next function is the one that solves the linear system. Here, too, no
743 * changes are necessary:
746 * template <int spacedim>
747 * void LaplaceBeltramiProblem<spacedim>::solve()
749 * SolverControl solver_control(solution.size(), 1e-7 * system_rhs.l2_norm());
750 * SolverCG<Vector<double>> cg(solver_control);
752 * PreconditionSSOR<SparseMatrix<double>> preconditioner;
753 * preconditioner.initialize(system_matrix, 1.2);
755 * cg.solve(system_matrix, solution, system_rhs, preconditioner);
763 * <a name="LaplaceBeltramiProblemoutput_result"></a>
764 * <h4>LaplaceBeltramiProblem::output_result</h4>
768 * This is the function that generates graphical output from the
769 * solution. Most of it is boilerplate code, but there are two points worth
774 * - The DataOut::add_data_vector() function can take two kinds of vectors:
775 * Either vectors that have one value per degree of freedom defined by the
776 * DoFHandler object previously attached via DataOut::attach_dof_handler();
777 * and vectors that have one value for each cell of the triangulation, for
778 * example to output estimated errors for each cell. Typically, the
779 * DataOut class knows to tell these two kinds of vectors apart: there are
780 * almost always more degrees of freedom than cells, so we can
781 * differentiate by the two kinds looking at the length of a vector. We
782 * could do the same here, but only because we got lucky: we use a half
783 * sphere. If we had used the whole sphere as domain and @f$Q_1@f$ elements,
784 * we would have the same number of cells as vertices and consequently the
785 * two kinds of vectors would have the same number of elements. To avoid
786 * the resulting confusion, we have to tell the DataOut::add_data_vector()
787 * function which kind of vector we have: DoF data. This is what the third
788 * argument to the function does.
789 * - The DataOut::build_patches() function can generate output that subdivides
790 * each cell so that visualization programs can resolve curved manifolds
791 * or higher polynomial degree shape functions better. We here subdivide
792 * each element in each coordinate direction as many times as the
793 * polynomial degree of the finite element in use.
796 * template <int spacedim>
797 * void LaplaceBeltramiProblem<spacedim>::output_results() const
799 * DataOut<dim, DoFHandler<dim, spacedim>> data_out;
800 * data_out.attach_dof_handler(dof_handler);
801 * data_out.add_data_vector(
804 * DataOut<dim, DoFHandler<dim, spacedim>>::type_dof_data);
805 * data_out.build_patches(mapping, mapping.get_degree());
807 * const std::string filename =
808 * "solution-" + std::to_string(spacedim) + "d.vtk";
809 * std::ofstream output(filename);
810 * data_out.write_vtk(output);
818 * <a name="LaplaceBeltramiProblemcompute_error"></a>
819 * <h4>LaplaceBeltramiProblem::compute_error</h4>
823 * This is the last piece of functionality: we want to compute the error in
824 * the numerical solution. It is a verbatim copy of the code previously
825 * shown and discussed in @ref step_7 "step-7". As mentioned in the introduction, the
826 * <code>Solution</code> class provides the (tangential) gradient of the
827 * solution. To avoid evaluating the error only a superconvergence points,
828 * we choose a quadrature rule of sufficiently high order.
831 * template <int spacedim>
832 * void LaplaceBeltramiProblem<spacedim>::compute_error() const
834 * Vector<float> difference_per_cell(triangulation.n_active_cells());
835 * VectorTools::integrate_difference(mapping,
838 * Solution<spacedim>(),
839 * difference_per_cell,
840 * QGauss<dim>(2 * fe.degree + 1),
841 * VectorTools::H1_norm);
843 * double h1_error = VectorTools::compute_global_error(triangulation,
844 * difference_per_cell,
845 * VectorTools::H1_norm);
846 * std::cout << "H1 error = " << h1_error << std::endl;
854 * <a name="LaplaceBeltramiProblemrun"></a>
855 * <h4>LaplaceBeltramiProblem::run</h4>
859 * The last function provides the top-level logic. Its contents are
863 * template <int spacedim>
864 * void LaplaceBeltramiProblem<spacedim>::run()
866 * make_grid_and_dofs();
872 * } // namespace Step38
878 * <a name="Themainfunction"></a>
879 * <h3>The main() function</h3>
883 * The remainder of the program is taken up by the <code>main()</code>
884 * function. It follows exactly the general layout first introduced in @ref step_6 "step-6"
885 * and used in all following tutorial programs:
892 * using namespace Step38;
894 * LaplaceBeltramiProblem<3> laplace_beltrami;
895 * laplace_beltrami.run();
897 * catch (std::exception &exc)
899 * std::cerr << std::endl
901 * << "----------------------------------------------------"
903 * std::cerr << "Exception on processing: " << std::endl
904 * << exc.what() << std::endl
905 * << "Aborting!" << std::endl
906 * << "----------------------------------------------------"
912 * std::cerr << std::endl
914 * << "----------------------------------------------------"
916 * std::cerr << "Unknown exception!" << std::endl
917 * << "Aborting!" << std::endl
918 * << "----------------------------------------------------"
926<a name="Results"></a><h1>Results</h1>
929When you run the program, the following output should be printed on screen:
932Surface mesh has 1280 cells.
933Surface mesh has 5185 degrees of freedom.
938By playing around with the number of global refinements in the
939<code>LaplaceBeltrami::make_grid_and_dofs</code> function you increase or decrease mesh
940refinement. For example, doing one more refinement and only running the 3d surface
941problem yields the following
945Surface mesh has 5120 cells.
946Surface mesh has 20609 degrees of freedom.
950This is what we expect: make the mesh size smaller by a factor of two and the
951error goes down by a factor of four (remember that we use bi-quadratic
952elements). The full sequence of errors from one to five refinements looks like
953this, neatly following the theoretically predicted pattern:
962Finally, the program produces graphical output that we can visualize. Here is
963a plot of the results:
965<img src="https://www.dealii.org/images/steps/developer/step-38.solution-3d.png" alt="">
967The program also works for 1d curves in 2d, not just 2d surfaces in 3d. You
968can test this by changing the template argument in <code>main()</code> like
971 LaplaceBeltramiProblem<2> laplace_beltrami;
973The domain is a curve in 2d, and we can visualize the solution by using the
974third dimension (and color) to denote the value of the function @f$u(x)@f$. This
975then looks like so (the white curve is the domain, the colored curve is the
976solution extruded into the third dimension, clearly showing the change in sign
977as the curve moves from one quadrant of the domain into the adjacent one):
979<img src="https://www.dealii.org/images/steps/developer/step-38.solution-2d.png" alt="">
982<a name="extensions"></a>
983<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
986Computing on surfaces only becomes interesting if the surface is more
987interesting than just a half sphere. To achieve this, deal.II can read
988meshes that describe surfaces through the usual GridIn class. Or, in case you
989have an analytic description, a simple mesh can sometimes be stretched and
990bent into a shape we are interested in.
992Let us consider a relatively simple example: we take the half sphere we used
993before, we stretch it by a factor of 10 in the z-direction, and then we jumble
994the x- and y-coordinates a bit. Let's show the computational domain and the
995solution
first before we go into details of the implementation below:
997<img src=
"https://www.dealii.org/images/steps/developer/step-38.warp-1.png" alt=
"">
999<img src=
"https://www.dealii.org/images/steps/developer/step-38.warp-2.png" alt=
"">
1002function. It needs a way to
transform each individual mesh
point to a
1003different position. Let us here use the following, rather simple function
1004(remember: stretch in
one direction, jumble in the other two):
1007template <
int spacedim>
1008Point<spacedim> warp(const
Point<spacedim> &p)
1011 q[spacedim-1] *= 10;
1022If we followed the <code>LaplaceBeltrami::make_grid_and_dofs</code> function, we would
1023extract the half spherical surface mesh as before, warp it into the shape we
1024want, and
refine as often as necessary. This is not quite as simple as we
'd
1025like here, though: refining requires that we have an appropriate manifold
1026object attached to the triangulation that describes where new vertices of the
1027mesh should be located upon refinement. I'm sure it
's possible to describe
1028this manifold in a not-too-complicated way by simply undoing the
1029transformation above (yielding the spherical surface again), finding the
1030location of a new point on the sphere, and then re-warping the result. But I'm
1031a lazy person, and since doing
this is not really the
point here, let
's just
1032make our lives a bit easier: we'll
extract the half sphere,
refine it as
1033often as necessary, get rid of the
object that describes the manifold since we
1034now no longer need it, and then
finally warp the mesh. With the function
1035above,
this would look as follows:
1038template <
int spacedim>
1039void LaplaceBeltrami<spacedim>::make_grid_and_dofs()
1047 std::set<types::boundary_id> boundary_ids;
1048 boundary_ids.insert(0);
1053 std::ofstream x(
"x"), y(
"y");
1058 std::cout <<
"Surface mesh has " <<
triangulation.n_active_cells()
1065Note that the only essential addition is the line marked with
1066asterisks. It is worth pointing out
one other thing here, though: because we
1067detach the manifold description from the surface mesh, whenever we use a
1068mapping
object in the rest of the program, it has no curves boundary
1069description to go on any more. Rather, it will have to use the implicit,
1070FlatManifold class that is used on all parts of the domain not
1071explicitly assigned a different manifold object. Consequently, whether we use
1073using a bilinear approximation.
1075All these drawbacks aside, the resulting pictures are still pretty. The only
1076other differences to what's in @ref step_38 "step-38" is that we changed the right hand side
1077to @f$f(\mathbf x)=
\sin x_3@f$ and the boundary
values (through the
1078<code>Solution</code>
class) to @f$u(\mathbf x)|_{\partial\Omega}=
\cos x_3@f$. Of
1079course, we now no longer know the exact solution, so the computation of the
1084<a name=
"PlainProg"></a>
1085<h1> The plain program</h1>
1086@include
"step-38.cc"
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
constexpr Number trace(const SymmetricTensor< 2, dim, Number > &d)
numbers::NumberTraits< Number >::real_type norm() const
void refine_global(const unsigned int times=1)
VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &x)
VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &x)
VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &x)
__global__ void set(Number *val, const Number s, const size_type N)
std::map< typename MeshType< dim - 1, spacedim >::cell_iterator, typename MeshType< dim, spacedim >::face_iterator > 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)
static const types::blas_int one
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)
void run(const Iterator &begin, const typename identity< Iterator >::type &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)
static constexpr double PI
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation