]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Attempted a first fix to step-34. Still not working. I'll finish it in a few hours.
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 22 Apr 2009 17:57:10 +0000 (17:57 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 22 Apr 2009 17:57:10 +0000 (17:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@18690 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-34/doc/intro.dox
deal.II/examples/step-34/step-34.cc

index 0d1da58479db027991dc333fbe4ca2142b96502d..35a6f5a52cefd2fef702f91c9a92aa63d51087af 100644 (file)
@@ -16,15 +16,23 @@ usually modeled by the Euler equations of fluid dynamics:
   \frac{\partial }{\partial t}\mathbf{v} + (\mathbf{v}\cdot\nabla)\mathbf{v}
   &=
   -\frac{1}{\rho}\nabla p + \mathbf{g} 
-  \qquad &\text{in } \mathbb{R}^n \backslash \Omega
+  \qquad &\text{in } \Omega
   \\
   \nabla \cdot \mathbf{v}&=0  
-  &\text{in } \mathbb{R}^n\backslash\Omega
+  &\text{in } \Omega
 \f}
+
 where the fluid density $\rho$ and the acceleration $\mathbf{g}$ due
 to external forces are given and the velocity $\mathbf{v}$ and the
-pressure $p$ are the unknowns. Here $\Omega$ is a closed bounded
-region representing the body around which the fluid moves. 
+pressure $p$ are the unknowns. Here $\bar \Omega =
+\mathbb{R}^n\backslash \Omega$ is a closed bounded region representing
+the body around which the fluid moves.  
+
+The boundary of our exterior domain is defined as $\Gamma \cup
+\Gamma_\infty = \partial \Omega$ where $\Gamma = \partial \bar \Omega$
+while $\Gamma_\infty = \lim_{r\to\infty} \partial B_r(0)$. Notice that
+the outer normal to the domain of integration is inner to the body
+around which the fluid moves.
 
 The above equations can be derived from Navier-Stokes equations
 assuming that the effects due to viscosity are negligible compared to
@@ -44,10 +52,10 @@ flow without external forces:
   (\mathbf{v}\cdot\nabla)\mathbf{v}
   &=
   -\frac{1}{\rho}\nabla p
-  \qquad &\text{in } \mathbb{R}^n \backslash \Omega
+  \qquad &\text{in } \Omega
   \\
   \nabla \cdot \mathbf{v}&=0  
-  &\text{in } \mathbb{R}^n\backslash\Omega
+  &\text{in } \Omega
 \f}
 
 
@@ -56,8 +64,8 @@ boundary conditions
 \f[
   \label{eq:boundary-conditions}
   \begin{aligned}
-    -\mathbf{n}\cdot\mathbf{v}& = 0 \qquad && \text{ on } \partial\Omega \\
-    \mathbf{v}& = \mathbf{v}_\infty && \text{ when } |\mathbf{x}| \to \infty,
+    \mathbf{n}\cdot\mathbf{v}& = 0 \qquad && \text{ on } \Gamma = \partial \bar \Omega \\
+    \mathbf{v}& = \mathbf{v}_\infty && \text{ on } \Gamma_\infty,
   \end{aligned}
 \f]
 
@@ -65,10 +73,7 @@ which is to say that the body is at rest in our coordinate systems and
 is not permeable, and that the fluid has (constant) velocity
 $\mathbf{v}_\infty$ at infinity. An alternative viewpoint is that our
 coordinate system moves along with the body whereas the background
-fluid is at rest at infinity. Notice that we define the normal
-$\mathbf{n}$ as the <i>outer</i> normal to the domain $\Omega$, which
-is the opposite of the outer normal to the integration domain. This
-explains the reason for the minus sign in the above equation.
+fluid is at rest at infinity.
 
 For both stationary and non stationary flow, the solution process
 starts by solving for the velocity in the second equation and
@@ -86,13 +91,13 @@ which we find that the same equations hold (because $\nabla\cdot
 \f[
   \label{eq:boundary-conditions-tilde}
   \begin{aligned}
-    -\mathbf{n}\cdot\mathbf{\tilde{v}}& = \mathbf{n}\cdot\mathbf{v}_\infty \qquad && \text{ on } \partial\Omega \\
-    \mathbf{\tilde{v}}& = 0 && \text{ when } |\mathbf{x}| \to \infty,
+    \mathbf{n}\cdot\mathbf{\tilde{v}}& = -\mathbf{n}\cdot\mathbf{v}_\infty \qquad && \text{ on } \Gamma \\
+    \mathbf{\tilde{v}}& = 0 && \text{ on } \Gamma_\infty,
   \end{aligned}
 \f]
 
 If we assume that the fluid is irrotational, i.e., $\nabla \times
-\mathbf{v}=0$ in $\mathbb{R}^n\backslash\Omega$, we can represent the
+\mathbf{v}=0$ in $\Omega$, we can represent the
 velocity, and consequently also the perturbation velocity, as the
 gradient of a scalar function:
 \f[
@@ -102,10 +107,12 @@ and so the second part of Euler equations above can be rewritten
 as the homogenous Laplace equation for the unknown $\phi$:
 \f{align*}
 \label{laplace}
-\Delta\phi &= 0 \qquad &&\text{in}\ \mathbb{R}^n\backslash\Omega,
+\Delta\phi &= 0 \qquad &&\text{ in } \Omega,
           \\
-          -\mathbf{n}\cdot\nabla\phi &= \mathbf{n}\cdot\mathbf{v}_\infty 
-          && \text{on}\ \partial\Omega
+          \mathbf{n}\cdot\nabla\phi &= -\mathbf{n}\cdot\mathbf{v}_\infty 
+          && \text{ on }\ \Gamma\\
+          \mathbf{n}\cdot\nabla\phi &= 0 
+          && \text{ on }\ \Gamma_\infty
 \f}
 while the momentum equation reduces to Bernoulli's equation that expresses the
 pressure $p$ as a function of the potential $\phi$:
@@ -117,119 +124,138 @@ So we can solve the problem by solving the Laplace equation for the potential.
 We will now reformulate this equation in integral form using the
  Green identity:
 \f[\label{green}
-  \int_{\mathbb{R}^n\backslash\Omega}
-  (\Delta u)v\,dx + \int_{\partial\Omega} u\frac{\partial v}{\partial \mathbf{n}} \,ds, 
+  \int_{\Omega}
+  (\Delta u)v\,dx - \int_{\partial\Omega} \frac{\partial u}{\partial \mathbf{n} }v \,ds
   = 
-  \int_{\mathbb{R}^n\backslash\Omega}
-  (\Delta v)u\,dx + \int_{\partial\Omega} \frac{\partial u}{\partial \mathbf{n} }v \,ds
-\f]
-where, again, $\mathbf{n}$ is the normal to the surface of $\Omega$ pointing
-towards the fluid (note that the equality is typically stated with
-$\mathbf{n}$ pointing <i>outward</i> from the domain of integration whereas
-in our case it is pointing <i>inward</i>; the difference is only in the sign
-of the surface integrals).
-We also recall that the following functions,
-called fundamental solutions of the Laplace equation,
-\f[
-\begin{aligned}
-  \label{eq:3}
-    G(\mathbf{y}-\mathbf{x}) = &
-    -\frac{1}{2\pi}\ln|\mathbf{y}-\mathbf{x}|  \qquad && \text{for } n=2
-    \\
-    G(\mathbf{y}-\mathbf{x}) = &
-    \frac{1}{4\pi}\frac{1}{|\mathbf{y}-\mathbf{x}|}&& \text{for } n=3,  
-\end{aligned}
+  \int_{\Omega}
+  (\Delta v)u\,dx - \int_{\partial\Omega} u\frac{\partial v}{\partial \mathbf{n}} \,ds,
 \f]
+
+where $\mathbf{n}$ is normal to the surface of $\Omega$ pointing
+outwards from fluid (note that this normal is pointing <i>inward</i>
+the body around which the fluid moves).  We also recall that the
+following functions, called fundamental solutions of the Laplace
+equation, 
+
+\f[ \begin{aligned} \label{eq:3} G(\mathbf{y}-\mathbf{x}) = &
+-\frac{1}{2\pi}\ln|\mathbf{y}-\mathbf{x}| \qquad && \text{for } n=2 \\
+G(\mathbf{y}-\mathbf{x}) = &
+\frac{1}{4\pi}\frac{1}{|\mathbf{y}-\mathbf{x}|}&& \text{for } n=3,
+\end{aligned} 
+\f] 
+
 satisfy in a distributional sense the equation: 
-\f[
-  \Delta_y G(\mathbf{y}-\mathbf{x}) = \delta(\mathbf{y}-\mathbf{x}),
-\f]
+
+\f[ 
+\Delta_y G(\mathbf{y}-\mathbf{x}) = \delta(\mathbf{y}-\mathbf{x}),
+\f] 
+
 where the derivative is done in the variable $\mathbf{y}$.
 
 If we substitute $u$ and $v$ in the Green identity with the solution
 $\phi$ and with the fundamental solution of the Laplace equation
-respectively, as long as $\mathbf{x}$ is chosen in the region
-$\mathbb{R}^n\backslash\Omega$, we obtain:
+respectively, as long as $\mathbf{x}$ is chosen inside $\Omega$, we
+obtain:
+
 \f[
-  \phi(\mathbf{x})  
+  \phi(\mathbf{x}) - 
   \int_{\partial \Omega}\frac{\partial G(\mathbf{y}-\mathbf{x})}{\partial \mathbf{n}_y}\phi(\mathbf{y})\,ds_y 
   =
-  \int_{\partial \Omega}G(\mathbf{y}-\mathbf{x})\frac{\partial \phi}{\partial \mathbf{n}_y}(\mathbf{y})\,ds_y
-  \qquad \forall\mathbf{x}\in \mathbb{R}^n\backslash\Omega.
+  \int_{\partial \Omega}G(\mathbf{y}-\mathbf{x})\frac{\partial \phi}{\partial \mathbf{n}_y}(\mathbf{y})\,ds_y
+  \qquad \forall\mathbf{x}\in \Omega.
 \f]
-We can write this more compactly using the so-called Single and Double
-Layer Potential operators:
+
+We notice here that on $\Gamma_\infty$ we have that $\nabla \phi$ is
+zero, implying that $\phi_\infty$ is constant. It is a simple fact to
+verify that this in turn implies that
+
+\f[
+   \int_{\Gamma_\infty}\frac{\partial G(\mathbf{y}-\mathbf{x})}
+   {\partial \mathbf{n}_y}\phi(\mathbf{y})\,ds_y =
+    \phi_\infty \qquad \forall \mathbf{x} \in \mathbb{R}^n,
+\f]
+
+both in two and in three dimensions.
+
+We can then write the final problem in a more compatc way using the
+so-called Single and Double Layer Potential operators, which are
+integral operators defined only on the finite part of the
+boundary of $\Omega$:
+
 \f[\label{integral}
-  \phi(\mathbf{x}) + (D\phi)(\mathbf{x}) =
-  \left(S \frac{\partial \phi}{\partial n_y}\right)(\mathbf{x})
-  \qquad \forall\mathbf{x}\in \mathbb{R}^n\backslash\Omega.
+  \phi(\mathbf{x}) -\phi_\infty - (D\phi)(\mathbf{x}) =
+  -\left(S \frac{\partial \phi}{\partial n_y}\right)(\mathbf{x})
+  \qquad \forall\mathbf{x}\in \Omega.
 \f]
+
 (The name of these operators comes from the fact that they describe the
 electric potential in $\mathbb{R}^n$ due to a single thin sheet of charges
 along a surface, and due to a double sheet of charges and anti-charges along
 the surface, respectively.)
 
-In our case, we know the Neumann values of $\phi$ on the boundary:
-$-\mathbf{n}\cdot\nabla\phi = \mathbf{n}\cdot\mathbf{v}_\infty$.
+In our case, we know the Neumann values of $\phi$ on the boundary $\Gamma$:
+$\mathbf{n}\cdot\nabla\phi = -\mathbf{n}\cdot\mathbf{v}_\infty$.
 Consequently,
+
 \f[
-  \phi(\mathbf{x}) + (D\phi)(\mathbf{x}) = 
+  \phi(\mathbf{x}) -\phi_\infty - (D\phi)(\mathbf{x}) = 
    \left(S[\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
-   \qquad \forall\mathbf{x} \in \mathbb{R}^n\backslash\Omega.
+     \qquad \forall\mathbf{x} \in \Omega.
 \f]
-If we take the limit for $\mathbf{x}$ tending to $\partial\Omega$ of
+
+If we take the limit for $\mathbf{x}$ going to $\Gamma$ of
 the above equation, using well known properties of the single and double layer
-operators, we obtain an equation for $\phi$ just on the boundary of
-$\Omega$:
+operators, we obtain an equation for $\phi$ just on $\Gamma$:
 
 \f[\label{SD}
-  \alpha(\mathbf{x})\phi(\mathbf{x}) + (D\phi)(\mathbf{x}) = 
+  \alpha(\mathbf{x})\phi(\mathbf{x}) -\phi_\infty - (D\phi)(\mathbf{x}) = 
   \left(S [\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
-  \quad \mathbf{x}\in \partial\Omega,
+    \quad \mathbf{x}\in \partial\Omega,
 \f]
 which is the integral formulation we were looking for, where the
-quantity $\alpha(\mathbf{x})$ is the fraction of angle or solid
+quantity $\alpha(\mathbf{x})$ is the fraction of exterior angle or solid
 angle by which the point $\mathbf{x}$ sees the domain $\Omega$. In particular,
-at points $\mathbf{x}$ where the boundary $\partial\Omega$ is differentiable
+at points $\mathbf{x}$ where the boundary $\Gamma$ is differentiable
 (i.e. smooth) we have $\alpha(\mathbf{x})=\frac 12$, but the value may be
 smaller or larger at points where the boundary has a corner or an edge.
 
 Substituting the single and double layer operators we get:
 \f[               
-  \alpha(\mathbf{x}) \phi(\mathbf{x}) 
-  - \frac{1}{2\pi}\int_{\partial \Omega}  \frac{
+  \alpha(\mathbf{x}) \phi(\mathbf{x}) -\phi_\infty
+  + \frac{1}{2\pi}\int_{\Gamma}  \frac{
   (\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y  }{ |\mathbf{y}-\mathbf{x}|^2 }
   \phi(\mathbf{x}) \,ds_y 
   =
-  -\frac{1}{2\pi}\int_{\partial \Omega}  \ln|\mathbf{y}-\mathbf{x}| \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
+   -\frac{1}{2\pi}\int_{\Gamma}  \ln|\mathbf{y}-\mathbf{x}| \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
 \f]                 
 for two dimensional flows and
 \f[               
-  \alpha(\mathbf{x}) \phi(\mathbf{x}) 
-   - \frac{1}{4\pi}\int_{\partial \Omega} \frac{ (\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y  }{ |\mathbf{y}-\mathbf{x}|^3 }\phi(\mathbf{y})\,ds_y
+  \alpha(\mathbf{x}) \phi(\mathbf{x})  -\phi_\infty
+   + \frac{1}{4\pi}\int_{\Gamma} \frac{ (\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y  }{ |\mathbf{y}-\mathbf{x}|^3 }\phi(\mathbf{y})\,ds_y
   =
-  \frac{1}{4\pi}\int_{\partial \Omega} \frac{1}{|\mathbf{y}-\mathbf{x}|} \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
+  \frac{1}{4\pi}\int_{\Gamma} \frac{1}{|\mathbf{y}-\mathbf{x}|} \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
 \f]                 
 for three dimensional flows, where the normal derivatives of the fundamental
 solutions have been written in a form that makes computation easier. In either
 case, $\phi$ is the solution of an integral equation posed entirely on the
-boundary since both $\mathbf{x},\mathbf{y}\in\partial\Omega$.
+boundary since both $\mathbf{x},\mathbf{y}\in\Gamma$.
 
 Notice that the fraction of angle (in 2d) or solid angle (in 3d)
 $\alpha(\mathbf{x})$ by which the point $\mathbf{x}$ sees the domain
 $\Omega$ can be defined using the double layer potential itself:
 \f[
-\alpha(\mathbf{x}) := 
-\frac{1}{2(n-1)\pi}\int_{\partial \Omega} \frac{ (\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y  }
+\alpha(\mathbf{x}) := 1 -
+\frac{1}{2(n-1)\pi}\int_{\Gamma} \frac{ (\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y  }
 { |\mathbf{y}-\mathbf{x}|^{n} }\phi(\mathbf{y})\,ds_y =
--\int_{\partial \Omega} \frac{ \partial G(\mathbf{y}-\mathbf{x}) }{\partial \mathbf{n}_y} \, ds_y.
+1 + \int_{\Gamma} \frac{ \partial G(\mathbf{y}-\mathbf{x}) }{\partial \mathbf{n}_y} \, ds_y.
 \f]
 
 The reason why this is possible can be understood if we consider the
 fact that the solution of a pure Neumann problem is known up to an
-arbitrary constant $c$, which means that, if we set the Neumann data
-to be zero, then any constant $\phi$ will be a solution, giving us an
-the explicit expression above for $\alpha(\mathbf{x})$.
+arbitrary constant $\phi_\infty$, which means that, if we set the
+Neumann data to be zero, then any constant $\phi = \phi_\infty$ will
+be a solution, giving us an the explicit expression above for
+$\alpha(\mathbf{x})$.
 
 While this example program is really only focused on the solution of the
 boundary integral equation, in a realistic setup one would still need to solve
@@ -238,9 +264,8 @@ $\phi(\mathbf{x})$ for all $\mathbf{x}\in\partial\Omega$. In the next step, we
 can compute (analytically, if we want) the solution $\phi(\mathbf{x})$ in all
 of $\mathbb{R}^n\backslash\Omega$. To this end, recall that we had
 \f[
-  \phi(\mathbf{x}) 
+  \phi(\mathbf{x}) - \phi_\infty
   =
-  -
   (D\phi)(\mathbf{x})
   +
   \left(S[\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
@@ -252,8 +277,8 @@ $\phi$ on the boundary we have just computed). Finally, we can then recover
 the velocity as $\mathbf{\tilde v}=\nabla \phi$. 
 
 Notice that the evaluation of the above formula for $\mathbf{x} \in
-\Omega$ should yield zero as a result, since the integration of the
-the Dirac delta $\delta(\mathbf{x})$ in the domain
+\Omega$ should yield $\phi_\infty$ as a result, since the integration
+of the the Dirac delta $\delta(\mathbf{x})$ in the domain
 $\mathbb{R}^n\backslash\Omega$ is always zero by definition.
 
 As a final test, let us verify that this velocity indeed satisfies the
@@ -413,7 +438,7 @@ This is exactly Bernoulli's law mentioned above.
 
 <h3>The numerical approximation</h3>
 
-Numerical approximations of %Boundary Integral Equations (BIE) are commonly
+Numerical approximations of Boundary Integral Equations (BIE) are commonly
 referred to as the boundary element method or panel method (the latter
 expression being used mostly in the computational fluid dynamics community).
 The goal of the following test problem is to solve the integral
@@ -424,7 +449,7 @@ treat boundary element problems almost as easily as finite element
 problems using the deal.II library.
 
 To this end, let $\mathcal{T}_h = \bigcup_i K_i$ be a subdivision of the
-manifold $\Gamma = \partial \Omega$ into $M$ line segments if $n=2$, or $M$
+manifold $\Gamma$ into $M$ line segments if $n=2$, or $M$
 quadrilaterals if $n=3$. We will call each individual segment or
 quadrilateral an <i>element</i> or <i>cell</i>, independently of the
 dimension $n$ of the surrounding space $\mathbb{R}^n$.
@@ -470,16 +495,16 @@ Given the datum $\mathbf{v}_\infty$, find a function $\phi_h$ in $V_h$
 such that the following $n\_dofs$ equations are satisfied:
 
 \f{align*}
-    \alpha(\mathbf{x}_i) \phi_h(\mathbf{x}_i)  
-    + \int_{\Gamma_y} \frac{ \partial G(\mathbf{y}-\mathbf{x}_i)}{\partial\mathbf{n}_y }
+    \alpha(\mathbf{x}_i) \phi_h(\mathbf{x}_i)
+    - \int_{\Gamma_y} \frac{ \partial G(\mathbf{y}-\mathbf{x}_i)}{\partial\mathbf{n}_y }
     \phi_h(\mathbf{y}) \,ds_y = 
     \int_{\Gamma_y} G(\mathbf{y}-\mathbf{x}_i) \, 
-    \mathbf{n}_y\cdot\mathbf{v_\infty} \,ds_y 
+      \mathbf{n}_y\cdot\mathbf{v_\infty} \,ds_y 
     ,
 \f}
 where the quantity $\alpha(\mathbf{x}_i)$ is the fraction of (solid)
 angle by which the point $\mathbf{x}_i$ sees the domain $\Omega$, as 
-explained above.
+explained above, and we arbitrarily set  $\phi_\infty$ to be zero.
 If the support points $\mathbf{x}_i$ are chosen appropriately, then the
 problem can be written as the following linear system:
 \f[
@@ -492,11 +517,11 @@ where
 \mathbf{A}_{ij}&= 
 \alpha(\mathbf{x}_i) \psi_j(\mathbf{x}_i)
 =
--\int_\Gamma 
-\frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y 
-\psi_j(\mathbf{x}_i)
-\\
-\mathbf{N}_{ij}&= \int_\Gamma
+ 1+\int_\Gamma 
+ \frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y 
+ \psi_j(\mathbf{x}_i)
+ \\
+ \mathbf{N}_{ij}&= -\int_\Gamma
   \frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}
   \psi_j(\mathbf{y}) \,ds_y 
 \\
@@ -516,10 +541,10 @@ $\mathbf{A}$ is diagonal with entries
 \f[
   \mathbf{A}_{ii}
   = 
-  -\int_\Gamma 
+  1+\int_\Gamma 
   \frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y 
   =
-  -\sum_j N_{ij},
+  1-\sum_j N_{ij},
 \f]
 where we have used that $\sum_j \psi_j(\mathbf{y})=1$ for the usual Lagrange
 elements. 
@@ -538,7 +563,7 @@ element $\hat K$.
 Before discussing specifics of this integration in the
 next section, let us point out that the matrix $\mathbf{A}+\mathbf{N}$
 is rank deficient. This is mostly easily seen by realizing that
-$\mathbf{A}=-(\mathbf{N}\mathbf{e})\mathbf{e}^T$ where $\mathbf{e}$ is a
+$\mathbf{A}=1-(\mathbf{N}\mathbf{e})\mathbf{e}^T$ where $\mathbf{e}$ is a
 vector of all ones. Consequently,
 $\mathbf{A}+\mathbf{N} =
 \mathbf{N}(\mathbf{I}-\mathbf{e}\mathbf{e}^T)$. Even if $\mathbf{N}$
index 071ed9b29bc1c59e211f4e30b38830b43c896006..c4e09fefad1b397449ba6074b39b35f74c606ae4 100644 (file)
 #include <base/utilities.h>
 
 #include <lac/full_matrix.h>
+#include <lac/matrix_lib.h>
 #include <lac/vector.h>
-#include <lac/sparse_direct.h>
+#include <lac/solver_control.h>
+#include <lac/solver_gmres.h>
+#include <lac/precondition.h>
 
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
@@ -85,15 +88,15 @@ namespace LaplaceKernel
   {
     switch(dim)
       {
-       case 2:
-           return (-std::log(R.norm()) / (2*numbers::PI) );
-
-       case 3:
-           return (1./( R.norm()*4*numbers::PI ) );
-
-       default:
-           Assert(false, ExcInternalError());
-           return 0.;
+      case 2:
+       return (-std::log(R.norm()) / (2*numbers::PI) );
+       
+      case 3:
+       return (1./( R.norm()*4*numbers::PI ) );
+       
+      default:
+       Assert(false, ExcInternalError());
+       return 0.;
       }
   }
         
@@ -104,14 +107,14 @@ namespace LaplaceKernel
   {
     switch(dim)
       {
-       case 2:
-           return R / ( -2*numbers::PI * R.square());
-       case 3:
-           return R / ( -4*numbers::PI * R.square()*R.norm() );
-
-       default:
-           Assert(false, ExcInternalError());
-           return Point<dim>();
+      case 2:
+       return R / ( -2*numbers::PI * R.square());
+      case 3:
+       return R / ( -4*numbers::PI * R.square() * R.norm() );
+       
+      default:
+       Assert(false, ExcInternalError());
+       return Point<dim>();
       }
   }
 }
@@ -140,6 +143,8 @@ class BEMProblem
 {
   public:
     BEMProblem();
+
+    ~BEMProblem();
     
     void run();
 
@@ -185,38 +190,51 @@ class BEMProblem
                                     // are used where necessary.
     void assemble_system();
 
-                                    // Notwithstanding the fact that
-                                    // the matrix is full, we use a
-                                    // SparseMatrix object and the
-                                    // SparseDirectUMFPACK solver,
-                                    // since in our experience it
-                                    // works better than using, for
-                                    // example, the LapackFullMatrix
-                                    // class. Of course, using a
-                                    // SparseMatrix object to store
-                                    // the matrix is wasteful, but at
-                                    // least for the moment that is
-                                    // all the SparseDirectUMFPACK
-                                    // class can deal with.
+                                    // There are two options for the
+                                    // solution of this problem. The
+                                    // first is to use a direct
+                                    // solver, and the second is to
+                                    // use an iterative solver. We
+                                    // opt for the second option.
                                     //
-                                    // An alternative approach would
-                                    // be the use of the GMRES
-                                    // method; however the
-                                    // construction of an efficient
-                                    // preconditioner for boundary
-                                    // element methods is not a
-                                    // trivial issue, and we won't
-                                    // treat this problem here.
+                                    // The matrix that we assemble is
+                                    // not symmetric, and we opt to
+                                    // use the GMRES method; however
+                                    // the construction of an
+                                    // efficient preconditioner for
+                                    // boundary element methods is
+                                    // not a trivial issue. Here we
+                                    // use a non preconditioned GMRES
+                                    // solver. The options for the
+                                    // iterative solver, such as the
+                                    // tolerance, the maximum number
+                                    // of iterations, are selected
+                                    // through the parameter file.
                                     //
-                                    // Moreover, we should notice
-                                    // that the solution we will
-                                    // obtain will only be unique up
-                                    // to an additive constant. This
+                                    // There is a catch, however. The
+                                    // iterative solver has some
+                                    // difficulties in treating our
+                                    // matrix, because of its special
+                                    // structure. The solution we are
+                                    // trying to obtain will only be
+                                    // unique up to an additive
+                                    // constant: unless we eliminate
+                                    // the constants from the
+                                    // computation of the residuals,
+                                    // the GMRES solver will not be
+                                    // able to find a solution. This
                                     // is taken care of in the
                                     // <code>solve_system()</code>
-                                    // method, which filters out the
-                                    // mean value of the solution at
-                                    // the end of the computation.
+                                    // method, which constructs a
+                                    // ProductMatrix between the
+                                    // system matrix and a
+                                    // MeanValueFilter class to
+                                    // eliminate the mean value at
+                                    // each iteration of GMRES. The
+                                    // solution we obtain is already
+                                    // with zero mean value, and the
+                                    // solver converges very quickly,
+                                    // even without a preconditioner.
     void solve_system();
 
                                     // Once we obtained the solution,
@@ -339,18 +357,35 @@ class BEMProblem
                                     // system might be solved by
                                     // direct LU decomposition, or by
                                     // iterative methods. In this
-                                    // example we use the
-                                    // SparseDirectUMFPACK solver,
-                                    // applied to a "fake" sparse
-                                    // matrix (a sparse matrix will
-                                    // all entries different from
-                                    // zero). We found that this
-                                    // method is faster than using a
-                                    // LapackFullMatrix object.
-
-    SparsityPattern             sparsity;
-    SparseMatrix<double>        system_matrix;    
-    Vector<double>              system_rhs;
+                                    // example we use an
+                                    // unpreconditioned GMRES
+                                    // method. Building a
+                                    // preconditioner for BEM method
+                                    // is non trivial, and we don't
+                                    // treat this subject here.
+                                    //
+                                    // Notice, moreover, that
+                                    // FullMatrix objects require
+                                    // their dimensions to be passed
+                                    // at construction time. We don't
+                                    // know yet what dimension the
+                                    // matrix will have, and we
+                                    // instantiate only a
+                                    // SmartPointer, which will later
+                                    // be initialized to a new
+                                    // FullMatrix. Since we construct
+                                    // the actual matrix using the
+                                    // new operator, we have to take
+                                    // care of freeing the used
+                                    // memory whenever we exit or
+                                    // whenever we refine the
+                                    // triangulation. This is done
+                                    // both in the distructor and in
+                                    // the refine_and_resize()
+                                    // function.
+
+    SmartPointer<FullMatrix<double> >        system_matrix;    
+    Vector<double>                           system_rhs;
 
                                     // The next two variables will
                                     // denote the solution $\phi$ as
@@ -424,6 +459,8 @@ class BEMProblem
 
     std_cxx0x::shared_ptr<Quadrature<dim-1> > quadrature;
     unsigned int singular_quadrature_order;
+    
+    SolverControl solver_control;
 
     unsigned int n_cycles;
     unsigned int external_refinement;
@@ -469,6 +506,17 @@ BEMProblem<dim>::BEMProblem()
                wind(dim)
 {}
 
+
+template <int dim>
+BEMProblem<dim>::~BEMProblem() {
+  if(system_matrix != 0) {
+      FullMatrix<double> * ptr = system_matrix;
+      system_matrix = 0;
+      delete ptr;
+  }
+}
+
+
 template <int dim> 
 void BEMProblem<dim>::read_parameters (const std::string &filename)
 {
@@ -538,6 +586,19 @@ void BEMProblem<dim>::read_parameters (const std::string &filename)
                                   // us to have only one parameter
                                   // file for both 2 and 3
                                   // dimensional problems.
+                                  //
+                                  // Notice that from a mathematical
+                                  // point of view, the wind function
+                                  // on the boundary should satisfy
+                                  // the condition
+                                  // $\int_{\partial\Omega}
+                                  // \mathbf{v}\cdot \mathbf{n} d
+                                  // \Gamma = 0$, for the problem to
+                                  // have a solution. If this
+                                  // condition is not satisfied, then
+                                  // no solution can be found, and
+                                  // the solver will answer
+                                  // erratically. 
   prm.enter_subsection("Wind function 2d");
   {
     Functions::ParsedFunction<2>::declare_parameters(prm, 2);
@@ -551,21 +612,30 @@ void BEMProblem<dim>::read_parameters (const std::string &filename)
     prm.set("Function expression", "1; 1; 1");
   }
   prm.leave_subsection();
-
+  
   prm.enter_subsection("Exact solution 2d");
   {
     Functions::ParsedFunction<2>::declare_parameters(prm);
-    prm.set("Function expression", "x+y");
+    prm.set("Function expression", "-x-y");
   }
   prm.leave_subsection();
 
   prm.enter_subsection("Exact solution 3d");
   {
     Functions::ParsedFunction<3>::declare_parameters(prm);
-    prm.set("Function expression", "x+y+z");
+    prm.set("Function expression", "-x-y-z");
   }
   prm.leave_subsection();
 
+
+                                  // In the solver section, we set
+                                  // all SolverControl
+                                  // parameters. The object will be
+                                  // fed to the GMRES solver.
+  prm.enter_subsection("Solver");
+  SolverControl::declare_parameters(prm);
+  prm.leave_subsection();
+
                                   // After declaring all these
                                   // parameters to the
                                   // ParameterHandler object, let's
@@ -604,6 +674,11 @@ void BEMProblem<dim>::read_parameters (const std::string &filename)
   }
   prm.leave_subsection();
 
+  prm.enter_subsection("Solver");
+  solver_control.parse_parameters(prm);
+  prm.leave_subsection();
+
+
                                   // Finally, here's another example
                                   // of how to use parameter files in
                                   // dimension independent
@@ -707,20 +782,10 @@ void BEMProblem<dim>::read_domain()
                                 // freedom, and resizes matrices and
                                 // vectors.
                                 //
-                                // Note that the matrix is a full
-                                // matrix and that consequently we
-                                // have to build a sparsity pattern
-                                // that contains every single
-                                // entry. Notwithstanding this fact,
-                                // the SparseMatrix class coupled
-                                // with the SparseDirectUMFPACK
-                                // solver are still faster than
-                                // Lapack solvers for full
-                                // matrices. The drawback is that we
-                                // need to assemble a full
-                                // SparsityPattern, which is not the
-                                // most efficient way to store a full
-                                // matrix.
+                                // Note that we need to handle the
+                                // FullMatrix system_matrix
+                                // carefully, since it requires the
+                                // size.
 
 template <int dim>
 void BEMProblem<dim>::refine_and_resize()
@@ -730,14 +795,13 @@ void BEMProblem<dim>::refine_and_resize()
   dh.distribute_dofs(fe);
     
   const unsigned int n_dofs =  dh.n_dofs();
-    
-  system_matrix.clear();
-  sparsity.reinit(n_dofs, n_dofs, n_dofs);
-  for(unsigned int i=0; i<n_dofs;++i)
-    for(unsigned int j=0; j<n_dofs; ++j)
-      sparsity.add(i,j);
-  sparsity.compress();
-  system_matrix.reinit(sparsity);
+  
+  if(system_matrix != 0) {
+      FullMatrix<double> * ptr = system_matrix;
+      system_matrix = 0;
+      delete ptr;
+  }
+  system_matrix = new FullMatrix<double> (n_dofs, n_dofs);
     
   system_rhs.reinit(n_dofs);
   phi.reinit(n_dofs);
@@ -1206,13 +1270,13 @@ void BEMProblem<dim>::assemble_system()
              delete singular_quadrature;
          }
             
-                                          // Finally, we need to add the
-                                          // contributions of the current
-                                          // cell to the global matrix:
+                                          // Finally, we need to add
+                                          // the contributions of the
+                                          // current cell to the
+                                          // global matrix.
          for(unsigned int j=0; j<fe.dofs_per_cell; ++j) 
-           system_matrix.add(i,
-                             local_dof_indices[j],
-                             local_matrix_row_i(j));
+             (*system_matrix)(i,local_dof_indices[j]) 
+                 += local_matrix_row_i(j);
        }
     }
 
@@ -1242,39 +1306,47 @@ void BEMProblem<dim>::assemble_system()
                                   // matrix:
   Vector<double> ones(dh.n_dofs());
   ones.add(-1.);
-
-  system_matrix.vmult(alpha, ones);
+  
+  system_matrix->vmult(alpha, ones);
+  alpha.add(1);
   for(unsigned int i = 0; i<dh.n_dofs(); ++i)
-    system_matrix.add(i,i,alpha(i));
+      (*system_matrix)(i,i) +=  alpha(i);
 }
 
 
                                 // @sect4{BEMProblem::solve_system}
 
                                 // The next function simply solves
-                                // the linear system. As described,
-                                // we use the SparseDirectUMFPACK
-                                // direct solver to compute the
-                                // inverse of the matrix (in reality
-                                // it only produces an LU
-                                // decomposition) and then apply this
-                                // inverse to the right hand side to
-                                // yield the solution.
+                                // the linear system.
                                 //
                                 // As mentioned in the introduction,
-                                // the solution is only known up to a
-                                // constant potential. We solve this
-                                // issue by subtracting the mean
-                                // value of the vector from each
-                                // vector entry to normalize it.
+                                // the fact that the system matrix
+                                // has a non null kernel, requires us
+                                // to be careful in case we wish to
+                                // use an iterative solver. To
+                                // address this issue, we use two new
+                                // instruments of the library: the
+                                // MeanValueFilter class, and the
+                                // ProductMatrix class. The
+                                // MeanValueFilter has the same
+                                // interface of a matrix, with the
+                                // effect of subtracting the mean
+                                // value to source vector. We cascade
+                                // this operator with the system
+                                // matrix, and we obtain a matrix
+                                // whose result is renormalized to a
+                                // zero mean value Vector. This
+                                // object is then passed to a GMRES
+                                // solver.
 template <int dim>
 void BEMProblem<dim>::solve_system()
 {
-  SparseDirectUMFPACK inverse_matrix;
-  inverse_matrix.initialize (system_matrix);
-  inverse_matrix.vmult (phi, system_rhs);
-
-  phi.add(-phi.mean_value());
+    PrimitiveVectorMemory<Vector<double> > mem;
+    MeanValueFilter filter;
+    ProductMatrix<Vector<double> > system(filter, *system_matrix, mem);
+       
+    SolverGMRES<Vector<double> > solver (solver_control);
+    solver.solve (system, phi, system_rhs, PreconditionIdentity());
 }
 
 

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.