The step-18 tutorial program

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

Introduction

This tutorial program is another one in the series on the elasticity problem 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 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 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 main documentation page in your installation of deal.II, or on the deal.II webpage.

Quasistatic elastic deformation

Motivation of the model

In general, time-dependent small elastic deformations are described by the elastic wave equation

\[ \rho \frac{\partial^2 \mathbf{u}}{\partial t^2} + c \frac{\partial \mathbf{u}}{\partial t} - \textrm{div}\ ( C \varepsilon(\mathbf{u})) = \mathbf{f} \qquad \textrm{in}\ \Omega, \]

where $\mathbf{u}=\mathbf{u} (\mathbf{x},t)$ is the deformation of the body, $\rho$ and $c$ the density and attenuation coefficient, and $\mathbf{f}$ external forces. In addition, initial conditions

\[ \mathbf{u}(\cdot, 0) = \mathbf{u}_0(\cdot) \qquad \textrm{on}\ \Omega, \]

and Dirichlet (displacement) or Neumann (traction) boundary conditions need to be specified for a unique solution:

\begin{eqnarray*} \mathbf{u}(\mathbf{x},t) &=& \mathbf{d}(\mathbf{x},t) \qquad \textrm{on}\ \Gamma_D\subset\partial\Omega, \\ \mathbf{n} \ C \varepsilon(\mathbf{u}(\mathbf{x},t)) &=& \mathbf{b}(\mathbf{x},t) \qquad \textrm{on}\ \Gamma_N=\partial\Omega\backslash\Gamma_D. \end{eqnarray*}

In above formulation, $\varepsilon(\mathbf{u})= \frac 12 (\nabla \mathbf{u} + \nabla \mathbf{u}^T)$ is the symmetric gradient of the displacement, also called the strain. $C$ is a tensor of rank 4, called the stress-strain tensor that contains knowledge of the elastic strength of the material; its symmetry properties make sure that it maps symmetric tensors of rank 2 (“matrices” of dimension $d$, where $d$ is the spatial dimensionality) onto symmetric tensors of the same rank. We will comment on the roles of the strain and stress tensors more below. For the moment it suffices to say that we interpret the term $\textrm{div}\ ( C \varepsilon(\mathbf{u}))$ as the vector with components $\frac \partial{\partial x_j} C_{ijkl} \varepsilon(\mathbf{u})_{kl}$, where summation over indices $j,k,l$ is implied.

The quasistatic limit of this equation is motivated as follows: each small perturbation of the body, for example by changes in boundary condition or the forcing function, will result in a corresponding change in the configuration of the body. In general, this will be in the form of waves radiating away from the location of the disturbance. Due to the presence of the damping term, these waves will be attenuated on a time scale of, say, $\tau$. Now, assume that all changes in external forcing happen on times scales that are much larger than $\tau$. In that case, the dynamic nature of the change is unimportant: we can consider the body to always be in static equilibrium, i.e. we can assume that at all times the body satisfies

\begin{eqnarray*} - \textrm{div}\ ( C \varepsilon(\mathbf{u})) &=& \mathbf{f}(\mathbf{x},t) \qquad \textrm{in}\ \Omega, \\ \mathbf{u}(\mathbf{x},t) &=& \mathbf{d}(\mathbf{x},t) \qquad \textrm{on}\ \Gamma_D, \\ \mathbf{n} \ C \varepsilon(\mathbf{u}(\mathbf{x},t)) &=& \mathbf{b}(\mathbf{x},t) \qquad \textrm{on}\ \Gamma_N. \end{eqnarray*}

Note that the differential equation does not contain any time derivatives any more -- all time dependence is introduced through boundary conditions and a possibly time-varying force function $\mathbf{f}(\mathbf{x},t)$. The changes in configuration can therefore be considered as being stationary instantaneously. An alternative view of this is that $t$ is not really a time variable, but only a time-like parameter that governs the evolution of the problem.

While these equations are sufficient to describe small deformations, computing large deformations is a little more complicated. To do so, let us first introduce a tensorial stress variable $\sigma$, and write the differential equations in terms of the stress:

\begin{eqnarray*} - \textrm{div}\ \sigma &=& \mathbf{f}(\mathbf{x},t) \qquad \textrm{in}\ \Omega(t), \\ \mathbf{u}(\mathbf{x},t) &=& \mathbf{d}(\mathbf{x},t) \qquad \textrm{on}\ \Gamma_D\subset\partial\Omega(t), \\ \mathbf{n} \ C \varepsilon(\mathbf{u}(\mathbf{x},t)) &=& \mathbf{b}(\mathbf{x},t) \qquad \textrm{on}\ \Gamma_N=\partial\Omega(t)\backslash\Gamma_D. \end{eqnarray*}

Note that these equations are posed on a domain $\Omega(t)$ that changes with time, with the boundary moving according to the displacements $\mathbf{u}(\mathbf{x},t)$ of the points on the boundary. To complete this system, we have to specify the incremental relationship between the stress and the strain, as follows:

\[ \dot\sigma = C \varepsilon (\dot{\mathbf{u}}), \qquad \qquad \textrm{[stress-strain]} \]

where a dot indicates a time derivative. Both the stress $\sigma$ and the strain $\varepsilon(\mathbf{u})$ are symmetric tensors of rank 2.

Time discretization

Numerically, this system is solved as follows: first, we discretize the time component using a backward Euler scheme. This leads to a discrete equilibrium of force at time step $n$:

\[ -\textrm{div}\ \sigma^n = f^n, \]

where

\[ \sigma^n = \sigma^{n-1} + C \varepsilon (\Delta \mathbf{u}^n), \]

and $\Delta \mathbf{u}^n$ the incremental displacement for time step $n$. In addition, we have to specify initial data $\mathbf{u}(\cdot,0)=\mathbf{u}_0$. This way, if we want to solve for the displacement increment, we have to solve the following system:

\begin{eqnarray*} - \textrm{div}\ C \varepsilon(\Delta\mathbf{u}^n) &=& \mathbf{f} + \textrm{div}\ \sigma^{n-1} \qquad \textrm{in}\ \Omega(t_{n-1}), \\ \Delta \mathbf{u}^n(\mathbf{x},t) &=& \mathbf{d}(\mathbf{x},t_n) - \mathbf{d}(\mathbf{x},t_{n-1}) \qquad \textrm{on}\ \Gamma_D\subset\partial\Omega(t_{n-1}), \\ \mathbf{n} \ C \varepsilon(\Delta \mathbf{u}^n(\mathbf{x},t)) &=& \mathbf{b}(\mathbf{x},t_n)-\mathbf{b}(\mathbf{x},t_{n-1}) \qquad \textrm{on}\ \Gamma_N=\partial\Omega(t_{n-1})\backslash\Gamma_D. \end{eqnarray*}

The weak form of this set of equations, which as usual is the basis for the finite element formulation, reads as follows: find $\Delta \mathbf{u}^n \in \{v\in H^1(\Omega(t_{n-1}))^d: v|_{\Gamma_D}=\mathbf{d}(\cdot,t_n) - \mathbf{d}(\cdot,t_{n-1})\}$ such that

\begin{eqnarray*} (C \varepsilon(\Delta\mathbf{u}^n), \varepsilon(\varphi) )_{\Omega(t_{n-1})} = (\mathbf{f}, \varphi)_{\Omega(t_{n-1})} -(\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})} \\ +(\mathbf{b}(\mathbf{x},t_n)-\mathbf{b}(\mathbf{x},t_{n-1}), \varphi)_{\Gamma_N} \\ \forall \varphi \in \{\mathbf{v}\in H^1(\Omega(t_{n-1}))^d: \mathbf{v}|_{\Gamma_D}=0\}. \qquad \qquad \textrm{[linear-system]} \end{eqnarray*}

We note that, for simplicity, in the program we will always assume that there are no boundary forces, i.e. $\mathbf{b} = 0$, and that the deformation of the body is driven by body forces $\mathbf{f}$ and prescribed boundary displacements $\mathbf{d}$ alone. It is also worth noting that when integrating by parts, we would get terms of the form $(C \varepsilon(\Delta\mathbf{u}^n), \nabla \varphi )_{\Omega(t_{n-1})}$, but that we replace it with the term involving the symmetric gradient $\varepsilon(\varphi)$ instead of $\nabla\varphi$. Due to the symmetry of $C$, the two terms are equivalent, but the symmetric version avoids a potential for round-off to render the resulting matrix slightly non-symmetric.

The system at time step $n$, to be solved on the old domain $\Omega(t_{n-1})$, has exactly the form of a stationary elastic problem, and is therefore similar to what we have already implemented in previous example programs. We will therefore not comment on the space discretization beyond saying that we again use lowest order continuous finite elements.

There are differences, however:

  1. We have to move (update) the mesh after each time step, in order to be able to solve the next time step on a new domain;

  2. We need to know $\sigma^{n-1}$ to compute the next incremental displacement, i.e. we need to compute it at the end of the time step to make sure it is available for the next time step. Essentially, the stress variable is our window to the history of deformation of the body.
These two operations are done in the functions move_mesh and update_quadrature_point_history in the program. While moving the mesh is only a technicality, updating the stress is a little more complicated and will be discussed in the next section.

Updating the stress variable

As indicated above, we need to have the stress variable $\sigma^n$ available when computing time step $n+1$, and we can compute it using

\[ \sigma^n = \sigma^{n-1} + C \varepsilon (\Delta \mathbf{u}^n). \qquad \qquad \textrm{[stress-update]} \]

There are, despite the apparent simplicity of this equation, two questions that we need to discuss. The first concerns the way we store $\sigma^n$: even if we compute the incremental updates $\Delta\mathbf{u}^n$ using lowest-order finite elements, then its symmetric gradient $\varepsilon(\Delta\mathbf{u}^n)$ is in general still a function that is not easy to describe. In particular, it is not a piecewise constant function, and on general meshes (with cells that are not rectangles parallel to the coordinate axes) or with non-constant stress-strain tensors $C$ it is not even a bi- or trilinear function. Thus, it is a priori not clear how to store $\sigma^n$ in a computer program.

To decide this, we have to see where it is used. The only place where we require the stress is in the term $(\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})}$. In practice, we of course replace this term by numerical quadrature:

\[ (\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})} = \sum_{K\subset {T}} (\sigma^{n-1},\varepsilon(\varphi))_K \approx \sum_{K\subset {T}} \sum_q w_q \ \sigma^{n-1}(\mathbf{x}_q) : \varepsilon(\varphi(\mathbf{x}_q), \]

where $w_q$ are the quadrature weights and $\mathbf{x}_q$ the quadrature points on cell $K$. This should make clear that what we really need is not the stress $\sigma^{n-1}$ in itself, but only the values of the stress in the quadrature points on all cells. This, however, is a simpler task: we only have to provide a data structure that is able to hold one symmetric tensor of rank 2 for each quadrature point on all cells (or, since we compute in parallel, all quadrature points of all cells that the present MPI process “owns”). At the end of each time step we then only have to evaluate $\varepsilon(\Delta \mathbf{u}^n(\mathbf{x}_q))$, multiply it by the stress-strain tensor $C$, and use the result to update the stress $\sigma^n(\mathbf{x}_q)$ at quadrature point $q$.

The second complication is not visible in our notation as chosen above. It is due to the fact that we compute $\Delta u^n$ on the domain $\Omega(t_{n-1})$, and then use this displacement increment to both update the stress as well as move the mesh nodes around to get to $\Omega(t_n)$ on which the next increment is computed. What we have to make sure, in this context, is that moving the mesh does not only involve moving around the nodes, but also making corresponding changes to the stress variable: the updated stress is a variable that is defined with respect to the coordinate system of the material in the old domain, and has to be transferred to the new domain. The reason for this can be understood as follows: locally, the incremental deformation $\Delta\mathbf{u}$ can be decomposed into three parts, a linear translation (the constant part of the displacement increment field in the neighborhood of a point), a dilational component (that part of the gradient of the displacement field that has a nonzero divergence), and a rotation. A linear translation of the material does not affect the stresses that are frozen into it -- the stress values are simply translated along. The dilational or compressional change produces a corresponding stress update. However, the rotational component does not necessarily induce a nonzero stress update (think, in 2d, for example of the situation where $\Delta\mathbf{u}=(y, -x)^T$, with which $\varepsilon(\Delta \mathbf{u})=0$). Nevertheless, if the the material was pre-stressed in a certain direction, then this direction will be rotated along with the material. To this end, we have to define a rotation matrix $R(\Delta \mathbf{u}^n)$ that describes, in each point the rotation due to the displacement increments. It is not hard to see that the actual dependence of $R$ on $\Delta \mathbf{u}^n$ can only be through the curl of the displacement, rather than the displacement itself or its full gradient (as mentioned above, the constant components of the increment describe translations, its divergence the dilational modes, and the curl the rotational modes). Since the exact form of $R$ is cumbersome, we only state it in the program code, and note that the correct updating formula for the stress variable is then

\[ \sigma^n = R(\Delta \mathbf{u}^n)^T [\sigma^{n-1} + C \varepsilon (\Delta \mathbf{u}^n)] R(\Delta \mathbf{u}^n). \qquad \qquad \textrm{[stress-update+rot]} \]

Both stress update and rotation are implemented in the function update_quadrature_point_history of the example program.

Parallel graphical output

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 problem were implemented in the library at the time, and the problem was deferred to a later version.

This functionality has been implemented in the meantime, and this is the time to explain its use. Basically, what we need to do is let every process generate graphical output for that subset of cells that it owns, write them into separate files and have a way to merge them later on. At this point, it should be noted that none of the graphical output formats known to the author of this program allows for a simple way to later re-read it and merge it with other files corresponding to the same simulation. What deal.II therefore offers is the following: When you call the DataOut::build_patches function, an intermediate format is generated that contains all the information for the data on each cell. Usually, this intermediate format is then further processed and converted into one of the graphical formats that we can presently write, such as gmv, eps, ucd, gnuplot, or a number of other ones. Once written in these formats, there is no way to reconstruct the necessary information to merge multiple blocks of output. However, the base classes of DataOut also allow to simply dump the intermediate format to a file, from which it can later be recovered without loss of information.

This has two advantages: first, simulations may just dump the intermediate format data during run-time, and the user may later decide which particular graphics format she wants to have. This way, she does not have to re-run the entire simulation if graphical output is requested in a different format. One typical case is that one would like to take a quick look at the data with gnuplot, and then create high-quality pictures using GMV or OpenDX. Since both can be generated out of the intermediate format without problem, there is no need to re-run the simulation.

In the present context, of more interest is the fact that in contrast to any of the other formats, it is simple to merge multiple files of intermediate format, if they belong to the same simulation. This is what we will do here: we will generate one output file in intermediate format for each processor that belongs to this computation (in the sequential case, this will simply be a single file). They may then later be read in and merged so that we can output a single file in whatever graphical format is requested.

The way to do this is to first instruct the DataOutBase class to write intermediate format rather than in gmv or any other graphical format. This is simple: just use data_out.write_deal_II_intermediate. We will write to a file called solution-TTTT.TTTT.d2 if there is only one processor, or files solution-TTTT.TTTT.NNN.d2 if this is really a parallel job. Here, TTTT.TTTT denotes the time for which this output has been generated, and NNN 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 step-19 tutorial program: for example, for the first time step, call it through

    ../step-19/step-19 solution-0001.0000.*.d2 solution-0001.0000.gmv
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 step-19 tutorial program.

Overall structure of the program

The overall structure of the program can be inferred from the run() function that first calls do_initial_timestep() for the first time step, and then do_timestep() on all subsequent time steps. The difference between these functions is only that in the first time step we start on a coarse mesh, solve on it, refine the mesh adaptively, and then start again with a clean state on that new mesh. This procedure gives us a better starting mesh, although we should of course keep adapting the mesh as iterations proceed -- this isn't done in this program, but commented on below.

The common part of the two functions treating time steps is the following sequence of operations on the present mesh:

With this general structure of the code, we only have to define what case we want to solve. For the present program, we have chosen to simulate the quasistatic deformation of a vertical cylinder for which the bottom boundary is fixed and the top boundary is pushed down at a prescribed vertical velocity. However, the horizontal velocity of the top boundary is left unspecified -- one can imagine this situation as a well-greased plate pushing from the top onto the cylinder, the points on the top boundary of the cylinder being allowed to slide horizontally along the surface of the plate, but forced to move downward by the plate. The inner and outer boundaries of the cylinder are free and not subject to any prescribed deflection or traction. In addition, gravity acts on the body.

The program text will reveal more about how to implement this situation, and the results section will show what displacement pattern comes out of this simulation.

Possible directions for extensions

The program as is does not really solve an equation that has many applications in practice: quasi-static material deformation based on a purely elastic law is almost boring. However, the program may serve as the starting point for more interesting experiments, and that indeed was the initial motivation for writing it. Here are some suggestions of what the program is missing and in what direction it may be extended:

Plasticity models

The most obvious extension is to use a more

realistic material model for large-scale quasistatic deformation. The natural choice for this would be plasticity, in which a nonlinear relationship between stress and strain replaces equation [stress-strain]. Plasticity models are usually rather complicated to program since the stress-strain dependence is generally non-smooth. The material can be thought of being able to withstand only a maximal stress (the yield stress) after which it starts to “flow”. A mathematical description to this can be given in the form of a variational inequality, which alternatively can be treated as minimizing the elastic energy

\[ E(\mathbf{u}) = (\varepsilon(\mathbf{u}), C\varepsilon(\mathbf{u}))_{\Omega} - (\mathbf{f}, \mathbf{u})_{\Omega} - (\mathbf{b}, \mathbf{u})_{\Gamma_N}, \]

subject to the constraint

\[ f(\sigma(\mathbf{u})) \le 0 \]

on the stress. This extension makes the problem to be solved in each time step nonlinear, so we need another loop within each time step.

Without going into further details of this model, we refer to the excellent book by Simo and Hughes on “Computational Inelasticity” for a comprehensive overview of computational strategies for solving plastic models. Alternatively, a brief but concise description of an algorithm for plasticity is given in an article by S. Commend, A. Truty, and Th. Zimmermann, titled “Stabilized finite elements applied to elastoplasticity: I. Mixed displacement-pressure formulation” (Computer Methods in Applied Mechanics and Engineering, vol. 193, pp. 3559-3586, 2004).

Stabilization issues

The formulation we have chosen, i.e. using

piecewise (bi-, tri-)linear elements for all components of the displacement vector, and treating the stress as a variable dependent on the displacement is appropriate for most materials. However, this so-called displacement-based formulation becomes unstable and exhibits spurious modes for incompressible or nearly-incompressible materials. While fluids are usually not elastic (in most cases, the stress depends on velocity gradients, not displacement gradients, although there are exceptions such as electro-rheologic fluids), there are a few solids that are nearly incompressible, for example rubber. Another case is that many plasticity models ultimately let the material become incompressible, although this is outside the scope of the present program.

Incompressibility is characterized by Poisson's ratio

\[ \nu = \frac{\lambda}{2(\lambda+\mu)}, \]

where $\lambda,\mu$ are the Lam\'e constants of the material. Physical constraints indicate that $-1\le \nu\le \frac 12$ (the condition also follows from mathematical stability considerations). If $\nu$ approaches $\frac 12$, then the material becomes incompressible. In that case, pure displacement-based formulations are no longer appropriate for the solution of such problems, and stabilization techniques have to be employed for a stable and accurate solution. The book and paper cited above give indications as to how to do this, but there is also a large volume of literature on this subject; a good start to get an overview of the topic can be found in the references of the paper by H.-Y. Duan and Q. Lin on “Mixed finite elements of least-squares type for elasticity” (Computer Methods in Applied Mechanics and Engineering, vol. 194, pp. 1093-1112, 2005).

Refinement during timesteps

In the present form, the program

only refines the initial mesh a number of times, but then never again. For any kind of realistic simulation, one would want to extend this so that the mesh is refined and coarsened every few time steps instead. This is not hard to do, in fact, but has been left for future tutorial programs or as an exercise, if you wish. The main complication one has to overcome is that one has to transfer the data that is stored in the quadrature points of the cells of the old mesh to the new mesh, preferably by some sort of projection scheme. This is only slightly messy in the sequential case; in fact, the functions FETools::get_projection_from_quadrature_points_matrix will do the projection, and the FiniteElement::get_restriction_matrix and FiniteElement::get_prolongation_matrix functions will do the transfer between mother and child cells. However, it becomes complicated once we run the program in parallel, since then each process only stores this data for the cells it owned on the old mesh, and it may need to know the 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 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.

Ensuring mesh regularity

At present, the program makes no attempt

to make sure that a cell, after moving its vertices at the end of the time step, still has a valid geometry (i.e. that its Jacobian determinant is positive and bounded away from zero everywhere). It is, in fact, not very hard to set boundary values and forcing terms in such a way that one gets distorted and inverted cells rather quickly. Certainly, in some cases of large deformation, this is unavoidable with a mesh of finite mesh size, but in some other cases this should be preventable by appropriate mesh refinement and/or a reduction of the time step size. The program does not do that, but a more sophisticated version definitely should employ some sort of heuristic defining what amount of deformation of cells is acceptable, and what isn't.

Compiling the program

Finally, just to remind everyone: the program runs in 3d (see the definition of the elastic_problem variable in main(), unlike almost all of the other example programs. While the compiler doesn't care what dimension it compiles for, the linker has to know which library to link with. And as explained in other places, this requires slight changes to the Makefile compared to the other tutorial programs. In particular, everywhere where the 2d versions of libraries are mentioned, one needs to change this to 3d, although this is already done in the distributed version of the Makefile. Conversely, if you want to run the program in 2d (after making the necessary changes to accommodate for a 2d geometry), you have to change the Makefile back to allow for 2d.

The commented program

First the usual list of header files that have already been used in previous example programs:

 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <base/logstream.h>
 #include <base/conditional_ostream.h>
 #include <base/utilities.h>
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/petsc_vector.h>
 #include <lac/petsc_parallel_vector.h>
 #include <lac/petsc_parallel_sparse_matrix.h>
 #include <lac/petsc_solver.h>
 #include <lac/petsc_precondition.h>
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <grid/grid_refinement.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <grid/tria_boundary_lib.h>
 #include <grid/grid_tools.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <dofs/dof_constraints.h>
 #include <dofs/dof_renumbering.h>
 #include <fe/fe_values.h>
 #include <fe/fe_system.h>
 #include <fe/fe_q.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
 #include <numerics/error_estimator.h>

And here the only two new things among the header files: an include file in which symmetric tensors of rank 2 and 4 are implemented, as introduced in the introduction:

 #include <base/symmetric_tensor.h>

And a header that implements filters for iterators looping over all cells. We will use this when selecting only those cells for output that are owned by the present process in a parallel program:

 #include <grid/filtered_iterator.h>

This is then simply C++ again:

 #include <fstream>
 #include <iostream>
 #include <sstream>
 #include <iomanip>

The last step is as in all previous programs:

 using namespace dealii;

So much for the header files. As a matter of good practice, I have started to put everything that corresponds to a certain project into a namespace of its own, named after the problem that we are solving:

 namespace QuasiStaticElasticity
 {

The PointHistory class

As was mentioned in the introduction, we have to store the old stress in quadrature point so that we can compute the residual forces at this point during the next time step. This alone would not warrant a structure with only one member, but in more complicated applications, we would have to store more information in quadrature points as well, such as the history variables of plasticity, etc. In essence, we have to store everything that affects the present state of the material here, which in plasticity is determined by the deformation history variables.

We will not give this class any meaningful functionality beyond being able to store data, i.e. there are no constructors, destructors, or other member functions. In such cases of `dumb' classes, we usually opt to declare them as struct rather than class, to indicate that they are closer to C-style structures than C++-style classes.

   template <int dim>
   struct PointHistory
   {
       SymmetricTensor<2,dim> old_stress;
   };

The stress-strain tensor

Next, we define the linear relationship between the stress and the strain in elasticity. It is given by a tensor of rank 4 that is usually written in the form $C_{ijkl} = \mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}) + \lambda \delta_{ij} \delta_{kl}$. This tensor maps symmetric tensor of rank 2 to symmetric tensors of rank 2. A function implementing its creation for given values of the Lame constants lambda and mu is straightforward:

   template <int dim>
   SymmetricTensor<4,dim>
   get_stress_strain_tensor (const double lambda, const double mu)
   {
     SymmetricTensor<4,dim> tmp;
     for (unsigned int i=0; i<dim; ++i)
       for (unsigned int j=0; j<dim; ++j)
         for (unsigned int k=0; k<dim; ++k)
           for (unsigned int l=0; l<dim; ++l)
             tmp[i][j][k][l] = (((i==k) && (j==l) ? mu : 0.0) +
                                ((i==l) && (j==k) ? mu : 0.0) +
                                ((i==j) && (k==l) ? lambda : 0.0));
     return tmp;
   }

With this function, we will define a static member variable of the main class below that will be used throughout the program as the stress-strain tensor. Note that in more elaborate programs, this will probably be a member variable of some class instead, or a function that returns the stress-strain relationship depending on other input. For example in damage theory models, the Lame constants are considered a function of the prior stress/strain history of a point. Conversely, in plasticity the form of the stress-strain tensor is modified if the material has reached the yield stress in a certain point, and possibly also depending on its prior history.

In the present program, however, we assume that the material is completely elastic and linear, and a constant stress-strain tensor is sufficient for our present purposes.

Auxiliary functions

Before the rest of the program, here are a few functions that we need as tools. These are small functions that are called in inner loops, so we mark them as inline.

The first one computes the symmetric strain tensor for shape function shape_func at quadrature point q_point by forming the symmetric gradient of this shape function. We need that when we want to form the matrix, for example.

We should note that in previous examples where we have treated vector-valued problems, we have always asked the finite element object in which of the vector component the shape function is actually non-zero, and thereby avoided to compute any terms that we could prove were zero anyway. For this, we used the fe.system_to_component_index function that returns in which component a shape function was zero, and also that the fe_values.shape_value and fe_values.shape_grad functions only returned the value and gradient of the single non-zero component of a shape function if this is a vector-valued element.

This was an optimization, and if it isn't terribly time critical, we can get away with a simpler technique: just ask the fe_values for the value or gradient of a given component of a given shape function at a given quadrature point. This is what the fe_values.shape_grad_component(shape_func,q_point,i) call does: return the full gradient of the ith component of shape function shape_func at quadrature point q_point. If a certain component of a certain shape function is always zero, then this will simply always return zero.

As mentioned, using fe_values.shape_grad_component instead of the combination of fe.system_to_component_index and fe_values.shape_grad may be less efficient, but its implementation is optimized for such cases and shouldn't be a big slowdown. We demonstrate the technique here since it is so much simpler and straightforward.

   template <int dim>
   inline
   SymmetricTensor<2,dim>
   get_strain (const FEValues<dim> &fe_values,
               const unsigned int   shape_func,
               const unsigned int   q_point)
   {

Declare a temporary that will hold the return value:

     SymmetricTensor<2,dim> tmp;

First, fill diagonal terms which are simply the derivatives in direction i of the i component of the vector-valued shape function:

     for (unsigned int i=0; i<dim; ++i)
       tmp[i][i] = fe_values.shape_grad_component (shape_func,q_point,i)[i];

Then fill the rest of the strain tensor. Note that since the tensor is symmetric, we only have to compute one half (here: the upper right corner) of the off-diagonal elements, and the implementation of the SymmetricTensor class makes sure that at least to the outside the symmetric entries are also filled (in practice, the class of course stores only one copy). Here, we have picked the upper right half of the tensor, but the lower left one would have been just as good:

     for (unsigned int i=0; i<dim; ++i)
       for (unsigned int j=i+1; j<dim; ++j)
         tmp[i][j]
           = (fe_values.shape_grad_component (shape_func,q_point,i)[j] +
              fe_values.shape_grad_component (shape_func,q_point,j)[i]) / 2;
   
     return tmp;
   }

The second function does something very similar (and therefore is given the same name): compute the symmetric strain tensor from the gradient of a vector-valued field. If you already have a solution field, the fe_values.get_function_grads function allows you to extract the gradients of each component of your solution field at a quadrature point. It returns this as a vector of rank-1 tensors: one rank-1 tensor (gradient) per vector component of the solution. From this we have to reconstruct the (symmetric) strain tensor by transforming the data storage format and symmetrization. We do this in the same way as above, i.e. we avoid a few computations by filling first the diagonal and then only one half of the symmetric tensor (the SymmetricTensor class makes sure that it is sufficient to write only one of the two symmetric components).

Before we do this, though, we make sure that the input has the kind of structure we expect: that is that there are dim vector components, i.e. one displacement component for each coordinate direction. We test this with the Assert macro that will simply abort our program if the condition is not met.

   template <int dim>
   inline
   SymmetricTensor<2,dim>
   get_strain (const std::vector<Tensor<1,dim> > &grad)
   {
     Assert (grad.size() == dim, ExcInternalError());
 
     SymmetricTensor<2,dim> strain;
     for (unsigned int i=0; i<dim; ++i)
       strain[i][i] = grad[i][i];
     
     for (unsigned int i=0; i<dim; ++i)
       for (unsigned int j=i+1; j<dim; ++j)
         strain[i][j] = (grad[i][j] + grad[j][i]) / 2;
     
     return strain;
   }

Finally, below we will need a function that computes the rotation matrix induced by a displacement at a given point. In fact, of course, the displacement at a single point only has a direction and a magnitude, it is the change in direction and magnitude that induces rotations. In effect, the rotation matrix can be computed from the gradients of a displacement, or, more specifically, from the curl.

The formulas by which the rotation matrices are determined are a little awkward, especially in 3d. For 2d, there is a simpler way, so we implement this function twice, once for 2d and once for 3d, so that we can compile and use the program in both space dimensions if so desired -- after all, deal.II is all about dimension independent programming and reuse of algorithm thoroughly tested with cheap computations in 2d, for the more expensive computations in 3d. Here is one case, where we have to implement different algorithms for 2d and 3d, but then can write the rest of the program in a way that is independent of the space dimension.

So, without further ado to the 2d implementation:

   Tensor<2,2>
   get_rotation_matrix (const std::vector<Tensor<1,2> > &grad_u)
   {

First, compute the curl of the velocity field from the gradients. Note that we are in 2d, so the rotation is a scalar:

     const double curl = (grad_u[1][0] - grad_u[0][1]);

From this, compute the angle of rotation:

     const double angle = std::atan (curl);

And from this, build the antisymmetric rotation matrix:

     const double t[2][2] = {{ cos(angle), sin(angle) },
                             {-sin(angle), cos(angle) }};
     return Tensor<2,2>(t);
   }

The 3d case is a little more contrived:

   Tensor<2,3>
   get_rotation_matrix (const std::vector<Tensor<1,3> > &grad_u)
   {

Again first compute the curl of the velocity field. This time, it is a real vector:

     const Point<3> curl (grad_u[2][1] - grad_u[1][2],
                          grad_u[0][2] - grad_u[2][0],
                          grad_u[1][0] - grad_u[0][1]);

From this vector, using its magnitude, compute the tangent of the angle of rotation, and from it the actual angle:

     const double tan_angle = std::sqrt(curl*curl);
     const double angle = std::atan (tan_angle);

Now, here's one problem: if the angle of rotation is too small, that means that there is no rotation going on (for example a translational motion). In that case, the rotation matrix is the identity matrix.

The reason why we stress that is that in this case we have that tan_angle==0. Further down, we need to divide by that number in the computation of the axis of rotation, and we would get into trouble when dividing doing so. Therefore, let's shortcut this and simply return the identity matrix if the angle of rotation is really small:

     if (angle < 1e-9)
       {
         static const double rotation[3][3]
           = {{ 1, 0, 0}, { 0, 1, 0 }, { 0, 0, 1 } };
         static const Tensor<2,3> rot(rotation);
         return rot;
       }

Otherwise compute the real rotation matrix. The algorithm for this is not exactly obvious, but can be found in a number of books, particularly on computer games where rotation is a very frequent operation. Online, you can find a description at http://www.makegames.com/3drotation/ and (this particular form, with the signs as here) at http://www.gamedev.net/reference/articles/article1199.asp:

     const double c = std::cos(angle);
     const double s = std::sin(angle);
     const double t = 1-c;
 
     const Point<3> axis = curl/tan_angle;
     const double rotation[3][3]
       = {{ t*axis[0]*axis[0]+c,
            t*axis[0]*axis[1]+s*axis[2],
            t*axis[0]*axis[2]-s*axis[1]},
          { t*axis[0]*axis[1]-s*axis[2],
            t*axis[1]*axis[1]+c,
            t*axis[1]*axis[2]+s*axis[0]},
          { t*axis[0]*axis[2]+s*axis[1],
            t*axis[1]*axis[1]-s*axis[0],
            t*axis[2]*axis[2]+c  } };
     return Tensor<2,3>(rotation);
   }

The TopLevel class

This is the main class of the program. Since the namespace already indicates what problem we are solving, let's call it by what it does: it directs the flow of the program, i.e. it is the toplevel driver.

The member variables of this class are essentially as before, i.e. it has to have a triangulation, a DoF handler and associated objects such as constraints, variables that describe the linear system, etc. There are a good number of more member functions now, which we will explain below.

The external interface of the class, however, is unchanged: it has a public constructor and desctructor, and it has a run function that initiated all the work.

   template <int dim>
   class TopLevel 
   {
     public:
       TopLevel ();
       ~TopLevel ();
       void run ();
     
     private:

The private interface is more extensive than in step-17. First, we obviously need functions that create the initial mesh, set up the variables that describe the linear system on the present mesh (i.e. matrices and vectors), and then functions that actually assemble the system, direct what has to be solved in each time step, a function that solves the linear system that arises in each timestep (and returns the number of iterations it took), and finally output the solution vector on the currect mesh:

       void create_coarse_grid ();
     
       void setup_system ();
       
       void assemble_system ();
       
       void solve_timestep ();
 
       unsigned int solve_linear_problem ();
 
       void output_results () const;

All, except for the first two, of these functions are called in each timestep. Since the first time step is a little special, we have separate functions that describe what has to happen in a timestep: one for the first, and one for all following timesteps:

       void do_initial_timestep ();
 
       void do_timestep ();

Then we need a whole bunch of functions that do various things. The first one refines the initial grid: we start on the coarse grid with a pristine state, solve the problem, then look at it and refine the mesh accordingly, and start the same process over again, again with a pristine state. Thus, refining the initial mesh is somewhat simpler than refining a grid between two successive time steps, since it does not involve transferring data from the old to the new triangulation, in particular the history data that is stored in each quadrature point.

       void refine_initial_grid ();

At the end of each time step, we want to move the mesh vertices around according to the incremental displacement computed in this time step. This is the function in which this is done:

       void move_mesh ();

Next are two functions that handle the history variables stored in each quadrature point. The first one is called before the first timestep to set up a pristine state for the history variables. It only works on those quadrature points on cells that belong to the present processor:

       void setup_quadrature_point_history ();

The second one updates the history variables at the end of each timestep:

       void update_quadrature_point_history ();

After the member functions, here are the member variables. The first ones have all been discussed in more detail in previous example programs:

       Triangulation<dim>   triangulation;
 
       FESystem<dim>        fe;
 
       DoFHandler<dim>      dof_handler;
 
       ConstraintMatrix     hanging_node_constraints;

One difference of this program is that we declare the quadrature formula in the class declaration. The reason is that in all the other programs, it didn't do much harm if we had used different quadrature formulas when computing the matrix and the righ hand side, for example. However, in the present case it does: we store information in the quadrature points, so we have to make sure all parts of the program agree on where they are and how many there are on each cell. Thus, let us first declare the quadrature formula that will be used throughout...

       const QGauss<dim>          quadrature_formula;

... and then also have a vector of history objects, one per quadrature point on those cells for which we are responsible (i.e. we don't store history data for quadrature points on cells that are owned by other processors).

       std::vector<PointHistory<dim> > quadrature_point_history;

The way this object is accessed is through a user pointer that each cell, face, or edge holds: it is a void* pointer that can be used by application programs to associate arbitrary data to cells, faces, or edges. What the program actually does with this data is within its own responsibility, the library just allocates some space for these pointers, and application programs can set and read the pointers for each of these objects.

Further: we need the objects of linear systems to be solved, i.e. matrix, right hand side vector, and the solution vector. Since we anticipate solving big problems, we use the same types as in step-17, i.e. distributed parallel matrices and vectors built on top of the PETSc library. Conveniently, they can also be used when running on only a single machine, in which case this machine happens to be the only one in our parallel universe.

However, as a difference to step-17, we do not store the solution vector -- which here is the incremental displacements computed in each time step -- in a distributed fashion. I.e., of course it must be a distributed vector when computing it, but immediately after that we make sure each processor has a complete copy. The reason is that we had already seen in step-17 that many functions needed a complete copy. While it is not hard to get it, this requires communication on the network, and is thus slow. In addition, these were repeatedly the same operations, which is certainly undesirable unless the gains of not always having to store the entire vector outweighs it. When writing this program, it turned out that we need a complete copy of the solution in so many places that it did not seem worthwhile to only get it when necessary. Instead, we opted to obtain the complete copy once and for all, and instead get rid of the distributed copy immediately. Thus, note that the declaration of inremental_displacement does not denote a distribute vector as would be indicated by the middle namespace MPI:

       PETScWrappers::MPI::SparseMatrix system_matrix;
 
       PETScWrappers::MPI::Vector       system_rhs;
 
       PETScWrappers::Vector            incremental_displacement;

The next block of variables is then related to the time dependent nature of the problem: they denote the length of the time interval which we want to simulate, the present time and number of time step, and length of present timestep:

       double       present_time;
       double       present_timestep;
       double       end_time;
       unsigned int timestep_no;

Then a few variables that have to do with parallel processing: first, a variable denoting the MPI communicator we use, and then two numbers telling us how many participating processors there are, and where in this world we are. Finally, a stream object that makes sure only one processor is actually generating output to the console. This is all the same as in step-17:

       MPI_Comm mpi_communicator;
 
       const unsigned int n_mpi_processes;
 
       const unsigned int this_mpi_process;
 
       ConditionalOStream pcout;

Here is a vector where each entry denotes the numbers of degrees of freedom that are stored on the processor with that particular number:

       std::vector<unsigned int> local_dofs_per_process;

Next, how many degrees of freedom the present processor stores. This is, of course, an abbreviation to local_dofs_per_process[this_mpi_process].

       unsigned int         n_local_dofs;

In the same direction, also cache how many cells the present processor owns. Note that the cells that belong to a processor are not necessarily contiguously numbered (when iterating over them using active_cell_iterator).

       unsigned int         n_local_cells;

Finally, we have a static variable that denotes the linear relationship between the stress and strain. Since it is a constant object that does not depend on any input (at least not in this program), we make it a static variable and will initialize it in the same place where we define the constructor of this class:

       static const SymmetricTensor<4,dim> stress_strain_tensor;
   };

The BodyForce class

Before we go on to the main functionality of this program, we have to define what forces will act on the body whose deformation we want to study. These may either be body forces or boundary forces. Body forces are generally mediated by one of the four basic physical types of forces: gravity, strong and weak interaction, and electromagnetism. Unless one wants to consider subatomic objects (for which quasistatic deformation is irrelevant and an inappropriate description anyway), only gravity and electromagnetic forces need to be considered. Let us, for simplicity assume that our body has a certain mass density, but is either non-magnetic and not electrically conducting or that there are no significant electromagnetic fields around. In that case, the body forces are simply rho g, where rho is the material density and g is a vector in negative z-direction with magnitude 9.81 m/s^2. Both the density and g are defined in the function, and we take as the density 7700 kg/m^3, a value commonly assumed for steel.

To be a little more general and to be able to do computations in 2d as well, we realize that the body force is always a function returning a dim dimensional vector. We assume that gravity acts along the negative direction of the last, i.e. dim-1th coordinate. The rest of the implementation of this function should be mostly self-explanatory given similar definitions in previous example programs. Note that the body force is independent of the location; to avoid compiler warnings about unused function arguments, we therefore comment out the name of the first argument of the vector_value function:

   template <int dim>
   class BodyForce :  public Function<dim> 
   {
     public:
       BodyForce ();
     
       virtual
       void
       vector_value (const Point<dim> &p,
                     Vector<double>   &values) const;
 
       virtual
       void
       vector_value_list (const std::vector<Point<dim> > &points,
                          std::vector<Vector<double> >   &value_list) const;
   };
 
 
   template <int dim>
   BodyForce<dim>::BodyForce ()
                   :
                   Function<dim> (dim)
   {}
 
 
   template <int dim>
   inline
   void
   BodyForce<dim>::vector_value (const Point<dim> &/ *p* /,
                                 Vector<double>   &values) const 
   {
     Assert (values.size() == dim, 
             ExcDimensionMismatch (values.size(), dim));
 
     const double g   = 9.81;
     const double rho = 7700;
     
     values = 0;
     values(dim-1) = -rho * g;
   }
 
 
 
   template <int dim>
   void
   BodyForce<dim>::vector_value_list (const std::vector<Point<dim> > &points,
                                      std::vector<Vector<double> >   &value_list) const 
   {
     const unsigned int n_points = points.size();
 
     Assert (value_list.size() == n_points, 
             ExcDimensionMismatch (value_list.size(), n_points));
 
     for (unsigned int p=0; p<n_points; ++p)
       BodyForce<dim>::vector_value (points[p],
                                     value_list[p]);
   }

The IncrementalBoundaryValue class

In addition to body forces, movement can be induced by boundary forces and forced boundary displacement. The latter case is equivalent to forces being chosen in such a way that they induce certain displacement.

For quasistatic displacement, typical boundary forces would be pressure on a body, or tangential friction against another body. We chose a somewhat simpler case here: we prescribe a certain movement of (parts of) the boundary, or at least of certain components of the displacement vector. We describe this by another vector-valued function that, for a given point on the boundary, returns the prescribed displacement.

Since we have a time-dependent problem, the displacement increment of the boundary equals the displacement accumulated during the length of the timestep. The class therefore has to know both the present time and the length of the present time step, and can then approximate the incremental displacement as the present velocity times the present timestep.

For the purposes of this program, we choose a simple form of boundary displacement: we displace the top boundary with constant velocity downwards. The rest of the boundary is either going to be fixed (and is then described using an object of type ZeroFunction) or free (Neumann-type, in which case nothing special has to be done). The implementation of the class describing the constant downward motion should then be obvious using the knowledge we gained through all the previous example programs:

   template <int dim>
   class IncrementalBoundaryValues :  public Function<dim> 
   {
     public:
       IncrementalBoundaryValues (const double present_time,
                                  const double present_timestep);
     
       virtual
       void
       vector_value (const Point<dim> &p,
                     Vector<double>   &values) const;
 
       virtual
       void
       vector_value_list (const std::vector<Point<dim> > &points,
                          std::vector<Vector<double> >   &value_list) const;
 
     private:
       const double velocity;
       const double present_time;
       const double present_timestep;
   };
 
 
   template <int dim>
   IncrementalBoundaryValues<dim>::
   IncrementalBoundaryValues (const double present_time,
                              const double present_timestep)
                   :
                   Function<dim> (dim),
                   velocity (.1),
                   present_time (present_time),
                   present_timestep (present_timestep)
   {}
 
 
   template <int dim>
   void
   IncrementalBoundaryValues<dim>::
   vector_value (const Point<dim> &/ *p* /,
                 Vector<double>   &values) const 
   {
     Assert (values.size() == dim, 
             ExcDimensionMismatch (values.size(), dim));
 
     values = 0;
     values(2) = -present_timestep * velocity;
   }
 
 
 
   template <int dim>
   void
   IncrementalBoundaryValues<dim>::
   vector_value_list (const std::vector<Point<dim> > &points,
                      std::vector<Vector<double> >   &value_list) const 
   {
     const unsigned int n_points = points.size();
 
     Assert (value_list.size() == n_points, 
             ExcDimensionMismatch (value_list.size(), n_points));
 
     for (unsigned int p=0; p<n_points; ++p)
       IncrementalBoundaryValues<dim>::vector_value (points[p],
                                                     value_list[p]);
   }

Implementation of the TopLevel class

Now for the implementation of the main class. First, we initialize the stress-strain tensor, which we have declared as a static const variable. We chose Lame constants that are appropriate for steel:

   template <int dim>
   const SymmetricTensor<4,dim>
   TopLevel<dim>::stress_strain_tensor
   = get_stress_strain_tensor<dim> (/ *lambda = * / 9.695e10,
                                    / *mu     = * / 7.617e10);

The public interface

The next step is the definition of constructors and descructors. There are no surprises here: we choose linear and continuous finite elements for each of the dim vector components of the solution, and a Gaussian quadrature formula with 2 points in each coordinate direction. The destructor should be obvious:

   template <int dim>
   TopLevel<dim>::TopLevel ()
                   :
                   fe (FE_Q<dim>(1), dim),
                   dof_handler (triangulation),
                   quadrature_formula (2),
                   mpi_communicator (MPI_COMM_WORLD),
                   n_mpi_processes (Utilities::System::get_n_mpi_processes(mpi_communicator)),
                   this_mpi_process (Utilities::System::get_this_mpi_process(mpi_communicator)),
                   pcout (std::cout, this_mpi_process == 0)
   {}
 
 
 
   template <int dim>
   TopLevel<dim>::~TopLevel () 
   {
     dof_handler.clear ();
   }

The last of the public functions is the one that directs all the work, run(). It initializes the variables that describe where in time we presently are, then runs the first time step, then loops over all the other time steps. Note that for simplicity we use a fixed time step, whereas a more sophisticated program would of course have to choose it in some more reasonable way adaptively:

   template <int dim>
   void TopLevel<dim>::run () 
   {
     present_time = 0;
     present_timestep = 1;
     end_time = 10;
     timestep_no = 0;
   
     do_initial_timestep ();
 
     while (present_time < end_time)
       do_timestep ();
   }

TopLevel::create_coarse_grid

The next function in the order in which they were declared above is the one that creates the coarse grid from which we start. For this example program, we want to compute the deformation of a cylinder under axial compression. The first step therefore is to generate a mesh for a cylinder of length 3 and with inner and outer radii of 0.8 and 1, respectively. Fortunately, there is a library function for such a mesh.

In a second step, we have to associated boundary conditions with the upper and lower faces of the cylinder. We choose a boundary indicator of 0 for the boundary faces that are characterized by their midpoints having z-coordinates of either 0 (bottom face), an indicator of 1 for z=3 (top face); finally, we use boundary indicator 2 for all faces on the inside of the cylinder shell, and 3 for the outside.

   template <int dim>
   void TopLevel<dim>::create_coarse_grid ()
   {
     const double inner_radius = 0.8,
                  outer_radius = 1;
     GridGenerator::cylinder_shell (triangulation,
                                    3, inner_radius, outer_radius);
     for (typename Triangulation<dim>::active_cell_iterator
            cell=triangulation.begin_active();
          cell!=triangulation.end(); ++cell)
       for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
         if (cell->face(f)->at_boundary())
           {
             const Point<dim> face_center = cell->face(f)->center();
             
             if (face_center[2] == 0)
               cell->face(f)->set_boundary_indicator (0);
             else if (face_center[2] == 3)
               cell->face(f)->set_boundary_indicator (1);
             else if (std::sqrt(face_center[0]*face_center[0] +
                                face_center[1]*face_center[1])
                      <
                      (inner_radius + outer_radius) / 2)
               cell->face(f)->set_boundary_indicator (2);
             else
               cell->face(f)->set_boundary_indicator (3);
           }

In order to make sure that new vertices are placed correctly on mesh refinement, we have to associate objects describing those parts of the boundary that do not consist of straight parts. Corresponding to the cylinder shell generator function used above, there are classes that can be used to describe the geometry of cylinders. We need to use different objects for the inner and outer parts of the cylinder, with different radii; the second argument to the constructor indicates the axis around which the cylinder revolves -- in this case the z-axis. Note that the boundary objects need to live as long as the triangulation does; we can achieve this by making the objects static, which means that they live as long as the program runs:

     static const CylinderBoundary<dim> inner_cylinder (inner_radius, 2);
     static const CylinderBoundary<dim> outer_cylinder (outer_radius, 2);

We then attach these two objects to the triangulation, and make them correspond to boundary indicators 2 and 3:

     triangulation.set_boundary (2, inner_cylinder);
     triangulation.set_boundary (3, outer_cylinder);

There's one more thing we have to take care of (we should have done so above already, but for didactic reasons it was more appropriate to handle it after discussing boundary objects). Boundary indicators in deal.II, for mostly historic reasons, serve a dual purpose: they describe the type of a boundary for other places in a program where different boundary conditions are implemented; and they describe which boundary object (as the ones associated above) should be queried when new boundary points need to be placed upon mesh refinement. In the prefix to this function, we have discussed the boundary condition issue, and the boundary geometry issue was mentioned just above. But there is a case where we have to be careful with geometry: what happens if a cell is refined that has two faces with different boundary indicators? For example one at the edges of the cylinder? In that case, the library wouldn't know where to put new points in the middle of edges (one of the twelve lines of a hexahedron). In fact, the library doesn't even care about the boundary indicator of adjacent faces when refining edges: it considers the boundary indicators associated with the edges themselves. So what do we want to happen with the edges of the cylinder shell: they sit on both faces with boundary indicators 2 or 3 (inner or outer shell) and 0 or 1 (for which no boundary objects have been specified, and for which the library therefore assumes straight lines). Obviously, we want these lines to follow the curved shells, so we have to assign all edges along faces with boundary indicators 2 or 3 these same boundary indicators to make sure they are refined using the appropriate geometry objects. This is easily done:

     for (typename Triangulation<dim>::active_face_iterator
            face=triangulation.begin_active_face();
          face!=triangulation.end_face(); ++face)
       if (face->at_boundary())
         if ((face->boundary_indicator() == 2)
             ||
             (face->boundary_indicator() == 3))
           for (unsigned int edge = 0; edge<GeometryInfo<dim>::lines_per_face;
                ++edge)
             face->line(edge)
               ->set_boundary_indicator (face->boundary_indicator());

Once all this is done, we can refine the mesh once globally:

     triangulation.refine_global (1);

As the final step, we need to set up a clean state of the data that we store in the quadrature points on all cells that are treated on the present processor. To do so, we also have to know which processors are ours in the first place. This is done in the following two function calls:

     GridTools::partition_triangulation (n_mpi_processes, triangulation);
     setup_quadrature_point_history ();  
   }

TopLevel::setup_system

The next function is the one that sets up the data structures for a given mesh. This is done in most the same way as in step-17: distribute the degrees of freedom, then sort these degrees of freedom in such a way that each processor gets a contiguous chunk of them. Note that subdivions into chunks for each processor is handled in the functions that create or refine grids, unlike in the previous example program (the point where this happens is mostly a matter of taste; here, we chose to do it when grids are created since in the do_initial_timestep and do_timestep functions we want to output the number of cells on each processor at a point where we haven't called the present function yet).

   template <int dim>
   void TopLevel<dim>::setup_system ()
   {
     dof_handler.distribute_dofs (fe);
     DoFRenumbering::subdomain_wise (dof_handler);

The next thing is to store some information for later use on how many cells or degrees of freedom the present processor, or any of the processors has to work on. First the cells local to this processor...

     n_local_cells
       = GridTools::count_cells_with_subdomain_association (triangulation,
                                                            this_mpi_process);

...and then a list of numbers of how many degrees of freedom each processor has to handle:

     local_dofs_per_process.resize (n_mpi_processes);
     for (unsigned int i=0; i<n_mpi_processes; ++i)
       local_dofs_per_process[i]
         = DoFTools::count_dofs_with_subdomain_association (dof_handler, i);

Finally, make it easier to denote how many degrees of freedom the present process has to deal with, by introducing an abbreviation:

     n_local_dofs = local_dofs_per_process[this_mpi_process];

The next step is to set up constraints due to hanging nodes. This has been handled many times before:

     hanging_node_constraints.clear ();
     DoFTools::make_hanging_node_constraints (dof_handler,
                                              hanging_node_constraints);
     hanging_node_constraints.close ();

And then we have to set up the matrix. Here we deviate from step-17, in which we simply used PETSc's ability to just know about the size of the matrix and later allocate those nonzero elements that are being written to. While this works just fine from a correctness viewpoint, it is not at all efficient: if we don't give PETSc a clue as to which elements are written to, it is (at least at the time of this writing) unbearably slow when we set the elements in the matrix for the first time (i.e. in the first timestep). Later on, when the elements have been allocated, everything is much faster. In experiments we made, the first timestep can be accelerated by almost two orders of magnitude if we instruct PETSc which elements will be used and which are not.

To do so, we first generate the sparsity pattern of the matrix we are going to work with, and make sure that the condensation of hanging node constraints add the necessary additional entries in the sparsity pattern:

     CompressedSparsityPattern sparsity_pattern (dof_handler.n_dofs(),
                                                 dof_handler.n_dofs());
     DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
     hanging_node_constraints.condense (sparsity_pattern);

Note that we have used the CompressedSparsityPattern class here that was already introduced in step-11, rather than the SparsityPattern class that we have used in all other cases. The reason for this is that for the latter class to work we have to give an initial upper bound for the number of entries in each row, a task that is traditionally done by DoFHandler::max_couplings_between_dofs(). However, this function suffers from a serious problem: it has to compute an upper bound to the number of nonzero entries in each row, and this is a rather complicated task, in particular in 3d. In effect, while it is quite accurate in 2d, it often comes up with much too large a number in 3d, and in that case the SparsityPattern allocates much too much memory at first, often several 100 MBs. This is later corrected when DoFTools::make_sparsity_pattern is called and we realize that we don't need all that much memory, but at time it is already too late: for large problems, the temporary allocation of too much memory can lead to out-of-memory situations.

In order to avoid this, we resort to the CompressedSparsityPattern class that is slower but does not require any up-front estimate on the number of nonzero entries per row. It therefore only ever allocates as much memory as it needs at any given time, and we can build it even for large 3d problems.

It is also worth noting that the sparsity pattern we construct is global, i.e. comprises all degrees of freedom whether they will be owned by the processor we are on or another one (in case this program is run in parallel via MPI). This of course is not optimal -- it limits the size of the problems we can solve, since storing the entire sparsity pattern (even if only for a short time) on each processor does not scale well. However, there are several more places in the program in which we do this, for example we always keep the global triangulation and DoF handler objects around, even if we only work on part of them. At present, deal.II does not have the necessary facilities to completely distribute these objects (a task that, indeed, is very hard to achieve with adaptive meshes, since well-balanced subdivisions of a domain tend to become unbalanced as the mesh is adaptively refined).

With this data structure, we can then go to the PETSc sparse matrix and tell it to pre-allocate all the entries we will later want to write to:

     system_matrix.reinit (mpi_communicator,
                           sparsity_pattern,
                           local_dofs_per_process,
                           local_dofs_per_process,
                           this_mpi_process);

After this point, no further explicit knowledge of the sparsity pattern is required any more and we can let the sparsity_pattern variable go out of scope without any problem.

The last task in this function is then only to reset the right hand side vector as well as the solution vector to its correct size; remember that the solution vector is a local one, unlike the right hand side that is a distributed parallel one and therefore needs to know the MPI communicator over which it is supposed to transmit messages:

     system_rhs.reinit (mpi_communicator, dof_handler.n_dofs(), n_local_dofs);
     incremental_displacement.reinit (dof_handler.n_dofs());
   }

TopLevel::assemble_system

Again, assembling the system matrix and right hand side follows the same structure as in many example programs before. In particular, it is mostly equivalent to step-17, except for the different right hand side that now only has to take into account internal stresses. In addition, assembling the matrix is made significantly more transparent by using the SymmetricTensor class: note the elegance of forming the scalar products of symmetric tensors of rank 2 and 4. The implementation is also more general since it is independent of the fact that we may or may not be using an isotropic elasticity tensor.

The first part of the assembly routine is as always:

   template <int dim>
   void TopLevel<dim>::assemble_system () 
   {
     system_rhs = 0;
     system_matrix = 0;
 
     FEValues<dim> fe_values (fe, quadrature_formula, 
                              update_values   | update_gradients |
                              update_quadrature_points | update_JxW_values);
 
     const unsigned int   dofs_per_cell = fe.dofs_per_cell;
     const unsigned int   n_q_points    = quadrature_formula.size();
 
     FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
     Vector<double>       cell_rhs (dofs_per_cell);
 
     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
     BodyForce<dim>      body_force;
     std::vector<Vector<double> > body_force_values (n_q_points,
                                                     Vector<double>(dim));

As in step-17, we only need to loop over all cells that belong to the present processor:

     typename DoFHandler<dim>::active_cell_iterator
       cell = dof_handler.begin_active(),
       endc = dof_handler.end();
     for (; cell!=endc; ++cell)
       if (cell->subdomain_id() == this_mpi_process)
         {
           cell_matrix = 0;
           cell_rhs = 0;
 
           fe_values.reinit (cell);

Then loop over all indices i,j and quadrature points and assemble the system matrix contributions from this cell. Note how we extract the symmetric gradients (strains) of the shape functions at a given quadrature point from the FEValues object, and the elegance with which we form the triple contraction eps_phi_i : C : eps_phi_j; the latter needs to be compared to the clumsy computations needed in step-17, both in the introduction as well as in the respective place in the program:

           for (unsigned int i=0; i<dofs_per_cell; ++i)
             for (unsigned int j=0; j<dofs_per_cell; ++j) 
               for (unsigned int q_point=0; q_point<n_q_points;
                    ++q_point)
                 {
                   const SymmetricTensor<2,dim>
                     eps_phi_i = get_strain (fe_values, i, q_point),
                     eps_phi_j = get_strain (fe_values, j, q_point);
 
                   cell_matrix(i,j) 
                     += (eps_phi_i * stress_strain_tensor * eps_phi_j
                         *
                         fe_values.JxW (q_point));
                 }

Then also assemble the local right hand side contributions. For this, we need to access the prior stress value in this quadrature point. To get it, we use the user pointer of this cell that points into the global array to the quadrature point data corresponding to the first quadrature point of the present cell, and then add an offset corresponding to the index of the quadrature point we presently consider:

           const PointHistory<dim> *local_quadrature_points_data
             = reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());

In addition, we need the values of the external body forces at the quadrature points on this cell:

           body_force.vector_value_list (fe_values.get_quadrature_points(),
                                         body_force_values);

Then we can loop over all degrees of freedom on this cell and compute local contributions to the right hand side:

           for (unsigned int i=0; i<dofs_per_cell; ++i)
             {
               const unsigned int 
                 component_i = fe.system_to_component_index(i).first;
           
               for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
                 {
                   const SymmetricTensor<2,dim> &old_stress
                     = local_quadrature_points_data[q_point].old_stress;
                 
                   cell_rhs(i) += (body_force_values[q_point](component_i) *
                                   fe_values.shape_value (i,q_point)
                                   -
                                   old_stress *
                                   get_strain (fe_values,i,q_point))
                                  *
                                  fe_values.JxW (q_point);
                 }
             }

Now that we have the local contributions to the linear system, we need to transfer it into the global objects. This is done exactly as in step-17:

           cell->get_dof_indices (local_dof_indices);
 
           hanging_node_constraints
             .distribute_local_to_global (cell_matrix,
                                          local_dof_indices,
                                          system_matrix);
 
           hanging_node_constraints
             .distribute_local_to_global (cell_rhs,
                                          local_dof_indices,
                                          system_rhs);
         }

The last step is to again fix up boundary values, just as we already did in previous programs. A slight complication is that the apply_boundary_values function wants to have a solution vector compatible with the matrix and right hand side (i.e. here a distributed parallel vector, rather than the sequential vector we use in this program) in order to preset the entries of the solution vector with the correct boundary values. We provide such a compatible vector in the form of a temporary vector which we then copy into the sequential one.

We make up for this complication by showing how boundary values can be used flexibly: following the way we create the triangulation, there are three distinct boundary indicators used to describe the domain, corresponding to the bottom and top faces, as well as the inner/outer surfaces. We would like to impose boundary conditions of the following type: The inner and outer cylinder surfaces are free of external forces, a fact that corresponds to natural (Neumann-type) boundary conditions for which we don't have to do anything. At the bottom, we want no movement at all, corresponding to the cylinder being clamped or cemented in at this part of the boundary. At the top, however, we want a prescribed vertical downward motion compressing the cylinder; in addition, we only want to restrict the vertical movement, but not the horizontal ones -- one can think of this situation as a well-greased plate sitting on top of the cylinder pushing it downwards: the atoms of the cylinder are forced to move downward, but they are free to slide horizontally along the plate.

The way to describe this is as follows: for boundary indicator zero (bottom face) we use a dim-dimensional zero function representing no motion in any coordinate direction. For the boundary with indicator 1 (top surface), we use the IncrementalBoundaryValues class, but we specify an additional argument to the VectorTools::interpolate_boundary_values function denoting which vector components it should apply to; this is a vector of bools for each vector component and because we only want to restrict vertical motion, it has only its last component set:

     std::vector<bool> z_component (dim, false);
     z_component[dim-1] = true;
     std::map<unsigned int,double> boundary_values;
     VectorTools::
       interpolate_boundary_values (dof_handler,
                                    0,
                                    ZeroFunction<dim> (dim),
                                    boundary_values);
     VectorTools::
       interpolate_boundary_values (dof_handler,
                                    1,
                                    IncrementalBoundaryValues<dim>(present_time,
                                                                   present_timestep),
                                    boundary_values,
                                    z_component);
     
     PETScWrappers::MPI::Vector tmp (mpi_communicator, dof_handler.n_dofs(),
                                     n_local_dofs);
     MatrixTools::apply_boundary_values (boundary_values,
                                         system_matrix, tmp,
                                         system_rhs, false);
     incremental_displacement = tmp;
   }

TopLevel::solve_timestep

The next function is the one that controls what all has to happen within a timestep. The order of things should be relatively self-explanatory from the function names:

   template <int dim>
   void TopLevel<dim>::solve_timestep ()
   {
     pcout << "    Assembling system..." << std::flush;
     assemble_system ();
     pcout << " norm of rhs is " << system_rhs.l2_norm()
           << std::endl;
       
     const unsigned int n_iterations = solve_linear_problem ();
   
     pcout << "    Solver converged in " << n_iterations
           &