Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
step-42.h
Go to the documentation of this file.
1 
1623  * {
1624  * TimerOutput::Scope t(computing_timer, "Setup: distribute DoFs");
1625  * dof_handler.distribute_dofs(fe);
1626  *
1627  * locally_owned_dofs = dof_handler.locally_owned_dofs();
1628  * locally_relevant_dofs.clear();
1630  * locally_relevant_dofs);
1631  * }
1632  *
1633  * /* setup hanging nodes and Dirichlet constraints */
1634  * {
1635  * TimerOutput::Scope t(computing_timer, "Setup: constraints");
1636  * constraints_hanging_nodes.reinit(locally_relevant_dofs);
1638  * constraints_hanging_nodes);
1639  * constraints_hanging_nodes.close();
1640  *
1641  * pcout << " Number of active cells: "
1642  * << triangulation.n_global_active_cells() << std::endl
1643  * << " Number of degrees of freedom: " << dof_handler.n_dofs()
1644  * << std::endl;
1645  *
1646  * compute_dirichlet_constraints();
1647  * }
1648  *
1649  * /* initialization of vectors and the active set */
1650  * {
1651  * TimerOutput::Scope t(computing_timer, "Setup: vectors");
1652  * solution.reinit(locally_relevant_dofs, mpi_communicator);
1653  * newton_rhs.reinit(locally_owned_dofs, mpi_communicator);
1654  * newton_rhs_uncondensed.reinit(locally_owned_dofs, mpi_communicator);
1655  * diag_mass_matrix_vector.reinit(locally_owned_dofs, mpi_communicator);
1656  * fraction_of_plastic_q_points_per_cell.reinit(
1657  * triangulation.n_active_cells());
1658  *
1659  * active_set.clear();
1660  * active_set.set_size(dof_handler.n_dofs());
1661  * }
1662  *
1663  * @endcode
1664  *
1665  * Finally, we set up sparsity patterns and matrices.
1666  * We temporarily (ab)use the system matrix to also build the (diagonal)
1667  * matrix that we use in eliminating degrees of freedom that are in contact
1668  * with the obstacle, but we then immediately set the Newton matrix back
1669  * to zero.
1670  *
1671  * @code
1672  * {
1673  * TimerOutput::Scope t(computing_timer, "Setup: matrix");
1674  * TrilinosWrappers::SparsityPattern sp(locally_owned_dofs,
1675  * mpi_communicator);
1676  *
1677  * DoFTools::make_sparsity_pattern(dof_handler,
1678  * sp,
1679  * constraints_dirichlet_and_hanging_nodes,
1680  * false,
1682  * mpi_communicator));
1683  * sp.compress();
1684  * newton_matrix.reinit(sp);
1685  *
1686  *
1687  * TrilinosWrappers::SparseMatrix &mass_matrix = newton_matrix;
1688  *
1689  * assemble_mass_matrix_diagonal(mass_matrix);
1690  *
1691  * const unsigned int start = (newton_rhs.local_range().first),
1692  * end = (newton_rhs.local_range().second);
1693  * for (unsigned int j = start; j < end; ++j)
1694  * diag_mass_matrix_vector(j) = mass_matrix.diag_element(j);
1695  * diag_mass_matrix_vector.compress(VectorOperation::insert);
1696  *
1697  * mass_matrix = 0;
1698  * }
1699  * }
1700  *
1701  *
1702  * @endcode
1703  *
1704  *
1705  * <a name="PlasticityContactProblemcompute_dirichlet_constraints"></a>
1706  * <h4>PlasticityContactProblem::compute_dirichlet_constraints</h4>
1707  *
1708 
1709  *
1710  * This function, broken out of the preceding one, computes the constraints
1711  * associated with Dirichlet-type boundary conditions and puts them into the
1712  * <code>constraints_dirichlet_and_hanging_nodes</code> variable by merging
1713  * with the constraints that come from hanging nodes.
1714  *
1715 
1716  *
1717  * As laid out in the introduction, we need to distinguish between two
1718  * cases:
1719  * - If the domain is a box, we set the displacement to zero at the bottom,
1720  * and allow vertical movement in z-direction along the sides. As
1721  * shown in the <code>make_grid()</code> function, the former corresponds
1722  * to boundary indicator 6, the latter to 8.
1723  * - If the domain is a half sphere, then we impose zero displacement along
1724  * the curved part of the boundary, associated with boundary indicator zero.
1725  *
1726  * @code
1727  * template <int dim>
1728  * void PlasticityContactProblem<dim>::compute_dirichlet_constraints()
1729  * {
1730  * constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs);
1731  * constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);
1732  *
1733  * if (base_mesh == "box")
1734  * {
1735  * @endcode
1736  *
1737  * interpolate all components of the solution
1738  *
1739  * @code
1741  * dof_handler,
1742  * 6,
1743  * EquationData::BoundaryValues<dim>(),
1744  * constraints_dirichlet_and_hanging_nodes,
1745  * ComponentMask());
1746  *
1747  * @endcode
1748  *
1749  * interpolate x- and y-components of the
1750  * solution (this is a bit mask, so apply
1751  * operator| )
1752  *
1753  * @code
1754  * const FEValuesExtractors::Scalar x_displacement(0);
1755  * const FEValuesExtractors::Scalar y_displacement(1);
1757  * dof_handler,
1758  * 8,
1759  * EquationData::BoundaryValues<dim>(),
1760  * constraints_dirichlet_and_hanging_nodes,
1761  * (fe.component_mask(x_displacement) |
1762  * fe.component_mask(y_displacement)));
1763  * }
1764  * else
1766  * dof_handler,
1767  * 0,
1768  * EquationData::BoundaryValues<dim>(),
1769  * constraints_dirichlet_and_hanging_nodes,
1770  * ComponentMask());
1771  *
1772  * constraints_dirichlet_and_hanging_nodes.close();
1773  * }
1774  *
1775  *
1776  *
1777  * @endcode
1778  *
1779  *
1780  * <a name="PlasticityContactProblemassemble_mass_matrix_diagonal"></a>
1781  * <h4>PlasticityContactProblem::assemble_mass_matrix_diagonal</h4>
1782  *
1783 
1784  *
1785  * The next helper function computes the (diagonal) mass matrix that
1786  * is used to determine the active set of the active set method we use in
1787  * the contact algorithm. This matrix is of mass matrix type, but unlike
1788  * the standard mass matrix, we can make it diagonal (even in the case of
1789  * higher order elements) by using a quadrature formula that has its
1790  * quadrature points at exactly the same locations as the interpolation points
1791  * for the finite element are located. We achieve this by using a
1792  * QGaussLobatto quadrature formula here, along with initializing the finite
1793  * element with a set of interpolation points derived from the same quadrature
1794  * formula. The remainder of the function is relatively straightforward: we
1795  * put the resulting matrix into the given argument; because we know the
1796  * matrix is diagonal, it is sufficient to have a loop over only @f$i@f$ and
1797  * not over @f$j@f$. Strictly speaking, we could even avoid multiplying the
1798  * shape function's values at quadrature point <code>q_point</code> by itself
1799  * because we know the shape value to be a vector with exactly one one which
1800  * when dotted with itself yields one. Since this function is not time
1801  * critical we add this term for clarity.
1802  *
1803  * @code
1804  * template <int dim>
1805  * void PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal(
1806  * TrilinosWrappers::SparseMatrix &mass_matrix)
1807  * {
1808  * QGaussLobatto<dim - 1> face_quadrature_formula(fe.degree + 1);
1809  *
1810  * FEFaceValues<dim> fe_values_face(fe,
1811  * face_quadrature_formula,
1812  * update_values | update_JxW_values);
1813  *
1814  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1815  * const unsigned int n_face_q_points = face_quadrature_formula.size();
1816  *
1817  * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1818  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1819  *
1820  * const FEValuesExtractors::Vector displacement(0);
1821  *
1822  * for (const auto &cell : dof_handler.active_cell_iterators())
1823  * if (cell->is_locally_owned())
1824  * for (const auto &face : cell->face_iterators())
1825  * if (face->at_boundary() && face->boundary_id() == 1)
1826  * {
1827  * fe_values_face.reinit(cell, face);
1828  * cell_matrix = 0;
1829  *
1830  * for (unsigned int q_point = 0; q_point < n_face_q_points;
1831  * ++q_point)
1832  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1833  * cell_matrix(i, i) +=
1834  * (fe_values_face[displacement].value(i, q_point) *
1835  * fe_values_face[displacement].value(i, q_point) *
1836  * fe_values_face.JxW(q_point));
1837  *
1838  * cell->get_dof_indices(local_dof_indices);
1839  *
1840  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1841  * mass_matrix.add(local_dof_indices[i],
1842  * local_dof_indices[i],
1843  * cell_matrix(i, i));
1844  * }
1845  * mass_matrix.compress(VectorOperation::add);
1846  * }
1847  *
1848  *
1849  * @endcode
1850  *
1851  *
1852  * <a name="PlasticityContactProblemupdate_solution_and_constraints"></a>
1853  * <h4>PlasticityContactProblem::update_solution_and_constraints</h4>
1854  *
1855 
1856  *
1857  * The following function is the first function we call in each Newton
1858  * iteration in the <code>solve_newton()</code> function. What it does is
1859  * to project the solution onto the feasible set and update the active set
1860  * for the degrees of freedom that touch or penetrate the obstacle.
1861  *
1862 
1863  *
1864  * In order to function, we first need to do some bookkeeping: We need
1865  * to write into the solution vector (which we can only do with fully
1866  * distributed vectors without ghost elements) and we need to read
1867  * the Lagrange multiplier and the elements of the diagonal mass matrix
1868  * from their respective vectors (which we can only do with vectors that
1869  * do have ghost elements), so we create the respective vectors. We then
1870  * also initialize the constraints object that will contain constraints
1871  * from contact and all other sources, as well as an object that contains
1872  * an index set of all locally owned degrees of freedom that are part of
1873  * the contact:
1874  *
1875  * @code
1876  * template <int dim>
1877  * void PlasticityContactProblem<dim>::update_solution_and_constraints()
1878  * {
1879  * std::vector<bool> dof_touched(dof_handler.n_dofs(), false);
1880  *
1881  * TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
1882  * mpi_communicator);
1883  * distributed_solution = solution;
1884  *
1885  * TrilinosWrappers::MPI::Vector lambda(locally_relevant_dofs,
1886  * mpi_communicator);
1887  * lambda = newton_rhs_uncondensed;
1888  *
1889  * TrilinosWrappers::MPI::Vector diag_mass_matrix_vector_relevant(
1890  * locally_relevant_dofs, mpi_communicator);
1891  * diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
1892  *
1893  *
1894  * all_constraints.reinit(locally_relevant_dofs);
1895  * active_set.clear();
1896  *
1897  * @endcode
1898  *
1899  * The second part is a loop over all cells in which we look at each
1900  * point where a degree of freedom is defined whether the active set
1901  * condition is true and we need to add this degree of freedom to
1902  * the active set of contact nodes. As we always do, if we want to
1903  * evaluate functions at individual points, we do this with an
1904  * FEValues object (or, here, an FEFaceValues object since we need to
1905  * check contact at the surface) with an appropriately chosen quadrature
1906  * object. We create this face quadrature object by choosing the
1907  * "support points" of the shape functions defined on the faces
1908  * of cells (for more on support points, see this
1909  * @ref GlossSupport "glossary entry"). As a consequence, we have as
1910  * many quadrature points as there are shape functions per face and
1911  * looping over quadrature points is equivalent to looping over shape
1912  * functions defined on a face. With this, the code looks as follows:
1913  *
1914  * @code
1915  * Quadrature<dim - 1> face_quadrature(fe.get_unit_face_support_points());
1916  * FEFaceValues<dim> fe_values_face(fe,
1917  * face_quadrature,
1918  * update_quadrature_points);
1919  *
1920  * const unsigned int dofs_per_face = fe.dofs_per_face;
1921  * const unsigned int n_face_q_points = face_quadrature.size();
1922  *
1923  * std::vector<types::global_dof_index> dof_indices(dofs_per_face);
1924  *
1925  * for (const auto &cell : dof_handler.active_cell_iterators())
1926  * if (!cell->is_artificial())
1927  * for (const auto &face : cell->face_iterators())
1928  * if (face->at_boundary() && face->boundary_id() == 1)
1929  * {
1930  * fe_values_face.reinit(cell, face);
1931  * face->get_dof_indices(dof_indices);
1932  *
1933  * for (unsigned int q_point = 0; q_point < n_face_q_points;
1934  * ++q_point)
1935  * {
1936  * @endcode
1937  *
1938  * At each quadrature point (i.e., at each support point of a
1939  * degree of freedom located on the contact boundary), we then
1940  * ask whether it is part of the z-displacement degrees of
1941  * freedom and if we haven't encountered this degree of
1942  * freedom yet (which can happen for those on the edges
1943  * between faces), we need to evaluate the gap between the
1944  * deformed object and the obstacle. If the active set
1945  * condition is true, then we add a constraint to the
1946  * AffineConstraints object that the next Newton update needs
1947  * to satisfy, set the solution vector's corresponding element
1948  * to the correct value, and add the index to the IndexSet
1949  * object that stores which degree of freedom is part of the
1950  * contact:
1951  *
1952  * @code
1953  * const unsigned int component =
1954  * fe.face_system_to_component_index(q_point).first;
1955  *
1956  * const unsigned int index_z = dof_indices[q_point];
1957  *
1958  * if ((component == 2) && (dof_touched[index_z] == false))
1959  * {
1960  * dof_touched[index_z] = true;
1961  *
1962  * const Point<dim> this_support_point =
1963  * fe_values_face.quadrature_point(q_point);
1964  *
1965  * const double obstacle_value =
1966  * obstacle->value(this_support_point, 2);
1967  * const double solution_here = solution(index_z);
1968  * const double undeformed_gap =
1969  * obstacle_value - this_support_point(2);
1970  *
1971  * const double c = 100.0 * e_modulus;
1972  * if ((lambda(index_z) /
1973  * diag_mass_matrix_vector_relevant(index_z) +
1974  * c * (solution_here - undeformed_gap) >
1975  * 0) &&
1976  * !constraints_hanging_nodes.is_constrained(index_z))
1977  * {
1978  * all_constraints.add_line(index_z);
1979  * all_constraints.set_inhomogeneity(index_z,
1980  * undeformed_gap);
1981  * distributed_solution(index_z) = undeformed_gap;
1982  *
1983  * active_set.add_index(index_z);
1984  * }
1985  * }
1986  * }
1987  * }
1988  *
1989  * @endcode
1990  *
1991  * At the end of this function, we exchange data between processors updating
1992  * those ghost elements in the <code>solution</code> variable that have been
1993  * written by other processors. We then merge the Dirichlet constraints and
1994  * those from hanging nodes into the AffineConstraints object that already
1995  * contains the active set. We finish the function by outputting the total
1996  * number of actively constrained degrees of freedom for which we sum over
1997  * the number of actively constrained degrees of freedom owned by each
1998  * of the processors. This number of locally owned constrained degrees of
1999  * freedom is of course the number of elements of the intersection of the
2000  * active set and the set of locally owned degrees of freedom, which
2001  * we can get by using <code>operator&</code> on two IndexSets:
2002  *
2003  * @code
2004  * distributed_solution.compress(VectorOperation::insert);
2005  * solution = distributed_solution;
2006  *
2007  * all_constraints.close();
2008  * all_constraints.merge(constraints_dirichlet_and_hanging_nodes);
2009  *
2010  * pcout << " Size of active set: "
2011  * << Utilities::MPI::sum((active_set & locally_owned_dofs).n_elements(),
2012  * mpi_communicator)
2013  * << std::endl;
2014  * }
2015  *
2016  *
2017  * @endcode
2018  *
2019  *
2020  * <a name="PlasticityContactProblemassemble_newton_system"></a>
2021  * <h4>PlasticityContactProblem::assemble_newton_system</h4>
2022  *
2023 
2024  *
2025  * Given the complexity of the problem, it may come as a bit of a surprise
2026  * that assembling the linear system we have to solve in each Newton iteration
2027  * is actually fairly straightforward. The following function builds the
2028  * Newton right hand side and Newton matrix. It looks fairly innocent because
2029  * the heavy lifting happens in the call to
2030  * <code>ConstitutiveLaw::get_linearized_stress_strain_tensors()</code> and in
2031  * particular in AffineConstraints::distribute_local_to_global(), using the
2032  * constraints we have previously computed.
2033  *
2034  * @code
2035  * template <int dim>
2036  * void PlasticityContactProblem<dim>::assemble_newton_system(
2037  * const TrilinosWrappers::MPI::Vector &linearization_point)
2038  * {
2039  * TimerOutput::Scope t(computing_timer, "Assembling");
2040  *
2041  * QGauss<dim> quadrature_formula(fe.degree + 1);
2042  * QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
2043  *
2044  * FEValues<dim> fe_values(fe,
2045  * quadrature_formula,
2046  * update_values | update_gradients |
2047  * update_JxW_values);
2048  *
2049  * FEFaceValues<dim> fe_values_face(fe,
2050  * face_quadrature_formula,
2051  * update_values | update_quadrature_points |
2052  * update_JxW_values);
2053  *
2054  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
2055  * const unsigned int n_q_points = quadrature_formula.size();
2056  * const unsigned int n_face_q_points = face_quadrature_formula.size();
2057  *
2058  * const EquationData::BoundaryForce<dim> boundary_force;
2059  * std::vector<Vector<double>> boundary_force_values(n_face_q_points,
2060  * Vector<double>(dim));
2061  *
2062  * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
2063  * Vector<double> cell_rhs(dofs_per_cell);
2064  *
2065  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2066  *
2067  * const FEValuesExtractors::Vector displacement(0);
2068  *
2069  * for (const auto &cell : dof_handler.active_cell_iterators())
2070  * if (cell->is_locally_owned())
2071  * {
2072  * fe_values.reinit(cell);
2073  * cell_matrix = 0;
2074  * cell_rhs = 0;
2075  *
2076  * std::vector<SymmetricTensor<2, dim>> strain_tensor(n_q_points);
2077  * fe_values[displacement].get_function_symmetric_gradients(
2078  * linearization_point, strain_tensor);
2079  *
2080  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2081  * {
2082  * SymmetricTensor<4, dim> stress_strain_tensor_linearized;
2083  * SymmetricTensor<4, dim> stress_strain_tensor;
2084  * constitutive_law.get_linearized_stress_strain_tensors(
2085  * strain_tensor[q_point],
2086  * stress_strain_tensor_linearized,
2087  * stress_strain_tensor);
2088  *
2089  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2090  * {
2091  * @endcode
2092  *
2093  * Having computed the stress-strain tensor and its
2094  * linearization, we can now put together the parts of the
2095  * matrix and right hand side. In both, we need the linearized
2096  * stress-strain tensor times the symmetric gradient of
2097  * @f$\varphi_i@f$, i.e. the term @f$I_\Pi\varepsilon(\varphi_i)@f$,
2098  * so we introduce an abbreviation of this term. Recall that
2099  * the matrix corresponds to the bilinear form
2100  * @f$A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))@f$
2101  * in the notation of the accompanying publication, whereas
2102  * the right hand side is @f$F_i=([I_\Pi-P_\Pi
2103  * C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))@f$ where @f$u@f$
2104  * is the current linearization points (typically the last
2105  * solution). This might suggest that the right hand side will
2106  * be zero if the material is completely elastic (where
2107  * @f$I_\Pi=P_\Pi@f$) but this ignores the fact that the right
2108  * hand side will also contain contributions from
2109  * non-homogeneous constraints due to the contact.
2110  *
2111 
2112  *
2113  * The code block that follows this adds contributions that
2114  * are due to boundary forces, should there be any.
2115  *
2116  * @code
2117  * const SymmetricTensor<2, dim> stress_phi_i =
2118  * stress_strain_tensor_linearized *
2119  * fe_values[displacement].symmetric_gradient(i, q_point);
2120  *
2121  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
2122  * cell_matrix(i, j) +=
2123  * (stress_phi_i *
2124  * fe_values[displacement].symmetric_gradient(j, q_point) *
2125  * fe_values.JxW(q_point));
2126  *
2127  * cell_rhs(i) +=
2128  * ((stress_phi_i -
2129  * stress_strain_tensor *
2130  * fe_values[displacement].symmetric_gradient(i,
2131  * q_point)) *
2132  * strain_tensor[q_point] * fe_values.JxW(q_point));
2133  * }
2134  * }
2135  *
2136  * for (const auto &face : cell->face_iterators())
2137  * if (face->at_boundary() && face->boundary_id() == 1)
2138  * {
2139  * fe_values_face.reinit(cell, face);
2140  *
2141  * boundary_force.vector_value_list(
2142  * fe_values_face.get_quadrature_points(),
2143  * boundary_force_values);
2144  *
2145  * for (unsigned int q_point = 0; q_point < n_face_q_points;
2146  * ++q_point)
2147  * {
2148  * Tensor<1, dim> rhs_values;
2149  * rhs_values[2] = boundary_force_values[q_point][2];
2150  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2151  * cell_rhs(i) +=
2152  * (fe_values_face[displacement].value(i, q_point) *
2153  * rhs_values * fe_values_face.JxW(q_point));
2154  * }
2155  * }
2156  *
2157  * cell->get_dof_indices(local_dof_indices);
2158  * all_constraints.distribute_local_to_global(cell_matrix,
2159  * cell_rhs,
2160  * local_dof_indices,
2161  * newton_matrix,
2162  * newton_rhs,
2163  * true);
2164  * }
2165  *
2166  * newton_matrix.compress(VectorOperation::add);
2167  * newton_rhs.compress(VectorOperation::add);
2168  * }
2169  *
2170  *
2171  *
2172  * @endcode
2173  *
2174  *
2175  * <a name="PlasticityContactProblemcompute_nonlinear_residual"></a>
2176  * <h4>PlasticityContactProblem::compute_nonlinear_residual</h4>
2177  *
2178 
2179  *
2180  * The following function computes the nonlinear residual of the equation
2181  * given the current solution (or any other linearization point). This
2182  * is needed in the linear search algorithm where we need to try various
2183  * linear combinations of previous and current (trial) solution to
2184  * compute the (real, globalized) solution of the current Newton step.
2185  *
2186 
2187  *
2188  * That said, in a slight abuse of the name of the function, it actually
2189  * does significantly more. For example, it also computes the vector
2190  * that corresponds to the Newton residual but without eliminating
2191  * constrained degrees of freedom. We need this vector to compute contact
2192  * forces and, ultimately, to compute the next active set. Likewise, by
2193  * keeping track of how many quadrature points we encounter on each cell
2194  * that show plastic yielding, we also compute the
2195  * <code>fraction_of_plastic_q_points_per_cell</code> vector that we
2196  * can later output to visualize the plastic zone. In both of these cases,
2197  * the results are not necessary as part of the line search, and so we may
2198  * be wasting a small amount of time computing them. At the same time, this
2199  * information appears as a natural by-product of what we need to do here
2200  * anyway, and we want to collect it once at the end of each Newton
2201  * step, so we may as well do it here.
2202  *
2203 
2204  *
2205  * The actual implementation of this function should be rather obvious:
2206  *
2207  * @code
2208  * template <int dim>
2209  * void PlasticityContactProblem<dim>::compute_nonlinear_residual(
2210  * const TrilinosWrappers::MPI::Vector &linearization_point)
2211  * {
2212  * QGauss<dim> quadrature_formula(fe.degree + 1);
2213  * QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
2214  *
2215  * FEValues<dim> fe_values(fe,
2216  * quadrature_formula,
2217  * update_values | update_gradients |
2218  * update_JxW_values);
2219  *
2220  * FEFaceValues<dim> fe_values_face(fe,
2221  * face_quadrature_formula,
2222  * update_values | update_quadrature_points |
2223  * update_JxW_values);
2224  *
2225  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
2226  * const unsigned int n_q_points = quadrature_formula.size();
2227  * const unsigned int n_face_q_points = face_quadrature_formula.size();
2228  *
2229  * const EquationData::BoundaryForce<dim> boundary_force;
2230  * std::vector<Vector<double>> boundary_force_values(n_face_q_points,
2231  * Vector<double>(dim));
2232  *
2233  * Vector<double> cell_rhs(dofs_per_cell);
2234  *
2235  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2236  *
2237  * const FEValuesExtractors::Vector displacement(0);
2238  *
2239  * newton_rhs = 0;
2240  * newton_rhs_uncondensed = 0;
2241  *
2242  * fraction_of_plastic_q_points_per_cell = 0;
2243  *
2244  * for (const auto &cell : dof_handler.active_cell_iterators())
2245  * if (cell->is_locally_owned())
2246  * {
2247  * fe_values.reinit(cell);
2248  * cell_rhs = 0;
2249  *
2250  * std::vector<SymmetricTensor<2, dim>> strain_tensors(n_q_points);
2251  * fe_values[displacement].get_function_symmetric_gradients(
2252  * linearization_point, strain_tensors);
2253  *
2254  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2255  * {
2256  * SymmetricTensor<4, dim> stress_strain_tensor;
2257  * const bool q_point_is_plastic =
2258  * constitutive_law.get_stress_strain_tensor(
2259  * strain_tensors[q_point], stress_strain_tensor);
2260  * if (q_point_is_plastic)
2261  * ++fraction_of_plastic_q_points_per_cell(
2262  * cell->active_cell_index());
2263  *
2264  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2265  * {
2266  * cell_rhs(i) -=
2267  * (strain_tensors[q_point] * stress_strain_tensor *
2268  * fe_values[displacement].symmetric_gradient(i, q_point) *
2269  * fe_values.JxW(q_point));
2270  *
2271  * Tensor<1, dim> rhs_values;
2272  * rhs_values = 0;
2273  * cell_rhs(i) += (fe_values[displacement].value(i, q_point) *
2274  * rhs_values * fe_values.JxW(q_point));
2275  * }
2276  * }
2277  *
2278  * for (const auto &face : cell->face_iterators())
2279  * if (face->at_boundary() && face->boundary_id() == 1)
2280  * {
2281  * fe_values_face.reinit(cell, face);
2282  *
2283  * boundary_force.vector_value_list(
2284  * fe_values_face.get_quadrature_points(),
2285  * boundary_force_values);
2286  *
2287  * for (unsigned int q_point = 0; q_point < n_face_q_points;
2288  * ++q_point)
2289  * {
2290  * Tensor<1, dim> rhs_values;
2291  * rhs_values[2] = boundary_force_values[q_point][2];
2292  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2293  * cell_rhs(i) +=
2294  * (fe_values_face[displacement].value(i, q_point) *
2295  * rhs_values * fe_values_face.JxW(q_point));
2296  * }
2297  * }
2298  *
2299  * cell->get_dof_indices(local_dof_indices);
2300  * constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(
2301  * cell_rhs, local_dof_indices, newton_rhs);
2302  *
2303  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2304  * newton_rhs_uncondensed(local_dof_indices[i]) += cell_rhs(i);
2305  * }
2306  *
2307  * fraction_of_plastic_q_points_per_cell /= quadrature_formula.size();
2308  * newton_rhs.compress(VectorOperation::add);
2309  * newton_rhs_uncondensed.compress(VectorOperation::add);
2310  * }
2311  *
2312  *
2313  *
2314  * @endcode
2315  *
2316  *
2317  * <a name="PlasticityContactProblemsolve_newton_system"></a>
2318  * <h4>PlasticityContactProblem::solve_newton_system</h4>
2319  *
2320 
2321  *
2322  * The last piece before we can discuss the actual Newton iteration
2323  * on a single mesh is the solver for the linear systems. There are
2324  * a couple of complications that slightly obscure the code, but
2325  * mostly it is just setup then solve. Among the complications are:
2326  *
2327 
2328  *
2329  * - For the hanging nodes we have to apply
2330  * the AffineConstraints::set_zero function to newton_rhs.
2331  * This is necessary if a hanging node with solution value @f$x_0@f$
2332  * has one neighbor with value @f$x_1@f$ which is in contact with the
2333  * obstacle and one neighbor @f$x_2@f$ which is not in contact. Because
2334  * the update for the former will be prescribed, the hanging node constraint
2335  * will have an inhomogeneity and will look like @f$x_0 = x_1/2 +
2336  * \text{gap}/2@f$. So the corresponding entries in the right-hand-side are
2337  * non-zero with a meaningless value. These values we have to set to zero.
2338  * - Like in @ref step_40 "step-40", we need to shuffle between vectors that do and do
2339  * not have ghost elements when solving or using the solution.
2340  *
2341 
2342  *
2343  * The rest of the function is similar to @ref step_40 "step-40" and
2344  * @ref step_41 "step-41" except that we use a BiCGStab solver
2345  * instead of CG. This is due to the fact that for very small hardening
2346  * parameters @f$\gamma@f$, the linear system becomes almost semidefinite though
2347  * still symmetric. BiCGStab appears to have an easier time with such linear
2348  * systems.
2349  *
2350  * @code
2351  * template <int dim>
2352  * void PlasticityContactProblem<dim>::solve_newton_system()
2353  * {
2354  * TimerOutput::Scope t(computing_timer, "Solve");
2355  *
2356  * TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
2357  * mpi_communicator);
2358  * distributed_solution = solution;
2359  *
2360  * constraints_hanging_nodes.set_zero(distributed_solution);
2361  * constraints_hanging_nodes.set_zero(newton_rhs);
2362  *
2363  * TrilinosWrappers::PreconditionAMG preconditioner;
2364  * {
2365  * TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner");
2366  *
2367  * std::vector<std::vector<bool>> constant_modes;
2368  * DoFTools::extract_constant_modes(dof_handler,
2369  * ComponentMask(),
2370  * constant_modes);
2371  *
2372  * TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
2373  * additional_data.constant_modes = constant_modes;
2374  * additional_data.elliptic = true;
2375  * additional_data.n_cycles = 1;
2376  * additional_data.w_cycle = false;
2377  * additional_data.output_details = false;
2378  * additional_data.smoother_sweeps = 2;
2379  * additional_data.aggregation_threshold = 1e-2;
2380  *
2381  * preconditioner.initialize(newton_matrix, additional_data);
2382  * }
2383  *
2384  * {
2385  * TimerOutput::Scope t(computing_timer, "Solve: iterate");
2386  *
2387  * TrilinosWrappers::MPI::Vector tmp(locally_owned_dofs, mpi_communicator);
2388  *
2389  * const double relative_accuracy = 1e-8;
2390  * const double solver_tolerance =
2391  * relative_accuracy *
2392  * newton_matrix.residual(tmp, distributed_solution, newton_rhs);
2393  *
2394  * SolverControl solver_control(newton_matrix.m(), solver_tolerance);
2395  * SolverBicgstab<TrilinosWrappers::MPI::Vector> solver(solver_control);
2396  * solver.solve(newton_matrix,
2397  * distributed_solution,
2398  * newton_rhs,
2399  * preconditioner);
2400  *
2401  * pcout << " Error: " << solver_control.initial_value() << " -> "
2402  * << solver_control.last_value() << " in "
2403  * << solver_control.last_step() << " Bicgstab iterations."
2404  * << std::endl;
2405  * }
2406  *
2407  * all_constraints.distribute(distributed_solution);
2408  *
2409  * solution = distributed_solution;
2410  * }
2411  *
2412  *
2413  * @endcode
2414  *
2415  *
2416  * <a name="PlasticityContactProblemsolve_newton"></a>
2417  * <h4>PlasticityContactProblem::solve_newton</h4>
2418  *
2419 
2420  *
2421  * This is, finally, the function that implements the damped Newton method
2422  * on the current mesh. There are two nested loops: the outer loop for the
2423  * Newton iteration and the inner loop for the line search which will be used
2424  * only if necessary. To obtain a good and reasonable starting value we solve
2425  * an elastic problem in the very first Newton step on each mesh (or only on
2426  * the first mesh if we transfer solutions between meshes). We do so by
2427  * setting the yield stress to an unreasonably large value in these iterations
2428  * and then setting it back to the correct value in subsequent iterations.
2429  *
2430 
2431  *
2432  * Other than this, the top part of this function should be
2433  * reasonably obvious. We initialize the variable
2434  * <code>previous_residual_norm</code> to the most negative value
2435  * representable with double precision numbers so that the
2436  * comparison whether the current residual is less than that of the
2437  * previous step will always fail in the first step.
2438  *
2439  * @code
2440  * template <int dim>
2441  * void PlasticityContactProblem<dim>::solve_newton()
2442  * {
2443  * TrilinosWrappers::MPI::Vector old_solution(locally_owned_dofs,
2444  * mpi_communicator);
2445  * TrilinosWrappers::MPI::Vector residual(locally_owned_dofs,
2446  * mpi_communicator);
2447  * TrilinosWrappers::MPI::Vector tmp_vector(locally_owned_dofs,
2448  * mpi_communicator);
2449  * TrilinosWrappers::MPI::Vector locally_relevant_tmp_vector(
2450  * locally_relevant_dofs, mpi_communicator);
2451  * TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
2452  * mpi_communicator);
2453  *
2454  * double residual_norm;
2455  * double previous_residual_norm = -std::numeric_limits<double>::max();
2456  *
2457  * const double correct_sigma = sigma_0;
2458  *
2459  * IndexSet old_active_set(active_set);
2460  *
2461  * for (unsigned int newton_step = 1; newton_step <= 100; ++newton_step)
2462  * {
2463  * if (newton_step == 1 &&
2464  * ((transfer_solution && current_refinement_cycle == 0) ||
2465  * !transfer_solution))
2466  * constitutive_law.set_sigma_0(1e+10);
2467  * else if (newton_step == 2 || current_refinement_cycle > 0 ||
2468  * !transfer_solution)
2469  * constitutive_law.set_sigma_0(correct_sigma);
2470  *
2471  * pcout << " " << std::endl;
2472  * pcout << " Newton iteration " << newton_step << std::endl;
2473  * pcout << " Updating active set..." << std::endl;
2474  *
2475  * {
2476  * TimerOutput::Scope t(computing_timer, "update active set");
2477  * update_solution_and_constraints();
2478  * }
2479  *
2480  * pcout << " Assembling system... " << std::endl;
2481  * newton_matrix = 0;
2482  * newton_rhs = 0;
2483  * assemble_newton_system(solution);
2484  *
2485  * pcout << " Solving system... " << std::endl;
2486  * solve_newton_system();
2487  *
2488  * @endcode
2489  *
2490  * It gets a bit more hairy after we have computed the
2491  * trial solution @f$\tilde{\mathbf u}@f$ of the current Newton step.
2492  * We handle a highly nonlinear problem so we have to damp
2493  * Newton's method using a line search. To understand how we do this,
2494  * recall that in our formulation, we compute a trial solution
2495  * in each Newton step and not the update between old and new solution.
2496  * Since the solution set is a convex set, we will use a line
2497  * search that tries linear combinations of the
2498  * previous and the trial solution to guarantee that the
2499  * damped solution is in our solution set again.
2500  * At most we apply 5 damping steps.
2501  *
2502 
2503  *
2504  * There are exceptions to when we use a line search. First,
2505  * if this is the first Newton step on any mesh, then we don't have
2506  * any point to compare the residual to, so we always accept a full
2507  * step. Likewise, if this is the second Newton step on the first mesh
2508  * (or the second on any mesh if we don't transfer solutions from mesh
2509  * to mesh), then we have computed the first of these steps using just
2510  * an elastic model (see how we set the yield stress sigma to an
2511  * unreasonably large value above). In this case, the first Newton
2512  * solution was a purely elastic one, the second one a plastic one,
2513  * and any linear combination would not necessarily be expected to
2514  * lie in the feasible set -- so we just accept the solution we just
2515  * got.
2516  *
2517 
2518  *
2519  * In either of these two cases, we bypass the line search and just
2520  * update residual and other vectors as necessary.
2521  *
2522  * @code
2523  * if ((newton_step == 1) ||
2524  * (transfer_solution && newton_step == 2 &&
2525  * current_refinement_cycle == 0) ||
2526  * (!transfer_solution && newton_step == 2))
2527  * {
2528  * compute_nonlinear_residual(solution);
2529  * old_solution = solution;
2530  *
2531  * residual = newton_rhs;
2532  * const unsigned int start_res = (residual.local_range().first),
2533  * end_res = (residual.local_range().second);
2534  * for (unsigned int n = start_res; n < end_res; ++n)
2535  * if (all_constraints.is_inhomogeneously_constrained(n))
2536  * residual(n) = 0;
2537  *
2538  * residual.compress(VectorOperation::insert);
2539  *
2540  * residual_norm = residual.l2_norm();
2541  *
2542  * pcout << " Accepting Newton solution with residual: "
2543  * << residual_norm << std::endl;
2544  * }
2545  * else
2546  * {
2547  * for (unsigned int i = 0; i < 5; ++i)
2548  * {
2549  * distributed_solution = solution;
2550  *
2551  * const double alpha = std::pow(0.5, static_cast<double>(i));
2552  * tmp_vector = old_solution;
2553  * tmp_vector.sadd(1 - alpha, alpha, distributed_solution);
2554  *
2555  * TimerOutput::Scope t(computing_timer, "Residual and lambda");
2556  *
2557  * locally_relevant_tmp_vector = tmp_vector;
2558  * compute_nonlinear_residual(locally_relevant_tmp_vector);
2559  * residual = newton_rhs;
2560  *
2561  * const unsigned int start_res = (residual.local_range().first),
2562  * end_res = (residual.local_range().second);
2563  * for (unsigned int n = start_res; n < end_res; ++n)
2564  * if (all_constraints.is_inhomogeneously_constrained(n))
2565  * residual(n) = 0;
2566  *
2567  * residual.compress(VectorOperation::insert);
2568  *
2569  * residual_norm = residual.l2_norm();
2570  *
2571  * pcout
2572  * << " Residual of the non-contact part of the system: "
2573  * << residual_norm << std::endl
2574  * << " with a damping parameter alpha = " << alpha
2575  * << std::endl;
2576  *
2577  * if (residual_norm < previous_residual_norm)
2578  * break;
2579  * }
2580  *
2581  * solution = tmp_vector;
2582  * old_solution = solution;
2583  * }
2584  *
2585  * previous_residual_norm = residual_norm;
2586  *
2587  *
2588  * @endcode
2589  *
2590  * The final step is to check for convergence. If the active set
2591  * has not changed across all processors and the residual is
2592  * less than a threshold of @f$10^{-10}@f$, then we terminate
2593  * the iteration on the current mesh:
2594  *
2595  * @code
2596  * if (Utilities::MPI::sum((active_set == old_active_set) ? 0 : 1,
2597  * mpi_communicator) == 0)
2598  * {
2599  * pcout << " Active set did not change!" << std::endl;
2600  * if (residual_norm < 1e-10)
2601  * break;
2602  * }
2603  *
2604  * old_active_set = active_set;
2605  * }
2606  * }
2607  *
2608  * @endcode
2609  *
2610  *
2611  * <a name="PlasticityContactProblemrefine_grid"></a>
2612  * <h4>PlasticityContactProblem::refine_grid</h4>
2613  *
2614 
2615  *
2616  * If you've made it this far into the deal.II tutorial, the following
2617  * function refining the mesh should not pose any challenges to you
2618  * any more. It refines the mesh, either globally or using the Kelly
2619  * error estimator, and if so asked also transfers the solution from
2620  * the previous to the next mesh. In the latter case, we also need
2621  * to compute the active set and other quantities again, for which we
2622  * need the information computed by <code>compute_nonlinear_residual()</code>.
2623  *
2624  * @code
2625  * template <int dim>
2626  * void PlasticityContactProblem<dim>::refine_grid()
2627  * {
2628  * if (refinement_strategy == RefinementStrategy::refine_global)
2629  * {
2630  * for (typename Triangulation<dim>::active_cell_iterator cell =
2631  * triangulation.begin_active();
2632  * cell != triangulation.end();
2633  * ++cell)
2634  * if (cell->is_locally_owned())
2635  * cell->set_refine_flag();
2636  * }
2637  * else
2638  * {
2639  * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
2640  * KellyErrorEstimator<dim>::estimate(
2641  * dof_handler,
2642  * QGauss<dim - 1>(fe.degree + 2),
2643  * std::map<types::boundary_id, const Function<dim> *>(),
2644  * solution,
2645  * estimated_error_per_cell);
2646  *
2647  * parallel::distributed::GridRefinement ::refine_and_coarsen_fixed_number(
2648  * triangulation, estimated_error_per_cell, 0.3, 0.03);
2649  * }
2650  *
2651  * triangulation.prepare_coarsening_and_refinement();
2652  *
2653  * parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::Vector>
2654  * solution_transfer(dof_handler);
2655  * if (transfer_solution)
2656  * solution_transfer.prepare_for_coarsening_and_refinement(solution);
2657  *
2658  * triangulation.execute_coarsening_and_refinement();
2659  *
2660  * setup_system();
2661  *
2662  * if (transfer_solution)
2663  * {
2664  * TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
2665  * mpi_communicator);
2666  * solution_transfer.interpolate(distributed_solution);
2667  *
2668  * @endcode
2669  *
2670  * enforce constraints to make the interpolated solution conforming on
2671  * the new mesh:
2672  *
2673  * @code
2674  * constraints_hanging_nodes.distribute(distributed_solution);
2675  *
2676  * solution = distributed_solution;
2677  * compute_nonlinear_residual(solution);
2678  * }
2679  * }
2680  *
2681  *
2682  * @endcode
2683  *
2684  *
2685  * <a name="PlasticityContactProblemmove_mesh"></a>
2686  * <h4>PlasticityContactProblem::move_mesh</h4>
2687  *
2688 
2689  *
2690  * The remaining three functions before we get to <code>run()</code>
2691  * have to do with generating output. The following one is an attempt
2692  * at showing the deformed body in its deformed configuration. To this
2693  * end, this function takes a displacement vector field and moves every
2694  * vertex of the (local part) of the mesh by the previously computed
2695  * displacement. We will call this function with the current
2696  * displacement field before we generate graphical output, and we will
2697  * call it again after generating graphical output with the negative
2698  * displacement field to undo the changes to the mesh so made.
2699  *
2700 
2701  *
2702  * The function itself is pretty straightforward. All we have to do
2703  * is keep track which vertices we have already touched, as we
2704  * encounter the same vertices multiple times as we loop over cells.
2705  *
2706  * @code
2707  * template <int dim>
2708  * void PlasticityContactProblem<dim>::move_mesh(
2709  * const TrilinosWrappers::MPI::Vector &displacement) const
2710  * {
2711  * std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
2712  *
2713  * for (const auto &cell : dof_handler.active_cell_iterators())
2714  * if (cell->is_locally_owned())
2715  * for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2716  * if (vertex_touched[cell->vertex_index(v)] == false)
2717  * {
2718  * vertex_touched[cell->vertex_index(v)] = true;
2719  *
2720  * Point<dim> vertex_displacement;
2721  * for (unsigned int d = 0; d < dim; ++d)
2722  * vertex_displacement[d] =
2723  * displacement(cell->vertex_dof_index(v, d));
2724  *
2725  * cell->vertex(v) += vertex_displacement;
2726  * }
2727  * }
2728  *
2729  *
2730  *
2731  * @endcode
2732  *
2733  *
2734  * <a name="PlasticityContactProblemoutput_results"></a>
2735  * <h4>PlasticityContactProblem::output_results</h4>
2736  *
2737 
2738  *
2739  * Next is the function we use to actually generate graphical output. The
2740  * function is a bit tedious, but not actually particularly complicated.
2741  * It moves the mesh at the top (and moves it back at the end), then
2742  * computes the contact forces along the contact surface. We can do
2743  * so (as shown in the accompanying paper) by taking the untreated
2744  * residual vector and identifying which degrees of freedom
2745  * correspond to those with contact by asking whether they have an
2746  * inhomogeneous constraints associated with them. As always, we need
2747  * to be mindful that we can only write into completely distributed
2748  * vectors (i.e., vectors without ghost elements) but that when we
2749  * want to generate output, we need vectors that do indeed have
2750  * ghost entries for all locally relevant degrees of freedom.
2751  *
2752  * @code
2753  * template <int dim>
2754  * void PlasticityContactProblem<dim>::output_results(
2755  * const unsigned int current_refinement_cycle)
2756  * {
2757  * TimerOutput::Scope t(computing_timer, "Graphical output");
2758  *
2759  * pcout << " Writing graphical output... " << std::flush;
2760  *
2761  * move_mesh(solution);
2762  *
2763  * @endcode
2764  *
2765  * Calculation of the contact forces
2766  *
2767  * @code
2768  * TrilinosWrappers::MPI::Vector distributed_lambda(locally_owned_dofs,
2769  * mpi_communicator);
2770  * const unsigned int start_res = (newton_rhs_uncondensed.local_range().first),
2771  * end_res = (newton_rhs_uncondensed.local_range().second);
2772  * for (unsigned int n = start_res; n < end_res; ++n)
2773  * if (all_constraints.is_inhomogeneously_constrained(n))
2774  * distributed_lambda(n) =
2775  * newton_rhs_uncondensed(n) / diag_mass_matrix_vector(n);
2776  * distributed_lambda.compress(VectorOperation::insert);
2777  * constraints_hanging_nodes.distribute(distributed_lambda);
2778  *
2779  * TrilinosWrappers::MPI::Vector lambda(locally_relevant_dofs,
2780  * mpi_communicator);
2781  * lambda = distributed_lambda;
2782  *
2783  * TrilinosWrappers::MPI::Vector distributed_active_set_vector(
2784  * locally_owned_dofs, mpi_communicator);
2785  * distributed_active_set_vector = 0.;
2786  * for (const auto index : active_set)
2787  * distributed_active_set_vector[index] = 1.;
2788  * distributed_lambda.compress(VectorOperation::insert);
2789  *
2790  * TrilinosWrappers::MPI::Vector active_set_vector(locally_relevant_dofs,
2791  * mpi_communicator);
2792  * active_set_vector = distributed_active_set_vector;
2793  *
2794  * DataOut<dim> data_out;
2795  *
2796  * data_out.attach_dof_handler(dof_handler);
2797  *
2798  * const std::vector<DataComponentInterpretation::DataComponentInterpretation>
2799  * data_component_interpretation(
2800  * dim, DataComponentInterpretation::component_is_part_of_vector);
2801  * data_out.add_data_vector(solution,
2802  * std::vector<std::string>(dim, "displacement"),
2803  * DataOut<dim>::type_dof_data,
2804  * data_component_interpretation);
2805  * data_out.add_data_vector(lambda,
2806  * std::vector<std::string>(dim, "contact_force"),
2807  * DataOut<dim>::type_dof_data,
2808  * data_component_interpretation);
2809  * data_out.add_data_vector(active_set_vector,
2810  * std::vector<std::string>(dim, "active_set"),
2811  * DataOut<dim>::type_dof_data,
2812  * data_component_interpretation);
2813  *
2814  * Vector<float> subdomain(triangulation.n_active_cells());
2815  * for (unsigned int i = 0; i < subdomain.size(); ++i)
2816  * subdomain(i) = triangulation.locally_owned_subdomain();
2817  * data_out.add_data_vector(subdomain, "subdomain");
2818  *
2819  * data_out.add_data_vector(fraction_of_plastic_q_points_per_cell,
2820  * "fraction_of_plastic_q_points");
2821  *
2822  * data_out.build_patches();
2823  *
2824  * @endcode
2825  *
2826  * In the remainder of the function, we generate one VTU file on
2827  * every processor, indexed by the subdomain id of this processor.
2828  * On the first processor, we then also create a <code>.pvtu</code>
2829  * file that indexes <i>all</i> of the VTU files so that the entire
2830  * set of output files can be read at once. These <code>.pvtu</code>
2831  * are used by Paraview to describe an entire parallel computation's
2832  * output files. We then do the same again for the competitor of
2833  * Paraview, the VisIt visualization program, by creating a matching
2834  * <code>.visit</code> file.
2835  *
2836  * @code
2837  * const std::string master_name = data_out.write_vtu_with_pvtu_record(
2838  * output_dir, "solution", current_refinement_cycle, mpi_communicator, 2);
2839  * pcout << master_name << std::endl;
2840  *
2841  * TrilinosWrappers::MPI::Vector tmp(solution);
2842  * tmp *= -1;
2843  * move_mesh(tmp);
2844  * }
2845  *
2846  *
2847  * @endcode
2848  *
2849  *
2850  * <a name="PlasticityContactProblemoutput_contact_force"></a>
2851  * <h4>PlasticityContactProblem::output_contact_force</h4>
2852  *
2853 
2854  *
2855  * This last auxiliary function computes the contact force by
2856  * calculating an integral over the contact pressure in z-direction
2857  * over the contact area. For this purpose we set the contact
2858  * pressure lambda to 0 for all inactive dofs (whether a degree
2859  * of freedom is part of the contact is determined just as
2860  * we did in the previous function). For all
2861  * active dofs, lambda contains the quotient of the nonlinear
2862  * residual (newton_rhs_uncondensed) and corresponding diagonal entry
2863  * of the mass matrix (diag_mass_matrix_vector). Because it is
2864  * not unlikely that hanging nodes show up in the contact area
2865  * it is important to apply constraints_hanging_nodes.distribute
2866  * to the distributed_lambda vector.
2867  *
2868  * @code
2869  * template <int dim>
2870  * void PlasticityContactProblem<dim>::output_contact_force() const
2871  * {
2872  * TrilinosWrappers::MPI::Vector distributed_lambda(locally_owned_dofs,
2873  * mpi_communicator);
2874  * const unsigned int start_res = (newton_rhs_uncondensed.local_range().first),
2875  * end_res = (newton_rhs_uncondensed.local_range().second);
2876  * for (unsigned int n = start_res; n < end_res; ++n)
2877  * if (all_constraints.is_inhomogeneously_constrained(n))
2878  * distributed_lambda(n) =
2879  * newton_rhs_uncondensed(n) / diag_mass_matrix_vector(n);
2880  * else
2881  * distributed_lambda(n) = 0;
2882  * distributed_lambda.compress(VectorOperation::insert);
2883  * constraints_hanging_nodes.distribute(distributed_lambda);
2884  *
2885  * TrilinosWrappers::MPI::Vector lambda(locally_relevant_dofs,
2886  * mpi_communicator);
2887  * lambda = distributed_lambda;
2888  *
2889  * double contact_force = 0.0;
2890  *
2891  * QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
2892  * FEFaceValues<dim> fe_values_face(fe,
2893  * face_quadrature_formula,
2895  *
2896  * const unsigned int n_face_q_points = face_quadrature_formula.size();
2897  *
2898  * const FEValuesExtractors::Vector displacement(0);
2899  *
2900  * for (const auto &cell : dof_handler.active_cell_iterators())
2901  * if (cell->is_locally_owned())
2902  * for (const auto &face : cell->face_iterators())
2903  * if (face->at_boundary() && face->boundary_id() == 1)
2904  * {
2905  * fe_values_face.reinit(cell, face);
2906  *
2907  * std::vector<Tensor<1, dim>> lambda_values(n_face_q_points);
2908  * fe_values_face[displacement].get_function_values(lambda,
2909  * lambda_values);
2910  *
2911  * for (unsigned int q_point = 0; q_point < n_face_q_points;
2912  * ++q_point)
2913  * contact_force +=
2914  * lambda_values[q_point][2] * fe_values_face.JxW(q_point);
2915  * }
2916  * contact_force = Utilities::MPI::sum(contact_force, MPI_COMM_WORLD);
2917  *
2918  * pcout << "Contact force = " << contact_force << std::endl;
2919  * }
2920  *
2921  *
2922  * @endcode
2923  *
2924  *
2925  * <a name="PlasticityContactProblemrun"></a>
2927  *
2928 
2929  *
2930  * As in all other tutorial programs, the <code>run()</code> function contains
2931  * the overall logic. There is not very much to it here: in essence, it
2932  * performs the loops over all mesh refinement cycles, and within each, hands
2933  * things over to the Newton solver in <code>solve_newton()</code> on the
2934  * current mesh and calls the function that creates graphical output for
2935  * the so-computed solution. It then outputs some statistics concerning both
2936  * run times and memory consumption that has been collected over the course of
2937  * computations on this mesh.
2938  *
2939  * @code
2940  * template <int dim>
2942  * {
2943  * computing_timer.reset();
2944  * for (; current_refinement_cycle < n_refinement_cycles;
2945  * ++current_refinement_cycle)
2946  * {
2947  * {
2948  * TimerOutput::Scope t(computing_timer, "Setup");
2949  *
2950  * pcout << std::endl;
2951  * pcout << "Cycle " << current_refinement_cycle << ':' << std::endl;
2952  *
2953  * if (current_refinement_cycle == 0)
2954  * {
2955  * make_grid();
2956  * setup_system();
2957  * }
2958  * else
2959  * {
2960  * TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
2961  * refine_grid();
2962  * }
2963  * }
2964  *
2965  * solve_newton();
2966  *
2967  * output_results(current_refinement_cycle);
2968  *
2969  * computing_timer.print_summary();
2970  * computing_timer.reset();
2971  *
2974  * pcout << "Peak virtual memory used, resident in kB: " << stats.VmSize
2975  * << " " << stats.VmRSS << std::endl;
2976  *
2977  * if (base_mesh == "box")
2978  * output_contact_force();
2979  * }
2980  * }
2981  * } // namespace Step42
2982  *
2983  * @endcode
2984  *
2985  *
2986  * <a name="Thecodemaincodefunction"></a>
2987  * <h3>The <code>main</code> function</h3>
2988  *
2989 
2990  *
2991  * There really isn't much to the <code>main()</code> function. It looks
2992  * like they always do:
2993  *
2994  * @code
2995  * int main(int argc, char *argv[])
2996  * {
2997  * using namespace dealii;
2998  * using namespace Step42;
2999  *
3000  * try
3001  * {
3002  * ParameterHandler prm;
3003  * PlasticityContactProblem<3>::declare_parameters(prm);
3004  * if (argc != 2)
3005  * {
3006  * std::cerr << "*** Call this program as <./step-42 input.prm>"
3007  * << std::endl;
3008  * return 1;
3009  * }
3010  *
3011  * prm.parse_input(argv[1]);
3012  * Utilities::MPI::MPI_InitFinalize mpi_initialization(
3013  * argc, argv, numbers::invalid_unsigned_int);
3014  * {
3015  * PlasticityContactProblem<3> problem(prm);
3016  * problem.run();
3017  * }
3018  * }
3019  * catch (std::exception &exc)
3020  * {
3021  * std::cerr << std::endl
3022  * << std::endl
3023  * << "----------------------------------------------------"
3024  * << std::endl;
3025  * std::cerr << "Exception on processing: " << std::endl
3026  * << exc.what() << std::endl
3027  * << "Aborting!" << std::endl
3028  * << "----------------------------------------------------"
3029  * << std::endl;
3030  *
3031  * return 1;
3032  * }
3033  * catch (...)
3034  * {
3035  * std::cerr << std::endl
3036  * << std::endl
3037  * << "----------------------------------------------------"
3038  * << std::endl;
3039  * std::cerr << "Unknown exception!" << std::endl
3040  * << "Aborting!" << std::endl
3041  * << "----------------------------------------------------"
3042  * << std::endl;
3043  * return 1;
3044  * }
3045  *
3046  * return 0;
3047  * }
3048  * @endcode
3049 <a name="Results"></a><h1>Results</h1>
3050 
3051 
3052 The directory that contains this program also contains a number of input
3053 parameter files that can be used to create various different
3054 simulations. For example, running the program with the
3055 <code>p1_adaptive.prm</code> parameter file (using a ball as obstacle and the
3056 box as domain) on 16 cores produces output like this:
3057 @code
3058  Using output directory 'p1adaptive/'
3059  FE degree 1
3060  transfer solution false
3061 
3062 Cycle 0:
3063  Number of active cells: 512
3064  Number of degrees of freedom: 2187
3065 
3066  Newton iteration 1
3067  Updating active set...
3068  Size of active set: 1
3069  Assembling system...
3070  Solving system...
3071  Error: 173.076 -> 1.64265e-06 in 7 Bicgstab iterations.
3072  Accepting Newton solution with residual: 1.64265e-06
3073 
3074  Newton iteration 2
3075  Updating active set...
3076  Size of active set: 1
3077  Assembling system...
3078  Solving system...
3079  Error: 57.3622 -> 3.23721e-07 in 8 Bicgstab iterations.
3080  Accepting Newton solution with residual: 24.9028
3081  Active set did not change!
3082 
3083  Newton iteration 3
3084  Updating active set...
3085  Size of active set: 1
3086  Assembling system...
3087  Solving system...
3088  Error: 24.9028 -> 9.94326e-08 in 7 Bicgstab iterations.
3089  Residual of the non-contact part of the system: 1.63333
3090  with a damping parameter alpha = 1
3091  Active set did not change!
3092 
3093 ...
3094 
3095  Newton iteration 6
3096  Updating active set...
3097  Size of active set: 1
3098  Assembling system...
3099  Solving system...
3100  Error: 1.43188e-07 -> 3.56218e-16 in 8 Bicgstab iterations.
3101  Residual of the non-contact part of the system: 4.298e-14
3102  with a damping parameter alpha = 1
3103  Active set did not change!
3104  Writing graphical output... p1_adaptive/solution-00.pvtu
3105 
3106 
3107 +---------------------------------------------+------------+------------+
3108 | Total wallclock time elapsed since start | 1.13s | |
3109 | | | |
3110 | Section | no. calls | wall time | % of total |
3111 +---------------------------------+-----------+------------+------------+
3112 | Assembling | 6 | 0.463s | 41% |
3113 | Graphical output | 1 | 0.0257s | 2.3% |
3114 | Residual and lambda | 4 | 0.0754s | 6.7% |
3115 | Setup | 1 | 0.227s | 20% |
3116 | Setup: constraints | 1 | 0.0347s | 3.1% |
3117 | Setup: distribute DoFs | 1 | 0.0441s | 3.9% |
3118 | Setup: matrix | 1 | 0.0119s | 1.1% |
3119 | Setup: vectors | 1 | 0.00155s | 0.14% |
3120 | Solve | 6 | 0.246s | 22% |
3121 | Solve: iterate | 6 | 0.0631s | 5.6% |
3122 | Solve: setup preconditioner | 6 | 0.167s | 15% |
3123 | update active set | 6 | 0.0401s | 3.6% |
3124 +---------------------------------+-----------+------------+------------+
3125 
3126 Peak virtual memory used, resident in kB: 541884 77464
3127 Contact force = 37.3058
3128 
3129 ...
3130 
3131 Cycle 3:
3132  Number of active cells: 14652
3133  Number of degrees of freedom: 52497
3134 
3135  Newton iteration 1
3136  Updating active set...
3137  Size of active set: 145
3138  Assembling system...
3139  Solving system...
3140  Error: 296.309 -> 2.72484e-06 in 10 Bicgstab iterations.
3141  Accepting Newton solution with residual: 2.72484e-06
3142 
3143 ...
3144 
3145  Newton iteration 10
3146  Updating active set...
3147  Size of active set: 145
3148  Assembling system...
3149  Solving system...
3150  Error: 2.71541e-07 -> 1.5428e-15 in 27 Bicgstab iterations.
3151  Residual of the non-contact part of the system: 1.89261e-13
3152  with a damping parameter alpha = 1
3153  Active set did not change!
3154  Writing graphical output... p1_adaptive/solution-03.pvtu
3155 
3156 
3157 +---------------------------------------------+------------+------------+
3158 | Total wallclock time elapsed since start | 38.4s | |
3159 | | | |
3160 | Section | no. calls | wall time | % of total |
3161 +---------------------------------+-----------+------------+------------+
3162 | Assembling | 10 | 22.5s | 58% |
3163 | Graphical output | 1 | 0.327s | 0.85% |
3164 | Residual and lambda | 9 | 3.75s | 9.8% |
3165 | Setup | 1 | 4.83s | 13% |
3166 | Setup: constraints | 1 | 0.578s | 1.5% |
3167 | Setup: distribute DoFs | 1 | 0.71s | 1.8% |
3168 | Setup: matrix | 1 | 0.111s | 0.29% |
3169 | Setup: refine mesh | 1 | 4.83s | 13% |
3170 | Setup: vectors | 1 | 0.00548s | 0.014% |
3171 | Solve | 10 | 5.49s | 14% |
3172 | Solve: iterate | 10 | 3.5s | 9.1% |
3173 | Solve: setup preconditioner | 10 | 1.84s | 4.8% |
3174 | update active set | 10 | 0.662s | 1.7% |
3175 +---------------------------------+-----------+------------+------------+
3176 
3177 Peak virtual memory used, resident in kB: 566052 105788
3178 Contact force = 56.794
3179 
3180 ...
3181 @endcode
3182 
3183 The tables at the end of each cycle show information about computing time
3184 (these numbers are of course specific to the machine on which this output
3185 was produced)
3186 and the number of calls of different parts of the program like assembly or
3187 calculating the residual, for the most recent mesh refinement cycle. Some of
3188 the numbers above can be improved by transferring the solution from one mesh to
3189 the next, an option we have not exercised here. Of course, you can also make
3190 the program run faster, especially on the later refinement cycles, by just
3191 using more processors: the accompanying paper shows good scaling to at least
3192 1000 cores.
3193 
3194 In a typical run, you can observe that for every refinement step, the active
3195 set - the contact points - are iterated out at first. After that the Newton
3196 method has only to resolve the plasticity. For the finer meshes,
3197 quadratic convergence can be observed for the last 4 or 5 Newton iterations.
3198 
3199 We will not discuss here in all detail what happens with each of the input
3200 files. Rather, let us just show pictures of the solution (the left half of the
3201 domain is omitted if cells have zero quadrature points at which the plastic
3202 inequality is active):
3203 
3204 <table align="center">
3205  <tr>
3206  <td>
3207  <img src="https://www.dealii.org/images/steps/developer/step-42.CellConstitutionColorbar.png">
3208  </td>
3209  <td>
3210  <img src="https://www.dealii.org/images/steps/developer/step-42.CellConstitutionBall2.png" alt="" width="70%">
3211  </td>
3212  <td valign="top">
3213  &nbsp;
3214  </td>
3215  <td>
3216  <img src="https://www.dealii.org/images/steps/developer/step-42.CellConstitutionLi2.png" alt="" alt="" width="70%">
3217  </td>
3218  </tr>
3219 </table>
3220 
3221 The picture shows the adaptive refinement and as well how much a cell is
3222 plastified during the contact with the ball. Remember that we consider the
3223 norm of the deviator part of the stress in each quadrature point to
3224 see if there is elastic or plastic behavior.
3225 The blue
3226 color means that this cell contains only elastic quadrature points in
3227 contrast to the red cells in which all quadrature points are plastified.
3228 In the middle of the top surface -
3229 where the mesh is finest - a very close look shows the dimple caused by the
3230 obstacle. This is the result of the <code>move_mesh()</code>
3231 function. However, because the indentation of the obstacles we consider here
3232 is so small, it is hard to discern this effect; one could play with displacing
3233 vertices of the mesh by a multiple of the computed displacement.
3234 
3235 Further discussion of results that can be obtained using this program is
3236 provided in the publication mentioned at the very top of this page.
3237 
3238 
3239 <a name="extensions"></a>
3240 <a name="Possibilitiesforextensions"></a><h1>Possibilities for extensions</h1>
3241 
3242 
3243 There are, as always, multiple possibilities for extending this program. From
3244 an algorithmic perspective, this program goes about as far as one can at the
3245 time of writing, using the best available algorithms for the contact
3246 inequality, the plastic nonlinearity, and the linear solvers. However, there
3247 are things one would like to do with this program as far as more realistic
3248 situations are concerned:
3249 <ul>
3250 <li> Extend the program from a static to a quasi-static situation, perhaps by
3251 choosing a backward-Euler-scheme for the time discretization. Some theoretical
3252 results can be found in the PhD thesis by Jörg Frohne, <i>FEM-Simulation
3253 der Umformtechnik metallischer Oberfl&auml;chen im Mikrokosmos</i>, University
3254 of Siegen, Germany, 2011.
3255 
3256 <li> It would also be an interesting advance to consider a contact problem
3257 with friction. In almost every mechanical process friction has a big
3258 influence. To model this situation, we have to take into account tangential
3259 stresses at the contact surface. Friction also adds another inequality to
3260 our problem since body and obstacle will typically stick together as long as
3261 the tangential stress does not exceed a certain limit, beyond which the two
3262 bodies slide past each other.
3263 
3264 <li> If we already simulate a frictional contact, the next step to consider
3265 is heat generation over the contact zone. The heat that is
3266 caused by friction between two bodies raises the temperature in the
3267 deformable body and entails an change of some material parameters.
3268 
3269 <li> It might be of interest to implement more accurate, problem-adapted error
3270 estimators for contact as well as for the plasticity.
3271 </ul>
3272  *
3273  *
3274 <a name="PlainProg"></a>
3275 <h1> The plain program</h1>
3276 @include "step-42.cc"
3277 */
LAPACKSupport::diagonal
@ diagonal
Matrix is diagonal.
Definition: lapack_support.h:121
TrilinosWrappers::MPI::Vector
Definition: trilinos_vector.h:400
LocalIntegrators::L2::mass_matrix
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: l2.h:63
TimerOutput::Scope
Definition: timer.h:554
TrilinosWrappers::SparseMatrix
Definition: trilinos_sparse_matrix.h:515
VectorOperation::insert
@ insert
Definition: vector_operation.h:49
Differentiation::SD::OptimizerType::lambda
@ lambda
Utilities::System::get_memory_stats
void get_memory_stats(MemoryStats &stats)
Definition: utilities.cc:971
FEValuesExtractors::Scalar
Definition: fe_values_extractors.h:95
ComponentMask
Definition: component_mask.h:83
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
second
Point< 2 > second
Definition: grid_out.cc:4353
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
TrilinosWrappers::SparsityPattern
Definition: trilinos_sparsity_pattern.h:279
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: dof_tools_sparsity.cc:63
VectorTools::interpolate_boundary_values
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
FEValuesExtractors::Vector
Definition: fe_values_extractors.h:150
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
std_cxx17::apply
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
Utilities::System::MemoryStats::VmSize
unsigned long int VmSize
Definition: utilities.h:861
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, AffineConstraints< number > &constraints)
Definition: dof_tools_constraints.cc:1725
QGaussLobatto
Definition: quadrature_lib.h:76
QGauss
Definition: quadrature_lib.h:40
value
static const bool value
Definition: dof_tools_constraints.cc:433
AffineConstraints
Definition: affine_constraints.h:180
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
DoFTools::extract_locally_relevant_dofs
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1173
std::pow
inline ::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &x, const Number p)
Definition: vectorization.h:5428
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
FETools::interpolate
void interpolate(const DoFHandlerType1< dim, spacedim > &dof1, const InVector &u1, const DoFHandlerType2< dim, spacedim > &dof2, OutVector &u2)
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
FEFaceValues< dim >
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
first
Point< 2 > first
Definition: grid_out.cc:4352
Utilities::System::MemoryStats::VmRSS
unsigned long int VmRSS
Definition: utilities.h:872
Utilities::System::MemoryStats
Definition: utilities.h:851