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. $\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.
+region representing 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
+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$):
+flow without external forces:
\f{align*}
(\mathbf{v}\cdot\nabla)\mathbf{v}
&=
\f[
\label{eq:boundary-conditions}
\begin{aligned}
- \mathbf{n}\cdot\mathbf{v}& = 0 \qquad && \text{ on } \partial\Omega \\
+ -\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 is to say that the body moves is at rest in our coordinate systems and
+
+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.
+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.
For both stationary and non stationary flow, the solution process
starts by solving for the velocity in the second equation and
\f[
\label{eq:boundary-conditions-tilde}
\begin{aligned}
- \mathbf{n}\cdot\mathbf{\tilde{v}}& = -\mathbf{v}_\infty \qquad && \text{ on } \partial\Omega \\
+ -\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,
\end{aligned}
\f]
\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
+ -\mathbf{n}\cdot\nabla\phi &= \mathbf{n}\cdot\mathbf{v}_\infty
+ && \text{on}\ \partial\Omega
\f}
while the momentum equation reduces to the Bernoulli's equation
\f[
-\frac{p}{\rho} + \frac{\partial \phi}{\partial t} +g z
-+\frac{1}{2} | \nabla \phi |^2 = 0 \in \Omega,
+\frac{p}{\rho} +\frac{1}{2} | \nabla \phi |^2 = 0 \in \Omega,
\f]
and the pressure and velocity are uncoupled.
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
+ (\Delta u)v\,dx + \int_{\partial\Omega} u\frac{\partial v}{\partial \mathbf{n}} \,ds,
=
\int_{\mathbb{R}^n\backslash\Omega}
- (-\Delta v)u\,dx - \int_{\partial\Omega} u\frac{\partial v}{\partial \mathbf{n}} \,ds,
+ (\Delta v)u\,dx + \int_{\partial\Omega} \frac{\partial u}{\partial \mathbf{n} }v \,ds
\f]
-where $\mathbf{n}$ is the normal to the surface of $\Omega$ pointing
+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
\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{y}-\mathbf{x}) = &
+ -\frac{1}{2\pi}\ln|\mathbf{y}-\mathbf{x}| \qquad && \text{for } n=2
\\
- G(\mathbf{x}-\mathbf{y}) = &
- -\frac{1}{4\pi}\frac{1}{|\mathbf{x}-\mathbf{y}|}&& \text{for } n=3,
+ 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 variational sense the equations:
+satisfy in a variational sense the equation:
\f[
- -\Delta_x G(\mathbf{x}-\mathbf{y}) = \delta(\mathbf{x}-\mathbf{y}),
+ \Delta_y G(\mathbf{y}-\mathbf{x}) = \delta(\mathbf{y}-\mathbf{x}),
\f]
-where the derivation is done in the variable $\mathbf{x}$.
+where the derivation 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, we obtain:
+respectively, as long as $\mathbf{x}$ is chosen in the region
+$\mathbb{R}^n\backslash\Omega$, 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
+ \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.
\f]
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})
+ \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.
\f]
(The name of these operators comes from the fact that they describe the
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$.
+$-\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.
+ \phi(\mathbf{x}) + (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.
\f]
If we take the limit for $\mathbf{x}$ tending to $\partial\Omega$ of
the above equation, using well known properties of the single and double layer
$\Omega$:
\f[\label{SD}
- \alpha(\mathbf{x})\phi(\mathbf{x}) =
- - \left(S [\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
- +
- (D\phi)(\mathbf{x})
+ \alpha(\mathbf{x})\phi(\mathbf{x}) + (D\phi)(\mathbf{x}) =
+ \left(S [\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
\quad \mathbf{x}\in \partial\Omega,
\f]
which is the integral formulation we were looking for, where the
-quantity $\alpha(\mathbf{x}_i)$ is the fraction of solid angle by
-which the point $\mathbf{x}_i$ sees the domain $\Omega$.
+quantity $\alpha(\mathbf{x}_i)$ is the fraction of angle or solid
+angle by which the point $\mathbf{x}_i$ sees the domain $\Omega$.
-Substituting the single
- and double layer operators we get:
+Substituting the single and double layer operators we get:
\f[
- \alpha(\mathbf{x}) \phi(\mathbf{x})=
- \frac{1}{2\pi}\int_{\partial \Omega} \ln|\mathbf{x}-\mathbf{y}| \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
- +
- \frac{1}{2\pi}\int_{\partial \Omega} \frac{ (\mathbf{x}-\mathbf{y})\cdot\mathbf{n}_y }{ |\mathbf{x}-\mathbf{y}|^2 }\,ds_y
+ \alpha(\mathbf{x}) \phi(\mathbf{x})
+ - \frac{1}{2\pi}\int_{\partial \Omega} \frac{ (\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y }{ |\mathbf{y}-\mathbf{x}|^2 }\,ds_y
+ =
+ -\frac{1}{2\pi}\int_{\partial \Omega} \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{1}{|\mathbf{x}-\mathbf{y}|} \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
- +
- \frac{1}{4\pi}\int_{\partial \Omega} \frac{ (\mathbf{x}-\mathbf{y})\cdot\mathbf{n}_y }{ |\mathbf{x}-\mathbf{y}|^3 }\phi(\mathbf{y})\,ds_y
+ \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
+ =
+ \frac{1}{4\pi}\int_{\partial \Omega} \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
$\Omega$ can be defined using the double layer potential itself:
\f[
\alpha(\mathbf{x}) :=
-\int_{\partial \Omega} \frac{ (\mathbf{x}-\mathbf{y})\cdot\mathbf{n}_y }
-{ |\mathbf{x}-\mathbf{y}|^{dim} }\phi(\mathbf{y})\,ds_y =
-\int_{\partial \Omega} \frac{ \partial G(\mathbf{x}-\mathbf{y}) }{\partial \mathbf{n}_y} \, ds_y
+\frac{1}{2(n-1)\pi}\int_{\partial \Omega} \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.
\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})$.
+
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
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})
- +
+ \phi(\mathbf{x}) +
(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.
\f]
where now we have everything that is on the right hand side ($S$ and $D$ are
$\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
+$\mathbb{R}^n\backslash\Omega$ is always zero by definition.
+
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
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
+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
\f[
\label{eq:definition-of-element}
\phi_h(\mathbf{x}) := \phi_i \psi_i(\mathbf{x}), \qquad
- \boldsymbol{\alpha} := \{ \phi_i \},
+ \boldsymbol{\phi} := \{ \phi_i \},
\f]
where summation is implied over repeated indexes. Note that we could use
discontinuous elements here — in fact, there is no real reason to use
<h3> Collocation boundary element method </h3>
-By far, the most common approximation of boundary integral equations,
+By far, the most common approximation of boundary integral equations
is by use of the collocation based boundary element method.
This method requires the evaluation of the boundary integral equation
at a number of collocation points which is equal to the number of
unknowns of the system. The choice of these points is a delicate
matter, that requires a careful study. Assume that these points are
-known for the moment, and call them $\mathbf x_i$ with $i=0...n_dofs$.
+known for the moment, and call them $\mathbf x_i$ with $i=0...n\_dofs$.
The problem then becomes:
Given the datum $\mathbf{v}_\infty$, find a function $\phi_h$ in $V_h$
-such that, the following $n_dofs$ equations are satisfied:
+such that, the following $n\_dofs$ equations are satisfied:
\f{align*}
- \alpha(mathbf{x}_i) \phi_h(\mathbf{x}_i) =
- & \int_{\Gamma_y} G(\mathbf{x}_i-\mathbf{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
- \\
- & + \int_{\Gamma_y} \frac{ \partial
- G(\mathbf{x}_i-\mathbf{y})}{\partial\mathbf{n}_y }
- \phi_h(\mathbf{y}) \,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
problem can be written as the following linear system:
\f[
\label{eq:linear-system}
-(\mathbf{A}-\mathbf{N})\boldsymbol\phi = \mathbf{b},
+(\mathbf{A}+\mathbf{N})\boldsymbol\phi = \mathbf{b},
\f]
where
\f[
\begin{aligned}
-\mathbf{A}_{ii}&= \int_\Gamma
-\frac{\partial G(\mathbf{x}_i-\mathbf{y})}{\partial \mathbf{n}_y}\,ds_y \\
+\mathbf{A}_{ii}&= -\int_\Gamma
+\frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y = \sum_{j=0}^{n\_dofs} N_{ij}\\
\mathbf{N}_{ij}&= \int_\Gamma
- \frac{\partial G(\mathbf{x}_i-\mathbf{y})}{\partial \mathbf{n}_y}
+ \frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}
\psi_j(\mathbf{y}) \,ds_y
\\
\mathbf{b}_i&= \int_\Gamma
- G(\mathbf{x}_i-\mathbf{y}) \, \mathbf{n}_y\cdot\mathbf{v_\infty}
+ G(\mathbf{y}-\mathbf{x}_i) \, \mathbf{n}_y\cdot\mathbf{v_\infty}
\psi_i(\mathbf{y}) ds_y.
\end{aligned}
\f]
-The computation of the entries of the matrices $\mathbf{A}$,
-$\mathbf{N}$ 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$.
+From a linear algebra point of view, the best possible choice of the
+collocation points is the one that renders the matrix
+$\mathbf{A}+\mathbf{N}$ the most diagonally dominant. A natural choice
+is then to select the $\mathbf{x}_i$ collocation points to be the
+support points of the nodal basis functions $\psi_i(\mathbf{x})$.
+
+With this choice of collocation points, the computation of the entries
+of the matrices $\mathbf{A}$, $\mathbf{N}$ and of the right hand side
+$\mathbf{b}$ requires the evaluation of singular integrals on the
+elements $K_i$ of the triangulation $\mathcal{T}_h$.
As usual in these cases, all integrations are performed on a reference
simple domain, i.e., we assume that each element $K_i$ of
change of variables from the real element $K_i$ to the reference
element $\hat K$.
-<h3> Singular integrals in two dimension. </h3>
+<h3> Treating the singular integrals. </h3>
In two dimensions it is not necessary to compute the diagonal elements
$\mathbf{N}_{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
+zero because $\mathbf{n}_y$ and $(\mathbf{y}-\mathbf{x})$ are
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}
- \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
-(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>
-
-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 $\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)$.
-
-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}
+ \int_{K_i}
+ \ln|\mathbf{y}-\mathbf{x}_i| \, \mathbf{n}_y\cdot\mathbf{v_\infty} \,ds_y.
\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.
-
+This can be easily treated by the QGaussLogR quadrature
+formula.
-<h3>What the program does</h3>
+Similarly, it is possible to use the QGaussOneOverR quadrature formula
+to perform the singular integrations in three dimensions. The
+interested reader will find detailed explanations on how these
+quadrature rules work in their documentation.
+The resulting matrix $\mathbf{A}+\mathbf{N}$ is full. Depending on its
+size, it might be convenient to use a direct solver or an iterative
+one. For the purpose of this example code, we chose to use only a
+direct solver, which limits the applicability of this method to
+relatively small problems. Remember however that it is possible to
+obtain very accurate results with relatively few surface elements.
// (e.g., Stokes or Maxwell).
//
// The most noticeable difference is the fact that the final
- // matrix is full, and that we have two nested loops on cells
- // instead of the usual one we have in finite element method.
- //
- // The reason for this is that while the basis functions have a
- // compact support, their convolution with the fundamental
- // solution of the laplace equation is global, and needs to be
- // integrated against all other basis functions.
+ // matrix is full, and that we have a nested loop inside the usual
+ // loop on cells that visits all support points of the degrees of
+ // freedom. Moreover, when the support point lyes inside the cell
+ // which we are visiting, then the integral we perform becomes
+ // singular.
//
// The practical consequence is that we have two sets of
// quadrature formulas, finite element values and temporary
- // elements, one for the inner integration and one for the outer
- // integration. We allow for different quadrature rules to be used
- // in the two integrations to preserve generality and to allow,
- // for example, the use of collocation method (by specifying midpoint
- // quadrature formula on the outer integration).
+ // elements, one for standard integration and one for the singular
+ // integration, which are used where necessary.
void assemble_system();
- // The only difference in the solution of the system is that the
- // matrix is a LAPACKFullMatrix, which requires a different
- // treatment with respect to what we saw in most of the other
- // examples. Besides from this detail, things proceeds pretty much
- // in the same way as usual.
+ // Notwithstanding the fact that the matrix is full, we use a
+ // SparseMatrix object and the SparseDirectUMFPACK solver, since
+ // in our experience it works better then using, for example, the
+ // LapackFullMatrix class. An alternative approach would be the
+ // use of 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.
+ //
+ // We should notice moreover that the solution we will obtain will
+ // only be unique up to an additive constant. This is taken care
+ // of in the solve_system method, which filters out the mean value
+ // of the solution at the end of the computation.
void solve_system();
// Once we obtained a solution on the codimension one domain, we
// argument, and produces the usual finite element classes that we
// saw in all previous examples.
- Triangulation<dim-1, dim> tria;
- FE_Q<dim-1,dim> fe;
- DoFHandler<dim-1,dim> dh;
+ Triangulation<dim-1, dim> tria;
+ FE_Q<dim-1,dim> fe;
+ DoFHandler<dim-1,dim> dh;
// In BEM methods, the matrix that is generated is
// dense. Depending on the size of the problem, the final system
// might be solved by direct LU decomposition, or by iterative
- // methods. Just for the purpose of illustrating the use of the
- // LAPACK classes, we opt for LU decomposition of the final
- // system. Note that this will be very inefficient when the number
- // of dofs grows, since it is of order $n^3$.
+ // 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 the LapackFullMatrix object.
- SparsityPattern sparsity;
- SparseMatrix<double> system_matrix;
- Vector<double> system_rhs;
- Vector<double> phi;
+ SparsityPattern sparsity;
+ SparseMatrix<double> system_matrix;
+ Vector<double> system_rhs;
+ Vector<double> phi;
+ Vector<double> alpha;
// The reconstruction of the solution in the entire space is done
// on a continuous finite element grid of dimension dim. These are
// the usual ones, and we don't comment any further on them.
- Triangulation<dim> external_tria;
- FE_Q<dim> external_fe;
- DoFHandler<dim> external_dh;
- Vector<double> external_phi;
+ Triangulation<dim> external_tria;
+ FE_Q<dim> external_fe;
+ DoFHandler<dim> external_dh;
+ Vector<double> external_phi;
// The following variables are the one that we fill through a
// parameter file. The new objects that we use in this example
// The QuadratureSelector class allows us to generate quadrature
// formulas based on an identifying string and on the possible
// degree of the formula itself. We used this to allow custom
- // selection of the quadrature formulas for the inner integration.
+ // selection of the quadrature formulas for the standard
+ // integration, and to define the order of the singular quadrature
+ // rule.
//
// Notice that the pointer given below for the quadrature rule is
// only used for non singular integrals. Whenever the integral is
// wanted to extend the solution to the entire domain.
Functions::ParsedFunction<dim> wind;
SmartPointer<Quadrature<dim-1> > quadrature_pointer;
+ unsigned int singular_quadrature_order;
unsigned int n_cycles;
unsigned int external_refinement;
+ bool run_in_this_dimension;
bool extend_solution;
};
public:
// The following two functions are the actual calculations of the
// single and double layer potential kernels, that is G and Grad
- // G. They are well defined only if the vector $R = x-y$ is
+ // G. They are well defined only if the vector $R = y-x$ is
// different from zero.
//
// Whenever the integration is performed with the singularity
// used that allows one to integrate arbitrary functions against a
// singular weight on the reference cell.
//
- // In order to do so, it is necessary to provide a
- // "desingularized" single and double layer potentials which can
- // then be integrated on the given cell. When the @p
- // factor_out_singularity parameter is set to true, then the
- // computed kernels do not conatain the singular factor, which is
- // included in the quadrature formulas as a weighting function.
+ // There are two options when the integral is singular. One could
+ // take into account the singularity inside the quadrature formula
+ // as a weigthing function, or one could use a quadrature formula
+ // that is taylored to integrate singular objects, but where the
+ // actual weighting function is one. The use of the first method
+ // requires the user to provide a "desingularized" single and
+ // double layer potentials which can then be integrated on the
+ // given cell. When the @p factor_out_singularity parameter is set
+ // to true, then the computed kernels do not conatain the singular
+ // factor, which is included in the quadrature formulas as a
+ // weighting function. This works best in two dimension, where the
+ // singular integrals are integrals along a segment of a
+ // logarithmic singularity.
+ //
+ // These integrals are somewhat delicate, because inserting a
+ // factor Jx in the variable of integration does not result only
+ // in a factor J appearing as a constant factor on the entire
+ // integral, but also on an additional integral to be added, that
+ // contains the logarithm of J. For this reason in two dimensions
+ // we opt for the desingularized kernel, and use the QGaussLogR
+ // quadrature formula, that takes care of integrating the correct
+ // weight for us.
+ //
+ // In the three dimensional case the singular integral is taken
+ // care of using the QGaussOneOverR quadrature formula. We could
+ // use the desingularized kernel here as well, but this would
+ // require us to be careful about the different scaling of r in
+ // the reference cell and in real space. The quadrature formula
+ // uses as weight 1/r in local coordinates, while we need to
+ // integrate 1/R in real coordinates. A factor of r/R has to be
+ // introduced in the quadrature formula. This can be done
+ // manually, or we simply calculate the standard kernels and then
+ // use a desingularized quadrature formula, i.e., one which is
+ // taylored for singular integrals, but whose weight is 1 instead
+ // of the singularity.
//
// Notice that the QGaussLog quadrature formula is made to
// integrate f(x)ln|x-x0|, but the kernel for two dimensional
// problems has the opposite sign. This is taken care of by
// switching the sign of the two dimensional desingularized
// kernel.
+ //
+ // The last argument to both funcitons is simply ignored in three
+ // dimensions.
double single_layer(const Point<dim> &R,
- bool factor_out_singularity = false);
+ bool factor_out_2d_singularity = false);
Point<dim> double_layer(const Point<dim> &R,
- bool factor_out_singularity = false);
+ bool factor_out_2d_singularity = false);
};
template <int dim>
void BEMProblem<dim>::read_parameters(std::string filename) {
+ deallog << std::endl << "Parsing parameter file " << filename << std::endl
+ << "for a " << dim << " dimensional simulation. " << std::endl;
+
ParameterHandler prm;
prm.declare_entry("Number of cycles", "4", Patterns::Integer());
prm.declare_entry("External refinement", "5", Patterns::Integer());
- prm.declare_entry("Extend solution on the -2,2 box", "false", Patterns::Bool());
+ prm.declare_entry("Extend solution on the -2,2 box", "true", Patterns::Bool());
+ prm.declare_entry("Run 2d simulation", "true", Patterns::Bool());
+ prm.declare_entry("Run 3d simulation", "true", Patterns::Bool());
- prm.enter_subsection("Quadrature rule");
+ prm.enter_subsection("Quadrature rules");
prm.declare_entry("Quadrature type", "gauss",
- Patterns::Selection(QuadratureSelector<(dim-1)>::get_quadrature_names()));
- prm.declare_entry("Quadrature order", "5", Patterns::Integer());
+ Patterns::Selection(QuadratureSelector<(dim-1)>::get_quadrature_names()));
+ prm.declare_entry("Quadrature order", "4", Patterns::Integer());
+ prm.declare_entry("Singular quadrature order", "5", Patterns::Integer());
prm.leave_subsection();
+ // For both two and three dimensions, we set the input data to be
+ // such that the solution is $x+y+c$ or $x+y+z+c$.
prm.enter_subsection("Wind function 2d");
Functions::ParsedFunction<2>::declare_parameters(prm, 2);
+ prm.set("Function expression","1; 1");
prm.leave_subsection();
prm.enter_subsection("Wind function 3d");
Functions::ParsedFunction<3>::declare_parameters(prm, 3);
+ prm.set("Function expression","1; 1; 1");
prm.leave_subsection();
prm.read_input(filename);
-
- n_cycles = prm.get_integer("Number of cycles");
+
+ n_cycles = prm.get_integer("Number of cycles");
external_refinement = prm.get_integer("External refinement");
extend_solution = prm.get_bool("Extend solution on the -2,2 box");
+ run_in_this_dimension = prm.get_bool("Run " +
+ Utilities::int_to_string(dim) +
+ "d simulation");
- prm.enter_subsection("Quadrature rule");
+ prm.enter_subsection("Quadrature rules");
static QuadratureSelector<dim-1> quadrature
- (prm.get("Quadrature type"),
- prm.get_integer("Quadrature order"));
+ (prm.get("Quadrature type"),
+ prm.get_integer("Quadrature order"));
+ singular_quadrature_order = prm.get_integer("Singular quadrature order");
prm.leave_subsection();
prm.enter_subsection(std::string("Wind function ")+
- Utilities::int_to_string(dim)+std::string("d"));
+ Utilities::int_to_string(dim)+std::string("d"));
wind.parse_parameters(prm);
prm.leave_subsection();
template <int dim>
double LaplaceKernel<dim>::single_layer(const Point<dim> &R,
- bool factor_out_singularity) {
- if(factor_out_singularity == true)
- return (dim == 2 ? -1. : 1.)/(2*(dim-1)*numbers::PI);
- else
- if(dim == 2)
- return (-std::log(R.norm()) / (2*numbers::PI) );
- else if(dim == 3)
- return (1./( R.norm()*4*numbers::PI ) );
- else {
- Assert(false, ExcInternalError());
- return 0.;
- }
+ bool factor_out_2d_singularity) {
+ switch(dim) {
+ case 2:
+ if(factor_out_2d_singularity == true)
+ return -1./(2*numbers::PI);
+ else
+ return (-std::log(R.norm()) / (2*numbers::PI) );
+ break;
+ case 3:
+ return (1./( R.norm()*4*numbers::PI ) );
+ break;
+ defualt:
+ Assert(false, ExcInternalError());
+ return 0.;
+ break;
+ }
return 0.;
}
-
+
template <int dim>
Point<dim> LaplaceKernel<dim>::double_layer(const Point<dim> &R,
- bool factor_out_singularity) {
+ bool factor_out_2d_singularity) {
Point<dim> D(R);
switch(dim) {
case 2:
- factor_out_singularity ? D *= 0 : D /= -2*numbers::PI * R.square();
- break;
+ factor_out_2d_singularity ? D *= 0 : D /= -2*numbers::PI * R.square();
+ break;
case 3:
- D /= ( -4*numbers::PI * R.square() *
- ( factor_out_singularity ? 1. : R.norm() ) );
- break;
+ D /= ( -4*numbers::PI * R.square()*R.norm() );
+ break;
default:
- Assert(false, ExcInternalError());
- break;
+ Assert(false, ExcInternalError());
+ break;
}
return D;
}
-
+
template <int dim>
void BEMProblem<dim>::read_domain() {
// formats which are compatible with the boundary element method
// capabilities of deal.II. In particular we can use either UCD or
// GMSH formats. In both cases, we have to be particularly careful
- // with
- // the orientation of the mesh, because, unlike in the standard
- // finite element case, no reordering or compatibility check is
- // performed here.
+ // with the orientation of the mesh, because, unlike in the
+ // standard finite element case, no reordering or compatibility
+ // check is performed here.
//
// All meshes are considered as oriented, because they are
// embedded in a higher dimensional space. See the documentation
GridIn<dim-1, dim> gi;
gi.attach_triangulation (tria);
if(dim == 3) {
- std::ifstream in ("coarse_sphere.inp");
- gi.read_ucd (in);
+ std::ifstream in ("coarse_sphere.inp");
+ gi.read_ucd (in);
} else if(dim == 2) {
- std::ifstream in ("coarse_circle.inp");
- gi.read_ucd (in);
+ std::ifstream in ("coarse_circle.inp");
+ gi.read_ucd (in);
}
tria.set_boundary(1, boundary);
}
const unsigned int ndofs = dh.n_dofs();
deallog << "Levels: " << tria.n_levels()
- << ", potential dofs: " << ndofs << endl;
+ << ", potential dofs: " << ndofs << endl;
// The matrix is a full matrix. Notwithstanding this fact, the
// SparseMatrix class coupled with the SparseDirectUMFPACK solver
system_matrix.clear();
sparsity.reinit(ndofs, ndofs, ndofs);
for(unsigned int i=0; i<ndofs;++i)
- for(unsigned int j=0; j<ndofs; ++j)
- sparsity.add(i,j);
+ for(unsigned int j=0; j<ndofs; ++j)
+ sparsity.add(i,j);
sparsity.compress();
system_matrix.reinit(sparsity);
system_rhs.reinit(ndofs);
phi.reinit(ndofs);
+ alpha.reinit(ndofs);
}
template <int dim>
void BEMProblem<dim>::assemble_system() {
typename DoFHandler<dim-1,dim>::active_cell_iterator
- celli = dh.begin_active(),
- cellj = dh.begin_active(),
- endc = dh.end();
+ cell = dh.begin_active(),
+ endc = dh.end();
// Quadrature formula for the integration of the kernel in non
// singular cells. This quadrature is selected with the parameter
Quadrature<dim-1> &quadrature = *quadrature_pointer;
// We create initially the singular quadratures for the
- // threedimensional problem, since in this case it is only
+ // threedimensional problem, since in this case they only
// dependent on the reference element. This quadrature is a
// standard Gauss quadrature formula reparametrized in such a way
// that allows one to integrate singularities of the kind 1/R
// centered at one of the vertices. Here we define a vector of
- // four such quadratures that will be used later on.
+ // four such quadratures that will be used later on, only in the
+ // three dimensional case.
vector<QGaussOneOverR<2> > sing_quadratures_3d;
- for(unsigned int i=0; i<4; ++i)
- sing_quadratures_3d.push_back(QGaussOneOverR<2>(quadrature.size(), i));
-
-
+ for(unsigned int i=0; i<4; ++i) {
+ sing_quadratures_3d.push_back
+ (QGaussOneOverR<2>(singular_quadrature_order, i, true));
+ }
+
FEValues<dim-1,dim> fe_v(fe, quadrature,
- update_values |
- update_cell_normal_vectors |
- update_quadrature_points |
- update_JxW_values);
+ update_values |
+ update_cell_normal_vectors |
+ update_quadrature_points |
+ update_JxW_values);
const unsigned int n_q_points = fe_v.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(fe.dofs_per_cell);
vector<Vector<double> > cell_wind(n_q_points, Vector<double>(dim) );
double normal_wind;
- Vector<double> local_rhs(fe.dofs_per_cell);
- FullMatrix<double> local_matrix(fe.dofs_per_cell, fe.dofs_per_cell);
+ // Unlike in finite element methods, if we use a collocation
+ // boundary element method, then in each assembly loop we only
+ // assemble the informations that refer to the coupling between
+ // one degree of freedom (the degree associated with support point
+ // i) and the current cell. This is done using a vector of
+ // fe.dofs_per_cell elements, which will then be distributed to
+ // the matrix in the global row i.
+ Vector<double> local_matrix_row_i(fe.dofs_per_cell);
// The kernel.
LaplaceKernel<dim> kernel;
// support of the ith basis function, while j runs on inner
// integration. We perform this check here to ensure that we are
// not trying to use this code for high order elements. It will
- // only work with Q1 elements, that is, for fe_dofs_per_cell =
+ // only work with Q1 elements, that is, for fe.dofs_per_cell ==
// GeometryInfo<dim>::vertices_per_cell.
AssertThrow(fe.dofs_per_cell == GeometryInfo<dim-1>::vertices_per_cell,
- ExcDimensionMismatch(fe.dofs_per_cell, GeometryInfo<dim-1>::vertices_per_cell));
-
- for(; celli != endc; ++celli) {
-
- // On the outer cell, we only need to know how to go from
- // local numbering to global numbering. Each degree of freedom
- // is associated with its support point, which is the ith
- // vertex of the cell.
- celli->get_dof_indices(dofs_i);
-
- for(cellj = dh.begin_active(); cellj != endc; ++cellj) {
-
- // If we are on the same cell, then the integrals we are
- // performing are singular, and they require a special
- // treatment, as explained in the introduction.
- //
- // In all other cases, standard Gauss quadrature rules can
- // be used.
- bool is_singular = (cellj->index() == celli->index());
-
- local_rhs = 0;
- local_matrix = 0;
-
- fe_v.reinit(cellj);
- cellj->get_dof_indices(dofs_j);
-
- const vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
- const vector<Point<dim> > &normals = fe_v.get_cell_normal_vectors();
- wind.vector_value_list(q_points, cell_wind);
-
- if(is_singular == false) {
- for(unsigned int q=0; q<n_q_points; ++q) {
- normal_wind = 0;
- for(unsigned int d=0; d<dim; ++d)
- normal_wind += normals[q][d]*cell_wind[q](d);
-
- for(unsigned int i=0; i<fe.dofs_per_cell; ++i) {
-
- // Distance between the external support point
- // and the quadrature point on the internal
- // cell.
- R = celli->vertex(i)-q_points[q];
-
- local_rhs(i) += ( kernel.single_layer(R) *
- normal_wind *
- fe_v.JxW(q) );
-
- for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
-
- local_matrix(i,j) += ( ( kernel.double_layer(R) *
- normals[q] ) *
- fe_v.shape_value(j,q) *
- fe_v.JxW(q) );
- }
- }
- }
- } else {
- // Now we treat the more delicate case. If we are
- // here, it means that the cell that runs on the j
- // index and the one that runs on the i index are the
- // same. In this case both the single and the double
- // layer potential are singular, and they require a
- // special treatment, as explained in the
- // introduction.
- //
- // In the two dimensional case we perform the
- // integration using a QGaussLogR quadrature formula,
- // which is specifically designed to integrate
- // logarithmic singularities on the unit interval,
- // while in three dimensions we use the
- // QGaussOneOverR, which allows us to integrate 1/R
- // singularities on the vertices of the reference
- // element. Since we don't want to rebuild the two
- // dimensional quadrature formula at each singular
- // integration, we built them outside the loop on the
- // cells, and we only use a pointer to that quadrature
- // here.
- //
- // Notice that in one dimensional integration this is
- // not possible, since we need to know the scaling
- // parameter for the quadrature, which is not known a
- // priori.
- //
- // Dimension independent programming here is a little
- // tricky, but can be achieved via dynamic casting. We
- // check that everything went ok with an assertion at
- // the end of this block. Notice that the dynamic cast
- // will only work when the dimension is the correct
- // one, in which case it is possible to cast a
- // QGaussLogR and QGaussOneOverR to a Quadrature<1>
- // and Quadrature<2> object.
- //
- // In the other cases this won't be called, and even
- // if it was, the dynamic_cast function would just
- // return a null pointer. We check that this is not
- // the case with the Assert at the end.
- //
- // Notice that in two dimensions the singular
- // quadrature rule depends also on the size of the
- // current cell. For this reason, it is necessary to
- // create a new quadrature for each singular
- // integration. Since we create it using the new
- // operator of C++, we also need to destroy it using
- // the dual of new: delete. This is done at the end,
- // and only if dim == 2.
- Quadrature<dim-1> * singular_quadrature;
- for(unsigned int i=0; i<fe.dofs_per_cell; ++i) {
- if(dim == 2) {
- singular_quadrature = dynamic_cast<Quadrature<dim-1> *>(
- new QGaussLogR<1>(quadrature.size(),
- Point<1>((double)i),
- 1./cellj->measure()));
- } else {
- singular_quadrature = dynamic_cast<Quadrature<dim-1> *>(
- & sing_quadratures_3d[i]);
- }
-
- Assert(singular_quadrature, ExcInternalError());
-
-
- FEValues<dim-1,dim> fe_v_singular(fe, *singular_quadrature,
- update_jacobians |
- update_values |
- update_cell_normal_vectors |
- update_quadrature_points );
- fe_v_singular.reinit(cellj);
-
- static vector<Vector<double> > singular_cell_wind( (*singular_quadrature).size(),
- Vector<double>(dim) );
-
- const vector<Point<dim> > &singular_normals = fe_v_singular.get_cell_normal_vectors();
- const vector<Point<dim> > &singular_q_points = fe_v_singular.get_quadrature_points();
-
- wind.vector_value_list(singular_q_points, singular_cell_wind);
-
- for(unsigned int q=0; q<singular_quadrature->size(); ++q) {
- R = celli->vertex(i)-singular_q_points[q];
- double normal_wind = 0;
- for(unsigned int d=0; d<dim; ++d)
- normal_wind += (singular_cell_wind[q](d)*
- singular_normals[q][d]);
-
- local_rhs(i) += ( kernel.single_layer(R, is_singular) *
- normal_wind *
- fe_v_singular.JxW(q) );
-
- for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
- local_matrix(i,j) += (( kernel.double_layer(R, is_singular) *
- singular_normals[q]) *
- fe_v_singular.shape_value(j,q) *
- fe_v_singular.JxW(q) );
- }
- }
- if(dim==2) delete singular_quadrature;
- }
- }
- // Move the local matrix and rhs to the global one.
- 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.add(dofs_i[i],dofs_j[j], local_matrix(i,j));
- }
- }
+ ExcDimensionMismatch(fe.dofs_per_cell,
+ GeometryInfo<dim-1>::vertices_per_cell));
+
+ // Now that we checked that the number of vertices is equal to the
+ // number of degrees of freedom, we construct a vector of support
+ // points which will be used in the local integrations.
+ std::vector<Point<dim> > support_points(dh.n_dofs());
+ DoFTools::map_dofs_to_support_points<dim-1, dim>( StaticMappingQ1<dim-1, dim>::mapping,
+ dh, support_points);
+
+ for(cell = dh.begin_active(); cell != endc; ++cell) {
+
+ fe_v.reinit(cell);
+ cell->get_dof_indices(dofs);
+
+ const vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
+ const vector<Point<dim> > &normals = fe_v.get_cell_normal_vectors();
+ wind.vector_value_list(q_points, cell_wind);
+
+
+ for(unsigned int i=0; i<dh.n_dofs() ; ++i) {
+
+ local_matrix_row_i = 0;
+
+ // The integral that we need to perform is singular if one
+ // of the local degrees of freedom is the same of the
+ // support point i. Here we check wether this is the case,
+ // and we store which one is the singular index.
+ bool is_singular = false;
+ unsigned int singular_index = numbers::invalid_unsigned_int;
+
+ for(unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ if(dofs[j] == i) {
+ singular_index = j;
+ is_singular = true;
+ }
+
+ if(is_singular == false) {
+ for(unsigned int q=0; q<n_q_points; ++q) {
+ normal_wind = 0;
+ for(unsigned int d=0; d<dim; ++d)
+ normal_wind += normals[q][d]*cell_wind[q](d);
+
+ // Distance between the external support point
+ // and the quadrature point on the internal
+ // cell.
+ R = q_points[q] - support_points[i];
+
+ system_rhs(i) += ( kernel.single_layer(R) *
+ normal_wind *
+ fe_v.JxW(q) );
+
+ for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
+
+ local_matrix_row_i(j) += ( ( kernel.double_layer(R) *
+ normals[q] ) *
+ fe_v.shape_value(j,q) *
+ fe_v.JxW(q) );
+ }
+ }
+ } else {
+ // Now we treat the more delicate case. If we are
+ // here, it means that the cell that runs on the j
+ // index contains the support_point[i]. In this case
+ // both the single and the double layer potential are
+ // singular, and they require a special treatment, as
+ // explained in the introduction.
+ //
+ // In the two dimensional case we perform the integration
+ // using a QGaussLogR quadrature formula, which is
+ // specifically designed to integrate logarithmic
+ // singularities on the unit interval, while in three
+ // dimensions we use the QGaussOneOverR, which allows us to
+ // integrate 1/R singularities on the vertices of the
+ // reference element. Since we don't want to rebuild the two
+ // dimensional quadrature formula at each singular
+ // integration, we built them outside the loop on the cells,
+ // and we only use a pointer to that quadrature here.
+ //
+ // Notice that in one dimensional integration this is not
+ // possible, since we need to know the scaling parameter for
+ // the quadrature, which is not known a priori.
+ //
+ // Dimension independent programming here is a little tricky,
+ // but can be achieved via dynamic casting. We check that
+ // everything went ok with an assertion at the end of this
+ // block. Notice that the dynamic cast will only work when the
+ // dimension is the correct one, in which case it is possible
+ // to cast a QGaussLogR and QGaussOneOverR to a Quadrature<1>
+ // and Quadrature<2> object.
+ //
+ // In the other cases this won't be called, and even if it
+ // was, the dynamic_cast function would just return a null
+ // pointer. We check that this is not the case with the Assert
+ // at the end.
+ //
+ // Notice that in two dimensions the singular quadrature rule
+ // depends also on the size of the current cell. For this
+ // reason, it is necessary to create a new quadrature for each
+ // singular integration. Since we create it using the new
+ // operator of C++, we also need to destroy it using the dual
+ // of new: delete. This is done at the end, and only if dim ==
+ // 2.
+ Assert(singular_index != numbers::invalid_unsigned_int,
+ ExcInternalError());
+
+ Quadrature<dim-1> * singular_quadrature;
+
+ if(dim == 2) {
+ singular_quadrature = dynamic_cast<Quadrature<dim-1> *>(
+ new QGaussLogR<1>(singular_quadrature_order,
+ Point<1>((double)singular_index),
+ 1./cell->measure()));
+ } else {
+ singular_quadrature = dynamic_cast<Quadrature<dim-1> *>(
+ & sing_quadratures_3d[singular_index]);
+ }
+
+ Assert(singular_quadrature, ExcInternalError());
+
+ FEValues<dim-1,dim> fe_v_singular (fe, *singular_quadrature,
+ update_jacobians |
+ update_values |
+ update_cell_normal_vectors |
+ update_quadrature_points );
+
+ fe_v_singular.reinit(cell);
+
+ static vector<Vector<double> > singular_cell_wind( (*singular_quadrature).size(),
+ Vector<double>(dim) );
+
+ const vector<Point<dim> > &singular_normals = fe_v_singular.get_cell_normal_vectors();
+ const vector<Point<dim> > &singular_q_points = fe_v_singular.get_quadrature_points();
+
+ wind.vector_value_list(singular_q_points, singular_cell_wind);
+
+ for(unsigned int q=0; q<singular_quadrature->size(); ++q) {
+ R = singular_q_points[q]- support_points[i];
+ double normal_wind = 0;
+ for(unsigned int d=0; d<dim; ++d)
+ normal_wind += (singular_cell_wind[q](d)*
+ singular_normals[q][d]);
+
+ system_rhs(i) += ( kernel.single_layer(R, is_singular) *
+ normal_wind *
+ fe_v_singular.JxW(q) );
+
+ for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
+ local_matrix_row_i(j) += (( kernel.double_layer(R, is_singular) *
+ singular_normals[q]) *
+ fe_v_singular.shape_value(j,q) *
+ fe_v_singular.JxW(q) );
+ }
+ }
+ if(dim==2) {
+ delete singular_quadrature;
+ }
+ }
+
+ // Move the local matrix row to the global matrix.
+ for(unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ system_matrix.add(i,dofs[j], local_matrix_row_i(j));
+ }
}
- // One quick way to compute the matrix of the solid angles, is to
- // use the Neumann matrix itself. It is enough to multiply the
- // matrix with the vector of ones, to get the diagonal matrix of
- // the alpha solid angles.
- Vector<double> ones(dh.n_dofs()), alpha(dh.n_dofs());
+ // One quick way to compute the diagonal matrix of the solid
+ // angles, is to use the Neumann matrix itself. It is enough to
+ // multiply the matrix with the vector of ones, to get the
+ // diagonal matrix of the alpha solid angles.
+ Vector<double> ones(dh.n_dofs());
for(unsigned int i=0; i<dh.n_dofs(); ++i)
- ones(i) = 1.;
+ ones(i) = -1.;
system_matrix.vmult(alpha, ones);
- for(unsigned int i=0; i<dh.n_dofs(); ++i)
- system_matrix.add(i,i,-alpha(i));
+ for(unsigned int i = 0; i<dh.n_dofs(); ++i) {
+ system_matrix.add(i,i,alpha(i));
+ }
}
template <int dim>
LU.vmult(phi, system_rhs);
// Since we are solving a purely Neumann problem, the solution is
- // only known up to a constant potential. We filter out the mean
- // value using the MeanValueFilter class.
- MeanValueFilter mean_filter;
- mean_filter.filter(phi);
+ // only known up to a constant potential. We solve this issue by
+ // subtracting the mean value of the vector from each vector
+ // entry.
+ double mean_value = phi.mean_value();
+ for(unsigned int i=0; i<phi.size(); ++i)
+ phi(i) -= mean_value;
}
external_phi.reinit(external_dh.n_dofs());
typename DoFHandler<dim-1,dim>::active_cell_iterator
- cell = dh.begin_active(),
- endc = dh.end();
+ cell = dh.begin_active(),
+ endc = dh.end();
+
-
Quadrature<dim-1> &quadrature = *quadrature_pointer;
FEValues<dim-1,dim> fe_v(fe, quadrature,
- update_values |
- update_cell_normal_vectors |
- update_quadrature_points |
- update_JxW_values);
+ update_values |
+ update_cell_normal_vectors |
+ update_quadrature_points |
+ update_JxW_values);
const unsigned int n_q_points = fe_v.n_quadrature_points;
vector<unsigned int> dofs(fe.dofs_per_cell);
vector<double> local_phi(n_q_points);
+ vector<double> normal_wind(n_q_points);
vector<Vector<double> > local_wind(n_q_points, Vector<double>(dim) );
- double normal_wind;
LaplaceKernel<dim> kernel;
Point<dim> R;
typename DoFHandler<dim>::active_cell_iterator
- external_cell = external_dh.begin_active(),
- external_endc = external_dh.end();
-
- vector<unsigned int> external_dofs(external_fe.dofs_per_cell);
- vector<bool> dof_is_treated(external_dh.n_dofs(), false);
-
-
- for(; external_cell != external_endc; ++external_cell) {
- external_cell->get_dof_indices(external_dofs);
-
- for(unsigned int i=0; i<external_fe.dofs_per_cell; ++i)
- if(dof_is_treated[external_dofs[i]] == false) {
-
- dof_is_treated[external_dofs[i]] = true;
-
- external_phi(external_dofs[i]) = 0;
-
- for(cell = dh.begin_active(); cell != endc; ++cell) {
- fe_v.reinit(cell);
-
- const vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
- const vector<Point<dim> > &normals = fe_v.get_cell_normal_vectors();
-
- cell->get_dof_indices(dofs);
- fe_v.get_function_values(phi, local_phi);
-
- wind.vector_value_list(q_points, local_wind);
-
- for(unsigned int q=0; q<n_q_points; ++q) {
- normal_wind = 0;
- for(unsigned int d=0; d<dim; ++d)
- normal_wind += normals[q][d]*local_wind[q](d);
-
- R = external_cell->vertex(i) - q_points[q];
-
- external_phi(external_dofs[i]) += ( ( kernel.single_layer(R) *
- normal_wind -
- //
- ( kernel.double_layer(R) *
- normals[q] ) *
- local_phi[q] ) *
- fe_v.JxW(q) );
- }
- }
- }
+ external_cell = external_dh.begin_active(),
+ external_endc = external_dh.end();
+
+ std::vector<Point<dim> > external_support_points(external_dh.n_dofs());
+ DoFTools::map_dofs_to_support_points<dim>( StaticMappingQ1<dim>::mapping,
+ external_dh, external_support_points);
+
+ for(cell = dh.begin_active(); cell != endc; ++cell) {
+ fe_v.reinit(cell);
+
+ const vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
+ const vector<Point<dim> > &normals = fe_v.get_cell_normal_vectors();
+
+ cell->get_dof_indices(dofs);
+ fe_v.get_function_values(phi, local_phi);
+
+ wind.vector_value_list(q_points, local_wind);
+
+ for(unsigned int q=0; q<n_q_points; ++q){
+ normal_wind[q] = 0;
+ for(unsigned int d=0; d<dim; ++d)
+ normal_wind[q] += normals[q][d]*local_wind[q](d);
+ }
+
+ for(unsigned int i=0; i<external_dh.n_dofs(); ++i) {
+
+ for(unsigned int q=0; q<n_q_points; ++q) {
+
+ R = q_points[q] - external_support_points[i];
+
+ external_phi(i) += ( ( kernel.single_layer(R) *
+ normal_wind[q] +
+ //
+ (- kernel.double_layer(R) *
+ normals[q] ) *
+ local_phi[q] ) *
+ fe_v.JxW(q) );
+ }
+ }
}
DataOut<dim, DoFHandler<dim> > dataout;
dataout.attach_dof_handler(dh);
dataout.add_data_vector(phi, "phi");
+ dataout.add_data_vector(alpha, "alpha");
dataout.build_patches();
std::string filename = ( Utilities::int_to_string(dim) +
- "d_boundary_solution_" +
- Utilities::int_to_string(cycle) +
- ".vtk" );
+ "d_boundary_solution_" +
+ Utilities::int_to_string(cycle) +
+ ".vtk" );
std::ofstream file(filename.c_str());
dataout.write_vtk(file);
void BEMProblem<dim>::run() {
read_parameters("parameters.prm");
- read_domain();
-
- for(unsigned int cycle=0; cycle<n_cycles; ++cycle) {
- refine_and_resize();
- assemble_system();
- solve_system();
- output_results(cycle);
+ if(run_in_this_dimension == true) {
+ read_domain();
+
+ for(unsigned int cycle=0; cycle<n_cycles; ++cycle) {
+ refine_and_resize();
+ assemble_system();
+ solve_system();
+ output_results(cycle);
+ }
+ if(extend_solution == true)
+ interpolate();
+ } else {
+ deallog << "Run in dimension " << dim
+ << " explicitly disabled in parameter file. "
+ << std::endl;
}
- if(extend_solution == true)
- interpolate();
}
projects / dealii-svn.git / commitdiff
