1621 * dof_handler.distribute_dofs(fe);
1623 * locally_owned_dofs = dof_handler.locally_owned_dofs();
1624 * locally_relevant_dofs =
1631 * constraints_hanging_nodes.reinit(locally_owned_dofs,
1632 * locally_relevant_dofs);
1634 * constraints_hanging_nodes);
1635 * constraints_hanging_nodes.close();
1637 * pcout <<
" Number of active cells: "
1639 * <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1642 * compute_dirichlet_constraints();
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(
1655 * active_set.clear();
1656 * active_set.set_size(dof_handler.n_dofs());
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
1671 * mpi_communicator);
1675 * constraints_dirichlet_and_hanging_nodes,
1678 * mpi_communicator));
1680 * newton_matrix.reinit(sp);
1685 * assemble_mass_matrix_diagonal(mass_matrix);
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);
1701 * <a name=
"step_42-PlasticityContactProblemcompute_dirichlet_constraints"></a>
1702 * <h4>PlasticityContactProblem::compute_dirichlet_constraints</h4>
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.
1713 * As laid out in the introduction, we need to distinguish between two
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.
1723 *
template <
int dim>
1724 *
void PlasticityContactProblem<dim>::compute_dirichlet_constraints()
1726 * constraints_dirichlet_and_hanging_nodes.reinit(locally_owned_dofs,
1727 * locally_relevant_dofs);
1728 * constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);
1730 *
if (base_mesh ==
"box")
1740 * EquationData::BoundaryValues<dim>(),
1741 * constraints_dirichlet_and_hanging_nodes,
1747 * solution (
this is a bit mask, so apply
1756 * EquationData::BoundaryValues<dim>(),
1757 * constraints_dirichlet_and_hanging_nodes,
1758 * (fe.component_mask(x_displacement) |
1759 * fe.component_mask(y_displacement)));
1765 * EquationData::BoundaryValues<dim>(),
1766 * constraints_dirichlet_and_hanging_nodes,
1769 * constraints_dirichlet_and_hanging_nodes.close();
1777 * <a name=
"step_42-PlasticityContactProblemassemble_mass_matrix_diagonal"></a>
1778 * <h4>PlasticityContactProblem::assemble_mass_matrix_diagonal</h4>
1782 * The next helper function computes the (diagonal) @ref GlossMassMatrix
"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
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.
1801 * template <int dim>
1802 * void PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal(
1803 * TrilinosWrappers::SparseMatrix &mass_matrix)
1805 * const QGaussLobatto<dim - 1> face_quadrature_formula(fe.degree + 1);
1807 * FEFaceValues<dim> fe_values_face(fe,
1808 * face_quadrature_formula,
1809 * update_values | update_JxW_values);
1811 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1812 * const unsigned int n_face_q_points = face_quadrature_formula.size();
1814 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1815 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1817 * const FEValuesExtractors::Vector displacement(0);
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)
1824 * fe_values_face.reinit(cell, face);
1827 * for (unsigned int q_point = 0; q_point < n_face_q_points;
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));
1835 * cell->get_dof_indices(local_dof_indices);
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));
1842 * mass_matrix.compress(VectorOperation::add);
1849 * <a name="step_42-PlasticityContactProblemupdate_solution_and_constraints"></a>
1850 * <h4>PlasticityContactProblem::update_solution_and_constraints</h4>
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.
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
1873 * template <int dim>
1874 * void PlasticityContactProblem<dim>::update_solution_and_constraints()
1876 * std::vector<bool> dof_touched(dof_handler.n_dofs(), false);
1878 * TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
1879 * mpi_communicator);
1880 * distributed_solution = solution;
1882 * TrilinosWrappers::MPI::Vector lambda(locally_relevant_dofs,
1883 * mpi_communicator);
1884 * lambda = newton_rhs_uncondensed;
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;
1891 * all_constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
1892 * active_set.clear();
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:
1912 * const Quadrature<dim - 1> face_quadrature(
1913 * fe.get_unit_face_support_points());
1914 * FEFaceValues<dim> fe_values_face(fe,
1916 * update_quadrature_points);
1918 * const unsigned int dofs_per_face = fe.n_dofs_per_face();
1919 * const unsigned int n_face_q_points = face_quadrature.size();
1921 * std::vector<types::global_dof_index> dof_indices(dofs_per_face);
1923 * for (const auto &cell : dof_handler.active_cell_iterators())
1924 * if (!cell->is_artificial())
1925 * for (const auto &face : cell->face_iterators())
1926 * if (face->at_boundary() && face->boundary_id() == 1)
1928 * fe_values_face.reinit(cell, face);
1929 * face->get_dof_indices(dof_indices);
1931 * for (unsigned int q_point = 0; q_point < n_face_q_points;
1936 * At each quadrature point (i.e., at each support point of a
1937 * degree of freedom located on the contact boundary), we then
1938 * ask whether it is part of the z-displacement degrees of
1939 * freedom and if we haven't encountered
this degree of
1940 * freedom yet (which can happen
for those on the edges
1941 * between faces), we need to evaluate the gap between the
1942 * deformed
object and the obstacle. If the active set
1943 * condition is
true, then we add a constraint to the
1945 * to satisfy, set the solution vector
's corresponding element
1946 * to the correct value, and add the index to the IndexSet
1947 * object that stores which degree of freedom is part of the
1951 * const unsigned int component =
1952 * fe.face_system_to_component_index(q_point).first;
1954 * const unsigned int index_z = dof_indices[q_point];
1956 * if ((component == 2) && (dof_touched[index_z] == false))
1958 * dof_touched[index_z] = true;
1960 * const Point<dim> this_support_point =
1961 * fe_values_face.quadrature_point(q_point);
1963 * const double obstacle_value =
1964 * obstacle->value(this_support_point, 2);
1965 * const double solution_here = solution(index_z);
1966 * const double undeformed_gap =
1967 * obstacle_value - this_support_point[2];
1969 * const double c = 100.0 * e_modulus;
1970 * if ((lambda(index_z) /
1971 * diag_mass_matrix_vector_relevant(index_z) +
1972 * c * (solution_here - undeformed_gap) >
1974 * !constraints_hanging_nodes.is_constrained(index_z))
1976 * all_constraints.add_constraint(index_z,
1979 * distributed_solution(index_z) = undeformed_gap;
1981 * active_set.add_index(index_z);
1989 * At the end of this function, we exchange data between processors updating
1990 * those ghost elements in the <code>solution</code> variable that have been
1991 * written by other processors. We then merge the Dirichlet constraints and
1992 * those from hanging nodes into the AffineConstraints object that already
1993 * contains the active set. We finish the function by outputting the total
1994 * number of actively constrained degrees of freedom for which we sum over
1995 * the number of actively constrained degrees of freedom owned by each
1996 * of the processors. This number of locally owned constrained degrees of
1997 * freedom is of course the number of elements of the intersection of the
1998 * active set and the set of locally owned degrees of freedom, which
1999 * we can get by using <code>operator&</code> on two IndexSets:
2002 * distributed_solution.compress(VectorOperation::insert);
2003 * solution = distributed_solution;
2005 * all_constraints.close();
2006 * all_constraints.merge(constraints_dirichlet_and_hanging_nodes);
2008 * pcout << " Size of active set: "
2009 * << Utilities::MPI::sum((active_set & locally_owned_dofs).n_elements(),
2018 * <a name="step_42-PlasticityContactProblemassemble_newton_system"></a>
2019 * <h4>PlasticityContactProblem::assemble_newton_system</h4>
2023 * Given the complexity of the problem, it may come as a bit of a surprise
2024 * that assembling the linear system we have to solve in each Newton iteration
2025 * is actually fairly straightforward. The following function builds the
2026 * Newton right hand side and Newton matrix. It looks fairly innocent because
2027 * the heavy lifting happens in the call to
2028 * <code>ConstitutiveLaw::get_linearized_stress_strain_tensors()</code> and in
2029 * particular in AffineConstraints::distribute_local_to_global(), using the
2030 * constraints we have previously computed.
2033 * template <int dim>
2034 * void PlasticityContactProblem<dim>::assemble_newton_system(
2035 * const TrilinosWrappers::MPI::Vector &linearization_point)
2037 * TimerOutput::Scope t(computing_timer, "Assembling");
2039 * const QGauss<dim> quadrature_formula(fe.degree + 1);
2040 * const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
2042 * FEValues<dim> fe_values(fe,
2043 * quadrature_formula,
2044 * update_values | update_gradients |
2045 * update_JxW_values);
2047 * FEFaceValues<dim> fe_values_face(fe,
2048 * face_quadrature_formula,
2049 * update_values | update_quadrature_points |
2050 * update_JxW_values);
2052 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
2053 * const unsigned int n_q_points = quadrature_formula.size();
2054 * const unsigned int n_face_q_points = face_quadrature_formula.size();
2056 * const EquationData::BoundaryForce<dim> boundary_force;
2057 * std::vector<Vector<double>> boundary_force_values(n_face_q_points,
2058 * Vector<double>(dim));
2060 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
2061 * Vector<double> cell_rhs(dofs_per_cell);
2063 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2065 * const FEValuesExtractors::Vector displacement(0);
2067 * for (const auto &cell : dof_handler.active_cell_iterators())
2068 * if (cell->is_locally_owned())
2070 * fe_values.reinit(cell);
2074 * std::vector<SymmetricTensor<2, dim>> strain_tensor(n_q_points);
2075 * fe_values[displacement].get_function_symmetric_gradients(
2076 * linearization_point, strain_tensor);
2078 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2080 * SymmetricTensor<4, dim> stress_strain_tensor_linearized;
2081 * SymmetricTensor<4, dim> stress_strain_tensor;
2082 * constitutive_law.get_linearized_stress_strain_tensors(
2083 * strain_tensor[q_point],
2084 * stress_strain_tensor_linearized,
2085 * stress_strain_tensor);
2087 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2091 * Having computed the stress-strain tensor and its
2092 * linearization, we can now put together the parts of the
2093 * matrix and right hand side. In both, we need the linearized
2094 * stress-strain tensor times the symmetric gradient of
2095 * @f$\varphi_i@f$, i.e. the term @f$I_\Pi\varepsilon(\varphi_i)@f$,
2096 * so we introduce an abbreviation of this term. Recall that
2097 * the matrix corresponds to the bilinear form
2098 * @f$A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))@f$
2099 * in the notation of the accompanying publication, whereas
2100 * the right hand side is @f$F_i=([I_\Pi-P_\Pi
2101 * C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))@f$ where @f$u@f$
2102 * is the current linearization points (typically the last
2103 * solution). This might suggest that the right hand side will
2104 * be zero if the material is completely elastic (where
2105 * @f$I_\Pi=P_\Pi@f$) but this ignores the fact that the right
2106 * hand side will also contain contributions from
2107 * non-homogeneous constraints due to the contact.
2111 * The code block that follows this adds contributions that
2112 * are due to boundary forces, should there be any.
2115 * const SymmetricTensor<2, dim> stress_phi_i =
2116 * stress_strain_tensor_linearized *
2117 * fe_values[displacement].symmetric_gradient(i, q_point);
2119 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
2120 * cell_matrix(i, j) +=
2122 * fe_values[displacement].symmetric_gradient(j, q_point) *
2123 * fe_values.JxW(q_point));
2127 * stress_strain_tensor *
2128 * fe_values[displacement].symmetric_gradient(i,
2130 * strain_tensor[q_point] * fe_values.JxW(q_point));
2134 * for (const auto &face : cell->face_iterators())
2135 * if (face->at_boundary() && face->boundary_id() == 1)
2137 * fe_values_face.reinit(cell, face);
2139 * boundary_force.vector_value_list(
2140 * fe_values_face.get_quadrature_points(),
2141 * boundary_force_values);
2143 * for (unsigned int q_point = 0; q_point < n_face_q_points;
2146 * Tensor<1, dim> rhs_values;
2147 * rhs_values[2] = boundary_force_values[q_point][2];
2148 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2150 * (fe_values_face[displacement].value(i, q_point) *
2151 * rhs_values * fe_values_face.JxW(q_point));
2155 * cell->get_dof_indices(local_dof_indices);
2156 * all_constraints.distribute_local_to_global(cell_matrix,
2158 * local_dof_indices,
2164 * newton_matrix.compress(VectorOperation::add);
2165 * newton_rhs.compress(VectorOperation::add);
2173 * <a name="step_42-PlasticityContactProblemcompute_nonlinear_residual"></a>
2174 * <h4>PlasticityContactProblem::compute_nonlinear_residual</h4>
2178 * The following function computes the nonlinear residual of the equation
2179 * given the current solution (or any other linearization point). This
2180 * is needed in the linear search algorithm where we need to try various
2181 * linear combinations of previous and current (trial) solution to
2182 * compute the (real, globalized) solution of the current Newton step.
2186 * That said, in a slight abuse of the name of the function, it actually
2187 * does significantly more. For example, it also computes the vector
2188 * that corresponds to the Newton residual but without eliminating
2189 * constrained degrees of freedom. We need this vector to compute contact
2190 * forces and, ultimately, to compute the next active set. Likewise, by
2191 * keeping track of how many quadrature points we encounter on each cell
2192 * that show plastic yielding, we also compute the
2193 * <code>fraction_of_plastic_q_points_per_cell</code> vector that we
2194 * can later output to visualize the plastic zone. In both of these cases,
2195 * the results are not necessary as part of the line search, and so we may
2196 * be wasting a small amount of time computing them. At the same time, this
2197 * information appears as a natural by-product of what we need to do here
2198 * anyway, and we want to collect it once at the end of each Newton
2199 * step, so we may as well do it here.
2203 * The actual implementation of this function should be rather obvious:
2206 * template <int dim>
2207 * void PlasticityContactProblem<dim>::compute_nonlinear_residual(
2208 * const TrilinosWrappers::MPI::Vector &linearization_point)
2210 * const QGauss<dim> quadrature_formula(fe.degree + 1);
2211 * const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
2213 * FEValues<dim> fe_values(fe,
2214 * quadrature_formula,
2215 * update_values | update_gradients |
2216 * update_JxW_values);
2218 * FEFaceValues<dim> fe_values_face(fe,
2219 * face_quadrature_formula,
2220 * update_values | update_quadrature_points |
2221 * update_JxW_values);
2223 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
2224 * const unsigned int n_q_points = quadrature_formula.size();
2225 * const unsigned int n_face_q_points = face_quadrature_formula.size();
2227 * const EquationData::BoundaryForce<dim> boundary_force;
2228 * std::vector<Vector<double>> boundary_force_values(n_face_q_points,
2229 * Vector<double>(dim));
2231 * Vector<double> cell_rhs(dofs_per_cell);
2233 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2235 * const FEValuesExtractors::Vector displacement(0);
2238 * newton_rhs_uncondensed = 0;
2240 * fraction_of_plastic_q_points_per_cell = 0;
2242 * for (const auto &cell : dof_handler.active_cell_iterators())
2243 * if (cell->is_locally_owned())
2245 * fe_values.reinit(cell);
2248 * std::vector<SymmetricTensor<2, dim>> strain_tensors(n_q_points);
2249 * fe_values[displacement].get_function_symmetric_gradients(
2250 * linearization_point, strain_tensors);
2252 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2254 * SymmetricTensor<4, dim> stress_strain_tensor;
2255 * const bool q_point_is_plastic =
2256 * constitutive_law.get_stress_strain_tensor(
2257 * strain_tensors[q_point], stress_strain_tensor);
2258 * if (q_point_is_plastic)
2259 * ++fraction_of_plastic_q_points_per_cell(
2260 * cell->active_cell_index());
2262 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2265 * (strain_tensors[q_point] * stress_strain_tensor *
2266 * fe_values[displacement].symmetric_gradient(i, q_point) *
2267 * fe_values.JxW(q_point));
2269 * Tensor<1, dim> rhs_values;
2271 * cell_rhs(i) += (fe_values[displacement].value(i, q_point) *
2272 * rhs_values * fe_values.JxW(q_point));
2276 * for (const auto &face : cell->face_iterators())
2277 * if (face->at_boundary() && face->boundary_id() == 1)
2279 * fe_values_face.reinit(cell, face);
2281 * boundary_force.vector_value_list(
2282 * fe_values_face.get_quadrature_points(),
2283 * boundary_force_values);
2285 * for (unsigned int q_point = 0; q_point < n_face_q_points;
2288 * Tensor<1, dim> rhs_values;
2289 * rhs_values[2] = boundary_force_values[q_point][2];
2290 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2292 * (fe_values_face[displacement].value(i, q_point) *
2293 * rhs_values * fe_values_face.JxW(q_point));
2297 * cell->get_dof_indices(local_dof_indices);
2298 * constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(
2299 * cell_rhs, local_dof_indices, newton_rhs);
2301 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2302 * newton_rhs_uncondensed(local_dof_indices[i]) += cell_rhs(i);
2305 * fraction_of_plastic_q_points_per_cell /= quadrature_formula.size();
2306 * newton_rhs.compress(VectorOperation::add);
2307 * newton_rhs_uncondensed.compress(VectorOperation::add);
2315 * <a name="step_42-PlasticityContactProblemsolve_newton_system"></a>
2316 * <h4>PlasticityContactProblem::solve_newton_system</h4>
2320 * The last piece before we can discuss the actual Newton iteration
2321 * on a single mesh is the solver for the linear systems. There are
2322 * a couple of complications that slightly obscure the code, but
2323 * mostly it is just setup then solve. Among the complications are:
2327 * - For the hanging nodes we have to apply
2328 * the AffineConstraints::set_zero function to newton_rhs.
2329 * This is necessary if a hanging node with solution value @f$x_0@f$
2330 * has one neighbor with value @f$x_1@f$ which is in contact with the
2331 * obstacle and one neighbor @f$x_2@f$ which is not in contact. Because
2332 * the update for the former will be prescribed, the hanging node constraint
2333 * will have an inhomogeneity and will look like @f$x_0 = x_1/2 +
2334 * \text{gap}/2@f$. So the corresponding entries in the right-hand-side are
2335 * non-zero with a meaningless value. These values we have to set to zero.
2336 * - Like in @ref step_40 "step-40", we need to shuffle between vectors that do and do
2337 * not have ghost elements when solving or using the solution.
2341 * The rest of the function is similar to @ref step_40 "step-40" and
2342 * @ref step_41 "step-41" except that we use a BiCGStab solver
2343 * instead of CG. This is due to the fact that for very small hardening
2344 * parameters @f$\gamma@f$, the linear system becomes almost semidefinite though
2345 * still symmetric. BiCGStab appears to have an easier time with such linear
2349 * template <int dim>
2350 * void PlasticityContactProblem<dim>::solve_newton_system()
2352 * TimerOutput::Scope t(computing_timer, "Solve");
2354 * TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
2355 * mpi_communicator);
2356 * distributed_solution = solution;
2358 * constraints_hanging_nodes.set_zero(distributed_solution);
2359 * constraints_hanging_nodes.set_zero(newton_rhs);
2361 * TrilinosWrappers::PreconditionAMG preconditioner;
2363 * TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner");
2365 * const std::vector<std::vector<bool>> constant_modes =
2366 * DoFTools::extract_constant_modes(dof_handler, ComponentMask());
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;
2377 * preconditioner.initialize(newton_matrix, additional_data);
2381 * TimerOutput::Scope t(computing_timer, "Solve: iterate");
2383 * TrilinosWrappers::MPI::Vector tmp(locally_owned_dofs, mpi_communicator);
2385 * const double relative_accuracy = 1e-8;
2386 * const double solver_tolerance =
2387 * relative_accuracy *
2388 * newton_matrix.residual(tmp, distributed_solution, newton_rhs);
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,
2397 * pcout << " Error: " << solver_control.initial_value() << " -> "
2398 * << solver_control.last_value() << " in "
2399 * << solver_control.last_step() << " Bicgstab iterations."
2403 * all_constraints.distribute(distributed_solution);
2405 * solution = distributed_solution;
2412 * <a name="step_42-PlasticityContactProblemsolve_newton"></a>
2413 * <h4>PlasticityContactProblem::solve_newton</h4>
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.
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.
2436 * template <int dim>
2437 * void PlasticityContactProblem<dim>::solve_newton()
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);
2450 * double residual_norm;
2451 * double previous_residual_norm = -std::numeric_limits<double>::max();
2453 * const double correct_sigma = sigma_0;
2455 * IndexSet old_active_set(active_set);
2457 * for (unsigned int newton_step = 1; newton_step <= 100; ++newton_step)
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);
2467 * pcout << ' ' << std::endl;
2468 * pcout << " Newton iteration " << newton_step << std::endl;
2469 * pcout << " Updating active set..." << std::endl;
2472 * TimerOutput::Scope t(computing_timer, "update active set");
2473 * update_solution_and_constraints();
2476 * pcout << " Assembling system... " << std::endl;
2477 * newton_matrix = 0;
2479 * assemble_newton_system(solution);
2481 * pcout << " Solving system... " << std::endl;
2482 * solve_newton_system();
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.
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
2515 * In either of these two cases, we bypass the line search and just
2516 * update residual and other vectors as necessary.
2519 *
if ((newton_step == 1) ||
2520 * (transfer_solution && newton_step == 2 &&
2521 * current_refinement_cycle == 0) ||
2522 * (!transfer_solution && newton_step == 2))
2524 * compute_nonlinear_residual(solution);
2525 * old_solution = solution;
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))
2536 * residual_norm = residual.l2_norm();
2538 * pcout <<
" Accepting Newton solution with residual: "
2539 * << residual_norm << std::endl;
2543 *
for (
unsigned int i = 0; i < 5; ++i)
2545 * distributed_solution = solution;
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);
2553 * locally_relevant_tmp_vector = tmp_vector;
2554 * compute_nonlinear_residual(locally_relevant_tmp_vector);
2555 * residual = newton_rhs;
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))
2565 * residual_norm = residual.l2_norm();
2568 * <<
" Residual of the non-contact part of the system: "
2569 * << residual_norm << std::endl
2570 * <<
" with a damping parameter alpha = " << alpha
2573 *
if (residual_norm < previous_residual_norm)
2577 * solution = tmp_vector;
2578 * old_solution = solution;
2581 * previous_residual_norm = residual_norm;
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:
2593 * mpi_communicator) == 0)
2595 * pcout <<
" Active set did not change!" <<
std::endl;
2596 *
if (residual_norm < 1e-10)
2600 * old_active_set = active_set;
2607 * <a name=
"step_42-PlasticityContactProblemrefine_grid"></a>
2608 * <h4>PlasticityContactProblem::refine_grid</h4>
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>.
2621 * template <int dim>
2622 * void PlasticityContactProblem<dim>::refine_grid()
2624 * if (refinement_strategy == RefinementStrategy::refine_global)
2626 * for (typename Triangulation<dim>::active_cell_iterator cell =
2627 * triangulation.begin_active();
2628 * cell != triangulation.end();
2630 * if (cell->is_locally_owned())
2631 * cell->set_refine_flag();
2635 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
2636 * KellyErrorEstimator<dim>::estimate(
2638 * QGauss<dim - 1>(fe.degree + 2),
2639 * std::map<types::boundary_id, const Function<dim> *>(),
2641 * estimated_error_per_cell);
2643 * parallel::distributed::GridRefinement ::refine_and_coarsen_fixed_number(
2644 * triangulation, estimated_error_per_cell, 0.3, 0.03);
2647 * triangulation.prepare_coarsening_and_refinement();
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);
2654 * triangulation.execute_coarsening_and_refinement();
2658 * if (transfer_solution)
2660 * TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
2661 * mpi_communicator);
2662 * solution_transfer.interpolate(distributed_solution);
2666 * enforce constraints to make the interpolated solution conforming on
2670 * constraints_hanging_nodes.distribute(distributed_solution);
2672 * solution = distributed_solution;
2673 * compute_nonlinear_residual(solution);
2681 * <a name="step_42-PlasticityContactProblemmove_mesh"></a>
2682 * <h4>PlasticityContactProblem::move_mesh</h4>
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.
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.
2703 * template <int dim>
2704 * void PlasticityContactProblem<dim>::move_mesh(
2705 * const TrilinosWrappers::MPI::Vector &displacement) const
2707 * std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
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)
2714 * vertex_touched[cell->vertex_index(v)] = true;
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));
2721 * cell->vertex(v) += vertex_displacement;
2730 * <a name="step_42-PlasticityContactProblemoutput_results"></a>
2731 * <h4>PlasticityContactProblem::output_results</h4>
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.
2749 * template <int dim>
2750 * void PlasticityContactProblem<dim>::output_results(
2751 * const unsigned int current_refinement_cycle)
2753 * TimerOutput::Scope t(computing_timer, "Graphical output");
2755 * pcout << " Writing graphical output... " << std::flush;
2757 * move_mesh(solution);
2761 * Calculation of the contact forces
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);
2775 * TrilinosWrappers::MPI::Vector lambda(locally_relevant_dofs,
2776 * mpi_communicator);
2777 * lambda = distributed_lambda;
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);
2786 * TrilinosWrappers::MPI::Vector active_set_vector(locally_relevant_dofs,
2787 * mpi_communicator);
2788 * active_set_vector = distributed_active_set_vector;
2790 * DataOut<dim> data_out;
2792 * data_out.attach_dof_handler(dof_handler);
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);
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");
2815 * data_out.add_data_vector(fraction_of_plastic_q_points_per_cell,
2816 * "fraction_of_plastic_q_points");
2818 * data_out.build_patches();
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.
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;
2846 * <a name=
"step_42-PlasticityContactProblemoutput_contact_force"></a>
2847 * <h4>PlasticityContactProblem::output_contact_force</h4>
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.
2865 *
template <
int dim>
2866 *
void PlasticityContactProblem<dim>::output_contact_force() const
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);
2877 * distributed_lambda(n) = 0;
2879 * constraints_hanging_nodes.distribute(distributed_lambda);
2882 * mpi_communicator);
2883 *
lambda = distributed_lambda;
2885 *
double contact_force = 0.0;
2887 *
const QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
2889 * face_quadrature_formula,
2892 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
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)
2901 * fe_values_face.
reinit(cell, face);
2903 * std::vector<Tensor<1, dim>> lambda_values(n_face_q_points);
2904 * fe_values_face[displacement].get_function_values(lambda,
2907 *
for (
unsigned int q_point = 0; q_point < n_face_q_points;
2910 * lambda_values[q_point][2] * fe_values_face.JxW(q_point);
2914 * pcout <<
"Contact force = " << contact_force << std::endl;
2921 * <a name=
"step_42-PlasticityContactProblemrun"></a>
2922 * <h4>PlasticityContactProblem::run</h4>
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.
2936 *
template <
int dim>
2937 *
void PlasticityContactProblem<dim>::run()
2939 * computing_timer.reset();
2940 *
for (; current_refinement_cycle < n_refinement_cycles;
2941 * ++current_refinement_cycle)
2946 * pcout << std::endl;
2947 * pcout <<
"Cycle " << current_refinement_cycle <<
':' << std::endl;
2949 *
if (current_refinement_cycle == 0)
2963 * output_results(current_refinement_cycle);
2965 * computing_timer.print_summary();
2966 * computing_timer.reset();
2970 * pcout <<
"Peak virtual memory used, resident in kB: " << stats.VmSize
2971 * <<
' ' << stats.VmRSS << std::endl;
2973 *
if (base_mesh ==
"box")
2974 * output_contact_force();
2982 * <a name=
"step_42-Thecodemaincodefunction"></a>
2983 * <h3>The <code>main</code> function</h3>
2987 * There really isn
't much to the <code>main()</code> function. It looks
2988 * like they always do:
2991 * int main(int argc, char *argv[])
2993 * using namespace dealii;
2994 * using namespace Step42;
2998 * ParameterHandler prm;
2999 * PlasticityContactProblem<3>::declare_parameters(prm);
3002 * std::cerr << "*** Call this program as <./step-42 input.prm>"
3007 * prm.parse_input(argv[1]);
3008 * Utilities::MPI::MPI_InitFinalize mpi_initialization(
3009 * argc, argv, numbers::invalid_unsigned_int);
3011 * PlasticityContactProblem<3> problem(prm);
3015 * catch (std::exception &exc)
3017 * std::cerr << std::endl
3019 * << "----------------------------------------------------"
3021 * std::cerr << "Exception on processing: " << std::endl
3022 * << exc.what() << std::endl
3023 * << "Aborting!" << std::endl
3024 * << "----------------------------------------------------"
3031 * std::cerr << std::endl
3033 * << "----------------------------------------------------"
3035 * std::cerr << "Unknown exception!" << std::endl
3036 * << "Aborting!" << std::endl
3037 * << "----------------------------------------------------"
3045<a name="step_42-Results"></a><h1>Results</h1>
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:
3054 Using output directory 'p1adaptive/
'
3056 transfer solution false
3059 Number of active cells: 512
3060 Number of degrees of freedom: 2187
3063 Updating active set...
3064 Size of active set: 1
3065 Assembling system...
3067 Error: 173.076 -> 1.64265e-06 in 7 Bicgstab iterations.
3068 Accepting Newton solution with residual: 1.64265e-06
3071 Updating active set...
3072 Size of active set: 1
3073 Assembling 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!
3080 Updating active set...
3081 Size of active set: 1
3082 Assembling 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!
3092 Updating active set...
3093 Size of active set: 1
3094 Assembling 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
3103+---------------------------------------------+------------+------------+
3104| Total wallclock time elapsed since start | 1.13s | |
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+---------------------------------+-----------+------------+------------+
3122Peak virtual memory used, resident in kB: 541884 77464
3123Contact force = 37.3058
3128 Number of active cells: 14652
3129 Number of degrees of freedom: 52497
3132 Updating active set...
3133 Size of active set: 145
3134 Assembling system...
3136 Error: 296.309 -> 2.72484e-06 in 10 Bicgstab iterations.
3137 Accepting Newton solution with residual: 2.72484e-06
3142 Updating active set...
3143 Size of active set: 145
3144 Assembling 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
3153+---------------------------------------------+------------+------------+
3154| Total wallclock time elapsed since start | 38.4s | |
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+---------------------------------+-----------+------------+------------+
3173Peak virtual memory used, resident in kB: 566052 105788
3174Contact force = 56.794
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
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
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.
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):
3200<table align="center">
3203 <img src="https://www.dealii.org/images/steps/developer/step-42.CellConstitutionColorbar.png">
3206 <img src="https://www.dealii.org/images/steps/developer/step-42.CellConstitutionBall2.png" alt="" width="70%">
3212 <img src="https://www.dealii.org/images/steps/developer/step-42.CellConstitutionLi2.png" alt="" alt="" width="70%">
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.
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.
3231Further discussion of results that can be obtained using this program is
3232provided in the publication mentioned at the very top of this page.
3235<a name="step-42-extensions"></a>
3236<a name="step_42-Possibilitiesforextensions"></a><h1>Possibilities for extensions</h1>
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:
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ächen im Mikrokosmos</i>, University
3250of Siegen, Germany, 2011.
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.
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.
3265<li> It might be of interest to implement more accurate, problem-adapted error
3266estimators for contact as well as for the plasticity.
3270<a name="step_42-PlainProg"></a>
3271<h1> The plain program</h1>
3272@include "step-42.cc"
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
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)
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)
void get_memory_stats(MemoryStats &stats)
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)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation