fields representing real and imaginary parts. We also revisit the
ParameterHandler class first used in step-19, which provides a
convenient way for reading parameters from a configuration file at runtime
-without the need to recompile the program code.
+without the need to recompile the program code.
The equations covered here fall into the class of vector-valued problems. A
-toplevel overview of this topic can be found in the @ref vector_valued module.
+top level overview of this topic can be found in the @ref vector_valued module.
<h3>Problem setting</h3>
-The original purpose of this program is to simulate the focussing properties
+The original purpose of this program is to simulate the focusing properties
of an ultrasound wave generated by a transducer lens with variable
geometry. Recent applications in medical imaging use ultrasound waves not only
for imaging purposes, but also to excite certain local effects in a
other imaging techniques. A vital ingredient for these methods is the ability
to focus the intensity of the ultrasound wave in a particular part of the
material, ideally in a point, to be able to examine the properties of the
-material at that particular location.
+material at that particular location.
To derive a model for this problem, we think of ultrasound as a pressure wave
-governed by the wave equation:
+governed by the wave equation:
@f[
\frac{\partial^2 U}{\partial t^2} - c^2 \Delta U = 0
@f]
transducer lens and $\Gamma_2$ an absorbing boundary (that is, we want to
choose boundary conditions on $\Gamma_2$ in such a way that they imitate a
larger domain). On $\Gamma_1$, the transducer generates a wave of constant
-frequency ${\omega}>0$ and constant amplitude (that we chose to be 1 here):
+frequency ${\omega}>0$ and constant amplitude (that we chose to be 1 here):
@f[
U(x,t) = \cos{\omega t}, \qquad x\in \Gamma_1
@f]
we see that for $u$ we have
@f{eqnarray*}
-\omega^2 u(x) - c^2\Delta u(x) &=& 0, \qquad x\in\Omega,\\
-u(x) &=& 1, \qquad x\in\Gamma_1.
+u(x) &=& 1, \qquad x\in\Gamma_1.
@f}
For finding suitable conditions on $\Gamma_2$ that model an absorbing
for $V$ to solve the wave equation, $|k|={\frac{\omega}{c}}$ must
hold. Suppose that this wave hits the boundary in $x_0\in\Gamma_2$ at a right
angle, i.e. $n=\frac{k}{|k|}$ with $n$ denoting the outer unit normal of
-$\Omega$ in $x_0$. Then at $x_0$, this wave satisfies the equation
+$\Omega$ in $x_0$. Then at $x_0$, this wave satisfies the equation
@f[
c (n\cdot\nabla V) + \frac{\partial V}{\partial t} = (i\, c\, |k| - i\, \omega) V = 0.
@f]
-Hence, by enforcing the boundary condition
+Hence, by enforcing the boundary condition
@f[
c (n\cdot\nabla U) + \frac{\partial U}{\partial t} = 0, \qquad x\in\Gamma_2,
@f]
will pass through the boundary as if it wasn't here whereas the remaining
fraction of the wave will be reflected back into the domain.
-If we are willing to accept this as a sufficient approximation to an absorbing boundary we finally arrive at the following problem for $u$:
+If we are willing to accept this as a sufficient approximation to an absorbing boundary we finally arrive at the following problem for $u$:
@f{eqnarray*}
-\omega^2 u - c^2\Delta u &=& 0, \qquad x\in\Omega,\\
c (n\cdot\nabla u) + i\,\omega\,u &=&0, \qquad x\in\Gamma_2,\\
-u &=& 1, \qquad x\in\Gamma_1.
+u &=& 1, \qquad x\in\Gamma_1.
@f}
This is a Helmholtz equation (similar to the one in step-7, but this time with ''the bad sign'') with Dirichlet data on $\Gamma_1$ and mixed boundary conditions on $\Gamma_2$. Because of the condition on $\Gamma_2$, we cannot just treat the equations for real and imaginary parts of $u$ separately. What we can do however is to view the PDE for $u$ as a system of two PDEs for the real and imaginary parts of $u$, with the boundary condition on $\Gamma_2$ representing the coupling terms between the two components of the system. This works along the following lines: Let $v=\textrm{Re}\;u,\; w=\textrm{Im}\;u$, then in terms of $v$ and $w$ we have the following system:
@f{eqnarray*}
\left.\begin{array}{ccc}
- -\omega^2 v - c^2\Delta v &=& 0 \quad\\
+ -\omega^2 v - c^2\Delta v &=& 0 \quad\\
-\omega^2 w - c^2\Delta w &=& 0 \quad
\end{array}\right\} &\;& x\in\Omega,
\\
\left.\begin{array}{ccc}
- c (n\cdot\nabla v) - \omega\,w &=& 0 \quad\\
+ c (n\cdot\nabla v) - \omega\,w &=& 0 \quad\\
c (n\cdot\nabla w) + \omega\,v &=& 0 \quad
\end{array}\right\} &\;& x\in\Gamma_2,
\\
\left.\begin{array}{ccc}
- v &=& 1 \quad\\
+ v &=& 1 \quad\\
w &=& 0 \quad
\end{array}\right\} &\;& x\in\Gamma_1.
@f}
-\omega^2 \langle \psi, w \rangle_{\mathrm{L}^2(\Omega)} + c^2 \langle \nabla \psi, \nabla w \rangle_{\mathrm{L}^2(\Omega)} + c \omega \langle \psi, v \rangle_{\mathrm{L}^2(\Gamma_2)} &=& 0.
@f}
-We choose finite element spaces $V_h$ and $W_h$ with bases $\{\phi_j\}_{j=1}^n, \{\psi_j\}_{j=1}^n$ and look for approximate solutions
+We choose finite element spaces $V_h$ and $W_h$ with bases $\{\phi_j\}_{j=1}^n, \{\psi_j\}_{j=1}^n$ and look for approximate solutions
@f[
-v_h = \sum_{j=1}^n \alpha_j \phi_j, \;\; w_h = \sum_{j=1}^n \beta_j \psi_j.
+v_h = \sum_{j=1}^n \alpha_j \phi_j, \;\; w_h = \sum_{j=1}^n \beta_j \psi_j.
@f]
Plugging into the variational form yields the equation system
@f[
\renewcommand{\arraystretch}{2.0}
\left.\begin{array}{ccc}
-\sum_{j=1}^n \left(-\omega^2 \langle \phi_i, \phi_j \rangle_{\mathrm{L}^2(\Omega)} +c^2 \langle \nabla \phi_i, \nabla \phi_j \rangle_{\mathrm{L}^2(\Omega)}\right)\alpha_j - \left(c\omega \langle \phi_i,\psi_j\rangle_{\mathrm{L}^2(\Gamma_2)}\right)\beta_j &=& 0 \\
+\sum_{j=1}^n \left(-\omega^2 \langle \phi_i, \phi_j \rangle_{\mathrm{L}^2(\Omega)} +c^2 \langle \nabla \phi_i, \nabla \phi_j \rangle_{\mathrm{L}^2(\Omega)}\right)\alpha_j - \left(c\omega \langle \phi_i,\psi_j\rangle_{\mathrm{L}^2(\Gamma_2)}\right)\beta_j &=& 0 \\
\sum_{j=1}^n \left(-\omega^2 \langle \psi_i, \psi_j \rangle_{\mathrm{L}^2(\Omega)} +c^2 \langle \nabla \psi_i, \nabla \psi_j \rangle_{\mathrm{L}^2(\Omega)}\right)\beta_j + \left(c\omega \langle \psi_i,\phi_j\rangle_{\mathrm{L}^2(\Gamma_2)}\right)\alpha_j &=& 0
\end{array}\right\}\;\;\forall\; i =1,\ldots,n.
@f]
-In matrix notation:
+In matrix notation:
@f[
\renewcommand{\arraystretch}{2.0}
\left(
\end{array}
\right)
@f]
-(One should not be fooled by the right hand side being zero here, that is
-because we haven't included the Dirichlet boundary data yet.)
-Because of the alternating sign in the off-diagonal blocks, we can already
+(One should not be fooled by the right hand side being zero here, that is
+because we haven't included the Dirichlet boundary data yet.)
+Because of the alternating sign in the off-diagonal blocks, we can already
see that this system is non-symmetric, in fact it is even indefinite.
Of course, there is no necessity to choose the spaces $V_h$ and $W_h$ to be
the same. However, we expect real and imaginary part of the solution to
-have similar properties and will therefore indeed take $V_h=W_h$ in the
-implementation, and also use the same basis functions $\phi_i = \psi_i$ for
-both spaces. The reason for the notation using different symbols is just that
-it allows us to distinguish between shape functions for $v$ and $w$, as this
-distinction plays an important role in the implementation.
+have similar properties and will therefore indeed take $V_h=W_h$ in the
+implementation, and also use the same basis functions $\phi_i = \psi_i$ for
+both spaces. The reason for the notation using different symbols is just that
+it allows us to distinguish between shape functions for $v$ and $w$, as this
+distinction plays an important role in the implementation.
<h3>The test case</h3>
-For the computations, we will consider wave propagation in the unit square,
-with ultrasound generated by a transducer lens that is shaped like a segment
-of the circle with center at $(0.5, d)$ and a
+For the computations, we will consider wave propagation in the unit square,
+with ultrasound generated by a transducer lens that is shaped like a segment
+of the circle with center at $(0.5, d)$ and a
radius slightly greater than $d$; this shape should lead to a focusing of the sound
wave at the center of the circle. Varying $d$ changes the "focus" of the lens
and affects the spatial distribution of the intensity of $u$, where our main
-concern is how well $|u|=\sqrt{v^2+w^2}$ is focused.
+concern is how well $|u|=\sqrt{v^2+w^2}$ is focused.
In the program below, we will implement the complex-valued Helmholtz equations
using the formulation with split real and imaginary parts. We will also
discuss how to generate a domain that looks like a square with a slight bulge
simulating the transducer (in the
<code>UltrasoundProblem<dim>::make_grid()</code> function), and how to
-generate graphical output that not only contains the solution components $v$ and
+generate graphical output that not only contains the solution components $v$ and
$w$, but also the magnitude $\sqrt{v^2+w^2}$ directly in the output file (in
-<code>UltrasoundProblem<dim>::output_results()</code>). Finally, we use the
-ParameterHandler class to easily read parameters like the focal distance $d$,
-wave speed $c$, frequency $\omega$, and a number of other parameters from an
-input file at run-time, rather than fixing those parameters in the source code
+<code>UltrasoundProblem<dim>::output_results()</code>). Finally, we use the
+ParameterHandler class to easily read parameters like the focal distance $d$,
+wave speed $c$, frequency $\omega$, and a number of other parameters from an
+input file at run-time, rather than fixing those parameters in the source code
where we would have to re-compile every time we want to change parameters.
Here's the console output of the program in debug mode:
@code
-examples/\step-29> make run
-============================ Running \step-29
-DEAL::Generating grid... done (1.11607s)
+> make run
+[ 66%] Built target step-29
+[100%] Run step-29 with Debug configuration
+DEAL::Generating grid... done (0.832682s)
DEAL:: Number of active cells: 25600
-DEAL::Setting up system... done (1.10807s)
+DEAL::Setting up system... done (0.614761s)
DEAL:: Number of degrees of freedom: 51842
-DEAL::Assembling system matrix... done (6.50841s)
-DEAL::Solving linear system... done (2.89218s)
-DEAL::Generating output... done (4.52428s)
+DEAL::Assembling system matrix... done (2.82536s)
+DEAL::Solving linear system... done (2.27971s)
+DEAL::Generating output... done (1.84418s)
+[100%] Built target run
@endcode
(Of course, execution times will differ if you run the program
locally.) The fact that most of the time is spent on assembling
-the system matrix and generating output is due to the many assertion
+the system matrix and generating output is due to the many assertions
that need to be checked in debug mode. In release mode these parts
of the program run much faster whereas solving the linear system is
hardly sped up at all:
@code
-============================ Running \step-29
-DEAL::Generating grid... done (0.0280020s)
+> make run
+[ 66%] Built target step-29
+Scanning dependencies of target run
+[100%] Run step-29 with Release configuration
+DEAL::Generating grid... done (0.0144960s)
DEAL:: Number of active cells: 25600
-DEAL::Setting up system... done (0.112007s)
+DEAL::Setting up system... done (0.0356880s)
DEAL:: Number of degrees of freedom: 51842
-DEAL::Assembling system matrix... done (0.160010s)
-DEAL::Solving linear system... done (2.21614s)
-DEAL::Generating output... done (0.880055s)
+DEAL::Assembling system matrix... done (0.0436570s)
+DEAL::Solving linear system... done (1.54733s)
+DEAL::Generating output... done (0.720528s)
+[100%] Built target run
@endcode
The graphical output of the program looks as follows:
global refinements:
@code
-DEAL::Generating grid... done (0.0320020s)
+> make run
+[ 66%] Built target step-29
+[100%] Run step-29 with Release configuration
+DEAL::Generating grid... done (0.0135260s)
DEAL:: Number of active cells: 25600
-DEAL::Setting up system... done (0.104006s)
+DEAL::Setting up system... done (0.0213910s)
DEAL:: Number of degrees of freedom: 51842
-DEAL::Assembling system matrix... done (0.164011s)
-DEAL::Solving linear system... done (2.23214s)
-DEAL::Generating output... done (0.900056s)
-
-DEAL::Generating grid... done (0.132009s)
+DEAL::Assembling system matrix... done (0.0414300s)
+DEAL::Solving linear system... done (1.56621s)
+DEAL::Generating output... done (0.729605s)
+[100%] Built target run
+
+> make run
+[ 66%] Built target step-29
+[100%] Run step-29 with Release configuration
+DEAL::Generating grid... done (0.0668490s)
DEAL:: Number of active cells: 102400
-DEAL::Setting up system... done (0.408025s)
+DEAL::Setting up system... done (0.109694s)
DEAL:: Number of degrees of freedom: 206082
-DEAL::Assembling system matrix... done (0.656041s)
-DEAL::Solving linear system... done (14.8849s)
-DEAL::Generating output... done (3.57222s)
-
-DEAL::Generating grid... done (0.504031s)
+DEAL::Assembling system matrix... done (0.160784s)
+DEAL::Solving linear system... done (7.86577s)
+DEAL::Generating output... done (2.89320s)
+[100%] Built target run
+
+> make run
+[ 66%] Built target step-29
+[100%] Run step-29 with Release configuration
+DEAL::Generating grid... done (0.293154s)
DEAL:: Number of active cells: 409600
-DEAL::Setting up system... done (1.72011s)
+DEAL::Setting up system... done (0.491301s)
DEAL:: Number of degrees of freedom: 821762
-DEAL::Assembling system matrix... done (2.63216s)
-DEAL::Solving linear system... done (117.811s)
-DEAL::Generating output... done (15.1489s)
+DEAL::Assembling system matrix... done (0.605386s)
+DEAL::Solving linear system... done (45.1989s)
+DEAL::Generating output... done (11.2292s)
+[100%] Built target run
@endcode
Each time we refine the mesh once, so the number of cells and degrees
// @sect3{Include files}
-// The following header files are unchanged from step-7 and have been
-// discussed before:
+// The following header files have all been discussed before:
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
: Function<dim>(2)
{}
- virtual void vector_value(const Point<dim> &p,
- Vector<double> & values) const override;
+ virtual void vector_value(const Point<dim> & /*p*/,
+ Vector<double> &values) const override
+ {
+ Assert(values.size() == 2, ExcDimensionMismatch(values.size(), 2));
+
+ values(0) = 1;
+ values(1) = 0;
+ }
virtual void
vector_value_list(const std::vector<Point<dim>> &points,
- std::vector<Vector<double>> & value_list) const override;
- };
-
-
- template <int dim>
- inline void
- DirichletBoundaryValues<dim>::vector_value(const Point<dim> & /*p*/,
- Vector<double> &values) const
- {
- Assert(values.size() == 2, ExcDimensionMismatch(values.size(), 2));
-
- values(0) = 1;
- values(1) = 0;
- }
-
-
- template <int dim>
- void DirichletBoundaryValues<dim>::vector_value_list(
- const std::vector<Point<dim>> &points,
- std::vector<Vector<double>> & value_list) const
- {
- Assert(value_list.size() == points.size(),
- ExcDimensionMismatch(value_list.size(), points.size()));
+ std::vector<Vector<double>> & value_list) const override
+ {
+ Assert(value_list.size() == points.size(),
+ ExcDimensionMismatch(value_list.size(), points.size()));
- for (unsigned int p = 0; p < points.size(); ++p)
- DirichletBoundaryValues<dim>::vector_value(points[p], value_list[p]);
- }
+ for (unsigned int p = 0; p < points.size(); ++p)
+ DirichletBoundaryValues<dim>::vector_value(points[p], value_list[p]);
+ }
+ };
// @sect3{The <code>ParameterReader</code> class}
{
public:
UltrasoundProblem(ParameterHandler &);
- ~UltrasoundProblem();
void run();
private:
, fe(FE_Q<dim>(1), 2)
{}
-
- template <int dim>
- UltrasoundProblem<dim>::~UltrasoundProblem()
- {
- dof_handler.clear();
- }
-
// @sect4{<code>UltrasoundProblem::make_grid</code>}
// Here we setup the grid for our domain. As mentioned in the exposition,
// boundary indicator.
GridGenerator::subdivided_hyper_cube(triangulation, 5, 0, 1);
- typename Triangulation<dim>::cell_iterator cell = triangulation.begin(),
- endc = triangulation.end();
-
- for (; cell != endc; ++cell)
+ for (auto &cell : triangulation.cell_iterators())
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
- typename DoFHandler<dim>::active_cell_iterator cell =
- dof_handler.begin_active(),
- endc = dof_handler.end();
-
- for (; cell != endc; ++cell)
+ for (const auto &cell : dof_handler.active_cell_iterators())
{
// On each cell, we first need to reset the local contribution matrix
// and request the FEValues object to compute the shape functions for
projects / dealii.git / commitdiff