This tutorial depends on step-41, step-40.
This program was contributed by Jörg Frohne (University of Siegen, Germany) while on a long-term visit to Texas A&M University, with significant contributions by Timo Heister and Wolfgang Bangerth.
The code described here provides the basis for the numerical experiments shown in the following paper:
J. Frohne, T. Heister, W. Bangerth: Efficient numerical methods for the large-scale, parallel solution of elastoplastic contact problems. Accepted for publication in International Journal for Numerical Methods in Engineering, 2015.
Introduction
This example is an extension of step-41, considering a 3d contact problem with an elasto-plastic material behavior with isotropic hardening in three dimensions. In other words, it considers how a three-dimensional body deforms if one pushes into it a rigid obstacle (the contact problem) where deformation is governed by an elasto-plastic material law (a material that can only accommodate a certain maximal stress) that hardens as deformation accumulates. To show what we intend to do before going into too many details, let us just show a picture of what the solution will look like (the deformable body is a cube - only half of which is actually shown -, the obstacle corresponds to a Chinese character that is discussed below):
This problem description implies that we have to take care of an additional nonlinearity compared to step-41: the material behavior. Since we consider a three dimensional problem here, we also have to account for the fact that the contact area is at the boundary of the deformable body now, rather than in the interior. Finally, compared to step-41, we also have to deal with hanging nodes in both the handling of the linear system as well as of the inequality constraints as we would like to use an adaptive mesh; in the latter case, we will have to deal with prioritizing whether the constraints from the hanging nodes or from the inequalities are more important.
Since you can very easily reach a few million degrees of freedom in three dimensions, even with adaptive mesh refinement, we decided to use Trilinos and p4est to run our code in parallel, building on the framework of step-40 for the parallelization. Additional pointers for parallelization can be found in step-32.
Classical formulation
The classical formulation of the problem possesses the following form:
\begin{align*} \varepsilon(\mathbf u) &= A\sigma + \varepsilon^p & &\quad\text{in } \Omega,\\ -\textrm{div}\ \sigma &= \mathbf f & &\quad\text{in } \Omega,\\ \varepsilon^p:(\tau - \sigma) &\geq 0\quad\forall\tau\text{ with }\mathcal{F}(\tau)\leq 0 & &\quad\text{in } \Omega,\\ \mathbf u &= 0 & &\quad\text{on }\Gamma_D,\\ \sigma \mathbf n - [\mathbf n \cdot(\sigma \mathbf n)]\mathbf n &= 0, \quad \mathbf n \cdot (\sigma \mathbf n) \leq 0 & &\quad\text{on }\Gamma_C,\\ (\mathbf n \cdot (\sigma \mathbf n))(\mathbf n \cdot \mathbf u - g) &= 0,\quad \mathbf n \cdot \mathbf u - g \leq 0 & &\quad\text{on } \Gamma_C. \end{align*}
Here, the first of these equations defines the relationship between strain \(\varepsilon(\mathbf u)=\frac{1}{2}\left(\nabla \mathbf u + \nabla \mathbf u^T\right)\) and stress \(\sigma\) via the fourth-order compliance tensor \(A\); \(\varepsilon^p\) provides the plastic component of the strain to ensure that the stress does not exceed the yield stress. We will only consider isotropic materials for which \(A\) can be expressed in terms of the Lam\'e moduli \(\lambda\) and \(\mu\) or alternatively in terms of the bulk modulus \(\kappa\) and \(\mu\). The second equation is the force balance; we will here not consider any body forces and henceforth assume that \(\mathbf f=0\). The complementarity condition in the third line implies that \(\varepsilon^p=0\) if \(\mathcal{F}(\sigma)< 0\) but that \(\varepsilon^p\) may be a nonzero tensor if and only if \(\mathcal{F}(\sigma) = 0\), and in particular that in this case \(\varepsilon^p\) must point in the direction \(\partial \mathcal{F}(\sigma)/\partial \sigma\). The inequality \(\mathcal{F}(\sigma)\le 0\) is a statement of the fact that plastic materials can only support a finite amount of stress; in other words, they react with plastic deformations \(\varepsilon^p\) if external forces would result in a stress \(\sigma\) for which \(\mathcal{F}(\sigma)> 0\) would result. A typical form for this yield function is \(\mathcal{F}(\sigma)=|\sigma^D|-\sigma_{\text{yield}}\) where \(\tau^D = \tau - \dfrac{1}{3}tr(\tau)I\) is the deviatoric part of a tensor and \(|\cdot|\) denotes the Frobenius norm.
Further equations describe a fixed, zero displacement on \(\Gamma_D\) and that on the surface \(\Gamma_C=\partial\Omega\backslash\Gamma_D\) where contact may appear, the normal force \(\sigma_n=\mathbf n \cdot (\sigma(\mathbf u) \mathbf n)\) exerted by the obstacle is inward (no "pull" by the obstacle on our body) and with zero tangential component \(\mathbf \sigma_t= \sigma \mathbf n - \mathbf \sigma_n \mathbf n = \sigma \mathbf n - [\mathbf n \cdot(\sigma \mathbf n)]\mathbf n\). The last condition is again a complementarity condition that implies that on \(\Gamma_C\), the normal force can only be nonzero if the body is in contact with the obstacle; the second part describes the impenetrability of the obstacle and the body. The last two equations are commonly referred to as the Signorini contact conditions.
Most materials - especially metals - have the property that they show some hardening as a result of deformation. In other words, \(\sigma_{\text{yield}}\) increases with deformation. In practice, it is not the elastic deformation that results in hardening, but the plastic component. There are different constitutive laws to describe those material behaviors. The simplest one is called linear isotropic hardening described by the flow function \(\mathcal{F}(\sigma,\varepsilon^p) = \vert\sigma^D\vert - (\sigma_0 + \gamma^{\text{iso}}|\varepsilon^p|)\).
Reformulation as a variational inequality
It is generally rather awkward to deal with inequalities. Here, we have to deal with two: plasticity and the contact problem. As described in more detail in the paper mentioned at the top of this page, one can at least reformulate the plasticity in a way that makes it look like a nonlinearity that we can then treat with Newton's method. This is slightly tricky mathematically since the nonlinearity is not just some smooth function but instead has kinks where the stress reaches the yield stress; however, it can be shown for such semismooth functions that Newton's method still converges.
Without going into details, we will also get rid of the stress as an independent variable and instead work exclusively with the displacements \(\mathbf u\). Ultimately, the goal of this reformulation is that we will want to end up with a symmetric, positive definite problem - such as a linearized elasticity problem with spatially variable coefficients resulting from the plastic behavior - that needs to be solved in each Newton step. We want this because there are efficient and scalable methods for the solution of such linear systems, such as CG preconditioned with an algebraic multigrid. This is opposed to the saddle point problem akin to the mixed Laplace (see step-20) we would get were we to continue with the mixed formulation containing both displacements and stresses, and for which step-20 already gives a hint at how difficult it is to construct good solvers and preconditioners.
With this said, let us simply state the problem we obtain after reformulation (again, details can be found in the paper): Find a displacement \(\mathbf u \in V^+\) so that
\begin{align*} \left(P_{\Pi}(C\varepsilon(\mathbf u)),\varepsilon(\varphi) - \varepsilon(\mathbf u)\right) \geq 0,\quad \forall \varphi\in V^+. \end{align*}
where the projector \(P_\Pi\) is defined as
\begin{align*} P_{\Pi}(\tau) \dealcoloneq \begin{cases} \tau, & \text{if }\vert\tau^D\vert \leq \sigma_0,\\ \left[ \dfrac{\gamma^{\text{iso}}}{2\mu + \gamma^{\text{iso}}} + \left(1-\dfrac{\gamma^{\text{iso}}}{2\mu + \gamma^{\text{iso}}}\right)\dfrac{\sigma_0}{\vert\tau^D\vert} \right]\tau^D + \dfrac{1}{3}\text{trace}(\tau) I, & \text{if }\vert\tau^D\vert > \sigma_0, \end{cases} \end{align*}
and the space \(V^+\) is the space of all displacements that satisfy the contact condition:
\begin{align*} V &= \left\{ \mathbf u\in \left[H^1(\Omega)\right]^{d}: \mathbf u = 0 \text{ on } \Gamma_D\right\}, \\ V^+ &= \left\{ \mathbf u\in V: \mathbf n \cdot \mathbf u\leq g \text{ on } \Gamma_C \right\}. \end{align*}
In the actual code, we will use the abbreviation \(\gamma=\dfrac{\gamma^{\text{iso}}}{2\mu + \gamma^{\text{iso}}}\).
Given this formulation, we will apply two techniques:
- Run a Newton method to iterate out the nonlinearity in the projector.
- Run an active set method for the contact condition, in much the same way as we did in step-41.
A strict approach would keep the active set fixed while we iterate the Newton method to convergence (or maybe the other way around: find the final active set before moving on to the next Newton iteration). In practice, it turns out that it is sufficient to do only a single Newton step per active set iteration, and so we will iterate over them concurrently. We will also, every once in a while, refine the mesh.
A Newton method for the plastic nonlinearity
As mentioned, we will treat the nonlinearity of the operator \(P_\Pi\) by applying a Newton method, despite the fact that the operator is not differentiable in the strict sense. However, it satisfies the conditions of slant differentiability and this turns out to be enough for Newton's method to work. The resulting method then goes by the name semi-smooth Newton method, which sounds impressive but is, in reality, just a Newton method applied to a semi-smooth function with an appropriately chosen "derivative".
In the current case, we will run our iteration by solving in each iteration \(i\) the following equation (still an inequality, but linearized):
\begin{align*} \label{eq:linearization} \left(I_{\Pi}\varepsilon(\tilde {\mathbf u}^{i}), \varepsilon(\varphi) - \varepsilon(\tilde {\mathbf u}^{i})\right) \geq \left(\left(I_{\Pi}\varepsilon({\mathbf u}^{i-1}), \varepsilon(\varphi) - \varepsilon(\tilde {\mathbf u}^{i})\right) - \left(P_{\Pi}(C\varepsilon({\mathbf u}^{i-1})), \varepsilon(\varphi) - \varepsilon(\tilde {\mathbf u}^{i})\right)\right), \quad \forall \varphi\in V^+, \end{align*}
where the rank-4 tensor \(I_\Pi=I_\Pi(\varepsilon^D(\mathbf u^{i-1}))\) given by
\begin{align} I_\Pi = \begin{cases} C_{\mu} + C_{\kappa}, & \hspace{-8em} \text{if } \vert C\varepsilon^D(\mathbf u^{i-1}) \vert \leq \sigma_0, \\ \frac{\gamma^{\text{iso}}}{2\mu + \gamma^{\text{iso}}} C_{\mu} + \frac{\left(1-\frac{\gamma^{\text{iso}}}{2\mu + \gamma^{\text{iso}}}\right)\sigma_0}{\vert C\varepsilon^D(\mathbf u^{i-1}) \vert}\left(C_{\mu} - 2\mu\dfrac{C\varepsilon^D(\mathbf u^{i-1})\otimes C\varepsilon^D(\mathbf u^{i-1})}{\vert C\varepsilon^D(\mathbf u^{i-1})\vert^2}\right) + C_{\kappa}, & \text{ else.} \end{cases} \end{align}
This tensor is the (formal) linearization of \(P_\Pi(C\cdot)\) around \(\varepsilon^D(\mathbf u^{i-1})\). For the linear isotropic material we consider here, the bulk and shear components of the projector are given by
\begin{gather*} C_{\kappa} = \kappa I\otimes I, \qquad\qquad\qquad\qquad C_{\mu} = 2\mu\left(\mathbb{I} - \dfrac{1}{3} I\otimes I\right), \end{gather*}
where \(I\) and \(\mathbb{I}\) are the identity tensors of rank 2 and 4, respectively.
Note that this problem corresponds to a linear elastic contact problem where \(I_\Pi\) plays the role of the elasticity tensor \(C=A^{-1}\). Indeed, if the material is not plastic at a point, then \(I_\Pi=C\). However, at places where the material is plastic, \(I_\Pi\) is a spatially varying function. In any case, the system we have to solve for the Newton iterate \(\tilde {\mathbf u}^{i}\) gets us closer to the goal of rewriting our problem in a way that allows us to use well-known solvers and preconditioners for elliptic systems.
As a final note about the Newton method let us mention that as is common with Newton methods we need to globalize it by controlling the step length. In other words, while the system above solves for \(\tilde {\mathbf u}^{i}\), the final iterate will rather be
\begin{align*} {\mathbf u}^{i} = {\mathbf u}^{i-1} + \alpha_i (\tilde {\mathbf u}^{i} - {\mathbf u}^{i-1}) \end{align*}
where the difference in parentheses on the right takes the role of the traditional Newton direction, \(\delta {\mathbf u}^{i}\). We will determine \(\alpha^i\) using a standard line search.
Active Set methods to solve the saddle point problem
This linearized problem to be solved in each Newton step is essentially like in step-41. The only difference consists in the fact that the contact area is at the boundary instead of in the domain. But this has no further consequence so that we refer to the documentation of step-41 with the only hint that \(\mathcal{S}\) contains all the vertices at the contact boundary \(\Gamma_C\) this time. As there, what we need to do is keep a subset of degrees of freedom fixed, leading to additional constraints that one can write as a saddle point problem. However, as discussed in the paper, by writing these constraints in an appropriate way that removes the coupling between degrees of freedom, we end up with a set of nodes that essentially just have Dirichlet values attached to them.
Overall algorithm
The algorithm outlined above combines the damped semismooth Newton-method, which we use for the nonlinear constitutive law, with the semismooth Newton method for the contact. It works as follows:
-
Initialize the active and inactive sets \(\mathcal{A}_i\) and \(\mathcal{F}_i\) such that \(\mathcal{S} = \mathcal{A}_i \cup \mathcal{F}_i\) and \(\mathcal{A}_i \cap \mathcal{F}_i = \emptyset\) and set \(i = 1\). Here, \(\mathcal{S}\) is the set of all degrees of freedom located at the surface of the domain where contact may happen. The start value \(\hat U^0 \dealcoloneq P_{\mathcal{A}_k}(0)\) fulfills our obstacle condition, i.e., we project an initial zero displacement onto the set of feasible displacements.
-
Assemble the Newton matrix \(A_{pq} \dealcoloneq a'( U^{i-1};\varphi_p,\varphi_q)\) and the right-hand-side \(F(\hat U^{i-1})\). These correspond to the linearized Newton step, ignoring for the moment the contact inequality.
-
Find the primal-dual pair \((\tilde U^i,\Lambda^i)\) that satisfies
\begin{align*} A\tilde U^i + B\Lambda^i & = F, &\\ \left[B^T\tilde U^i\right]_p & = G_p & \forall p\in\mathcal{A}_i,\\ \Lambda^i_p & = 0 & \forall p\in\mathcal{F}_i. \end{align*}
As in step-41, we can obtain the solution to this problem by eliminating those degrees of freedom in \({\cal A}_i\) from the first equation and obtain a linear system \(\hat {\hat A}(U^{i-1}) \tilde U^i = \hat {\hat H}(U^{i-1})\).
-
Damp the Newton iteration for \(i>2\) by applying a line search and calculating a linear combination of \(U^{i-1}\) and \(\tilde U^i\). This requires finding an \(\alpha^i_l \dealcoloneq 2^{-l},(l=0,\ldots,10)\) so that
\begin{gather*}U^i \dealcoloneq \alpha^i_l\bar U^i + (1-\alpha^i_l)U^{i-1}\end{gather*}
satisfies
\begin{gather*} \vert {\hat R}\left({\mathbf u}^{i}\right) \vert < \vert {\hat R}\left({\mathbf u}^{i-1}\right) \vert. \end{gather*}
with \({\hat R}\left({\mathbf u}\right)=\left(P_{Pi}(C\varepsilon(u)),\varepsilon(\varphi^{i}_p\right)\) with the exceptions of (i) elements \(p\in\mathcal{A}_i\) where we set \({\hat R}\left({\mathbf u}\right)=0\), and (ii) elements that correspond to hanging nodes, which we eliminate in the usual manner.
-
Define the new active and inactive sets by
\begin{gather*}\mathcal{A}_{i+1} \dealcoloneq \lbrace p\in\mathcal{S}:\Lambda^i_p + c\left(\left[B^TU^i\right]_p - G_p\right) > 0\rbrace,\end{gather*}
\begin{gather*}\mathcal{F}_{i+1} \dealcoloneq \lbrace p\in\mathcal{S}:\Lambda^i_p + c\left(\left[B^TU^i\right]_p - G_p\right) \leq 0\rbrace.\end{gather*}
-
Project \(U^i\) so that it satisfies the contact inequality,
\begin{gather*}\hat U^i \dealcoloneq P_{\mathcal{A}_{i+1}}(U^i).\end{gather*}
Here, \(P_{\mathcal{A}}(U)\) is the projection of the active components in \(\mathcal{A}\) to the gap
\begin{gather*}P_{\mathcal{A}}(U)_p \dealcoloneq \begin{cases} U_p, & \textrm{if}\quad p\notin\mathcal{A}\\ g_{h,p}, & \textrm{if}\quad p\in\mathcal{A}, \end{cases}\end{gather*}
where \(g_{h,p}\) is the gap denoting the distance of the obstacle from the undisplaced configuration of the body.
-
If \(\mathcal{A}_{i+1} = \mathcal{A}_k\) and \(\left\| {\hat R}\left({\mathbf u}^{i}\right) \right\|_{\ell_2} < \delta\) then stop, else set \(i=i+1\) and go to step (1). This step ensures that we only stop iterations if both the correct active set has been found and the plasticity has been iterated to sufficient accuracy.
In step 3 of this algorithm, the matrix \(B\in\mathbb{R}^{n\times m}\), \(n>m\) describes the coupling of the bases for the displacements and Lagrange multiplier (contact forces) and it is not quadratic in our situation since \(\Lambda^k\) is only defined on \(\Gamma_C\), i.e., the surface where contact may happen. As shown in the paper, we can choose \(B\) to be a matrix that has only one entry per row, (see also Hüeber, Wohlmuth: A primal-dual active set strategy for non-linear multibody contact problems, Comput. Methods Appl. Mech. Engrg. 194, 2005, pp. 3147-3166). The vector \(G\) is defined by a suitable approximation \(g_h\) of the gap \(g\)
\begin{gather*}G_p = \begin{cases} g_{h,p}, & \text{if}\quad p\in\mathcal{S}\\ 0, & \text{if}\quad p\notin\mathcal{S}. \end{cases}\end{gather*}
Adaptive mesh refinement
Since we run our program in 3d, the computations the program performs are expensive. Consequently using adaptive mesh refinement is an important step towards staying within acceptable run-times. To make our lives easier we simply choose the KellyErrorEstimator that is already implemented in deal.II. We hand the solution vector to it which contains the displacement \(u\). As we will see in the results it yields a quite reasonable adaptive mesh for the contact zone as well as for plasticity.
Implementation
This tutorial is essentially a mixture of step-40 and step-41 but instead of PETSc we let the Trilinos library deal with parallelizing the linear algebra (like in step-32). Since we are trying to solve a similar problem like in step-41 we will use the same methods but now in parallel.
A difficulty is handling of the constraints from the Dirichlet conditions, hanging nodes and the inequality condition that arises from the contact. For this purpose we create three objects of type AffineConstraints that describe the various constraints and that we will combine as appropriate in each iteration.
Compared to step-41, the programs has a few new classes:
-
ConstitutiveLaw describes the plastic behavior of the material
-
SphereObstacle describes a sphere that serves as the obstacle that is pushed into the deformable, elastoplastic body. Whether this or the next class is used to describe the obstacle is determined from the input parameter file.
-
ChineseObstacle (and a helper class) is a class that allows us to read in an obstacle from a file. In the example we will show in the results section, this file will be 'obstacle_file.dat' and will correspond to data that shows the Chinese, Japanese or Korean symbol for force or power (see http://www.orientaloutpost.com/ : "This word can be used for motivation - it can also mean power/motion/propulsion/force. It can be anything internal or external that keeps you going. This is the safest way to express motivation in Chinese. If your audience is Japanese, please see the other entry for motivation. This is a word in Japanese and Korean, but it means "motive
power" or "kinetic energy" (without the motivation meaning that you are probably looking for)"). In essence, we will pretend that we have a stamp (i.e., a mask that corresponds to a flat bottomed obstacle with no pieces of intermediate height) that we press into the body. The symbol in question looks as follows (see also the picture at the top of this section on how the end result looks like):
Other than that, let us comment only on the following aspects:
The commented program
Include files
The set of include files is not much of a surprise any more at this time:
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/fe_field_function.h>
#include <fstream>
#include <iostream>
Finally, we include two system headers that let us create a directory for output files. The first header provides the mkdir function and the second lets us determine what happened if mkdir fails.
#include <sys/stat.h>
#include <cerrno>
namespace Step42
{
The ConstitutiveLaw class template
This class provides an interface for a constitutive law, i.e., for the relationship between strain \(\varepsilon(\mathbf u)\) and stress \(\sigma\). In this example we are using an elastoplastic material behavior with linear, isotropic hardening. Such materials are characterized by Young's modulus \(E\), Poisson's ratio \(\nu\), the initial yield stress \(\sigma_0\) and the isotropic hardening parameter \(\gamma\). For \(\gamma = 0\) we obtain perfect elastoplastic behavior.
As explained in the paper that describes this program, the first Newton steps are solved with a completely elastic material model to avoid having to deal with both nonlinearities (plasticity and contact) at once. To this end, this class has a function set_sigma_0() that we use later on to simply set \(\sigma_0\) to a very large value – essentially guaranteeing that the actual stress will not exceed it, and thereby producing an elastic material. When we are ready to use a plastic model, we set \(\sigma_0\) back to its proper value, using the same function. As a result of this approach, we need to leave sigma_0 as the only non-const member variable of this class.
template <int dim>
class ConstitutiveLaw
{
public:
ConstitutiveLaw(const double E,
const double nu,
const double sigma_0,
const double gamma);
void set_sigma_0(double sigma_zero);
bool get_stress_strain_tensor(
void get_linearized_stress_strain_tensors(
private:
const double kappa;
double sigma_0;
const double gamma;
};
The constructor of the ConstitutiveLaw class sets the required material parameter for our deformable body. Material parameters for elastic isotropic media can be defined in a variety of ways, such as the pair \(E, \nu\) (elastic modulus and Poisson's number), using the Lame parameters \(\lambda,mu\) or several other commonly used conventions. Here, the constructor takes a description of material parameters in the form of \(E,\nu\), but since this turns out to these are not the coefficients that appear in the equations of the plastic projector, we immediately convert them into the more suitable set \(\kappa,\mu\) of bulk and shear moduli. In addition, the constructor takes \(\sigma_0\) (the yield stress absent any plastic strain) and \(\gamma\) (the hardening parameter) as arguments. In this constructor, we also compute the two principal components of the stress-strain relation and its linearization.
template <int dim>
ConstitutiveLaw<dim>::ConstitutiveLaw(double E,
double nu,
double sigma_0,
double gamma)
: kappa(
E / (3 * (1 - 2 * nu)))
, sigma_0(sigma_0)
, gamma(gamma)
, stress_strain_tensor_kappa(kappa *
, stress_strain_tensor_mu(
3.0))
{}
template <int dim>
void ConstitutiveLaw<dim>::set_sigma_0(double sigma_zero)
{
sigma_0 = sigma_zero;
}
ConstitutiveLaw::get_stress_strain_tensor
This is the principal component of the constitutive law. It computes the fourth order symmetric tensor that relates the strain to the stress according to the projection given above, when evaluated at a particular strain point. We need this function to calculate the nonlinear residual in PlasticityContactProblem::residual_nl_system() where we multiply this tensor with the strain given in a quadrature point. The computations follow the formulas laid out in the introduction. In comparing the formulas there with the implementation below, recall that \(C_\mu : \varepsilon = \tau_D\) and that \(C_\kappa : \varepsilon = \kappa \text{trace}(\varepsilon) I = \frac 13 \text{trace}(\tau) I\).
The function returns whether the quadrature point is plastic to allow for some statistics downstream on how many of the quadrature points are plastic and how many are elastic.
template <int dim>
bool ConstitutiveLaw<dim>::get_stress_strain_tensor(
{
stress_tensor =
(stress_strain_tensor_kappa + stress_strain_tensor_mu) * strain_tensor;
const double deviator_stress_tensor_norm = deviator_stress_tensor.
norm();
stress_strain_tensor = stress_strain_tensor_mu;
if (deviator_stress_tensor_norm > sigma_0)
{
const double beta = sigma_0 / deviator_stress_tensor_norm;
stress_strain_tensor *= (gamma + (1 - gamma) * beta);
}
stress_strain_tensor += stress_strain_tensor_kappa;
return (deviator_stress_tensor_norm > sigma_0);
}
ConstitutiveLaw::get_linearized_stress_strain_tensors
This function returns the linearized stress strain tensor, linearized around the solution \(u^{i-1}\) of the previous Newton step \(i-1\). The parameter strain_tensor (commonly denoted \(\varepsilon(u^{i-1})\)) must be passed as an argument, and serves as the linearization point. The function returns the derivative of the nonlinear constitutive law in the variable stress_strain_tensor, as well as the stress-strain tensor of the linearized problem in stress_strain_tensor_linearized. See PlasticityContactProblem::assemble_nl_system where this function is used.
template <int dim>
void ConstitutiveLaw<dim>::get_linearized_stress_strain_tensors(
{
stress_tensor =
(stress_strain_tensor_kappa + stress_strain_tensor_mu) * strain_tensor;
stress_strain_tensor = stress_strain_tensor_mu;
stress_strain_tensor_linearized = stress_strain_tensor_mu;
const double deviator_stress_tensor_norm = deviator_stress_tensor.
norm();
if (deviator_stress_tensor_norm > sigma_0)
{
const double beta = sigma_0 / deviator_stress_tensor_norm;
stress_strain_tensor *= (gamma + (1 - gamma) * beta);
stress_strain_tensor_linearized *= (gamma + (1 - gamma) * beta);
deviator_stress_tensor /= deviator_stress_tensor_norm;
stress_strain_tensor_linearized -=
(1 - gamma) * beta * 2 *
mu *
}
stress_strain_tensor += stress_strain_tensor_kappa;
stress_strain_tensor_linearized += stress_strain_tensor_kappa;
}
Equation data: boundary forces, boundary values, obstacles
The following should be relatively standard. We need classes for the boundary forcing term (which we here choose to be zero) and boundary values on those part of the boundary that are not part of the contact surface (also chosen to be zero here).
namespace EquationData
{
template <int dim>
class BoundaryForce :
public Function<dim>
{
public:
BoundaryForce();
const unsigned int component = 0) const override;
};
template <int dim>
BoundaryForce<dim>::BoundaryForce()
{}
template <int dim>
double BoundaryForce<dim>::value(
const Point<dim> &,
const unsigned int) const
{
return 0.;
}
template <int dim>
void BoundaryForce<dim>::vector_value(
const Point<dim> &p,
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = BoundaryForce<dim>::value(p, c);
}
template <int dim>
class BoundaryValues :
public Function<dim>
{
public:
BoundaryValues();
const unsigned int component = 0) const override;
};
template <int dim>
BoundaryValues<dim>::BoundaryValues()
{}
template <int dim>
double BoundaryValues<dim>::value(
const Point<dim> &,
const unsigned int) const
{
return 0.;
}
The SphereObstacle class
The following class is the first of two obstacles that can be selected from the input file. It describes a sphere centered at position \(x=y=0.5, z=z_{\text{surface}}+0.59\) and radius \(r=0.6\), where \(z_{\text{surface}}\) is the vertical position of the (flat) surface of the deformable body. The function's value returns the location of the obstacle for a given \(x,y\) value if the point actually lies below the sphere, or a large positive value that can't possibly interfere with the deformation if it lies outside the "shadow" of the sphere.
template <int dim>
class SphereObstacle :
public Function<dim>
{
public:
SphereObstacle(const double z_surface);
const unsigned int component = 0) const override;
private:
const double z_surface;
};
template <int dim>
SphereObstacle<dim>::SphereObstacle(const double z_surface)
, z_surface(z_surface)
{}
template <int dim>
double SphereObstacle<dim>::value(
const Point<dim> & p,
const unsigned int component) const
{
if (component == 0)
return p(0);
else if (component == 1)
return p(1);
else if (component == 2)
{
if ((p(0) - 0.5) * (p(0) - 0.5) + (p(1) - 0.5) * (p(1) - 0.5) < 0.36)
return (-std::sqrt(0.36 - (p(0) - 0.5) * (p(0) - 0.5) -
(p(1) - 0.5) * (p(1) - 0.5)) +
z_surface + 0.59);
else
return 1000;
}
return 1e9;
preceding Assert
}
template <int dim>
void SphereObstacle<dim>::vector_value(
const Point<dim> &p,
{
values(c) = SphereObstacle<dim>::value(p, c);
}
The BitmapFile and ChineseObstacle classes
The following two classes describe the obstacle outlined in the introduction, i.e., the Chinese character. The first of the two, BitmapFile is responsible for reading in data from a picture file stored in pbm ascii format. This data will be bilinearly interpolated and thereby provides a function that describes the obstacle. (The code below shows how one can construct a function by interpolating between given data points. One could use the Functions::InterpolatedUniformGridData, introduced after this tutorial program was written, which does exactly what we want here, but it is instructive to see how to do it by hand.)
The data which we read from the file will be stored in a double std::vector named obstacle_data. This vector composes the base to calculate a piecewise bilinear function as a polynomial interpolation. The data we will read from a file consists of zeros (white) and ones (black).
The hx,hy variables denote the spacing between pixels in \(x\) and \(y\) directions. nx,ny are the numbers of pixels in each of these directions. get_value() returns the value of the image at a given location, interpolated from the adjacent pixel values.
template <int dim>
class BitmapFile
{
public:
BitmapFile(const std::string &name);
double get_value(const double x, const double y) const;
private:
std::vector<double> obstacle_data;
double hx, hy;
int nx, ny;
double get_pixel_value(const int i, const int j) const;
};
The constructor of this class reads in the data that describes the obstacle from the given file name.
template <int dim>
BitmapFile<dim>::BitmapFile(const std::string &name)
: obstacle_data(0)
, hx(0)
, hy(0)
, nx(0)
, ny(0)
{
std::ifstream f(name);
ExcMessage(std::string(
"Can't read from file <") + name +
">!"));
std::string temp;
f >> temp >> nx >> ny;
for (int k = 0; k < nx * ny; ++k)
{
double val;
f >> val;
obstacle_data.push_back(val);
}
hx = 1.0 / (nx - 1);
hy = 1.0 / (ny - 1);
std::cout << "Read obstacle from file <" << name << ">" << std::endl
<< "Resolution of the scanned obstacle picture: " << nx
<< " x " << ny << std::endl;
}
The following two functions return the value of a given pixel with coordinates \(i,j\), which we identify with the values of a function defined at positions i*hx, j*hy, and at arbitrary coordinates \(x,y\) where we do a bilinear interpolation between point values returned by the first of the two functions. In the second function, for each \(x,y\), we first compute the (integer) location of the nearest pixel coordinate to the bottom left of \(x,y\), and then compute the coordinates \(\xi,\eta\) within this pixel. We truncate both kinds of variables from both below and above to avoid problems when evaluating the function outside of its defined range as may happen due to roundoff errors.
template <int dim>
double BitmapFile<dim>::get_pixel_value(const int i, const int j) const
{
assert(i >= 0 && i < nx);
assert(j >= 0 && j < ny);
return obstacle_data[nx * (ny - 1 - j) + i];
}
template <int dim>
double BitmapFile<dim>::get_value(const double x, const double y) const
{
const int ix = std::min(std::max(static_cast<int>(x / hx), 0), nx - 2);
const int iy = std::min(std::max(static_cast<int>(y / hy), 0), ny - 2);
const double xi = std::min(std::max((x - ix * hx) / hx, 1.), 0.);
const double eta = std::min(std::max((y - iy * hy) / hy, 1.), 0.);
return ((1 - xi) * (1 - eta) * get_pixel_value(ix, iy) +
xi * (1 - eta) * get_pixel_value(ix + 1, iy) +
(1 - xi) * eta * get_pixel_value(ix, iy + 1) +
xi * eta * get_pixel_value(ix + 1, iy + 1));
}
Finally, this is the class that actually uses the class above. It has a BitmapFile object as a member that describes the height of the obstacle. As mentioned above, the BitmapFile class will provide us with a mask, i.e., values that are either zero or one (and, if you ask for locations between pixels, values that are interpolated between zero and one). This class translates this to heights that are either 0.001 below the surface of the deformable body (if the BitmapFile class reports a one at this location) or 0.999 above the obstacle (if the BitmapFile class reports a zero). The following function should then be self-explanatory.
template <int dim>
class ChineseObstacle :
public Function<dim>
{
public:
ChineseObstacle(const std::string &filename, const double z_surface);
const unsigned int component = 0) const override;
private:
const BitmapFile<dim> input_obstacle;
double z_surface;
};
template <int dim>
ChineseObstacle<dim>::ChineseObstacle(const std::string &filename,
const double z_surface)
, input_obstacle(filename)
, z_surface(z_surface)
{}
template <int dim>
double ChineseObstacle<dim>::value(
const Point<dim> & p,
const unsigned int component) const
{
if (component == 0)
return p(0);
if (component == 1)
return p(1);
else if (component == 2)
{
if (p(0) >= 0.0 && p(0) <= 1.0 && p(1) >= 0.0 && p(1) <= 1.0)
return z_surface + 0.999 - input_obstacle.get_value(p(0), p(1));
}
return 1e9;
preceding Assert
}
template <int dim>
void ChineseObstacle<dim>::vector_value(
const Point<dim> &p,
{
values(c) = ChineseObstacle<dim>::value(p, c);
}
}
The PlasticityContactProblem class template
This is the main class of this program and supplies all functions and variables needed to describe the nonlinear contact problem. It is close to step-41 but with some additional features like handling hanging nodes, a Newton method, using Trilinos and p4est for parallel distributed computing. To deal with hanging nodes makes life a bit more complicated since we need another AffineConstraints object now. We create a Newton method for the active set method for the contact situation and to handle the nonlinear operator for the constitutive law.
The general layout of this class is very much like for most other tutorial programs. To make our life a bit easier, this class reads a set of input parameters from an input file. These parameters, using the ParameterHandler class, are declared in the declare_parameters function (which is static so that it can be called before we even create an object of the current type), and a ParameterHandler object that has been used to read an input file will then be passed to the constructor of this class.
The remaining member functions are by and large as we have seen in several of the other tutorial programs, though with additions for the current nonlinear system. We will comment on their purpose as we get to them further below.
template <int dim>
class PlasticityContactProblem
{
public:
private:
void make_grid();
void setup_system();
void compute_dirichlet_constraints();
void update_solution_and_constraints();
void
void assemble_newton_system(
void compute_nonlinear_residual(
void solve_newton_system();
void solve_newton();
void refine_grid();
void output_results(const std::string &filename_base);
void output_contact_force() const;
As far as member variables are concerned, we start with ones that we use to indicate the MPI universe this program runs on, a stream we use to let exactly one processor produce output to the console (see step-17) and a variable that is used to time the various sections of the program:
MPI_Comm mpi_communicator;
The next group describes the mesh and the finite element space. In particular, for this parallel program, the finite element space has associated with it variables that indicate which degrees of freedom live on the current processor (the index sets, see also step-40 and the Parallel computing with multiple processors using distributed memory documentation module) as well as a variety of constraints: those imposed by hanging nodes, by Dirichlet boundary conditions, and by the active set of contact nodes. Of the three AffineConstraints variables defined here, the first only contains hanging node constraints, the second also those associated with Dirichlet boundary conditions, and the third these plus the contact constraints.
The variable active_set consists of those degrees of freedom constrained by the contact, and we use fraction_of_plastic_q_points_per_cell to keep track of the fraction of quadrature points on each cell where the stress equals the yield stress. The latter is only used to create graphical output showing the plastic zone, but not for any further computation; the variable is a member variable of this class since the information is computed as a by-product of computing the residual, but is used only much later. (Note that the vector is a vector of length equal to the number of active cells on the local mesh; it is never used to exchange information between processors and can therefore be a regular deal.II vector.)
const unsigned int n_initial_global_refinements;
const unsigned int fe_degree;
The next block of variables corresponds to the solution and the linear systems we need to form. In particular, this includes the Newton matrix and right hand side; the vector that corresponds to the residual (i.e., the Newton right hand side) but from which we have not eliminated the various constraints and that is used to determine which degrees of freedom need to be constrained in the next iteration; and a vector that corresponds to the diagonal of the \(B\) matrix briefly mentioned in the introduction and discussed in the accompanying paper.
The next block contains the variables that describe the material response:
const double e_modulus, nu, gamma, sigma_0;
ConstitutiveLaw<dim> constitutive_law;
And then there is an assortment of other variables that are used to identify the mesh we are asked to build as selected by the parameter file, the obstacle that is being pushed into the deformable body, the mesh refinement strategy, whether to transfer the solution from one mesh to the next, and how many mesh refinement cycles to perform. As possible, we mark these kinds of variables as const to help the reader identify which ones may or may not be modified later on (the output directory being an exception – it is never modified outside the constructor but it is awkward to initialize in the member-initializer-list following the colon in the constructor since there we have only one shot at setting it; the same is true for the mesh refinement criterion):
const std::string base_mesh;
const std::shared_ptr<const Function<dim>> obstacle;
struct RefinementStrategy
{
enum value
{
refine_global,
refine_percentage,
refine_fix_dofs
};
};
typename RefinementStrategy::value refinement_strategy;
const bool transfer_solution;
std::string output_dir;
const unsigned int n_refinement_cycles;
unsigned int current_refinement_cycle;
};
Implementation of the PlasticityContactProblem class
PlasticityContactProblem::declare_parameters
Let us start with the declaration of run-time parameters that can be selected in the input file. These values will be read back in the constructor of this class to initialize the member variables of this class:
template <int dim>
{
"polynomial degree",
"1",
"Polynomial degree of the FE_Q finite element space, typically 1 or 2.");
"2",
"Number of initial global mesh refinement steps before "
"the first computation.");
"refinement strategy",
"percentage",
"Mesh refinement strategy:\n"
" global: one global refinement\n"
" percentage: a fixed percentage of cells gets refined using the Kelly estimator.");
"5",
"Number of adaptive mesh refinement cycles to run.");
"obstacle",
"sphere",
"The name of the obstacle to use. This may either be 'sphere' if we should "
"use a spherical obstacle, or 'read from file' in which case the obstacle "
"will be read from a file named 'obstacle.pbm' that is supposed to be in "
"ASCII PBM format.");
"output directory",
"",
"Directory for output files (graphical output and benchmark "
"statistics). If empty, use the current directory.");
"transfer solution",
"false",
"Whether the solution should be used as a starting guess "
"for the next finer mesh. If false, then the iteration starts at "
"zero on every mesh.");
"box",
"Select the shape of the domain: 'box' or 'half sphere'");
}
The PlasticityContactProblem constructor
Given the declarations of member variables as well as the declarations of run-time parameters that are read from the input file, there is nothing surprising in this constructor. In the body we initialize the mesh refinement strategy and the output directory, creating such a directory if necessary.
template <int dim>
PlasticityContactProblem<dim>::PlasticityContactProblem(
: mpi_communicator(MPI_COMM_WORLD)
, computing_timer(MPI_COMM_WORLD,
pcout,
, n_initial_global_refinements(
prm.get_integer("number of initial refinements"))
, triangulation(mpi_communicator)
, fe_degree(prm.get_integer("polynomial degree"))
, dof_handler(triangulation)
, e_modulus(200000)
, nu(0.3)
, gamma(0.01)
, sigma_0(400.0)
, constitutive_law(e_modulus, nu, sigma_0, gamma)
, base_mesh(prm.get("base mesh"))
, obstacle(prm.get("obstacle") == "read from file" ?
new EquationData::ChineseObstacle<dim>(
"obstacle.pbm",
(base_mesh == "box" ? 1.0 : 0.5))) :
new EquationData::SphereObstacle<dim>(
base_mesh == "box" ? 1.0 : 0.5)))
, transfer_solution(prm.get_bool("transfer solution"))
, n_refinement_cycles(prm.get_integer("number of cycles"))
, current_refinement_cycle(0)
{
std::string strat = prm.
get(
"refinement strategy");
if (strat == "global")
refinement_strategy = RefinementStrategy::refine_global;
else if (strat == "percentage")
refinement_strategy = RefinementStrategy::refine_percentage;
else
output_dir = prm.
get(
"output directory");
if (output_dir != "" && *(output_dir.rbegin()) != '/')
output_dir += "/";
If necessary, create a new directory for the output.
{
const int ierr = mkdir(output_dir.c_str(), 0777);
}
pcout << " Using output directory '" << output_dir << "'" << std::endl;
pcout << " FE degree " << fe_degree << std::endl;
pcout << " transfer solution " << (transfer_solution ? "true" : "false")
<< std::endl;
}
PlasticityContactProblem::make_grid
The next block deals with constructing the starting mesh. We will use the following helper function and the first block of the make_grid() to construct a mesh that corresponds to a half sphere. deal.II has a function that creates such a mesh, but it is in the wrong location and facing the wrong direction, so we need to shift and rotate it a bit before using it.
For later reference, as described in the documentation of GridGenerator::half_hyper_ball(), the flat surface of the halfsphere has boundary indicator zero, while the remainder has boundary indicator one.
{
return {in(2), in(1), -in(0)};
}
template <int dim>
void PlasticityContactProblem<dim>::make_grid()
{
if (base_mesh == "half sphere")
{
const double radius = 0.8;
Since we will attach a different manifold below, we immediately clear the default manifold description:
triangulation.reset_all_manifolds();
triangulation.set_manifold(0, manifold_description);
}
Alternatively, create a hypercube mesh. After creating it, assign boundary indicators as follows:
In other words, the boundary indicators of the sides of the cube are 8. The boundary indicator of the bottom is 6 and the top has indicator 1. We set these by looping over all cells of all faces and looking at coordinate values of the cell center, and will make use of these indicators later when evaluating which boundary will carry Dirichlet boundary conditions or will be subject to potential contact. (In the current case, the mesh contains only a single cell, and all of its faces are on the boundary, so both the loop over all cells and the query whether a face is on the boundary are, strictly speaking, unnecessary; we retain them simply out of habit: this kind of code can be found in many programs in essentially this form.)
else
{
for (const auto &cell : triangulation.active_cell_iterators())
for (unsigned int face_no = 0;
face_no < GeometryInfo<dim>::faces_per_cell;
++face_no)
if (cell->face(face_no)->at_boundary())
{
if (std::fabs(cell->face(face_no)->center()[2] - p2[2]) < 1e-12)
cell->face(face_no)->set_boundary_id(1);
if (std::fabs(cell->face(face_no)->center()[0] - p1[0]) <
1e-12 ||
std::fabs(cell->face(face_no)->center()[0] - p2[0]) <
1e-12 ||
std::fabs(cell->face(face_no)->center()[1] - p1[1]) <
1e-12 ||
std::fabs(cell->face(face_no)->center()[1] - p2[1]) < 1e-12)
cell->face(face_no)->set_boundary_id(8);
if (std::fabs(cell->face(face_no)->center()[2] - p1[2]) < 1e-12)
cell->face(face_no)->set_boundary_id(6);
}
}
triangulation.refine_global(n_initial_global_refinements);
}
PlasticityContactProblem::setup_system
The next piece in the puzzle is to set up the DoFHandler, resize vectors and take care of various other status variables such as index sets and constraint matrices.
In the following, each group of operations is put into a brace-enclosed block that is being timed by the variable declared at the top of the block (the constructor of the TimerOutput::Scope variable starts the timed section, the destructor that is called at the end of the block stops it again).
template <int dim>
void PlasticityContactProblem<dim>::setup_system()
{
/ * setup dofs and get index sets for locally owned and relevant dofs * /
{
dof_handler.distribute_dofs(fe);
locally_owned_dofs = dof_handler.locally_owned_dofs();
locally_relevant_dofs.
clear();
locally_relevant_dofs);
}
/ * setup hanging nodes and Dirichlet constraints * /
{
constraints_hanging_nodes.reinit(locally_relevant_dofs);
constraints_hanging_nodes);
constraints_hanging_nodes.close();
pcout << " Number of active cells: "
<< triangulation.n_global_active_cells() << std::endl
<< " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
compute_dirichlet_constraints();
}
/ * initialization of vectors and the active set * /
{
solution.
reinit(locally_relevant_dofs, mpi_communicator);
newton_rhs.reinit(locally_owned_dofs, mpi_communicator);
newton_rhs_uncondensed.reinit(locally_owned_dofs, mpi_communicator);
diag_mass_matrix_vector.reinit(locally_owned_dofs, mpi_communicator);
fraction_of_plastic_q_points_per_cell.reinit(
triangulation.n_active_cells());
active_set.clear();
active_set.set_size(dof_handler.n_dofs());
}
Finally, we set up sparsity patterns and matrices. We temporarily (ab)use the system matrix to also build the (diagonal) matrix that we use in eliminating degrees of freedom that are in contact with the obstacle, but we then immediately set the Newton matrix back to zero.
{
mpi_communicator);
sp,
constraints_dirichlet_and_hanging_nodes,
false,
mpi_communicator));
sp.compress();
newton_matrix.reinit(sp);
assemble_mass_matrix_diagonal(mass_matrix);
const unsigned int start = (newton_rhs.local_range().first),
end = (newton_rhs.local_range().second);
for (unsigned int j = start; j < end; ++j)
diag_mass_matrix_vector(j) =
mass_matrix.diag_element(j);
}
}
PlasticityContactProblem::compute_dirichlet_constraints
This function, broken out of the preceding one, computes the constraints associated with Dirichlet-type boundary conditions and puts them into the constraints_dirichlet_and_hanging_nodes variable by merging with the constraints that come from hanging nodes.
As laid out in the introduction, we need to distinguish between two cases:
- If the domain is a box, we set the displacement to zero at the bottom, and allow vertical movement in z-direction along the sides. As shown in the
make_grid() function, the former corresponds to boundary indicator 6, the latter to 8.
- If the domain is a half sphere, then we impose zero displacement along the curved part of the boundary, associated with boundary indicator zero.
template <int dim>
void PlasticityContactProblem<dim>::compute_dirichlet_constraints()
{
constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs);
constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);
if (base_mesh == "box")
{
interpolate all components of the solution
dof_handler,
6,
EquationData::BoundaryValues<dim>(),
constraints_dirichlet_and_hanging_nodes,
interpolate x- and y-components of the solution (this is a bit mask, so apply operator| )
dof_handler,
8,
EquationData::BoundaryValues<dim>(),
constraints_dirichlet_and_hanging_nodes,
(fe.component_mask(x_displacement) |
fe.component_mask(y_displacement)));
}
else
dof_handler,
0,
EquationData::BoundaryValues<dim>(),
constraints_dirichlet_and_hanging_nodes,
constraints_dirichlet_and_hanging_nodes.close();
}
PlasticityContactProblem::assemble_mass_matrix_diagonal
The next helper function computes the (diagonal) mass matrix that is used to determine the active set of the active set method we use in the contact algorithm. This matrix is of mass matrix type, but unlike the standard mass matrix, we can make it diagonal (even in the case of higher order elements) by using a quadrature formula that has its quadrature points at exactly the same locations as the interpolation points for the finite element are located. We achieve this by using a QGaussLobatto quadrature formula here, along with initializing the finite element with a set of interpolation points derived from the same quadrature formula. The remainder of the function is relatively straightforward: we put the resulting matrix into the given argument; because we know the matrix is diagonal, it is sufficient to have a loop over only \(i\) and not over \(j\). Strictly speaking, we could even avoid multiplying the shape function's values at quadrature point q_point by itself because we know the shape value to be a vector with exactly one one which when dotted with itself yields one. Since this function is not time critical we add this term for clarity.
template <int dim>
void PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal(
{
face_quadrature_formula,
const unsigned int n_face_q_points = face_quadrature_formula.size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
cell->face(face)->boundary_id() == 1)
{
fe_values_face.reinit(cell, face);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
(fe_values_face[displacement].value(i, q_point) *
fe_values_face[displacement].value(i, q_point) *
fe_values_face.JxW(q_point));
cell->get_dof_indices(local_dof_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
local_dof_indices[i],
}
}
PlasticityContactProblem::update_solution_and_constraints
The following function is the first function we call in each Newton iteration in the solve_newton() function. What it does is to project the solution onto the feasible set and update the active set for the degrees of freedom that touch or penetrate the obstacle.
In order to function, we first need to do some bookkeeping: We need to write into the solution vector (which we can only do with fully distributed vectors without ghost elements) and we need to read the Lagrange multiplier and the elements of the diagonal mass matrix from their respective vectors (which we can only do with vectors that do have ghost elements), so we create the respective vectors. We then also initialize the constraints object that will contain constraints from contact and all other sources, as well as an object that contains an index set of all locally owned degrees of freedom that are part of the contact:
template <int dim>
void PlasticityContactProblem<dim>::update_solution_and_constraints()
{
std::vector<bool> dof_touched(dof_handler.n_dofs(), false);
mpi_communicator);
distributed_solution = solution;
mpi_communicator);
lambda = newton_rhs_uncondensed;
locally_relevant_dofs, mpi_communicator);
diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
all_constraints.reinit(locally_relevant_dofs);
active_set.clear();
The second part is a loop over all cells in which we look at each point where a degree of freedom is defined whether the active set condition is true and we need to add this degree of freedom to the active set of contact nodes. As we always do, if we want to evaluate functions at individual points, we do this with an FEValues object (or, here, an FEFaceValues object since we need to check contact at the surface) with an appropriately chosen quadrature object. We create this face quadrature object by choosing the "support points" of the shape functions defined on the faces of cells (for more on support points, see this glossary entry). As a consequence, we have as many quadrature points as there are shape functions per face and looping over quadrature points is equivalent to looping over shape functions defined on a face. With this, the code looks as follows:
Quadrature<dim - 1> face_quadrature(fe.get_unit_face_support_points());
face_quadrature,
const unsigned int dofs_per_face = fe.dofs_per_face;
const unsigned int n_face_q_points = face_quadrature.size();
std::vector<types::global_dof_index> dof_indices(dofs_per_face);
for (const auto &cell : dof_handler.active_cell_iterators())
if (!cell->is_artificial())
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
cell->face(face)->boundary_id() == 1)
{
fe_values_face.reinit(cell, face);
cell->face(face)->get_dof_indices(dof_indices);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
At each quadrature point (i.e., at each support point of a degree of freedom located on the contact boundary), we then ask whether it is part of the z-displacement degrees of freedom and if we haven't encountered this degree of freedom yet (which can happen for those on the edges between faces), we need to evaluate the gap between the deformed object and the obstacle. If the active set condition is true, then we add a constraint to the AffineConstraints object that the next Newton update needs to satisfy, set the solution vector's corresponding element to the correct value, and add the index to the IndexSet object that stores which degree of freedom is part of the contact:
const unsigned int component =
fe.face_system_to_component_index(q_point).first;
const unsigned int index_z = dof_indices[q_point];
if ((component == 2) && (dof_touched[index_z] == false))
{
dof_touched[index_z] = true;
fe_values_face.quadrature_point(q_point);
const double obstacle_value =
obstacle->value(this_support_point, 2);
const double solution_here = solution(index_z);
const double undeformed_gap =
obstacle_value - this_support_point(2);
const double c = 100.0 * e_modulus;
if ((lambda(index_z) /
diag_mass_matrix_vector_relevant(index_z) +
c * (solution_here - undeformed_gap) >
0) &&
!constraints_hanging_nodes.is_constrained(index_z))
{
all_constraints.add_line(index_z);
all_constraints.set_inhomogeneity(index_z,
undeformed_gap);
distributed_solution(index_z) = undeformed_gap;
active_set.add_index(index_z);
}
}
}
}
At the end of this function, we exchange data between processors updating those ghost elements in the solution variable that have been written by other processors. We then merge the Dirichlet constraints and those from hanging nodes into the AffineConstraints object that already contains the active set. We finish the function by outputting the total number of actively constrained degrees of freedom for which we sum over the number of actively constrained degrees of freedom owned by each of the processors. This number of locally owned constrained degrees of freedom is of course the number of elements of the intersection of the active set and the set of locally owned degrees of freedom, which we can get by using operator& on two IndexSets:
solution = distributed_solution;
all_constraints.close();
all_constraints.merge(constraints_dirichlet_and_hanging_nodes);
pcout << " Size of active set: "
mpi_communicator)
<< std::endl;
}
PlasticityContactProblem::assemble_newton_system
Given the complexity of the problem, it may come as a bit of a surprise that assembling the linear system we have to solve in each Newton iteration is actually fairly straightforward. The following function builds the Newton right hand side and Newton matrix. It looks fairly innocent because the heavy lifting happens in the call to ConstitutiveLaw::get_linearized_stress_strain_tensors() and in particular in AffineConstraints::distribute_local_to_global(), using the constraints we have previously computed.
template <int dim>
void PlasticityContactProblem<dim>::assemble_newton_system(
{
QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
quadrature_formula,
face_quadrature_formula,
const unsigned int n_q_points = quadrature_formula.
size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
const EquationData::BoundaryForce<dim> boundary_force;
std::vector<Vector<double>> boundary_force_values(n_face_q_points,
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_rhs = 0;
std::vector<SymmetricTensor<2, dim>> strain_tensor(n_q_points);
fe_values[displacement].get_function_symmetric_gradients(
linearization_point, strain_tensor);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
constitutive_law.get_linearized_stress_strain_tensors(
strain_tensor[q_point],
stress_strain_tensor_linearized,
stress_strain_tensor);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
Having computed the stress-strain tensor and its linearization, we can now put together the parts of the matrix and right hand side. In both, we need the linearized stress-strain tensor times the symmetric gradient of \(\varphi_i\), i.e. the term \(I_\Pi\varepsilon(\varphi_i)\), so we introduce an abbreviation of this term. Recall that the matrix corresponds to the bilinear form \(A_{ij}=(I_\Pi\varepsilon(\varphi_i),\varepsilon(\varphi_j))\) in the notation of the accompanying publication, whereas the right hand side is \(F_i=([I_\Pi-P_\Pi C]\varepsilon(\varphi_i),\varepsilon(\mathbf u))\) where \(u\) is the current linearization points (typically the last solution). This might suggest that the right hand side will be zero if the material is completely elastic (where \(I_\Pi=P_\Pi\)) but this ignores the fact that the right hand side will also contain contributions from non-homogeneous constraints due to the contact.
The code block that follows this adds contributions that are due to boundary forces, should there be any.
stress_strain_tensor_linearized *
fe_values[displacement].symmetric_gradient(i, q_point);
for (unsigned int j = 0; j < dofs_per_cell; ++j)
(stress_phi_i *
fe_values[displacement].symmetric_gradient(j, q_point) *
cell_rhs(i) +=
((stress_phi_i -
stress_strain_tensor *
fe_values[displacement].symmetric_gradient(i,
q_point)) *
strain_tensor[q_point] * fe_values.
JxW(q_point));
}
}
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
cell->face(face)->boundary_id() == 1)
{
fe_values_face.
reinit(cell, face);
boundary_force.vector_value_list(
fe_values_face.get_quadrature_points(),
boundary_force_values);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
rhs_values[2] = boundary_force_values[q_point][2];
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) +=
(fe_values_face[displacement].value(i, q_point) *
rhs_values * fe_values_face.JxW(q_point));
}
}
cell->get_dof_indices(local_dof_indices);
all_constraints.distribute_local_to_global(cell_matrix,
cell_rhs,
local_dof_indices,
newton_matrix,
newton_rhs,
true);
}
}
PlasticityContactProblem::compute_nonlinear_residual
The following function computes the nonlinear residual of the equation given the current solution (or any other linearization point). This is needed in the linear search algorithm where we need to try various linear combinations of previous and current (trial) solution to compute the (real, globalized) solution of the current Newton step.
That said, in a slight abuse of the name of the function, it actually does significantly more. For example, it also computes the vector that corresponds to the Newton residual but without eliminating constrained degrees of freedom. We need this vector to compute contact forces and, ultimately, to compute the next active set. Likewise, by keeping track of how many quadrature points we encounter on each cell that show plastic yielding, we also compute the fraction_of_plastic_q_points_per_cell vector that we can later output to visualize the plastic zone. In both of these cases, the results are not necessary as part of the line search, and so we may be wasting a small amount of time computing them. At the same time, this information appears as a natural by-product of what we need to do here anyway, and we want to collect it once at the end of each Newton step, so we may as well do it here.
The actual implementation of this function should be rather obvious:
template <int dim>
void PlasticityContactProblem<dim>::compute_nonlinear_residual(
{
QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
quadrature_formula,
face_quadrature_formula,
const unsigned int n_q_points = quadrature_formula.
size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
const EquationData::BoundaryForce<dim> boundary_force;
std::vector<Vector<double>> boundary_force_values(n_face_q_points,
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
newton_rhs = 0;
newton_rhs_uncondensed = 0;
fraction_of_plastic_q_points_per_cell = 0;
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_rhs = 0;
std::vector<SymmetricTensor<2, dim>> strain_tensors(n_q_points);
fe_values[displacement].get_function_symmetric_gradients(
linearization_point, strain_tensors);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const bool q_point_is_plastic =
constitutive_law.get_stress_strain_tensor(
strain_tensors[q_point], stress_strain_tensor);
if (q_point_is_plastic)
++fraction_of_plastic_q_points_per_cell(
cell->active_cell_index());
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
cell_rhs(i) -=
(strain_tensors[q_point] * stress_strain_tensor *
fe_values[displacement].symmetric_gradient(i, q_point) *
rhs_values = 0;
cell_rhs(i) += (fe_values[displacement].value(i, q_point) *
rhs_values * fe_values.
JxW(q_point));
}
}
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
cell->face(face)->boundary_id() == 1)
{
fe_values_face.
reinit(cell, face);
boundary_force.vector_value_list(
fe_values_face.get_quadrature_points(),
boundary_force_values);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
rhs_values[2] = boundary_force_values[q_point][2];
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) +=
(fe_values_face[displacement].value(i, q_point) *
rhs_values * fe_values_face.JxW(q_point));
}
}
cell->get_dof_indices(local_dof_indices);
constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(
cell_rhs, local_dof_indices, newton_rhs);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
newton_rhs_uncondensed(local_dof_indices[i]) += cell_rhs(i);
}
fraction_of_plastic_q_points_per_cell /= quadrature_formula.
size();
}
PlasticityContactProblem::solve_newton_system
The last piece before we can discuss the actual Newton iteration on a single mesh is the solver for the linear systems. There are a couple of complications that slightly obscure the code, but mostly it is just setup then solve. Among the complications are:
- For the hanging nodes we have to apply the AffineConstraints::set_zero function to newton_rhs. This is necessary if a hanging node with solution value \(x_0\) has one neighbor with value \(x_1\) which is in contact with the obstacle and one neighbor \(x_2\) which is not in contact. Because the update for the former will be prescribed, the hanging node constraint will have an inhomogeneity and will look like \(x_0 = x_1/2 + \text{gap}/2\). So the corresponding entries in the right-hand-side are non-zero with a meaningless value. These values we have to set to zero.
- Like in step-40, we need to shuffle between vectors that do and do not have ghost elements when solving or using the solution.
The rest of the function is similar to step-40 and step-41 except that we use a BiCGStab solver instead of CG. This is due to the fact that for very small hardening parameters \(\gamma\), the linear system becomes almost semidefinite though still symmetric. BiCGStab appears to have an easier time with such linear systems.
template <int dim>
void PlasticityContactProblem<dim>::solve_newton_system()
{
mpi_communicator);
distributed_solution = solution;
constraints_hanging_nodes.set_zero(distributed_solution);
constraints_hanging_nodes.set_zero(newton_rhs);
{
std::vector<std::vector<bool>> constant_modes;
constant_modes);
preconditioner.
initialize(newton_matrix, additional_data);
}
{
const double relative_accuracy = 1
e-8;
const double solver_tolerance =
relative_accuracy *
newton_matrix.residual(tmp, distributed_solution, newton_rhs);
SolverControl solver_control(newton_matrix.m(), solver_tolerance);
solver.solve(newton_matrix,
distributed_solution,
newton_rhs,
preconditioner);
pcout << " Error: " << solver_control.initial_value() << " -> "
<< solver_control.last_value() << " in "
<< solver_control.last_step() << " Bicgstab iterations."
<< std::endl;
}
all_constraints.distribute(distributed_solution);
solution = distributed_solution;
}
PlasticityContactProblem::solve_newton
This is, finally, the function that implements the damped Newton method on the current mesh. There are two nested loops: the outer loop for the Newton iteration and the inner loop for the line search which will be used only if necessary. To obtain a good and reasonable starting value we solve an elastic problem in the very first Newton step on each mesh (or only on the first mesh if we transfer solutions between meshes). We do so by setting the yield stress to an unreasonably large value in these iterations and then setting it back to the correct value in subsequent iterations.
Other than this, the top part of this function should be reasonably obvious. We initialize the variable previous_residual_norm to the most negative value representable with double precision numbers so that the comparison whether the current residual is less than that of the previous step will always fail in the first step.
template <int dim>
void PlasticityContactProblem<dim>::solve_newton()
{
mpi_communicator);
mpi_communicator);
mpi_communicator);
locally_relevant_dofs, mpi_communicator);
mpi_communicator);
double residual_norm;
double previous_residual_norm = -std::numeric_limits<double>::max();
const double correct_sigma = sigma_0;
for (unsigned int newton_step = 1; newton_step <= 100; ++newton_step)
{
if (newton_step == 1 &&
((transfer_solution && current_refinement_cycle == 0) ||
!transfer_solution))
constitutive_law.set_sigma_0(1e+10);
else if (newton_step == 2 || current_refinement_cycle > 0 ||
!transfer_solution)
constitutive_law.set_sigma_0(correct_sigma);
pcout << " " << std::endl;
pcout << " Newton iteration " << newton_step << std::endl;
pcout << " Updating active set..." << std::endl;
{
update_solution_and_constraints();
}
pcout << " Assembling system... " << std::endl;
newton_matrix = 0;
newton_rhs = 0;
assemble_newton_system(solution);
pcout << " Solving system... " << std::endl;
solve_newton_system();
It gets a bit more hairy after we have computed the trial solution \(\tilde{\mathbf u}\) of the current Newton step. We handle a highly nonlinear problem so we have to damp Newton's method using a line search. To understand how we do this, recall that in our formulation, we compute a trial solution in each Newton step and not the update between old and new solution. Since the solution set is a convex set, we will use a line search that tries linear combinations of the previous and the trial solution to guarantee that the damped solution is in our solution set again. At most we apply 5 damping steps.
There are exceptions to when we use a line search. First, if this is the first Newton step on any mesh, then we don't have any point to compare the residual to, so we always accept a full step. Likewise, if this is the second Newton step on the first mesh (or the second on any mesh if we don't transfer solutions from mesh to mesh), then we have computed the first of these steps using just an elastic model (see how we set the yield stress sigma to an unreasonably large value above). In this case, the first Newton solution was a purely elastic one, the second one a plastic one, and any linear combination would not necessarily be expected to lie in the feasible set – so we just accept the solution we just got.
In either of these two cases, we bypass the line search and just update residual and other vectors as necessary.
if ((newton_step == 1) ||
(transfer_solution && newton_step == 2 &&
current_refinement_cycle == 0) ||
(!transfer_solution && newton_step == 2))
{
compute_nonlinear_residual(solution);
old_solution = solution;
residual = newton_rhs;
const unsigned int start_res = (residual.local_range().first),
end_res = (residual.local_range().second);
for (unsigned int n = start_res; n < end_res; ++n)
if (all_constraints.is_inhomogeneously_constrained(n))
residual(n) = 0;
residual_norm = residual.l2_norm();
pcout << " Accepting Newton solution with residual: "
<< residual_norm << std::endl;
}
else
{
for (unsigned int i = 0; i < 5; ++i)
{
distributed_solution = solution;
const double alpha = std::pow(0.5, static_cast<double>(i));
tmp_vector = old_solution;
tmp_vector.sadd(1 - alpha, alpha, distributed_solution);
locally_relevant_tmp_vector = tmp_vector;
compute_nonlinear_residual(locally_relevant_tmp_vector);
residual = newton_rhs;
const unsigned int start_res = (residual.local_range().first),
end_res = (residual.local_range().second);
for (unsigned int n = start_res; n < end_res; ++n)
if (all_constraints.is_inhomogeneously_constrained(n))
residual(n) = 0;
residual_norm = residual.l2_norm();
pcout
<< " Residual of the non-contact part of the system: "
<< residual_norm << std::endl
<< " with a damping parameter alpha = " << alpha
<< std::endl;
if (residual_norm < previous_residual_norm)
break;
}
solution = tmp_vector;
old_solution = solution;
}
previous_residual_norm = residual_norm;
The final step is to check for convergence. If the active set has not changed across all processors and the residual is less than a threshold of \(10^{-10}\), then we terminate the iteration on the current mesh:
mpi_communicator) == 0)
{
pcout << " Active set did not change!" << std::endl;
if (residual_norm < 1e-10)
break;
}
old_active_set = active_set;
}
}
PlasticityContactProblem::refine_grid
If you've made it this far into the deal.II tutorial, the following function refining the mesh should not pose any challenges to you any more. It refines the mesh, either globally or using the Kelly error estimator, and if so asked also transfers the solution from the previous to the next mesh. In the latter case, we also need to compute the active set and other quantities again, for which we need the information computed by compute_nonlinear_residual().
template <int dim>
void PlasticityContactProblem<dim>::refine_grid()
{
if (refinement_strategy == RefinementStrategy::refine_global)
{
triangulation.begin_active();
cell != triangulation.end();
++cell)
if (cell->is_locally_owned())
cell->set_refine_flag();
}
else
{
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
dof_handler,
solution,
estimated_error_per_cell);
parallel::distributed::GridRefinement ::refine_and_coarsen_fixed_number(
triangulation, estimated_error_per_cell, 0.3, 0.03);
}
triangulation.prepare_coarsening_and_refinement();
solution_transfer(dof_handler);
if (transfer_solution)
solution_transfer.prepare_for_coarsening_and_refinement(solution);
triangulation.execute_coarsening_and_refinement();
setup_system();
if (transfer_solution)
{
mpi_communicator);
solution_transfer.interpolate(distributed_solution);
enforce constraints to make the interpolated solution conforming on the new mesh:
constraints_hanging_nodes.distribute(distributed_solution);
solution = distributed_solution;
compute_nonlinear_residual(solution);
}
}
PlasticityContactProblem::move_mesh
The remaining three functions before we get to run() have to do with generating output. The following one is an attempt at showing the deformed body in its deformed configuration. To this end, this function takes a displacement vector field and moves every vertex of the (local part) of the mesh by the previously computed displacement. We will call this function with the current displacement field before we generate graphical output, and we will call it again after generating graphical output with the negative displacement field to undo the changes to the mesh so made.
The function itself is pretty straightforward. All we have to do is keep track which vertices we have already touched, as we encounter the same vertices multiple times as we loop over cells.
template <int dim>
void PlasticityContactProblem<dim>::move_mesh(
{
std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
if (vertex_touched[cell->vertex_index(v)] == false)
{
vertex_touched[cell->vertex_index(v)] = true;
for (
unsigned int d = 0;
d < dim; ++
d)
vertex_displacement[d] =
displacement(cell->vertex_dof_index(v, d));
cell->vertex(v) += vertex_displacement;
}
}
PlasticityContactProblem::output_results
Next is the function we use to actually generate graphical output. The function is a bit tedious, but not actually particularly complicated. It moves the mesh at the top (and moves it back at the end), then computes the contact forces along the contact surface. We can do so (as shown in the accompanying paper) by taking the untreated residual vector and identifying which degrees of freedom correspond to those with contact by asking whether they have an inhomogeneous constraints associated with them. As always, we need to be mindful that we can only write into completely distributed vectors (i.e., vectors without ghost elements) but that when we want to generate output, we need vectors that do indeed have ghost entries for all locally relevant degrees of freedom.
template <int dim>
void PlasticityContactProblem<dim>::output_results(
const std::string &filename_base)
{
pcout << " Writing graphical output... " << std::flush;
move_mesh(solution);
Calculation of the contact forces
mpi_communicator);
const unsigned int start_res = (newton_rhs_uncondensed.local_range().first),
end_res = (newton_rhs_uncondensed.local_range().second);
for (unsigned int n = start_res; n < end_res; ++n)
if (all_constraints.is_inhomogeneously_constrained(n))
distributed_lambda(n) =
newton_rhs_uncondensed(n) / diag_mass_matrix_vector(n);
constraints_hanging_nodes.distribute(distributed_lambda);
mpi_communicator);
lambda = distributed_lambda;
locally_owned_dofs, mpi_communicator);
distributed_active_set_vector = 0.;
for (const auto index : active_set)
distributed_active_set_vector[index] = 1.;
mpi_communicator);
active_set_vector = distributed_active_set_vector;
const std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(
std::vector<std::string>(dim, "displacement"),
data_component_interpretation);
std::vector<std::string>(dim, "contact_force"),
data_component_interpretation);
std::vector<std::string>(dim, "active_set"),
data_component_interpretation);
for (unsigned int i = 0; i < subdomain.size(); ++i)
subdomain(i) = triangulation.locally_owned_subdomain();
"fraction_of_plastic_q_points");
In the remainder of the function, we generate one VTU file on every processor, indexed by the subdomain id of this processor. On the first processor, we then also create a .pvtu file that indexes all of the VTU files so that the entire set of output files can be read at once. These .pvtu are used by Paraview to describe an entire parallel computation's output files. We then do the same again for the competitor of Paraview, the Visit visualization program, by creating a matching .visit file.
const std::string filename =
(output_dir + filename_base + "-" +
std::ofstream output_vtu((filename + ".vtu"));
pcout << output_dir + filename_base << ".pvtu" << std::endl;
{
std::vector<std::string> filenames;
for (unsigned int i = 0;
++i)
filenames.push_back(filename_base + "-" +
std::ofstream pvtu_master_output(
(output_dir + filename_base + ".pvtu"));
std::ofstream visit_master_output(
(output_dir + filename_base + ".visit"));
}
tmp *= -1;
move_mesh(tmp);
}
PlasticityContactProblem::output_contact_force
This last auxiliary function computes the contact force by calculating an integral over the contact pressure in z-direction over the contact area. For this purpose we set the contact pressure lambda to 0 for all inactive dofs (whether a degree of freedom is part of the contact is determined just as we did in the previous function). For all active dofs, lambda contains the quotient of the nonlinear residual (newton_rhs_uncondensed) and corresponding diagonal entry of the mass matrix (diag_mass_matrix_vector). Because it is not unlikely that hanging nodes show up in the contact area it is important to apply constraints_hanging_nodes.distribute to the distributed_lambda vector.
template <int dim>
void PlasticityContactProblem<dim>::output_contact_force() const
{
mpi_communicator);
const unsigned int start_res = (newton_rhs_uncondensed.local_range().first),
end_res = (newton_rhs_uncondensed.local_range().second);
for (unsigned int n = start_res; n < end_res; ++n)
if (all_constraints.is_inhomogeneously_constrained(n))
distributed_lambda(n) =
newton_rhs_uncondensed(n) / diag_mass_matrix_vector(n);
else
distributed_lambda(n) = 0;
constraints_hanging_nodes.distribute(distributed_lambda);
mpi_communicator);
lambda = distributed_lambda;
double contact_force = 0.0;
QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
face_quadrature_formula,
const unsigned int n_face_q_points = face_quadrature_formula.size();
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
cell->face(face)->boundary_id() == 1)
{
fe_values_face.reinit(cell, face);
std::vector<Tensor<1, dim>> lambda_values(n_face_q_points);
fe_values_face[displacement].get_function_values(lambda,
lambda_values);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
contact_force +=
lambda_values[q_point][2] * fe_values_face.JxW(q_point);
}
pcout << "Contact force = " << contact_force << std::endl;
}
PlasticityContactProblem::run
As in all other tutorial programs, the run() function contains the overall logic. There is not very much to it here: in essence, it performs the loops over all mesh refinement cycles, and within each, hands things over to the Newton solver in solve_newton() on the current mesh and calls the function that creates graphical output for the so-computed solution. It then outputs some statistics concerning both run times and memory consumption that has been collected over the course of computations on this mesh.
template <int dim>
void PlasticityContactProblem<dim>::run()
{
computing_timer.reset();
for (; current_refinement_cycle < n_refinement_cycles;
++current_refinement_cycle)
{
{
pcout << std::endl;
pcout << "Cycle " << current_refinement_cycle << ':' << std::endl;
if (current_refinement_cycle == 0)
{
make_grid();
setup_system();
}
else
{
refine_grid();
}
}
solve_newton();
output_results(std::string("solution-") +
computing_timer.print_summary();
computing_timer.reset();
pcout <<
"Peak virtual memory used, resident in kB: " << stats.
VmSize <<
" " << stats.
VmRSS << std::endl;
if (base_mesh == "box")
output_contact_force();
}
}
}
The main function
There really isn't much to the main() function. It looks like they always do:
int main(int argc, char *argv[])
{
using namespace Step42;
try
{
PlasticityContactProblem<3>::declare_parameters(prm);
if (argc != 2)
{
std::cerr << "*** Call this program as <./step-42 input.prm>"
<< std::endl;
return 1;
}
{
PlasticityContactProblem<3> problem(prm);
problem.run();
}
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
Results
The directory that contains this program also contains a number of input parameter files that can be used to create various different simulations. For example, running the program with the p1_adaptive.prm parameter file (using a ball as obstacle and the box as domain) on 16 cores produces output like this:
Using output directory 'p1adaptive/'
FE degree 1
transfer solution false
Cycle 0:
Number of active cells: 512
Number of degrees of freedom: 2187
Newton iteration 1
Updating active set...
Size of active set: 1
Assembling system...
Solving system...
Error: 173.076 -> 1.64265e-06 in 7 Bicgstab iterations.
Accepting Newton solution with residual: 1.64265e-06
Newton iteration 2
Updating active set...
Size of active set: 1
Assembling system...
Solving system...
Error: 57.3622 -> 3.23721e-07 in 8 Bicgstab iterations.
Accepting Newton solution with residual: 24.9028
Active set did not change!
Newton iteration 3
Updating active set...
Size of active set: 1
Assembling system...
Solving system...
Error: 24.9028 -> 9.94326e-08 in 7 Bicgstab iterations.
Residual of the non-contact part of the system: 1.63333
with a damping parameter alpha = 1
Active set did not change!
...
Newton iteration 6
Updating active set...
Size of active set: 1
Assembling system...
Solving system...
Error: 1.43188e-07 -> 3.56218e-16 in 8 Bicgstab iterations.
Residual of the non-contact part of the system: 4.298e-14
with a damping parameter alpha = 1
Active set did not change!
Writing graphical output... p1_adaptive/solution-00.pvtu
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 1.13s | |
| | | |
| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assembling | 6 | 0.463s | 41% |
| Graphical output | 1 | 0.0257s | 2.3% |
| Residual and lambda | 4 | 0.0754s | 6.7% |
| Setup | 1 | 0.227s | 20% |
| Setup: constraints | 1 | 0.0347s | 3.1% |
| Setup: distribute DoFs | 1 | 0.0441s | 3.9% |
| Setup:
matrix | 1 | 0.0119s | 1.1% |
| Setup: vectors | 1 | 0.00155s | 0.14% |
| Solve | 6 | 0.246s | 22% |
| Solve: iterate | 6 | 0.0631s | 5.6% |
| Solve: setup preconditioner | 6 | 0.167s | 15% |
| update active set | 6 | 0.0401s | 3.6% |
+---------------------------------+-----------+------------+------------+
Peak virtual memory used, resident in kB: 541884 77464
Contact force = 37.3058
...
Cycle 3:
Number of active cells: 14652
Number of degrees of freedom: 52497
Newton iteration 1
Updating active set...
Size of active set: 145
Assembling system...
Solving system...
Error: 296.309 -> 2.72484e-06 in 10 Bicgstab iterations.
Accepting Newton solution with residual: 2.72484e-06
...
Newton iteration 10
Updating active set...
Size of active set: 145
Assembling system...
Solving system...
Error: 2.71541e-07 -> 1.5428e-15 in 27 Bicgstab iterations.
Residual of the non-contact part of the system: 1.89261e-13
with a damping parameter alpha = 1
Active set did not change!
Writing graphical output... p1_adaptive/solution-03.pvtu
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 38.4s | |
| | | |
| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assembling | 10 | 22.5s | 58% |
| Graphical output | 1 | 0.327s | 0.85% |
| Residual and lambda | 9 | 3.75s | 9.8% |
| Setup | 1 | 4.83s | 13% |
| Setup: constraints | 1 | 0.578s | 1.5% |
| Setup: distribute DoFs | 1 | 0.71s | 1.8% |
| Setup:
matrix | 1 | 0.111s | 0.29% |
| Setup:
refine mesh | 1 | 4.83s | 13% |
| Setup: vectors | 1 | 0.00548s | 0.014% |
| Solve | 10 | 5.49s | 14% |
| Solve: iterate | 10 | 3.5s | 9.1% |
| Solve: setup preconditioner | 10 | 1.84s | 4.8% |
| update active set | 10 | 0.662s | 1.7% |
+---------------------------------+-----------+------------+------------+
Peak virtual memory used, resident in kB: 566052 105788
Contact force = 56.794
...
The tables at the end of each cycle show information about computing time (these numbers are of course specific to the machine on which this output was produced) and the number of calls of different parts of the program like assembly or calculating the residual, for the most recent mesh refinement cycle. Some of the numbers above can be improved by transferring the solution from one mesh to the next, an option we have not exercised here. Of course, you can also make the program run faster, especially on the later refinement cycles, by just using more processors: the accompanying paper shows good scaling to at least 1000 cores.
In a typical run, you can observe that for every refinement step, the active set - the contact points - are iterated out at first. After that the Newton method has only to resolve the plasticity. For the finer meshes, quadratic convergence can be observed for the last 4 or 5 Newton iterations.
We will not discuss here in all detail what happens with each of the input files. Rather, let us just show pictures of the solution (the left half of the domain is omitted if cells have zero quadrature points at which the plastic inequality is active):
The picture shows the adaptive refinement and as well how much a cell is plastified during the contact with the ball. Remember that we consider the norm of the deviator part of the stress in each quadrature point to see if there is elastic or plastic behavior. The blue color means that this cell contains only elastic quadrature points in contrast to the red cells in which all quadrature points are plastified. In the middle of the top surface - where the mesh is finest - a very close look shows the dimple caused by the obstacle. This is the result of the move_mesh() function. However, because the indentation of the obstacles we consider here is so small, it is hard to discern this effect; one could play with displacing vertices of the mesh by a multiple of the computed displacement.
Further discussion of results that can be obtained using this program is provided in the publication mentioned at the very top of this page.
Possibilities for extensions
There are, as always, multiple possibilities for extending this program. From an algorithmic perspective, this program goes about as far as one can at the time of writing, using the best available algorithms for the contact inequality, the plastic nonlinearity, and the linear solvers. However, there are things one would like to do with this program as far as more realistic situations are concerned:
-
Extend the program from a static to a quasi-static situation, perhaps by choosing a backward-Euler-scheme for the time discretization. Some theoretical results can be found in the PhD thesis by Jörg Frohne, FEM-Simulation der Umformtechnik metallischer Oberflächen im Mikrokosmos, University of Siegen, Germany, 2011.
-
It would also be an interesting advance to consider a contact problem with friction. In almost every mechanical process friction has a big influence. To model this situation, we have to take into account tangential stresses at the contact surface. Friction also adds another inequality to our problem since body and obstacle will typically stick together as long as the tangential stress does not exceed a certain limit, beyond which the two bodies slide past each other.
-
If we already simulate a frictional contact, the next step to consider is heat generation over the contact zone. The heat that is caused by friction between two bodies raises the temperature in the deformable body and entails an change of some material parameters.
-
It might be of interest to implement more accurate, problem-adapted error estimators for contact as well as for the plasticity.
The plain program
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/fe_field_function.h>
#include <fstream>
#include <iostream>
#include <sys/stat.h>
#include <cerrno>
namespace Step42
{
template <int dim>
class ConstitutiveLaw
{
public:
ConstitutiveLaw(const double E,
const double nu,
const double sigma_0,
const double gamma);
void set_sigma_0(double sigma_zero);
bool get_stress_strain_tensor(
void get_linearized_stress_strain_tensors(
private:
const double kappa;
double sigma_0;
const double gamma;
};
template <int dim>
ConstitutiveLaw<dim>::ConstitutiveLaw(double E,
double nu,
double sigma_0,
double gamma)
: kappa(
E / (3 * (1 - 2 * nu)))
, sigma_0(sigma_0)
, gamma(gamma)
, stress_strain_tensor_kappa(kappa *
, stress_strain_tensor_mu(
3.0))
{}
template <int dim>
void ConstitutiveLaw<dim>::set_sigma_0(double sigma_zero)
{
sigma_0 = sigma_zero;
}
template <int dim>
bool ConstitutiveLaw<dim>::get_stress_strain_tensor(
{
stress_tensor =
(stress_strain_tensor_kappa + stress_strain_tensor_mu) * strain_tensor;
const double deviator_stress_tensor_norm = deviator_stress_tensor.
norm();
stress_strain_tensor = stress_strain_tensor_mu;
if (deviator_stress_tensor_norm > sigma_0)
{
const double beta = sigma_0 / deviator_stress_tensor_norm;
stress_strain_tensor *= (gamma + (1 - gamma) * beta);
}
stress_strain_tensor += stress_strain_tensor_kappa;
return (deviator_stress_tensor_norm > sigma_0);
}
template <int dim>
void ConstitutiveLaw<dim>::get_linearized_stress_strain_tensors(
{
stress_tensor =
(stress_strain_tensor_kappa + stress_strain_tensor_mu) * strain_tensor;
stress_strain_tensor = stress_strain_tensor_mu;
stress_strain_tensor_linearized = stress_strain_tensor_mu;
const double deviator_stress_tensor_norm = deviator_stress_tensor.
norm();
if (deviator_stress_tensor_norm > sigma_0)
{
const double beta = sigma_0 / deviator_stress_tensor_norm;
stress_strain_tensor *= (gamma + (1 - gamma) * beta);
stress_strain_tensor_linearized *= (gamma + (1 - gamma) * beta);
deviator_stress_tensor /= deviator_stress_tensor_norm;
stress_strain_tensor_linearized -=
(1 - gamma) * beta * 2 *
mu *
}
stress_strain_tensor += stress_strain_tensor_kappa;
stress_strain_tensor_linearized += stress_strain_tensor_kappa;
}
namespace EquationData
{
template <int dim>
class BoundaryForce :
public Function<dim>
{
public:
BoundaryForce();
const unsigned int component = 0) const override;
};
template <int dim>
BoundaryForce<dim>::BoundaryForce()
{}
template <int dim>
double BoundaryForce<dim>::value(
const Point<dim> &,
const unsigned int) const
{
return 0.;
}
template <int dim>
void BoundaryForce<dim>::vector_value(
const Point<dim> &p,
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = BoundaryForce<dim>::value(p, c);
}
template <int dim>
class BoundaryValues :
public Function<dim>
{
public:
BoundaryValues();
const unsigned int component = 0) const override;
};
template <int dim>
BoundaryValues<dim>::BoundaryValues()
{}
template <int dim>
double BoundaryValues<dim>::value(
const Point<dim> &,
const unsigned int) const
{
return 0.;
}
template <int dim>
class SphereObstacle :
public Function<dim>
{
public:
SphereObstacle(const double z_surface);
const unsigned int component = 0) const override;
private:
const double z_surface;
};
template <int dim>
SphereObstacle<dim>::SphereObstacle(const double z_surface)
, z_surface(z_surface)
{}
template <int dim>
double SphereObstacle<dim>::value(
const Point<dim> & p,
const unsigned int component) const
{
if (component == 0)
return p(0);
else if (component == 1)
return p(1);
else if (component == 2)
{
if ((p(0) - 0.5) * (p(0) - 0.5) + (p(1) - 0.5) * (p(1) - 0.5) < 0.36)
return (-std::sqrt(0.36 - (p(0) - 0.5) * (p(0) - 0.5) -
(p(1) - 0.5) * (p(1) - 0.5)) +
z_surface + 0.59);
else
return 1000;
}
return 1e9;
}
template <int dim>
void SphereObstacle<dim>::vector_value(
const Point<dim> &p,
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = SphereObstacle<dim>::value(p, c);
}
template <int dim>
class BitmapFile
{
public:
BitmapFile(const std::string &name);
double get_value(const double x, const double y) const;
private:
std::vector<double> obstacle_data;
double hx, hy;
int nx, ny;
double get_pixel_value(const int i, const int j) const;
};
template <int dim>
BitmapFile<dim>::BitmapFile(const std::string &name)
: obstacle_data(0)
, hx(0)
, hy(0)
, nx(0)
, ny(0)
{
std::ifstream f(name);
ExcMessage(std::string(
"Can't read from file <") + name +
">!"));
std::string temp;
f >> temp >> nx >> ny;
for (int k = 0; k < nx * ny; ++k)
{
double val;
f >> val;
obstacle_data.push_back(val);
}
hx = 1.0 / (nx - 1);
hy = 1.0 / (ny - 1);
std::cout << "Read obstacle from file <" << name << ">" << std::endl
<< "Resolution of the scanned obstacle picture: " << nx
<< " x " << ny << std::endl;
}
template <int dim>
double BitmapFile<dim>::get_pixel_value(const int i, const int j) const
{
assert(i >= 0 && i < nx);
assert(j >= 0 && j < ny);
return obstacle_data[nx * (ny - 1 - j) + i];
}
template <int dim>
double BitmapFile<dim>::get_value(const double x, const double y) const
{
const int ix = std::min(std::max(static_cast<int>(x / hx), 0), nx - 2);
const int iy = std::min(std::max(static_cast<int>(y / hy), 0), ny - 2);
const double xi = std::min(std::max((x - ix * hx) / hx, 1.), 0.);
const double eta = std::min(std::max((y - iy * hy) / hy, 1.), 0.);
return ((1 - xi) * (1 - eta) * get_pixel_value(ix, iy) +
xi * (1 - eta) * get_pixel_value(ix + 1, iy) +
(1 - xi) * eta * get_pixel_value(ix, iy + 1) +
xi * eta * get_pixel_value(ix + 1, iy + 1));
}
template <int dim>
class ChineseObstacle :
public Function<dim>
{
public:
ChineseObstacle(const std::string &filename, const double z_surface);
const unsigned int component = 0) const override;
private:
const BitmapFile<dim> input_obstacle;
double z_surface;
};
template <int dim>
ChineseObstacle<dim>::ChineseObstacle(const std::string &filename,
const double z_surface)
, input_obstacle(filename)
, z_surface(z_surface)
{}
template <int dim>
double ChineseObstacle<dim>::value(
const Point<dim> & p,
const unsigned int component) const
{
if (component == 0)
return p(0);
if (component == 1)
return p(1);
else if (component == 2)
{
if (p(0) >= 0.0 && p(0) <= 1.0 && p(1) >= 0.0 && p(1) <= 1.0)
return z_surface + 0.999 - input_obstacle.get_value(p(0), p(1));
}
return 1e9;
}
template <int dim>
void ChineseObstacle<dim>::vector_value(
const Point<dim> &p,
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = ChineseObstacle<dim>::value(p, c);
}
}
template <int dim>
class PlasticityContactProblem
{
public:
private:
void make_grid();
void setup_system();
void compute_dirichlet_constraints();
void update_solution_and_constraints();
void
void assemble_newton_system(
void compute_nonlinear_residual(
void solve_newton_system();
void solve_newton();
void refine_grid();
void output_results(const std::string &filename_base);
void output_contact_force() const;
MPI_Comm mpi_communicator;
const unsigned int n_initial_global_refinements;
const unsigned int fe_degree;
Vector<float> fraction_of_plastic_q_points_per_cell;
const double e_modulus, nu, gamma, sigma_0;
ConstitutiveLaw<dim> constitutive_law;
const std::string base_mesh;
const std::shared_ptr<const Function<dim>> obstacle;
struct RefinementStrategy
{
enum value
{
refine_global,
refine_percentage,
refine_fix_dofs
};
};
typename RefinementStrategy::value refinement_strategy;
const bool transfer_solution;
std::string output_dir;
const unsigned int n_refinement_cycles;
unsigned int current_refinement_cycle;
};
template <int dim>
{
"polynomial degree",
"1",
"Polynomial degree of the FE_Q finite element space, typically 1 or 2.");
"2",
"Number of initial global mesh refinement steps before "
"the first computation.");
"refinement strategy",
"percentage",
"Mesh refinement strategy:\n"
" global: one global refinement\n"
" percentage: a fixed percentage of cells gets refined using the Kelly estimator.");
"5",
"Number of adaptive mesh refinement cycles to run.");
"obstacle",
"sphere",
"The name of the obstacle to use. This may either be 'sphere' if we should "
"use a spherical obstacle, or 'read from file' in which case the obstacle "
"will be read from a file named 'obstacle.pbm' that is supposed to be in "
"ASCII PBM format.");
"output directory",
"",
"Directory for output files (graphical output and benchmark "
"statistics). If empty, use the current directory.");
"transfer solution",
"false",
"Whether the solution should be used as a starting guess "
"for the next finer mesh. If false, then the iteration starts at "
"zero on every mesh.");
"box",
"Select the shape of the domain: 'box' or 'half sphere'");
}
template <int dim>
PlasticityContactProblem<dim>::PlasticityContactProblem(
: mpi_communicator(MPI_COMM_WORLD)
, computing_timer(MPI_COMM_WORLD,
pcout,
, n_initial_global_refinements(
prm.get_integer("number of initial refinements"))
, triangulation(mpi_communicator)
, fe_degree(prm.get_integer("polynomial degree"))
, dof_handler(triangulation)
, e_modulus(200000)
, nu(0.3)
, gamma(0.01)
, sigma_0(400.0)
, constitutive_law(e_modulus, nu, sigma_0, gamma)
, base_mesh(prm.get("base mesh"))
, obstacle(prm.get("obstacle") == "read from file" ?
new EquationData::ChineseObstacle<dim>(
"obstacle.pbm",
(base_mesh == "box" ? 1.0 : 0.5))) :
new EquationData::SphereObstacle<dim>(
base_mesh == "box" ? 1.0 : 0.5)))
, transfer_solution(prm.get_bool("transfer solution"))
, n_refinement_cycles(prm.get_integer("number of cycles"))
, current_refinement_cycle(0)
{
std::string strat = prm.
get(
"refinement strategy");
if (strat == "global")
refinement_strategy = RefinementStrategy::refine_global;
else if (strat == "percentage")
refinement_strategy = RefinementStrategy::refine_percentage;
else
output_dir = prm.
get(
"output directory");
if (output_dir != "" && *(output_dir.rbegin()) != '/')
output_dir += "/";
{
const int ierr = mkdir(output_dir.c_str(), 0777);
}
pcout << " Using output directory '" << output_dir << "'" << std::endl;
pcout << " FE degree " << fe_degree << std::endl;
pcout << " transfer solution " << (transfer_solution ? "true" : "false")
<< std::endl;
}
{
return {in(2), in(1), -in(0)};
}
template <int dim>
void PlasticityContactProblem<dim>::make_grid()
{
if (base_mesh == "half sphere")
{
const double radius = 0.8;
triangulation.reset_all_manifolds();
triangulation.set_manifold(0, manifold_description);
}
else
{
for (const auto &cell : triangulation.active_cell_iterators())
for (unsigned int face_no = 0;
face_no < GeometryInfo<dim>::faces_per_cell;
++face_no)
if (cell->face(face_no)->at_boundary())
{
if (std::fabs(cell->face(face_no)->center()[2] - p2[2]) < 1e-12)
cell->face(face_no)->set_boundary_id(1);
if (std::fabs(cell->face(face_no)->center()[0] - p1[0]) <
1e-12 ||
std::fabs(cell->face(face_no)->center()[0] - p2[0]) <
1e-12 ||
std::fabs(cell->face(face_no)->center()[1] - p1[1]) <
1e-12 ||
std::fabs(cell->face(face_no)->center()[1] - p2[1]) < 1e-12)
cell->face(face_no)->set_boundary_id(8);
if (std::fabs(cell->face(face_no)->center()[2] - p1[2]) < 1e-12)
cell->face(face_no)->set_boundary_id(6);
}
}
triangulation.refine_global(n_initial_global_refinements);
}
template <int dim>
void PlasticityContactProblem<dim>::setup_system()
{
{
dof_handler.distribute_dofs(fe);
locally_owned_dofs = dof_handler.locally_owned_dofs();
locally_relevant_dofs.
clear();
locally_relevant_dofs);
}
{
constraints_hanging_nodes.reinit(locally_relevant_dofs);
constraints_hanging_nodes);
constraints_hanging_nodes.close();
pcout << " Number of active cells: "
<< triangulation.n_global_active_cells() << std::endl
<< " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
compute_dirichlet_constraints();
}
{
solution.
reinit(locally_relevant_dofs, mpi_communicator);
newton_rhs.reinit(locally_owned_dofs, mpi_communicator);
newton_rhs_uncondensed.reinit(locally_owned_dofs, mpi_communicator);
diag_mass_matrix_vector.reinit(locally_owned_dofs, mpi_communicator);
fraction_of_plastic_q_points_per_cell.reinit(
triangulation.n_active_cells());
active_set.clear();
active_set.set_size(dof_handler.n_dofs());
}
{
mpi_communicator);
sp,
constraints_dirichlet_and_hanging_nodes,
false,
mpi_communicator));
sp.compress();
newton_matrix.reinit(sp);
assemble_mass_matrix_diagonal(mass_matrix);
const unsigned int start = (newton_rhs.local_range().first),
end = (newton_rhs.local_range().second);
for (unsigned int j = start; j < end; ++j)
diag_mass_matrix_vector(j) =
mass_matrix.diag_element(j);
}
}
template <int dim>
void PlasticityContactProblem<dim>::compute_dirichlet_constraints()
{
constraints_dirichlet_and_hanging_nodes.reinit(locally_relevant_dofs);
constraints_dirichlet_and_hanging_nodes.merge(constraints_hanging_nodes);
if (base_mesh == "box")
{
dof_handler,
6,
EquationData::BoundaryValues<dim>(),
constraints_dirichlet_and_hanging_nodes,
dof_handler,
8,
EquationData::BoundaryValues<dim>(),
constraints_dirichlet_and_hanging_nodes,
(fe.component_mask(x_displacement) |
fe.component_mask(y_displacement)));
}
else
dof_handler,
0,
EquationData::BoundaryValues<dim>(),
constraints_dirichlet_and_hanging_nodes,
constraints_dirichlet_and_hanging_nodes.close();
}
template <int dim>
void PlasticityContactProblem<dim>::assemble_mass_matrix_diagonal(
{
face_quadrature_formula,
const unsigned int n_face_q_points = face_quadrature_formula.size();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
cell->face(face)->boundary_id() == 1)
{
fe_values_face.reinit(cell, face);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
for (unsigned int i = 0; i < dofs_per_cell; ++i)
(fe_values_face[displacement].value(i, q_point) *
fe_values_face[displacement].value(i, q_point) *
fe_values_face.JxW(q_point));
cell->get_dof_indices(local_dof_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
local_dof_indices[i],
}
}
template <int dim>
void PlasticityContactProblem<dim>::update_solution_and_constraints()
{
std::vector<bool> dof_touched(dof_handler.n_dofs(), false);
mpi_communicator);
distributed_solution = solution;
mpi_communicator);
lambda = newton_rhs_uncondensed;
locally_relevant_dofs, mpi_communicator);
diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
all_constraints.reinit(locally_relevant_dofs);
active_set.clear();
Quadrature<dim - 1> face_quadrature(fe.get_unit_face_support_points());
face_quadrature,
const unsigned int dofs_per_face = fe.dofs_per_face;
const unsigned int n_face_q_points = face_quadrature.size();
std::vector<types::global_dof_index> dof_indices(dofs_per_face);
for (const auto &cell : dof_handler.active_cell_iterators())
if (!cell->is_artificial())
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
cell->face(face)->boundary_id() == 1)
{
fe_values_face.reinit(cell, face);
cell->face(face)->get_dof_indices(dof_indices);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
const unsigned int component =
fe.face_system_to_component_index(q_point).first;
const unsigned int index_z = dof_indices[q_point];
if ((component == 2) && (dof_touched[index_z] == false))
{
dof_touched[index_z] = true;
fe_values_face.quadrature_point(q_point);
const double obstacle_value =
obstacle->value(this_support_point, 2);
const double solution_here = solution(index_z);
const double undeformed_gap =
obstacle_value - this_support_point(2);
const double c = 100.0 * e_modulus;
if ((lambda(index_z) /
diag_mass_matrix_vector_relevant(index_z) +
c * (solution_here - undeformed_gap) >
0) &&
!constraints_hanging_nodes.is_constrained(index_z))
{
all_constraints.add_line(index_z);
all_constraints.set_inhomogeneity(index_z,
undeformed_gap);
distributed_solution(index_z) = undeformed_gap;
active_set.add_index(index_z);
}
}
}
}
solution = distributed_solution;
all_constraints.close();
all_constraints.merge(constraints_dirichlet_and_hanging_nodes);
pcout << " Size of active set: "
mpi_communicator)
<< std::endl;
}
template <int dim>
void PlasticityContactProblem<dim>::assemble_newton_system(
{
QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
quadrature_formula,
face_quadrature_formula,
const unsigned int n_q_points = quadrature_formula.
size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
const EquationData::BoundaryForce<dim> boundary_force;
std::vector<Vector<double>> boundary_force_values(n_face_q_points,
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_rhs = 0;
std::vector<SymmetricTensor<2, dim>> strain_tensor(n_q_points);
fe_values[displacement].get_function_symmetric_gradients(
linearization_point, strain_tensor);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
constitutive_law.get_linearized_stress_strain_tensors(
strain_tensor[q_point],
stress_strain_tensor_linearized,
stress_strain_tensor);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
stress_strain_tensor_linearized *
fe_values[displacement].symmetric_gradient(i, q_point);
for (unsigned int j = 0; j < dofs_per_cell; ++j)
(stress_phi_i *
fe_values[displacement].symmetric_gradient(j, q_point) *
cell_rhs(i) +=
((stress_phi_i -
stress_strain_tensor *
fe_values[displacement].symmetric_gradient(i,
q_point)) *
strain_tensor[q_point] * fe_values.
JxW(q_point));
}
}
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
cell->face(face)->boundary_id() == 1)
{
fe_values_face.
reinit(cell, face);
boundary_force.vector_value_list(
fe_values_face.get_quadrature_points(),
boundary_force_values);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
rhs_values[2] = boundary_force_values[q_point][2];
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) +=
(fe_values_face[displacement].value(i, q_point) *
rhs_values * fe_values_face.JxW(q_point));
}
}
cell->get_dof_indices(local_dof_indices);
all_constraints.distribute_local_to_global(cell_matrix,
cell_rhs,
local_dof_indices,
newton_matrix,
newton_rhs,
true);
}
}
template <int dim>
void PlasticityContactProblem<dim>::compute_nonlinear_residual(
{
QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
quadrature_formula,
face_quadrature_formula,
const unsigned int n_q_points = quadrature_formula.
size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
const EquationData::BoundaryForce<dim> boundary_force;
std::vector<Vector<double>> boundary_force_values(n_face_q_points,
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
newton_rhs = 0;
newton_rhs_uncondensed = 0;
fraction_of_plastic_q_points_per_cell = 0;
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_rhs = 0;
std::vector<SymmetricTensor<2, dim>> strain_tensors(n_q_points);
fe_values[displacement].get_function_symmetric_gradients(
linearization_point, strain_tensors);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const bool q_point_is_plastic =
constitutive_law.get_stress_strain_tensor(
strain_tensors[q_point], stress_strain_tensor);
if (q_point_is_plastic)
++fraction_of_plastic_q_points_per_cell(
cell->active_cell_index());
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
cell_rhs(i) -=
(strain_tensors[q_point] * stress_strain_tensor *
fe_values[displacement].symmetric_gradient(i, q_point) *
rhs_values = 0;
cell_rhs(i) += (fe_values[displacement].value(i, q_point) *
rhs_values * fe_values.
JxW(q_point));
}
}
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
cell->face(face)->boundary_id() == 1)
{
fe_values_face.
reinit(cell, face);
boundary_force.vector_value_list(
fe_values_face.get_quadrature_points(),
boundary_force_values);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
{
rhs_values[2] = boundary_force_values[q_point][2];
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) +=
(fe_values_face[displacement].value(i, q_point) *
rhs_values * fe_values_face.JxW(q_point));
}
}
cell->get_dof_indices(local_dof_indices);
constraints_dirichlet_and_hanging_nodes.distribute_local_to_global(
cell_rhs, local_dof_indices, newton_rhs);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
newton_rhs_uncondensed(local_dof_indices[i]) += cell_rhs(i);
}
fraction_of_plastic_q_points_per_cell /= quadrature_formula.
size();
}
template <int dim>
void PlasticityContactProblem<dim>::solve_newton_system()
{
mpi_communicator);
distributed_solution = solution;
constraints_hanging_nodes.set_zero(distributed_solution);
constraints_hanging_nodes.set_zero(newton_rhs);
{
std::vector<std::vector<bool>> constant_modes;
constant_modes);
preconditioner.
initialize(newton_matrix, additional_data);
}
{
const double relative_accuracy = 1
e-8;
const double solver_tolerance =
relative_accuracy *
newton_matrix.residual(tmp, distributed_solution, newton_rhs);
SolverControl solver_control(newton_matrix.m(), solver_tolerance);
solver.solve(newton_matrix,
distributed_solution,
newton_rhs,
preconditioner);
pcout << " Error: " << solver_control.initial_value() << " -> "
<< solver_control.last_value() << " in "
<< solver_control.last_step() << " Bicgstab iterations."
<< std::endl;
}
all_constraints.distribute(distributed_solution);
solution = distributed_solution;
}
template <int dim>
void PlasticityContactProblem<dim>::solve_newton()
{
mpi_communicator);
mpi_communicator);
mpi_communicator);
locally_relevant_dofs, mpi_communicator);
mpi_communicator);
double residual_norm;
double previous_residual_norm = -std::numeric_limits<double>::max();
const double correct_sigma = sigma_0;
for (unsigned int newton_step = 1; newton_step <= 100; ++newton_step)
{
if (newton_step == 1 &&
((transfer_solution && current_refinement_cycle == 0) ||
!transfer_solution))
constitutive_law.set_sigma_0(1e+10);
else if (newton_step == 2 || current_refinement_cycle > 0 ||
!transfer_solution)
constitutive_law.set_sigma_0(correct_sigma);
pcout << " " << std::endl;
pcout << " Newton iteration " << newton_step << std::endl;
pcout << " Updating active set..." << std::endl;
{
update_solution_and_constraints();
}
pcout << " Assembling system... " << std::endl;
newton_matrix = 0;
newton_rhs = 0;
assemble_newton_system(solution);
pcout << " Solving system... " << std::endl;
solve_newton_system();
if ((newton_step == 1) ||
(transfer_solution && newton_step == 2 &&
current_refinement_cycle == 0) ||
(!transfer_solution && newton_step == 2))
{
compute_nonlinear_residual(solution);
old_solution = solution;
residual = newton_rhs;
const unsigned int start_res = (residual.local_range().first),
end_res = (residual.local_range().second);
for (unsigned int n = start_res; n < end_res; ++n)
if (all_constraints.is_inhomogeneously_constrained(n))
residual(n) = 0;
residual_norm = residual.l2_norm();
pcout << " Accepting Newton solution with residual: "
<< residual_norm << std::endl;
}
else
{
for (unsigned int i = 0; i < 5; ++i)
{
distributed_solution = solution;
const double alpha = std::pow(0.5, static_cast<double>(i));
tmp_vector = old_solution;
tmp_vector.sadd(1 - alpha, alpha, distributed_solution);
locally_relevant_tmp_vector = tmp_vector;
compute_nonlinear_residual(locally_relevant_tmp_vector);
residual = newton_rhs;
const unsigned int start_res = (residual.local_range().first),
end_res = (residual.local_range().second);
for (unsigned int n = start_res; n < end_res; ++n)
if (all_constraints.is_inhomogeneously_constrained(n))
residual(n) = 0;
residual_norm = residual.l2_norm();
pcout
<< " Residual of the non-contact part of the system: "
<< residual_norm << std::endl
<< " with a damping parameter alpha = " << alpha
<< std::endl;
if (residual_norm < previous_residual_norm)
break;
}
solution = tmp_vector;
old_solution = solution;
}
previous_residual_norm = residual_norm;
mpi_communicator) == 0)
{
pcout << " Active set did not change!" << std::endl;
if (residual_norm < 1e-10)
break;
}
old_active_set = active_set;
}
}
template <int dim>
void PlasticityContactProblem<dim>::refine_grid()
{
if (refinement_strategy == RefinementStrategy::refine_global)
{
triangulation.begin_active();
cell != triangulation.end();
++cell)
if (cell->is_locally_owned())
cell->set_refine_flag();
}
else
{
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
dof_handler,
solution,
estimated_error_per_cell);
parallel::distributed::GridRefinement ::refine_and_coarsen_fixed_number(
triangulation, estimated_error_per_cell, 0.3, 0.03);
}
triangulation.prepare_coarsening_and_refinement();
solution_transfer(dof_handler);
if (transfer_solution)
solution_transfer.prepare_for_coarsening_and_refinement(solution);
triangulation.execute_coarsening_and_refinement();
setup_system();
if (transfer_solution)
{
mpi_communicator);
solution_transfer.interpolate(distributed_solution);
constraints_hanging_nodes.distribute(distributed_solution);
solution = distributed_solution;
compute_nonlinear_residual(solution);
}
}
template <int dim>
void PlasticityContactProblem<dim>::move_mesh(
{
std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
if (vertex_touched[cell->vertex_index(v)] == false)
{
vertex_touched[cell->vertex_index(v)] = true;
for (
unsigned int d = 0;
d < dim; ++
d)
vertex_displacement[d] =
displacement(cell->vertex_dof_index(v, d));
cell->vertex(v) += vertex_displacement;
}
}
template <int dim>
void PlasticityContactProblem<dim>::output_results(
const std::string &filename_base)
{
pcout << " Writing graphical output... " << std::flush;
move_mesh(solution);
mpi_communicator);
const unsigned int start_res = (newton_rhs_uncondensed.local_range().first),
end_res = (newton_rhs_uncondensed.local_range().second);
for (unsigned int n = start_res; n < end_res; ++n)
if (all_constraints.is_inhomogeneously_constrained(n))
distributed_lambda(n) =
newton_rhs_uncondensed(n) / diag_mass_matrix_vector(n);
constraints_hanging_nodes.distribute(distributed_lambda);
mpi_communicator);
lambda = distributed_lambda;
locally_owned_dofs, mpi_communicator);
distributed_active_set_vector = 0.;
for (const auto index : active_set)
distributed_active_set_vector[index] = 1.;
mpi_communicator);
active_set_vector = distributed_active_set_vector;
const std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(
std::vector<std::string>(dim, "displacement"),
data_component_interpretation);
std::vector<std::string>(dim, "contact_force"),
data_component_interpretation);
std::vector<std::string>(dim, "active_set"),
data_component_interpretation);
Vector<float> subdomain(triangulation.n_active_cells());
for (unsigned int i = 0; i < subdomain.size(); ++i)
subdomain(i) = triangulation.locally_owned_subdomain();
"fraction_of_plastic_q_points");
const std::string filename =
(output_dir + filename_base + "-" +
std::ofstream output_vtu((filename + ".vtu"));
pcout << output_dir + filename_base << ".pvtu" << std::endl;
{
std::vector<std::string> filenames;
for (unsigned int i = 0;
++i)
filenames.push_back(filename_base + "-" +
std::ofstream pvtu_master_output(
(output_dir + filename_base + ".pvtu"));
std::ofstream visit_master_output(
(output_dir + filename_base + ".visit"));
}
tmp *= -1;
move_mesh(tmp);
}
template <int dim>
void PlasticityContactProblem<dim>::output_contact_force() const
{
mpi_communicator);
const unsigned int start_res = (newton_rhs_uncondensed.local_range().first),
end_res = (newton_rhs_uncondensed.local_range().second);
for (unsigned int n = start_res; n < end_res; ++n)
if (all_constraints.is_inhomogeneously_constrained(n))
distributed_lambda(n) =
newton_rhs_uncondensed(n) / diag_mass_matrix_vector(n);
else
distributed_lambda(n) = 0;
constraints_hanging_nodes.distribute(distributed_lambda);
mpi_communicator);
lambda = distributed_lambda;
double contact_force = 0.0;
QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
face_quadrature_formula,
const unsigned int n_face_q_points = face_quadrature_formula.size();
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
if (cell->face(face)->at_boundary() &&
cell->face(face)->boundary_id() == 1)
{
fe_values_face.reinit(cell, face);
std::vector<Tensor<1, dim>> lambda_values(n_face_q_points);
fe_values_face[displacement].get_function_values(lambda,
lambda_values);
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
contact_force +=
lambda_values[q_point][2] * fe_values_face.JxW(q_point);
}
pcout << "Contact force = " << contact_force << std::endl;
}
template <int dim>
void PlasticityContactProblem<dim>::run()
{
computing_timer.reset();
for (; current_refinement_cycle < n_refinement_cycles;
++current_refinement_cycle)
{
{
pcout << std::endl;
pcout << "Cycle " << current_refinement_cycle << ':' << std::endl;
if (current_refinement_cycle == 0)
{
make_grid();
setup_system();
}
else
{
refine_grid();
}
}
solve_newton();
output_results(std::string("solution-") +
computing_timer.print_summary();
computing_timer.reset();
pcout <<
"Peak virtual memory used, resident in kB: " << stats.
VmSize <<
" " << stats.
VmRSS << std::endl;
if (base_mesh == "box")
output_contact_force();
}
}
}
int main(int argc, char *argv[])
{
using namespace Step42;
try
{
PlasticityContactProblem<3>::declare_parameters(prm);
if (argc != 2)
{
std::cerr << "*** Call this program as <./step-42 input.prm>"
<< std::endl;
return 1;
}
{
PlasticityContactProblem<3> problem(prm);
problem.run();
}
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
1.8.14