]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
More text
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 9 Sep 2006 21:25:45 +0000 (21:25 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 9 Sep 2006 21:25:45 +0000 (21:25 +0000)
git-svn-id: https://svn.dealii.org/trunk@13874 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-24/doc/intro.dox
deal.II/examples/step-24/step-24.cc

index 939112b4fa085beb25d142947174f0e897ecb145..96dd500d918d1d2685b3417289fbec6ee89f7299 100644 (file)
 <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&amp;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}    \\
@@ -136,7 +162,7 @@ c_0^2\,\delta t\,\theta\,A+c_0\,B  &  M   \\
               \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]
 
@@ -147,16 +173,16 @@ G_1 \\
 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}
index fed5cdf9e2c96dba16095da9f28dc669283c24ab..ec5ac84e3c105029d723807097ecf404f6de4413 100644 (file)
@@ -1,73 +1,62 @@
-
-/*    $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
 {
@@ -76,223 +65,168 @@ 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),
@@ -300,188 +234,254 @@ TATForwardProblem<dim>::TATForwardProblem () :
                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
+                                // &mdash; 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");
@@ -490,66 +490,61 @@ void TATForwardProblem<dim>::output_results (const unsigned int timestep_number)
 
   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);
@@ -558,14 +553,14 @@ void TATForwardProblem<dim>::run ()
   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;
@@ -574,13 +569,13 @@ void TATForwardProblem<dim>::run ()
       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;
@@ -589,13 +584,13 @@ void TATForwardProblem<dim>::run ()
       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;
@@ -613,12 +608,12 @@ void TATForwardProblem<dim>::run ()
 
       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;
@@ -628,7 +623,7 @@ void TATForwardProblem<dim>::run ()
       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;
@@ -644,54 +639,82 @@ void TATForwardProblem<dim>::run ()
       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;
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.