step-42.h

Go to the documentation of this file.

 *     {
 *       TimerOutput::Scope t(computing_timer, "Setup: distribute DoFs");
 *       dof_handler.distribute_dofs(fe);
 * 
 *       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 and the active set */
 *     {
 *       TimerOutput::Scope t(computing_timer, "Setup: vectors");
 *       solution.reinit(locally_relevant_dofs, mpi_communicator);
 *       newton_rhs.reinit(locally_owned_dofs, mpi_communicator);
 *       newton_rhs_uncondensed.reinit(locally_owned_dofs, mpi_communicator);
 *       diag_mass_matrix_vector.reinit(locally_owned_dofs, mpi_communicator);
 *       fraction_of_plastic_q_points_per_cell.reinit(
 *         triangulation.n_active_cells());
 * 
 *       active_set.clear();
 *       active_set.set_size(dof_handler.n_dofs());
 *     }
 * 
 * @endcode
 * 
 * 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.
 * 
 * @code
 *     {
 *       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,
 *                                       Utilities::MPI::this_mpi_process(
 *                                         mpi_communicator));
 *       sp.compress();
 *       newton_matrix.reinit(sp);
 * 
 * 
 *       TrilinosWrappers::SparseMatrix &mass_matrix = newton_matrix;
 * 
 *       assemble_mass_matrix_diagonal(mass_matrix);
 * 
 *       const unsigned int start = (newton_rhs.local_range().first),
 *                          end   = (newton_rhs.local_range().second);
 *       for (unsigned int j = start; j < end; ++j)
 *         diag_mass_matrix_vector(j) = mass_matrix.diag_element(j);
 *       diag_mass_matrix_vector.compress(VectorOperation::insert);
 * 
 *       mass_matrix = 0;
 *     }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemcompute_dirichlet_constraints"></a> 
 * <h4>PlasticityContactProblem::compute_dirichlet_constraints</h4>
 * 

 * 
 * 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.
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::compute_dirichlet_constraints()
 *   {
 *     constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs);
 *     constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);
 * 
 *     if (base_mesh == "box")
 *       {
 * @endcode
 * 
 * interpolate all components of the solution
 * 
 * @code
 *         VectorTools::interpolate_boundary_values(
 *           dof_handler,
 *           6,
 *           EquationData::BoundaryValues<dim>(),
 *           constraints_dirichlet_and_hanging_nodes,
 *           ComponentMask());
 * 
 * @endcode
 * 
 * interpolate x- and y-components of the
 * solution (this is a bit mask, so apply
 * operator| )
 * 
 * @code
 *         const FEValuesExtractors::Scalar x_displacement(0);
 *         const FEValuesExtractors::Scalar y_displacement(1);
 *         VectorTools::interpolate_boundary_values(
 *           dof_handler,
 *           8,
 *           EquationData::BoundaryValues<dim>(),
 *           constraints_dirichlet_and_hanging_nodes,
 *           (fe.component_mask(x_displacement) |
 *            fe.component_mask(y_displacement)));
 *       }
 *     else
 *       VectorTools::interpolate_boundary_values(
 *         dof_handler,
 *         0,
 *         EquationData::BoundaryValues<dim>(),
 *         constraints_dirichlet_and_hanging_nodes,
 *         ComponentMask());
 * 
 *     constraints_dirichlet_and_hanging_nodes.close();
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemassemble_mass_matrix_diagonal"></a> 
 * <h4>PlasticityContactProblem::assemble_mass_matrix_diagonal</h4>
 * 

 * 
 * The next helper function computes the (diagonal) mass matrix that
 * is used to determine the active set of the active set method we use in
 * the contact algorithm. This matrix is of mass matrix type, but unlike
 * the standard mass matrix, we can make it diagonal (even in the case of
 * higher order elements) by using a quadrature formula that has its
 * quadrature points at exactly the same locations as the interpolation points
 * for the finite element are located. We achieve this by using a
 * QGaussLobatto quadrature formula here, along with initializing the finite
 * element with a set of interpolation points derived from the same quadrature
 * formula. The remainder of the function is relatively straightforward: we
 * put the resulting matrix into the given argument; because we know the
 * matrix is diagonal, it is sufficient to have a loop over only @f$i@f$ and
 * not over @f$j@f$. Strictly speaking, we could even avoid multiplying the
 * shape function's values at quadrature point <code>q_point</code> by itself
 * because we know the shape value to be a vector with exactly one one which
 * when dotted with itself yields one. Since this function is not time
 * critical we add this term for clarity.
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal(
 *     TrilinosWrappers::SparseMatrix &mass_matrix)
 *   {
 *     QGaussLobatto<dim - 1> face_quadrature_formula(fe.degree + 1);
 * 
 *     FEFaceValues<dim> fe_values_face(fe,
 *                                      face_quadrature_formula,
 *                                      update_values | update_JxW_values);
 * 
 *     const unsigned int dofs_per_cell   = fe.dofs_per_cell;
 *     const unsigned int n_face_q_points = face_quadrature_formula.size();
 * 
 *     FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 * 
 *     const FEValuesExtractors::Vector displacement(0);
 * 
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       if (cell->is_locally_owned())
 *         for (const auto &face : cell->face_iterators())
 *           if (face->at_boundary() && face->boundary_id() == 1)
 *             {
 *               fe_values_face.reinit(cell, face);
 *               cell_matrix = 0;
 * 
 *               for (unsigned int q_point = 0; q_point < n_face_q_points;
 *                    ++q_point)
 *                 for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                   cell_matrix(i, i) +=
 *                     (fe_values_face[displacement].value(i, q_point) *
 *                      fe_values_face[displacement].value(i, q_point) *
 *                      fe_values_face.JxW(q_point));
 * 
 *               cell->get_dof_indices(local_dof_indices);
 * 
 *               for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                 mass_matrix.add(local_dof_indices[i],
 *                                 local_dof_indices[i],
 *                                 cell_matrix(i, i));
 *             }
 *     mass_matrix.compress(VectorOperation::add);
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemupdate_solution_and_constraints"></a> 
 * <h4>PlasticityContactProblem::update_solution_and_constraints</h4>
 * 

 * 
 * The following function is the first function we call in each Newton
 * iteration in the <code>solve_newton()</code> function. What it does is
 * to project the solution onto the feasible set and update the active set
 * for the degrees of freedom that touch or penetrate the obstacle.
 *   

 * 
 * In order to function, we first need to do some bookkeeping: We need
 * to write into the solution vector (which we can only do with fully
 * distributed vectors without ghost elements) and we need to read
 * the Lagrange multiplier and the elements of the diagonal mass matrix
 * from their respective vectors (which we can only do with vectors that
 * do have ghost elements), so we create the respective vectors. We then
 * also initialize the constraints object that will contain constraints
 * from contact and all other sources, as well as an object that contains
 * an index set of all locally owned degrees of freedom that are part of
 * the contact:
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::update_solution_and_constraints()
 *   {
 *     std::vector<bool> dof_touched(dof_handler.n_dofs(), false);
 * 
 *     TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
 *                                                        mpi_communicator);
 *     distributed_solution = solution;
 * 
 *     TrilinosWrappers::MPI::Vector lambda(locally_relevant_dofs,
 *                                          mpi_communicator);
 *     lambda = newton_rhs_uncondensed;
 * 
 *     TrilinosWrappers::MPI::Vector diag_mass_matrix_vector_relevant(
 *       locally_relevant_dofs, mpi_communicator);
 *     diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
 * 
 * 
 *     all_constraints.reinit(locally_relevant_dofs);
 *     active_set.clear();
 * 
 * @endcode
 * 
 * The second part is a loop over all cells in which we look at each
 * point where a degree of freedom is defined whether the active set
 * condition is true and we need to add this degree of freedom to
 * the active set of contact nodes. As we always do, if we want to
 * evaluate functions at individual points, we do this with an
 * FEValues object (or, here, an FEFaceValues object since we need to
 * check contact at the surface) with an appropriately chosen quadrature
 * object. We create this face quadrature object by choosing the
 * "support points" of the shape functions defined on the faces
 * of cells (for more on support points, see this
 * @ref GlossSupport "glossary entry"). As a consequence, we have as
 * many quadrature points as there are shape functions per face and
 * looping over quadrature points is equivalent to looping over shape
 * functions defined on a face. With this, the code looks as follows:
 * 
 * @code
 *     Quadrature<dim - 1> face_quadrature(fe.get_unit_face_support_points());
 *     FEFaceValues<dim>   fe_values_face(fe,
 *                                      face_quadrature,
 *                                      update_quadrature_points);
 * 
 *     const unsigned int dofs_per_face   = fe.dofs_per_face;
 *     const unsigned int n_face_q_points = face_quadrature.size();
 * 
 *     std::vector<types::global_dof_index> dof_indices(dofs_per_face);
 * 
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       if (!cell->is_artificial())
 *         for (const auto &face : cell->face_iterators())
 *           if (face->at_boundary() && face->boundary_id() == 1)
 *             {
 *               fe_values_face.reinit(cell, face);
 *               face->get_dof_indices(dof_indices);
 * 
 *               for (unsigned int q_point = 0; q_point < n_face_q_points;
 *                    ++q_point)
 *                 {
 * @endcode
 * 
 * At each quadrature point (i.e., at each support point of a
 * degree of freedom located on the contact boundary), we then
 * ask whether it is part of the z-displacement degrees of
 * freedom and if we haven't encountered this degree of
 * freedom yet (which can happen for those on the edges
 * between faces), we need to evaluate the gap between the
 * deformed object and the obstacle. If the active set
 * condition is true, then we add a constraint to the
 * AffineConstraints object that the next Newton update needs
 * to satisfy, set the solution vector's corresponding element
 * to the correct value, and add the index to the IndexSet
 * object that stores which degree of freedom is part of the
 * contact:
 * 
 * @code
 *                   const unsigned int component =
 *                     fe.face_system_to_component_index(q_point).first;
 * 
 *                   const unsigned int index_z = dof_indices[q_point];
 * 
 *                   if ((component == 2) && (dof_touched[index_z] == false))
 *                     {
 *                       dof_touched[index_z] = true;
 * 
 *                       const Point<dim> this_support_point =
 *                         fe_values_face.quadrature_point(q_point);
 * 
 *                       const double obstacle_value =
 *                         obstacle->value(this_support_point, 2);
 *                       const double solution_here = solution(index_z);
 *                       const double undeformed_gap =
 *                         obstacle_value - this_support_point(2);
 * 
 *                       const double c = 100.0 * e_modulus;
 *                       if ((lambda(index_z) /
 *                                diag_mass_matrix_vector_relevant(index_z) +
 *                              c * (solution_here - undeformed_gap) >
 *                            0) &&
 *                           !constraints_hanging_nodes.is_constrained(index_z))
 *                         {
 *                           all_constraints.add_line(index_z);
 *                           all_constraints.set_inhomogeneity(index_z,
 *                                                             undeformed_gap);
 *                           distributed_solution(index_z) = undeformed_gap;
 * 
 *                           active_set.add_index(index_z);
 *                         }
 *                     }
 *                 }
 *             }
 * 
 * @endcode
 * 
 * At the end of this function, we exchange data between processors updating
 * those ghost elements in the <code>solution</code> variable that have been
 * written by other processors. We then merge the Dirichlet constraints and
 * those from hanging nodes into the AffineConstraints object that already
 * contains the active set. We finish the function by outputting the total
 * number of actively constrained degrees of freedom for which we sum over
 * the number of actively constrained degrees of freedom owned by each
 * of the processors. This number of locally owned constrained degrees of
 * freedom is of course the number of elements of the intersection of the
 * active set and the set of locally owned degrees of freedom, which
 * we can get by using <code>operator&</code> on two IndexSets:
 * 
 * @code
 *     distributed_solution.compress(VectorOperation::insert);
 *     solution = distributed_solution;
 * 
 *     all_constraints.close();
 *     all_constraints.merge(constraints_dirichlet_and_hanging_nodes);
 * 
 *     pcout << "         Size of active set: "
 *           << Utilities::MPI::sum((active_set & locally_owned_dofs).n_elements(),
 *                                  mpi_communicator)
 *           << std::endl;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemassemble_newton_system"></a> 
 * <h4>PlasticityContactProblem::assemble_newton_system</h4>
 * 

 * 
 * 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::distribute_local_to_global(), using the
 * constraints we have previously computed.
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::assemble_newton_system(
 *     const TrilinosWrappers::MPI::Vector &linearization_point)
 *   {
 *     TimerOutput::Scope t(computing_timer, "Assembling");
 * 
 *     QGauss<dim>     quadrature_formula(fe.degree + 1);
 *     QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
 * 
 *     FEValues<dim> fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients |
 *                               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::BoundaryForce<dim> boundary_force;
 *     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);
 * 
 *     const FEValuesExtractors::Vector displacement(0);
 * 
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       if (cell->is_locally_owned())
 *         {
 *           fe_values.reinit(cell);
 *           cell_matrix = 0;
 *           cell_rhs    = 0;
 * 
 *           std::vector<SymmetricTensor<2, dim>> strain_tensor(n_q_points);
 *           fe_values[displacement].get_function_symmetric_gradients(
 *             linearization_point, 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; ++i)
 *                 {
 * @endcode
 * 
 * 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
 * @f$\varphi_i@f$, i.e. the term @f$I_\Pi\varepsilon(\varphi_i)@f$,
 * so we introduce an abbreviation of this term. Recall that
 * the matrix corresponds to the bilinear form
 * @f$A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))@f$
 * in the notation of the accompanying publication, whereas
 * the right hand side is @f$F_i=([I_\Pi-P_\Pi
 * C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))@f$ where @f$u@f$
 * 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
 * @f$I_\Pi=P_\Pi@f$) 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.
 * 
 * @code
 *                   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 -
 *                       stress_strain_tensor *
 *                         fe_values[displacement].symmetric_gradient(i,
 *                                                                    q_point)) *
 *                      strain_tensor[q_point] * fe_values.JxW(q_point));
 *                 }
 *             }
 * 
 *           for (const auto &face : cell->face_iterators())
 *             if (face->at_boundary() && face->boundary_id() == 1)
 *               {
 *                 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;
 *                     rhs_values[2] = boundary_force_values[q_point][2];
 *                     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);
 *           all_constraints.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);
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemcompute_nonlinear_residual"></a> 
 * <h4>PlasticityContactProblem::compute_nonlinear_residual</h4>
 * 

 * 
 * 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:
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::compute_nonlinear_residual(
 *     const TrilinosWrappers::MPI::Vector &linearization_point)
 *   {
 *     QGauss<dim>     quadrature_formula(fe.degree + 1);
 *     QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
 * 
 *     FEValues<dim> fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients |
 *                               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::BoundaryForce<dim> boundary_force;
 *     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             = 0;
 *     newton_rhs_uncondensed = 0;
 * 
 *     fraction_of_plastic_q_points_per_cell = 0;
 * 
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       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);
 * 
 *           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->active_cell_index());
 * 
 *               for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                 {
 *                   cell_rhs(i) -=
 *                     (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 (const auto &face : cell->face_iterators())
 *             if (face->at_boundary() && face->boundary_id() == 1)
 *               {
 *                 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;
 *                     rhs_values[2] = boundary_force_values[q_point][2];
 *                     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);
 * 
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             newton_rhs_uncondensed(local_dof_indices[i]) += cell_rhs(i);
 *         }
 * 
 *     fraction_of_plastic_q_points_per_cell /= quadrature_formula.size();
 *     newton_rhs.compress(VectorOperation::add);
 *     newton_rhs_uncondensed.compress(VectorOperation::add);
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemsolve_newton_system"></a> 
 * <h4>PlasticityContactProblem::solve_newton_system</h4>
 * 

 * 
 * 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::set_zero function to newton_rhs.
 * This is necessary if a hanging node with solution value @f$x_0@f$
 * has one neighbor with value @f$x_1@f$ which is in contact with the
 * obstacle and one neighbor @f$x_2@f$ 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 @f$x_0 = x_1/2 +
 * \text{gap}/2@f$. So the corresponding entries in the right-hand-side are
 * non-zero with a meaningless value. These values we have to set to zero.
 * - Like in @ref step_40 "step-40", we need to shuffle between vectors that do and do
 * not have ghost elements when solving or using the solution.
 *   

 * 
 * The rest of the function is similar to @ref step_40 "step-40" and
 * @ref step_41 "step-41" except that we use a BiCGStab solver
 * instead of CG. This is due to the fact that for very small hardening
 * parameters @f$\gamma@f$, the linear system becomes almost semidefinite though
 * still symmetric. BiCGStab appears to have an easier time with such linear
 * systems.
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::solve_newton_system()
 *   {
 *     TimerOutput::Scope t(computing_timer, "Solve");
 * 
 *     TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
 *                                                        mpi_communicator);
 *     distributed_solution = solution;
 * 
 *     constraints_hanging_nodes.set_zero(distributed_solution);
 *     constraints_hanging_nodes.set_zero(newton_rhs);
 * 
 *     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 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;
 *     }
 * 
 *     all_constraints.distribute(distributed_solution);
 * 
 *     solution = distributed_solution;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemsolve_newton"></a> 
 * <h4>PlasticityContactProblem::solve_newton</h4>
 * 

 * 
 * 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 the 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. We initialize the variable
 * <code>previous_residual_norm</code> to the most negative value
 * representable with double precision numbers so that the
 * comparison whether the current residual is less than that of the
 * previous step will always fail in the first step.
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<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);
 * 
 *     double residual_norm;
 *     double previous_residual_norm = -std::numeric_limits<double>::max();
 * 
 *     const double correct_sigma = sigma_0;
 * 
 *     IndexSet old_active_set(active_set);
 * 
 *     for (unsigned int newton_step = 1; newton_step <= 100; ++newton_step)
 *       {
 *         if (newton_step == 1 &&
 *             ((transfer_solution && current_refinement_cycle == 0) ||
 *              !transfer_solution))
 *           constitutive_law.set_sigma_0(1e+10);
 *         else if (newton_step == 2 || current_refinement_cycle > 0 ||
 *                  !transfer_solution)
 *           constitutive_law.set_sigma_0(correct_sigma);
 * 
 *         pcout << " " << std::endl;
 *         pcout << "   Newton iteration " << newton_step << std::endl;
 *         pcout << "      Updating active set..." << std::endl;
 * 
 *         {
 *           TimerOutput::Scope t(computing_timer, "update active set");
 *           update_solution_and_constraints();
 *         }
 * 
 *         pcout << "      Assembling system... " << std::endl;
 *         newton_matrix = 0;
 *         newton_rhs    = 0;
 *         assemble_newton_system(solution);
 * 
 *         pcout << "      Solving system... " << std::endl;
 *         solve_newton_system();
 * 
 * @endcode
 * 
 * It gets a bit more hairy after we have computed the
 * trial solution @f$\tilde{\mathbf u}@f$ 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.
 * 
 * @code
 *         if ((newton_step == 1) ||
 *             (transfer_solution && newton_step == 2 &&
 *              current_refinement_cycle == 0) ||
 *             (!transfer_solution && newton_step == 2))
 *           {
 *             compute_nonlinear_residual(solution);
 *             old_solution = solution;
 * 
 *             residual                     = newton_rhs;
 *             const unsigned int start_res = (residual.local_range().first),
 *                                end_res   = (residual.local_range().second);
 *             for (unsigned int n = start_res; n < end_res; ++n)
 *               if (all_constraints.is_inhomogeneously_constrained(n))
 *                 residual(n) = 0;
 * 
 *             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 = solution;
 * 
 *                 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;
 *                 compute_nonlinear_residual(locally_relevant_tmp_vector);
 *                 residual = newton_rhs;
 * 
 *                 const unsigned int start_res = (residual.local_range().first),
 *                                    end_res   = (residual.local_range().second);
 *                 for (unsigned int n = start_res; n < end_res; ++n)
 *                   if (all_constraints.is_inhomogeneously_constrained(n))
 *                     residual(n) = 0;
 * 
 *                 residual.compress(VectorOperation::insert);
 * 
 *                 residual_norm = residual.l2_norm();
 * 
 *                 pcout
 *                   << "      Residual of the non-contact part of the system: "
 *                   << residual_norm << std::endl
 *                   << "         with a damping parameter alpha = " << alpha
 *                   << std::endl;
 * 
 *                 if (residual_norm < previous_residual_norm)
 *                   break;
 *               }
 * 
 *             solution     = tmp_vector;
 *             old_solution = solution;
 *           }
 * 
 *         previous_residual_norm = residual_norm;
 * 
 * 
 * @endcode
 * 
 * The final step is to check for convergence. If the active set
 * has not changed across all processors and the residual is
 * less than a threshold of @f$10^{-10}@f$, then we terminate
 * the iteration on the current mesh:
 * 
 * @code
 *         if (Utilities::MPI::sum((active_set == old_active_set) ? 0 : 1,
 *                                 mpi_communicator) == 0)
 *           {
 *             pcout << "      Active set did not change!" << std::endl;
 *             if (residual_norm < 1e-10)
 *               break;
 *           }
 * 
 *         old_active_set = active_set;
 *       }
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemrefine_grid"></a> 
 * <h4>PlasticityContactProblem::refine_grid</h4>
 * 

 * 
 * 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>.
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::refine_grid()
 *   {
 *     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
 *       {
 *         Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 *         KellyErrorEstimator<dim>::estimate(
 *           dof_handler,
 *           QGauss<dim - 1>(fe.degree + 2),
 *           std::map<types::boundary_id, const Function<dim> *>(),
 *           solution,
 *           estimated_error_per_cell);
 * 
 *         parallel::distributed::GridRefinement ::refine_and_coarsen_fixed_number(
 *           triangulation, estimated_error_per_cell, 0.3, 0.03);
 *       }
 * 
 *     triangulation.prepare_coarsening_and_refinement();
 * 
 *     parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::Vector>
 *       solution_transfer(dof_handler);
 *     if (transfer_solution)
 *       solution_transfer.prepare_for_coarsening_and_refinement(solution);
 * 
 *     triangulation.execute_coarsening_and_refinement();
 * 
 *     setup_system();
 * 
 *     if (transfer_solution)
 *       {
 *         TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
 *                                                            mpi_communicator);
 *         solution_transfer.interpolate(distributed_solution);
 * 
 * @endcode
 * 
 * enforce constraints to make the interpolated solution conforming on
 * the new mesh:
 * 
 * @code
 *         constraints_hanging_nodes.distribute(distributed_solution);
 * 
 *         solution = distributed_solution;
 *         compute_nonlinear_residual(solution);
 *       }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemmove_mesh"></a> 
 * <h4>PlasticityContactProblem::move_mesh</h4>
 * 

 * 
 * 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.
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::move_mesh(
 *     const TrilinosWrappers::MPI::Vector &displacement) const
 *   {
 *     std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
 * 
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       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;
 *             }
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemoutput_results"></a> 
 * <h4>PlasticityContactProblem::output_results</h4>
 * 

 * 
 * 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.
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::output_results(
 *     const unsigned int current_refinement_cycle)
 *   {
 *     TimerOutput::Scope t(computing_timer, "Graphical output");
 * 
 *     pcout << "      Writing graphical output... " << std::flush;
 * 
 *     move_mesh(solution);
 * 
 * @endcode
 * 
 * Calculation of the contact forces
 * 
 * @code
 *     TrilinosWrappers::MPI::Vector distributed_lambda(locally_owned_dofs,
 *                                                      mpi_communicator);
 *     const unsigned int start_res = (newton_rhs_uncondensed.local_range().first),
 *                        end_res = (newton_rhs_uncondensed.local_range().second);
 *     for (unsigned int n = start_res; n < end_res; ++n)
 *       if (all_constraints.is_inhomogeneously_constrained(n))
 *         distributed_lambda(n) =
 *           newton_rhs_uncondensed(n) / diag_mass_matrix_vector(n);
 *     distributed_lambda.compress(VectorOperation::insert);
 *     constraints_hanging_nodes.distribute(distributed_lambda);
 * 
 *     TrilinosWrappers::MPI::Vector lambda(locally_relevant_dofs,
 *                                          mpi_communicator);
 *     lambda = distributed_lambda;
 * 
 *     TrilinosWrappers::MPI::Vector distributed_active_set_vector(
 *       locally_owned_dofs, mpi_communicator);
 *     distributed_active_set_vector = 0.;
 *     for (const auto index : active_set)
 *       distributed_active_set_vector[index] = 1.;
 *     distributed_lambda.compress(VectorOperation::insert);
 * 
 *     TrilinosWrappers::MPI::Vector active_set_vector(locally_relevant_dofs,
 *                                                     mpi_communicator);
 *     active_set_vector = distributed_active_set_vector;
 * 
 *     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);
 *     data_out.add_data_vector(lambda,
 *                              std::vector<std::string>(dim, "contact_force"),
 *                              DataOut<dim>::type_dof_data,
 *                              data_component_interpretation);
 *     data_out.add_data_vector(active_set_vector,
 *                              std::vector<std::string>(dim, "active_set"),
 *                              DataOut<dim>::type_dof_data,
 *                              data_component_interpretation);
 * 
 *     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();
 * 
 * @endcode
 * 
 * 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.
 * 
 * @code
 *     const std::string master_name = data_out.write_vtu_with_pvtu_record(
 *       output_dir, "solution", current_refinement_cycle, mpi_communicator, 2);
 *     pcout << master_name << std::endl;
 * 
 *     TrilinosWrappers::MPI::Vector tmp(solution);
 *     tmp *= -1;
 *     move_mesh(tmp);
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemoutput_contact_force"></a> 
 * <h4>PlasticityContactProblem::output_contact_force</h4>
 * 

 * 
 * This last auxiliary function computes the contact force by
 * calculating an integral over the contact pressure in z-direction
 * over the contact area. For this purpose we set the contact
 * pressure lambda to 0 for all inactive dofs (whether a degree
 * of freedom is part of the contact is determined just as
 * we did in the previous function). For all
 * active dofs, lambda contains the quotient of the nonlinear
 * residual (newton_rhs_uncondensed) and corresponding diagonal entry
 * of the mass matrix (diag_mass_matrix_vector). Because it is
 * not unlikely that hanging nodes show up in the contact area
 * it is important to apply constraints_hanging_nodes.distribute
 * to the distributed_lambda vector.
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::output_contact_force() const
 *   {
 *     TrilinosWrappers::MPI::Vector distributed_lambda(locally_owned_dofs,
 *                                                      mpi_communicator);
 *     const unsigned int start_res = (newton_rhs_uncondensed.local_range().first),
 *                        end_res = (newton_rhs_uncondensed.local_range().second);
 *     for (unsigned int n = start_res; n < end_res; ++n)
 *       if (all_constraints.is_inhomogeneously_constrained(n))
 *         distributed_lambda(n) =
 *           newton_rhs_uncondensed(n) / diag_mass_matrix_vector(n);
 *       else
 *         distributed_lambda(n) = 0;
 *     distributed_lambda.compress(VectorOperation::insert);
 *     constraints_hanging_nodes.distribute(distributed_lambda);
 * 
 *     TrilinosWrappers::MPI::Vector lambda(locally_relevant_dofs,
 *                                          mpi_communicator);
 *     lambda = distributed_lambda;
 * 
 *     double contact_force = 0.0;
 * 
 *     QGauss<dim - 1>   face_quadrature_formula(fe.degree + 1);
 *     FEFaceValues<dim> fe_values_face(fe,
 *                                      face_quadrature_formula,
 *                                      update_values | update_JxW_values);
 * 
 *     const unsigned int n_face_q_points = face_quadrature_formula.size();
 * 
 *     const FEValuesExtractors::Vector displacement(0);
 * 
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       if (cell->is_locally_owned())
 *         for (const auto &face : cell->face_iterators())
 *           if (face->at_boundary() && face->boundary_id() == 1)
 *             {
 *               fe_values_face.reinit(cell, face);
 * 
 *               std::vector<Tensor<1, dim>> lambda_values(n_face_q_points);
 *               fe_values_face[displacement].get_function_values(lambda,
 *                                                                lambda_values);
 * 
 *               for (unsigned int q_point = 0; q_point < n_face_q_points;
 *                    ++q_point)
 *                 contact_force +=
 *                   lambda_values[q_point][2] * fe_values_face.JxW(q_point);
 *             }
 *     contact_force = Utilities::MPI::sum(contact_force, MPI_COMM_WORLD);
 * 
 *     pcout << "Contact force = " << contact_force << std::endl;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="PlasticityContactProblemrun"></a> 
 * <h4>PlasticityContactProblem::run</h4>
 * 

 * 
 * 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.
 * 
 * @code
 *   template <int dim>
 *   void PlasticityContactProblem<dim>::run()
 *   {
 *     computing_timer.reset();
 *     for (; current_refinement_cycle < n_refinement_cycles;
 *          ++current_refinement_cycle)
 *       {
 *         {
 *           TimerOutput::Scope t(computing_timer, "Setup");
 * 
 *           pcout << std::endl;
 *           pcout << "Cycle " << current_refinement_cycle << ':' << std::endl;
 * 
 *           if (current_refinement_cycle == 0)
 *             {
 *               make_grid();
 *               setup_system();
 *             }
 *           else
 *             {
 *               TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
 *               refine_grid();
 *             }
 *         }
 * 
 *         solve_newton();
 * 
 *         output_results(current_refinement_cycle);
 * 
 *         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 (base_mesh == "box")
 *           output_contact_force();
 *       }
 *   }
 * } // namespace Step42
 * 
 * @endcode
 * 
 * 
 * <a name="Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 

 * 
 * There really isn't much to the <code>main()</code> function. It looks
 * like they always do:
 * 
 * @code
 * int main(int argc, char *argv[])
 * {
 *   using namespace dealii;
 *   using namespace Step42;
 * 
 *   try
 *     {
 *       ParameterHandler prm;
 *       PlasticityContactProblem<3>::declare_parameters(prm);
 *       if (argc != 2)
 *         {
 *           std::cerr << "*** Call this program as <./step-42 input.prm>"
 *                     << std::endl;
 *           return 1;
 *         }
 * 
 *       prm.parse_input(argv[1]);
 *       Utilities::MPI::MPI_InitFinalize mpi_initialization(
 *         argc, argv, numbers::invalid_unsigned_int);
 *       {
 *         PlasticityContactProblem<3> 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;
 * }
 * @endcode
<a name="Results"></a><h1>Results</h1>


The directory that contains this program also contains a number of input
parameter files that can be used to create various different
simulations. For example, running the program with the
<code>p1_adaptive.prm</code> parameter file (using a ball as obstacle and the
box as domain) on 16 cores produces output like this:
@code
    Using output directory 'p1adaptive/'
    FE degree 1
    transfer solution false

Cycle 0:
   Number of active cells: 512
   Number of degrees of freedom: 2187

  Newton iteration 1
      Updating active set...
         Size of active set: 1
      Assembling system...
      Solving system...
         Error: 173.076 -> 1.64265e-06 in 7 Bicgstab iterations.
      Accepting Newton solution with residual: 1.64265e-06

   Newton iteration 2
      Updating active set...
         Size of active set: 1
      Assembling system...
      Solving system...
         Error: 57.3622 -> 3.23721e-07 in 8 Bicgstab iterations.
      Accepting Newton solution with residual: 24.9028
      Active set did not change!

   Newton iteration 3
      Updating active set...
         Size of active set: 1
      Assembling system...
      Solving system...
         Error: 24.9028 -> 9.94326e-08 in 7 Bicgstab iterations.
      Residual of the non-contact part of the system: 1.63333
         with a damping parameter alpha = 1
      Active set did not change!

...

  Newton iteration 6
      Updating active set...
         Size of active set: 1
      Assembling system...
      Solving system...
         Error: 1.43188e-07 -> 3.56218e-16 in 8 Bicgstab iterations.
      Residual of the non-contact part of the system: 4.298e-14
         with a damping parameter alpha = 1
      Active set did not change!
      Writing graphical output... p1_adaptive/solution-00.pvtu


+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |      1.13s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assembling                      |         6 |     0.463s |        41% |
| Graphical output                |         1 |    0.0257s |       2.3% |
| Residual and lambda             |         4 |    0.0754s |       6.7% |
| Setup                           |         1 |     0.227s |        20% |
| Setup: constraints              |         1 |    0.0347s |       3.1% |
| Setup: distribute DoFs          |         1 |    0.0441s |       3.9% |
| Setup: matrix                   |         1 |    0.0119s |       1.1% |
| Setup: vectors                  |         1 |   0.00155s |      0.14% |
| Solve                           |         6 |     0.246s |        22% |
| Solve: iterate                  |         6 |    0.0631s |       5.6% |
| Solve: setup preconditioner     |         6 |     0.167s |        15% |
| update active set               |         6 |    0.0401s |       3.6% |
+---------------------------------+-----------+------------+------------+

Peak virtual memory used, resident in kB: 541884 77464
Contact force = 37.3058

...

Cycle 3:
   Number of active cells: 14652
   Number of degrees of freedom: 52497

   Newton iteration 1
      Updating active set...
         Size of active set: 145
      Assembling system...
      Solving system...
         Error: 296.309 -> 2.72484e-06 in 10 Bicgstab iterations.
      Accepting Newton solution with residual: 2.72484e-06

...

   Newton iteration 10
      Updating active set...
         Size of active set: 145
      Assembling system...
      Solving system...
         Error: 2.71541e-07 -> 1.5428e-15 in 27 Bicgstab iterations.
      Residual of the non-contact part of the system: 1.89261e-13
         with a damping parameter alpha = 1
      Active set did not change!
      Writing graphical output... p1_adaptive/solution-03.pvtu


+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |      38.4s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assembling                      |        10 |      22.5s |        58% |
| Graphical output                |         1 |     0.327s |      0.85% |
| Residual and lambda             |         9 |      3.75s |       9.8% |
| Setup                           |         1 |      4.83s |        13% |
| Setup: constraints              |         1 |     0.578s |       1.5% |
| Setup: distribute DoFs          |         1 |      0.71s |       1.8% |
| Setup: matrix                   |         1 |     0.111s |      0.29% |
| Setup: refine mesh              |         1 |      4.83s |        13% |
| Setup: vectors                  |         1 |   0.00548s |     0.014% |
| Solve                           |        10 |      5.49s |        14% |
| Solve: iterate                  |        10 |       3.5s |       9.1% |
| Solve: setup preconditioner     |        10 |      1.84s |       4.8% |
| update active set               |        10 |     0.662s |       1.7% |
+---------------------------------+-----------+------------+------------+

Peak virtual memory used, resident in kB: 566052 105788
Contact force = 56.794

...
@endcode

The tables at the end of each cycle show information about computing time
(these numbers are of course specific to the machine on which this output
was produced)
and the number of calls of different parts of the program like assembly or
calculating the residual, for the most recent mesh refinement cycle. Some of
the numbers above can be improved by transferring the solution from one mesh to
the next, an option we have not exercised here. Of course, you can also make
the program run faster, especially on the later refinement cycles, by just
using more processors: the accompanying paper shows good scaling to at least
1000 cores.

In a typical run, you can observe that for every refinement step, the active
set - the contact points - are iterated out at first. After that the Newton
method has only to resolve the plasticity. For the finer meshes,
quadratic convergence can be observed for the last 4 or 5 Newton iterations.

We will not discuss here in all detail what happens with each of the input
files. Rather, let us just show pictures of the solution (the left half of the
domain is omitted if cells have zero quadrature points at which the plastic
inequality is active):

<table align="center">
  <tr>
    <td>
    <img src="https://www.dealii.org/images/steps/developer/step-42.CellConstitutionColorbar.png">
    </td>
    <td>
    <img src="https://www.dealii.org/images/steps/developer/step-42.CellConstitutionBall2.png" alt="" width="70%">
    </td>
    <td valign="top">
      &nbsp;
    </td>
    <td>
    <img src="https://www.dealii.org/images/steps/developer/step-42.CellConstitutionLi2.png" alt="" alt="" width="70%">
    </td>
  </tr>
</table>

The picture shows the adaptive refinement and as well how much a cell is
plastified during the contact with the ball. Remember that we consider the
norm of the deviator part of the stress in each quadrature point to
see if there is elastic or plastic behavior.
The blue
color means that this cell contains only elastic quadrature points in
contrast to the red cells in which all quadrature points are plastified.
In the middle of the top surface -
where the mesh is finest - a very close look shows the dimple caused by the
obstacle. This is the result of the <code>move_mesh()</code>
function. However, because the indentation of the obstacles we consider here
is so small, it is hard to discern this effect; one could play with displacing
vertices of the mesh by a multiple of the computed displacement.

Further discussion of results that can be obtained using this program is
provided in the publication mentioned at the very top of this page.


<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h1>Possibilities for extensions</h1>


There are, as always, multiple possibilities for extending this program. From
an algorithmic perspective, this program goes about as far as one can at the
time of writing, using the best available algorithms for the contact
inequality, the plastic nonlinearity, and the linear solvers. However, there
are things one would like to do with this program as far as more realistic
situations are concerned:
<ul>
<li> Extend the program from a static to a quasi-static situation, perhaps by
choosing a backward-Euler-scheme for the time discretization. Some theoretical
results can be found in the PhD thesis by Jörg Frohne, <i>FEM-Simulation
der Umformtechnik metallischer Oberfl&auml;chen im Mikrokosmos</i>, University
of Siegen, Germany, 2011.

<li> It would also be an interesting advance to consider a contact problem
with friction. In almost every mechanical process friction has a big
influence.  To model this situation, we have to take into account tangential
stresses at the contact surface. Friction also adds another inequality to
our problem since body and obstacle will typically stick together as long as
the tangential stress does not exceed a certain limit, beyond which the two
bodies slide past each other.

<li> If we already simulate a frictional contact, the next step to consider
is heat generation over the contact zone. The heat that is
caused by friction between two bodies raises the temperature in the
deformable body and entails an change of some material parameters.

<li> It might be of interest to implement more accurate, problem-adapted error
estimators for contact as well as for the plasticity.
</ul>
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-42.cc"
*/