+<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
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:
(-\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}
-\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
\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>
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 — 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>
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
\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}
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$.
$\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}
\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>
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)$.
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.
\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]
\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]
projects / dealii-svn.git / commitdiff