]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Convert intro from latex to doxygen format.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 25 Feb 2012 10:33:25 +0000 (10:33 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 25 Feb 2012 10:33:25 +0000 (10:33 +0000)
git-svn-id: https://svn.dealii.org/trunk@25170 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-41/doc/intro.dox
deal.II/examples/step-41/doc/step-41-doc.tex [deleted file]

index 3438be49a4d6c6d715a8ff44f2685cb5d026bbe7..2b5601442f993b019ba423a4469f31dca8c4b60a 100644 (file)
@@ -6,3 +6,413 @@ This material is based upon work partly supported by the ...
 </i>
 
 
+<h3>Introduction</h3>
+
+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:
+
+<table align="center" class="tutorial" cellspacing="3" cellpadding="3">
+  <tr>
+    <td align="center">
+        @image html "step-41.displacement.png"
+    </td>
+    <td align="center">
+        @image html "step-41.active-set.png"
+    </td>
+  </tr>
+</table>
+
+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.
+
+
+<h3>Classical formulation</h3>
+
+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).
+
+
+<h3>Derivation of the variational inequality</h3>
+
+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:
+
+<i>Find a function $u\in G$ with</i>
+@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.
+
+
+
+<h3>Formulation as a saddle point problem</h3>
+
+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:
+
+<i>Find $u\in V$ and $\lambda\in K$ such that</i>
+@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}
+<i>with</i>
+@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.
+
+
+
+<h3>Active Set methods to solve the saddle point problem</h3>
+
+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.
+
+<h3>The primal-dual active set algorithm</h3>
+
+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):
+  <br>
+  Then either $[BU^k]_i<G_i$ and $\Lambda^k_i=0$ (penetration) or $\Lambda^k_i<0$ and $[BU^k]_i=G_i$ (pressing load).
+- [2.] $\Lambda^k_i + c([BU^k]_i - G_i)\geq 0$ (i inactive):
+  <br>
+  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.
+
+
+<h3>Implementation</h3>
+
+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, <code>assemble_mass_matrix_diagonal</code> and
+  <code>update_solution_and_constraints</code>.
+- 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 (file)
index e43cdaa..0000000
+++ /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<G_i$ and $\Lambda^k_i=0$ (penetration) or $\Lambda^k_i<0$ and $[BU^k]_i=G_i$ (pressing load).
- \item [2.] $\Lambda^k_i + c([BU^k]_i - G_i)\geq 0$ (i inactive):
-
-  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).
-\end{itemize}
-
-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
-\begin{eqnarray*}
- &A U + B\Lambda = F,&\\
- &C([BU-G]_i, \Lambda_i) = 0,
-\qquad \forall i.&
-\end{eqnarray*}
-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
-\begin{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}
-\end{equation*}
-\begin{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}
-\end{equation*}
-This suggest a semismooth Newton step of the form
-\begin{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},
-\end{equation*}
-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
-\begin{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}.
-\end{equation*}
-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$:
-\begin{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}.
-\end{equation*}
-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
-\begin{equation*}
- \Lambda^k_S = B^{-1}\left(f_{\mathcal{S}} - A_{\mathcal{S}}U^k_{\mathcal{S}}\right).
-\end{equation*}
-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:
-\begin{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}.
-\end{equation*}
-Fortunately, this form is easy to arrive at: we simply build the usual Laplace
-linear system
-\begin{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},
-\end{equation*}
-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.
-
-
-\section{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 are, 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.
-
-\end{document}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.