415 * vector_value_list(
const std::vector<
Point<dim>> &points,
420 *
for (
unsigned int p = 0; p < points.size(); ++p)
421 * DirichletBoundaryValues<dim>::vector_value(points[p], value_list[p]);
428 * <a name=
"step_29-ThecodeParameterReadercodeclass"></a>
429 * <h3>The <code>ParameterReader</code>
class</h3>
434 * and reading parameters from an input file. It includes a function
435 * <code>declare_parameters</code> that declares all the necessary
436 * parameters and a <code>read_parameters</code> function that is called
437 * from
outside to initiate the parameter reading process.
444 *
void read_parameters(
const std::string &);
447 *
void declare_parameters();
458 * : prm(paramhandler)
464 * <a name=
"step_29-codeParameterReaderdeclare_parameterscode"></a>
465 * <h4><code>ParameterReader::declare_parameters</code></h4>
469 * The <code>declare_parameters</code> function declares all the parameters
471 * along with their
types, range conditions and the subsections they appear
472 * in. We will wrap all the entries that go into a section in a pair of
473 * braces to force the editor to indent them by one
level, making it simpler
474 * to read which entries together form a section:
477 *
void ParameterReader::declare_parameters()
481 * Parameters
for mesh and geometry include the number of global
482 * refinement steps that are applied to the
initial coarse mesh and the
483 * focal distance @f$d@f$ of the transducer lens. For the number of refinement
484 * steps, we allow integer
values in the range @f$[0,\infty)@f$, where the
486 * half-open interval. For the focal distance any number greater than
490 * prm.enter_subsection(
"Mesh & geometry parameters");
492 * prm.declare_entry(
"Number of refinements",
495 *
"Number of global mesh refinement steps "
496 *
"applied to initial coarse grid");
498 * prm.declare_entry(
"Focal distance",
501 *
"Distance of the focal point of the lens "
504 * prm.leave_subsection();
508 * The next subsection is devoted to the physical parameters appearing in
509 * the equation, which are the frequency @f$\omega@f$ and wave speed
510 * @f$c@f$. Again, both need to lie in the half-open interval @f$[0,\infty)@f$
512 * end-point as argument:
515 * prm.enter_subsection(
"Physical constants");
521 * prm.leave_subsection();
526 * Last but not least we would like to be able to change some properties
527 * of the output, like filename and format, through entries in the
528 * configuration file, which is the purpose of the last subsection:
531 * prm.enter_subsection(
"Output parameters");
533 * prm.declare_entry(
"Output filename",
536 *
"Name of the output file (without extension)");
540 * Since different output formats may require different parameters
for
541 * generating output (like
for example, postscript output needs
542 * viewpoint angles, line widths, colors etc), it would be cumbersome
if
543 * we had to declare all these parameters by hand
for every possible
544 * output format supported in the library. Instead, each output format
545 * has a <code>FormatFlags::declare_parameters</code> function, which
546 * declares all the parameters specific to that format in an own
547 * subsection. The following call of
549 * <code>declare_parameters</code>
for all available output formats, so
550 * that
for each format an own subsection will be created with
551 * parameters declared
for that particular output format. (The actual
552 *
value of the
template parameter in the call, <code>@<1@></code>
553 * above, does not matter here: the function does the same work
554 *
independent of the dimension, but happens to be in a
555 *
template-parameter-dependent
class.) To find out what parameters
556 * there are
for which output format, you can either consult the
557 * documentation of the
DataOutBase class, or simply
run this program
558 * without a parameter file present. It will then create a file with all
559 * declared parameters set to their
default values, which can
560 * conveniently serve as a starting
point for setting the parameters to
566 * prm.leave_subsection();
572 * <a name=
"step_29-codeParameterReaderread_parameterscode"></a>
573 * <h4><code>ParameterReader::read_parameters</code></h4>
577 * This is the main function in the ParameterReader
class. It gets called
578 * from
outside,
first declares all the parameters, and then reads them from
579 * the input file whose filename is provided by the caller. After the call
580 * to
this function is complete, the <code>prm</code>
object can be used to
581 * retrieve the
values of the parameters read in from the file :
584 *
void ParameterReader::read_parameters(
const std::string ¶meter_file)
586 * declare_parameters();
588 * prm.parse_input(parameter_file);
596 * <a name=
"step_29-ThecodeComputeIntensitycodeclass"></a>
597 * <h3>The <code>ComputeIntensity</code>
class</h3>
601 * As mentioned in the introduction, the quantity that we are really after
602 * is the spatial distribution of the intensity of the ultrasound wave,
603 * which corresponds to @f$|u|=\
sqrt{v^2+
w^2}@f$. Now we could just be content
604 * with having @f$v@f$ and @f$w@f$ in our output, and use a suitable visualization
605 * or postprocessing tool to derive @f$|u|@f$ from the solution we
606 * computed. However, there is also a way to output data derived from the
607 * solution in deal.II, and we are going to make use of
this mechanism here.
612 * vectors containing output data to a
DataOut object. There is a special
613 * version of
this function that in addition to the data vector has an
615 * function is used
for output is that at each
point where output data is to
619 * quantities from the values, the gradients and the
second derivatives of
620 * the finite element function represented by the data vector (in the case
621 * of face related data, normal vectors are available as well). Hence, this
622 * allows us to output any quantity that can locally be derived from the
623 * values of the solution and its derivatives. Of course, the ultrasound
624 * intensity @f$|u|@f$ is such a quantity and its computation doesn't even
625 * involve any derivatives of @f$v@f$ or @f$w@f$.
630 * this functionality, and we need to derive our own class from it in order
631 * to implement the functions specified by the interface. In the most
632 * general case one has to implement several member functions but if the
633 * output quantity is a single scalar then some of this boilerplate code can
635 * can derive from that one instead. This is what the
636 * <code>ComputeIntensity</code> class does:
643 * ComputeIntensity();
647 * std::vector<
Vector<double>> &computed_quantities)
const override;
652 * In the constructor, we need to call the constructor of the base
class
653 * with two arguments. The
first denotes the name by which the single
scalar
654 * quantity computed by
this class should be represented in output files. In
655 * our
case, the postprocessor has @f$|u|@f$ as output, so we use
"Intensity".
659 * The
second argument is a set of flags that indicate which data is needed
660 * by the postprocessor in order to compute the output quantities. This can
663 * documented in
UpdateFlags. Of course, computation of the derivatives
664 *
requires additional resources, so only the flags
for data that are really
665 * needed should be given here, just as we
do when we use
FEValues objects.
666 * In our
case, only the function values of @f$v@f$ and @f$w@f$ are needed to
667 * compute @f$|u|@f$, so we
're good with the update_values flag.
671 * ComputeIntensity<dim>::ComputeIntensity()
672 * : DataPostprocessorScalar<dim>("Intensity", update_values)
678 * The actual postprocessing happens in the following function. Its input is
679 * an object that stores values of the function (which is here vector-valued)
680 * representing the data vector given to DataOut::add_data_vector, evaluated
681 * at all evaluation points where we generate output, and some tensor objects
682 * representing derivatives (that we don't use here since @f$|u|@f$ is computed
683 * from just @f$v@f$ and @f$w@f$). The derived quantities are returned in the
684 * <code>computed_quantities</code> vector. Remember that
this function may
685 * only use data
for which the respective update flag is specified by
686 * <code>get_needed_update_flags</code>. For example, we may not use the
687 * derivatives here, since our implementation of
688 * <code>get_needed_update_flags</code> requests that only function
values
693 *
void ComputeIntensity<dim>::evaluate_vector_field(
697 *
AssertDimension(computed_quantities.size(), inputs.solution_values.size());
701 * The computation itself is straightforward: We iterate over each
702 * entry in the output vector and compute @f$|u|@f$ from the
703 * corresponding
values of @f$v@f$ and @f$w@f$. We
do this by creating a
704 * complex number @f$u@f$ and then calling `
std::abs()` on the
705 * result. (One may be tempted to call `std::norm()`, but in a
706 * historical quirk, the
C++ committee decided that `std::norm()`
707 * should
return the <i>square</i> of the absolute
value --
708 * thereby not satisfying the properties mathematicians require of
709 * something called a
"norm".)
712 *
for (
unsigned int p = 0; p < computed_quantities.size(); ++p)
717 *
const std::complex<double> u(inputs.solution_values[p](0),
718 * inputs.solution_values[p](1));
720 * computed_quantities[p](0) =
std::abs(u);
728 * <a name=
"step_29-ThecodeUltrasoundProblemcodeclass"></a>
729 * <h3>The <code>UltrasoundProblem</code>
class</h3>
733 * Finally here is the main
class of this program. It's member
functions
734 * are very similar to the previous examples, in particular @ref step_4
"step-4", and the
735 * list of member variables does not contain any major surprises either.
737 * as a reference to allow easy access to the parameters from all
functions
738 * of the
class. Since we are working with vector valued finite elements,
739 * the FE
object we are
using is of type
FESystem.
743 *
class UltrasoundProblem
751 *
void setup_system();
752 *
void assemble_system();
754 *
void output_results()
const;
772 * reference. It also initializes the DoF-Handler and the finite element
773 * system, which consists of two copies of the
scalar @f$Q_1@f$ field, one
for
774 * @f$v@f$ and one
for @f$w@f$. In other words, we want the finite element space
775 * @f$Q_1\times Q_1 = Q_1^2@f$, which is easily constructed and passed as the
776 * constructor argument to the
FESystem class (i.e., the type of the `fe`
777 * member being initialized here):
784 * , fe(
FE_Q<dim>(1) ^ 2)
790 * <a name=
"step_29-codeUltrasoundProblemmake_gridcode"></a>
791 * <h4><code>UltrasoundProblem::make_grid</code></h4>
795 * Here we set up the grid for our domain. As mentioned in the exposition,
796 * the geometry is just a unit square (in 2
d) with the part of the boundary
797 * that represents the transducer lens replaced by a sector of a circle.
801 * void UltrasoundProblem<dim>::make_grid()
805 * First we generate some logging output and start a timer so we can
806 * compute execution time when this function is done:
809 *
std::cout <<
"Generating grid... ";
814 * Then we query the
values for the focal distance of the transducer lens
821 *
const double focal_distance = prm.get_double(
"Focal distance");
822 *
const unsigned int n_refinements = prm.get_integer(
"Number of refinements");
824 * prm.leave_subsection();
828 * Next, two points are defined
for position and focal
point of the
829 * transducer lens, which is the
center of the circle whose segment will
830 * form the transducer part of the boundary. Notice that
this is the only
831 *
point in the program where things are slightly different in 2
d and 3
d.
832 * Even though
this tutorial only deals with the 2
d case, the necessary
833 * additions to make
this program functional in 3
d are so minimal that we
834 * opt
for including them:
846 * As
initial coarse grid we take a simple unit square with 5 subdivisions
847 * in each direction. The number of subdivisions is chosen so that the
848 * line segment @f$[0.4,0.6]@f$ that we want to designate as the transducer
849 * boundary is spanned by a single face. Then we step through all cells to
850 * find the faces where the transducer is to be located, which in fact is
851 * just the single edge from 0.4 to 0.6 on the x-axis. This is where we
852 * want the refinements to be made according to a circle shaped boundary,
853 * so we mark
this edge with a different manifold indicator. Since we will
854 * Dirichlet boundary conditions on the transducer, we also change its
855 * boundary indicator.
861 * for (const auto &face : cell->face_iterators())
862 * if (face->at_boundary() &&
863 * ((face->
center() - transducer).norm_square() < 0.01))
865 * face->set_boundary_id(1);
866 * face->set_manifold_id(1);
871 * is used (which, of course, in 2d just represents a circle), with
center
879 * Now global refinement is executed. Cells near the transducer location
880 * will be automatically refined according to the circle shaped boundary
881 * of the transducer lens:
888 * Lastly, we generate some more logging output. We stop the timer and
889 * query the number of CPU seconds elapsed since the beginning of the
894 * std::cout <<
"done (" << timer.cpu_time() <<
"s)" << std::endl;
896 * std::cout <<
" Number of active cells: " <<
triangulation.n_active_cells()
904 * <a name=
"step_29-codeUltrasoundProblemsetup_systemcode"></a>
905 * <h4><code>UltrasoundProblem::setup_system</code></h4>
909 * Initialization of the system
matrix, sparsity patterns and vectors are
910 * the same as in previous examples and therefore
do not need further
911 * comment. As in the previous function, we also output the
run time of what
916 *
void UltrasoundProblem<dim>::setup_system()
918 * std::cout <<
"Setting up system... ";
921 * dof_handler.distribute_dofs(fe);
925 * sparsity_pattern.copy_from(dsp);
927 * system_matrix.reinit(sparsity_pattern);
928 * system_rhs.reinit(dof_handler.n_dofs());
929 * solution.reinit(dof_handler.n_dofs());
932 * std::cout <<
"done (" << timer.cpu_time() <<
"s)" << std::endl;
934 * std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
942 * <a name=
"step_29-codeUltrasoundProblemassemble_systemcode"></a>
943 * <h4><code>UltrasoundProblem::assemble_system</code></h4>
947 * As before,
this function takes care of assembling the system
matrix and
948 * right hand side vector:
952 *
void UltrasoundProblem<dim>::assemble_system()
954 * std::cout <<
"Assembling system matrix... ";
960 * and store them in local variables, as they will be used frequently
961 * throughout
this function.
969 *
const double omega = prm.get_double(
"omega"), c = prm.get_double(
"c");
971 * prm.leave_subsection();
975 * As usual,
for computing integrals ordinary Gauss quadrature rule is
976 * used. Since our bilinear form involves boundary integrals on
977 * @f$\Gamma_2@f$, we also need a quadrature rule
for surface integration on
978 * the faces, which are @f$dim-1@f$ dimensional:
981 *
const QGauss<dim> quadrature_formula(fe.degree + 1);
982 *
const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
984 *
const unsigned int n_q_points = quadrature_formula.size(),
985 * n_face_q_points = face_quadrature_formula.size(),
986 * dofs_per_cell = fe.n_dofs_per_cell();
991 * part of the bilinear form that involves integration on @f$\Omega@f$, we
'll
992 * need the values and gradients of the shape functions, and of course the
993 * quadrature weights. For the terms involving the boundary integrals,
994 * only shape function values and the quadrature weights are necessary.
997 * FEValues<dim> fe_values(fe,
998 * quadrature_formula,
999 * update_values | update_gradients |
1000 * update_JxW_values);
1002 * FEFaceValues<dim> fe_face_values(fe,
1003 * face_quadrature_formula,
1004 * update_values | update_JxW_values);
1008 * As usual, the system matrix is assembled cell by cell, and we need a
1009 * matrix for storing the local cell contributions as well as an index
1010 * vector to transfer the cell contributions to the appropriate location
1011 * in the global system matrix after.
1014 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1015 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1017 * for (const auto &cell : dof_handler.active_cell_iterators())
1021 * On each cell, we first need to reset the local contribution matrix
1022 * and request the FEValues object to compute the shape functions for
1027 * fe_values.reinit(cell);
1029 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1031 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1035 * At this point, it is important to keep in mind that we are
1036 * dealing with a finite element system with two
1037 * components. Due to the way we constructed this FESystem,
1038 * namely as the Cartesian product of two scalar finite
1039 * element fields, each shape function has only a single
1040 * nonzero component (they are, in deal.II lingo, @ref
1041 * GlossPrimitive "primitive"). Hence, each shape function
1042 * can be viewed as one of the @f$\phi@f$'s or @f$\psi@f$
's from the
1043 * introduction, and similarly the corresponding degrees of
1044 * freedom can be attributed to either @f$\alpha@f$ or @f$\beta@f$.
1045 * As we iterate through all the degrees of freedom on the
1046 * current cell however, they do not come in any particular
1047 * order, and so we cannot decide right away whether the DoFs
1048 * with index @f$i@f$ and @f$j@f$ belong to the real or imaginary part
1049 * of our solution. On the other hand, if you look at the
1050 * form of the system matrix in the introduction, this
1051 * distinction is crucial since it will determine to which
1052 * block in the system matrix the contribution of the current
1053 * pair of DoFs will go and hence which quantity we need to
1054 * compute from the given two shape functions. Fortunately,
1055 * the FESystem object can provide us with this information,
1056 * namely it has a function
1057 * FESystem::system_to_component_index(), that for each local
1058 * DoF index returns a pair of integers of which the first
1059 * indicates to which component of the system the DoF
1060 * belongs. The second integer of the pair indicates which
1061 * index the DoF has in the scalar base finite element field,
1062 * but this information is not relevant here. If you want to
1063 * know more about this function and the underlying scheme
1064 * behind primitive vector valued elements, take a look at
1065 * @ref step_8 "step-8" or the @ref vector_valued topic, where these topics
1066 * are explained in depth.
1069 * if (fe.system_to_component_index(i).first ==
1070 * fe.system_to_component_index(j).first)
1074 * If both DoFs @f$i@f$ and @f$j@f$ belong to same component,
1075 * i.e. their shape functions are both @f$\phi@f$'s or both
1076 * @f$\psi@f$
's, the contribution will end up in one of the
1077 * diagonal blocks in our system matrix, and since the
1078 * corresponding entries are computed by the same formula,
1079 * we do not bother if they actually are @f$\phi@f$ or @f$\psi@f$
1080 * shape functions. We can simply compute the entry by
1081 * iterating over all quadrature points and adding up
1082 * their contributions, where values and gradients of the
1083 * shape functions are supplied by our FEValues object.
1089 * for (unsigned int q_point = 0; q_point < n_q_points;
1091 * cell_matrix(i, j) +=
1092 * (((fe_values.shape_value(i, q_point) *
1093 * fe_values.shape_value(j, q_point)) *
1094 * (-omega * omega) +
1095 * (fe_values.shape_grad(i, q_point) *
1096 * fe_values.shape_grad(j, q_point)) *
1098 * fe_values.JxW(q_point));
1102 * You might think that we would have to specify which
1103 * component of the shape function we'd like to evaluate
1106 * are primitive, they have only one
nonzero component,
1107 * and the
FEValues class is smart enough to figure out
1108 * that we are definitely interested in
this one
nonzero
1119 * We also have to add contributions due to boundary terms. To
this
1120 *
end, we
loop over all faces of the current cell and see
if first it
1121 * is at the boundary, and
second has the correct boundary indicator
1122 * associated with @f$\Gamma_2@f$, the part of the boundary where we have
1123 * absorbing boundary conditions:
1126 *
for (
const auto face_no : cell->face_indices())
1127 * if (cell->face(face_no)->at_boundary() &&
1132 * These faces will certainly contribute to the off-
diagonal
1134 * object to provide us with the shape function
values on this
1138 * fe_face_values.
reinit(cell, face_no);
1143 * Next, we
loop through all DoFs of the current cell to find
1144 * pairs that belong to different components and both have
1145 * support on the current face_no:
1148 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1149 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
1150 *
if ((fe.system_to_component_index(i).first !=
1151 * fe.system_to_component_index(j).first) &&
1152 * fe.has_support_on_face(i, face_no) &&
1153 * fe.has_support_on_face(j, face_no))
1157 * face is not strictly necessary:
if we don
't check for
1158 * it we would simply add up terms to the local cell
1159 * matrix that happen to be zero because at least one of
1160 * the shape functions happens to be zero. However, we can
1161 * save that work by adding the checks above.
1165 * In either case, these DoFs will contribute to the
1166 * boundary integrals in the off-diagonal blocks of the
1167 * system matrix. To compute the integral, we loop over
1168 * all the quadrature points on the face and sum up the
1169 * contribution weighted with the quadrature weights that
1170 * the face quadrature rule provides. In contrast to the
1171 * entries on the diagonal blocks, here it does matter
1172 * which one of the shape functions is a @f$\psi@f$ and which
1173 * one is a @f$\phi@f$, since that will determine the sign of
1174 * the entry. We account for this by a simple conditional
1175 * statement that determines the correct sign. Since we
1176 * already checked that DoF @f$i@f$ and @f$j@f$ belong to
1177 * different components, it suffices here to test for one
1178 * of them to which component it belongs.
1181 * for (unsigned int q_point = 0; q_point < n_face_q_points;
1183 * cell_matrix(i, j) +=
1184 * ((fe.system_to_component_index(i).first == 0) ? -1 :
1186 * fe_face_values.shape_value(i, q_point) *
1187 * fe_face_values.shape_value(j, q_point) * c * omega *
1188 * fe_face_values.JxW(q_point);
1193 * Now we are done with this cell and have to transfer its
1194 * contributions from the local to the global system matrix. To this
1195 * end, we first get a list of the global indices of the this cells
1199 * cell->get_dof_indices(local_dof_indices);
1204 * ...and then add the entries to the system matrix one by one:
1207 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1208 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1209 * system_matrix.add(local_dof_indices[i],
1210 * local_dof_indices[j],
1211 * cell_matrix(i, j));
1217 * The only thing left are the Dirichlet boundary values on @f$\Gamma_1@f$,
1218 * which is characterized by the boundary indicator 1. The Dirichlet
1219 * values are provided by the <code>DirichletBoundaryValues</code> class
1223 * std::map<types::global_dof_index, double> boundary_values;
1224 * VectorTools::interpolate_boundary_values(dof_handler,
1226 * DirichletBoundaryValues<dim>(),
1229 * MatrixTools::apply_boundary_values(boundary_values,
1235 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
1243 * <a name="step_29-codeUltrasoundProblemsolvecode"></a>
1244 * <h4><code>UltrasoundProblem::solve</code></h4>
1248 * As already mentioned in the introduction, the system matrix is neither
1249 * symmetric nor definite, and so it is not quite obvious how to come up
1250 * with an iterative solver and a preconditioner that do a good job on this
1251 * matrix. (For more on this topic, see also the
1252 * <a href="#extensions">Possibilities for extensions</a> section below.)
1253 * We chose instead to go a different way and solve the linear
1254 * system with the sparse LU decomposition provided by UMFPACK. This is
1255 * often a good first choice for 2d problems and works reasonably well even
1256 * for moderately large numbers of DoFs. The deal.II interface to UMFPACK
1257 * is implemented in the SparseDirectUMFPACK class, which is very easy to
1258 * use and allows us to solve our linear system with just 3 lines of code.
1262 * Note again that for compiling this example program, you need to have the
1263 * deal.II library built with UMFPACK support.
1266 * template <int dim>
1267 * void UltrasoundProblem<dim>::solve()
1269 * std::cout << "Solving linear system... ";
1274 * The code to solve the linear system is short: First, we allocate an
1275 * object of the right type. The following call to
1276 * SparseDirectUMFPACK::solve() takes as argument the matrix to decompose,
1277 * and a vector that upon input equals the right hand side of the linear
1278 * system to be solved, and upon output contains the solution of the linear
1279 * system. To satisfy this input/output requirement, we first assign the
1280 * right hand side vector to the `solution` variable.
1283 * SparseDirectUMFPACK A_direct;
1285 * solution = system_rhs;
1286 * A_direct.solve(system_matrix, solution);
1289 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
1297 * <a name="step_29-codeUltrasoundProblemoutput_resultscode"></a>
1298 * <h4><code>UltrasoundProblem::output_results</code></h4>
1302 * Here we output our solution @f$v@f$ and @f$w@f$ as well as the derived quantity
1303 * @f$|u|@f$ in the format specified in the parameter file. Most of the work for
1304 * deriving @f$|u|@f$ from @f$v@f$ and @f$w@f$ was already done in the implementation of
1305 * the <code>ComputeIntensity</code> class, so that the output routine is
1306 * rather straightforward and very similar to what is done in the previous
1310 * template <int dim>
1311 * void UltrasoundProblem<dim>::output_results() const
1313 * std::cout << "Generating output... ";
1318 * Define objects of our <code>ComputeIntensity</code> class and a DataOut
1322 * ComputeIntensity<dim> intensities;
1323 * DataOut<dim> data_out;
1325 * data_out.attach_dof_handler(dof_handler);
1329 * Next we query the output-related parameters from the ParameterHandler.
1330 * The DataOut::parse_parameters call acts as a counterpart to the
1331 * DataOutInterface<1>::declare_parameters call in
1332 * <code>ParameterReader::declare_parameters</code>. It collects all the
1333 * output format related parameters from the ParameterHandler and sets the
1334 * corresponding properties of the DataOut object accordingly.
1337 * prm.enter_subsection("Output parameters");
1339 * const std::string output_filename = prm.get("Output filename");
1340 * data_out.parse_parameters(prm);
1342 * prm.leave_subsection();
1346 * Now we put together the filename from the base name provided by the
1347 * ParameterHandler and the suffix which is provided by the DataOut class
1348 * (the default suffix is set to the right type that matches the one set
1349 * in the .prm file through parse_parameters()):
1352 * const std::string filename = output_filename + data_out.default_suffix();
1354 * std::ofstream output(filename);
1358 * The solution vectors @f$v@f$ and @f$w@f$ are added to the DataOut object in the
1362 * std::vector<std::string> solution_names;
1363 * solution_names.emplace_back("Re_u");
1364 * solution_names.emplace_back("Im_u");
1366 * data_out.add_data_vector(solution, solution_names);
1370 * For the intensity, we just call <code>add_data_vector</code> again, but
1371 * this with our <code>ComputeIntensity</code> object as the second
1372 * argument, which effectively adds @f$|u|@f$ to the output data:
1375 * data_out.add_data_vector(solution, intensities);
1379 * The last steps are as before. Note that the actual output format is now
1380 * determined by what is stated in the input file, i.e. one can change the
1381 * output format without having to re-compile this program:
1384 * data_out.build_patches();
1385 * data_out.write(output);
1388 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
1396 * <a name="step_29-codeUltrasoundProblemruncode"></a>
1397 * <h4><code>UltrasoundProblem::run</code></h4>
1401 * Here we simply execute our functions one after the other:
1404 * template <int dim>
1405 * void UltrasoundProblem<dim>::run()
1409 * assemble_system();
1413 * } // namespace Step29
1419 * <a name="step_29-Thecodemaincodefunction"></a>
1420 * <h4>The <code>main</code> function</h4>
1424 * Finally the <code>main</code> function of the program. It has the same
1425 * structure as in almost all of the other tutorial programs. The only
1426 * exception is that we define ParameterHandler and
1427 * <code>ParameterReader</code> objects, and let the latter read in the
1428 * parameter values from a textfile called <code>@ref step_29 "step-29".prm</code>. The
1429 * values so read are then handed over to an instance of the UltrasoundProblem
1437 * using namespace dealii;
1438 * using namespace Step29;
1440 * ParameterHandler prm;
1441 * ParameterReader param(prm);
1442 * param.read_parameters("step-29.prm");
1444 * UltrasoundProblem<2> ultrasound_problem(prm);
1445 * ultrasound_problem.run();
1447 * catch (std::exception &exc)
1449 * std::cerr << std::endl
1451 * << "----------------------------------------------------"
1453 * std::cerr << "Exception on processing: " << std::endl
1454 * << exc.what() << std::endl
1455 * << "Aborting!" << std::endl
1456 * << "----------------------------------------------------"
1462 * std::cerr << std::endl
1464 * << "----------------------------------------------------"
1466 * std::cerr << "Unknown exception!" << std::endl
1467 * << "Aborting!" << std::endl
1468 * << "----------------------------------------------------"
1475<a name="step_29-Results"></a><h1>Results</h1>
1478The current program reads its run-time parameters from an input file
1479called <code>step-29.prm</code> that looks like this:
1481subsection Mesh & geometry parameters
1482 # Distance of the focal point of the lens to the x-axis
1483 set Focal distance = 0.3
1485 # Number of global mesh refinement steps applied to initial coarse grid
1486 set Number of refinements = 5
1490subsection Physical constants
1499subsection Output parameters
1500 # Name of the output file (without extension)
1501 set Output file = solution
1503 # A name for the output format to be used
1504 set Output format = vtu
1508As can be seen, we set
1509@f$d=0.3@f$, which amounts to a focus of the transducer lens
1510at @f$x=0.5@f$, @f$y=0.3@f$. The coarse mesh is refined 5 times,
1511resulting in 160x160 cells, and the output is written in vtu
1512format. The parameter reader understands many more parameters
1513pertaining in particular to the generation of output, but we
1514need none of these parameters here and therefore stick with
1515their default values.
1517Here's the console output of the program in debug mode:
1521[ 33%] Building CXX
object CMakeFiles/step-29.dir/step-29.cc.o
1522[ 66%] Linking CXX executable @ref step_29
"step-29"
1523[ 66%] Built target @ref step_29
"step-29"
1524[100%] Run @ref step_29
"step-29" with Debug configuration
1525Generating grid... done (0.675234s)
1526 Number of active cells: 25600
1527Setting up system... done (0.375661s)
1528 Number of degrees of freedom: 51842
1529Assembling system
matrix... done (1.96864s)
1530Solving linear system... done (0.861862s)
1531Generating output... done (1.66913s)
1532[100%] Built target
run
1535(Of course, execution times will differ
if you run the program
1536locally.) The fact that most of the time is spent on assembling
1537the system
matrix and generating output is due to the many assertions
1538that need to be checked in debug mode. In release mode these parts
1539of the program
run much faster whereas solving the linear system is
1540hardly sped up at all:
1544[ 33%] Building CXX
object CMakeFiles/step-29.dir/step-29.cc.o
1545[ 66%] Linking CXX executable @ref step_29
"step-29"
1546[ 66%] Built target @ref step_29
"step-29"
1547[100%] Run @ref step_29
"step-29" with Release configuration
1548Generating grid... done (0.005207s)
1549 Number of active cells: 25600
1550Setting up system... done (0.018111s)
1551 Number of degrees of freedom: 51842
1552Assembling system
matrix... done (0.021501s)
1553Solving linear system... done (0.590303s)
1554Generating output... done (0.090787s)
1555[100%] Built target
run
1558The graphical output of the program looks as follows:
1561<table align=
"center" class=
"doxtable">
1564 <img src=
"https://www.dealii.org/images/steps/developer/step-29.v.png" alt=
"v = Re(u)">
1567 <img src=
"https://www.dealii.org/images/steps/developer/step-29.w.png" alt=
"w = Im(u)">
1572 <img src=
"https://www.dealii.org/images/steps/developer/step-29.intensity.png" alt=
"|u|">
1577The
first two pictures show the real and imaginary parts of
1578@f$u@f$, whereas the last shows the intensity @f$|u|@f$. One can clearly
1579see that the intensity is focused around the focal
point of the
1580lens (0.5, 0.3), and that the focus
1581is rather sharp in @f$x@f$-direction but more blurred in @f$y@f$-direction, which is a
1582consequence of the geometry of the focusing lens, its finite aperture,
1583and the wave nature of the problem.
1585Because colorful graphics are
always fun, and to stress the focusing
1586effects some more, here is another set of images highlighting how well
1587the intensity is actually focused in @f$x@f$-direction:
1589<table align=
"center" class=
"doxtable">
1592 <img src=
"https://www.dealii.org/images/steps/developer/step-29.surface.png" alt=
"|u|">
1595 <img src=
"https://www.dealii.org/images/steps/developer/step-29.contours.png" alt=
"|u|">
1601As a
final note, the structure of the program makes it easy to
1602determine which parts of the program
scale nicely as the mesh is
1603refined and which parts don
't. Here are the run times for 5, 6, and 7
1608[ 66%] Built target @ref step_29 "step-29"
1609[100%] Run @ref step_29 "step-29" with Release configuration
1610Generating grid... done (0.008767s)
1611 Number of active cells: 25600
1612Setting up system... done (0.007043s)
1613 Number of degrees of freedom: 51842
1614Assembling system matrix... done (0.020942s)
1615Solving linear system... done (0.585853s)
1616Generating output... done (0.092681s)
1617[100%] Built target run
1620[ 66%] Built target @ref step_29 "step-29"
1621[100%] Run @ref step_29 "step-29" with Release configuration
1622Generating grid... done (0.027986s)
1623 Number of active cells: 102400
1624Setting up system... done (0.070022s)
1625 Number of degrees of freedom: 206082
1626Assembling system matrix... done (0.075252s)
1627Solving linear system... done (4.13989s)
1628Generating output... done (0.369556s)
1629[100%] Built target run
1632[ 66%] Built target @ref step_29 "step-29"
1633[100%] Run @ref step_29 "step-29" with Release configuration
1634Generating grid... done (0.186439s)
1635 Number of active cells: 409600
1636Setting up system... done (0.272828s)
1637 Number of degrees of freedom: 821762
1638Assembling system matrix... done (0.316252s)
1639Solving linear system... done (33.1424s)
1640Generating output... done (1.36328s)
1641[100%] Built target run
1644Each time we refine the mesh once more, so the number of cells and degrees
1645of freedom roughly quadruples from each step to the next. As can be seen,
1646generating the grid, setting up degrees of freedom, assembling the
1647linear system, and generating output scale pretty closely to linear,
1648whereas solving the linear system is an operation that requires 8
1649times more time each time the number of degrees of freedom is
1650increased by a factor of 4, i.e. it is @f${\cal O}(N^{3/2})@f$. This can
1651be explained by the fact that (using optimal ordering) the
1652bandwidth of a finite element matrix is @f$B={\cal O}(N^{(dim-1)/dim})@f$,
1653and the effort to solve a banded linear system using LU decomposition
1654is @f${\cal O}(BN)@f$. This also explains why the program does run in 3d
1655as well (after changing the dimension on the
1656<code>UltrasoundProblem</code> object), but scales very badly and
1657takes extraordinary patience before it finishes solving the linear
1658system on a mesh with appreciable resolution, even though all the
1659other parts of the program scale very nicely.
1663<a name="step-29-extensions"></a>
1664<a name="step_29-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1667An obvious possible extension for this program is to run it in 3d
1668— after all, the world around us is three-dimensional, and
1669ultrasound beams propagate in three-dimensional media. You can try
1670this by simply changing the template parameter of the principal class
1671in <code>main()</code> and running it. This won't get you very far,
1672though: certainly not
if you
do 5 global refinement steps as set in
1673the parameter file. You
'll simply run out of memory as both the mesh
1674(with its @f$(2^5)^3 \cdot 5^3=2^{15}\cdot 125 \approx 4\cdot 10^6@f$ cells)
1675and in particular the sparse direct solver take too much memory. You
1676can solve with 3 global refinement steps, however, if you have a bit
1677of time: in early 2011, the direct solve takes about half an
1678hour. What you'll notice, however, is that the solution is completely
1679wrong: the mesh size is simply not small enough to resolve the
1680solution
's waves accurately, and you can see this in plots of the
1681solution. Consequently, this is one of the cases where adaptivity is
1682indispensable if you don't just want to
throw a bigger (presumably
1686<a name=
"step_29-PlainProg"></a>
1687<h1> The plain program</h1>
1688@include
"step-29.cc"
static void declare_parameters(ParameterHandler &prm)
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={})
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
#define AssertDimension(dim1, dim2)
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())
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_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
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)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation