]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add a few things.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 15 Mar 2009 07:03:27 +0000 (07:03 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 15 Mar 2009 07:03:27 +0000 (07:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@18485 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 366fe42b2e2b8fd8ec39685f387f1ce2dcafbb18..98881f9775a229507d54fb8225ec2c5c07a10ff2 100644 (file)
@@ -1,63 +1,91 @@
+<br>
+
+<i>This program was contributed by Luca Heltai. 
+</i>
+
 <a name="Intro"></a>
 
 <h1>Introduction</h1>
 
 <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:
+The incompressible 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}
+\f{align*}
   \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 } \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]
+  \nabla \cdot \mathbf{v}&=0  
+  &\text{in } \mathbb{R}^n\backslash\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. 
+region representing the body around which the fluid moves. $\rho$ is
+independent of the pressure: due to the second equation the flow field is
+incompressible and so the density can not depend on the pressure. 
+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, inertial forces and the external
+forces. This is the opposite case of the Stokes equations
+discussed in @ref step_22 "step-22" which are the limit case of dominant
+viscosity, i.e. where the velocity is so small that inertia forces can be
+neglected. On the other hand, owing to the assumed incompressibility, the
+equations are not suited for very high speed gas flows where compressibility
+and the equation of state of the gas have to be taken into account, leading
+to the Euler equations of gas dynamics, a hyperbolic system.
+
+For the purpose of this tutorial program, we will consider only stationary
+flow without external forces (i.e. $\mathbf{g}=0$, though we will consider 
+forces due to the obstacle $\Omega$):
+\f{align*}
+  (\mathbf{v}\cdot\nabla)\mathbf{v}
+  &=
+  -\frac{1}{\rho}\nabla p
+  \qquad &\text{in } \mathbb{R}^n \backslash \Omega
+  \\
+  \nabla \cdot \mathbf{v}&=0  
+  &\text{in } \mathbb{R}^n\backslash\Omega
+\f}
 
-Uniqueness of the solution is ensured by adding the boundary conditions
+
+Uniqueness of the solution of the Euler equations 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,
+    \mathbf{n}\cdot\mathbf{v}& = 0 \qquad && \text{ on } \partial\Omega \\
+    \mathbf{v}& = \mathbf{v}_\infty && \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. 
+which is to say that the body moves 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. 
 
 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:
+The first step in this process is to change the frame of reference from a
+coordinate system moving along with the body to one in which the body moves
+through a fluid that is at rest at infinity. This can be expressed by
+introducing a new velocity $\mathbf{\tilde{v}}=\mathbf{v}-\mathbf{v}_\infty$ for
+which we find that the same equations hold (because $\nabla\cdot
+\mathbf{v}_\infty=0$) and we have boundary conditions
 \f[
-\nabla\cdot\mathbf{v} = 
-\nabla\cdot(\mathbf{v_\infty} + \mathbf{v_p}) = 
-\nabla\cdot\mathbf{v_p} = 0.
+  \label{eq:boundary-conditions-tilde}
+  \begin{aligned}
+    \mathbf{n}\cdot\mathbf{\tilde{v}}& = -\mathbf{v}_\infty \qquad && \text{ on } \partial\Omega \\
+    \mathbf{\tilde{v}}& = 0 && \text{ when } |\mathbf{x}| \to \infty,
+  \end{aligned}
 \f]
 
 If we assume that the fluid is irrotational, i.e., $\nabla \times
@@ -65,13 +93,17 @@ If we assume that the fluid is irrotational, i.e., $\nabla \times
 velocity, and consequently also the perturbation velocity, as the
 gradient of a scalar function:
 \f[
-  \mathbf{v_p}=\nabla\phi,
+  \mathbf{\tilde{v}}=\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]
+\f{align*}
+\label{laplace}
+\Delta\phi &= 0 \qquad &&\text{in}\ \mathbb{R}^n\backslash\Omega,
+          \\
+\mathbf{n}\cdot\nabla\phi &= -\mathbf{n}\cdot\mathbf{v}_\infty 
+&& \text{on}\ \partial\Omega
+\f}
 
 We will now reformulate this equation in integral form using the
  Green identity:
@@ -83,10 +115,12 @@ We will now reformulate this equation in integral form using the
   (-\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:
+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}
@@ -97,55 +131,56 @@ called fundamental solutions of Laplace equation:
     -\frac{1}{4\pi}\frac{1}{|\mathbf{x}-\mathbf{y}|}&& \text{for } n=3,  
 \end{aligned}
 \f]
-satisfy in variational sense
+satisfy in a variational sense the equations
 \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
+If we substitute $u$ and $v$ in the Green identity with the solution
+$\phi$ and with the fundamental solution of the 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,
+  \qquad \forall\mathbf{x}\in \mathbb{R}^n\backslash\Omega.
 \f]
-that we can write more compactly using the Single and Double
+We can write this more compactly using the so-called 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.
+  \qquad \forall\mathbf{x}\in \mathbb{R}^n\backslash\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$.
+Consequently,
+\f[
+  \phi(\mathbf{x}) = -\left(S[\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
+  + 
+  (D\phi)(\mathbf{x})
+  \qquad \forall\mathbf{x}\in \mathbb{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
+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$:
 
 \f[\label{SD}
-  \frac{1}{2}\phi(\mathbf{x}) = \left(S \frac{\partial \phi}{\partial \mathbf{n}_y}\right)(\mathbf{x})
+  \frac{1}{2}\phi(\mathbf{x}) = - \left(S [\mathbf{n}\cdot\mathbf{v}_\infty]\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:
+which is the integral formulation we were looking for. Substituting the single
+  and double layer operators we get:
 \f[               
   \pi\phi(\mathbf{x})=
   \int_{\partial \Omega}  \ln|\mathbf{x}-\mathbf{y}| \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
@@ -159,7 +194,109 @@ for two dimensional flows and
   \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.
+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$.
+
+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
+for the velocities. To this end, note that we have just computed
+$\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}) = -\left(S[\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
+  + 
+  (D\phi)(\mathbf{x})
+  \qquad \forall\mathbf{x}\in \mathbb{R}^n\backslash\Omega.
+\f]
+where now we have everything that is on the right hand side ($S$ and $D$ are
+integrals we can evaluate, the normal velocity on the boundary is given, and
+$\phi$ on the boundary we have just computed). Finally, we can then recover
+the velocity as $\mathbf{\tilde v}=\nabla \phi$. 
+
+As a final test, let us verify that this velocity indeed satisfies the
+momentum balance equation for a stationary flow field, i.e. whether
+$\mathbf{v}\cdot\nabla\mathbf{v} = -\frac 1\rho \nabla p$ where
+$\mathbf{v}=\mathbf{\tilde
+v}+\mathbf{v}_\infty=\nabla\phi+\mathbf{v}_\infty$ for some (unknown) pressure
+$p$ and given $\rho,\mathbf{g}$. The left hand side of this
+equation equates to
+@f{align*}
+  \mathbf{v}\cdot\nabla\mathbf{v}
+  &=
+  [(\nabla\phi+\mathbf{v}_\infty)\cdot\nabla] (\nabla\phi+\mathbf{v}_\infty)
+  \\
+  &=
+  [(\nabla\phi+\mathbf{v}_\infty)\cdot\nabla] (\nabla\phi)
+@f}
+where we have used that $\mathbf{v}_\infty$ is constant. We would like to
+write this expression as the gradient of something (remember that $\rho$ is a
+constant). The next step is more
+convenient if we consider the components of the equation individually
+(summation over indices that appear twice is implied):
+@f{align*}
+  [\mathbf{v}\cdot\nabla\mathbf{v}]_i
+  &=
+  (\partial_j\phi+v_{\infty,j}) \partial_j \partial_i\phi
+  \\
+  &=
+  \partial_j [(\partial_j\phi+v_{\infty,j}) \partial_i\phi]
+  -
+  \partial_j [(\partial_j\phi+v_{\infty,j})] \partial_i\phi  
+  \\
+  &=
+  \partial_j [(\partial_j\phi+v_{\infty,j}) \partial_i\phi]
+@f}
+because $\partial_j \partial_j\phi = \Delta \phi = 0$ and $\textrm{div}
+\ \mathbf{v}_\infty=0$. Next,
+@f{align*}
+  [\mathbf{v}\cdot\nabla\mathbf{v}]_i
+  &=
+  \partial_j [(\partial_j\phi+v_{\infty,j}) \partial_i\phi]
+  \\
+  &=
+  \partial_j [(\partial_j\phi) (\partial_i\phi)]
+  +
+  \partial_j [v_{\infty,j} \partial_i\phi]
+  \\
+  &=
+  \partial_j [(\partial_j\phi) (\partial_i\phi)]
+  +
+  \partial_j [v_{\infty,j}] \partial_i\phi
+  +
+  v_{\infty,j} \partial_j \partial_i\phi
+  \\
+  &=
+  \partial_j [(\partial_j\phi) (\partial_i\phi)]
+  +
+  v_{\infty,j} \partial_j \partial_i\phi
+  \\
+  &=
+  \partial_i \partial_j [(\partial_j\phi) \phi]
+  -
+  \partial_j [\partial_i (\partial_j\phi) \phi]
+  +
+  \partial_i [v_{\infty,j} \partial_j \phi]
+  -
+  \partial_i [v_{\infty,j}] \partial_j \phi
+@f}
+Again, the last term disappears because $\mathbf{v}_\infty$ is constant and we
+can merge the first and third term into one:
+@f{align*}
+  [\mathbf{v}\cdot\nabla\mathbf{v}]_i
+  &=
+  \partial_i (\partial_j [(\partial_j\phi) \phi + v_{\infty,j} \partial_j \phi])
+  -
+  \partial_j [\partial_i (\partial_j\phi) \phi]
+  \\
+  &=
+  \partial_i [(\partial_j\phi)(\partial_j \phi) + v_{\infty,j} \partial_j \phi]
+  -
+  \partial_j [\partial_i (\partial_j\phi) \phi]
+@f}
+
+
 
 <h3>The numerical approximation</h3>
 
@@ -168,44 +305,38 @@ 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,
+formulation of the 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:
+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
+To this end, let $\mathcal{T}_h = \bigcup_i 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 <i>element</i> or <i>cell</i>, independently on the
-dimension $n$ of the surrounding space $R^n$.
-
+quadrilateral an <i>element</i> or <i>cell</i>, independently of the
+dimension $n$ of the surrounding space $\mathbb{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\},
+  \forall i\},
 \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
+i.e., one if $\mathbf{x}$ belongs to $K_i$, and zero otherwise.  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,
+  \boldsymbol{\alpha} := \{ \phi_i \},
 \f]
-where summation  is implied over repeated indexes.
+where summation  is implied over repeated indexes. Note that we use
+discontinuous elements here &mdash; in fact, there is no real reason to use
+continuous ones in the first place since the integral formulation does not
+imply any derivatives on our trial functions so continuity is unnecessary.
+
 
 <h3> Galerkin boundary element method </h3>
 
@@ -215,20 +346,17 @@ 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[
+\f{align*}
   \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}
+    \int_{\Gamma_x} \phi_h(\mathbf{x}) \eta(\mathbf{x})\,ds_x =
+    & 2\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 +
+    \eta(\mathbf{x}) \,ds_x\,ds_y 
     \\
-    & \frac{1}{\pi} \int_{\Gamma_x}\int_{\Gamma_y} \frac{ \partial
+    & + 2\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]
+\f}
 
 The linearity of the integral operator makes this problem equivalent 
 to solving the linear system
@@ -240,11 +368,11 @@ where
 \f[
 \begin{aligned}
 \mathbf{M}_{ij}&= |K_i|\delta_{ij}\\
-\mathbf{A}_{ij}&= \frac{1}{\pi} \int_{K_i}\int_{K_j}
+\mathbf{A}_{ij}&= 2\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}}
+\mathbf{b}_i&= 2\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}
@@ -254,11 +382,11 @@ 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
+As usual in these 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
+boundary element $\hat K := [0,1]^{n-1}$, and we perform the integrations after a
 change of variables from the real element $K_i$ to the reference
 element $\hat K$.
 
@@ -268,7 +396,8 @@ 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
+orthogonal (on our polygonal approximation of the boundary of $\Omega$), 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}
@@ -304,7 +433,10 @@ we end up with:
   \ln|\alpha| \, \mathbf{n}_y\cdot\mathbf{v_\infty}
    \, ,
 \f]
-which can be computed analytically with Gauss-log quadrature formulae.
+which can be computed analytically with Gauss-log quadrature formulae
+(i.e. Gauss formula derived for integrals with a logarithmic weight; deal.II
+conveniently has the QGaussLog class for this).
+
 
 <h3>Singular integrals in three dimensions</h3>
 
@@ -326,7 +458,7 @@ We are interested in calculating the integrals
     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
+where $\mathbf{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)$.
@@ -340,7 +472,7 @@ Introducing the tangent vectors to the surface of the panels,
 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{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.
@@ -374,7 +506,7 @@ take the form above, where the terms $I_S$ and $I_D$ are given by
       \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)}%
+    \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]
@@ -382,7 +514,7 @@ 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)}%
+  \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]
 

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.