+++ /dev/null
-\documentclass{article}
-
-\usepackage{amsmath}
-\usepackage{amssymb}
-
-\title{Documentation of step-42, An obstacle problem for elasto-plastic material behavior in three dimensions}
-\author{Joerg Frohne}
-\date{Juni, 2012}
-
-\begin{document}
-
-\section{Introduction}
-
-This example is an extension of step-41, considering a 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
-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
-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.\\
-
-\begin{huge}
-{distributed}
-\end{huge}
-
-
-\section{Classical formulation}
-
-The classical formulation of the problem possesses the following form:
-\begin{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
-\end{align*}
-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
-$$\mathcal{F}(\tau) = \vert\tau^D\vert - \sigma_0¸\quad\text{with}\quad \tau^D
-= \tau - \dfrac{1}{3}tr(\tau)I,$$
-$\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
-$$ \sigma_n = \sigma n\cdot n$$
-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. We refer that you have to ensure that the obstacle does not hit
-the boundary of $\Gamma_C$
-
-\section{Derivation of the variational inequality}
-
-As a starting point to derive the equations above, let us imagine that we want
-to minimise an energy functional:
-$$E(\tau) := \dfrac{1}{2}\int\limits_{\Omega}\tau A \tau d\tau,\quad \tau\in \Pi W^{\textrm{div}}$$
-with
-$$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$$ and
-$$\Pi \Sigma:=\lbrace \tau\in \Sigma, \mathcal{F}(\tau)\leq 0\rbrace$$
-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:
-$$L(\tau,\varphi) := E(\tau) + (\varphi, \textrm{div}(\tau)),\quad \lbrace\tau,\varphi\rbrace\in\Pi W^{\textrm{div}}\times V^+$$
-with
-$$V^+ := \lbrace u\in V: u_n\leq g \text{ on } \Gamma_C \rbrace$$
-$$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$$
-By building the Fr\'echet derivatives of $L$ for both components we obtain the
-dual formulation for the stationary case which is known as \textbf{Hencky-Type-Model}:\\
-Find a pair $\lbrace\sigma,u\rbrace\in \Pi W\times V^+$ with
-$$\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}}$$
-$$-\left(\textrm{div}(\sigma),\varphi - u\right) \geq 0,\quad \forall \varphi\in V^+.$$
-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
-$$\left(\sigma,\tau - \sigma\right) - \left(C\varepsilon(u), \tau - \sigma\right) \geq 0,\quad \forall \tau\in \Pi W$$
-$$\left(\sigma,\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in V^+.$$
-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.
-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
-$$\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}))$$
-$$\left(\sigma,\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in V^+,$$
-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\\
-$$\left(\sigma - C\varepsilon(u), \tau - \sigma\right) \geq 0,\quad \forall \tau\in \Pi W,$$
-which yields with the second inequality:\\
-Find the displacement $u\in V^+$ with
-$$\left(P_{\Pi}(C\varepsilon(u)),\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in V^+,$$
-with the projection:
-$$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}$$
-with the radius
-$$\hat\alpha := \sigma_0 + \gamma^{\text{iso}}\xi .$$
-With the relation $\xi = \vert\varepsilon(u) - A\sigma\vert$ it is possible to eliminate $\xi$ inside the projection $P_{\Pi}$:\\
-$$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}$$
-$$\alpha := \sigma_0 + \dfrac{\gamma^{\text{iso}}}{2\mu+\gamma^{\text{iso}}}\left(\vert\tau^D\vert - \sigma_0\right) ,$$
-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
-$$\gamma := \dfrac{\gamma^{\text{iso}}}{2\mu +
-\gamma^{\text{iso}}}\in[0,1)\text{ with }\gamma^{\text{iso}}\in[0,\infty),$$
-$$\beta :=\dfrac{\sigma_0}{\vert\tau^D\vert}.$$ 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 $$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}.$$
-For further details 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 $$\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)$$
-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.
-
-\section{Linearization of the constitutive law for the Newton method}
-
-To apply Newton method we have to linearise the following semi-linearform
-$$a(\psi;\varphi) :=
-\left(P_{\Pi}(C\varepsilon(\psi)),\varepsilon(\varphi)\right).$$ 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
-\begin{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^+.
-\end{eqnarray*}
-Neglecting the constant terms $ a(u^i;u^i)$ and $ a'(u^i;u^i,u^i)$ we obtain the
-following minimisation problem $$\dfrac{1}{2} a'(u^i;u^{i+1},u^{i+1}) -
-F(u^{i+1})\rightarrow \textrm{min},\quad u^{i+1}\in V^+$$ with $$F(\varphi) :=
-\left(a'(u^{i};u^{i},\varphi) - a(u^{i},\varphi) \right).$$
-In the case of our constitutive law the Fr\'echet derivative of the
-semi-linearform $a(.;.)$ at the point $u^i$ is
-
-$$a'(u^i;\psi,\varphi) =
-(I(x)\varepsilon(\psi),\varepsilon(\varphi)),\quad x\in\Omega,$$ $$
-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}
-$$
-with
-$$C_{\mu} := 2\mu\left(\mathbb{I} - \dfrac{1}{3} I\otimes
-I\right)\quad\text{(shear part of the stress strain tensor)},$$
-$$C_{\kappa} := \kappa I\otimes I\quad\text{(bulk part of the stress strain
-tensor)},$$
-$$\tau^D := C\varepsilon^D(u^i).$$
-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.
-
-\section{Formulation as a saddle point problem}
-
-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$, $$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 \}$$
-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:
-
-\textit{Find $u\in V$ and $\lambda\in K$ such that}
-\begin{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,
-\end{align*}
-\textit{with}
-\begin{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'.
-\end{align*}
-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.\\
-
-\noindent
-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.\\
-
-\noindent
-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.
-
-\section{Active Set methods to solve the saddle point problem}
-
-The linearized problem is essentially like a pure elastic problem with contact 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.
-
-\section{The primal-dual active set algorithm combined with the inexact semi smooth
-Newton method}
-
-Now we describe an algorithm that 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:
-\begin{itemize}
- \item[(0)] 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.
- \item[(1)] Assemble the Newton matrix $A := a'(\hat
- U^{k-1};\varphi_p,\varphi_q)$ and the right-hand-side $F(\hat U^{k-1})$.
- \item[(2)] Find the primal-dual pair $(\bar U^k,\Lambda^k)$ that satisfies
- \begin{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.
- \end{align*}
- \item[(3)] 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 $$U^k := \alpha_i\bar U^k +
- (1-\alpha_i)U^{k-1}$$ yields $$\vert
- F\left(U^{k}\right) \vert < \vert F\left(U^{k-1}\right) \vert.$$
- \item[(4)] Define the new active and inactive sets by
- $$\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,$$
- $$\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.$$
- Projection $U^k$ so that it holds the second equation in (2)
- $$\hat U^K := P_{\mathcal{A}_{k+1}}(U^k).$$
- \item[(5)] 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).
-\end{itemize}
-\noindent
-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
-$$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}$$\\
-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
-set strategy for non-linear multibody contact problems, Comput. Methods Appl. Mech. Engrg.
-194, 2005, pp. 3147-3166)
-$$ \int\limits_{\Gamma_C}\psi_i(x)\varphi_j(x)dx =
-\delta_{ij}\int\limits_{\Gamma_C}\varphi_j(x)dx$$
-this yields
-$$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}$$
-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}$.
-The vector $G$ is defined by a suitable approximation $g_h$ of the gap $g$
-$$G_p = \begin{cases}
-g_{h,p}\int\limits_{\Gamma_C}\varphi_p(x)dx, & \text{if}\quad p\in\mathcal{S}\\
-0, & \text{if}\quad p\notin\mathcal{S}.
-\end{cases}$$\\
-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, Schaefer, 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.
-
-\section{Adaptive mesh refinement}
-
-Since we run our program in 3d, there is a good reason to use adaptive
-mesh refined. To make life a bit easier we are choosing 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.
-
-\section{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.
-
-Another difficulty is the handling of the different constraints from
-(the Dirichlet conditions), the 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 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)'').
-
-\section{Results}
-
-
-
-\end{document}
projects / dealii-svn.git / commitdiff