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-24.h
Go to the documentation of this file.
1
528 = 0) const override
529 *   {
530 *   static const std::array<Source, 5> sources{
531 *   {Source(Point<dim>(0, 0), 0.025),
532 *   Source(Point<dim>(-0.135, 0), 0.05),
533 *   Source(Point<dim>(0.17, 0), 0.03),
534 *   Source(Point<dim>(-0.25, 0), 0.02),
535 *   Source(Point<dim>(-0.05, -0.15), 0.015)}};
536 *  
537 *   for (const auto &source : sources)
538 *   if (p.distance(source.location) < source.radius)
539 *   return 1;
540 *  
541 *   return 0;
542 *   }
543 *  
544 *   private:
545 *   struct Source
546 *   {
547 *   Source(const Point<dim> &l, const double r)
548 *   : location(l)
549 *   , radius(r)
550 *   {}
551 *  
552 *   const Point<dim> location;
553 *   const double radius;
554 *   };
555 *   };
556 *  
557 *  
558 * @endcode
559 *
560 *
561 * <a name="step_24-ImplementationofthecodeTATForwardProblemcodeclass"></a>
562 * <h3>Implementation of the <code>TATForwardProblem</code> class</h3>
563 *
564
565 *
566 * Let's start again with the constructor. Setting the member variables is
567 * straightforward. We use the acoustic wave speed of mineral oil (in
568 * millimeters per microsecond, a common unit in experimental biomedical
569 * imaging) since this is where many of the experiments we want to compare
570 * the output with are made in. The Crank-Nicolson scheme is used again,
571 * i.e. theta is set to 0.5. The time step is later selected to satisfy @f$k =
572 * \frac hc@f$: here we initialize it to an invalid number.
573 *
574 * @code
575 *   template <int dim>
576 *   TATForwardProblem<dim>::TATForwardProblem()
577 *   : fe(1)
578 *   , dof_handler(triangulation)
579 *   , time_step(std::numeric_limits<double>::quiet_NaN())
580 *   , time(time_step)
581 *   , timestep_number(1)
582 *   , theta(0.5)
583 *   , wave_speed(1.437)
584 *   {
585 * @endcode
586 *
587 * The second task in the constructor is to initialize the array that
588 * holds the detector locations. The results of this program were compared
589 * with experiments in which the step size of the detector spacing is 2.25
590 * degree, corresponding to 160 detector locations. The radius of the
591 * scanning circle is selected to be half way between the center and the
592 * boundary to avoid that the remaining reflections from the imperfect
593 * boundary condition spoils our numerical results.
594 *
595
596 *
597 * The locations of the detectors are then calculated in clockwise
598 * order. Note that the following of course only works if we are computing
599 * in 2d, a condition that we guard with an assertion. If we later wanted
600 * to run the same program in 3d, we would have to add code here for the
601 * initialization of detector locations in 3d. Due to the assertion, there
602 * is no way we can forget to do this.
603 *
604 * @code
605 *   Assert(dim == 2, ExcNotImplemented());
606 *  
607 *   const double detector_step_angle = 2.25;
608 *   const double detector_radius = 0.5;
609 *  
610 *   for (double detector_angle = 2 * numbers::PI; detector_angle >= 0;
611 *   detector_angle -= detector_step_angle / 360 * 2 * numbers::PI)
612 *   detector_locations.push_back(
613 *   Point<dim>(std::cos(detector_angle), std::sin(detector_angle)) *
614 *   detector_radius);
615 *   }
616 *  
617 *  
618 *  
619 * @endcode
620 *
621 *
622 * <a name="step_24-TATForwardProblemsetup_system"></a>
623 * <h4>TATForwardProblem::setup_system</h4>
624 *
625
626 *
627 * The following system is pretty much what we've already done in @ref step_23 "step-23",
628 * but with two important differences. First, we have to create a circular
629 * (or spherical) mesh around the origin, with a radius of 1. This nothing
630 * new: we've done so before in @ref step_6 "step-6" and @ref step_10 "step-10", where we also explain
631 * how the PolarManifold or SphericalManifold object places new points on
632 * concentric circles when a cell is refined, which we will use here as
633 * well.
634 *
635
636 *
637 * One thing we had to make sure is that the time step satisfies the CFL
638 * condition discussed in the introduction of @ref step_23 "step-23". Back in that program,
639 * we ensured this by hand by setting a timestep that matches the mesh
640 * width, but that was error prone because if we refined the mesh once more
641 * we would also have to make sure the time step is changed. Here, we do
642 * that automatically: we ask a library function for the minimal diameter of
643 * any cell. Then we set @f$k=\frac h{c_0}@f$. The only problem is: what exactly
644 * is @f$h@f$? The point is that there is really no good theory on this question
645 * for the wave equation. It is known that for uniformly refined meshes
646 * consisting of rectangles, @f$h@f$ is the minimal edge length. But for meshes
647 * on general quadrilaterals, the exact relationship appears to be unknown,
648 * i.e. it is unknown what properties of cells are relevant for the CFL
649 * condition. The problem is that the CFL condition follows from knowledge
650 * of the smallest eigenvalue of the Laplace matrix, and that can only be
651 * computed analytically for simply structured meshes.
652 *
653
654 *
655 * The upshot of all this is that we're not quite sure what exactly we
656 * should take for @f$h@f$. The function GridTools::minimal_cell_diameter
657 * computes the minimal diameter of all cells. If the cells were all squares
658 * or cubes, then the minimal edge length would be the minimal diameter
659 * divided by <code>std::sqrt(dim)</code>. We simply generalize this,
660 * without theoretical justification, to the case of non-uniform meshes.
661 *
662
663 *
664 * The only other significant change is that we need to build the boundary
665 * mass matrix. We will comment on this further down below.
666 *
667 * @code
668 *   template <int dim>
669 *   void TATForwardProblem<dim>::setup_system()
670 *   {
671 *   const Point<dim> center;
674 *  
675 *   time_step = GridTools::minimal_cell_diameter(triangulation) / wave_speed /
676 *   std::sqrt(1. * dim);
677 *  
678 *   std::cout << "Number of active cells: " << triangulation.n_active_cells()
679 *   << std::endl;
680 *  
681 *   dof_handler.distribute_dofs(fe);
682 *  
683 *   std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
684 *   << std::endl
685 *   << std::endl;
686 *  
687 *   DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
688 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
689 *   sparsity_pattern.copy_from(dsp);
690 *  
691 *   system_matrix.reinit(sparsity_pattern);
692 *   mass_matrix.reinit(sparsity_pattern);
693 *   laplace_matrix.reinit(sparsity_pattern);
694 *  
695 *   MatrixCreator::create_mass_matrix(dof_handler,
696 *   QGauss<dim>(fe.degree + 1),
697 *   mass_matrix);
699 *   QGauss<dim>(fe.degree + 1),
700 *   laplace_matrix);
701 *  
702 * @endcode
703 *
704 * The second difference, as mentioned, to @ref step_23 "step-23" is that we need to
705 * build the boundary mass matrix that grew out of the absorbing boundary
706 * conditions.
707 *
708
709 *
710 * A first observation would be that this matrix is much sparser than the
711 * regular mass matrix, since none of the shape functions with purely
712 * interior support contribute to this matrix. We could therefore
713 * optimize the storage pattern to this situation and build up a second
714 * sparsity pattern that only contains the nonzero entries that we
715 * need. There is a trade-off to make here: first, we would have to have a
716 * second sparsity pattern object, so that costs memory. Secondly, the
717 * matrix attached to this sparsity pattern is going to be smaller and
718 * therefore requires less memory; it would also be faster to perform
719 * matrix-vector multiplications with it. The final argument, however, is
720 * the one that tips the scale: we are not primarily interested in
721 * performing matrix-vector with the boundary matrix alone (though we need
722 * to do that for the right hand side vector once per time step), but
723 * mostly wish to add it up to the other matrices used in the first of the
724 * two equations since this is the one that is going to be multiplied with
725 * once per iteration of the CG method, i.e. significantly more often. It
726 * is now the case that the SparseMatrix::add class allows to add one
727 * matrix to another, but only if they use the same sparsity pattern (the
728 * reason being that we can't add nonzero entries to a matrix after the
729 * sparsity pattern has been created, so we simply require that the two
730 * matrices have the same sparsity pattern).
731 *
732
733 *
734 * So let's go with that:
735 *
736 * @code
737 *   boundary_matrix.reinit(sparsity_pattern);
738 *  
739 * @endcode
740 *
741 * The second thing to do is to actually build the matrix. Here, we need
742 * to integrate over faces of cells, so first we need a quadrature object
743 * that works on <code>dim-1</code> dimensional objects. Secondly, the
744 * FEFaceValues variant of FEValues that works on faces, as its name
745 * suggest. And finally, the other variables that are part of the assembly
746 * machinery. All of this we put between curly braces to limit the scope
747 * of these variables to where we actually need them.
748 *
749
750 *
751 * The actual act of assembling the matrix is then fairly straightforward:
752 * we loop over all cells, over all faces of each of these cells, and then
753 * do something only if that particular face is at the boundary of the
754 * domain. Like this:
755 *
756 * @code
757 *   {
758 *   const QGauss<dim - 1> quadrature_formula(fe.degree + 1);
759 *   FEFaceValues<dim> fe_values(fe,
760 *   quadrature_formula,
762 *  
763 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
764 *   const unsigned int n_q_points = quadrature_formula.size();
765 *  
766 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
767 *  
768 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
769 *  
770 *   for (const auto &cell : dof_handler.active_cell_iterators())
771 *   for (const auto &face : cell->face_iterators())
772 *   if (face->at_boundary())
773 *   {
774 *   cell_matrix = 0;
775 *  
776 *   fe_values.reinit(cell, face);
777 *  
778 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
779 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
780 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
781 *   cell_matrix(i, j) += (fe_values.shape_value(i, q_point) *
782 *   fe_values.shape_value(j, q_point) *
783 *   fe_values.JxW(q_point));
784 *  
785 *   cell->get_dof_indices(local_dof_indices);
786 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
787 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
788 *   boundary_matrix.add(local_dof_indices[i],
789 *   local_dof_indices[j],
790 *   cell_matrix(i, j));
791 *   }
792 *   }
793 *  
794 *   system_matrix.copy_from(mass_matrix);
795 *   system_matrix.add(time_step * time_step * theta * theta * wave_speed *
796 *   wave_speed,
797 *   laplace_matrix);
798 *   system_matrix.add(wave_speed * theta * time_step, boundary_matrix);
799 *  
800 *  
801 *   solution_p.reinit(dof_handler.n_dofs());
802 *   old_solution_p.reinit(dof_handler.n_dofs());
803 *   system_rhs_p.reinit(dof_handler.n_dofs());
804 *  
805 *   solution_v.reinit(dof_handler.n_dofs());
806 *   old_solution_v.reinit(dof_handler.n_dofs());
807 *   system_rhs_v.reinit(dof_handler.n_dofs());
808 *  
809 *   constraints.close();
810 *   }
811 *  
812 *  
813 * @endcode
814 *
815 *
816 * <a name="step_24-TATForwardProblemsolve_pandTATForwardProblemsolve_v"></a>
817 * <h4>TATForwardProblem::solve_p and TATForwardProblem::solve_v</h4>
818 *
819
820 *
821 * The following two functions, solving the linear systems for the pressure
822 * and the velocity variable, are taken pretty much verbatim (with the
823 * exception of the change of name from @f$u@f$ to @f$p@f$ of the primary variable)
824 * from @ref step_23 "step-23":
825 *
826 * @code
827 *   template <int dim>
828 *   void TATForwardProblem<dim>::solve_p()
829 *   {
830 *   SolverControl solver_control(1000, 1e-8 * system_rhs_p.l2_norm());
831 *   SolverCG<Vector<double>> cg(solver_control);
832 *  
833 *   cg.solve(system_matrix, solution_p, system_rhs_p, PreconditionIdentity());
834 *  
835 *   std::cout << " p-equation: " << solver_control.last_step()
836 *   << " CG iterations." << std::endl;
837 *   }
838 *  
839 *  
840 *  
841 *   template <int dim>
842 *   void TATForwardProblem<dim>::solve_v()
843 *   {
844 *   SolverControl solver_control(1000, 1e-8 * system_rhs_v.l2_norm());
845 *   SolverCG<Vector<double>> cg(solver_control);
846 *  
847 *   cg.solve(mass_matrix, solution_v, system_rhs_v, PreconditionIdentity());
848 *  
849 *   std::cout << " v-equation: " << solver_control.last_step()
850 *   << " CG iterations." << std::endl;
851 *   }
852 *  
853 *  
854 *  
855 * @endcode
856 *
857 *
858 * <a name="step_24-TATForwardProblemoutput_results"></a>
859 * <h4>TATForwardProblem::output_results</h4>
860 *
861
862 *
863 * The same holds here: the function is from @ref step_23 "step-23".
864 *
865 * @code
866 *   template <int dim>
867 *   void TATForwardProblem<dim>::output_results() const
868 *   {
869 *   DataOut<dim> data_out;
870 *  
871 *   data_out.attach_dof_handler(dof_handler);
872 *   data_out.add_data_vector(solution_p, "P");
873 *   data_out.add_data_vector(solution_v, "V");
874 *  
875 *   data_out.build_patches();
876 *  
877 *   const std::string filename =
878 *   "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu";
879 *   DataOutBase::VtkFlags vtk_flags;
881 *   std::ofstream output(filename);
882 *   data_out.write_vtu(output);
883 *   }
884 *  
885 *  
886 *  
887 * @endcode
888 *
889 *
890 * <a name="step_24-TATForwardProblemrun"></a>
891 * <h4>TATForwardProblem::run</h4>
892 *
893
894 *
895 * This function that does most of the work is pretty much again like in
896 * @ref step_23 "step-23", though we make things a bit clearer by using the vectors G1 and
897 * G2 mentioned in the introduction. Compared to the overall memory
898 * consumption of the program, the introduction of a few temporary vectors
899 * isn't doing much harm.
900 *
901
902 *
903 * The only changes to this function are: first, that we do not have to
904 * project initial values for the velocity @f$v@f$, since we know that it is
905 * zero. And second that we evaluate the solution at the detector locations
906 * computed in the constructor. This is done using the
907 * VectorTools::point_value function. These values are then written to a
908 * file that we open at the beginning of the function.
909 *
910 * @code
911 *   template <int dim>
912 *   void TATForwardProblem<dim>::run()
913 *   {
914 *   setup_system();
915 *  
916 *   VectorTools::project(dof_handler,
917 *   constraints,
918 *   QGauss<dim>(fe.degree + 1),
919 *   InitialValuesP<dim>(),
920 *   old_solution_p);
921 *   old_solution_v = 0;
922 *  
923 *  
924 *   std::ofstream detector_data("detectors.dat");
925 *  
926 *   Vector<double> tmp(solution_p.size());
927 *   Vector<double> G1(solution_p.size());
928 *   Vector<double> G2(solution_v.size());
929 *  
930 *   const double end_time = 0.7;
931 *   for (time = time_step; time <= end_time;
932 *   time += time_step, ++timestep_number)
933 *   {
934 *   std::cout << std::endl;
935 *   std::cout << "time_step " << timestep_number << " @ t=" << time
936 *   << std::endl;
937 *  
938 *   mass_matrix.vmult(G1, old_solution_p);
939 *   mass_matrix.vmult(tmp, old_solution_v);
940 *   G1.add(time_step * (1 - theta), tmp);
941 *  
942 *   mass_matrix.vmult(G2, old_solution_v);
943 *   laplace_matrix.vmult(tmp, old_solution_p);
944 *   G2.add(-wave_speed * wave_speed * time_step * (1 - theta), tmp);
945 *  
946 *   boundary_matrix.vmult(tmp, old_solution_p);
947 *   G2.add(wave_speed, tmp);
948 *  
949 *   system_rhs_p = G1;
950 *   system_rhs_p.add(time_step * theta, G2);
951 *  
952 *   solve_p();
953 *  
954 *   system_rhs_v = G2;
955 *   laplace_matrix.vmult(tmp, solution_p);
956 *   system_rhs_v.add(-time_step * theta * wave_speed * wave_speed, tmp);
957 *  
958 *   boundary_matrix.vmult(tmp, solution_p);
959 *   system_rhs_v.add(-wave_speed, tmp);
960 *  
961 *   solve_v();
962 *  
963 *   output_results();
964 *  
965 *   detector_data << time;
966 *   for (unsigned int i = 0; i < detector_locations.size(); ++i)
967 *   detector_data << ' '
968 *   << VectorTools::point_value(dof_handler,
969 *   solution_p,
970 *   detector_locations[i])
971 *   << ' ';
972 *   detector_data << std::endl;
973 *  
974 *   old_solution_p = solution_p;
975 *   old_solution_v = solution_v;
976 *   }
977 *   }
978 *   } // namespace Step24
979 *  
980 *  
981 *  
982 * @endcode
983 *
984 *
985 * <a name="step_24-Thecodemaincodefunction"></a>
986 * <h3>The <code>main</code> function</h3>
987 *
988
989 *
990 * What remains is the main function of the program. There is nothing here
991 * that hasn't been shown in several of the previous programs:
992 *
993 * @code
994 *   int main()
995 *   {
996 *   try
997 *   {
998 *   using namespace Step24;
999 *  
1000 *   TATForwardProblem<2> forward_problem_solver;
1001 *   forward_problem_solver.run();
1002 *   }
1003 *   catch (std::exception &exc)
1004 *   {
1005 *   std::cerr << std::endl
1006 *   << std::endl
1007 *   << "----------------------------------------------------"
1008 *   << std::endl;
1009 *   std::cerr << "Exception on processing: " << std::endl
1010 *   << exc.what() << std::endl
1011 *   << "Aborting!" << std::endl
1012 *   << "----------------------------------------------------"
1013 *   << std::endl;
1014 *  
1015 *   return 1;
1016 *   }
1017 *   catch (...)
1018 *   {
1019 *   std::cerr << std::endl
1020 *   << std::endl
1021 *   << "----------------------------------------------------"
1022 *   << std::endl;
1023 *   std::cerr << "Unknown exception!" << std::endl
1024 *   << "Aborting!" << std::endl
1025 *   << "----------------------------------------------------"
1026 *   << std::endl;
1027 *   return 1;
1028 *   }
1029 *  
1030 *   return 0;
1031 *   }
1032 * @endcode
1033<a name="step_24-Results"></a><h1>Results</h1>
1034
1035
1036The program writes both graphical data for each time step as well as the
1037values evaluated at each detector location to disk. We then
1038draw them in plots. Experimental data were also collected for comparison.
1039Currently our experiments have only been done in two dimensions by
1040circularly scanning a single detector. The tissue sample here is a thin slice
1041in the @f$X-Y@f$ plane (@f$Z=0@f$), and we assume that signals from other @f$Z@f$
1042directions won't contribute to the data. Consequently, we only have to compare
1043our experimental data with two dimensional simulated data.
1044
1045<a name="step_24-Oneabsorber"></a><h3> One absorber </h3>
1046
1047
1048This movie shows the thermoacoustic waves generated by a single small absorber
1049propagating in the medium (in our simulation, we assume the medium is mineral
1050oil, which has a acoustic speed of 1.437 @f$\frac{mm}{\mu s}@f$):
1051
1052<img src="https://www.dealii.org/images/steps/developer/step-24.one_movie.gif" alt="">
1053
1054For a single absorber, we of course have to change the
1055<code>InitialValuesP</code> class accordingly.
1056
1057Next, let us compare experimental and computational results. The visualization
1058uses a technique long used in seismology, where the data of each detector is
1059plotted all in one graph. The way this is done is by offsetting each
1060detector's signal a bit compared to the previous one. For example, here is a
1061plot of the first four detectors (from bottom to top, with time in
1062microseconds running from left to right) using the source setup used in the
1063program, to make things a bit more interesting compared to the present case of
1064only a single source:
1065
1066<img src="https://www.dealii.org/images/steps/developer/step-24.traces.png" alt="">
1067
1068One thing that can be seen, for example, is that the arrival of the second and
1069fourth signals shifts to earlier times for greater detector numbers (i.e. the
1070topmost ones), but not the first and the third; this can be interpreted to
1071mean that the origin of these signals must be closer to the latter detectors
1072than to the former ones.
1073
1074If we stack not only 4, but all 160 detectors in one graph, the individual
1075lines blur, but where they run together they create a pattern of darker or
1076lighter grayscales. The following two figures show the results obtained at
1077the detector locations stacked in that way. The left figure is obtained from
1078experiments, and the right is the simulated data.
1079In the experiment, a single small strong absorber was embedded in
1080weaker absorbing tissue:
1081
1082<table width="100%">
1083<tr>
1084<td>
1085<img src="https://www.dealii.org/images/steps/developer/step-24.one.png" alt="">
1086</td>
1087<td>
1088<img src="https://www.dealii.org/images/steps/developer/step-24.one_s.png" alt="">
1089</td>
1090</tr>
1091</table>
1092
1093It is obvious that the source location is closer to the detectors at angle
1094@f$180^\circ@f$. All the other signals that can be seen in the experimental data
1095result from the fact that there are weak absorbers also in the rest of the
1096tissue, which surrounds the signals generated by the small strong absorber in
1097the center. On the other hand, in the simulated data, we only simulate the
1098small strong absorber.
1099
1100In reality, detectors have limited bandwidth. The thermoacoustic waves passing
1101through the detector will therefore be filtered. By using a high-pass filter
1102(implemented in MATLAB and run against the data file produced by this program),
1103the simulated results can be made to look closer to the experimental
1104data:
1105
1106<img src="https://www.dealii.org/images/steps/developer/step-24.one_sf.png" alt="">
1107
1108In our simulations, we see spurious signals behind the main wave that
1109result from numerical artifacts. This problem can be alleviated by using finer
1110mesh, resulting in the following plot:
1111
1112<img src="https://www.dealii.org/images/steps/developer/step-24.one_s2.png" alt="">
1113
1114
1115
1116<a name="step_24-Multipleabsorbers"></a><h3>Multiple absorbers</h3>
1117
1118
1119To further verify the program, we will also show simulation results for
1120multiple absorbers. This corresponds to the case that is actually implemented
1121in the program. The following movie shows the propagation of the generated
1122thermoacoustic waves in the medium by multiple absorbers:
1123
1124<img src="https://www.dealii.org/images/steps/developer/step-24.multi_movie.gif" alt="">
1125
1126Experimental data and our simulated data are compared in the following two
1127figures:
1128<table width="100%">
1129<tr>
1130<td>
1131<img src="https://www.dealii.org/images/steps/developer/step-24.multi.png" alt="">
1132</td>
1133<td>
1134<img src="https://www.dealii.org/images/steps/developer/step-24.multi_s.png" alt="">
1135</td>
1136</tr>
1137</table>
1138
1139Note that in the experimental data, the first signal (i.e. the left-most dark
1140line) results from absorption at the tissue boundary, and therefore reaches
1141the detectors first and before any of the signals from the interior. This
1142signal is also faintly visible at the end of the traces, around 30 @f$\mu s@f$,
1143which indicates that the signal traveled through the entire tissue to reach
1144detectors at the other side, after all the signals originating from the
1145interior have reached them.
1146
1147As before, the numerical result better matches experimental ones by applying a
1148bandwidth filter that matches the actual behavior of detectors (left) and by
1149choosing a finer mesh (right):
1150
1151<table width="100%">
1152<tr>
1153<td>
1154<img src="https://www.dealii.org/images/steps/developer/step-24.multi_sf.png" alt="">
1155</td>
1156<td>
1157<img src="https://www.dealii.org/images/steps/developer/step-24.multi_s2.png" alt="">
1158</td>
1159</tr>
1160</table>
1161
1162One of the important differences between the left and the right figure is that
1163the curves look much less "angular" at the right. The angularity comes from
1164the fact that while waves in the continuous equation travel equally fast in
1165all directions, this isn't the case after discretization: there, waves that
1166travel diagonal to cells move at slightly different speeds to those that move
1167parallel to mesh lines. This anisotropy leads to wave fronts that aren't
1168perfectly circular (and would produce sinusoidal signals in the stacked
1169plots), but are bulged out in certain directions. To make things worse, the
1170circular mesh we use (see for example @ref step_6 "step-6" for a view of the
1171coarse mesh) is not isotropic either. The net result is that the signal fronts
1172are not sinusoidal unless the mesh is sufficiently fine. The right image is a
1173lot better in this respect, though artifacts in the form of trailing spurious
1174waves can still be seen.
1175 *
1176 *
1177<a name="step_24-PlainProg"></a>
1178<h1> The plain program</h1>
1179@include "step-24.cc"
1180*/
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
Definition point.h:111
void add(const size_type i, const size_type j, const number value)
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
Point< 3 > center
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
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.
void hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), 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.
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:470
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
DataOutBase::CompressionLevel compression_level