* depending on the ConditionalOStream object being active (default)
* or not. The condition of this object can be changed by
* set_condition() and in the constructor. This class is used in the
- * @ref step_17 "step-17", @ref step_18 "step-18", @ref step_32 "step-32",
- * @ref step_33 "step-33", and @ref step_35 "step-35"
+ * step-17, step-18, step-32,
+ * step-33, and step-35
* tutorial programs.
*
* This class is mostly useful in parallel computations. Ordinarily, you would
* that are not signed to be omitted.
* </ol>
*
- * A detailed discussion of this class can also be found in the @ref step_7 "step-7" and
- * @ref step_13 "step-13" example programs.
+ * A detailed discussion of this class can also be found in the step-7 and
+ * step-13 example programs.
*
* @ingroup textoutput
* @author Ralf Hartmann, 1999
* generate a single output file from it, which may be again in
* intermediate format or any of the final formats. This latter option
* is most helpful for parallel programs: as demonstrated in the
- * @ref step_17 "step-17" example program, it is possible to let only one processor
+ * step-17 example program, it is possible to let only one processor
* generate the graphical output for the entire parallel program, but
* this can become vastly inefficient if many processors are involved,
* because the load is no longer balanced. The way out is to let each
* be merged, while this is almost impossible once the data has been
* written out in any of the supported established graphics formats.
*
- * This second use scenario is explained in some detail in the @ref step_18 "step-18"
+ * This second use scenario is explained in some detail in the step-18
* example program.
*
- * Both these applications are implemented in the @ref step_19 "step-19" example program.
+ * Both these applications are implemented in the step-19 example program.
* There, a slight complication is also explained: in order to read data back
* into this object, you have to know the template parameters for the space
* dimension which were used when writing the data. If this knowledge is
* (such as in a simple format converter), then it needs to be figured out at
* run time, even though the compiler already needs it at compile time. A way
* around using the DataOutBase::determine_intermediate_format_dimensions()
- * function is explained in @ref step_19 "step-19".
+ * function is explained in step-19.
*
* Note that the intermediate format is what its name suggests: a
* direct representation of internal data. It isn't standardized and
* patches.
*
* The use of this function is
- * demonstrated in @ref step_19 "step-19".
+ * demonstrated in step-19.
*/
void merge (const DataOutReader<dim,spacedim> &other);
* components are zero. It is especially useful as a weight function
* for <tt>VectorTools::integrate_difference</tt>, where it allows to
* integrate only one or a few vector components, rather than the
- * entire vector-valued solution. See the @ref step_20 "step-20"
+ * entire vector-valued solution. See the step-20
* tutorial program for a detailed explanation.
*
* @ingroup functions
* from the fparser.txt file that comes with the fparser library. Please refer
* also to that file both for clarifications on how this wrapper works, as
* well as for any issue regarding the licence that applies to this class.
- * This class is used in the @ref step_33 "step-33" and @ref step_36 "step-36"
+ * This class is used in the step-33 and step-36
* tutorial programs (the latter being much simpler to understand).
*
* By using this class you indicate that you accept the terms of the licence
*
* The ParameterHandler class is discussed in detail in the @ref
* step_19 "step-19" example program, and is used in more realistic
- * situations in @ref step_29 "step-29", @ref step_33 "step-33"
- * and @ref step_34 "step-34".
+ * situations in step-29, step-33
+ * and step-34.
*
* <h3>Declaring entries</h3>
*
/**
* Friendly interface to the FunctionParser class. This class is
* meant as a wrapper for the FunctionParser class. It is used in the
- * @ref step_34 "step-34" tutorial program.
+ * step-34 tutorial program.
*
* It provides two
* methods to declare and parse a ParameterHandler object and creates
* to prescribe the format and the captions the columns are written
* with in tex mode.
*
- * A detailed explanation of this class is also given in the @ref step_13 "step-13" tutorial
+ * A detailed explanation of this class is also given in the step-13 tutorial
* program.
*
*
* Alternatively, you can also restart the timer instead of resetting
* it. The times between successive calls to start()/ stop() will then be
* accumulated. The usage of this class is also explained in the
- * @ref step_12 "step-12" and @ref step_29 "step-29" tutorial programs.
+ * step-12 and step-29 tutorial programs.
*
* @note Implementation of this class is system
* dependent. Unfortunately, it does not work with multithreading
* object. In this case, we did a lot of other stuff, so that the time
* proportions of the functions we measured are far away from 100 precent.
*
- * See the @ref step_32 "step-32" tutorial program for usage of this class.
+ * See the step-32 tutorial program for usage of this class.
*
* @ingroup utilities
* @author M. Kronbichler, 2009.
* implementations of Trilinos, sets up the MPI communicator in case the
* programs are run in %parallel, and correctly terminates all processes
* when the destructor is called. An example usage of this class is shown
- * in the tutorial program @ref step_32 "step-32".
+ * in the tutorial program step-32.
*/
#ifdef DEAL_II_USE_TRILINOS
namespace Trilinos
* problem is to use an intermediate
* compressed sparsity pattern that only
* allocates memory on demand. Refer to
- * the @ref step_11 "step-11" example
+ * the step-11 example
* program on how to do this. The problem
* is also discussed in the documentation
* of the module on @ref Sparsity.
* sparsity patterns produced by the various algorithms when using the
* $Q_2^d\times Q_1$ element combination typically employed in the
* discretization of Stokes equations, when used on the mesh obtained
- * in @ref step_22 "step-22" after one adaptive mesh refinement in
+ * in step-22 after one adaptive mesh refinement in
* 3d. The space dimension together with the coupled finite element
* leads to a rather dense system matrix with, on average around 180
* nonzero entries per row. After applying each of the reordering
* right and bottom left.
*
* The goal of reordering strategies is to improve the
- * preconditioner. In @ref step_22 "step-22" we use a SparseILU to
+ * preconditioner. In step-22 we use a SparseILU to
* preconditioner for the velocity-velocity block at the top left. The
* quality of the preconditioner can then be measured by the number of
* CG iterations required to solve a linear system with this
* of iterations is only slightly less than with BOOST's Cuthill-McKee
* implementation, runtime is significantly less. This, again, may be due
* to cache effects. As a consequence, this is the algorithm best suited
- * to the testcase, and is in fact used in @ref step_22 "step-22".
+ * to the testcase, and is in fact used in step-22.
* </td>
* <td>
* Minimum degree enumeration as produced by calling the BOOST implementation
* DoFRenumbering namespace.
*
* This algorithm is used in
- * @ref step_22 "step-22".
+ * step-22.
*/
template <class DH>
void
* possible and results in a
* blocking of several components
* into one. This is discussed in
- * @ref step_22 "step-22". If you
+ * step-22. If you
* omit this argument, the same
* order as given by the finite
* element is used.
* function. No call to
* ConstraintMatrix::condense() is then
* necessary. This process is explained
- * in @ref step_27 "step-27".
+ * in step-27.
*
* In case the constraints are
* already taken care of in this
* pattern that corresponds to the
* subdomain_id for which it is
* responsible. This feature is
- * @ref step_32 "step-32".
+ * step-32.
*/
template <class DH, class SparsityPattern>
static
* pattern that corresponds to the
* subdomain_id for which it is
* responsible. This feature is
- * @ref step_32 "step-32".
+ * step-32.
*/
template <class DH, class SparsityPattern>
static
* case, an assertion is thrown.
*
* This function is used in the
- * @ref step_22 "step-22",
- * @ref step_31 "step-31", and
- * @ref step_32 "step-32" tutorial
+ * step-22,
+ * step-31, and
+ * step-32 tutorial
* programs.
*/
template <int dim, int spacedim>
* Make a constraint matrix for the
* constraints that result from zero
* boundary values. This function is used
- * in @ref step_36 "step-36", for
+ * in step-36, for
* example.
*/
template <int dim, int spacedim, template <int, int> class DH>
*
* The use of this function is
* explained extensively in the
- * @ref step_8 "step-8" and @ref
+ * step-8 and @ref
* step_20 "step-20" tutorial
* programs as well as in the
* @ref vector_valued module.
* finite element object, which is composed of several other finite elements
* that are possibly of different type. The result is then a vector-valued
* finite element. Vector valued elements are discussed in a number of
- * tutorial programs, for example @ref step_8 "step-8", @ref step_20
- * "step-20", @ref step_21 "step-21", and in particular in the @ref vector_valued
+ * tutorial programs, for example step-8, @ref step_20
+ * "step-20", step-21, and in particular in the @ref vector_valued
* module.
*
* The overall numbering of degrees of freedom is as follows: for each
* quadrature points.
*
* A use for this function is described
- * in the introduction to the @ref step_18 "step-18"
+ * in the introduction to the step-18
* example program.
*
* The opposite of this function,
* called with a pair of iterators denoting a range on which they
* shall act, by choosing a filtered iterator instead of usual ones.
*
- * This class is used in @ref step_18 "step-18" and
- * @ref step_32 "step-32".
+ * This class is used in step-18 and
+ * step-32.
*
*
* <h3>Predicates</h3>
*
* @note For an example of the
* use of this function see the
- * @ref step_28 "step-28"
+ * step-28
* tutorial program.
*/
template <int dim>
/**
* Return the diamater of the smallest
* active cell of a triangulation. See
- * @ref step_24 "step-24" for an example
+ * step-24 for an example
* of use of this function.
*/
template <int dim, int spacedim>
* surprising fact may lead to programs that do the wrong thing
* since they are not written for the additional cases that can
* happen with anisotropic meshes, see the discussion in the
- * introduction to @ref step_30 "step-30".
+ * introduction to step-30.
*
* <li> @p none:
* Select no smoothing at all.
* This material_id may be set upon construction of a
* triangulation (through the CellData data structure), or later
* through use of cell iterators. For a typical use of this
- * functionality, see the @ref step_28 "step-28" tutorial
+ * functionality, see the step-28 tutorial
* program. The functions of the GridGenerator namespace typically
* set the material ID of all cells to zero. When reading a
* triangulation, the material id must be specified in the input
* <tt>cell->neighbor(neighbor)->face(face_no)->child(subface_no)==cell->face(neighbor)</tt>.
* This case in 3d and how it can happen
* is discussed in the introduction of the
- * @ref step_30 "step-30" tutorial program.
+ * step-30 tutorial program.
*
* This function is impossible
* for <tt>dim==1</tt>.
{
/**
- * An hp equivalent of the ::FEValues class. See the @ref step_27 "step-27"
+ * An hp equivalent of the ::FEValues class. See the step-27
* tutorial program for examples of use.
*
* The idea of this class is as follows: when one assembled matrices in the hp
* DataOut_DoFData::add_data_vector
* functions, this can be achieved.
*
- * See the @ref step_22 "step-22" tutorial
+ * See the step-22 tutorial
* program for an example on how this
* information can be used in practice.
*
* it allows to attach DataPostprocessor arguments to compute derived
* information from a data vector at each quadrature point (for example, the
* Mach number in hypersonic flow can be computed from density and velocities;
- * @ref step_29 "step-29" also shows an example); another piece of information
+ * step-29 also shows an example); another piece of information
* specified through arguments with default values is how certain output
* components should be interpreted, i.e. whether each component of the data
* is logically an independent scalar field, or whether some of them together
* By default, this class produces patches for all active cells. Sometimes,
* this is not what you want, maybe because they are simply too many (and too
* small to be seen individually) or because you only want to see a certain
- * region of the domain (for example in parallel programs such as the @ref step_18 "step-18"
+ * region of the domain (for example in parallel programs such as the step-18
* example program), or for some other reason.
*
* For this, internally build_patches() does not generate
* and possibly the first and second derivates of the solution. Examples are
* the calculation Mach numbers from velocity and density in supersonic flow
* computations, or the computation of the magnitude of a complex-valued
- * solution as demonstrated in @ref step_29 "step-29". Other uses are shown in
- * @ref step_33 "step-33". This class offers the interface to perform such
+ * solution as demonstrated in step-29. Other uses are shown in
+ * step-33. This class offers the interface to perform such
* postprocessing. Given the values and derivatives of the solution at those
* points where we want to generated output, the functions of this class can
* be overloaded to compute new quantities.
* an MPI system.
*
* This function is used in
- * @ref step_17 "step-17" and
- * @ref step_18 "step-18".
+ * step-17 and
+ * step-18.
*/
static void
apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
*
* The integration_loop() function loop takes most of the information
* that it needs to pass to loop() from an IntegrationInfoBox
- * object. Its use is explained in @ref step_38 "step-38", but in
+ * object. Its use is explained in step-38, but in
* short it requires functions that do the local integration on a
* cell, interior or boundary face, and it needs an object (called
* "assembler") that copies these local contributions into the global
* Transfers a discrete FE function (like a solution vector) by interpolation
* while refining and/or coarsening a grid. During interpolation the
* vector is reinitialized to the new size and filled with the interpolated
- * values. This class is used in the @ref step_15 "step-15",
- * @ref step_31 "step-31", and @ref step_33 "step-33" tutorial programs.
+ * values. This class is used in the step-15,
+ * step-31, and step-33 tutorial programs.
*
* <h3>Usage</h3>
*
* soltrans.refine_interpolate(solutions_old, solutions);
* @endverbatim
* This is used in several of the tutorial programs, for example
- * @ref step_31 "step-31".
+ * step-31.
*
* <li> If the grid will be refined AND coarsened
* then use @p SolutionTransfer as follows
* \phi_i(x) \phi_j(x) dx$ and $f_i = \int_\Omega f(x) \phi_i(x)
* dx$. The solution vector $v$ then is the nodal representation of
* the projection <i>f<sub>h</sub></i>. The project() functions are
- * used in the @ref step_21 "step-21" and @ref step_23 "step-23"
+ * used in the step-21 and step-23
* tutorial programs.
*
* In order to get proper results, it be may necessary to treat
*
* The use of this function is
* explained in more detail in
- * @ref step_31 "step-31". It
+ * step-31. It
* doesn't make much sense in 1d,
* so the function throws an
* exception in that case.
# however, undo that change if the dollar sign was escaped with a backslash
s/\\\@f\$/\$/g;
-# we don't let doxygen put everything into a namespace
+# We don't let doxygen put everything into a namespace
# dealii. consequently, doxygen can't link references that contain an
# explicit dealii:: qualification. remove it and replace it by the
# global scope
# that's not our fault. This is reported here:
# http://bugzilla.gnome.org/show_bug.cgi?id=365053
s/dealii::/::/g;
+
+# Replace all occurrences of something like step-xx by
+# @ref step_xx "step-xx"
+# so that doxygen properly cross references them. Before we had
+# this rule, we actually had to write this sequence out in our
+# documentation. Unfortunately, as a consequence, there are vestiges
+# of this style, so we can't substitute things that look like
+# "step-xx".
+s/(?<!\")step-(\d+)(?!\")/\@ref step_\1 \"step-\1\"/gi;
*
* Rather, one has to check them every time during execution of a
* program. Within deal.II, this is done using the AssertThrow macro
- * introduced in @ref step_9 "step-9", @ref step_13 "step-13", and
+ * introduced in step-9, step-13, and
* following tutorial programs. The macro checks a condition, and if
* violated throws an exception of one of the types declared in this
* module, using the C++ <code>throw</code> mechanism. Since these
* @code
* Assert (false, ExcInternalError());
* @endcode
- * See the @ref step_7 "step-7" and several other of the tutorial programs for
+ * See the step-7 and several other of the tutorial programs for
* a use of this construct.
*
* As mentioned above, the program is terminated once a call to
Currently, the geodynamics testsuite contains the following
programs:
-- @ref step_8 "step-8": Elasticity
-- @ref step_16 "step-16": A %parallel elasticity solver
-- @ref step_20 "step-20": Porous media flow
-- @ref step_21 "step-21": Multiphase flow through porous media
-- @ref step_22 "step-22": Stokes flow
-- @ref step_31 "step-31": Thermal convection (Boussinesq flow)
-- @ref step_32 "step-32": A %parallel Boussinesq solver for mantle convection
+- step-8: Elasticity
+- step-16: A %parallel elasticity solver
+- step-20: Porous media flow
+- step-21: Multiphase flow through porous media
+- step-22: Stokes flow
+- step-31: Thermal convection (Boussinesq flow)
+- step-32: A %parallel Boussinesq solver for mantle convection
Some of these programs were developed under contract from the California
Institute of Technology with support by the National Science Foundation
run on larger scale machines. For example, a small mantle convection code
built on deal.II may be used to determine whether second order elements are
useful for this purpose (see, for example, the results shown in
- @ref step_31 "step-31"). If so, then one may use this result to implement
+ step-31). If so, then one may use this result to implement
second, rather than first, order elements in dedicated, large-scale mantle
convection codes such as that which
Ghattas and Zhong are building and that may run on 10,000s of processors, a
* by creating one block for each vector component. However, this is not
* always so, and the DoFRenumbering::component_wise function allows to group
* several vector components into the same block (see, for example, the @ref
- * step_22 "step-22" or @ref step_31 "step-31" tutorial programs, as opposed
- * to @ref step_20 "step-20"). </dd>
+ * step_22 "step-22" or step-31 tutorial programs, as opposed
+ * to step-20). </dd>
*
* <dt class="glossary">@anchor GlossComponent <b>Component</b></dt>
*
*
* The function Triangulation::create_triangulation, which is called
* by the various functions in GridGenerator and GridIn (but can also
- * be called from user code, see @ref step_14 "step-14" will signal
+ * be called from user code, see step-14 will signal
* the creation of coarse meshes with distorted cells by throwing an
* exception of type Triangulation::DistortedCellList. There are
* legitimate cases for creating meshes with distorted cells (in
* from http://www.math.tamu.edu/~bangerth/publications.html .
*
* The numerical examples shown in that paper are generated with a slightly
- * modified version of @ref step_27 "step-27". The main difference to that
+ * modified version of step-27. The main difference to that
* tutorial program is that various operations in the program were timed for
* the paper to compare different options and show that $hp$ methods are
* really not all that expensive.
* "Adaptive multilevel methods with local smoothing", that
* describes many of the algorithms and data structures used in the implementation
* of the multigrid framework of deal.II. It underlies the implementation of
- * the classes that are used in @ref step_16 "step-16" for multigrid
+ * the classes that are used in step-16 for multigrid
* methods.
*
* The full reference for this paper is as follows:
* element is primitive. This includes, in particular, all scalar
* elements as well as vector-valued elements assembled via the
* FESystem class from other primitive (for example scalar) elements
- * as shown in @ref step_8 "step-8", @ref step_29 "step_29", @ref
+ * as shown in step-8, @ref step_29 "step_29", @ref
* step_22 "step-22" and several others. On the other hand,
* the FE_RaviartThomas class used
- * in @ref step_20 "step-20" and @ref step_21 "step-21", or the
+ * in step-20 and step-21, or the
* FE_Nedelec class provide non-primitive finite elements because
* there, each vector-value shape function may have several non-zero
* components.</dd>
* cell is associated with.
*
* For programs that are parallelized based on MPI but where each processor
- * stores the entire triangulation (as in, for example, @ref step_18 "step-18"
- * or @ref step_32 "step-32", subdomain ids are assigned to cells by
+ * stores the entire triangulation (as in, for example, step-18
+ * or step-32, subdomain ids are assigned to cells by
* partitioning a mesh, and each MPI process then only works on those cells it
* "owns", i.e. that belong to a subdomain that it is associated with
* (traditionally, this is the case for the subdomain id whose numerical value
* @defgroup hp hp finite element support
*
* Classes and functions that have to do with hp finite elements. The
- * @ref step_27 "step-27" tutorial program gives an overview of how to
+ * step-27 tutorial program gives an overview of how to
* use the classes in this namespace.
*
* The hp namespace implements the algorithms and data structures used for
* a common strategy is to first build and intermediate sparsity pattern that
* uses a less efficient storage scheme during construction of the sparsity
* pattern and later copy it directly into the static, compressed form (see,
- * for example the @ref step_11 "step-11", @ref step_18 "step-18", and @ref
+ * for example the step-11, step-18, and @ref
* step_27 "step-27" tutorial programs).
*
*
* up to twice the memory in the worst case) for additional speed which is
* noticeable in cases with many nonzero entries. CompressedSetSparsityPattern
* on the other hand uses a lot more memory but may perform better in cases with
- * many nonzero entries per row. See for example the @ref step_27 "step-27"
- * and @ref step_22 "step-22" tutorial programs.
+ * many nonzero entries per row. See for example the step-27
+ * and step-22 tutorial programs.
*
* As a guideline you should start using CompressedSparsityPattern and try the
* other variants if you run into problems. Switching between them should be as
*
* These classes inherit the same tradeoffs regarding their efficiency from
* their non-block classes (see above). See their documentation and
- * @ref step_22 "step-22" for more information.
+ * step-22 for more information.
*
* @ingroup Matrices
*/
* Finally, we have several auxiliary classes, namely MGLevelObject,
* which stores an object on each level
*
- * See the @ref step_16 "step-16" example program on how to use this
+ * See the step-16 example program on how to use this
* functionality.
*/
* Interfaces to Trilinos exist in the TrilinosWrappers namespace,
* making matrices, vectors, and solvers look like the corresponding
* deal.II classes. Their use is explained in the @ref step_31
- * "step-31", @ref step_32 "step-32", and @ref step_33 "step-33"
+ * "step-31", step-32, and step-33
* tutorial programs. The <a
* href="../../readme-petsc-trilinos.html">PETSc and Trilinos
* ReadMe</a> file explains how to configure deal.II to use this
* vectors.
*
* Some discussion on this topic can be found in the discussion of the
- * InverseMatrix class in @ref step_20 "step-20".
+ * InverseMatrix class in step-20.
*
* @ingroup LAC
*/
* solution variable is not a scalar function, but a vector-valued function or
* a set of functions. This includes, for example:
* <ul>
- * <li>The elasticity equation discussed in @ref step_8 "step-8",
- * @ref step_17 "step-17", and @ref step_18 "step-18" in which the
+ * <li>The elasticity equation discussed in step-8,
+ * step-17, and step-18 in which the
* solution is the vector-valued displacement at each point.
* <li>The mixed Laplace equation and its extensions discussed in
- * @ref step_20 "step-20", and @ref step_21 "step-21" in which the
+ * step-20, and step-21 in which the
* solution is the scalar pressure and the vector-valued velocity
* at each point.
* <li>The Stokes equation and its extensions discussed in
- * @ref step_22 "step-22", and @ref step_31 "step-31" in which again
+ * step-22, and step-31 in which again
* the solution is the scalar pressure and the vector-valued velocity
* at each point.
* <li>Complex-valued solutions consisting of real and imaginary parts, as
- * discussed for example in @ref step_29 "step-29".
+ * discussed for example in step-29.
* </ul>
*
* This page gives an overview of how to implement such vector-valued problems
* first with a weak (variational) formulation of the problem that takes into
* account all the solution variables. To understand how to do that, let us
* consider a simple example, the mixed Laplace equations discussed in
- * @ref step_20 "step-20":
+ * step-20:
@f{eqnarray*}
\textbf{u} + \nabla p &=& 0,
\\
* component.
*
* As it turns out these (equivalent) choices do not lead to a stable scheme
- * for the mixed Laplace equation. In @ref step_20 "step-20", we therefore use
+ * for the mixed Laplace equation. In step-20, we therefore use
* a Raviart-Thomas element for the velocities. What exactly this means may be
* of less concern for now except that the FE_RaviartThomas class describes
* elements that already have <code>dim</code> components. For the pressure,
*
* The next step is to assemble the linear system. How to do this for the
* simple case of a scalar problem has been shown in many tutorial programs,
- * starting with @ref step_3 "step-3". Here we will show how to do it for
+ * starting with step-3. Here we will show how to do it for
* vector problems.
*
* How to do this is possibly best explained by showing an example
* divergence $\sum_{d=0}^{dim-1} \frac{\partial\phi_d}{\partial x_d}$.
* </ul>
* Other examples of using extractors and views are shown in tutorial programs
- * @ref step_21 "step-21",
- * @ref step_22 "step-22",
- * @ref step_31 "step-31" and a few other programs.
+ * step-21,
+ * step-22,
+ * step-31 and a few other programs.
*
*
* @anchor VVAlternative
* There are situations in which one can optimize the assembly of a matrix or
* right hand side vector a bit, using knowledge of the finite element in
* use. Consider, for example, the bilinear form of the elasticity equations
- * which we are concerned with first in @ref step_8 "step-8":
+ * which we are concerned with first in step-8:
*
@f[
a({\mathbf u}, {\mathbf v}) =
}
* @endcode
*
- * Now, this is not the code used in @ref step_8 "step-8". In fact,
+ * Now, this is not the code used in step-8. In fact,
* if one used the above code over the one implemented in that program,
* it would run about 8 per cent slower. It can be improved (bringing
* down the penalty to about 4 per cent) by taking a close look at the
}
* @endcode
*
- * So if, again, this is not the code we use in @ref step_8 "step-8", what do
+ * So if, again, this is not the code we use in step-8, what do
* we do there? The answer rests on the finite element we use. There, we use the
* following element:
* @code
* component is going to be nonzero. This information is provided by the
* FiniteElement::system_to_component_index function. What can be done with
* it, using the example above, is explained in detail in
- * @ref step_8 "step-8".
+ * step-8.
*
*
* @anchor VVBlockSolvers
* vector-valued problem. However, then it also has to be solved. This is more
* complicated. Naively, one could just consider the matrix as a whole. For
* most problems, this matrix is not going to be definite (except for special
- * cases like the elasticity equations covered in @ref step_8 "step-8" and
- * @ref step_17 "step-17"). It will, often, also not be symmetric. This rather
+ * cases like the elasticity equations covered in step-8 and
+ * step-17). It will, often, also not be symmetric. This rather
* general class of matrices presents problems for iterative solvers: the lack
* of structural properties prevents the use of most efficient methods and
* preconditioners. While it can be done, the solution process will therefore
* What is needed is to re-enumerate degrees of freedom so that velocities
* come first and pressures last. This can be done using the
* DoFRenumbering::component_wise function, as explained in @ref step_20
- * "step-20", @ref step_21 "step-21", @ref step_22 "step-22", and @ref step_31
+ * "step-20", step-21, step-22, and @ref step_31
* "step-31". After this, at least the degrees of freedom are partitioned
* properly.
*
* But then we still have to make use of it, i.e. we have to come up with a
- * solver that uses the structure. For example, in @ref step_20 "step-20", we
+ * solver that uses the structure. For example, in step-20, we
* do a block elimination of the linear system
@f{eqnarray*}
\left(
* large whole matrix we had before.
*
* How a solver like this is implemented is explained in more detail in @ref
- * step_20 "step-20", @ref step_31 "step-31", and a few other tutorial
+ * step_20 "step-20", step-31, and a few other tutorial
* programs. What we would like to point out here is that we now need a way to
* extract certain parts of a matrix or vector: if we are to multiply, say,
* the $U$ part of the solution vector by the $M$ part of the global matrix,
*
* To show how to do this, let us consider the second equation $MU=F-BP$ to be
* solved above. This can be achieved using the following sequence similar to
- * what we have in @ref step_20 "step-20":
+ * what we have in step-20:
* @code
Vector<double> tmp (solution.block(0).size());
system_matrix.block(0,1).vmult (tmp, solution.block(1));
* allocation for the outer vector as well as for all the inner vectors.
* Secondly, maybe we are only interested in the velocities,
* for example to solve an advection problem in a second stage (as, for
- * example, in @ref step_21 "step-21" or @ref step_31 "step-31"). In that
+ * example, in step-21 or step-31). In that
* case, one would have to hand-extract these values like so:
* @code
for (unsigned int q=0; q<n_q_points; ++q)
# $Id$
# Version: $Name$
#
-# Copyright (C) 2006, 2007 by the deal.II authors
+# Copyright (C) 2006, 2007, 2010 by the deal.II authors
#
# This file is subject to QPL and may not be distributed
# without copyright and license information. Please refer
s/\t/ /g;
if (($state == $program_mode) && m!^\s*//!)
- {
+ {
$state = $comment_mode;
print " * \@endcode\n";
print " * \n";
# if in comment mode and no comment line: toggle state.
# don't do so, if only a blank line
elsif (($state == $comment_mode) && !m!^\s*//! && !m!^\s*$!)
- {
+ {
$state = $program_mode;
print " * \n";
print " * \@code\n";
}
-
- if ($state == $comment_mode)
+
+ if ($state == $comment_mode)
{
# in comment mode: first skip leading whitespace and
# comment // signs
$_ = "\n * <a name=\"$sect_name\"></a> \n * $_";
}
- # then replace references to other example programs automatically:
- s!step-(\d+)!\@ref step_\1 \"step-\1\"!g;
-
# finally print this line
print " * $_\n";
{
# in program mode, output the program line. the only thing we need
# to do is to avoid $ signs because that confuses doxygen. since
- # we don't want formulas rendered in the program text anyway,
+ # we don't want formulas rendered in the program text anyway,
# simply replace them by spaces (it would be nice to suppress their
# meaning somehow, but I don't know how...)
s/\$//g;
# $Id$
# Version: $Name$
#
-# Copyright (C) 2006, 2007 by the deal.II authors
+# Copyright (C) 2006, 2007, 2010 by the deal.II authors
#
# This file is subject to QPL and may not be distributed
# without copyright and license information. Please refer
s/\t/ /g;
if (($state == $program_mode) && m!^\s*//!)
- {
+ {
$state = $comment_mode;
print "</code></pre>\n";
print "\n";
# if in comment mode and no comment line: toggle state.
# don't do so, if only a blank line
elsif (($state == $comment_mode) && !m!^\s*//! && !m!^\s*$!)
- {
+ {
$state = $program_mode;
print "</p>\n";
print "\n";
print "<pre><code>\n";
}
-
- if ($state == $comment_mode)
+
+ if ($state == $comment_mode)
{
# in comment mode: first skip leading whitespace and
# comment // signs
# replace double dashes in comments by —
s!--!—!g;
- # then replace references to other example programs automatically:
- s!(step-\d+)!<a href="\1.html" target="body">\1</a>!g;
-
# finally print this line
print $_, "\n";
<h3>General</h3>
<ol>
-<li> <p>New: @ref step_33 "step-33", written mainly by David Neckels,
+<li> <p>New: step-33, written mainly by David Neckels,
illustrates how the Euler equations of compressible inviscid gas
dynamics can be solved. More generally, the program also provides
a framework for the solution of arbitrary hyperbolic, nonlinear
(WB, 2008/05/23)
</p>
-<li> <p>Improved: @ref step_16 "step-16" has been updated to show the
+<li> <p>Improved: step-16 has been updated to show the
use of Multigrid with Dirichlet boundary conditions.
<br>
(GK, 2007/05/13)
</p>
-<li> <p>New: @ref step_22 "step-22", written by Martin Kronbichler and Wolfgang
+<li> <p>New: step-22, written by Martin Kronbichler and Wolfgang
Bangerth, shows how to solve, and in particular how to solve
efficiently, the Stokes equations that describe viscous
incompressible fluid flow.
(WB, 2008/03/28)
</p>
-<li> <p>New: @ref step_29 "step-29" demonstrates how problems involving
+<li> <p>New: step-29 demonstrates how problems involving
complex %numbers can be implemented by viewing real and imaginary parts
of a complex-valued solution as the two components of a vector-valued
function.
<li> <p>New: A significantly simpler way to code the assembly of linear
systems for problems with more than one solution variable has been
implemented. This is explained in detail in the report on @ref vector_valued
- and tutorial programs @ref step_20 "step-20" and @ref step_21 "step-21"
+ and tutorial programs step-20 and step-21
have been converted as well.
<br>
(WB, 2008/02/23)
only the VTK reader in DataOutBase::write_vtk as well as the deal.II
intermediate format writer make use of it.
<p>
- The use of this ability is explained in the @ref step_22 "step-22" and
- @ref step_33 "step-33" tutorial programs.
+ The use of this ability is explained in the step-22 and
+ step-33 tutorial programs.
<br>
(WB 2007/10/11)
</p> </li>
that may be used to indicate that certain components of the data
logically form a vector field. The output formats for which this
is presently implemented then indicate this fact in the output file. The
- mechanism is shown in action in @ref step_22 "step-22" and
- @ref step_33 "step-33".
+ mechanism is shown in action in step-22 and
+ step-33.
<br>
(WB 2007/10/14)
</p></li>
<ol>
<li>
<p>
- Updated: The @ref step_22 "step-22" tutorial program now uses some
+ Updated: The step-22 tutorial program now uses some
advanced features of the ConstraintMatrix such as imhomogeneities
for implementing Dirichlet boundary condition and condensation of
constraints on the fly.
<li>
<p>
- New: A new tutorial program @ref step_34 "step-34" was added to the
+ New: A new tutorial program step-34 was added to the
library that shows the usage of the new codimension one functionality
recently added to the library. In this tutorial we show the use of
bondary element methods on piecewise constant functions defined over
<li>
<p>
- Updated: The @ref step_33 "step-33" tutorial program now uses the
+ Updated: The step-33 tutorial program now uses the
Trilinos wrapper
classes for accessing Trilinos linear algebra routines instead of
implementing that by hand.
if you intend two-dimensional triangulations of two-dimensional
domains.
<br>
- The new functionality is described in the @ref step_34 "step-34"
+ The new functionality is described in the step-34
tutorial program.
<br>
(Cataldo Manigrasso, Luca Heltai 2008/12/05)
<li>
<p>
- New: A new tutorial program, @ref step_31 "step-31", is contributed by
+ New: A new tutorial program, step-31, is contributed by
Martin Kronbichler and Wolfgang Bangerth and shows the
solution of the Boussinesq approximation to thermally driven
convection in an incompressible fluid.
<li>
<p>
- Updated: In the @ref step_22 "step-22" tutorial program the generation of
+ Updated: In the step-22 tutorial program the generation of
the sparsity pattern using the class BlockCompressedSetSparsityPattern has
been replaced by a pattern of the class BlockCompressedSimpleSparsityPattern,
which uses far less memory and is slightly faster.
<li>
<p>
- Fixed: In the @ref step_23 "step-23" tutorial program the terms corresponding
+ Fixed: In the step-23 tutorial program the terms corresponding
to the external force were computed incorrectly. The error wasn't visible
right away since in the program the right hand side is assumed to be zero.
This has now been fixed.
<li>
<p>
- New: A new tutorial program, @ref step_28 "step-28", is contributed by
+ New: A new tutorial program, step-28, is contributed by
Yaqi Wang. It illustrates the solution
of a coupled system of diffusion equations relevant to
nuclear reactor physics where we use different meshes for
<p>
New: Contributed by Tobias Leicht as well as Ralf Hartmann, deal.II
now supports anisotropic refinement of meshes. This functionality
- is explained in the @ref step_30 "step-30" tutorial program.
+ is explained in the step-30 tutorial program.
<br>
(Tobias Leicht, RH 2008/07/08)
</p>
<li>
<p>
- Fixed: In the new @ref step_33 "step-33" tutorial program there was
+ Fixed: In the new step-33 tutorial program there was
a place where we incorrectly passed the diameter of a cell, rather
than a face, to a function. This is now fixed.
<br>
<li>
<p>
- Fixed: The new @ref step_33 "step-33" tutorial program had a place where
+ Fixed: The new step-33 tutorial program had a place where
we didn't release some memory. This is now fixed.
<br>
(Chih-Che Chueh, WB 2008/06/02)
New: FEValues objects can now directly calculate the Laplacians of
solution functions at quadrature points. This is equivalent to the
trace of the Hessians, but faster and more convenient in user code.
- @ref step_31 "step-31" program has also been updated to use this.
+ step-31 program has also been updated to use this.
<br>
(Martin Kronbichler 2008/12/15)
</p>
<code>fe_values[velocities].get_function_values (global_solution,
local_velocity_values)</code>. See the new section at the end of the
@ref vector_valued module for an overview of this facility. The
- @ref step_31 "step-31" program has also been updated to use this.
+ step-31 program has also been updated to use this.
<br>
(WB 2008/12/13)
</p>
<li>
<p>
- When running in parallel, the @ref step_18 "step-18" tutorial program
+ When running in parallel, the step-18 tutorial program
produced an error indicating that resetting user pointers was not
possible. This is now fixed.
</p>
support multigrid also for continuous finite elements on adaptively
refined meshes (uniformly refined meshes and discontinuous elements
have worked for a long time). As part of this process,
- @ref step_16 "step-16" has
+ step-16 has
been rewritten and now solves the same problem
- @ref step_6 "step-6" solves, just
+ step-6 solves, just
with a multigrid solver.
<br>
(Bärbel Janssen, WB 2010/02/13)
<li>
<p>
- New: There is now a new tutorial program, @ref step_35 "step-35",
+ New: There is now a new tutorial program, step-35,
contributed by Abner Salgado-Gonzalez, that implements a solver
for the Navier-Stokes equations using a decoupled projection
scheme.
<li>
<p>
- Fixed: @ref step_31 "step-31" had a bug in the computation of the
+ Fixed: step-31 had a bug in the computation of the
global scaling parameter in the function that evaluates the artificial
viscosity: we computed
$c(\mathbf{u},T) =
\ |\mathrm{diam}(\Omega)|^{\alpha-2}$. This didn't matter much in this
program because $\mathrm{diam}(\Omega)=2^{1/\textrm{dim}}$ and so is close
to one. It would matter, however, if the domain had been different, as
- it is, for example, in @ref step_32 "step-32".
+ it is, for example, in step-32.
<br>
(WB 2009/08/19)
</p>
<li>
<p>
- New: There is now a new tutorial program, @ref step_36 "step-36",
+ New: There is now a new tutorial program, step-36,
contributed by Toby D. Young and Wolfgang Bangerth, that demonstrates
solving eigenvalue problems.
<br>
<h3>The software</h3>
The step-14 example program builds heavily on the techniques already used in
-the @ref step_13 "step-13" program. Its implementation of the dual weighted residual error
+the step-13 program. Its implementation of the dual weighted residual error
estimator explained above is done by deriving a second class, properly called
<code>DualSolver</code>, from the <code>Solver</code> base class, and having a class
(<code>WeightedResidual</code>) that joins the two again and controls the solution
This example shows the basic usage of the multilevel functions in
-deal.II. It solves the same problem as used in @ref step_6 "step-6",
+deal.II. It solves the same problem as used in step-6,
but demonstrating the things one has to provide when using multigrid
as a preconditioner. In particular, this requires that we define a
hierarchy of levels, provide transfer operators from one level to the
<h3>The testcase</h3>
The problem we solve here is exactly the same as in
-@ref step_6 "step-6", the only difference being the solver we use
+step-6, the only difference being the solver we use
here. You may want to look there for a definition of what we solve,
right hand side and boundary conditions. Obviously, the program would
also work if we changed the geometry and other pieces of data that
<h1>Results</h1>
The output that this program generates is, of course, the same as that
-of @ref step_6 "step-6", so you may see there for more results. On the
+of step-6, so you may see there for more results. On the
other hand, since no tutorial program is a good one unless it has at
least one colorful picture, here is, again, the solution:
<tr><td>5120</td> <td>94</td> <td>11</td> </tr>
<tr><td>20480</td><td>180</td><td>13</td></tr>
</table>
-This isn't only fewer iterations than in @ref step_5 "step-5" (each of which
+This isn't only fewer iterations than in step-5 (each of which
is, however, much more expensive) but more importantly, the number of
iterations also grows much more slowly under mesh refinement (again,
it would be almost constant if the coefficient was constant rather
<h3> Possible extensions </h3>
A close inspection of this program's performance shows that it is mostly
-dominated by matrix-vector operations. @ref step_37 "step-37" shows one way
+dominated by matrix-vector operations. step-37 shows one way
how this can be avoided by working with matrix-free methods.
Another avenue would be to use algebraic multigrid methods. The
sort of multigrid hierarchy for good performance but can figure out
level matrices and similar things out by itself. Algebraic multigrid
methods do exactly this, and we will use them in
-@ref step_31 "step-31" for the solution of a Stokes problem.
+step-31 for the solution of a Stokes problem.
This program does not introduce any new mathematical ideas; in fact, all it
-does is to do the exact same computations that @ref step_8 "step-8"
+does is to do the exact same computations that step-8
already does, but it does so in a different manner: instead of using deal.II's
own linear algebra classes, we build everything on top of classes deal.II
provides that wrap around the linear algebra implementation of the <a
As one would hope for, the x- (left) and y-displacements (right) shown here
-closely match what we already saw in @ref step_8 "step-8". As shown
-there and in @ref step_22 "step-22", we could as well have produced a
+closely match what we already saw in step-8. As shown
+there and in step-22, we could as well have produced a
vector plot of the displacement field, rather than plotting it as two
separate scalar fields. What may be more interesting,
though, is to look at the mesh and partition at this step:
This tutorial program is another one in the series on the elasticity problem
-that we have already started with @ref step_8 "step-8" and @ref step_17 "step-17". It extends it into two
+that we have already started with step-8 and step-17. It extends it into two
different directions: first, it solves the quasistatic but time dependent
elasticity problem for large deformations with a Lagrangian mesh movement
approach. Secondly, it shows some more techniques for solving such problems
using %parallel processing with PETSc's linear algebra. In addition to this, we
-show how to work around the main bottleneck of @ref step_17 "step-17", namely that we
+show how to work around the main bottleneck of step-17, namely that we
generated graphical output from only one process, and that this scaled very
badly with larger numbers of processes and on large problems. Finally, a good
number of assorted improvements and techniques are demonstrated that have not
been shown yet in previous programs.
-As before in @ref step_17 "step-17", the program runs just as fine on a single sequential
+As before in step-17, the program runs just as fine on a single sequential
machine as long as you have PETSc installed. Information on how to tell
deal.II about a PETSc installation on your system can be found in the deal.II
README file, which is linked to from the <a href="../../index.html">main
<h3>Parallel graphical output</h3>
-In the @ref step_17 "step-17" example program, the main bottleneck for %parallel computations
+In the step-17 example program, the main bottleneck for %parallel computations
was that only the first processor generated output for the entire domain.
Since generating graphical output is expensive, this did not scale well when
large numbers of processors were involved. However, no viable ways around this
been generated, and <code>NNN</code> the number of the MPI process that did this.
The next step is to convert this file or these files into whatever
-format you like. The program that does this is the @ref step_19 "step-19" tutorial program:
+format you like. The program that does this is the step-19 tutorial program:
for example, for the first time step, call it through
@code
../step-19/step-19 solution-0001.0000.*.d2 solution-0001.0000.gmv
@endcode
to merge all the intermediate format files into a single file in GMV
format. More details on the parameters of this program and what it can do for
-you can be found in the documentation of the @ref step_19 "step-19" tutorial program.
+you can be found in the documentation of the step-19 tutorial program.
\varepsilon(\varphi_j(\mathbf{x}_q))],
@f]
with quadrature points $\mathbf{x}_q$ and weights $w_q$. We have built these
- contributions before, in @ref step_8 "step-8" and @ref step_17 "step-17", but in both of these cases we
+ contributions before, in step-8 and step-17, but in both of these cases we
have done so rather clumsily by using knowledge of how the rank-4 tensor $C$
is composed, and considering individual elements of the strain tensors
$\varepsilon(\varphi_i),\varepsilon(\varphi_j)$. This is not really
solving the linear system using the CG method, using an incomplete LU
decomposition as a preconditioner (in the %parallel case, it uses an ILU of
each processor's block separately). It is virtually unchanged
- from @ref step_17 "step-17".
+ from step-17.
<li> <code>update_quadrature_point_history ()</code> [via
<code>solve_timestep ()</code>]: Based on the displacement field $\Delta \mathbf{u}^n$ computed before, we update the stress values in all quadrature points
values of the quadrature point data on other cells if the corresponding cells
on the new mesh are assigned to this process after subdividing the new mesh. A
global communication of these data elements is therefore necessary, making the
-entire process a little more unpleasant. The @ref step_28 "step-28" tutorial
+entire process a little more unpleasant. The step-28 tutorial
program shows how to work with different meshes at the same time, albeit for a
different kind of problem, giving indications on how to approach the problem
for time-dependent adaptivity as well.
Let us convert these files in deal.II intermediate format to gmv
format (this assumes that you have already compiled the
-@ref step_19 "step-19" example program):
+step-19 example program):
@code
examples/step-18> ../step-19/step-19
<h1>Introduction</h1>
-In @ref step_18 "step-18", we saw a need to write
+In step-18, we saw a need to write
output files in an intermediate format: in a parallel program, it doesn't scale
well if all processors participate in computing a result, and then only a
single processor generates the graphical output. Rather, each of them should
Let us try to run this program on a set of input files generated by a modified
-@ref step_18 "step-18" run on 32 nodes of a
+step-18 run on 32 nodes of a
cluster. The computation was rather big, with more
than 350,000 cells and some 1.2M unknowns. That makes for 32 rather big
intermediate files that we will try to merge using the present program. Here is
If visualized, the output looks very much like that shown for
-@ref step_18 "step-18". But that's not quite as
+step-18. But that's not quite as
important for the moment, rather we are interested in showing how to use the
parameter file. To this end, remember that if no parameter file is given, or if
it is empty, all the default values listed above are used. However, whatever we
definition of solvers and preconditioners, and finally the actual test case we
are going to solve.
-We are going to extend this tutorial program in @ref step_21 "step-21" to
+We are going to extend this tutorial program in step-21 to
solve not only the mixed Laplace equation, but add another equation that
describes the transport of a mixture of two fluids.
The necessary question then is: how do we do this in a program?
Vector-valued elements have already been discussed in previous tutorial
-programs, the first time and in detail in @ref step_8 "step-8". The main difference there
+programs, the first time and in detail in step-8. The main difference there
was that the vector-valued space $V_h$ is uniform in all its components: the
$dim$ components of the displacement vector are all equal and from the same
function space. What we could therefore do was to build $V_h$ as the outer
product of the $dim$ times the usual $Q(1)$ finite element space, and by this
make sure that all our shape functions have only a single non-zero vector
component. Instead of dealing with vector-valued shape functions, all we did
-in @ref step_8 "step-8" was therefore to look at the (scalar) only non-zero component and
+in step-8 was therefore to look at the (scalar) only non-zero component and
use the <code>fe.system_to_component_index(i).first</code> call to figure out
which component this actually is.
gradients of shape functions at quadrature points. For example, we would call
<code>fe_values.shape_value(i,q_point)</code> to obtain the value of the
<code>i</code>th shape function on the quadrature point with number
-<code>q_point</code>. Later, in @ref step_8 "step-8" and other tutorial programs, we learned
+<code>q_point</code>. Later, in step-8 and other tutorial programs, we learned
that this function call also works for vector-valued shape functions (of
primitive finite elements), and that it returned the value of the only
non-zero component of shape function <code>i</code> at quadrature point
store the result, as this saves us a few otherwise repeated computations (it is
possible to save even more repeated operations by calculating all relevant
quantities in advance and then only inserting the results in the actual loop,
-see @ref step_22 "step-22" for a realization of that approach).
+see step-22 for a realization of that approach).
The final result then looks like this, working in every space dimension:
@code
There is one final term that we have to take care of: the right hand side
contained the term $(g,{\mathbf v}\cdot {\mathbf n})_{\partial\Omega}$, constituting the
weak enforcement of pressure boundary conditions. We have already seen in
-@ref step_7 "step-7" how to deal with face integrals: essentially exactly the same as with
+step-7 how to deal with face integrals: essentially exactly the same as with
domain integrals, except that we have to use the FEFaceValues class
instead of <code>FEValues</code>. To compute the boundary term we then simply have
to loop over all boundary faces and integrate there. The mechanism works in
In this code, the constructor takes a reference to a block sparse matrix for
the entire system, and a reference to an object representing the inverse of
the mass matrix. It stores these using <code>SmartPointer</code> objects (see
-@ref step_7 "step-7"), and additionally allocates two temporary vectors <code>tmp1</code> and
+step-7), and additionally allocates two temporary vectors <code>tmp1</code> and
<code>tmp2</code> for the vectors labeled $w,y$ in the list above.
In the matrix-vector multiplication function, the product $Sv$ is performed in
@image html step-20.u-random.png
@image html step-20.p-random.png
-We will use these permeability fields again in @ref step_21 "step-21".
+We will use these permeability fields again in step-21.
"step-18").
The equations covered here are an extension of the material already covered in
-@ref step_20 "step-20". In particular, they fall into the class of
+step-20. In particular, they fall into the class of
vector-valued problems. A toplevel overview of this topic can be found in the
@ref vector_valued module.
So far, this looks like an ordinary stationary, Poisson-like equation that we
can solve right away with the techniques of the first few tutorial programs
-(take a look at @ref step_6 "step-6", for example, for something very
+(take a look at step-6, for example, for something very
similar). However, we have not said anything yet about the saturation, which
of course is going to change as the fluids move around.
constraint is that the divergence of the electric displacement field
equals the charge density, $\textrm{div}\ \mathbf D = \rho$ and that the
divergence of the magnetic flux density is zero: $\textrm{div}\ \mathbf
-B = 0$); even the quasistatic model of @ref step_18 "step-18" falls into this
+B = 0$); even the quasistatic model of step-18 falls into this
category. We will see that the different character of the two equations
will inform our discretization strategy for the two equations.
total velocity $\mathbf u$ appearing in the transport equation as a function
of the pressure, but can rather take the primary variable for it. Given the
saddle point structure of the first two equations and their similarity to the
-mixed Laplace formulation we have introduced in @ref step_20 "step-20", it
+mixed Laplace formulation we have introduced in step-20, it
will come as no surprise that we will use a mixed discretization again.
But let's postpone this for a moment. The first business we have with these
variables on the present cell, whereas $S_-,\mathbf{u}_-$ (needed on the
inflow part of the boundary of $K$) are quantities taken from the neighboring
cell. Some more context on discontinuous element techniques and evaluation of
-fluxes can also be found in @ref step_12 "step-12".
+fluxes can also be found in step-12.
<h3>Linear solvers</h3>
The linear solvers used in this program are a straightforward extension of the
-ones used in @ref step_20 "step-20". Essentially, we simply have to extend
+ones used in step-20. Essentially, we simply have to extend
everything from
two to three solution components. If we use the discrete spaces
mentioned above and put shape functions into the bilinear forms, we
The solution scheme then involves the following steps:
<ol>
<li>Solve for the pressure $p^{n+1}$ using the Schur complement
- technique introduced in @ref step_20 "step-20".
+ technique introduced in step-20.
<li>Solve for the velocity $\mathbf u^{n+1}$ as also discussed in
- @ref step_20 "step-20".
+ step-20.
<li>Compute the term $F_3-\triangle t\; H \mathbf u^{n+1}$, using
the just computed velocities.
anything useful once we have negative saturations or ones larger than one.
Note that we will have similar restrictions on the time step also in @ref
-step_23 "step-23" and @ref step_24 "step-24" where we solve the time dependent
+step_23 "step-23" and step-24 where we solve the time dependent
wave equation, another hyperbolic problem. We will also come back to the issue
of time step choice below in the section on <a href="#extensions">possible
extensions to this program</a>.
"step-20":
<ul>
<li>A function that models a single, winding crack that snakes through the
- domain. In analgy to @ref step_20 "step-20", but taking care of the slightly
+ domain. In analgy to step-20, but taking care of the slightly
different geometry we have here, we describe this by the following function:
@f[
k(\mathbf x)
replacing the continuous variables and test functions by finite element
functions in finite dimensional spaces $\textbf V_{g,h}\subset \textbf V_g,
Q_h\subset Q$, we have to make sure that $\textbf V_h,Q_h$ also satisfy the LBB
-conditions. This is similar to what we had to do in @ref step_20 "step-20".
+conditions. This is similar to what we had to do in step-20.
For the Stokes equations, there are a number of possible choices to ensure
that the finite element spaces are compatible with the LBB condition. A simple
@f}
for all test functions $\textbf v_h, q_h$. Assembling the linear system
associated with this problem follows the same lines used in @ref step_20
-"step-20", @ref step_21 "step-21", and explained in detail in the @ref
+"step-20", step-21, and explained in detail in the @ref
vector_valued module.
F \\ G
\end{array}\right),
@f}
-Like in @ref step_20 "step-20" and @ref step_21 "step-21", we will solve this
+Like in step-20 and step-21, we will solve this
system of equations by forming the Schur complement, i.e. we will first find
the solution $P$ of
@f{eqnarray*}
<li>
While the outer preconditioner has become simpler compared to the
-mixed Laplace case discussed in @ref step_20 "step-20", the issue of
+mixed Laplace case discussed in step-20, the issue of
the inner solver has become more complicated. In the mixed Laplace
discretization, the Schur complement has the form $B^TM^{-1}B$. Thus,
every time we multiplied with the Schur complement, we had to solve a
system_matrix, system_rhs);
@endcode
-This technique is further discussed in the @ref step_27 "step-27" tutorial
+This technique is further discussed in the step-27 tutorial
program. All we need to know here is that this functions does three things
at once: it writes the local data into the global matrix and right hand
side, it distributes the hanging node constraints and additionally
unsatisfactory, if one compares with an elliptic problem of similar
size. Either way, we will introduce below a number of improvements to the
linear solver, a discussion that we will re-consider again with additional
-options in the @ref step_31 "step-31" program.
+options in the step-31 program.
<a name="improved-ilu">
<h5>Better ILU decomposition by smart reordering</h5>
Another idea to improve the situation even more would be to choose a
preconditioner that makes CG for the (0,0) matrix $A$ converge in a
mesh-independent number of iterations, say 10 to 30. We have seen such a
-canditate in @ref step_16 "step-16": multigrid.
+canditate in step-16: multigrid.
<h5>Block Schur complement preconditioner</h5>
Even with a good preconditioner for $A$, we still
An ultimate linear solver for this problem could be imagined as a
combination of an optimal
preconditioner for $A$ (e.g. multigrid) and the block preconditioner
-described above, which is the approach taken in the @ref step_31 "step-31"
+described above, which is the approach taken in the step-31
tutorial program.
<h5>No block matrices and vectors</h5>
This is the first of a number of tutorial programs that will finally
cover "real" time-dependent problems, not the slightly odd form of time
-dependence found in @ref step_18 "step-18" or the DAE model of @ref step_21 "step-21". In particular, this program introduces
-the wave equation in a bounded domain. Later, @ref step_24 "step-24"
+dependence found in step-18 or the DAE model of step-21. In particular, this program introduces
+the wave equation in a bounded domain. Later, step-24
will consider an example of absorbing boundary conditions, and @ref
step_25 "step-25" a kind of nonlinear wave equation producing
solutions called solitons.
used different shape functions, then we would have to compute
integrals that contain shape functions defined on two meshes. This is a
somewhat messy process that we omit here, but that is treated in some
-detail in @ref step_28 "step-28".
+detail in step-28.
Under these conditions (i.e. a mesh that doesn't change), one can optimize the
solution procedure a bit by basically eliminating the solution of the second
<li>More complicated domains or more refined meshes. Remember that the time
step needs to be bounded by the mesh width, so changing the mesh should
always involve also changing the time step. We will come back to this issue
- in @ref step_24 "step-24".
+ in step-24.
<li>Variable coefficients: In real media, the wave speed is often
variable. In particular, the "real" wave equation in realistic media would
<h3>Weak form and discretization</h3>
-As in @ref step_23 "step-23", one first introduces a second variable, which is
+As in step-23, one first introduces a second variable, which is
defined as the derivative of the pressure potential:
@f[
v = \frac{\partial\bar{p}}{\partial t}
to the thermoacoustic problem $f=0$.
The semi-discretized, weak version of this model, using the general $\theta$ scheme
-introduced in @ref step_23 "step-23" is then:
+introduced in step-23 is then:
@f{eqnarray*}
\left(\frac{\bar{p}^n-\bar{p}^{n-1}}{k},\phi\right)_\Omega-
\left(\theta v^{n}+(1-\theta)v^{n-1},\phi\right)_\Omega & = & 0 \\
(-c_0^2k \theta A-c_0 B)\bar{p}^n-Mv^{n} & = &
(c_0^2k(1-\theta)A-c_0B)\bar{p}^{n-1}-Mv^{n-1}+c_0^2k(\theta F^{n}+(1-\theta)F^{n-1}).
@f}
-The matrices $M$ and $A$ are here as in @ref step_23 "step-23", and the
+The matrices $M$ and $A$ are here as in step-23, and the
boundary mass matrix
@f[
B_{ij} = \left(\varphi_i,\varphi_j\right)_{\partial\Omega}
<h3>What the program does</h3>
-Compared to @ref step_23 "step-23", this programs adds the treatment of a
+Compared to step-23, this programs adds the treatment of a
simple absorbing boundary conditions. In addition, it deals with data obtained
from actual experimental measurements. To this end, we need to evaluate the
solution at points at which the experiment also evaluates a real pressure
parallel to mesh lines. This anisotropy leads to wave fronts that aren't
perfectly circular (and would produce sinusoidal signals in the stacked
plots), but are bulged out in certain directions. To make things worse, the
-circular mesh we use (see for example @ref step_6 "step-6" for a view of the
+circular mesh we use (see for example step-6 for a view of the
coarse mesh) is not isotropic either. The net result is that the signal fronts
are not sinusoidal unless the mesh is sufficiently fine. The right image is a
lot better in this respect, though artificts in the form of trailing spurious
u(x,t_0) &=& u_0(x).
\f}
It is a nonlinear equation similar to the wave equation we
-discussed in @ref step_23 "step-23" and @ref step_24 "step-24".
+discussed in step-23 and step-24.
We have chosen to enforce zero Neumann boundary conditions in order for waves
to reflect off the boundaries of our domain. It should be noted, however, that
Dirichlet boundary conditions are not appropriate for this problem. Even
derivative directly and accurately, then we would need a large
stencil (i.e., several time steps would need to be kept in the
memory), which could become expensive. Therefore, in complete analogy
-to what we did in @ref step_23 "step-23" and @ref step_24 "step-24",
+to what we did in step-23 and step-24,
we split the
second-order (in time) sine-Gordon equation into a system of two
first-order (in time) equations, which we call the split, or velocity,
that utilizes the symmetry, such as, for example, SymmLQ, which is also
implemented in deal.II.
-This program uses a clever optimization over @ref step_23 "step-23" and @ref
+This program uses a clever optimization over step-23 and @ref
step_24 "step-24": If you read the above formulas closely, it becomes clear
that the velocity $V$ only ever appears in products with the mass matrix. In
-@ref step_23 "step-23" and @ref step_24 "step-24", we were, therefore, a bit
+step-23 and step-24, we were, therefore, a bit
wasteful: in each time step, we would solve a linear system with the mass
matrix, only to multiply the solution of that system by $M$ again in the next
time step. This can, of course, be avoided, and we do so in this program.
This tutorial program attempts to show how to use $hp$ finite element methods
with deal.II. It solves the Laplace equation and so builds only on the first
-few tutorial programs, in particular on @ref step_4 "step-4" for dimension
-independent programming and @ref step_6 "step-6" for adaptive mesh
+few tutorial programs, in particular on step-4 for dimension
+independent programming and step-6 for adaptive mesh
refinement.
The $hp$ finite element method was proposed in the early 1980s by
The next task we have to consider is what to do with the list of
finite element objects we want to use. In previous tutorial programs,
-starting with @ref step_2 "step-2", we have seen that the DoFHandler
+starting with step-2, we have seen that the DoFHandler
class is responsible for making the connection between a mesh
(described by a Triangulation object) and a finite element, by
allocating the correct number of degrees of freedom for each vertex,
The main idea is of course unchanged: we have to loop over all cells, assemble
local contributions, and then copy them into the global objects. As discussed
-in some detail first in @ref step_3 "step-3", deal.II has the FEValues class
+in some detail first in step-3, deal.II has the FEValues class
that pulls finite element description, mapping, and quadrature formula
together and aids in evaluating values and gradients of shape functions as
well as other information on each of the quadrature points mapped to the real
which are the cells where the error is largest, and then refine them. In many
of the other tutorial programs, we use the KellyErrorEstimator class to get an
indication of the size of the error on a cell, although we also discuss more
-complicated strategies in some programs, most importantly in @ref step_14 "step-14".
+complicated strategies in some programs, most importantly in step-14.
In any case, as long as the decision is only "refine this cell" or "do not
refine this cell", the actual refinement step is not particularly
lead to a program running out of memory when in fact there would be plenty of
memory for all computations.
-A solution to the problem has already been discussed in @ref step_11 "step-11"
-and @ref step_18 "step-18". It used an intermediate object of type
+A solution to the problem has already been discussed in step-11
+and step-18. It used an intermediate object of type
CompressedSparsityPattern. This class uses a different memory storage scheme
that is optimized to <i>creating</i> a sparsity pattern when maximal numbers
of entries per row are not accurately available, but is unsuitable for use as
time of creating the sparsity pattern, or while copying local contributions
into the global matrix object (or global vector).
-So, instead of the code snippet (taken from @ref step_6 "step-6")
+So, instead of the code snippet (taken from step-6)
@code
DoFTools::make_hanging_node_constraints (dof_handler,
hanging_node_constraints);
hanging_node_constraints.condense (system_rhs);
@endcode
We now replace copying and later condensing into one step, as already shown in
-@ref step_17 "step-17", @ref step_18 "step-18", and @ref step_22 "step-22":
+step-17, step-18, and step-22:
@code
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
interpolates the Dirichlet boundary conditions already in the setup phase in
order to tell the ConstraintMatrix object about them, and then do the
transfer from local to global data on matrix and vector simultaneously. This
-is exactly what we've shown in the @ref step_22 "step-22" tutorial program.
+is exactly what we've shown in the step-22 tutorial program.
<h3>The test case</h3>
The test case we will solve with this program is a re-take of the one we
-already look at in @ref step_14 "step-14": we solve the Laplace equation
+already look at in step-14: we solve the Laplace equation
@f[
-\Delta u = f
@f]
$u$. We do so on the domain $[-1,1]^2\backslash[-\frac 12,\frac 12]^2$, i.e. a
square with a square hole in the middle.
-The difference to @ref step_14 "step-14" is of course that we use $hp$ finite
+The difference to step-14 is of course that we use $hp$ finite
elements for the solution. The testcase is of interest because it has
re-entrant corners in the corners of the hole, at which the solution has
singularities. We therefore expect that the solution will be smooth in the
freedom, at least on the later grids when we have elements of relatively
high order (in 3d, the fraction of constrained degrees of freedom can be up
to 30%). This is, in fact, on the same order of magnitude as for non-$hp$
-discretizations. For example, in the last step of the @ref step_6 "step-6"
+discretizations. For example, in the last step of the step-6
program, we have 18401 degrees of freedom, 4104 of which are
constrained. The difference is that in the latter program, each constrained
hanging node is constrained against only the two adjacent degrees of
equations. This short example demonstrates how this can be implemented in
deal.II by using an <code>FE_system</code> object to stack two finite element
fields representing real and imaginary parts. We also revisit the
-ParameterHandler class first used in @ref step_19 "step-19", which provides a
+ParameterHandler class first used in step-19, which provides a
convenient way for reading parameters from a configuration file at runtime
without the need to recompile the program code.
c (n\cdot\nabla u) + i\,\omega\,u &=&0, \qquad x\in\Gamma_2,\\
u &=& 1, \qquad x\in\Gamma_1.
@f}
-This is a Helmholtz equation (similar to the one in @ref step_7 "step-7", but this time with ''the bad sign'') with Dirichlet data on $\Gamma_1$ and mixed boundary conditions on $\Gamma_2$. Because of the condition on $\Gamma_2$, we cannot just treat the equations for real and imaginary parts of $u$ separately. What we can do however is to view the PDE for $u$ as a system of two PDEs for the real and imaginary parts of $u$, with the boundary condition on $\Gamma_2$ representing the coupling terms between the two components of the system. This works along the following lines: Let $v=\textrm{Re}\;u,\; w=\textrm{Im}\;u$, then in terms of $v$ and $w$ we have the following system:
+This is a Helmholtz equation (similar to the one in step-7, but this time with ''the bad sign'') with Dirichlet data on $\Gamma_1$ and mixed boundary conditions on $\Gamma_2$. Because of the condition on $\Gamma_2$, we cannot just treat the equations for real and imaginary parts of $u$ separately. What we can do however is to view the PDE for $u$ as a system of two PDEs for the real and imaginary parts of $u$, with the boundary condition on $\Gamma_2$ representing the coupling terms between the two components of the system. This works along the following lines: Let $v=\textrm{Re}\;u,\; w=\textrm{Im}\;u$, then in terms of $v$ and $w$ we have the following system:
@f{eqnarray*}
\left.\begin{array}{ccc}
-\omega^2 v - c^2\Delta v &=& 0 \quad\\
resulting in 160x160 cells, and the output is written in gmv
format. The parameter reader understands many more parameters
pertaining in particular to the generation of output, see the
-explanation in @ref step_19 "step-19", but we need none of these
+explanation in step-19, but we need none of these
parameters here and therefore stick with their default values.
Here's the console output of the program in debug mode:
<li>
Observe convergence: We will only discuss computing errors in norms in
- @ref step_7 "step-7", but it is easy to check that computations converge
+ step-7, but it is easy to check that computations converge
already here. For example, we could evaluate the value of the solution in a
single point and compare the value for different %numbers of global
refinement (the number of global refinement steps is set in
This example is devoted to <em>anisotropic refinement</em>, which extends to
possibilities of local refinement. In most parts, this is a modification of the
-@ref step_12 "step-12" tutorial program, we use the same DG method for a linear transport
+step-12 tutorial program, we use the same DG method for a linear transport
equation. This program will cover the following topics:
<ol>
<li> <em>Anisotropic refinement</em>: What is the meaning of anisotropic refinement?
</ol>
The discretization itself will not be discussed, and neither will
implementation techniques not specific to anisotropic refinement used
-here. Please refer to @ref step_12 "step-12" for this.
+here. Please refer to step-12 for this.
Please note, at the moment of writing this tutorial program, anisotropic
refinement is only fully implemented for discontinuous Galerkin Finite
Another important aspect, and the most important one in this tutorial, is
the treatment of neighbor-relations when assembling jump terms on the
faces between cells. Looking at the documentation of the
-assemble_system functions in @ref step_12 "step-12" we notice, that we need to decide if a
+assemble_system functions in step-12 we notice, that we need to decide if a
neighboring cell is coarser, finer or on the same (refinement) level as our
current cell. These decisions do not work in the same way for anisotropic
refinement as the information given by the <em>level</em> of a cell is not
use <code>else if (!cell->coarser_neighbor(face_no))</code>. The code inside this
block can be left untouched. However, there is one thing to mention here: If
we want to use a rule, which cell should assemble certain terms on a given
- face we might think of the rule presented in @ref step_12 "step-12". We know that we have to
+ face we might think of the rule presented in step-12. We know that we have to
leave out the part about comparing our cell's level with that of the neighbor
and replace it with the test for a coarser neighbor presented above. However,
we also have to consider the possibility that neighboring cells of same
<li> Call Triangulation<dim>::execute_coarsening_and_refinement to perform the
requested refinement, using the requested isotropic and anisotropic flags.
</ol>
-This approach is similar to the one we have used in @ref step_27 "step-27"
+This approach is similar to the one we have used in step-27
for hp refinement and
has the great advantage of flexibility: Any error indicator can be
used in the anisotropic process, i.e. if you have quite involved a posteriori
<h3>The problem</h3>
-We solve the linear transport equation presented in @ref step_12 "step-12". The domain is
+We solve the linear transport equation presented in step-12. The domain is
extended to cover $[-1,1]\times[0,1]$ in 2D, where the flow field $\beta$ describes a
counterclockwise quarter circle around the origin in the right half of the
domain and is parallel to the x-axis in the left part of the domain. The inflow
boundary is again located at $x=1$ and along the positive part of the x-axis,
-and the boundary conditions are chosen as in @ref step_12 "step-12". Compared to @ref step_12 "step-12" we
+and the boundary conditions are chosen as in step-12. Compared to step-12 we
only use the more effective second assembling technique. In order to make
comparisons more effective, we decided to keep function names like @p
assemble_system2 even if there is only one of these routines in this tutorial
[(\nabla{\mathbf u}) + (\nabla {\mathbf u})^T]$ is the symmetric
gradient of the velocity. As can be seen, velocity and pressure
solve a Stokes equation describing the motion of an incompressible
-fluid, an equation we have previously considered in @ref step_22 "step-22"; we
+fluid, an equation we have previously considered in step-22; we
will draw extensively on the experience we have gained in that program, in
particular with regard to efficient linear Stokes solvers.
coefficient $\beta$,
the temperature <i>T</i> and the gravity vector <b>g</b> pointing
downward. (A derivation of why the right hand side looks like it looks
-is given in the introduction of @ref step_32 "step-32".)
+is given in the introduction of step-32.)
While the first two equations describe how the fluid reacts to
temperature differences by moving around, the third equation states
how the fluid motion affects the temperature field: it is an advection
<h3>Solution approach</h3>
-Like the equations solved in @ref step_21 "step-21", we here have a
+Like the equations solved in step-21, we here have a
system of differential-algebraic equations (DAE): with respect to the time
variable, only the temperature equation is a differential equation
whereas the Stokes system for <b>u</b> and <i>p</i> has no
time-derivatives and is therefore of the sort of an algebraic
constraint that has to hold at each time instant. The main difference
-to @ref step_21 "step-21" is that the algebraic constraint there was a
+to step-21 is that the algebraic constraint there was a
mixed Laplace system of the form
@f{eqnarray*}
\mathbf u + {\mathbf K}\lambda \nabla p &=& 0, \\
where $\nabla \cdot \eta \varepsilon (\cdot)$ is an operator similar to the
Laplacian $\Delta$ applied to a vector field.
-Given the similarity to what we have done in @ref step_21 "step-21",
+Given the similarity to what we have done in step-21,
it may not come as a surprise that we choose a similar approach,
although we will have to make adjustments for the change in operator
in the top-left corner of the differential operator.
<h4>Time stepping</h4>
The structure of the problem as a DAE allows us to use the same strategy as
-we have already used in @ref step_21 "step-21", i.e. we use a time lag
+we have already used in step-21, i.e. we use a time lag
scheme: we first solve the temperature equation (using an extrapolated
velocity field), and then insert the new temperature solution into the right
hand side of the velocity equation. The way we implement this in our code
and then the temperature equation with an extrapolated velocity field to
time <i>n</i>.
-In contrast to @ref step_21 "step-21", we'll use a higher order time
+In contrast to step-21, we'll use a higher order time
stepping scheme here, namely the <a
href="http://en.wikipedia.org/wiki/Backward_differentiation_formula">Backward
Differentiation Formula scheme of order 2 (BDF-2 in short)</a> that replaces
Like solving the mixed Laplace equations, solving the Stokes equations
requires us to choose particular pairs of finite elements for
velocities and pressure variables. Because this has already been discussed in
-@ref step_22 "step-22", we only cover this topic briefly:
+step-22, we only cover this topic briefly:
Here, we use the
stable pair $Q_{p+1}^d \times Q_p, p\ge 1$. These are continuous
elements, so we can form the weak form of the Stokes equation without
@f}
This is exactly the same as what we already discussed in
-@ref step_22 "step-22" and there is not much more to say about this here.
+step-22 and there is not much more to say about this here.
<h4>Stabilization, weak form and space discretization for the temperature equation</h4>
upwinding, stabilization, or all of this. One way to achieve this is
to use discontinuous elements (i.e. the FE_DGQ class that we used, for
example, in the discretization of the transport equation in
-@ref step_12 "step-12", or in discretizing the pressure in
-@ref step_20 "step-20" and @ref step_21 "step-21") and to define a
+step-12, or in discretizing the pressure in
+step-20 and step-21) and to define a
flux at the interface between cells that takes into account
upwinding. If we had a pure advection problem this would probably be
the simplest way to go. However, here we have some diffusion as well,
<h5>Linear solvers for the Stokes problem</h5>
Solving the linear equations coming from the Stokes system has been
-discussed in great detail in @ref step_22 "step-22". In particular, in
+discussed in great detail in step-22. In particular, in
the results section of that program, we have discussed a number of
alternative linear solver strategies that turned out to be more
efficient than the original approach. The best alternative
matrices. In particular, it turned out that <i>S</i> is spectrally
equivalent to the mass matrix and consequently replacing $\tilde
S^{-1}$ by a CG solver applied to the mass matrix on the pressure
-space was a good choice. In a small deviation from @ref step_22 "step-22", we
+space was a good choice. In a small deviation from step-22, we
here have a coefficient $\eta$ in the momentum equation, and by the same
derivation as there we should arrive at the conclusion that it is the weighted
mass matrix with entries $\tilde S_{ij}=(\eta^{-1}\varphi_i,\varphi_j)$ that
the vector-valued velocity field, i.e.
$A_{ij} = (\varepsilon {\mathbf v}_i, 2\eta \varepsilon ({\mathbf
v}_j))$.
-In @ref step_22 "step-22" we used a sparse LU decomposition (using the
+In step-22 we used a sparse LU decomposition (using the
SparseDirectUMFPACK class) of <i>A</i> for $\tilde A^{-1}$ — the
perfect preconditioner — in 2d, but for 3d memory and compute
time is not usually sufficient to actually compute this decomposition;
To sum this whole story up, we can observe:
<ul>
<li> Compared to building a preconditioner from the original matrix <i>A</i>
- resulting from the symmetric gradient as we did in @ref step_22 "step-22",
+ resulting from the symmetric gradient as we did in step-22,
we have to expect that the preconditioner based on the Laplace bilinear form
performs worse since it does not take into account the coupling between
vector components.
finite element with several vector components, and a single DoFHandler object.
Sometimes, we have substructured the resulting matrix into blocks to
facilitate particular solver schemes; this was, for example, the case in the
-@ref step_22 "step-22" program for the Stokes equations upon which the current
+step-22 program for the Stokes equations upon which the current
program is based.
We could of course do the same here. The linear system that we would get would
<h4>Using Trilinos</h4>
In much the same way as we used PETSc to support our linear algebra needs in
-@ref step_17 "step-17" and @ref step_18 "step-18", we use interfaces to the <a
+step-17 and step-18, we use interfaces to the <a
href="http://trilinos.sandia.gov">Trilinos</a> library (see @ref
SoftwareTrilinos) in this program. Trilinos is a very large collection of
everything that has to do with linear and nonlinear algebra, as well as all
the Trilinos ML package that implements an Algebraic Multigrid (AMG)
method. We will use this preconditioner to precondition the second order
operator part of the momentum equation. The ability to solve problems in
-%parallel will be explored in @ref step_32 "step-32", using the same problem as
+%parallel will be explored in step-32, using the same problem as
discussed here.
boundary, we will require no-normal flux for the velocity
($\mathrm{n}\cdot\mathrm{u}=0$) and for the temperature
($\mathrm{n}\cdot\nabla T=0$). This is one of the cases discussed in the
-introduction of @ref step_22 "step-22" and fixes one component of the velocity
+introduction of step-22 and fixes one component of the velocity
while allowing flow to be %parallel to the boundary. There remain
<code>dim-1</code> components to be fixed, namely the tangential components of
the normal stress; for these, we choose homogenous conditions which means that
this value for $c_R$ appears to work just fine for the current
program, we corrected the formula in the program and set $c_R$ to a
value that reproduces exactly the results we had before. We will,
-however, revisit this issue again in @ref step_32 "step-32".
+however, revisit this issue again in step-32.
Now, however, back to the discussion of what values of $c_k$ and
$\beta$ to choose:
There are various ways to extend the current program. Of particular interest
is, of course, to make it faster and/or increase the resolution of the
-program, in particular in 3d. This is the topic of the @ref step_32 "step-32"
+program, in particular in 3d. This is the topic of the step-32
tutorial program which will implement strategies to solve this problem in
%parallel on a cluster.
<a name="Intro"></a>
<h1>Introduction</h1>
-This program does pretty much exactly what @ref step_31 "step-31" already
+This program does pretty much exactly what step-31 already
does: it solves the Boussinesq equations that describe the motion of a fluid
whose temperature is not in equilibrium. As such, all the equations we have
-described in @ref step_31 "step-31" still hold: we solve the same partial
+described in step-31 still hold: we solve the same partial
differential equation, using the same finite element scheme, the same time
stepping algorithm, and the same stabilization method for the temperature
advection-diffusion equation. As a consequence, you may first want to
understand that program before you work on the current one.
-The difference between @ref step_31 "step-31" and the current program is that
+The difference between step-31 and the current program is that
here we want to do things in %parallel, using both the availability of many
machines in a cluster (with parallelization based on MPI) as well as many
processor cores within a single machine (with parallelization based on
<h3> Using the "right" pressure </h3>
-In @ref step_31 "step-31", we used the following Stokes model for the
+In step-31, we used the following Stokes model for the
velocity and pressure field:
@f{eqnarray*}
-\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=&
\\
\nabla \cdot {\mathbf u} &=& 0.
@f}
-This is exactly the form we used in @ref step_31 "step-31", and it was
+This is exactly the form we used in step-31, and it was
appropriate to do so because all changes in the fluid flow are only
driven by the dynamic pressure that results from temperature
differences.
\nabla \cdot \kappa \nabla T &=& \gamma,
@f}
augmented by appropriate boundary and initial conditions. As discussed
-in @ref step_31 "step-31", we will solve this set of equations by
+in step-31, we will solve this set of equations by
solving for a Stokes problem first in each time step, and then moving
the temperature equation forward by one time interval.
quite obviously, does not make sense, and we should not be surprised
that doing so is eventually going to come back hurting us.
-So why is this an issue here, but not in @ref step_31 "step-31"? The
+So why is this an issue here, but not in step-31? The
reason back there is that everything was nicely balanced: velocities
were on the order of one, the pressure likewise, the viscosity was
one, and the domain had a diameter of $\sqrt{2}$. As a result, while
<h3> Changes to the Stokes preconditioner </h3>
In this tutorial program, we apply a variant of the preconditioner used in
-@ref step_31 "step-31". That preconditioner was built to operate on the
+step-31. That preconditioner was built to operate on the
system matrix <i>M</i> in block form such that the product matrix
@f{eqnarray*}
P^{-1} M
equations. BiCGStab harmonizes nicely with the ILU preconditioner on a
pressure mass matrix as approximation for $S^{-1}$, so it is at least as
good as GMRES in this example. Keep in mind the discussion in the results
-section of the @ref step_22 "step-22" tutorial program, where we observed
+section of the step-22 tutorial program, where we observed
that BiCGStab does <i>not</i> like inner solves with CG, which made us
prefer GMRES in step-31.
-As a final note, let us remark that in @ref step_31 "step-31" we computed the
+As a final note, let us remark that in step-31 we computed the
Schur complement $S=B A^{-1} B^T$ by approximating
$-\text{div}(-\eta\Delta)^{-1}\nabla \approx \frac 1{\eta} \mathbf{1}$. Now,
however, we have re-scaled the $B$ and $B^T$ operators. So $S$ should now
<h3> Changes to the artificial viscosity stabilization </h3>
-As in @ref step_31 "step-31", we will use an artificial viscosity of
+As in step-31, we will use an artificial viscosity of
the form
@f{eqnarray*}
\nu_\alpha(T)|_K
$c(\mathbf{u},T) =
c_R\ \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T)
\ |\mathrm{diam}(\Omega)|^{\alpha-2}$ (for the meaning of the various
-terms in these formulas, see @ref step_31 "step-31". In the results
+terms in these formulas, see step-31. In the results
section of that program, we have discussed our choice for $c_R$ and
how we arrived at the value used there mostly by accident, and in more
detail how $\beta$ was chosen. For the
Parallelization of scientific codes across multiple machines in a cluster of
computers is almost always done using the Message Passing Interface
(MPI). This program is no exception to that, and it follows the
-@ref step_17 "step-17" and @ref step_18 "step-18" programs in this.
+step-17 and step-18 programs in this.
MPI is a rather awkward interface to program with, and so we usually try to
not use it directly but through an interface layer that abstracts most of the
of the same functionality. It is, however, a project that is several years
younger, is written in C++ and by people who generally have put a significant
emphasis on software design. We have already used Trilinos in
-@ref step_31 "step-31", and will do so again here, with the difference that we
+step-31, and will do so again here, with the difference that we
will use its %parallel capabilities.
deal.II's Trilinos interfaces encapsulate pretty much everything Trilinos
<h3> The testcase </h3>
The setup for this program is mildly reminiscent of the problem we wanted to
-solve in the first place (see the introduction of @ref step_31 "step-31"):
+solve in the first place (see the introduction of step-31):
convection in the earth mantle. As a consequence, we choose the following
data, all of which appears in the program in units of meters and seconds (the
SI system) even if we list them here in other units.
<h1>Results</h1>
When run, the program simulates convection in 3d in much the same way
-as @ref step_31 "step-31" did, though with an entirely different testcase.
+as step-31 did, though with an entirely different testcase.
<h3>Comparison of results with step-31</h3>
<h3>Possibilities for extensions</h3>
Apart from the various possibilities for extensions already outlined
-in the @ref step_31 "step-31", here are a few more ideas:
+in the step-31, here are a few more ideas:
<ul>
<li> The temperature field we get in our simulations after a while
Discretization happens in the usual way, taking into account that this
is a hyperbolic problem in the same style as the simple one discussed
-in @ref step_12 "step-12":
+in step-12:
We choose a finite element space $V_h$, and integrate our conservation law against
our (vector-valued) test function $\mathbf{z} \in V_h$. We then integrate by parts and approximate the
boundary flux with a <i> numerical </i> flux $\mathbf{H}$,
assuming that the effects due to viscosity are negligible compared to
those due to the pressure gradient, inertial forces and the external
forces. This is the opposite case of the Stokes equations discussed in
-@ref step_22 "step-22" which are the limit case of dominant viscosity,
+step-22 which are the limit case of dominant viscosity,
i.e. where the velocity is so small that inertia forces can be
neglected. On the other hand, owing to the assumed incompressibility,
the equations are not suited for very high speed gas flows where
where $n$ is the outer unit normal. The boundary conditions on $\Gamma_2$ are often
used to model outflow conditions.
-In previous tutorial programs (see for instance @ref step_20 "step-20" and
-@ref step_22 "step-22") we have seen
+In previous tutorial programs (see for instance step-20 and
+step-22) we have seen
how to solve the time-independent Stokes equations using a Schur complement approach. For the
time-dependent case, after time discretization, we would arrive at a system like
@f{align*}
<h3>Implementation details</h3>
The program below is essentially just a slightly modified version of
-@ref step_4 "step-4". The things that are different are the following:
+step-4. The things that are different are the following:
<ul>
longer has a single solution vector, but a whole set of vectors for
the various eigenfunctions we want to compute.</li>
-<li>We use PETSc matrices and vectors as in @ref step_17 "step-17" and
-@ref step_18 "step-18" since that is what the SLEPc eigenvalue solvers
+<li>We use PETSc matrices and vectors as in step-17 and
+step-18 since that is what the SLEPc eigenvalue solvers
require.</li>
<li>The function <code>EigenvalueProblem::solve</code> is entirely
further with this but stop at making only a few of the things that one
could select at run time actual input file parameters. In order to see
what could be done in this regard, take a look at @ref step_29
-"step-29", @ref step_33 "step-33", and in particular @ref step_19
+"step-29", step-33, and in particular @ref step_19
"step-19".</li>
<li>We use the FunctionParser class to make the potential $V(\mathbf
</ul>
The rest of the program follows in a pretty straightforward way from
-@ref step_4 "step-4".
+step-4.
an input file, but we can also do so right away when creating the output
file. The facility to do that is the DataPostprocessor class that can
be used in conjunction with the DataOut class. Examples of how this
-can be done can be found in @ref step_29 "step-29" and
-@ref step_33 "step-33".
+can be done can be found in step-29 and
+step-33.
<li> What happens if the particle in the box has %internal degrees of
freedom? For example, if the particle were a spin-$1/2$ particle? In
decreases efficiency.
The implementation into the multigrid framework is then straightforward. We
-will only need some minor modifications compared to @ref step_16 "step-16".
+will only need some minor modifications compared to step-16.
<h3>The test case</h3>
<h3>Program output</h3>
-Since this example solves the same problem as @ref step_5 "step-5" (except for
+Since this example solves the same problem as step-5 (except for
a different coefficient), we refer to the graphical output there. Here, we
evaluate some aspects of the multigrid solver.
the inflow part of the boundary of the domain and ${\bf n}$ denotes
the unit outward normal to the boundary $\Gamma$. This equation is the
conservative version of the advection equation already considered in
-@ref step_9 "step-9" of this tutorial.
+step-9 of this tutorial.
In particular, we solve the advection equation on
$\Omega=[0,1]^2$ with ${\mathbf \beta}=\frac{1}{|x|}(-x_2, x_1)$
representing a circular counterclockwise flow field, and $g=1$
number of CG iterations needed (available through the
<code>solver_control.last_step()</code> call, see @ref step_4
"step-4") as well as CPU time needed (using the Timer class,
-discussed, for example, in @ref step_12 "step-12") and get the
+discussed, for example, in step-12) and get the
following results (left: iterations; right: CPU time):
<TABLE WIDTH="60%" ALIGN="center">
The message to take away from this is not that simplicity in
preconditioners is always best. While this may be true for the current
problem, it definitely is not once we move to more complicated
-problems (elasticity or Stokes, for examples @ref step_8 "step-8" or
-@ref step_22 "step-22"). Secondly, all of these preconditioners still
+problems (elasticity or Stokes, for examples step-8 or
+step-22). Secondly, all of these preconditioners still
lead to an increase in the number of iterations as the number $N$ of
degrees of freedom grows, for example ${\cal O}(N^\alpha)$; this, in
turn, leads to a total growth in effort as ${\cal O}(N^{1+\alpha})$
undesirable, we would really like to solve linear systems with $N$
unknowns in a total of ${\cal O}(N)$ work; there is a class
of preconditioners that can achieve this, namely geometric (@ref
-step_16 "step-16") or algebraic (@ref step_31 "step-31") multigrid
+step_16 "step-16") or algebraic (step-31) multigrid
preconditioners. They are, however, significantly more complex than
the preconditioners outlined above. Finally, the last message to take
home is that today (in 2008), linear systems with 100,000 unknowns are
This is the form in which we will implement the local stiffness matrix and
right hand side vectors.
-As a final note: in the @ref step_17 "step-17" example program, we will
+As a final note: in the step-17 example program, we will
revisit the elastic problem laid out here, and will show how to solve it in
%parallel on a cluster of computers. The resulting program will thus be able to
solve this problem to significantly higher accuracy, and more efficiently if
-this is required. In addition, in @ref step_20 "step-20", @ref step_21
+this is required. In addition, in step-20, @ref step_21
"step-21", as well as a few other of the later tutorial programs, we will
revisit some vector-valued problems and show a few techniques that may make it
simpler to actually go through all the stuff shown above, with
fact not two independent scalar fields, but the components of a vector
field. Depending on what graphical file format we use, we could as
well have specified this in the program already and then written into
-the output file. The @ref step_22 "step-22" shows how to do this in
+the output file. The step-22 shows how to do this in
the output_results() function.
The next image shows the final grid after eight steps of refinement:
* There are several, exchangeable variations of this class, see @ref Sparsity,
* section 'Dynamic block sparsity patterns' for more information.
*
- * This class is used in @ref step_22 "step-22" and @ref step_31 "step-31".
+ * This class is used in step-22 and step-31.
*
* @author Timo Heister, 2008
*/
* (in the sense that the method compress() needs to be called before the
* pattern can be used).
*
- * This class is used in @ref step_32 "step-32".
+ * This class is used in step-32.
*
* @author Martin Kronbichler, 2008, 2009
*/
* module.
*
* This class is an example of the "dynamic" type of @ref Sparsity. It
- * is discussed in the @ref step_27 "step-27" and @ref step_22
+ * is discussed in the step-27 and @ref step_22
* "step-22" tutorial programs.
*
* <h3>Interface</h3>
* sp.copy_from (compressed_pattern);
* @endverbatim
*
- * See also @ref step_11 "step-11" and @ref step_18 "step-18" for usage
+ * See also step-11 and step-18 for usage
* patterns of the related CompressedSparsityPattern class, and
- * @ref step_27 "step-27" of the current class.
+ * step-27 of the current class.
*
* <h3>Notes</h3>
*
* container is used. This solution might not be the fastest in all
* situations, but seems to work much better than the
* CompressedSparsityPattern in the context of hp-adaptivity (see for
- * example @ref step_27 "step-27"), or generally when there are many
+ * example step-27), or generally when there are many
* nonzero entries in each row of a matrix (see @ref step_22
* "step-22"). On the other hand, a benchmark where nonzero entries
* were randomly inserted into the sparsity pattern revealed that this
* sp.copy_from (compressed_pattern);
* @endverbatim
*
- * See also @ref step_11 "step-11" and @ref step_18 "step-18" for usage
+ * See also step-11 and step-18 for usage
* patterns.
*
* <h3>Notes</h3>
* CompressedSetSparsityPattern class, as
* this turned out to be more efficient
* for hp finite element programs such as
- * @ref step_27 "step-27"). Using the
+ * step-27). Using the
* present scheme, we only need to
* allocate memory once for every 8 added
* entries, and we waste a lot less
* constraints occur as hanging nodes even if all used finite elements are
* identical. While they are most frequent for hanging nodes, constraints of
* the given form appear also in other contexts, see for example the
- * application the @ref step_11 "step-11" tutorial program.
+ * application the step-11 tutorial program.
*
* Homogenous constraints of this form also arise in the context of vector-valued
* fields, for example if one wants to enforce boundary conditions of the form
* the constraint. This is exactly the form we need for hanging node
* constraints, where we need to constrain one degree of freedom in terms of
* others. There are other conditions of this form possible, for example for
- * implementing mean value conditions as is done in the @ref step_11 "step-11"
+ * implementing mean value conditions as is done in the step-11
* tutorial program. The name of the class stems from the fact that these
* constraints can be represented in matrix form as $X x = b$, and this object
* then describes the matrix $X$ (as well as, incidentally, the vector $b$ --
* constraints where $b=0$, thus the name). The most frequent way to
* create/fill objects of this type is using the
* DoFTools::make_hanging_node_constraints() function. The use of these
- * objects is first explained in @ref step_6 "step-6".
+ * objects is first explained in step-6.
*
* Matrices of the present type are organized in lines (rows), but only those
* lines are stored where constraints are present. New constraints are added
* before).
*
* The use of ConstraintMatrix for implementing Dirichlet boundary conditions
- * is discussed in the @ref step_22 "step-22" tutorial program.
+ * is discussed in the step-22 tutorial program.
*
*
* <h3>Avoiding explicit condensation</h3>
* is typically the case for 3d, or high polynomial degree computations, as
* well as for hp finite element methods, see for example the @ref hp_paper
* "hp paper". This is the case discussed in the hp tutorial program, @ref
- * step_27 "step-27", as well as in @ref step_22 "step-22" and @ref step_31
+ * step_27 "step-27", as well as in step-22 and @ref step_31
* "step-31".
*
* <li>
* is, for example, the case for the PETSc and Trilinos wrapper classes,
* where we have no access to the underlying representation of the matrix,
* and therefore cannot efficiently implement the condense()
- * operation). This is the case discussed in @ref step_17 "step-17", @ref
- * step_18 "step-18", @ref step_31 "step-31", and @ref step_32 "step-32".
+ * operation). This is the case discussed in step-17, @ref
+ * step_18 "step-18", step-31, and step-32.
* </ul>
*
* In this case, one possibility is to distribute local entries to the final
*
* As explained in the
* @ref hp_paper "hp paper"
- * and in @ref step_27 "step-27",
+ * and in step-27,
* first allocating a sparsity pattern
* and later coming back and allocating
* additional entries for those matrix
PreconditionIdentity());
@endverbatim
*
- * See the @ref step_3 "step-3" tutorial program for an example and
+ * See the step-3 tutorial program for an example and
* additional explanations.
*
* Alternatively, the IdentityMatrix class can be used to precondition
* for the generalized eigenvalue problem $Ax=B\lambda x$, where the
* variable <code>const unsigned int size_of_spectrum</code> tells
* SLEPc the number of eigenvector/eigenvalue pairs to solve for: See
- * also @ref step_36 "step-36" for a hands-on example.
+ * also step-36 for a hands-on example.
*
* An alternative implementation to the one above is to use the API
* internals directly within the application code. In this way the
* preconditioner to be based on. A call of the respective
* <code>vmult</code> function does call the respective operation in the
* Trilinos package, where it is called <code>ApplyInverse</code>. Use of
- * this class is explained in the @ref step_31 "step-31" tutorial program.
+ * this class is explained in the step-31 tutorial program.
*
* Since the Trilinos objects we want to use are heavily dependent on Epetra
* objects, we recommend using this class in conjunction with Trilinos