From: bangerth Date: Sat, 9 Sep 2006 21:25:45 +0000 (+0000) Subject: More text X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b218b3e6f47dd898c4a4200129fff0d93d4d25c5;p=dealii-svn.git More text git-svn-id: https://svn.dealii.org/trunk@13874 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-24/doc/intro.dox b/deal.II/examples/step-24/doc/intro.dox index 939112b4fa..96dd500d91 100644 --- a/deal.II/examples/step-24/doc/intro.dox +++ b/deal.II/examples/step-24/doc/intro.dox @@ -1,134 +1,160 @@

Introduction

-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. - - -

Problem

- -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. + + +

The problem

+ +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. +

Weak form and discretization

-

Weak form and Discretization

- -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} diff --git a/deal.II/examples/step-24/step-24.cc b/deal.II/examples/step-24/step-24.cc index fed5cdf9e2..ec5ac84e3c 100644 --- a/deal.II/examples/step-24/step-24.cc +++ b/deal.II/examples/step-24/step-24.cc @@ -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 -#include - // 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 -#include - // 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 -#include -#include -#include -#include -#include -#include - -#include + // @sect3{Include files} + + // The following have all been covered + // previously: #include #include -#include -#include +#include +#include + #include #include #include #include #include + +#include +#include +#include +#include +#include + +#include +#include +#include #include +#include +#include + +#include #include #include -#include - // These are for c++ #include #include #include - -#include - -#include - // @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 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 triangulation; FE_Q fe; DoFHandler dof_handler; + ConstraintMatrix constraints; + SparsityPattern sparsity_pattern; - SparseMatrix system_matrix_p; - SparseMatrix system_matrix_v; + SparseMatrix system_matrix; SparseMatrix mass_matrix; SparseMatrix laplace_matrix; + + Vector solution_p, solution_v; + Vector old_solution_p, old_solution_v; + Vector system_rhs_p, system_rhs_v; + + double time, time_step; + unsigned int timestep_number; + const double theta; + + // SparseMatrix 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 solution_p; - Vector old_solution_p; - Vector system_rhs_p; - - Vector solution_v; - Vector old_solution_v; - Vector 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 -class RightHandSide_p : public Function +class RightHandSideP : public Function { public: - RightHandSide_p () : Function() {}; + RightHandSideP () : Function() {}; virtual double value (const Point &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 -class RightHandSide_v : public Function +class RightHandSideV : public Function { public: - RightHandSide_v () : Function() {}; + RightHandSideV () : Function() {}; virtual double value (const Point &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 -class InitialValues_p : public Function +class InitialValuesP : public Function { public: - InitialValues_p () : Function() {}; + InitialValuesP () : Function() {}; - virtual double value (const Point &p, + virtual double value (const Point &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 -class InitialValues_v : public Function +class InitialValuesV : public Function { public: - InitialValues_v () : Function() {}; + InitialValuesV () : Function() {}; virtual double value (const Point &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 -double RightHandSide_p::value (const Point &/*p*/, - const unsigned int /*component*/) const +double RightHandSideP::value (const Point &/*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 -double RightHandSide_v::value (const Point &/*p*/, - const unsigned int /*component*/) const +double RightHandSideV::value (const Point &/*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 -double InitialValues_p::value (const Point &p, - const unsigned int /*component*/) const +double InitialValuesP::value (const Point &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(-0.135,0))<0.05) - return 1; + return 1; if (p.distance(Point(0.17,0))<0.03) - return 1; + return 1; if (p.distance(Point(-0.25,0))<0.02) - return 1; + return 1; if (p.distance(Point(-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 -double InitialValues_v::value (const Point &/*p*/, - const unsigned int /*component*/) const +double InitialValuesV::value (const Point &/*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 -double point_value (const DoFHandler &dof, - const Vector &fe_function, - const Point &point) -{ - // Define a map that maps the unit cell to a - // a general grid cell with straight lines in - // dim dimensions - static const MappingQ1 mapping; - const FiniteElement& 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::active_cell_iterator cell = GridTools::find_active_cell_around_point (dof, point); - - const Point unit_point = mapping.transform_real_to_unit_cell(cell, point); - Assert (GeometryInfo::is_inside_unit_cell (unit_point), - ExcInternalError()); - - const Quadrature quadrature (unit_point); - FEValues 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 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 TATForwardProblem::TATForwardProblem () : fe (1), @@ -300,188 +234,254 @@ TATForwardProblem::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 n_refinements 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 -void TATForwardProblem::make_grid_and_dofs () +void TATForwardProblem::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(), 1.); static const HyperBallBoundary 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 -void TATForwardProblem::assemble_system () -{ MatrixCreator::create_mass_matrix (dof_handler, QGauss(3), mass_matrix); - MatrixCreator::create_mass_matrix (dof_handler, QTrapez(), - mass_matrix); - mass_matrix /= 2; - MatrixCreator::create_laplace_matrix (dof_handler, QGauss(3), laplace_matrix); - MatrixCreator::create_laplace_matrix (dof_handler, QTrapez(), - laplace_matrix); - laplace_matrix /= 2; - const QGauss 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 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 dim-1 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 quadrature_formula(3); + FEFaceValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - for (unsigned int f=0; f::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::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->at_boundary(f)) + { + cell_matrix = 0; + + fe_values.reinit (cell, f); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); for (unsigned int i=0; iget_dof_indices (local_dof_indices); - for (unsigned int i=0; i void TATForwardProblem::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 void TATForwardProblem::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 -void TATForwardProblem::output_results (const unsigned int timestep_number) const +void TATForwardProblem::output_results () const { - - DataOut data_out; - + DataOut 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::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 void TATForwardProblem::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(3), InitialValues_p(), - old_solution_p); + QGauss(3), InitialValuesP(), + old_solution_p); VectorTools::project (dof_handler,constraints, - QGauss(3), InitialValues_v(), - old_solution_v); + QGauss(3), InitialValuesV(), + 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(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(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(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 detector_x (n_detectors+1); Vector 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 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::run () std::cout<< "Total number of time steps = "<< n_steps <main 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; }