From 8f1782d82b1098688272fdbb739008f7bc182ae4 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Sun, 5 Apr 2009 18:55:27 +0000 Subject: [PATCH] Step-34: Rewrote almost all formulas in a consistent manner (almost all signs were wrong). Removed all references to Galerkin BEM. Now it is a purely collocation BEM, that uses advanced integration rules just added to deal.II. Once MappingQ of arbitrary order are implemented, this should allow for very advanced BEM solutions, with very few unknowns. Added a couple of pictures of solutions, and the parameter file at the end. Still missing is a completion of Wolfgang addition to the introduction. git-svn-id: https://svn.dealii.org/trunk@18552 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-34/doc/intro.dox | 327 ++++------ deal.II/examples/step-34/doc/results.dox | 80 ++- deal.II/examples/step-34/parameters.prm | 19 +- deal.II/examples/step-34/step-34.cc | 752 ++++++++++++----------- 4 files changed, 618 insertions(+), 560 deletions(-) diff --git a/deal.II/examples/step-34/doc/intro.dox b/deal.II/examples/step-34/doc/intro.dox index 1d559028a0..ab75eb9db4 100644 --- a/deal.II/examples/step-34/doc/intro.dox +++ b/deal.II/examples/step-34/doc/intro.dox @@ -24,23 +24,22 @@ usually modeled by the Euler equations of fluid dynamics: 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} &= @@ -57,15 +56,19 @@ boundary conditions \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 outer 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 @@ -83,7 +86,7 @@ which we find that the same equations hold (because $\nabla\cdot \f[ \label{eq:boundary-conditions-tilde} \begin{aligned} - \mathbf{n}\cdot\mathbf{\tilde{v}}& = -\mathbf{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] @@ -101,13 +104,12 @@ as the homogenous Laplace equation for the unknown $\phi$: \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. @@ -115,12 +117,12 @@ We will now reformulate this equation in integral form using the Green identity: \f[\label{green} \int_{\mathbb{R}^n\backslash\Omega} - (-\Delta u)v\,dx - \int_{\partial\Omega} \frac{\partial u}{\partial \mathbf{n} }v \,ds + (\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 outward from the domain of integration whereas in our case it is pointing inward; the difference is only in the sign @@ -130,34 +132,35 @@ called fundamental solutions of the Laplace equation, \f[ \begin{aligned} \label{eq:3} - G(\mathbf{x}-\mathbf{y}) = & - -\frac{1}{2\pi}\ln|\mathbf{x}-\mathbf{y}| \qquad && \text{for } n=2 + G(\mathbf{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 @@ -166,13 +169,12 @@ 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$. +$-\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 @@ -180,30 +182,27 @@ operators, we obtain an equation for $\phi$ just on the boundary of $\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 @@ -215,11 +214,17 @@ $\alpha(\mathbf{x})$ by which the point $\mathbf{x}$ sees the domain $\Omega$ can be defined using the double layer potential itself: \f[ \alpha(\mathbf{x}) := -\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 @@ -227,9 +232,10 @@ $\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}) - + + \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 @@ -237,6 +243,11 @@ 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$. +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 @@ -330,7 +341,7 @@ The goal of the following test problem is to solve the integral 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 @@ -352,7 +363,7 @@ $\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 \}, + \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 @@ -362,28 +373,27 @@ and often in the literature only piecewise constant elements are used.

Collocation boundary element method

-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 @@ -393,27 +403,33 @@ If the $\mathbf{x}_i$ support points are chosen correctly, then the 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 @@ -423,137 +439,32 @@ boundary element $\hat K := [0,1]^{n-1}$, and we perform the integrations after change of variables from the real element $K_i$ to the reference element $\hat K$. -

Singular integrals in two dimension.

+

Treating the singular integrals.

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). - - -

Singular integrals in three dimensions

- -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. -

What the program does

+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. diff --git a/deal.II/examples/step-34/doc/results.dox b/deal.II/examples/step-34/doc/results.dox index a1e1e35249..067230e652 100644 --- a/deal.II/examples/step-34/doc/results.dox +++ b/deal.II/examples/step-34/doc/results.dox @@ -1,3 +1,81 @@

Results

-

Simple Example

+We ran the program using the following parameters.prm file: +@verbatim +# Listing of Parameters +# --------------------- +set Extend solution on the -2,2 box = true +set External refinement = 5 +set Number of cycles = 4 + + +subsection Quadrature rules + set Quadrature order = 4 + set Quadrature type = gauss + set Singular quadrature order = 5 +end + + +subsection Wind function 2d + # Any constant used inside the function which is not a variable name. + set Function constants = + + # Separate vector valued expressions by ';' as ',' is used internally by the + # function parser. + set Function expression = 1; 1 + + # The name of the variables as they will be used in the function, separated + # by ','. + set Variable names = x,y,t +end + + +subsection Wind function 3d + # Any constant used inside the function which is not a variable name. + set Function constants = + + # Separate vector valued expressions by ';' as ',' is used internally by the + # function parser. + set Function expression = 1; 1; 1 + + # The name of the variables as they will be used in the function, separated + # by ','. + set Variable names = x,y,z,t +end +@endverbatim + +When we run the program, the following is printed on screen: +@verbatim + +DEAL::Parsing parameter file parameters.prm +DEAL::for a 2 dimensional simulation. +DEAL::Levels: 2, potential dofs: 20 +DEAL::Levels: 3, potential dofs: 40 +DEAL::Levels: 4, potential dofs: 80 +DEAL::Levels: 5, potential dofs: 160 + +DEAL::Parsing parameter file parameters.prm +DEAL::for a 3 dimensional simulation. +DEAL::Levels: 2, potential dofs: 26 +DEAL::Levels: 3, potential dofs: 98 +DEAL::Levels: 4, potential dofs: 386 +DEAL::Levels: 5, potential dofs: 1538 +@endverbatim + +The result of running these computations is a bunch of output files that we +can pass to our visualization program of choice. + +The output files are of two kind: the potential on the boundary +element surface, and the potential extended to the outer and inner +domain. The combination of the two for the two dimensional case look +like + +@image html step-34_2d.png + +while in three dimensions we only show the potential on the surface, +together with a countur plot: + +@image html step-34_3d.png + +The exact solution for this model problem should be $x+y$ or +$x+y+z$ on the surface, which is precisely what we obtain. \ No newline at end of file diff --git a/deal.II/examples/step-34/parameters.prm b/deal.II/examples/step-34/parameters.prm index b61ad55275..aaf72cb5b3 100644 --- a/deal.II/examples/step-34/parameters.prm +++ b/deal.II/examples/step-34/parameters.prm @@ -1,13 +1,16 @@ # Listing of Parameters # --------------------- -set Extend solution on the -2,2 box = false -set External refinement = 5 -set Number of cycles = 3 +set Extend solution on the -2,2 box = true +set External refinement = 8 +set Number of cycles = 4 +set Run 2d simulation = true +set Run 3d simulation = false -subsection Quadrature rule - set Quadrature order = 5 - set Quadrature type = gauss +subsection Quadrature rules + set Quadrature order = 4 + set Quadrature type = gauss + set Singular quadrature order = 5 end @@ -17,7 +20,7 @@ subsection Wind function 2d # Separate vector valued expressions by ';' as ',' is used internally by the # function parser. - set Function expression = 1; 0 + set Function expression = 1; 1 # default: 0; 0 # The name of the variables as they will be used in the function, separated # by ','. @@ -31,7 +34,7 @@ subsection Wind function 3d # Separate vector valued expressions by ';' as ',' is used internally by the # function parser. - set Function expression = 1; 0; 0 + set Function expression = 1; 1; 1 # default: 0; 0; 0 # The name of the variables as they will be used in the function, separated # by ','. diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index 0b7b0924b9..77e7f27166 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -91,28 +91,30 @@ public: // (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 @@ -143,31 +145,32 @@ private: // argument, and produces the usual finite element classes that we // saw in all previous examples. - Triangulation tria; - FE_Q fe; - DoFHandler dh; + Triangulation tria; + FE_Q fe; + DoFHandler 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 system_matrix; - Vector system_rhs; - Vector phi; + SparsityPattern sparsity; + SparseMatrix system_matrix; + Vector system_rhs; + Vector phi; + Vector 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 external_tria; - FE_Q external_fe; - DoFHandler external_dh; - Vector external_phi; + Triangulation external_tria; + FE_Q external_fe; + DoFHandler external_dh; + Vector external_phi; // The following variables are the one that we fill through a // parameter file. The new objects that we use in this example @@ -182,7 +185,9 @@ private: // 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 @@ -194,8 +199,10 @@ private: // wanted to extend the solution to the entire domain. Functions::ParsedFunction wind; SmartPointer > quadrature_pointer; + unsigned int singular_quadrature_order; unsigned int n_cycles; unsigned int external_refinement; + bool run_in_this_dimension; bool extend_solution; }; @@ -207,7 +214,7 @@ class LaplaceKernel 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 @@ -215,22 +222,54 @@ public: // 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 &R, - bool factor_out_singularity = false); + bool factor_out_2d_singularity = false); Point double_layer(const Point &R, - bool factor_out_singularity = false); + bool factor_out_2d_singularity = false); }; @@ -246,40 +285,54 @@ BEMProblem::BEMProblem() : template void BEMProblem::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 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(); @@ -289,42 +342,45 @@ void BEMProblem::read_parameters(std::string filename) { template double LaplaceKernel::single_layer(const Point &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 Point LaplaceKernel::double_layer(const Point &R, - bool factor_out_singularity) { + bool factor_out_2d_singularity) { Point 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 void BEMProblem::read_domain() { @@ -337,10 +393,9 @@ void BEMProblem::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 @@ -362,11 +417,11 @@ void BEMProblem::read_domain() { GridIn 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); } @@ -382,7 +437,7 @@ void BEMProblem::refine_and_resize() { 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 @@ -391,22 +446,22 @@ void BEMProblem::refine_and_resize() { system_matrix.clear(); sparsity.reinit(ndofs, ndofs, ndofs); for(unsigned int i=0; i void BEMProblem::assemble_system() { typename DoFHandler::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 @@ -415,33 +470,40 @@ void BEMProblem::assemble_system() { Quadrature &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 > 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 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 dofs_i(fe.dofs_per_cell); - vector dofs_j(fe.dofs_per_cell); + vector dofs(fe.dofs_per_cell); vector > cell_wind(n_q_points, Vector(dim) ); double normal_wind; - Vector local_rhs(fe.dofs_per_cell); - FullMatrix 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 local_matrix_row_i(fe.dofs_per_cell); // The kernel. LaplaceKernel kernel; @@ -452,183 +514,184 @@ void BEMProblem::assemble_system() { // 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::vertices_per_cell. AssertThrow(fe.dofs_per_cell == GeometryInfo::vertices_per_cell, - ExcDimensionMismatch(fe.dofs_per_cell, GeometryInfo::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 > &q_points = fe_v.get_quadrature_points(); - const vector > &normals = fe_v.get_cell_normal_vectors(); - wind.vector_value_list(q_points, cell_wind); - - if(is_singular == false) { - for(unsigned int q=0; qvertex(i)-q_points[q]; - - local_rhs(i) += ( kernel.single_layer(R) * - normal_wind * - fe_v.JxW(q) ); - - for(unsigned int j=0; j - // 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 * singular_quadrature; - for(unsigned int i=0; i *>( - new QGaussLogR<1>(quadrature.size(), - Point<1>((double)i), - 1./cellj->measure())); - } else { - singular_quadrature = dynamic_cast *>( - & sing_quadratures_3d[i]); - } - - Assert(singular_quadrature, ExcInternalError()); - - - FEValues fe_v_singular(fe, *singular_quadrature, - update_jacobians | - update_values | - update_cell_normal_vectors | - update_quadrature_points ); - fe_v_singular.reinit(cellj); - - static vector > singular_cell_wind( (*singular_quadrature).size(), - Vector(dim) ); - - const vector > &singular_normals = fe_v_singular.get_cell_normal_vectors(); - const vector > &singular_q_points = fe_v_singular.get_quadrature_points(); - - wind.vector_value_list(singular_q_points, singular_cell_wind); - - for(unsigned int q=0; qsize(); ++q) { - R = celli->vertex(i)-singular_q_points[q]; - double normal_wind = 0; - for(unsigned int d=0; d::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 > support_points(dh.n_dofs()); + DoFTools::map_dofs_to_support_points( StaticMappingQ1::mapping, + dh, support_points); + + for(cell = dh.begin_active(); cell != endc; ++cell) { + + fe_v.reinit(cell); + cell->get_dof_indices(dofs); + + const vector > &q_points = fe_v.get_quadrature_points(); + const vector > &normals = fe_v.get_cell_normal_vectors(); + wind.vector_value_list(q_points, cell_wind); + + + for(unsigned int i=0; i + // 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 * singular_quadrature; + + if(dim == 2) { + singular_quadrature = dynamic_cast *>( + new QGaussLogR<1>(singular_quadrature_order, + Point<1>((double)singular_index), + 1./cell->measure())); + } else { + singular_quadrature = dynamic_cast *>( + & sing_quadratures_3d[singular_index]); + } + + Assert(singular_quadrature, ExcInternalError()); + + FEValues fe_v_singular (fe, *singular_quadrature, + update_jacobians | + update_values | + update_cell_normal_vectors | + update_quadrature_points ); + + fe_v_singular.reinit(cell); + + static vector > singular_cell_wind( (*singular_quadrature).size(), + Vector(dim) ); + + const vector > &singular_normals = fe_v_singular.get_cell_normal_vectors(); + const vector > &singular_q_points = fe_v_singular.get_quadrature_points(); + + wind.vector_value_list(singular_q_points, singular_cell_wind); + + for(unsigned int q=0; qsize(); ++q) { + R = singular_q_points[q]- support_points[i]; + double normal_wind = 0; + for(unsigned int d=0; d 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 ones(dh.n_dofs()); for(unsigned int i=0; i @@ -638,10 +701,12 @@ void BEMProblem::solve_system() { 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::interpolate() { external_phi.reinit(external_dh.n_dofs()); typename DoFHandler::active_cell_iterator - cell = dh.begin_active(), - endc = dh.end(); + cell = dh.begin_active(), + endc = dh.end(); + - Quadrature &quadrature = *quadrature_pointer; FEValues 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 dofs(fe.dofs_per_cell); vector local_phi(n_q_points); + vector normal_wind(n_q_points); vector > local_wind(n_q_points, Vector(dim) ); - double normal_wind; LaplaceKernel kernel; Point R; typename DoFHandler::active_cell_iterator - external_cell = external_dh.begin_active(), - external_endc = external_dh.end(); - - vector external_dofs(external_fe.dofs_per_cell); - vector 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 > &q_points = fe_v.get_quadrature_points(); - const vector > &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; qvertex(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 > external_support_points(external_dh.n_dofs()); + DoFTools::map_dofs_to_support_points( StaticMappingQ1::mapping, + external_dh, external_support_points); + + for(cell = dh.begin_active(); cell != endc; ++cell) { + fe_v.reinit(cell); + + const vector > &q_points = fe_v.get_quadrature_points(); + const vector > &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 > dataout; @@ -747,12 +806,13 @@ void BEMProblem::output_results(unsigned int cycle) { 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); @@ -762,16 +822,22 @@ template void BEMProblem::run() { read_parameters("parameters.prm"); - read_domain(); - - for(unsigned int cycle=0; cycle