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