<a name="Intro"></a>
<h1>Introduction</h1>
-This project is to simulate the thermoacoustic tomography imaging. In thermoacoustic
-tomography,pulsed electromagnetic energy is delivered into biological issues.
-Tissues absorbe the energy and then generate thermoacoustic waves through
-thermoelastic expansion. The forward problem is a wave propagation problem.
-
-
-<h3>Problem</h3>
-
-Thermal equation without considering thermal diffusion is
+This program grew out of a student project by Xing Jin at Texas A&M
+University. Most of the work for this program is by her.
+
+The program is part of a project that aims to simulate thermoacoustic
+tomography imaging. In thermoacoustic tomography, pulsed electromagnetic
+energy is delivered into biological issues. Tissues absorb some of this energy
+and those parts of the tissue that absorb the most energy generate
+thermoacoustic waves through thermoelastic expansion. For imaging, one uses
+that different kinds of tissue, most importantly healthy and diseased tissue,
+absorb different amounts of energy and therefore expand at different
+rates. The experimental setup is to measure the amplitude of the pressure
+waves generated by these sources on the surface of the tissue and try to
+reconstruct the source distributions, which is indicative for the distribution
+of absorbers and therefore of different kinds of tissue. Part of this project
+is to compare simulated data with actual measurements, so one has to solve the
+"forward problem", i.e. the wave equation that describes the propagation of
+pressure waves in tissue. This program is therefore a continuation of @ref
+step_23 "step-23", where the wave equation was first introduced.
+
+
+<h3>The problem</h3>
+
+The temperature at a given location, neglecting thermal diffusion, can be
+stated as
@f[
\rho C_p \frac{\partial}{\partial t}T(t,\mathbf r) = H(t,\mathbf r)
@f]
-Where $\rho (\mathbf r) $ is the density; $C_p (\mathbf r) $ is the specific heat;
-$ T(t,\mathbf r)$ is the temperature rise due to the delivered microwave
-energy; and $H(t,\mathbf r)$ is the heating function defined as the thermal energy
-per time and volume transformed from deposited microwave energy.
-
-Assume tissues have heterogeneous dielectric properties but homogeneous acoustic
-properties. The basic acoustic generation equation in an acoustically homogeneous
-medium is the linear inviscid force equation
-
+Here $\rho (\mathbf r) $ is the density; $C_p (\mathbf r) $ is the specific
+heat; $\frac{\partial T}{\partial t}(t,\mathbf r)$ is the temperature rise due
+to the delivered microwave energy; and $H(t,\mathbf r)$ is the heating
+function defined as the thermal energy per time and volume transformed from
+deposited microwave energy.
+
+Let us assume that tissues have heterogeneous dielectric properties but
+homogeneous acoustic properties. The basic acoustic generation equation in an
+acoustically homogeneous medium can be described as follows: if $u$ is the
+vector-valued displacement, then tissue certainly reacts to changes in
+pressure by accelleration:
@f[
\rho \frac{\partial^2}{\partial t^2}u(t,\mathbf r) =
--\nabla p(t,\mathbf r)
+-\nabla p(t,\mathbf r).
@f]
-
-and the expansion equation:
-
+Furthermore, it expands based on changes in temperature:
@f[
\nabla \cdot u(t,\mathbf r) = -\frac{p(t,\mathbf r)}{\rho c_0^2}+\beta T(t,\mathbf r)
@f]
-The original problem can be described as:
-
+If we combine these equations and assume that heating only happens on a time
+scale much shorter than wave propagation through tissue (i.e. the temporal
+length of the microwave pulse that heats the tissue is much shorter than the
+time it takes a wave to cross the domain), then we can rewrite the above
+equations as follows:
@f[
-\Delta p-\frac{1}{c_0^2} \frac{\partial^2 p }{\partial^2 t} = \lambda \delta(t)a(\mathbf r)
+\Delta p-\frac{1}{c_0^2} \frac{\partial^2 p}{\partial^2 t} = \lambda \delta(t)a(\mathbf r)
@f]
-
-where $\lambda = - \frac{\beta}{C_p}$.
-
-The forward propogation problem can be changed to solve a wave equation with
-initial conditions as follows:
-
+where $\lambda = - \frac{\beta}{C_p}$. This corresponds to a wave equation
+with initial conditions as follows:
@f{eqnarray*}
-\Delta \bar{p}- \frac{1}{c_0^2} \frac{\partial^2 \bar{p}}{\partial^2 t} & = & f(t,\mathbf r) \\
-
-\bar{p}(0,\mathbf r)=\lambda a(\mathbf r) & = & b(\mathbf r)
+\Delta \bar{p}- \frac{1}{c_0^2} \frac{\partial^2 \bar{p}}{\partial^2 t} & = &
+f(t,\mathbf r) \\
+\bar{p}(0,\mathbf r) &=&\lambda a(\mathbf r) = b(\mathbf r)
@f}
+In the inverse problem, it is this right hand side $\lambda a(\mathbf r)$ that
+one would like to recover, since it is a map of absorption strengths for
+microwave energy, and therefore presumably an indicator to discern healthy
+from diseased tissue.
+
+In real application, the thermoacoustic source is very small as compared to
+the medium. The propagation path of the thermoacoustic waves can then be
+approximated as from the source to the infinity. Furthermore, detectors are
+only a limited distance from the source. One only needs to evaluate the values
+when the thermoacoustic waves pass through the detectors, although they do
+continue beyond. This is therefore a problem where we are only interested in a
+small part of an infinite medium, and we do not want waves generated somewhere
+to be reflected at the boundary of the domain which we consider
+interesting. Rather, we would like to simulate only that part of the wave
+field that is contained inside the domain of interest, and waves that hit the
+boundary of that domain to simply pass undisturbed through the boundary. In
+other words, we would like the boundary to absorb any waves that hit it.
+
+In general, this is a hard problem: Good absorbing boundary conditions are
+nonlinear and/or numerically very expensive. We therefore opt for a simple
+first order approximation to absorbing boundary conditions that reads
+@f[
+\frac{\partial\bar{p}}{\partial\mathbf n} =
+-\frac{1}{c_0} \frac{\partial\bar{p}}{\partial t}
+@f]
+Here, $\frac{\partial\bar{p}}{\partial\mathbf n}$ is the normal derivative at
+the boundary. It should be noted that this is not a particularly good boundary
+condition, but it is one of the very few that are reasonably simple to implement.
+<h3>Weak form and discretization</h3>
-<h3>Weak form and Discretization</h3>
-
-One first introduces a second variable, which is defined as the derivative of
-the pressure potential.
-
+As in @ref step_23 "step-23", one first introduces a second variable, which is
+defined as the derivative of the pressure potential:
@f[
v = \frac{\partial\bar{p}}{\partial t}
@f]
With the second variables, one then transform the forward problem into
two seperate equations:
-
@f{eqnarray*}
\bar{p}_{t} - v & = & 0 \\
\Delta\bar{p} - \frac{1}{c_0^2}\,v_{t} & = & f
@f}
-
with initial conditions:
-
@f{eqnarray*}
\bar{p}(0,\mathbf r) & = & b(r) \\
v(0,\mathbf r)=\bar{p}_t(0,\mathbf r) & = & 1
@f}
-In real application, the thermoacoustic source is very small as compared to the medium.
-The propagation path of the thermoacoustic waves can be approximated as from the source
-to the infinity. And the detector is in limited distance from the source. One only needs to
-evaluate the values when the thermoacoustic waves pass through the detectors. For this specific
-detection geometry, One then chooses the absorbing boundary condition for the simulation.
-
-@f[
-\frac{\partial\bar{p}}{\partial\mathbf n} =
--\frac{1}{c_0} \frac{\partial\bar{p}}{\partial t}
-@f]
-
-$\frac{\partial\bar{p}}{\partial\mathbf n}$ is normal derivative at the boundary. This is the
-time-varying FEM model. To implement FEM for time dependent problem , one discretizes according
-to $t$ and obtains:
-
+The semi-discretized, weak version of this model, using the general $\theta$ scheme
+introduced in @ref step_23 "step-23" is then:
@f{eqnarray*}
-(\frac{\bar{p}^n-\bar{p}^{n-1}}{\delta t},\phi)_\Omega-
-(\theta v^{n}+(1-\theta)v^{n-1},\phi)_\Omega & = & 0 \\
--(\Delta((\theta\bar{p}^n+(1-\theta)\bar{p}^{n-1})),\nabla)_\Omega-
-\frac{1}{c_0}(\frac{\bar{p}^n-\bar{p}^{n-1}}{\delta t},\phi)_\partial\Omega -
-\frac{1}{c_0^2}(\frac{v^n-v^{n-1}}{\delta t},\phi)_\Omega & =
-& \theta f^{n}+(1-\theta)f^{n-1}
+\left(\frac{\bar{p}^n-\bar{p}^{n-1}}{k},\phi\right)_\Omega-
+\left(\theta v^{n}+(1-\theta)v^{n-1},\phi\right)_\Omega & = & 0 \\
+-\left(\nabla((\theta\bar{p}^n+(1-\theta)\bar{p}^{n-1})),\nabla\phi\right)_\Omega-
+\frac{1}{c_0}\left(\frac{\bar{p}^n-\bar{p}^{n-1}}{k},\phi\right)_{\partial\Omega} -
+\frac{1}{c_0^2}\left(\frac{v^n-v^{n-1}}{k},\phi\right)_\Omega & =
+& \theta f^{n}+(1-\theta)f^{n-1},
@f}
+where $\phi$ is an arbitrary test function, and where we have used the
+absorbing boundary condition to integrate by parts:
+absoring boundary conditions are incorporated into the weak form by using
+@f[
+\int_\Omega\varphi \, \Delta p\; dx =
+-\int_\Omega\nabla \varphi \cdot \nabla p dx +
+\int_{\partial\Omega}\varphi \frac{\partial p}{\partial t}ds.
+@f]
-The weak formulation of the problem is obtained by multiplying the above two equations
-with test functions and integrating some terms by parts:
-
+From this we obtain the discrete model by introducing a finite number of shape
+functions, and get
@f{eqnarray*}
-M\bar{p}^{n}-(\delta t \theta)M v^{n-1} & = & M\bar{p}^{n-1}+\delta t (1-\theta)Mv^{n-1}\\
+M\bar{p}^{n}-k \theta M v^{n-1} & = & M\bar{p}^{n-1}+k (1-\theta)Mv^{n-1},\\
-(-c_0^2\delta t \theta A-c_0 B)\bar{p}^n-Mv^{n} & = &
-(c_0^2\delta t(1-\theta)A-c_0B)\bar{p}^{n-1}-Mv^{n-1}+c_0^2\delta t(\theta F^{n}+(1-\theta)F^{n-1})
+(-c_0^2k \theta A-c_0 B)\bar{p}^n-Mv^{n} & = &
+(c_0^2k(1-\theta)A-c_0B)\bar{p}^{n-1}-Mv^{n-1}+c_0^2k(\theta F^{n}+(1-\theta)F^{n-1}).
@f}
-
-Here, the absoring boundary conditions are incorporated into the weak form by using
-
-@f[
-\int_\Omega\varphi(\nabla\cdot(\nabla p))dx =
--\int_\Omega\nabla \varphi \cdot \nabla p dx +
-\int_{\partial\Omega}\varphi \frac{\partial p}{\partial t}ds
+The matrices $M$ and $A$ are here as in @ref step_23 "step-23", and the
+boundary mass matrix
+@f[
+ B_{ij} = \left(\varphi_i,\varphi_j\right)_{\partial\Omega}
@f]
+results from the use of absorbing boundary conditions.
-Where $\varphi$ is the test function.
-
-Pressure and its derivative are the two variables one is interested in the above equations,
-One can write the above two equations as a matrix form with the pressure and its derivative as
+Above two equations can be rewritten in a matrix form with the pressure and its derivative as
an unknown vecotor:
@f[
\left(\begin{array}{cc}
- M & -\delta t\theta M \\
-c_0^2\,\delta t\,\theta\,A+c_0\,B & M \\
+ M & -k\theta M \\
+c_0^2\,k\,\theta\,A+c_0\,B & M \\
\end{array} \right)\\
\left(\begin{array}{c}
\bar{p}^{n} \\
\end{array}\right)=\\
\left(\begin{array}{l}
G_1 \\
- G_2 -(\theta F^{n}+(1-\theta)F ^{n-1})c_{0}^{2}\delta t \\
+ G_2 -(\theta F^{n}+(1-\theta)F ^{n-1})c_{0}^{2}k \\
\end{array}\right)
@f]
G_2 \\
\end{array} \right)=\\
\left(\begin{array}{l}
- M\bar{p}^{n-1}+\delta t(1-\theta)Mv^{n-1}\\
- (-c_{0}^{2}\delta t (1-\theta)A+c_0 B)\bar{p}^{n-1} +Mv^{n-1}
+ M\bar{p}^{n-1}+k(1-\theta)Mv^{n-1}\\
+ (-c_{0}^{2}k (1-\theta)A+c_0 B)\bar{p}^{n-1} +Mv^{n-1}
\end{array}\right)
@f]
-By some simply transformation, one obtains two iterative equations for
-the pressure potential and its derivative:
+By simple transformations, one then obtains two equations for
+the pressure potential and its derivative, just as in the previous tutorial program:
@f{eqnarray*}
-(M+(\delta t\,\theta\,c_{0})^{2}A+c_0\delta t\theta B)\bar{p}^{n} & = &
-G_{1}+(\delta t\, \theta)G_{2}-(c_0\delta t)^2\theta (\theta F^{n}+(1-\theta)F^{n-1}) \\
-Mv^n & = & -(c_0^2\,\delta t\, \theta\, A+c_0B)\bar{p}^{n}+ G_2 -
-c_0^2\delta t(\theta F^{n}+(1-\theta)F^{n-1})
+(M+(k\,\theta\,c_{0})^{2}A+c_0k\theta B)\bar{p}^{n} & = &
+G_{1}+(k\, \theta)G_{2}-(c_0k)^2\theta (\theta F^{n}+(1-\theta)F^{n-1}) \\
+Mv^n & = & -(c_0^2\,k\, \theta\, A+c_0B)\bar{p}^{n}+ G_2 -
+c_0^2k(\theta F^{n}+(1-\theta)F^{n-1})
@f}
-
-/* $Id: project.cc modified from heat-equation.cc 2006/03/05 $ */
-/* Author: Xing Jin */
-/* */
-/* $Id: step-4.cc,v 1.34 2006/02/06 21:33:10 wolf Exp $ */
+/* $Id: step-4.cc,v 1.34 2006/02/06 21:33:10 wolf Exp $ */
/* Version: $Name: $ */
/* */
-/* Copyright (c) 1999,2000,2001,2002,2003,2004,2005,2006 */
-/* by the deal.II authors. */
+/* Copyright (C) 2006 by the deal.II authors */
+/* Author: Xing Jin, Wolfgang Bangerth, Texas A&M University, 2006 */
+/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
/* to the file deal.II/doc/license.html for the text and */
/* further information on this license. */
- // @sect3{Include files}
- // Most include files have been covered in
- // step-6 and will not be further commented on
-#include <grid/tria.h>
-#include <dofs/dof_handler.h>
- // We will need to read the value at a specific
- // location. This including file is needed for
- // finding a cell that contains a given point
-#include <grid/grid_tools.h>
-#include <grid/grid_generator.h>
- // Because the scanning geometry is on a circle,
- // the boundaries are not straight lines, so
- // we need some classes to predefine some
- // boundary description
-#include <grid/tria_boundary_lib.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_iterator.h>
-#include <dofs/dof_accessor.h>
-#include <fe/fe_q.h>
-#include <dofs/dof_tools.h>
-#include <fe/fe_values.h>
-
-#include <fe/mapping_q1.h>
+ // @sect3{Include files}
+
+ // The following have all been covered
+ // previously:
#include <base/quadrature_lib.h>
#include <base/function.h>
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
+#include <base/logstream.h>
+#include <base/utilities.h>
+
#include <lac/vector.h>
#include <lac/full_matrix.h>
#include <lac/sparse_matrix.h>
#include <lac/solver_cg.h>
#include <lac/precondition.h>
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
#include <dofs/dof_constraints.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+
+#include <numerics/data_out.h>
#include <numerics/matrices.h>
#include <numerics/vectors.h>
-#include <numerics/data_out.h>
- // These are for c++
#include <fstream>
#include <iostream>
#include <sstream>
-
-#include <base/logstream.h>
-
-#include <base/point.h>
- // @sect3{"The forward problem" class template}
- // The main class is similar to the wave equation.
- // The difference is that we add an absorbing
- // boundary condition. Because we are only interested
- // in values at specific locations, we define some
- // parameters to obtain the coordinates of those
- // locations.
+ // @sect3{The "forward problem" class template}
+
+ // The main class is similar to the wave
+ // equation. The difference is that we add
+ // an absorbing boundary condition. Because
+ // we are only interested in values at
+ // specific locations, we define some
+ // parameters to obtain the coordinates of
+ // those locations.
template <int dim>
class TATForwardProblem
{
void run ();
private:
- void make_grid_and_dofs ();
- void assemble_system ();
+ void setup_system ();
void solve_p ();
void solve_v ();
- void output_results (const unsigned int timestep_number) const;
+ void output_results () const;
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
+ ConstraintMatrix constraints;
+
SparsityPattern sparsity_pattern;
- SparseMatrix<double> system_matrix_p;
- SparseMatrix<double> system_matrix_v;
+ SparseMatrix<double> system_matrix;
SparseMatrix<double> mass_matrix;
SparseMatrix<double> laplace_matrix;
+
+ Vector<double> solution_p, solution_v;
+ Vector<double> old_solution_p, old_solution_v;
+ Vector<double> system_rhs_p, system_rhs_v;
+
+ double time, time_step;
+ unsigned int timestep_number;
+ const double theta;
+
+ //
SparseMatrix<double> boundary_matrix;
- // Number of refinement
- unsigned int n_refinements;
- // The acoustic speed in the medium $c_0$
- double acoustic_speed;
- // This parameter is needed for discritizing
- // time-dependent problem
- double theta;
-
- // The total data collection time
- double total_time;
-
- // The size of the time step
- double time_step;
- // The detector circullarly scan the target region.
- // The step size of the detector is in angles
- double step_angle;
- // The scanning radius
- double radius;
-
-
- Vector<double> solution_p;
- Vector<double> old_solution_p;
- Vector<double> system_rhs_p;
-
- Vector<double> solution_v;
- Vector<double> old_solution_v;
- Vector<double> system_rhs_v;
-
-
+ // Number of refinement
+ const unsigned int n_refinements;
+ // The acoustic speed in the medium $c_0$
+ const double acoustic_speed;
+
+ // The detector circullarly scan the target region.
+ // The step size of the detector is in angles
+ const double step_angle;
+ // The scanning radius
+ const double radius;
+
+ const double end_time;
};
- // Declare a class template for the right hand side
- // of the pressure potential
+ // Declare a class template for the right hand side
+ // of the pressure potential
template <int dim>
-class RightHandSide_p : public Function<dim>
+class RightHandSideP : public Function<dim>
{
public:
- RightHandSide_p () : Function<dim>() {};
+ RightHandSideP () : Function<dim>() {};
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
- // Declare a class template for the right hand side
- // of the derivative of the pressure potential
+ // Declare a class template for the right hand side
+ // of the derivative of the pressure potential
template <int dim>
-class RightHandSide_v : public Function<dim>
+class RightHandSideV : public Function<dim>
{
public:
- RightHandSide_v () : Function<dim>() {};
+ RightHandSideV () : Function<dim>() {};
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
- // Declare a class template for the initial values
- // of the pressure potential
+ // Declare a class template for the initial values
+ // of the pressure potential
template <int dim>
-class InitialValues_p : public Function<dim>
+class InitialValuesP : public Function<dim>
{
public:
- InitialValues_p () : Function<dim>() {};
+ InitialValuesP () : Function<dim>() {};
- virtual double value (const Point<dim> &p,
+ virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
- // Declare a class template for the initial values
- // of the derivative of the pressure potential
+ // Declare a class template for the initial values
+ // of the derivative of the pressure potential
template <int dim>
-class InitialValues_v : public Function<dim>
+class InitialValuesV : public Function<dim>
{
public:
- InitialValues_v () : Function<dim>() {};
+ InitialValuesV () : Function<dim>() {};
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
- // Here is the function to set the right hand side
- // values to be zero for pressure potential
+ // Here is the function to set the right hand side
+ // values to be zero for pressure potential
template <int dim>
-double RightHandSide_p<dim>::value (const Point<dim> &/*p*/,
- const unsigned int /*component*/) const
+double RightHandSideP<dim>::value (const Point<dim> &/*p*/,
+ const unsigned int /*component*/) const
{
return 0;
}
- // Similarly we set the right-hand size of the
- // derivative of the pressure potential to be
- // zero
+ // Similarly we set the right-hand size of the
+ // derivative of the pressure potential to be
+ // zero
template <int dim>
-double RightHandSide_v<dim>::value (const Point<dim> &/*p*/,
- const unsigned int /*component*/) const
+double RightHandSideV<dim>::value (const Point<dim> &/*p*/,
+ const unsigned int /*component*/) const
{
return 0;
}
- // The sources of the thermoacoustic waves
- // are small absorbers. We will compare the
- // simulation results with the experimental
- // data.
+ // The sources of the thermoacoustic waves
+ // are small absorbers. We will compare the
+ // simulation results with the experimental
+ // data.
template <int dim>
-double InitialValues_p<dim>::value (const Point<dim> &p,
- const unsigned int /*component*/) const
+double InitialValuesP<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
{
if (std::sqrt(p.square())< 0.025 )
- return 1;
- // The "distance" function is used to compute
- // the Euclidian distance between two points.
+ return 1;
+ // The "distance" function is used to compute
+ // the Euclidian distance between two points.
if (p.distance(Point<dim>(-0.135,0))<0.05)
- return 1;
+ return 1;
if (p.distance(Point<dim>(0.17,0))<0.03)
- return 1;
+ return 1;
if (p.distance(Point<dim>(-0.25,0))<0.02)
- return 1;
+ return 1;
if (p.distance(Point<dim>(-0.05,-0.15))<0.015)
- return 1;
+ return 1;
- return 0;
+ return 0;
}
- // Initial value for the derivative of
- // pressure potential is set to zero
+ // Initial value for the derivative of
+ // pressure potential is set to zero
template <int dim>
-double InitialValues_v<dim>::value (const Point<dim> &/*p*/,
- const unsigned int /*component*/) const
+double InitialValuesV<dim>::value (const Point<dim> &/*p*/,
+ const unsigned int /*component*/) const
{
- return 0;
-}
-
- // Evaluate point values at arbitrary locations
- // In real situation, we collect data by placing
- // a detector in the medium. By scanning the detector,
- // we obtain a series projections of the target
- // from different viewing angles. By using a
- // circular radon transform, we can reconstruct
- // the energy distribution in the target area from
- // the measurements obtained by the detectors.
-
-template <int dim>
-double point_value (const DoFHandler<dim> &dof,
- const Vector<double> &fe_function,
- const Point<dim> &point)
-{
- // Define a map that maps the unit cell to a
- // a general grid cell with straight lines in
- // dim dimensions
- static const MappingQ1<dim> mapping;
- const FiniteElement<dim>& fe = dof.get_fe();
-
- Assert(fe.n_components() == 1,
- ExcMessage ("Finite element is not scalar as is necessary for this function"));
-
- // First find the cell in which this point
- // is, initialize a quadrature rule with
- // it, and then a FEValues object
- // The algorithm will first look for the
- // surrounding cell on a coarse grid, and
- // then recersively checking its sibling
- // cells.
- const typename DoFHandler<dim>::active_cell_iterator cell = GridTools::find_active_cell_around_point (dof, point);
-
- const Point<dim> unit_point = mapping.transform_real_to_unit_cell(cell, point);
- Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
- ExcInternalError());
-
- const Quadrature<dim> quadrature (unit_point);
- FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
- fe_values.reinit(cell);
-
- // Then use this to get at the values of
- // the given fe_function at this point
- std::vector<double> u_value(1);
- fe_values.get_function_values(fe_function, u_value);
-
- return u_value[0];
+ return 0;
}
- // @sect4{Initialize the problem}
- // Acoustic_speed here is the acoustic speed
- // in the medium. Specifically we use acoustic speed
- // in mineral oil. We use Crank-Nicolson scheme
- // for our time-dependent problem, therefore theta is
- // set to be 0.5. The step size of the detector
- // is 2.25 degree, which means we need 160 steps
- // in order to finish a circular scan. The radius of the
- // scanning circle is select to be half way between
- // the center and the boundary to avoid the reflections
- // from the the boundary, so as to miminize the
- // interference brought by the inperfect absorbing
- // boundary condition. The time step is selected
- // to satisfy $k = h/c$
+ // @sect4{Initialize the problem}
+ // Acoustic_speed here is the acoustic speed
+ // in the medium. Specifically we use acoustic speed
+ // in mineral oil. We use Crank-Nicolson scheme
+ // for our time-dependent problem, therefore theta is
+ // set to be 0.5. The step size of the detector
+ // is 2.25 degree, which means we need 160 steps
+ // in order to finish a circular scan. The radius of the
+ // scanning circle is select to be half way between
+ // the center and the boundary to avoid the reflections
+ // from the the boundary, so as to miminize the
+ // interference brought by the inperfect absorbing
+ // boundary condition. The time step is selected
+ // to satisfy $k = h/c$
template <int dim>
TATForwardProblem<dim>::TATForwardProblem () :
fe (1),
n_refinements (7),
acoustic_speed (1.437),
theta (0.5),
- total_time (0.7),
+ end_time (0.7),
time_step (0.5/std::pow(2.,1.0*n_refinements)/acoustic_speed),
step_angle (2.25),
radius (0.5)
{}
- // This is similar to step-6 except that
- // the mesh generated is a hyper_ball. We select
- // hyper_ball instead of hyper_cube because of
- // our data collection geometry is on a circular in
- // 2-D, and on a sphere in 3-D.
+
+
+ // @sect4{TATForwardProblem::setup_system}
+
+ // The following system is pretty much what
+ // we've already done in @ref step_23
+ // "step-23", but with two important
+ // differences. First, we have to create a
+ // circular (or spherical) mesh around the
+ // origin, with a radius of 1. This nothing
+ // new: we've done so before in @ref step_6
+ // "step-6", @ref step_10 "step-10", and @ref
+ // step_11 "step-11", where we also explain
+ // how to attach a boundary object to a
+ // triangulation to be used whenever the
+ // triangulation needs to know where new
+ // boundary points lie when a cell is
+ // refined. Following this, the mesh is
+ // refined <code>n_refinements</code> times
+ // — this variable was introduced to
+ // make sure the time step size is always
+ // compatible with the cell size, and
+ // therefore satisfies the CFL condition that
+ // was talked about in the introduction of
+ // @ref step_23 "step-23".
+ //
+ // The only other significant change is that
+ // we need to build the boundary mass
+ // matrix. We will comment on this further
+ // down below.
template <int dim>
-void TATForwardProblem<dim>::make_grid_and_dofs ()
+void TATForwardProblem<dim>::setup_system ()
{
- // In two dimensional domain. The center of the
- // circle shall be the point (0,0) and the radius
- // is 1
- const Point<2> center (0,0);
- GridGenerator::hyper_ball (triangulation, center, 1);
- // Accordingly, we use hyper ball boundary
- // instead of hyper cube.
+ GridGenerator::hyper_ball (triangulation, Point<dim>(), 1.);
static const HyperBallBoundary<dim> boundary_description(center);
triangulation.set_boundary (0,boundary_description);
- // The mesh is refined n_refinements times
triangulation.refine_global (n_refinements);
- std::cout << " Number of active cells: "
+ std::cout << "Number of active cells: "
<< triangulation.n_active_cells()
- << std::endl
- << " Total number of cells: "
- << triangulation.n_cells()
- << std::endl;
+ << std::endl;
dof_handler.distribute_dofs (fe);
- std::cout << " Number of degrees of freedom: "
+ std::cout << "Number of degrees of freedom: "
<< dof_handler.n_dofs()
+ << std::endl
<< std::endl;
-
sparsity_pattern.reinit (dof_handler.n_dofs(),
dof_handler.n_dofs(),
dof_handler.max_couplings_between_dofs());
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
sparsity_pattern.compress();
- // We will do the following for both
- // the pressure potential and its derivative
- system_matrix_p.reinit (sparsity_pattern);
- system_matrix_v.reinit (sparsity_pattern);
+ system_matrix.reinit (sparsity_pattern);
mass_matrix.reinit (sparsity_pattern);
laplace_matrix.reinit (sparsity_pattern);
- boundary_matrix.reinit (sparsity_pattern);
- solution_p.reinit (dof_handler.n_dofs());
- old_solution_p.reinit (dof_handler.n_dofs());
- system_rhs_p.reinit (dof_handler.n_dofs());
-
- solution_v.reinit (dof_handler.n_dofs());
- old_solution_v.reinit (dof_handler.n_dofs());
- system_rhs_v.reinit (dof_handler.n_dofs());
-
-}
-
-
- // @sect3{ Assemble system}
- // Because we used absorbing boundary condition in the
- // simulation, a new boundary matrix is introduced.
- // We need to assemble boundary matrix. The detailed
- // description for assembling matrix is discussed in
- // step-3.
-template <int dim>
-void TATForwardProblem<dim>::assemble_system ()
-{
MatrixCreator::create_mass_matrix (dof_handler, QGauss<dim>(3),
mass_matrix);
- MatrixCreator::create_mass_matrix (dof_handler, QTrapez<dim>(),
- mass_matrix);
- mass_matrix /= 2;
-
MatrixCreator::create_laplace_matrix (dof_handler, QGauss<dim>(3),
laplace_matrix);
- MatrixCreator::create_laplace_matrix (dof_handler, QTrapez<dim>(),
- laplace_matrix);
- laplace_matrix /= 2;
- const QGauss<dim-1> quadrature_formula(3);
+ // The second difference, as mentioned, to
+ // @ref step_23 "step-23" is that we need
+ // to build the boundary mass matrix that
+ // grew out of the absorbing boundary
+ // conditions.
+ //
+ // A first observation would be that this
+ // matrix is much sparser than the regular
+ // mass matrix, since none of the shape
+ // functions with purely interior support
+ // contributes to this matrix. We could
+ // therefore optimize the storage pattern
+ // to this situation and build up a second
+ // sparsity pattern that only contains the
+ // nonzero entries that we need. There is a
+ // trade-off to make here: first, we would
+ // have to have a second sparsity pattern
+ // object, so that costs memory. Secondly,
+ // the matrix attached to this sparsity
+ // pattern is going to be smaller and
+ // therefore requires less memore; it would
+ // also be faster to perform matrix-vector
+ // multiplications with it. The final
+ // argument, however, is the one that tips
+ // the scale: we are not primarily
+ // interested in performing matrix-vector
+ // with the boundary matrix alone (though
+ // we need to do that for the right hand
+ // side vector once per time step), but
+ // mostly wish to add it up to the other
+ // matrices used in the first of the two
+ // equations since this is the one that is
+ // going to be multiplied with once per
+ // iteration of the CG method,
+ // i.e. significantly more often. It is now
+ // the case that the SparseMatrix::add
+ // class allows to add one matrix to
+ // another, but only if they use the same
+ // sparsity pattern (the reason being that
+ // we can't add nonzero entries to a matrix
+ // after the sparsity pattern has been
+ // created, so we simply require that the
+ // two matrices have the same sparsity
+ // pattern.
+ //
+ // So let's go with that:
+ boundary_matrix.reinit (sparsity_pattern);
- FEFaceValues<dim> fe_values (fe, quadrature_formula,
- update_values | update_JxW_values);
+ // The second thing to do is to actually
+ // build the matrix. Here, we need to
+ // integrate over faces of cells, so first
+ // we need a quadrature object that works
+ // on <code>dim-1</code> dimensional
+ // objects. Secondly, the FEFaceValues
+ // variant of FEValues that works on faces,
+ // as its name suggest. And finally, the
+ // other variables that are part of the
+ // assembly machinery. All of this we put
+ // between curly braces to limit the scope
+ // of these variables to where we actually
+ // need them.
+ //
+ // The actual act of assembling the matrix
+ // is then fairly straightforward: we loop
+ // over all cells, over all faces of each
+ // of these cells, and then do something
+ // only if that particular face is at the
+ // boundary of the domain. Like this:
+ {
+ const QGauss<dim-1> quadrature_formula(3);
+ FEFaceValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_JxW_values);
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
- const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ std::vector<unsigned int> 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 (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (cell->at_boundary(f))
- {
- cell_matrix = 0;
-
- fe_values.reinit (cell, f);
-
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->at_boundary(f))
+ {
+ cell_matrix = 0;
+
+ fe_values.reinit (cell, f);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (fe_values.shape_value(i,q_point) *
+ fe_values.shape_value(j,q_point) *
+ fe_values.JxW(q_point));
+
+ cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- cell_matrix(i,j) += (fe_values.shape_value(i,q_point) *
- fe_values.shape_value(j,q_point) *
- fe_values.JxW(q_point));
-
- cell->get_dof_indices (local_dof_indices);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- boundary_matrix.add (local_dof_indices[i],
- local_dof_indices[j],
- cell_matrix(i,j));
- }
-
- // The system matrix of pressure potential
- // is shown in the introduction
- system_matrix_p = 0;
- system_matrix_p.copy_from (mass_matrix);
- system_matrix_p.add (time_step*time_step*theta*theta*acoustic_speed*acoustic_speed, laplace_matrix);
- system_matrix_p.add (acoustic_speed*theta*time_step, boundary_matrix);
- // The system matrix of the derivative
- // of the pressure potential is same as
- // the mass matrix
- system_matrix_v = 0;
- system_matrix_v.copy_from (mass_matrix);
+ boundary_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+ }
+
+ }
+
+ system_matrix.copy_from (mass_matrix);
+ system_matrix.add (time_step * time_step * theta * theta *
+ acoustic_speed * acoustic_speed,
+ laplace_matrix);
+ system_matrix.add (acoustic_speed * theta * time_step, boundary_matrix);
+
+
+ solution_p.reinit (dof_handler.n_dofs());
+ old_solution_p.reinit (dof_handler.n_dofs());
+ system_rhs_p.reinit (dof_handler.n_dofs());
+ solution_v.reinit (dof_handler.n_dofs());
+ old_solution_v.reinit (dof_handler.n_dofs());
+ system_rhs_v.reinit (dof_handler.n_dofs());
+
+ constraints.close ();
}
- // We will solve two equations.
- // We first solve for pressure potential
- // at a time step, then solve the derivative
- // of the pressure potential.
+ // @sect4{TATForwardProblem::solve_p and TATForwardProblem::solve_v}
+
+ // The following two functions, solving the
+ // linear systems for the pressure and the
+ // velocity variable, are taken pretty much
+ // verbatim (with the exception of the change
+ // of name from $u$ to $p$ of the primary
+ // variable) from @ref step_23 "step-23":
template <int dim>
void TATForwardProblem<dim>::solve_p ()
{
- SolverControl solver_control (1000, 1e-10);
+ SolverControl solver_control (1000, 1e-8*system_rhs_p.l2_norm());
SolverCG<> cg (solver_control);
- cg.solve (system_matrix_p, solution_p, system_rhs_p,
+
+ cg.solve (system_matrix, solution_p, system_rhs_p,
PreconditionIdentity());
- std::cout << " " << solver_control.last_step()
- << " CG iterations needed to obtain convergence."
+ std::cout << " p-equation: " << solver_control.last_step()
+ << " CG iterations."
<< std::endl;
}
-
- // To solve the derivative of the pressure potential
+
+
template <int dim>
void TATForwardProblem<dim>::solve_v ()
{
- SolverControl solver_control (1000, 1e-10);
+ SolverControl solver_control (1000, 1e-8*system_rhs_v.l2_norm());
SolverCG<> cg (solver_control);
-
- cg.solve (system_matrix_v, solution_v, system_rhs_v,
+
+ cg.solve (mass_matrix, solution_v, system_rhs_v,
PreconditionIdentity());
- std::cout << " " << solver_control.last_step()
- << " CG iterations needed to obtain convergence."
+ std::cout << " v-equation: " << solver_control.last_step()
+ << " CG iterations."
<< std::endl;
}
- // We output the solution for pressure potential
- // at each time step in "vtk" format.
+
+
+ // @sect4{TATForwardProblem::output_results}
+
+ // The same holds here: the function is from
+ // @ref step_23 "step-23".
template <int dim>
-void TATForwardProblem<dim>::output_results (const unsigned int timestep_number) const
+void TATForwardProblem<dim>::output_results () const
{
-
- DataOut<dim> data_out;
-
+ DataOut<dim> data_out;
+
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution_p, "P");
data_out.add_data_vector (solution_v, "V");
std::ostringstream filename;
filename << "solution-"
- << timestep_number<<".vtk";
+ << Utilities::int_to_string (timestep_number, 3)
+ << ".gnuplot";
std::ofstream output (filename.str().c_str());
- data_out.write_vtk (output);
-
-
+ data_out.write_gnuplot (output);
}
+
+//XXX
// This is the main function
//
template <int dim>
void TATForwardProblem<dim>::run ()
{
- std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
-
- make_grid_and_dofs();
- assemble_system ();
-
- ConstraintMatrix constraints;
- constraints.close();
+ setup_system();
VectorTools::project (dof_handler,constraints,
- QGauss<dim>(3), InitialValues_p<dim>(),
- old_solution_p);
+ QGauss<dim>(3), InitialValuesP<dim>(),
+ old_solution_p);
VectorTools::project (dof_handler,constraints,
- QGauss<dim>(3), InitialValues_v<dim>(),
- old_solution_v);
+ QGauss<dim>(3), InitialValuesV<dim>(),
+ old_solution_v);
- unsigned int timestep_number = 1;
+ timestep_number = 1;
unsigned int n_steps;
unsigned int n_detectors;
double scanning_angle;
- // Number of time steps is defined as the
- // ratio of the total time to the time step
- n_steps=static_cast<unsigned int>(std::floor(total_time/time_step));
- // Number of detector positions is defined
- // as the ratio of 360 degrees to the step
- // angle
+ // Number of time steps is defined as the
+ // ratio of the total time to the time step
+ n_steps=static_cast<unsigned int>(std::floor(end_time/time_step));
+ // Number of detector positions is defined
+ // as the ratio of 360 degrees to the step
+ // angle
n_detectors=static_cast<unsigned int>(std::ceil(360/step_angle));
- // Define two vectors to hold the coordinates
- // of the detectors in the scanning
- // geometry
+ // Define two vectors to hold the coordinates
+ // of the detectors in the scanning
+ // geometry
Vector<double> detector_x (n_detectors+1);
Vector<double> detector_y (n_detectors+1);
- // Define a vector to hold the value obtained
- // by the detector
+ // Define a vector to hold the value obtained
+ // by the detector
Vector<double> project_dat (n_steps * n_detectors +1);
- // Get the coordinates of the detector on the
- // different locations of the circle.
- // Scanning angle is viewing angle at
- // current position. The coordinates of
- // the detectors are calculated from the radius
- // and scanning angle.
+ // Get the coordinates of the detector on the
+ // different locations of the circle.
+ // Scanning angle is viewing angle at
+ // current position. The coordinates of
+ // the detectors are calculated from the radius
+ // and scanning angle.
scanning_angle=0;
for (unsigned int i=1; i<=n_detectors; i++){
- // Scanning clockwisely. We need to change the angles
- // into radians because std::cos and std:sin accept
- // values in radian only
+ // Scanning clockwisely. We need to change the angles
+ // into radians because std::cos and std:sin accept
+ // values in radian only
scanning_angle -= step_angle/180 * 3.14159265;
detector_x(i) = radius * std::cos(scanning_angle);
detector_y(i) = radius * std::sin(scanning_angle);
std::cout<< "Total number of time steps = "<< n_steps <<std::endl;
std::cout<< "Total number of detectors = "<< n_detectors << std::endl;
- // Open a file to write the data
- // obtained by the detectors
+ // Open a file to write the data
+ // obtained by the detectors
std::ofstream proj_out;
proj_out.open("proj.dat");
- for (double time = time_step; time<=total_time; time+=time_step, ++timestep_number)
+ for (double time = time_step; time<=end_time; time+=time_step, ++timestep_number)
{
std::cout << std::endl;
std::cout<< "time_step " << timestep_number << " @ t=" << time << std::endl;
Vector<double> tmp2 (solution_v.size());
Vector<double> F1 (solution_p.size());
Vector<double> F2 (solution_v.size());
- // Calculate G1 as defined in the introduction section
+ // Calculate G1 as defined in the introduction section
mass_matrix.vmult (tmp1, old_solution_p);
mass_matrix.vmult (tmp2, old_solution_v);
F1 = tmp1;
F1.add(time_step * (1-theta), tmp2);
- // Calculate G2 as defined in the introduction section
+ // Calculate G2 as defined in the introduction section
mass_matrix.vmult (tmp1, old_solution_v);
laplace_matrix.vmult (tmp2, old_solution_p);
F2 = tmp1;
boundary_matrix.vmult (tmp1,old_solution_p);
F2.add(acoustic_speed,tmp1);
- // Compute the pressure potential p, the formula
- // has been presented in the introduction section
+ // Compute the pressure potential p, the formula
+ // has been presented in the introduction section
system_rhs_p = F1;
system_rhs_p.add(time_step * theta , F2);
- RightHandSide_p<dim> rhs_function_p;
+ RightHandSideP<dim> rhs_function_p;
rhs_function_p.set_time (time);
tmp1=0;
solve_p ();
- // Compute the derivative potential pressure.
- // The formula has been presented in the introduction
- // section. The potential derivative is calculated
- // after the potential pressure because the calculation
- // depends on the current value of the potential
- // pressure
+ // Compute the derivative potential pressure.
+ // The formula has been presented in the introduction
+ // section. The potential derivative is calculated
+ // after the potential pressure because the calculation
+ // depends on the current value of the potential
+ // pressure
system_rhs_v = F2;
tmp1 = 0;
boundary_matrix.vmult(tmp1, solution_p);
system_rhs_v.add(-acoustic_speed,tmp1);
- RightHandSide_v<dim> rhs_function_v;
+ RightHandSideV<dim> rhs_function_v;
rhs_function_v.set_time (time);
tmp2 = 0;
system_rhs_p.add(-(1-theta)*time_step*acoustic_speed*acoustic_speed,tmp2);
solve_v ();
- // Compute the energy in the system.By checking
- // energy change in the system, we can verify
- // the correctness of the code.
+ // Compute the energy in the system.By checking
+ // energy change in the system, we can verify
+ // the correctness of the code.
double energy = (mass_matrix.matrix_scalar_product(solution_v,solution_v)+
- acoustic_speed*acoustic_speed*laplace_matrix.matrix_scalar_product(solution_p,solution_p))/2;
+ acoustic_speed*acoustic_speed*laplace_matrix.matrix_scalar_product(solution_p,solution_p))/2;
std::cout << "energy= " << energy << std::endl;
- // output_results (timestep_number);
+ output_results ();
- // Evaluate the value at specific locations.
- // For 2-D, it is on a circle. For 1-D,
- // it is a point detector.
+ // Evaluate the value at specific locations.
+ // For 2-D, it is on a circle. For 1-D,
+ // it is a point detector.
proj_out << time ;
for (unsigned i=1 ; i<=n_detectors; i++){
- project_dat((timestep_number-1)*n_detectors+i)=point_value (dof_handler,solution_p,
- Point<2>(detector_x(i),detector_y(i)));
+ project_dat((timestep_number-1)*n_detectors+i)
+ = VectorTools::point_value (dof_handler,solution_p,
+ Point<2>(detector_x(i),detector_y(i)));
proj_out << " "<< project_dat((timestep_number-1)*n_detectors+i)<<" " ;
}
proj_out<<std::endl;
- // Update the values for the pressure potential
- // and its derivative.
+ // Update the values for the pressure potential
+ // and its derivative.
old_solution_p = solution_p;
solution_p = 0;
old_solution_v = solution_v;
solution_v = 0;
}
- proj_out.close();
-
-
}
- // @sect3{The "main" function}
- // The main function calls the above functions
- // in the order of their appearances.
-int main ()
+
+
+ // @sect3{The <code>main</code> function}
+
+ // What remains is the main function of the
+ // program. There is nothing here that hasn't
+ // been shown in several of the previous
+ // programs:
+int main ()
{
- deallog.depth_console (0);
- {
- TATForwardProblem<2> TAT_forward_2d;
- TAT_forward_2d.run ();
- }
+ try
+ {
+ deallog.depth_console (0);
+ TATForwardProblem<2> forward_problem_solver;
+ forward_problem_solver.run ();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
return 0;
}