1 /* ---------------------------------------------------------------------
2 *
3 * Copyright (C) 2007 - 2020 by the deal.II authors
4 *
5 * This file is part of the deal.II library.
6 *
7 * The deal.II library is free software; you can use it, redistribute
8 * it, and/or modify it under the terms of the GNU Lesser General
9 * Public License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 * The full text of the license can be found in the file LICENSE.md at
12 * the top level directory of deal.II.
13 *
14 * ---------------------------------------------------------------------
15
16 *
17 * Author: Moritz Allmaras, Texas A&M University, 2007
18 */
19
20
21 // @sect3{Include files}
22
23 // The following header files have all been discussed before:
24
25 #include <deal.II/base/quadrature_lib.h>
26 #include <deal.II/base/function.h>
27 #include <deal.II/base/logstream.h>
28
29 #include <deal.II/lac/vector.h>
30 #include <deal.II/lac/full_matrix.h>
31 #include <deal.II/lac/sparse_matrix.h>
32 #include <deal.II/lac/dynamic_sparsity_pattern.h>
33
34 #include <deal.II/grid/tria.h>
35 #include <deal.II/grid/grid_generator.h>
36 #include <deal.II/grid/tria_accessor.h>
37 #include <deal.II/grid/tria_iterator.h>
38 #include <deal.II/grid/manifold_lib.h>
39
40 #include <deal.II/dofs/dof_handler.h>
41 #include <deal.II/dofs/dof_accessor.h>
42 #include <deal.II/dofs/dof_tools.h>
43
44 #include <deal.II/fe/fe_q.h>
45 #include <deal.II/fe/fe_values.h>
46
47 #include <deal.II/numerics/matrix_tools.h>
48 #include <deal.II/numerics/data_out.h>
49 #include <deal.II/numerics/vector_tools.h>
50
51 #include <iostream>
52 #include <fstream>
53
54 // This header file contains the necessary declarations for the
55 // ParameterHandler class that we will use to read our parameters from a
56 // configuration file:
57 #include <deal.II/base/parameter_handler.h>
58
59 // For solving the linear system, we'll use the sparse LU-decomposition
60 // provided by UMFPACK (see the SparseDirectUMFPACK class), for which the
61 // following header file is needed. Note that in order to compile this
62 // tutorial program, the deal.II-library needs to be built with UMFPACK
63 // support, which is enabled by default:
64 #include <deal.II/lac/sparse_direct.h>
65
66 // The FESystem class allows us to stack several FE-objects to one compound,
67 // vector-valued finite element field. The necessary declarations for this
68 // class are provided in this header file:
69 #include <deal.II/fe/fe_system.h>
70
71 // Finally, include the header file that declares the Timer class that we will
72 // use to determine how much time each of the operations of our program takes:
73 #include <deal.II/base/timer.h>
74
75 // As the last step at the beginning of this program, we put everything that
76 // is in this program into its namespace and, within it, make everything that
77 // is in the deal.II namespace globally available, without the need to prefix
78 // everything with <code>dealii</code><code>::</code>:
79 namespace Step29
80 {
81 using namespace dealii;
82
83
84 // @sect3{The <code>DirichletBoundaryValues</code> class}
85
86 // First we define a class for the function representing the Dirichlet
87 // boundary values. This has been done many times before and therefore does
88 // not need much explanation.
89 //
90 // Since there are two values $v$ and $w$ that need to be prescribed at the
91 // boundary, we have to tell the base class that this is a vector-valued
92 // function with two components, and the <code>vector_value</code> function
93 // and its cousin <code>vector_value_list</code> must return vectors with
94 // two entries. In our case the function is very simple, it just returns 1
95 // for the real part $v$ and 0 for the imaginary part $w$ regardless of the
96 // point where it is evaluated.
97 template <int dim>
98 class DirichletBoundaryValues : public Function<dim>
99 {
100 public:
DirichletBoundaryValues()101 DirichletBoundaryValues()
102 : Function<dim>(2)
103 {}
104
vector_value(const Point<dim> &,Vector<double> & values) const105 virtual void vector_value(const Point<dim> & /*p*/,
106 Vector<double> &values) const override
107 {
108 Assert(values.size() == 2, ExcDimensionMismatch(values.size(), 2));
109
110 values(0) = 1;
111 values(1) = 0;
112 }
113
114 virtual void
vector_value_list(const std::vector<Point<dim>> & points,std::vector<Vector<double>> & value_list) const115 vector_value_list(const std::vector<Point<dim>> &points,
116 std::vector<Vector<double>> & value_list) const override
117 {
118 Assert(value_list.size() == points.size(),
119 ExcDimensionMismatch(value_list.size(), points.size()));
120
121 for (unsigned int p = 0; p < points.size(); ++p)
122 DirichletBoundaryValues<dim>::vector_value(points[p], value_list[p]);
123 }
124 };
125
126 // @sect3{The <code>ParameterReader</code> class}
127
128 // The next class is responsible for preparing the ParameterHandler object
129 // and reading parameters from an input file. It includes a function
130 // <code>declare_parameters</code> that declares all the necessary
131 // parameters and a <code>read_parameters</code> function that is called
132 // from outside to initiate the parameter reading process.
133 class ParameterReader : public Subscriptor
134 {
135 public:
136 ParameterReader(ParameterHandler &);
137 void read_parameters(const std::string &);
138
139 private:
140 void declare_parameters();
141 ParameterHandler &prm;
142 };
143
144 // The constructor stores a reference to the ParameterHandler object that is
145 // passed to it:
ParameterReader(ParameterHandler & paramhandler)146 ParameterReader::ParameterReader(ParameterHandler ¶mhandler)
147 : prm(paramhandler)
148 {}
149
150 // @sect4{<code>ParameterReader::declare_parameters</code>}
151
152 // The <code>declare_parameters</code> function declares all the parameters
153 // that our ParameterHandler object will be able to read from input files,
154 // along with their types, range conditions and the subsections they appear
155 // in. We will wrap all the entries that go into a section in a pair of
156 // braces to force the editor to indent them by one level, making it simpler
157 // to read which entries together form a section:
declare_parameters()158 void ParameterReader::declare_parameters()
159 {
160 // Parameters for mesh and geometry include the number of global
161 // refinement steps that are applied to the initial coarse mesh and the
162 // focal distance $d$ of the transducer lens. For the number of refinement
163 // steps, we allow integer values in the range $[0,\infty)$, where the
164 // omitted second argument to the Patterns::Integer object denotes the
165 // half-open interval. For the focal distance any number greater than
166 // zero is accepted:
167 prm.enter_subsection("Mesh & geometry parameters");
168 {
169 prm.declare_entry("Number of refinements",
170 "6",
171 Patterns::Integer(0),
172 "Number of global mesh refinement steps "
173 "applied to initial coarse grid");
174
175 prm.declare_entry("Focal distance",
176 "0.3",
177 Patterns::Double(0),
178 "Distance of the focal point of the lens "
179 "to the x-axis");
180 }
181 prm.leave_subsection();
182
183 // The next subsection is devoted to the physical parameters appearing in
184 // the equation, which are the frequency $\omega$ and wave speed
185 // $c$. Again, both need to lie in the half-open interval $[0,\infty)$
186 // represented by calling the Patterns::Double class with only the left
187 // end-point as argument:
188 prm.enter_subsection("Physical constants");
189 {
190 prm.declare_entry("c", "1.5e5", Patterns::Double(0), "Wave speed");
191
192 prm.declare_entry("omega", "5.0e7", Patterns::Double(0), "Frequency");
193 }
194 prm.leave_subsection();
195
196
197 // Last but not least we would like to be able to change some properties
198 // of the output, like filename and format, through entries in the
199 // configuration file, which is the purpose of the last subsection:
200 prm.enter_subsection("Output parameters");
201 {
202 prm.declare_entry("Output filename",
203 "solution",
204 Patterns::Anything(),
205 "Name of the output file (without extension)");
206
207 // Since different output formats may require different parameters for
208 // generating output (like for example, postscript output needs
209 // viewpoint angles, line widths, colors etc), it would be cumbersome if
210 // we had to declare all these parameters by hand for every possible
211 // output format supported in the library. Instead, each output format
212 // has a <code>FormatFlags::declare_parameters</code> function, which
213 // declares all the parameters specific to that format in an own
214 // subsection. The following call of
215 // DataOutInterface<1>::declare_parameters executes
216 // <code>declare_parameters</code> for all available output formats, so
217 // that for each format an own subsection will be created with
218 // parameters declared for that particular output format. (The actual
219 // value of the template parameter in the call, <code>@<1@></code>
220 // above, does not matter here: the function does the same work
221 // independent of the dimension, but happens to be in a
222 // template-parameter-dependent class.) To find out what parameters
223 // there are for which output format, you can either consult the
224 // documentation of the DataOutBase class, or simply run this program
225 // without a parameter file present. It will then create a file with all
226 // declared parameters set to their default values, which can
227 // conveniently serve as a starting point for setting the parameters to
228 // the values you desire.
229 DataOutInterface<1>::declare_parameters(prm);
230 }
231 prm.leave_subsection();
232 }
233
234 // @sect4{<code>ParameterReader::read_parameters</code>}
235
236 // This is the main function in the ParameterReader class. It gets called
237 // from outside, first declares all the parameters, and then reads them from
238 // the input file whose filename is provided by the caller. After the call
239 // to this function is complete, the <code>prm</code> object can be used to
240 // retrieve the values of the parameters read in from the file:
read_parameters(const std::string & parameter_file)241 void ParameterReader::read_parameters(const std::string ¶meter_file)
242 {
243 declare_parameters();
244
245 prm.parse_input(parameter_file);
246 }
247
248
249
250 // @sect3{The <code>ComputeIntensity</code> class}
251
252 // As mentioned in the introduction, the quantity that we are really after
253 // is the spatial distribution of the intensity of the ultrasound wave,
254 // which corresponds to $|u|=\sqrt{v^2+w^2}$. Now we could just be content
255 // with having $v$ and $w$ in our output, and use a suitable visualization
256 // or postprocessing tool to derive $|u|$ from the solution we
257 // computed. However, there is also a way to output data derived from the
258 // solution in deal.II, and we are going to make use of this mechanism here.
259
260 // So far we have always used the DataOut::add_data_vector function to add
261 // vectors containing output data to a DataOut object. There is a special
262 // version of this function that in addition to the data vector has an
263 // additional argument of type DataPostprocessor. What happens when this
264 // function is used for output is that at each point where output data is to
265 // be generated, the DataPostprocessor::evaluate_scalar_field() or
266 // DataPostprocessor::evaluate_vector_field() function of the
267 // specified DataPostprocessor object is invoked to compute the output
268 // quantities from the values, the gradients and the second derivatives of
269 // the finite element function represented by the data vector (in the case
270 // of face related data, normal vectors are available as well). Hence, this
271 // allows us to output any quantity that can locally be derived from the
272 // values of the solution and its derivatives. Of course, the ultrasound
273 // intensity $|u|$ is such a quantity and its computation doesn't even
274 // involve any derivatives of $v$ or $w$.
275
276 // In practice, the DataPostprocessor class only provides an interface to
277 // this functionality, and we need to derive our own class from it in order
278 // to implement the functions specified by the interface. In the most
279 // general case one has to implement several member functions but if the
280 // output quantity is a single scalar then some of this boilerplate code can
281 // be handled by a more specialized class, DataPostprocessorScalar and we
282 // can derive from that one instead. This is what the
283 // <code>ComputeIntensity</code> class does:
284 template <int dim>
285 class ComputeIntensity : public DataPostprocessorScalar<dim>
286 {
287 public:
288 ComputeIntensity();
289
290 virtual void evaluate_vector_field(
291 const DataPostprocessorInputs::Vector<dim> &inputs,
292 std::vector<Vector<double>> &computed_quantities) const override;
293 };
294
295 // In the constructor, we need to call the constructor of the base class
296 // with two arguments. The first denotes the name by which the single scalar
297 // quantity computed by this class should be represented in output files. In
298 // our case, the postprocessor has $|u|$ as output, so we use "Intensity".
299 //
300 // The second argument is a set of flags that indicate which data is needed
301 // by the postprocessor in order to compute the output quantities. This can
302 // be any subset of update_values, update_gradients and update_hessians
303 // (and, in the case of face data, also update_normal_vectors), which are
304 // documented in UpdateFlags. Of course, computation of the derivatives
305 // requires additional resources, so only the flags for data that are really
306 // needed should be given here, just as we do when we use FEValues objects.
307 // In our case, only the function values of $v$ and $w$ are needed to
308 // compute $|u|$, so we're good with the update_values flag.
309 template <int dim>
ComputeIntensity()310 ComputeIntensity<dim>::ComputeIntensity()
311 : DataPostprocessorScalar<dim>("Intensity", update_values)
312 {}
313
314
315 // The actual postprocessing happens in the following function. Its input is
316 // an object that stores values of the function (which is here vector-valued)
317 // representing the data vector given to DataOut::add_data_vector, evaluated
318 // at all evaluation points where we generate output, and some tensor objects
319 // representing derivatives (that we don't use here since $|u|$ is computed
320 // from just $v$ and $w$). The derived quantities are returned in the
321 // <code>computed_quantities</code> vector. Remember that this function may
322 // only use data for which the respective update flag is specified by
323 // <code>get_needed_update_flags</code>. For example, we may not use the
324 // derivatives here, since our implementation of
325 // <code>get_needed_update_flags</code> requests that only function values
326 // are provided.
327 template <int dim>
evaluate_vector_field(const DataPostprocessorInputs::Vector<dim> & inputs,std::vector<Vector<double>> & computed_quantities) const328 void ComputeIntensity<dim>::evaluate_vector_field(
329 const DataPostprocessorInputs::Vector<dim> &inputs,
330 std::vector<Vector<double>> & computed_quantities) const
331 {
332 Assert(computed_quantities.size() == inputs.solution_values.size(),
333 ExcDimensionMismatch(computed_quantities.size(),
334 inputs.solution_values.size()));
335
336 // The computation itself is straightforward: We iterate over each
337 // entry in the output vector and compute $|u|$ from the
338 // corresponding values of $v$ and $w$. We do this by creating a
339 // complex number $u$ and then calling `std::abs()` on the
340 // result. (One may be tempted to call `std::norm()`, but in a
341 // historical quirk, the C++ committee decided that `std::norm()`
342 // should return the <i>square</i> of the absolute value --
343 // thereby not satisfying the properties mathematicians require of
344 // something called a "norm".)
345 for (unsigned int i = 0; i < computed_quantities.size(); i++)
346 {
347 Assert(computed_quantities[i].size() == 1,
348 ExcDimensionMismatch(computed_quantities[i].size(), 1));
349 Assert(inputs.solution_values[i].size() == 2,
350 ExcDimensionMismatch(inputs.solution_values[i].size(), 2));
351
352 const std::complex<double> u(inputs.solution_values[i](0),
353 inputs.solution_values[i](1));
354
355 computed_quantities[i](0) = std::abs(u);
356 }
357 }
358
359
360 // @sect3{The <code>UltrasoundProblem</code> class}
361
362 // Finally here is the main class of this program. It's member functions
363 // are very similar to the previous examples, in particular step-4, and the
364 // list of member variables does not contain any major surprises either.
365 // The ParameterHandler object that is passed to the constructor is stored
366 // as a reference to allow easy access to the parameters from all functions
367 // of the class. Since we are working with vector valued finite elements,
368 // the FE object we are using is of type FESystem.
369 template <int dim>
370 class UltrasoundProblem
371 {
372 public:
373 UltrasoundProblem(ParameterHandler &);
374 void run();
375
376 private:
377 void make_grid();
378 void setup_system();
379 void assemble_system();
380 void solve();
381 void output_results() const;
382
383 ParameterHandler &prm;
384
385 Triangulation<dim> triangulation;
386 DoFHandler<dim> dof_handler;
387 FESystem<dim> fe;
388
389 SparsityPattern sparsity_pattern;
390 SparseMatrix<double> system_matrix;
391 Vector<double> solution, system_rhs;
392 };
393
394
395
396 // The constructor takes the ParameterHandler object and stores it in a
397 // reference. It also initializes the DoF-Handler and the finite element
398 // system, which consists of two copies of the scalar Q1 field, one for $v$
399 // and one for $w$:
400 template <int dim>
UltrasoundProblem(ParameterHandler & param)401 UltrasoundProblem<dim>::UltrasoundProblem(ParameterHandler ¶m)
402 : prm(param)
403 , dof_handler(triangulation)
404 , fe(FE_Q<dim>(1), 2)
405 {}
406
407 // @sect4{<code>UltrasoundProblem::make_grid</code>}
408
409 // Here we setup the grid for our domain. As mentioned in the exposition,
410 // the geometry is just a unit square (in 2d) with the part of the boundary
411 // that represents the transducer lens replaced by a sector of a circle.
412 template <int dim>
make_grid()413 void UltrasoundProblem<dim>::make_grid()
414 {
415 // First we generate some logging output and start a timer so we can
416 // compute execution time when this function is done:
417 std::cout << "Generating grid... ";
418 Timer timer;
419
420 // Then we query the values for the focal distance of the transducer lens
421 // and the number of mesh refinement steps from our ParameterHandler
422 // object:
423 prm.enter_subsection("Mesh & geometry parameters");
424
425 const double focal_distance = prm.get_double("Focal distance");
426 const unsigned int n_refinements = prm.get_integer("Number of refinements");
427
428 prm.leave_subsection();
429
430 // Next, two points are defined for position and focal point of the
431 // transducer lens, which is the center of the circle whose segment will
432 // form the transducer part of the boundary. Notice that this is the only
433 // point in the program where things are slightly different in 2D and 3D.
434 // Even though this tutorial only deals with the 2D case, the necessary
435 // additions to make this program functional in 3D are so minimal that we
436 // opt for including them:
437 const Point<dim> transducer =
438 (dim == 2) ? Point<dim>(0.5, 0.0) : Point<dim>(0.5, 0.5, 0.0);
439 const Point<dim> focal_point = (dim == 2) ?
440 Point<dim>(0.5, focal_distance) :
441 Point<dim>(0.5, 0.5, focal_distance);
442
443
444 // As initial coarse grid we take a simple unit square with 5 subdivisions
445 // in each direction. The number of subdivisions is chosen so that the
446 // line segment $[0.4,0.6]$ that we want to designate as the transducer
447 // boundary is spanned by a single face. Then we step through all cells to
448 // find the faces where the transducer is to be located, which in fact is
449 // just the single edge from 0.4 to 0.6 on the x-axis. This is where we
450 // want the refinements to be made according to a circle shaped boundary,
451 // so we mark this edge with a different manifold indicator. Since we will
452 // Dirichlet boundary conditions on the transducer, we also change its
453 // boundary indicator.
454 GridGenerator::subdivided_hyper_cube(triangulation, 5, 0, 1);
455
456 for (auto &cell : triangulation.cell_iterators())
457 for (const auto &face : cell->face_iterators())
458 if (face->at_boundary() &&
459 ((face->center() - transducer).norm_square() < 0.01))
460 {
461 face->set_boundary_id(1);
462 face->set_manifold_id(1);
463 }
464 // For the circle part of the transducer lens, a SphericalManifold object
465 // is used (which, of course, in 2D just represents a circle), with center
466 // computed as above.
467 triangulation.set_manifold(1, SphericalManifold<dim>(focal_point));
468
469 // Now global refinement is executed. Cells near the transducer location
470 // will be automatically refined according to the circle shaped boundary
471 // of the transducer lens:
472 triangulation.refine_global(n_refinements);
473
474 // Lastly, we generate some more logging output. We stop the timer and
475 // query the number of CPU seconds elapsed since the beginning of the
476 // function:
477 timer.stop();
478 std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
479
480 std::cout << " Number of active cells: " << triangulation.n_active_cells()
481 << std::endl;
482 }
483
484
485 // @sect4{<code>UltrasoundProblem::setup_system</code>}
486 //
487 // Initialization of the system matrix, sparsity patterns and vectors are
488 // the same as in previous examples and therefore do not need further
489 // comment. As in the previous function, we also output the run time of what
490 // we do here:
491 template <int dim>
setup_system()492 void UltrasoundProblem<dim>::setup_system()
493 {
494 std::cout << "Setting up system... ";
495 Timer timer;
496
497 dof_handler.distribute_dofs(fe);
498
499 DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
500 DoFTools::make_sparsity_pattern(dof_handler, dsp);
501 sparsity_pattern.copy_from(dsp);
502
503 system_matrix.reinit(sparsity_pattern);
504 system_rhs.reinit(dof_handler.n_dofs());
505 solution.reinit(dof_handler.n_dofs());
506
507 timer.stop();
508 std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
509
510 std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
511 << std::endl;
512 }
513
514
515 // @sect4{<code>UltrasoundProblem::assemble_system</code>}
516
517 // As before, this function takes care of assembling the system matrix and
518 // right hand side vector:
519 template <int dim>
assemble_system()520 void UltrasoundProblem<dim>::assemble_system()
521 {
522 std::cout << "Assembling system matrix... ";
523 Timer timer;
524
525 // First we query wavespeed and frequency from the ParameterHandler object
526 // and store them in local variables, as they will be used frequently
527 // throughout this function.
528
529 prm.enter_subsection("Physical constants");
530
531 const double omega = prm.get_double("omega"), c = prm.get_double("c");
532
533 prm.leave_subsection();
534
535 // As usual, for computing integrals ordinary Gauss quadrature rule is
536 // used. Since our bilinear form involves boundary integrals on
537 // $\Gamma_2$, we also need a quadrature rule for surface integration on
538 // the faces, which are $dim-1$ dimensional:
539 QGauss<dim> quadrature_formula(fe.degree + 1);
540 QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
541
542 const unsigned int n_q_points = quadrature_formula.size(),
543 n_face_q_points = face_quadrature_formula.size(),
544 dofs_per_cell = fe.n_dofs_per_cell();
545
546 // The FEValues objects will evaluate the shape functions for us. For the
547 // part of the bilinear form that involves integration on $\Omega$, we'll
548 // need the values and gradients of the shape functions, and of course the
549 // quadrature weights. For the terms involving the boundary integrals,
550 // only shape function values and the quadrature weights are necessary.
551 FEValues<dim> fe_values(fe,
552 quadrature_formula,
553 update_values | update_gradients |
554 update_JxW_values);
555
556 FEFaceValues<dim> fe_face_values(fe,
557 face_quadrature_formula,
558 update_values | update_JxW_values);
559
560 // As usual, the system matrix is assembled cell by cell, and we need a
561 // matrix for storing the local cell contributions as well as an index
562 // vector to transfer the cell contributions to the appropriate location
563 // in the global system matrix after.
564 FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
565 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
566
567 for (const auto &cell : dof_handler.active_cell_iterators())
568 {
569 // On each cell, we first need to reset the local contribution matrix
570 // and request the FEValues object to compute the shape functions for
571 // the current cell:
572 cell_matrix = 0;
573 fe_values.reinit(cell);
574
575 for (unsigned int i = 0; i < dofs_per_cell; ++i)
576 {
577 for (unsigned int j = 0; j < dofs_per_cell; ++j)
578 {
579 // At this point, it is important to keep in mind that we are
580 // dealing with a finite element system with two
581 // components. Due to the way we constructed this FESystem,
582 // namely as the Cartesian product of two scalar finite
583 // element fields, each shape function has only a single
584 // nonzero component (they are, in deal.II lingo, @ref
585 // GlossPrimitive "primitive"). Hence, each shape function
586 // can be viewed as one of the $\phi$'s or $\psi$'s from the
587 // introduction, and similarly the corresponding degrees of
588 // freedom can be attributed to either $\alpha$ or $\beta$.
589 // As we iterate through all the degrees of freedom on the
590 // current cell however, they do not come in any particular
591 // order, and so we cannot decide right away whether the DoFs
592 // with index $i$ and $j$ belong to the real or imaginary part
593 // of our solution. On the other hand, if you look at the
594 // form of the system matrix in the introduction, this
595 // distinction is crucial since it will determine to which
596 // block in the system matrix the contribution of the current
597 // pair of DoFs will go and hence which quantity we need to
598 // compute from the given two shape functions. Fortunately,
599 // the FESystem object can provide us with this information,
600 // namely it has a function
601 // FESystem::system_to_component_index, that for each local
602 // DoF index returns a pair of integers of which the first
603 // indicates to which component of the system the DoF
604 // belongs. The second integer of the pair indicates which
605 // index the DoF has in the scalar base finite element field,
606 // but this information is not relevant here. If you want to
607 // know more about this function and the underlying scheme
608 // behind primitive vector valued elements, take a look at
609 // step-8 or the @ref vector_valued module, where these topics
610 // are explained in depth.
611 if (fe.system_to_component_index(i).first ==
612 fe.system_to_component_index(j).first)
613 {
614 // If both DoFs $i$ and $j$ belong to same component,
615 // i.e. their shape functions are both $\phi$'s or both
616 // $\psi$'s, the contribution will end up in one of the
617 // diagonal blocks in our system matrix, and since the
618 // corresponding entries are computed by the same formula,
619 // we do not bother if they actually are $\phi$ or $\psi$
620 // shape functions. We can simply compute the entry by
621 // iterating over all quadrature points and adding up
622 // their contributions, where values and gradients of the
623 // shape functions are supplied by our FEValues object.
624
625 for (unsigned int q_point = 0; q_point < n_q_points;
626 ++q_point)
627 cell_matrix(i, j) +=
628 (((fe_values.shape_value(i, q_point) *
629 fe_values.shape_value(j, q_point)) *
630 (-omega * omega) +
631 (fe_values.shape_grad(i, q_point) *
632 fe_values.shape_grad(j, q_point)) *
633 c * c) *
634 fe_values.JxW(q_point));
635
636 // You might think that we would have to specify which
637 // component of the shape function we'd like to evaluate
638 // when requesting shape function values or gradients from
639 // the FEValues object. However, as the shape functions
640 // are primitive, they have only one nonzero component,
641 // and the FEValues class is smart enough to figure out
642 // that we are definitely interested in this one nonzero
643 // component.
644 }
645 }
646 }
647
648
649 // We also have to add contributions due to boundary terms. To this
650 // end, we loop over all faces of the current cell and see if first it
651 // is at the boundary, and second has the correct boundary indicator
652 // associated with $\Gamma_2$, the part of the boundary where we have
653 // absorbing boundary conditions:
654 for (const auto face_no : cell->face_indices())
655 if (cell->face(face_no)->at_boundary() &&
656 (cell->face(face_no)->boundary_id() == 0))
657 {
658 // These faces will certainly contribute to the off-diagonal
659 // blocks of the system matrix, so we ask the FEFaceValues
660 // object to provide us with the shape function values on this
661 // face:
662 fe_face_values.reinit(cell, face_no);
663
664
665 // Next, we loop through all DoFs of the current cell to find
666 // pairs that belong to different components and both have
667 // support on the current face_no:
668 for (unsigned int i = 0; i < dofs_per_cell; ++i)
669 for (unsigned int j = 0; j < dofs_per_cell; ++j)
670 if ((fe.system_to_component_index(i).first !=
671 fe.system_to_component_index(j).first) &&
672 fe.has_support_on_face(i, face_no) &&
673 fe.has_support_on_face(j, face_no))
674 // The check whether shape functions have support on a
675 // face is not strictly necessary: if we don't check for
676 // it we would simply add up terms to the local cell
677 // matrix that happen to be zero because at least one of
678 // the shape functions happens to be zero. However, we can
679 // save that work by adding the checks above.
680
681 // In either case, these DoFs will contribute to the
682 // boundary integrals in the off-diagonal blocks of the
683 // system matrix. To compute the integral, we loop over
684 // all the quadrature points on the face and sum up the
685 // contribution weighted with the quadrature weights that
686 // the face quadrature rule provides. In contrast to the
687 // entries on the diagonal blocks, here it does matter
688 // which one of the shape functions is a $\psi$ and which
689 // one is a $\phi$, since that will determine the sign of
690 // the entry. We account for this by a simple conditional
691 // statement that determines the correct sign. Since we
692 // already checked that DoF $i$ and $j$ belong to
693 // different components, it suffices here to test for one
694 // of them to which component it belongs.
695 for (unsigned int q_point = 0; q_point < n_face_q_points;
696 ++q_point)
697 cell_matrix(i, j) +=
698 ((fe.system_to_component_index(i).first == 0) ? -1 :
699 1) *
700 fe_face_values.shape_value(i, q_point) *
701 fe_face_values.shape_value(j, q_point) * c * omega *
702 fe_face_values.JxW(q_point);
703 }
704
705 // Now we are done with this cell and have to transfer its
706 // contributions from the local to the global system matrix. To this
707 // end, we first get a list of the global indices of the this cells
708 // DoFs...
709 cell->get_dof_indices(local_dof_indices);
710
711
712 // ...and then add the entries to the system matrix one by one:
713 for (unsigned int i = 0; i < dofs_per_cell; ++i)
714 for (unsigned int j = 0; j < dofs_per_cell; ++j)
715 system_matrix.add(local_dof_indices[i],
716 local_dof_indices[j],
717 cell_matrix(i, j));
718 }
719
720
721 // The only thing left are the Dirichlet boundary values on $\Gamma_1$,
722 // which is characterized by the boundary indicator 1. The Dirichlet
723 // values are provided by the <code>DirichletBoundaryValues</code> class
724 // we defined above:
725 std::map<types::global_dof_index, double> boundary_values;
726 VectorTools::interpolate_boundary_values(dof_handler,
727 1,
728 DirichletBoundaryValues<dim>(),
729 boundary_values);
730
731 MatrixTools::apply_boundary_values(boundary_values,
732 system_matrix,
733 solution,
734 system_rhs);
735
736 timer.stop();
737 std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
738 }
739
740
741
742 // @sect4{<code>UltrasoundProblem::solve</code>}
743
744 // As already mentioned in the introduction, the system matrix is neither
745 // symmetric nor definite, and so it is not quite obvious how to come up
746 // with an iterative solver and a preconditioner that do a good job on this
747 // matrix. We chose instead to go a different way and solve the linear
748 // system with the sparse LU decomposition provided by UMFPACK. This is
749 // often a good first choice for 2D problems and works reasonably well even
750 // for a large number of DoFs. The deal.II interface to UMFPACK is given by
751 // the SparseDirectUMFPACK class, which is very easy to use and allows us to
752 // solve our linear system with just 3 lines of code.
753
754 // Note again that for compiling this example program, you need to have the
755 // deal.II library built with UMFPACK support.
756 template <int dim>
solve()757 void UltrasoundProblem<dim>::solve()
758 {
759 std::cout << "Solving linear system... ";
760 Timer timer;
761
762 // The code to solve the linear system is short: First, we allocate an
763 // object of the right type. The following <code>initialize</code> call
764 // provides the matrix that we would like to invert to the
765 // SparseDirectUMFPACK object, and at the same time kicks off the
766 // LU-decomposition. Hence, this is also the point where most of the
767 // computational work in this program happens.
768 SparseDirectUMFPACK A_direct;
769 A_direct.initialize(system_matrix);
770
771 // After the decomposition, we can use <code>A_direct</code> like a matrix
772 // representing the inverse of our system matrix, so to compute the
773 // solution we just have to multiply with the right hand side vector:
774 A_direct.vmult(solution, system_rhs);
775
776 timer.stop();
777 std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
778 }
779
780
781
782 // @sect4{<code>UltrasoundProblem::output_results</code>}
783
784 // Here we output our solution $v$ and $w$ as well as the derived quantity
785 // $|u|$ in the format specified in the parameter file. Most of the work for
786 // deriving $|u|$ from $v$ and $w$ was already done in the implementation of
787 // the <code>ComputeIntensity</code> class, so that the output routine is
788 // rather straightforward and very similar to what is done in the previous
789 // tutorials.
790 template <int dim>
output_results() const791 void UltrasoundProblem<dim>::output_results() const
792 {
793 std::cout << "Generating output... ";
794 Timer timer;
795
796 // Define objects of our <code>ComputeIntensity</code> class and a DataOut
797 // object:
798 ComputeIntensity<dim> intensities;
799 DataOut<dim> data_out;
800
801 data_out.attach_dof_handler(dof_handler);
802
803 // Next we query the output-related parameters from the ParameterHandler.
804 // The DataOut::parse_parameters call acts as a counterpart to the
805 // DataOutInterface<1>::declare_parameters call in
806 // <code>ParameterReader::declare_parameters</code>. It collects all the
807 // output format related parameters from the ParameterHandler and sets the
808 // corresponding properties of the DataOut object accordingly.
809 prm.enter_subsection("Output parameters");
810
811 const std::string output_filename = prm.get("Output filename");
812 data_out.parse_parameters(prm);
813
814 prm.leave_subsection();
815
816 // Now we put together the filename from the base name provided by the
817 // ParameterHandler and the suffix which is provided by the DataOut class
818 // (the default suffix is set to the right type that matches the one set
819 // in the .prm file through parse_parameters()):
820 const std::string filename = output_filename + data_out.default_suffix();
821
822 std::ofstream output(filename);
823
824 // The solution vectors $v$ and $w$ are added to the DataOut object in the
825 // usual way:
826 std::vector<std::string> solution_names;
827 solution_names.emplace_back("Re_u");
828 solution_names.emplace_back("Im_u");
829
830 data_out.add_data_vector(solution, solution_names);
831
832 // For the intensity, we just call <code>add_data_vector</code> again, but
833 // this with our <code>ComputeIntensity</code> object as the second
834 // argument, which effectively adds $|u|$ to the output data:
835 data_out.add_data_vector(solution, intensities);
836
837 // The last steps are as before. Note that the actual output format is now
838 // determined by what is stated in the input file, i.e. one can change the
839 // output format without having to re-compile this program:
840 data_out.build_patches();
841 data_out.write(output);
842
843 timer.stop();
844 std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
845 }
846
847
848
849 // @sect4{<code>UltrasoundProblem::run</code>}
850
851 // Here we simply execute our functions one after the other:
852 template <int dim>
run()853 void UltrasoundProblem<dim>::run()
854 {
855 make_grid();
856 setup_system();
857 assemble_system();
858 solve();
859 output_results();
860 }
861 } // namespace Step29
862
863
864 // @sect4{The <code>main</code> function}
865
866 // Finally the <code>main</code> function of the program. It has the same
867 // structure as in almost all of the other tutorial programs. The only
868 // exception is that we define ParameterHandler and
869 // <code>ParameterReader</code> objects, and let the latter read in the
870 // parameter values from a textfile called <code>step-29.prm</code>. The
871 // values so read are then handed over to an instance of the UltrasoundProblem
872 // class:
main()873 int main()
874 {
875 try
876 {
877 using namespace dealii;
878 using namespace Step29;
879
880 ParameterHandler prm;
881 ParameterReader param(prm);
882 param.read_parameters("step-29.prm");
883
884 UltrasoundProblem<2> ultrasound_problem(prm);
885 ultrasound_problem.run();
886 }
887 catch (std::exception &exc)
888 {
889 std::cerr << std::endl
890 << std::endl
891 << "----------------------------------------------------"
892 << std::endl;
893 std::cerr << "Exception on processing: " << std::endl
894 << exc.what() << std::endl
895 << "Aborting!" << std::endl
896 << "----------------------------------------------------"
897 << std::endl;
898 return 1;
899 }
900 catch (...)
901 {
902 std::cerr << std::endl
903 << std::endl
904 << "----------------------------------------------------"
905 << std::endl;
906 std::cerr << "Unknown exception!" << std::endl
907 << "Aborting!" << std::endl
908 << "----------------------------------------------------"
909 << std::endl;
910 return 1;
911 }
912 return 0;
913 }
914