/* ---------------------------------------------------------------------
 * $Id: elastoplastic.cc 31592 2013-11-08 16:47:28Z Ghorashi $
 *
 * Copyright (C) 2012 - 2013 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
 * Public License as published by the Free Software Foundation; either
 * 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.
 *
 * ---------------------------------------------------------------------

 *
 * Authors: Seyed Shahram Ghorashi, Bauhaus-Universit\"at Weimar, 2014
 *          Joerg Frohne, Texas A&M University and
 *                        University of Siegen, 2012, 2013
 *          Wolfgang Bangerth, Texas A&M University, 2012, 2013
 *          Timo Heister, Texas A&M University, 2013
 */

// @sect3{Include files}
// The set of include files is not much of a surprise any more at this time:
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/table_handler.h>

#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/sparse_direct.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/grid_tools.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/manifold_lib.h>

#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>

#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>

#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_tools.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/numerics/fe_field_function.h>
#include <deal.II/numerics/solution_transfer.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 <deal.II/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 <deal.II/grid/filtered_iterator.h>

#include <fstream>
#include <iostream>

// This final include file provides the <code>mkdir</code> function
// that we will use to create a directory for output files, if necessary:
#include <sys/stat.h>

namespace ElastoPlastic
{
  using namespace dealii;

  void
  extrude_triangulation(const Triangulation<2, 2> &input,
                        const unsigned int n_slices,
                        const double height,
                        Triangulation<3,3> &result)
  {
    //  Assert (input.n_levels() == 1,
    //          ExcMessage ("The input triangulations must be coarse meshes."));
    Assert(result.n_cells()==0, ExcMessage("resultin Triangulation need to be empty upon calling extrude_triangulation."));
    Assert(height>0, ExcMessage("The height in extrude_triangulation needs to be positive."));
    Assert(n_slices>=2, ExcMessage("The number of slices in extrude_triangulation needs to be at least 2."));

    std::vector<Point<3> > points(n_slices*input.n_used_vertices());
    std::vector<CellData<3> > cells;
    cells.reserve((n_slices-1)*input.n_active_cells());

    for (unsigned int slice=0; slice<n_slices; ++slice)
      {
        for (unsigned int i=0; i<input.n_vertices(); ++i)

          {
            if (input.get_used_vertices()[i])
              {
                const Point<2> &v = input.get_vertices()[i];
                points[i+slice*input.n_vertices()](0) = v(0);
                points[i+slice*input.n_vertices()](1) = v(1);
                points[i+slice*input.n_vertices()](2) = height * slice / (n_slices-1);
              }
          }
      }

    for (Triangulation<2,2>::cell_iterator
         cell = input.begin_active(); cell != input.end(); ++cell)
      {
        for (unsigned int slice=0; slice<n_slices-1; ++slice)
          {
            CellData<3> this_cell;
            for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
              {
                this_cell.vertices[v]
                  = cell->vertex_index(v)+slice*input.n_used_vertices();
                this_cell.vertices[v+GeometryInfo<2>::vertices_per_cell]
                  = cell->vertex_index(v)+(slice+1)*input.n_used_vertices();
              }

            this_cell.material_id = cell->material_id();
            cells.push_back(this_cell);
          }
      }

    SubCellData s;
    types::boundary_id bid=0;
    s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2);
    for (Triangulation<2,2>::cell_iterator
         cell = input.begin_active(); cell != input.end(); ++cell)
      {
        CellData<2> quad;
        for (unsigned int f=0; f<4; ++f)
          if (cell->at_boundary(f))
            {
              quad.boundary_id = cell->face(f)->boundary_id();
              bid = std::max(bid, quad.boundary_id);
              for (unsigned int slice=0; slice<n_slices-1; ++slice)
                {
                  quad.vertices[0] = cell->face(f)->vertex_index(0)+slice*input.n_used_vertices();
                  quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_used_vertices();
                  quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_used_vertices();
                  quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_used_vertices();
                  s.boundary_quads.push_back(quad);
                }
            }
      }

    for (Triangulation<2,2>::cell_iterator
         cell = input.begin_active(); cell != input.end(); ++cell)
      {
        CellData<2> quad;
        quad.boundary_id = bid + 1;
        quad.vertices[0] = cell->vertex_index(0);
        quad.vertices[1] = cell->vertex_index(1);
        quad.vertices[2] = cell->vertex_index(2);
        quad.vertices[3] = cell->vertex_index(3);
        s.boundary_quads.push_back(quad);

        quad.boundary_id = bid + 2;
        for (int i=0; i<4; ++i)
          quad.vertices[i] += (n_slices-1)*input.n_used_vertices();
        s.boundary_quads.push_back(quad);
      }

    result.create_triangulation (points,
                                 cells,
                                 s);
  }

  namespace Evaluation
  {


    template <int dim>
    double get_von_Mises_stress(const SymmetricTensor<2, dim> &stress)
    {

      //      if (dim == 2)
      //      {
      //        von_Mises_stress = std::sqrt(  stress[0][0]*stress[0][0]
      //                                                         + stress[1][1]*stress[1][1]
      //                                                         - stress[0][0]*stress[1][1]
      //                                                         + 3*stress[0][1]*stress[0][1]);
      //      }else if (dim == 3)
      //      {
      //        von_Mises_stress = std::sqrt(  stress[0][0]*stress[0][0]
      //                                                       + stress[1][1]*stress[1][1]
      //                                                       + stress[2][2]*stress[2][2]
      //                                                       - stress[0][0]*stress[1][1]
      //                                                       - stress[1][1]*stress[2][2]
      //                                                       - stress[0][0]*stress[2][2]
      //                                                         + 3*( stress[0][1]*stress[0][1]
      //                                                              +stress[1][2]*stress[1][2]
      //                                                              +stress[0][2]*stress[0][2]) );
      //      }

      // -----------------------------------------------
      // "Perforated_strip_tension"
      // plane stress
//      const double von_Mises_stress = std::sqrt(  stress[0][0]*stress[0][0]
//                                                + stress[1][1]*stress[1][1]
//                                                - stress[0][0]*stress[1][1]
//                                                + 3*stress[0][1]*stress[0][1]);
      // -----------------------------------------------
      // otherwise
      // plane strain / 3d case
      const double von_Mises_stress = std::sqrt(1.5) * (deviator(stress)).norm();
      // -----------------------------------------------



      return von_Mises_stress;
    }


    template <int dim>
    class PointValuesEvaluation
    {
    public:
      PointValuesEvaluation (const Point<dim>  &evaluation_point);

      void compute (const DoFHandler<dim>  &dof_handler,
                    const Vector<double>   &solution,
                    Vector<double>         &point_values);

      DeclException1 (ExcEvaluationPointNotFound,
                      Point<dim>,
                      << "The evaluation point " << arg1
                      << " was not found among the vertices of the present grid.");
    private:
      const Point<dim>  evaluation_point;
    };


    template <int dim>
    PointValuesEvaluation<dim>::
    PointValuesEvaluation (const Point<dim>  &evaluation_point)
      :
      evaluation_point (evaluation_point)
    {}



    template <int dim>
    void
    PointValuesEvaluation<dim>::
    compute (const DoFHandler<dim>  &dof_handler,
             const Vector<double>   &solution,
             Vector<double>         &point_values)
    {
      const unsigned int dofs_per_vertex = dof_handler.get_fe().dofs_per_vertex;
      AssertThrow (point_values.size() == dofs_per_vertex,
                   ExcDimensionMismatch (point_values.size(), dofs_per_vertex));
      point_values = 1e20;

      typename DoFHandler<dim>::active_cell_iterator
      cell = dof_handler.begin_active(),
      endc = dof_handler.end();
      bool evaluation_point_found = false;
      for (; (cell!=endc) && !evaluation_point_found; ++cell)
        {
          if (cell->is_locally_owned() && !evaluation_point_found)
            for (unsigned int vertex=0;
                 vertex<GeometryInfo<dim>::vertices_per_cell;
                 ++vertex)
              {
                if (cell->vertex(vertex).distance (evaluation_point)
                    <
                    cell->diameter() * 1e-8)
                  {
                    for (unsigned int id=0; id!=dofs_per_vertex; ++id)
                      {
                        point_values[id] = solution(cell->vertex_dof_index(vertex,id));
                      }

                    evaluation_point_found = true;
                    break;
                  }
              }
        }

      AssertThrow (evaluation_point_found,
                   ExcEvaluationPointNotFound(evaluation_point));
    }


  }

  // @sect3{The <code>PointHistory</code> 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 <code>struct</code> rather than <code>class</code>, to
  // indicate that they are closer to C-style structures than C++-style
  // classes.
  template <int dim>
  struct PointHistory
  {
    SymmetricTensor<2,dim> old_stress;
    SymmetricTensor<2,dim> old_strain;
    Point<dim> point;
  };


  // @sect3{The <code>ConstitutiveLaw</code> class template}

  // This class provides an interface for a constitutive law, i.e., for the
  // relationship between strain $\varepsilon(\mathbf u)$ and stress
  // $\sigma$. In this example we are using an elastoplastic material behavior
  // with linear, isotropic hardening. Such materials are characterized by
  // Young's modulus $E$, Poisson's ratio $\nu$, the initial yield stress
  // $\sigma_0$ and the isotropic hardening parameter $\gamma$.  For $\gamma =
  // 0$ we obtain perfect elastoplastic behavior.
  //
  // As explained in the paper that describes this program, the first Newton
  // steps are solved with a completely elastic material model to avoid having
  // to deal with both nonlinearities (plasticity and contact) at once. To this
  // end, this class has a function <code>set_sigma_0()</code> that we use later
  // on to simply set $\sigma_0$ to a very large value -- essentially
  // guaranteeing that the actual stress will not exceed it, and thereby
  // producing an elastic material. When we are ready to use a plastic model, we
  // set $\sigma_0$ back to its proper value, using the same function.  As a
  // result of this approach, we need to leave <code>sigma_0</code> as the only
  // non-const member variable of this class.
  template <int dim>
  class ConstitutiveLaw
  {
  public:
    ConstitutiveLaw (const double E,
                     const double nu,
                     const double sigma_0,
                     const double gamma);

    void
    set_sigma_0 (double sigma_zero);

    bool
    get_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
                              SymmetricTensor<4, dim> &stress_strain_tensor) const;

    bool
    get_grad_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
                                   const std::vector<Tensor<2, dim> > &point_hessian,
                                   Tensor<5, dim> &stress_strain_tensor_grad) const;

    void
    get_linearized_stress_strain_tensors (const SymmetricTensor<2, dim> &strain_tensor,
                                          SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
                                          SymmetricTensor<4, dim> &stress_strain_tensor) const;

  private:
    const double kappa;
    const double mu;
    double       sigma_0;
    const double gamma;

    const SymmetricTensor<4, dim> stress_strain_tensor_kappa;
    const SymmetricTensor<4, dim> stress_strain_tensor_mu;
  };

  // The constructor of the ConstitutiveLaw class sets the required material
  // parameter for our deformable body. Material parameters for elastic
  // isotropic media can be defined in a variety of ways, such as the pair $E,
  // \nu$ (elastic modulus and Poisson's number), using the Lame parameters
  // $\lambda,mu$ or several other commonly used conventions. Here, the
  // constructor takes a description of material parameters in the form of
  // $E,\nu$, but since this turns out to these are not the coefficients that
  // appear in the equations of the plastic projector, we immediately convert
  // them into the more suitable set $\kappa,\mu$ of bulk and shear moduli.  In
  // addition, the constructor takes $\sigma_0$ (the yield stress absent any
  // plastic strain) and $\gamma$ (the hardening parameter) as arguments. In
  // this constructor, we also compute the two principal components of the
  // stress-strain relation and its linearization.
  template <int dim>
  ConstitutiveLaw<dim>::ConstitutiveLaw (double E,
                                         double nu,
                                         double sigma_0,
                                         double gamma)
    :
    //--------------------
    // Plane stress
//    kappa (((E*(1+2*nu)) / (std::pow((1+nu),2))) / (3 * (1 - 2 * (nu / (1+nu))))),
//    mu (((E*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu))))),
    //--------------------
    // 3d and plane strain
    kappa (E / (3 * (1 - 2 * nu))),
    mu (E / (2 * (1 + nu))),
    //--------------------
    sigma_0(sigma_0),
    gamma(gamma),
    stress_strain_tensor_kappa (kappa
                                * outer_product(unit_symmetric_tensor<dim>(),
                                                unit_symmetric_tensor<dim>())),
    stress_strain_tensor_mu (2 * mu
                             * (identity_tensor<dim>()
                                - outer_product(unit_symmetric_tensor<dim>(),
                                                unit_symmetric_tensor<dim>()) / 3.0))
  {}


  template <int dim>
  void
  ConstitutiveLaw<dim>::set_sigma_0 (double sigma_zero)
  {
    sigma_0 = sigma_zero;
  }


  // @sect4{ConstitutiveLaw::get_stress_strain_tensor}

  // This is the principal component of the constitutive law. It projects the
  // deviatoric part of the stresses in a quadrature point back to the yield
  // stress (i.e., the original yield stress $\sigma_0$ plus the term that
  // describes linear isotropic hardening).  We need this function to calculate
  // the nonlinear residual in PlasticityContactProblem::residual_nl_system. The
  // computations follow the formulas laid out in the introduction.
  //
  // The function returns whether the quadrature point is plastic to allow for
  // some statistics downstream on how many of the quadrature points are
  // plastic and how many are elastic.
  template <int dim>
  bool
  ConstitutiveLaw<dim>::
  get_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
                            SymmetricTensor<4, dim> &stress_strain_tensor) const
  {
    SymmetricTensor<2, dim> stress_tensor;
    stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
                    * strain_tensor;

//    const SymmetricTensor<2, dim> deviator_stress_tensor = deviator(stress_tensor);
//    const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
    const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);

    stress_strain_tensor = stress_strain_tensor_mu;
    if (von_Mises_stress > sigma_0)
      {
        const double beta = sigma_0 / von_Mises_stress;
        stress_strain_tensor *= (gamma + (1 - gamma) * beta);
      }

    stress_strain_tensor += stress_strain_tensor_kappa;

    return (von_Mises_stress > sigma_0);
  }


  template <int dim>
  bool
  ConstitutiveLaw<dim>::
  get_grad_stress_strain_tensor (const SymmetricTensor<2, dim> &strain_tensor,
                                 const std::vector<Tensor<2, dim> > &point_hessian,
                                 Tensor<5, dim> &stress_strain_tensor_grad) const
  {
    SymmetricTensor<2, dim> stress_tensor;
    stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
                    * strain_tensor;

    const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);

    if (von_Mises_stress > sigma_0)
      {
        const SymmetricTensor<2, dim> deviator_strain_tensor = deviator(strain_tensor);
        const double deviator_strain_tensor_norm = deviator_strain_tensor.norm();
        const double multiplier = -(1-gamma)*sigma_0/(2*mu*std::pow(deviator_strain_tensor_norm,3));

        Vector<double> multiplier_vector(dim);
        multiplier_vector = 0;

        for (unsigned int i=0; i!=dim; ++i)
          for (unsigned int m=0; m!=dim; ++m)
            for (unsigned int n=0; n!=dim; ++n)
              {
                multiplier_vector(i) += deviator_strain_tensor[m][n] *
                                        ( 0.5*( point_hessian[m][n][i] + point_hessian[n][m][i] )
                                          + ( m==n && dim==2 ? -1/dim*(point_hessian[0][0][i]
                                                                       + point_hessian[1][1][i]) : 0 )
                                          + ( m==n && dim==3 ? -1/dim*(point_hessian[0][0][i]
                                                                       + point_hessian[1][1][i]
                                                                       + point_hessian[2][2][i]) : 0 ) );
              }

        // -----------------------------------------------
        // "Perforated_strip_tension"
        // plane stress
//      const double VM_factor = std::sqrt(2);
        // -----------------------------------------------
        // otherwise
        // plane strain / 3d case
        const double VM_factor = std::sqrt(1.5);
        // -----------------------------------------------

        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)
                for (unsigned int m=0; m!=dim; ++m)
                  {
                    stress_strain_tensor_grad[i][j][k][l][m] = 1/VM_factor
                                                               * multiplier
                                                               * stress_strain_tensor_mu[i][j][k][l]
                                                               * multiplier_vector(m);
                  }

      }
    else
      {
        stress_strain_tensor_grad = 0;
      }

    return (von_Mises_stress > sigma_0);
  }


  // @sect4{ConstitutiveLaw::get_linearized_stress_strain_tensors}

  // This function returns the linearized stress strain tensor, linearized
  // around the solution $u^{i-1}$ of the previous Newton step $i-1$.  The
  // parameter <code>strain_tensor</code> (commonly denoted
  // $\varepsilon(u^{i-1})$) must be passed as an argument, and serves as the
  // linearization point. The function returns the derivative of the nonlinear
  // constitutive law in the variable stress_strain_tensor, as well as the
  // stress-strain tensor of the linearized problem in
  // stress_strain_tensor_linearized.  See
  // PlasticityContactProblem::assemble_nl_system where this function is used.
  template <int dim>
  void
  ConstitutiveLaw<dim>::
  get_linearized_stress_strain_tensors (const SymmetricTensor<2, dim> &strain_tensor,
                                        SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
                                        SymmetricTensor<4, dim> &stress_strain_tensor) const
  {
    SymmetricTensor<2, dim> stress_tensor;
    stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
                    * strain_tensor;

    stress_strain_tensor = stress_strain_tensor_mu;
    stress_strain_tensor_linearized = stress_strain_tensor_mu;

    SymmetricTensor<2, dim> deviator_stress_tensor = deviator(stress_tensor);
    const double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
    const double von_Mises_stress = Evaluation::get_von_Mises_stress(stress_tensor);

    if (von_Mises_stress > sigma_0)
      {
        const double beta = sigma_0 / von_Mises_stress;
        stress_strain_tensor *= (gamma + (1 - gamma) * beta);
        stress_strain_tensor_linearized *= (gamma + (1 - gamma) * beta);
        deviator_stress_tensor /= deviator_stress_tensor_norm;
        stress_strain_tensor_linearized -= (1 - gamma) * beta * 2 * mu
                                           * outer_product(deviator_stress_tensor,
                                                           deviator_stress_tensor);
      }

    stress_strain_tensor += stress_strain_tensor_kappa;
    stress_strain_tensor_linearized += stress_strain_tensor_kappa;
  }

  // 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
    // <code>tan_angle==0</code>. Further down, we need to divide by that
    // number in the computation of the axis of rotation, and we would get
    // into trouble when 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 } };
        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);
  }


  // <h3>Equation data: Body forces, boundary forces,
  // incremental boundary values</h3>
  //
  // The following should be relatively standard. We need classes for
  // the boundary forcing term (which we here choose to be zero)
  // and incremental boundary values.
  namespace EquationData
  {

    /*
    template <int dim>
    class BoundaryForce : public Function<dim>
    {
    public:
      BoundaryForce ();

      virtual
      double value (const Point<dim> &p,
                    const unsigned int component = 0) const;

      virtual
      void vector_value (const Point<dim> &p,
                         Vector<double> &values) const;
    };

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


    template <int dim>
    double
    BoundaryForce<dim>::value (const Point<dim> &,
                               const unsigned int) const
    {
      return 0.;
    }

    template <int dim>
    void
    BoundaryForce<dim>::vector_value (const Point<dim> &p,
                                      Vector<double> &values) const
    {
      for (unsigned int c = 0; c < this->n_components; ++c)
        values(c) = BoundaryForce<dim>::value(p, c);
    }

    // @sect3{The <code>BodyForce</code> class}
    // 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 <code>rho g</code>, where <code>rho</code> is the material
    // density and <code>g</code> is a vector in negative z-direction with
    // magnitude 9.81 m/s^2.  Both the density and <code>g</code> are defined in
    // the function, and we take as the density 7700 kg/m^3, a value commonly
    // assumed for steel.
    //
    // To be a little more general and to be able to do computations in 2d as
    // well, we realize that the body force is always a function returning a
    // <code>dim</code> dimensional vector. We assume that gravity acts along
    // the negative direction of the last, i.e. <code>dim-1</code>th
    // coordinate. The rest of the implementation of this function should be
    // mostly self-explanatory given similar 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
    // <code>vector_value</code> 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]);
    }

    // @sect3{The <code>IncrementalBoundaryValue</code> 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 <code>ZeroFunction</code>) 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]);
    }
    */

    // ----------------------------- TimoshenkoBeam ---------------------------------------
    /*
    template <int dim>
    class IncrementalBoundaryForce : public Function<dim>
    {
    public:
      IncrementalBoundaryForce (const double present_time,
                                const double end_time);

      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 present_time,
                   end_time,
                   shear_force,
                   length,
                   depth,
                   thickness;
    };

    template <int dim>
    IncrementalBoundaryForce<dim>::
    IncrementalBoundaryForce (const double present_time,
                                const double end_time)
    :
    Function<dim>(dim),
    present_time (present_time),
    end_time (end_time),
    shear_force (2e4),
    length (.48),
    depth (.12),
    thickness (.01)
    {}

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

      // compute traction on the right face of Timoshenko beam problem, t_bar
      double inertia_moment = (thickness*std::pow(depth,3)) / 12;

      double x = p(0);
      double y = p(1);

      AssertThrow(std::fabs(x-length)<1e-12, ExcNotImplemented());

      values(0) = 0;
      values(1) = - shear_force/(2*inertia_moment) * ( depth*depth/4-y*y );

      // compute the fraction of imposed force
      const double frac = present_time/end_time;

      values *= frac;
    }

    template <int dim>
    void
    IncrementalBoundaryForce<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)
        IncrementalBoundaryForce<dim>::vector_value (points[p],
            value_list[p]);
    }


    template <int dim>
    class BodyForce :  public ZeroFunction<dim>
    {
    public:
      BodyForce () : ZeroFunction<dim> (dim) {}
    };

    template <int dim>
    class IncrementalBoundaryValues :  public Function<dim>
    {
    public:
      IncrementalBoundaryValues (const double present_time,
                                 const double end_time);

      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 present_time,
                   end_time,
                   shear_force,
                   Youngs_modulus,
                   Poissons_ratio,
                   length,
                   depth,
                   thickness;
    };


    template <int dim>
    IncrementalBoundaryValues<dim>::
    IncrementalBoundaryValues (const double present_time,
                               const double end_time)
    :
    Function<dim> (dim),
    present_time (present_time),
    end_time (end_time),
    shear_force (2e4),
    Youngs_modulus (2.e11),
    Poissons_ratio (.3),
    length (.48),
    depth (.12),
    thickness (.01)
    {}


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


      // compute exact displacement of Timoshenko beam problem, u_bar
      double inertia_moment = (thickness*std::pow(depth,3)) / 12;

      double x = p(0);
      double y = p(1);

      double fac = shear_force / (6*Youngs_modulus*inertia_moment);

      values(0) =  fac * y * ( (6*length-3*x)*x + (2+Poissons_ratio)*(y*y-depth*depth/4) );
      values(1) = -fac* ( 3*Poissons_ratio*y*y*(length-x) + 0.25*(4+5*Poissons_ratio)*depth*depth*x + (3*length-x)*x*x );

      // compute the fraction of imposed force
      const double frac = present_time/end_time;

      values *= frac;
    }



    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]);
    }
    */

    // ------------------------- Thick_tube_internal_pressure ----------------------------------
    /*
    template <int dim>
      class IncrementalBoundaryForce : public Function<dim>
      {
      public:
        IncrementalBoundaryForce (const double present_time,
                                  const double end_time);

        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 present_time,
                     end_time,
                     pressure,
                     inner_radius;
      };

      template <int dim>
      IncrementalBoundaryForce<dim>::
      IncrementalBoundaryForce (const double present_time,
                                const double end_time)
      :
      Function<dim>(dim),
      present_time (present_time),
      end_time (end_time),
      pressure (0.6*2.4e8),
    //    pressure (1.94e8),
      inner_radius(.1)
      {}

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

        const double eps = 1.e-7 * inner_radius,
                     radius = p.norm();
        // compute traction on the inner boundary, t_bar
        AssertThrow(radius < (eps+inner_radius), ExcInternalError());

        const double theta = std::atan2(p(1),p(0));

        values(0) = pressure * std::cos(theta);
        values(1) = pressure * std::sin(theta);

        // compute the fraction of imposed force
        const double frac = present_time/end_time;

        values *= frac;
      }

      template <int dim>
      void
      IncrementalBoundaryForce<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)
          IncrementalBoundaryForce<dim>::vector_value (points[p],
              value_list[p]);
      }


      template <int dim>
      class BodyForce :  public ZeroFunction<dim>
      {
      public:
        BodyForce () : ZeroFunction<dim> (dim) {}
      };


      template <int dim>
      class IncrementalBoundaryValues :  public Function<dim>
      {
      public:
        IncrementalBoundaryValues (const double present_time,
                                   const double end_time);

        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 present_time,
                     end_time;
      };


      template <int dim>
      IncrementalBoundaryValues<dim>::
      IncrementalBoundaryValues (const double present_time,
                                 const double end_time)
      :
      Function<dim> (dim),
      present_time (present_time),
      end_time (end_time)
      {}


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

        values = 0.;
      }



      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]);
      }
      */

    // ------------------------- Perforated_strip_tension ----------------------------------
    /*
    template <int dim>
    class IncrementalBoundaryForce : public Function<dim>
    {
    public:
      IncrementalBoundaryForce (const double present_time,
                                const double end_time);

      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 present_time,
                   end_time;
    };

    template <int dim>
    IncrementalBoundaryForce<dim>::
    IncrementalBoundaryForce (const double present_time,
                              const double end_time)
    :
    Function<dim>(dim),
    present_time (present_time),
    end_time (end_time)
    {}

    template <int dim>
    void
    IncrementalBoundaryForce<dim>::vector_value (const Point<dim> &p,
                                                 Vector<double> &values) const
    {
      AssertThrow (values.size() == dim,
                   ExcDimensionMismatch (values.size(), dim));

      values = 0;

      // compute the fraction of imposed force
      const double frac = present_time/end_time;

      values *= frac;
    }

    template <int dim>
    void
    IncrementalBoundaryForce<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)
        IncrementalBoundaryForce<dim>::vector_value (points[p],
            value_list[p]);
    }


    template <int dim>
    class BodyForce :  public ZeroFunction<dim>
    {
    public:
      BodyForce () : ZeroFunction<dim> (dim) {}
    };


    template <int dim>
    class IncrementalBoundaryValues :  public Function<dim>
    {
    public:
      IncrementalBoundaryValues (const double present_time,
                                 const double end_time);

      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 present_time,
                   end_time,
                   imposed_displacement,
                   height;
    };


    template <int dim>
    IncrementalBoundaryValues<dim>::
    IncrementalBoundaryValues (const double present_time,
                               const double end_time)
    :
    Function<dim> (dim),
    present_time (present_time),
    end_time (end_time),
    imposed_displacement (0.00055),
    height (0.18)
    {}


    template <int dim>
    void
    IncrementalBoundaryValues<dim>::
    vector_value (const Point<dim> &p,
                  Vector<double>   &values) const
    {
      AssertThrow (values.size() == dim,
                   ExcDimensionMismatch (values.size(), dim));

      const double eps = 1.e-8 * height;

      values = 0.;

      // impose displacement only on the top edge
      if (std::abs(p[1]-height) < eps)
      {
        // compute the fraction of imposed displacement
        const double inc_frac = 1/end_time;

        values(1) = inc_frac*imposed_displacement;
      }

    }



    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]);
    }
    */

    // ------------------------- Cantiliver_beam_3d ----------------------------------
    template <int dim>
    class IncrementalBoundaryForce : public Function<dim>
    {
    public:
      IncrementalBoundaryForce (const double present_time,
                                const double end_time);

      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 present_time,
            end_time,
            pressure,
            height;
    };

    template <int dim>
    IncrementalBoundaryForce<dim>::
    IncrementalBoundaryForce (const double present_time,
                              const double end_time)
      :
      Function<dim>(dim),
      present_time (present_time),
      end_time (end_time),
      pressure (6e6),
      height (200e-3)
    {}

    template <int dim>
    void
    IncrementalBoundaryForce<dim>::vector_value (const Point<dim> &p,
                                                 Vector<double> &values) const
    {
      AssertThrow (dim == 3, ExcNotImplemented());
      AssertThrow (values.size() == dim,
                   ExcDimensionMismatch (values.size(), dim));

      const double eps = 1.e-7 * height;

      // pressure should be imposed on the top surface, y = height
      AssertThrow(std::abs(p[1]-(height/2)) < eps, ExcInternalError());

      values = 0;

      values(1) = -pressure;

      // compute the fraction of imposed force
      const double frac = present_time/end_time;

      values *= frac;
    }

    template <int dim>
    void
    IncrementalBoundaryForce<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)
        IncrementalBoundaryForce<dim>::vector_value (points[p], value_list[p]);
    }


    template <int dim>
    class BodyForce :  public ZeroFunction<dim>
    {
    public:
      BodyForce () : ZeroFunction<dim> (dim) {}
    };


    template <int dim>
    class IncrementalBoundaryValues :  public Function<dim>
    {
    public:
      IncrementalBoundaryValues (const double present_time,
                                 const double end_time);

      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 present_time,
            end_time;
    };


    template <int dim>
    IncrementalBoundaryValues<dim>::
    IncrementalBoundaryValues (const double present_time,
                               const double end_time)
      :
      Function<dim> (dim),
      present_time (present_time),
      end_time (end_time)
    {}


    template <int dim>
    void
    IncrementalBoundaryValues<dim>::
    vector_value (const Point<dim> &/*p*/,
                  Vector<double>   &values) const
    {
      AssertThrow (values.size() == dim,
                   ExcDimensionMismatch (values.size(), dim));
      AssertThrow (dim == 3, ExcNotImplemented());

      values = 0.;
    }


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

    // -------------------------------------------------------------------------------
  }


  namespace DualFunctional
  {

    template <int dim>
    class DualFunctionalBase : public Subscriptor
    {
    public:
      virtual
      void
      assemble_rhs (const DoFHandler<dim>      &dof_handler,
                    const Vector<double>       &solution,
                    const ConstitutiveLaw<dim> &constitutive_law,
                    const DoFHandler<dim>      &dof_handler_dual,
                    Vector<double>             &rhs_dual) const = 0;
    };


    template <int dim>
    class PointValuesEvaluation : public DualFunctionalBase<dim>
    {
    public:
      PointValuesEvaluation (const Point<dim> &evaluation_point);

      virtual
      void
      assemble_rhs (const DoFHandler<dim>      &dof_handler,
                    const Vector<double>       &solution,
                    const ConstitutiveLaw<dim> &constitutive_law,
                    const DoFHandler<dim>      &dof_handler_dual,
                    Vector<double>             &rhs_dual) const;

      DeclException1 (ExcEvaluationPointNotFound,
                      Point<dim>,
                      << "The evaluation point " << arg1
                      << " was not found among the vertices of the present grid.");

    protected:
      const Point<dim> evaluation_point;
    };


    template <int dim>
    PointValuesEvaluation<dim>::
    PointValuesEvaluation (const Point<dim> &evaluation_point)
      :
      evaluation_point (evaluation_point)
    {}


    template <int dim>
    void
    PointValuesEvaluation<dim>::
    assemble_rhs (const DoFHandler<dim>      &/*dof_handler*/,
                  const Vector<double>       &/*solution*/,
                  const ConstitutiveLaw<dim> &/*constitutive_law*/,
                  const DoFHandler<dim>      &dof_handler_dual,
                  Vector<double>             &rhs_dual) const
    {
      rhs_dual.reinit (dof_handler_dual.n_dofs());
      const unsigned int dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex;

      typename DoFHandler<dim>::active_cell_iterator
      cell_dual = dof_handler_dual.begin_active(),
      endc_dual = dof_handler_dual.end();
      for (; cell_dual!=endc_dual; ++cell_dual)
        for (unsigned int vertex=0;
             vertex<GeometryInfo<dim>::vertices_per_cell;
             ++vertex)
          if (cell_dual->vertex(vertex).distance(evaluation_point)
              < cell_dual->diameter()*1e-8)
            {
              for (unsigned int id=0; id!=dofs_per_vertex; ++id)
                {
                  rhs_dual(cell_dual->vertex_dof_index(vertex,id)) = 1;
                }
              return;
            }

      AssertThrow (false, ExcEvaluationPointNotFound(evaluation_point));
    }


    template <int dim>
    class PointXDerivativesEvaluation : public DualFunctionalBase<dim>
    {
    public:
      PointXDerivativesEvaluation (const Point<dim> &evaluation_point);

      virtual
      void
      assemble_rhs (const DoFHandler<dim>      &dof_handler,
                    const Vector<double>       &solution,
                    const ConstitutiveLaw<dim> &constitutive_law,
                    const DoFHandler<dim>      &dof_handler_dual,
                    Vector<double>             &rhs_dual) const;

      DeclException1 (ExcEvaluationPointNotFound,
                      Point<dim>,
                      << "The evaluation point " << arg1
                      << " was not found among the vertices of the present grid.");

    protected:
      const Point<dim> evaluation_point;
    };


    template <int dim>
    PointXDerivativesEvaluation<dim>::
    PointXDerivativesEvaluation (const Point<dim> &evaluation_point)
      :
      evaluation_point (evaluation_point)
    {}


    template <int dim>
    void
    PointXDerivativesEvaluation<dim>::
    assemble_rhs (const DoFHandler<dim>      &/*dof_handler*/,
                  const Vector<double>       &/*solution*/,
                  const ConstitutiveLaw<dim> &/*constitutive_law*/,
                  const DoFHandler<dim>      &dof_handler_dual,
                  Vector<double>             &rhs_dual) const
    {
      rhs_dual.reinit (dof_handler_dual.n_dofs());

      QGauss<dim> quadrature(4);
      FEValues<dim>  fe_values (dof_handler_dual.get_fe(), quadrature,
                                update_gradients |
                                update_quadrature_points  |
                                update_JxW_values);
      const unsigned int n_q_points = fe_values.n_quadrature_points;
      Assert ( n_q_points==quadrature.size() , ExcInternalError() );
      const unsigned int dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell;

      Vector<double> cell_rhs (dofs_per_cell);
      std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

      double total_volume = 0;

      typename DoFHandler<dim>::active_cell_iterator
      cell = dof_handler_dual.begin_active(),
      endc = dof_handler_dual.end();
      for (; cell!=endc; ++cell)
        if (cell->center().distance(evaluation_point) <=
            cell->diameter())
          {
            fe_values.reinit (cell);
            cell_rhs = 0;

            for (unsigned int q=0; q<n_q_points; ++q)
              {
                for (unsigned int i=0; i<dofs_per_cell; ++i)
                  {
                    cell_rhs(i) += fe_values.shape_grad(i,q)[0] *
                                   fe_values.JxW (q);
                  }

                total_volume += fe_values.JxW (q);
              }

            cell->get_dof_indices (local_dof_indices);
            for (unsigned int i=0; i<dofs_per_cell; ++i)
              {
                rhs_dual(local_dof_indices[i]) += cell_rhs(i);
              }
          }

      AssertThrow (total_volume > 0,
                   ExcEvaluationPointNotFound(evaluation_point));

      rhs_dual *= 1./total_volume;
    }



    template <int dim>
    class MeanDisplacementFace : public DualFunctionalBase<dim>
    {
    public:
      MeanDisplacementFace (const unsigned int face_id,
                            const std::vector<bool> comp_mask);

      virtual
      void
      assemble_rhs (const DoFHandler<dim>      &dof_handler,
                    const Vector<double>       &solution,
                    const ConstitutiveLaw<dim> &constitutive_law,
                    const DoFHandler<dim>      &dof_handler_dual,
                    Vector<double>             &rhs_dual) const;

    protected:
      const unsigned int face_id;
      const std::vector<bool> comp_mask;
    };


    template <int dim>
    MeanDisplacementFace<dim>::
    MeanDisplacementFace (const unsigned int face_id,
                          const std::vector<bool> comp_mask )
      :
      face_id (face_id),
      comp_mask (comp_mask)
    {
      AssertThrow(comp_mask.size() == dim,
                  ExcDimensionMismatch (comp_mask.size(), dim) );
    }


    template <int dim>
    void
    MeanDisplacementFace<dim>::
    assemble_rhs (const DoFHandler<dim>      &/*dof_handler*/,
                  const Vector<double>       &/*solution*/,
                  const ConstitutiveLaw<dim> &/*constitutive_law*/,
                  const DoFHandler<dim>      &dof_handler_dual,
                  Vector<double>             &rhs_dual) const
    {
      AssertThrow (dim >= 2, ExcNotImplemented());

      rhs_dual.reinit (dof_handler_dual.n_dofs());

      const QGauss<dim-1> face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1);
      FEFaceValues<dim> fe_face_values (dof_handler_dual.get_fe(), face_quadrature,
                                        update_values | update_JxW_values);

      const unsigned int  dofs_per_vertex = dof_handler_dual.get_fe().dofs_per_vertex;
      const unsigned int  dofs_per_cell = dof_handler_dual.get_fe().dofs_per_cell;
      const unsigned int  n_face_q_points = face_quadrature.size();

      AssertThrow(dofs_per_vertex == dim,
                  ExcDimensionMismatch (dofs_per_vertex, dim) );

      std::vector<unsigned int> comp_vector(dofs_per_vertex);
      for (unsigned int i=0; i!=dofs_per_vertex; ++i)
        {
          if (comp_mask[i])
            {
              comp_vector[i] = 1;
            }
        }

      Vector<double>       cell_rhs (dofs_per_cell);

      std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

      // bound_size : size of the boundary, in 2d is the length
      //              and in the 3d case, area
      double bound_size = 0.;

      typename DoFHandler<dim>::active_cell_iterator
      cell = dof_handler_dual.begin_active(),
      endc = dof_handler_dual.end();
      bool evaluation_face_found = false;
      for (; cell!=endc; ++cell)
        {
          cell_rhs = 0;
          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (cell->face(face)->at_boundary()
                  &&
                  cell->face(face)->boundary_id() == face_id)
                {
                  if (!evaluation_face_found)
                    {
                      evaluation_face_found = true;
                    }
                  fe_face_values.reinit (cell, face);

                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      bound_size += fe_face_values.JxW(q_point);

                      for (unsigned int i=0; i<dofs_per_cell; ++i)
                        {
                          const unsigned int
                          component_i = dof_handler_dual.get_fe().system_to_component_index(i).first;

                          cell_rhs(i) += (fe_face_values.shape_value(i,q_point) *
                                          comp_vector[component_i] *
                                          fe_face_values.JxW(q_point));
                        }

                    }

                }
            }

          cell->get_dof_indices (local_dof_indices);
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            {
              rhs_dual(local_dof_indices[i]) += cell_rhs(i);
            }

        }

      AssertThrow(evaluation_face_found, ExcInternalError());

      rhs_dual /= bound_size;
    }



    template <int dim>
    class MeanStressFace : public DualFunctionalBase<dim>
    {
    public:
      MeanStressFace (const unsigned int face_id,
                      const std::vector<std::vector<unsigned int> > &comp_stress);

      virtual
      void
      assemble_rhs (const DoFHandler<dim>      &dof_handler,
                    const Vector<double>       &solution,
                    const ConstitutiveLaw<dim> &constitutive_law,
                    const DoFHandler<dim>      &dof_handler_dual,
                    Vector<double>             &rhs_dual) const;

    protected:
      const unsigned int face_id;
      const std::vector<std::vector<unsigned int> >  comp_stress;
    };


    template <int dim>
    MeanStressFace<dim>::
    MeanStressFace (const unsigned int face_id,
                    const std::vector<std::vector<unsigned int> > &comp_stress )
      :
      face_id (face_id),
      comp_stress (comp_stress)
    {
      AssertThrow(comp_stress.size() == dim,
                  ExcDimensionMismatch (comp_stress.size(), dim) );
    }


    template <int dim>
    void
    MeanStressFace<dim>::
    assemble_rhs (const DoFHandler<dim>      &dof_handler,
                  const Vector<double>       &solution,
                  const ConstitutiveLaw<dim> &constitutive_law,
                  const DoFHandler<dim>      &dof_handler_dual,
                  Vector<double>             &rhs_dual) const
    {
      AssertThrow (dim >= 2, ExcNotImplemented());

      rhs_dual.reinit (dof_handler_dual.n_dofs());

      const QGauss<dim-1> face_quadrature(dof_handler_dual.get_fe().tensor_degree()+1);

      FEFaceValues<dim> fe_face_values (dof_handler.get_fe(), face_quadrature,
                                        update_gradients);
      FEFaceValues<dim> fe_face_values_dual (dof_handler_dual.get_fe(), face_quadrature,
                                             update_gradients | update_JxW_values);

      const unsigned int  dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell;
      const unsigned int  n_face_q_points = face_quadrature.size();

      std::vector<SymmetricTensor<2, dim> > strain_tensor(n_face_q_points);
      SymmetricTensor<4, dim> stress_strain_tensor;

      Vector<double>      cell_rhs (dofs_per_cell_dual);

      std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_dual);

      // bound_size : size of the boundary, in 2d is the length
      //              and in the 3d case, area
      double bound_size = 0.;

      bool evaluation_face_found = false;

      typename DoFHandler<dim>::active_cell_iterator
      cell_dual = dof_handler_dual.begin_active(),
      endc_dual = dof_handler_dual.end(),
      cell = dof_handler.begin_active();

      const FEValuesExtractors::Vector displacement(0);

      for (; cell_dual!=endc_dual; ++cell_dual, ++cell)
        {
          cell_rhs = 0;
          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (cell_dual->face(face)->at_boundary()
                  &&
                  cell_dual->face(face)->boundary_id() == face_id)
                {
                  if (!evaluation_face_found)
                    {
                      evaluation_face_found = true;
                    }

                  fe_face_values.reinit (cell, face);
                  fe_face_values_dual.reinit (cell_dual, face);

                  fe_face_values[displacement].get_function_symmetric_gradients(solution,
                      strain_tensor);

                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      bound_size += fe_face_values_dual.JxW(q_point);

                      constitutive_law.get_stress_strain_tensor(strain_tensor[q_point],
                                                                stress_strain_tensor);

                      for (unsigned int i=0; i<dofs_per_cell_dual; ++i)
                        {
                          const SymmetricTensor<2, dim>
                          stress_phi_i = stress_strain_tensor
                                         * fe_face_values_dual[displacement].symmetric_gradient(i, q_point);

                          for (unsigned int k=0; k!=dim; ++k)
                            {
                              for (unsigned int l=0; l!=dim; ++l)
                                {
                                  if ( comp_stress[k][l] == 1 )
                                    {
                                      cell_rhs(i) += stress_phi_i[k][l]
                                                     *
                                                     fe_face_values_dual.JxW(q_point);
                                    }

                                }
                            }

                        }

                    }

                }
            }

          cell_dual->get_dof_indices (local_dof_indices);
          for (unsigned int i=0; i<dofs_per_cell_dual; ++i)
            {
              rhs_dual(local_dof_indices[i]) += cell_rhs(i);
            }

        }

      AssertThrow(evaluation_face_found, ExcInternalError());

      rhs_dual /= bound_size;

    }


    template <int dim>
    class MeanStressDomain : public DualFunctionalBase<dim>
    {
    public:
      MeanStressDomain (const std::string &base_mesh,
                        const std::vector<std::vector<unsigned int> > &comp_stress);

      virtual
      void
      assemble_rhs (const DoFHandler<dim>      &dof_handler,
                    const Vector<double>       &solution,
                    const ConstitutiveLaw<dim> &constitutive_law,
                    const DoFHandler<dim>      &dof_handler_dual,
                    Vector<double>             &rhs_dual) const;

    protected:
      const std::string base_mesh;
      const std::vector<std::vector<unsigned int> >  comp_stress;
    };


    template <int dim>
    MeanStressDomain<dim>::
    MeanStressDomain (const std::string &base_mesh,
                      const std::vector<std::vector<unsigned int> > &comp_stress )
      :
      base_mesh (base_mesh),
      comp_stress (comp_stress)
    {
      AssertThrow(comp_stress.size() == dim,
                  ExcDimensionMismatch (comp_stress.size(), dim) );
    }


    template <int dim>
    void
    MeanStressDomain<dim>::
    assemble_rhs (const DoFHandler<dim>      &dof_handler,
                  const Vector<double>       &solution,
                  const ConstitutiveLaw<dim> &constitutive_law,
                  const DoFHandler<dim>      &dof_handler_dual,
                  Vector<double>             &rhs_dual) const
    {
      AssertThrow (base_mesh == "Cantiliver_beam_3d", ExcNotImplemented());
      AssertThrow (dim == 3, ExcNotImplemented());

      // Mean stress at the specified domain is of interest.
      // The interest domains are located on the bottom and top of the flanges
      // close to the clamped face, z = 0
      // top domain: height/2 - thickness_flange <= y <= height/2
      //             0 <= z <= 2 * thickness_flange
      // bottom domain: -height/2 <= y <= -height/2 + thickness_flange
      //             0 <= z <= 2 * thickness_flange

      const double height = 200e-3,
                   thickness_flange = 10e-3;

      rhs_dual.reinit (dof_handler_dual.n_dofs());

      const QGauss<dim> quadrature_formula(dof_handler_dual.get_fe().tensor_degree()+1);

      FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula,
                               update_gradients);
      FEValues<dim> fe_values_dual (dof_handler_dual.get_fe(), quadrature_formula,
                                    update_gradients | update_JxW_values);

      const unsigned int  dofs_per_cell_dual = dof_handler_dual.get_fe().dofs_per_cell;
      const unsigned int  n_q_points = quadrature_formula.size();

      std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
      SymmetricTensor<4, dim> stress_strain_tensor;

      Vector<double>      cell_rhs (dofs_per_cell_dual);

      std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_dual);

      // domain_size : size of the interested domain, in 2d is the area
      //              and in the 3d case, volume
      double domain_size = 0.;

      bool evaluation_domain_found = false;

      typename DoFHandler<dim>::active_cell_iterator
      cell_dual = dof_handler_dual.begin_active(),
      endc_dual = dof_handler_dual.end(),
      cell = dof_handler.begin_active();

      const FEValuesExtractors::Vector displacement(0);

      for (; cell_dual!=endc_dual; ++cell_dual, ++cell)
        {
          const double y = cell->center()[1],
                       z = cell->center()[2];
          // top domain: height/2 - thickness_flange <= y <= height/2
          //             0 <= z <= 2 * thickness_flange
          // bottom domain: -height/2 <= y <= -height/2 + thickness_flange
          //             0 <= z <= 2 * thickness_flange
          if ( ((z > 0) && (z < 2*thickness_flange)) &&
               ( ((y > height/2 - thickness_flange) && (y < height/2)) ||
                 ((y > -height/2) && (y < -height/2 + thickness_flange)) ) )
            {
              cell_rhs = 0;

              if (!evaluation_domain_found)
                {
                  evaluation_domain_found = true;
                }

              fe_values.reinit(cell);
              fe_values_dual.reinit(cell_dual);

              fe_values[displacement].get_function_symmetric_gradients(solution,
                                                                       strain_tensor);

              for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
                {
                  domain_size += fe_values_dual.JxW(q_point);

                  constitutive_law.get_stress_strain_tensor(strain_tensor[q_point],
                                                            stress_strain_tensor);

                  for (unsigned int i=0; i<dofs_per_cell_dual; ++i)
                    {
                      const SymmetricTensor<2, dim>
                      stress_phi_i = stress_strain_tensor
                                     * fe_values_dual[displacement].symmetric_gradient(i, q_point);

                      for (unsigned int k=0; k!=dim; ++k)
                        {
                          for (unsigned int l=0; l!=dim; ++l)
                            {
                              if ( comp_stress[k][l] == 1 )
                                {
                                  cell_rhs(i) += stress_phi_i[k][l]
                                                 *
                                                 fe_values_dual.JxW(q_point);
                                }

                            }
                        }

                    }

                }

            }

          cell_dual->get_dof_indices (local_dof_indices);
          for (unsigned int i=0; i<dofs_per_cell_dual; ++i)
            {
              rhs_dual(local_dof_indices[i]) += cell_rhs(i);
            }

        }

      AssertThrow(evaluation_domain_found, ExcInternalError());

      rhs_dual /= domain_size;

    }


    template <int dim>
    class MeanStrainEnergyFace : public DualFunctionalBase<dim>
    {
    public:
      MeanStrainEnergyFace (const unsigned int face_id,
                            const Function<dim>      &lambda_function,
                            const Function<dim>      &mu_function );

      void assemble_rhs_nonlinear (const DoFHandler<dim> &primal_dof_handler,
                                   const Vector<double>  &primal_solution,
                                   const DoFHandler<dim> &dof_handler,
                                   Vector<double>        &rhs) const;

    protected:
      const unsigned int face_id;
      const SmartPointer<const Function<dim> >       lambda_function;
      const SmartPointer<const Function<dim> >       mu_function;
    };


    template <int dim>
    MeanStrainEnergyFace<dim>::
    MeanStrainEnergyFace (const unsigned int face_id,
                          const Function<dim>      &lambda_function,
                          const Function<dim>      &mu_function )
      :
      face_id (face_id),
      lambda_function (&lambda_function),
      mu_function (&mu_function)
    {}


    template <int dim>
    void
    MeanStrainEnergyFace<dim>::
    assemble_rhs_nonlinear (const DoFHandler<dim> &primal_dof_handler,
                            const Vector<double>  &primal_solution,
                            const DoFHandler<dim> &dof_handler,
                            Vector<double>        &rhs) const
    {
      // Assemble right hand side of the dual problem when the quantity of interest is
      // a nonlinear functinoal. In this case, the QoI should be linearized which depends
      // on the solution of the primal problem.
      // The extracter of the linearized QoI functional is the gradient of the the original
      // QoI functional with the primal solution values.

      AssertThrow (dim >= 2, ExcNotImplemented());

      rhs.reinit (dof_handler.n_dofs());

      const QGauss<dim-1> face_quadrature(dof_handler.get_fe().tensor_degree()+1);
      FEFaceValues<dim> primal_fe_face_values (primal_dof_handler.get_fe(), face_quadrature,
                                               update_quadrature_points |
                                               update_gradients | update_hessians |
                                               update_JxW_values);

      FEFaceValues<dim> fe_face_values (dof_handler.get_fe(), face_quadrature,
                                        update_values);

      const unsigned int  dofs_per_vertex = primal_dof_handler.get_fe().dofs_per_vertex;
      const unsigned int  n_face_q_points = face_quadrature.size();
      const unsigned int  dofs_per_cell = dof_handler.get_fe().dofs_per_cell;

      AssertThrow(dofs_per_vertex == dim,
                  ExcDimensionMismatch (dofs_per_vertex, dim) );

      std::vector< std::vector< Tensor<1,dim> > > primal_solution_gradients;
      primal_solution_gradients.resize(n_face_q_points);

      std::vector<std::vector<Tensor<2,dim> > >   primal_solution_hessians;
      primal_solution_hessians.resize (n_face_q_points);

      for (unsigned int i=0; i!=n_face_q_points; ++i)
        {
          primal_solution_gradients[i].resize (dofs_per_vertex);
          primal_solution_hessians[i].resize  (dofs_per_vertex);
        }

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

      Vector<double>      cell_rhs (dofs_per_cell);

      std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

      // bound_size : size of the boundary, in 2d is the length
      //              and in the 3d case, area
      double bound_size  = 0.;

      bool evaluation_face_found = false;

      typename DoFHandler<dim>::active_cell_iterator
      primal_cell = primal_dof_handler.begin_active(),
      primal_endc = primal_dof_handler.end();

      typename DoFHandler<dim>::active_cell_iterator
      cell = dof_handler.begin_active(),
      endc = dof_handler.end();

      for (; cell!=endc; ++cell, ++primal_cell)
        {
          cell_rhs = 0;
          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
            {
              if (cell->face(face)->at_boundary()
                  &&
                  cell->face(face)->boundary_id() == face_id)
                {
                  if (!evaluation_face_found)
                    {
                      evaluation_face_found = true;
                    }
                  primal_fe_face_values.reinit (primal_cell, face);

                  primal_fe_face_values.get_function_gradients (primal_solution,
                                                                primal_solution_gradients);

                  primal_fe_face_values.get_function_hessians (primal_solution,
                                                               primal_solution_hessians);

                  lambda_function->value_list (primal_fe_face_values.get_quadrature_points(), lambda_values);
                  mu_function->value_list     (primal_fe_face_values.get_quadrature_points(), mu_values);

                  fe_face_values.reinit (cell, face);

                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      bound_size += primal_fe_face_values.JxW(q_point);

                      for (unsigned int m=0; m<dofs_per_cell; ++m)
                        {
                          const unsigned int
                          component_m = dof_handler.get_fe().system_to_component_index(m).first;

                          for (unsigned int i=0; i!=dofs_per_vertex; ++i)
                            {
                              for (unsigned int j=0; j!=dofs_per_vertex; ++j)
                                {
                                  cell_rhs(m) += fe_face_values.shape_value(m,q_point) *
                                                 (
                                                   lambda_values[q_point] *
                                                   (
                                                     primal_solution_hessians[q_point][i][i][component_m] * primal_solution_gradients[q_point][j][j]
                                                     +
                                                     primal_solution_gradients[q_point][i][i] * primal_solution_hessians[q_point][j][j][component_m]
                                                   )
                                                   +
                                                   mu_values[q_point] *
                                                   (
                                                     2*primal_solution_hessians[q_point][j][i][component_m] * primal_solution_gradients[q_point][j][i]
                                                     +
                                                     primal_solution_hessians[q_point][i][j][component_m] * primal_solution_gradients[q_point][j][i]
                                                     +
                                                     primal_solution_gradients[q_point][i][j] * primal_solution_hessians[q_point][j][i][component_m]
                                                   )
                                                 ) *
                                                 primal_fe_face_values.JxW(q_point);

                                }
                            }

                        } // end loop DoFs


                    }  // end loop Gauss points

                }  // end if face
            }  // end loop face

          cell->get_dof_indices (local_dof_indices);
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            {
              rhs(local_dof_indices[i]) += cell_rhs(i);
            }

        }  // end loop cell

      AssertThrow(evaluation_face_found, ExcInternalError());

      rhs *= 1./(2*bound_size);

    }


  }


  // DualSolver class
  template <int dim>
  class DualSolver
  {
  public:
    DualSolver (const Triangulation<dim>                        &triangulation,
                const FESystem<dim>                             &fe,
                const Vector<double>                            &solution,
                const ConstitutiveLaw<dim>                      &constitutive_law,
                const DualFunctional::DualFunctionalBase<dim>   &dual_functional,
                const unsigned int                              &timestep_no,
                const std::string                               &output_dir,
                const std::string                               &base_mesh,
                const double                                    &present_time,
                const double                                    &end_time);

    void compute_error_DWR (Vector<float> &estimated_error_per_cell);

    ~DualSolver ();

  private:
    void setup_system ();
    void compute_dirichlet_constraints ();
    void assemble_matrix ();
    void assemble_rhs ();
    void solve ();
    void output_results ();

    const FESystem<dim>       &fe;
    DoFHandler<dim>            dof_handler;
    const Vector<double>       solution;

    const unsigned int         fe_degree;


    const unsigned int         fe_degree_dual;
    FESystem<dim>              fe_dual;
    DoFHandler<dim>            dof_handler_dual;

    const QGauss<dim>          quadrature_formula;
    const QGauss<dim - 1>      face_quadrature_formula;

    AffineConstraints<double>  constraints_hanging_nodes_dual;
    AffineConstraints<double>  constraints_dirichlet_and_hanging_nodes_dual;

    SparsityPattern            sparsity_pattern_dual;
    SparseMatrix<double>       system_matrix_dual;
    Vector<double>             system_rhs_dual;
    Vector<double>             solution_dual;

    const ConstitutiveLaw<dim> constitutive_law;

    const SmartPointer<const Triangulation<dim> > triangulation;
    const SmartPointer<const DualFunctional::DualFunctionalBase<dim> > dual_functional;

    unsigned int               timestep_no;
    std::string                output_dir;
    const std::string          base_mesh;
    double                     present_time;
    double                     end_time;
  };


  template<int dim>
  DualSolver<dim>::
  DualSolver (const Triangulation<dim>                        &triangulation,
              const FESystem<dim>                             &fe,
              const Vector<double>                            &solution,
              const ConstitutiveLaw<dim>                      &constitutive_law,
              const DualFunctional::DualFunctionalBase<dim>   &dual_functional,
              const unsigned int                              &timestep_no,
              const std::string                               &output_dir,
              const std::string                               &base_mesh,
              const double                                    &present_time,
              const double                                    &end_time)
    :
    fe (fe),
    dof_handler (triangulation),
    solution(solution),
    fe_degree(fe.tensor_degree()),
    fe_degree_dual(fe_degree + 1),
    fe_dual(FE_Q<dim>(fe_degree_dual), dim),
    dof_handler_dual (triangulation),
    quadrature_formula (fe_degree_dual + 1),
    face_quadrature_formula (fe_degree_dual + 1),
    constitutive_law (constitutive_law),
    triangulation (&triangulation),
    dual_functional (&dual_functional),
    timestep_no (timestep_no),
    output_dir (output_dir),
    base_mesh (base_mesh),
    present_time (present_time),
    end_time (end_time)
  {}


  template<int dim>
  DualSolver<dim>::~DualSolver()
  {
    dof_handler_dual.clear ();
  }


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

    dof_handler_dual.distribute_dofs (fe_dual);
    std::cout << "    Number of degrees of freedom in dual problem:  "
              << dof_handler_dual.n_dofs()
              << std::endl;

    constraints_hanging_nodes_dual.clear ();
    DoFTools::make_hanging_node_constraints (dof_handler_dual,
                                             constraints_hanging_nodes_dual);
    constraints_hanging_nodes_dual.close ();

    compute_dirichlet_constraints();

    sparsity_pattern_dual.reinit (dof_handler_dual.n_dofs(),
                                  dof_handler_dual.n_dofs(),
                                  dof_handler_dual.max_couplings_between_dofs());
    DoFTools::make_sparsity_pattern (dof_handler_dual, sparsity_pattern_dual);

//    constraints_hanging_nodes_dual.condense (sparsity_pattern_dual);
    constraints_dirichlet_and_hanging_nodes_dual.condense (sparsity_pattern_dual);

    sparsity_pattern_dual.compress();

    system_matrix_dual.reinit (sparsity_pattern_dual);

    solution_dual.reinit (dof_handler_dual.n_dofs());
    system_rhs_dual.reinit (dof_handler_dual.n_dofs());

  }

  template<int dim>
  void DualSolver<dim>::compute_dirichlet_constraints()
  {
    constraints_dirichlet_and_hanging_nodes_dual.clear ();
    constraints_dirichlet_and_hanging_nodes_dual.merge(constraints_hanging_nodes_dual);

    std::vector<bool> component_mask(dim);

    if (base_mesh == "Timoshenko beam")
      {
        VectorTools::interpolate_boundary_values(dof_handler_dual,
                                                 0,
                                                 EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                 constraints_dirichlet_and_hanging_nodes_dual,
                                                 ComponentMask());
      }
    else if (base_mesh == "Thick_tube_internal_pressure")
      {
        // the boundary x = 0
        component_mask[0] = true;
        component_mask[1] = false;
        VectorTools::interpolate_boundary_values (dof_handler_dual,
                                                  2,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes_dual,
                                                  component_mask);
        // the boundary y = 0
        component_mask[0] = false;
        component_mask[1] = true;
        VectorTools::interpolate_boundary_values (dof_handler_dual,
                                                  3,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes_dual,
                                                  component_mask);
      }
    else if (base_mesh == "Perforated_strip_tension")
      {
        // the boundary x = 0
        component_mask[0] = true;
        component_mask[1] = false;
        component_mask[2] = false;
        VectorTools::interpolate_boundary_values (dof_handler_dual,
                                                  4,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes_dual,
                                                  component_mask);
        // the boundary y = 0
        component_mask[0] = false;
        component_mask[1] = true;
        component_mask[2] = false;
        VectorTools::interpolate_boundary_values (dof_handler_dual,
                                                  1,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes_dual,
                                                  component_mask);
        // the boundary y = imposed incremental displacement
        component_mask[0] = false;
        component_mask[1] = true;
        component_mask[2] = false;
        VectorTools::interpolate_boundary_values (dof_handler_dual,
                                                  3,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes_dual,
                                                  component_mask);
      }
    else if (base_mesh == "Cantiliver_beam_3d")
      {
        // the boundary x = y = z = 0
        component_mask[0] = true;
        component_mask[1] = true;
        component_mask[2] = true;
        VectorTools::interpolate_boundary_values (dof_handler_dual,
                                                  1,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes_dual,
                                                  component_mask);
      }
    else
      {
        AssertThrow(false, ExcNotImplemented());
      }

    constraints_dirichlet_and_hanging_nodes_dual.close();
  }


  template<int dim>
  void DualSolver<dim>::assemble_matrix()
  {
    FEValues<dim> fe_values(fe, quadrature_formula, update_gradients);

    FEValues<dim> fe_values_dual(fe_dual, quadrature_formula,
                                 update_values | update_gradients | update_JxW_values);

    const unsigned int dofs_per_cell_dual = fe_dual.dofs_per_cell;
    const unsigned int n_q_points         = quadrature_formula.size();

    FullMatrix<double> cell_matrix (dofs_per_cell_dual, dofs_per_cell_dual);

    std::vector<types::global_dof_index>   local_dof_indices(dofs_per_cell_dual);

    typename DoFHandler<dim>::active_cell_iterator
    cell_dual = dof_handler_dual.begin_active(),
    endc_dual = dof_handler_dual.end(),
    cell = dof_handler.begin_active();

    const FEValuesExtractors::Vector displacement(0);

    for (; cell_dual != endc_dual; ++cell_dual, ++cell)
      if (cell_dual->is_locally_owned())
        {
          fe_values.reinit(cell);

          fe_values_dual.reinit(cell_dual);
          cell_matrix = 0;

          std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
          fe_values[displacement].get_function_symmetric_gradients(solution,
                                                                   strain_tensor);

          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
            {
              SymmetricTensor<4, dim> stress_strain_tensor_linearized;
              SymmetricTensor<4, dim> stress_strain_tensor;
              constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
                                                                    stress_strain_tensor_linearized,
                                                                    stress_strain_tensor);

              for (unsigned int i = 0; i < dofs_per_cell_dual; ++i)
                {
                  const SymmetricTensor<2, dim>
                  stress_phi_i = stress_strain_tensor_linearized
                                 * fe_values_dual[displacement].symmetric_gradient(i, q_point);

                  for (unsigned int j = 0; j < dofs_per_cell_dual; ++j)
                    cell_matrix(i, j) += (stress_phi_i
                                          * fe_values_dual[displacement].symmetric_gradient(j, q_point)
                                          * fe_values_dual.JxW(q_point));

                }

            }

          cell_dual->get_dof_indices(local_dof_indices);
          constraints_dirichlet_and_hanging_nodes_dual.distribute_local_to_global(cell_matrix,
              local_dof_indices,
              system_matrix_dual);

        }

  }


  template<int dim>
  void DualSolver<dim>::assemble_rhs()
  {
    dual_functional->assemble_rhs (dof_handler, solution, constitutive_law,
                                   dof_handler_dual, system_rhs_dual);
    constraints_dirichlet_and_hanging_nodes_dual.condense (system_rhs_dual);
  }


  template<int dim>
  void DualSolver<dim>::solve()
  {
    // +++  direct solver +++++++++
    SparseDirectUMFPACK   A_direct;
    A_direct.initialize(system_matrix_dual);

    // After the decomposition, we can use A_direct like a matrix representing
    // the inverse of our system matrix, so to compute the solution we just
    // have to multiply with the right hand side vector:
    A_direct.vmult(solution_dual, system_rhs_dual);

    // ++++  iterative solver ++ CG ++++ doesn't work
//    SolverControl solver_control (5000, 1e-12);
//    SolverCG<> cg (solver_control);
//
//    PreconditionSSOR<> preconditioner;
//    preconditioner.initialize(system_matrix_dual, 1.2);
//
//    cg.solve (system_matrix_dual, solution_dual, system_rhs_dual,
//              preconditioner);

    // ++++  iterative solver ++ BiCGStab ++++++ doesn't work
//    SolverControl solver_control (5000, 1e-12);
//    SolverBicgstab<> bicgstab (solver_control);
//
//    PreconditionJacobi<> preconditioner;
//    preconditioner.initialize(system_matrix_dual, 1.0);
//
//    bicgstab.solve (system_matrix_dual, solution_dual, system_rhs_dual,
//                    preconditioner);

    // +++++++++++++++++++++++++++++++++++++++++++++++++

    constraints_dirichlet_and_hanging_nodes_dual.distribute (solution_dual);
  }

  template<int dim>
  void DualSolver<dim>::output_results()
  {
    std::string filename = (output_dir + "dual-solution-" +
                            Utilities::int_to_string(timestep_no, 4) + ".vtk");
    std::ofstream output (filename.c_str());
    DataOut<dim> data_out;
    data_out.attach_dof_handler (dof_handler_dual);
    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.add_data_vector (solution_dual, solution_names);
    data_out.build_patches ();
    data_out.write_vtk (output);
  }

  template<int dim>
  void DualSolver<dim>::compute_error_DWR (Vector<float> &estimated_error_per_cell)
  {
    Assert (estimated_error_per_cell.size() == triangulation->n_global_active_cells(),
            ExcDimensionMismatch (estimated_error_per_cell.size(), triangulation->n_global_active_cells()));

    // solve the dual problem
    setup_system ();
    assemble_matrix ();
    assemble_rhs ();
    solve ();
    output_results ();

    // compuate the dual weights
    Vector<double> primal_solution (dof_handler_dual.n_dofs());
    FETools::interpolate (dof_handler,
                          solution,
                          dof_handler_dual,
                          constraints_dirichlet_and_hanging_nodes_dual,
                          primal_solution);

    AffineConstraints<double> constraints_hanging_nodes;
    DoFTools::make_hanging_node_constraints (dof_handler,
                                             constraints_hanging_nodes);
    constraints_hanging_nodes.close();
    Vector<double> dual_weights (dof_handler_dual.n_dofs());
    FETools::interpolation_difference (dof_handler_dual,
                                       constraints_dirichlet_and_hanging_nodes_dual,
                                       solution_dual,
                                       dof_handler,
                                       constraints_hanging_nodes,
                                       dual_weights);

    // estimate the error
    FEValues<dim> fe_values(fe_dual, quadrature_formula,
                            update_values    |
                            update_gradients |
                            update_hessians  |
                            update_quadrature_points |
                            update_JxW_values);

    const unsigned int n_q_points      = quadrature_formula.size();
    std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
    SymmetricTensor<4, dim> stress_strain_tensor_linearized;
    SymmetricTensor<4, dim> stress_strain_tensor;
    Tensor<5, dim>          stress_strain_tensor_grad;
    std::vector<std::vector<Tensor<2,dim> > > cell_hessians (n_q_points);
    for (unsigned int i=0; i!=n_q_points; ++i)
      {
        cell_hessians[i].resize (dim);
      }
    std::vector<Vector<double> > dual_weights_cell_values (n_q_points, Vector<double>(dim));

    const EquationData::BodyForce<dim> body_force;
    std::vector<Vector<double> > body_force_values (n_q_points, Vector<double>(dim));
    const FEValuesExtractors::Vector displacement(0);


    FEFaceValues<dim> fe_face_values_cell(fe_dual, face_quadrature_formula,
                                          update_values           |
                                          update_quadrature_points|
                                          update_gradients        |
                                          update_JxW_values       |
                                          update_normal_vectors),
                                          fe_face_values_neighbor (fe_dual, face_quadrature_formula,
                                              update_values     |
                                              update_gradients  |
                                              update_JxW_values |
                                              update_normal_vectors);
    FESubfaceValues<dim> fe_subface_values_cell (fe_dual, face_quadrature_formula,
                                                 update_gradients);

    const unsigned int n_face_q_points = face_quadrature_formula.size();
    std::vector<Vector<double> > jump_residual (n_face_q_points, Vector<double>(dim));
    std::vector<Vector<double> > dual_weights_face_values (n_face_q_points, Vector<double>(dim));

    std::vector<std::vector<Tensor<1,dim> > > cell_grads(n_face_q_points);
    for (unsigned int i=0; i!=n_face_q_points; ++i)
      {
        cell_grads[i].resize (dim);
      }
    std::vector<std::vector<Tensor<1,dim> > > neighbor_grads(n_face_q_points);
    for (unsigned int i=0; i!=n_face_q_points; ++i)
      {
        neighbor_grads[i].resize (dim);
      }
    SymmetricTensor<2, dim> q_cell_strain_tensor;
    SymmetricTensor<2, dim> q_neighbor_strain_tensor;
    SymmetricTensor<4, dim> cell_stress_strain_tensor;
    SymmetricTensor<4, dim> neighbor_stress_strain_tensor;


    typename std::map<typename DoFHandler<dim>::face_iterator, Vector<double> >
    face_integrals;
    typename DoFHandler<dim>::active_cell_iterator
    cell = dof_handler_dual.begin_active(),
    endc = dof_handler_dual.end();
    for (; cell!=endc; ++cell)
      if (cell->is_locally_owned())
        {
          for (unsigned int face_no=0;
               face_no<GeometryInfo<dim>::faces_per_cell;
               ++face_no)
            {
              face_integrals[cell->face(face_no)].reinit (dim);
              face_integrals[cell->face(face_no)] = -1e20;
            }
        }

    std::vector<Vector<float> > error_indicators_vector;
    error_indicators_vector.resize( triangulation->n_active_cells(),
                                    Vector<float>(dim) );

    // ----------------- estimate_some -------------------------
    cell = dof_handler_dual.begin_active();
    unsigned int present_cell = 0;
    for (; cell!=endc; ++cell, ++present_cell)
      if (cell->is_locally_owned())
        {
          // --------------- integrate_over_cell -------------------
          fe_values.reinit(cell);
          body_force.vector_value_list(fe_values.get_quadrature_points(),
                                       body_force_values);
          fe_values[displacement].get_function_symmetric_gradients(primal_solution,
                                                                   strain_tensor);
          fe_values.get_function_hessians(primal_solution, cell_hessians);

          fe_values.get_function_values(dual_weights,
                                        dual_weights_cell_values);

          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
            {
              constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
                                                                    stress_strain_tensor_linearized,
                                                                    stress_strain_tensor);
              constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point],
                                                             cell_hessians[q_point],
                                                             stress_strain_tensor_grad);

              for (unsigned int i=0; i!=dim; ++i)
                {
                  error_indicators_vector[present_cell](i) +=
                    body_force_values[q_point](i)*
                    dual_weights_cell_values[q_point](i)*
                    fe_values.JxW(q_point);
                  for (unsigned int j=0; j!=dim; ++j)
                    {
                      for (unsigned int k=0; k!=dim; ++k)
                        {
                          for (unsigned int l=0; l!=dim; ++l)
                            {
                              error_indicators_vector[present_cell](i) +=
                                ( stress_strain_tensor[i][j][k][l]*
                                  0.5*(cell_hessians[q_point][k][l][j]
                                       +
                                       cell_hessians[q_point][l][k][j])
                                  + stress_strain_tensor_grad[i][j][k][l][j] * strain_tensor[q_point][k][l]
                                ) *
                                dual_weights_cell_values[q_point](i) *
                                fe_values.JxW(q_point);
                            }
                        }
                    }

                }

            }
          // -------------------------------------------------------
          // compute face_integrals
          for (unsigned int face_no=0;
               face_no<GeometryInfo<dim>::faces_per_cell;
               ++face_no)
            {
              if (cell->face(face_no)->at_boundary())
                {
                  for (unsigned int id=0; id!=dim; ++id)
                    {
                      face_integrals[cell->face(face_no)](id) = 0;
                    }
                  continue;
                }

              if ((cell->neighbor(face_no)->has_children() == false) &&
                  (cell->neighbor(face_no)->level() == cell->level()) &&
                  (cell->neighbor(face_no)->index() < cell->index()))
                continue;

              if (cell->at_boundary(face_no) == false)
                if (cell->neighbor(face_no)->level() < cell->level())
                  continue;


              if (cell->face(face_no)->has_children() == false)
                {
                  // ------------- integrate_over_regular_face -----------
                  fe_face_values_cell.reinit(cell, face_no);
                  fe_face_values_cell.get_function_gradients (primal_solution,
                                                              cell_grads);

                  Assert (cell->neighbor(face_no).state() == IteratorState::valid,
                          ExcInternalError());
                  const unsigned int
                  neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
                  const typename DoFHandler<dim>::active_cell_iterator
                  neighbor = cell->neighbor(face_no);

                  fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
                  fe_face_values_neighbor.get_function_gradients (primal_solution,
                                                                  neighbor_grads);

                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      q_cell_strain_tensor = 0.;
                      q_neighbor_strain_tensor = 0.;
                      for (unsigned int i=0; i!=dim; ++i)
                        {
                          for (unsigned int j=0; j!=dim; ++j)
                            {
                              q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
                                                                cell_grads[q_point][j][i] );
                              q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
                                                                    neighbor_grads[q_point][j][i] );
                            }
                        }

                      constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
                                                                 cell_stress_strain_tensor);
                      constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
                                                                 neighbor_stress_strain_tensor);

                      jump_residual[q_point] = 0.;
                      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)
                                    {
                                      jump_residual[q_point](i) += (cell_stress_strain_tensor[i][j][k][l]*
                                                                    q_cell_strain_tensor[k][l]
                                                                    -
                                                                    neighbor_stress_strain_tensor[i][j][k][l]*
                                                                    q_neighbor_strain_tensor[k][l] )*
                                                                   fe_face_values_cell.normal_vector(q_point)[j];
                                    }
                                }
                            }
                        }

                    }

                  fe_face_values_cell.get_function_values (dual_weights,
                                                           dual_weights_face_values);

                  Vector<double> face_integral_vector(dim);
                  face_integral_vector = 0;
                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      for (unsigned int i=0; i!=dim; ++i)
                        {
                          face_integral_vector(i) += jump_residual[q_point](i) *
                                                     dual_weights_face_values[q_point](i) *
                                                     fe_face_values_cell.JxW(q_point);
                        }
                    }

                  Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
                          ExcInternalError());

                  for (unsigned int i=0; i!=dim; ++i)
                    {
                      Assert (face_integrals[cell->face(face_no)](i) == -1e20,
                              ExcInternalError());
                      face_integrals[cell->face(face_no)](i) = face_integral_vector(i);

                    }

                  // -----------------------------------------------------
                }
              else
                {
                  // ------------- integrate_over_irregular_face ---------
                  const typename DoFHandler<dim>::face_iterator
                  face = cell->face(face_no);
                  const typename DoFHandler<dim>::cell_iterator
                  neighbor = cell->neighbor(face_no);
                  Assert (neighbor.state() == IteratorState::valid,
                          ExcInternalError());
                  Assert (neighbor->has_children(),
                          ExcInternalError());

                  const unsigned int
                  neighbor_neighbor = cell->neighbor_of_neighbor (face_no);

                  for (unsigned int subface_no=0;
                       subface_no<face->n_children(); ++subface_no)
                    {
                      const typename DoFHandler<dim>::active_cell_iterator
                      neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no);
                      Assert (neighbor_child->face(neighbor_neighbor) ==
                              cell->face(face_no)->child(subface_no),
                              ExcInternalError());

                      fe_subface_values_cell.reinit (cell, face_no, subface_no);
                      fe_subface_values_cell.get_function_gradients (primal_solution,
                                                                     cell_grads);
                      fe_face_values_neighbor.reinit (neighbor_child,
                                                      neighbor_neighbor);
                      fe_face_values_neighbor.get_function_gradients (primal_solution,
                                                                      neighbor_grads);

                      for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                        {
                          q_cell_strain_tensor = 0.;
                          q_neighbor_strain_tensor = 0.;
                          for (unsigned int i=0; i!=dim; ++i)
                            {
                              for (unsigned int j=0; j!=dim; ++j)
                                {
                                  q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
                                                                    cell_grads[q_point][j][i] );
                                  q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
                                                                        neighbor_grads[q_point][j][i] );
                                }
                            }

                          constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
                                                                     cell_stress_strain_tensor);
                          constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
                                                                     neighbor_stress_strain_tensor);

                          jump_residual[q_point] = 0.;
                          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)
                                        {
                                          jump_residual[q_point](i) += (-cell_stress_strain_tensor[i][j][k][l]*
                                                                        q_cell_strain_tensor[k][l]
                                                                        +
                                                                        neighbor_stress_strain_tensor[i][j][k][l]*
                                                                        q_neighbor_strain_tensor[k][l] )*
                                                                       fe_face_values_neighbor.normal_vector(q_point)[j];
                                        }
                                    }
                                }
                            }

                        }

                      fe_face_values_neighbor.get_function_values (dual_weights,
                                                                   dual_weights_face_values);

                      Vector<double> face_integral_vector(dim);
                      face_integral_vector = 0;
                      for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                        {
                          for (unsigned int i=0; i!=dim; ++i)
                            {
                              face_integral_vector(i) += jump_residual[q_point](i) *
                                                         dual_weights_face_values[q_point](i) *
                                                         fe_face_values_neighbor.JxW(q_point);
                            }
                        }

                      for (unsigned int i=0; i!=dim; ++i)
                        {
                          face_integrals[neighbor_child->face(neighbor_neighbor)](i) = face_integral_vector(i);
                        }

                    }

                  Vector<double> sum (dim);
                  sum = 0;
                  for (unsigned int subface_no=0;
                       subface_no<face->n_children(); ++subface_no)
                    {
                      Assert (face_integrals.find(face->child(subface_no)) !=
                              face_integrals.end(),
                              ExcInternalError());
                      for (unsigned int i=0; i!=dim; ++i)
                        {
                          Assert (face_integrals[face->child(subface_no)](i) != -1e20,
                                  ExcInternalError());
                          sum(i) += face_integrals[face->child(subface_no)](i);
                        }
                    }
                  for (unsigned int i=0; i!=dim; ++i)
                    {
                      face_integrals[face](i) = sum(i);
                    }


                  // -----------------------------------------------------
                }


            }
        }
    // ----------------------------------------------------------

    present_cell=0;
    cell = dof_handler_dual.begin_active();
    for (; cell!=endc; ++cell, ++present_cell)
      if (cell->is_locally_owned())
        {
          for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
               ++face_no)
            {
              Assert(face_integrals.find(cell->face(face_no)) !=
                     face_integrals.end(),
                     ExcInternalError());

              for (unsigned int id=0; id!=dim; ++id)
                {
                  error_indicators_vector[present_cell](id)
                  -= 0.5*face_integrals[cell->face(face_no)](id);
                }

            }

          estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm();

        }
  }



  // @sect3{The <code>PlasticityContactProblem</code> class template}

  // This is the main class of this program and supplies all functions
  // and variables needed to describe
  // the nonlinear contact problem. It is
  // close to step-41 but with some additional
  // features like handling hanging nodes,
  // a Newton method, using Trilinos and p4est
  // for parallel distributed computing.
  // To deal with hanging nodes makes
  // life a bit more complicated since
  // we need another AffineConstraints object now.
  // We create a Newton method for the
  // active set method for the contact
  // situation and to handle the nonlinear
  // operator for the constitutive law.
  //
  // The general layout of this class is very much like for most other tutorial programs.
  // To make our life a bit easier, this class reads a set of input parameters from an input file. These
  // parameters, using the ParameterHandler class, are declared in the <code>declare_parameters</code>
  // function (which is static so that it can be called before we even create an object of the current
  // type), and a ParameterHandler object that has been used to read an input file will then be passed
  // to the constructor of this class.
  //
  // The remaining member functions are by and large as we have seen in several of the other tutorial
  // programs, though with additions for the current nonlinear system. We will comment on their purpose
  // as we get to them further below.
  template <int dim>
  class ElastoPlasticProblem
  {
  public:
    ElastoPlasticProblem (const ParameterHandler &prm);

    void run ();

    static void declare_parameters (ParameterHandler &prm);

  private:
    void make_grid ();
    void setup_system ();
    void compute_dirichlet_constraints ();
    void assemble_newton_system (const TrilinosWrappers::MPI::Vector &linearization_point,
                                 const TrilinosWrappers::MPI::Vector &delta_linearization_point);
    void compute_nonlinear_residual (const TrilinosWrappers::MPI::Vector &linearization_point);
    void solve_newton_system ();
    void solve_newton ();
    void compute_error ();
    void compute_error_residual (const TrilinosWrappers::MPI::Vector &tmp_solution);
    void refine_grid ();
    void move_mesh (const TrilinosWrappers::MPI::Vector &displacement) const;
    void output_results (const std::string &filename_base);

    // Next are three 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 ();

    // As far as member variables are concerned, we start with ones that we use to
    // indicate the MPI universe this program runs on, and then two numbers
    // telling us how many participating processors there are, and where in
    // this world we are., a stream we use to let
    // exactly one processor produce output to the console (see step-17) and
    // a variable that is used to time the various sections of the program:
    MPI_Comm           mpi_communicator;
    const unsigned int n_mpi_processes;
    const unsigned int this_mpi_process;
    ConditionalOStream pcout;
    TimerOutput        computing_timer;

    // The next group describes the mesh and the finite element space.
    // In particular, for this parallel program, the finite element
    // space has associated with it variables that indicate which degrees
    // of freedom live on the current processor (the index sets, see
    // also step-40 and the @ref distributed documentation module) as
    // well as a variety of constraints: those imposed by hanging nodes,
    // by Dirichlet boundary conditions, and by the active set of
    // contact nodes. Of the three AffineConstraints objects defined
    // here, the first only contains hanging node constraints, the
    // second also those associated with Dirichlet boundary conditions,
    // and the third these plus the contact constraints.
    //
    // The variable <code>active_set</code> consists of those degrees
    // of freedom constrained by the contact, and we use
    // <code>fraction_of_plastic_q_points_per_cell</code> to keep
    // track of the fraction of quadrature points on each cell where
    // the stress equals the yield stress. The latter is only used to
    // create graphical output showing the plastic zone, but not for
    // any further computation; the variable is a member variable of
    // this class since the information is computed as a by-product
    // of computing the residual, but is used only much later. (Note
    // that the vector is a vector of length equal to the number of
    // active cells on the <i>local mesh</i>; it is never used to
    // exchange information between processors and can therefore be
    // a regular deal.II vector.)
    const unsigned int                        n_initial_global_refinements;
    parallel::distributed::Triangulation<dim> triangulation;

    const unsigned int fe_degree;
    FESystem<dim>      fe;
    DoFHandler<dim>    dof_handler;

    IndexSet           locally_owned_dofs;
    IndexSet           locally_relevant_dofs;

    AffineConstraints<double> constraints_hanging_nodes;
    AffineConstraints<double> constraints_dirichlet_and_hanging_nodes;

    Vector<float>      fraction_of_plastic_q_points_per_cell;

    // 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 right 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;
    const QGauss<dim - 1>      face_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 <code>user pointer</code>
    // that each cell, face, or edge holds: it is a <code>void*</code> 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.


    // The next block of variables corresponds to the solution
    // and the linear systems we need to form. In particular, this
    // includes the Newton matrix and right hand side; the vector
    // that corresponds to the residual (i.e., the Newton right hand
    // side) but from which we have not eliminated the various
    // constraints and that is used to determine which degrees of
    // freedom need to be constrained in the next iteration; and
    // a vector that corresponds to the diagonal of the $B$ matrix
    // briefly mentioned in the introduction and discussed in the
    // accompanying paper.
    TrilinosWrappers::SparseMatrix    newton_matrix;

    TrilinosWrappers::MPI::Vector     solution;
    TrilinosWrappers::MPI::Vector     incremental_displacement;
    TrilinosWrappers::MPI::Vector     newton_rhs;
    TrilinosWrappers::MPI::Vector     newton_rhs_residual;

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

    // The next block contains the variables that describe the material
    // response:
    const double         e_modulus, nu, sigma_0, gamma;
    ConstitutiveLaw<dim> constitutive_law;

    // And then there is an assortment of other variables that are used
    // to identify the mesh we are asked to build as selected by the
    // parameter file, the obstacle that is being pushed into the
    // deformable body, the mesh refinement strategy, whether to transfer
    // the solution from one mesh to the next, and how many mesh
    // refinement cycles to perform. As possible, we mark these kinds
    // of variables as <code>const</code> to help the reader identify
    // which ones may or may not be modified later on (the output directory
    // being an exception -- it is never modified outside the constructor
    // but it is awkward to initialize in the member-initializer-list
    // following the colon in the constructor since there we have only
    // one shot at setting it; the same is true for the mesh refinement
    // criterion):
    const std::string                                  base_mesh;

    struct RefinementStrategy
    {
      enum value
      {
        refine_global,
        refine_percentage,
        refine_fix_dofs
      };
    };
    typename RefinementStrategy::value                 refinement_strategy;

    struct ErrorEstimationStrategy
    {
      enum value
      {
        kelly_error,
        residual_error,
        weighted_residual_error,
        weighted_kelly_error
      };
    };
    typename ErrorEstimationStrategy::value            error_estimation_strategy;

    Vector<float>                                      estimated_error_per_cell;

    const bool                                         transfer_solution;
    std::string                                        output_dir;
    TableHandler                                       table_results,
                                                       table_results_2,
                                                       table_results_3;

    unsigned int                                       current_refinement_cycle;

    const double                                       max_relative_error;
    float                                              relative_error;

    const bool                                         show_stresses;
  };


  // @sect3{Implementation of the <code>PlasticityContactProblem</code> class}

  // @sect4{PlasticityContactProblem::declare_parameters}

  // Let us start with the declaration of run-time parameters that can be
  // selected in the input file. These values will be read back in the
  // constructor of this class to initialize the member variables of this
  // class:
  template <int dim>
  void
  ElastoPlasticProblem<dim>::declare_parameters (ParameterHandler &prm)
  {
    prm.declare_entry("polynomial degree", "1",
                      Patterns::Integer(),
                      "Polynomial degree of the FE_Q finite element space, typically 1 or 2.");
    prm.declare_entry("number of initial refinements", "2",
                      Patterns::Integer(),
                      "Number of initial global mesh refinement steps before "
                      "the first computation.");
    prm.declare_entry("refinement strategy", "percentage",
                      Patterns::Selection("global|percentage"),
                      "Mesh refinement strategy:\n"
                      " global: one global refinement\n"
                      " percentage: a fixed percentage of cells gets refined using the selected error estimator.");
    prm.declare_entry("error estimation strategy", "kelly_error",
                      Patterns::Selection("kelly_error|residual_error|weighted_residual_error"),
                      "Error estimation strategy:\n"
                      " kelly_error: Kelly error estimator\n"
                      " residual_error: residual-based error estimator\n"
                      " weighted_residual_error: dual weighted residual (Goal-oriented) error estimator.\n");
    prm.declare_entry("maximum relative error","0.05",
                      Patterns::Double(),
                      "maximum relative error which plays the role of a criteria for refinement.");
    prm.declare_entry("number of cycles", "5",
                      Patterns::Integer(),
                      "Number of adaptive mesh refinement cycles to run.");
    prm.declare_entry("output directory", "",
                      Patterns::Anything(),
                      "Directory for output files (graphical output and benchmark "
                      "statistics). If empty, use the current directory.");
    prm.declare_entry("transfer solution", "true",
                      Patterns::Bool(),
                      "Whether the solution should be used as a starting guess "
                      "for the next finer mesh. If false, then the iteration starts at "
                      "zero on every mesh.");
    prm.declare_entry("base mesh", "Thick_tube_internal_pressure",
                      Patterns::Selection("Timoshenko beam|Thick_tube_internal_pressure|"
                                          "Perforated_strip_tension|Cantiliver_beam_3d"),
                      "Select the shape of the domain: 'box' or 'half sphere'");
    prm.declare_entry("elasticity modulus","2.e11",
                      Patterns::Double(),
                      "Elasticity modulus of the material in MPa (N/mm2)");
    prm.declare_entry("Poissons ratio","0.3",
                      Patterns::Double(),
                      "Poisson's ratio of the material");
    prm.declare_entry("yield stress","2.e11",
                      Patterns::Double(),
                      "Yield stress of the material in MPa (N/mm2)");
    prm.declare_entry("isotropic hardening parameter","0.",
                      Patterns::Double(),
                      "Isotropic hardening parameter of the material");
    prm.declare_entry("show stresses", "false",
                      Patterns::Bool(),
                      "Whether illustrates the stresses and von Mises stresses or not.");


  }


  // @sect4{The <code>PlasticityContactProblem</code> constructor}

  // Given the declarations of member variables as well as the
  // declarations of run-time parameters that are read from the input
  // file, there is nothing surprising in this constructor. In the body
  // we initialize the mesh refinement strategy and the output directory,
  // creating such a directory if necessary.
  template <int dim>
  ElastoPlasticProblem<dim>::
  ElastoPlasticProblem (const ParameterHandler &prm)
    :
    mpi_communicator(MPI_COMM_WORLD),
    n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)),
    this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),
    pcout(std::cout, this_mpi_process == 0),
    computing_timer(MPI_COMM_WORLD, pcout, TimerOutput::never,
                    TimerOutput::wall_times),

    n_initial_global_refinements (prm.get_integer("number of initial refinements")),
    triangulation(mpi_communicator),
    fe_degree (prm.get_integer("polynomial degree")),
    fe(FE_Q<dim>(QGaussLobatto<1>(fe_degree+1)), dim),
    dof_handler(triangulation),
    quadrature_formula (fe_degree + 1),
    face_quadrature_formula (fe_degree + 1),

    e_modulus (prm.get_double("elasticity modulus")),
    nu (prm.get_double("Poissons ratio")),
    sigma_0(prm.get_double("yield stress")),
    gamma (prm.get_double("isotropic hardening parameter")),
    constitutive_law (e_modulus,
                      nu,
                      sigma_0,
                      gamma),

    base_mesh (prm.get("base mesh")),

    transfer_solution (prm.get_bool("transfer solution")),
    table_results(),
    table_results_2(),
    table_results_3(),
    max_relative_error (prm.get_double("maximum relative error")),
    show_stresses (prm.get_bool("show stresses"))
  {
    std::string strat = prm.get("refinement strategy");
    if (strat == "global")
      refinement_strategy = RefinementStrategy::refine_global;
    else if (strat == "percentage")
      refinement_strategy = RefinementStrategy::refine_percentage;
    else
      AssertThrow (false, ExcNotImplemented());

    strat = prm.get("error estimation strategy");
    if (strat == "kelly_error")
      error_estimation_strategy = ErrorEstimationStrategy::kelly_error;
    else if (strat == "residual_error")
      error_estimation_strategy = ErrorEstimationStrategy::residual_error;
    else if (strat == "weighted_residual_error")
      error_estimation_strategy = ErrorEstimationStrategy::weighted_residual_error;
    else
      AssertThrow(false, ExcNotImplemented());

    output_dir = prm.get("output directory");
    if (output_dir != "" && *(output_dir.rbegin()) != '/')
      output_dir += "/";
    mkdir(output_dir.c_str(), 0777);

    pcout << "    Using output directory '" << output_dir << "'" << std::endl;
    pcout << "    FE degree " << fe_degree << std::endl;
    pcout << "    transfer solution "
          << (transfer_solution ? "true" : "false") << std::endl;
  }



  // @sect4{PlasticityContactProblem::make_grid}

  // The next block deals with constructing the starting mesh.
  // We will use the following helper function and the first
  // block of the <code>make_grid()</code> to construct a
  // mesh that corresponds to a half sphere. deal.II has a function
  // that creates such a mesh, but it is in the wrong location
  // and facing the wrong direction, so we need to shift and rotate
  // it a bit before using it.
  //
  // For later reference, as described in the documentation of
  // GridGenerator::half_hyper_ball(), the flat surface of the halfsphere
  // has boundary indicator zero, while the remainder has boundary
  // indicator one.
  Point<3>
  rotate_half_sphere (const Point<3> &in)
  {
    return Point<3>(in(2), in(1), -in(0));
  }

  template <int dim>
  void
  ElastoPlasticProblem<dim>::make_grid ()
  {
    if (base_mesh == "Timoshenko beam")
      {
        AssertThrow (dim == 2, ExcNotImplemented());

        const double length = .48,
                     depth  = .12;

        const Point<dim> point_1(0, -depth/2),
              point_2(length, depth/2);

        std::vector<unsigned int> repetitions(2);
        repetitions[0] = 4;
        repetitions[1] = 1;
        GridGenerator::subdivided_hyper_rectangle(triangulation, repetitions, point_1, point_2);


        // give the indicators to boundaries for specification,
        //
        //     ________100______
        //     |                |
        //   0 |                | 5
        //     |________________|
        //             100
        // 0 to essential boundary conditions (left edge) which are as default
        // 100 to the null boundaries (upper and lower edges) where we do not need to take care of them
        // 5 to the natural boundaries (right edge) for imposing the traction force
        typename Triangulation<dim>::cell_iterator
        cell = triangulation.begin_active(),
        endc = triangulation.end();
        for (; cell!=endc; ++cell)
          {
            for (unsigned int face=0; face!=GeometryInfo<dim>::faces_per_cell; ++face)
              {
                if ( std::fabs(cell->face(face)->center()(0)-length) < 1e-12 )
                  {
                    cell->face(face)->set_manifold_id(5);
                  }
                else if ( ( std::fabs(cell->face(face)->center()(1)-(depth/2)) < 1e-12 )
                          ||
                          ( std::fabs(cell->face(face)->center()(1)-(-depth/2)) < 1e-12 ) )
                  {
                    cell->face(face)->set_manifold_id(100);
                  }

              }
          }

        triangulation.refine_global(n_initial_global_refinements);

      }
    else if (base_mesh == "Thick_tube_internal_pressure")
      {
        // Example 1 from the paper: Zhong Z., .... A new numerical method for determining
        // collapse load-carrying capacity of structure made of elasto-plastic material,
        // J. Cent. South Univ. (2014) 21: 398-404
        AssertThrow (dim == 2, ExcNotImplemented());

        const Point<dim> center(0, 0);
        const double inner_radius = .1,
                     outer_radius = .2;
        GridGenerator::quarter_hyper_shell(triangulation,
                                           center, inner_radius, outer_radius,
                                           0, true);

        // give the indicators to boundaries for specification,

        /*    _____
             |     \
             |       \
           2 |         \ 1
             |_          \
               \          \
              0 \         |
                 |________|
                     3
        */
        // 0 - inner boundary  - natural boundary condition - impose the traction force
        // 1 - outer boundary  - free boundary - we do not need to take care of them
        // 2 - left boundary   - essential boundary condition - constrained to move along the x direction
        // 3 - bottom boundary - essential boundary condition - constrained to move along the y direction

        const SphericalManifold<dim> inner_boundary_description(center);
        triangulation.set_manifold (0, inner_boundary_description);

        const SphericalManifold<dim> outer_boundary_description(center);
        triangulation.set_manifold (1, outer_boundary_description);

        triangulation.refine_global(n_initial_global_refinements);

        triangulation.reset_manifold (0);
        triangulation.reset_manifold (1);

      }
    else if (base_mesh == "Perforated_strip_tension")
      {
        // Example 2 from the paper: Zhong Z., .... A new numerical method for determining
        // collapse load-carrying capacity of structure made of elasto-plastic material,
        // J. Cent. South Univ. (2014) 21: 398-404
        AssertThrow (dim == 3, ExcNotImplemented());

        const int dim_2d = 2;
        const Point<dim_2d> center_2d(0, 0);
        const double inner_radius = 0.05,
                     outer_radius = 0.1,
                     height = 0.18,
                     thickness = 0.004;
//                   thickness = 0.01;

        Triangulation<dim_2d> triangulation_1,
                      triangulation_2,
                      triangulation_2d;

        const double eps = 1e-7 * inner_radius;
        {
          Point<dim_2d> point;

          GridGenerator::quarter_hyper_shell(triangulation_1,
                                             center_2d, inner_radius, outer_radius,
                                             2);

          // Modify the triangulation_1
          typename Triangulation<dim_2d>::active_cell_iterator
          cell = triangulation_1.begin_active(),
          endc = triangulation_1.end();
          std::vector<bool> treated_vertices(triangulation_1.n_vertices(), false);
          for (; cell != endc; ++cell)
            {
              for (unsigned int f=0; f<GeometryInfo<dim_2d>::faces_per_cell; ++f)
                if (cell->face(f)->at_boundary() && cell->face(f)->center()(0)>eps &&
                    cell->face(f)->center()(1)>eps )
                  {
                    // distance of the face center from the center
                    point(0) = cell->face(f)->center()(0) - center_2d(0);
                    point(1) = cell->face(f)->center()(1) - center_2d(1);
                    if ( point.norm() > (inner_radius + eps) )
                      {
                        for (unsigned int v=0; v < GeometryInfo<dim_2d>::vertices_per_face; ++v)
                          {
                            unsigned int vv = cell->face(f)->vertex_index(v);
                            if (treated_vertices[vv] == false)
                              {
                                treated_vertices[vv] = true;
                                if (vv==1)
                                  {
                                    cell->face(f)->vertex(v) = center_2d+Point<dim_2d>(outer_radius,outer_radius);
                                  }
                              }
                          }
                      }

                  }
            }

        }

        // Make the triangulation_2, a rectangular above the triangulation_1
        {
          const Point<dim_2d> point1 (0, outer_radius),
                point2 (outer_radius, height);

          GridGenerator::hyper_rectangle(triangulation_2, point1, point2);

        }

        // make the triangulation_2d and refine it
        {
          // Merge the two triangulation_1 and triangulation_2
          GridGenerator::merge_triangulations(triangulation_1, triangulation_2, triangulation_2d);

          // Assign boundary indicators to the boundary faces
          /*
           *
           *    /\ y
           *     |
           *      _____3_____
           *     |          |
           *     |          |
           *   4 |          |
           *     |          |
           *     |          | 2
           *     |_         |
           *        \       |
           *      10 \      |
           *         |______|   ____________\  x
           *            1                   /
           */
          {
            typename Triangulation<dim_2d>::active_cell_iterator
            cell = triangulation_2d.begin_active(),
            endc = triangulation_2d.end();
            for (; cell != endc; ++cell)
              {
                for (unsigned int f=0; f<GeometryInfo<dim_2d>::faces_per_cell; ++f)
                  {
                    if (cell->face(f)->at_boundary())
                      {
                        if ( std::fabs(cell->face(f)->center()(1)) < eps )
                          {
                            cell->face(f)->set_manifold_id(1);
                          }
                        else if ( std::fabs(cell->face(f)->center()(0)-outer_radius) < eps )
                          {
                            cell->face(f)->set_manifold_id(2);
                          }
                        else if ( std::fabs(cell->face(f)->center()(1)-height) < eps )
                          {
                            cell->face(f)->set_manifold_id(3);
                          }
                        else if ( std::fabs(cell->face(f)->center()(0)) < eps )
                          {
                            cell->face(f)->set_manifold_id(4);
                          }
                        else
                          {
                            cell->face(f)->set_all_boundary_ids(10);
                          }

                      }
                  }
              }

          }

          const SphericalManifold<dim_2d> inner_boundary_description(center_2d);
          triangulation_2d.set_manifold (10, inner_boundary_description);

          triangulation_2d.refine_global(3);

          triangulation_2d.reset_manifold (10);
        }

        // Extrude the triangulation_2d and make it 3d
//      GridGenerator::extrude_triangulation(triangulation_2d,
//                                           2, thickness, triangulation);
        extrude_triangulation(triangulation_2d,
                              2, thickness, triangulation);

        // Assign boundary indicators to the boundary faces
        /*
         *
         *    /\ y
         *     |
         *      _____3_____
         *     |          |
         *     |          |
         *   4 |          |
         *     |    5|6   |
         *     |          | 2
         *     |_         |
         *        \       |
         *      10 \      |
         *         |______|   ____________\  x
         *            1                   /
         */
        {
          Tensor<1,dim> dist_vector;
          Point<dim> center(center_2d(0), center_2d(1), 0);

          typename Triangulation<dim>::active_cell_iterator
          cell = triangulation.begin_active(),
          endc = triangulation.end();
          for (; cell != endc; ++cell)
            {
              for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
                {
                  if (cell->face(f)->at_boundary())
                    {
                      dist_vector = cell->face(f)->center() - center;

                      if ( std::fabs(dist_vector[1]) < eps )
                        {
                          cell->face(f)->set_manifold_id(1);
                        }
                      else if ( std::fabs(dist_vector[0]-outer_radius) < eps )
                        {
                          cell->face(f)->set_manifold_id(2);
                        }
                      else if ( std::fabs(dist_vector[1]-height) < eps )
                        {
                          cell->face(f)->set_manifold_id(3);
                        }
                      else if ( std::fabs(dist_vector[0]) < eps )
                        {
                          cell->face(f)->set_manifold_id(4);
                        }
                      else if ( std::fabs(dist_vector[2]) < eps )
                        {
                          cell->face(f)->set_manifold_id(5);
                        }
                      else if ( std::fabs(dist_vector[2]-thickness) < eps )
                        {
                          cell->face(f)->set_manifold_id(6);
                        }
                      else
                        {
                          cell->face(f)->set_all_boundary_ids(10);
                        }

                    }
                }
            }

        }

        const CylindricalManifold<dim> inner_boundary_description(2);
        triangulation.set_manifold (10, inner_boundary_description);

        triangulation.refine_global(n_initial_global_refinements);

        triangulation.reset_manifold (10);

      }
    else if (base_mesh == "Cantiliver_beam_3d")
      {
        // A rectangular tube made of Aluminium
        // http://www.google.de/imgres?imgurl=http%3A%2F%2Fwww.americanaluminum.com%2Fimages%2Fstockshape-rectangletube.gif&imgrefurl=http%3A%2F%2Fwww.americanaluminum.com%2Fstandard%2FrectangleTube&h=280&w=300&tbnid=VPDNh4-DJz4wyM%3A&zoom=1&docid=9DoGJCkOeFqiSM&ei=L1AuVfG5GMvtO7DggdAF&tbm=isch&client=ubuntu&iact=rc&uact=3&dur=419&page=1&start=0&ndsp=33&ved=0CGYQrQMwFQ
        // approximation of beam 17250
        // units are in meter

        AssertThrow (dim == 3, ExcNotImplemented());

        const int dim_2d = 2;

        const double length = .7,
                     width = 80e-3,
                     height = 200e-3,
                     thickness_web = 10e-3,
                     thickness_flange = 10e-3;

        Triangulation<dim_2d> triangulation_b,
                      triangulation_t,
                      triangulation_l,
                      triangulation_r,
                      triangulation_2d;

        const double eps = 1e-7 * width;
        // Make the triangulation_b, a rectangular at the bottom of rectangular tube
        {
          const Point<dim_2d> point1 (-width/2, -height/2),
                point2 (width/2, -(height/2)+thickness_flange);

          std::vector<unsigned int> repetitions(dim_2d);
          repetitions[0] = 8;
          repetitions[1] = 1;

          GridGenerator::subdivided_hyper_rectangle(triangulation_b, repetitions, point1, point2);
        }

        // Make the triangulation_t, a rectangular at the top of rectangular tube
        {
          const Point<dim_2d> point1 (-width/2, (height/2)-thickness_flange),
                point2 (width/2, height/2);

          std::vector<unsigned int> repetitions(dim_2d);
          repetitions[0] = 8;
          repetitions[1] = 1;

          GridGenerator::subdivided_hyper_rectangle(triangulation_t, repetitions, point1, point2);
        }

        // Make the triangulation_l, a rectangular at the left of rectangular tube
        {
          const Point<dim_2d> point1 (-width/2, -(height/2)+thickness_flange),
                point2 (-(width/2)+thickness_web, (height/2)-thickness_flange);

          std::vector<unsigned int> repetitions(dim_2d);
          repetitions[0] = 1;
          repetitions[1] = 18;

          GridGenerator::subdivided_hyper_rectangle(triangulation_l, repetitions, point1, point2);
        }

        // Make the triangulation_r, a rectangular at the right of rectangular tube
        {
          const Point<dim_2d> point1 ((width/2)-thickness_web, -(height/2)+thickness_flange),
                point2 (width/2, (height/2)-thickness_flange);

          std::vector<unsigned int> repetitions(dim_2d);
          repetitions[0] = 1;
          repetitions[1] = 18;

          GridGenerator::subdivided_hyper_rectangle(triangulation_r, repetitions, point1, point2);
        }

        // make the triangulation_2d
        {
          // merging every two triangles to make triangulation_2d
          Triangulation<dim_2d> triangulation_bl,
                        triangulation_blr;

          GridGenerator::merge_triangulations(triangulation_b, triangulation_l, triangulation_bl);
          GridGenerator::merge_triangulations(triangulation_bl, triangulation_r, triangulation_blr);
          GridGenerator::merge_triangulations(triangulation_blr, triangulation_t, triangulation_2d);
        }

        // Extrude the triangulation_2d and make it 3d
        const unsigned int n_slices = static_cast<int>(length*1000/20) + 1;
        extrude_triangulation(triangulation_2d,
                              n_slices, length, triangulation);

        // Assign boundary indicators to the boundary faces
        /*
         *
         *                     A
         *            ---------*----------
         *           /                   /|
         *          /                   / |
         *         /                   /  |
         *        /       2    length /   |
         *       /                   /    |
         *      /                   /     |
         *     /                   /      |
         *    /        width      /       |
         *    --------------------        |
         *    | --------1-------. |       |
         *    | :               : |       |
         *    | :               : |h      |
         *    | :      y   z    : |e      |
         *    | :       | /     : |i     /
         *    |1:       |___ x  :1|g    /
         *    | :               : |h   /
         *    | :               : |t  /
         *    | :               : |  /
         *    | :               : | /
         *    | ----------------- |/
         *    ---------1----------/
         *
         *   face id:
         *   Essential boundary condition:
         *   1: z = 0: clamped, fixed in x, y and z directions
         *   Natural/Newmann boundary condition:
         *   2: y = height/2: traction face: pressure on the surface
         *   Quantity of interest:
         *   displacement at Point A (x=0, y=height/2, z=length)
         */
        {
          Tensor<1,dim> dist_vector;
          Point<dim> center(0, 0, 0);

          typename Triangulation<dim>::active_cell_iterator
          cell = triangulation.begin_active(),
          endc = triangulation.end();
          for (; cell != endc; ++cell)
            {
              for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
                {
                  if (cell->face(f)->at_boundary())
                    {
                      dist_vector = cell->face(f)->center() - center;

                      if ( std::fabs(dist_vector[2]) < eps )
                        {
                          cell->face(f)->set_manifold_id(1);
                        }
                      else if ( std::fabs(dist_vector[1]-(height/2)) < eps )
                        {
                          cell->face(f)->set_manifold_id(2);
                        }
                      else
                        {
                          cell->face(f)->set_all_boundary_ids(0);
                        }

                    }
                }
            }

        }

        triangulation.refine_global(n_initial_global_refinements);

      }
    else
      {
        AssertThrow(false, ExcNotImplemented());
      }

    pcout << "    Number of active cells:       "
          << triangulation.n_active_cells()
          << std::endl;
  }



  // @sect4{PlasticityContactProblem::setup_system}

  // The next piece in the puzzle is to set up the DoFHandler, resize
  // vectors and take care of various other status variables such as
  // index sets and constraint matrices.
  //
  // In the following, each group of operations is put into a brace-enclosed
  // block that is being timed by the variable declared at the top of the
  // block (the constructor of the TimerOutput::Scope variable starts the
  // timed section, the destructor that is called at the end of the block
  // stops it again).
  template <int dim>
  void
  ElastoPlasticProblem<dim>::setup_system ()
  {
    /* setup dofs and get index sets for locally owned and relevant dofs */
    TimerOutput::Scope t(computing_timer, "Setup");
    {
      TimerOutput::Scope t(computing_timer, "Setup: distribute DoFs");
      dof_handler.distribute_dofs(fe);
      pcout << "    Number of degrees of freedom: "
            << dof_handler.n_dofs()
            << std::endl;

      locally_owned_dofs = dof_handler.locally_owned_dofs();
      locally_relevant_dofs.clear();
      DoFTools::extract_locally_relevant_dofs(dof_handler,
                                              locally_relevant_dofs);
    }

    /* setup hanging nodes and Dirichlet constraints */
    {
      TimerOutput::Scope t(computing_timer, "Setup: constraints");
      constraints_hanging_nodes.reinit(locally_relevant_dofs);
      DoFTools::make_hanging_node_constraints(dof_handler,
                                              constraints_hanging_nodes);
      constraints_hanging_nodes.close();

      pcout << "   Number of active cells: "
            << triangulation.n_global_active_cells() << std::endl
            << "   Number of degrees of freedom: " << dof_handler.n_dofs()
            << std::endl;

      compute_dirichlet_constraints();
    }

    /* initialization of vectors*/
    {
      TimerOutput::Scope t(computing_timer, "Setup: vectors");
      if (timestep_no==1 || current_refinement_cycle!=0)
        {
          solution.reinit(locally_relevant_dofs, mpi_communicator);
        }
      incremental_displacement.reinit(locally_relevant_dofs, mpi_communicator);
      newton_rhs.reinit(locally_owned_dofs, mpi_communicator);
      newton_rhs_residual.reinit(locally_owned_dofs, mpi_communicator);
      fraction_of_plastic_q_points_per_cell.reinit(triangulation.n_active_cells());
    }

    // Finally, we set up sparsity patterns and matrices.
    // We temporarily (ab)use the system matrix to also build the (diagonal)
    // matrix that we use in eliminating degrees of freedom that are in contact
    // with the obstacle, but we then immediately set the Newton matrix back
    // to zero.
    {
      TimerOutput::Scope t(computing_timer, "Setup: matrix");
      TrilinosWrappers::SparsityPattern sp(locally_owned_dofs,
                                           mpi_communicator);

      DoFTools::make_sparsity_pattern(dof_handler, sp,
                                      constraints_dirichlet_and_hanging_nodes, false,
                                      this_mpi_process);
      sp.compress();
      newton_matrix.reinit(sp);
    }
  }


  // @sect4{PlasticityContactProblem::compute_dirichlet_constraints}

  // This function, broken out of the preceding one, computes the constraints
  // associated with Dirichlet-type boundary conditions and puts them into the
  // <code>constraints_dirichlet_and_hanging_nodes</code> variable by merging
  // with the constraints that come from hanging nodes.
  //
  // As laid out in the introduction, we need to distinguish between two
  // cases:
  // - If the domain is a box, we set the displacement to zero at the bottom,
  //   and allow vertical movement in z-direction along the sides. As
  //   shown in the <code>make_grid()</code> function, the former corresponds
  //   to boundary indicator 6, the latter to 8.
  // - If the domain is a half sphere, then we impose zero displacement along
  //   the curved part of the boundary, associated with boundary indicator zero.
  template <int dim>
  void
  ElastoPlasticProblem<dim>::compute_dirichlet_constraints ()
  {
    constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs);
    constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);

    std::vector<bool> component_mask(dim);

    if (base_mesh == "Timoshenko beam")
      {
        VectorTools::interpolate_boundary_values(dof_handler,
                                                 0,
                                                 EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                 constraints_dirichlet_and_hanging_nodes,
                                                 ComponentMask());
      }
    else if (base_mesh == "Thick_tube_internal_pressure")
      {
        // the boundary x = 0
        component_mask[0] = true;
        component_mask[1] = false;
        VectorTools::interpolate_boundary_values (dof_handler,
                                                  2,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes,
                                                  component_mask);
        // the boundary y = 0
        component_mask[0] = false;
        component_mask[1] = true;
        VectorTools::interpolate_boundary_values (dof_handler,
                                                  3,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes,
                                                  component_mask);
      }
    else if (base_mesh == "Perforated_strip_tension")
      {
        // the boundary x = 0
        component_mask[0] = true;
        component_mask[1] = false;
        component_mask[2] = false;
        VectorTools::interpolate_boundary_values (dof_handler,
                                                  4,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes,
                                                  component_mask);
        // the boundary y = 0
        component_mask[0] = false;
        component_mask[1] = true;
        component_mask[2] = false;
        VectorTools::interpolate_boundary_values (dof_handler,
                                                  1,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes,
                                                  component_mask);
        // the boundary y = imposed incremental displacement
        component_mask[0] = false;
        component_mask[1] = true;
        component_mask[2] = false;
        VectorTools::interpolate_boundary_values (dof_handler,
                                                  3,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes,
                                                  component_mask);
      }
    else if (base_mesh == "Cantiliver_beam_3d")
      {
        // the boundary x = y = z = 0
        component_mask[0] = true;
        component_mask[1] = true;
        component_mask[2] = true;
        VectorTools::interpolate_boundary_values (dof_handler,
                                                  1,
                                                  EquationData::IncrementalBoundaryValues<dim>(present_time, end_time),
                                                  constraints_dirichlet_and_hanging_nodes,
                                                  component_mask);
      }
    else
      {
        AssertThrow(false, ExcNotImplemented());
      }


    constraints_dirichlet_and_hanging_nodes.close();
  }


  // @sect4{PlasticityContactProblem::assemble_newton_system}

  // Given the complexity of the problem, it may come as a bit of a surprise
  // that assembling the linear system we have to solve in each Newton iteration
  // is actually fairly straightforward. The following function builds the Newton
  // right hand side and Newton matrix. It looks fairly innocent because the
  // heavy lifting happens in the call to
  // <code>ConstitutiveLaw::get_linearized_stress_strain_tensors()</code> and in
  // particular in AffineConstraints<double>::distribute_local_to_global(), using the
  // constraints we have previously computed.
  template <int dim>
  void
  ElastoPlasticProblem<dim>::
  assemble_newton_system (const TrilinosWrappers::MPI::Vector &/*linearization_point*/,
                          const TrilinosWrappers::MPI::Vector &delta_linearization_point)
  {
    TimerOutput::Scope t(computing_timer, "Assembling");

    types::boundary_id traction_surface_id = numbers::invalid_boundary_id;
    if (base_mesh == "Timoshenko beam")
      {
        traction_surface_id = 5;
      }
    else if (base_mesh == "Thick_tube_internal_pressure")
      {
        traction_surface_id = 0;
      }
    else if (base_mesh == "Cantiliver_beam_3d")
      {
        traction_surface_id = 2;
      }
    else
      {
        AssertThrow(false, ExcNotImplemented());
      }

    FEValues<dim> fe_values(fe, quadrature_formula,
                            update_values | update_gradients |
                            update_quadrature_points | update_JxW_values);

    FEFaceValues<dim> fe_values_face(fe, face_quadrature_formula,
                                     update_values | 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();
    const unsigned int n_face_q_points = face_quadrature_formula.size();


    const EquationData::BodyForce<dim>     body_force;
    std::vector<Vector<double> >           body_force_values(n_q_points,
                                                             Vector<double>(dim));

    const EquationData::
    IncrementalBoundaryForce<dim>          boundary_force(present_time, end_time);
    std::vector<Vector<double> >           boundary_force_values(n_face_q_points,
        Vector<double>(dim));

    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<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
    std::vector<SymmetricTensor<2, dim> > incremental_strain_tensor(n_q_points);

    typename DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();

    const FEValuesExtractors::Vector displacement(0);

    for (; cell != endc; ++cell)
      if (cell->is_locally_owned())
        {
          fe_values.reinit(cell);
          cell_matrix = 0;
          cell_rhs = 0;

          fe_values[displacement].get_function_symmetric_gradients(delta_linearization_point,
                                                                   incremental_strain_tensor);

          // For assembling the local right hand side contributions, we need
          // to access the prior linearized 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_history
            = reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
          Assert (local_quadrature_points_history >=
                  &quadrature_point_history.front(),
                  ExcInternalError());
          Assert (local_quadrature_points_history <
                  &quadrature_point_history.back(),
                  ExcInternalError());

          // 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);

          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
            {
              SymmetricTensor<2, dim> tmp_strain_tensor_qpoint;
              tmp_strain_tensor_qpoint = local_quadrature_points_history[q_point].old_strain
                                         + incremental_strain_tensor[q_point];

              SymmetricTensor<4, dim> stress_strain_tensor_linearized;
              SymmetricTensor<4, dim> stress_strain_tensor;
              constitutive_law.get_linearized_stress_strain_tensors(tmp_strain_tensor_qpoint,
                                                                    stress_strain_tensor_linearized,
                                                                    stress_strain_tensor);

              Tensor<1, dim> rhs_values_body_force;
              for (unsigned int i = 0; i < dim; ++i)
                {
                  rhs_values_body_force[i] = body_force_values[q_point][i];
                }

              for (unsigned int i = 0; i < dofs_per_cell; ++i)
                {
                  // Having computed the stress-strain tensor and its linearization,
                  // we can now put together the parts of the matrix and right hand side.
                  // In both, we need the linearized stress-strain tensor times the
                  // symmetric gradient of $\varphi_i$, i.e. the term $I_\Pi\varepsilon(\varphi_i)$,
                  // so we introduce an abbreviation of this term. Recall that the
                  // matrix corresponds to the bilinear form
                  // $A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))$ in the
                  // notation of the accompanying publication, whereas the right
                  // hand side is $F_i=([I_\Pi-P_\Pi C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))$
                  // where $u$ is the current linearization points (typically the last solution).
                  // This might suggest that the right hand side will be zero if the material
                  // is completely elastic (where $I_\Pi=P_\Pi$) but this ignores the fact
                  // that the right hand side will also contain contributions from
                  // non-homogeneous constraints due to the contact.
                  //
                  // The code block that follows this adds contributions that are due to
                  // boundary forces, should there be any.
                  const SymmetricTensor<2, dim>
                  stress_phi_i = stress_strain_tensor_linearized
                                 * fe_values[displacement].symmetric_gradient(i, q_point);

                  for (unsigned int j = 0; j < dofs_per_cell; ++j)
                    cell_matrix(i, j) += (stress_phi_i
                                          * fe_values[displacement].symmetric_gradient(j, q_point)
                                          * fe_values.JxW(q_point));

                  cell_rhs(i) += (
                                   ( stress_phi_i
                                     * incremental_strain_tensor[q_point] )
                                   -
                                   ( ( stress_strain_tensor
                                       * fe_values[displacement].symmetric_gradient(i, q_point))
                                     * tmp_strain_tensor_qpoint )
                                   +
                                   ( fe_values[displacement].value(i, q_point)
                                     * rhs_values_body_force )
                                 ) * fe_values.JxW(q_point);

                }
            }

          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
            if (cell->face(face)->at_boundary()
                &&
                cell->face(face)->boundary_id() == traction_surface_id)
              {
                fe_values_face.reinit(cell, face);

                boundary_force.vector_value_list(fe_values_face.get_quadrature_points(),
                                                 boundary_force_values);

                for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                  {
                    Tensor<1, dim> rhs_values;
                    for (unsigned int i = 0; i < dim; ++i)
                      {
                        rhs_values[i] = boundary_force_values[q_point][i];
                      }
                    for (unsigned int i = 0; i < dofs_per_cell; ++i)
                      cell_rhs(i) += (fe_values_face[displacement].value(i, q_point)
                                      * rhs_values
                                      * fe_values_face.JxW(q_point));
                  }
              }

          cell->get_dof_indices(local_dof_indices);
          constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_matrix, cell_rhs,
              local_dof_indices,
              newton_matrix,
              newton_rhs,
              true);

        }

    newton_matrix.compress(VectorOperation::add);
    newton_rhs.compress(VectorOperation::add);
  }



  // @sect4{PlasticityContactProblem::compute_nonlinear_residual}

  // The following function computes the nonlinear residual of the equation
  // given the current solution (or any other linearization point). This
  // is needed in the linear search algorithm where we need to try various
  // linear combinations of previous and current (trial) solution to
  // compute the (real, globalized) solution of the current Newton step.
  //
  // That said, in a slight abuse of the name of the function, it actually
  // does significantly more. For example, it also computes the vector
  // that corresponds to the Newton residual but without eliminating
  // constrained degrees of freedom. We need this vector to compute contact
  // forces and, ultimately, to compute the next active set. Likewise, by
  // keeping track of how many quadrature points we encounter on each cell
  // that show plastic yielding, we also compute the
  // <code>fraction_of_plastic_q_points_per_cell</code> vector that we
  // can later output to visualize the plastic zone. In both of these cases,
  // the results are not necessary as part of the line search, and so we may
  // be wasting a small amount of time computing them. At the same time, this
  // information appears as a natural by-product of what we need to do here
  // anyway, and we want to collect it once at the end of each Newton
  // step, so we may as well do it here.
  //
  // The actual implementation of this function should be rather obvious:
  template <int dim>
  void
  ElastoPlasticProblem<dim>::
  compute_nonlinear_residual (const TrilinosWrappers::MPI::Vector &linearization_point)
  {
    types::boundary_id traction_surface_id = numbers::invalid_boundary_id;
    if (base_mesh == "Timoshenko beam")
      {
        traction_surface_id = 5;
      }
    else if (base_mesh == "Thick_tube_internal_pressure")
      {
        traction_surface_id = 0;
      }
    else if (base_mesh == "Cantiliver_beam_3d")
      {
        traction_surface_id = 2;
      }
    else
      {
        AssertThrow(false, ExcNotImplemented());
      }

    FEValues<dim> fe_values(fe, quadrature_formula,
                            update_values | update_gradients | update_quadrature_points |
                            update_JxW_values);

    FEFaceValues<dim> fe_values_face(fe, face_quadrature_formula,
                                     update_values | 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();
    const unsigned int n_face_q_points = face_quadrature_formula.size();

    const EquationData::BodyForce<dim>     body_force;
    std::vector<Vector<double> >           body_force_values(n_q_points,
                                                             Vector<double>(dim));

    const EquationData::
    IncrementalBoundaryForce<dim>          boundary_force(present_time, end_time);
    std::vector<Vector<double> >           boundary_force_values(n_face_q_points,
        Vector<double>(dim));

    Vector<double> cell_rhs(dofs_per_cell);

    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

    const FEValuesExtractors::Vector displacement(0);

    newton_rhs_residual = 0;

    fraction_of_plastic_q_points_per_cell = 0;

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

          std::vector<SymmetricTensor<2, dim> > strain_tensors(n_q_points);
          fe_values[displacement].get_function_symmetric_gradients(linearization_point,
                                                                   strain_tensors);

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

          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
            {
              SymmetricTensor<4, dim> stress_strain_tensor;
              const bool q_point_is_plastic
                = constitutive_law.get_stress_strain_tensor(strain_tensors[q_point],
                                                            stress_strain_tensor);
              if (q_point_is_plastic)
                ++fraction_of_plastic_q_points_per_cell(cell_number);

              Tensor<1, dim> rhs_values_body_force;
              for (unsigned int i = 0; i < dim; ++i)
                {
                  rhs_values_body_force[i] = body_force_values[q_point][i];
                }

              for (unsigned int i = 0; i < dofs_per_cell; ++i)
                {
                  cell_rhs(i) += (fe_values[displacement].value(i, q_point)
                                  * rhs_values_body_force
                                  -
                                  strain_tensors[q_point]
                                  * stress_strain_tensor
                                  * fe_values[displacement].symmetric_gradient(i, q_point)
                                 )
                                 * fe_values.JxW(q_point);

                  Tensor<1, dim> rhs_values;
                  rhs_values = 0;
                  cell_rhs(i) += (fe_values[displacement].value(i, q_point)
                                  * rhs_values
                                  * fe_values.JxW(q_point));
                }
            }

          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
            if (cell->face(face)->at_boundary()
                && cell->face(face)->boundary_id() == traction_surface_id)
              {
                fe_values_face.reinit(cell, face);

                boundary_force.vector_value_list(fe_values_face.get_quadrature_points(),
                                                 boundary_force_values);

                for (unsigned int q_point = 0; q_point < n_face_q_points;
                     ++q_point)
                  {
                    Tensor<1, dim> rhs_values;
                    for (unsigned int i = 0; i < dim; ++i)
                      {
                        rhs_values[i] = boundary_force_values[q_point][i];
                      }
                    for (unsigned int i = 0; i < dofs_per_cell; ++i)
                      cell_rhs(i) += (fe_values_face[displacement].value(i, q_point) * rhs_values
                                      * fe_values_face.JxW(q_point));
                  }
              }

          cell->get_dof_indices(local_dof_indices);
          constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(cell_rhs,
              local_dof_indices,
              newton_rhs_residual);

        }

    fraction_of_plastic_q_points_per_cell /= quadrature_formula.size();
    newton_rhs_residual.compress(VectorOperation::add);

  }





  // @sect4{PlasticityContactProblem::solve_newton_system}

  // The last piece before we can discuss the actual Newton iteration
  // on a single mesh is the solver for the linear systems. There are
  // a couple of complications that slightly obscure the code, but
  // mostly it is just setup then solve. Among the complications are:
  //
  // - For the hanging nodes we have to apply
  //   the AffineConstraints<double>::set_zero function to newton_rhs.
  //   This is necessary if a hanging node with solution value $x_0$
  //   has one neighbor with value $x_1$ which is in contact with the
  //   obstacle and one neighbor $x_2$ which is not in contact. Because
  //   the update for the former will be prescribed, the hanging node constraint
  //   will have an inhomogeneity and will look like $x_0 = x_1/2 + \text{gap}/2$.
  //   So the corresponding entries in the
  //   ride-hang-side are non-zero with a
  //   meaningless value. These values we have to
  //   to set to zero.
  // - Like in step-40, we need to shuffle between vectors that do and do
  //   do not have ghost elements when solving or using the solution.
  //
  // The rest of the function is similar to step-40 and
  // step-41 except that we use a BiCGStab solver
  // instead of CG. This is due to the fact that for very small hardening
  // parameters $\gamma$, the linear system becomes almost semidefinite though
  // still symmetric. BiCGStab appears to have an easier time with such linear
  // systems.
  template <int dim>
  void
  ElastoPlasticProblem<dim>::solve_newton_system ()
  {
    TimerOutput::Scope t(computing_timer, "Solve");

    TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs, mpi_communicator);
    distributed_solution = incremental_displacement;

    constraints_hanging_nodes.set_zero(distributed_solution);
    constraints_hanging_nodes.set_zero(newton_rhs);

    // ------- Solver Bicgstab --- Preconditioner AMG -------------------
//    TrilinosWrappers::PreconditionAMG preconditioner;
//    {
//      TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner");
//
//      std::vector<std::vector<bool> > constant_modes;
//      DoFTools::extract_constant_modes(dof_handler, ComponentMask(),
//                                       constant_modes);
//
//      TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
//      additional_data.constant_modes = constant_modes;
//      additional_data.elliptic = true;
//      additional_data.n_cycles = 1;
//      additional_data.w_cycle = false;
//      additional_data.output_details = false;
//      additional_data.smoother_sweeps = 2;
//      additional_data.aggregation_threshold = 1e-2;
//
//      preconditioner.initialize(newton_matrix, additional_data);
//    }

//    {
//      TimerOutput::Scope t(computing_timer, "Solve: iterate");
//
//      TrilinosWrappers::MPI::Vector tmp(locally_owned_dofs, mpi_communicator);
//
////      const double relative_accuracy = 1e-8;
//      const double relative_accuracy = 1e-2;
//      const double solver_tolerance  = relative_accuracy
//                                       * newton_matrix.residual(tmp, distributed_solution,
//                                                                newton_rhs);
//
//      SolverControl solver_control(newton_matrix.m(),
//                                   solver_tolerance);
//      SolverBicgstab<TrilinosWrappers::MPI::Vector> solver(solver_control);
//      solver.solve(newton_matrix, distributed_solution,
//                   newton_rhs, preconditioner);
//
//      pcout << "         Error: " << solver_control.initial_value()
//            << " -> " << solver_control.last_value() << " in "
//            << solver_control.last_step() << " Bicgstab iterations."
//            << std::endl;
//    }

    // ------- Solver CG --- Preconditioner SSOR -------------------
    TrilinosWrappers::PreconditionSSOR preconditioner;
    {
      TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner");

      TrilinosWrappers::PreconditionSSOR::AdditionalData additional_data;
      preconditioner.initialize(newton_matrix, additional_data);
    }

    {
      TimerOutput::Scope t(computing_timer, "Solve: iterate");

      TrilinosWrappers::MPI::Vector tmp(locally_owned_dofs, mpi_communicator);

//      const double relative_accuracy = 1e-8;
      const double relative_accuracy = 1e-2;
      const double solver_tolerance  = relative_accuracy
                                       * newton_matrix.residual(tmp, distributed_solution,
                                                                newton_rhs);

//      SolverControl solver_control(newton_matrix.m(),
//                                   solver_tolerance);
      SolverControl solver_control(10*newton_matrix.m(),
                                   solver_tolerance);
      SolverCG<TrilinosWrappers::MPI::Vector> solver(solver_control);
      solver.solve(newton_matrix, distributed_solution,
                   newton_rhs, preconditioner);

      pcout << "         Error: " << solver_control.initial_value()
            << " -> " << solver_control.last_value() << " in "
            << solver_control.last_step() << " CG iterations."
            << std::endl;
    }
    // ........................................................

    constraints_dirichlet_and_hanging_nodes.distribute(distributed_solution);

    incremental_displacement = distributed_solution;
  }


  // @sect4{PlasticityContactProblem::solve_newton}

  // This is, finally, the function that implements the damped Newton method
  // on the current mesh. There are two nested loops: the outer loop for the Newton
  // iteration and the inner loop for the line search which
  // will be used only if necessary. To obtain a good and reasonable
  // starting value we solve an elastic problem in very first Newton step on each
  // mesh (or only on the first mesh if we transfer solutions between meshes). We
  // do so by setting the yield stress to an unreasonably large value in these
  // iterations and then setting it back to the correct value in subsequent
  // iterations.
  //
  // Other than this, the top part of this function should be reasonably
  // obvious:
  template <int dim>
  void
  ElastoPlasticProblem<dim>::solve_newton ()
  {
    TrilinosWrappers::MPI::Vector old_solution(locally_owned_dofs, mpi_communicator);
    TrilinosWrappers::MPI::Vector residual(locally_owned_dofs, mpi_communicator);
    TrilinosWrappers::MPI::Vector tmp_vector(locally_owned_dofs, mpi_communicator);
    TrilinosWrappers::MPI::Vector locally_relevant_tmp_vector(locally_relevant_dofs, mpi_communicator);
    TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs, mpi_communicator);
    TrilinosWrappers::MPI::Vector tmp_solution(locally_owned_dofs, mpi_communicator);

    double residual_norm;
    double previous_residual_norm = -std::numeric_limits<double>::max();

    double disp_norm,
           previous_disp_norm = 0;

    const double correct_sigma = sigma_0;

    const unsigned int max_newton_iter = 100;

    for (unsigned int newton_step = 1; newton_step <= max_newton_iter; ++newton_step)
      {
        if (newton_step == 1
            &&
            ((transfer_solution && timestep_no == 1)
             ||
             !transfer_solution))
          constitutive_law.set_sigma_0(1e+10);
        else
          constitutive_law.set_sigma_0(correct_sigma);

        pcout << " " << std::endl;
        pcout << "   Newton iteration " << newton_step << std::endl;

        pcout << "      Assembling system... " << std::endl;
        newton_matrix       = 0;
        newton_rhs          = 0;
        newton_rhs_residual = 0;

        tmp_solution = solution;
        tmp_solution += incremental_displacement;
        assemble_newton_system(tmp_solution,
                               incremental_displacement);

        pcout << "      Solving system... " << std::endl;
        solve_newton_system();

        // It gets a bit more hairy after we have computed the
        // trial solution $\tilde{\mathbf u}$ of the current Newton step.
        // We handle a highly nonlinear problem so we have to damp
        // Newton's method using a line search. To understand how we do this,
        // recall that in our formulation, we compute a trial solution
        // in each Newton step and not the update between old and new solution.
        // Since the solution set is a convex set, we will use a line
        // search that tries linear combinations of the
        // previous and the trial solution to guarantee that the
        // damped solution is in our solution set again.
        // At most we apply 5 damping steps.
        //
        // There are exceptions to when we use a line search. First,
        // if this is the first Newton step on any mesh, then we don't have
        // any point to compare the residual to, so we always accept a full
        // step. Likewise, if this is the second Newton step on the first mesh (or
        // the second on any mesh if we don't transfer solutions from
        // mesh to mesh), then we have computed the first of these steps using
        // just an elastic model (see how we set the yield stress sigma to
        // an unreasonably large value above). In this case, the first Newton
        // solution was a purely elastic one, the second one a plastic one,
        // and any linear combination would not necessarily be expected to
        // lie in the feasible set -- so we just accept the solution we just
        // got.
        //
        // In either of these two cases, we bypass the line search and just
        // update residual and other vectors as necessary.
        if ((newton_step==1)
            ||
            (transfer_solution && newton_step == 2 && current_refinement_cycle == 0)
            ||
            (!transfer_solution && newton_step == 2))
          {
            tmp_solution = solution;
            tmp_solution += incremental_displacement;
            compute_nonlinear_residual(tmp_solution);
            old_solution = incremental_displacement;

            residual = newton_rhs_residual;

            residual.compress(VectorOperation::insert);

            residual_norm = residual.l2_norm();

            pcout << "      Accepting Newton solution with residual: "
                  << residual_norm << std::endl;
          }
        else
          {
            for (unsigned int i = 0; i < 5; ++i)
              {
                distributed_solution = incremental_displacement;

                const double alpha = std::pow(0.5, static_cast<double>(i));
                tmp_vector = old_solution;
                tmp_vector.sadd(1 - alpha, alpha, distributed_solution);

                TimerOutput::Scope t(computing_timer, "Residual and lambda");

                locally_relevant_tmp_vector = tmp_vector;
                tmp_solution = solution;
                tmp_solution += locally_relevant_tmp_vector;
                compute_nonlinear_residual(tmp_solution);
                residual = newton_rhs_residual;

                residual.compress(VectorOperation::insert);

                residual_norm = residual.l2_norm();

                pcout << "      Residual of the system: "
                      << residual_norm << std::endl
                      << "         with a damping parameter alpha = " << alpha
                      << std::endl;

                if (residual_norm < previous_residual_norm)
                  break;
              }

            incremental_displacement = tmp_vector;
            old_solution = incremental_displacement;
          }

        disp_norm = incremental_displacement.l2_norm();


        // The final step is to check for convergence. If the residual is
        // less than a threshold of $10^{-10}$, then we terminate
        // the iteration on the current mesh:
//        if (residual_norm < 1e-10)
        if (residual_norm < 1e-7)
          break;

        pcout << "    difference of two consecutive incremental displacement l2 norm : "
              << std::abs(disp_norm - previous_disp_norm) << std::endl;
        if ( std::abs(disp_norm - previous_disp_norm) < 1e-10 &&
             (residual_norm < 1e-5 || std::abs(residual_norm - previous_residual_norm)<1e-9) )
          {
            pcout << " Convergence by difference of two consecutive solution! " << std::endl;
            break;
          }


        previous_residual_norm = residual_norm;
        previous_disp_norm = disp_norm;
      }
  }

  // @sect4{PlasticityContactProblem::compute_error}

  template <int dim>
  void
  ElastoPlasticProblem<dim>::compute_error ()
  {
    TrilinosWrappers::MPI::Vector   tmp_solution(locally_owned_dofs, mpi_communicator);
    tmp_solution = solution;
    tmp_solution += incremental_displacement;

    estimated_error_per_cell.reinit (triangulation.n_active_cells());
    if (error_estimation_strategy == ErrorEstimationStrategy::kelly_error)
      {
        using FunctionMap = std::map<types::boundary_id, const Function<dim> *>;

        KellyErrorEstimator<dim>::estimate(dof_handler,
                                           QGauss<dim - 1>(fe.degree + 2),
                                           std::map<types::boundary_id, const Function<dim> *>(),
                                           tmp_solution,
                                           estimated_error_per_cell);

      }
    else if (error_estimation_strategy == ErrorEstimationStrategy::residual_error)
      {
        compute_error_residual(tmp_solution);

      }
    else if (error_estimation_strategy == ErrorEstimationStrategy::weighted_residual_error)
      {
        // make a non-parallel copy of tmp_solution
        Vector<double> copy_solution(tmp_solution);

        // the dual function definition (it should be defined previously, e.g. input file)
        if (base_mesh == "Timoshenko beam")
          {
            double length = .48,
                   depth  = .12;

            const Point<dim> evaluation_point(length, -depth/2);

            DualFunctional::PointValuesEvaluation<dim> dual_functional(evaluation_point);

            DualSolver<dim> dual_solver(triangulation, fe,
                                        copy_solution,
                                        constitutive_law, dual_functional,
                                        timestep_no, output_dir, base_mesh,
                                        present_time, end_time);

            dual_solver.compute_error_DWR (estimated_error_per_cell);

          }
        else if (base_mesh == "Thick_tube_internal_pressure")
          {
            const unsigned int face_id = 0;
            std::vector<std::vector<unsigned int> > comp_stress(dim);
            for (unsigned int i=0; i!=dim; ++i)
              {
                comp_stress[i].resize(dim);
                for (unsigned int j=0; j!=dim; ++j)
                  {
                    comp_stress[i][j] = 1;
                  }
              }

            DualFunctional::MeanStressFace<dim> dual_functional(face_id, comp_stress);

            DualSolver<dim> dual_solver(triangulation, fe,
                                        copy_solution,
                                        constitutive_law, dual_functional,
                                        timestep_no, output_dir, base_mesh,
                                        present_time, end_time);

            dual_solver.compute_error_DWR (estimated_error_per_cell);

          }
        else if (base_mesh == "Perforated_strip_tension")
          {
            // .........................................
            // Mean stress_yy over the bottom boundary
            const unsigned int face_id = 1;
            std::vector<std::vector<unsigned int> > comp_stress(dim);
            for (unsigned int i=0; i!=dim; ++i)
              {
                comp_stress[i].resize(dim);
                for (unsigned int j=0; j!=dim; ++j)
                  {
                    comp_stress[i][j] = 0;
                  }
              }
            comp_stress[1][1] = 1;

            DualFunctional::MeanStressFace<dim> dual_functional(face_id, comp_stress);

            // .........................................

            DualSolver<dim> dual_solver(triangulation, fe,
                                        copy_solution,
                                        constitutive_law, dual_functional,
                                        timestep_no, output_dir, base_mesh,
                                        present_time, end_time);

            dual_solver.compute_error_DWR (estimated_error_per_cell);

          }
        else if (base_mesh == "Cantiliver_beam_3d")
          {
            // Quantity of interest:
            // -----------------------------------------------------------
            // displacement at Point A (x=0, y=height/2, z=length)
            /*
            const double length = .7,
                         height = 200e-3;

            const Point<dim> evaluation_point(0, height/2, length);

            DualFunctional::PointValuesEvaluation<dim> dual_functional(evaluation_point);
            */

            // -----------------------------------------------------------
            // Mean stress at the specified domain is of interest.
            // The interest domains are located on the bottom and top of the flanges
            // close to the clamped face, z = 0
            // top domain: height/2 - thickness_flange <= y <= height/2
            //             0 <= z <= 2 * thickness_flange
            // bottom domain: -height/2 <= y <= -height/2 + thickness_flange
            //             0 <= z <= 2 * thickness_flange

            std::vector<std::vector<unsigned int> > comp_stress(dim);
            for (unsigned int i=0; i!=dim; ++i)
              {
                comp_stress[i].resize(dim);
                for (unsigned int j=0; j!=dim; ++j)
                  {
                    comp_stress[i][j] = 1;
                  }
              }
            DualFunctional::MeanStressDomain<dim> dual_functional(base_mesh, comp_stress);

            // -----------------------------------------------------------

            DualSolver<dim> dual_solver(triangulation, fe,
                                        copy_solution,
                                        constitutive_law, dual_functional,
                                        timestep_no, output_dir, base_mesh,
                                        present_time, end_time);

            dual_solver.compute_error_DWR (estimated_error_per_cell);

          }
        else
          {
            AssertThrow(false, ExcNotImplemented());
          }


      }
    else
      {
        AssertThrow(false, ExcNotImplemented());
      }


    relative_error = estimated_error_per_cell.l2_norm() / tmp_solution.l2_norm();

    pcout << "Estimated relative error = " << relative_error << std::endl;

  }

  template <int dim>
  void
  ElastoPlasticProblem<dim>::compute_error_residual (const TrilinosWrappers::MPI::Vector &tmp_solution)
  {
    FEValues<dim> fe_values(fe, quadrature_formula,
                            update_values    |
                            update_gradients |
                            update_hessians  |
                            update_quadrature_points |
                            update_JxW_values);

    const unsigned int n_q_points      = quadrature_formula.size();
    std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
    SymmetricTensor<4, dim> stress_strain_tensor_linearized;
    SymmetricTensor<4, dim> stress_strain_tensor;
    Tensor<5, dim>          stress_strain_tensor_grad;
    std::vector<std::vector<Tensor<2,dim> > > cell_hessians (n_q_points);
    for (unsigned int i=0; i!=n_q_points; ++i)
      {
        cell_hessians[i].resize (dim);
      }
    const EquationData::BodyForce<dim> body_force;

    std::vector<Vector<double> > body_force_values (n_q_points, Vector<double>(dim));
    const FEValuesExtractors::Vector displacement(0);


    FEFaceValues<dim> fe_face_values_cell(fe, face_quadrature_formula,
                                          update_values           |
                                          update_quadrature_points|
                                          update_gradients        |
                                          update_JxW_values       |
                                          update_normal_vectors),
                                          fe_face_values_neighbor (fe, face_quadrature_formula,
                                              update_values     |
                                              update_gradients  |
                                              update_JxW_values |
                                              update_normal_vectors);
    FESubfaceValues<dim> fe_subface_values_cell (fe, face_quadrature_formula,
                                                 update_gradients);

    const unsigned int n_face_q_points = face_quadrature_formula.size();
    std::vector<Vector<double> > jump_residual (n_face_q_points, Vector<double>(dim));
    std::vector<std::vector<Tensor<1,dim> > > cell_grads(n_face_q_points);
    for (unsigned int i=0; i!=n_face_q_points; ++i)
      {
        cell_grads[i].resize (dim);
      }
    std::vector<std::vector<Tensor<1,dim> > > neighbor_grads(n_face_q_points);
    for (unsigned int i=0; i!=n_face_q_points; ++i)
      {
        neighbor_grads[i].resize (dim);
      }
    SymmetricTensor<2, dim> q_cell_strain_tensor;
    SymmetricTensor<2, dim> q_neighbor_strain_tensor;
    SymmetricTensor<4, dim> cell_stress_strain_tensor;
    SymmetricTensor<4, dim> neighbor_stress_strain_tensor;


    typename std::map<typename DoFHandler<dim>::face_iterator, Vector<double> >
    face_integrals;
    typename DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();
    for (; cell!=endc; ++cell)
      if (cell->is_locally_owned())
        {
          for (unsigned int face_no=0;
               face_no<GeometryInfo<dim>::faces_per_cell;
               ++face_no)
            {
              face_integrals[cell->face(face_no)].reinit (dim);
              face_integrals[cell->face(face_no)] = -1e20;
            }
        }

    std::vector<Vector<float> > error_indicators_vector;
    error_indicators_vector.resize( triangulation.n_active_cells(),
                                    Vector<float>(dim) );

    // ----------------- estimate_some -------------------------
    cell = dof_handler.begin_active();
    unsigned int present_cell = 0;
    for (; cell!=endc; ++cell, ++present_cell)
      if (cell->is_locally_owned())
        {
          // --------------- integrate_over_cell -------------------
          fe_values.reinit(cell);
          body_force.vector_value_list(fe_values.get_quadrature_points(),
                                       body_force_values);
          fe_values[displacement].get_function_symmetric_gradients(tmp_solution,
                                                                   strain_tensor);
          fe_values.get_function_hessians(tmp_solution, cell_hessians);

          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
            {
              constitutive_law.get_linearized_stress_strain_tensors(strain_tensor[q_point],
                                                                    stress_strain_tensor_linearized,
                                                                    stress_strain_tensor);
              constitutive_law.get_grad_stress_strain_tensor(strain_tensor[q_point],
                                                             cell_hessians[q_point],
                                                             stress_strain_tensor_grad);

              for (unsigned int i=0; i!=dim; ++i)
                {
                  error_indicators_vector[present_cell](i) +=
                    body_force_values[q_point](i)*fe_values.JxW(q_point);
                  for (unsigned int j=0; j!=dim; ++j)
                    {
                      for (unsigned int k=0; k!=dim; ++k)
                        {
                          for (unsigned int l=0; l!=dim; ++l)
                            {
                              error_indicators_vector[present_cell](i) +=
                                ( stress_strain_tensor[i][j][k][l]*
                                  0.5*(cell_hessians[q_point][k][l][j]
                                       +
                                       cell_hessians[q_point][l][k][j])
                                  + stress_strain_tensor_grad[i][j][k][l][j] * strain_tensor[q_point][k][l]
                                ) *
                                fe_values.JxW(q_point);
                            }
                        }
                    }

                }

            }
          // -------------------------------------------------------
          // compute face_integrals
          for (unsigned int face_no=0;
               face_no<GeometryInfo<dim>::faces_per_cell;
               ++face_no)
            {
              if (cell->face(face_no)->at_boundary())
                {
                  for (unsigned int id=0; id!=dim; ++id)
                    {
                      face_integrals[cell->face(face_no)](id) = 0;
                    }
                  continue;
                }

              if ((cell->neighbor(face_no)->has_children() == false) &&
                  (cell->neighbor(face_no)->level() == cell->level()) &&
                  (cell->neighbor(face_no)->index() < cell->index()))
                continue;

              if (cell->at_boundary(face_no) == false)
                if (cell->neighbor(face_no)->level() < cell->level())
                  continue;


              if (cell->face(face_no)->has_children() == false)
                {
                  // ------------- integrate_over_regular_face -----------
                  fe_face_values_cell.reinit(cell, face_no);
                  fe_face_values_cell.get_function_gradients (tmp_solution,
                                                              cell_grads);

                  Assert (cell->neighbor(face_no).state() == IteratorState::valid,
                          ExcInternalError());
                  const unsigned int
                  neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
                  const typename DoFHandler<dim>::active_cell_iterator
                  neighbor = cell->neighbor(face_no);

                  fe_face_values_neighbor.reinit(neighbor, neighbor_neighbor);
                  fe_face_values_neighbor.get_function_gradients (tmp_solution,
                                                                  neighbor_grads);

                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      q_cell_strain_tensor = 0.;
                      q_neighbor_strain_tensor = 0.;
                      for (unsigned int i=0; i!=dim; ++i)
                        {
                          for (unsigned int j=0; j!=dim; ++j)
                            {
                              q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
                                                                cell_grads[q_point][j][i] );
                              q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
                                                                    neighbor_grads[q_point][j][i] );
                            }
                        }

                      constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
                                                                 cell_stress_strain_tensor);
                      constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
                                                                 neighbor_stress_strain_tensor);

                      jump_residual[q_point] = 0.;
                      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)
                                    {
                                      jump_residual[q_point](i) += (cell_stress_strain_tensor[i][j][k][l]*
                                                                    q_cell_strain_tensor[k][l]
                                                                    -
                                                                    neighbor_stress_strain_tensor[i][j][k][l]*
                                                                    q_neighbor_strain_tensor[k][l] )*
                                                                   fe_face_values_cell.normal_vector(q_point)[j];
                                    }
                                }
                            }
                        }

                    }

                  Vector<double> face_integral_vector(dim);
                  face_integral_vector = 0;
                  for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                    {
                      for (unsigned int i=0; i!=dim; ++i)
                        {
                          face_integral_vector(i) += jump_residual[q_point](i) *
                                                     fe_face_values_cell.JxW(q_point);
                        }
                    }

                  Assert (face_integrals.find (cell->face(face_no)) != face_integrals.end(),
                          ExcInternalError());

                  for (unsigned int i=0; i!=dim; ++i)
                    {
                      Assert (face_integrals[cell->face(face_no)](i) == -1e20,
                              ExcInternalError());
                      face_integrals[cell->face(face_no)](i) = face_integral_vector(i);

                    }

                  // -----------------------------------------------------
                }
              else
                {
                  // ------------- integrate_over_irregular_face ---------
                  const typename DoFHandler<dim>::face_iterator
                  face = cell->face(face_no);
                  const typename DoFHandler<dim>::cell_iterator
                  neighbor = cell->neighbor(face_no);
                  Assert (neighbor.state() == IteratorState::valid,
                          ExcInternalError());
                  Assert (neighbor->has_children(),
                          ExcInternalError());

                  const unsigned int
                  neighbor_neighbor = cell->neighbor_of_neighbor (face_no);

                  for (unsigned int subface_no=0;
                       subface_no<face->n_children(); ++subface_no)
                    {
                      const typename DoFHandler<dim>::active_cell_iterator
                      neighbor_child = cell->neighbor_child_on_subface (face_no, subface_no);
                      Assert (neighbor_child->face(neighbor_neighbor) ==
                              cell->face(face_no)->child(subface_no),
                              ExcInternalError());

                      fe_subface_values_cell.reinit (cell, face_no, subface_no);
                      fe_subface_values_cell.get_function_gradients (tmp_solution,
                                                                     cell_grads);
                      fe_face_values_neighbor.reinit (neighbor_child,
                                                      neighbor_neighbor);
                      fe_face_values_neighbor.get_function_gradients (tmp_solution,
                                                                      neighbor_grads);

                      for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                        {
                          q_cell_strain_tensor = 0.;
                          q_neighbor_strain_tensor = 0.;
                          for (unsigned int i=0; i!=dim; ++i)
                            {
                              for (unsigned int j=0; j!=dim; ++j)
                                {
                                  q_cell_strain_tensor[i][j] = 0.5*(cell_grads[q_point][i][j] +
                                                                    cell_grads[q_point][j][i] );
                                  q_neighbor_strain_tensor[i][j] = 0.5*(neighbor_grads[q_point][i][j] +
                                                                        neighbor_grads[q_point][j][i] );
                                }
                            }

                          constitutive_law.get_stress_strain_tensor (q_cell_strain_tensor,
                                                                     cell_stress_strain_tensor);
                          constitutive_law.get_stress_strain_tensor (q_neighbor_strain_tensor,
                                                                     neighbor_stress_strain_tensor);

                          jump_residual[q_point] = 0.;
                          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)
                                        {
                                          jump_residual[q_point](i) += (-cell_stress_strain_tensor[i][j][k][l]*
                                                                        q_cell_strain_tensor[k][l]
                                                                        +
                                                                        neighbor_stress_strain_tensor[i][j][k][l]*
                                                                        q_neighbor_strain_tensor[k][l] )*
                                                                       fe_face_values_neighbor.normal_vector(q_point)[j];
                                        }
                                    }
                                }
                            }

                        }

                      Vector<double> face_integral_vector(dim);
                      face_integral_vector = 0;
                      for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
                        {
                          for (unsigned int i=0; i!=dim; ++i)
                            {
                              face_integral_vector(i) += jump_residual[q_point](i) *
                                                         fe_face_values_neighbor.JxW(q_point);
                            }
                        }

                      for (unsigned int i=0; i!=dim; ++i)
                        {
                          face_integrals[neighbor_child->face(neighbor_neighbor)](i) = face_integral_vector(i);
                        }

                    }

                  Vector<double> sum (dim);
                  sum = 0;
                  for (unsigned int subface_no=0;
                       subface_no<face->n_children(); ++subface_no)
                    {
                      Assert (face_integrals.find(face->child(subface_no)) !=
                              face_integrals.end(),
                              ExcInternalError());
                      for (unsigned int i=0; i!=dim; ++i)
                        {
                          Assert (face_integrals[face->child(subface_no)](i) != -1e20,
                                  ExcInternalError());
                          sum(i) += face_integrals[face->child(subface_no)](i);
                        }
                    }
                  for (unsigned int i=0; i!=dim; ++i)
                    {
                      face_integrals[face](i) = sum(i);
                    }


                  // -----------------------------------------------------
                }


            }
        }
    // ----------------------------------------------------------

    present_cell=0;
    cell = dof_handler.begin_active();
    for (; cell!=endc; ++cell, ++present_cell)
      if (cell->is_locally_owned())
        {
          for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
               ++face_no)
            {
              Assert(face_integrals.find(cell->face(face_no)) !=
                     face_integrals.end(),
                     ExcInternalError());

              for (unsigned int id=0; id!=dim; ++id)
                {
                  error_indicators_vector[present_cell](id)
                  -= 0.5*face_integrals[cell->face(face_no)](id);
                }

            }

          estimated_error_per_cell(present_cell) = error_indicators_vector[present_cell].l2_norm();

        }

  }


  // @sect4{PlasticityContactProblem::refine_grid}

  // If you've made it this far into the deal.II tutorial, the following
  // function refining the mesh should not pose any challenges to you
  // any more. It refines the mesh, either globally or using the Kelly
  // error estimator, and if so asked also transfers the solution from
  // the previous to the next mesh. In the latter case, we also need
  // to compute the active set and other quantities again, for which we
  // need the information computed by <code>compute_nonlinear_residual()</code>.
  template <int dim>
  void
  ElastoPlasticProblem<dim>::refine_grid ()
  {
    // ---------------------------------------------------------------
    // Make a field variable for history varibales to be able to
    // transfer the data to the quadrature points of the new mesh
    FE_DGQ<dim> history_fe (1);
    DoFHandler<dim> history_dof_handler (triangulation);
    history_dof_handler.distribute_dofs (history_fe);
    std::vector< std::vector< Vector<double> > >
    history_stress_field (dim, std::vector< Vector<double> >(dim)),
                         local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
                         local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));


    std::vector< std::vector< Vector<double> > >
    history_strain_field (dim, std::vector< Vector<double> >(dim)),
                         local_history_strain_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
                         local_history_strain_fe_values (dim, std::vector< Vector<double> >(dim));

    for (unsigned int i=0; i<dim; ++i)
      for (unsigned int j=0; j<dim; ++j)
        {
          history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
          local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
          local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);

          history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
          local_history_strain_values_at_qpoints[i][j].reinit(quadrature_formula.size());
          local_history_strain_fe_values[i][j].reinit(history_fe.dofs_per_cell);
        }
    FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
                                             quadrature_formula.size());
    FETools::compute_projection_from_quadrature_points_matrix
    (history_fe,
     quadrature_formula, quadrature_formula,
     qpoint_to_dof_matrix);
    typename DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end(),
    dg_cell = history_dof_handler.begin_active();
    for (; cell!=endc; ++cell, ++dg_cell)
      if (cell->is_locally_owned())
        {
          PointHistory<dim> *local_quadrature_points_history
            = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
          Assert (local_quadrature_points_history >=
                  &quadrature_point_history.front(),
                  ExcInternalError());
          Assert (local_quadrature_points_history <
                  &quadrature_point_history.back(),
                  ExcInternalError());
          for (unsigned int i=0; i<dim; ++i)
            for (unsigned int j=0; j<dim; ++j)
              {
                for (unsigned int q=0; q<quadrature_formula.size(); ++q)
                  {
                    local_history_stress_values_at_qpoints[i][j](q)
                      = local_quadrature_points_history[q].old_stress[i][j];

                    local_history_strain_values_at_qpoints[i][j](q)
                      = local_quadrature_points_history[q].old_strain[i][j];
                  }
                qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
                                            local_history_stress_values_at_qpoints[i][j]);
                dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
                                         history_stress_field[i][j]);

                qpoint_to_dof_matrix.vmult (local_history_strain_fe_values[i][j],
                                            local_history_strain_values_at_qpoints[i][j]);
                dg_cell->set_dof_values (local_history_strain_fe_values[i][j],
                                         history_strain_field[i][j]);
              }
        }


    // ---------------------------------------------------------------
    // Refine the mesh
    if (refinement_strategy == RefinementStrategy::refine_global)
      {
        for (typename Triangulation<dim>::active_cell_iterator
             cell = triangulation.begin_active();
             cell != triangulation.end(); ++cell)
          if (cell->is_locally_owned())
            cell->set_refine_flag ();
      }
    else
      {
        const double refine_fraction_cells = .3,
                     coarsen_fraction_cells = .03;
//        const double refine_fraction_cells = .1,
//                     coarsen_fraction_cells = .3;

        parallel::distributed::GridRefinement
        ::refine_and_coarsen_fixed_number(triangulation,
                                          estimated_error_per_cell,
                                          refine_fraction_cells, coarsen_fraction_cells);
      }

    triangulation.prepare_coarsening_and_refinement();

    parallel::distributed::SolutionTransfer<dim,
             TrilinosWrappers::MPI::Vector> solution_transfer(dof_handler);
    solution_transfer.prepare_for_coarsening_and_refinement(solution);


    parallel::distributed::SolutionTransfer<dim,
             TrilinosWrappers::MPI::Vector> incremental_displacement_transfer(dof_handler);
    if (transfer_solution)
      incremental_displacement_transfer.prepare_for_coarsening_and_refinement(incremental_displacement);

    SolutionTransfer<dim, Vector<double> > history_stress_field_transfer0(history_dof_handler),
                     history_stress_field_transfer1(history_dof_handler),
                     history_stress_field_transfer2(history_dof_handler);
    history_stress_field_transfer0.prepare_for_coarsening_and_refinement(history_stress_field[0]);
    if ( dim > 1)
      {
        history_stress_field_transfer1.prepare_for_coarsening_and_refinement(history_stress_field[1]);
      }
    if ( dim == 3)
      {
        history_stress_field_transfer2.prepare_for_coarsening_and_refinement(history_stress_field[2]);
      }

    SolutionTransfer<dim, Vector<double> > history_strain_field_transfer0(history_dof_handler),
                     history_strain_field_transfer1(history_dof_handler),
                     history_strain_field_transfer2(history_dof_handler);
    history_strain_field_transfer0.prepare_for_coarsening_and_refinement(history_strain_field[0]);
    if ( dim > 1)
      {
        history_strain_field_transfer1.prepare_for_coarsening_and_refinement(history_strain_field[1]);
      }
    if ( dim == 3)
      {
        history_strain_field_transfer2.prepare_for_coarsening_and_refinement(history_strain_field[2]);
      }

    triangulation.execute_coarsening_and_refinement();
    pcout << "    Number of active cells:       "
          << triangulation.n_active_cells()
          << std::endl;

    setup_system();
    setup_quadrature_point_history ();


    TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs, mpi_communicator);
//    distributed_solution = solution;
    solution_transfer.interpolate(distributed_solution);
    solution = distributed_solution;

    if (transfer_solution)
      {
        TrilinosWrappers::MPI::Vector distributed_incremental_displacement(locally_owned_dofs, mpi_communicator);
//        distributed_incremental_displacement = incremental_displacement;
        incremental_displacement_transfer.interpolate(distributed_incremental_displacement);
        incremental_displacement = distributed_incremental_displacement;
//        compute_nonlinear_residual(incremental_displacement);
      }

    // ---------------------------------------------------
    history_dof_handler.distribute_dofs (history_fe);
    // stress
    std::vector< std::vector< Vector<double> > >
    distributed_history_stress_field (dim, std::vector< Vector<double> >(dim));
    for (unsigned int i=0; i<dim; ++i)
      for (unsigned int j=0; j<dim; ++j)
        {
          distributed_history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
        }

    history_stress_field_transfer0.interpolate(history_stress_field[0], distributed_history_stress_field[0]);
    if ( dim > 1)
      {
        history_stress_field_transfer1.interpolate(history_stress_field[1], distributed_history_stress_field[1]);
      }
    if ( dim == 3)
      {
        history_stress_field_transfer2.interpolate(history_stress_field[2], distributed_history_stress_field[2]);
      }

    history_stress_field = distributed_history_stress_field;

    // strain
    std::vector< std::vector< Vector<double> > >
    distributed_history_strain_field (dim, std::vector< Vector<double> >(dim));
    for (unsigned int i=0; i<dim; ++i)
      for (unsigned int j=0; j<dim; ++j)
        {
          distributed_history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
        }

    history_strain_field_transfer0.interpolate(history_strain_field[0], distributed_history_strain_field[0]);
    if ( dim > 1)
      {
        history_strain_field_transfer1.interpolate(history_strain_field[1], distributed_history_strain_field[1]);
      }
    if ( dim == 3)
      {
        history_strain_field_transfer2.interpolate(history_strain_field[2], distributed_history_strain_field[2]);
      }

    history_strain_field = distributed_history_strain_field;

    // ---------------------------------------------------------------
    // Transfer the history data to the quadrature points of the new mesh
    // In a final step, we have to get the data back from the now
    // interpolated global field to the quadrature points on the
    // new mesh. The following code will do that:

    FullMatrix<double> dof_to_qpoint_matrix (quadrature_formula.size(),
                                             history_fe.dofs_per_cell);
    FETools::compute_interpolation_to_quadrature_points_matrix
    (history_fe,
     quadrature_formula,
     dof_to_qpoint_matrix);
    cell = dof_handler.begin_active();
    endc = dof_handler.end();
    dg_cell = history_dof_handler.begin_active();
    for (; cell != endc; ++cell, ++dg_cell)
      if (cell->is_locally_owned())
        {
          PointHistory<dim> *local_quadrature_points_history
            = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
          Assert (local_quadrature_points_history >=
                  &quadrature_point_history.front(),
                  ExcInternalError());
          Assert (local_quadrature_points_history <
                  &quadrature_point_history.back(),
                  ExcInternalError());
          for (unsigned int i=0; i<dim; ++i)
            for (unsigned int j=0; j<dim; ++j)
              {
                dg_cell->get_dof_values (history_stress_field[i][j],
                                         local_history_stress_fe_values[i][j]);
                dof_to_qpoint_matrix.vmult (local_history_stress_values_at_qpoints[i][j],
                                            local_history_stress_fe_values[i][j]);

                dg_cell->get_dof_values (history_strain_field[i][j],
                                         local_history_strain_fe_values[i][j]);
                dof_to_qpoint_matrix.vmult (local_history_strain_values_at_qpoints[i][j],
                                            local_history_strain_fe_values[i][j]);
                for (unsigned int q=0; q<quadrature_formula.size(); ++q)
                  {
                    local_quadrature_points_history[q].old_stress[i][j]
                      = local_history_stress_values_at_qpoints[i][j](q);

                    local_quadrature_points_history[q].old_strain[i][j]
                      = local_history_strain_values_at_qpoints[i][j](q);
                  }
              }


        }
  }

  // @sect4{ElastoPlasticProblem::setup_quadrature_point_history}

  // At the beginning of our computations, we needed to set up initial values
  // of the history variables, such as the existing stresses in the material,
  // that we store in each quadrature point. As mentioned above, we use the
  // <code>user_pointer</code> for this that is available in each cell.
  //
  // To put this into larger perspective, we note that if we had previously
  // available stresses in our model (which we assume do not exist for the
  // purpose of this program), then we would need to interpolate the field of
  // preexisting stresses to the quadrature points. Likewise, if we were to
  // simulate elasto-plastic materials with hardening/softening, then we would
  // have to store additional history variables like the present yield stress
  // of the accumulated plastic strains in each quadrature
  // points. Pre-existing hardening or weakening would then be implemented by
  // interpolating these variables in the present function as well.
  template <int dim>
  void ElastoPlasticProblem<dim>::setup_quadrature_point_history ()
  {
    // What we need to do here is to first count how many quadrature points
    // are within the responsibility of this processor. This, of course,
    // equals the number of cells that belong to this processor times the
    // number of quadrature points our quadrature formula has on each cell.
    //
    // For good measure, we also set all user pointers of all cells, whether
    // ours of not, to the null pointer. This way, if we ever access the user
    // pointer of a cell which we should not have accessed, a segmentation
    // fault will let us know that this should not have happened:
    unsigned int our_cells = 0;
    for (typename Triangulation<dim>::active_cell_iterator
         cell = triangulation.begin_active();
         cell != triangulation.end(); ++cell)
      if (cell->is_locally_owned())
        ++our_cells;

    triangulation.clear_user_data();

    // Next, allocate as many quadrature objects as we need. Since the
    // <code>resize</code> function does not actually shrink the amount of
    // allocated memory if the requested new size is smaller than the old
    // size, we resort to a trick to first free all memory, and then
    // reallocate it: we declare an empty vector as a temporary variable and
    // then swap the contents of the old vector and this temporary
    // variable. This makes sure that the
    // <code>quadrature_point_history</code> is now really empty, and we can
    // let the temporary variable that now holds the previous contents of the
    // vector go out of scope and be destroyed. In the next step. we can then
    // re-allocate as many elements as we need, with the vector
    // default-initializing the <code>PointHistory</code> objects, which
    // includes setting the stress variables to zero.
    {
      std::vector<PointHistory<dim> > tmp;
      tmp.swap (quadrature_point_history);
    }
    quadrature_point_history.resize (our_cells *
                                     quadrature_formula.size());

    // Finally loop over all cells again and set the user pointers from the
    // cells that belong to the present processor to point to the first
    // quadrature point objects corresponding to this cell in the vector of
    // such objects:
    unsigned int history_index = 0;
    for (typename Triangulation<dim>::active_cell_iterator
         cell = triangulation.begin_active();
         cell != triangulation.end(); ++cell)
      if (cell->is_locally_owned())
        {
          cell->set_user_pointer (&quadrature_point_history[history_index]);
          history_index += quadrature_formula.size();
        }

    // At the end, for good measure make sure that our count of elements was
    // correct and that we have both used up all objects we allocated
    // previously, and not point to any objects beyond the end of the
    // vector. Such defensive programming strategies are always good checks to
    // avoid accidental errors and to guard against future changes to this
    // function that forget to update all uses of a variable at the same
    // time. Recall that constructs using the <code>Assert</code> macro are
    // optimized away in optimized mode, so do not affect the run time of
    // optimized runs:
    Assert (history_index == quadrature_point_history.size(),
            ExcInternalError());
  }

  // @sect4{ElastoPlasticProblem::update_quadrature_point_history}

  // At the end of each time step, we should have computed an incremental
  // displacement update so that the material in its new configuration
  // accommodates for the difference between the external body and boundary
  // forces applied during this time step minus the forces exerted through
  // preexisting internal stresses. In order to have the preexisting
  // stresses available at the next time step, we therefore have to update the
  // preexisting stresses with the stresses due to the incremental
  // displacement computed during the present time step. Ideally, the
  // resulting sum of internal stresses would exactly counter all external
  // forces. Indeed, a simple experiment can make sure that this is so: if we
  // choose boundary conditions and body forces to be time independent, then
  // the forcing terms (the sum of external forces and internal stresses)
  // should be exactly zero. If you make this experiment, you will realize
  // from the output of the norm of the right hand side in each time step that
  // this is almost the case: it is not exactly zero, since in the first time
  // step the incremental displacement and stress updates were computed
  // relative to the undeformed mesh, which was then deformed. In the second
  // time step, we again compute displacement and stress updates, but this
  // time in the deformed mesh -- there, the resulting updates are very small
  // but not quite zero. This can be iterated, and in each such iteration the
  // residual, i.e. the norm of the right hand side vector, is reduced; if one
  // makes this little experiment, one realizes that the norm of this residual
  // decays exponentially with the number of iterations, and after an initial
  // very rapid decline is reduced by roughly a factor of about 3.5 in each
  // iteration (for one testcase I looked at, other testcases, and other
  // numbers of unknowns change the factor, but not the exponential decay).

  // In a sense, this can then be considered as a quasi-timestepping scheme to
  // resolve the nonlinear problem of solving large-deformation elasticity on
  // a mesh that is moved along in a Lagrangian manner.
  //
  // Another complication is that the existing (old) stresses are defined on
  // the old mesh, which we will move around after updating the stresses. If
  // this mesh update involves rotations of the cell, then we need to also
  // rotate the updated stress, since it was computed relative to the
  // coordinate system of the old cell.
  //
  // Thus, what we need is the following: on each cell which the present
  // processor owns, we need to extract the old stress from the data stored
  // with each quadrature point, compute the stress update, add the two
  // together, and then rotate the result together with the incremental
  // rotation computed from the incremental displacement at the present
  // quadrature point. We will detail these steps below:
  template <int dim>
  void ElastoPlasticProblem<dim>::
  update_quadrature_point_history ()
  {
    // First, set up an <code>FEValues</code> object by which we will evaluate
    // the displacements and the gradients thereof at the
    // quadrature points, together with a vector that will hold this
    // information:
    FEValues<dim> fe_values (fe, quadrature_formula,
                             update_values | update_gradients |
                             update_quadrature_points);

    const unsigned int n_q_points = quadrature_formula.size();

    std::vector<SymmetricTensor<2, dim> > incremental_strain_tensor(n_q_points);
    SymmetricTensor<4, dim> stress_strain_tensor;


    // Then loop over all cells and do the job in the cells that belong to our
    // subdomain:

    typename DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();

    const FEValuesExtractors::Vector displacement(0);

    for (;  cell != endc; ++cell)
      if (cell->is_locally_owned())
        {
          // Next, get a pointer to the quadrature point history data local to
          // the present cell, and, as a defensive measure, make sure that
          // this pointer is within the bounds of the global array:
          PointHistory<dim> *local_quadrature_points_history
            = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
          Assert (local_quadrature_points_history >=
                  &quadrature_point_history.front(),
                  ExcInternalError());
          Assert (local_quadrature_points_history <
                  &quadrature_point_history.back(),
                  ExcInternalError());

          // Then initialize the <code>FEValues</code> object on the present
          // cell, and extract the strains of the displacement at the
          // quadrature points
          fe_values.reinit (cell);
          fe_values[displacement].get_function_symmetric_gradients(incremental_displacement,
                                                                   incremental_strain_tensor);

          // Then loop over the quadrature points of this cell:
          for (unsigned int q=0; q<quadrature_formula.size(); ++q)
            {
              local_quadrature_points_history[q].old_strain +=
                incremental_strain_tensor[q];

              constitutive_law.get_stress_strain_tensor(local_quadrature_points_history[q].old_strain,
                                                        stress_strain_tensor);

              // The result of these operations is then written back into
              // the original place:
              local_quadrature_points_history[q].old_stress
                = stress_strain_tensor *  local_quadrature_points_history[q].old_strain;

              local_quadrature_points_history[q].point
                = fe_values.get_quadrature_points ()[q];
            }
        }
  }


  // @sect4{PlasticityContactProblem::move_mesh}

  // The remaining three functions before we get to <code>run()</code>
  // have to do with generating output. The following one is an attempt
  // at showing the deformed body in its deformed configuration. To this
  // end, this function takes a displacement vector field and moves every
  // vertex of the (local part) of the mesh by the previously computed
  // displacement. We will call this function with the current
  // displacement field before we generate graphical output, and we will
  // call it again after generating graphical output with the negative
  // displacement field to undo the changes to the mesh so made.
  //
  // The function itself is pretty straightforward. All we have to do
  // is keep track which vertices we have already touched, as we
  // encounter the same vertices multiple times as we loop over cells.
  template <int dim>
  void
  ElastoPlasticProblem<dim>::
  move_mesh (const TrilinosWrappers::MPI::Vector &displacement) const
  {
    std::vector<bool> vertex_touched(triangulation.n_vertices(), false);

    for (typename DoFHandler<dim>::active_cell_iterator cell =
           dof_handler.begin_active();
         cell != dof_handler.end(); ++cell)
      if (cell->is_locally_owned())
        for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
          if (vertex_touched[cell->vertex_index(v)] == false)
            {
              vertex_touched[cell->vertex_index(v)] = true;

              Point<dim> vertex_displacement;
              for (unsigned int d = 0; d < dim; ++d)
                vertex_displacement[d] = displacement(cell->vertex_dof_index(v, d));

              cell->vertex(v) += vertex_displacement;
            }
  }



  // @sect4{PlasticityContactProblem::output_results}

  // Next is the function we use to actually generate graphical output. The
  // function is a bit tedious, but not actually particularly complicated.
  // It moves the mesh at the top (and moves it back at the end), then
  // computes the contact forces along the contact surface. We can do
  // so (as shown in the accompanying paper) by taking the untreated
  // residual vector and identifying which degrees of freedom
  // correspond to those with contact by asking whether they have an
  // inhomogeneous constraints associated with them. As always, we need
  // to be mindful that we can only write into completely distributed
  // vectors (i.e., vectors without ghost elements) but that when we
  // want to generate output, we need vectors that do indeed have
  // ghost entries for all locally relevant degrees of freedom.
  template <int dim>
  void
  ElastoPlasticProblem<dim>::output_results (const std::string &filename_base)
  {
    TimerOutput::Scope t(computing_timer, "Graphical output");

    pcout << "      Writing graphical output... " << std::flush;

    TrilinosWrappers::MPI::Vector magnified_solution(solution);

    const double magnified_factor = 3;
    magnified_solution *= magnified_factor;

    move_mesh(magnified_solution);

    DataOut<dim> data_out;

    data_out.attach_dof_handler(dof_handler);

    //
    const std::vector<DataComponentInterpretation::DataComponentInterpretation>
    data_component_interpretation(dim, DataComponentInterpretation::component_is_part_of_vector);
    data_out.add_data_vector(solution,
                             std::vector<std::string> (dim, "displacement"),
                             DataOut<dim>::type_dof_data, data_component_interpretation);

    //
    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:
        AssertThrow (false, ExcNotImplemented());
      }

    data_out.add_data_vector (solution, solution_names);


    //
    Vector<float> subdomain(triangulation.n_active_cells());
    for (unsigned int i = 0; i < subdomain.size(); ++i)
      subdomain(i) = triangulation.locally_owned_subdomain();
    data_out.add_data_vector(subdomain, "subdomain");

    //
    data_out.add_data_vector(fraction_of_plastic_q_points_per_cell,
                             "fraction_of_plastic_q_points");

    //
    data_out.build_patches();

    // In the remainder of the function, we generate one VTU file on
    // every processor, indexed by the subdomain id of this processor.
    // On the first processor, we then also create a <code>.pvtu</code>
    // file that indexes <i>all</i> of the VTU files so that the entire
    // set of output files can be read at once. These <code>.pvtu</code>
    // are used by Paraview to describe an entire parallel computation's
    // output files. We then do the same again for the competitor of
    // Paraview, the Visit visualization program, by creating a matching
    // <code>.visit</code> file.
    const std::string filename =
      (output_dir + filename_base + "-"
       + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));

    std::ofstream output_vtu((filename + ".vtu").c_str());
    data_out.write_vtu(output_vtu);
    pcout << output_dir + filename_base << ".pvtu" << std::endl;


    if (this_mpi_process == 0)
      {
        std::vector<std::string> filenames;
        for (unsigned int i = 0; i < n_mpi_processes; ++i)
          filenames.push_back(filename_base + "-" +
                              Utilities::int_to_string(i, 4) +
                              ".vtu");

        std::ofstream pvtu_master_output((output_dir + filename_base + ".pvtu").c_str());
        data_out.write_pvtu_record(pvtu_master_output, filenames);

        std::ofstream visit_master_output((output_dir + filename_base + ".visit").c_str());
        data_out.write_pvtu_record(visit_master_output, filenames);

        // produce eps files for mesh illustration
        std::ofstream output_eps((filename + ".eps").c_str());
        GridOut grid_out;
        grid_out.write_eps(triangulation, output_eps);
      }

    // Extrapolate the stresses from Gauss point to the nodes
    SymmetricTensor<2, dim> stress_at_qpoint;

    FE_DGQ<dim> history_fe (1);
    DoFHandler<dim> history_dof_handler (triangulation);
    history_dof_handler.distribute_dofs (history_fe);
    std::vector< std::vector< Vector<double> > >
    history_stress_field (dim, std::vector< Vector<double> >(dim)),
                         local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
                         local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
    for (unsigned int i=0; i<dim; ++i)
      for (unsigned int j=0; j<dim; ++j)
        {
          history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
          local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
          local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
        }

    Vector<double>  VM_stress_field (history_dof_handler.n_dofs()),
           local_VM_stress_values_at_qpoints (quadrature_formula.size()),
           local_VM_stress_fe_values (history_fe.dofs_per_cell);

    FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
                                             quadrature_formula.size());
    FETools::compute_projection_from_quadrature_points_matrix
    (history_fe,
     quadrature_formula, quadrature_formula,
     qpoint_to_dof_matrix);

    typename DoFHandler<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end(),
    dg_cell = history_dof_handler.begin_active();

    const FEValuesExtractors::Vector displacement(0);

    for (; cell!=endc; ++cell, ++dg_cell)
      if (cell->is_locally_owned())
        {
          PointHistory<dim> *local_quadrature_points_history
            = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
          Assert (local_quadrature_points_history >=
                  &quadrature_point_history.front(),
                  ExcInternalError());
          Assert (local_quadrature_points_history <
                  &quadrature_point_history.back(),
                  ExcInternalError());

          // Then loop over the quadrature points of this cell:
          for (unsigned int q=0; q<quadrature_formula.size(); ++q)
            {
              stress_at_qpoint = local_quadrature_points_history[q].old_stress;

              for (unsigned int i=0; i<dim; ++i)
                for (unsigned int j=i; j<dim; ++j)
                  {
                    local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
                  }

              local_VM_stress_values_at_qpoints(q) = Evaluation::get_von_Mises_stress(stress_at_qpoint);

            }


          for (unsigned int i=0; i<dim; ++i)
            for (unsigned int j=i; j<dim; ++j)
              {
                qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
                                            local_history_stress_values_at_qpoints[i][j]);
                dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
                                         history_stress_field[i][j]);
              }

          qpoint_to_dof_matrix.vmult (local_VM_stress_fe_values,
                                      local_VM_stress_values_at_qpoints);
          dg_cell->set_dof_values (local_VM_stress_fe_values,
                                   VM_stress_field);


        }

    // Save stresses on nodes by nodal averaging
    // construct a DoFHandler object based on FE_Q with 1 degree of freedom
    // in order to compute stresses on nodes (by applying nodal averaging)
    // Therefore, each vertex has one degree of freedom
    FE_Q<dim>          fe_1 (1);
    DoFHandler<dim>    dof_handler_1 (triangulation);
    dof_handler_1.distribute_dofs (fe_1);

    AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(),
                ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices()));

    std::vector< std::vector< Vector<double> > >
    history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
    for (unsigned int i=0; i<dim; ++i)
      for (unsigned int j=0; j<dim; ++j)
        {
          history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
        }

    Vector<double>  VM_stress_on_vertices (dof_handler_1.n_dofs()),
           counter_on_vertices (dof_handler_1.n_dofs());
    VM_stress_on_vertices = 0;
    counter_on_vertices = 0;

    cell = dof_handler.begin_active();
    dg_cell = history_dof_handler.begin_active();
    typename DoFHandler<dim>::active_cell_iterator
    cell_1 = dof_handler_1.begin_active();
    for (; cell!=endc; ++cell, ++dg_cell, ++cell_1)
      if (cell->is_locally_owned())
        {
          dg_cell->get_dof_values (VM_stress_field,
                                   local_VM_stress_fe_values);

          for (unsigned int i=0; i<dim; ++i)
            for (unsigned int j=0; j<dim; ++j)
              {
                dg_cell->get_dof_values (history_stress_field[i][j],
                                         local_history_stress_fe_values[i][j]);
              }

          for  (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
            {
              types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0);

              // begin check
              //            Point<dim> point1, point2;
              //            point1 = cell_1->vertex(v);
              //            point2 = dg_cell->vertex(v);
              //            AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError());
              // end check

              counter_on_vertices (dof_1_vertex) += 1;

              VM_stress_on_vertices (dof_1_vertex) += local_VM_stress_fe_values (v);

              for (unsigned int i=0; i<dim; ++i)
                for (unsigned int j=0; j<dim; ++j)
                  {
                    history_stress_on_vertices[i][j](dof_1_vertex) +=
                      local_history_stress_fe_values[i][j](v);
                  }

            }
        }

    for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
      {
        VM_stress_on_vertices(id) /= counter_on_vertices(id);

        for (unsigned int i=0; i<dim; ++i)
          for (unsigned int j=0; j<dim; ++j)
            {
              history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
            }
      }

    // Save figures of stresses
    if (show_stresses)
      {
        {
          DataOut<dim>  data_out;
          data_out.attach_dof_handler (history_dof_handler);


          data_out.add_data_vector (history_stress_field[0][0], "stress_xx");
          data_out.add_data_vector (history_stress_field[1][1], "stress_yy");
          data_out.add_data_vector (history_stress_field[0][1], "stress_xy");
          data_out.add_data_vector (VM_stress_field, "Von_Mises_stress");

          if (dim == 3)
            {
              data_out.add_data_vector (history_stress_field[0][2], "stress_xz");
              data_out.add_data_vector (history_stress_field[1][2], "stress_yz");
              data_out.add_data_vector (history_stress_field[2][2], "stress_zz");
            }

          data_out.build_patches ();

          const std::string filename_base_stress = ("stress-" + filename_base);

          const std::string filename =
            (output_dir + filename_base_stress + "-"
             + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));

          std::ofstream output_vtu((filename + ".vtu").c_str());
          data_out.write_vtu(output_vtu);
          pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;

          if (this_mpi_process == 0)
            {
              std::vector<std::string> filenames;
              for (unsigned int i = 0; i < n_mpi_processes; ++i)
                filenames.push_back(filename_base_stress + "-" +
                                    Utilities::int_to_string(i, 4) +
                                    ".vtu");

              std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
              data_out.write_pvtu_record(pvtu_master_output, filenames);

              std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
              data_out.write_pvtu_record(visit_master_output, filenames);
            }


        }

        {
          DataOut<dim>  data_out;
          data_out.attach_dof_handler (dof_handler_1);


          data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_xx_averaged");
          data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_yy_averaged");
          data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_xy_averaged");
          data_out.add_data_vector (VM_stress_on_vertices, "Von_Mises_stress_averaged");

          if (dim == 3)
            {
              data_out.add_data_vector (history_stress_on_vertices[0][2], "stress_xz_averaged");
              data_out.add_data_vector (history_stress_on_vertices[1][2], "stress_yz_averaged");
              data_out.add_data_vector (history_stress_on_vertices[2][2], "stress_zz_averaged");
            }

          data_out.build_patches ();

          const std::string filename_base_stress = ("averaged-stress-" + filename_base);

          const std::string filename =
            (output_dir + filename_base_stress + "-"
             + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));

          std::ofstream output_vtu((filename + ".vtu").c_str());
          data_out.write_vtu(output_vtu);
          pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;

          if (this_mpi_process == 0)
            {
              std::vector<std::string> filenames;
              for (unsigned int i = 0; i < n_mpi_processes; ++i)
                filenames.push_back(filename_base_stress + "-" +
                                    Utilities::int_to_string(i, 4) +
                                    ".vtu");

              std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
              data_out.write_pvtu_record(pvtu_master_output, filenames);

              std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
              data_out.write_pvtu_record(visit_master_output, filenames);
            }


        }
        // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      }

    magnified_solution *= -1;
    move_mesh(magnified_solution);

    // Timoshenko beam
    if (base_mesh == "Timoshenko beam")
      {
        const double length = .48,
                     depth  = .12;

        Point<dim> intersted_point(length, -depth/2);
        Point<dim> vertex_displacement;
        bool vertex_found = false;

        for (typename DoFHandler<dim>::active_cell_iterator cell =
               dof_handler.begin_active();
             cell != dof_handler.end(); ++cell)
          if (cell->is_locally_owned() && !vertex_found)
            for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
              if ( std::fabs(cell->vertex(v)[0] - intersted_point[0])<1e-6 &&
                   std::fabs(cell->vertex(v)[1] - intersted_point[1])<1e-6)
                {
                  vertex_found = true;

                  for (unsigned int d = 0; d < dim; ++d)
                    vertex_displacement[d] = solution(cell->vertex_dof_index(v, d));

                  break;
                }

        pcout << "   Number of active cells: "
              << triangulation.n_global_active_cells() << std::endl
              << "   Number of degrees of freedom: " << dof_handler.n_dofs()
              << std::endl;

        AssertThrow(vertex_found, ExcInternalError());
        std::cout << "Displacement at the point (" << intersted_point[0]
                  << ", " << intersted_point[1] << ") is "
                  << "(" << vertex_displacement[0]
                  << ", " << vertex_displacement[1] << ").\n";

        Vector<double> vertex_exact_displacement(dim);
        EquationData::IncrementalBoundaryValues<dim> incremental_boundary_values(present_time, end_time);
        incremental_boundary_values.vector_value (intersted_point, vertex_exact_displacement);

        std::cout << "Exact displacement at the point (" << intersted_point[0]
                  << ", " << intersted_point[1] << ") is "
                  << "(" << vertex_exact_displacement[0]
                  << ", " << vertex_exact_displacement[1] << ").\n\n";

      }
    else if (base_mesh == "Thick_tube_internal_pressure")
      {
        const double pressure (0.6*2.4e8),
              inner_radius (.1);
//      const double pressure (1.94e8),
//                   inner_radius (.1);


        // Plane stress
//      const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu)))));
        // 3d and plane strain
        const double mu (e_modulus / (2 * (1 + nu)));

        const Point<dim> point_A(inner_radius, 0.);
        Vector<double>   disp_A(dim);

        // make a non-parallel copy of solution
        Vector<double> copy_solution(solution);

        Evaluation::PointValuesEvaluation<dim> point_values_evaluation(point_A);

        point_values_evaluation.compute (dof_handler, copy_solution, disp_A);

        table_results.add_value("time step", timestep_no);
        table_results.add_value("Cells", triangulation.n_global_active_cells());
        table_results.add_value("DoFs", dof_handler.n_dofs());
        table_results.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
        table_results.add_value("4*mu*u_A/(sigma_0*a)", 4*mu*disp_A(0)/(sigma_0*inner_radius));

        // Compute stresses in the POLAR coordinates, 1- save it on Gauss points,
        // 2- extrapolate them to nodes and taking their avarages (nodal avaraging)
        AssertThrow (dim == 2, ExcNotImplemented());

        // we define a rotation matrix to be able to transform the stress
        // from the Cartesian coordinate to the polar coordinate
        Tensor<2, dim> rotation_matrix; // [cos sin; -sin cos]    , sigma_r = rot * sigma * rot^T

        FEValues<dim> fe_values (fe, quadrature_formula, update_quadrature_points |
                                 update_values | update_gradients);

        const unsigned int n_q_points = quadrature_formula.size();

        std::vector<SymmetricTensor<2, dim> > strain_tensor(n_q_points);
        SymmetricTensor<4, dim> stress_strain_tensor;
        Tensor<2, dim>  stress_at_qpoint;

        FE_DGQ<dim> history_fe (1);
        DoFHandler<dim> history_dof_handler (triangulation);
        history_dof_handler.distribute_dofs (history_fe);
        std::vector< std::vector< Vector<double> > >
        history_stress_field (dim, std::vector< Vector<double> >(dim)),
                             local_history_stress_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
                             local_history_stress_fe_values (dim, std::vector< Vector<double> >(dim));
        for (unsigned int i=0; i<dim; ++i)
          for (unsigned int j=0; j<dim; ++j)
            {
              history_stress_field[i][j].reinit(history_dof_handler.n_dofs());
              local_history_stress_values_at_qpoints[i][j].reinit(quadrature_formula.size());
              local_history_stress_fe_values[i][j].reinit(history_fe.dofs_per_cell);
            }

        FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
                                                 quadrature_formula.size());
        FETools::compute_projection_from_quadrature_points_matrix
        (history_fe,
         quadrature_formula, quadrature_formula,
         qpoint_to_dof_matrix);

        typename DoFHandler<dim>::active_cell_iterator
        cell = dof_handler.begin_active(),
        endc = dof_handler.end(),
        dg_cell = history_dof_handler.begin_active();

        const FEValuesExtractors::Vector displacement(0);

        for (; cell!=endc; ++cell, ++dg_cell)
          if (cell->is_locally_owned())
            {
              PointHistory<dim> *local_quadrature_points_history
                = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
              Assert (local_quadrature_points_history >=
                      &quadrature_point_history.front(),
                      ExcInternalError());
              Assert (local_quadrature_points_history <
                      &quadrature_point_history.back(),
                      ExcInternalError());

              // Then loop over the quadrature points of this cell:
              for (unsigned int q=0; q<quadrature_formula.size(); ++q)
                {
                  stress_at_qpoint = local_quadrature_points_history[q].old_stress;

                  // transform the stress from the Cartesian coordinate to the polar coordinate
                  const Point<dim> point = local_quadrature_points_history[q].point;
                  const double theta = std::atan2(point(1),point(0));

                  // rotation matrix
                  rotation_matrix[0][0] = std::cos(theta);
                  rotation_matrix[0][1] = std::sin(theta);
                  rotation_matrix[1][0] = -std::sin(theta);
                  rotation_matrix[1][1] = std::cos(theta);

                  // stress in polar coordinate
                  stress_at_qpoint = rotation_matrix * stress_at_qpoint * transpose(rotation_matrix);

                  for (unsigned int i=0; i<dim; ++i)
                    for (unsigned int j=i; j<dim; ++j)
                      {
                        local_history_stress_values_at_qpoints[i][j](q) = stress_at_qpoint[i][j];
                      }

                }


              for (unsigned int i=0; i<dim; ++i)
                for (unsigned int j=i; j<dim; ++j)
                  {
                    qpoint_to_dof_matrix.vmult (local_history_stress_fe_values[i][j],
                                                local_history_stress_values_at_qpoints[i][j]);
                    dg_cell->set_dof_values (local_history_stress_fe_values[i][j],
                                             history_stress_field[i][j]);
                  }

            }

        {
          DataOut<dim>  data_out;
          data_out.attach_dof_handler (history_dof_handler);


          data_out.add_data_vector (history_stress_field[0][0], "stress_rr");
          data_out.add_data_vector (history_stress_field[1][1], "stress_tt");
          data_out.add_data_vector (history_stress_field[0][1], "stress_rt");

          data_out.build_patches ();

          const std::string filename_base_stress = ("stress-polar-" + filename_base);

          const std::string filename =
            (output_dir + filename_base_stress + "-"
             + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));

          std::ofstream output_vtu((filename + ".vtu").c_str());
          data_out.write_vtu(output_vtu);
          pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;

          if (this_mpi_process == 0)
            {
              std::vector<std::string> filenames;
              for (unsigned int i = 0; i < n_mpi_processes; ++i)
                filenames.push_back(filename_base_stress + "-" +
                                    Utilities::int_to_string(i, 4) +
                                    ".vtu");

              std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
              data_out.write_pvtu_record(pvtu_master_output, filenames);

              std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
              data_out.write_pvtu_record(visit_master_output, filenames);
            }


        }

        // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        // construct a DoFHandler object based on FE_Q with 1 degree of freedom
        // in order to compute stresses on nodes (by applying nodal averaging)
        // Therefore, each vertex has one degree of freedom
        FE_Q<dim>          fe_1 (1);
        DoFHandler<dim>    dof_handler_1 (triangulation);
        dof_handler_1.distribute_dofs (fe_1);

        AssertThrow(dof_handler_1.n_dofs() == triangulation.n_vertices(),
                    ExcDimensionMismatch(dof_handler_1.n_dofs(),triangulation.n_vertices()));

        std::vector< std::vector< Vector<double> > >
        history_stress_on_vertices (dim, std::vector< Vector<double> >(dim));
        for (unsigned int i=0; i<dim; ++i)
          for (unsigned int j=0; j<dim; ++j)
            {
              history_stress_on_vertices[i][j].reinit(dof_handler_1.n_dofs());
            }

        Vector<double>  counter_on_vertices (dof_handler_1.n_dofs());
        counter_on_vertices = 0;

        cell = dof_handler.begin_active();
        dg_cell = history_dof_handler.begin_active();
        typename DoFHandler<dim>::active_cell_iterator
        cell_1 = dof_handler_1.begin_active();
        for (; cell!=endc; ++cell, ++dg_cell, ++cell_1)
          if (cell->is_locally_owned())
            {

              for (unsigned int i=0; i<dim; ++i)
                for (unsigned int j=0; j<dim; ++j)
                  {
                    dg_cell->get_dof_values (history_stress_field[i][j],
                                             local_history_stress_fe_values[i][j]);
                  }

              for  (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
                {
                  types::global_dof_index dof_1_vertex = cell_1->vertex_dof_index(v, 0);

                  // begin check
//            Point<dim> point1, point2;
//            point1 = cell_1->vertex(v);
//            point2 = dg_cell->vertex(v);
//            AssertThrow(point1.distance(point2) < cell->diameter()*1e-8, ExcInternalError());
                  // end check

                  counter_on_vertices (dof_1_vertex) += 1;

                  for (unsigned int i=0; i<dim; ++i)
                    for (unsigned int j=0; j<dim; ++j)
                      {
                        history_stress_on_vertices[i][j](dof_1_vertex) +=
                          local_history_stress_fe_values[i][j](v);
                      }

                }
            }

        for (unsigned int id=0; id<dof_handler_1.n_dofs(); ++id)
          {
            for (unsigned int i=0; i<dim; ++i)
              for (unsigned int j=0; j<dim; ++j)
                {
                  history_stress_on_vertices[i][j](id) /= counter_on_vertices(id);
                }
          }


        {
          DataOut<dim>  data_out;
          data_out.attach_dof_handler (dof_handler_1);


          data_out.add_data_vector (history_stress_on_vertices[0][0], "stress_rr_averaged");
          data_out.add_data_vector (history_stress_on_vertices[1][1], "stress_tt_averaged");
          data_out.add_data_vector (history_stress_on_vertices[0][1], "stress_rt_averaged");

          data_out.build_patches ();

          const std::string filename_base_stress = ("averaged-stress-polar-" + filename_base);

          const std::string filename =
            (output_dir + filename_base_stress + "-"
             + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));

          std::ofstream output_vtu((filename + ".vtu").c_str());
          data_out.write_vtu(output_vtu);
          pcout << output_dir + filename_base_stress << ".pvtu" << std::endl;

          if (this_mpi_process == 0)
            {
              std::vector<std::string> filenames;
              for (unsigned int i = 0; i < n_mpi_processes; ++i)
                filenames.push_back(filename_base_stress + "-" +
                                    Utilities::int_to_string(i, 4) +
                                    ".vtu");

              std::ofstream pvtu_master_output((output_dir + filename_base_stress + ".pvtu").c_str());
              data_out.write_pvtu_record(pvtu_master_output, filenames);

              std::ofstream visit_master_output((output_dir + filename_base_stress + ".visit").c_str());
              data_out.write_pvtu_record(visit_master_output, filenames);
            }


        }
        // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

        if ( std::abs( (present_time/end_time)*(pressure/sigma_0) - 0.6 ) <
             .501*(present_timestep/end_time)*(pressure/sigma_0) )
          {

            // table_results_2: presenting the stress_rr and stress_tt on the nodes of bottom edge
            const unsigned int face_id = 3;

            std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);

            bool evaluation_face_found = false;

            typename DoFHandler<dim>::active_cell_iterator
            cell = dof_handler.begin_active(),
            endc = dof_handler.end(),
            cell_1 = dof_handler_1.begin_active();
            for (; cell!=endc; ++cell, ++cell_1)
              if (cell->is_locally_owned())
                {
                  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
                    {
                      if (cell->face(face)->at_boundary()
                          &&
                          cell->face(face)->boundary_id() == face_id)
                        {
                          if (!evaluation_face_found)
                            {
                              evaluation_face_found = true;
                            }


                          for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
                            {
                              types::global_dof_index dof_1_vertex =
                                cell_1->face(face)->vertex_dof_index(v, 0);
                              if (!vertices_found[dof_1_vertex])
                                {

                                  const Point<dim> vertex_coordinate = cell_1->face(face)->vertex(v);

                                  table_results_2.add_value("x coordinate", vertex_coordinate[0]);
                                  table_results_2.add_value("stress_rr", history_stress_on_vertices[0][0](dof_1_vertex));
                                  table_results_2.add_value("stress_tt", history_stress_on_vertices[1][1](dof_1_vertex));
                                  table_results_2.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);

                                  vertices_found[dof_1_vertex] = true;
                                }
                            }

                        }
                    }

                }

            AssertThrow(evaluation_face_found, ExcInternalError());

            // table_results_3: presenting the mean stress_rr of the nodes on the inner radius
            const unsigned int face_id_2 = 0;

            Tensor<2, dim> stress_node,
                   mean_stress_polar;
            mean_stress_polar = 0;

            std::vector<bool> vertices_found_2 (dof_handler_1.n_dofs(), false);
            unsigned int no_vertices_found = 0;

            evaluation_face_found = false;

            cell = dof_handler.begin_active(),
            endc = dof_handler.end(),
            cell_1 = dof_handler_1.begin_active();
            for (; cell!=endc; ++cell, ++cell_1)
              if (cell->is_locally_owned())
                {
                  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
                    {
                      if (cell->face(face)->at_boundary()
                          &&
                          cell->face(face)->boundary_id() == face_id_2)
                        {
                          if (!evaluation_face_found)
                            {
                              evaluation_face_found = true;
                            }


                          for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
                            {
                              types::global_dof_index dof_1_vertex =
                                cell_1->face(face)->vertex_dof_index(v, 0);
                              if (!vertices_found_2[dof_1_vertex])
                                {
                                  for (unsigned int ir=0; ir<dim; ++ir)
                                    for (unsigned int ic=0; ic<dim; ++ic)
                                      stress_node[ir][ic] = history_stress_on_vertices[ir][ic](dof_1_vertex);

                                  mean_stress_polar += stress_node;

                                  vertices_found_2[dof_1_vertex] = true;
                                  ++no_vertices_found;
                                }
                            }

                        }
                    }

                }

            AssertThrow(evaluation_face_found, ExcInternalError());

            mean_stress_polar /= no_vertices_found;

            table_results_3.add_value("time step", timestep_no);
            table_results_3.add_value("pressure/sigma_0", (pressure*present_time/end_time)/sigma_0);
            table_results_3.add_value("Cells", triangulation.n_global_active_cells());
            table_results_3.add_value("DoFs", dof_handler.n_dofs());
            table_results_3.add_value("radius", inner_radius);
            table_results_3.add_value("mean stress_rr", mean_stress_polar[0][0]);
            table_results_3.add_value("mean stress_tt", mean_stress_polar[1][1]);


          }


      }
    else if (base_mesh == "Perforated_strip_tension")
      {
        const double imposed_displacement (0.00055),
              inner_radius (0.05);

        // Plane stress
//      const double mu (((e_modulus*(1+2*nu)) / (std::pow((1+nu),2))) / (2 * (1 + (nu / (1+nu)))));
        // 3d and plane strain

        // table_results: Demonstrates the result of displacement at the top left corner versus imposed tension
        /*
        {
          const Point<dim> point_C(0., height);
          Vector<double>   disp_C(dim);

          // make a non-parallel copy of solution
          Vector<double> copy_solution(solution);

          typename Evaluation::PointValuesEvaluation<dim>::
          PointValuesEvaluation point_values_evaluation(point_C);

          point_values_evaluation.compute (dof_handler, copy_solution, disp_C);

          table_results.add_value("time step", timestep_no);
          table_results.add_value("Cells", triangulation.n_global_active_cells());
          table_results.add_value("DoFs", dof_handler.n_dofs());
          table_results.add_value("4*mu*u_C/(sigma_0*r)", 4*mu*disp_C(1)/(sigma_0*inner_radius));
        }
        */

        // compute average sigma_yy on the bottom edge
        double stress_yy_av;
        {
          stress_yy_av = 0;
          const unsigned int face_id = 1;

          std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);
          unsigned int no_vertices_in_face = 0;

          bool evaluation_face_found = false;

          typename DoFHandler<dim>::active_cell_iterator
          cell = dof_handler.begin_active(),
          endc = dof_handler.end(),
          cell_1 = dof_handler_1.begin_active();
          for (; cell!=endc; ++cell, ++cell_1)
            if (cell->is_locally_owned())
              {
                for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
                  {
                    if (cell->face(face)->at_boundary()
                        &&
                        cell->face(face)->boundary_id() == face_id)
                      {
                        if (!evaluation_face_found)
                          {
                            evaluation_face_found = true;
                          }


                        for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
                          {
                            types::global_dof_index dof_1_vertex =
                              cell_1->face(face)->vertex_dof_index(v, 0);
                            if (!vertices_found[dof_1_vertex])
                              {
                                stress_yy_av += history_stress_on_vertices[1][1](dof_1_vertex);
                                ++no_vertices_in_face;

                                vertices_found[dof_1_vertex] = true;
                              }
                          }

                      }
                  }

              }

          AssertThrow(evaluation_face_found, ExcInternalError());

          stress_yy_av /= no_vertices_in_face;

        }

        // table_results_2: Demonstrate the stress_yy on the nodes of bottom edge

//      if ( std::abs( (stress_yy_av/sigma_0) - .91 ) < .2 )
        if ( (timestep_no) % 19 == 0 )
//      if ( true )
          {
            const unsigned int face_id = 1;

            std::vector<bool> vertices_found (dof_handler_1.n_dofs(), false);

            bool evaluation_face_found = false;

            typename DoFHandler<dim>::active_cell_iterator
            cell = dof_handler.begin_active(),
            endc = dof_handler.end(),
            cell_1 = dof_handler_1.begin_active();
            for (; cell!=endc; ++cell, ++cell_1)
              if (cell->is_locally_owned())
                {
                  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
                    {
                      if (cell->face(face)->at_boundary()
                          &&
                          cell->face(face)->boundary_id() == face_id)
                        {
                          if (!evaluation_face_found)
                            {
                              evaluation_face_found = true;
                            }


                          for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
                            {
                              types::global_dof_index dof_1_vertex =
                                cell_1->face(face)->vertex_dof_index(v, 0);

                              const Point<dim> vertex_coordinate = cell_1->face(face)->vertex(v);

                              if (!vertices_found[dof_1_vertex] && std::abs(vertex_coordinate[2])<1.e-8)
                                {
                                  table_results_2.add_value("x", vertex_coordinate[0]);
                                  table_results_2.add_value("x/r", vertex_coordinate[0]/inner_radius);
                                  table_results_2.add_value("stress_xx/sigma_0", history_stress_on_vertices[0][0](dof_1_vertex)/sigma_0);
                                  table_results_2.add_value("stress_yy/sigma_0", history_stress_on_vertices[1][1](dof_1_vertex)/sigma_0);
                                  table_results_2.add_value("stress_yy_av/sigma_0", stress_yy_av/sigma_0);
                                  table_results_2.add_value("Imposed u_y", (imposed_displacement*present_time/end_time));

                                  vertices_found[dof_1_vertex] = true;
                                }
                            }

                        }
                    }

                }

            AssertThrow(evaluation_face_found, ExcInternalError());

          }

        // table_results_3: Demonstrate the Stress_mean (average tensile stress)
        //  on the bottom edge versus epsilon_yy on the bottom left corner
        {
          double strain_yy_A = 0.;

          // compute strain_yy_A
          // Since the point A is the node on the bottom left corner,
          // we need to work just with one element
          {
            const Point<dim> point_A(inner_radius, 0, 0);

            Vector<double>  local_strain_yy_values_at_qpoints (quadrature_formula.size()),
                   local_strain_yy_fe_values (history_fe.dofs_per_cell);

            SymmetricTensor<2, dim> strain_at_qpoint;

            typename DoFHandler<dim>::active_cell_iterator
            cell = dof_handler.begin_active(),
            endc = dof_handler.end(),
            dg_cell = history_dof_handler.begin_active();

            bool cell_found = false;

            for (; cell!=endc; ++cell, ++dg_cell)
              if (cell->is_locally_owned() && !cell_found)
                {
                  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
                    if ( std::fabs(cell->vertex(v)[0] - point_A[0])<1e-6 &&
                         std::fabs(cell->vertex(v)[1] - point_A[1])<1e-6 &&
                         std::fabs(cell->vertex(v)[2] - point_A[2])<1e-6)
                      {
                        PointHistory<dim> *local_quadrature_points_history
                          = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
                        Assert (local_quadrature_points_history >=
                                &quadrature_point_history.front(),
                                ExcInternalError());
                        Assert (local_quadrature_points_history <
                                &quadrature_point_history.back(),
                                ExcInternalError());

                        // Then loop over the quadrature points of this cell:
                        for (unsigned int q=0; q<quadrature_formula.size(); ++q)
                          {
                            strain_at_qpoint = local_quadrature_points_history[q].old_strain;

                            local_strain_yy_values_at_qpoints(q) = strain_at_qpoint[1][1];
                          }

                        qpoint_to_dof_matrix.vmult (local_strain_yy_fe_values,
                                                    local_strain_yy_values_at_qpoints);

                        strain_yy_A = local_strain_yy_fe_values (v);

                        cell_found = true;
                        break;
                      }

                }

          }

          table_results_3.add_value("time step", timestep_no);
          table_results_3.add_value("Cells", triangulation.n_global_active_cells());
          table_results_3.add_value("DoFs", dof_handler.n_dofs());
          table_results_3.add_value("Imposed u_y", (imposed_displacement*present_time/end_time));
          table_results_3.add_value("mean_tensile_stress/sigma_0", stress_yy_av/sigma_0);
          table_results_3.add_value("E*strain_yy-A/sigma_0", e_modulus*strain_yy_A/sigma_0);

        }


        if (std::abs(present_time-end_time) < 1.e-7)
          {
            table_results_2.set_precision("Imposed u_y", 6);
            table_results_3.set_precision("Imposed u_y", 6);
          }

      }
    else if (base_mesh == "Cantiliver_beam_3d")
      {
        const double pressure (6e6),
              length (.7),
              height (200e-3);

        // table_results: Demonstrates the result of displacement at the top front point, Point A
        {
          // Quantity of interest:
          // displacement at Point A (x=0, y=height/2, z=length)

          const Point<dim> point_A(0, height/2, length);
          Vector<double>   disp_A(dim);

          // make a non-parallel copy of solution
          Vector<double> copy_solution(solution);

          Evaluation::PointValuesEvaluation<dim> point_values_evaluation(point_A);

          point_values_evaluation.compute (dof_handler, copy_solution, disp_A);

          table_results.add_value("time step", timestep_no);
          table_results.add_value("Cells", triangulation.n_global_active_cells());
          table_results.add_value("DoFs", dof_handler.n_dofs());
          table_results.add_value("pressure", pressure*present_time/end_time);
          table_results.add_value("u_A", disp_A(1));
        }

        {
          // demonstrate the location and maximum von-Mises stress in the
          // specified domain close to the clamped face, z = 0
          // top domain: height/2 - thickness_flange <= y <= height/2
          //             0 <= z <= 2 * thickness_flange
          // bottom domain: -height/2 <= y <= -height/2 + thickness_flange
          //             0 <= z <= 2 * thickness_flange

          double VM_stress_max (0);
          Point<dim> point_max;

          SymmetricTensor<2, dim> stress_at_qpoint;

          typename DoFHandler<dim>::active_cell_iterator
          cell = dof_handler.begin_active(),
          endc = dof_handler.end();

          const FEValuesExtractors::Vector displacement(0);

          for (; cell!=endc; ++cell)
            if (cell->is_locally_owned())
              {
                PointHistory<dim> *local_quadrature_points_history
                  = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
                Assert (local_quadrature_points_history >=
                        &quadrature_point_history.front(),
                        ExcInternalError());
                Assert (local_quadrature_points_history <
                        &quadrature_point_history.back(),
                        ExcInternalError());

                // Then loop over the quadrature points of this cell:
                for (unsigned int q=0; q<quadrature_formula.size(); ++q)
                  {
                    stress_at_qpoint = local_quadrature_points_history[q].old_stress;

                    const double VM_stress = Evaluation::get_von_Mises_stress(stress_at_qpoint);
                    if (VM_stress > VM_stress_max)
                      {
                        VM_stress_max = VM_stress;
                        point_max = local_quadrature_points_history[q].point;
                      }

                  }
              }

          table_results.add_value("maximum von_Mises stress", VM_stress_max);
          table_results.add_value("x", point_max[0]);
          table_results.add_value("y", point_max[1]);
          table_results.add_value("z", point_max[2]);

        }

      }


  }


  // @sect4{PlasticityContactProblem::run}

  // As in all other tutorial programs, the <code>run()</code> function contains
  // the overall logic. There is not very much to it here: in essence, it
  // performs the loops over all mesh refinement cycles, and within each, hands
  // things over to the Newton solver in <code>solve_newton()</code> on the
  // current mesh and calls the function that creates graphical output for
  // the so-computed solution. It then outputs some statistics concerning both
  // run times and memory consumption that has been collected over the course of
  // computations on this mesh.
  template <int dim>
  void
  ElastoPlasticProblem<dim>::run ()
  {
    computing_timer.reset();

    present_time = 0;
    present_timestep = 1;
    end_time = 10;
    timestep_no = 0;

    make_grid();

    // ----------------------------------------------------------------
    //    base_mesh == "Thick_tube_internal_pressure"
    /*
    const Point<dim> center(0, 0);
    const double inner_radius = .1,
        outer_radius = .2;

    const SphericalManifold<dim> inner_boundary_description(center, inner_radius);
    triangulation.set_manifold (0, inner_boundary_description);

    const SphericalManifold<dim> outer_boundary_description(center, outer_radius);
    triangulation.set_manifold (1, outer_boundary_description);
    */
    // ----------------------------------------------------------------
    //    base_mesh == "Perforated_strip_tension"
    /*
    const double inner_radius = 0.05;

    const CylinderBoundary<dim> inner_boundary_description(inner_radius, 2);
    triangulation.set_manifold (10, inner_boundary_description);
    */
    // ----------------------------------------------------------------

    setup_quadrature_point_history ();

    while (present_time < end_time)
      {
        present_time += present_timestep;
        ++timestep_no;

        if (present_time > end_time)
          {
            present_timestep -= (present_time - end_time);
            present_time = end_time;
          }
        pcout << std::endl;
        pcout << "Time step " << timestep_no << " at time " << present_time
              << std::endl;

        relative_error = max_relative_error * 10;
        current_refinement_cycle = 0;

        setup_system();


        // ------------------------ Refinement based on the relative error -------------------------------

        while (relative_error >= max_relative_error)
          {
            solve_newton();
            compute_error();

            if ( (timestep_no > 1) && (current_refinement_cycle>0) && (relative_error >= max_relative_error) )
              {
                pcout << "The relative error, " << relative_error
                      << " , is still more than maximum relative error, "
                      << max_relative_error << ", but we move to the next increment.\n";
                relative_error = .1 * max_relative_error;
              }

            if (relative_error >= max_relative_error)
              {
                TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
                ++current_refinement_cycle;
                refine_grid();
              }

          }

        // ------------------------ Refinement based on the number of refinement --------------------------
        /*
        bool continue_loop = true;
        while (continue_loop)
        {
          solve_newton();
          compute_error();

          if ( (timestep_no == 1) && (current_refinement_cycle < 1) )
          {
            TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
            ++current_refinement_cycle;
            refine_grid();
          }else
          {
            continue_loop = false;
          }

        }
        */

        // -------------------------------------------------------------------------------------------------

        solution += incremental_displacement;

        update_quadrature_point_history ();

        output_results((std::string("solution-") +
                        Utilities::int_to_string(timestep_no, 4)).c_str());

        computing_timer.print_summary();
        computing_timer.reset();

        Utilities::System::MemoryStats stats;
        Utilities::System::get_memory_stats(stats);
        pcout << "Peak virtual memory used, resident in kB: " << stats.VmSize << " "
              << stats.VmRSS << std::endl;


        if (std::abs(present_time-end_time) < 1.e-7)
          {
            const std::string filename = (output_dir + "Results");

            std::ofstream output_txt((filename + ".txt").c_str());

            pcout << std::endl;
            table_results.write_text(output_txt);
            pcout << std::endl;
            table_results_2.write_text(output_txt);
            pcout << std::endl;
            table_results_3.write_text(output_txt);
            pcout << std::endl;
          }

      }

    if (base_mesh == "Thick_tube_internal_pressure")
      {
        triangulation.reset_manifold (0);
        triangulation.reset_manifold (1);
      }
    else if (base_mesh == "Perforated_strip_tension")
      {
        triangulation.reset_manifold (10);
      }

  }
}

// @sect3{The <code>main</code> function}

// There really isn't much to the <code>main()</code> function. It looks
// like they always do:
int main (int argc, char *argv[])
{
  using namespace dealii;
  using namespace ElastoPlastic;

  try
    {
      deallog.depth_console(0);
      ParameterHandler prm;
      const int dim = 3;
      ElastoPlasticProblem<dim>::declare_parameters(prm);
      if (argc != 2)
        {
          std::cerr << "*** Call this program as <./elastoplastic input.prm>" << std::endl;
          return 1;
        }

      prm.parse_input(argv[1]);
      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
      {
        ElastoPlasticProblem<dim> problem(prm);
        problem.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;
}