From: heltai Date: Tue, 3 Mar 2009 10:58:40 +0000 (+0000) Subject: Enanced comments for step-34, rewrote assembly routine, added external parameter... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=be2763d125f2e827edad17c1a6a0a66397633f88;p=dealii-svn.git Enanced comments for step-34, rewrote assembly routine, added external parameter file. git-svn-id: https://svn.dealii.org/trunk@18443 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-34/Makefile b/deal.II/examples/step-34/Makefile index 96fe808f62..3dc5c73fbf 100644 --- a/deal.II/examples/step-34/Makefile +++ b/deal.II/examples/step-34/Makefile @@ -153,4 +153,3 @@ Makefile.dep: $(target).cc Makefile \ # To make the dependencies known to `make', we finally have to include # them: include Makefile.dep - diff --git a/deal.II/examples/step-34/doc/intro.dox b/deal.II/examples/step-34/doc/intro.dox index 0d85ffb09c..9eddab7663 100644 --- a/deal.II/examples/step-34/doc/intro.dox +++ b/deal.II/examples/step-34/doc/intro.dox @@ -1,9 +1,394 @@ +

Introduction

-

The problem

+

Irrotational Flow

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

The numerical approximation

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

Galerkin boundary element method

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

Singular Integrals in two dimension.

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

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

Galerkin BEM

What the program does

diff --git a/deal.II/examples/step-34/parameters.prm b/deal.II/examples/step-34/parameters.prm new file mode 100644 index 0000000000..b4995437f6 --- /dev/null +++ b/deal.II/examples/step-34/parameters.prm @@ -0,0 +1,25 @@ +# Listing of Parameters +# --------------------- + +set Number of cycles = 3 + +subsection Outer quadrature rule + set Quadrature order = 0 + set Quadrature type = midpoint +end + + +subsection Wind function + # Any constant used inside the function which is not a variable name. + set Function constants = angle=0 + + # Separate vector valued expressions by ';' as ',' is used internally by the + # function parser. + set Function expression = cos(angle); sin(angle); 0 + + # The name of the variables as they will be used in the function, separated + # by ','. + set Variable names = x,y,z,t +end + + diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index abb217697d..c39ccb761d 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -1,4 +1,4 @@ -//---------------------------- bem_integration.cc --------------------------- +//---------------------------- step-34.cc --------------------------- // $Id: testsuite.html 13373 2006-07-13 13:12:08Z manigrasso $ // Version: $Name$ // @@ -11,7 +11,7 @@ // // Authors: Luca Heltai, Cataldo Manigrasso // -//---------------------------- bem_integration.cc --------------------------- +//---------------------------- step-34.cc --------------------------- // @@ -21,10 +21,12 @@ #include #include +#include #include #include #include #include +#include #include #include #include @@ -59,27 +61,53 @@ class LaplaceKernelIntegration { public: - LaplaceKernelIntegration(); + LaplaceKernelIntegration(FiniteElement &fe); ~LaplaceKernelIntegration(); void run(); - void compute_SD_integral_on_cell(vector &dst, - typename DoFHandler::active_cell_iterator &cell, - const Point &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 > > &dst, + typename DoFHandler::active_cell_iterator &cell, + const vector > &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 > fe_values; + SmartPointer > fe_values; }; template @@ -89,71 +117,74 @@ public: 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 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 fe; + FESystem fev; - /** The potential degrees of freedom. */ + // Finite element space used to smoothen the potential solution + // (from piecewise constant to continuous piecewise quadratic) + FE_Q fe_q; + + // And the relevant degrees of freedom. DoFHandler dh; - - /** The velocity degrees of freedom. */ DoFHandler dhv; + DoFHandler 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 > system_matrix; - /** Single layer potential matrix. */ - FullMatrix single_layer_matrix; - - /** Normal component of the wind field. */ - Vector Vn; - - /** System rhs. */ - Vector rhs; - - /** Potential. */ + // The right hand side, the potential and its smoothed version + Vector system_rhs; Vector phi; - - /** Something else. */ - Vector vs; + Vector 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 wind; + SmartPointer > quadrature_pointer; + unsigned int n_cycles; }; template BEMProblem::BEMProblem() : fe(0), fev(FE_DGP(0), dim), + fe_q(FE_Q(2)), dh(tria), - dhv(tria) + dhv(tria), + dhq(tria), + wind(dim) {} - template BEMProblem::~BEMProblem() { LAPACKFullMatrix * p = system_matrix; @@ -162,10 +193,49 @@ BEMProblem::~BEMProblem() { } +template +void BEMProblem::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::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 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 > qps(5); qps[0] = Point<2>(0,0); qps[1] = Point<2>(0,1); @@ -182,50 +252,24 @@ LaplaceKernelIntegration<2>::LaplaceKernelIntegration() template LaplaceKernelIntegration::~LaplaceKernelIntegration() { - FEValues *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 *fp = fe_values; fe_values = 0; delete fp; } template <> -void -LaplaceKernelIntegration<2>::compute_SD_integral_on_cell(vector &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 > jacobians = fe_values->get_jacobians(); - vector > 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 i_S(4); - vector 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 double -LaplaceKernelIntegration::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; @@ -245,9 +289,9 @@ LaplaceKernelIntegration::term_S (const Point<3> &r, } -template +template <> double -LaplaceKernelIntegration::term_D (const Point<3> &r, +LaplaceKernelIntegration<3>::term_D (const Point<3> &r, const Point<3> &a1, const Point<3> &a2) { @@ -264,6 +308,55 @@ LaplaceKernelIntegration::term_D (const Point<3> &r, } +template <> +void +LaplaceKernelIntegration<3>::compute_SD_integral_on_cell(vector > > &dstvv, + DoFHandler<2,3>::active_cell_iterator &cell, + const vector > &q_points) +{ + fe_values->reinit(cell); + vector > jacobians = fe_values->get_jacobians(); + vector > 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; idofs_per_cell; ++i) { + vector > & 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 &point = q_points[outer_q]; + vector &dst = dstv[outer_q]; + r_c = point-cell->center(); + double rn_c = r_c*n_c; + vector i_S(4); + vector 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 void BEMProblem::read_domain() { // Center of the ball. It is the origin by default. @@ -307,75 +400,132 @@ void BEMProblem::refine_and_resize() { system_matrix = new LAPACKFullMatrix(ndofs, ndofs); - Vn.reinit(ndofs); - rhs.reinit(ndofs); + system_rhs.reinit(ndofs); phi.reinit(ndofs); - vs.reinit(nvdofs); } template void BEMProblem::assemble_system() { - Point wind; - wind[0] = 1.; - typename DoFHandler::active_cell_iterator celli = dh.begin_active(), cellj = dh.begin_active(), - cellv = dhv.begin_active(), + cellvi = dhv.begin_active(), endc = dh.end(); - QMidpoint midpoint; - FEValues 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 &quadrature_outer = *quadrature_pointer; + QMidpoint quadrature_mid; + + + FEValues fe_outer(fe, quadrature_outer, + update_values | + update_cell_normal_vectors | + update_quadrature_points); + + FEValues 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 dofs_i(fe.dofs_per_cell); + vector dofs_j(fe.dofs_per_cell); + vector dofs_v_i(fev.dofs_per_cell); - vector dofsi(fe_mid.n_quadrature_points); - vector dofsj(fe_mid.n_quadrature_points); - vector dofsv(dim); + vector > > single_double_layer_potentials + (fe.dofs_per_cell, vector > + (n_q_points_outer, vector (2, 0.) ) ); + + vector > cell_wind(n_q_points_inner, Vector(dim) ); + vector normal_wind(n_q_points_inner); - // Temporary matrix - FullMatrix _B(dh.n_dofs(), dh.n_dofs()); + Vector local_rhs(fe.dofs_per_cell); + FullMatrix local_matrix(fe.dofs_per_cell, fe.dofs_per_cell); // The kernel. - LaplaceKernelIntegration kernel; - - // i runs on points, j runs on cells. - for(; celli != endc; ++celli, ++cellv) { - fe_mid.reinit(celli); - Point ci = celli->center(); - Point 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 kernel(fe); + + // i runs on outer integration, j runs on inner integration. + for(; celli != endc; ++celli, ++cellvi) { + fe_outer.reinit(celli); + + const vector > &q_points_outer = fe_outer.get_quadrature_points(); + const vector > &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 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 > &q_points_inner = fe_inner.get_quadrature_points(); + const vector > &normals_inner = fe_inner.get_cell_normal_vectors(); + wind.vector_value_list(q_points_inner, cell_wind); + + + for(unsigned int q=0; qindex() != 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 void BEMProblem::solve_system() { - phi.swap(rhs); + phi.swap(system_rhs); system_matrix->compute_lu_factorization(); system_matrix->apply_lu_factorization(phi, false); } @@ -399,11 +549,11 @@ void BEMProblem::output_results(unsigned int cycle) { template void BEMProblem::run() { - read_domain(); - const unsigned int number_of_cycles = 4; + read_parameters("parameters.prm"); + read_domain(); - for(unsigned int cycle=0; cycle