From: bangerth Date: Sat, 25 Feb 2012 10:33:25 +0000 (+0000) Subject: Convert intro from latex to doxygen format. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=15f4efef78d810136a25d5ce8f36a308347e45cf;p=dealii-svn.git Convert intro from latex to doxygen format. git-svn-id: https://svn.dealii.org/trunk@25170 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-41/doc/intro.dox b/deal.II/examples/step-41/doc/intro.dox index 3438be49a4..2b5601442f 100644 --- a/deal.II/examples/step-41/doc/intro.dox +++ b/deal.II/examples/step-41/doc/intro.dox @@ -6,3 +6,413 @@ This material is based upon work partly supported by the ... +

Introduction

+ +This example is based on the Laplace equation in 2d and deals with the +question what happens if a membrane is deflected by some external force but is +also constrained by an obstacle. In other words, think of a elastic membrane +clamped at the boundary to a rectangular frame (we choose $\Omega = +\left[-1,1\right]^2$) and that sags through due to gravity acting on it. What +happens now if there is an obstacle under the membrane that prevents it from +reaching its equilibrium position if gravity was the only existing force? In +the current example program, we will consider that under the membrane is a +stair step obstacle against which gravity pushes the membrane. + +This problem is typically called the "obstacle problem", and it results in a +variational inequality, rather than a variational equation when put into the +weak form. We will below derive it from the classical formulation, but before we +go on to discuss the mathematics let us show how the solution of the problem we +will consider in this tutorial program looks to gain some intuition of what +we should expect: + + + + + + +
+ @image html "step-41.displacement.png" + + @image html "step-41.active-set.png" +
+ +Here, at the left, we see the displacement of the membrane. The shape +of the obstacle underneath is clearly visible. On the right, we overlay which +parts of the membrane are in contact with the obstacle. We will later call +this set of points the "active set" to indicate that an inequality constraint +is active there. + + +

Classical formulation

+ +The classical formulation of the problem possesses the following form: +@f{align*} + -\textrm{div}\ \sigma &\geq f & &\quad\text{in } \Omega,\\ + \sigma &= \nabla u & &\quad\text{in } \Omega,\\ + u(\mathbf x) &= 0 & &\quad\text{on }\partial\Omega,\\ +(-\Delta u - f)(u - g) &= 0 & &\quad\text{in } \Omega,\\ + u(\mathbf x) &\geq g(\mathbf x) & &\quad\text{in } \Omega +@f} +with $u\in H^2(\Omega)$. $u$ is a scalar valued function that denotes the +vertical displacement of the membrane. The first equation is called equilibrium +condition with a force of areal density $f$. Here, we will consider this force +to be gravity. The second one is known as Hooke's Law that says that the stresses +$\sigma$ are proportional to the gradient of the displacements $u$ (the +proportionality constant, often denoted by $E$, has been set to one here, +without loss of generality; if it is constant, it can be put into the right +hand side function). At the boundary we have zero Dirichlet +conditions. Obviously, the first two equations can be combined to yield +$-\Delta u \ge f$. + +Intuitively, gravity acts downward and so $f(\mathbf x)$ is a negative +function (we choose $f=-10$ in this program). The first condition then means +that the total force acting on the membrane is gravity plus something +positive: namely the upward force that the obstacle exerts on the membrane at +those places where the two of them are in contact. How big is this additional +force? We don't know yet (and neither do we know "where" it actually acts) but +it must be so that the membrane doesn't penetrate the obstacle. + +The fourth equality above together with the last inequality forms the obstacle +condition which has to hold at every point of the whole domain. The latter of +these two means that the membrane must be above the obstacle $g(\mathbf x)$ +everywhere. The second to last equation, often called the "complementarity +condition" says that where the membrane is not in contact with the obstacle +(i.e., those $\mathbf x$ where $u(\mathbf x) - g(\mathbf x) \neq 0$), then +$-\Delta u=f$ at these locations; in other words, no additional forces act +there, as expected. On the other hand, where $u=g$ we can have $-\Delta u-f +\neq 0$, i.e., there can be additional forces (though there don't have to be: +it is possible for the membrane to just touch, not press against, the +obstacle). + + +

Derivation of the variational inequality

+ +An obvious way to obtain the variational formulation of the obstacle problem is to consider the total potential energy: +@f{equation*} + E(u):=\dfrac{1}{2}\int\limits_{\Omega} \nabla u \cdot \nabla - \int\limits_{\Omega} fu. +@f} +We have to find a solution $u\in G$ of the following minimization problem: +@f{equation*} + E(u)\leq E(v)\quad \forall v\in G, +@f} +with the convex set of admissble displacements: +@f{equation*} + G:=\lbrace v\in V: v\geq g \text{ a.e. in } \Omega\rbrace,\quad V:=H^1_0(\Omega). +@f} +This set takes care of the third and fifth conditions above (the boundary +values and the complementarity condition). + +Consider now the minimizer $u\in G$ of $E$ and any other function $v\in +G$. Then the function +@f{equation*} + F(\varepsilon) := E(u+\varepsilon(v-u)),\quad\varepsilon\in\left[0,1\right], +@f} +takes its minimum at $\varepsilon = 0$ (because $u$ is a minimizer of the +energy functional $E(\cdot)$), so that $F'(0)\geq 0$ for any choice +of $v$. Note that +$u+\varepsilon(v-u) = (1-\varepsilon)u+\varepsilon v\in G$ because of the +convexity of $G$. If we compute $F'(\varepsilon)\vert_{\varepsilon=0}$ it +yields the variational formulation we are searching for: + +Find a function $u\in G$ with +@f{equation*} + \left(\nabla u, \nabla(v-u)\right) \geq \left(f,v-u\right) \quad \forall v\in G. +@f} + +This is the typical form of variational inequalities, where not just $v$ +appears in the bilinear form but in fact $v-u$. The reason is this: if $u$ is +not constrained, then we can find test functions $v$ in $G$ so that $v-u$ can have +any sign. By choosing test functions $v_1,v_2$ so that $v_1-u = -(v_2-u)$ it +follows that the inequality can only hold for both $v_1$ and $v_2$ if the two +sides are in fact equal, i.e., we obtain a variational equality. + +On the other hand, if $u=g$ then $G$ only allows test functions $v$ so that in fact +$v-u\ge 0$. This means that we can't test the equation with both $v-u$ and +$-(v-u)$ as above, and so we can no longer conclude that the two sides are in +fact equal. Thus, this mimicks the way we have discussed the complementarity +condition above. + + + +

Formulation as a saddle point problem

+ +The variational inequality above is awkward to work with. We would therefore +like to reformulate it as an equivalent saddle point problem. Set $V:=H^1_0(\Omega)$. +We introduce a Lagrange multiplier $\lambda$ and the convex cone $K\subset V'$, $V'$ +dual space of $V$, $K:=\{\mu\in V': \langle\mu,v\rangle\geq 0,\quad \forall +v\in V, v \le 0 \}$ of +Lagrange multipliers, where $\langle\cdot,\cdot\rangle$ denotes the duality +pairing between $V'$ and $V$. Intuitively, $K$ is the cone of all "non-positive +functions", except that $K\subset (H_0^1)'$ and so contains other objects +besides regular functions as well. +This yields: + +Find $u\in V$ and $\lambda\in K$ such that +@f{align*} + a(u,v) + b(v,\lambda) &= f(v),\quad &&v\in V\\ + b(u,\mu - \lambda) &\leq \langle g,\mu - \lambda\rangle,\quad&&\mu\in K, +@f} +with +@f{align*} + a(u,v) &:= \left(\nabla u, \nabla v\right),\quad &&u,v\in V\\ + b(u,\mu) &:= \langle g-u,\mu\rangle,\quad &&u\in V,\quad\mu\in V'. +@f} +In other words, we can consider $\lambda$ as the negative of the additional, positive force that the +obstacle exerts on the membrane. The inequality in the second line of the +statement above only appears to have the wrong sign because we have +$\mu-\lambda<0$ at points where $\lambda=0$, given the definition of $K$. + +The existence and uniqueness of $(u,\lambda)\in V\times K$ of this saddle +point problem has been stated in Grossmann and Roos: Numerical treatment of +partial differential equations, Springer-Verlag, Heidelberg-Berlin, 2007, 596 +pages, ISBN 978-3-540-71582-5. + + + +

Active Set methods to solve the saddle point problem

+ +There are different methods to solve the variational inequality. As one +possibility you can understand the saddle point problem as a convex quadratic program (QP) with +inequality constraints. + +To get there, let us assume that we discretize both $u$ and $\lambda$ with the +same finite element space, for example the usual $Q_k$ spaces. We would then +get the equations +@f{eqnarray*} + &A U + B\Lambda = F,&\\ + &[BU-G]_i \geq 0, \quad \Lambda_i \leq 0,\quad \Lambda_i[BU-G]_i = 0 +\qquad \forall i.& +@f} +where $B$ is the mass matrix on the chosen finite element space and the +indices $i$ above are for all degrees of freedom in the set $\cal S$ of degrees of +freedom located in the interior of the domain +(we have Dirichlet conditions on the perimeter). However, we +can make our life simpler if we use a particular quadrature rule when +assembling all terms that yield this mass matrix, namely a quadrature formula +where quadrature points are only located at the interpolation points at +which shape functions are defined; since all but one shape function are zero +at these locations, we get a diagonal mass matrix with +@f{align*} + B_{ii} = \int_\Omega \varphi_i(\mathbf x)^2\ \textrm{d}x, + \qquad + B_{ij}=0 \ \text{for } i\neq j. +@f} +To define $G$ we use the same technique as for $B$. In other words, we +define +@f{align*} + G_{i} = \int_\Omega g_h(x) \varphi_i(\mathbf x)\ \textrm{d}x. +@f} +where $g_h$ is a +suitable approximation of $g$ and $\mathbf x_i$ is the support point of the +$i$the shape function. The integral in the definition of $B_{ii}$ and $G_i$ +are then approximated by the trapezoidal rule. +With this, the equations above can be restated as +@f{eqnarray*} + &A U + B\Lambda = F,&\\ + &U_i-B_{ii}^{-1}G_i \ge 0, \quad \Lambda_i \leq 0,\quad \Lambda_i[U_i-B_{ii}^{-1}G_i] = 0 +\qquad \forall i\in{\cal S}.& +@f} + +Now we define for each degree of freedom $i$ the function +@f{equation*} + C([BU]_i,\Lambda_i):=-\Lambda_i + \min\lbrace 0, \Lambda_i + c([BU]_i - G_i) \rbrace, +@f} +with some $c>0$. (In this program we choose $c = 100$. It is a kind of a +penalty parameter which depends on the problem itself and needs to be chosen +large enough; for example there is no convergence for $c = 1$ using the +current program if we use 7 global refinements.) + +After some headscratching one can then convince oneself that the inequalities +above can equivalently be rewritten as +@f{equation*} + C([BU]_i,\Lambda_i) = 0, \qquad \forall i\in{\cal S}. +@f} +The primal-dual active set strategy we will use here is an iterative scheme which is based on +this condition to predict the next active and inactive sets $\mathcal{A}_k$ and +$\mathcal{F}_k$ (that is, those complementary sets of indices $i$ for which +$U_i$ is either equal to or not equal to the value of the obstacle +$B^{-1}G$). For a more in depth treatment of this approach, see 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. + +

The primal-dual active set algorithm

+ +The algorithm for the primal-dual active set method works as follows: + +- [(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$. +- [(1)] Find the primal-dual pair $(U^k,\Lambda^k)$ that satisfies + @f{align*} + AU^k + B\Lambda^k &= F,\\ + [BU^k]_i &= G\quad&&\forall i\in\mathcal{A}_k,\\ + \Lambda_i^k &= 0\quad&&\forall i\in\mathcal{F}_k. + @f} + Note that the second and third conditions imply that exactly $|S|$ unknowns + are fixed, with the first condition yielding the remaining $|S|$ equations + necessary to determine both $U$ and $\Lambda$. +- [(2)] Define the new active and inactive sets by + @f{equation*} + \begin{split} + \mathcal{A}_{k+1}:=\lbrace i\in\mathcal{S}:\Lambda^k_i + c([BU^k]_i - G_i)< 0\rbrace,\\ + \mathcal{F}_{k+1}:=\lbrace i\in\mathcal{S}:\Lambda^k_i + c([BU^k]_i - G_i)\geq 0\rbrace. + \end{split} + @f} +- [(3)] If $\mathcal{A}_{k+1}=\mathcal{A}_k$ (and then, obviously, also + $\mathcal{F}_{k+1}=\mathcal{F}_k$) then stop, else set $k=k+1$ and go to step + (1). + +The method is called "primal-dual" because it uses both primal (the +displacement $U$) as well as dual variables (the Lagrange multiplier +$\Lambda$) to determine the next active set. + +At the end of this section, let us add two observations. First, +for any primal-dual pair $(U^k,\Lambda^k)$ that satisfies these +condition, we can distinguish the following cases: + +- [1.] $\Lambda^k_i + c([BU^k]_i - G_i) < 0$ (i active): +
+ Then either $[BU^k]_i + Then either $[BU^k]_i\geq G_i$ and $\Lambda^k_i=0$ (no contact) or $\Lambda^k_i\geq0$ and $[BU^k]_i=G_i$ (unpressing load). + +Second, the method above appears untuitively correct and useful but a bit ad +hoc. However, it can be derived in a concisely in the following way. To this +end, note that we'd like to solve the nonlinear system +@f{eqnarray*} + &A U + B\Lambda = F,&\\ + &C([BU-G]_i, \Lambda_i) = 0, +\qquad \forall i.& +@f} +We can iteratively solve this by always linearizing around the previous +iterate (i.e., applying a Newton method), but for this we need to linearize +the function $C(\cdot,\cdot)$ that is not differentiable. That said, it is +slantly differentiable, and in fact we have +@f{equation*} + \dfrac{\partial}{\partial U^k_i}C([BU^k]_i,\Lambda^k_i) = \begin{cases} + cB_{ii},& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)< 0\\ + 0,& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)\geq 0. + \end{cases} +@f} +@f{equation*} + \dfrac{\partial}{\partial\Lambda^k_i}C([BU^k]_i,\Lambda^k_i) = \begin{cases} + 0,& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)< 0\\ + -1,& \text{if}\ \Lambda^k_i + c([BU^k]_i - G_i)\geq 0. + \end{cases} +@f} +This suggest a semismooth Newton step of the form +@f{equation*} + \begin{pmatrix} + A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} & B_{\mathcal{F}_k} & 0\\ + A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} & 0 & B_{\mathcal{A}_k}\\ + 0 & 0 & -Id_{\mathcal{F}_k} & 0\\ + 0 & cB_{\mathcal{A}_k} & 0 & 0 +\end{pmatrix} +\begin{pmatrix} + \delta U^k_{\mathcal{F}_k}\\ \delta U^k_{\mathcal{A}_k}\\ \delta \Lambda^k_{\mathcal{F}_k}\\ \delta \Lambda^k_{\mathcal{A}_k} +\end{pmatrix} += +-\begin{pmatrix} + (AU^k + \Lambda^k - F)_{\mathcal{F}_k}\\ (AU^k + \Lambda^k - F)_{\mathcal{A}_k}\\ -\Lambda^k_{\mathcal{F}_k}\\ c(B_{\mathcal{A}_k} U^k - G)_{\mathcal{A}_k} +\end{pmatrix}, +@f} +where we have split matrices $A,B$ as well as vectors in the natural way into +rows and columns whose indices belong to either the active set +${\mathcal{A}_k}$ or the inactive set ${\mathcal{F}_k}$. + +Rather than solving for updates $\delta U, \delta \Lambda$, we can also solve +for the variables we are interested in right away by setting $\delta U^k := +U^{k+1} - U^k$ and $\delta \Lambda^k := \Lambda^{k+1} - \Lambda^k$ and +bringing all known terms to the right hand side. This yields +@f{equation*} +\begin{pmatrix} + A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} & B_{\mathcal{F}_k} & 0\\ + A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} & 0 & B_{\mathcal{A}_k}\\ + 0 & 0 & Id_{\mathcal{F}_k} & 0\\ + 0 & B_{\mathcal{A}_k} & 0 & 0 +\end{pmatrix} +\begin{pmatrix} + U^k_{\mathcal{F}_k}\\ U^k_{\mathcal{A}_k}\\ \Lambda^k_{\mathcal{F}_k}\\ \Lambda^k_{\mathcal{A}_k} +\end{pmatrix} += +\begin{pmatrix} + F_{\mathcal{F}_k}\\ F_{\mathcal{A}_k}\\ 0\\ G_{\mathcal{A}_k} +\end{pmatrix}. +@f} +These are the equations outlined above in the description of the basic algorithm. + +We could even drive this a bit further. +It's easy to see that we can eliminate the third row and the third column +because it implies $\Lambda_{\mathcal{F}_k} = 0$: +@f{equation*} +\begin{pmatrix} + A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} & 0\\ + A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} & B_{\mathcal{A}_k}\\ + 0 & B_{\mathcal{A}_k} & 0 +\end{pmatrix} +\begin{pmatrix} + U^k_{\mathcal{F}_k}\\ U^k_{\mathcal{A}_k}\\ \Lambda^k_{\mathcal{A}_k} +\end{pmatrix} += +\begin{pmatrix} + F_{\mathcal{F}_k}\\ F_{\mathcal{A}_k}\\ G_{\mathcal{A}_k} +\end{pmatrix}. +@f} +This shows that one in fact only needs to solve for the Lagrange multipliers +located on the active set. By considering the second row one would then recover +the full Lagrange multiplier vector through +@f{equation*} + \Lambda^k_S = B^{-1}\left(f_{\mathcal{S}} - A_{\mathcal{S}}U^k_{\mathcal{S}}\right). +@f} +Because of the third row and the fact that $B_{\mathcal{A}_k}$ is a diagonal matrix we are able +to calculate $U^k_{\mathcal{A}_k}=B^{-1}_{\mathcal{A}_k}G_{\mathcal{A}_k}$ directly. We can therefore also write the +linear system as follows: +@f{equation*} +\begin{pmatrix} + A_{\mathcal{F}_k\mathcal{F}_k} & 0\\ + 0 & Id_{\mathcal{A}_k} \\ +\end{pmatrix} +\begin{pmatrix} + U^k_{\mathcal{F}_k}\\ U^k_{\mathcal{A}_k} +\end{pmatrix} += +\begin{pmatrix} + F_{\mathcal{F}_k} - A_{\mathcal{F}_k\mathcal{A}_k}B^{-1}_{\mathcal{A}_k}G_{\mathcal{A}_k} + \\ + B_{\mathcal{A}_k}^{-1}G_{\mathcal{A}_k} +\end{pmatrix}. +@f} +Fortunately, this form is easy to arrive at: we simply build the usual Laplace +linear system +@f{equation*} +\begin{pmatrix} + A_{\mathcal{F}_k\mathcal{F}_k} & A_{\mathcal{F}_k\mathcal{A}_k} \\ + A_{\mathcal{A}_k\mathcal{F}_k} & A_{\mathcal{A}_k\mathcal{A}_k} +\end{pmatrix} +\begin{pmatrix} + U^k_{\mathcal{F}_k}\\ U^k_{\mathcal{A}_k} +\end{pmatrix} += +\begin{pmatrix} + F_{\mathcal{F}_k}\\ F_{\mathcal{A}_k} +\end{pmatrix}, +@f} +and then let the ConstraintMatrix class eliminate all constrained degrees of +freedom, namely $U^k_{\mathcal{A}_k}=B^{-1}_{\mathcal{A}_k}G_{\mathcal{A}_k}$, +in the same way as if the dofs in $\mathcal{A}_k$ were Dirichlet data. The +result linear system (the second to last one above) is symmetric and positive +definite and we solve it with a CG-method +and the AMG preconditioner from Trilinos. + + +

Implementation

+ +This tutorial is quite similar to step-4. The general structure of the program +follows step-4 with minor differences: +- We need two new methods, assemble_mass_matrix_diagonal and + update_solution_and_constraints. +- We need new member variables that denote the constraints we have here. +- We change the preconditioner for the solver. + + diff --git a/deal.II/examples/step-41/doc/step-41-doc.tex b/deal.II/examples/step-41/doc/step-41-doc.tex deleted file mode 100644 index e43cdaa282..0000000000 --- a/deal.II/examples/step-41/doc/step-41-doc.tex +++ /dev/null @@ -1,411 +0,0 @@ -\documentclass{article} - -\usepackage{amsmath} -\usepackage{amssymb} -\usepackage{a4wide} -\usepackage{graphicx} - -\title{Documentation of step-41, The obstacle problem} -\author{Joerg Frohne} -\date{November 11th, 2011} -\begin{document} -\maketitle - -\section{Introduction} - -This example is based on the Laplace equation in 2d and deals with the -question what happens if a membrane is deflected by some external force but is -also constrained by an obstacle. In other words, think of a elastic membrane -clamped at the boundary to a rectangular frame (we choose $\Omega = -\left[-1,1\right]^2$) and that sags through due to gravity acting on it. What -happens now if there is an obstacle under the membrane that prevents it from -reaching its equilibrium position if gravity was the only existing force? In -the current example program, we will consider that under the membrane is a -stair step obstacle against which gravity pushes the membrane. - -This problem is typically called the "obstacle problem", and it results in a -variational inequality, rather than a variational equation when put into the -weak form. We will below derive it from the classical formulation, but before we -go on to discuss the mathematics let us show how the solution of the problem we -will consider in this tutorial program looks to gain some intuition of what -we should expect: - -XXX (see files step-41.*.png) XXX - -Here, at the left, we see the displacement of the membrane. The shape -of the obstacle underneath is clearly visible. On the right, we overlay which -parts of the membrane are in contact with the obstacle. We will later call -this set of points the "active set" to indicate that an inequality constraint -is active there. - - -\section{Classical formulation} - -The classical formulation of the problem possesses the following form: -\begin{align*} - -\textrm{div}\ \sigma &\geq f & &\quad\text{in } \Omega,\\ - \sigma &= \nabla u & &\quad\text{in } \Omega,\\ - u(\mathbf x) &= 0 & &\quad\text{on }\partial\Omega,\\ -(-\Delta u - f)(u - g) &= 0 & &\quad\text{in } \Omega,\\ - u(\mathbf x) &\geq g(\mathbf x) & &\quad\text{in } \Omega -\end{align*} -with $u\in H^2(\Omega)$. $u$ is a scalar valued function that denotes the -vertical displacement of the membrane. The first equation is called equilibrium -condition with a force of areal density $f$. Here, we will consider this force -to be gravity. The second one is known as Hooke's Law that says that the stresses -$\sigma$ are proportional to the gradient of the displacements $u$ (the -proportionality constant, often denoted by $E$, has been set to one here, -without loss of generality; if it is constant, it can be put into the right -hand side function). At the boundary we have zero Dirichlet -conditions. Obviously, the first two equations can be combined to yield -$-\Delta u \ge f$. - -Intuitively, gravity acts downward and so $f(\mathbf x)$ is a negative -function (we choose $f=-10$ in this program). The first condition then means -that the total force acting on the membrane is gravity plus something -positive: namely the upward force that the obstacle exerts on the membrane at -those places where the two of them are in contact. How big is this additional -force? We don't know yet (and neither do we know "where" it actually acts) but -it must be so that the membrane doesn't penetrate the obstacle. - -The fourth equality above together with the last inequality forms the obstacle -condition which has to hold at every point of the whole domain. The latter of -these two means that the membrane must be above the obstacle $g(\mathbf x)$ -everywhere. The second to last equation, often called the "complementarity -condition" says that where the membrane is not in contact with the obstacle -(i.e., those $\mathbf x$ where $u(\mathbf x) - g(\mathbf x) \neq 0$), then -$-\Delta u=f$ at these locations; in other words, no additional forces act -there, as expected. On the other hand, where $u=g$ we can have $-\Delta u-f -\neq 0$, i.e., there can be additional forces (though there don't have to be: -it is possible for the membrane to just touch, not press against, the -obstacle). - - -\section{Derivation of the variational inequality} - -An obvious way to obtain the variational formulation of the obstacle problem is to consider the total potential energy: -\begin{equation*} - E(u):=\dfrac{1}{2}\int\limits_{\Omega} \nabla u \cdot \nabla - \int\limits_{\Omega} fu. -\end{equation*} -We have to find a solution $u\in G$ of the following minimization problem: -\begin{equation*} - E(u)\leq E(v)\quad \forall v\in G, -\end{equation*} -with the convex set of admissble displacements: -\begin{equation*} - G:=\lbrace v\in V: v\geq g \text{ a.e. in } \Omega\rbrace,\quad V:=H^1_0(\Omega). -\end{equation*} -This set takes care of the third and fifth conditions above (the boundary -values and the complementarity condition). - -Consider now the minimizer $u\in G$ of $E$ and any other function $v\in -G$. Then the function -\begin{equation*} - F(\varepsilon) := E(u+\varepsilon(v-u)),\quad\varepsilon\in\left[0,1\right], -\end{equation*} -takes its minimum at $\varepsilon = 0$ (because $u$ is a minimizer of the -energy functional $E(\cdot)$), so that $F'(0)\geq 0$ for any choice -of $v$. Note that -$u+\varepsilon(v-u) = (1-\varepsilon)u+\varepsilon v\in G$ because of the -convexity of $G$. If we compute $F'(\varepsilon)\vert_{\varepsilon=0}$ it -yields the variational formulation we are searching for: - -\textit{Find a function $u\in G$ with} -\begin{equation*} - \left(\nabla u, \nabla(v-u)\right) \geq \left(f,v-u\right) \quad \forall v\in G. -\end{equation*} - -This is the typical form of variational inequalities, where not just $v$ -appears in the bilinear form but in fact $v-u$. The reason is this: if $u$ is -not constrained, then we can find test functions $v$ in $G$ so that $v-u$ can have -any sign. By choosing test functions $v_1,v_2$ so that $v_1-u = -(v_2-u)$ it -follows that the inequality can only hold for both $v_1$ and $v_2$ if the two -sides are in fact equal, i.e., we obtain a variational equality. - -On the other hand, if $u=g$ then $G$ only allows test functions $v$ so that in fact -$v-u\ge 0$. This means that we can't test the equation with both $v-u$ and -$-(v-u)$ as above, and so we can no longer conclude that the two sides are in -fact equal. Thus, this mimicks the way we have discussed the complementarity -condition above. - - - -\section{Formulation as a saddle point problem} - -The variational inequality above is awkward to work with. We would therefore -like to reformulate it as an equivalent saddle point problem. Set $V:=H^1_0(\Omega)$. -We introduce a Lagrange multiplier $\lambda$ and the convex cone $K\subset V'$, $V'$ -dual space of $V$, $K:=\{\mu\in V': \langle\mu,v\rangle\geq 0,\quad \forall -v\in V, v \le 0 \}$ of -Lagrange multipliers, where $\langle\cdot,\cdot\rangle$ denotes the duality -pairing between $V'$ and $V$. Intuitively, $K$ is the cone of all "non-positive -functions", except that $K\subset (H_0^1)'$ 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*} - a(u,v) + b(v,\lambda) &= f(v),\quad &&v\in V\\ - b(u,\mu - \lambda) &\leq \langle g,\mu - \lambda\rangle,\quad&&\mu\in K, -\end{align*} -\textit{with} -\begin{align*} - a(u,v) &:= \left(\nabla u, \nabla v\right),\quad &&u,v\in V\\ - b(u,\mu) &:= \langle g-u,\mu\rangle,\quad &&u\in V,\quad\mu\in V'. -\end{align*} -In other words, we can consider $\lambda$ as the negative of the additional, positive force that the -obstacle exerts on the membrane. The inequality in the second line of the -statement above only appears to have the wrong sign because we have -$\mu-\lambda<0$ at points where $\lambda=0$, given the definition of $K$. - -The existence and uniqueness of $(u,\lambda)\in V\times K$ of this saddle -point problem has been stated in Grossmann and Roos: Numerical treatment of -partial differential equations, Springer-Verlag, Heidelberg-Berlin, 2007, 596 -pages, ISBN 978-3-540-71582-5. - - - -\section{Active Set methods to solve the saddle point problem} - -There are different methods to solve the variational inequality. As one -possibility you can understand the saddle point problem as a convex quadratic program (QP) with -inequality constraints. - -To get there, let us assume that we discretize both $u$ and $\lambda$ with the -same finite element space, for example the usual $Q_k$ spaces. We would then -get the equations -\begin{eqnarray*} - &A U + B\Lambda = F,&\\ - &[BU-G]_i \geq 0, \quad \Lambda_i \leq 0,\quad \Lambda_i[BU-G]_i = 0 -\qquad \forall i.& -\end{eqnarray*} -where $B$ is the mass matrix on the chosen finite element space and the -indices $i$ above are for all degrees of freedom in the set $\cal S$ of degrees of -freedom located in the interior of the domain -(we have Dirichlet conditions on the perimeter). However, we -can make our life simpler if we use a particular quadrature rule when -assembling all terms that yield this mass matrix, namely a quadrature formula -where quadrature points are only located at the interpolation points at -which shape functions are defined; since all but one shape function are zero -at these locations, we get a diagonal mass matrix with -\begin{align*} - B_{ii} = \int_\Omega \varphi_i(\mathbf x)^2\ \textrm{d}x, - \qquad - B_{ij}=0 \ \text{for } i\neq j. -\end{align*} -To define $G$ we use the same technique as for $B$. In other words, we -define -\begin{align*} - G_{i} = \int_\Omega g_h(x) \varphi_i(\mathbf x)\ \textrm{d}x. -\end{align*} -where $g_h$ is a -suitable approximation of $g$ and $\mathbf x_i$ is the support point of the -$i$the shape function. The integral in the definition of $B_{ii}$ and $G_i$ -are then approximated by the trapezoidal rule. -With this, the equations above can be restated as -\begin{eqnarray*} - &A U + B\Lambda = F,&\\ - &U_i-B_{ii}^{-1}G_i \ge 0, \quad \Lambda_i \leq 0,\quad \Lambda_i[U_i-B_{ii}^{-1}G_i] = 0 -\qquad \forall i\in{\cal S}.& -\end{eqnarray*} - -Now we define for each degree of freedom $i$ the function -\begin{equation*} - C([BU]_i,\Lambda_i):=-\Lambda_i + \min\lbrace 0, \Lambda_i + c([BU]_i - G_i) \rbrace, -\end{equation*} -with some $c>0$. (In this program we choose $c = 100$. It is a kind of a -penalty parameter which depends on the problem itself and needs to be chosen -large enough; for example there is no convergence for $c = 1$ using the -current program if we use 7 global refinements.) - -After some headscratching one can then convince oneself that the inequalities -above can equivalently be rewritten as -\begin{equation*} - C([BU]_i,\Lambda_i) = 0, \qquad \forall i\in{\cal S}. -\end{equation*} -The primal-dual active set strategy we will use here is an iterative scheme which is based on -this condition to predict the next active and inactive sets $\mathcal{A}_k$ and -$\mathcal{F}_k$ (that is, those complementary sets of indices $i$ for which -$U_i$ is either equal to or not equal to the value of the obstacle -$B^{-1}G$). For a more in depth treatment of this approach, see 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. - -\section{The primal-dual active set algorithm} - -The algorithm for the primal-dual active set method 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$. - \item [(1)] Find the primal-dual pair $(U^k,\Lambda^k)$ that satisfies - \begin{align*} - AU^k + B\Lambda^k &= F,\\ - [BU^k]_i &= G\quad&&\forall i\in\mathcal{A}_k,\\ - \Lambda_i^k &= 0\quad&&\forall i\in\mathcal{F}_k. - \end{align*} - Note that the second and third conditions imply that exactly $|S|$ unknowns - are fixed, with the first condition yielding the remaining $|S|$ equations - necessary to determine both $U$ and $\Lambda$. - \item [(2)] Define the new active and inactive sets by - \begin{equation*} - \begin{split} - \mathcal{A}_{k+1}:=\lbrace i\in\mathcal{S}:\Lambda^k_i + c([BU^k]_i - G_i)< 0\rbrace,\\ - \mathcal{F}_{k+1}:=\lbrace i\in\mathcal{S}:\Lambda^k_i + c([BU^k]_i - G_i)\geq 0\rbrace. - \end{split} - \end{equation*} - \item [(3)] If $\mathcal{A}_{k+1}=\mathcal{A}_k$ (and then, obviously, also $\mathcal{F}_{k+1}=\mathcal{F}_k$) then stop, else set $k=k+1$ and go to step (1). -\end{itemize} -The method is called "primal-dual" because it uses both primal (the -displacement $U$) as well as dual variables (the Lagrange multiplier -$\Lambda$) to determine the next active set. - -At the end of this section, let us add two observations. First, -for any primal-dual pair $(U^k,\Lambda^k)$ that satisfies these -condition, we can distinguish the following cases: -\begin{itemize} - \item [1.] $\Lambda^k_i + c([BU^k]_i - G_i) < 0$ (i active): - - Then either $[BU^k]_i