deal.II version GIT relicensing-1721-g8100761196 2024-08-31 12:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-38.h
Go to the documentation of this file.
1) const
487 *   {
488 *   return (-8. * p[0] * p[1]);
489 *   }
490 *  
491 *  
492 *   template <>
493 *   double RightHandSide<3>::value(const Point<3> &p,
494 *   const unsigned int /*component*/) const
495 *   {
496 *   using numbers::PI;
497 *  
499 *  
500 *   hessian[0][0] = -PI * PI * sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
501 *   hessian[1][1] = -PI * PI * sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
502 *   hessian[2][2] = sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
503 *  
504 *   hessian[0][1] = -PI * PI * cos(PI * p[0]) * sin(PI * p[1]) * exp(p[2]);
505 *   hessian[1][0] = -PI * PI * cos(PI * p[0]) * sin(PI * p[1]) * exp(p[2]);
506 *  
507 *   hessian[0][2] = PI * cos(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
508 *   hessian[2][0] = PI * cos(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
509 *  
510 *   hessian[1][2] = -PI * sin(PI * p[0]) * sin(PI * p[1]) * exp(p[2]);
511 *   hessian[2][1] = -PI * sin(PI * p[0]) * sin(PI * p[1]) * exp(p[2]);
512 *  
514 *   gradient[0] = PI * cos(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
515 *   gradient[1] = -PI * sin(PI * p[0]) * sin(PI * p[1]) * exp(p[2]);
516 *   gradient[2] = sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
517 *  
518 *   Point<3> normal = p;
519 *   normal /= p.norm();
520 *  
521 *   return (-trace(hessian) + 2 * (gradient * normal) +
522 *   (hessian * normal) * normal);
523 *   }
524 *  
525 *  
526 * @endcode
527 *
528 *
529 * <a name="step_38-ImplementationofthecodeLaplaceBeltramiProblemcodeclass"></a>
530 * <h3>Implementation of the <code>LaplaceBeltramiProblem</code> class</h3>
531 *
532
533 *
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
537 * DoF handler to the triangulation:
538 *
539 * @code
540 *   template <int spacedim>
541 *   LaplaceBeltramiProblem<spacedim>::LaplaceBeltramiProblem(
542 *   const unsigned degree)
543 *   : fe(degree)
544 *   , dof_handler(triangulation)
545 *   , mapping(degree)
546 *   {}
547 *  
548 *  
549 * @endcode
550 *
551 *
552 * <a name="step_38-LaplaceBeltramiProblemmake_grid_and_dofs"></a>
553 * <h4>LaplaceBeltramiProblem::make_grid_and_dofs</h4>
554 *
555
556 *
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
563 * using the facilities in the GridGenerator namespace.
564 *
565
566 *
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.
581 *
582
583 *
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.
595 *
596
597 *
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
600 * programs.
601 *
602 * @code
603 *   template <int spacedim>
604 *   void LaplaceBeltramiProblem<spacedim>::make_grid_and_dofs()
605 *   {
606 *   {
607 *   Triangulation<spacedim> volume_mesh;
608 *   GridGenerator::half_hyper_ball(volume_mesh);
609 *  
610 *   const std::set<types::boundary_id> boundary_ids = {0};
611 *  
612 *   GridGenerator::extract_boundary_mesh(volume_mesh,
613 *   triangulation,
614 *   boundary_ids);
615 *   }
616 *   triangulation.set_all_manifold_ids(0);
617 *   triangulation.set_manifold(0, SphericalManifold<dim, spacedim>());
618 *  
619 *   triangulation.refine_global(4);
620 *  
621 *   std::cout << "Surface mesh has " << triangulation.n_active_cells()
622 *   << " cells." << std::endl;
623 *  
624 *   dof_handler.distribute_dofs(fe);
625 *  
626 *   std::cout << "Surface mesh has " << dof_handler.n_dofs()
627 *   << " degrees of freedom." << std::endl;
628 *  
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);
632 *  
633 *   system_matrix.reinit(sparsity_pattern);
634 *  
635 *   solution.reinit(dof_handler.n_dofs());
636 *   system_rhs.reinit(dof_handler.n_dofs());
637 *   }
638 *  
639 *  
640 * @endcode
641 *
642 *
643 * <a name="step_38-LaplaceBeltramiProblemassemble_system"></a>
644 * <h4>LaplaceBeltramiProblem::assemble_system</h4>
645 *
646
647 *
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
655 * either:
656 *
657 * @code
658 *   template <int spacedim>
659 *   void LaplaceBeltramiProblem<spacedim>::assemble_system()
660 *   {
661 *   system_matrix = 0;
662 *   system_rhs = 0;
663 *  
664 *   const QGauss<dim> quadrature_formula(2 * fe.degree);
665 *   FEValues<dim, spacedim> fe_values(mapping,
666 *   fe,
667 *   quadrature_formula,
668 *   update_values | update_gradients |
669 *   update_quadrature_points |
670 *   update_JxW_values);
671 *  
672 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
673 *   const unsigned int n_q_points = quadrature_formula.size();
674 *  
675 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
676 *   Vector<double> cell_rhs(dofs_per_cell);
677 *  
678 *   std::vector<double> rhs_values(n_q_points);
679 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
680 *  
681 *   RightHandSide<spacedim> rhs;
682 *  
683 *   for (const auto &cell : dof_handler.active_cell_iterators())
684 *   {
685 *   cell_matrix = 0;
686 *   cell_rhs = 0;
687 *  
688 *   fe_values.reinit(cell);
689 *  
690 *   rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
691 *  
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);
698 *  
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);
703 *  
704 *   cell->get_dof_indices(local_dof_indices);
705 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
706 *   {
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));
711 *  
712 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
713 *   }
714 *   }
715 *  
716 *   std::map<types::global_dof_index, double> boundary_values;
717 *   VectorTools::interpolate_boundary_values(
718 *   mapping, dof_handler, 0, Solution<spacedim>(), boundary_values);
719 *  
720 *   MatrixTools::apply_boundary_values(
721 *   boundary_values, system_matrix, solution, system_rhs, false);
722 *   }
723 *  
724 *  
725 *  
726 * @endcode
727 *
728 *
729 * <a name="step_38-LaplaceBeltramiProblemsolve"></a>
730 * <h4>LaplaceBeltramiProblem::solve</h4>
731 *
732
733 *
734 * The next function is the one that solves the linear system. Here, too, no
735 * changes are necessary:
736 *
737 * @code
738 *   template <int spacedim>
739 *   void LaplaceBeltramiProblem<spacedim>::solve()
740 *   {
741 *   SolverControl solver_control(solution.size(), 1e-7 * system_rhs.l2_norm());
742 *   SolverCG<Vector<double>> cg(solver_control);
743 *  
744 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
745 *   preconditioner.initialize(system_matrix, 1.2);
746 *  
747 *   cg.solve(system_matrix, solution, system_rhs, preconditioner);
748 *   }
749 *  
750 *  
751 *  
752 * @endcode
753 *
754 *
755 * <a name="step_38-LaplaceBeltramiProblemoutput_result"></a>
756 * <h4>LaplaceBeltramiProblem::output_result</h4>
757 *
758
759 *
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
762 * pointing out:
763 *
764
765 *
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.
786 *
787 * @code
788 *   template <int spacedim>
789 *   void LaplaceBeltramiProblem<spacedim>::output_results() const
790 *   {
791 *   DataOut<dim, spacedim> data_out;
792 *   data_out.attach_dof_handler(dof_handler);
793 *   data_out.add_data_vector(solution,
794 *   "solution",
795 *   DataOut<dim, spacedim>::type_dof_data);
796 *   data_out.build_patches(mapping, mapping.get_degree());
797 *  
798 *   const std::string filename =
799 *   "solution-" + std::to_string(spacedim) + "d.vtk";
800 *   std::ofstream output(filename);
801 *   data_out.write_vtk(output);
802 *   }
803 *  
804 *  
805 *  
806 * @endcode
807 *
808 *
809 * <a name="step_38-LaplaceBeltramiProblemcompute_error"></a>
810 * <h4>LaplaceBeltramiProblem::compute_error</h4>
811 *
812
813 *
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.
820 *
821 * @code
822 *   template <int spacedim>
823 *   void LaplaceBeltramiProblem<spacedim>::compute_error() const
824 *   {
825 *   Vector<float> difference_per_cell(triangulation.n_active_cells());
826 *   VectorTools::integrate_difference(mapping,
827 *   dof_handler,
828 *   solution,
829 *   Solution<spacedim>(),
830 *   difference_per_cell,
831 *   QGauss<dim>(2 * fe.degree + 1),
832 *   VectorTools::H1_norm);
833 *  
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;
838 *   }
839 *  
840 *  
841 *  
842 * @endcode
843 *
844 *
845 * <a name="step_38-LaplaceBeltramiProblemrun"></a>
846 * <h4>LaplaceBeltramiProblem::run</h4>
847 *
848
849 *
850 * The last function provides the top-level logic. Its contents are
851 * self-explanatory:
852 *
853 * @code
854 *   template <int spacedim>
855 *   void LaplaceBeltramiProblem<spacedim>::run()
856 *   {
857 *   make_grid_and_dofs();
858 *   assemble_system();
859 *   solve();
860 *   output_results();
861 *   compute_error();
862 *   }
863 *   } // namespace Step38
864 *  
865 *  
866 * @endcode
867 *
868 *
869 * <a name="step_38-Themainfunction"></a>
870 * <h3>The main() function</h3>
871 *
872
873 *
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:
877 *
878 * @code
879 *   int main()
880 *   {
881 *   try
882 *   {
883 *   using namespace Step38;
884 *  
885 *   LaplaceBeltramiProblem<3> laplace_beltrami;
886 *   laplace_beltrami.run();
887 *   }
888 *   catch (std::exception &exc)
889 *   {
890 *   std::cerr << std::endl
891 *   << std::endl
892 *   << "----------------------------------------------------"
893 *   << std::endl;
894 *   std::cerr << "Exception on processing: " << std::endl
895 *   << exc.what() << std::endl
896 *   << "Aborting!" << std::endl
897 *   << "----------------------------------------------------"
898 *   << std::endl;
899 *   return 1;
900 *   }
901 *   catch (...)
902 *   {
903 *   std::cerr << std::endl
904 *   << std::endl
905 *   << "----------------------------------------------------"
906 *   << std::endl;
907 *   std::cerr << "Unknown exception!" << std::endl
908 *   << "Aborting!" << std::endl
909 *   << "----------------------------------------------------"
910 *   << std::endl;
911 *   return 1;
912 *   }
913 *  
914 *   return 0;
915 *   }
916 * @endcode
917<a name="step_38-Results"></a><h1>Results</h1>
918
919
920When you run the program, the following output should be printed on screen:
921
922@verbatim
923Surface mesh has 1280 cells.
924Surface mesh has 5185 degrees of freedom.
925H1 error = 0.0217136
926@endverbatim
927
928
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
933output:
934
935@verbatim
936Surface mesh has 5120 cells.
937Surface mesh has 20609 degrees of freedom.
938H1 error = 0.00543481
939@endverbatim
940
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:
945@verbatim
9460.339438
9470.0864385
9480.0217136
9490.00543481
9500.00135913
951@endverbatim
952
953Finally, the program produces graphical output that we can visualize. Here is
954a plot of the results:
955
956<img src="https://www.dealii.org/images/steps/developer/step-38.solution-3d.png" alt="">
957
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
960so:
961@code
962 LaplaceBeltramiProblem<2> laplace_beltrami;
963@endcode
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):
969
970<img src="https://www.dealii.org/images/steps/developer/step-38.solution-2d.png" alt="">
971
972
973<a name="step-38-extensions"></a>
974<a name="step_38-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
975
976
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.
982
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:
987
988<img src="https://www.dealii.org/images/steps/developer/step-38.warp-1.png" alt="">
989
990<img src="https://www.dealii.org/images/steps/developer/step-38.warp-2.png" alt="">
991
992The way to produce such a mesh is by using the GridTools::transform()
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):
996
997@code
998template <int spacedim>
999Point<spacedim> warp(const Point<spacedim> &p)
1000{
1001 Point<spacedim> q = p;
1002 q[spacedim-1] *= 10;
1003
1004 if (spacedim >= 2)
1005 q[0] += 2*std::sin(q[spacedim-1]);
1006 if (spacedim >= 3)
1007 q[1] += 2*std::cos(q[spacedim-1]);
1008
1009 return q;
1010}
1011@endcode
1012
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:
1027
1028@code
1029template <int spacedim>
1030void LaplaceBeltrami<spacedim>::make_grid_and_dofs()
1031{
1032 {
1033 Triangulation<spacedim> volume_mesh;
1034 GridGenerator::half_hyper_ball(volume_mesh);
1035
1036 volume_mesh.refine_global(4);
1037
1038 const std::set<types::boundary_id> boundary_ids = {0};
1039
1041 boundary_ids);
1042 GridTools::transform(&warp<spacedim>, triangulation); /* ** */
1043 std::ofstream x("x"), y("y");
1044 GridOut().write_gnuplot(volume_mesh, x);
1046 }
1047
1048 std::cout << "Surface mesh has " << triangulation.n_active_cells()
1049 << " cells."
1050 << std::endl;
1051 ...
1052}
1053@endcode
1054
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
1062MappingQ(2), MappingQ(15) or MappingQ1, each cell of our mesh will be mapped
1063using a bilinear approximation.
1064
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
1071number.
1072 *
1073 *
1074<a name="step_38-PlainProg"></a>
1075<h1> The plain program</h1>
1076@include "step-38.cc"
1077*/
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition grid_out.cc:4608
Definition point.h:111
numbers::NumberTraits< Number >::real_type norm() const
void refine_global(const unsigned int times=1)
Point< 3 > vertices[4]
Point< 2 > first
Definition grid_out.cc:4623
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 > &center=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)
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
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
Definition numbers.h:254
::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 > &)