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