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