<br>
<i>This program was contributed by Jörg Frohne (University of Siegen,
-Germany) while on a long-term visit to Texas A&M University.
+Germany) while on a long-term visit to Texas A&M University, with significant
+contributions by Timo Heister and Wolfgang Bangerth.
<br>
+<br>
+The code described here provides the basis for the numerical experiments shown
+in the following paper:
+<br>
+ J. Frohne, T. Heister, W. Bangerth: <b>Efficient numerical methods for the large-scale, parallel
+ solution of elastoplastic contact problems</b>. Submitted, 2013.
+</em>
</i>
<a name="Intro"></a>
<h3>Introduction</h3>
-This example is an extension of step-41, considering a contact problem with an
+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.
-That means that we have to take care of an additional nonlinearity: the
-material behavior. Since we consider a three dimensional problem here, a
-separate difference to step-41 is that the contact area is at the boundary of
+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 we 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):
+
+@image html step-42-CellConstitutionLi2.png
+
+
+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 because of the adaptive mesh
-in both the handling of the linear
-system as well as of the inequality constraints; in the latter case, we will
+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.
-
-@f{huge}
-{distributed}
-@f}
+the parallelization. Additional pointers for parallelization can be found in
+step-32.
<h3>Classical formulation</h3>
The classical formulation of the problem possesses the following form:
@f{align*}
- \varepsilon(u) &= A\sigma + \lambda & &\quad\text{in } \Omega,\\
- \lambda(\tau - \sigma) &\geq 0\quad\forall\tau\text{ with
- }\mathcal{F}(\tau)\leq 0 & &\quad\text{in } \Omega,\\
- -\textrm{\textrm{div}}\ \sigma &= f & &\quad\text{in } \Omega,\\
- u(\mathbf x) &= 0 & &\quad\text{on }\Gamma_D,\\
- \sigma_t(u) &= 0,\quad\sigma_n(u)\leq 0 & &\quad\text{on }\Gamma_C,\\
-\sigma_n(u)(u_n - g) &= 0,\quad u_n(\mathbf x) - g(\mathbf x) \leq 0 & &\quad\text{on } \Gamma_C
+ \varepsilon(\mathbf u) &= A\sigma + \varepsilon^p & &\quad\text{in } \Omega,\\
+ -\textrm{\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 \cdot \mathbf n - [\mathbf n \cdot(\sigma \cdot \mathbf n)]\mathbf n &= 0,
+ \quad \mathbf n \cdot (\sigma \cdot
+ \mathbf n) \leq 0 & &\quad\text{on }\Gamma_C,\\
+ (\mathbf n \cdot (\sigma \cdot
+ \mathbf n))(\mathbf n \cdot \mathbf u - g) &= 0,\quad \mathbf n
+ \cdot \mathbf u - g \leq 0 & &\quad\text{on } \Gamma_C.
@f}
-with $u\in H^2(\Omega),\Omega\subset\mathbb{R}^3$. The vector valued
-function $u$ denotes the displacement in the deformable body. The first two lines describe the
-elasto-plastic material behavior. Therein the equation shows the
-strain of the deformation $\varepsilon (u)$ as the additive decomposition of the
-elastic part $A\sigma$ and the plastic part $\lambda$. $A$ is defined as the compliance tensor of fourth order which contains some material constants and $\sigma$ as the
-symmetric stress tensor of second order. So we have to consider
-the inequality in the second row in a pointwise sense where $\lambda(\tau -
-\sigma)$ is the inner product of two symmetric tensors of second order. Furthermore we have to
-distinguish two cases.
-
-The continuous and convex function $\mathcal{F}$ denotes the von Mises flow function
-@f{gather*}\mathcal{F}(\tau) = \vert\tau^D\vert - \sigma_0�\quad\text{with}\quad \tau^D
-= \tau - \dfrac{1}{3}tr(\tau)I,@f}
-$\sigma_0$ as yield stress and $\vert .\vert$ as the Frobenius norm. If there
-are no plastic deformations in a particular point - that is $\lambda=0$ - this yields $\vert\sigma^D\vert <
-\sigma_0$ and otherwise if $\lambda > 0$ it follows that $\vert\sigma^D\vert = \sigma_0$.
-That means if the stress is smaller than the yield stress there are only elastic
-deformations in that point.
-
-To consider it the other way around if the deviator stress $\sigma^D$ is in a
-norm bigger than the yield stress then $\sigma^D$ has to be projected back to the yield surface and there are plastic deformations which means $\lambda$
-would be positive for that particular point. We refer that the stresses are
-computed by Hooke's law for isotropic materials. You can find the description at the end of section 3. Else if the norm of the deviator stress tensor is smaller or equal the yield stress then $\lambda$ is zero and there are no plastic deformations in
-that point.
-
-There the index $D$ denotes the deviator part of for example the stress where
-$tr(.)$ is the trace of a tensor. The definition shows an additive decomposition
-of the stress $\sigma$ into a hydrostatic part (or volumetric part) $\dfrac{1}{3}tr(\tau)I$ and the deviator
-part $\sigma^D$. For metal the deviator stress composes the main indicator for
-plastic deformations.
-
-The third equation is called equilibrium condition with a force of volume
-density $f$ which we will neglect in our example.
-The boundary of $\Omega$ separates as follows $\Gamma=\Gamma_D\bigcup\Gamma_C$ and $\Gamma_D\bigcap\Gamma_C=\emptyset$.
-At the boundary $\Gamma_D$ we have zero Dirichlet conditions. $\Gamma_C$ denotes the potential contact boundary.
-
-The last two lines describe the so-called Signorini contact conditions. If there is no contact the normal stress
-@f{gather*} \sigma_n = \sigma n\cdot n@f}
-is zero with the outward normal $n$. If there is contact ($u_n = g$) the tangential stress $\sigma_t = \sigma\cdot n - \sigma_n n$
-vanishes, because we consider a frictionless situation and the normal stress is
-negative. The gap $g$ comes with the start configuration of the obstacle and the
-deformable body.
-
-
-<h3>Derivation of the variational inequality</h3>
-
-As a starting point to derive the equations above, let us imagine that we want
-to minimise an energy functional:
-@f{gather*}E(\tau) := \dfrac{1}{2}\int\limits_{\Omega}\tau A \tau d\tau,\quad \tau\in \Pi W^{\textrm{div}}@f}
-with
-@f{gather*}W^{\textrm{div}}:=\lbrace \tau\in
-L^2(\Omega,\mathbb{R}^{\textrm{dim}\times\textrm{dim}}_{\textrm{sym}}):\textrm{div}(\tau)\in L^2(\Omega,\mathbb{R}^{\textrm{dim}})\rbrace@f} and
-@f{gather*}\Pi \Sigma:=\lbrace \tau\in \Sigma, \mathcal{F}(\tau)\leq 0\rbrace@f}
-as the set of admissible stresses which is defined
-by a continuous, convex flow function $\mathcal{F}$.
-
-With the goal of deriving the dual formulation of the minimisation
-problem, we define a Lagrange function:
-@f{gather*}L(\tau,\varphi) := E(\tau) + (\varphi, \textrm{div}(\tau)),\quad \lbrace\tau,\varphi\rbrace\in\Pi W^{\textrm{div}}\times V^+@f}
-with
-@f{gather*}V^+ := \lbrace u\in V: u_n\leq g \text{ on } \Gamma_C \rbrace@f}
-@f{gather*}V:=\left[ H_0^1 \right]^{\textrm{dim}}:=\lbrace u\in \left[H^1(\Omega)\right]^{\textrm{dim}}: u
-= 0 \text{ on } \Gamma_D\rbrace@f}
-By building the Fréchet derivatives of $L$ for both components we obtain the
-dual formulation for the stationary case which is known as <i>Hencky-Type-Model</i>:\\
-Find a pair $\lbrace\sigma,u\rbrace\in \Pi W\times V^+$ with
-@f{gather*}\left(A\sigma,\tau - \sigma\right) + \left(u, \textrm{div}(\tau) - \textrm{div}(\sigma)\right) \geq 0,\quad \forall \tau\in \Pi W^{\textrm{div}}@f}
-@f{gather*}-\left(\textrm{div}(\sigma),\varphi - u\right) \geq 0,\quad \forall \varphi\in V^+.@f}
-By integrating by parts and multiplying the first inequality by the elastic
-tensor $C=A^{-1}$ we achieve the primal-mixed version of our problem:
-Find a pair $\lbrace\sigma,u\rbrace\in \Pi W\times V^+$ with
-@f{gather*}\left(\sigma,\tau - \sigma\right) - \left(C\varepsilon(u), \tau - \sigma\right) \geq 0,\quad \forall \tau\in \Pi W@f}
-@f{gather*}\left(\sigma,\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in V^+.@f}
-Therein $\varepsilon$ denotes the linearised deformation tensor with $\varepsilon(u) := \dfrac{1}{2}\left(\nabla u + \nabla u^T\right)$ for small deformations.
-
-Most materials - especially metals - have the property that they show some hardening effects during the forming process.
+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 <i>yield function</i> 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) \cdot
+ \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 \cdot \mathbf n - \mathbf \sigma_n \mathbf n
+= \sigma \cdot \mathbf n - [\mathbf n \cdot(\sigma \cdot \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}(\tau,\eta) = \vert\tau^D\vert - (\sigma_0 + \gamma^{\text{iso}}\eta)$ where
-$\eta$ is the norm of the plastic strain $\eta = \vert \varepsilon -
-A\sigma\vert$.
-It can be considered by establishing an additional term in our primal-mixed formulation:
-Find a pair $\lbrace(\sigma,\xi),u\rbrace\in \Pi (W\times L^2(\Omega,\mathbb{R}))\times V^+$ with
-@f{gather*}\left(\sigma,\tau - \sigma\right) - \left(C\varepsilon(u), \tau - \sigma\right) + \gamma^{\text{iso}}\left( \xi, \eta - \xi\right) \geq 0,\quad \forall (\tau,\eta)\in \Pi (W,L^2(\Omega,\mathbb{R}))@f}
-@f{gather*}\left(\sigma,\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in V^+,@f}
-with the hardening parameter $\gamma^{\text{iso}} > 0$.
-
-Now we want to derive a primal problem which only depends on the displacement $u$. For that purpose we
-set $\eta = \xi$ and eliminate the stress $\sigma$ by applying the projection
-theorem (see Grossmann, Roos: Numerical Treatment of Partial Differential
-Equations, Springer-Verlag Berlin Heidelberg, 2007 and Frohne: FEM-Simulation
-der Umformtechnik metallischer Oberflächen im Mikrokosmos, Ph.D. thesis,
-University of Siegen, Germany, 2011) on
-@f{gather*}\left(\sigma - C\varepsilon(u), \tau - \sigma\right) \geq 0,\quad \forall \tau\in \Pi W,@f}
-which yields with the second inequality:\\
-Find the displacement $u\in V^+$ with
-@f{gather*}\left(P_{\Pi}(C\varepsilon(u)),\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in V^+,@f}
-with the projection:
-@f{gather*}P_{\Pi}(\tau):=\begin{cases}
- \tau, & \text{if }\vert\tau^D\vert \leq \sigma_0 + \gamma^{\text{iso}}\xi,\\
- \hat\alpha\dfrac{\tau^D}{\vert\tau^D\vert} + \dfrac{1}{3}tr(\tau), & \text{if }\vert\tau^D\vert > \sigma_0 + \gamma^{\text{iso}}\xi,
- \end{cases}@f}
-with the radius
-@f{gather*}\hat\alpha := \sigma_0 + \gamma^{\text{iso}}\xi .@f}
-With the relation $\xi = \vert\varepsilon(u) - A\sigma\vert$ it is possible to eliminate $\xi$ inside the projection $P_{\Pi}$:\\
-@f{gather*}P_{\Pi}(\tau):=\begin{cases}
- \tau, & \text{if }\vert\tau^D\vert \leq \sigma_0,\\
- \alpha\dfrac{\tau^D}{\vert\tau^D\vert} + \dfrac{1}{3}tr(\tau), & \text{if }\vert\tau^D\vert > \sigma_0,
- \end{cases}@f}
-@f{gather*}\alpha := \sigma_0 + \dfrac{\gamma^{\text{iso}}}{2\mu+\gamma^{\text{iso}}}\left(\vert\tau^D\vert - \sigma_0\right) ,@f}
-with a further material parameter $\mu>0$ called shear modulus. We refer that
-this only possible for isotropic plasticity.
-
-To make things a bit easier from now on we denote
-@f{gather*}\gamma := \dfrac{\gamma^{\text{iso}}}{2\mu +
-\gamma^{\text{iso}}}\in[0,1)\text{ with }\gamma^{\text{iso}}\in[0,\infty),@f}
-@f{gather*}\beta :=\dfrac{\sigma_0}{\vert\tau^D\vert}.@f}
-If $\gamma^{\text{iso}}$ tends to zero $\gamma$ tends also to zero. And if $\gamma^{\text{iso}}$ tends to
-infinity $\gamma$ tends to one. This allows us to reformulate our problem as
-follows
-@f{gather*}P_{\Pi}(\tau):=\begin{cases}
- \tau, & \text{if }\vert\tau^D\vert \leq \sigma_0,\\
- \gamma\tau^D + (1-\gamma)\beta\tau^D
- + \dfrac{1}{3}tr(\tau), & \text{if }\vert\tau^D\vert >
- \sigma_0, \end{cases}.@f}
-For further details e.g., see Suttmeier: On Plasticity with Hardening:
-An Adaptive Finite Element Discretisation, International Mathematical Forum, 5,
-2010, no. 52, 2591-2601.
-
-So what we do is to calculate the stresses by using Hooke's law for linear elastic, isotropic materials
-@f{gather*}\sigma = C \varepsilon(u) = 2\mu \varepsilon^D(u) + \kappa tr(\varepsilon(u))I = \left[2\mu\left(\mathbb{I} -\dfrac{1}{3} I\otimes I\right) + \kappa I\otimes I\right]\varepsilon(u)@f}
-with the material parameter $\kappa>0$ (bulk modulus). The variables $I$ and
-$\mathbb{I}$ denote the identity tensors of second and forth order. In that
-notation $2\mu \varepsilon^D(u)$ is the deviatoric part and $\kappa
-tr(\varepsilon(u))$ the volumetric part of the stress tensor.
-
-In the next step we test in a pointwise sense where the deviator part of the
-stress in a norm is bigger than the yield stress. If there are such points we
-project the deviator stress in those points back to the yield surface. Methods of this kind are called projections algorithm or radial-return-algorithm.
-
-Now we have a primal formulation of our elasto-plastic contact problem which only depends on the displacement $u$.
-It consists of a nonlinear variational inequality and has a unique solution as
-it satisfies the theorem of Lions and Stampaccia. A proof can be found in
-Rodrigues: Obstacle Problems in Mathematical Physics, North-Holland, Amsterdam,
-1987.
-
-To handle the nonlinearity of the constitutive law we use a Newton method and to deal with the contact we apply an
-active set method like in step-41. To be more concrete we combine both methods to an inexact semi smooth Newton
-method - inexact since we use an iterative solver for the linearised problems in each Newton step.
-
-
-<h3>Linearisation of the constitutive law for the Newton method</h3>
-
-For the Newton method we have to linearise the following semi-linearform
-@f{gather*}a(\psi;\varphi) := \left(P_{\Pi}(C\varepsilon(\psi)),\varepsilon(\varphi)\right).@f}
-Because we have to find the solution $u$ in the convex set $V^+$, we have to
-apply an SQP-method (SQP: sequential quadratic programming). That means we have
-to solve a minimisation problem for a known $u^i$ in every SQP-step of the form
-@f{eqnarray*}
- & & a(u^{i};u^{i+1} - u^i) + \dfrac{1}{2}a'(u^i;u^{i+1} - u^i,u^{i+1} - u^i)\\
- &=& a(u^i;u^{i+1}) - a(u^i;u^i) +\\
- & & \dfrac{1}{2}\left( a'(u^i;u^{i+1},u^{i+1}) - 2a'(u^i;u^i,u^{i+1}) - a'(u^i;u^i,u^i)\right)\\
- &\rightarrow& \textrm{min},\quad u^{i+1}\in V^+.
-@f}
-Neglecting the constant terms $ a(u^i;u^i)$ and $ a'(u^i;u^i,u^i)$ we obtain the
-following minimisation problem @f{gather*}\dfrac{1}{2} a'(u^i;u^{i+1},u^{i+1}) - F(u^i)\rightarrow \textrm{min},\quad u^{i+1}\in V^+@f} with
-@f{gather*}F(\varphi) := \left(a'(\varphi;\varphi,u^{i+1}) - a(\varphi;u^{i+1}) \right).@f}
-In the case of our constitutive law the Fréchet derivative of the
-semi-linearform $a(.;.)$ at the point $u^i$ is
-
-@f{gather*}a'(u^i;\psi,\varphi) =
-(I(x)\varepsilon(\psi),\varepsilon(\varphi)),\quad x\in\Omega,@f}
-@f{gather*}
-I(x) := \begin{cases}
-C_{\mu} + C_{\kappa}, &
-\quad \vert \tau^D \vert \leq \sigma_0\\
-\gamma C_{\mu} + (1-\gamma)\beta\left(C_{\mu} -
-2\mu\dfrac{\tau^D\otimes\tau^D}{\vert\tau^D\vert^2}\right) + C_{\kappa}, &\quad
-\vert \tau^D \vert > \sigma_0
-\end{cases}
+$\mathcal{F}(\sigma,\varepsilon^p) = \vert\sigma^D\vert - (\sigma_0 +
+\gamma^{\text{iso}}|\varepsilon^p|)$.
+
+
+<h3>Reformulation as a variational inequality</h3>
+
+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 <i>semismooth</i> 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
+@f{align*}
+\left(P_{\Pi}(C\varepsilon(\mathbf u)),\varepsilon(\varphi) - \varepsilon(\mathbf u)\right) \geq 0,\quad \forall \varphi\in V^+.
@f}
-with
-@f{gather*}C_{\mu} := 2\mu\left(\mathbb{I} - \dfrac{1}{3} I\otimes
-I\right)\quad\text{(shear part of the stress strain tensor)},@f}
-@f{gather*}C_{\kappa} := \kappa I\otimes I\quad\text{(bulk part of the stress strain
-tensor)},@f}
-@f{gather*}\tau^D := C\varepsilon^D(u^i).@f}
-Remark that $a(.;.)$ is not differentiable in the common sense but it is
-slantly differentiable like the function for the contact problem and again we refer to
-Hintermueller, Ito, Kunisch: The primal-dual active set strategy as a semismooth Newton method, SIAM J. OPTIM., 2003, Vol. 13, No. 3, pp. 865-888.
-Again the first case is for elastic and the second for plastic deformation.
-
-
-<h3>Formulation as a saddle point problem</h3>
-
-Just as in step-41 we compose a saddle point problem out of the minimisation
-problem. Again we do so to gain a formulation that allows us to solve a linear
-system of equations finally.
-
-We introduce a Lagrange multiplier $\lambda$ and the convex cone $K\subset W'$,
-$W'$ dual space of the trace space $W:=\left[ H_0^{\frac{1}{2}}(\Gamma_C)
-\right]^{\textrm{dim}}$ of $V$ restricted to $\Gamma_C$,
-@f{gather*}K:=\{\mu\in W':\mu_T = 0,\quad\langle\mu n,v\rangle_{\Gamma_C}\geq 0,\quad
-\forall v\in H_0^{\frac{1}{2}}(\Gamma_C), v \ge 0\text{ on }\Gamma_C \}@f}
-of Lagrange multipliers, where $\langle\cdot,\cdot\rangle$
-denotes the duality pairing, i.e. a boundary integral, between $W'$ and $W$.
-Intuitively, $K$ is the cone of all "non-positive functions", except that $ K\subset
-\left( \left[ H_0^{\frac{1}{2}}(\Gamma_C) \right]^{\textrm{dim}} \right)' $ and so contains other
-objects besides regular functions as well. This yields:
-
-Find $u\in V$ and $\lambda\in K$ such that
+where the projector $P_\Pi$ is defined as
@f{align*}
- \hat{a}(u,v) + b(v,\lambda) &= f(v),\quad &&v\in V\\
- b(u,\mu - \lambda) &\leq \langle g,(\mu -
- \lambda)n\rangle_{\Gamma_C},\quad&&\mu\in K,
+ P_{\Pi}(\tau):=\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}
@f}
-with
+and the space $V^+$ is the space of all displacements that satisfy the contact
+condition:
@f{align*}
- \hat{a}(u,v) &:= a'(u^i;u,v)\\
- b(u,\mu) &:= \langle un,\mu n\rangle_{\Gamma_C},\quad &&u\in V,\quad\mu\in W'.
+ 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\}.
@f}
-As in the section before $u^i$ denotes the linearization point for the
-semi-linearform $a(.;.)$. In contrast to step-41 we directly consider $\lambda$
-as the additional, positive force $\sigma(u)n$ that the obstacle
-exerts on the boundary $\Gamma_C$ of the body.
-The existence and uniqueness of the analytical solution $(u,\lambda)\in V\times
-K$ of this saddle point problem has been stated in Glowinski, Lions and Tr\'{e}moli\`{e}res: Numerical
-Analysis of Variational Inequalities, North-Holland, 1981.
-
-NOTE: In this example as well as in the further documentation we make the
-assumption that the normal vector $n$ equals to $(0,0,1)$. This comes up with
-the starting condition of our deformable body.
+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.
+
+
+<h3>A Newton method for the plastic nonlinearity</h3>
+
+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 <i>slant</i>
+differentiability and this turns out to be enough for Newton's method to work.
+The resulting method then goes by the name <i>semi-smooth Newton method</i>,
+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):
+@f{align*}
+ \label{eq:linearization}
+ \left(I_{\Pi}\varepsilon(\tilde {\mathbf u}^{i}),
+ \varepsilon(\varphi) - \varepsilon(\tilde {\mathbf u}^{i})\right) \geq 0,
+ \quad \forall \varphi\in V^+,
+@f}
+where the rank-4 tensor $I_\Pi=I_\Pi(\varepsilon^D(\mathbf u^{i-1}))$ given by
+@f{align}
+ I_\Pi = \begin{cases}
+ C_{\mu} + C_{\kappa}, & \hspace{-8em}\hfill \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}
+@f}
+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
+@f{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),
+@f}
+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
+@f{align*}
+ {\mathbf u}^{i} = {\mathbf u}^{i-1} + \alpha_i (\tilde {\mathbf u}^{i} - {\mathbf u}^{i-1})
+@f}
+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.
<h3>Active Set methods to solve the saddle point problem</h3>
-The linearized problem is essentially like a pure elastic problem with contact like
+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.
+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.
-<h3>The primal-dual active set algorithm combined with the inexact semi smooth
-Newton method</h3>
+<h3>Overall algorithm</h3>
-Now we describe an algorithm that combines the damped semismooth Newton-method,
+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 sums up the results of the sections before and works as follows:
+method for the contact. It works as follows:
<ol>
- <li> Initialize $\mathcal{A}_k$ and $\mathcal{F}_k$, such that
- $\mathcal{S} = \mathcal{A}_k \cup \mathcal{F}_k$ and $\mathcal{A}_k \cap
- \mathcal{F}_k = \emptyset$ and set $k = 1$. The start value $\hat U^0 :=
- P_{\mathcal{A}_k}(0)$ fulfills our obstacle condition.
- <li> Assemble the Newton matrix $A := a'(\hat
- U^{k-1};\varphi_p,\varphi_q)$ and the right-hand-side $F(\hat U^{k-1})$.
- <li> Find the primal-dual pair $(\bar U^k,\Lambda^k)$ that satisfies
+ <li> 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 :=
+ P_{\mathcal{A}_k}(0)$ fulfills our obstacle condition, i.e., we project an
+ initial zero displacement onto the set of feasible displacements.
+
+ <li> Assemble the Newton matrix $A_{pq} := 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.
+
+ <li> Find the primal-dual pair $(\tilde U^i,\Lambda^i)$ that satisfies
@f{align*}
- A\bar U^k + B\Lambda^k & = F, &\\
- \left[B^T\bar U^k\right]_p & = G_p & \forall p\in\mathcal{A}_k,\\
- \Lambda^k_p & = 0 & \forall p\in\mathcal{F}_k.
+ 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.
@f}
- <li> Damping for $k>2$ by applying a line search and calculating a linear
- combination of $U^{k-1}$ and $\bar U^k$. Find an
- $\alpha_i:=2^{-i},(i=0,\ldots,10)$ so that
- @f{gather*}U^k := \alpha_i\bar U^k +
- (1-\alpha_i)U^{k-1}@f}
- yields
- @f{gather*}\vert F\left(U^{k}\right) \vert < \vert F\left(U^{k-1}\right) \vert.\f}
+ 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})$.
+
+
+
+ <li> 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:=2^{-l},(l=0,\ldots,10)$ so that
+ @f{gather*}U^i := \alpha^i_l\bar U^i +
+ (1-\alpha^i_l)U^{i-1}@f}
+ satisfies
+ @f{gather*}
+ \vert F\left(U^{i}\right) \vert < \vert F\left(U^{i-1}\right) \vert.
+ \f}
+
<li> Define the new active and inactive sets by
- @f{gather*}\mathcal{A}_{k+1}:=\lbrace p\in\mathcal{S}:\Lambda^k_p +
- c\left(\left[B^TU^k\right]_p - G_p\right) > 0\rbrace,@f}
- @f{gather*}\mathcal{F}_{k+1}:=\lbrace p\in\mathcal{S}:\Lambda^k_p +
- c\left(\left[B^TU^k\right]_p - G_p\right) \leq 0\rbrace.@f}
- Projection $U^k$ so that it holds the second equation in (2)
- @f{gather*}\hat U^K := P_{\mathcal{A}_{k+1}}(U^k).@f}
- <li> If $\mathcal{A}_{k+1} = \mathcal{A}_k$ and $\vert
- F\left(U^{k}\right) \vert < \delta$ then stop, else set $k=k+1$ and go to
- step (1).
+ @f{gather*}\mathcal{A}_{i+1}:=\lbrace p\in\mathcal{S}:\Lambda^i_p +
+ c\left(\left[B^TU^i\right]_p - G_p\right) > 0\rbrace,@f}
+ @f{gather*}\mathcal{F}_{i+1}:=\lbrace p\in\mathcal{S}:\Lambda^i_p +
+ c\left(\left[B^TU^i\right]_p - G_p\right) \leq 0\rbrace.@f}
+
+ <li>Project $U^i$ so that it satisfies the contact inequality,
+ @f{gather*}\hat U^i := P_{\mathcal{A}_{i+1}}(U^i).@f}
+ Here,
+ $P_{\mathcal{A}}(U)$ is the projection of the active
+ components in $\mathcal{A}$ to the gap
+ @f{gather*}P_{\mathcal{A}}(U)_p:=\begin{cases}
+ U_p, & \textrm{if}\quad p\notin\mathcal{A}\\
+ g_{h,p}, & \textrm{if}\quad
+ p\in\mathcal{A},
+ \end{cases}@f}
+ where $g_{h,p}$ is the <i>gap</i> denoting the distance of the obstacle
+ from the undisplaced configuration of the body.
+
+ <li> If $\mathcal{A}_{i+1} = \mathcal{A}_k$ and $\vert
+ F\left(U^{i}\right) \vert < \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.
</ol>
-The subscript $p$ denotes a vertex and the meaning of the
-decorated and none decorated $U$ is as follows: $\bar U$ denotes the solution of the linear system of equations in (2), $U$ is the
-damped solution and equals to $\bar U$ if the damping parameter $\alpha_0 =
-1$ and $\hat U := P_{\mathcal{A}}(U)$ is the projection of the active
-components in $\mathcal{A}$ to the gap
-
-@f{gather*}P_{\mathcal{A}}(U):=\begin{cases}
-U_p, & \textrm{if}\quad p\notin\mathcal{A}\\
-g_{h,p}, & \textrm{if}\quad
-p\in\mathcal{A}.
-\end{cases}@f}\\
-The matrix $B\in\mathbb{R}^{n\times m}$, $n>m$ describes the coupling of the
+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$. Due to the ansatz functions $\psi_i$ (scalar valued) of the
-Lagrange multiplier are fulfilling the following biorthogonal condition (see Hüeber, Wohlmuth: A primal–dual active
+$\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)
-@f{gather} \int\limits_{\Gamma_C}\psi_i(x)\varphi_j(x)dx =
-\delta_{ij}\int\limits_{\Gamma_C}\varphi_j(x)dx@f}
-this yields
-@f{gather*}B_{pq} = \begin{cases}
-\int\limits_{\Gamma_C}\varphi_p(x)dxI_3, & \text{if}\quad p=q,\quad
-p,q\in\mathcal{S}\\
-0I_3, & \text{if}\quad p\neq q,\quad p\textrm{ or }q\notin\mathcal{S}.
-\end{cases}@f}
-Here $I_3$ denotes the three-dimensional identity matrix.
-In our program we use the structure of a quadratic sparse for
-$B\in\mathbb{R}^{n\times n}$ and for $\Lambda^k$ a vector with length $n$ where
-$\Lambda^k_p = 0$ for $p\notin \mathcal{S}$.
+194, 2005, pp. 3147-3166).
The vector $G$ is defined by a suitable approximation $g_h$ of the gap $g$
@f{gather*}G_p = \begin{cases}
-g_{h,p}\int\limits_{\Gamma_C}\varphi_p(x)dx, & \text{if}\quad p\in\mathcal{S}\\
+g_{h,p}, & \text{if}\quad p\in\mathcal{S}\\
0, & \text{if}\quad p\notin\mathcal{S}.
\end{cases}@f}
-Note that $G_p$ is a three-dimensional vector and that again we applied the
-biorthogonal property of the Lagrange multiplier ansatz functions to the
-integral $\int\limits_{\Gamma_C}g_h(x)\varphi_p(x)dx$ with $g_h(x)=\sum\limits_i
-g_{h,p}\varphi_p(x)$ (see the reference mentioned above).
-
-Compared to step-41, step (1) is added but it should be clear
-from the sections above that we only linearize the problem. In step (2) we have
-to solve a linear system of equations again. And now the solution has to fulfill two stopping
-criteria. $\mathcal{A}_{k+1} = \mathcal{A}_k$ makes sure that the contact zones
-are iterated out and the second ensures an accurate enough residual which means
-that the plastic zones are also iterated out.
-
-A similar method can also be found in Brunssen, Schmid, Schäfer,
-Wohlmuth: A fast and robust iterative solver for nonlinear contact problems
-using a primal-dual active set strategy and algebraic multigrid, Int. J. Numer.
-Meth. Engng, 2007, 69, pp. 524-543. But in advance we apply a line search to
-obtain a more robust method regarding the start value. Solving an elastic
-problem in the very first step ($k=1$) we get a reasonable start value but as
-you can see in the results damping is important if we choose for example a ball as
-obstacle.
-
-Damping our Newton method is more important for the nonlinearity cause by the
-constitutive law as for the contact. For this reason we start to damp our method
-for $k>2$ when we have two plastic iterations. Note that $U^1$ is a solution of
-an elastic problem and $U^2$ is the first plastic solution. A linear combination
-between these both results in stresses which are not in the convex set
-of the feasible stresses.
-
<h3>Adaptive mesh refinement</h3>
-Since we run our program in 3d, there is a good reason to use adaptive
-mesh refinement. To make things a bit easier we are choosing the
+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!
+as for plasticity.
<h3>Implementation</h3>
(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.
-Another difficulty is the handling of the different constraints from
-(the Dirichlet conditions), the hanging nodes and the inequality condition that
+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
-ConstraintMatrix.
-
-Beside the ConstitutiveLaw class there is another new class called Input. This
-class allows us to read in an obstacle from a file. In our example the file
-'obstacle_file.dat' contains data which describe an Chinese, Japanese or
-Korean symbol for force or power. (See http://www.orientaloutpost.com/ :
+ConstraintMatrix that describe the various constraints and that we will
+combine as appropriate in each iteration.
+
+Compared to step-41, the programs has two new classes:
+<code>ConstitutiveLaw</code> describes the plastic behavior of the material,
+and <code>Input</code> is responsible for the parameters that describe
+the overall setup. In particular, the latter
+class 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)".)
+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 - see the picture at
+the top of this section.
projects / dealii-svn.git / commitdiff