Reference documentation for deal.II version 9.6.0
\(\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
487) const
488 *   {
489 *   return (-8. * p[0] * p[1]);
490 *   }
491 *  
492 *  
493 *   template <>
494 *   double RightHandSide<3>::value(const Point<3> &p,
495 *   const unsigned int /*component*/) const
496 *   {
497 *   using numbers::PI;
498 *  
500 *  
501 *   hessian[0][0] = -PI * PI * sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
502 *   hessian[1][1] = -PI * PI * sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
503 *   hessian[2][2] = sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
504 *  
505 *   hessian[0][1] = -PI * PI * cos(PI * p[0]) * sin(PI * p[1]) * exp(p[2]);
506 *   hessian[1][0] = -PI * PI * cos(PI * p[0]) * sin(PI * p[1]) * exp(p[2]);
507 *  
508 *   hessian[0][2] = PI * cos(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
509 *   hessian[2][0] = PI * cos(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
510 *  
511 *   hessian[1][2] = -PI * sin(PI * p[0]) * sin(PI * p[1]) * exp(p[2]);
512 *   hessian[2][1] = -PI * sin(PI * p[0]) * sin(PI * p[1]) * exp(p[2]);
513 *  
515 *   gradient[0] = PI * cos(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
516 *   gradient[1] = -PI * sin(PI * p[0]) * sin(PI * p[1]) * exp(p[2]);
517 *   gradient[2] = sin(PI * p[0]) * cos(PI * p[1]) * exp(p[2]);
518 *  
519 *   Point<3> normal = p;
520 *   normal /= p.norm();
521 *  
522 *   return (-trace(hessian) + 2 * (gradient * normal) +
523 *   (hessian * normal) * normal);
524 *   }
525 *  
526 *  
527 * @endcode
528 *
529 *
530 * <a name="step_38-ImplementationofthecodeLaplaceBeltramiProblemcodeclass"></a>
531 * <h3>Implementation of the <code>LaplaceBeltramiProblem</code> class</h3>
532 *
533
534 *
535 * The rest of the program is actually quite unspectacular if you know
536 * @ref step_4 "step-4". Our first step is to define the constructor, setting the
537 * polynomial degree of the finite element and mapping, and associating the
538 * DoF handler to the triangulation:
539 *
540 * @code
541 *   template <int spacedim>
542 *   LaplaceBeltramiProblem<spacedim>::LaplaceBeltramiProblem(
543 *   const unsigned degree)
544 *   : fe(degree)
545 *   , dof_handler(triangulation)
546 *   , mapping(degree)
547 *   {}
548 *  
549 *  
550 * @endcode
551 *
552 *
553 * <a name="step_38-LaplaceBeltramiProblemmake_grid_and_dofs"></a>
554 * <h4>LaplaceBeltramiProblem::make_grid_and_dofs</h4>
555 *
556
557 *
558 * The next step is to create the mesh, distribute degrees of freedom, and
559 * set up the various variables that describe the linear system. All of
560 * these steps are standard with the exception of how to create a mesh that
561 * describes a surface. We could generate a mesh for the domain we are
562 * interested in, generate a triangulation using a mesh generator, and read
563 * it in using the GridIn class. Or, as we do here, we generate the mesh
564 * using the facilities in the GridGenerator namespace.
565 *
566
567 *
568 * In particular, what we're going to do is this (enclosed between the set
569 * of braces below): we generate a <code>spacedim</code> dimensional mesh
570 * for the half disk (in 2d) or half ball (in 3d), using the
571 * GridGenerator::half_hyper_ball function. This function sets the boundary
572 * indicators of all faces on the outside of the boundary to zero for the
573 * ones located on the perimeter of the disk/ball, and one on the straight
574 * part that splits the full disk/ball into two halves. The next step is the
575 * main point: The GridGenerator::extract_boundary_mesh function creates a
576 * mesh that consists of those cells that are the faces of the previous mesh,
577 * i.e. it describes the <i>surface</i> cells of the original (volume)
578 * mesh. However, we do not want all faces: only those on the perimeter of
579 * the disk or ball which carry boundary indicator zero; we can select these
580 * cells using a set of boundary indicators that we pass to
581 * GridGenerator::extract_boundary_mesh.
582 *
583
584 *
585 * There is one point that needs to be mentioned. In order to refine a
586 * surface mesh appropriately if the manifold is curved (similarly to
587 * refining the faces of cells that are adjacent to a curved boundary), the
588 * triangulation has to have an object attached to it that describes where
589 * new vertices should be located. If you don't attach such a boundary
590 * object, they will be located halfway between existing vertices; this is
591 * appropriate if you have a domain with straight boundaries (e.g. a
592 * polygon) but not when, as here, the manifold has curvature. So for things
593 * to work properly, we need to attach a manifold object to our (surface)
594 * triangulation, in much the same way as we've already done in 1d for the
595 * boundary. We create such an object and attach it to the triangulation.
596 *
597
598 *
599 * The final step in creating the mesh is to refine it a number of
600 * times. The rest of the function is the same as in previous tutorial
601 * programs.
602 *
603 * @code
604 *   template <int spacedim>
605 *   void LaplaceBeltramiProblem<spacedim>::make_grid_and_dofs()
606 *   {
607 *   {
608 *   Triangulation<spacedim> volume_mesh;
609 *   GridGenerator::half_hyper_ball(volume_mesh);
610 *  
611 *   const std::set<types::boundary_id> boundary_ids = {0};
612 *  
613 *   GridGenerator::extract_boundary_mesh(volume_mesh,
614 *   triangulation,
615 *   boundary_ids);
616 *   }
617 *   triangulation.set_all_manifold_ids(0);
618 *   triangulation.set_manifold(0, SphericalManifold<dim, spacedim>());
619 *  
620 *   triangulation.refine_global(4);
621 *  
622 *   std::cout << "Surface mesh has " << triangulation.n_active_cells()
623 *   << " cells." << std::endl;
624 *  
625 *   dof_handler.distribute_dofs(fe);
626 *  
627 *   std::cout << "Surface mesh has " << dof_handler.n_dofs()
628 *   << " degrees of freedom." << std::endl;
629 *  
630 *   DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
631 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
632 *   sparsity_pattern.copy_from(dsp);
633 *  
634 *   system_matrix.reinit(sparsity_pattern);
635 *  
636 *   solution.reinit(dof_handler.n_dofs());
637 *   system_rhs.reinit(dof_handler.n_dofs());
638 *   }
639 *  
640 *  
641 * @endcode
642 *
643 *
644 * <a name="step_38-LaplaceBeltramiProblemassemble_system"></a>
645 * <h4>LaplaceBeltramiProblem::assemble_system</h4>
646 *
647
648 *
649 * The following is the central function of this program, assembling the
650 * matrix that corresponds to the surface Laplacian (Laplace-Beltrami
651 * operator). Maybe surprisingly, it actually looks exactly the same as for
652 * the regular Laplace operator discussed in, for example, @ref step_4 "step-4". The key
653 * is that the FEValues::shape_grad() function does the magic: It returns
654 * the surface gradient @f$\nabla_K \phi_i(x_q)@f$ of the @f$i@f$th shape function
655 * at the @f$q@f$th quadrature point. The rest then does not need any changes
656 * either:
657 *
658 * @code
659 *   template <int spacedim>
660 *   void LaplaceBeltramiProblem<spacedim>::assemble_system()
661 *   {
662 *   system_matrix = 0;
663 *   system_rhs = 0;
664 *  
665 *   const QGauss<dim> quadrature_formula(2 * fe.degree);
666 *   FEValues<dim, spacedim> fe_values(mapping,
667 *   fe,
668 *   quadrature_formula,
669 *   update_values | update_gradients |
670 *   update_quadrature_points |
671 *   update_JxW_values);
672 *  
673 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
674 *   const unsigned int n_q_points = quadrature_formula.size();
675 *  
676 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
677 *   Vector<double> cell_rhs(dofs_per_cell);
678 *  
679 *   std::vector<double> rhs_values(n_q_points);
680 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
681 *  
682 *   RightHandSide<spacedim> rhs;
683 *  
684 *   for (const auto &cell : dof_handler.active_cell_iterators())
685 *   {
686 *   cell_matrix = 0;
687 *   cell_rhs = 0;
688 *  
689 *   fe_values.reinit(cell);
690 *  
691 *   rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
692 *  
693 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
694 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
695 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
696 *   cell_matrix(i, j) += fe_values.shape_grad(i, q_point) *
697 *   fe_values.shape_grad(j, q_point) *
698 *   fe_values.JxW(q_point);
699 *  
700 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
701 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
702 *   cell_rhs(i) += fe_values.shape_value(i, q_point) *
703 *   rhs_values[q_point] * fe_values.JxW(q_point);
704 *  
705 *   cell->get_dof_indices(local_dof_indices);
706 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
707 *   {
708 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
709 *   system_matrix.add(local_dof_indices[i],
710 *   local_dof_indices[j],
711 *   cell_matrix(i, j));
712 *  
713 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
714 *   }
715 *   }
716 *  
717 *   std::map<types::global_dof_index, double> boundary_values;
718 *   VectorTools::interpolate_boundary_values(
719 *   mapping, dof_handler, 0, Solution<spacedim>(), boundary_values);
720 *  
721 *   MatrixTools::apply_boundary_values(
722 *   boundary_values, system_matrix, solution, system_rhs, false);
723 *   }
724 *  
725 *  
726 *  
727 * @endcode
728 *
729 *
730 * <a name="step_38-LaplaceBeltramiProblemsolve"></a>
731 * <h4>LaplaceBeltramiProblem::solve</h4>
732 *
733
734 *
735 * The next function is the one that solves the linear system. Here, too, no
736 * changes are necessary:
737 *
738 * @code
739 *   template <int spacedim>
740 *   void LaplaceBeltramiProblem<spacedim>::solve()
741 *   {
742 *   SolverControl solver_control(solution.size(), 1e-7 * system_rhs.l2_norm());
743 *   SolverCG<Vector<double>> cg(solver_control);
744 *  
745 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
746 *   preconditioner.initialize(system_matrix, 1.2);
747 *  
748 *   cg.solve(system_matrix, solution, system_rhs, preconditioner);
749 *   }
750 *  
751 *  
752 *  
753 * @endcode
754 *
755 *
756 * <a name="step_38-LaplaceBeltramiProblemoutput_result"></a>
757 * <h4>LaplaceBeltramiProblem::output_result</h4>
758 *
759
760 *
761 * This is the function that generates graphical output from the
762 * solution. Most of it is boilerplate code, but there are two points worth
763 * pointing out:
764 *
765
766 *
767 * - The DataOut::add_data_vector() function can take two kinds of vectors:
768 * Either vectors that have one value per degree of freedom defined by the
769 * DoFHandler object previously attached via DataOut::attach_dof_handler();
770 * and vectors that have one value for each cell of the triangulation, for
771 * example to output estimated errors for each cell. Typically, the
772 * DataOut class knows to tell these two kinds of vectors apart: there are
773 * almost always more degrees of freedom than cells, so we can
774 * differentiate by the two kinds looking at the length of a vector. We
775 * could do the same here, but only because we got lucky: we use a half
776 * sphere. If we had used the whole sphere as domain and @f$Q_1@f$ elements,
777 * we would have the same number of cells as vertices and consequently the
778 * two kinds of vectors would have the same number of elements. To avoid
779 * the resulting confusion, we have to tell the DataOut::add_data_vector()
780 * function which kind of vector we have: DoF data. This is what the third
781 * argument to the function does.
782 * - The DataOut::build_patches() function can generate output that subdivides
783 * each cell so that visualization programs can resolve curved manifolds
784 * or higher polynomial degree shape functions better. We here subdivide
785 * each element in each coordinate direction as many times as the
786 * polynomial degree of the finite element in use.
787 *
788 * @code
789 *   template <int spacedim>
790 *   void LaplaceBeltramiProblem<spacedim>::output_results() const
791 *   {
792 *   DataOut<dim, spacedim> data_out;
793 *   data_out.attach_dof_handler(dof_handler);
794 *   data_out.add_data_vector(solution,
795 *   "solution",
796 *   DataOut<dim, spacedim>::type_dof_data);
797 *   data_out.build_patches(mapping, mapping.get_degree());
798 *  
799 *   const std::string filename =
800 *   "solution-" + std::to_string(spacedim) + "d.vtk";
801 *   std::ofstream output(filename);
802 *   data_out.write_vtk(output);
803 *   }
804 *  
805 *  
806 *  
807 * @endcode
808 *
809 *
810 * <a name="step_38-LaplaceBeltramiProblemcompute_error"></a>
811 * <h4>LaplaceBeltramiProblem::compute_error</h4>
812 *
813
814 *
815 * This is the last piece of functionality: we want to compute the error in
816 * the numerical solution. It is a verbatim copy of the code previously
817 * shown and discussed in @ref step_7 "step-7". As mentioned in the introduction, the
818 * <code>Solution</code> class provides the (tangential) gradient of the
819 * solution. To avoid evaluating the error only a superconvergence points,
820 * we choose a quadrature rule of sufficiently high order.
821 *
822 * @code
823 *   template <int spacedim>
824 *   void LaplaceBeltramiProblem<spacedim>::compute_error() const
825 *   {
826 *   Vector<float> difference_per_cell(triangulation.n_active_cells());
827 *   VectorTools::integrate_difference(mapping,
828 *   dof_handler,
829 *   solution,
830 *   Solution<spacedim>(),
831 *   difference_per_cell,
832 *   QGauss<dim>(2 * fe.degree + 1),
833 *   VectorTools::H1_norm);
834 *  
835 *   double h1_error = VectorTools::compute_global_error(triangulation,
836 *   difference_per_cell,
837 *   VectorTools::H1_norm);
838 *   std::cout << "H1 error = " << h1_error << std::endl;
839 *   }
840 *  
841 *  
842 *  
843 * @endcode
844 *
845 *
846 * <a name="step_38-LaplaceBeltramiProblemrun"></a>
847 * <h4>LaplaceBeltramiProblem::run</h4>
848 *
849
850 *
851 * The last function provides the top-level logic. Its contents are
852 * self-explanatory:
853 *
854 * @code
855 *   template <int spacedim>
856 *   void LaplaceBeltramiProblem<spacedim>::run()
857 *   {
858 *   make_grid_and_dofs();
859 *   assemble_system();
860 *   solve();
861 *   output_results();
862 *   compute_error();
863 *   }
864 *   } // namespace Step38
865 *  
866 *  
867 * @endcode
868 *
869 *
870 * <a name="step_38-Themainfunction"></a>
871 * <h3>The main() function</h3>
872 *
873
874 *
875 * The remainder of the program is taken up by the <code>main()</code>
876 * function. It follows exactly the general layout first introduced in @ref step_6 "step-6"
877 * and used in all following tutorial programs:
878 *
879 * @code
880 *   int main()
881 *   {
882 *   try
883 *   {
884 *   using namespace Step38;
885 *  
886 *   LaplaceBeltramiProblem<3> laplace_beltrami;
887 *   laplace_beltrami.run();
888 *   }
889 *   catch (std::exception &exc)
890 *   {
891 *   std::cerr << std::endl
892 *   << std::endl
893 *   << "----------------------------------------------------"
894 *   << std::endl;
895 *   std::cerr << "Exception on processing: " << std::endl
896 *   << exc.what() << std::endl
897 *   << "Aborting!" << std::endl
898 *   << "----------------------------------------------------"
899 *   << std::endl;
900 *   return 1;
901 *   }
902 *   catch (...)
903 *   {
904 *   std::cerr << std::endl
905 *   << std::endl
906 *   << "----------------------------------------------------"
907 *   << std::endl;
908 *   std::cerr << "Unknown exception!" << std::endl
909 *   << "Aborting!" << std::endl
910 *   << "----------------------------------------------------"
911 *   << std::endl;
912 *   return 1;
913 *   }
914 *  
915 *   return 0;
916 *   }
917 * @endcode
918<a name="step_38-Results"></a><h1>Results</h1>
919
920
921When you run the program, the following output should be printed on screen:
922
923@verbatim
924Surface mesh has 1280 cells.
925Surface mesh has 5185 degrees of freedom.
926H1 error = 0.0217136
927@endverbatim
928
929
930By playing around with the number of global refinements in the
931<code>LaplaceBeltrami::make_grid_and_dofs</code> function you increase or decrease mesh
932refinement. For example, doing one more refinement and only running the 3d surface
933problem yields the following
934output:
935
936@verbatim
937Surface mesh has 5120 cells.
938Surface mesh has 20609 degrees of freedom.
939H1 error = 0.00543481
940@endverbatim
941
942This is what we expect: make the mesh size smaller by a factor of two and the
943error goes down by a factor of four (remember that we use bi-quadratic
944elements). The full sequence of errors from one to five refinements looks like
945this, neatly following the theoretically predicted pattern:
946@verbatim
9470.339438
9480.0864385
9490.0217136
9500.00543481
9510.00135913
952@endverbatim
953
954Finally, the program produces graphical output that we can visualize. Here is
955a plot of the results:
956
957<img src="https://www.dealii.org/images/steps/developer/step-38.solution-3d.png" alt="">
958
959The program also works for 1d curves in 2d, not just 2d surfaces in 3d. You
960can test this by changing the template argument in <code>main()</code> like
961so:
962@code
963 LaplaceBeltramiProblem<2> laplace_beltrami;
964@endcode
965The domain is a curve in 2d, and we can visualize the solution by using the
966third dimension (and color) to denote the value of the function @f$u(x)@f$. This
967then looks like so (the white curve is the domain, the colored curve is the
968solution extruded into the third dimension, clearly showing the change in sign
969as the curve moves from one quadrant of the domain into the adjacent one):
970
971<img src="https://www.dealii.org/images/steps/developer/step-38.solution-2d.png" alt="">
972
973
974<a name="step-38-extensions"></a>
975<a name="step_38-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
976
977
978Computing on surfaces only becomes interesting if the surface is more
979interesting than just a half sphere. To achieve this, deal.II can read
980meshes that describe surfaces through the usual GridIn class. Or, in case you
981have an analytic description, a simple mesh can sometimes be stretched and
982bent into a shape we are interested in.
983
984Let us consider a relatively simple example: we take the half sphere we used
985before, we stretch it by a factor of 10 in the z-direction, and then we jumble
986the x- and y-coordinates a bit. Let's show the computational domain and the
987solution first before we go into details of the implementation below:
988
989<img src="https://www.dealii.org/images/steps/developer/step-38.warp-1.png" alt="">
990
991<img src="https://www.dealii.org/images/steps/developer/step-38.warp-2.png" alt="">
992
993The way to produce such a mesh is by using the GridTools::transform()
994function. It needs a way to transform each individual mesh point to a
995different position. Let us here use the following, rather simple function
996(remember: stretch in one direction, jumble in the other two):
997
998@code
999template <int spacedim>
1000Point<spacedim> warp(const Point<spacedim> &p)
1001{
1002 Point<spacedim> q = p;
1003 q[spacedim-1] *= 10;
1004
1005 if (spacedim >= 2)
1006 q[0] += 2*std::sin(q[spacedim-1]);
1007 if (spacedim >= 3)
1008 q[1] += 2*std::cos(q[spacedim-1]);
1009
1010 return q;
1011}
1012@endcode
1013
1014If we followed the <code>LaplaceBeltrami::make_grid_and_dofs</code> function, we would
1015extract the half spherical surface mesh as before, warp it into the shape we
1016want, and refine as often as necessary. This is not quite as simple as we'd
1017like here, though: refining requires that we have an appropriate manifold
1018object attached to the triangulation that describes where new vertices of the
1019mesh should be located upon refinement. I'm sure it's possible to describe
1020this manifold in a not-too-complicated way by simply undoing the
1021transformation above (yielding the spherical surface again), finding the
1022location of a new point on the sphere, and then re-warping the result. But I'm
1023a lazy person, and since doing this is not really the point here, let's just
1024make our lives a bit easier: we'll extract the half sphere, refine it as
1025often as necessary, get rid of the object that describes the manifold since we
1026now no longer need it, and then finally warp the mesh. With the function
1027above, this would look as follows:
1028
1029@code
1030template <int spacedim>
1031void LaplaceBeltrami<spacedim>::make_grid_and_dofs()
1032{
1033 {
1034 Triangulation<spacedim> volume_mesh;
1035 GridGenerator::half_hyper_ball(volume_mesh);
1036
1037 volume_mesh.refine_global(4);
1038
1039 const std::set<types::boundary_id> boundary_ids = {0};
1040
1042 boundary_ids);
1043 GridTools::transform(&warp<spacedim>, triangulation); /* ** */
1044 std::ofstream x("x"), y("y");
1045 GridOut().write_gnuplot(volume_mesh, x);
1047 }
1048
1049 std::cout << "Surface mesh has " << triangulation.n_active_cells()
1050 << " cells."
1051 << std::endl;
1052 ...
1053}
1054@endcode
1055
1056Note that the only essential addition is the line marked with
1057asterisks. It is worth pointing out one other thing here, though: because we
1058detach the manifold description from the surface mesh, whenever we use a
1059mapping object in the rest of the program, it has no curves boundary
1060description to go on any more. Rather, it will have to use the implicit,
1061FlatManifold class that is used on all parts of the domain not
1062explicitly assigned a different manifold object. Consequently, whether we use
1063MappingQ(2), MappingQ(15) or MappingQ1, each cell of our mesh will be mapped
1064using a bilinear approximation.
1065
1066All these drawbacks aside, the resulting pictures are still pretty. The only
1067other differences to what's in @ref step_38 "step-38" is that we changed the right hand side
1068to @f$f(\mathbf x)=\sin x_3@f$ and the boundary values (through the
1069<code>Solution</code> class) to @f$u(\mathbf x)|_{\partial\Omega}=\cos x_3@f$. Of
1070course, we now no longer know the exact solution, so the computation of the
1071error at the end of <code>LaplaceBeltrami::run</code> will yield a meaningless
1072number.
1073 *
1074 *
1075<a name="step_38-PlainProg"></a>
1076<h1> The plain program</h1>
1077@include "step-38.cc"
1078*/
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
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
Point< 3 > vertices[4]
Point< 2 > first
Definition grid_out.cc:4623
__global__ void set(Number *val, const Number s, const size_type N)
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:259
::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 > &)