 Reference documentation for deal.II version 8.4.1
The step-8 tutorial program

1. Introduction
2. The commented program

# Introduction

In real life, most partial differential equations are really systems of equations. Accordingly, the solutions are usually vector-valued. The deal.II library supports such problems (see the extensive documentation in the Handling vector valued problems module), and we will show that that is mostly rather simple. The only more complicated problems are in assembling matrix and right hand side, but these are easily understood as well.

Note
The material presented here is also discussed in video lecture 19. (All video lectures are also available here.)

In this tutorial program we will want to solve the elastic equations. They are an extension to Laplace's equation with a vector-valued solution that describes the displacement in each space direction of a rigid body which is subject to a force. Of course, the force is also vector-valued, meaning that in each point it has a direction and an absolute value. The elastic equations are the following:

$- \partial_j (c_{ijkl} \partial_k u_l) = f_i, \qquad i=1\ldots d,$

where the values $$c_{ijkl}$$ are the stiffness coefficients and will usually depend on the space coordinates. In many cases, one knows that the material under consideration is isotropic, in which case by introduction of the two coefficients $$\lambda$$ and $$\mu$$ the coefficient tensor reduces to

$c_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}).$

The elastic equations can then be rewritten in much simpler a form:

$- \nabla \lambda (\nabla\cdot {\mathbf u}) - (\nabla \cdot \mu \nabla) {\mathbf u} - \nabla\cdot \mu (\nabla {\mathbf u})^T = {\mathbf f},$

and the respective bilinear form is then

$a({\mathbf u}, {\mathbf v}) = \left( \lambda \nabla\cdot {\mathbf u}, \nabla\cdot {\mathbf v} \right)_\Omega + \sum_{k,l} \left( \mu \partial_k u_l, \partial_k v_l \right)_\Omega + \sum_{k,l} \left( \mu \partial_k u_l, \partial_l v_k \right)_\Omega,$

or also writing the first term a sum over components:

$a({\mathbf u}, {\mathbf v}) = \sum_{k,l} \left( \lambda \partial_l u_l, \partial_k v_k \right)_\Omega + \sum_{k,l} \left( \mu \partial_k u_l, \partial_k v_l \right)_\Omega + \sum_{k,l} \left( \mu \partial_k u_l, \partial_l v_k \right)_\Omega.$

Note
As written, the equations above are generally considered to be the right description for the displacement of three-dimensional objects if the displacement is small and we can assume that Hooke's law is valid. In that case, the indices $$i,j,k,l$$ above all run over the set $$\{1,2,3\}$$ (or, in the C++ source, over $$\{0,1,2\}$$). However, as is, the program runs in 2d, and while the equations above also make mathematical sense in that case, they would only describe a truly two-dimensional solid. In particular, they are not the appropriate description of an $$x-y$$ cross-section of a body infinite in the $$z$$ direction; this is in contrast to many other two-dimensional equations that can be obtained by assuming that the body has infinite extent in $$z$$-direction and that the solution function does not depend on the $$z$$ coordinate. On the other hand, there are equations for two-dimensional models of elasticity; see for example the Wikipedia article on plane strain, antiplane shear and plan stress.

But let's get back to the original problem. How do we assemble the matrix for such an equation? A very long answer with a number of different alternatives is given in the documentation of the Handling vector valued problems module. Historically, the solution shown below was the only one available in the early years of the library. It turns out to also be the fastest. On the other hand, if a few per cent of compute time do not matter, there are simpler and probably more intuitive ways to assemble the linear system than the one discussed below but that weren't available until several years after this tutorial program was first written; if you are interested in them, take a look at the Handling vector valued problems module.

Let us go back to the question of how to assemble the linear system. The first thing we need is some knowledge about how the shape functions work in the case of vector-valued finite elements. Basically, this comes down to the following: let $$n$$ be the number of shape functions for the scalar finite element of which we build the vector element (for example, we will use bilinear functions for each component of the vector-valued finite element, so the scalar finite element is the FE_Q(1) element which we have used in previous examples already, and $$n=4$$ in two space dimensions). Further, let $$N$$ be the number of shape functions for the vector element; in two space dimensions, we need $$n$$ shape functions for each component of the vector, so $$N=2n$$. Then, the $$i$$th shape function of the vector element has the form

$\Phi_i({\mathbf x}) = \varphi_{\text{base}(i)}({\mathbf x})\ {\mathbf e}_{\text{comp}(i)},$

where $$e_l$$ is the $$l$$th unit vector, $$\text{comp}(i)$$ is the function that tells us which component of $$\Phi_i$$ is the one that is nonzero (for each vector shape function, only one component is nonzero, and all others are zero). $$\varphi_{\text{base}(i)}(x)$$ describes the space dependence of the shape function, which is taken to be the $$\text{base}(i)$$-th shape function of the scalar element. Of course, while $$i$$ is in the range $$0,\ldots,N-1$$, the functions $$\text{comp}(i)$$ and $$\text{base}(i)$$ have the ranges $$0,1$$ (in 2D) and $$0,\ldots,n-1$$, respectively.

For example (though this sequence of shape functions is not guaranteed, and you should not rely on it), the following layout could be used by the library:

\begin{eqnarray*} \Phi_0({\mathbf x}) &=& \left(\begin{array}{c} \varphi_0({\mathbf x}) \\ 0 \end{array}\right), \\ \Phi_1({\mathbf x}) &=& \left(\begin{array}{c} 0 \\ \varphi_0({\mathbf x}) \end{array}\right), \\ \Phi_2({\mathbf x}) &=& \left(\begin{array}{c} \varphi_1({\mathbf x}) \\ 0 \end{array}\right), \\ \Phi_3({\mathbf x}) &=& \left(\begin{array}{c} 0 \\ \varphi_1({\mathbf x}) \end{array}\right), \ldots \end{eqnarray*}

where here

$\text{comp}(0)=0, \quad \text{comp}(1)=1, \quad \text{comp}(2)=0, \quad \text{comp}(3)=1, \quad \ldots$

$\text{base}(0)=0, \quad \text{base}(1)=0, \quad \text{base}(2)=1, \quad \text{base}(3)=1, \quad \ldots$

In all but very rare cases, you will not need to know which shape function $$\varphi_{\text{base}(i)}$$ of the scalar element belongs to a shape function $$\Phi_i$$ of the vector element. Let us therefore define

$\phi_i = \varphi_{\text{base}(i)}$

by which we can write the vector shape function as

$\Phi_i({\mathbf x}) = \phi_{i}({\mathbf x})\ {\mathbf e}_{\text{comp}(i)}.$

You can now safely forget about the function $$\text{base}(i)$$, at least for the rest of this example program.

Now using this vector shape functions, we can write the discrete finite element solution as

${\mathbf u}_h({\mathbf x}) = \sum_i \Phi_i({\mathbf x})\ U_i$

with scalar coefficients $$U_i$$. If we define an analog function $${\mathbf v}_h$$ as test function, we can write the discrete problem as follows: Find coefficients $$U_i$$ such that

$a({\mathbf u}_h, {\mathbf v}_h) = ({\mathbf f}, {\mathbf v}_h) \qquad \forall {\mathbf v}_h.$

If we insert the definition of the bilinear form and the representation of $${\mathbf u}_h$$ and $${\mathbf v}_h$$ into this formula:

\begin{eqnarray*} \sum_{i,j} U_i V_j \sum_{k,l} \left\{ \left( \lambda \partial_l (\Phi_i)_l, \partial_k (\Phi_j)_k \right)_\Omega + \left( \mu \partial_l (\Phi_i)_k, \partial_l (\Phi_j)_k \right)_\Omega + \left( \mu \partial_l (\Phi_i)_k, \partial_k (\Phi_j)_l \right)_\Omega \right\} \\ = \sum_j V_j \sum_l \left( f_l, (\Phi_j)_l \right)_\Omega. \end{eqnarray*}

We note that here and in the following, the indices $$k,l$$ run over spatial directions, i.e. $$0\le k,l < d$$, and that indices $$i,j$$ run over degrees of freedoms.

The local stiffness matrix on cell $$K$$ therefore has the following entries:

$A^K_{ij} = \sum_{k,l} \left\{ \left( \lambda \partial_l (\Phi_i)_l, \partial_k (\Phi_j)_k \right)_K + \left( \mu \partial_l (\Phi_i)_k, \partial_l (\Phi_j)_k \right)_K + \left( \mu \partial_l (\Phi_i)_k, \partial_k (\Phi_j)_l \right)_K \right\},$

where $$i,j$$ now are local degrees of freedom and therefore $$0\le i,j < N$$. In these formulas, we always take some component of the vector shape functions $$\Phi_i$$, which are of course given as follows (see their definition):

$(\Phi_i)_l = \phi_i \delta_{l,\text{comp}(i)},$

with the Kronecker symbol $$\delta_{nm}$$. Due to this, we can delete some of the sums over $$k$$ and $$l$$:

\begin{eqnarray*} A^K_{ij} &=& \sum_{k,l} \Bigl\{ \left( \lambda \partial_l \phi_i\ \delta_{l,\text{comp}(i)}, \partial_k \phi_j\ \delta_{k,\text{comp}(j)} \right)_K \\ &\qquad\qquad& + \left( \mu \partial_l \phi_i\ \delta_{k,\text{comp}(i)}, \partial_l \phi_j\ \delta_{k,\text{comp}(j)} \right)_K + \left( \mu \partial_l \phi_i\ \delta_{k,\text{comp}(i)}, \partial_k \phi_j\ \delta_{l,\text{comp}(j)} \right)_K \Bigr\} \\ &=& \left( \lambda \partial_{\text{comp}(i)} \phi_i, \partial_{\text{comp}(j)} \phi_j \right)_K + \sum_l \left( \mu \partial_l \phi_i, \partial_l \phi_j \right)_K \ \delta_{\text{comp}(i),\text{comp}(j)} + \left( \mu \partial_{\text{comp}(j)} \phi_i, \partial_{\text{comp}(i)} \phi_j \right)_K \\ &=& \left( \lambda \partial_{\text{comp}(i)} \phi_i, \partial_{\text{comp}(j)} \phi_j \right)_K + \left( \mu \nabla \phi_i, \nabla \phi_j \right)_K \ \delta_{\text{comp}(i),\text{comp}(j)} + \left( \mu \partial_{\text{comp}(j)} \phi_i, \partial_{\text{comp}(i)} \phi_j \right)_K. \end{eqnarray*}

Likewise, the contribution of cell $$K$$ to the right hand side vector is

\begin{eqnarray*} f^K_j &=& \sum_l \left( f_l, (\Phi_j)_l \right)_K \\ &=& \sum_l \left( f_l, \phi_j \delta_{l,\text{comp}(j)} \right)_K \\ &=& \left( f_{\text{comp}(j)}, \phi_j \right)_K. \end{eqnarray*}

This is the form in which we will implement the local stiffness matrix and right hand side vectors.

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 step-20, 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 FiniteElement::system_to_component_index etc.

# The commented program

### Include files

As usual, the first few include files are already known, so we will not comment on them further.

#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>

In this example, we need vector-valued finite elements. The support for these can be found in the following include file :

#include <deal.II/fe/fe_system.h>

We will compose the vector-valued finite elements from regular Q1 elements which can be found here, as usual:

#include <deal.II/fe/fe_q.h>

This again is C++:

#include <fstream>
#include <iostream>

The last step is as in previous programs. In particular, just like in step-7, we pack everything that's specific to this program into a namespace of its own.

namespace Step8
{
using namespace dealii;

### The ElasticProblem class template

The main class is, except for its name, almost unchanged with respect to the step-6 example.

The only change is the use of a different class for the fe variable: Instead of a concrete finite element class such as FE_Q, we now use a more generic one, FESystem. In fact, FESystem is not really a finite element itself in that it does not implement shape functions of its own. Rather, it is a class that can be used to stack several other elements together to form one vector-valued finite element. In our case, we will compose the vector-valued element of FE_Q(1) objects, as shown below in the constructor of this class.

template <int dim>
class ElasticProblem
{
public:
ElasticProblem ();
~ElasticProblem ();
void run ();
private:
void setup_system ();
void assemble_system ();
void solve ();
void refine_grid ();
void output_results (const unsigned int cycle) const;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
ConstraintMatrix hanging_node_constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};

### Right hand side values

Before going over to the implementation of the main class, we declare and define the class which describes the right hand side. This time, the right hand side is vector-valued, as is the solution, so we will describe the changes required for this in some more detail.

The first thing is that vector-valued functions have to have a constructor, since they need to pass down to the base class of how many components the function consists. The default value in the constructor of the base class is one (i.e.: a scalar function), which is why we did not need not define a constructor for the scalar function used in previous programs.

template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide ();

The next change is that we want a replacement for the value function of the previous examples. There, a second parameter component was given, which denoted which component was requested. Here, we implement a function that returns the whole vector of values at the given place at once, in the second argument of the function. The obvious name for such a replacement function is vector_value.

Secondly, in analogy to the value_list function, there is a function vector_value_list, which returns the values of the vector-valued function at several points at once:

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;
};

This is the constructor of the right hand side class. As said above, it only passes down to the base class the number of components, which is dim in the present case (one force component in each of the dim space directions).

Some people would have moved the definition of such a short function right into the class declaration. We do not do that, as a matter of style: the deal.II style guides require that class declarations contain only declarations, and that definitions are always to be found outside. This is, obviously, as much as matter of taste as indentation, but we try to be consistent in this direction.

template <int dim>
RightHandSide<dim>::RightHandSide ()
:
Function<dim> (dim)
{}

Next the function that returns the whole vector of values at the point p at once.

To prevent cases where the return vector has not previously been set to the right size we test for this case and otherwise throw an exception at the beginning of the function. Note that enforcing that output arguments already have the correct size is a convention in deal.II, and enforced almost everywhere. The reason is that we would otherwise have to check at the beginning of the function and possibly change the size of the output vector. This is expensive, and would almost always be unnecessary (the first call to the function would set the vector to the right size, and subsequent calls would only have to do redundant checks). In addition, checking and possibly resizing the vector is an operation that can not be removed if we can't rely on the assumption that the vector already has the correct size; this is in contract to the Assert call that is completely removed if the program is compiled in optimized mode.

Likewise, if by some accident someone tried to compile and run the program in only one space dimension (in which the elastic equations do not make much sense since they reduce to the ordinary Laplace equation), we terminate the program in the second assertion. The program will work just fine in 3d, however.

template <int dim>
inline
void RightHandSide<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
Assert (values.size() == dim,
ExcDimensionMismatch (values.size(), dim));
Assert (dim >= 2, ExcNotImplemented());

The rest of the function implements computing force values. We will use a constant (unit) force in x-direction located in two little circles (or spheres, in 3d) around points (0.5,0) and (-0.5,0), and y-force in an area around the origin; in 3d, the z-component of these centers is zero as well.

For this, let us first define two objects that denote the centers of these areas. Note that upon construction of the Point objects, all components are set to zero.

Point<dim> point_1, point_2;
point_1(0) = 0.5;
point_2(0) = -0.5;

If now the point p is in a circle (sphere) of radius 0.2 around one of these points, then set the force in x-direction to one, otherwise to zero:

if (((p-point_1).norm_square() < 0.2*0.2) ||
((p-point_2).norm_square() < 0.2*0.2))
values(0) = 1;
else
values(0) = 0;

Likewise, if p is in the vicinity of the origin, then set the y-force to 1, otherwise to zero:

if (p.norm_square() < 0.2*0.2)
values(1) = 1;
else
values(1) = 0;
}

Now, this is the function of the right hand side class that returns the values at several points at once. The function starts out with checking that the number of input and output arguments is equal (the sizes of the individual output vectors will be checked in the function that we call further down below). Next, we define an abbreviation for the number of points which we shall work on, to make some things simpler below.

template <int dim>
void RightHandSide<dim>::vector_value_list (const std::vector<Point<dim> > &points,
std::vector<Vector<double> > &value_list) const
{
Assert (value_list.size() == points.size(),
ExcDimensionMismatch (value_list.size(), points.size()));
const unsigned int n_points = points.size();

Finally we treat each of the points. In one of the previous examples, we have explained why the value_list/vector_value_list function had been introduced: to prevent us from calling virtual functions too frequently. On the other hand, we now need to implement the same function twice, which can lead to confusion if one function is changed but the other is not.

We can prevent this situation by calling RightHandSide::vector_value on each point in the input list. Note that by giving the full name of the function, including the class name, we instruct the compiler to explicitly call this function, and not to use the virtual function call mechanism that would be used if we had just called vector_value. This is important, since the compiler generally can't make any assumptions which function is called when using virtual functions, and it therefore can't inline the called function into the site of the call. On the contrary, here we give the fully qualified name, which bypasses the virtual function call, and consequently the compiler knows exactly which function is called and will inline above function into the present location. (Note that we have declared the vector_value function above inline, though modern compilers are also able to inline functions even if they have not been declared as inline).

It is worth noting why we go to such length explaining what we do. Using this construct, we manage to avoid any inconsistency: if we want to change the right hand side function, it would be difficult to always remember that we always have to change two functions in the same way. Using this forwarding mechanism, we only have to change a single place (the vector_value function), and the second place (the vector_value_list function) will always be consistent with it. At the same time, using virtual function call bypassing, the code is no less efficient than if we had written it twice in the first place:

for (unsigned int p=0; p<n_points; ++p)
RightHandSide<dim>::vector_value (points[p],
value_list[p]);
}

### The ElasticProblem class implementation

#### ElasticProblem::ElasticProblem

Following is the constructor of the main class. As said before, we would like to construct a vector-valued finite element that is composed of several scalar finite elements (i.e., we want to build the vector-valued element so that each of its vector components consists of the shape functions of a scalar element). Of course, the number of scalar finite elements we would like to stack together equals the number of components the solution function has, which is dim since we consider displacement in each space direction. The FESystem class can handle this: we pass it the finite element of which we would like to compose the system of, and how often it shall be repeated:

template <int dim>
ElasticProblem<dim>::ElasticProblem ()
:
dof_handler (triangulation),
fe (FE_Q<dim>(1), dim)
{}

In fact, the FESystem class has several more constructors which can perform more complex operations than just stacking together several scalar finite elements of the same type into one; we will get to know these possibilities in later examples.

#### ElasticProblem::~ElasticProblem

The destructor, on the other hand, is exactly as in step-6:

template <int dim>
ElasticProblem<dim>::~ElasticProblem ()
{
dof_handler.clear ();
}

#### ElasticProblem::setup_system

Setting up the system of equations is identical to the function used in the step-6 example. The DoFHandler class and all other classes used here are fully aware that the finite element we want to use is vector-valued, and take care of the vector-valuedness of the finite element themselves. (In fact, they do not, but this does not need to bother you: since they only need to know how many degrees of freedom there are per vertex, line and cell, and they do not ask what they represent, i.e. whether the finite element under consideration is vector-valued or whether it is, for example, a scalar Hermite element with several degrees of freedom on each vertex).

template <int dim>
void ElasticProblem<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
hanging_node_constraints.clear ();
hanging_node_constraints);
hanging_node_constraints.close ();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
dsp,
hanging_node_constraints,
/ *keep_constrained_dofs = * / true);
sparsity_pattern.copy_from (dsp);
system_matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
}

#### ElasticProblem::assemble_system

The big changes in this program are in the creation of matrix and right hand side, since they are problem-dependent. We will go through that process step-by-step, since it is a bit more complicated than in previous examples.

The first parts of this function are the same as before, however: setting up a suitable quadrature formula, initializing an FEValues object for the (vector-valued) finite element we use as well as the quadrature object, and declaring a number of auxiliary arrays. In addition, we declare the ever same two abbreviations: n_q_points and dofs_per_cell. The number of degrees of freedom per cell we now obviously ask from the composed finite element rather than from the underlying scalar Q1 element. Here, it is dim times the number of degrees of freedom per cell of the Q1 element, though this is not explicit knowledge we need to care about:

template <int dim>
void ElasticProblem<dim>::assemble_system ()
{
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<types::global_dof_index> local_dof_indices (dofs_per_cell);

As was shown in previous examples as well, we need a place where to store the values of the coefficients at all the quadrature points on a cell. In the present situation, we have two coefficients, lambda and mu.

std::vector<double> lambda_values (n_q_points);
std::vector<double> mu_values (n_q_points);

Well, we could as well have omitted the above two arrays since we will use constant coefficients for both lambda and mu, which can be declared like this. They both represent functions always returning the constant value 1.0. Although we could omit the respective factors in the assemblage of the matrix, we use them here for purpose of demonstration.

ConstantFunction<dim> lambda(1.), mu(1.);

Then again, we need to have the same for the right hand side. This is exactly as before in previous examples. However, we now have a vector-valued right hand side, which is why the data type of the rhs_values array is changed. We initialize it by n_q_points elements, each of which is a Vector<double> with dim elements.

RightHandSide<dim> right_hand_side;
std::vector<Vector<double> > rhs_values (n_q_points,

Now we can begin with the loop over all cells:

endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
cell_rhs = 0;
fe_values.reinit (cell);

Next we get the values of the coefficients at the quadrature points. Likewise for the right hand side:

rhs_values);

Then assemble the entries of the local stiffness matrix and right hand side vector. This follows almost one-to-one the pattern described in the introduction of this example. One of the few comments in place is that we can compute the number comp(i), i.e. the index of the only nonzero vector component of shape function i using the fe.system_to_component_index(i).first function call below.

(By accessing the first variable of the return value of the system_to_component_index function, you might already have guessed that there is more in it. In fact, the function returns a std::pair<unsigned int, unsigned int>, of which the first element is comp(i) and the second is the value base(i) also noted in the introduction, i.e. the index of this shape function within all the shape functions that are nonzero in this component, i.e. base(i) in the diction of the introduction. This is not a number that we are usually interested in, however.)

With this knowledge, we can assemble the local matrix contributions:

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 j=0; j<dofs_per_cell; ++j)
{
const unsigned int
component_j = fe.system_to_component_index(j).first;
for (unsigned int q_point=0; q_point<n_q_points;
++q_point)
{
+=

The first term is (lambda d_i u_i, d_j v_j) + (mu d_i u_j, d_j v_i). Note that shape_grad(i,q_point) returns the gradient of the only nonzero component of the i-th shape function at quadrature point q_point. The component comp(i) of the gradient, which is the derivative of this only nonzero vector component of the i-th shape function with respect to the comp(i)th coordinate is accessed by the appended brackets.

(
lambda_values[q_point])
+
mu_values[q_point])
+

The second term is (mu nabla u_i, nabla v_j). We need not access a specific component of the gradient, since we only have to compute the scalar product of the two gradients, of which an overloaded version of the operator* takes care, as in previous examples.

Note that by using the ?: operator, we only do this if comp(i) equals comp(j), otherwise a zero is added (which will be optimized away by the compiler).

((component_i == component_j) ?
mu_values[q_point]) :
0)
)
*
fe_values.JxW(q_point);
}
}
}

Assembling the right hand side is also just as discussed in the introduction:

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)
cell_rhs(i) += fe_values.shape_value(i,q_point) *
rhs_values[q_point](component_i) *
fe_values.JxW(q_point);
}

The transfer from local degrees of freedom into the global matrix and right hand side vector does not depend on the equation under consideration, and is thus the same as in all previous examples. The same holds for the elimination of hanging nodes from the matrix and right hand side, once we are done with assembling the entire linear system:

cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
local_dof_indices[j],
cell_matrix(i,j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
hanging_node_constraints.condense (system_matrix);
hanging_node_constraints.condense (system_rhs);

The interpolation of the boundary values needs a small modification: since the solution function is vector-valued, so need to be the boundary values. The ZeroFunction constructor accepts a parameter that tells it that it shall represent a vector valued, constant zero function with that many components. By default, this parameter is equal to one, in which case the ZeroFunction object would represent a scalar function. Since the solution vector has dim components, we need to pass dim as number of components to the zero function as well.

std::map<types::global_dof_index,double> boundary_values;
0,
boundary_values);
system_matrix,
solution,
system_rhs);
}

#### ElasticProblem::solve

The solver does not care about where the system of equations comes, as long as it stays positive definite and symmetric (which are the requirements for the use of the CG solver), which the system indeed is. Therefore, we need not change anything.

template <int dim>
void ElasticProblem<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve (system_matrix, solution, system_rhs,
preconditioner);
hanging_node_constraints.distribute (solution);
}

#### ElasticProblem::refine_grid

The function that does the refinement of the grid is the same as in the step-6 example. The quadrature formula is adapted to the linear elements again. Note that the error estimator by default adds up the estimated obtained from all components of the finite element solution, i.e., it uses the displacement in all directions with the same weight. If we would like the grid to be adapted to the x-displacement only, we could pass the function an additional parameter which tells it to do so and do not consider the displacements in all other directions for the error indicators. However, for the current problem, it seems appropriate to consider all displacement components with equal weight.

template <int dim>
void ElasticProblem<dim>::refine_grid ()
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
solution,
estimated_error_per_cell);
estimated_error_per_cell,
0.3, 0.03);
}

#### ElasticProblem::output_results

The output happens mostly as has been shown in previous examples already. The only difference is that the solution function is vector valued. The DataOut class takes care of this automatically, but we have to give each component of the solution vector a different name.

template <int dim>
void ElasticProblem<dim>::output_results (const unsigned int cycle) const
{
std::string filename = "solution-";
filename += ('0' + cycle);
Assert (cycle < 10, ExcInternalError());
filename += ".vtk";
std::ofstream output (filename.c_str());
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);

As said above, we need a different name for each component of the solution function. To pass one name for each component, a vector of strings is used. Since the number of components is the same as the number of dimensions we are working in, the following switch statement is used.

We note that some graphics programs have restriction as to what characters are allowed in the names of variables. The library therefore supports only the minimal subset of these characters that is supported by all programs. Basically, these are letters, numbers, underscores, and some other characters, but in particular no whitespace and minus/hyphen. The library will throw an exception otherwise, at least if in debug mode.

After listing the 1d, 2d, and 3d case, it is good style to let the program die if we run upon a case which we did not consider. Remember that the Assert macro generates an exception if the condition in the first parameter is not satisfied. Of course, the condition false can never be satisfied, so the program will always abort whenever it gets to the default statement:

std::vector<std::string> solution_names;
switch (dim)
{
case 1:
solution_names.push_back ("displacement");
break;
case 2:
solution_names.push_back ("x_displacement");
solution_names.push_back ("y_displacement");
break;
case 3:
solution_names.push_back ("x_displacement");
solution_names.push_back ("y_displacement");
solution_names.push_back ("z_displacement");
break;
default:
Assert (false, ExcNotImplemented());
}

After setting up the names for the different components of the solution vector, we can add the solution vector to the list of data vectors scheduled for output. Note that the following function takes a vector of strings as second argument, whereas the one which we have used in all previous examples accepted a string there. In fact, the latter function is only a shortcut for the function which we call here: it puts the single string that is passed to it into a vector of strings with only one element and forwards that to the other function.

data_out.build_patches ();
data_out.write_vtk (output);
}

#### ElasticProblem::run

The run function does the same things as in step-6, for example. This time, we use the square [-1,1]^d as domain, and we refine it twice globally before starting the first iteration.

The reason for refining twice is a bit accidental: we use the QGauss quadrature formula with two points in each direction for integration of the right hand side; that means that there are four quadrature points on each cell (in 2D). If we only refine the initial grid once globally, then there will be only four quadrature points in each direction on the domain. However, the right hand side function was chosen to be rather localized and in that case, by pure chance, it happens that all quadrature points lie at points where the the right hand side function is zero (in mathematical terms, the quadrature points happen to be at points outside the support of the right hand side function). The right hand side vector computed with quadrature will then contain only zeroes (even though it would of course be nonzero if we had computed the right hand side vector exactly using the integral) and the solution of the system of equations is the zero vector, i.e., a finite element function that is zero everywhere. In a sense, we should not be surprised that this is happening since we have chosen an initial grid that is totally unsuitable for the problem at hand.

The unfortunate thing is that if the discrete solution is constant, then the error indicators computed by the KellyErrorEstimator class are zero for each cell as well, and the call to refine_and_coarsen_fixed_number on the triangulation object will not flag any cells for refinement (why should it if the indicated error is zero for each cell?). The grid in the next iteration will therefore consist of four cells only as well, and the same problem occurs again.

The conclusion needs to be: while of course we will not choose the initial grid to be well-suited for the accurate solution of the problem, we must at least choose it such that it has the chance to capture the important features of the solution. In this case, it needs to be able to see the right hand side. Thus, we refine twice globally. (Any larger number of global refinement steps would of course also work.)

template <int dim>
void ElasticProblem<dim>::run ()
{
for (unsigned int cycle=0; cycle<8; ++cycle)
{
std::cout << "Cycle " << cycle << ':' << std::endl;
if (cycle == 0)
{
GridGenerator::hyper_cube (triangulation, -1, 1);
triangulation.refine_global (2);
}
else
refine_grid ();
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl;
setup_system ();
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
assemble_system ();
solve ();
output_results (cycle);
}
}
}

### The main function

After closing the Step8 namespace in the last line above, the following is the main function of the program and is again exactly like in step-6 (apart from the changed class names, of course).

int main ()
{
try
{
Step8::ElasticProblem<2> elastic_problem_2d;
elastic_problem_2d.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}

# Results

There is not much to be said about the results of this program, other than that they look nice. All images were made using Visit from the output files that the program wrote to disk. The first two pictures show the $$x$$- and $$y$$-displacements as a scalar components:  You can clearly see the sources of $$x$$-displacement around $$x=0.5$$ and $$x=-0.5$$, and of $$y$$-displacement at the origin. The next image shows the final grid after eight steps of refinement: What one frequently would like to do is to show the displacement as a vector field, i.e., show vectors that for each point show the direction and magnitude of displacement. Unfortunately, that's a bit more involved. To understand why this is so, remember that we have just defined our finite element as a collection of two components (in dim=2 dimensions). Nowhere have we said that this is not just a pressure and a concentration (two scalar quantities) but that the two components actually are the parts of a vector-valued quantity, namely the displacement. Absent this knowledge, the DataOut class assumes that all individual variables we print are separate scalars, and Visit then faithfully assumes that this is indeed what it is. In other words, once we have written the data as scalars, there is nothing in Visit that allows us to paste these two scalar fields back together as a vector field. Where we would have to attack this problem is at the root, namely in ElasticProblem::output_results. We won't do so here but instead refer the reader to the step-22 program where we show how to do this for a more general situation. That said, we couldn't help generating the data anyway that would show how this would look if implemented as discussed in step-22. The vector field then looks like this (Visit randomly selects a few hundred vertices from which to draw the vectors; drawing them from each individual vertex would make the picture unreadable): We note that one may have intuitively expected the solution to be symmetric about the $$x$$- and $$y$$-axes since the $$x$$- and $$y$$-forces are symmetric with respect to these axes. However, the force considered as a vector is not symmetric and consequently neither is the solution.

# The plain program

/* ---------------------------------------------------------------------
*
* Copyright (C) 2000 - 2015 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, University of Heidelberg, 2000
*/
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <fstream>
#include <iostream>
namespace Step8
{
using namespace dealii;
template <int dim>
class ElasticProblem
{
public:
ElasticProblem ();
~ElasticProblem ();
void run ();
private:
void setup_system ();
void assemble_system ();
void solve ();
void refine_grid ();
void output_results (const unsigned int cycle) const;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
ConstraintMatrix hanging_node_constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide ();
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>
RightHandSide<dim>::RightHandSide ()
:
Function<dim> (dim)
{}
template <int dim>
inline
void RightHandSide<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
Assert (values.size() == dim,
ExcDimensionMismatch (values.size(), dim));
Assert (dim >= 2, ExcNotImplemented());
Point<dim> point_1, point_2;
point_1(0) = 0.5;
point_2(0) = -0.5;
if (((p-point_1).norm_square() < 0.2*0.2) ||
((p-point_2).norm_square() < 0.2*0.2))
values(0) = 1;
else
values(0) = 0;
if (p.norm_square() < 0.2*0.2)
values(1) = 1;
else
values(1) = 0;
}
template <int dim>
void RightHandSide<dim>::vector_value_list (const std::vector<Point<dim> > &points,
std::vector<Vector<double> > &value_list) const
{
Assert (value_list.size() == points.size(),
ExcDimensionMismatch (value_list.size(), points.size()));
const unsigned int n_points = points.size();
for (unsigned int p=0; p<n_points; ++p)
RightHandSide<dim>::vector_value (points[p],
value_list[p]);
}
template <int dim>
ElasticProblem<dim>::ElasticProblem ()
:
dof_handler (triangulation),
fe (FE_Q<dim>(1), dim)
{}
template <int dim>
ElasticProblem<dim>::~ElasticProblem ()
{
dof_handler.clear ();
}
template <int dim>
void ElasticProblem<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
hanging_node_constraints.clear ();
hanging_node_constraints);
hanging_node_constraints.close ();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
dsp,
hanging_node_constraints,
/*keep_constrained_dofs = */ true);
sparsity_pattern.copy_from (dsp);
system_matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
}
template <int dim>
void ElasticProblem<dim>::assemble_system ()
{
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<types::global_dof_index> local_dof_indices (dofs_per_cell);
std::vector<double> lambda_values (n_q_points);
std::vector<double> mu_values (n_q_points);
ConstantFunction<dim> lambda(1.), mu(1.);
RightHandSide<dim> right_hand_side;
std::vector<Vector<double> > rhs_values (n_q_points,
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
cell_rhs = 0;
fe_values.reinit (cell);
rhs_values);
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 j=0; j<dofs_per_cell; ++j)
{
const unsigned int
component_j = fe.system_to_component_index(j).first;
for (unsigned int q_point=0; q_point<n_q_points;
++q_point)
{
+=
(
lambda_values[q_point])
+
mu_values[q_point])
+
((component_i == component_j) ?
mu_values[q_point]) :
0)
)
*
fe_values.JxW(q_point);
}
}
}
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)
cell_rhs(i) += fe_values.shape_value(i,q_point) *
rhs_values[q_point](component_i) *
fe_values.JxW(q_point);
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
local_dof_indices[j],
cell_matrix(i,j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
hanging_node_constraints.condense (system_matrix);
hanging_node_constraints.condense (system_rhs);
std::map<types::global_dof_index,double> boundary_values;
0,
boundary_values);
system_matrix,
solution,
system_rhs);
}
template <int dim>
void ElasticProblem<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve (system_matrix, solution, system_rhs,
preconditioner);
hanging_node_constraints.distribute (solution);
}
template <int dim>
void ElasticProblem<dim>::refine_grid ()
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
solution,
estimated_error_per_cell);
estimated_error_per_cell,
0.3, 0.03);
}
template <int dim>
void ElasticProblem<dim>::output_results (const unsigned int cycle) const
{
std::string filename = "solution-";
filename += ('0' + cycle);
Assert (cycle < 10, ExcInternalError());
filename += ".vtk";
std::ofstream output (filename.c_str());
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
std::vector<std::string> solution_names;
switch (dim)
{
case 1:
solution_names.push_back ("displacement");
break;
case 2:
solution_names.push_back ("x_displacement");
solution_names.push_back ("y_displacement");
break;
case 3:
solution_names.push_back ("x_displacement");
solution_names.push_back ("y_displacement");
solution_names.push_back ("z_displacement");
break;
default:
Assert (false, ExcNotImplemented());
}
data_out.build_patches ();
data_out.write_vtk (output);
}
template <int dim>
void ElasticProblem<dim>::run ()
{
for (unsigned int cycle=0; cycle<8; ++cycle)
{
std::cout << "Cycle " << cycle << ':' << std::endl;
if (cycle == 0)
{
GridGenerator::hyper_cube (triangulation, -1, 1);
triangulation.refine_global (2);
}
else
refine_grid ();
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl;
setup_system ();
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
assemble_system ();
solve ();
output_results (cycle);
}
}
}
int main ()
{
try
{
Step8::ElasticProblem<2> elastic_problem_2d;
elastic_problem_2d.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}