deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+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-24.h
Go to the documentation of this file.
1 = 0) const override
531 *   {
532 *   static const std::array<Source, 5> sources{
533 *   {Source(Point<dim>(0, 0), 0.025),
534 *   Source(Point<dim>(-0.135, 0), 0.05),
535 *   Source(Point<dim>(0.17, 0), 0.03),
536 *   Source(Point<dim>(-0.25, 0), 0.02),
537 *   Source(Point<dim>(-0.05, -0.15), 0.015)}};
538 *  
539 *   for (const auto &source : sources)
540 *   if (p.distance(source.location) < source.radius)
541 *   return 1;
542 *  
543 *   return 0;
544 *   }
545 *  
546 *   private:
547 *   struct Source
548 *   {
549 *   Source(const Point<dim> &l, const double r)
550 *   : location(l)
551 *   , radius(r)
552 *   {}
553 *  
554 *   const Point<dim> location;
555 *   const double radius;
556 *   };
557 *   };
558 *  
559 *  
560 * @endcode
561 *
562 *
563 * <a name="step_24-ImplementationofthecodeTATForwardProblemcodeclass"></a>
564 * <h3>Implementation of the <code>TATForwardProblem</code> class</h3>
565 *
566
567 *
568 * Let's start again with the constructor. Setting the member variables is
569 * straightforward. We use the acoustic wave speed of mineral oil (in
570 * millimeters per microsecond, a common unit in experimental biomedical
571 * imaging) since this is where many of the experiments we want to compare
572 * the output with are made in. The Crank-Nicolson scheme is used again,
573 * i.e. theta is set to 0.5. The time step is later selected to satisfy @f$k =
574 * \frac hc@f$: here we initialize it to an invalid number.
575 *
576 * @code
577 *   template <int dim>
578 *   TATForwardProblem<dim>::TATForwardProblem()
579 *   : fe(1)
580 *   , dof_handler(triangulation)
581 *   , time_step(std::numeric_limits<double>::quiet_NaN())
582 *   , time(time_step)
583 *   , timestep_number(1)
584 *   , theta(0.5)
585 *   , wave_speed(1.437)
586 *   {
587 * @endcode
588 *
589 * The second task in the constructor is to initialize the array that
590 * holds the detector locations. The results of this program were compared
591 * with experiments in which the step size of the detector spacing is 2.25
592 * degree, corresponding to 160 detector locations. The radius of the
593 * scanning circle is selected to be half way between the center and the
594 * boundary to avoid that the remaining reflections from the imperfect
595 * boundary condition spoils our numerical results.
596 *
597
598 *
599 * The locations of the detectors are then calculated in clockwise
600 * order. Note that the following of course only works if we are computing
601 * in 2d, a condition that we guard with an assertion. If we later wanted
602 * to run the same program in 3d, we would have to add code here for the
603 * initialization of detector locations in 3d. Due to the assertion, there
604 * is no way we can forget to do this.
605 *
606 * @code
607 *   Assert(dim == 2, ExcNotImplemented());
608 *  
609 *   const double detector_step_angle = 2.25;
610 *   const double detector_radius = 0.5;
611 *  
612 *   for (double detector_angle = 2 * numbers::PI; detector_angle >= 0;
613 *   detector_angle -= detector_step_angle / 360 * 2 * numbers::PI)
614 *   detector_locations.push_back(
615 *   Point<dim>(std::cos(detector_angle), std::sin(detector_angle)) *
616 *   detector_radius);
617 *   }
618 *  
619 *  
620 *  
621 * @endcode
622 *
623 *
624 * <a name="step_24-TATForwardProblemsetup_system"></a>
625 * <h4>TATForwardProblem::setup_system</h4>
626 *
627
628 *
629 * The following system is pretty much what we've already done in @ref step_23 "step-23",
631 * (or spherical) mesh around the origin, with a radius of 1. This nothing
632 * new: we've done so before in @ref step_6 "step-6" and @ref step_10 "step-10", where we also explain
633 * how the PolarManifold or SphericalManifold object places new points on
634 * concentric circles when a cell is refined, which we will use here as
635 * well.
636 *
637
638 *
639 * One thing we had to make sure is that the time step satisfies the CFL
640 * condition discussed in the introduction of @ref step_23 "step-23". Back in that program,
641 * we ensured this by hand by setting a timestep that matches the mesh
642 * width, but that was error prone because if we refined the mesh once more
643 * we would also have to make sure the time step is changed. Here, we do
644 * that automatically: we ask a library function for the minimal diameter of
645 * any cell. Then we set @f$k=\frac h{c_0}@f$. The only problem is: what exactly
646 * is @f$h@f$? The point is that there is really no good theory on this question
647 * for the wave equation. It is known that for uniformly refined meshes
648 * consisting of rectangles, @f$h@f$ is the minimal edge length. But for meshes
649 * on general quadrilaterals, the exact relationship appears to be unknown,
650 * i.e. it is unknown what properties of cells are relevant for the CFL
651 * condition. The problem is that the CFL condition follows from knowledge
652 * of the smallest eigenvalue of the Laplace matrix, and that can only be
653 * computed analytically for simply structured meshes.
654 *
655
656 *
657 * The upshot of all this is that we're not quite sure what exactly we
659 * computes the minimal diameter of all cells. If the cells were all squares
663 *
664
665 *
666 * The only other significant change is that we need to build the boundary
668 *
669 * @code
670 *   template <int dim>
672 *   {
673 *   const Point<dim> center;
675 *   triangulation.refine_global(7);
676 *  
678 *   std::sqrt(1. * dim);
679 *  
680 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
681 *   << std::endl;
682 *  
683 *   dof_handler.distribute_dofs(fe);
684 *  
685 *   std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
686 *   << std::endl
687 *   << std::endl;
688 *  
689 *   DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
690 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
691 *   sparsity_pattern.copy_from(dsp);
692 *  
693 *   system_matrix.reinit(sparsity_pattern);
694 *   mass_matrix.reinit(sparsity_pattern);
695 *   laplace_matrix.reinit(sparsity_pattern);
696 *  
698 *   QGauss<dim>(fe.degree + 1),
699 *   mass_matrix);
701 *   QGauss<dim>(fe.degree + 1),
703 *  
704 * @endcode
705 *
706 * The second difference, as mentioned, to @ref step_23 "step-23" is that we need to
707 * build the boundary mass matrix that grew out of the absorbing boundary
708 * conditions.
709 *
710
711 *
714 * interior support contribute to this matrix. We could therefore
715 * optimize the storage pattern to this situation and build up a second
716 * sparsity pattern that only contains the nonzero entries that we
718 * second sparsity pattern object, so that costs memory. Secondly, the
719 * matrix attached to this sparsity pattern is going to be smaller and
720 * therefore requires less memory; it would also be faster to perform
723 * performing matrix-vector with the boundary matrix alone (though we need
724 * to do that for the right hand side vector once per time step), but
725 * mostly wish to add it up to the other matrices used in the first of the
727 * once per iteration of the CG method, i.e. significantly more often. It
728 * is now the case that the SparseMatrix::add class allows to add one
729 * matrix to another, but only if they use the same sparsity pattern (the
730 * reason being that we can't add nonzero entries to a matrix after the
731 * sparsity pattern has been created, so we simply require that the two
732 * matrices have the same sparsity pattern).
733 *
734
735 *
736 * So let's go with that:
737 *
738 * @code
739 *   boundary_matrix.reinit(sparsity_pattern);
740 *  
741 * @endcode
742 *
744 * to integrate over faces of cells, so first we need a quadrature object
745 * that works on <code>dim-1</code> dimensional objects. Secondly, the
746 * FEFaceValues variant of FEValues that works on faces, as its name
750 *
751
752 *
754 * we loop over all cells, over all faces of each of these cells, and then
755 * do something only if that particular face is at the boundary of the
756 * domain. Like this:
757 *
758 * @code
759 *   {
760 *   const QGauss<dim - 1> quadrature_formula(fe.degree + 1);
761 *   FEFaceValues<dim> fe_values(fe,
764 *  
765 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
766 *   const unsigned int n_q_points = quadrature_formula.size();
767 *  
768 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
769 *  
770 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
771 *  
772 *   for (const auto &cell : dof_handler.active_cell_iterators())
773 *   for (const auto &face : cell->face_iterators())
774 *   if (face->at_boundary())
775 *   {
776 *   cell_matrix = 0;
777 *  
778 *   fe_values.reinit(cell, face);
779 *  
780 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
781 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
782 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
783 *   cell_matrix(i, j) += (fe_values.shape_value(i, q_point) *
784 *   fe_values.shape_value(j, q_point) *
785 *   fe_values.JxW(q_point));
786 *  
787 *   cell->get_dof_indices(local_dof_indices);
788 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
789 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
790 *   boundary_matrix.add(local_dof_indices[i],
791 *   local_dof_indices[j],
792 *   cell_matrix(i, j));
793 *   }
794 *   }
795 *  
796 *   system_matrix.copy_from(mass_matrix);
797 *   system_matrix.add(time_step * time_step * theta * theta * wave_speed *
798 *   wave_speed,
800 *   system_matrix.add(wave_speed * theta * time_step, boundary_matrix);
801 *  
802 *  
803 *   solution_p.reinit(dof_handler.n_dofs());
804 *   old_solution_p.reinit(dof_handler.n_dofs());
805 *   system_rhs_p.reinit(dof_handler.n_dofs());
806 *  
807 *   solution_v.reinit(dof_handler.n_dofs());
808 *   old_solution_v.reinit(dof_handler.n_dofs());
809 *   system_rhs_v.reinit(dof_handler.n_dofs());
810 *  
811 *   constraints.close();
812 *   }
813 *  
814 *  
815 * @endcode
816 *
817 *
818 * <a name="step_24-TATForwardProblemsolve_pandTATForwardProblemsolve_v"></a>
820 *
821
822 *
825 * exception of the change of name from @f$u@f$ to @f$p@f$ of the primary variable)
826 * from @ref step_23 "step-23":
827 *
828 * @code
829 *   template <int dim>
831 *   {
832 *   SolverControl solver_control(1000, 1e-8 * system_rhs_p.l2_norm());
833 *   SolverCG<Vector<double>> cg(solver_control);
834 *  
835 *   cg.solve(system_matrix, solution_p, system_rhs_p, PreconditionIdentity());
836 *  
837 *   std::cout << " p-equation: " << solver_control.last_step()
838 *   << " CG iterations." << std::endl;
839 *   }
840 *  
841 *  
842 *  
843 *   template <int dim>
845 *   {
846 *   SolverControl solver_control(1000, 1e-8 * system_rhs_v.l2_norm());
847 *   SolverCG<Vector<double>> cg(solver_control);
848 *  
849 *   cg.solve(mass_matrix, solution_v, system_rhs_v, PreconditionIdentity());
850 *  
851 *   std::cout << " v-equation: " << solver_control.last_step()
852 *   << " CG iterations." << std::endl;
853 *   }
854 *  
855 *  
856 *  
857 * @endcode
858 *
859 *
860 * <a name="step_24-TATForwardProblemoutput_results"></a>
862 *
863
864 *
865 * The same holds here: the function is from @ref step_23 "step-23".
866 *
867 * @code
868 *   template <int dim>
870 *   {
871 *   DataOut<dim> data_out;
872 *  
873 *   data_out.attach_dof_handler(dof_handler);
874 *   data_out.add_data_vector(solution_p, "P");
875 *   data_out.add_data_vector(solution_v, "V");
876 *  
877 *   data_out.build_patches();
878 *  
879 *   const std::string filename =
880 *   "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu";
881 *   DataOutBase::VtkFlags vtk_flags;
882 *   vtk_flags.compression_level = DataOutBase::CompressionLevel::best_speed;
883 *   std::ofstream output(filename);
884 *   data_out.write_vtu(output);
885 *   }
886 *  
887 *  
888 *  
889 * @endcode
890 *
891 *
892 * <a name="step_24-TATForwardProblemrun"></a>
893 * <h4>TATForwardProblem::run</h4>
894 *
895
896 *
897 * This function that does most of the work is pretty much again like in
898 * @ref step_23 "step-23", though we make things a bit clearer by using the vectors G1 and
901 * isn't doing much harm.
902 *
903
904 *
905 * The only changes to this function are: first, that we do not have to
906 * project initial values for the velocity @f$v@f$, since we know that it is
907 * zero. And second that we evaluate the solution at the detector locations
908 * computed in the constructor. This is done using the
909 * VectorTools::point_value function. These values are then written to a
910 * file that we open at the beginning of the function.
911 *
912 * @code
913 *   template <int dim>
914 *   void TATForwardProblem<dim>::run()
915 *   {
916 *   setup_system();
917 *  
918 *   VectorTools::project(dof_handler,
919 *   constraints,
920 *   QGauss<dim>(fe.degree + 1),
921 *   InitialValuesP<dim>(),
922 *   old_solution_p);
923 *   old_solution_v = 0;
924 *  
925 *  
926 *   std::ofstream detector_data("detectors.dat");
927 *  
928 *   Vector<double> tmp(solution_p.size());
929 *   Vector<double> G1(solution_p.size());
930 *   Vector<double> G2(solution_v.size());
931 *  
932 *   const double end_time = 0.7;
933 *   for (time = time_step; time <= end_time;
934 *   time += time_step, ++timestep_number)
935 *   {
936 *   std::cout << std::endl;
937 *   std::cout << "time_step " << timestep_number << " @ t=" << time
938 *   << std::endl;
939 *  
940 *   mass_matrix.vmult(G1, old_solution_p);
941 *   mass_matrix.vmult(tmp, old_solution_v);
942 *   G1.add(time_step * (1 - theta), tmp);
943 *  
944 *   mass_matrix.vmult(G2, old_solution_v);
945 *   laplace_matrix.vmult(tmp, old_solution_p);
946 *   G2.add(-wave_speed * wave_speed * time_step * (1 - theta), tmp);
947 *  
948 *   boundary_matrix.vmult(tmp, old_solution_p);
949 *   G2.add(wave_speed, tmp);
950 *  
951 *   system_rhs_p = G1;
952 *   system_rhs_p.add(time_step * theta, G2);
953 *  
954 *   solve_p();
955 *  
956 *   system_rhs_v = G2;
957 *   laplace_matrix.vmult(tmp, solution_p);
958 *   system_rhs_v.add(-time_step * theta * wave_speed * wave_speed, tmp);
959 *  
960 *   boundary_matrix.vmult(tmp, solution_p);
961 *   system_rhs_v.add(-wave_speed, tmp);
962 *  
963 *   solve_v();
964 *  
965 *   output_results();
966 *  
967 *   detector_data << time;
968 *   for (unsigned int i = 0; i < detector_locations.size(); ++i)
969 *   detector_data << ' '
970 *   << VectorTools::point_value(dof_handler,
971 *   solution_p,
972 *   detector_locations[i])
973 *   << ' ';
974 *   detector_data << std::endl;
975 *  
976 *   old_solution_p = solution_p;
977 *   old_solution_v = solution_v;
978 *   }
979 *   }
980 *   } // namespace Step24
981 *  
982 *  
983 *  
984 * @endcode
985 *
986 *
987 * <a name="step_24-Thecodemaincodefunction"></a>
988 * <h3>The <code>main</code> function</h3>
989 *
990
991 *
992 * What remains is the main function of the program. There is nothing here
993 * that hasn't been shown in several of the previous programs:
994 *
995 * @code
996 *   int main()
997 *   {
998 *   try
999 *   {
1000 *   using namespace Step24;
1001 *  
1003 *   forward_problem_solver.run();
1004 *   }
1005 *   catch (std::exception &exc)
1006 *   {
1007 *   std::cerr << std::endl
1008 *   << std::endl
1009 *   << "----------------------------------------------------"
1010 *   << std::endl;
1011 *   std::cerr << "Exception on processing: " << std::endl
1012 *   << exc.what() << std::endl
1013 *   << "Aborting!" << std::endl
1014 *   << "----------------------------------------------------"
1015 *   << std::endl;
1016 *  
1017 *   return 1;
1018 *   }
1019 *   catch (...)
1020 *   {
1021 *   std::cerr << std::endl
1022 *   << std::endl
1023 *   << "----------------------------------------------------"
1024 *   << std::endl;
1025 *   std::cerr << "Unknown exception!" << std::endl
1026 *   << "Aborting!" << std::endl
1027 *   << "----------------------------------------------------"
1028 *   << std::endl;
1029 *   return 1;
1030 *   }
1031 *  
1032 *   return 0;
1033 *   }
1034 * @endcode
1035<a name="step_24-Results"></a><h1>Results</h1>
1036
1037
1038The program writes both graphical data for each time step as well as the
1041Currently our experiments have only been done in two dimensions by
1043in the @f$X-Y@f$ plane (@f$Z=0@f$), and we assume that signals from other @f$Z@f$
1044directions won't contribute to the data. Consequently, we only have to compare
1045our experimental data with two dimensional simulated data.
1046
1047<a name="step_24-Oneabsorber"></a><h3> One absorber </h3>
1048
1049
1050This movie shows the thermoacoustic waves generated by a single small absorber
1051propagating in the medium (in our simulation, we assume the medium is mineral
1052oil, which has a acoustic speed of 1.437 @f$\frac{mm}{\mu s}@f$):
1053
1054<img src="https://www.dealii.org/images/steps/developer/step-24.one_movie.gif" alt="">
1055
1056For a single absorber, we of course have to change the
1057<code>InitialValuesP</code> class accordingly.
1058
1059Next, let us compare experimental and computational results. The visualization
1060uses a technique long used in seismology, where the data of each detector is
1061plotted all in one graph. The way this is done is by offsetting each
1064microseconds running from left to right) using the source setup used in the
1067
1068<img src="https://www.dealii.org/images/steps/developer/step-24.traces.png" alt="">
1069
1071fourth signals shifts to earlier times for greater detector numbers (i.e. the
1075
1076If we stack not only 4, but all 160 detectors in one graph, the individual
1077lines blur, but where they run together they create a pattern of darker or
1079the detector locations stacked in that way. The left figure is obtained from
1083
1084<table width="100%">
1085<tr>
1086<td>
1087<img src="https://www.dealii.org/images/steps/developer/step-24.one.png" alt="">
1088</td>
1089<td>
1090<img src="https://www.dealii.org/images/steps/developer/step-24.one_s.png" alt="">
1091</td>
1092</tr>
1093</table>
1094
1096@f$180^\circ@f$. All the other signals that can be seen in the experimental data
1101
1103through the detector will therefore be filtered. By using a high-pass filter
1105the simulated results can be made to look closer to the experimental
1106data:
1107
1108<img src="https://www.dealii.org/images/steps/developer/step-24.one_sf.png" alt="">
1109
1113
1114<img src="https://www.dealii.org/images/steps/developer/step-24.one_s2.png" alt="">
1115
1116
1117
1118<a name="step_24-Multipleabsorbers"></a><h3>Multiple absorbers</h3>
1119
1120
1121To further verify the program, we will also show simulation results for
1125
1126<img src="https://www.dealii.org/images/steps/developer/step-24.multi_movie.gif" alt="">
1127
1129figures:
1130<table width="100%">
1131<tr>
1132<td>
1133<img src="https://www.dealii.org/images/steps/developer/step-24.multi.png" alt="">
1134</td>
1135<td>
1136<img src="https://www.dealii.org/images/steps/developer/step-24.multi_s.png" alt="">
1137</td>
1138</tr>
1139</table>
1140
1141Note that in the experimental data, the first signal (i.e. the left-most dark
1142line) results from absorption at the tissue boundary, and therefore reaches
1143the detectors first and before any of the signals from the interior. This
1146detectors at the other side, after all the signals originating from the
1147interior have reached them.
1148
1150bandwidth filter that matches the actual behavior of detectors (left) and by
1151choosing a finer mesh (right):
1152
1153<table width="100%">
1154<tr>
1155<td>
1156<img src="https://www.dealii.org/images/steps/developer/step-24.multi_sf.png" alt="">
1157</td>
1158<td>
1159<img src="https://www.dealii.org/images/steps/developer/step-24.multi_s2.png" alt="">
1160</td>
1161</tr>
1162</table>
1163
1165the curves look much less "angular" at the right. The angularity comes from
1167all directions, this isn't the case after discretization: there, waves that
1168travel diagonal to cells move at slightly different speeds to those that move
1169parallel to mesh lines. This anisotropy leads to wave fronts that aren't
1172circular mesh we use (see for example @ref step_6 "step-6" for a view of the
1177 *
1178 *
1179<a name="step_24-PlainProg"></a>
1180<h1> The plain program</h1>
1181@include "step-24.cc"
1182*/
Definition point.h:113
void add(const size_type i, const size_type j, const number value)
friend class Tensor
Definition tensor.h:865
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
Definition tensor.h:851
Point< 2 > second
Definition grid_out.cc:4630
Point< 2 > first
Definition grid_out.cc:4629
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
std::vector< index_type > data
Definition mpi.cc:740
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > &center={}, const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
double diameter(const Triangulation< dim, spacedim > &tria)
@ matrix
Contents is actually a matrix.
constexpr char A
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57
void create_mass_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, MatrixType &matrix, const Function< spacedim, typename MatrixType::value_type > *const a=nullptr, const AffineConstraints< typename MatrixType::value_type > &constraints=AffineConstraints< typename MatrixType::value_type >())
void create_laplace_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, MatrixType &matrix, const Function< spacedim, typename MatrixType::value_type > *const a=nullptr, const AffineConstraints< typename MatrixType::value_type > &constraints=AffineConstraints< typename MatrixType::value_type >())
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
VectorType::value_type * end(VectorType &V)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:474
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &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)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation