]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Enanced comments for step-34, rewrote assembly routine, added external parameter...
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 3 Mar 2009 10:58:40 +0000 (10:58 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 3 Mar 2009 10:58:40 +0000 (10:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@18443 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-34/Makefile
deal.II/examples/step-34/doc/intro.dox
deal.II/examples/step-34/parameters.prm [new file with mode: 0644]
deal.II/examples/step-34/step-34.cc

index 96fe808f621be0e7c9e42a0e5f39785bdd62c116..3dc5c73fbf9e1fb2f5921e7295ea93503379a9d4 100644 (file)
@@ -153,4 +153,3 @@ Makefile.dep: $(target).cc Makefile \
 # To make the dependencies known to `make', we finally have to include
 # them:
 include Makefile.dep
-
index 0d85ffb09c4c7bac2a8972c0e057f6d3f367574c..9eddab766370a10ef4de4c25870234180381f969 100644 (file)
@@ -1,9 +1,394 @@
 <a name="Intro"></a>
+
 <h1>Introduction</h1>
 
-<h3>The problem</h3>
+<h3> Irrotational Flow </h3>
+The motion of an inviscid fluid past a body (for example air past an
+airplane wing, or air or water past a propeller) is usually modeled by
+the Euler equations of fluid dynamics:
+
+\f[
+  \left\{
+  \begin{array}{rclr}
+  \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
+  \\
+  \nabla \cdot \mathbf{v}&=&0  
+  &\qquad \text{in } \mathbb{R}^n\backslash\Omega
+%  \\
+%  \mathbf{n}\cdot\mathbf{v}&=& 0
+%  &\qquad \text{on } \partial\Omega,
+  \end{array}\right.,
+\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. 
+
+Uniqueness of the solution is ensured by adding the boundary conditions
+\f[
+  \label{eq:boundary-conditions}
+  \begin{aligned}
+    \mathbf{n}\cdot\mathbf{v}& =0 \qquad && \text{ on } \partial\Omega \\
+    \mathbf{v}& =0 && \text{ when } |\mathbf{x}| \to \infty,
+  \end{aligned}
+\f]
+which mean that the body is not permeable, and that the fluid is
+assumed at rest at infinity.
+
+The above equations can be derived from Navier-Stokes
+equations assuming that the effects due to viscosity are negligible
+compared to those due to the pressure gradient and to the external
+forces. 
+
+For both stationary and non stationary flow, the solution process
+starts by solving for the velocity in the second equation and
+substituting in the first equation in order to find the pressure.
+
+The solution of the stationary Euler equations is typically performed
+in order to understand the behavior of the given (possibly complex)
+geometry when a prescribed motion is enforced on the system. 
+
+The first step in this process is to change the frame of reference,
+putting it on the body $\Omega$, which sees a uniform background
+velocity field, and a perturbation due to its own presence:
+\f[
+\nabla\cdot\mathbf{v} = 
+\nabla\cdot(\mathbf{v_\infty} + \mathbf{v_p}) = 
+\nabla\cdot\mathbf{v_p} = 0.
+\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
+velocity, and consequently also the perturbation velocity, as the
+gradient of a scalar function:
+\f[
+  \mathbf{v_p}=\nabla\phi,
+\f]
+and so the second part of Euler equations above can be rewritten
+as the homogenous Laplace equation for the unknown $\phi$:
+\f[\label{laplace}
+\Delta\phi = 0.
+\f]
+
+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} \frac{\partial u}{\partial \mathbf{n} }v \,ds
+  = 
+  \int_{\mathbb{R}^n\backslash\Omega}
+  (-\Delta v)u\,dx - \int_{\partial\Omega} u\frac{\partial v}{\partial \mathbf{n}} \,ds,
+\f]
+where $\mathbf{n}$ is the normal to the surface of $\Omega$ pointing
+towards the fluid. 
+
+We also remind the reader that the following functions,
+called fundamental solutions of Laplace equation:
+\f[
+\begin{aligned}
+  \label{eq:3}
+    G(\mathbf{x}-\mathbf{y}) = &
+    -\frac{1}{2\pi}\ln|\mathbf{x}-\mathbf{y}|  \qquad && \text{for } n=2
+    \\
+    G(\mathbf{x}-\mathbf{y}) = &
+    -\frac{1}{4\pi}\frac{1}{|\mathbf{x}-\mathbf{y}|}&& \text{for } n=3,  
+\end{aligned}
+\f]
+satisfy in variational sense: 
+\f[
+  -\Delta_x G(\mathbf{x}-\mathbf{y}) = \delta(\mathbf{x}-\mathbf{y}),
+\f]
+where the derivation is done in the variable $\mathbf{x}$.
+
+If we substitute $u$ and $v$ in the green identity with the solution
+$\phi$ and with the fundamental solution of Laplace equation
+respectively, we obtain:
+\f[
+  \phi(\mathbf{x})=\int_{\partial \Omega}G(\mathbf{x}-\mathbf{y})\frac{\partial \phi}{\partial \mathbf{n}_y}(\mathbf{y})\,ds_y
+  +
+  \int_{\partial \Omega}\frac{\partial G(\mathbf{x}-\mathbf{y})}{\partial \mathbf{n}_y}\phi(\mathbf{y})\,ds_y 
+  \quad \mathbf{x}\in \mathbf{R}^n\backslash\Omega,
+\f]
+that we can write more compactly using the Single and Double
+Layer Potential operators:
+\f[\label{integral}
+  \phi(\mathbf{x}) = \left(S \frac{\partial \phi}{\partial n_y}\right)(\mathbf{x})
+  + 
+  (D\phi)(\mathbf{x})
+  \quad \mathbf{x}\in \mathbf{R}^n\backslash\Omega.
+\f]
+
+It can be shown that this is equivalent to solving the homogeneous
+laplace equation with the given Neumann boundary values. Notice also
+that the last equation lets one calculate $\phi$ in any point of the
+domain once its expression on the boundary $\partial\Omega$ is known.
+If we take the limit for $\mathbf{x}$ tending to $\partial\Omega$ of
+the above equation, using well known properties of the SLP and DLP
+operators, we obtain an equation for $\phi$ just on the boundary of
+$\Omega$:
+
+\f[\label{SD}
+  \frac{1}{2}\phi(\mathbf{x}) = \left(S \frac{\partial \phi}{\partial \mathbf{n}_y}\right)(\mathbf{x})
+  + 
+  (D\phi)(\mathbf{x})
+  \quad \mathbf{x}\in \partial\Omega,
+\f]
+which is the integral formulation we were looking for. Imposing the 
+boundary conditions we obtain:
+\f[
+\mathbf{n} \cdot( \mathbf{v}_\infty + \mathbf{v}_p)=0
+\quad \Rightarrow \quad
+\mathbf{n}\cdot\mathbf{v_p}=-\mathbf{n}\cdot\mathbf{v_\infty}
+\quad \Rightarrow \quad
+\frac{\partial \phi}{\partial\mathbf{n}} = -\mathbf{n}\cdot\mathbf{v_\infty}
+\f]
+which can be substituted in the single layer potential equation to obtain:
+\f[               
+  \pi\phi(\mathbf{x})=
+  \int_{\partial \Omega}  \ln|\mathbf{x}-\mathbf{y}| \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
+  +
+  \int_{\partial \Omega}  \frac{ (\mathbf{x}-\mathbf{y})\cdot\mathbf{n}_y  }{ |\mathbf{x}-\mathbf{y}|^2 }\,ds_y 
+\f]                 
+for two dimensional flows and
+\f[               
+  2\pi\phi(\mathbf{x})=\int_{\partial \Omega} -\frac{1}{|\mathbf{x}-\mathbf{y}|} \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
+  +
+  \int_{\partial \Omega} \frac{ (\mathbf{x}-\mathbf{y})\cdot\mathbf{n}_y  }{ |\mathbf{x}-\mathbf{y}|^3 }\phi(\mathbf{y})\,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.
+
+<h3>The numerical approximation</h3>
+
+Numerical approximations of Boundary Integral Equations are commonly
+referred to as the boundary element method or panel method (the latter
+is used mostly in the computational fluid dynamics community).
+
+The goal of the following test problem is to solve the integral
+formulation of Laplace equation with Neumann boundary conditions,
+using a circle and a sphere respectively in two and three space
+dimensions, illustrating along the way the features that allow one to
+treat boundary value problems almost as easily as finite element
+problems using the deal.II library:
+
+\f[
+  \label{eq:strong-continuous}
+  \frac{1}{2}\phi(\mathbf{x})=
+  \int_{\partial \Omega}G(\mathbf{x}-\mathbf{y})\frac{\partial \phi}{\partial \mathbf{n}_y}(\mathbf{y})\,ds_y
+  +
+  \int_{\partial \Omega}\frac{\partial G(\mathbf{x}-\mathbf{y})}{\partial \mathbf{n}_y}\phi(\mathbf{y})\,ds_y 
+  \quad \mathbf{x}\in \partial \Omega.
+\f]
+
+Let $\mathcal{T}_h = \cup_{i=0}^M K_i$ be a subdivision of the
+manifold $\Gamma = \partial \Omega$ into $M$ segments if $n=2$ or $M$
+quadrilaterals if $n=3$. We will call each individual segment or
+quadrilateral an \emph{element} or \emph{cell}, independently on the
+dimension $n$ of the surrounding space $R^n$.
+
+We define the finite dimensional space $V_h$ as 
+\f[
+  \label{eq:definition-Vh}
+  V_h := \{ v \in L^2(\Gamma) \text{ s.t. } v|_{K_i} \in \mathcal{P}^0(K_i), 
+  i = 0, \dots, M\},
+\f]
+with basis functions $\psi_i(\mathbf{x}) = \chi_{K_i}(\mathbf{x})$,
+i.e., one if $\mathbf{x}$ belongs to $K_i$, and zero otherwise
+(section~\ref{sec:dofs}).  An element $\phi_h$ of $V_h$ is uniquely
+identified by the vector $\boldsymbol{\alpha}$ of its coefficients
+$\phi_i$, that is:
+\f[
+  \label{eq:definition-of-element}
+  \phi_h(\mathbf{x}) := \phi_i \psi_i(\mathbf{x}), \qquad 
+  \boldsymbol{\alpha} := \{ \phi_i \}_{i=0}^M,
+\f]
+where summation  is implied over repeated indexes.
+
+<h3> Galerkin boundary element method </h3>
+
+The usual Galerkin approach for the discretization of the above
+problem gives us the following variational formulation:
+
+Given the datum $\mathbf{v}_\infty$, find a function $\phi_h$ in $V_h$
+such that, for any $\eta$ in $V_h$ the following equation is
+satisfied:
+\f[
+  \label{eq:galerkin-continuous}
+  \begin{split}
+    \int_{\Gamma_x} \phi_h(\mathbf{x}) \eta(\mathbf{x})\,ds_x &=
+    \\
+    & \frac{1}{\pi} \int_{\Gamma_x} \int_{\Gamma_y}
+    G(\mathbf{x}-\mathbf{y}) \, \mathbf{n}_y\cdot\mathbf{v_\infty}
+    \eta(\mathbf{x}) \,ds_x\,ds_y +
+    \\
+    & \frac{1}{\pi} \int_{\Gamma_x}\int_{\Gamma_y} \frac{ \partial
+      G(\mathbf{x}-\mathbf{y})}{\partial\mathbf{n}_y }
+    \phi_h(\mathbf{y})\eta(\mathbf{x}) \,ds_x\,ds_y.
+  \end{split}
+\f]
+
+The linearity of the integral operator makes this problem equivalent 
+to solving the linear system
+\f[
+\label{eq:linear-system}
+(\mathbf{M}-\mathbf{A})\boldsymbol\alpha = \mathbf{b},
+\f]
+where
+\f[
+\begin{aligned}
+\mathbf{M}_{ij}&= |K_i|\delta_{ij}\\
+\mathbf{A}_{ij}&= \frac{1}{\pi} \int_{K_i}\int_{K_j}
+  \frac{\partial G(\mathbf{x}-\mathbf{y})}{\partial \mathbf{n}_y}
+  \psi_i(\mathbf{x})\psi_j(\mathbf{y}) \,ds_x\,ds_y 
+\\
+\mathbf{b}_i&= \frac{1}{\pi}\int_{K_i}  \int_{\Gamma_{h,y}}
+   G(\mathbf{x}-\mathbf{y})  \, \mathbf{n}_y\cdot\mathbf{v_\infty}
+  \psi_i(\mathbf{y}) \,ds_x\,ds_y.
+\end{aligned}
+\f]
+
+The computation of the entries of the matrix $\mathbf{A}$ and of the
+right hand side $\mathbf{b}$ require the evaluation of singular
+integrals on the elements $K_i$ of the triangulation $\mathcal{T}_h$.
+
+As usual in this cases, all integrations are performed on a reference
+simple domain, i.e., we assume that each element $K_i$ of
+$\mathcal{T}_h$ can be expressed as a linear (in two dimensions) or
+bi-linear (in three dimensions) transformation of the reference
+element $\hat K := [0,1]^n$, and we perform the integrations after a
+change of variables from the real element $K_i$ to the reference
+element $\hat K$.
+
+<h3> Singular Integrals in two dimension. </h3>
+
+In two dimensions it is not necessary to compute the diagonal elements
+$\mathbf{A}_{ii}$ of the system matrix, since, even if the denominator
+goes to zero when $\mathbf{x}=\mathbf{y}$, the numerator is always
+zero because $\mathbf{n}_y$ and $(\mathbf{x}-\mathbf{y})$ are
+orthogonal, and the only singular integral arises in the computation
+of $\mathbf{b}_i$ on the i-th element of $\mathcal{T}_h$:
+\f[
+  \frac{1}{\pi}
+  \int_{K_i}  \int_{K_i}
+  \ln|\mathbf{x}-\mathbf{y}| \, \mathbf{n}_y\cdot\mathbf{v_\infty}
+  \,ds_x\,ds_y \, .
+\f]
+
+Using the linear transformations above and defining
+$J=|\mathbf{V}^i_1-\mathbf{V}_0^i|$, we obtain
+\f[
+  \frac{1}{\pi}
+  \int_0^1 d\eta  \int_0^1 d\lambda \, J^2
+  (\ln|\lambda-\eta|+\ln(J)) \, 
+   \mathbf{n}_y\cdot\mathbf{v_\infty}
+  \, .
+\f]
+
+After another change of variables
+\f[
+\begin{aligned}
+\lambda-\eta=&\alpha
+\\
+\lambda+\eta=&\beta \, ,
+\end{aligned}
+\f]
+we end up with:
+\f[
+  \frac{J^2}{\pi}\ln(J) \mathbf{n}_y\cdot\mathbf{v_\infty}
+  \, + \,
+  \frac{2J^2}{\pi}
+  \int_{-1}^1 d\alpha  \int_0^2 d\beta
+  \ln|\alpha| \, \mathbf{n}_y\cdot\mathbf{v_\infty}
+   \, ,
+\f]
+which can be computed analytically with Gauss-log quadrature formulae.
+
+<h3>Singular integrals in three dimensions</h3>
+
+In three dimensions the computation of the integrals is somewhat more
+complicated. Some results in this direction can be found
+in~\cite{newman} and~\cite{morino-chen-suciu}, and even if the two
+methods are identical for low order elements, in~\cite{newman} the
+method is extended to higher order approximations, using multipole
+expansion.
+
+In this work we started with the implementation of a low order method,
+inspired by~\cite{morino-chen-suciu}.
+
+We are interested in calculating the integrals 
+\f[
+  \label{eq:slp-dlp-on-panel}
+  \begin{split}
+    S_i(x,y,z) =  & - \frac{1}{2\pi} \int_{K_i} \frac{1}{R} \mathtt{d}S \\
+    D_i(x,y,z) =  & - \frac{1}{2\pi} \int_{K_i} \frac{\mathbf{R}\cdot\mathbf{n}}{R^3} \mathtt{d}S,
+  \end{split}
+\f]
+where $R(x,y,z,\eta,\xi)$ is the distance from the point
+$\mathbf{x}=(x,y,z)$ and the point $\mathbf{y}(\eta, \xi)$ on the
+surface of the panel identified by the local coordinates $(\eta,
+\xi)$.
+
+Introducing the tangent vectors to the surface of the panels,
+\f[
+  \label{eq:tangential}
+  \mathbf{a}_1 = \frac{\partial \mathbf{y}}{\partial \eta}, 
+  \qquad \mathbf{a}_2  = \frac{\partial \mathbf{y}}{\partial \xi},
+\f]
+we can define the normal vector and the element differential area as
+\f[
+  \label{eq:normal-diff-area}
+  \mathbf{n}(\eta, \xi) = \frac{\mathbf{a}_1\times \mathbf{a}_2}% 
+  {|\mathbf{a}_1\times \mathbf{a}_2|}, \qquad
+  \mathtt{d} S = |\mathbf{a}_1\times \mathbf{a}_2| 
+  \mathtt{d} \eta \mathtt{d} \xi.
+\f]
+
+The single and double layer potential integrals on a panel can then be
+expressed analytically by rewriting them as
+\f[
+  \label{eq:kernels-by-parts}
+  \int_0^1 \int_0^1 \frac{\partial^2 I_X(\eta, \xi)}{\partial \eta \partial \xi} \mathtt{d} S, 
+\f]
+which implies
+\f[
+  \label{eq:integral-form-of-slp-and-dlp}
+  \begin{split}
+    S_i = & I_S(1,1)-I_S(1,0) + I_S(0,1) - I_S(0,0) \\
+    D_i = & I_D(1,1)-I_D(1,0) + I_D(0,1) - I_D(0,0).
+  \end{split}
+\f]
+
+The single and double layer integrals on a quadrilateral panel then
+take the form above, where the terms $I_S$ and $I_D$ are given by
+\f[
+  \label{eq:def-I-S}
+  \begin{split}
+    I_S(\eta,\xi) = - \frac{1}{2\pi} \Bigg( & -
+    \frac{(\mathbf{R}\times \mathbf{a}_1)\cdot
+      \mathbf{n}}{|\mathbf{a}_1|}
+    \sinh^{-1}\left( \frac{\mathbf{R}\cdot\mathbf{a}_1}{|\mathbf{R}\times \mathbf{a}_1|} \right) \\
+    & + \frac{(\mathbf{R}\times \mathbf{a}_2)\cdot
+      \mathbf{n}}{|\mathbf{a}_2|}
+    \sinh^{-1}\left( \frac{\mathbf{R}\cdot\mathbf{a}_2}{|\mathbf{R}\times \mathbf{a}_2|} \right) \\
+    & + \mathbf{R}\cdot \mathbf{n} \tan^{-1}
+    \left( \frac{(\mathbf{R}\times\mathbf{a}_1)\cdot(\mathbf{R}\times\mathbf{a}_2)}%
+      {R\mathbf{R}\cdot(\mathbf{a}_1\times\mathbf{a}_2)}\right) \Bigg),
+  \end{split}
+\f]
+and 
+\f[
+  \label{eq:def-I-D}
+  I_D(\eta,\xi) = \frac{1}{2\pi}  \tan^{-1}
+  \left( \frac{(\mathbf{R}\times\mathbf{a}_1)\cdot(\mathbf{R}\times\mathbf{a}_2)}%
+    {R\mathbf{R}\cdot(\mathbf{a}_1\times\mathbf{a}_2)}\right).
+\f]
+
+The resulting matrix $\mathbf{A}$ is full. Depending on its size, it
+might be convenient to use a direct solver or an iterative one.
 
-<h3>Galerkin BEM</h3>
 
 <h3>What the program does</h3>
 
diff --git a/deal.II/examples/step-34/parameters.prm b/deal.II/examples/step-34/parameters.prm
new file mode 100644 (file)
index 0000000..b499543
--- /dev/null
@@ -0,0 +1,25 @@
+# Listing of Parameters
+# ---------------------
+
+set Number of cycles = 3
+
+subsection Outer quadrature rule
+  set Quadrature order = 0
+  set Quadrature type  = midpoint
+end
+
+
+subsection Wind function
+  # Any constant used inside the function which is not a variable name.
+  set Function constants  = angle=0
+
+  # Separate vector valued expressions by ';' as ',' is used internally by the
+  # function parser.
+  set Function expression = cos(angle); sin(angle); 0
+
+  # The name of the variables as they will be used in the function, separated
+  # by ','.
+  set Variable names      = x,y,z,t
+end
+
+
index abb217697dc73e750625212ab1ea4421b298f1ca..c39ccb761d2663b6ecd3f8b838a9cba6f56c2f12 100644 (file)
@@ -1,4 +1,4 @@
-//----------------------------  bem_integration.cc  ---------------------------
+//----------------------------  step-34.cc  ---------------------------
 //    $Id: testsuite.html 13373 2006-07-13 13:12:08Z manigrasso $
 //    Version: $Name$ 
 //
@@ -11,7 +11,7 @@
 //
 //    Authors: Luca Heltai, Cataldo Manigrasso
 //
-//----------------------------  bem_integration.cc  ---------------------------
+//----------------------------  step-34.cc  ---------------------------
 
 
 // 
 
 #include <base/convergence_table.h>
 #include <base/quadrature_lib.h>
+#include <base/quadrature_selector.h>
 #include <base/table.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <dofs/dof_renumbering.h>
+#include <base/parsed_function.h>
 #include <fe/fe_dgp.h>
 #include <fe/fe_system.h>
 #include <fe/fe_tools.h>
@@ -59,27 +61,53 @@ class LaplaceKernelIntegration
 {
 public:
 
-    LaplaceKernelIntegration();
+    LaplaceKernelIntegration(FiniteElement<dim-1,dim> &fe);
     ~LaplaceKernelIntegration();
   
     void run();
     
-    void compute_SD_integral_on_cell(vector<double> &dst,
-                                    typename DoFHandler<dim,dim+1>::active_cell_iterator &cell,
-                                    const Point<dim+1> &point);
+    // This functions computes the integral of the single and double
+    // layer potentials on the cell given as a parameter, at the
+    // quadrature points @p q. In practice this function produces the objects
+    // 
+    // \f[
+    // \text{dst}_{ik0} := \int_{\text{cell}} G(y - \text[q]_k) \phi_i dy
+    // \f]
+    // 
+    // and 
+    //
+    // \f[
+    // \text{dst}_{ik1} := \int_{\text{cell}} \frac{\partial
+    // G}{\partial \textbf{n}} (y - \text[q]_k) \phi_i dy
+    // \f]
+    void compute_SD_integral_on_cell(vector<vector<vector<double> > > &dst,
+                                    typename DoFHandler<dim-1,dim>::active_cell_iterator &cell,
+                                    const vector<Point<dim> > &q);
 
 private:
+    // The following two helper functions should only be called when
+    // dim=3. If this is not the case, the default implementation is
+    // to throw an exception. When the dimension is equal to two, it
+    // is possible to compute the singular integrals using the
+    // GaussLog quadrature formulas.
+
     double term_S(const Point<3> &r,
                  const Point<3> &a1,
                  const Point<3> &a2,
                  const Point<3> &n,
-                 const double &rn_c);
+                 const double &rn_c) {
+       AssertThrow(false, ExcImpossibleInDim());
+       return 0;
+    };
 
     double term_D(const Point<3> &r,
                  const Point<3> &a1,
-                 const Point<3> &a2);
+                 const Point<3> &a2) {
+       AssertThrow(false, ExcImpossibleInDim());
+       return 0;
+    };
 
-    SmartPointer<FEValues<dim,dim+1> > fe_values;
+    SmartPointer<FEValues<dim-1,dim> > fe_values;
 };
 
 template <int dim> 
@@ -89,71 +117,74 @@ public:
     BEMProblem();
     ~BEMProblem();
     
-    /** Starts the Boundary Element Method Computation. */
+    // Read parameters.
+    void read_parameters(std::string filename);
+    
+    // Starts the Boundary Element Method Computation.
     void run();
     
-    /** Initialize mesh and vector space. */
+    // Initialize mesh and vector space. 
     void read_domain();
 
-    /** Refine and resize all vectors for the active step. */
+    // Refine and resize all vectors for the active step. 
     void refine_and_resize();
     
-    /** Assemble the two system matrices as well as the system right
-     * hands side. */
+    // Assemble the two system matrices as well as the system right
+    // hands side.
     void assemble_system();
 
-    /** Solve the system. */
+    // Solve the system. 
     void solve_system();
 
-    /** Output results for the given cycle. */
+    // Output results for the given cycle. 
     void output_results(unsigned int cycle);
     
 private:
-    /** The boundary element method triangulation. */
+    // The boundary element method triangulation. 
     Triangulation<dim-1, dim> tria;
 
-    /** The finite element space for the potential. */
-    FE_DGP<2,3> fe;
-    
-    /** The finite element space for the velocity. */
-    FESystem<2,3> fev;
+    // The finite element spaces for the potential and the velocity. 
+    FE_DGP<dim-1,dim> fe;
+    FESystem<dim-1,dim> fev;
     
-    /** The potential degrees of freedom. */
+    // Finite element space used to smoothen the potential solution
+    // (from piecewise constant to continuous piecewise quadratic)
+    FE_Q<dim-1, dim>  fe_q;
+
+    // And the relevant degrees of freedom.
     DoFHandler<dim-1,dim> dh;
-    
-    /** The velocity degrees of freedom. */
     DoFHandler<dim-1,dim> dhv;
+    DoFHandler<dim-1,dim> dhq;
 
-    /** The system matrix. This is I-C. Since the LAPACKFullMatrix
-     * does not have a reinit method, we need to work around this a
-     * little. */
+    // The system matrix. This is I-C. Since the LAPACKFullMatrix does
+    // not have a reinit method, we need to work around this a little.
     SmartPointer<LAPACKFullMatrix<double> > system_matrix;
     
-    /** Single layer potential matrix. */
-    FullMatrix<double> single_layer_matrix;
-    
-    /** Normal component of the wind field. */
-    Vector<double> Vn;
-    
-    /** System rhs. */
-    Vector<double> rhs;
-    
-    /** Potential. */
+    // The right hand side, the potential and its smoothed version
+    Vector<double> system_rhs;
     Vector<double> phi;
-    
-    /** Something else. */
-    Vector<double> vs;
+    Vector<double> smooth_phi;
+    
+    // These are the parameters that we read in from a parameter file.
+    // In particular we define the wind function and the outer
+    // quadrature.  We use a parsed function, for its ease of
+    // definition, and the quadrature formula
+    Functions::ParsedFunction<dim> wind;
+    SmartPointer<Quadrature<dim-1> > quadrature_pointer;
+    unsigned int n_cycles;
 };
 
 template <int dim>
 BEMProblem<dim>::BEMProblem() :
     fe(0),
     fev(FE_DGP<dim-1,dim>(0), dim),
+    fe_q(FE_Q<dim-1,dim>(2)),
     dh(tria),
-    dhv(tria)
+    dhv(tria),
+    dhq(tria),
+    wind(dim)
 {}
 
-
 template <int dim>
 BEMProblem<dim>::~BEMProblem() {
     LAPACKFullMatrix<double> * p = system_matrix;
@@ -162,10 +193,49 @@ BEMProblem<dim>::~BEMProblem() {
 }
 
 
+template <int dim> 
+void BEMProblem<dim>::read_parameters(std::string filename) {
+    ParameterHandler prm;
+    
+    prm.declare_entry("Number of cycles", "4", Patterns::Integer());
+    
+    prm.enter_subsection("Outer quadrature rule");
+    prm.declare_entry("Quadrature type", "midpoint", 
+                     Patterns::Selection(QuadratureSelector<(dim-1)>::get_quadrature_names()));
+    prm.declare_entry("Quadrature order", "0", Patterns::Integer());
+    prm.leave_subsection();
+    
+    prm.enter_subsection("Wind function");
+    Functions::ParsedFunction<dim>::declare_parameters(prm, dim);
+    prm.leave_subsection();
+    
+    prm.read_input(filename);
+    
+    n_cycles = prm.get_integer("Number of cycles");                  
+                     
+    prm.enter_subsection("Outer quadrature rule");
+    static QuadratureSelector<dim-1> quadrature
+                     (prm.get("Quadrature type"),
+                      prm.get_integer("Quadrature order"));
+    prm.leave_subsection();
+    
+    prm.enter_subsection("Wind function");
+    wind.parse_parameters(prm);
+    prm.leave_subsection();
+
+    quadrature_pointer = &quadrature;
+}
+
 template <>
-LaplaceKernelIntegration<2>::LaplaceKernelIntegration()
+LaplaceKernelIntegration<3>::LaplaceKernelIntegration(FiniteElement<2,3> &fe)
 {
-    static FE_DGP<2,3> fe(0);
+    // In order to perform the two dimensional singular integration on
+    // the given cell, we use standard formulas derived by Morino and
+    // Chu, as explained in the introduction. In order to do so, we
+    // generate a custom quadrature point with the four vertices and
+    // the middle point. We won't use the weights, and we set them to
+    // 1.
+
     vector<Point<2> > qps(5);
     qps[0] = Point<2>(0,0);
     qps[1] = Point<2>(0,1);
@@ -182,50 +252,24 @@ LaplaceKernelIntegration<2>::LaplaceKernelIntegration()
 
 template <int dim>
 LaplaceKernelIntegration<dim>::~LaplaceKernelIntegration() {
-    FEValues<dim,dim+1> *fp = fe_values;
+    // We delete the pointer. Since this was created via the new
+    // operator, we need to destroy it using delete. But delete does
+    // not take smart pointers, which implies we need to first remove
+    // detach the smart pointer from the fe_values object, and then
+    // delete it.
+    FEValues<dim-1,dim> *fp = fe_values;
     fe_values = 0;
     delete fp;
 }
 
 
 template <>
-void
-LaplaceKernelIntegration<2>::compute_SD_integral_on_cell(vector<double> &dst,
-                                                        DoFHandler<2,3>::active_cell_iterator &cell,
-                                                        const Point<3> &point)
-{
-    Assert(dst.size() == 2,
-          ExcDimensionMismatch(dst.size(), 2));
-    fe_values->reinit(cell);
-    vector<Tensor<2,3> > jacobians = fe_values->get_jacobians();
-    vector<Point<3> > quad_points = fe_values->get_quadrature_points();
-    Point<3> r,a1,a2,n,r_c,n_c;
-    r_c = point-cell->center();
-    n_c = jacobians[4][2];
-    double rn_c = r_c*n_c;
-    vector<double> i_S(4);
-    vector<double> i_D(4);
-    for (unsigned int q_point=0; q_point < 4; ++q_point)
-    {
-       r = point-quad_points[q_point];
-       a1 = jacobians[q_point][0];
-       a2 = jacobians[q_point][1];
-       n =  jacobians[q_point][2];
-       i_S[q_point]=term_S(r,a1,a2,n,rn_c);
-       i_D[q_point]=term_D(r,a1,a2);
-    }
-    dst[0] = (i_S[3]-i_S[1]-i_S[2]+i_S[0]);
-    dst[1] = (i_D[3]-i_D[1]-i_D[2]+i_D[0]);
-
-}
-
-template <int dim>
 double
-LaplaceKernelIntegration<dim>::term_S (const Point<3> &r,
-                                      const Point<3> &a1,
-                                      const Point<3> &a2,
-                                      const Point<3> &n,
-                                      const double &rn_c)
+LaplaceKernelIntegration<3>::term_S (const Point<3> &r,
+                                    const Point<3> &a1,
+                                    const Point<3> &a2,
+                                    const Point<3> &n,
+                                    const double &rn_c)
 {
     Point<3> ra1, ra2, a12;
 
@@ -245,9 +289,9 @@ LaplaceKernelIntegration<dim>::term_S (const Point<3> &r,
 
 }
 
-template <int dim>
+template <>
 double
-LaplaceKernelIntegration<dim>::term_D (const Point<3> &r,
+LaplaceKernelIntegration<3>::term_D (const Point<3> &r,
                                       const Point<3> &a1,
                                       const Point<3> &a2)    
 {
@@ -264,6 +308,55 @@ LaplaceKernelIntegration<dim>::term_D (const Point<3> &r,
 
 }
 
+template <>
+void
+LaplaceKernelIntegration<3>::compute_SD_integral_on_cell(vector<vector<vector<double> > > &dstvv,
+                                                        DoFHandler<2,3>::active_cell_iterator &cell,
+                                                        const vector<Point<3> > &q_points)
+{
+    fe_values->reinit(cell);
+    vector<Tensor<2,3> > jacobians = fe_values->get_jacobians();
+    vector<Point<3> > quad_points = fe_values->get_quadrature_points();
+
+    Point<3> r,a1,a2,n,r_c,n_c;
+    
+    Assert(dstvv.size() == fe_values->dofs_per_cell,
+          ExcDimensionMismatch(dstvv.size(), fe_values->dofs_per_cell));
+          
+    for(unsigned int i=0; i<fe_values->dofs_per_cell; ++i) {
+       vector<vector<double> > & dstv = dstvv[i];
+       Assert(dstv.size() == q_points.size(),
+              ExcDimensionMismatch(dstv.size(), q_points.size()));
+    
+       /* Check only the first size. */
+       Assert(dstv[0].size() == 2,
+              ExcDimensionMismatch(dstv[0].size(), 2));
+       
+    
+       n_c = jacobians[4][2];
+       
+       for(unsigned int outer_q=0; outer_q<q_points.size(); ++outer_q) {
+           const Point<3> &point = q_points[outer_q];
+           vector<double> &dst = dstv[outer_q];
+           r_c = point-cell->center();
+           double rn_c = r_c*n_c;
+           vector<double> i_S(4);
+           vector<double> i_D(4);
+           for (unsigned int inner_q_point=0; inner_q_point < 4; ++inner_q_point)
+           {
+               r = point-quad_points[inner_q_point];
+               a1 = jacobians[inner_q_point][0];
+               a2 = jacobians[inner_q_point][1];
+               n =  jacobians[inner_q_point][2];
+               i_S[inner_q_point]= term_S(r,a1,a2,n,rn_c) * fe_values->shape_value(i,inner_q_point);
+               i_D[inner_q_point]= term_D(r,a1,a2) * fe_values->shape_value(i,inner_q_point);
+           }
+           dst[0] = (i_S[3]-i_S[1]-i_S[2]+i_S[0]);
+           dst[1] = (i_D[3]-i_D[1]-i_D[2]+i_D[0]);
+       }
+    }
+}
+
 template <int dim>
 void BEMProblem<dim>::read_domain() {
     // Center of the ball. It is the origin by default.
@@ -307,75 +400,132 @@ void BEMProblem<dim>::refine_and_resize() {
     
     system_matrix = new LAPACKFullMatrix<double>(ndofs, ndofs);
 
-    Vn.reinit(ndofs);
-    rhs.reinit(ndofs);
+    system_rhs.reinit(ndofs);
     phi.reinit(ndofs);
-    vs.reinit(nvdofs);
 }    
 
 template <int dim>
 void BEMProblem<dim>::assemble_system() {
     
-    Point<dim> wind;
-    wind[0] = 1.;
-    
     typename DoFHandler<dim-1,dim>::active_cell_iterator
        celli = dh.begin_active(),
        cellj = dh.begin_active(),
-       cellv = dhv.begin_active(),
+       cellvi = dhv.begin_active(),
        endc = dh.end();
     
-    QMidpoint<dim-1> midpoint;
-    FEValues<dim-1,dim> fe_mid(fe, midpoint,
-                              update_values |
-                              update_cell_normal_vectors |
-                              update_quadrature_points);
+    // Outer quadrature rule. If we choose midpoint quadrature rule,
+    // then this is a collocation method. If we choose any other
+    // Quadrature rule, then this is Galerkin method.
+    Quadrature<dim-1> &quadrature_outer = *quadrature_pointer;
+    QMidpoint<dim-1> quadrature_mid;
+
+    
+    FEValues<dim-1,dim> fe_outer(fe, quadrature_outer,
+                                update_values |
+                                update_cell_normal_vectors |
+                                update_quadrature_points);
+    
+    FEValues<dim-1,dim> fe_inner(fe, quadrature_mid,
+                                update_values |
+                                update_cell_normal_vectors |
+                                update_quadrature_points);
+    
+    const unsigned int n_q_points_outer = fe_outer.n_quadrature_points;
+    const unsigned int n_q_points_inner = fe_inner.n_quadrature_points;
+    
+    vector<unsigned int> dofs_i(fe.dofs_per_cell);
+    vector<unsigned int> dofs_j(fe.dofs_per_cell);
+    vector<unsigned int> dofs_v_i(fev.dofs_per_cell);
     
-    vector<unsigned int> dofsi(fe_mid.n_quadrature_points);
-    vector<unsigned int> dofsj(fe_mid.n_quadrature_points);
-    vector<unsigned int> dofsv(dim);
+    vector<vector<vector<double> > > single_double_layer_potentials
+       (fe.dofs_per_cell, vector<vector<double> >
+        (n_q_points_outer, vector<double> (2, 0.) ) ); 
+
+    vector<Vector<double> > cell_wind(n_q_points_inner, Vector<double>(dim) );
+    vector<double> normal_wind(n_q_points_inner);
     
-    // Temporary matrix 
-    FullMatrix<double> _B(dh.n_dofs(), dh.n_dofs());
+    Vector<double>     local_rhs(fe.dofs_per_cell);
+    FullMatrix<double>  local_matrix(fe.dofs_per_cell, fe.dofs_per_cell);
     
     // The kernel.
-    LaplaceKernelIntegration<dim-1> kernel;
-    
-    // i runs on points, j runs on cells.
-    for(; celli != endc; ++celli, ++cellv) {
-           fe_mid.reinit(celli);
-           Point<dim> ci = celli->center();
-           Point<dim> ni = fe_mid.cell_normal_vector(0);
-
-           celli->get_dof_indices(dofsi);
-           cellv->get_dof_indices(dofsv);
-
-           // Vn vector:
-           Vn(dofsi[0]) = wind*ni;
-           vs(dofsv[0]) = ni[0];
-           vs(dofsv[1]) = ni[1];
-           vs(dofsv[2]) = ni[2];
+    LaplaceKernelIntegration<dim> kernel(fe);
+    
+    // i runs on outer integration, j runs on inner integration.
+    for(; celli != endc; ++celli, ++cellvi) {
+       fe_outer.reinit(celli);
+       
+       const vector<Point<dim> > &q_points_outer = fe_outer.get_quadrature_points();
+       const vector<Point<dim> > &normals_outer = fe_outer.get_cell_normal_vectors();
+
+       celli->get_dof_indices(dofs_i);
+       cellvi->get_dof_indices(dofs_v_i);
         
-           // Now the two matrices.
-           for(cellj = dh.begin_active(); cellj != endc; ++cellj) {
-               vector<double> SD(2,0.); // Single and Double layers.
-               cellj->get_dof_indices(dofsj);
-               kernel.compute_SD_integral_on_cell(SD,
-                                                  cellj,
-                                                  ci);
-               _B (dofsi[0], dofsj[0]) += -SD[0];
-               if(dofsi[0] != dofsj[0])
-                   (*system_matrix)(dofsi[0], dofsj[0]) += -SD[1];
-               if(dofsi[0] == dofsj[0])
-                   (*system_matrix)(dofsi[0], dofsj[0]) += 1.;
+       // Now the mass matrix and the single and double layer
+       // potentials. Notice that instead of integrating the single
+       // layer potential against the normal velocity, we integrate
+       // it against the average value of the velocity in the given
+       // cell.
+       //
+       // The reason for proceeding in this way is that in the
+       // current three-dimensional formulation we can only integrate
+       // the constants against the single layer potential. While
+       // this is certainly a rough approximation, it suffices for
+       // the purpose of this example.
+       for(cellj = dh.begin_active(); cellj != endc; ++cellj) {
+           local_rhs = 0;
+           local_matrix = 0;
+           
+           fe_inner.reinit(cellj);
+           cellj->get_dof_indices(dofs_j);
+           
+           const vector<Point<dim> > &q_points_inner = fe_inner.get_quadrature_points();
+           const vector<Point<dim> > &normals_inner = fe_inner.get_cell_normal_vectors();
+           wind.vector_value_list(q_points_inner, cell_wind);
+           
+           
+           for(unsigned int q=0; q<n_q_points_inner; ++q) {
+               normal_wind[q] = 0;
+               for(unsigned int d=0; d<dim; ++d)
+                   normal_wind[q] += normals_inner[q][d] * cell_wind[q](d);
+           }
+           
+           kernel.compute_SD_integral_on_cell(single_double_layer_potentials, 
+                                              cellj, q_points_outer);
+           
+           for(unsigned int i=0; i<fe.dofs_per_cell; ++i) {
+               for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
+                   for(unsigned int qo=0; qo<n_q_points_outer; ++qo) {
+                       local_rhs(i) += ( - single_double_layer_potentials[j][qo][0]  * 
+                                         normal_wind[qo]                       *
+                                         fe_outer.shape_value(i,qo)           *
+                                         fe_outer.JxW(qo) );
+                       
+                       if(celli->index() != cellj->index())
+                           local_matrix(i,j) += ( -single_double_layer_potentials[j][qo][1] * 
+                                                  fe_outer.shape_value(i,qo)                *
+                                                  fe_outer.JxW(qo) );
+                       // When the indices are the same, we assemble
+                       // also the mass matrix.
+                       if(celli->index() == cellj->index())
+                           local_matrix(i,j) += ( fe_outer.shape_value(i,qo)        * 
+                                                  fe_outer.shape_value(j,qo)        *
+                                                  fe_outer.JxW(qo) );
+                   }
+               }
+           }
+           
+           for(unsigned int i=0; i<fe.dofs_per_cell; ++i) {
+               system_rhs(dofs_i[i]) += local_rhs(i);
+               for(unsigned int j=0; j<fe.dofs_per_cell; ++j) 
+                   (*system_matrix)(dofs_i[i],dofs_j[j]) += local_matrix(i,j);
            }
        }
-    _B.vmult(rhs, Vn);
+    }
 }
 
 template <int dim>
 void BEMProblem<dim>::solve_system() {
-    phi.swap(rhs);
+    phi.swap(system_rhs);
     system_matrix->compute_lu_factorization();
     system_matrix->apply_lu_factorization(phi, false);
 }
@@ -399,11 +549,11 @@ void BEMProblem<dim>::output_results(unsigned int cycle) {
 
 template <int dim>
 void BEMProblem<dim>::run() {
-    read_domain();
     
-    const unsigned int number_of_cycles = 4;
+    read_parameters("parameters.prm");
+    read_domain();
     
-    for(unsigned int cycle=0; cycle<number_of_cycles; ++cycle) {
+    for(unsigned int cycle=0; cycle<n_cycles; ++cycle) {
        refine_and_resize();
        assemble_system();
        solve_system();

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.