+/* */
+/* Copyright (C) 2007 by the deal.II authors and M. Allmaras */
+/* */
+/* 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}
// The following header files are unchanged
#include <fstream>
- // This header file is needed for
- // the ParameterHandler class that we will
- // use to read parameters from a
- // configuration file during runtime:
+ // This header file contains the necessary
+ // declarations for the ParameterHandler class
+ // that we will use to read our parameters from
+ // a configuration file:
#include <base/parameter_handler.h>
+ // For solving the linear system, we'll use
+ // the sparse LU-decomposition provided by UMFPACK,
+ // for which the following header file is needed.
+ // Note that in order to compile this tutorial program,
+ // the deal.II-library needs to be
+ // built with UMFPACK support, which can be most easily
+ // achieved by giving the <code> --with-umfpack</code>
+ // switch when configuring the library:
#include <lac/sparse_direct.h>
+
+ // The FESystem class allows us to stack
+ // several FE-objects to one compound, vector-valued finite
+ // element field. The necessary declarations for this class
+ // are provided in this header file:
#include <fe/fe_system.h>
+ // The last header is from the C++ standard library and provides
+ // functions that will allow us to measure execution time
+ // of various parts of our program:
+#include <ctime>
+
+ // Although we'll follow good deal.ii practice and keep
+ // all of the code dimension independent, we will
+ // really only consider the 2D problem here:
#define DIM 2
using namespace dealii;
+ // The next line provides a shorthand to the CLOCKS_PER_SEC
+ // constant, which is defined in the <code>ctime</code> header
+ // and contains the number of processor ticks in each second.
+ // Henceforth we can compute actual time from ticks by simply dividing
+ // by tps.
+const double tps = CLOCKS_PER_SEC;
+
+ // @sect3{The <code>DirichletBoundaryValues</code> class}
+
+ // First we define a class for the function representing
+ // the Dirichlet boundary values. This has been done many times before
+ // and therefore does not need much explanation.
template <int dim>
class DirichletBoundaryValues : public Function<dim>
{
};
+ // Since there are two values $v$ and $w$ that need to be prescribed
+ // at the boundary, the boundary value function must return a vector
+ // with two entries. In our case the function is very simple,
+ // it just returns 1 for the real part $v$ and 0 for the imaginary
+ // part $w$ regardless of the point where it is evaluated.
template <int dim>
inline
void DirichletBoundaryValues<dim>::vector_value (const Point<dim> &/*p*/,
DirichletBoundaryValues<dim>::vector_value (points[p], value_list[p]);
}
+ // @sect3{The <code>ParameterReader</code> class}
+ // The next classis responsible for preparing the
+ // ParameterHandler object and reading parameters from
+ // an input file.
+ // It includes a function <code>declare_parameters</code>
+ // that declares all the necessary parameters
+ // and a <code>read_parameters</code>
+ // function that is called from outside to initiate
+ // the parameter reading process.
class ParameterReader : public Subscriptor
{
public:
ParameterReader(ParameterHandler &);
- void read_parameters();
+ void read_parameters(const std::string);
private:
void declare_parameters();
ParameterHandler &prm;
};
-
+ // The constructor stores a reference to
+ // the ParameterHandler object that is passed to it:
ParameterReader::ParameterReader(ParameterHandler ¶mhandler)
:
prm(paramhandler)
{}
+ // @sect4{<code>ParameterReader::declare_parameters</code>}
-void ParameterReader::read_parameters()
-{
- declare_parameters();
-
- const std::string parameter_file = "step-29.prm";
- prm.read_input (parameter_file);
-}
-
-
+ // The <code>declare_parameters</code> function declares all the parameters
+ // that our ParameterHandler object will discover in the input file,
+ // along with their types, range conditions and the subsections
+ // they appear in:
void ParameterReader::declare_parameters()
{
+ // Parameters for mesh and geometry include the number
+ // of global refinement steps that are applied to the initial
+ // coarse mesh and the focal distance $d$ of the transducer lens. For the number
+ // of refinement steps, we allow integer values between 1 and 10,
+ // and for the focal distance any number greater than zero:
prm.enter_subsection ("Mesh & geometry parameters");
prm.declare_entry("Number of refinements", "6",
prm.declare_entry("Focal distance", "0.3",
Patterns::Double(0),
"Distance of the focal point of the lens "
- "to the x-axis (or xy-plane in 3D)");
+ "to the x-axis");
prm.leave_subsection ();
+ // The next subsection is devoted to the physical parameters appearing
+ // in the equation, which are the frequency $\omega$
+ // and wave speed $c$:
prm.enter_subsection ("Physical constants");
prm.declare_entry("c", "1.5e5",
- Patterns::Double(),
+ Patterns::Double(0),
"Wave speed");
prm.declare_entry("omega", "5.0e7",
- Patterns::Double(),
+ Patterns::Double(0),
"Frequency");
prm.leave_subsection ();
+
+ // Last but not least we would like to be able to change
+ // some properties of the output, like filename and format,
+ // through entries in the configuration file, which is the
+ // purpose of the last subsection:
prm.enter_subsection ("Output parameters");
prm.declare_entry("Output file", "solution",
Patterns::Anything(),
"Name of the output file (without extension)");
+ // Since different output formats may require different
+ // parameters for generating output (like for example,
+ // postscript output needs viewpoint angles, line widths, colors
+ // etc), it would be cumbersome if we had to declare all these parameters
+ // by hand for every possible output format supported in the library. Instead,
+ // each output format has a <code>FormatFlags::declare_parameters</code>
+ // function, which declares all the parameters specific to that format in
+ // an own subsection. The following call of
+ // DataOutInterface<1>::declare_parameters executes
+ // <code>declare_parameters</code> for all available output formats, so that
+ // for each format an own subsection will be created with parameters declared
+ // for that particular output format.
+ // To find out what parameters there are for which output format, you can either
+ // consult the documentation of the DataOutBase class, or simply run this
+ // program without a parameter file present. It will then create a file with all
+ // declared parameters set to their default values, which can conveniently serve
+ // as a starting point for setting the parameters to the values you desire.
DataOutInterface<1>::declare_parameters (prm);
prm.leave_subsection ();
}
+ // @sect4{<code>ParameterReader::read_parameters</code>}
+
+ // This is the main function in the ParameterReader class.
+ // It gets called from outside and first initiates declaration of
+ // the parameters, and then tries to read them from the input file whose
+ // filename is provided by the caller.
+void ParameterReader::read_parameters(const std::string parameter_file)
+{
+ declare_parameters();
+
+ prm.read_input (parameter_file);
+}
+
+
+ // @sect3{The <code>ComputeIntensity</code> class}
+
+ // As mentioned in the introduction, the quantitiy that we
+ // are really after is the spatial distribution of
+ // the intensity of the ultrasound wave, which corresponds
+ // to $|u|=\sqrt{v^2+w^2}$. Now we could just be content with
+ // having $v$ and $w$ in our output, and use a suitable
+ // visualization or postprocessing tool to derive $|u|$ from the
+ // solution we computed. However, there is also a way to output
+ // data derived from the solution in deal.II, and we are going
+ // to make use of this mechanism here.
+
+ // So far we have always used the DataOut::add_data_vector function
+ // to add vectors contaning output data to a DataOut object.
+ // There is a special version of this function
+ // that in addition to the data vector has an additional argument of
+ // type DataPostprocessor. What happens when this function
+ // is used for output is that at each point where output data
+ // is to be generated, the compute_derived_quantities function
+ // of the specified DataPostprocessor object is invoked to compute
+ // the output quantities from the values, the gradients and the
+ // second derivatives of the finite element function represented
+ // by the data vector (in the case of face related data, normal vectors
+ // are available as well). Hence, this allows us to output any quantity
+ // that can locally be derived from the values of the solution and
+ // its derivatives.
+ // Of course, the ultrasound intensity $|u|$ is such a quantity and
+ // its computation doesn't even involve any derivatives of $v$ or $w$.
+
+ // In practice, the DataPostprocessor class only provides an
+ // interface to this functionality, and we need to derive our own
+ // class from it in order to
+ // implement the functions specified by the interface.
+ // This is what the <code>ComputeIntensity</code> class is about.
+ // Notice that all its member functions are implementations of
+ // virtual functions defined by the interface class DataPostprocessor.
template <int dim>
-class Postprocessor : public DataPostprocessor<dim>
+class ComputeIntensity : public DataPostprocessor<dim>
{
public:
const std::vector< Vector< double > > &,
const std::vector< std::vector< Tensor< 1, dim > > > &,
const std::vector< std::vector< Tensor< 2, dim > > > &,
- const std::vector< Point< dim > > &,
+ const std::vector< Point< dim > > &,
std::vector< Vector< double > > &
) const;
unsigned int n_output_variables () const;
};
-
+ // The <code>get_names</code> function returns a vector of strings
+ // representing the names we assign to the individual
+ // quantities that our postprocessor outputs. In our
+ // case, the postprocessor has only $|u|$ as an output, so we
+ // need to provide just that one name:
template <int dim>
std::vector<std::string>
-Postprocessor<dim>::get_names() const
+ComputeIntensity<dim>::get_names() const
{
std::vector<std::string> field_names;
-
- field_names.push_back("Re_u");
- field_names.push_back("Im_u");
field_names.push_back("Intensity");
return field_names;
}
-
+ // The next function returns a set of flags that indicate
+ // which data is needed by the postprocessor in order to
+ // compute the output quantities.
+ // This can be any subset of update_values,
+ // update_gradients and update_hessians
+ // (and, in the case of face data, also
+ // update_normal_vectors).
+ // Of course, computation of the derivatives requires additional
+ // resources, so only the flags for data that is really needed
+ // should be given here. In our case, only the function values
+ // of $v$ and $w$ are needed to compute $|u|$, so we're good
+ // with the update_values flag.
template <int dim>
UpdateFlags
-Postprocessor<dim>::get_needed_update_flags () const
+ComputeIntensity<dim>::get_needed_update_flags () const
{
return update_values;
}
-
+ // To allow the caller to find out how many derived quantities
+ // are returned by the postprocessor, the
+ // <code>n_output_variables</code> function is used. Since
+ // we compute only $|u|$, the correct value to return
+ // in our case is just 1:
template <int dim>
unsigned int
-Postprocessor<dim>::n_output_variables () const
+ComputeIntensity<dim>::n_output_variables () const
{
- return 3;
+ return 1;
}
+ // The actual prostprocessing happens in the following function.
+ // Its inputs are a vector representing point values of the function
+ // and some tensor objects representing derivatives (that we don't
+ // use here since $|u|$ is computed from just $v$ and $w$).
+ // The derived quantities are returned in the
+ // <code>computed_quantities</code> vector.
+ // Remember that this function may only use data for which the
+ // respective update flag is specified by
+ // <code>get_needed_update_flags</code>. For example, we may
+ // not use the derivatives here,
+ // since our implementation of <code>get_needed_update_flags</code>
+ // requests that only function values are provided.
template <int dim>
void
-Postprocessor<dim>::compute_derived_quantities_vector (
- const std::vector< Vector< double > > &uh,
- const std::vector< std::vector< Tensor< 1, dim > > > &/*duh*/,
- const std::vector< std::vector< Tensor< 2, dim > > > &/*dduh*/,
- const std::vector< Point< dim > > &/*normals*/,
- std::vector< Vector< double > > &computed_quantities
+ComputeIntensity<dim>::compute_derived_quantities_vector (
+ const std::vector< Vector< double > > & uh,
+ const std::vector< std::vector< Tensor< 1, dim > > > & /*duh*/,
+ const std::vector< std::vector< Tensor< 2, dim > > > & /*dduh*/,
+ const std::vector< Point< dim > > & /*normals*/,
+ std::vector< Vector< double > > & computed_quantities
) const
{
Assert(computed_quantities.size() == uh.size(),
ExcDimensionMismatch (computed_quantities.size(), uh.size()));
+ // The computation itself is straightforward: We iterate
+ // over each entry in the output vector and compute
+ // $|u|$ from the corresponding values of $v$ and $w$:
for (unsigned int i=0; i<computed_quantities.size(); i++)
{
- Assert(computed_quantities[i].size() == 3,
- ExcDimensionMismatch (computed_quantities[i].size(), 3));
+ Assert(computed_quantities[i].size() == 1,
+ ExcDimensionMismatch (computed_quantities[i].size(), 1));
Assert(uh[i].size() == 2, ExcDimensionMismatch (uh[i].size(), 2));
- computed_quantities[i](0) = uh[i](0);
- computed_quantities[i](1) = uh[i](1);
- computed_quantities[i](2) = sqrt(uh[i](0)*uh[i](0) + uh[i](1)*uh[i](1));
+ computed_quantities[i](0) = sqrt(uh[i](0)*uh[i](0) + uh[i](1)*uh[i](1));
}
}
+ // @sect3{The <code>UltrasoundProblem</code> class}
+
+ // Finally here is the main class of this program.
+ // It's member functions are very similar to the previous
+ // examples and the list of member variables does not contain
+ // any major surprises either.
+ // The ParameterHandler object that is passed
+ // to the constructor is stored as a reference to allow
+ // easy access to the parameters from all functions of the class.
+ // Since we are working with vector valued finite elements, the
+ // FE object we are using is of type FESystem.
template <int dim>
class UltrasoundProblem
{
void run ();
private:
- static double get_omega (ParameterHandler &);
- static double get_c (ParameterHandler &);
-
void make_grid ();
void setup_system ();
void assemble_system ();
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution, system_rhs;
-
- const double c, omega;
};
-template <int dim>
-double
-UltrasoundProblem<dim>::get_omega (ParameterHandler &prm)
-{
- prm.enter_subsection ("Physical constants");
- double omega_tmp = prm.get_double("omega");
- prm.leave_subsection ();
-
- return omega_tmp;
-}
-
-
-template <int dim>
-double
-UltrasoundProblem<dim>::get_c (ParameterHandler &prm)
-{
- prm.enter_subsection ("Physical constants");
- double c_tmp = prm.get_double("c");
- prm.leave_subsection ();
-
- return c_tmp;
-}
-
+ // The constructor takes the ParameterHandler object and stores
+ // it in a reference. It also initializes the DoF-Handler and
+ // the finite element system, which consists of two copies
+ // of the scalar Q1 field, one for $v$ and one for $w$:
template <int dim>
UltrasoundProblem<dim>::UltrasoundProblem (ParameterHandler& param)
:
prm(param),
dof_handler(triangulation),
- fe(FE_Q<dim>(1), 2),
- c(get_c(prm)),
- omega(get_omega(prm))
+ fe(FE_Q<dim>(1), 2)
{}
dof_handler.clear();
}
+ // @sect4{<code>UltrasoundProblem::make_grid</code>}
+ // Here we setup the grid for our domain.
+ // As mentioned in the exposition, the geometry is just a unit square
+ // with the part of the boundary that represents the transducer
+ // lens replaced by a sector of a circle.
template <int dim>
void UltrasoundProblem<dim>::make_grid ()
{
+ // First we generate some logging output
+ // and store the current number of ticks to be able to
+ // compute execution time when this function is done:
+ deallog << "Generating grid... ";
+ clock_t start = clock();
+
+ // Then we query the values for the focal distance of the
+ // transducer lens and the number of mesh refinement steps
+ // from our ParameterHandler object:
prm.enter_subsection ("Mesh & geometry parameters");
const double focal_distance = prm.get_double("Focal distance");
prm.leave_subsection ();
- GridGenerator::subdivided_hyper_cube (triangulation, 5, 0, 1);
-
+ // Next, two points are defined for position and focal point
+ // of the transducer lens, which is the center of the circle
+ // whose segment will form the transducer part of the boundary. We
+ // compute the radius of this circle in such a way that the
+ // segment fits in the interval [0.4,0.6] on the x-axis.
+ // Notice that this is the only point in the program where things
+ // are slightly different in 2D and 3D.
+ // Even though this tutorial only deals with the 2D case,
+ // the necessary additions to make this program functional
+ // in 3D are so minimal that we opt for including them:
const Point<dim> transducer = (dim == 2) ?
Point<dim> (0.5, 0.0) :
Point<dim> (0.5, 0.5, 0.0),
focal_point.distance(transducer)) +
((dim==2) ? 0.01 : 0.02));
+
+ // As initial coarse grid we take a simple unit square with 5 subdivisions
+ // in each direction. Then we step through all cells to find the
+ // faces where the transducer is to be located, which in fact is just
+ // the single edge from 0.4 to 0.6 on the x-axis. This is where we want
+ // the refinements to be made according to a circle shaped boundary,
+ // so we mark this edge with a different boundary indicator.
+ GridGenerator::subdivided_hyper_cube (triangulation, 5, 0, 1);
+
typename Triangulation<dim>::cell_iterator
cell = triangulation.begin (),
endc = triangulation.end();
cell->face(face)->set_boundary_indicator (1);
+ // For the circle part of the transducer lens, a hyper-ball object is used
+ // (which, of course, in 2D just represents a circle),
+ // with radius and center as computed above. Then we assign this boundary-object
+ // to the part of the boundary with boundary indicator 1:
const HyperBallBoundary<dim> boundary(focal_point, radius);
triangulation.set_boundary(1, boundary);
+ // Now the global refinement is executed. Cells near the transducer
+ // location will be automatically refined according to the
+ // circle shaped boundary of the transducer lens:
triangulation.refine_global (N_ref);
- deallog << " Number of active cells: "
- << triangulation.n_active_cells()
+ // The next line releases the triangulation's
+ // pointer to the boundary object that we just created, which
+ // is necessary since the boundary object will be destructed
+ // as we leave this function
+ // and we don't want the triangulation to keep a hanging pointer.
+ triangulation.set_boundary(1);
+
+ // Lastly, we generate some more logging output. By querying
+ // the present number of ticks again and comparing to
+ // what we had at the beginning of the function, we can
+ // calculate execution time by dividing by tps.
+ // Note that the resolution of the <code>clock()</code> function
+ // is implementation depended, and also the <code>clock_t</code> values
+ // it returns may overflow, so this way of measuring execution
+ // time should be taken with a grain of salt as it may not
+ // be very accurate and even completely wrong for longer timespans:
+ clock_t end = clock();
+ deallog << "done ("
+ << (end - start) / tps
+ << "s)"
<< std::endl;
- triangulation.set_boundary(1);
+ deallog << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl;
}
+ // @sect4{<code>UltrasoundProblem::setup_system</code>}
+ // Initialization of the system matrix, sparsity patterns
+ // and vectors are the same as in previous examples
+ // and therefore do not need further comment:
template <int dim>
void UltrasoundProblem<dim>::setup_system ()
{
- dof_handler.distribute_dofs (fe);
+ deallog << "Setting up system... ";
+ clock_t start = clock();
- deallog << " Number of degrees of freedom: "
- << dof_handler.n_dofs()
- << std::endl;
+ dof_handler.distribute_dofs (fe);
sparsity_pattern.reinit (dof_handler.n_dofs(),
dof_handler.n_dofs(),
system_matrix.reinit (sparsity_pattern);
system_rhs.reinit (dof_handler.n_dofs());
solution.reinit (dof_handler.n_dofs());
+
+ clock_t end = clock();
+ deallog << "done ("
+ << (end - start) / tps
+ << "s)"
+ << std::endl;
+
+ deallog << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
}
+ // @sect4{<code>UltrasoundProblem::assemble_system</code>}
+ // As before, this function takes care of assembling the
+ // system matrix and right hand side vector:
template <int dim>
void UltrasoundProblem<dim>::assemble_system ()
{
- const double om2 = omega * omega;
- const double c2 = c * c;
+ deallog << "Assembling system matrix... ";
+ clock_t start = clock();
+
+ // First we query wavespeed and frequency from the
+ // ParameterHandler object and store them in local variables,
+ // as they will be used frequently throughout this
+ // function.
+
+ prm.enter_subsection ("Physical constants");
+
+ const double omega = prm.get_double("omega"),
+ c = prm.get_double("c");
+
+ prm.leave_subsection ();
+ // As usual, for computing integrals ordinary Gauss quadrature
+ // rule is used. Since our bilinear form involves boundary integrals
+ // on $\Gamma_2$, we also need a quadrature rule for surface
+ // integration on the faces, which are dim-1 dimensional:
QGauss<dim> quadrature_formula(2);
QGauss<dim-1> face_quadrature_formula(2);
n_face_q_points = face_quadrature_formula.n_quadrature_points,
dofs_per_cell = fe.dofs_per_cell;
- FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
-
- std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
+ // The FEValues objects will evaluate the shape functions for us.
+ // For the part of the bilinear form that involves integration on
+ // $\Omega$, we'll need the values and gradients
+ // of the shape functions, and of course the quadrature weights.
+ // For the terms involving the boundary integrals, only shape function
+ // values and the quadrature weights are necessary.
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_JxW_values);
FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
update_values | update_JxW_values);
+ // As usual, the system matrix is assembled cell by cell,
+ // and we need a matrix for storing the local cell contributions
+ // as well as an index vector to transfer the cell contributions to the
+ // appropriate location in the global system matrix after.
+ FullMatrix<double> cell_matrix (dofs_per_cell, 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)
{
+
+ // On each cell, we first need to reset the local contribution
+ // matrix and request the FEValues object to compute the shape
+ // functions for the current cell:
cell_matrix = 0;
fe_values.reinit (cell);
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
+
+ // At this point, it is important to keep in mind that we are dealing with a
+ // finite element system with two components. Due
+ // to the way we constructed this FESystem, namely as the cartesian product of
+ // two scalar finite element fields, each shape function
+ // has only a single nonzero component (they are, in deal.II lingo, primitive).
+ // Hence, each shape function can be viewed as one of the $\phi$'s or $\psi$'s
+ // from the introduction, and similarly
+ // the corresponding degrees of freedom can be attributed to either $\alpha$ or $\beta$.
+ // As we iterate through all the degrees of freedom on the current cell however,
+ // they do not come in any particular order, and so we cannot decide right away
+ // whether the DoFs with index i and j belong to the real or imaginary part of our solution.
+ // But if you look at the form of the system matrix in the introduction, this disctinction
+ // is crucial since it will determine to which block in the system matrix the
+ // contribution of the current pair of DoFs will go and hence which quantity we need to
+ // compute from the given two shape functions.
+ // Fortunately, the FESystem object can provide us with this information, namely it
+ // has a function FESystem::system_to_component_index, that for each local DoF index
+ // returns a pair of integers of which the first indicates to which component of the
+ // system the DoF belongs. The second integer of the pair indicates
+ // which index the DoF has in the scalar base finite element field, but this information
+ // is not relevant here. If you want to know more about this function and the underlying
+ // scheme behind primitive vector valued elements, take a look at step-8,
+ // where these topics are explained in depth.
if (fe.system_to_component_index(i).first ==
fe.system_to_component_index(j).first)
{
+
+ // If both DoFs i and j belong to same component, i.e. their shape functions are
+ // both $\phi$'s or both $\psi$'s, the contribution will end up in one of the diagonal
+ // blocks in our system matrix, and since the corresponding entries are computed
+ // by the same formula, we do not bother if they actually are
+ // $\phi$ or $\psi$ shape functions. We can simply compute the entry
+ // by iterating over all quadrature points and adding up their contributions,
+ // where values and gradients of the shape functions are supplied by our
+ // FEValues object.
+
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_matrix(i,j) += (((fe_values.shape_value(i,q_point) *
fe_values.shape_value(j,q_point)) *
- (- om2)
+ (- omega * omega)
+
(fe_values.shape_grad(i,q_point) *
fe_values.shape_grad(j,q_point)) *
- c2) *
+ c * c) *
fe_values.JxW(q_point));
+
+ // You might think that we would have to specify which
+ // component of the shape function we'd like to evaluate when requesting shape
+ // function values or gradients from the FEValues object. However, as the shape
+ // functions are primitive, they have only one nonzero component, and the
+ // FEValues class is smart enough to figure out that we are definitely interested in
+ // this one nonzero component.
}
}
}
+
+ // For DoFs that belong to different components of the system, i.e. one DoF
+ // representing a $\phi$ and the other a $\psi$, a contribution is only
+ // possible in the off-diagonal blocks of the system matrix. The entries
+ // in these blocks consist of a boundary integral on $\Gamma_2$, so we
+ // should first check if the current cell is on the boundary at all, since
+ // if it is not, its shape functions will certainly not have support on the boundary.
if (cell->at_boundary())
+
+ // If the current cell is at the boundary, we look through its
+ // faces to identify the ones that lie on $\Gamma_2$:
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary() &&
(cell->face(face)->boundary_indicator() == 0) )
{
+
+
+ // These faces will certainly contribute to the off-diagonal blocks of the
+ // system matrix, so we ask the FEFaceValues object to provide us with the
+ // shape function values on this face:
fe_face_values.reinit (cell, face);
+
+ // Next, we loop through all DoFs of the current cell to find pairs that
+ // belong to different components and both have support on the current
+ // face:
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
if ((fe.system_to_component_index(i).first !=
fe.has_support_on_face(i, face) &&
fe.has_support_on_face(j, face))
+
+ // These DoFs will then contribute to the boundary integrals
+ // in the off-diagonal blocks of the system matrix. To compute the
+ // integral, we loop over all the quadrature points on the face and
+ // sum up the contribution weighted with the quadrature weights that
+ // the face quadrature rule provides.
+ // In contrast to the entries on the diagonal blocks, here it does
+ // matter which one of the shape functions is a $\psi$ and which one
+ // is a $\phi$, since that will determine the sign of the entry.
+ // We account for this by a simple conditional statement
+ // that determines the correct sign. Since we already checked
+ // that DoF i and j belong to different components, so it suffices here
+ // to test for one of them to which component it belongs.
for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
- cell_matrix(i,j) += ((fe.system_to_component_index(i).first) ? 1 : (-1)) *
+ cell_matrix(i,j) += ((fe.system_to_component_index(i).first) ? -1 : 1) *
fe_face_values.shape_value(i,q_point) *
fe_face_values.shape_value(j,q_point) *
c *
fe_face_values.JxW(q_point);
}
+ // Now we are done with this cell and have to transfer its contributions
+ // from the local to the global system matrix. To this end,
+ // we first get a list of the global indices of the this cells DoFs:
cell->get_dof_indices (local_dof_indices);
+
+ // and then add the entries to the system matrix one by one:
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add (local_dof_indices[i],
cell_matrix(i,j));
}
+
+ // The only thing left are the Dirichlet boundary values on
+ // $\Gamma_1$, which is characterized by the boundary
+ // indicator 1. The Dirichlet values are provided by
+ // the <code>DirichletBoundaryValues</code> class we defined above:
std::map<unsigned int,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
1,
system_matrix,
solution,
system_rhs);
+
+ clock_t end = clock();
+ deallog << "done ("
+ << (end - start) / tps
+ << "s)"
+ << std::endl;
}
+
+ // @sect4{<code>UltrasoundProblem::solve</code>}
+
template <int dim>
void UltrasoundProblem<dim>::solve ()
{
+ deallog << "Solving linear system... ";
+ clock_t start = clock();
+
+ // As already mentioned in the introduction, the system matrix
+ // is neither symmetric nor definite, and so it is not
+ // quite obvious how to come up with an iterative solver
+ // and a preconditioner that do a good job on this matrix.
+ // We chose instead to go a different way and solve the linear
+ // system with the sparse LU decomposition provided by
+ // UMFPACK. This is often a good first choice for 2D problems
+ // and works reasonably well even for a large number of DoFs.
+ // The deal.II interface to UMFPACK is given by the SparseDirectUMFPACK
+ // class, which is very easy to use and allows us to solve our
+ // linear system with just 3 lines of code.
+
+ // Note again that for compiling this example program, you need
+ // to have the deal.II library built with UMFPACK support, which
+ // can be achieved by providing the <code> --with-umfpack</code>
+ // switch to the configure script prior to compilation of the library.
SparseDirectUMFPACK A_direct;
+ // The <code>initialize</code> call provides the matrix that we would like to invert
+ // to the SparseDirectUMFPACK object, and at the same
+ // time kicks off the LU-decomposition. Hence, this is also the point
+ // where most of the computational work in this program happens.
A_direct.initialize(system_matrix);
+
+ // After the decomposition, we can use <code>A_direct</code> like a matrix representing
+ // the inverse of our system matrix, so to compute the solution we just have
+ // to multiply with the right hand side vector:
A_direct.vmult(solution,system_rhs);
+
+ clock_t end = clock();
+ deallog << "done ("
+ << (end - start) / tps
+ << "s)"
+ << std::endl;
}
+
+ // @sect4{<code>UltrasoundProblem::output_results</code>}
+
+ // Here we output our solution $v$ and $w$ as well as the
+ // derived quantity $|u|$ in the
+ // format specified in the parameter file. Most of the
+ // work for deriving $|u|$ from $v$ and $w$ was already
+ // done in the implementation of the <code>ComputeIntensity</code> class,
+ // so that the output routine is rather straightforward and very similar
+ // to what is done in the previous tutorials.
template <int dim>
void UltrasoundProblem<dim>::output_results () const
{
- Postprocessor<dim> pproc;
+ deallog << "Generating output... ";
+ clock_t start = clock();
+
+ // Define objects of our <code>ComputeIntensity</code> class and a DataOut
+ // object:
+ ComputeIntensity<dim> intensities;
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
+ // Next we query the output-related parameters from the ParameterHandler:
prm.enter_subsection("Output parameters");
const std::string output_file = prm.get("Output file"),
output_format = prm.get("Output format");
+ // The DataOut::parse_parameters call acts as a counterpart to the
+ // DataOutInterface<1>::declare_parameters call in
+ // <code>ParameterReader::declare_parameters</code>. It collects all
+ // the output format related parameters from the ParameterHandler
+ // and sets the corresponding properties of the
+ // DataOut object accordingly.
data_out.parse_parameters(prm);
prm.leave_subsection ();
+ // Since the ParameterHandler provides the output format
+ // parameter as a string, we need to convert it to
+ // a format flag that can be understood by the DataOut object.
+ // The following function takes care of this:
DataOutBase::OutputFormat format = DataOutBase::parse_output_format(output_format);
+ // Now we put together the filename from the base name provided
+ // by the ParameterHandler and the suffix which is derived
+ // from the format by the DataOutBase::default_suffix function:
const std::string filename = output_file +
DataOutBase::default_suffix(format);
std::ofstream output (filename.c_str());
- data_out.add_data_vector (solution, pproc);
+ // The solution vectors $v$ and $w$ are added to the DataOut
+ // object in the usual way:
+ std::vector<std::string> solution_names;
+ solution_names.push_back ("Re_u");
+ solution_names.push_back ("Im_u");
+
+ data_out.add_data_vector (solution, solution_names);
+
+ // For the intensity, we just call <code>add_data_vector</code> again,
+ // but this with our <code>ComputeIntensity</code> object as the second argument,
+ // which effectively adds $|u|$ to the output data:
+ data_out.add_data_vector (solution, intensities);
+
+ // The last steps are as before:
data_out.build_patches ();
data_out.write (output, format);
+
+ clock_t end = clock();
+ deallog << "done ("
+ << (end - start) / tps
+ << "s)"
+ << std::endl;
}
+
+ // @sect4{<code>UltrasoundProblem::run</code>}
+ // Here we simply execute our functions one after the other:
template <int dim>
void UltrasoundProblem<dim>::run ()
{
}
+ // @sect4{The <code>main</code> function}
+
+ // Finally the <code>main</code> function of the program:
int main ()
{
try
{
- ParameterHandler prm;
+ // In 1D, the description of the domain
+ // and the boundary conditions is not very sensible, so
+ // exclude this case:
+ Assert (DIM > 1, ExcNotImplemented());
+ // Next define ParameterHandler and <code>ParameterReader</code> objects,
+ // and let the latter read in the parameter values from
+ // a textfile called <code>step-29.prm</code>:
+ ParameterHandler prm;
ParameterReader param(prm);
- param.read_parameters();
-
- Assert (DIM > 1, ExcNotImplemented());
+ param.read_parameters("step-29.prm");
+ // Lastly, we instantiate our main class with the ParameterHandler
+ // object and start the computations:
UltrasoundProblem<DIM> ultrasound_problem (prm);
ultrasound_problem.run ();
}
<< std::endl;
return 1;
}
-
return 0;
}