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