// @sect3{Include files}
// The most fundamental class in the
- // library is the ``Triangulation''
+ // library is the <code>Triangulation</code>
// class, which is declared here:
#include <grid/tria.h>
// We need the following two includes
// Now we want to write a graphical
// representation of the mesh to an
- // output file. The ``GridOut''
+ // output file. The <code>GridOut</code>
// class of deal.II can do that in
// a number of different output
// formats; here, we choose
// marked for further
// refinement, obviously). By
// convention, we almost always
- // use the names ``cell'' and
- // ``endc'' for the iterator
+ // use the names <code>cell</code> and
+ // <code>endc</code> for the iterator
// pointing to the present cell
// and to the
- // ``one-past-the-end''
+ // <code>one-past-the-end</code>
// iterator:
Triangulation<2>::active_cell_iterator
cell = triangulation.begin_active(),
// by which we find out about
// the number of vertices of
// a cell. Using the
- // ``GeometryInfo'' class, we
+ // <code>GeometryInfo</code> class, we
// will later have an easier
// time getting the program
// to also run in 3d: we only
// have to change all
- // occurrences of ``@<2@>'' to
- // ``@<3@>'', and do not have
+ // occurrences of <code>@<2@></code> to
+ // <code>@<3@></code>, and do not have
// to audit our code for the
// hidden appearance of magic
// numbers like a 4 that
#include <fe/fe_values.h>
// This is the only new one: in it,
- // we declare the ``MappingQ'' class
+ // we declare the <code>MappingQ</code> class
// which we will use for polynomial
// mappings of arbitrary order:
#include <fe/mapping_q.h>
// range of the number of digits
// which a double variable can hold,
// we rather declare the reference
- // value as a ``long double'' and
+ // value as a <code>long double</code> and
// give it a number of extra digits:
const long double pi = 3.141592653589793238462643;
// generates a triangulation of a
// circle (hyperball) and outputs the
// Qp mapping of its cells for
- // different values of ``p''. Then,
+ // different values of <code>p</code>. Then,
// we refine the grid once and do so
// again.
template <int dim>
// reasonable character sets
// nowadays), but also assumes
// that the increment
- // ``refinement'' is less than
+ // <code>refinement</code> is less than
// ten. This is therefore more
// a quick hack if we know
// exactly the values which the
// increment can assume. A
// better implementation would
// use the
- // ``std::istringstream''
+ // <code>std::istringstream</code>
// class to generate a name.
std::string filename_base = "ball";
filename_base += '0'+refinement;
// For this, first set up
// an object describing the
// mapping. This is done
- // using the ``MappingQ''
+ // using the <code>MappingQ</code>
// class, which takes as
// argument to the
// constructor the
// fact: if you want a
// piecewise linear
// mapping, then you could
- // give a value of ``1'' to
+ // give a value of <code>1</code> to
// the
// constructor. However,
// for linear mappings, so
// many things can be
// generated simpler that
// there is another class,
- // called ``MappingQ1''
+ // called <code>MappingQ1</code>
// which does exactly the
// same is if you gave an
- // degree of ``1'' to the
- // ``MappingQ'' class, but
+ // degree of <code>1</code> to the
+ // <code>MappingQ</code> class, but
// does so significantly
- // faster. ``MappingQ1'' is
+ // faster. <code>MappingQ1</code> is
// also the class that is
// implicitly used
// throughout the library
// object. This argument
// has a default value, and
// if no value is given a
- // simple ``MappingQ1''
+ // simple <code>MappingQ1</code>
// object is taken, which
// we briefly described
// above. This would then
// provides the corresponding `JxW'
// values of each cell. (Note that
// `JxW' is meant to abbreviate
- // ``Jacobian determinant times
- // weight''; since in numerical
+ // <code>Jacobian determinant times
+ // weight</code>; since in numerical
// quadrature the two factors always
// occur at the same places, we only
// offer the combined quantity,
// tells the FEValues object
// that it needs not compute
// other quantities upon
- // calling the ``reinit''
+ // calling the <code>reinit</code>
// function, thus saving
// computation time.
//
// long double) function
// implemented. Note that
// this also concerns the
- // second call as the ``fabs''
- // function in the ``std''
+ // second call as the <code>fabs</code>
+ // function in the <code>std</code>
// namespace is overloaded on
// its argument types, so there
// exists a version taking
- // and returning a ``long double'',
+ // and returning a <code>long double</code>,
// in contrast to the global
// namespace where only one such
// function is declared (which
// Just this one is new: it declares
// a class
- // ``CompressedSparsityPattern'',
+ // <code>CompressedSparsityPattern</code>,
// which we will use and explain
// further down below.
#include <lac/compressed_sparsity_pattern.h>
// class looks rather the same, with
// the sole structural difference
// that the functions
- // ``assemble_system'' now calls
- // ``solve'' itself, and is thus
- // called ``assemble_and_solve'', and
+ // <code>assemble_system</code> now calls
+ // <code>solve</code> itself, and is thus
+ // called <code>assemble_and_solve</code>, and
// that the output function was
// dropped since the solution
// function is so boring that it is
// Construct such an object, by
// initializing the variables. Here,
// we use linear finite elements (the
- // argument to the ``fe'' variable
+ // argument to the <code>fe</code> variable
// denotes the polynomial degree),
// and mappings of given order. Print
// to screen what we are about to do.
// The first task is to set up the
// variables for this problem. This
// includes generating a valid
- // ``DoFHandler'' object, as well as
+ // <code>DoFHandler</code> object, as well as
// the sparsity patterns for the
// matrix, and the object
// representing the constraints that
// this, we first want a list of
// those nodes which are actually
// at the boundary. The
- // ``DoFTools'' class has a
+ // <code>DoFTools</code> class has a
// function that returns an array
- // of boolean values where ``true''
+ // of boolean values where <code>true</code>
// indicates that the node is at
// the boundary. The second
// argument denotes a mask
// have a scalar finite element
// anyway, this mask consists of
// only one entry, and its value
- // must be ``true''.
+ // must be <code>true</code>.
std::vector<bool> boundary_dofs (dof_handler.n_dofs(), false);
DoFTools::extract_boundary_dofs (dof_handler, std::vector<bool>(1,true),
boundary_dofs);
// first pick out the first
// boundary node from this list. We
// do that by searching for the
- // first ``true'' value in the
- // array (note that ``std::find''
+ // first <code>true</code> value in the
+ // array (note that <code>std::find</code>
// returns an iterator to this
// element), and computing its
// distance to the overall first
// computation on a once coarser
// grid), then add this one line
// constraining the
- // ``first_boundary_dof'' to the
+ // <code>first_boundary_dof</code> to the
// sum of other boundary DoFs each
// with weight -1. Finally, close
// the constraints object, i.e. do
// sparsity pattern. This is indeed
// a tricky task here. Usually, we
// just call
- // ``DoFTools::make_sparsity_pattern''
+ // <code>DoFTools::make_sparsity_pattern</code>
// and condense the result using
// the hanging node constraints. We
// have no hanging node constraints
// we have this global constraint
// on the boundary. This poses one
// severe problem in this context:
- // the ``SparsityPattern'' class
+ // the <code>SparsityPattern</code> class
// wants us to state beforehand the
// maximal number of entries per
// row, either for all rows or for
// can tell you this number in case
// you just have hanging node
// constraints (namely
- // ``DoFHandler::max_coupling_between_dofs''),
+ // <code>DoFHandler::max_coupling_between_dofs</code>),
// but how is this for the present
// case? The difficulty arises
// because the elimination of the
// given that allows allocation of
// only a reasonable amount of
// memory, there is a class
- // ``CompressedSparsityPattern'',
+ // <code>CompressedSparsityPattern</code>,
// that can help us out here. It
// does not require that we know in
// advance how many entries rows
// initializing it with the
// dimensions of the matrix and
// calling another function
- // ``DoFTools::make_sparsity_pattern''
+ // <code>DoFTools::make_sparsity_pattern</code>
// to get the sparsity pattern due
// to the differential operator,
// then condense it with the
// Finally, once we have the full
// pattern, we can initialize an
// object of type
- // ``SparsityPattern'' from it and
+ // <code>SparsityPattern</code> from it and
// in turn initialize the matrix
// with it. Note that this is
// actually necessary, since the
- // ``CompressedSparsityPattern'' is
+ // <code>CompressedSparsityPattern</code> is
// so inefficient compared to the
- // ``SparsityPattern'' class due to
+ // <code>SparsityPattern</code> class due to
// the more flexible data
// structures it has to use, that
// we can impossibly base the
// sparse matrix class on it, but
// rather need an object of type
- // ``SparsityPattern'', which we
+ // <code>SparsityPattern</code>, which we
// generate by copying from the
// intermediate object.
//
// As a further sidenote, you will
// notice that we do not explicitly
- // have to ``compress'' the
+ // have to <code>compress</code> the
// sparsity pattern here. This, of
// course, is due to the fact that
- // the ``copy_from'' function
+ // the <code>copy_from</code> function
// generates a compressed object
// right from the start, to which
// you cannot add new entries
- // anymore. The ``compress'' call
+ // anymore. The <code>compress</code> call
// is therefore implicit in the
- // ``copy_from'' call.
+ // <code>copy_from</code> call.
sparsity_pattern.copy_from (csp);
system_matrix.reinit (sparsity_pattern);
}
// of a right hand side vector from
// body or boundary forces. They
// take the mapping object, the
- // ``DoFHandler'' object
+ // <code>DoFHandler</code> object
// representing the degrees of
// freedom and the finite element
// in use, a quadrature formula to
// operator. For this reason, there
// are quite a large number of
// variants of these functions in
- // the ``MatrixCreator'' and
- // ``MatrixTools''
+ // the <code>MatrixCreator</code> and
+ // <code>MatrixTools</code>
// classes. Whenever you need a
// slightly different version of
// these functions than the ones
// mentioned above) has to be at
// least 2, this makes up for the
// formula above computing
- // ``gauss_degree''.
+ // <code>gauss_degree</code>.
//
// Since the generation of the body
// force contributions to the right
// later add them together. The
// reason we had to do so is that
// the
- // ``VectorTools::create_right_hand_side''
+ // <code>VectorTools::create_right_hand_side</code>
// and
- // ``VectorTools::create_boundary_right_hand_side''
+ // <code>VectorTools::create_boundary_right_hand_side</code>
// functions first clear the output
// vector, rather than adding up
// their results to previous
// a function in the library that
// does this, although in a
// slightly non-obvious way: the
- // ``VectorTools::integrate_difference''
+ // <code>VectorTools::integrate_difference</code>
// function integrates the norm of
// the difference between a finite
// element function and a
// (which we make us of here), and
// the one which we have used in
// previous examples which
- // implicitly uses ``MappingQ1''.
+ // implicitly uses <code>MappingQ1</code>.
// Also note that we take a
// quadrature formula of one degree
// higher, in order to avoid
// mappings of linear through
// cubic mappings. Note that
// since we need the object of
- // type ``LaplaceProblem@<2@>''
+ // type <code>LaplaceProblem@<2@></code>
// only once, we do not even
// name it, but create an
// unnamed such object and call
- // the ``run'' function of it,
+ // the <code>run</code> function of it,
// subsequent to which it is
// immediately destroyed again.
for (unsigned int mapping_degree=1; mapping_degree<=3; ++mapping_degree)
#include <numerics/data_out.h>
// This is the first new file. It
- // declares the ``MappingQ1'' class
+ // declares the <code>MappingQ1</code> class
// that gives the standard bilinear
// mapping. For bilinear mappings use
// an object of this class rather
// than an object of the
- // ``MappingQ(1)'' class, as the
- // ``MappingQ1'' class is optimized
+ // <code>MappingQ(1)</code> class, as the
+ // <code>MappingQ1</code> class is optimized
// due to the pre-knowledge of the
// actual polynomial degree 1.
#include <fe/mapping_q1.h>
// programs -- there isn't much user
// interaction with finite element
// classes at all: the are passed to
- // ``DoFHandler'' and ``FEValues''
+ // <code>DoFHandler</code> and <code>FEValues</code>
// objects, and that is about it.
#include <fe/fe_dgq.h>
// We are going to use the simplest
// refinement indicator.
#include <numerics/derivative_approximation.h>
// Finally we do some time comparison
- // using the ``Timer'' class.
+ // using the <code>Timer</code> class.
#include <base/timer.h>
// And this again is C++:
//
// First we define the classes
// representing the equation-specific
- // functions. Both classes, ``RHS''
- // and ``BoundaryValues'', are
- // derived from the ``Function''
- // class. Only the ``value_list''
+ // functions. Both classes, <code>RHS</code>
+ // and <code>BoundaryValues</code>, are
+ // derived from the <code>Function</code>
+ // class. Only the <code>value_list</code>
// function are implemented because
// only lists of function values are
// computed rather than single
};
- // The class ``Beta'' represents the
+ // The class <code>Beta</code> represents the
// vector valued flow field of the
// linear transport equation and is
- // not derived from the ``Function''
+ // not derived from the <code>Function</code>
// class as we prefer to get function
- // values of type ``Point'' rather
+ // values of type <code>Point</code> rather
// than of type
- // ``Vector@<double@>''. This, because
+ // <code>Vector@<double@></code>. This, because
// there exist scalar products
- // between ``Point'' and ``Point'' as
- // well as between ``Point'' and
- // ``Tensor'', simplifying terms like
+ // between <code>Point</code> and <code>Point</code> as
+ // well as between <code>Point</code> and
+ // <code>Tensor</code>, simplifying terms like
// $\beta\cdot n$ and
// $\beta\cdot\nabla v$.
//
// The implementation of the
- // ``value_list'' functions of these
+ // <code>value_list</code> functions of these
// classes are rather simple. For
// simplicity the right hand side is
// set to be zero but will be
// Next we define the
// equation-dependent and
// DG-method-dependent class
- // ``DGTransportEquation''. Its
+ // <code>DGTransportEquation</code>. Its
// member functions were already
// mentioned in the Introduction and
// will be explained
// below. Furthermore it includes
// objects of the previously defined
- // ``Beta'', ``RHS'' and
- // ``BoundaryValues'' function
+ // <code>Beta</code>, <code>RHS</code> and
+ // <code>BoundaryValues</code> function
// classes.
template <int dim>
class DGTransportEquation
// @sect4{Function: assemble_cell_term}
//
- // The ``assemble_cell_term''
+ // The <code>assemble_cell_term</code>
// function assembles the cell terms
// of the discretization.
- // ``ui_vi_matrix'' is a cell matrix,
+ // <code>ui_vi_matrix</code> is a cell matrix,
// i.e. for a DG method of degree 1,
// it is of size 4 times 4, and
- // ``cell_vector'' is of size 4.
+ // <code>cell_vector</code> is of size 4.
// When this function is invoked,
- // ``fe_v'' is already reinit'ed with the
+ // <code>fe_v</code> is already reinit'ed with the
// current cell before and includes
// all shape values needed.
template <int dim>
FullMatrix<double> &ui_vi_matrix,
Vector<double> &cell_vector) const
{
- // First we ask ``fe_v'' for the
+ // First we ask <code>fe_v</code> for the
// quadrature weights,
const std::vector<double> &JxW = fe_v.get_JxW_values ();
// Then the flow field beta and the
- // ``rhs_function'' are evaluated at
+ // <code>rhs_function</code> are evaluated at
// the quadrature points,
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
std::vector<double> rhs (fe_v.n_quadrature_points);
// @sect4{Function: assemble_boundary_term}
//
- // The ``assemble_boundary_term''
+ // The <code>assemble_boundary_term</code>
// function assembles the face terms
// at boundary faces. When this
- // function is invoked, ``fe_v'' is
+ // function is invoked, <code>fe_v</code> is
// already reinit'ed with the current
// cell and current face. Hence it
// provides the shape values on that
{
// Again, as in the previous
// function, we ask the
- // ``FEValues'' object for the
+ // <code>FEValues</code> object for the
// quadrature weights
const std::vector<double> &JxW = fe_v.get_JxW_values ();
// but here also for the normals.
// @sect4{Function: assemble_face_term1}
//
- // The ``assemble_face_term1''
+ // The <code>assemble_face_term1</code>
// function assembles the face terms
// corresponding to the first version
// of the DG method, cf. above. For
// all cell boundaries.
//
// When this function is invoked,
- // ``fe_v'' and ``fe_v_neighbor'' are
+ // <code>fe_v</code> and <code>fe_v_neighbor</code> are
// already reinit'ed with the current
// cell and the neighoring cell,
// respectively, as well as with the
// on the face.
//
// In addition to the cell matrix
- // ``ui_vi_matrix'' this function
+ // <code>ui_vi_matrix</code> this function
// gets a new argument
- // ``ue_vi_matrix'', that stores
+ // <code>ue_vi_matrix</code>, that stores
// contributions to the system matrix
// that are based on exterior values
// of $u$ and interior values of
- // $v$. Here we note that ``ue'' is
- // the short notation for ``u
- // exterior'' and represents $u_h^-$,
+ // $v$. Here we note that <code>ue</code> is
+ // the short notation for <code>u
+ // exterior</code> and represents $u_h^-$,
// see the introduction.
template <int dim>
void DGTransportEquation<dim>::assemble_face_term1(
// @sect4{Function: assemble_face_term2}
//
// Now we look at the
- // ``assemble_face_term2'' function
+ // <code>assemble_face_term2</code> function
// that assembles the face terms
// corresponding to the second
// version of the DG method,
// terms are given as a sum of
// integrals over all faces. Here we
// need two additional cell matrices
- // ``ui_ve_matrix'' and
- // ``ue_ve_matrix'' that will store
+ // <code>ui_ve_matrix</code> and
+ // <code>ue_ve_matrix</code> that will store
// contributions due to terms
// involving ui and ve as well as ue
// and ve.
// After these preparations, we
// proceed with the main part of this
// program. The main class, here
- // called ``DGMethod'' is basically
+ // called <code>DGMethod</code> is basically
// the main class of step-6. One of
// the differences is that there's no
// ConstraintMatrix object. This is,
// solutions to the problems
// corresponding to the two
// different assembling routines
- // ``assemble_system1'' and
- // ``assemble_system2'';
+ // <code>assemble_system1</code> and
+ // <code>assemble_system2</code>;
Vector<double> solution1;
Vector<double> solution2;
Vector<double> right_hand_side;
// @sect4{Function: assemble_system1}
//
// We proceed with the
- // ``assemble_system1'' function that
+ // <code>assemble_system1</code> function that
// implements the DG discretization
// in its first version. This
// function repeatedly calls the
- // ``assemble_cell_term'',
- // ``assemble_boundary_term'' and
- // ``assemble_face_term1'' functions
- // of the ``DGTransportEquation''
+ // <code>assemble_cell_term</code>,
+ // <code>assemble_boundary_term</code> and
+ // <code>assemble_face_term1</code> functions
+ // of the <code>DGTransportEquation</code>
// object. The
- // ``assemble_boundary_term'' covers
+ // <code>assemble_boundary_term</code> covers
// the first case mentioned in the
// introduction.
//
// 1. face is at boundary
//
// This function takes a
- // ``FEFaceValues'' object as
+ // <code>FEFaceValues</code> object as
// argument. In contrast to that
- // ``assemble_face_term1''
- // takes two ``FEFaceValuesBase''
+ // <code>assemble_face_term1</code>
+ // takes two <code>FEFaceValuesBase</code>
// objects; one for the shape
// functions on the current cell and
// the other for shape functions on
// the neighboring cell under
// consideration. Both objects are
- // either of class ``FEFaceValues''
- // or of class ``FESubfaceValues''
+ // either of class <code>FEFaceValues</code>
+ // or of class <code>FESubfaceValues</code>
// (both derived from
- // ``FEFaceValuesBase'') according to
+ // <code>FEFaceValuesBase</code>) according to
// the remaining cases mentioned
// in the introduction:
//
// 2. neighboring cell is finer
- // (current cell: ``FESubfaceValues'',
- // neighboring cell: ``FEFaceValues'');
+ // (current cell: <code>FESubfaceValues</code>,
+ // neighboring cell: <code>FEFaceValues</code>);
//
// 3. neighboring cell is of the same
// refinement level (both, current
// and neighboring cell:
- // ``FEFaceValues'');
+ // <code>FEFaceValues</code>);
//
// 4. neighboring cell is coarser
- // (current cell: ``FEFaceValues'',
+ // (current cell: <code>FEFaceValues</code>,
// neighboring cell:
- // ``FESubfaceValues'').
+ // <code>FESubfaceValues</code>).
//
// If we considered globally refined
// meshes then only case 3 would
std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
// First we create the
- // ``update_flags'' for the
- // ``FEValues'' and the
- // ``FEFaceValues'' objects.
+ // <code>update_flags</code> for the
+ // <code>FEValues</code> and the
+ // <code>FEFaceValues</code> objects.
const UpdateFlags update_flags = update_values
| update_gradients
| update_q_points
// vectors of the current cell.
const UpdateFlags neighbor_face_update_flags = update_values;
- // Then we create the ``FEValues''
+ // Then we create the <code>FEValues</code>
// object. Note, that since version
// 3.2.0 of deal.II the constructor
// of this class takes a
- // ``Mapping'' object as first
+ // <code>Mapping</code> object as first
// argument. Although the
- // constructor without ``Mapping''
+ // constructor without <code>Mapping</code>
// argument is still supported it
// is recommended to use the new
// constructor. This reduces the
// effect of `hidden magic' (the
// old constructor implicitely
- // assumes a ``MappingQ1'' mapping)
+ // assumes a <code>MappingQ1</code> mapping)
// and makes it easier to change
// the mapping object later.
FEValues<dim> fe_v (
mapping, fe, quadrature, update_flags);
// Similarly we create the
- // ``FEFaceValues'' and
- // ``FESubfaceValues'' objects for
+ // <code>FEFaceValues</code> and
+ // <code>FESubfaceValues</code> objects for
// both, the current and the
// neighboring cell. Within the
// following nested loop over all
// and vectors. Here we need two
// cell matrices, both for face
// terms that include test
- // functions ``vi'' (internal shape
+ // functions <code>vi</code> (internal shape
// functions, i.e. shape functions
// of the current cell). To be more
// precise, the first matrix will
for (;cell!=endc; ++cell)
{
// In the
- // ``assemble_face_term1''
+ // <code>assemble_face_term1</code>
// function contributions to
// the cell matrices and the
// cell vector are only
// ADDED. Therefore on each
// cell we need to reset the
- // ``ui_vi_matrix'' and
- // ``cell_vector'' to zero,
+ // <code>ui_vi_matrix</code> and
+ // <code>cell_vector</code> to zero,
// before assembling the cell terms.
ui_vi_matrix = 0;
cell_vector = 0;
- // Now we reinit the ``FEValues''
+ // Now we reinit the <code>FEValues</code>
// object for the current cell
fe_v.reinit (cell);
// and call the function
// that assembles the cell
// terms. The first argument is
- // the ``FEValues'' that was
+ // the <code>FEValues</code> that was
// previously reinit'ed on the
// current cell.
dg.assemble_cell_term(fe_v,
typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
// and clear the
- // ``ue_vi_matrix'' on each
+ // <code>ue_vi_matrix</code> on each
// face.
ue_vi_matrix = 0;
if (face->at_boundary())
{
// We reinit the
- // ``FEFaceValues''
+ // <code>FEFaceValues</code>
// object to the
// current face
fe_v_face.reinit (cell, face_no);
// note that the
// following part of
// code will not work
- // for ``dim==1''.
+ // for <code>dim==1</code>.
if (face->has_children())
{
// First we store
// neighbor-@>neighbor(neighbor2)
// equals the
// current cell
- // ``cell''.
+ // <code>cell</code>.
const unsigned int neighbor2=
cell->neighbor_of_neighbor(face_no);
// and set the
// cell
// iterator
- // ``neighbor_child''
+ // <code>neighbor_child</code>
// to the cell
// placed
// `behind' the
// We need to
// reset the
- // ``ue_vi_matrix''
+ // <code>ue_vi_matrix</code>
// on each
// subface
// because on
// each subface
- // the ``un''
+ // the <code>un</code>
// belong to
// different
// neighboring
// case (case
// 2) we employ
// the
- // ``FESubfaceValues''
+ // <code>FESubfaceValues</code>
// of the
// current
// cell (here
// and
// distribute
- // ``ue_vi_matrix''
+ // <code>ue_vi_matrix</code>
// to the
// system_matrix
for (unsigned int i=0; i<dofs_per_cell; ++i)
system_matrix.add(dofs[i], dofs_neighbor[k],
ue_vi_matrix(i,k));
}
- // End of ``if
- // (face-@>has_children())''
+ // End of <code>if
+ // (face-@>has_children())</code>
}
else
{
// We reinit
// the
- // ``FEFaceValues''
+ // <code>FEFaceValues</code>
// of the
// current and
// neighboring
fe_v_face_neighbor,
ui_vi_matrix,
ue_vi_matrix);
- // End of ``if
+ // End of <code>if
// (neighbor-@>level()
// ==
- // cell-@>level())''
+ // cell-@>level())</code>
}
else
{
// Reinit the
// appropriate
- // ``FEFaceValues''
+ // <code>FEFaceValues</code>
// and assemble
// the face
// terms.
// Now we get the
// dof indices of
// the
- // ``neighbor_child''
+ // <code>neighbor_child</code>
// cell,
neighbor->get_dof_indices (dofs_neighbor);
// and distribute the
- // ``ue_vi_matrix''.
+ // <code>ue_vi_matrix</code>.
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int k=0; k<dofs_per_cell; ++k)
system_matrix.add(dofs[i], dofs_neighbor[k],
ue_vi_matrix(i,k));
}
- // End of ``face not at boundary'':
+ // End of <code>face not at boundary</code>:
}
// End of loop over all faces:
}
// Finally we distribute the
- // ``ui_vi_matrix''
+ // <code>ui_vi_matrix</code>
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));
// @sect4{Function: assemble_system2}
//
// We proceed with the
- // ``assemble_system2'' function that
+ // <code>assemble_system2</code> function that
// implements the DG discretization
// in its second version. This
// function is very similar to the
- // ``assemble_system1''
+ // <code>assemble_system1</code>
// function. Therefore, here we only
// discuss the differences between
// the two functions. This function
// repeatedly calls the
- // ``assemble_face_term2'' function
+ // <code>assemble_face_term2</code> function
// of the DGTransportEquation object,
// that assembles the face terms
// written as a sum of integrals over
const UpdateFlags neighbor_face_update_flags = update_values;
// Here we do not need
- // ``fe_v_face_neighbor'' as case 4
+ // <code>fe_v_face_neighbor</code> as case 4
// does not occur.
FEValues<dim> fe_v (
mapping, fe, quadrature, update_flags);
// Additionally we need the
// following two cell matrices,
// both for face term that include
- // test function ``ve'' (external
+ // test function <code>ve</code> (external
// shape functions, i.e. shape
// functions of the neighboring
// cell). To be more precise, the
// difference quotients including the
// cell under consideration and its
// neighbors. This is done by the
- // ``DerivativeApproximation'' class
+ // <code>DerivativeApproximation</code> class
// that computes the approximate
// gradients in a way similar to the
- // ``GradientEstimation'' described
+ // <code>GradientEstimation</code> described
// in step-9 of this tutorial. In
// fact, the
- // ``DerivativeApproximation'' class
+ // <code>DerivativeApproximation</code> class
// was developed following the
- // ``GradientEstimation'' class of
+ // <code>GradientEstimation</code> class of
// step-9. Relating to the
// discussion in step-9, here we
// consider $h^{1+d/2}|\nabla_h
template <int dim>
void DGMethod<dim>::refine_grid ()
{
- // The ``DerivativeApproximation''
+ // The <code>DerivativeApproximation</code>
// class computes the gradients to
// float precision. This is
// sufficient as they are
}
- // The following ``run'' function is
+ // The following <code>run</code> function is
// similar to previous examples. The
// only difference is that the
// problem is assembled and solved
// twice on each refinement step;
- // first by ``assemble_system1'' that
+ // first by <code>assemble_system1</code> that
// implements the first version and
- // then by ``assemble_system2'' that
+ // then by <code>assemble_system2</code> that
// implements the second version of
// writing the DG
// discretization. Furthermore the
}
}
- // The following ``main'' function is
+ // The following <code>main</code> function is
// similar to previous examples and
// need not to be commented on.
int main ()
// a list of include files from the
// library, and as usual they are in
// the standard order which is
- // ``base'' -- ``lac'' -- ``grid'' --
- // ``dofs'' -- ``fe'' -- ``numerics''
+ // <code>base</code> -- <code>lac</code> -- <code>grid</code> --
+ // <code>dofs</code> -- <code>fe</code> -- <code>numerics</code>
// (as each of these categories
// roughly builds upon previous
// ones), then C++ standard headers:
// equation. In fact, they can
// evaluate every kind of solution,
// as long as it is described by a
- // ``DoFHandler'' object, and a
+ // <code>DoFHandler</code> object, and a
// solution vector. We define them
// here first, even before the
// classes that actually generate the
// From an abstract point of view, we
// declare a pure base class
// that provides an evaluation
- // operator ``operator()'' which will
+ // operator <code>operator()</code> which will
// do the evaluation of the solution
// (whatever derived classes might
- // consider an ``evaluation''). Since
+ // consider an <code>evaluation</code>). Since
// this is the only real function of
// this base class (except for some
// bookkeeping machinery), one
// usually terms such a class that
- // only has an ``operator()'' a
- // ``functor'' in C++ terminology,
+ // only has an <code>operator()</code> a
+ // <code>functor</code> in C++ terminology,
// since it is used just like a
// function object.
//
// Now for the abstract base class
// of evaluation classes: its main
// purpose is to declare a pure
- // virtual function ``operator()''
- // taking a ``DoFHandler'' object,
+ // virtual function <code>operator()</code>
+ // taking a <code>DoFHandler</code> object,
// and the solution vector. In
// order to be able to use pointers
// to this base class only, it also
// like to extract a point value
// from the solution, so the first
// class does this in its
- // ``operator()''. The actual point
+ // <code>operator()</code>. The actual point
// is given to this class through
// the constructor, as well as a
// table object into which it will
// In the step-9 example program,
// we have already seen how such an
// exception class can be declared,
- // using the ``DeclExceptionN''
+ // using the <code>DeclExceptionN</code>
// macros. We use this mechanism
// here again.
//
// beginning of the
// function, for example
// by a statement like
- // ``Assert
+ // <code>Assert
// (dof_handler.get_fe().dofs_per_vertex
// @> 0,
- // ExcNotImplemented())'',
+ // ExcNotImplemented())</code>,
// which should make it
// quite clear what is
// going wrong when the
// that that does not
// hurt here, since the
// statement
- // ``cell-@>vertex_dof_index(vertex,0)''
+ // <code>cell-@>vertex_dof_index(vertex,0)</code>
// would fail if we asked
// it to give us the DoF
// index of a vertex if
// solution there and the rest of
// the computations were useless
// anyway. So make sure through
- // the ``AssertThrow'' macro
+ // the <code>AssertThrow</code> macro
// already used in the step-9
// program that we have indeed
// found this point. If this is
// exception of the type that is
// given to it as second
// argument, but compared to a
- // straightforward ``throw''
+ // straightforward <code>throw</code>
// statement, it fills the
// exception object with a set of
// additional information, for
// line number where the
// exception was generated, and
// the condition that failed. If
- // you have a ``catch'' clause in
+ // you have a <code>catch</code> clause in
// your main function (as this
// program has), you will catch
// all exceptions that are not
AssertThrow (evaluation_point_found,
ExcEvaluationPointNotFound(evaluation_point));
// Note that we have used the
- // ``Assert'' macro in other
+ // <code>Assert</code> macro in other
// example programs as well. It
// differed from the
- // ``AssertThrow'' macro used
+ // <code>AssertThrow</code> macro used
// here in that it simply aborts
// the program, rather than
// throwing an exception, and
// program in debug mode, but
// should be checked always, also
// in production runs. Thus the
- // use of the ``AssertThrow''
+ // use of the <code>AssertThrow</code>
// macro here.
// Now, if we are sure that we
// @sect4{Generating output}
// A different, maybe slightly odd
- // kind of ``evaluation'' of a
+ // kind of <code>evaluation</code> of a
// solution is to output it to a
// file in a graphical
// format. Since in the evaluation
// functions we are given a
- // ``DoFHandler'' object and the
+ // <code>DoFHandler</code> object and the
// solution vector, we have all we
// need to do this, so we can do it
// in an evaluation class. The
//
// Since this class which generates
// the output is derived from the
- // common ``EvaluationBase'' base
+ // common <code>EvaluationBase</code> base
// class, its main interface is the
- // ``operator()''
+ // <code>operator()</code>
// function. Furthermore, it has a
// constructor taking a string that
// will be used as the base part of
// we write).
//
// Regarding the output format, the
- // ``DataOutInterface'' class
+ // <code>DataOutInterface</code> class
// (which is a base class of
- // ``DataOut'' through which we
+ // <code>DataOut</code> through which we
// will access its fields) provides
// an enumeration field
- // ``OutputFormat'', which lists
+ // <code>OutputFormat</code>, which lists
// names for all supported output
// formats. At the time of writing
// of this program, the supported
// graphics formats are represented
- // by the enum values ``ucd'',
- // ``gnuplot'', ``povray'',
- // ``eps'', ``gmv'', ``tecplot'',
- // ``tecplot_binary'', ``dx'', and
- // ``vtk'', but this list will
+ // by the enum values <code>ucd</code>,
+ // <code>gnuplot</code>, <code>povray</code>,
+ // <code>eps</code>, <code>gmv</code>, <code>tecplot</code>,
+ // <code>tecplot_binary</code>, <code>dx</code>, and
+ // <code>vtk</code>, but this list will
// certainly grow over time. Now,
// within various functions of that
// base class, you can use values
// (for example the default suffix
// used for files of each format),
// and you can call a generic
- // ``write'' function, which then
+ // <code>write</code> function, which then
// branches to the
- // ``write_gnuplot'',
- // ``write_ucd'', etc functions
+ // <code>write_gnuplot</code>,
+ // <code>write_ucd</code>, etc functions
// which we have used in previous
// examples already, based on the
// value of a second argument given
// particularly interesting feature
// over previous example programs
// is the use of the
- // ``DataOut::default_suffix''
+ // <code>DataOut::default_suffix</code>
// function, returning the usual
// suffix for files of a given
// format (e.g. ".eps" for
// encapsulated postscript files,
// ".gnuplot" for Gnuplot files),
// and of the generic
- // ``DataOut::write'' function with
+ // <code>DataOut::write</code> function with
// a second argument, which
// branches to the actual output
// functions for the different
// passed as second argument.
//
// Also note that we have to prefix
- // ``this-@>'' to access a member
+ // <code>this-@></code> to access a member
// variable of the template
// dependent base class. The reason
// here, and further down in the
// program is the same as the one
// described in the step-7 example
- // program (look for ``two-stage
- // name lookup'' there).
+ // program (look for <code>two-stage
+ // name lookup</code> there).
template <int dim>
void
SolutionOutput<dim>::operator () (const DoFHandler<dim> &dof_handler,
// style used in Smalltalk or Java
// programs, where all classes are
// derived from entirely abstract
- // classes ``Object'', even number
+ // classes <code>Object</code>, even number
// representations. The author
// admits that he does not
// particularly like the use of
// have to make sure that the
// triangulation exists until the
// destructor exits. We do this by
- // keeping a ``SmartPointer'' to
+ // keeping a <code>SmartPointer</code> to
// this triangulation, which uses a
// counter in the triangulation
// class to denote the fact that
// by this we allow that derived
// classes refine or coarsen the
// triangulation within the
- // ``refine_grid'' function.
+ // <code>refine_grid</code> function.
//
// Finally, we have a function
- // ``n_dofs'' is only a tool for
+ // <code>n_dofs</code> is only a tool for
// the driver functions to decide
// whether we want to go on with
// mesh refinement or not. It
// solving it, and calling the
// postprocessor objects on the
// solution. It implements the
- // ``solve_problem'' and
- // ``postprocess'' functions
+ // <code>solve_problem</code> and
+ // <code>postprocess</code> functions
// declared in the base class. It
// does not, however, implement the
- // ``refine_grid'' method, as mesh
+ // <code>refine_grid</code> method, as mesh
// refinement will be implemented
// in a number of derived classes.
//
// It also declares a new abstract
// virtual function,
- // ``assemble_rhs'', that needs to
+ // <code>assemble_rhs</code>, that needs to
// be overloaded in subclasses. The
// reason is that we will implement
// two different classes that will
// down to the base class's
// constructor, or are stored and
// used to generate a
- // ``DoFHandler'' object
+ // <code>DoFHandler</code> object
// later. Since finite elements and
// quadrature formula should match,
// it is also passed a quadrature
// object.
//
- // The ``solve_problem'' sets up
+ // The <code>solve_problem</code> sets up
// the data structures for the
// actual solution, calls the
// functions to assemble the linear
// system, and solves it.
//
- // The ``postprocess'' function
+ // The <code>postprocess</code> function
// finally takes an evaluation
// object and applies it to the
// computed solution.
//
- // The ``n_dofs'' function finally
+ // The <code>n_dofs</code> function finally
// implements the pure virtual
// function of the base class.
template <int dim>
// of the class. It does not do
// much except store pointers to
// the objects given, and generate
- // ``DoFHandler'' object
+ // <code>DoFHandler</code> object
// initialized with the given
// pointer to a triangulation. This
// causes the DoF handler to store
// already generate a finite
// element numbering (we only ask
// for that in the
- // ``solve_problem'' function).
+ // <code>solve_problem</code> function).
template <int dim>
Solver<dim>::Solver (Triangulation<dim> &triangulation,
const FiniteElement<dim> &fe,
// As stated above, the
- // ``postprocess'' function takes
+ // <code>postprocess</code> function takes
// an evaluation object, and
// applies it to the computed
// solution. This function may be
}
- // The ``n_dofs'' function should
+ // The <code>n_dofs</code> function should
// be self-explanatory:
template <int dim>
unsigned int
// overwriting their
// respective
// work. Previously, we have
- // used the ``acquire'' and
- // ``release'' functions of
+ // used the <code>acquire</code> and
+ // <code>release</code> functions of
// the mutex to lock and
// unlock the mutex,
// respectively. While this
// in the middle of the
// locked block, and forgets
// that before we call
- // ``return'', we also have
+ // <code>return</code>, we also have
// to unlock the mutex. This
// all is not be a problem
// here, but we want to show
// mutex, and on running the
// destructor unlocks it
// again. This is called the
- // ``scoped lock'' pattern
+ // <code>scoped lock</code> pattern
// (apparently invented by
// Doug Schmidt originally),
// and it works because
// objects are also run when
// we exit the function
// either through a
- // ``return'' statement, or
+ // <code>return</code> statement, or
// when an exception is
// raised. Thus, it is
// guaranteed that the mutex
cell_matrix(i,j));
// Here, at the brace, the
// current scope ends, so the
- // ``lock'' variable goes out
+ // <code>lock</code> variable goes out
// of existence and its
// destructor the mutex is
// unlocked.
// and do the second action in the
// main thread. Since only one
// thread is generated, we don't
- // use the ``Threads::ThreadGroup''
+ // use the <code>Threads::ThreadGroup</code>
// class here, but rather use the
// one created thread object
// directly to wait for this
//
// Note that taking up the address
// of the
- // ``DoFTools::make_hanging_node_constraints''
+ // <code>DoFTools::make_hanging_node_constraints</code>
// function is a little tricky,
// since there are actually three
// of them, one for each supported
// addresses of overloaded
// functions is somewhat
// complicated in C++, since the
- // address-of operator ``&'' in
+ // address-of operator <code>&</code> in
// that case returns more like a
// set of values (the addresses of
// all functions with that name),
// would like to have; for this, we
// could use a cast, but for more
// clarity, we assign it to a
- // temporary ``mhnc_p'' (short for
- // ``pointer to
- // make_hanging_node_constraints'')
+ // temporary <code>mhnc_p</code> (short for
+ // <code>pointer to
+ // make_hanging_node_constraints</code>)
// with the right type, and using
// this pointer instead.
template <int dim>
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
// Wait until the
- // ``hanging_node_constraints''
+ // <code>hanging_node_constraints</code>
// object is fully set up, then
// close it and use it to
// condense the sparsity pattern:
// that denotes the right hand side
// of the problem. A pointer to
// this object is stored (again as
- // a ``SmartPointer'', in order to
+ // a <code>SmartPointer</code>, in order to
// make sure that the function
// object is not deleted as long as
// it is still used by this class).
//
// The only functional part of this
- // class is the ``assemble_rhs''
+ // class is the <code>assemble_rhs</code>
// method that does what its name
// suggests.
template <int dim>
- // ... as does the ``assemble_rhs''
+ // ... as does the <code>assemble_rhs</code>
// function. Since this is
// explained in several of the
// previous example programs, we
// By now, all functions of the
// abstract base class except for
- // the ``refine_grid'' function
+ // the <code>refine_grid</code> function
// have been implemented. We will
// now have two classes that
// implement this function for the
- // ``PrimalSolver'' class, one
+ // <code>PrimalSolver</code> class, one
// doing global refinement, one a
// form of local refinement.
//
// The first, doing global
// refinement, is rather simple:
// its main function just calls
- // ``triangulation-@>refine_global
- // (1);'', which does all the work.
+ // <code>triangulation-@>refine_global
+ // (1);</code>, which does all the work.
//
- // Note that since the ``Base''
- // base class of the ``Solver''
+ // Note that since the <code>Base</code>
+ // base class of the <code>Solver</code>
// class is virtual, we have to
// declare a constructor that
// initializes the immediate base
// exact solution the function
// $u(x,y)=exp(x+sin(10y+5x^2))$. In more
// than two dimensions, simply repeat
- // the sine-factor with ``y''
- // replaced by ``z'' and so on. Given
+ // the sine-factor with <code>y</code>
+ // replaced by <code>z</code> and so on. Given
// this, the following two classes
// are probably straightforward from
// the previous examples.
// your program).
for (unsigned int step=0; true; ++step)
{
- // Then give the ``alive''
+ // Then give the <code>alive</code>
// indication for this
// iteration. Note that the
- // ``std::flush'' is needed to
+ // <code>std::flush</code> is needed to
// have the text actually
// appear on the screen, rather
// than only in some buffer
// from all adjacent cells.
//
// Given the interface of the
- // ``PointValueEvaluation'' class,
+ // <code>PointValueEvaluation</code> class,
// the declaration of this class
// provides little surprise, and
// neither does the constructor:
// element fields at
// certain points is done
// through the
- // ``FEValues'' class, so
+ // <code>FEValues</code> class, so
// we use that. The
// question is: the
- // ``FEValues'' object
+ // <code>FEValues</code> object
// needs to be a given a
// quadrature formula and
// can then compute the
// above.
//
// Thus: initialize the
- // ``FEValues'' object on
+ // <code>FEValues</code> object on
// this cell,
fe_values.reinit (cell);
// and extract the
// the grids generated. This again
// can be done with one such
// class. Its structure is analog
- // to the ``SolutionOutput'' class
+ // to the <code>SolutionOutput</code> class
// of the previous example program,
// so we do not discuss it here in
// more detail. Furthermore,
// This class is almost unchanged,
// with the exception that it
// declares two more functions:
- // ``output_solution'' will be used
+ // <code>output_solution</code> will be used
// to generate output files from
// the actual solutions computed by
// derived classes, and the
- // ``set_refinement_cycle''
+ // <code>set_refinement_cycle</code>
// function by which the testing
// framework sets the number of the
// refinement cycle to a local
// @sect4{The Laplace Solver class}
- // Likewise, the ``Solver'' class
+ // Likewise, the <code>Solver</code> class
// is entirely unchanged and will
// thus not be discussed.
template <int dim>
// @sect4{The PrimalSolver class}
- // The ``PrimalSolver'' class is
+ // The <code>PrimalSolver</code> class is
// also mostly unchanged except for
// overloading the functions
- // ``solve_problem'', ``n_dofs'',
- // and ``postprocess'' of the base
+ // <code>solve_problem</code>, <code>n_dofs</code>,
+ // and <code>postprocess</code> of the base
// class, and implementing the
- // ``output_solution''
+ // <code>output_solution</code>
// function. These overloaded
// functions do nothing particular
// besides calling the functions of
// requires us to write such
// functions for the following
// scenario: Besides the
- // ``PrimalSolver'' class, we will
- // have a ``DualSolver'', both
- // derived from ``Solver''. We will
+ // <code>PrimalSolver</code> class, we will
+ // have a <code>DualSolver</code>, both
+ // derived from <code>Solver</code>. We will
// then have a final classes which
// derived from these two, which
// will then have two instances of
- // the ``Solver'' class as its base
+ // the <code>Solver</code> class as its base
// classes. If we want, for
// example, the number of degrees
// of freedom of the primal solver,
// we would have to indicate this
// like so:
- // ``PrimalSolver@<dim@>::n_dofs()''.
+ // <code>PrimalSolver@<dim@>::n_dofs()</code>.
// However, the compiler does not
- // accept this since the ``n_dofs''
+ // accept this since the <code>n_dofs</code>
// function is actually from a base
- // class of the ``PrimalSolver''
+ // class of the <code>PrimalSolver</code>
// class, so we have to inject the
// name from the base to the
// derived class using these
// additional functions.
//
// Regarding the implementation of
- // the ``output_solution''
+ // the <code>output_solution</code>
// function, we keep the
- // ``GlobalRefinement'' and
- // ``RefinementKelly'' classes in
+ // <code>GlobalRefinement</code> and
+ // <code>RefinementKelly</code> classes in
// this program, and they can then
// rely on the default
// implementation of this function
// everything that describes a test
// case: here, these are two
// subclasses, one called
- // ``BoundaryValues'' for the
+ // <code>BoundaryValues</code> for the
// boundary values of the exact
// solution, and one called
- // ``RightHandSide'', and then a way
+ // <code>RightHandSide</code>, and then a way
// to generate the coarse grid. Since
// the solution of the previous
// example program looked like curved
// right hand side by simply giving
// the name of the outer class as a
// template argument to a class which
- // we call here ``Data::SetUp'', and
+ // we call here <code>Data::SetUp</code>, and
// it then creates objects for the
// inner classes. In this case, to
// get all that characterizes the
// curved ridge solution, we would
// simply generate an instance of
- // ``Data::SetUp@<Data::CurvedRidge@>'',
+ // <code>Data::SetUp@<Data::CurvedRidge@></code>,
// and everything we need to know
// about the solution would be static
// member variables and functions of
// addition by material properties,
// Neumann values, different boundary
// descriptors, etc. In that case,
- // the ``SetUp'' class might consist
+ // the <code>SetUp</code> class might consist
// of a dozen or more objects, and
// each descriptor class (like the
- // ``CurvedRidges'' class below)
+ // <code>CurvedRidges</code> class below)
// would have to provide them. Then,
// you will be happy to be able to
// change from one set of data to
// another by only changing the
- // template argument to the ``SetUp''
+ // template argument to the <code>SetUp</code>
// class at one place, rather than at
// many.
//
// obvious way, see below, with
// virtual abstract functions. It
// forces us to introduce a second
- // template parameter ``dim'' which
+ // template parameter <code>dim</code> which
// we need for the base class (which
// could be avoided using some
// template magic, but we omit that),
// simple, you don't have to touch
// the framework classes, only a
// structure like the
- // ``CurvedRidges'' one is needed.
+ // <code>CurvedRidges</code> one is needed.
namespace Data
{
// @sect4{The SetUpBase and SetUp classes}
// Based on the above description,
- // the ``SetUpBase'' class then
+ // the <code>SetUpBase</code> class then
// looks as follows. To allow using
- // the ``SmartPointer'' class with
+ // the <code>SmartPointer</code> class with
// this class, we derived from the
- // ``Subscriptor'' class.
+ // <code>Subscriptor</code> class.
template <int dim>
struct SetUpBase : public Subscriptor
{
// The class that is used to
// describe the boundary values and
- // right hand side of the ``curved
- // ridge'' problem already used in
+ // right hand side of the <code>curved
+ // ridge</code> problem already used in
// the step-13 example program is
// then like so:
template <int dim>
// why we have not chosen to
// implement the classes implementing
// a certain setup (like the
- // ``CurvedRidges'' class) directly
+ // <code>CurvedRidges</code> class) directly
// as classes derived from
- // ``Data::SetUpBase''. Indeed, we
+ // <code>Data::SetUpBase</code>. Indeed, we
// could have done very well so. The
// only reason is that then we would
// have to have member variables for
// the solution and right hand side
- // classes in the ``CurvedRidges''
+ // classes in the <code>CurvedRidges</code>
// class, as well as member functions
// overloading the abstract functions
// of the base class giving access to
// these member variables. The
- // ``SetUp'' class has the sole
+ // <code>SetUp</code> class has the sole
// reason to relieve us from the need
// to reiterate these member
// variables and functions that would
// However, there might be good
// reasons to actually implement
// classes derived from
- // ``Data::SetUpBase'', for example
+ // <code>Data::SetUpBase</code>, for example
// if the solution or right hand side
// classes require constructors that
// take arguments, which the
- // ``Data::SetUpBase'' class cannot
+ // <code>Data::SetUpBase</code> class cannot
// provide. In that case, subclassing
// is a worthwhile strategy. Other
// possibilities for special cases
// are to derive from
- // ``Data::SetUp@<SomeSetUp@>'' where
- // ``SomeSetUp'' denotes a class, or
+ // <code>Data::SetUp@<SomeSetUp@></code> where
+ // <code>SomeSetUp</code> denotes a class, or
// even to explicitly specialize
- // ``Data::SetUp@<SomeSetUp@>''. The
+ // <code>Data::SetUp@<SomeSetUp@></code>. The
// latter allows to transparently use
- // the way the ``SetUp'' class is
+ // the way the <code>SetUp</code> class is
// used for other set-ups, but with
// special actions taken for special
// arguments.
// here) was small, and the number of
// test cases was small as well. One
// then starts out by handcoding them
- // into a number of ``switch''
+ // into a number of <code>switch</code>
// statements. Over time, projects
// grow, and so does the number of
// test cases. The number of
- // ``switch'' statements grows with
+ // <code>switch</code> statements grows with
// that, and their length as well,
// and one starts to find ways to
// consider impossible examples where
// average over all cells that
// surround this point. The
// question which cells
- // ``surrounds'' the evaluation
+ // <code>surrounds</code> the evaluation
// point is made dependent on the
// mesh width by including those
// cells for which the distance of
// to zero:
rhs.reinit (dof_handler.n_dofs());
- // Initialize a ``FEValues''
+ // Initialize a <code>FEValues</code>
// object with a quadrature
// formula, have abbreviations
// for the number of quadrature
{
// If we have found such a
// cell, then initialize
- // the ``FEValues'' object
+ // the <code>FEValues</code> object
// and integrate the
// x-component of the
// gradient of each shape
// @sect4{The DualSolver class}
// In the same way as the
- // ``PrimalSolver'' class above, we
+ // <code>PrimalSolver</code> class above, we
// now implement a
- // ``DualSolver''. It has all the
+ // <code>DualSolver</code>. It has all the
// same features, the only
// difference is that it does not
// take a function object denoting
// a right hand side object, but
// now takes a
- // ``DualFunctionalBase'' object
+ // <code>DualFunctionalBase</code> object
// that will assemble the right
// hand side vector of the dual
// problem. The rest of the class
// triangulation, but different
// discretizations, it now becomes
// clear why we have made the
- // ``Base'' class a virtual one:
+ // <code>Base</code> class a virtual one:
// since the final class will be
// derived from both
- // ``PrimalSolver'' as well as
- // ``DualSolver'', it would have
- // two ``Base'' instances, would we
+ // <code>PrimalSolver</code> as well as
+ // <code>DualSolver</code>, it would have
+ // two <code>Base</code> instances, would we
// not have marked the inheritance
// as virtual. Since in many
// applications the base class
// In the private section, we
// have two functions that are
// used to call the
- // ``solve_problem'' functions
+ // <code>solve_problem</code> functions
// of the primal and dual base
// classes. These two functions
// will be called in parallel
- // by the ``solve_problem''
+ // by the <code>solve_problem</code>
// function of this class.
void solve_primal_problem ();
void solve_dual_problem ();
// error estimates on cells and
// faces, we need a number of
// helper objects, such as
- // ``FEValues'' and
- // ``FEFaceValues'' functions,
+ // <code>FEValues</code> and
+ // <code>FEFaceValues</code> functions,
// but also temporary objects
// storing the values and
// gradients of primal and dual
// In the implementation of this
// class, we first have the
- // constructors of the ``CellData''
- // and ``FaceData'' member classes,
- // and the ``WeightedResidual''
+ // constructors of the <code>CellData</code>
+ // and <code>FaceData</code> member classes,
+ // and the <code>WeightedResidual</code>
// constructor. They only
// initialize fields to their
// correct lengths, so we do not
// Now, it is becoming more
- // interesting: the ``refine_grid''
+ // interesting: the <code>refine_grid</code>
// function asks the error
// estimator to compute the
// cell-wise error indicators, then
// Since we want to output both the
// primal and the dual solution, we
- // overload the ``output_solution''
+ // overload the <code>output_solution</code>
// function. The only interesting
// feature of this function is that
// the primal and dual solutions
// are defined on different finite
// element spaces, which is not the
- // format the ``DataOut'' class
+ // format the <code>DataOut</code> class
// expects. Thus, we have to
// transfer them to a common finite
// element space. Since we want the
// primal space. For the
// interpolation, there is a
// library function, that takes a
- // ``ConstraintMatrix'' object
+ // <code>ConstraintMatrix</code> object
// including the hanging node
// constraints. The rest is
// standard.
// work-around worth mentioning: in
// this function, as in a couple of
// following ones, we have to
- // access the ``DoFHandler''
+ // access the <code>DoFHandler</code>
// objects and solutions of both
// the primal as well as of the
// dual solver. Since these are
- // members of the ``Solver'' base
+ // members of the <code>Solver</code> base
// class which exists twice in the
// class hierarchy leading to the
// present class (once as base
- // class of the ``PrimalSolver''
+ // class of the <code>PrimalSolver</code>
// class, once as base class of the
- // ``DualSolver'' class), we have
+ // <code>DualSolver</code> class), we have
// to disambiguate accesses to them
// by telling the compiler a member
// of which of these two instances
// through the class hierarchy
// which disambiguates the base
// class, for example writing
- // ``PrimalSolver::dof_handler'' to
+ // <code>PrimalSolver::dof_handler</code> to
// denote the member variable
- // ``dof_handler'' from the
- // ``Solver'' base class of the
- // ``PrimalSolver''
+ // <code>dof_handler</code> from the
+ // <code>Solver</code> base class of the
+ // <code>PrimalSolver</code>
// class. Unfortunately, this
// confuses gcc's version 2.96 (a
// version that was intended as a
// Thus, we have to work around
// this problem. We do this by
// introducing references to the
- // ``PrimalSolver'' and
- // ``DualSolver'' components of the
- // ``WeightedResidual'' object at
+ // <code>PrimalSolver</code> and
+ // <code>DualSolver</code> components of the
+ // <code>WeightedResidual</code> object at
// the beginning of the
// function. Since each of these
// has an unambiguous base class
- // ``Solver'', we can access the
+ // <code>Solver</code>, we can access the
// member variables we want through
// these references. However, we
// are now accessing protected
// member variables of these
// classes through a pointer other
- // than the ``this'' pointer (in
+ // than the <code>this</code> pointer (in
// fact, this is of course the
- // ``this'' pointer, but not
+ // <code>this</code> pointer, but not
// explicitly). This finally is the
// reason why we had to declare the
// present class a friend of the
// Add the data vectors for which
// we want output. Add them both,
- // the ``DataOut'' functions can
+ // the <code>DataOut</code> functions can
// handle as many data vectors as
// you wish to write to output:
data_out.add_data_vector (primal_solver.solution,
// the finite element space in
// which we have solved the dual
// problem: But, again as in the
- // ``WeightedResidual::output_solution''
+ // <code>WeightedResidual::output_solution</code>
// function we first need to
// create a ConstraintMatrix
// including the hanging node
// element space of the primal
// solution and subtracting it
// from z: use the
- // ``interpolate_difference''
+ // <code>interpolate_difference</code>
// function, that gives (z-I_hz)
// in the element space of the
// dual solution.
// residual contributions of
// this cell, and put them
// into the
- // ``error_indicators''
+ // <code>error_indicators</code>
// variable:
integrate_over_cell (cell, cell_index,
primal_solution,
// finite element field on the
// present cell. For this,
// initialize the
- // ``FEFaceValues'' object
+ // <code>FEFaceValues</code> object
// corresponding to this side of
// the face, and extract the
// gradients using that
// to find out with which face of
// the neighboring cell we have
// to work, i.e. the
- // ``home-many''the neighbor the
+ // <code>home-many</code>the neighbor the
// present cell is of the cell
// behind the present face. For
// this, there is a function, and
// we put the result into a
// variable with the name
- // ``neighbor_neighbor'':
+ // <code>neighbor_neighbor</code>:
const unsigned int
neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
// Then define an abbreviation
// for the neigbor cell,
// initialize the
- // ``FEFaceValues'' object on
+ // <code>FEFaceValues</code> object on
// that cell, and extract the
// gradients on that cell:
const active_cell_iterator neighbor = cell->neighbor(face_no);
// were not the case, then
// there would be either a
// bug in the
- // ``neighbor_neighbor''
+ // <code>neighbor_neighbor</code>
// function called above, or
// -- worse -- some function
// in the library did not
// values, and subsequently a
// variable of that type. It
// will default to
- // ``dual_weighted_error_estimator''.
+ // <code>dual_weighted_error_estimator</code>.
enum RefinementCriterion {
dual_weighted_error_estimator,
global_refinement,
// Next to last, a function
// that is used as a weight
// to the
- // ``RefinementWeightedKelly''
+ // <code>RefinementWeightedKelly</code>
// class. The default value
// of this pointer is zero,
// but you have to set it to
// some other value if you
// want to use the
- // ``weighted_kelly_indicator''
+ // <code>weighted_kelly_indicator</code>
// refinement criterion.
SmartPointer<const Function<dim> > kelly_weight;
descriptor.refinement_criterion
= Framework<dim>::ProblemDescription::dual_weighted_error_estimator;
// Here, we could as well have
- // used ``global_refinement''
+ // used <code>global_refinement</code>
// or
- // ``weighted_kelly_indicator''. Note
+ // <code>weighted_kelly_indicator</code>. Note
// that the information given
// about dual finite elements,
// dual functional, etc is only
// hand side. These are
// prepackaged in classes. We
// take here the description of
- // ``Exercise_2_3'', but you
+ // <code>Exercise_2_3</code>, but you
// can also use
- // ``CurvedRidges@<dim@>'':
+ // <code>CurvedRidges@<dim@></code>:
descriptor.data = new Data::SetUp<Data::Exercise_2_3<dim>,dim> ();
// Next set first a dual
// value at an
// evaluation point,
// represented by the classes
- // ``PointValueEvaluation''
+ // <code>PointValueEvaluation</code>
// in the namespaces of
// evaluation and dual
// functional classes. You can
// also set the
- // ``PointXDerivativeEvaluation''
+ // <code>PointXDerivativeEvaluation</code>
// classes for the x-derivative
// instead of the value
// at the evaluation point.
// The first thing we have here is a helper
// function that computes an even power $|v|^n$
// of a vector $v$, by evaluating
- // $(v\cdot v)^{n/2}. We need this in the
+ // $(v\cdot v)^{n/2}$. We need this in the
// computations below where we do not want to
// dwell on the fact that the gradient of the
// solution is actually a scalar in the 1d
// is obvious, note the assertion at the
// beginning of the function body, which
// makes sure that the exponent is indeed an
- // even number (here, we use that ``n/2'' is
+ // even number (here, we use that <code>n/2</code> is
// computed in integer arithmetic, i.e. any
// remainder of the division is
- // lost). ``ExcMessage'' is a pre-defined
+ // lost). <code>ExcMessage</code> is a pre-defined
// exception class that takes a string
// argument explaining what goes wrong. It is
// a simpler way to declare exceptions than
// exception class, we lose the ability to
// attach additional information at run-time
// to the exception message, such as the
- // value of the variable ``n''. By following
+ // value of the variable <code>n</code>. By following
// the way explained in above example
// programs, adding this feature is simple,
// though.
// So here comes the function that implements
- // the function object. The ``base'' value is
- // $x^{1/3}$, while ``random'' is a random
+ // the function object. The <code>base</code> value is
+ // $x^{1/3}$, while <code>random</code> is a random
// number between -1 and 1 (note that
- // ``rand()'' returns a random integer value
- // between zero and ``RAND_MAX''; to convert
+ // <code>rand()</code> returns a random integer value
+ // between zero and <code>RAND_MAX</code>; to convert
// it to a floating point value between 0 and
- // 2, we have to divide by ``RAND_MAX'' and
+ // 2, we have to divide by <code>RAND_MAX</code> and
// multiply by two -- note that the first
// multiplication has to happen in floating
// point arithmetic, so that the division is
// class. As in most of the previous example
// programs, the public interface of the
// class consists only of a constructor and a
- // ``run'' function that does the actual
+ // <code>run</code> function that does the actual
// work. The constructor takes an additional
// argument that indicates the number of the
// run we are presently performing. This
// the computations, doing one nonlinear
// step, refineming the mesh, doing a line
// search for step length computations,
- // etc. The ``energy'' function computes the
+ // etc. The <code>energy</code> function computes the
// value of the optimization functional on an
// arbitrary finite element function with
- // nodal values given on the ``DoFHandler''
+ // nodal values given on the <code>DoFHandler</code>
// given as an argument. Since it does not
// depend on the state of this object, we
- // declare this function as ``static''.
+ // declare this function as <code>static</code>.
//
// The member variables of this class are
// what we have seen before, and the
// size to the vector, and use library
// function that takes a function object,
// and interpolates the given vector living
- // on a ``DoFHandler'' to this function
+ // on a <code>DoFHandler</code> to this function
// object:
present_solution.reinit (dof_handler.n_dofs());
VectorTools::interpolate (dof_handler,
// Then we still have to make sure that we
// get the boundary values right. This
// could have been done inside the
- // ``InitializationValues'' class, but it
+ // <code>InitializationValues</code> class, but it
// is instructive to see how it can also be
// done, in particular since it is so
// simple in 1d. First, start out with an
// cell to zero. Note that the zeroth
// vertex is the left one, and that zero is
// the only valid second argument to the
- // call to ``vertex_dof_index'', since we
+ // call to <code>vertex_dof_index</code>, since we
// have a scalar finite element; thus,
// there is only a single component.
present_solution(cell->vertex_dof_index(0,0)) = 0;
matrix.reinit (sparsity_pattern);
residual.reinit (dof_handler.n_dofs());
- // Then we initialize a ``FEValues'' object
+ // Then we initialize a <code>FEValues</code> object
// with a 4-point Gauss quadrature
// formula. This object will be used to
// compute the values and gradients of the
// nonlinear step as outlined in the
// introduction to this example program. In
// order to compute values and gradients,
- // we need to pass the ``update_values''
- // and ``update_gradients'' flags to the
+ // we need to pass the <code>update_values</code>
+ // and <code>update_gradients</code> flags to the
// constructor, and the
- // ``update_JxW_values'' flag for the
+ // <code>update_JxW_values</code> flag for the
// Jacobian times the weight at a
// quadrature point. In addition, we need
// to have the coordinate values of each
// quadrature point in real space for the
// $x-u^3$ terms; to get these from the
- // ``FEValues'' object, we need to pass it
- // the ``update_q_points'' flag.
+ // <code>FEValues</code> object, we need to pass it
+ // the <code>update_q_points</code> flag.
//
// It is a simple calculation to figure out
// that for linear elements, the integrals
// therefore need to have the values and
// gradients of the previous solution at
// the quadrature points. We will get them
- // from the ``FEValues'' object above, and
+ // from the <code>FEValues</code> object above, and
// will put them into the following two
// variables:
std::vector<double> local_solution_values (n_q_points);
// the previous solution at the
// quadrature points. To get them, we
// don't actually have to do much,
- // except for giving the ``FEValues''
+ // except for giving the <code>FEValues</code>
// object the global node vector from
// which to compute this data, and a
// reference to the objects into which
// to put them. After the calls, the
- // ``local_solution_values'' and
- // ``local_solution_values'' variables
+ // <code>local_solution_values</code> and
+ // <code>local_solution_values</code> variables
// will contain values and gradients
// for each of the quadrature points on
// this cell.
// this program is ever going to be run in
// higher dimensions, then we should only
// evaluate for indicator zero, which is
- // why we have placed the ``if'' statement
+ // why we have placed the <code>if</code> statement
// in front of the second function call.
//
// Note that we need zero boundary
for (unsigned int step=0; step<5; ++step)
{
// At the present location, which is
- // ``present_solution+alpha*update'',
+ // <code>present_solution+alpha*update</code>,
// evaluate the energy
tmp = present_solution;
tmp.add (alpha, update);
const double f_a = energy (dof_handler, tmp);
// Then determine a finite difference
- // step length ``dalpha'', and also
+ // step length <code>dalpha</code>, and also
// evaluate the energy functional at
- // positions ``alpha+dalpha'' and
- // ``alpha-dalpha'' along the search
+ // positions <code>alpha+dalpha</code> and
+ // <code>alpha-dalpha</code> along the search
// direction:
const double dalpha = (alpha != 0 ? alpha/100 : 0.01);
// works in 1d. However, to make later
// extension to higher space dimensions
// simpler, we define a constant integer
- // ``dim'' at the beginning of the function;
+ // <code>dim</code> at the beginning of the function;
// by using this constant as template
// argument in all places, we are actually
// able to write most of the code as if it
// need to evaluate the gradient on the
// neighbor cells. To avoid some of the
// work needed to reinitialize a
- // ``FEValues'' object on a cell, we define
+ // <code>FEValues</code> object on a cell, we define
// another such object here that we will
// only use for the neighbor cells. The
// data we need from the side of the
// over all cells. Since we need to write
// the result for each cell into
// consecutive elements of a vector, we
- // also keep a running index ``cell_index''
+ // also keep a running index <code>cell_index</code>
// that we increase with each cell treated.
DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active (),
endc = dof_handler.end ();
for (unsigned int cell_index = 0; cell!=endc; ++cell, ++cell_index)
{
- // After initializing the ``FEValues''
+ // After initializing the <code>FEValues</code>
// object on each cell, use it to
// evaluate solution and first and
// second derivatives of it at the
// The next step is to evaluate the
// jump terms. To make computations
// somewhat simpler (and to free up the
- // ``local_*'' variables for use on
+ // <code>local_*</code> variables for use on
// neighboring elements), we define
// some convenience variables for the
// positions of the left and right cell
// actually check for this. If this
// would not be the case, an exception
// of the (predefined) class
- // ``ExcInternalError'' would be
+ // <code>ExcInternalError</code> would be
// thrown. Of course, this does not
// happen in this program, but it shows
// a way of defensive coding: if you
// the library: the quadrature classes
// do not promise any particular order
// of their quadrature points, so the
- // ``QTrapez'' class could in principle
+ // <code>QTrapez</code> class could in principle
// change the order of its two
// evaluation points. In that case,
// your code would tell you that
// really does what it is hoped to do.)
//
// Given that we are now sure that
- // ``x_left'' and ``x_right'',
+ // <code>x_left</code> and <code>x_right</code>,
// extracted from the zeroth and first
// quadrature point, are indeed the
// left and right vertex of the cell,
// we can also be sure that the values
- // we extract for ``u_left'' et al. are
+ // we extract for <code>u_left</code> et al. are
// the ones we expect them to be, since
// the order of these values must of
// course match the order of the
// cells may have totally
// independent refinement
// levels. Thus, we really need the
- // ``while'' loop, not only an
- // ``if'' clause.
+ // <code>while</code> loop, not only an
+ // <code>if</code> clause.
DoFHandler<dim>::cell_iterator left_neighbor = cell->neighbor(0);
while (left_neighbor->has_children())
left_neighbor = left_neighbor->child(1);
// With the so-found neighbor,
// initialize the second
- // ``FEValues'' object to it,
+ // <code>FEValues</code> object to it,
// extract the gradients of the
// solution there, and from this
// get the gradient at the
// interface (this is the first
- // element of ``local_gradients'',
+ // element of <code>local_gradients</code>,
// since the right end point of the
// neighbor cell has index 1) as a
// scalar value (this is the zeroth
// component of
- // ``local_gradients[1]''.
+ // <code>local_gradients[1]</code>.
neighbor_fe_values.reinit (left_neighbor);
neighbor_fe_values.get_function_grads (present_solution, local_gradients);
// examples, however, we would like to
// transfer the solution vector from the
// old to the new grid. This is what the
- // ``SolutionTransfer'' class is good for,
+ // <code>SolutionTransfer</code> class is good for,
// but it requires some preliminary
// work. First, we need to tag the cells
// that we want to refine or coarsen, as
// situations, the library will silently
// also have to refine the neighbor cell
// once. It does so by calling the
- // ``Triangulation@<dim@>::prepare_coarsening_and_refinement''
+ // <code>Triangulation@<dim@>::prepare_coarsening_and_refinement</code>
// function before actually doing the
// refinement and coarsening. This function
// flags a set of additional cells for
// this function are exactly the ones that
// will actually be refined or
// coarsened. Since the
- // ``SolutionTransfer'' class needs this
+ // <code>SolutionTransfer</code> class needs this
// information in order to store the data
// from the old mesh and transfer to the
// new one.
triangulation.prepare_coarsening_and_refinement();
// With this out of the way, we initialize
- // a ``SolutionTransfer'' object with the
- // present ``DoFHandler'' and attach the
+ // a <code>SolutionTransfer</code> object with the
+ // present <code>DoFHandler</code> and attach the
// solution vector to it:
SolutionTransfer<dim,double> solution_transfer(dof_handler);
solution_transfer.prepare_for_coarsening_and_refinement (present_solution);
// Finally, we retrieve the old solution
// interpolated to the new mesh. Since the
- // ``SolutionTransfer'' function does not
+ // <code>SolutionTransfer</code> function does not
// actually store the values of the old
// solution, but rather indices, we need to
// preserve the old solution vector until
// actually unnecessary in 1d, but
// necessary for higher space dimensions,
// so we show it anyway: the result of what
- // the ``SolutionTransfer'' class provides
+ // the <code>SolutionTransfer</code> class provides
// is a vector that is interpolated from
// the old to the new mesh. Unfortunately,
// it does not necessarily have the right
hanging_node_constraints.close ();
hanging_node_constraints.distribute (present_solution);
// This is wasteful, since we create a
- // ``ConstraintMatrix'' object that will be
+ // <code>ConstraintMatrix</code> object that will be
// recreated again in the next call to
- // ``setup_system_on_mesh'' immediately
+ // <code>setup_system_on_mesh</code> immediately
// afterwards. A more efficient
// implementation would make sure that it
// is created only once. We don't care so
// computes the energy of a nodal vector in
// the functional considered in this example
// program. Its idea is simple: take a nodal
- // vector and the ``DoFHandler'' object it is
+ // vector and the <code>DoFHandler</code> object it is
// living on, then loop over all cells and
// add up the local contributions to the
// energy:
const Vector<double> &function)
{
// First define the quadrature formula and
- // a ``FEValues'' object with which to
+ // a <code>FEValues</code> object with which to
// compute the values of the input function
// at the quadrature points. Note again
// that the integrand is a polynomial of
for (; cell!=endc; ++cell)
{
// On each cell, initialize the
- // ``FEValues'' object, and extract
+ // <code>FEValues</code> object, and extract
// values and gradients of the given
// function:
fe_values.reinit (cell);
// So here is the driver function,
- // ``run()''. It generate a coarse mesh,
+ // <code>run()</code>. It generate a coarse mesh,
// refines it a couple of times, and
// initializes the starting values. It then
// goes into a loop in which we first set up
}
- // Finally: ``main()''. This function does
+ // Finally: <code>main()</code>. This function does
// what all its counterparts in previous
// examples already did, i.e. create an
// object of the main class, and hand off
// need particularly for this example
// program and that weren't in
// step-8. First, we replace the
- // standard output ``std::cout'' by a
- // new stream ``pcout'' which is used
+ // standard output <code>std::cout</code> by a
+ // new stream <code>pcout</code> which is used
// in parallel computations for
// generating output only on one of
// the MPI processes.
// for partitioning our meshes so that they
// can be efficiently distributed across an
// MPI network. The partitioning algorithm is
- // implemented in the ``GridTools'' class,
+ // implemented in the <code>GridTools</code> class,
// and we need an additional include file for
- // a function in ``DoFRenumbering'' that
+ // a function in <code>DoFRenumbering</code> that
// allows to sort the indices associated with
// degrees of freedom so that they are
// numbered according to the subdomain they
// copied verbatim from step-8, so we only
// comment on the few things that are
// different. There is one (cosmetic) change
- // in that we let ``solve'' return a value,
+ // in that we let <code>solve</code> return a value,
// namely the number of iterations it took to
// converge, so that we can output this to
// the screen at the appropriate place. In
// addition, we introduce a stream-like
- // variable ``pcout'', explained below:
+ // variable <code>pcout</code>, explained below:
template <int dim>
class ElasticProblem
{
// to only have one process output
// everything once, for example the one
// with process number
- // zero. ``ConditionalOStream'' does
+ // zero. <code>ConditionalOStream</code> does
// exactly this: it acts as if it were a
// stream, but only forwards to a real,
// underlying stream if a flag is set. By
// setting this condition to
- // ``this_mpi_process==0'', we make sure
+ // <code>this_mpi_process==0</code>, we make sure
// that output is only generated from the
// first process and that we don't get
// the same lines of output over and over
//
// With this simple trick, we make sure
// that we don't have to guard each and
- // every write to ``std::cout'' by a
- // prefixed ``if(this_mpi_process==0)''.
+ // every write to <code>std::cout</code> by a
+ // prefixed <code>if(this_mpi_process==0)</code>.
ConditionalOStream pcout;
// The next few variables are taken
// fact that we use the parallel versions
// is denoted the fact that we use the
// classes from the
- // ``PETScWrappers::MPI'' namespace;
+ // <code>PETScWrappers::MPI</code> namespace;
// sequential versions of these classes
- // are in the ``PETScWrappers''
- // namespace, i.e. without the ``MPI''
+ // are in the <code>PETScWrappers</code>
+ // namespace, i.e. without the <code>MPI</code>
// part). Note also that we do not use a
// separate sparsity pattern, since PETSc
// manages that as part of its matrix
// computations. Note that if this is a
// sequential job without support by MPI,
// then PETSc provides some dummy type
- // for ``MPI_Comm'', so we do not have to
+ // for <code>MPI_Comm</code>, so we do not have to
// care here whether the job is really a
// parallel one:
MPI_Comm mpi_communicator;
// Then we have two variables that tell
// us where in the parallel world we
// are. The first of the following
- // variables, ``n_mpi_processes'' tells
+ // variables, <code>n_mpi_processes</code> tells
// us how many MPI processes there exist
// in total, while the second one,
- // ``this_mpi_process'', indicates which
+ // <code>this_mpi_process</code>, indicates which
// is the number of the present process
// within this space of processes. The
// latter variable will have a unique
// value for each process between zero
// and (less than)
- // ``n_mpi_processes''. If this program
+ // <code>n_mpi_processes</code>. If this program
// is run on a single machine without MPI
- // support, then their values are ``1''
- // and ``0'', respectively.
+ // support, then their values are <code>1</code>
+ // and <code>0</code>, respectively.
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
};
// system, there is one thing to do for a
// parallel program: we need to assign
// cells to each of the processes. We do
- // this by splitting (``partitioning'') the
+ // this by splitting (<code>partitioning</code>) the
// mesh cells into as many chunks
- // (``subdomains'') as there are processes
+ // (<code>subdomains</code>) as there are processes
// in this MPI job (if this is a sequential
// job, then there is only one job and all
// cells will get a zero as subdomain
// to the number of degrees of freedom),
// and also how many rows out of this
// global size are to be stored locally
- // (``n_local_dofs''). In addition, PETSc
+ // (<code>n_local_dofs</code>). In addition, PETSc
// needs to know how to partition the
// columns in the chunk of the matrix that
// is stored locally; for square matrices,
// the columns should be partitioned in the
// same way as the rows (indicated by the
- // second ``n_local_dofs'' in the call) but
+ // second <code>n_local_dofs</code> in the call) but
// in the case of rectangular matrices one
// has to partition the columns in the same
// way as vectors are partitioned with
// freedom are split in a way such that all
// DoFs in the interior of cells and between
// cells belonging to the same subdomain
- // belong to the process that ``owns'' the
+ // belong to the process that <code>owns</code> the
// cell. However, even then we sometimes need
// to assemble on a cell with a neighbor that
// belongs to a different process, and in
// this by hand, PETSc does all this for us
// by caching these elements locally, and
// sending them to the other processes as
- // necessary when we call the ``compress()''
+ // necessary when we call the <code>compress()</code>
// functions on the matrix and vector at the
// end of this function.
//
// to do this first part: instead of
// copying elements by hand into the
// global matrix, we use the
- // ``distribute_local_to_global''
+ // <code>distribute_local_to_global</code>
// functions below to take care of
// hanging nodes at the same
// time. The second step, elimination
// generality, the subdomain id is used to
// split a domain into several parts (we do
// this above, at the beginning of
- // ``setup_system''), and which allows to
+ // <code>setup_system</code>), and which allows to
// identify which subdomain a cell is
// living on. In this application, we have
// each process handle exactly one
// subdomain, so we identify the terms
- // ``subdomain'' and ``MPI process'' with
+ // <code>subdomain</code> and <code>MPI process</code> with
// each other.
//
// Apart from this, assembling the local
// whether we should also delete the
// column corresponding to a boundary
// node, or keep it (and passing
- // ``true'' as above means: yes, do
+ // <code>true</code> as above means: yes, do
// eliminate the column). If we do,
// then the resulting matrix will be
// symmetric again if it was before;
// Now all processes have computed error
// indicators for their own cells and
// stored them in the respective elements
- // of the ``local_error_per_cell''
+ // of the <code>local_error_per_cell</code>
// vector. The elements of this vector for
// cells not on the present process are
// zero. However, since all processes have
// process. They will subsequently have to
// be copied into another process's memory
// space then, an operation that PETSc does
- // for us when we call the ``compress''
+ // for us when we call the <code>compress</code>
// function. This inefficiency could be
// avoided with some more code, but we
// refrain from it since it is not a major
// Lastly, here is the driver function. It is
// almost unchanged from step-8, with the
- // exception that we replace ``std::cout'' by
- // the ``pcout'' stream. Apart from this, the
+ // exception that we replace <code>std::cout</code> by
+ // the <code>pcout</code> stream. Apart from this, the
// only other cosmetic change is that we
// output how many degrees of freedom there
// are per process, and how many iterations
}
- // So that's it, almost. ``main()'' works the
+ // So that's it, almost. <code>main()</code> works the
// same way as most of the main functions in
// the other example programs, i.e. it
- // delegates work to the ``run'' function of
+ // delegates work to the <code>run</code> function of
// a master object, and only wraps everything
// into some code to catch exceptions:
int main (int argc, char **argv)
// PETSc requires that we initialize it
// at the beginning of the program, and
// un-initialize it at the end. So we
- // call ``PetscInitialize'' and
- // ``PetscFinalize''. The original code
+ // call <code>PetscInitialize</code> and
+ // <code>PetscFinalize</code>. The original code
// sits in between, enclosed in braces
// to make sure that the
- // ``elastic_problem'' variable goes
+ // <code>elastic_problem</code> variable goes
// out of scope (and is destroyed)
// before we call
- // ``PetscFinalize''. (If we wouldn't
+ // <code>PetscFinalize</code>. (If we wouldn't
// use braces, the destructor of
- // ``elastic_problem'' would run after
- // ``PetscFinalize''; since the
+ // <code>elastic_problem</code> would run after
+ // <code>PetscFinalize</code>; since the
// destructor involves calls to PETSc
// functions, we would get strange
// error messages from PETSc.)
namespace QuasiStaticElasticity
{
- // @sect3{The ``PointHistory'' class}
+ // @sect3{The <code>PointHistory</code> class}
// As was mentioned in the introduction, we
// have to store the old stress in
// constructors, destructors, or other
// member functions. In such cases of
// `dumb' classes, we usually opt to
- // declare them as ``struct'' rather than
- // ``class'', to indicate that they are
+ // declare them as <code>struct</code> rather than
+ // <code>class</code>, to indicate that they are
// closer to C-style structures than
// C++-style classes.
template <int dim>
// need as tools. These are small
// functions that are called in
// inner loops, so we mark them as
- // ``inline''.
+ // <code>inline</code>.
//
// The first one computes the
// symmetric strain tensor for
- // shape function ``shape_func'' at
- // quadrature point ``q_point'' by
+ // shape function <code>shape_func</code> at
+ // quadrature point <code>q_point</code> by
// forming the symmetric gradient
// of this shape function. We need
// that when we want to form the
// avoided to compute any terms
// that we could prove were zero
// anyway. For this, we used the
- // ``fe.system_to_component_index''
+ // <code>fe.system_to_component_index</code>
// function that returns in which
// component a shape function was
// zero, and also that the
- // ``fe_values.shape_value'' and
- // ``fe_values.shape_grad''
+ // <code>fe_values.shape_value</code> and
+ // <code>fe_values.shape_grad</code>
// functions only returned the
// value and gradient of the single
// non-zero component of a shape
// it isn't terribly time critical,
// we can get away with a simpler
// technique: just ask the
- // ``fe_values'' for the value or
+ // <code>fe_values</code> for the value or
// gradient of a given component of
// a given shape function at a
// given quadrature point. This is
// what the
- // ``fe_values.shape_grad_component(shape_func,q_point,i)''
+ // <code>fe_values.shape_grad_component(shape_func,q_point,i)</code>
// call does: return the full
- // gradient of the ``i''th
+ // gradient of the <code>i</code>th
// component of shape function
- // ``shape_func'' at quadrature
- // point ``q_point''. If a certain
+ // <code>shape_func</code> at quadrature
+ // point <code>q_point</code>. If a certain
// component of a certain shape
// function is always zero, then
// this will simply always return
// zero.
//
// As mentioned, using
- // ``fe_values.shape_grad_component''
+ // <code>fe_values.shape_grad_component</code>
// instead of the combination of
- // ``fe.system_to_component_index''
- // and ``fe_values.shape_grad'' may
+ // <code>fe.system_to_component_index</code>
+ // and <code>fe_values.shape_grad</code> may
// be less efficient, but its
// implementation is optimized for
// such cases and shouldn't be a
// First, fill diagonal terms
// which are simply the
- // derivatives in direction ``i''
- // of the ``i'' component of the
+ // derivatives in direction <code>i</code>
+ // of the <code>i</code> component of the
// vector-valued shape
// function:
for (unsigned int i=0; i<dim; ++i)
// (here: the upper right corner)
// of the off-diagonal elements,
// and the implementation of the
- // ``SymmetricTensor'' class
+ // <code>SymmetricTensor</code> class
// makes sure that at least to
// the outside the symmetric
// entries are also filled (in
// of a vector-valued field. If you
// already have a solution field,
// the
- // ``fe_values.get_function_grads''
+ // <code>fe_values.get_function_grads</code>
// function allows you to extract
// the gradients of each component
// of your solution field at a
// by filling first the diagonal
// and then only one half of the
// symmetric tensor (the
- // ``SymmetricTensor'' class makes
+ // <code>SymmetricTensor</code> class makes
// sure that it is sufficient to
// write only one of the two
// symmetric components).
// Before we do this, though, we
// make sure that the input has the
// kind of structure we expect:
- // that is that there are ``dim''
+ // that is that there are <code>dim</code>
// vector components, i.e. one
// displacement component for each
// coordinate direction. We test
- // this with the ``Assert'' macro
+ // this with the <code>Assert</code> macro
// that will simply abort our
// program if the condition is not
// met.
//
// The reason why we stress that is that
// in this case we have that
- // ``tan_angle==0''. Further down, we
+ // <code>tan_angle==0</code>. Further down, we
// need to divide by that number in the
// computation of the axis of rotation,
// and we would get into trouble when
- // @sect3{The ``TopLevel'' class}
+ // @sect3{The <code>TopLevel</code> class}
// This is the main class of the
// program. Since the namespace already
// The external interface of the class,
// however, is unchanged: it has a public
// constructor and desctructor, and it has
- // a ``run'' function that initiated all
+ // a <code>run</code> function that initiated all
// the work.
template <int dim>
class TopLevel
std::vector<PointHistory<dim> > quadrature_point_history;
// The way this object is accessed is
- // through a ``user pointer'' that each
+ // through a <code>user pointer</code> that each
// cell, face, or edge holds: it is a
- // ``void*'' pointer that can be used
+ // <code>void*</code> pointer that can be used
// by application programs to associate
// arbitrary data to cells, faces, or
// edges. What the program actually
// for all, and instead get rid of the
// distributed copy immediately. Thus,
// note that the declaration of
- // ``inremental_displacement'' does not
+ // <code>inremental_displacement</code> does not
// denote a distribute vector as would
// be indicated by the middle namespace
- // ``MPI'':
+ // <code>MPI</code>:
PETScWrappers::MPI::SparseMatrix system_matrix;
PETScWrappers::MPI::Vector system_rhs;
// Next, how many degrees of freedom
// the present processor stores. This
// is, of course, an abbreviation to
- // ``local_dofs_per_process[this_mpi_process]''.
+ // <code>local_dofs_per_process[this_mpi_process]</code>.
unsigned int n_local_dofs;
// In the same direction, also
// necessarily contiguously
// numbered (when iterating
// over them using
- // ``active_cell_iterator'').
+ // <code>active_cell_iterator</code>).
unsigned int n_local_cells;
// Finally, we have a
};
- // @sect3{The ``BodyForce'' class}
+ // @sect3{The <code>BodyForce</code> class}
// Before we go on to the main
// functionality of this program, we have
// not electrically conducting or that
// there are no significant electromagnetic
// fields around. In that case, the body
- // forces are simply ``rho g'', where
- // ``rho'' is the material density and
- // ``g'' is a vector in negative
+ // forces are simply <code>rho g</code>, where
+ // <code>rho</code> is the material density and
+ // <code>g</code> is a vector in negative
// z-direction with magnitude 9.81 m/s^2.
- // Both the density and ``g'' are defined
+ // Both the density and <code>g</code> are defined
// in the function, and we take as the
// density 7700 kg/m^3, a value commonly
// assumed for steel.
// To be a little more general and to be
// able to do computations in 2d as well,
// we realize that the body force is always
- // a function returning a ``dim''
+ // a function returning a <code>dim</code>
// dimensional vector. We assume that
// gravity acts along the negative
- // direction of the last, i.e. ``dim-1''th
+ // direction of the last, i.e. <code>dim-1</code>th
// coordinate. The rest of the
// implementation of this function should
// be mostly self-explanatory given similar
// compiler warnings about unused function
// arguments, we therefore comment out the
// name of the first argument of the
- // ``vector_value'' function:
+ // <code>vector_value</code> function:
template <int dim>
class BodyForce : public Function<dim>
{
- // @sect3{The ``IncrementalBoundaryValue'' class}
+ // @sect3{The <code>IncrementalBoundaryValue</code> class}
// In addition to body forces, movement can
// be induced by boundary forces and forced
// rest of the boundary is either
// going to be fixed (and is then
// described using an object of
- // type ``ZeroFunction'') or free
+ // type <code>ZeroFunction</code>) or free
// (Neumann-type, in which case
// nothing special has to be done).
// The implementation of the
- // @sect3{Implementation of the ``TopLevel'' class}
+ // @sect3{Implementation of the <code>TopLevel</code> class}
// Now for the implementation of the main
// class. First, we initialize the
// constructors and descructors. There are
// no surprises here: we choose linear and
// continuous finite elements for each of
- // the ``dim'' vector components of the
+ // the <code>dim</code> vector components of the
// solution, and a Gaussian quadrature
// formula with 2 points in each coordinate
// direction. The destructor should be
// The last of the public functions is the
// one that directs all the work,
- // ``run()''. It initializes the variables
+ // <code>run()</code>. It initializes the variables
// that describe where in time we presently
// are, then runs the first time step, then
// loops over all the other time
// this happens is mostly a matter
// of taste; here, we chose to do
// it when grids are created since
- // in the ``do_initial_timestep''
- // and ``do_timestep'' functions we
+ // in the <code>do_initial_timestep</code>
+ // and <code>do_timestep</code> functions we
// want to output the number of
// cells on each processor at a
// point where we haven't called
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
hanging_node_constraints.condense (sparsity_pattern);
// Note that we have used the
- // ``CompressedSparsityPattern'' class
+ // <code>CompressedSparsityPattern</code> class
// here that was already introduced in
// step-11, rather than the
- // ``SparsityPattern'' class that we have
+ // <code>SparsityPattern</code> class that we have
// used in all other cases. The reason
// for this is that for the latter class
// to work we have to give an initial
// upper bound for the number of entries
// in each row, a task that is
// traditionally done by
- // ``DoFHandler::max_couplings_between_dofs()''. However,
+ // <code>DoFHandler::max_couplings_between_dofs()</code>. However,
// this function suffers from a serious
// problem: it has to compute an upper
// bound to the number of nonzero entries
// 3d. In effect, while it is quite
// accurate in 2d, it often comes up with
// much too large a number in 3d, and in
- // that case the ``SparsityPattern''
+ // that case the <code>SparsityPattern</code>
// allocates much too much memory at
// first, often several 100 MBs. This is
// later corrected when
- // ``DoFTools::make_sparsity_pattern'' is
+ // <code>DoFTools::make_sparsity_pattern</code> is
// called and we realize that we don't
// need all that much memory, but at time
// it is already too late: for large
// out-of-memory situations.
//
// In order to avoid this, we resort to
- // the ``CompressedSparsityPattern''
+ // the <code>CompressedSparsityPattern</code>
// class that is slower but does not
// require any up-front estimate on the
// number of nonzero entries per row. It
// After this point, no further explicit
// knowledge of the sparsity pattern is
// required any more and we can let the
- // ``sparsity_pattern'' variable go out
+ // <code>sparsity_pattern</code> variable go out
// of scope without any problem.
// The last task in this function
// stresses. In addition,
// assembling the matrix is made
// significantly more transparent
- // by using the ``SymmetricTensor''
+ // by using the <code>SymmetricTensor</code>
// class: note the elegance of
// forming the scalar products of
// symmetric tensors of rank 2 and
// symmetric gradients (strains) of
// the shape functions at a given
// quadrature point from the
- // ``FEValues'' object, and the
+ // <code>FEValues</code> object, and the
// elegance with which we form the
- // triple contraction ``eps_phi_i :
- // C : eps_phi_j''; the latter
+ // triple contraction <code>eps_phi_i :
+ // C : eps_phi_j</code>; the latter
// needs to be compared to the
// clumsy computations needed in
// step-17, both in the
// already did in previous
// programs. A slight
// complication is that the
- // ``apply_boundary_values''
+ // <code>apply_boundary_values</code>
// function wants to have a
// solution vector compatible
// with the matrix and right hand
// direction. For the boundary
// with indicator 1 (top
// surface), we use the
- // ``IncrementalBoundaryValues''
+ // <code>IncrementalBoundaryValues</code>
// class, but we specify an
// additional argument to the
- // ``VectorTools::interpolate_boundary_values''
+ // <code>VectorTools::interpolate_boundary_values</code>
// function denoting which vector
// components it should apply to;
// this is a vector of bools for
// vector and initialize it with
// the contents of the local
// variable (remember that the
- // ``apply_boundary_values''
+ // <code>apply_boundary_values</code>
// function called in
- // ``assemble_system'' preset the
+ // <code>assemble_system</code> preset the
// values of boundary nodes in this
// vector), solve with it, and at
// the end of the function copy it
// the introduction.
//
// The crucial part of this function is to
- // give the ``DataOut'' class a way to only
+ // give the <code>DataOut</code> class a way to only
// work on the cells that the present
// process owns. This class is already
// well-equipped for that: it has two
- // virtual functions ``first_cell'' and
- // ``next_cell'' that return the first cell
+ // virtual functions <code>first_cell</code> and
+ // <code>next_cell</code> that return the first cell
// to be worked on, and given one cell
// return the next cell to be worked on. By
// default, these functions return the
// first active cell (i.e. the first one
// that has no children) and the next
// active cell. What we have to do here is
- // derive a class from ``DataOut'' that
+ // derive a class from <code>DataOut</code> that
// overloads these two functions to only
// iterate over those cells with the right
// subdomain indicator.
//
// We do this at the beginning of this
- // function. The ``first_cell'' function
+ // function. The <code>first_cell</code> function
// just starts with the first active cell,
// and then iterates to the next cells
// while the cell presently under
// we don't try to keep iterating when we
// have hit the end iterator.
//
- // The ``next_cell'' function could be
+ // The <code>next_cell</code> function could be
// implemented in a similar way. However,
// we use this occasion as a pretext to
// introduce one more thing that the
// other words, it seems as if we can't
// compute the average stresses for all
// cells. However, remember that our
- // class derived from ``DataOut'' only
+ // class derived from <code>DataOut</code> only
// iterates over those cells that
// actually do belong to the present
// processor, i.e. we don't have to
// determine the name of the file
// we will want to write it
// to. We compose it of the
- // prefix ``solution-'', followed
+ // prefix <code>solution-</code>, followed
// by a representation of the
// present time written as a
// fixed point number so that
// would overflow if there were
// 1000 processes or more. Note
// that we choose to use
- // ``AssertThrow'' rather than
- // ``Assert'' since the number of
+ // <code>AssertThrow</code> rather than
+ // <code>Assert</code> since the number of
// processes is a variable that
// depends on input files or the
// way the process is started,
// in the program
// code. Therefore, it is
// inappropriate to use
- // ``Assert'' that is optimized
+ // <code>Assert</code> that is optimized
// away in optimized mode,
// whereas here we actually can
// assume that users will run the
// output to the console to update the
// person watching the screen on what is
// going on. As in step-17, the use of
- // ``pcout'' instead of ``std::cout'' makes
+ // <code>pcout</code> instead of <code>std::cout</code> makes
// sure that only one of the parallel
// processes is actually writing to the
// console, without having to explicitly
// function. First, how we get the
// displacement field at a given vertex
// using the
- // ``cell-@>vertex_dof_index(v,d)'' function
- // that returns the index of the ``d''th
- // degree of freedom at vertex ``v'' of the
+ // <code>cell-@>vertex_dof_index(v,d)</code> function
+ // that returns the index of the <code>d</code>th
+ // degree of freedom at vertex <code>v</code> of the
// given cell. In the present case,
// displacement in the k-th coordinate
// direction corresonds to the kth
// function like this bears a certain risk,
// because it uses knowledge of the order
// of elements that we have taken together
- // for this program in the ``FESystem''
+ // for this program in the <code>FESystem</code>
// element. If we decided to add an
// additional variable, for example a
// pressure variable for stabilization, and
// associated with vertices. This is indeed
// the case for the present Q1 element, as
// would be for all Qp elements of
- // polynomial order ``p''. However, it
+ // polynomial order <code>p</code>. However, it
// would not hold for discontinuous
// elements, or elements for mixed
// formulations. Secondly, it also rests on
// be. For general finite elements, the way
// to go would be to take a quadrature
// formula with the quadrature points in
- // the vertices of a cell. The ``QTrapez''
+ // the vertices of a cell. The <code>QTrapez</code>
// formula for the trapezoidal rule does
// exactly this. With this quadrature
// formula, we would then initialize an
- // ``FEValues'' object in each cell, and
+ // <code>FEValues</code> object in each cell, and
// use the
- // ``FEValues::get_function_values''
+ // <code>FEValues::get_function_values</code>
// function to obtain the values of the
// solution function in the quadrature
// points, i.e. the vertices of the
// cell. These are the only values that we
// really need, i.e. we are not at all
// interested in the weights (or the
- // ``JxW'' values) associated with this
+ // <code>JxW</code> values) associated with this
// particular quadrature formula, and this
// can be specified as the last argument in
- // the constructor to ``FEValues''. The
+ // the constructor to <code>FEValues</code>. The
// only point of minor inconvenience in
// this scheme is that we have to figure
// out which quadrature point corresponds
// this short function is the way in which
// the triangulation class exports
// information about its vertices: through
- // the ``Triangulation::n_vertices''
+ // the <code>Triangulation::n_vertices</code>
// function, it advertises how many
// vertices there are in the
// triangulation. Not all of them are
// the number of a vertex once it has come
// into existence, even if vertices with
// lower number go away. Secondly, the
- // location returned by ``cell-@>vertex(v)''
+ // location returned by <code>cell-@>vertex(v)</code>
// is not only a read-only object of type
- // ``Point@<dim@>'', but in fact a reference
+ // <code>Point@<dim@></code>, but in fact a reference
// that can be written to. This allows to
// move around the nodes of a mesh with
// relative ease, but it is worth pointing
// history variables, such as the existing
// stresses in the material, that we store
// in each quadrature point. As mentioned
- // above, we use the ``user_pointer'' for
+ // above, we use the <code>user_pointer</code> for
// this that is available in each cell.
//
// To put this into larger perspective, we
// Next, allocate as many quadrature
// objects as we need. Since the
- // ``resize'' function does not actually
+ // <code>resize</code> function does not actually
// shrink the amount of allocated memory
// if the requested new size is smaller
// than the old size, we resort to a
// and then swap the contents of the old
// vector and this temporary
// variable. This makes sure that the
- // ``quadrature_point_history'' is now
+ // <code>quadrature_point_history</code> is now
// really empty, and we can let the
// temporary variable that now holds the
// previous contents of the vector go out
// step. we can then re-allocate as many
// elements as we need, with the vector
// default-initializing the
- // ``PointHistory'' objects, which
+ // <code>PointHistory</code> objects, which
// includes setting the stress variables
// to zero.
{
// this function that forget to update
// all uses of a variable at the same
// time. Recall that constructs using the
- // ``Assert'' macro are optimized away in
+ // <code>Assert</code> macro are optimized away in
// optimized mode, so do not affect the
// run time of optimized runs:
Assert (history_index == quadrature_point_history.size(),
template <int dim>
void TopLevel<dim>::update_quadrature_point_history ()
{
- // First, set up an ``FEValues'' object
+ // First, set up an <code>FEValues</code> object
// by which we will evaluate the
// incremental displacements and the
// gradients thereof at the quadrature
&quadrature_point_history.back(),
ExcInternalError());
- // Then initialize the ``FEValues''
+ // Then initialize the <code>FEValues</code>
// object on the present cell, and
// extract the gradients of the
// displacement at the quadrature
// tensor by contraction from
// the left and right, after we
// expand the symmetric tensor
- // ``new_stress'' into a full
+ // <code>new_stress</code> into a full
// tensor:
const SymmetricTensor<2,dim> rotated_new_stress
= symmetrize(transpose(rotation) *
// result. When
// assigning the result
// to a
- // ``SymmetricTensor'',
+ // <code>SymmetricTensor</code>,
// the constuctor of
// that class checks
// the symmetry and
// This ends the project specific
// namespace
- // ``QuasiStaticElasticity''. The
+ // <code>QuasiStaticElasticity</code>. The
// rest is as usual and as already
- // shown in step-17: A ``main()''
+ // shown in step-17: A <code>main()</code>
// function that initializes and
// terminates PETSc, calls the
// classes that do the actual work,
// prints a general message, and then goes on
// to list the parameters that are allowed in
// the parameter file (the
- // ``ParameterHandler'' class has a function
+ // <code>ParameterHandler</code> class has a function
// to do exactly this; see the results
// section for what it prints):
void
// format, we nevertheless want to show how
// to work with parameter files.
//
- // In short, the ``ParameterHandler'' class
+ // In short, the <code>ParameterHandler</code> class
// works as follows: one declares the entries
// of parameters that can be given in input
// files together, and later on one can read
// value specified in the declaration of that
// parameter is used. After that, the program
// can query the values assigned to certain
- // parameters from the ``ParameterHandler''
+ // parameters from the <code>ParameterHandler</code>
// object.
//
// Declaring parameters can be done using the
- // ``ParameterHandler::declare_entry''
+ // <code>ParameterHandler::declare_entry</code>
// function. It's arguments are the name of a
// parameter, a default value (given as a
// string, even if the parameter is numeric
// describes constraints on values that may
// be passed to this parameter. In the
// example below, we use an object of type
- // ``Patterns::Anything'' to denote that
+ // <code>Patterns::Anything</code> to denote that
// there are no constraints on file names
// (this is, of course, not true -- the
// operating system does have constraints,
// but from an application standpoint, almost
// all names are valid). In other cases, one
// may, for example, use
- // ``Patterns::Integer'' to make sure that
+ // <code>Patterns::Integer</code> to make sure that
// only parameters are accepted that can be
// interpreted as integer values (it is also
// possible to specify bounds for integer
// values, and all values outside this range
- // are rejected), ``Patterns::Double'' for
+ // are rejected), <code>Patterns::Double</code> for
// floating point values, classes that make
// sure that the given parameter value is a
// comma separated list of things, etc. Take
- // a look at the ``Patterns'' namespace to
+ // a look at the <code>Patterns</code> namespace to
// see what is possible.
//
- // The fourth argument to ``declare_entry''
+ // The fourth argument to <code>declare_entry</code>
// is a help string that can be printed to
// document what this parameter is meant to
// be used for and other information you may
// fourth argument is the empty string.
//
// I always wanted to have an example program
- // describing the ``ParameterHandler'' class,
+ // describing the <code>ParameterHandler</code> class,
// because it is so particularly useful. It
// would have been useful in a number of
// previous example programs (for example, in
// out: declaring and querying these
// parameters was already done centralized in
// one place of the libray, namely the
- // ``DataOutInterface'' class that handles
+ // <code>DataOutInterface</code> class that handles
// exactly this -- managing parameters for
// input and output.
//
// So the second function call in this
// function is to let the
- // ``DataOutInterface'' declare a good number
+ // <code>DataOutInterface</code> declare a good number
// of parameters that control everything from
// the output format to what kind of output
// should be generated if output is written
// options up front, when output is
// generated, rather than playing around with
// them later on. The call to
- // ``DataOutInterface::declare_parameters''
+ // <code>DataOutInterface::declare_parameters</code>
// declares entries that allow to specify
// them in the parameter input file during
// run-time. If the parameter file does not
// contain entries for them, defaults are
// taken.
//
- // As a final note: ``DataOutInterface'' is a
+ // As a final note: <code>DataOutInterface</code> is a
// template, because it is usually used to
// write output for a specific space
// dimension. However, this program is
// parameter. Fortunately, declaring
// parameters is something that is space
// dimension independent, so we can just pick
- // one arbitrarily. We pick ``1'', but it
+ // one arbitrarily. We pick <code>1</code>, but it
// could have been any other number as well.
void declare_parameters ()
{
// the list of input files can't, so at least
// one parameter needs to be there. Together
// with the name of the program (the zeroth
- // parameter), ``argc'' must therefore be at
+ // parameter), <code>argc</code> must therefore be at
// least 2. If this is not the case, we print
// an error message and exit:
void
// Next, collect all parameters in a list
// that will be somewhat simpler to handle
- // than the ``argc''/``argv'' mechanism. We
+ // than the <code>argc</code>/<code>argv</code> mechanism. We
// omit the name of the executable at the
// zeroth index:
std::list<std::string> args;
// Then process all these
// parameters. If the parameter is
- // ``-p'', then there must be a
+ // <code>-p</code>, then there must be a
// parameter file following (which
// we should then read), in case of
- // ``-x'' it is the name of an
+ // <code>-x</code> it is the name of an
// output format. Finally, for
- // ``-o'' it is the name of the
+ // <code>-o</code> it is the name of the
// output file. In all cases, once
// we've treated a parameter, we
// remove it from the list of
// the dummy subsection, we would
// write something like this to
// extract the value of the boolean
- // flag (the ``prm.get'' function
+ // flag (the <code>prm.get</code> function
// returns the value of a parameter
// as a string, whereas the
- // ``prm.get_X'' functions return a
+ // <code>prm.get_X</code> functions return a
// value already converted to a
// different type):
prm.enter_subsection ("Dummy subsection");
// from all the input file, and read in the
// first file through a stream. Note that
// every time we open a file, we use the
- // ``AssertThrow'' macro to check whether the
+ // <code>AssertThrow</code> macro to check whether the
// file is really readable -- if it isn't
// then this will trigger an exception and
// corresponding output will be generated
- // from the exception handler in ``main()'':
+ // from the exception handler in <code>main()</code>:
template <int dim, int spacedim>
void do_convert ()
{
// stream, and parse what we got as the
// name of the output format into an
// identifier. Fortunately, the
- // ``DataOutBase'' class has a function
+ // <code>DataOutBase</code> class has a function
// that does this parsing for us, i.e. it
// knows about all the presently supported
// output formats and makes sure that they
//
// So here is what we do: from the first
// input file, we determine (using a function
- // in ``DataOutBase'' that exists for this
+ // in <code>DataOutBase</code> that exists for this
// purpose) these dimensions. We then have a
// series of switches that dispatch,
- // statically, to the ``do_convert''
+ // statically, to the <code>do_convert</code>
// functions with different template
// arguments. Not pretty, but works. Apart
// from this, the function does nothing --
// except making sure that it covered the
// dimensions for which it was called, using
- // the ``AssertThrow'' macro at places in the
+ // the <code>AssertThrow</code> macro at places in the
// code that shouldn't be reached:
void convert ()
{
// field on the triangulation.
//
// This function shows how to do this. The
- // object to consider is the ``DoFHandler''
+ // object to consider is the <code>DoFHandler</code>
// class template. Before we do so, however,
// we first need something that describes how
// many degrees of freedom are to be
// finite element space, the finite element
// base class stores this information. In the
// present context, we therefore create an
- // object of the derived class ``FE_Q'' that
+ // object of the derived class <code>FE_Q</code> that
// describes Lagrange elements. Its
// constructor takes one argument that states
// the polynomial degree of the element,
// given to the constructor would instead
// give us a bi-cubic element with one degree
// of freedom per vertex, two per line, and
- // four inside the cell. In general, ``FE_Q''
+ // four inside the cell. In general, <code>FE_Q</code>
// denotes the family of continuous elements
// with complete polynomials
// (i.e. tensor-product polynomials) up to
//
// We first need to create an object of this
// class and then pass it on to the
- // ``DoFHandler'' object to allocate storage
+ // <code>DoFHandler</code> object to allocate storage
// for the degrees of freedom (in deal.II
- // lingo: we ``distribute degrees of
- // freedom''). Note that the DoFHandler
+ // lingo: we <code>distribute degrees of
+ // freedom</code>). Note that the DoFHandler
// object will store a reference to this
// finite element object, so we need have to
// make sure its lifetime is at least as long
- // as that of the ``DoFHandler''; one way to
+ // as that of the <code>DoFHandler</code>; one way to
// make sure this is so is to make it static
// as well, in order to prevent its
// preemptive destruction. (However, the
sparsity_pattern.compress ();
// What actually happens in this call is
// the following: upon creation of a
- // ``SparsityPattern'' object, memory is
+ // <code>SparsityPattern</code> object, memory is
// allocated for a maximum number of
// entries per row (20 in our case). The
// call to
- // ``DoFTools::make_sparsity_pattern'' then
+ // <code>DoFTools::make_sparsity_pattern</code> then
// actually allocates entries as necessary,
// but the number of nonzero entries in any
// given row may be less than the 20 we
// have allocated memory for. To save
// memory and to simplify some other
// operations, one then needs to
- // ``compress'' the sparsity pattern before
+ // <code>compress</code> the sparsity pattern before
// anything else.
// With this, we can now write the results
// the sparsity pattern is symmetric. This
// should not come as a surprise, since we
// have not given the
- // ``DoFTools::make_sparsity_pattern'' any
+ // <code>DoFTools::make_sparsity_pattern</code> any
// information that would indicate that our
// bilinear form may couple shape functions
// in a non-symmetric way. You will also
// sparsity pattern is more localized around
// the diagonal. The only interesting part of
// the function is the first call to
- // ``DoFRenumbering::Cuthill_McKee'', the
+ // <code>DoFRenumbering::Cuthill_McKee</code>, the
// rest is essentially as before:
void renumber_dofs (DoFHandler<2> &dof_handler)
{
// 100,000s).
// It is worth noting that the
- // ``DoFRenumbering'' class offers a number
+ // <code>DoFRenumbering</code> class offers a number
// of other algorithms as well to renumber
// degrees of freedom. For example, it would
// of course be ideal if all couplings were
// possible by enumerating degrees of freedom
// from the inflow boundary along streamlines
// to the outflow boundary. Not surprisingly,
- // ``DoFRenumbering'' also has algorithms for
+ // <code>DoFRenumbering</code> also has algorithms for
// this.
// Finally, this is the main program. The
// only thing it does is to allocate and
// create the triangulation, then create a
- // ``DoFHandler'' object and associate it to
+ // <code>DoFHandler</code> object and associate it to
// the triangulation, and finally call above
// two functions on it:
int main ()
// spatial dependence, we consider it
// a tensor-valued function. The
// following include file provides
- // the ``TensorFunction'' class that
+ // the <code>TensorFunction</code> class that
// offers such functionality:
#include <base/tensor_function.h>
- // @sect3{The ``MixedLaplaceProblem'' class template}
+ // @sect3{The <code>MixedLaplaceProblem</code> class template}
// Again, since this is an adaptation
// of step-6, the main class is
// argument (and that there is a
// corresponding member variable to
// store this value) and the addition
- // of the ``compute_error'' function
+ // of the <code>compute_error</code> function
// in which, no surprise, we will
// compute the difference between the
// exact and the numerical solution
// exact solution for later
// computations of the error. Note
// that these functions have one,
- // one, and ``dim+1'' components,
+ // one, and <code>dim+1</code> components,
// respectively, and that we pass the
// number of components down to the
- // ``Function@<dim@>'' base class. For
+ // <code>Function@<dim@></code> base class. For
// the exact solution, we only
// declare the function that actually
// returns the entire solution vector
// because this is all that appears
// in the weak form -- the inverse of
// the permeability tensor,
- // ``KInverse''. For the purpose of
+ // <code>KInverse</code>. For the purpose of
// verifying the exactness of the
// solution and determining
// convergence orders, this tensor is
// Possibly unsurprising, deal.II
// also has a base class not only for
// scalar and generally vector-valued
- // functions (the ``Function'' base
+ // functions (the <code>Function</code> base
// class) but also for functions that
// return tensors of fixed dimension
- // and rank, the ``TensorFunction''
+ // and rank, the <code>TensorFunction</code>
// template. Here, the function under
// consideration returns a dim-by-dim
// matrix, i.e. a tensor of rank 2
- // and dimension ``dim''. We then
+ // and dimension <code>dim</code>. We then
// choose the template arguments of
// the base class appropriately.
//
// The interface that the
- // ``TensorFunction'' class provides
+ // <code>TensorFunction</code> class provides
// is essentially equivalent to the
- // ``Function'' class. In particular,
- // there exists a ``value_list''
+ // <code>Function</code> class. In particular,
+ // there exists a <code>value_list</code>
// function that takes a list of
// points at which to evaluate the
// function, and returns the values
//
// The only thing worth describing
// here is the constructor call of
- // the ``fe'' variable. The
- // ``FESystem'' class to which this
+ // the <code>fe</code> variable. The
+ // <code>FESystem</code> class to which this
// variable belongs has a number of
// different constructors that all
// refer to binding simpler elements
// want to couple a single RT(degree)
// element with a single DQ(degree)
// element. The constructor to
- // ``FESystem'' that does this
+ // <code>FESystem</code> that does this
// requires us to specity first the
// first base element (the
- // ``FE_RaviartThomas'' object of
+ // <code>FE_RaviartThomas</code> object of
// given degree) and then the number
// of copies for this base element,
// and then similarly the kind and
- // number of ``FE_DGQ''
+ // number of <code>FE_DGQ</code>
// elements. Note that the Raviart
- // Thomas element already has ``dim''
+ // Thomas element already has <code>dim</code>
// vector components, so that the
// coupled element will have
- // ``dim+1'' vector components, the
- // first ``dim'' of which correspond
+ // <code>dim+1</code> vector components, the
+ // first <code>dim</code> of which correspond
// to the velocity variable whereas the
// last one corresponds to the
// pressure.
// we constructed this element from
// its base elements, with the way we
// have done so in step-8: there, we
- // have built it as ``fe
- // (FE_Q@<dim@>(1), dim)'', i.e. we
- // have simply used ``dim'' copies of
- // the ``FE_Q(1)'' element, one copy
+ // have built it as <code>fe
+ // (FE_Q@<dim@>(1), dim)</code>, i.e. we
+ // have simply used <code>dim</code> copies of
+ // the <code>FE_Q(1)</code> element, one copy
// for the displacement in each
// coordinate direction.
template <int dim>
// blocks, so that we can allocate
// an appropriate amount of
// space. To this end, we call the
- // ``DoFTools::count_dofs_per_component''
+ // <code>DoFTools::count_dofs_per_component</code>
// function that counts how many
// shape functions are non-zero for
// a particular vector
- // component. We have ``dim+1''
+ // component. We have <code>dim+1</code>
// vector components, and we have
// to use the knowledge that for
// Raviart-Thomas elements all
// step, we allocate a 2x2 block
// pattern and then reinitialize
// each of the blocks to its
- // correct size using the ``n_u''
- // and ``n_p'' variables defined
+ // correct size using the <code>n_u</code>
+ // and <code>n_p</code> variables defined
// above that hold the number of
// velocity and pressure
// variables. In this second step,
// its knowledge about the sizes of
// the blocks it manages; this
// happens with the
- // ``sparsity_pattern.collect_sizes()''
+ // <code>sparsity_pattern.collect_sizes()</code>
// call:
const unsigned int
n_couplings = dof_handler.max_couplings_between_dofs();
// are all the usual steps, with the
// addition that we do not only
// allocate quadrature and
- // ``FEValues'' objects for the cell
+ // <code>FEValues</code> objects for the cell
// terms, but also for face
// terms. After that, we define the
// usual abbreviations for variables,
// rather only comment on
// implementational aspects.
- // @sect4{The ``InverseMatrix'' class template}
+ // @sect4{The <code>InverseMatrix</code> class template}
// The first component of our linear
// solver scheme was the creation of
// a class that acts like the inverse
// of a matrix, i.e. which has a
- // ``vmult'' function that multiplies
+ // <code>vmult</code> function that multiplies
// a vector with an inverse matrix by
// solving a linear system.
//
// purpose of this class, two
// comments are in order. First, the
// class is derived from the
- // ``Subscriptor'' class so that we
- // can use the ``SmartPointer'' class
+ // <code>Subscriptor</code> class so that we
+ // can use the <code>SmartPointer</code> class
// with inverse matrix objects. The
- // use of the ``Subscriptor'' class
+ // use of the <code>Subscriptor</code> class
// has been explained before in
// step-7 and step-20. The present
// class also sits on the receiving
// end of this
- // ``Subscriptor''/``SmartPointer''
+ // <code>Subscriptor</code>/<code>SmartPointer</code>
// pair: it holds its pointer to the
// matrix it is supposed to be the
// inverse of through a
- // ``SmartPointer'' to make sure that
+ // <code>SmartPointer</code> to make sure that
// this matrix is not destroyed while
// we still have a pointer to it.
//
// vectors that it will release again
// at the end of its operation. What
// this means is that through
- // repeated calls to the ``vmult''
+ // repeated calls to the <code>vmult</code>
// function of this class we have to
// allocate and release vectors over
// and over again.
// only once? In fact, deal.II offers
// a way to do exactly this. What all
// the linear solvers do is not to
- // allocate memory using ``new'' and
- // ``delete'', but rather to allocate
+ // allocate memory using <code>new</code> and
+ // <code>delete</code>, but rather to allocate
// them from an object derived from
- // the ``VectorMemory'' class (see
+ // the <code>VectorMemory</code> class (see
// the module on Vector memory
// management in the API reference
// manual). By default, the linear
// solvers use a derived class
- // ``PrimitiveVectorMemory'' that,
+ // <code>PrimitiveVectorMemory</code> that,
// ever time a vector is requested,
- // allocates one using ``new'', and
- // calls ``delete'' on it again once
+ // allocates one using <code>new</code>, and
+ // calls <code>delete</code> on it again once
// the solver returns it to the
- // ``PrimitiveVectorMemory''
+ // <code>PrimitiveVectorMemory</code>
// object. This is the appropriate
// thing to do if we do not
// anticipate that the vectors may be
// vector memory object holds on to
// them for later requests by linear
// solvers. The
- // ``GrowingVectorMemory'' class does
+ // <code>GrowingVectorMemory</code> class does
// exactly this: when asked by a
// linear solver for a vector, it
// first looks whether it has unused
// simply grows its pool. Vectors are
// only returned to the C++ runtime
// memory system once the
- // ``GrowingVectorMemory'' object is
+ // <code>GrowingVectorMemory</code> object is
// destroyed itself.
//
// What we therefore need to do is
// have the present matrix have an
// object of type
- // ``GrowingVectorMemory'' as a
+ // <code>GrowingVectorMemory</code> as a
// member variable and use it
// whenever we create a linear solver
// object. There is a slight
// complication here: Since the
- // ``vmult'' function is marked as
- // ``const'' (it doesn't change the
+ // <code>vmult</code> function is marked as
+ // <code>const</code> (it doesn't change the
// state of the object, after all,
// and simply operates on its
// arguments), it can only pass an
// such attempt as an error, if we
// didn't make use of a rarely used
// feature of C++: we mark the
- // variable as ``mutable''. What this
+ // variable as <code>mutable</code>. What this
// does is to allow us to change a
// member variable even from a
- // ``const'' member function.
+ // <code>const</code> member function.
template <class Matrix>
class InverseMatrix : public Subscriptor
{
}
- // @sect4{The ``SchurComplement'' class template}
+ // @sect4{The <code>SchurComplement</code> class template}
// The next class is the Schur
// complement class. Its rationale
// in the introduction. The only
// things we would like to note is
// that the class, too, is derived
- // from the ``Subscriptor'' class and
+ // from the <code>Subscriptor</code> class and
// that as mentioned above it stores
// pointers to the entire block
// matrix and the inverse of the mass
// matrix block using
- // ``SmartPointer'' objects.
+ // <code>SmartPointer</code> objects.
//
- // The ``vmult'' function requires
+ // The <code>vmult</code> function requires
// two temporary vectors that we do
// not want to re-allocate and free
// every time we call this
// function. Since here, we have full
// control over the use of these
// vectors (unlike above, where a
- // class called by the ``vmult''
+ // class called by the <code>vmult</code>
// function required these vectors,
- // not the ``vmult'' function
+ // not the <code>vmult</code> function
// itself), we allocate them
// directly, rather than going
- // through the ``VectorMemory''
+ // through the <code>VectorMemory</code>
// mechanism. However, again, these
// member variables do not carry any
// state between successive calls to
// (i.e., we never care what values
// they were set to the last time a
// member function was called), we
- // mark these vectors as ``mutable''.
+ // mark these vectors as <code>mutable</code>.
//
// The rest of the (short)
// implementation of this class is
// straightforward if you know the
// order of matrix-vector
// multiplications performed by the
- // ``vmult'' function:
+ // <code>vmult</code> function:
class SchurComplement : public Subscriptor
{
public:
}
- // @sect4{The ``ApproximateSchurComplement'' class template}
+ // @sect4{The <code>ApproximateSchurComplement</code> class template}
// The third component of our solver
// and preconditioner system is the
// class that approximates the Schur
// complement so we can form a
- // ``InverseMatrix@<ApproximateSchurComplement@>''
+ // <code>InverseMatrix@<ApproximateSchurComplement@></code>
// object that approximates the
// inverse of the Schur
// complement. It follows the same
// pattern as the Schur complement
// class, with the only exception
// that we do not multiply with the
- // inverse mass matrix in ``vmult'',
+ // inverse mass matrix in <code>vmult</code>,
// but rather just do a single Jacobi
// step. Consequently, the class also
// does not have to store a pointer
//
// To compute errors in the solution,
// we have already introduced the
- // ``VectorTools::integrate_difference''
+ // <code>VectorTools::integrate_difference</code>
// function in step-7 and
// step-11. However, there we only
// dealt with scalar solutions,
// have to do is to `mask' the
// components that we are interested
// in. This is easily done: the
- // ``VectorTools::integrate_difference''
+ // <code>VectorTools::integrate_difference</code>
// function takes as its last
// argument a pointer to a weight
// function (the parameter defaults
// should pass a function that
// represents the constant vector
// with a unit value in component
- // ``dim'', whereas for the velocity
+ // <code>dim</code>, whereas for the velocity
// the constant vector should be one
- // in the first ``dim'' components,
+ // in the first <code>dim</code> components,
// and zero in the location of the
// pressure.
//
// In deal.II, the
- // ``ComponentSelectFunction'' does
+ // <code>ComponentSelectFunction</code> does
// exactly this: it wants to know how
// many vector components the
// function it is to represent should
// have (in our case this would be
- // ``dim+1'', for the joint
+ // <code>dim+1</code>, for the joint
// velocity-pressure space) and which
// individual or range of components
// should be equal to one. We
// and a vector in which we will
// store the cellwise errors as
// computed by
- // ``integrate_difference'':
+ // <code>integrate_difference</code>:
template <int dim>
void MixedLaplaceProblem<dim>::compute_errors () const
{
// quadrature. This actually
// presents a slight twist here: if
// we naively chose an object of
- // type ``QGauss@<dim@>(degree+1)''
+ // type <code>QGauss@<dim@>(degree+1)</code>
// as one may be inclined to do
// (this is what we used for
// integrating the linear system),
// ingration. To avoid this
// problem, we simply use a
// trapezoidal rule and iterate it
- // ``degree+2'' times in each
+ // <code>degree+2</code> times in each
// coordinate direction (again as
// explained in step-7):
QTrapez<1> q_trapez;
}
- // @sect3{The ``main'' function}
+ // @sect3{The <code>main</code> function}
// The main function we stole from
// step-6 instead of step-4. It is
// refinement indicator.
#include <numerics/derivative_approximation.h>
// Finally we do some time comparison
- // using the ``Timer'' class.
+ // using the <code>Timer</code> class.
#include <base/timer.h>
// And this again is C++:
//
// First we define the classes
// representing the equation-specific
- // functions. Both classes, ``RHS''
- // and ``BoundaryValues'', are
- // derived from the ``Function''
- // class. Only the ``value_list''
+ // functions. Both classes, <code>RHS</code>
+ // and <code>BoundaryValues</code>, are
+ // derived from the <code>Function</code>
+ // class. Only the <code>value_list</code>
// function are implemented because
// only lists of function values are
// computed rather than single
};
- // The class ``Beta'' represents the
+ // The class <code>Beta</code> represents the
// vector valued flow field of the
// linear transport equation and is
- // not derived from the ``Function''
+ // not derived from the <code>Function</code>
// class as we prefer to get function
- // values of type ``Point'' rather
+ // values of type <code>Point</code> rather
// than of type
- // ``Vector@<double@>''. This, because
+ // <code>Vector@<double@></code>. This, because
// there exist scalar products
- // between ``Point'' and ``Point'' as
- // well as between ``Point'' and
- // ``Tensor'', simplifying terms like
+ // between <code>Point</code> and <code>Point</code> as
+ // well as between <code>Point</code> and
+ // <code>Tensor</code>, simplifying terms like
// $\beta\cdot n$ and
// $\beta\cdot\nabla v$.
//
// The implementation of the
- // ``value_list'' functions of these
+ // <code>value_list</code> functions of these
// classes are rather simple. For
// simplicity the right hand side is
// set to be zero but will be
// Next we define the
// equation-dependent and
// DG-method-dependent class
- // ``DGTransportEquation''. Its
+ // <code>DGTransportEquation</code>. Its
// member functions were already
// mentioned in the Introduction and
// will be explained
// below. Furthermore it includes
// objects of the previously defined
- // ``Beta'', ``RHS'' and
- // ``BoundaryValues'' function
+ // <code>Beta</code>, <code>RHS</code> and
+ // <code>BoundaryValues</code> function
// classes.
template <int dim>
class DGTransportEquation
// @sect4{Function: assemble_cell_term}
//
- // The ``assemble_cell_term''
+ // The <code>assemble_cell_term</code>
// function assembles the cell terms
// of the discretization.
- // ``u_v_matrix'' is a cell matrix,
+ // <code>u_v_matrix</code> is a cell matrix,
// i.e. for a DG method of degree 1,
// it is of size 4 times 4, and
- // ``cell_vector'' is of size 4.
+ // <code>cell_vector</code> is of size 4.
// When this function is invoked,
- // ``fe_v'' is already reinit'ed with the
+ // <code>fe_v</code> is already reinit'ed with the
// current cell before and includes
// all shape values needed.
template <int dim>
FullMatrix<double> &u_v_matrix,
Vector<double> &cell_vector) const
{
- // First we ask ``fe_v'' for the
+ // First we ask <code>fe_v</code> for the
// shape gradients, shape values and
// quadrature weights,
const std::vector<double> &JxW = fe_v.get_JxW_values ();
// Then the flow field beta and the
- // ``rhs_function'' are evaluated at
+ // <code>rhs_function</code> are evaluated at
// the quadrature points,
std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
std::vector<double> rhs (fe_v.n_quadrature_points);
// @sect4{Function: assemble_boundary_term}
//
- // The ``assemble_boundary_term''
+ // The <code>assemble_boundary_term</code>
// function assembles the face terms
// at boundary faces. When this
- // function is invoked, ``fe_v'' is
+ // function is invoked, <code>fe_v</code> is
// already reinit'ed with the current
// cell and current face. Hence it
// provides the shape values on that
Vector<double> &cell_vector) const
{
// Again, as in the previous
- // function, we ask the ``FEValues''
+ // function, we ask the <code>FEValues</code>
// object for the shape values and
// the quadrature weights
const std::vector<double> &JxW = fe_v.get_JxW_values ();
// @sect4{Function: assemble_face_term1}
//
- // The ``assemble_face_term1''
+ // The <code>assemble_face_term1</code>
// function assembles the face terms
// corresponding to the first version
// of the DG method, cf. above. For
// all cell boundaries.
//
// When this function is invoked,
- // ``fe_v'' and ``fe_v_neighbor'' are
+ // <code>fe_v</code> and <code>fe_v_neighbor</code> are
// already reinit'ed with the current
// cell and the neighoring cell,
// respectively, as well as with the
// on the face.
//
// In addition to the cell matrix
- // ``u_v_matrix'' this function has
+ // <code>u_v_matrix</code> this function has
// got a new argument
- // ``un_v_matrix'', that stores
+ // <code>un_v_matrix</code>, that stores
// contributions to the system matrix
// that are based on outer values of
// u, see $\hat u_h$ in the
// introduction, and inner values of
// v, see $v_h$. Here we note that
- // ``un'' is the short notation for
- // ``u_neighbor'' and represents
+ // <code>un</code> is the short notation for
+ // <code>u_neighbor</code> and represents
// $\hat u_h$.
template <int dim>
void DGTransportEquation<dim>::assemble_face_term1(
// @sect4{Function: assemble_face_term2}
//
// Now we look at the
- // ``assemble_face_term2'' function
+ // <code>assemble_face_term2</code> function
// that assembles the face terms
// corresponding to the second
// version of the DG method,
// terms are given as a sum of
// integrals over all faces. Here we
// need two additional cell matrices
- // ``u_vn_matrix'' and
- // ``un_vn_matrix'' that will store
+ // <code>u_vn_matrix</code> and
+ // <code>un_vn_matrix</code> that will store
// contributions due to terms
// involving u and vn as well as un
// and vn.
// After these preparations, we
// proceed with the main part of this
// program. The main class, here
- // called ``DGMethod'' is basically
+ // called <code>DGMethod</code> is basically
// the main class of step 6. One of
// the differences is that there's no
// ConstraintMatrix object. This is,
// solutions to the problems
// corresponding to the two
// different assembling routines
- // ``assemble_system1'' and
- // ``assemble_system2'';
+ // <code>assemble_system1</code> and
+ // <code>assemble_system2</code>;
Vector<double> solution1;
Vector<double> solution2;
Vector<double> right_hand_side;
// @sect4{Function: assemble_system1}
//
// We proceed with the
- // ``assemble_system1'' function that
+ // <code>assemble_system1</code> function that
// implements the DG discretization
// in its first version. This
// function repeatedly calls the
- // ``assemble_cell_term'',
- // ``assemble_boundary_term'' and
- // ``assemble_face_term1'' functions
- // of the ``DGTransportEquation''
+ // <code>assemble_cell_term</code>,
+ // <code>assemble_boundary_term</code> and
+ // <code>assemble_face_term1</code> functions
+ // of the <code>DGTransportEquation</code>
// object. The
- // ``assemble_boundary_term'' covers
+ // <code>assemble_boundary_term</code> covers
// the first case mentioned in the
// introduction.
//
// 1. face is at boundary
//
// This function takes a
- // ``FEFaceValues'' object as
+ // <code>FEFaceValues</code> object as
// argument. In contrast to that
- // ``assemble_face_term1''
- // takes two ``FEFaceValuesBase''
+ // <code>assemble_face_term1</code>
+ // takes two <code>FEFaceValuesBase</code>
// objects; one for the shape
// functions on the current cell and
// the other for shape functions on
// the neighboring cell under
// consideration. Both objects are
- // either of class ``FEFaceValues''
- // or of class ``FESubfaceValues''
+ // either of class <code>FEFaceValues</code>
+ // or of class <code>FESubfaceValues</code>
// (both derived from
- // ``FEFaceValuesBase'') according to
+ // <code>FEFaceValuesBase</code>) according to
// the remaining cases mentioned
// in the introduction:
//
// 2. neighboring cell is finer
- // (current cell: ``FESubfaceValues'',
- // neighboring cell: ``FEFaceValues'');
+ // (current cell: <code>FESubfaceValues</code>,
+ // neighboring cell: <code>FEFaceValues</code>);
//
// 3. neighboring cell is of the same
// refinement level (both, current
// and neighboring cell:
- // ``FEFaceValues'');
+ // <code>FEFaceValues</code>);
//
// 4. neighboring cell is coarser
- // (current cell: ``FEFaceValues'',
+ // (current cell: <code>FEFaceValues</code>,
// neighboring cell:
- // ``FESubfaceValues'').
+ // <code>FESubfaceValues</code>).
//
// If we considered globally refined
// meshes then only case 3 would
void DGMethod<dim>::assemble_system1 ()
{
// First we create the
- // ``UpdateFlags'' for the
- // ``FEValues'' and the
- // ``FEFaceValues'' objects.
+ // <code>UpdateFlags</code> for the
+ // <code>FEValues</code> and the
+ // <code>FEFaceValues</code> objects.
const UpdateFlags update_flags = update_values
| update_gradients
| update_q_points
// vectors of the current cell.
const UpdateFlags neighbor_face_update_flags = update_values;
- // Then we create the ``FEValues''
+ // Then we create the <code>FEValues</code>
// object. Here, we use the default
// MappingQ1. different mapping
// create a MappingCollection first
hp::FEValues<dim> fe_v_x (fe_collection, quadratures, update_flags);
// Similarly we create the
- // ``FEFaceValues'' and
- // ``FESubfaceValues'' objects for
+ // <code>FEFaceValues</code> and
+ // <code>FESubfaceValues</code> objects for
// both, the current and the
// neighboring cell. Within the
// following nested loop over all
// and vectors. Here we need two
// cell matrices, both for face
// terms that include test
- // functions ``v'' (shape functions
+ // functions <code>v</code> (shape functions
// of the current cell). To be more
// precise, the first matrix will
// include the `u and v terms' and
std::vector<unsigned int> dofs_neighbor;
// In the
- // ``assemble_face_term1''
+ // <code>assemble_face_term1</code>
// function contributions to
// the cell matrices and the
// cell vector are only
// ADDED. Therefore on each
// cell we need to reset the
- // ``u_v_matrix'' and
- // ``cell_vector'' to zero,
+ // <code>u_v_matrix</code> and
+ // <code>cell_vector</code> to zero,
// before assembling the cell terms.
u_v_matrix = 0;
cell_vector = 0;
- // Now we reinit the ``FEValues''
+ // Now we reinit the <code>FEValues</code>
// object for the current cell
fe_v_x.reinit (cell);
// and call the function
// that assembles the cell
// terms. The first argument is
- // the ``FEValues'' that was
+ // the <code>FEValues</code> that was
// previously reinit'ed on the
// current cell.
dg.assemble_cell_term(fe_v_x.get_present_fe_values (),
typename hp::DoFHandler<dim>::face_iterator face=cell->face(face_no);
// and clear the
- // ``un_v_matrix'' on each
+ // <code>un_v_matrix</code> on each
// face.
un_v_matrix = 0;
if (face->at_boundary())
{
// We reinit the
- // ``FEFaceValues''
+ // <code>FEFaceValues</code>
// object to the
// current face
fe_v_face_x.reinit (cell, face_no);
// note that the
// following part of
// code will not work
- // for ``dim==1''.
+ // for <code>dim==1</code>.
if (face->has_children())
{
// First we store
// neighbor-@>neighbor(neighbor2)
// equals the
// current cell
- // ``cell''.
+ // <code>cell</code>.
const unsigned int neighbor2=
cell->neighbor_of_neighbor(face_no);
// and set the
// cell
// iterator
- // ``neighbor_child''
+ // <code>neighbor_child</code>
// to the cell
// placed
// `behind' the
// We need to
// reset the
- // ``un_v_matrix''
+ // <code>un_v_matrix</code>
// on each
// subface
// because on
// each subface
- // the ``un''
+ // the <code>un</code>
// belong to
// different
// neighboring
// case (case
// 2) we employ
// the
- // ``FESubfaceValues''
+ // <code>FESubfaceValues</code>
// of the
// current
// cell (here
// and
// distribute
- // ``un_v_matrix''
+ // <code>un_v_matrix</code>
// to the
// system_matrix
for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
system_matrix.add(dofs[i], dofs_neighbor[k],
un_v_matrix(i,k));
}
- // End of ``if
- // (face-@>has_children())''
+ // End of <code>if
+ // (face-@>has_children())</code>
}
else
{
// We reinit
// the
- // ``FEFaceValues''
+ // <code>FEFaceValues</code>
// of the
// current and
// neighboring
fe_v_face_neighbor_x.get_present_fe_values (),
u_v_matrix,
un_v_matrix);
- // End of ``if
+ // End of <code>if
// (neighbor-@>level()
// ==
- // cell-@>level())''
+ // cell-@>level())</code>
}
else
{
// Reinit the
// appropriate
- // ``FEFaceValues''
+ // <code>FEFaceValues</code>
// and assemble
// the face
// terms.
// Now we get the
// dof indices of
// the
- // ``neighbor_child''
+ // <code>neighbor_child</code>
// cell,
dofs_neighbor.resize (neighbor->get_fe().dofs_per_cell);
neighbor->get_dof_indices (dofs_neighbor);
// and distribute the
- // ``un_v_matrix''.
+ // <code>un_v_matrix</code>.
for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
for (unsigned int k=0; k<neighbor->get_fe().dofs_per_cell; ++k)
system_matrix.add(dofs[i], dofs_neighbor[k],
un_v_matrix(i,k));
}
- // End of ``face not at boundary'':
+ // End of <code>face not at boundary</code>:
}
// End of loop over all faces:
}
// Finally we distribute the
- // ``u_v_matrix''
+ // <code>u_v_matrix</code>
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add(dofs[i], dofs[j], u_v_matrix(i,j));
// @sect4{Function: assemble_system2}
//
// We proceed with the
- // ``assemble_system2'' function that
+ // <code>assemble_system2</code> function that
// implements the DG discretization
// in its second version. This
// function is very similar to the
- // ``assemble_system1''
+ // <code>assemble_system1</code>
// function. Therefore, here we only
// discuss the differences between
// the two functions. This function
// repeatedly calls the
- // ``assemble_face_term2'' function
+ // <code>assemble_face_term2</code> function
// of the DGTransportEquation object,
// that assembles the face terms
// written as a sum of integrals over
const UpdateFlags neighbor_face_update_flags = update_values;
// Here we do not need
- // ``fe_v_face_neighbor'' as case 4
+ // <code>fe_v_face_neighbor</code> as case 4
// does not occur.
hp::FEValues<dim> fe_v_x (
fe_collection, quadratures, update_flags);
// Additionally we need the
// following two cell matrices,
// both for face term that include
- // test function ``vn'' (shape
+ // test function <code>vn</code> (shape
// functions of the neighboring
// cell). To be more precise, the
// first matrix will include the `u
// difference quotients including the
// cell under consideration and its
// neighbors. This is done by the
- // ``DerivativeApproximation'' class
+ // <code>DerivativeApproximation</code> class
// that computes the approximate
// gradients in a way similar to the
- // ``GradientEstimation'' described
+ // <code>GradientEstimation</code> described
// in Step 9 of this tutorial. In
// fact, the
- // ``DerivativeApproximation'' class
+ // <code>DerivativeApproximation</code> class
// was developed following the
- // ``GradientEstimation'' class of
+ // <code>GradientEstimation</code> class of
// Step 9. Relating to the
// discussion in Step 9, here we
// consider $h^{1+d/2}|\nabla_h
template <int dim>
void DGMethod<dim>::refine_grid ()
{
- // The ``DerivativeApproximation''
+ // The <code>DerivativeApproximation</code>
// class computes the gradients to
// float precision. This is
// sufficient as they are
}
- // The following ``run'' function is
+ // The following <code>run</code> function is
// similar to previous examples. The
// only difference is that the
// problem is assembled and solved
// twice on each refinement step;
- // first by ``assemble_system1'' that
+ // first by <code>assemble_system1</code> that
// implements the first version and
- // then by ``assemble_system2'' that
+ // then by <code>assemble_system2</code> that
// implements the second version of
// writing the DG
// discretization. Furthermore the
}
}
- // The following ``main'' function is
+ // The following <code>main</code> function is
// similar to previous examples and
// need not to be commented on.
int main ()
#include <iostream>
- // @sect3{The ``LaplaceProblem'' class}
+ // @sect3{The <code>LaplaceProblem</code> class}
// Instead of the procedural programming of
// previous examples, we encapsulate
// freedom. This is done by using the
// distribute_dofs function, as we have
// seen in the step-2 example. Since we use
- // the ``FE_Q'' class with a polynomial
+ // the <code>FE_Q</code> class with a polynomial
// degree of 1, i.e. bilinear elements,
// this associates one degree of freedom
// with each vertex. While we're at
// instead of giving a magically obtained
// maximal number of nonzero entries per
// row, we now use a function in the
- // ``DoFHandler'' class that can compute
+ // <code>DoFHandler</code> class that can compute
// this number for us:
sparsity_pattern.reinit (dof_handler.n_dofs(),
dof_handler.n_dofs(),
// objects. That's too much, so there is one
// type of class that orchestrates
// information exchange between these three:
- // the ``FEValues'' class. If given one
+ // the <code>FEValues</code> class. If given one
// instance of each three of these objects,
// it will be able to provide you with
// information about values and gradients of
// actually need is given as a bitwise
// connection of flags as the third
// argument to the constructor of
- // ``FEValues''. Since these values have to
+ // <code>FEValues</code>. Since these values have to
// be recomputed, or updated, every time we
// go to a new cell, all of these flags
- // start with the prefix ``update_'' and
+ // start with the prefix <code>update_</code> and
// then indicate what it actually is that
// we want updated. The flag to give if we
// want the values of the shape functions
- // computed is ``update_values''; for the
+ // computed is <code>update_values</code>; for the
// gradients it is
- // ``update_gradients''. The determinants
+ // <code>update_gradients</code>. The determinants
// of the Jacobians and the quadrature
// weights are always used together, so
// only the products (Jacobians times
- // weights, or short ``JxW'') are computed;
+ // weights, or short <code>JxW</code>) are computed;
// since we need them, we have to list
- // ``update_JxW_values'' as well:
+ // <code>update_JxW_values</code> as well:
FEValues<2> fe_values (fe, quadrature_formula,
update_values | update_gradients | update_JxW_values);
// The advantage of this proceeding is that
// determinant and the quadrature point
// weight (that one gets together by
// the call to
- // ``fe_values.JxW''). Finally, this is
+ // <code>fe_values.JxW</code>). Finally, this is
// repeated for all shape functions
// phi_i and phi_j:
for (unsigned int i=0; i<dofs_per_cell; ++i)
// rather than projecting it onto the
// boundary. There is a function in the
// library which does exactly this:
- // ``VectorTools::interpolate_boundary_values''. Its
+ // <code>VectorTools::interpolate_boundary_values</code>. Its
// parameters are (omitting parameters for
// which default values exist and that we
// don't care about): the DoFHandler object
// the boundary.
//
// The function describing the boundary
- // values is an object of type ``Function''
+ // values is an object of type <code>Function</code>
// or of a derived class. One of the
- // derived classes is ``ZeroFunction'',
+ // derived classes is <code>ZeroFunction</code>,
// which describes (not unexpectedly) a
// function which is zero everywhere. We
// create such an object in-place and pass
// here for all entries). This
// mapping of DoF numbers to
// boundary values is done by the
- // ``std::map'' class.
+ // <code>std::map</code> class.
std::map<unsigned int,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
// First, we need to have an object that
// knows how to tell the CG algorithm when
// to stop. This is done by using a
- // ``SolverControl'' object, and as
+ // <code>SolverControl</code> object, and as
// stopping criterion we say: stop after a
// maximum of 1000 iterations (which is far
// more than is needed for 1089 variables;
// the one which stops the iteration:
SolverControl solver_control (1000, 1e-12);
// Then we need the solver itself. The
- // template parameters to the ``SolverCG''
+ // template parameters to the <code>SolverCG</code>
// class are the matrix type and the type
// of the vectors, but the empty angle
// brackets indicate that we simply take
// the default arguments (which are
- // ``SparseMatrix@<double@>'' and
- // ``Vector@<double@>''):
+ // <code>SparseMatrix@<double@></code> and
+ // <code>Vector@<double@></code>):
SolverCG<> cg (solver_control);
// Now solve the system of equations. The
// To write the output to a file,
// we need an object which knows
// about output formats and the
- // like. This is the ``DataOut'' class,
+ // like. This is the <code>DataOut</code> class,
// and we need an object of that
// type:
DataOut<2> data_out;
// Now we have to tell it where to take the
// values from which it shall write. We
- // tell it which ``DoFHandler'' object to
+ // tell it which <code>DoFHandler</code> object to
// use, and the solution vector (and
// the name by which the solution variable
// shall appear in the output file). If
// handle. The reason is that we
// have separated the frontend
// (which knows about how to treat
- // ``DoFHandler'' objects and data
+ // <code>DoFHandler</code> objects and data
// vectors) from the back end (which
// knows many different output formats)
// and use an intermediate data
// Finally, the last function of this class
// is the main function which calls all the
- // other functions of the ``LaplaceProblem''
+ // other functions of the <code>LaplaceProblem</code>
// class. The order in which this is done
// resembles the order in which most finite
// element programs work. Since the names are
}
- // @sect3{The ``main'' function}
+ // @sect3{The <code>main</code> function}
// This is the main function of the
// program. Since the concept of a
#include <base/logstream.h>
- // @sect3{The ``LaplaceProblem'' class template}
+ // @sect3{The <code>LaplaceProblem</code> class template}
// This is again the same
- // ``LaplaceProblem'' class as in the
+ // <code>LaplaceProblem</code> class as in the
// previous example. The only
// difference is that we have now
// declared it as a class with a
// knows the size of the loop at compile time
// (remember that at the time when you define
// the template, the compiler doesn't know
- // the value of ``dim'', but when it later
+ // the value of <code>dim</code>, but when it later
// encounters a statement or declaration
- // ``RightHandSide@<2@>'', it will take the
+ // <code>RightHandSide@<2@></code>, it will take the
// template, replace all occurrences of dim
// by 2 and compile the resulting function);
// in other words, at the time of compiling
// above right away.
//
// The last thing to note is that a
- // ``Point@<dim@>'' denotes a point in
+ // <code>Point@<dim@></code> denotes a point in
// dim-dimensionsal space, and its individual
// components (i.e. `x', `y',
// ... coordinates) can be accessed using the
- // @sect3{Implementation of the ``LaplaceProblem'' class}
+ // @sect3{Implementation of the <code>LaplaceProblem</code> class}
// Next for the implementation of the class
// template that makes use of the functions
// above. As before, we will write everything
// as templates that have a formal parameter
- // ``dim'' that we assume unknown at the time
+ // <code>dim</code> that we assume unknown at the time
// we define the template functions. Only
// later, the compiler will find a
- // declaration of ``LaplaceProblem@<2@>'' (in
- // the ``main'' function, actually) and
- // compile the entire class with ``dim''
+ // declaration of <code>LaplaceProblem@<2@></code> (in
+ // the <code>main</code> function, actually) and
+ // compile the entire class with <code>dim</code>
// replaced by 2, a process referred to as
// `instantiation of a template'. When doing
// so, it will also replace instances of
- // ``RightHandSide@<dim@>'' by
- // ``RightHandSide@<2@>'' and instantiate the
+ // <code>RightHandSide@<dim@></code> by
+ // <code>RightHandSide@<2@></code> and instantiate the
// latter class from the class template.
//
// In fact, the compiler will also find a
- // declaration ``LaplaceProblem@<3@>'' in
- // ``main()''. This will cause it to again go
+ // declaration <code>LaplaceProblem@<3@></code> in
+ // <code>main()</code>. This will cause it to again go
// back to the general
- // ``LaplaceProblem@<dim@>'' template, replace
- // all occurrences of ``dim'', this time by
+ // <code>LaplaceProblem@<dim@></code> template, replace
+ // all occurrences of <code>dim</code>, this time by
// 3, and compile the class a second
// time. Note that the two instantiations
- // ``LaplaceProblem@<2@>'' and
- // ``LaplaceProblem@<3@>'' are completely
+ // <code>LaplaceProblem@<2@></code> and
+ // <code>LaplaceProblem@<3@></code> are completely
// independent classes; their only common
// feature is that they are both instantiated
// from the same general template, but they
// @sect4{LaplaceProblem::LaplaceProblem}
// After this introduction, here is the
- // constructor of the ``LaplaceProblem''
+ // constructor of the <code>LaplaceProblem</code>
// class. It specifies the desired polynomial
// degree of the finite elements and
// associates the DoFHandler to the
// square [-1,1]x[-1,1] in 2D, or on
// the cube [-1,1]x[-1,1]x[-1,1] in
// 3D; both can be termed
- // ``hyper_cube'', so we may use the
+ // <code>hyper_cube</code>, so we may use the
// same function in whatever
// dimension we are. Of course, the
// functions that create a hypercube
// either. This function therefore looks
// exactly like in the previous example,
// although it performs actions that in their
- // details are quite different if ``dim''
+ // details are quite different if <code>dim</code>
// happens to be 3. The only significant
// difference from a user's perspective is
// the number of cells resulting, which is
// quadrature points on the cell we are
// presently on (previously, we only
// required values and gradients of the
- // shape function from the ``FEValues''
+ // shape function from the <code>FEValues</code>
// object, as well as the quadrature
- // weights, ``JxW''). We can tell the
- // ``FEValues'' object to do for us by also
- // giving it the ``update_q_points'' flag:
+ // weights, <code>JxW</code>). We can tell the
+ // <code>FEValues</code> object to do for us by also
+ // giving it the <code>update_q_points</code> flag:
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_q_points | update_JxW_values);
// Note, that a cell is a quadrilateral in
// two space dimensions, but a hexahedron
// in 3D. In fact, the
- // ``active_cell_iterator'' data type is
+ // <code>active_cell_iterator</code> data type is
// something different, depending on the
// dimension we are in, but to the outside
// world they look alike and you will
// and j at point q_point and multiply
// it with the scalar weights JxW. This
// is actually what happens:
- // ``fe_values.shape_grad(i,q_point)''
- // returns a ``dim'' dimensional
+ // <code>fe_values.shape_grad(i,q_point)</code>
+ // returns a <code>dim</code> dimensional
// vector, represented by a
- // ``Tensor@<1,dim@>'' object, and the
+ // <code>Tensor@<1,dim@></code> object, and the
// operator* that multiplies it with
// the result of
- // ``fe_values.shape_grad(j,q_point)''
- // makes sure that the ``dim''
+ // <code>fe_values.shape_grad(j,q_point)</code>
+ // makes sure that the <code>dim</code>
// components of the two vectors are
// properly contracted, and the result
// is a scalar floating point number
// that then is multiplied with the
// weights. Internally, this operator*
// makes sure that this happens
- // correctly for all ``dim'' components
- // of the vectors, whether ``dim'' be
+ // correctly for all <code>dim</code> components
+ // of the vectors, whether <code>dim</code> be
// 2, 3, or any other space dimension;
// from a user's perspective, this is
// not something worth bothering with,
// values in this example, contrary to the
// one before. This is a simple task, we
// only have to replace the
- // ``ZeroFunction'' used there by an object
+ // <code>ZeroFunction</code> used there by an object
// of the class which describes the
// boundary values we would like to use
- // (i.e. the ``BoundaryValues'' class
+ // (i.e. the <code>BoundaryValues</code> class
// declared above):
std::map<unsigned int,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
// obtain it see the ReadMe file of
// deal.II). To write data in this format, we
// simply replace the
- // ``data_out.write_gnuplot'' call by
- // ``data_out.write_gmv''.
+ // <code>data_out.write_gnuplot</code> call by
+ // <code>data_out.write_gmv</code>.
//
// Since the program will run both 2d and 3d
// versions of the laplace solver, we use the
}
- // @sect3{The ``main'' function}
+ // @sect3{The <code>main</code> function}
// And this is the main function. It also
// looks mostly like in step-3, but if you
// look at the code below, note how we first
// create a variable of type
- // ``LaplaceProblem@<2@>'' (forcing the
+ // <code>LaplaceProblem@<2@></code> (forcing the
// compiler to compile the class template
- // with ``dim'' replaced by ``2'') and run a
+ // with <code>dim</code> replaced by <code>2</code>) and run a
// 2d simulation, and then we do the whole
// thing over in 3d.
//
//
// Each of the two blocks is enclosed in
// braces to make sure that the
- // ``laplace_problem_2d'' variable goes out
+ // <code>laplace_problem_2d</code> variable goes out
// of scope (and releases the memory it
// holds) before we move on to allocate
// memory for the 3d case. Without the
// additional braces, the
- // ``laplace_problem_2d'' variable would only
+ // <code>laplace_problem_2d</code> variable would only
// be destroyed at the end of the function,
// i.e. after running the 3d problem, and
// would needlessly hog memory while the 3d
// starting residual and the number of the
// iteration where convergence was
// detected. This can be suppressed through
- // the ``deallog.depth_console(0)'' call.
+ // the <code>deallog.depth_console(0)</code> call.
//
// The rationale here is the following: the
// deallog (i.e. deal-log, not de-allog)
// ... and this is too: We will
// convert integers to strings using
// the C++ stringstream class
- // ``ostringstream'':
+ // <code>ostringstream</code>:
#include <sstream>
- // @sect3{The ``LaplaceProblem'' class template}
+ // @sect3{The <code>LaplaceProblem</code> class template}
// The main class is mostly as in the
// previous example. The most visible
// change is that the function
- // ``make_grid_and_dofs'' has been
+ // <code>make_grid_and_dofs</code> has been
// removed, since creating the grid
- // is now done in the ``run''
+ // is now done in the <code>run</code>
// function and the rest of its
// functionality is now in
- // ``setup_system''. Apart from this,
+ // <code>setup_system</code>. Apart from this,
// everything is as before.
template <int dim>
class LaplaceProblem
};
- // @sect3{Nonconstant coefficients, using ``Assert''}
+ // @sect3{Nonconstant coefficients, using <code>Assert</code>}
// In step-4, we showed how to use
// non-constant boundary values and
// we want to use a variable
// coefficient in the elliptic
// operator instead. Of course, the
- // suitable object is a ``Function'',
+ // suitable object is a <code>Function</code>,
// as we have used for the right hand
// side and boundary values in the
// last example. We will use it
// again, but we implement another
- // function ``value_list'' which
+ // function <code>value_list</code> which
// takes a list of points and returns
// the values of the function at
// these points as a list. The reason
// why such a function is reasonable
// although we can get all the
- // information from the ``value''
+ // information from the <code>value</code>
// function as well will be explained
// below when assembling the matrix.
//
// points at once. Of course, we need
// to make sure that the values are
// the same as if we would ask the
- // ``value'' function for each point
+ // <code>value</code> function for each point
// individually.
//
// This method takes three
// component that should be zero here
// since we only have a single scalar
// function. Now, of course the size
- // of the output array (``values'')
+ // of the output array (<code>values</code>)
// must be the same as that of the
- // input array (``points''), and we
+ // input array (<code>points</code>), and we
// could simply assume that. However,
// in practice, it turns out that
// more than 90 per cent of
// invalid array sizes, etc, so we
// should try to make sure that the
// parameters are valid. For this,
- // the ``Assert'' macro is a good means,
+ // the <code>Assert</code> macro is a good means,
// since it makes sure that the
// condition which is given as first
// argument is valid, and if not
// should not slow down the program
// too much if you want to do large
// computations. To this end, the
- // ``Assert'' macro is only used in
+ // <code>Assert</code> macro is only used in
// debug mode and expands to nothing
// if in optimized mode. Therefore,
// while you test your program on
// mode to optimized mode is to go
// edit the Makefile in this
// directory. It should have a line
- // ``debug-mode = on''; simply
- // replace it by ``debug-mode = off''
+ // <code>debug-mode = on</code>; simply
+ // replace it by <code>debug-mode = off</code>
// and recompile your program. The
- // output of the ``make'' program
+ // output of the <code>make</code> program
// should already indicate to you
// that the program is now compiled
// in optimized mode, and it will
// two arrays is one of the most
// frequent checks, which is why
// there is already an exception
- // class ``ExcDimensionMismatch''
+ // class <code>ExcDimensionMismatch</code>
// that takes the sizes of two
// vectors and prints some output in
// case the condition is violated:
// trigger this exception at the
// end of the main program, and
// what output results from this
- // (see the ``Results'' section of
+ // (see the <code>Results</code> section of
// this example program). You will
// certainly notice that the output
// is quite well suited to quickly
// While we're at it, we can do
// another check: the coefficient
// is a scalar, but the
- // ``Function'' class also
+ // <code>Function</code> class also
// represents vector-valued
// function. A scalar function must
// therefore be considered as a
// ask is zero (we always count
// from zero). The following
// assertion checks this. If the
- // condition in the ``Assert'' call
+ // condition in the <code>Assert</code> call
// is violated, an exception of
- // type ``ExcRange'' will be
+ // type <code>ExcRange</code> will be
// triggered; that class takes the
// violating index as first
// argument, and the second and
// zero, of course. (The interval
// is half open since we also want
// to write exceptions like
- // ``ExcRange(i,0,v.size())'',
+ // <code>ExcRange(i,0,v.size())</code>,
// where an index must be between
// zero but less than the size of
// an array. To save us the effort
- // of writing ``v.size()-1'' in
+ // of writing <code>v.size()-1</code> in
// many places, the range is
// defined as half-open.)
Assert (component == 0,
// The rest of the function is
// uneventful: we define
- // ``n_q_points'' as an
+ // <code>n_q_points</code> as an
// abbreviation for the number of
// points for which function values
// are requested, and then simply
}
- // @sect3{The ``LaplaceProblem'' class implementation}
+ // @sect3{The <code>LaplaceProblem</code> class implementation}
// @sect4{LaplaceProblem::LaplaceProblem}
// @sect4{LaplaceProblem::setup_system}
// This is the function
- // ``make_grid_and_dofs'' from the
+ // <code>make_grid_and_dofs</code> from the
// previous example, minus the
// generation of the grid. Everything
// else is unchanged:
// here.
//
// Then, below, we will ask the
- // ``coefficient'' function object
+ // <code>coefficient</code> function object
// to compute the values of the
// coefficient at all quadrature
// points on one cell at once. The
// use it in the computation of the
// local contributions. This is
// what we do in the call to
- // ``coefficient.value_list'' in
+ // <code>coefficient.value_list</code> in
// the fourth line of the loop.
//
// The second change is how we make
// (symmetric successive
// overrelaxation), with a relaxation
// factor of 1.2. For this purpose,
- // the ``SparseMatrix'' class has a
+ // the <code>SparseMatrix</code> class has a
// function which does one SSOR step,
// and we need to package the address
// of this function together with the
// (which is the matrix to be
// inverted) and the relaxation
// factor into one object. The
- // ``PreconditionSSOR'' class does
- // this for us. (``PreconditionSSOR''
+ // <code>PreconditionSSOR</code> class does
+ // this for us. (<code>PreconditionSSOR</code>
// class takes a template argument
// denoting the matrix type it is
// supposed to work on. The default
- // value is ``SparseMatrix@<double@>'',
+ // value is <code>SparseMatrix@<double@></code>,
// which is exactly what we need
// here, so we simply stick with the
// default and do not specify
//
// With this, the rest of the
// function is trivial: instead of
- // the ``PreconditionIdentity''
+ // the <code>PreconditionIdentity</code>
// object we have created before, we
// now use the preconditioner we have
// declared, and the CG solver will
// which the results are to be
// written. We would like to have
// it of the form
- // ``solution-N.eps'', where N is
+ // <code>solution-N.eps</code>, where N is
// the number of the refinement
// cycle. Thus, we have to convert
// an integer to a part of a
// string; this can be done using
- // the ``sprintf'' function, but in
+ // the <code>sprintf</code> function, but in
// C++ there is a more elegant way:
// write everything into a special
// stream (just like writing into a
// conversions from integer to
// strings, and one could as well
// use stream modifiers such as
- // ``setw'', ``setprecision'', and
+ // <code>setw</code>, <code>setprecision</code>, and
// so on. In C++, you can do this
// by using the so-called stringstream
// classes:
<< ".eps";
// We can get whatever we wrote to the
- // stream using the ``str()'' function. The
+ // stream using the <code>str()</code> function. The
// result is a string which we have to
- // convert to a char* using the ``c_str()''
+ // convert to a char* using the <code>c_str()</code>
// function. Use that as filename for the
// output stream and then write the data to
// the file:
// The second to last thing in this
// program is the definition of the
- // ``run()'' function. In contrast to
+ // <code>run()</code> function. In contrast to
// the previous programs, we will
// compute on a sequence of meshes
// that after each iteration is
// previous examples, we
// have already used some
// of the functions from
- // the ``GridGenerator''
+ // the <code>GridGenerator</code>
// class. Here we would
// like to read a grid from
// a file where the cells
// particular in this case
// one would of course try
// to do something else if
- // ``dim'' is not equal to
+ // <code>dim</code> is not equal to
// two, e.g. create a grid
// using library
// functions. Aborting a
// grid. It is in UCD
// (unstructured cell data)
// format (but the ending
- // of the ``UCD''-file is
- // ``inp''), as supported
+ // of the <code>UCD</code>-file is
+ // <code>inp</code>), as supported
// as input format by the
// AVS Explorer (a
// visualization program),
// If you like to use
// another input format,
// you have to use an other
- // ``grid_in.read_xxx''
+ // <code>grid_in.read_xxx</code>
// function. (See the
// documentation of the
- // ``GridIn'' class to find
+ // <code>GridIn</code> class to find
// out what input formats
// are presently
// supported.)
}
- // @sect3{The ``main'' function}
+ // @sect3{The <code>main</code> function}
// The main function looks mostly
// like the one in the previous
// Finally, we have promised to
// trigger an exception in the
- // ``Coefficient'' class through
- // the ``Assert'' macro we have
+ // <code>Coefficient</code> class through
+ // the <code>Assert</code> macro we have
// introduced there. For this, we
- // have to call its ``value_list''
+ // have to call its <code>value_list</code>
// function with two arrays of
// different size (the number in
// parentheses behind the
// will import the declaration of
// H1-conforming finite element shape
// functions. This family of finite
- // elements is called ``FE_Q'', and
+ // elements is called <code>FE_Q</code>, and
// was used in all examples before
// already to define the usual bi- or
// tri-linear elements, but we will
// refined grids (just the grid, not
// the solution) in each step, so we
// need the following include file
- // instead of ``grid_in.h'':
+ // instead of <code>grid_in.h</code>:
#include <grid/grid_out.h>
// When using locally refined grids,
- // we will get so-called ``hanging
- // nodes''. However, the standard
+ // we will get so-called <code>hanging
+ // nodes</code>. However, the standard
// finite element methods assumes
// that the discrete solution spaces
// be continuous, so we need to make
- // @sect3{The ``LaplaceProblem'' class template}
+ // @sect3{The <code>LaplaceProblem</code> class template}
// The main class is again almost
// unchanged. Two additions, however,
// are made: we have added the
- // ``refine_grid'' function, which is
+ // <code>refine_grid</code> function, which is
// used to adaptively refine the grid
// (instead of the global refinement
// in the previous examples), and a
}
- // @sect3{The ``LaplaceProblem'' class implementation}
+ // @sect3{The <code>LaplaceProblem</code> class implementation}
// @sect4{LaplaceProblem::LaplaceProblem}
// quadratic element. To do so, we
// only have to replace the
// constructor argument (which was
- // ``1'' in all previous examples) by
+ // <code>1</code> in all previous examples) by
// the desired polynomial degree
- // (here ``2''):
+ // (here <code>2</code>):
template <int dim>
LaplaceProblem<dim>::LaplaceProblem () :
dof_handler (triangulation),
// to add it is a subtle change in
// the order of data elements in the
// class as compared to all previous
- // examples: the ``dof_handler''
+ // examples: the <code>dof_handler</code>
// object was defined before and not
- // after the ``fe'' object. Of course
+ // after the <code>fe</code> object. Of course
// we could have left this order
// unchanged, but we would like to
// show what happens if the order is
// following: when we distribute the
// degrees of freedom using the
// function call
- // ``dof_handler.distribute_dofs()'',
- // the ``dof_handler'' also stores a
+ // <code>dof_handler.distribute_dofs()</code>,
+ // the <code>dof_handler</code> also stores a
// pointer to the finite element in
// use. Since this pointer is used
// every now and then until either
// the degrees of freedom are
// re-distributed using another
// finite element object or until the
- // ``dof_handler'' object is
+ // <code>dof_handler</code> object is
// destroyed, it would be unwise if
// we would allow the finite element
// object to be deleted before the
- // ``dof_handler'' object. To
+ // <code>dof_handler</code> object. To
// disallow this, the DoF handler
// increases a counter inside the
// finite element object which counts
// how many objects use that finite
// element (this is what the
- // ``Subscriptor''/``SmartPointer''
+ // <code>Subscriptor</code>/<code>SmartPointer</code>
// class pair is used for, in case
// you want something like this for
// your own programs; see step-7 for
// exactly the behavior sketched
// above. The reason is that member
// variables of the
- // ``LaplaceProblem'' class are
+ // <code>LaplaceProblem</code> class are
// destructed bottom-up (i.e. in
// reverse order of their declaration
// in the class), as always in
// declaration is below the one of
// the DoF handler. This triggers the
// situation above, and an exception
- // will be raised when the ``fe''
+ // will be raised when the <code>fe</code>
// object is destructed. What needs
// to be done is to tell the
- // ``dof_handler'' object to release
+ // <code>dof_handler</code> object to release
// its lock to the finite element. Of
- // course, the ``dof_handler'' will
+ // course, the <code>dof_handler</code> will
// only release its lock if it really
// does not need the finite element
// any more, i.e. when all finite
// element related data is deleted
// from it. For this purpose, the
- // ``DoFHandler'' class has a
- // function ``clear'' which deletes
+ // <code>DoFHandler</code> class has a
+ // function <code>clear</code> which deletes
// all degrees of freedom, and
// releases its lock to the finite
// element. After this, you can
// there are now 9 degrees of freedom
// per cell, not only four, that can
// couple with each other. The
- // ``dof_Handler.max_couplings_between_dofs()''
+ // <code>dof_Handler.max_couplings_between_dofs()</code>
// call will take care of this,
// however:
template <int dim>
// hanging nodes. In the class
// desclaration, we have already
// allocated space for an object
- // ``hanging_node_constraints''
+ // <code>hanging_node_constraints</code>
// that will hold a list of these
// constraints (they form a matrix,
// which is reflected in the name
DoFTools::make_hanging_node_constraints (dof_handler,
hanging_node_constraints);
- // The next step is ``closing''
+ // The next step is <code>closing</code>
// this object. For this note that,
// in principle, the
- // ``ConstraintMatrix'' class can
+ // <code>ConstraintMatrix</code> class can
// hold other constraints as well,
// i.e. constraints that do not
// stem from hanging
// nodes. Sometimes, it is useful
// to use such constraints, in
// which case they may be added to
- // the ``ConstraintMatrix'' object
+ // the <code>ConstraintMatrix</code> object
// after the hanging node
// constraints were computed. After
// all constraints have been added,
// rearranged to perform some
// actions more efficiently. This
// postprocessing is done using the
- // ``close()'' function, after which
+ // <code>close()</code> function, after which
// no further constraints may be
// added any more:
hanging_node_constraints.close ();
// some space for them here. Since
// the process of elimination of
// these constrained nodes is
- // called ``condensation'', the
+ // called <code>condensation</code>, the
// functions that eliminate them
- // are called ``condense'' for both
+ // are called <code>condense</code> for both
// the system matrix and right hand
// side, as well as for the
// sparsity pattern.
// polynomial degree in the finite
// element shape functions. This is
// easy to change: the constructor of
- // the ``QGauss'' class takes the
+ // the <code>QGauss</code> class takes the
// number of quadrature points in
// each space direction. Previously,
// we had two points for bilinear
// worth noting, however, that under
// the hood several things are
// different than before. First, the
- // variables ``dofs_per_cell'' and
- // ``n_q_points'' now are 9 each,
+ // variables <code>dofs_per_cell</code> and
+ // <code>n_q_points</code> now are 9 each,
// where they were 4
// before. Introducing such variables
// as abbreviations is a good
// strategy to make code work with
// different elements without having
// to change too much code. Secondly,
- // the ``fe_values'' object of course
+ // the <code>fe_values</code> object of course
// needs to do other things as well,
// since the shape functions are now
// quadratic, rather than linear, in
// of the matrix and all other
// entries for this line are set to
// zero) but the computed values
- // are invalid (the ``condense''
+ // are invalid (the <code>condense</code>
// function modifies the system so
// that the values in the solution
// corresponding to constrained
// system still has a well-defined
// solution; we compute the correct
// values for these nodes at the
- // end of the ``solve'' function).
+ // end of the <code>solve</code> function).
// As almost all the stuff before,
// the interpolation of boundary
// have to do is to use the
// constraints to assign to them the
// values that they should have. This
- // process, called ``distributing''
+ // process, called <code>distributing</code>
// hanging nodes, computes the values
// of constrained nodes from the
// values of the unconstrained ones,
// Instead of global refinement, we
// now use a slightly more elaborate
// scheme. We will use the
- // ``KellyErrorEstimator'' class
+ // <code>KellyErrorEstimator</code> class
// which implements an error
// estimator for the Laplace
// equation; it can in principle
// way to test an adaptive program.
//
// The way the estimator works is to
- // take a ``DoFHandler'' object
+ // take a <code>DoFHandler</code> object
// describing the degrees of freedom
// and a vector of values for each
// degree of freedom as input and
// for each active cell of the
// triangulation (i.e. one value for
// each of the
- // ``triangulation.n_active_cells()''
+ // <code>triangulation.n_active_cells()</code>
// cells). To do so, it needs two
// additional pieces of information:
// a quadrature formula on the faces
// (i.e. quadrature formula on
- // ``dim-1'' dimensional objects. We
+ // <code>dim-1</code> dimensional objects. We
// use a 3-point Gauss rule again, a
// pick that is consistent and
// appropriate with the choice
// corresponding Neumann values. This
// information is represented by an
// object of type
- // ``FunctionMap@<dim@>::type'' that is
+ // <code>FunctionMap@<dim@>::type</code> that is
// essentially a map from boundary
// indicators to function objects
// describing Neumann boundary values
// The above function returned one
// error indicator value for each
// cell in the
- // ``estimated_error_per_cell''
+ // <code>estimated_error_per_cell</code>
// array. Refinement is now done as
// follows: refine those 30 per
// cent of the cells with the
// hack with an explicit assertion at
// the beginning of the function. If
// this assertion is triggered,
- // i.e. when ``cycle'' is larger than
+ // i.e. when <code>cycle</code> is larger than
// or equal to 10, an exception of
- // type ``ExcNotImplemented'' is
+ // type <code>ExcNotImplemented</code> is
// raised, indicating that some
// functionality is not implemented
// for this case (the functionality
// @sect4{LaplaceProblem::run}
// The final function before
- // ``main()'' is again the main
- // driver of the class, ``run()''. It
+ // <code>main()</code> is again the main
+ // driver of the class, <code>run()</code>. It
// is similar to the one of step-5,
// except that we generate a file in
// the program again instead of
}
- // @sect3{The ``main'' function}
+ // @sect3{The <code>main</code> function}
// The main function is unaltered in
// its functionality from the
// for all, this kind of exceptions
// is not switched off in optimized
// mode, in contrast to the
- // ``Assert'' macro which we have
+ // <code>Assert</code> macro which we have
// used to test against programming
// errors. If uncaught, these
// exceptions propagate the call tree
- // up to the ``main'' function, and
+ // up to the <code>main</code> function, and
// if they are not caught there
// either, the program is aborted. In
// many cases, like if there is not
// useful to write any larger program
// in this way, and you can do so by
// more or less copying this function
- // except for the ``try'' block that
+ // except for the <code>try</code> block that
// actually encodes the functionality
// particular to the present
// application.
// exception that was thrown is an
// object of a class that is
// derived from the C++ standard
- // class ``exception'', then we can
- // use the ``what'' member function
+ // class <code>exception</code>, then we can
+ // use the <code>what</code> member function
// to get a string which describes
// the reason why the exception was
// thrown.
// The deal.II exception classes
// are all derived from the
// standard class, and in
- // particular, the ``exc.what()''
+ // particular, the <code>exc.what()</code>
// function will return
// approximately the same string as
// would be generated if the
// exception was thrown using the
- // ``Assert'' macro. You have seen
+ // <code>Assert</code> macro. You have seen
// the output of such an exception
// in the previous example, and you
// then know that it contains the
// much that we can do except
// exiting the program with an
// error code (this is what the
- // ``return 1;'' does):
+ // <code>return 1;</code> does):
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
// If the exception that was thrown
// somewhere was not an object of a
// class derived from the standard
- // ``exception'' class, then we
+ // <code>exception</code> class, then we
// can't do anything at all. We
// then simply print an error
// message and exit.
// In this example, we will not use the
// numeration scheme which is used per
- // default by the ``DoFHandler'' class, but
+ // default by the <code>DoFHandler</code> class, but
// will renumber them using the Cuthill-McKee
// algorithm. As has already been explained
// in step-2, the necessary functions are
// how we can make sure that objects
// are not deleted while they are
// still in use. For this purpose,
- // deal.II has the ``SmartPointer''
+ // deal.II has the <code>SmartPointer</code>
// helper class, which is declared in
// this file:
#include <base/smartpointer.h>
// Next, we will want to use the
- // ``integrate_difference'' function
+ // <code>integrate_difference</code> function
// mentioned in the introduction, and we are
- // going to use a ``ConvergenceTable'' that
+ // going to use a <code>ConvergenceTable</code> that
// collects all important data during a run
// and prints it at the end as a table. These
// comes from the following two files:
#include <numerics/vectors.h>
#include <base/convergence_table.h>
// And finally, we need to use the
- // ``FEFaceValues'' class, which is
+ // <code>FEFaceValues</code> class, which is
// declared in the same file as the
- // ``FEValues'' class:
+ // <code>FEValues</code> class:
#include <fe/fe_values.h>
// We need one more include from standard
// First we assign values to the centers for
// the 1d case, where we place the centers
// equidistantly at -1/3, 0, and 1/3. The
- // ``template @<@>'' header for this definition
+ // <code>template @<@></code> header for this definition
// indicates an explicit specialization. This
// means, that the variable belongs to a
// template, but that instead of providing
// the compiler with a template from which it
// can specialize a concrete variable by
- // substituting ``dim'' with some concrete
+ // substituting <code>dim</code> with some concrete
// value, we provide a specialization
- // ourselves, in this case for ``dim=1''. If
+ // ourselves, in this case for <code>dim=1</code>. If
// the compiler then sees a reference to this
// variable in a place where the template
// argument equals one, it knows that it
// doesn't have to generate the variable from
- // a template by substituting ``dim'', but
+ // a template by substituting <code>dim</code>, but
// can immediately use the following
// definition:
template <>
Point<1>(+1.0 / 3.0) };
// Likewise, we can provide an explicit
- // specialization for ``dim=2''. We place the
+ // specialization for <code>dim=2</code>. We place the
// centers for the 2d case as follows:
template <>
const Point<2>
// dimensions. In this case, we simply
// provide the compiler with a template from
// which it can generate a concrete
- // instantiation by substituting ``dim'' with
+ // instantiation by substituting <code>dim</code> with
// a concrete value:
template <int dim>
const double SolutionBase<dim>::width = 1./3.;
// the classes representing these
// two. They both represent
// continuous functions, so they are
- // derived from the ``Function@<dim@>''
+ // derived from the <code>Function@<dim@></code>
// base class, and they also inherit
// the characteristics defined in the
- // ``SolutionBase'' class.
+ // <code>SolutionBase</code> class.
//
// The actual classes are declared in the
// following. Note that in order to compute
// more than we have done in previous
// examples, where all we provided was the
// value at one or a list of
- // points. Fortunately, the ``Function''
+ // points. Fortunately, the <code>Function</code>
// class also has virtual functions for the
// gradient, so we can simply overload the
// respective virtual member functions in the
- // ``Function'' base class. Note that the
- // gradient of a function in ``dim'' space
- // dimensions is a vector of size ``dim'',
+ // <code>Function</code> base class. Note that the
+ // gradient of a function in <code>dim</code> space
+ // dimensions is a vector of size <code>dim</code>,
// i.e. a tensor of rank 1 and dimension
- // ``dim''. As for so many other things, the
+ // <code>dim</code>. As for so many other things, the
// library provides a suitable class for
// this.
//
// elements of a base class that is
// template dependent (in this case
// the elements of
- // ``SolutionBase@<dim@>''), then the
+ // <code>SolutionBase@<dim@></code>), then the
// C++ language forces us to write
- // ``this-@>n_source_centers'' (for
- // example). Note that the ``this-@>''
+ // <code>this-@>n_source_centers</code> (for
+ // example). Note that the <code>this-@></code>
// qualification is not necessary if
// the base class is not template
// dependent, and also that the gcc
// is necessary is complicated; some
// books on C++ may explain it, so if
// you are interested you can look it
- // up under the phrase ``two-stage
- // (name) lookup''.
+ // up under the phrase <code>two-stage
+ // (name) lookup</code>.
template <int dim>
double Solution<dim>::value (const Point<dim> &p,
const unsigned int) const
// gradient of the solution. In order to
// accumulate the gradient from the
// contributions of the exponentials, we
- // allocate an object ``return_value'' that
+ // allocate an object <code>return_value</code> that
// denotes the mathematical quantity of a
- // tensor of rank ``1'' and dimension
- // ``dim''. Its default constructor sets it
+ // tensor of rank <code>1</code> and dimension
+ // <code>dim</code>. Its default constructor sets it
// to the vector containing only zeroes, so
// we need not explicitly care for its
// initialization.
//
// Note that we could as well have taken the
- // type of the object to be ``Point@<dim@>''
- // instead of ``Tensor@<1,dim@>''. Tensors of
+ // type of the object to be <code>Point@<dim@></code>
+ // instead of <code>Tensor@<1,dim@></code>. Tensors of
// rank 1 and points are almost exchangeable,
// and have only very slightly different
// mathematical meanings. In fact, the
- // ``Point@<dim@>'' class is derived from the
- // ``Tensor@<1,dim@>'' class, which makes up
+ // <code>Point@<dim@></code> class is derived from the
+ // <code>Tensor@<1,dim@></code> class, which makes up
// for their mutual exchange ability. Their
// main difference is in what they logically
// mean: points are points in space, such as
// mode as arguments.
//
// The rest of the member functions are as
- // before except for the ``process_solution''
+ // before except for the <code>process_solution</code>
// function: After the solution has been
// computed, we perform some analysis on it,
// such as computing the error in various
// a triangulation object, and we
// have a finite element object,
// and we also have an object of
- // type ``DoFHandler'' that uses
+ // type <code>DoFHandler</code> that uses
// both of the first two. These
// three objects all have a
// lifetime that is rather long
// they are destroyed at the very
// end. The question is: can we
// guarantee that the two objects
- // which the ``DoFHandler'' uses,
+ // which the <code>DoFHandler</code> uses,
// live at least as long as they
// are in use? This means that
- // the ``DoFHandler'' must have some
+ // the <code>DoFHandler</code> must have some
// kind of lock on the
// destruction of the other
// objects, and it can only
// to such potentially dangerous
// pointers are derived from a
// class called
- // ``Subscriptor''. For example,
- // the ``Triangulation'',
- // ``DoFHandler'', and a base
- // class of the ``FiniteElement''
+ // <code>Subscriptor</code>. For example,
+ // the <code>Triangulation</code>,
+ // <code>DoFHandler</code>, and a base
+ // class of the <code>FiniteElement</code>
// class are derived from
- // ``Subscriptor''. This latter
+ // <code>Subscriptor</code>. This latter
// class does not offer much
// functionality, but it has a
// built-in counter which we can
//
// On the other hand, if an object of a
// class that is derived from the
- // ``Subscriptor'' class is destroyed, it
+ // <code>Subscriptor</code> class is destroyed, it
// also has to call the destructor of the
- // ``Subscriptor'' class. In this
+ // <code>Subscriptor</code> class. In this
// destructor, there
// will then be a check whether the
// counter is really zero. If
// for the programmer to do so
// herself. The class that
// actually does all this is
- // called ``SmartPointer'' and
+ // called <code>SmartPointer</code> and
// takes as template parameter
// the data type of the object
// which it shall point to. The
// latter type may be any class,
// as long as it is derived from
- // the ``Subscriptor'' class.
+ // the <code>Subscriptor</code> class.
//
// In the present example program, we
// want to protect the finite element
// object from the situation that for
// some reason the finite element pointed
// to is destroyed while still in use. We
- // therefore use a ``SmartPointer'' to
+ // therefore use a <code>SmartPointer</code> to
// the finite element object; since the
// finite element object is actually
// never changed in our computations, we
- // pass a ``const FiniteElement@<dim@>'' as
+ // pass a <code>const FiniteElement@<dim@></code> as
// template argument to the
- // ``SmartPointer'' class. Note that the
+ // <code>SmartPointer</code> class. Note that the
// pointer so declared is assigned at
// construction time of the solve object,
// and destroyed upon destruction, so the
// lock on the destruction of the finite
// element object extends throughout the
- // lifetime of this ``HelmholtzProblem''
+ // lifetime of this <code>HelmholtzProblem</code>
// object.
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
// (like the number of cells, or the L2
// error of the numerical solution) will
// be generated and later printed. The
- // ``TableHandler'' can be used to
+ // <code>TableHandler</code> can be used to
// collect all this data and to output it
// at the end of the run as a table in a
// simple text or in LaTeX
// format. Here we don't only use the
- // ``TableHandler'' but we use the
- // derived class ``ConvergenceTable''
+ // <code>TableHandler</code> but we use the
+ // derived class <code>ConvergenceTable</code>
// that additionally evaluates rates of
// convergence:
ConvergenceTable convergence_table;
};
- // @sect3{The ``HelmholtzProblem'' class implementation}
+ // @sect3{The <code>HelmholtzProblem</code> class implementation}
// @sect4{HelmholtzProblem::HelmholtzProblem}
// bi-quadratic elements and therefore have
// to use sufficiently accurate quadrature
// formula. In addition, we need to compute
- // integrals over faces, i.e. ``dim-1''
+ // integrals over faces, i.e. <code>dim-1</code>
// dimensional objects. The declaration of a
// face quadrature formula is then
// straightforward:
// cell) to evaluate the right hand
// side function. The object we use
// to get at this information is
- // the ``FEValues'' class discussed
+ // the <code>FEValues</code> class discussed
// previously.
//
// For the face integrals, we only
// from the exact solution object
// (see below). The class that gives
// us this information is called
- // ``FEFaceValues'':
+ // <code>FEFaceValues</code>:
FEValues<dim> fe_values (*fe, quadrature_formula,
update_values | update_gradients |
update_q_points | update_JxW_values);
// the right hand side object are only
// querying data, never changing the
// object. We can therefore declare it
- // ``const'':
+ // <code>const</code>:
const RightHandSide<dim> right_hand_side;
std::vector<double> rhs_values (n_q_points);
// nonzero. To this end, we
// loop over all faces and
// check whether its boundary
- // indicator equals ``1'',
+ // indicator equals <code>1</code>,
// which is the value that we
// have assigned to that
// portions of the boundary
// composing Gamma2 in the
- // ``run()'' function further
+ // <code>run()</code> function further
// below. (The
// default value of boundary
- // indicators is ``0'', so faces
+ // indicators is <code>0</code>, so faces
// can only have an indicator
- // equal to ``1'' if we have
+ // equal to <code>1</code> if we have
// explicitly set it.)
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary()
// computation of the
// contour integral. This
// is done using the
- // ``reinit'' function
+ // <code>reinit</code> function
// which we already know
- // from the ``FEValue''
+ // from the <code>FEValue</code>
// class:
fe_face_values.reinit (cell, face);
// vector to the face at the
// present quadrature point
// obtained from the
- // ``fe_face_values''
+ // <code>fe_face_values</code>
// object. This is then used to
// compute the additional
// contribution of this face to
// we interpolate boundary values
// (denoted by the second parameter
// to
- // ``interpolate_boundary_values'')
+ // <code>interpolate_boundary_values</code>)
// does not represent the whole
// boundary any more. Rather, it is
// that portion of the boundary
//
// At the end of the switch, we have a
// default case that looks slightly strange:
- // an ``Assert'' statement with a ``false''
- // condition. Since the ``Assert'' macro
+ // an <code>Assert</code> statement with a <code>false</code>
+ // condition. Since the <code>Assert</code> macro
// raises an error whenever the condition is
// false, this means that whenever we hit
// this statement the program will be
// the difference between computed
// numerical solution and the
// continuous solution (described
- // by the ``Solution'' class
+ // by the <code>Solution</code> class
// defined at the top of this
// file), we first need a vector
// that will hold the norm of the
// accuracy with 16 digits is not
// so important for these
// quantities, we save some memory
- // by using ``float'' instead of
- // ``double'' values.
+ // by using <code>float</code> instead of
+ // <code>double</code> values.
//
// The next step is to use a function
// from the library which computes the
// cell, and taking the square root
// of that value. This is
// equivalent to taking the l2
- // (lower case ``l'') norm of the
+ // (lower case <code>l</code>) norm of the
// vector of norms on each cell:
Vector<float> difference_per_cell (triangulation.n_active_cells());
VectorTools::integrate_difference (dof_handler,
// By same procedure we get the H1
// semi-norm. We re-use the
- // ``difference_per_cell'' vector since it
+ // <code>difference_per_cell</code> vector since it
// is no longer used after computing the
- // ``L2_error'' variable above.
+ // <code>L2_error</code> variable above.
VectorTools::integrate_difference (dof_handler,
solution,
Solution<dim>(),
// by iterating the trapezoidal
// rule five times in each space
// direction. Note that the
- // constructor of the ``QIterated''
+ // constructor of the <code>QIterated</code>
// class takes a one-dimensional
// quadrature rule and a number
// that tells it how often it shall
// maximum value over all cell-wise
// entries, an operation that is
// conveniently done using the
- // ``Vector@<float@>::linfty'' function:
+ // <code>Vector@<float@>::linfty</code> function:
const QTrapez<1> q_trapez;
const QIterated<dim> q_iterated (q_trapez, 5);
VectorTools::integrate_difference (dof_handler,
// computed, we finally write some
// output. In addition, we add the
// important data to the
- // ``TableHandler'' by specifying
+ // <code>TableHandler</code> by specifying
// the key of the column and the value.
// Note that it is not necessary to
// define column keys beforehand -- it is
// @sect4{HelmholtzProblem::run}
// As in previous example programs,
- // the ``run'' function controls the
+ // the <code>run</code> function controls the
// flow of execution. The basic
// layout is as in previous examples:
// an outer loop over successively
// For this, we will use the
// following convention: Faces
// belonging to Gamma1 will have the
- // boundary indicator ``0'' (which is
+ // boundary indicator <code>0</code> (which is
// the default, so we don't have to
// set it explicitely), and faces
- // belonging to Gamma2 will use ``1''
+ // belonging to Gamma2 will use <code>1</code>
// as boundary indicator. To set
// these values, we loop over all
// cells, then over all faces of a
// part of the boundary that we want
// to denote by Gamma2, and if so set
// its boundary indicator to
- // ``1''. For the present program, we
+ // <code>1</code>. For the present program, we
// consider the left and bottom
// boundaries as Gamma2. We determine
// whether a face is part of that
// statements which we have already
// discussed in previous examples. The
// first step is to generate a suitable
- // filename (called ``gmv_filename'' here,
+ // filename (called <code>gmv_filename</code> here,
// since we want to output data in GMV
// format; we add the prefix to distinguish
// the filename from that used for other
// end, the finite element base class
// stores the maximal polynomial degree of
// shape functions in each coordinate
- // variable as a variable ``degree'', and
+ // variable as a variable <code>degree</code>, and
// we use for the switch statement (note
// that the polynomial degree of bilinear
// shape functions is really 2, since they
- // contain the term ``x*y''; however, the
+ // contain the term <code>x*y</code>; however, the
// polynomial degree in each coordinate
// variable is still only 1). We again use
// the same defensive programming technique
// to safeguard against the case that the
// polynomial degree has an unexpected
- // value, using the ``Assert (false,
- // ExcNotImplemented())'' idiom in the
+ // value, using the <code>Assert (false,
+ // ExcNotImplemented())</code> idiom in the
// default branch of the switch statement:
switch (fe->degree)
{
//
// In order to allow writing more
// than one sub-cell per actual
- // cell, the ``build_patches''
+ // cell, the <code>build_patches</code>
// function accepts a parameter
- // (the default is ``1'', which is
+ // (the default is <code>1</code>, which is
// why you haven't seen this
// parameter in previous
// examples). This parameter
// per space direction each cell
// shall be subdivided for
// output. For example, if you give
- // ``2'', this leads to 4 cells in
+ // <code>2</code>, this leads to 4 cells in
// 2D and 8 cells in 3D. For
// quadratic elements, two
// sub-cells per space direction is
// obviously the right choice, so
// this is what we choose. In
// general, for elements of
- // polynomial order ``q'', we use
- // ``q'' subdivisions, and the
+ // polynomial order <code>q</code>, we use
+ // <code>q</code> subdivisions, and the
// order of the elements is
// determined in the same way as
// above.
// After graphical output, we would also
// like to generate tables from the error
// computations we have done in
- // ``process_solution''. There, we have
+ // <code>process_solution</code>. There, we have
// filled a table object with the number of
// cells for each refinement step as well
// as the errors in different norms.
// fixed point notation. However, for
// columns one would like to see in
// scientific notation another function
- // call sets the ``scientific_flag'' to
- // ``true'', leading to floating point
+ // call sets the <code>scientific_flag</code> to
+ // <code>true</code>, leading to floating point
// representation of numbers.
convergence_table.set_precision("L2", 3);
convergence_table.set_precision("H1", 3);
// For the output of a table into a LaTeX
// file, the default captions of the
// columns are the keys given as argument
- // to the ``add_value'' functions. To have
+ // to the <code>add_value</code> functions. To have
// TeX captions that differ from the
// default ones you can specify them by the
// following function calls.
// After this, we can finally write the
// table to the standard output stream
- // ``std::cout'' (after one extra empty
+ // <code>std::cout</code> (after one extra empty
// line, to make things look
// prettier). Note, that the output in text
// format is quite simple and that
// output the convergence
// rates. This may be done by the
// functionality the
- // ``ConvergenceTable'' offers over
+ // <code>ConvergenceTable</code> offers over
// the regular
- // ``TableHandler''. However, we do
+ // <code>TableHandler</code>. However, we do
// it only for global refinement,
// since for adaptive refinement
// the determination of something
convergence_table.set_column_order (new_order);
// For everything that happened
- // to the ``ConvergenceTable''
+ // to the <code>ConvergenceTable</code>
// until this point, it would
// have been sufficient to use
// a simple
- // ``TableHandler''. Indeed, the
- // ``ConvergenceTable'' is
+ // <code>TableHandler</code>. Indeed, the
+ // <code>ConvergenceTable</code> is
// derived from the
- // ``TableHandler'' but it offers
+ // <code>TableHandler</code> but it offers
// the additional functionality
// of automatically evaluating
// convergence rates. For
// order to destroy the
// respective objects (i.e. the
// finite element and the
- // ``HelmholtzProblem'' object)
+ // <code>HelmholtzProblem</code> object)
// at the end of the block and
// before we go to the next
// run. This avoids conflicts
// is released immediately
// after one of the three runs
// has finished, and not only
- // at the end of the ``try''
+ // at the end of the <code>try</code>
// block.
{
std::cout << "Solving with Q1 elements, adaptive refinement" << std::endl
#include <iostream>
- // @sect3{The ``ElasticProblem'' class template}
+ // @sect3{The <code>ElasticProblem</code> class template}
// The main class is, except for its
// name, almost unchanged with
// respect to the step-6 example.
//
// The only change is the use of a
- // different class for the ``fe''
+ // different class for the <code>fe</code>
// variable: Instead of a concrete
// finite element class such as
- // ``FE_Q'', we now use a more
- // generic one, ``FESystem''. In
- // fact, ``FESystem'' is not really a
+ // <code>FE_Q</code>, we now use a more
+ // generic one, <code>FESystem</code>. In
+ // fact, <code>FESystem</code> is not really a
// finite element itself in that it
// does not implement shape functions
// of its own. Rather, it is a class
// one vector-valued finite
// element. In our case, we will
// compose the vector-valued element
- // of ``FE_Q(1)'' objects, as shown
+ // of <code>FE_Q(1)</code> objects, as shown
// below in the constructor of this
// class.
template <int dim>
// The next change is that we
// want a replacement for the
- // ``value'' function of the
+ // <code>value</code> function of the
// previous examples. There, a
- // second parameter ``component''
+ // second parameter <code>component</code>
// was given, which denoted which
// component was requested. Here,
// we implement a function that
// once, in the second argument
// of the function. The obvious
// name for such a replacement
- // function is ``vector_value''.
+ // function is <code>vector_value</code>.
//
// Secondly, in analogy to the
- // ``value_list'' function, there
+ // <code>value_list</code> function, there
// is a function
- // ``vector_value_list'', which
+ // <code>vector_value_list</code>, which
// returns the values of the
// vector-valued function at
// several points at once:
// right hand side class. As said
// above, it only passes down to the
// base class the number of
- // components, which is ``dim'' in
+ // components, which is <code>dim</code> in
// the present case (one force
- // component in each of the ``dim''
+ // component in each of the <code>dim</code>
// space directions).
//
// Some people would have moved the
// Next the function that returns
// the whole vector of values at the
- // point ``p'' at once.
+ // point <code>p</code> at once.
//
// To prevent cases where the return
// vector has not previously been set
// not be removed if we can't rely on
// the assumption that the vector
// already has the correct size; this
- // is in contract to the ``Assert''
+ // is in contract to the <code>Assert</code>
// call that is completely removed if
// the program is compiled in
// optimized mode.
// two objects that denote the
// centers of these areas. Note
// that upon construction of the
- // ``Point'' objects, all
+ // <code>Point</code> objects, all
// components are set to zero.
Point<dim> point_1, point_2;
point_1(0) = 0.5;
point_2(0) = -0.5;
- // If now the point ``p'' is in a
+ // If now the point <code>p</code> is in a
// circle (sphere) of radius 0.2
// around one of these points, then
// set the force in x-direction to
else
values(0) = 0;
- // Likewise, if ``p'' is in the
+ // Likewise, if <code>p</code> is in the
// vicinity of the origin, then set
// the y-force to 1, otherwise to
// zero:
// points. In one of the previous
// examples, we have explained why
// the
- // ``value_list''/``vector_value_list''
+ // <code>value_list</code>/<code>vector_value_list</code>
// function had been introduced: to
// prevent us from calling virtual
// functions too frequently. On the
//
// We can prevent this situation by
// calling
- // ``RightHandSide@<dim@>::vector_valued''
+ // <code>RightHandSide@<dim@>::vector_valued</code>
// on each point in the input
// list. Note that by giving the
// full name of the function,
// and not to use the virtual
// function call mechanism that
// would be used if we had just
- // called ``vector_value''. This is
+ // called <code>vector_value</code>. This is
// important, since the compiler
// generally can't make any
// assumptions which function is
// inline above function into the
// present location. (Note that we
// have declared the
- // ``vector_value'' function above
- // ``inline'', though modern
+ // <code>vector_value</code> function above
+ // <code>inline</code>, though modern
// compilers are also able to
// inline functions even if they
// have not been declared as
// functions in the same way. Using
// this forwarding mechanism, we
// only have to change a single
- // place (the ``vector_value''
+ // place (the <code>vector_value</code>
// function), and the second place
- // (the ``vector_value_list''
+ // (the <code>vector_value_list</code>
// function) will always be
// consistent with it. At the same
// time, using virtual function
- // @sect3{The ``ElasticProblem'' class implementation}
+ // @sect3{The <code>ElasticProblem</code> class implementation}
// @sect4{ElasticProblem::ElasticProblem}
// would like to stack together
// equals the number of components
// the solution function has, which
- // is ``dim'' since we consider
+ // is <code>dim</code> since we consider
// displacement in each space
- // direction. The ``FESystem'' class
+ // direction. The <code>FESystem</code> class
// can handle this: we pass it the
// finite element of which we would
// like to compose the system of, and
dof_handler (triangulation),
fe (FE_Q<dim>(1), dim)
{}
- // In fact, the ``FESystem'' class
+ // In fact, the <code>FESystem</code> class
// has several more constructors
// which can perform more complex
// operations than just stacking
// Setting up the system of equations
// is identitical to the function
// used in the step-6 example. The
- // ``DoFHandler'' class and all other
+ // <code>DoFHandler</code> class and all other
// classes used here are fully aware
// that the finite element we want to
// use is vector-valued, and take
// are the same as before, however:
// setting up a suitable quadrature
// formula, initializing an
- // ``FEValues'' object for the
+ // <code>FEValues</code> object for the
// (vector-valued) finite element we
// use as well as the quadrature
// object, and declaring a number of
// auxiliary arrays. In addition, we
// declare the ever same two
- // abbreviations: ``n_q_points'' and
- // ``dofs_per_cell''. The number of
+ // abbreviations: <code>n_q_points</code> and
+ // <code>dofs_per_cell</code>. The number of
// degrees of freedom per cell we now
// obviously ask from the composed
// finite element rather than from
// the underlying scalar Q1
- // element. Here, it is ``dim'' times
+ // element. Here, it is <code>dim</code> times
// the number of degrees of freedom
// per cell of the Q1 element, though
// this is not explicit knowledge we
// we now have a vector-valued
// right hand side, which is why
// the data type of the
- // ``rhs_values'' array is
+ // <code>rhs_values</code> array is
// changed. We initialize it by
- // ``n_q_points'' elements, each of
- // which is a ``Vector@<double@>''
- // with ``dim'' elements.
+ // <code>n_q_points</code> elements, each of
+ // which is a <code>Vector@<double@></code>
+ // with <code>dim</code> elements.
RightHandSide<dim> right_hand_side;
std::vector<Vector<double> > rhs_values (n_q_points,
Vector<double>(dim));
// example. One of the few
// comments in place is that we
// can compute the number
- // ``comp(i)'', i.e. the index
+ // <code>comp(i)</code>, i.e. the index
// of the only nonzero vector
// component of shape function
- // ``i'' using the
- // ``fe.system_to_component_index(i).first''
+ // <code>i</code> using the
+ // <code>fe.system_to_component_index(i).first</code>
// function call below.
//
// (By accessing the
- // ``first'' variable of
+ // <code>first</code> variable of
// the return value of the
- // ``system_to_component_index''
+ // <code>system_to_component_index</code>
// function, you might
// already have guessed
// that there is more in
// it. In fact, the
// function returns a
- // ``std::pair@<unsigned int,
- // unsigned int@>'', of
+ // <code>std::pair@<unsigned int,
+ // unsigned int@></code>, of
// which the first element
- // is ``comp(i)'' and the
+ // is <code>comp(i)</code> and the
// second is the value
- // ``base(i)'' also noted
+ // <code>base(i)</code> also noted
// in the introduction, i.e.
// the index
// of this shape function
// within all the shape
// functions that are nonzero
// in this component,
- // i.e. ``base(i)'' in the
+ // i.e. <code>base(i)</code> in the
// diction of the
// introduction. This is not a
// number that we are usually
// + (mu d_i u_j,
// d_j v_i).
// Note that
- // ``shape_grad(i,q_point)''
+ // <code>shape_grad(i,q_point)</code>
// returns the
// gradient of
// the only
// point
// q_point. The
// component
- // ``comp(i)'' of
+ // <code>comp(i)</code> of
// the gradient,
// which is the
// derivative of
// modification: since the solution
// function is vector-valued, so
// need to be the boundary
- // values. The ``ZeroFunction''
+ // values. The <code>ZeroFunction</code>
// constructor accepts a parameter
// that tells it that it shall
// represent a vector valued,
// many components. By default,
// this parameter is equal to one,
// in which case the
- // ``ZeroFunction'' object would
+ // <code>ZeroFunction</code> object would
// represent a scalar
// function. Since the solution
- // vector has ``dim'' components,
- // we need to pass ``dim'' as
+ // vector has <code>dim</code> components,
+ // we need to pass <code>dim</code> as
// number of components to the zero
// function as well.
std::map<unsigned int,double> boundary_values;
// been shown in previous examples
// already. The only difference is
// that the solution function is
- // vector valued. The ``DataOut''
+ // vector valued. The <code>DataOut</code>
// class takes care of this
// automatically, but we have to give
// each component of the solution
// number of components is the same
// as the number of dimensions we
// are working in, the following
- // ``switch'' statement is used.
+ // <code>switch</code> statement is used.
//
// We note that some graphics
// programs have restriction as to
// the program die if we run upon a
// case which we did not
// consider. Remember that the
- // ``Assert'' macro generates an
+ // <code>Assert</code> macro generates an
// exception if the condition in
// the first parameter is not
// satisfied. Of course, the
- // condition ``false'' can never be
+ // condition <code>false</code> can never be
// satisfied, so the program will
// always abort whenever it gets to
// the default statement:
// @sect4{ElasticProblem::run}
- // The ``run'' function does the same
+ // The <code>run</code> function does the same
// things as in step-6, for
// example. This time, we use the
// square [-1,1]^d as domain, and we
// starting the first iteration.
//
// The reason is the following: we
- // use the ``Gauss'' quadrature
+ // use the <code>Gauss</code> quadrature
// formula with two points in each
// direction for integration of the
// right hand side; that means that
// The unfortunate thing is that if
// the discrete solution is constant,
// then the error indicators computed
- // by the ``KellyErrorEstimator''
+ // by the <code>KellyErrorEstimator</code>
// class are zero for each cell as
// well, and the call to
- // ``refine_and_coarsen_fixed_number''
- // on the ``triangulation'' object
+ // <code>refine_and_coarsen_fixed_number</code>
+ // on the <code>triangulation</code> object
// will not flag any cells for
// refinement (why should it if the
// indicated error is zero for each
// needs to be able to see the right
// hand side. Thus, we refine twice
// globally. (Note that the
- // ``refine_global'' function is not
- // part of the ``GridRefinement''
+ // <code>refine_global</code> function is not
+ // part of the <code>GridRefinement</code>
// class in which
- // ``refine_and_coarsen_fixed_number''
+ // <code>refine_and_coarsen_fixed_number</code>
// is declared, for example. The
// reason is first that it is not an
// algorithm that computed refinement
// importantly that it actually
// performs the refinement, in
// contrast to the functions in
- // ``GridRefinement'' that only flag
+ // <code>GridRefinement</code> that only flag
// cells without actually refining
// the grid.)
template <int dim>
}
}
- // @sect3{The ``main'' function}
+ // @sect3{The <code>main</code> function}
// The main function is again exactly
// like in step-6 (apart from the
// functions are declared which we
// need to start new threads and to
// wait for threads to return
- // (i.e. the ``Thread'' class
- // and the ``spawn'' functions). The
+ // (i.e. the <code>Thread</code> class
+ // and the <code>spawn</code> functions). The
// second file has a class
- // ``MultithreadInfo'' (and a global
- // object ``multithread_info'' of
+ // <code>MultithreadInfo</code> (and a global
+ // object <code>multithread_info</code> of
// that type) which can be used to
// query the number of processors in
// your system, which is often useful
#include <base/multithread_info.h>
// The next new include file declares
- // a base class ``TensorFunction''
- // not unlike the ``Function'' class,
+ // a base class <code>TensorFunction</code>
+ // not unlike the <code>Function</code> class,
// but with the difference that the
// return value is tensor-valued
// rather than scalar of
// local contributions of a cell
// to the global matrix at the
// same time. This is done using
- // a ``Mutex'', which is an
+ // a <code>Mutex</code>, which is an
// object that can be owned by
// only one thread at a time. If
// a thread wants to write to the
// was not compiled to support
// multithreading (which you have
// to specify at the time you
- // call the ``./configure''
+ // call the <code>./configure</code>
// script in the top-level
// directory), then a dummy the
// actual data type of the
// typedef
- // ``Threads::ThreadMutex'' is a
+ // <code>Threads::ThreadMutex</code> is a
// class that provides all the
// functions needed for a mutex,
// but does nothing when they are
// vector field with as many compents
// as there are space dimensions. One
// could now use a class derived from
- // the ``Function'' base class, as we
+ // the <code>Function</code> base class, as we
// have done for boundary values and
// coefficients in previous examples,
// but there is another possibility
// in the library, namely a base
// class that describes tensor valued
// functions. In contrast to the
- // usual ``Function'' objects, we
+ // usual <code>Function</code> objects, we
// provide the compiler with
// knowledge on the size of the
// objects of the return type. This
// simple for usual vector-valued
// functions where memory has to be
// allocated on the heap (thus, the
- // ``Function::vector_value''
+ // <code>Function::vector_value</code>
// function has to be given the
// address of an object into which
// the result is to be written, in
// applications, to be honest...
//
// The interface of the
- // ``TensorFunction'' class is
+ // <code>TensorFunction</code> class is
// relatively close to that of the
- // ``Function'' class, so there is
+ // <code>Function</code> class, so there is
// probably no need to comment in
// detail the following declaration:
template <int dim>
// reasonable. The format is
// basically as follows: use the
// name of one of the macros
- // ``DeclExceptionN'', where
- // ``N'' denotes the number of
+ // <code>DeclExceptionN</code>, where
+ // <code>N</code> denotes the number of
// additional parameters which
// the exception object shall
// take. In this case, as we want
// the sizes of two vectors
// differ, we need two arguments,
// so we use
- // ``DeclException2''. The first
+ // <code>DeclException2</code>. The first
// parameter then describes the
// name of the exception, while
// the following declare the data
// types of the parameters. The
// last argument is a sequence of
// output directives that will be
- // piped into the ``std::cerr''
+ // piped into the <code>std::cerr</code>
// object, thus the strange
- // format with the leading ``@<@<''
+ // format with the leading <code>@<@<</code>
// operator and the like. Note
// that we can access the
// parameters which are passed to
// the exception upon
// construction (i.e. within the
- // ``Assert'' call) by using the
- // names ``arg1'' through
- // ``argN'', where ``N'' is the
+ // <code>Assert</code> call) by using the
+ // names <code>arg1</code> through
+ // <code>argN</code>, where <code>N</code> is the
// number of arguments as defined
// by the use of the respective
- // macro ``DeclExceptionN''.
+ // macro <code>DeclExceptionN</code>.
//
// To learn how the preprocessor
// expands this macro into actual
// Besides the advection field, we
// need two functions describing the
- // source terms (``right hand side'')
+ // source terms (<code>right hand side</code>)
// and the boundary values. First for
// the right hand side, which follows
// the same pattern as in previous
// constant function in the vicinity
// of a source point, which we denote
// by the constant static variable
- // ``center_point''. We set the
+ // <code>center_point</code>. We set the
// values of this center using the
// same template tricks as we have
// shown in the step-7 example
// has been shown previously,
// including the way to avoid virtual
// function calls in the
- // ``value_list'' function.
+ // <code>value_list</code> function.
template <int dim>
class RightHandSide : public Function<dim>
{
// The only new thing here is that we
// check for the value of the
- // ``component'' parameter. As this
+ // <code>component</code> parameter. As this
// is a scalar function, it is
// obvious that it only makes sense
// if the desired component has the
// index zero, so we assert that this
// is indeed the
- // case. ``ExcIndexRange'' is a
+ // case. <code>ExcIndexRange</code> is a
// global predefined exception
// (probably the one most often used,
// we therefore made it global
// Finally for the boundary values,
// which is just another class
- // derived from the ``Function'' base
+ // derived from the <code>Function</code> base
// class:
template <int dim>
class BoundaryValues : public Function<dim>
// power of the mesh size, as
// described in the introduction.
// This class is a simple version of
- // the ``DerivativeApproximation''
+ // the <code>DerivativeApproximation</code>
// class in the library, that uses
// similar techniques to obtain
// finite difference approximations
//
// The
// class has one public static
- // function ``estimate'' that is
+ // function <code>estimate</code> that is
// called to compute a vector of
// error indicators, and one private
// function that does the actual work
// functions or variables, so this is
// not really a class, but rather
// serves the purpose of a
- // ``namespace'' in C++. The reason
+ // <code>namespace</code> in C++. The reason
// that we chose a class over a
// namespace is that this way we can
// declare functions that are
// argument.
//
// Finally note that the
- // ``IndexInterval'' typedef is
+ // <code>IndexInterval</code> typedef is
// introduced as a convenient
// abbreviation for an otherwise
// lengthy type name.
// Now for the implementation of the
// main class. Constructor,
// destructor and the function
- // ``setup_system'' follow the same
+ // <code>setup_system</code> follow the same
// pattern that was used previously,
// so we need not comment on these
// three function:
// we were to use this information,
// we could use the value of the
// global variable
- // ``multithread_info.n_cpus'',
+ // <code>multithread_info.n_cpus</code>,
// which is determined at start-up
// time of your program
// automatically. (Note that if the
// systems assign roughly the same
// CPU ressources to all threads
// presently running. For this
- // reason, the ``MultithreadInfo''
+ // reason, the <code>MultithreadInfo</code>
// class contains a read-write
- // variable ``n_default_threads''
- // which is set to ``n_cpus'' by
+ // variable <code>n_default_threads</code>
+ // which is set to <code>n_cpus</code> by
// default, but can be set to
// another value. This variable is
// also queried by functions inside
// capable of keeping track of the
// threads we created, and allows
// us to wait until they all have
- // finished (to ``join'' them in
+ // finished (to <code>join</code> them in
// the language of threads). The
- // ``Threads::ThreadGroup'' class
+ // <code>Threads::ThreadGroup</code> class
// does this, which is basically
// just a container for objects of
- // type ``Threads::Thread'' that
+ // type <code>Threads::Thread</code> that
// represent a single thread;
- // ``Threads::Thread'' is what the
- // ``spawn'' function below will
+ // <code>Threads::Thread</code> is what the
+ // <code>spawn</code> function below will
// return when we start a new
// thread.
//
- // Note that both ``ThreadGroup''
- // and ``Thread'' have a template
+ // Note that both <code>ThreadGroup</code>
+ // and <code>Thread</code> have a template
// argument that represents the
// return type of the function
// being called on a separate
// thread. Since most of the
// functions that we will call on
// different threads have return
- // type ``void'', the template
+ // type <code>void</code>, the template
// argument has a default value
- // ``void'', so that in that case
+ // <code>void</code>, so that in that case
// it can be omitted. (However, you
// still need to write the angle
// brackets, even if they are
//
// If you did not configure for
// multi-threading, then the
- // ``spawn'' function that is
+ // <code>spawn</code> function that is
// supposed to start a new thread
// in parallel only executes the
// function which should be run in
// parallel, waits for it to return
// (i.e. the function is executed
// sequentially), and puts the
- // return value into the ``Thread''
+ // return value into the <code>Thread</code>
// object. Likewise, the function
- // ``join'' that is supposed to
+ // <code>join</code> that is supposed to
// wait for all spawned threads to
// return, returns immediately, as
// there can't be threads running.
// splitting a range of cells is a
// rather common task when using
// multi-threading, there is a
- // function in the ``Threads''
+ // function in the <code>Threads</code>
// namespace that does exactly
// this. In fact, it does this not
// only for a range of cell
// iterators, but for iterators in
// general, so you could use it for
- // ``std::vector@<T@>::iterator'' or
+ // <code>std::vector@<T@>::iterator</code> or
// usual pointers as well.
//
// The function returns a vector of
// present case, however, the data
// types of the two first
// parameters differ
- // (``begin_active'' returns an
- // ``active_iterator'', while
- // ``end'' returns a
- // ``raw_iterator''), and in this
+ // (<code>begin_active</code> returns an
+ // <code>active_iterator</code>, while
+ // <code>end</code> returns a
+ // <code>raw_iterator</code>), and in this
// case the C++ language requires
// us to specify the template type
// explicitely. For brevity, we
// which is available online as
// well. Suffice it to say that we
// spawn a new thread that calls
- // the ``assemble_system_interval''
+ // the <code>assemble_system_interval</code>
// function on the present object
- // (the ``this'' pointer), with the
+ // (the <code>this</code> pointer), with the
// arguments following in the
// second set of parentheses passed
- // as parameters. The ``spawn''
+ // as parameters. The <code>spawn</code>
// function return an object of
- // type ``Threads::Thread'', which
- // we put into the ``threads''
+ // type <code>Threads::Thread</code>, which
+ // we put into the <code>threads</code>
// container. If a thread exits,
// the return value of the function
// being called is put into a place
// such that the thread objects can
// access it using their
- // ``return_value'' function; since
+ // <code>return_value</code> function; since
// the function we call doesn't
// have a return value, this does
// not apply here. Note that you
// right hand side are
// assemblesd. Waiting for all the
// threads to finish can be done
- // using the ``joint_all'' function
- // in the ``ThreadGroup''
+ // using the <code>joint_all</code> function
+ // in the <code>ThreadGroup</code>
// container, which just calls
- // ``join'' on each of the thread
+ // <code>join</code> on each of the thread
// objects it stores.
//
// Again, if the library was not
// Now, this is the function that
// does the actual work. It is not
// very different from the
- // ``assemble_system'' functions of
+ // <code>assemble_system</code> functions of
// previous example programs, so we
// will again only comment on the
// differences. The mathematical
QGauss<dim-1> face_quadrature_formula(2);
// Finally, we need objects of type
- // ``FEValues'' and
- // ``FEFaceValues''. For the cell
+ // <code>FEValues</code> and
+ // <code>FEFaceValues</code>. For the cell
// terms we need the values and
// gradients of the shape
// functions, the quadrature points
cell_rhs = 0;
// ... then initialize
- // the ``FEValues'' object...
+ // the <code>FEValues</code> object...
fe_values.reinit (cell);
// ... obtain the values of
// as well. Of course, the
// bilinear form only contains
// contributions from the
- // ``inflow'' part of the
+ // <code>inflow</code> part of the
// boundary, but to find out
// whether a certain part of a
// face of the present cell is
// thread operates on these
// objects at a time, we have
// to lock it. This is done
- // using a ``Mutex'', which is
- // short for ``mutually
- // exclusive'': a thread that
+ // using a <code>Mutex</code>, which is
+ // short for <code>mutually
+ // exclusive</code>: a thread that
// wants to write to the global
// objects acquires this lock,
// but has to wait if it is
// can't be parallel threads
// and there is no need to
// synchronize. Thus, the
- // ``lock'' and ``release''
+ // <code>lock</code> and <code>release</code>
// functions are no-ops,
// i.e. they return without
// doing anything.
// all threads execute member
// functions of the same
// object, they have the same
- // ``this'' pointer and
+ // <code>this</code> pointer and
// therefore also operate on
- // the same ``lock''.
+ // the same <code>lock</code>.
};
}
// described in the introduction. The
// respective computations are made
// in the class
- // ``GradientEstimation''. The only
+ // <code>GradientEstimation</code>. The only
// difference to previous examples is
// that we refine a little more
// aggressively (0.5 instead of 0.3
// @sect3{GradientEstimation class implementation}
// Now for the implementation of the
- // ``GradientEstimation'' class. The
+ // <code>GradientEstimation</code> class. The
// first function does not much
// except for delegating work to the
// other function:
// cells into chunks of equal
// size. Just as we have used the
// function
- // ``Threads::split_range'' when
+ // <code>Threads::split_range</code> when
// assembling above, there is a
// function that computes intervals
// of roughly equal size from a
n_threads);
// In the same way as before, we
- // use a ``Threads::ThreadGroup''
+ // use a <code>Threads::ThreadGroup</code>
// object to collect the descriptor
// objects of different
// threads. Note that as the
// function called is not a member
// function, but rather a static
// function, we need not (and can
- // not) pass a ``this'' pointer to
- // the ``spawn'' function in this
+ // not) pass a <code>this</code> pointer to
+ // the <code>spawn</code> function in this
// case.
//
// Taking pointers to templated
// quite frequently that we can't
// directly insert taking the
// address of a function in the
- // call to ``encapsulate'' for one
+ // call to <code>encapsulate</code> for one
// or the other compiler, but have
// to take a temporary variable for
// that purpose. Here, in this
- // case, Compaq's ``cxx'' compiler
+ // case, Compaq's <code>cxx</code> compiler
// choked on the code so we use
// this workaround with the
// function pointer:
threads.join_all ();
// Note that if the value of the
// variable
- // ``multithread_info.n_default_threads''
+ // <code>multithread_info.n_default_threads</code>
// was one, or if the library was
// not configured to use threads,
// then the sequence of commands
// above reduced to a complicated
// way to simply call the
- // ``estimate_interval'' function
+ // <code>estimate_interval</code> function
// with the whole range of cells to
// work on. However, using the way
// above, we are able to write the
// of the cells. As usual with
// values of finite element
// functions, we use an object of
- // type ``FEValues'', and we use
+ // type <code>FEValues</code>, and we use
// (or mis-use in this case) the
// midpoint quadrature rule to get
// at the values at the
// center. Note that the
- // ``FEValues'' object only needs
+ // <code>FEValues</code> object only needs
// to compute the values at the
// centers, and the location of the
// quadrature points in real space
// in order to get at the vectors
- // ``y''.
+ // <code>y</code>.
QMidpoint<dim> midpoint_rule;
FEValues<dim> fe_midpoint_value (dof_handler.get_fe(),
midpoint_rule,
update_values | update_q_points);
// Then we need space foe the
- // tensor ``Y'', which is the sum
+ // tensor <code>Y</code>, which is the sum
// of outer products of the
// y-vectors.
Tensor<2,dim> Y;
// cell and advancing them using
// the given start and end
// index. Note that we can use the
- // ``advance'' function of the
+ // <code>advance</code> function of the
// standard C++ library, but that
// we have to cast the distance by
// which the iterator is to be
for (; cell!=endc; ++cell, ++error_on_this_cell)
{
// First initialize the
- // ``FEValues'' object, as well
- // as the ``Y'' tensor:
+ // <code>FEValues</code> object, as well
+ // as the <code>Y</code> tensor:
fe_midpoint_value.reinit (cell);
Y.clear ();
// to
// involuntarily
// exchange
- // ``n==1'' for
- // ``n==0'' or
+ // <code>n==1</code> for
+ // <code>n==0</code> or
// the like) and
// in the library
// (the
// Now loop over all active neighbors
// and collect the data we
// need. Allocate a vector just like
- // ``this_midpoint_value'' which we
+ // <code>this_midpoint_value</code> which we
// will use to store the value of the
// solution in the midpoint of the
// neighbor cell. We allocate it here
// thereon. Note that for
// this information we
// have to reinitialize the
- // ``FEValues'' object for
+ // <code>FEValues</code> object for
// the neighbor cell.
fe_midpoint_value.reinit (neighbor);
const Point<dim> neighbor_center = fe_midpoint_value.quadrature_point(0);
fe_midpoint_value.get_function_values (solution,
neighbor_midpoint_value);
- // Compute the vector ``y''
+ // Compute the vector <code>y</code>
// connecting the centers
// of the two cells. Note
// that as opposed to the
// introduction, we denote
- // by ``y'' the normalized
+ // by <code>y</code> the normalized
// difference vector, as
// this is the quantity
// used everywhere in the
// an approximation of the
// gradient for the present
// cell, then we need to have
- // passed over vectors ``y''
+ // passed over vectors <code>y</code>
// which span the whole space,
// otherwise we would not have
// all components of the
// reasonable to try to catch
// this error also in optimized
// mode. For this case, there
- // is the ``AssertThrow''
+ // is the <code>AssertThrow</code>
// macro: it checks the
// condition like the
- // ``Assert'' macro, but not
+ // <code>Assert</code> macro, but not
// only in debug mode; it then
// outputs an error message,
// but instead of terminating
// the program as in the case
- // of the ``Assert'' macro, the
+ // of the <code>Assert</code> macro, the
// exception is thrown using
- // the ``throw'' command of
+ // the <code>throw</code> command of
// C++. This way, one has the
// possibility to catch this
// error and take reasonable
// @sect3{Main function}
- // The ``main'' function is exactly
+ // The <code>main</code> function is exactly
// like in previous examples, with
// the only difference in the name of
// the main class that actually does
projects / dealii.git / commitdiff