step-41.h

Go to the documentation of this file.

,
 *                            const unsigned int component = 0) const override
 *       {
 *         (void)component;
 *         AssertIndexRange(component, 1);
 *   
 *         return -10;
 *       }
 *     };
 *   
 *   
 *   
 *     template <int dim>
 *     class BoundaryValues : public Function<dim>
 *     {
 *     public:
 *       virtual double value(const Point<dim> & /*p*/,
 *                            const unsigned int component = 0) const override
 *       {
 *         (void)component;
 *         AssertIndexRange(component, 1);
 *   
 *         return 0;
 *       }
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * We describe the obstacle function by a cascaded barrier (think: stair
 * steps):
 * 
 * @code
 *     template <int dim>
 *     class Obstacle : public Function<dim>
 *     {
 *     public:
 *       virtual double value(const Point<dim>  &p,
 *                            const unsigned int component = 0) const override
 *       {
 *         (void)component;
 *         Assert(component == 0, ExcIndexRange(component, 0, 1));
 *   
 *         if (p[0] < -0.5)
 *           return -0.2;
 *         else if (p[0] >= -0.5 && p[0] < 0.0)
 *           return -0.4;
 *         else if (p[0] >= 0.0 && p[0] < 0.5)
 *           return -0.6;
 *         else
 *           return -0.8;
 *       }
 *     };
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_41-ImplementationofthecodeObstacleProblemcodeclass"></a> 
 * <h3>Implementation of the <code>ObstacleProblem</code> class</h3>
 * 
 
 * 
 * 
 
 * 
 * 
 * <a name="step_41-ObstacleProblemObstacleProblem"></a> 
 * <h4>ObstacleProblem::ObstacleProblem</h4>
 * 
 
 * 
 * To everyone who has taken a look at the first few tutorial programs, the
 * constructor is completely obvious:
 * 
 * @code
 *     template <int dim>
 *     ObstacleProblem<dim>::ObstacleProblem()
 *       : fe(1)
 *       , dof_handler(triangulation)
 *     {}
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_41-ObstacleProblemmake_grid"></a> 
 * <h4>ObstacleProblem::make_grid</h4>
 * 
 
 * 
 * We solve our obstacle problem on the square @f$[-1,1]\times [-1,1]@f$ in
 * 2d. This function therefore just sets up one of the simplest possible
 * meshes.
 * 
 * @code
 *     template <int dim>
 *     void ObstacleProblem<dim>::make_grid()
 *     {
 *       GridGenerator::hyper_cube(triangulation, -1, 1);
 *       triangulation.refine_global(7);
 *   
 *       std::cout << "Number of active cells: " << triangulation.n_active_cells()
 *                 << std::endl
 *                 << "Total number of cells: " << triangulation.n_cells()
 *                 << std::endl;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_41-ObstacleProblemsetup_system"></a> 
 * <h4>ObstacleProblem::setup_system</h4>
 * 
 
 * 
 * In this first function of note, we set up the degrees of freedom handler,
 * resize vectors and matrices, and deal with the constraints. Initially,
 * the constraints are, of course, only given by boundary values, so we
 * interpolate them towards the top of the function.
 * 
 * @code
 *     template <int dim>
 *     void ObstacleProblem<dim>::setup_system()
 *     {
 *       dof_handler.distribute_dofs(fe);
 *       active_set.set_size(dof_handler.n_dofs());
 *   
 *       std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
 *                 << std::endl
 *                 << std::endl;
 *   
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                0,
 *                                                BoundaryValues<dim>(),
 *                                                constraints);
 *       constraints.close();
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
 *   
 *       system_matrix.reinit(dsp);
 *       complete_system_matrix.reinit(dsp);
 *   
 *       IndexSet solution_index_set = dof_handler.locally_owned_dofs();
 *       solution.reinit(solution_index_set, MPI_COMM_WORLD);
 *       system_rhs.reinit(solution_index_set, MPI_COMM_WORLD);
 *       complete_system_rhs.reinit(solution_index_set, MPI_COMM_WORLD);
 *       contact_force.reinit(solution_index_set, MPI_COMM_WORLD);
 *   
 * @endcode
 * 
 * The only other thing to do here is to compute the factors in the @f$B@f$
 * matrix which is used to scale the residual. As discussed in the
 * introduction, we'll use a little trick to make this mass matrix
 * diagonal, and in the following then first compute all of this as a
 * matrix and then extract the diagonal elements for later use:
 * 
 * @code
 *       TrilinosWrappers::SparseMatrix mass_matrix;
 *       mass_matrix.reinit(dsp);
 *       assemble_mass_matrix_diagonal(mass_matrix);
 *       diagonal_of_mass_matrix.reinit(solution_index_set);
 *       for (unsigned int j = 0; j < solution.size(); ++j)
 *         diagonal_of_mass_matrix(j) = mass_matrix.diag_element(j);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_41-ObstacleProblemassemble_system"></a> 
 * <h4>ObstacleProblem::assemble_system</h4>
 * 
 
 * 
 * This function at once assembles the system matrix and right-hand-side and
 * applied the constraints (both due to the active set as well as from
 * boundary values) to our system. Otherwise, it is functionally equivalent
 * to the corresponding function in, for example, @ref step_4 "step-4".
 * 
 * @code
 *     template <int dim>
 *     void ObstacleProblem<dim>::assemble_system()
 *     {
 *       std::cout << "   Assembling system..." << std::endl;
 *   
 *       system_matrix = 0;
 *       system_rhs    = 0;
 *   
 *       const QGauss<dim>  quadrature_formula(fe.degree + 1);
 *       RightHandSide<dim> right_hand_side;
 *   
 *       FEValues<dim> fe_values(fe,
 *                               quadrature_formula,
 *                               update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *   
 *       FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *       Vector<double>     cell_rhs(dofs_per_cell);
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *           cell_matrix = 0;
 *           cell_rhs    = 0;
 *   
 *           for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *             for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               {
 *                 for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                   cell_matrix(i, j) +=
 *                     (fe_values.shape_grad(i, q_point) *
 *                      fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
 *   
 *                 cell_rhs(i) +=
 *                   (fe_values.shape_value(i, q_point) *
 *                    right_hand_side.value(fe_values.quadrature_point(q_point)) *
 *                    fe_values.JxW(q_point));
 *               }
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *   
 *           constraints.distribute_local_to_global(cell_matrix,
 *                                                  cell_rhs,
 *                                                  local_dof_indices,
 *                                                  system_matrix,
 *                                                  system_rhs,
 *                                                  true);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_41-ObstacleProblemassemble_mass_matrix_diagonal"></a> 
 * <h4>ObstacleProblem::assemble_mass_matrix_diagonal</h4>
 * 
 
 * 
 * The next function is used in the computation of the diagonal mass matrix
 * @f$B@f$ used to scale variables in the active set method. As discussed in the
 * introduction, we get the mass matrix to be diagonal by choosing the
 * trapezoidal rule for quadrature. Doing so we don't really need the triple
 * loop over quadrature points, indices @f$i@f$ and indices @f$j@f$ any more and
 * can, instead, just use a double loop. The rest of the function is obvious
 * given what we have discussed in many of the previous tutorial programs.
 *   
 
 * 
 * Note that at the time this function is called, the constraints object
 * only contains boundary value constraints; we therefore do not have to pay
 * attention in the last copy-local-to-global step to preserve the values of
 * matrix entries that may later on be constrained by the active set.
 *   
 
 * 
 * Note also that the trick with the trapezoidal rule only works if we have
 * in fact @f$Q_1@f$ elements. For higher order elements, one would need to use
 * a quadrature formula that has quadrature points at all the support points
 * of the finite element. Constructing such a quadrature formula isn't
 * really difficult, but not the point here, and so we simply assert at the
 * top of the function that our implicit assumption about the finite element
 * is in fact satisfied.
 * 
 * @code
 *     template <int dim>
 *     void ObstacleProblem<dim>::assemble_mass_matrix_diagonal(
 *       TrilinosWrappers::SparseMatrix &mass_matrix)
 *     {
 *       Assert(fe.degree == 1, ExcNotImplemented());
 *   
 *       const QTrapezoid<dim> quadrature_formula;
 *       FEValues<dim>         fe_values(fe,
 *                               quadrature_formula,
 *                               update_values | update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points    = 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);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *           cell_matrix = 0;
 *   
 *           for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *             for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               cell_matrix(i, i) +=
 *                 (fe_values.shape_value(i, q_point) *
 *                  fe_values.shape_value(i, q_point) * fe_values.JxW(q_point));
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *   
 *           constraints.distribute_local_to_global(cell_matrix,
 *                                                  local_dof_indices,
 *                                                  mass_matrix);
 *         }
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_41-ObstacleProblemupdate_solution_and_constraints"></a> 
 * <h4>ObstacleProblem::update_solution_and_constraints</h4>
 * 
 
 * 
 * In a sense, this is the central function of this program.  It updates the
 * active set of constrained degrees of freedom as discussed in the
 * introduction and computes an AffineConstraints object from it that can then
 * be used to eliminate constrained degrees of freedom from the solution of
 * the next iteration. At the same time we set the constrained degrees of
 * freedom of the solution to the correct value, namely the height of the
 * obstacle.
 *   
 
 * 
 * Fundamentally, the function is rather simple: We have to loop over all
 * degrees of freedom and check the sign of the function @f$\Lambda^k_i +
 * c([BU^k]_i - G_i) = \Lambda^k_i + cB_i(U^k_i - [g_h]_i)@f$ because in our
 * case @f$G_i = B_i[g_h]_i@f$. To this end, we use the formula given in the
 * introduction by which we can compute the Lagrange multiplier as the
 * residual of the original linear system (given via the variables
 * <code>complete_system_matrix</code> and <code>complete_system_rhs</code>.
 * At the top of this function, we compute this residual using a function
 * that is part of the matrix classes.
 * 
 * @code
 *     template <int dim>
 *     void ObstacleProblem<dim>::update_solution_and_constraints()
 *     {
 *       std::cout << "   Updating active set..." << std::endl;
 *   
 *       const double penalty_parameter = 100.0;
 *   
 *       TrilinosWrappers::MPI::Vector lambda(
 *         complete_index_set(dof_handler.n_dofs()));
 *       complete_system_matrix.residual(lambda, solution, complete_system_rhs);
 *   
 * @endcode
 * 
 * compute contact_force[i] = - lambda[i] * diagonal_of_mass_matrix[i]
 * 
 * @code
 *       contact_force = lambda;
 *       contact_force.scale(diagonal_of_mass_matrix);
 *       contact_force *= -1;
 *   
 * @endcode
 * 
 * The next step is to reset the active set and constraints objects and to
 * start the loop over all degrees of freedom. This is made slightly more
 * complicated by the fact that we can't just loop over all elements of
 * the solution vector since there is no way for us then to find out what
 * location a DoF is associated with; however, we need this location to
 * test whether the displacement of a DoF is larger or smaller than the
 * height of the obstacle at this location.
 *     
 
 * 
 * We work around this by looping over all cells and DoFs defined on each
 * of these cells. We use here that the displacement is described using a
 * @f$Q_1@f$ function for which degrees of freedom are always located on the
 * vertices of the cell; thus, we can get the index of each degree of
 * freedom and its location by asking the vertex for this information. On
 * the other hand, this clearly wouldn't work for higher order elements,
 * and so we add an assertion that makes sure that we only deal with
 * elements for which all degrees of freedom are located in vertices to
 * avoid tripping ourselves with non-functional code in case someone wants
 * to play with increasing the polynomial degree of the solution.
 *     
 
 * 
 * The price to pay for having to loop over cells rather than DoFs is that
 * we may encounter some degrees of freedom more than once, namely each
 * time we visit one of the cells adjacent to a given vertex. We will
 * therefore have to keep track which vertices we have already touched and
 * which we haven't so far. We do so by using an array of flags
 * <code>dof_touched</code>:
 * 
 * @code
 *       constraints.clear();
 *       active_set.clear();
 *   
 *       const Obstacle<dim> obstacle;
 *       std::vector<bool>   dof_touched(dof_handler.n_dofs(), false);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         for (const auto v : cell->vertex_indices())
 *           {
 *             Assert(dof_handler.get_fe().n_dofs_per_cell() == cell->n_vertices(),
 *                    ExcNotImplemented());
 *   
 *             const unsigned int dof_index = cell->vertex_dof_index(v, 0);
 *   
 *             if (dof_touched[dof_index] == false)
 *               dof_touched[dof_index] = true;
 *             else
 *               continue;
 *   
 * @endcode
 * 
 * Now that we know that we haven't touched this DoF yet, let's get
 * the value of the displacement function there as well as the value
 * of the obstacle function and use this to decide whether the
 * current DoF belongs to the active set. For that we use the
 * function given above and in the introduction.
 *           
 
 * 
 * If we decide that the DoF should be part of the active set, we
 * add its index to the active set, introduce an inhomogeneous
 * equality constraint in the AffineConstraints object, and reset the
 * solution value to the height of the obstacle. Finally, the
 * residual of the non-contact part of the system serves as an
 * additional control (the residual equals the remaining,
 * unaccounted forces, and should be zero outside the contact zone),
 * so we zero out the components of the residual vector (i.e., the
 * Lagrange multiplier lambda) that correspond to the area where the
 * body is in contact; at the end of the loop over all cells, the
 * residual will therefore only consist of the residual in the
 * non-contact zone. We output the norm of this residual along with
 * the size of the active set after the loop.
 * 
 * @code
 *             const double obstacle_value = obstacle.value(cell->vertex(v));
 *             const double solution_value = solution(dof_index);
 *   
 *             if (lambda(dof_index) + penalty_parameter *
 *                                       diagonal_of_mass_matrix(dof_index) *
 *                                       (solution_value - obstacle_value) <
 *                 0)
 *               {
 *                 active_set.add_index(dof_index);
 *                 constraints.add_constraint(dof_index, {}, obstacle_value);
 *   
 *                 solution(dof_index) = obstacle_value;
 *   
 *                 lambda(dof_index) = 0;
 *               }
 *           }
 *       std::cout << "      Size of active set: " << active_set.n_elements()
 *                 << std::endl;
 *   
 *       std::cout << "   Residual of the non-contact part of the system: "
 *                 << lambda.l2_norm() << std::endl;
 *   
 * @endcode
 * 
 * In a final step, we add to the set of constraints on DoFs we have so
 * far from the active set those that result from Dirichlet boundary
 * values, and close the constraints object:
 * 
 * @code
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                0,
 *                                                BoundaryValues<dim>(),
 *                                                constraints);
 *       constraints.close();
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="step_41-ObstacleProblemsolve"></a> 
 * <h4>ObstacleProblem::solve</h4>
 * 
 
 * 
 * There is nothing to say really about the solve function. In the context
 * of a Newton method, we are not typically interested in very high accuracy
 * (why ask for a highly accurate solution of a linear problem that we know
 * only gives us an approximation of the solution of the nonlinear problem),
 * and so we use the ReductionControl class that stops iterations when
 * either an absolute tolerance is reached (for which we choose @f$10^{-12}@f$)
 * or when the residual is reduced by a certain factor (here, @f$10^{-3}@f$).
 * 
 * @code
 *     template <int dim>
 *     void ObstacleProblem<dim>::solve()
 *     {
 *       std::cout << "   Solving system..." << std::endl;
 *   
 *       ReductionControl                        reduction_control(100, 1e-12, 1e-3);
 *       SolverCG<TrilinosWrappers::MPI::Vector> solver(reduction_control);
 *       TrilinosWrappers::PreconditionAMG       precondition;
 *       precondition.initialize(system_matrix);
 *   
 *       solver.solve(system_matrix, solution, system_rhs, precondition);
 *       constraints.distribute(solution);
 *   
 *       std::cout << "      Error: " << reduction_control.initial_value() << " -> "
 *                 << reduction_control.last_value() << " in "
 *                 << reduction_control.last_step() << " CG iterations."
 *                 << std::endl;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_41-ObstacleProblemoutput_results"></a> 
 * <h4>ObstacleProblem::output_results</h4>
 * 
 
 * 
 * We use the vtk-format for the output.  The file contains the displacement
 * and a numerical representation of the active set.
 * 
 * @code
 *     template <int dim>
 *     void ObstacleProblem<dim>::output_results(const unsigned int iteration) const
 *     {
 *       std::cout << "   Writing graphical output..." << std::endl;
 *   
 *       TrilinosWrappers::MPI::Vector active_set_vector(
 *         dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
 *       for (const auto index : active_set)
 *         active_set_vector[index] = 1.;
 *   
 *       DataOut<dim> data_out;
 *   
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution, "displacement");
 *       data_out.add_data_vector(active_set_vector, "active_set");
 *       data_out.add_data_vector(contact_force, "lambda");
 *   
 *       data_out.build_patches();
 *   
 *       std::ofstream output_vtk("output_" +
 *                                Utilities::int_to_string(iteration, 3) + ".vtk");
 *       data_out.write_vtk(output_vtk);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_41-ObstacleProblemrun"></a> 
 * <h4>ObstacleProblem::run</h4>
 * 
 
 * 
 * This is the function which has the top-level control over everything.  It
 * is not very long, and in fact rather straightforward: in every iteration
 * of the active set method, we assemble the linear system, solve it, update
 * the active set and project the solution back to the feasible set, and
 * then output the results. The iteration is terminated whenever the active
 * set has not changed in the previous iteration.
 *   
 
 * 
 * The only trickier part is that we have to save the linear system (i.e.,
 * the matrix and right hand side) after assembling it in the first
 * iteration. The reason is that this is the only step where we can access
 * the linear system as built without any of the contact constraints
 * active. We need this to compute the residual of the solution at other
 * iterations, but in other iterations that linear system we form has the
 * rows and columns that correspond to constrained degrees of freedom
 * eliminated, and so we can no longer access the full residual of the
 * original equation.
 * 
 * @code
 *     template <int dim>
 *     void ObstacleProblem<dim>::run()
 *     {
 *       make_grid();
 *       setup_system();
 *   
 *       IndexSet active_set_old(active_set);
 *       for (unsigned int iteration = 0; iteration <= solution.size(); ++iteration)
 *         {
 *           std::cout << "Newton iteration " << iteration << std::endl;
 *   
 *           assemble_system();
 *   
 *           if (iteration == 0)
 *             {
 *               complete_system_matrix.copy_from(system_matrix);
 *               complete_system_rhs = system_rhs;
 *             }
 *   
 *           solve();
 *           update_solution_and_constraints();
 *           output_results(iteration);
 *   
 *           if (active_set == active_set_old)
 *             break;
 *   
 *           active_set_old = active_set;
 *   
 *           std::cout << std::endl;
 *         }
 *     }
 *   } // namespace Step41
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_41-Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * And this is the main function. It follows the pattern of all other main
 * functions. The call to initialize MPI exists because the Trilinos library
 * upon which we build our linear solvers in this program requires it.
 * 
 * @code
 *   int main(int argc, char *argv[])
 *   {
 *     try
 *       {
 *         using namespace dealii;
 *         using namespace Step41;
 *   
 *         Utilities::MPI::MPI_InitFinalize mpi_initialization(
 *           argc, argv, numbers::invalid_unsigned_int);
 *   
 * @endcode
 * 
 * This program can only be run in serial. Otherwise, throw an exception.
 * 
 * @code
 *         AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
 *                     ExcMessage(
 *                       "This program can only be run in serial, use ./step-41"));
 *   
 *         ObstacleProblem<2> obstacle_problem;
 *         obstacle_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="step_41-Results"></a><h1>Results</h1>
 
 
Running the program produces output like this:
@code
Number of active cells: 16384
Total number of cells: 21845
Number of degrees of freedom: 16641
 
Newton iteration 0
   Assembling system...
   Solving system...
      Error: 0.310059 -> 5.16619e-05 in 5 CG iterations.
   Updating active set...
      Size of active set: 13164
   Residual of the non-contact part of the system: 1.61863e-05
   Writing graphical output...
 
Newton iteration 1
   Assembling system...
   Solving system...
      Error: 1.11987 -> 0.00109377 in 6 CG iterations.
   Updating active set...
      Size of active set: 12363
   Residual of the non-contact part of the system: 3.9373
   Writing graphical output...
 
...
 
Newton iteration 17
   Assembling system...
   Solving system...
      Error: 0.00713308 -> 2.29249e-06 in 4 CG iterations.
   Updating active set...
      Size of active set: 5399
   Residual of the non-contact part of the system: 0.000957525
   Writing graphical output...
 
Newton iteration 18
   Assembling system...
   Solving system...
      Error: 0.000957525 -> 2.8033e-07 in 4 CG iterations.
   Updating active set...
      Size of active set: 5399
   Residual of the non-contact part of the system: 2.8033e-07
   Writing graphical output...
@endcode
 
The iterations end once the active set doesn't change any more (it has
5,399 constrained degrees of freedom at that point). The algebraic
precondition is apparently working nicely since we only need 4-6 CG
iterations to solve the linear system (although this also has a lot to
do with the fact that we are not asking for very high accuracy of the
linear solver).
 
More revealing is to look at a sequence of graphical output files
(every third step is shown, with the number of the iteration in the
leftmost column):
 
<table align="center">
  <tr>
    <td valign="top">
      0 &nbsp;
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.00.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.00.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.00.png" alt="">
    </td>
  </tr>
  <tr>
    <td valign="top">
      3 &nbsp;
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.03.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.03.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.03.png" alt="">
    </td>
  </tr>
  <tr>
    <td valign="top">
      6 &nbsp;
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.06.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.06.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.06.png" alt="">
    </td>
  </tr>
  <tr>
    <td valign="top">
      9 &nbsp;
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.09.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.09.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.09.png" alt="">
    </td>
  </tr>
  <tr>
    <td valign="top">
      12 &nbsp;
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.12.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.12.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.12.png" alt="">
    </td>
  </tr>
  <tr>
    <td valign="top">
      15 &nbsp;
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.15.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.15.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.15.png" alt="">
    </td>
  </tr>
  <tr>
    <td valign="top">
      18 &nbsp;
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.18.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.active-set.18.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.displacement.3d.18.png" alt="">
    </td>
  </tr>
</table>
 
The pictures show that in the first step, the solution (which has been
computed without any of the constraints active) bends through so much
that pretty much every interior point has to be bounced back to the
stairstep function, producing a discontinuous solution. Over the
course of the active set iterations, this unphysical membrane shape is
smoothed out, the contact with the lower-most stair step disappears,
and the solution stabilizes.
 
In addition to this, the program also outputs the values of the
Lagrange multipliers. Remember that these are the contact forces and
so should only be positive on the contact set, and zero outside. If,
on the other hand, a Lagrange multiplier is negative in the active
set, then this degree of freedom must be removed from the active
set. The following pictures show the multipliers in iterations 1, 9
and 18, where we use red and browns to indicate positive values, and
blue for negative values.
 
<table align="center">
  <tr>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.forces.01.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.forces.09.png" alt="">
    </td>
    <td valign="top">
      <img src="https://www.dealii.org/images/steps/developer/step-41.forces.18.png" alt="">
    </td>
  </tr>
  <tr>
    <td align="center">
      Iteration 1
    </td>
    <td align="center">
      Iteration 9
    </td>
    <td align="center">
      Iteration 18
    </td>
  </tr>
</table>
 
It is easy to see that the positive values converge nicely to moderate
values in the interior of the contact set and large upward forces at
the edges of the steps, as one would expect (to support the large
curvature of the membrane there); at the fringes of the active set,
multipliers are initially negative, causing the set to shrink until,
in iteration 18, there are no more negative multipliers and the
algorithm has converged.
 
 
 
<a name="step-41-extensions"></a>
<a name="step_41-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
As with any of the programs of this tutorial, there are a number of
obvious possibilities for extensions and experiments. The first one is
clear: introduce adaptivity. Contact problems are prime candidates for
adaptive meshes because the solution has lines along which it is less
regular (the places where contact is established between membrane and
obstacle) and other areas where the solution is very smooth (or, in
the present context, constant wherever it is in contact with the
obstacle). Adding this to the current program should not pose too many
difficulties, but it is not trivial to find a good error estimator for
that purpose.
 
A more challenging task would be an extension to 3d. The problem here
is not so much to simply make everything run in 3d. Rather, it is that
when a 3d body is deformed and gets into contact with an obstacle,
then the obstacle does not act as a constraining body force within the
domain as is the case here. Rather, the contact force only acts on the
boundary of the object. The inequality then is not in the differential
equation but in fact in the (Neumann-type) boundary conditions, though
this leads to a similar kind of variational
inequality. Mathematically, this means that the Lagrange multiplier
only lives on the surface, though it can of course be extended by zero
into the domain if that is convenient. As in the current program, one
does not need to form and store this Lagrange multiplier explicitly.
 
A further interesting problem for the 3d case is to consider contact problems
with friction. In almost every mechanical process friction has a big influence.
For the modelling we have to take into account tangential stresses at the contact
surface. Also we have to observe that friction adds another nonlinearity to
our problem.
 
Another nontrivial modification is to implement a more complex constitutive
law like nonlinear elasticity or elasto-plastic  material behavior.
The difficulty here is to handle the additional nonlinearity arising
through the nonlinear constitutive law.
 *
 *
<a name="step_41-PlainProg"></a>
<h1> The plain program</h1>
@include "step-41.cc"
*/