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