From: heltai Date: Mon, 23 Mar 2009 18:00:48 +0000 (+0000) Subject: Modified step-34 to implement exclusively the collocation boundary element X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cce158ed5fc6a9837db337d695cfef08b78d42d7;p=dealii-svn.git Modified step-34 to implement exclusively the collocation boundary element method. This is done now through the QGaussLog and QGaussOneOverR classes that were previously added to the library. Extended a little the comments, but there is still a lot of documentation to do. The solution is the correct one (phi(x) = x). TODO: Add convergence table for both the potential and the velocity. git-svn-id: https://svn.dealii.org/trunk@18503 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-34/doc/intro.dox b/deal.II/examples/step-34/doc/intro.dox index 98881f9775..1d559028a0 100644 --- a/deal.II/examples/step-34/doc/intro.dox +++ b/deal.II/examples/step-34/doc/intro.dox @@ -104,6 +104,12 @@ as the homogenous Laplace equation for the unknown $\phi$: \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, +\f] +and the pressure and velocity are uncoupled. We will now reformulate this equation in integral form using the Green identity: @@ -174,30 +180,46 @@ operators, we obtain an equation for $\phi$ just on the boundary of $\Omega$: \f[\label{SD} - \frac{1}{2}\phi(\mathbf{x}) = - \left(S [\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x}) + \alpha(\mathbf{x})\phi(\mathbf{x}) = + - \left(S [\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x}) + (D\phi)(\mathbf{x}) \quad \mathbf{x}\in \partial\Omega, \f] -which is the integral formulation we were looking for. Substituting the single +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$. + +Substituting the single and double layer operators we get: \f[ - \pi\phi(\mathbf{x})= - \int_{\partial \Omega} \ln|\mathbf{x}-\mathbf{y}| \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y + \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 + - \int_{\partial \Omega} \frac{ (\mathbf{x}-\mathbf{y})\cdot\mathbf{n}_y }{ |\mathbf{x}-\mathbf{y}|^2 }\,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 \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 + \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 + - \int_{\partial \Omega} \frac{ (\mathbf{x}-\mathbf{y})\cdot\mathbf{n}_y }{ |\mathbf{x}-\mathbf{y}|^3 }\phi(\mathbf{y})\,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 \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 case, $\phi$ is the solution of an integral equation posed entirely on the boundary since both $\mathbf{x},\mathbf{y}\in\partial\Omega$. +Notice that the fraction of angle (in 2d) or solid angle (in 3d) +$\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 +\f] + 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 @@ -216,7 +238,7 @@ $\phi$ on the boundary we have just computed). Finally, we can then recover the velocity as $\mathbf{\tilde v}=\nabla \phi$. As a final test, let us verify that this velocity indeed satisfies the -momentum balance equation for a stationary flow field, i.e. whether +momentum balance equation for a stationary flow field, i.e., whether $\mathbf{v}\cdot\nabla\mathbf{v} = -\frac 1\rho \nabla p$ where $\mathbf{v}=\mathbf{\tilde v}+\mathbf{v}_\infty=\nabla\phi+\mathbf{v}_\infty$ for some (unknown) pressure @@ -319,68 +341,79 @@ dimension $n$ of the surrounding space $\mathbb{R}^n$. We define the finite dimensional space $V_h$ as \f[ \label{eq:definition-Vh} - V_h := \{ v \in L^2(\Gamma) \text{ s.t. } v|_{K_i} \in \mathcal{P}^0(K_i), + V_h := \{ v \in C^0(\Gamma) \text{ s.t. } v|_{K_i} \in \mathcal{Q}^1(K_i), \forall i\}, \f] -with basis functions $\psi_i(\mathbf{x}) = \chi_{K_i}(\mathbf{x})$, -i.e., one if $\mathbf{x}$ belongs to $K_i$, and zero otherwise. An element -$\phi_h$ of $V_h$ is uniquely -identified by the vector $\boldsymbol{\alpha}$ of its coefficients +with basis functions $\psi_i(\mathbf{x})$ (the usual FE_Q finite element, +with the catch that this time it is defined on a manifold of codimension one). +An element $\phi_h$ of $V_h$ is uniquely +identified by the vector $\boldsymbol{\phi}$ 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 \}, \f] -where summation is implied over repeated indexes. Note that we use +where summation is implied over repeated indexes. Note that we could use discontinuous elements here — in fact, there is no real reason to use -continuous ones in the first place since the integral formulation does not -imply any derivatives on our trial functions so continuity is unnecessary. +continuous ones since the integral formulation does not +imply any derivatives on our trial functions so continuity is unnecessary, +and often in the literature only piecewise constant elements are used. +

Collocation boundary element method

-

Galerkin boundary element method

+By far, the most common approximation of boundary integral equations, +is by use of the collocation based boundary element method. -The usual Galerkin approach for the discretization of the above -problem gives us the following variational formulation: +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$. + +The problem then becomes: 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: +such that, the following $n_dofs$ equations are satisfied: + \f{align*} - \label{eq:galerkin-continuous} - \int_{\Gamma_x} \phi_h(\mathbf{x}) \eta(\mathbf{x})\,ds_x = - & 2\int_{\Gamma_x} \int_{\Gamma_y} - G(\mathbf{x}-\mathbf{y}) \, \mathbf{n}_y\cdot\mathbf{v_\infty} - \eta(\mathbf{x}) \,ds_x\,ds_y + \alpha(mathbf{x}_i) \phi_h(\mathbf{x}_i) = + & \int_{\Gamma_y} G(\mathbf{x}_i-\mathbf{y}) \, + \mathbf{n}_y\cdot\mathbf{v_\infty} \,ds_y \\ - & + 2\int_{\Gamma_x}\int_{\Gamma_y} \frac{ \partial - G(\mathbf{x}-\mathbf{y})}{\partial\mathbf{n}_y } - \phi_h(\mathbf{y})\eta(\mathbf{x}) \,ds_x\,ds_y. + & + \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 +explained above. -The linearity of the integral operator makes this problem equivalent -to solving the linear system +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{M}-\mathbf{A})\boldsymbol\alpha = \mathbf{b}, +(\mathbf{A}-\mathbf{N})\boldsymbol\phi = \mathbf{b}, \f] where \f[ \begin{aligned} -\mathbf{M}_{ij}&= |K_i|\delta_{ij}\\ -\mathbf{A}_{ij}&= 2\int_{K_i}\int_{K_j} - \frac{\partial G(\mathbf{x}-\mathbf{y})}{\partial \mathbf{n}_y} - \psi_i(\mathbf{x})\psi_j(\mathbf{y}) \,ds_x\,ds_y +\mathbf{A}_{ii}&= \int_\Gamma +\frac{\partial G(\mathbf{x}_i-\mathbf{y})}{\partial \mathbf{n}_y}\,ds_y \\ +\mathbf{N}_{ij}&= \int_\Gamma + \frac{\partial G(\mathbf{x}_i-\mathbf{y})}{\partial \mathbf{n}_y} + \psi_j(\mathbf{y}) \,ds_y \\ -\mathbf{b}_i&= 2\int_{K_i} \int_{\Gamma_{h,y}} - G(\mathbf{x}-\mathbf{y}) \, \mathbf{n}_y\cdot\mathbf{v_\infty} - \psi_i(\mathbf{y}) \,ds_x\,ds_y. +\mathbf{b}_i&= \int_\Gamma + G(\mathbf{x}_i-\mathbf{y}) \, \mathbf{n}_y\cdot\mathbf{v_\infty} + \psi_i(\mathbf{y}) 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$. +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$. As usual in these cases, all integrations are performed on a reference simple domain, i.e., we assume that each element $K_i$ of @@ -393,7 +426,7 @@ 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 +$\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 orthogonal (on our polygonal approximation of the boundary of $\Omega$), and diff --git a/deal.II/examples/step-34/parameters.prm b/deal.II/examples/step-34/parameters.prm index 9bd40dc837..b61ad55275 100644 --- a/deal.II/examples/step-34/parameters.prm +++ b/deal.II/examples/step-34/parameters.prm @@ -1,20 +1,16 @@ # Listing of Parameters # --------------------- -set Number of cycles = 3 -set External refinement = 5 +set Extend solution on the -2,2 box = false +set External refinement = 5 +set Number of cycles = 3 -subsection Inner quadrature rule - set Quadrature order = 2 +subsection Quadrature rule + set Quadrature order = 5 set Quadrature type = gauss end -subsection Outer quadrature rule - set Quadrature order = 0 - set Quadrature type = midpoint -end - subsection Wind function 2d # Any constant used inside the function which is not a variable name. set Function constants = @@ -28,6 +24,7 @@ subsection Wind function 2d 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 = diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index 5fdf13e53f..0b7b0924b9 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -45,6 +45,7 @@ #include #include #include +#include #include @@ -59,15 +60,14 @@ using namespace dealii; template -class LaplaceKernelIntegration; +class LaplaceKernel; template class BEMProblem { public: - BEMProblem(const unsigned int degree = 0); - ~BEMProblem(); + BEMProblem(); // The structure of a boundary element method code is very similar // to the structure of a finite element code. By now you should be @@ -144,7 +144,7 @@ private: // saw in all previous examples. Triangulation tria; - FE_DGP fe; + FE_Q fe; DoFHandler dh; // In BEM methods, the matrix that is generated is @@ -155,9 +155,10 @@ private: // system. Note that this will be very inefficient when the number // of dofs grows, since it is of order $n^3$. - SmartPointer > system_matrix; - Vector system_rhs; - Vector phi; + SparsityPattern sparsity; + SparseMatrix system_matrix; + Vector system_rhs; + Vector phi; // The reconstruction of the solution in the entire space is done // on a continuous finite element grid of dimension dim. These are @@ -169,134 +170,92 @@ private: 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 are the - // ParsedFunction object and the QuadratureSelector object. + // parameter file. The new objects that we use in this example + // are the ParsedFunction object and the QuadratureSelector + // object. // // The ParsedFunction class allows us to easily and quickly define // new function objects via parameter files, with custom - // definitions which can be very - // complex (see the documentation of that class for all the - // available options). + // definitions which can be very complex (see the documentation of + // that class for all the available options). // // 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 quadrature formulas for the inner as well as the - // outer integration in the calculation of the boundary element - // matrix. + // selection of the quadrature formulas for the inner integration. // - // Notice that selecting the midpoint rule as the outer - // integration formula on uniformly refined meshes is equivalent - // (up to a scaling factor) to solving the boundary element method - // via collocation instead of Galerkin technique. + // Notice that the pointer given below for the quadrature rule is + // only used for non singular integrals. Whenever the integral is + // singular, then only the degree of the quadrature pointer is + // used, and the integration is a special one (see the + // assemble_matrix below for further details). + // + // We also define a couple of parameters which are used in case we + // wanted to extend the solution to the entire domain. Functions::ParsedFunction wind; - SmartPointer > outer_quadrature_pointer; - SmartPointer > inner_quadrature_pointer; + SmartPointer > quadrature_pointer; unsigned int n_cycles; unsigned int external_refinement; + bool extend_solution; }; template -class LaplaceKernelIntegration +class LaplaceKernel { public: - - LaplaceKernelIntegration(const FiniteElement &fe); - ~LaplaceKernelIntegration(); - - // 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) rhs(y) dy - // \f] + // 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 + // different from zero. // - // and + // Whenever the integration is performed with the singularity + // inside the given cell, then a special quadrature formula is + // used that allows one to integrate arbitrary functions against a + // singular weight on the reference cell. // - // \f[ - // \text{dst}_{ik1} := \int_{\text{cell}} \frac{\partial - // G}{\partial \textbf{n}} (y - \text[q]_k) \phi_i(y) dy - // \f] - void compute_SD_integral_on_cell(vector > > &dst, - typename DoFHandler::active_cell_iterator &cell, - const vector > &q, - const Function &rhs); - - // The following two functions are the actual calculations of the - // single and double layer potential kernels, with a minus sign in - // front of them. They are well defined only if the vector $R = - // x-y$ is different from zero. - double nS(const Point &R); - Point nD(const Point &R); - -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) { - AssertThrow(false, ExcImpossibleInDim()); - return 0; - }; - - double term_D(const Point<3> &r, - const Point<3> &a1, - const Point<3> &a2) { - AssertThrow(false, ExcImpossibleInDim()); - return 0; - }; - - SmartPointer > fe; - SmartPointer > fe_values; + // 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. + // + // 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. + double single_layer(const Point &R, + bool factor_out_singularity = false); + Point double_layer(const Point &R, + bool factor_out_singularity = false); }; template -BEMProblem::BEMProblem(const unsigned int degree) : - fe(degree), +BEMProblem::BEMProblem() : + fe(1), dh(tria), external_fe(1), external_dh(external_tria), wind(dim) {} -template -BEMProblem::~BEMProblem() { - LAPACKFullMatrix * p = system_matrix; - system_matrix = 0; - delete p; -} - - template void BEMProblem::read_parameters(std::string filename) { 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.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("Inner quadrature rule"); + prm.enter_subsection("Quadrature rule"); prm.declare_entry("Quadrature type", "gauss", Patterns::Selection(QuadratureSelector<(dim-1)>::get_quadrature_names())); - prm.declare_entry("Quadrature order", "2", Patterns::Integer()); + prm.declare_entry("Quadrature order", "5", Patterns::Integer()); prm.leave_subsection(); prm.enter_subsection("Wind function 2d"); @@ -310,17 +269,11 @@ void BEMProblem::read_parameters(std::string filename) { prm.read_input(filename); n_cycles = prm.get_integer("Number of cycles"); - external_refinement = prm.get_integer("External refinement"); - - prm.enter_subsection("Outer quadrature rule"); - static QuadratureSelector outer_quadrature - (prm.get("Quadrature type"), - prm.get_integer("Quadrature order")); - prm.leave_subsection(); + external_refinement = prm.get_integer("External refinement"); + extend_solution = prm.get_bool("Extend solution on the -2,2 box"); - - prm.enter_subsection("Inner quadrature rule"); - static QuadratureSelector inner_quadrature + prm.enter_subsection("Quadrature rule"); + static QuadratureSelector quadrature (prm.get("Quadrature type"), prm.get_integer("Quadrature order")); prm.leave_subsection(); @@ -330,200 +283,48 @@ void BEMProblem::read_parameters(std::string filename) { wind.parse_parameters(prm); prm.leave_subsection(); - outer_quadrature_pointer = &outer_quadrature; - inner_quadrature_pointer = &inner_quadrature; + quadrature_pointer = &quadrature; } - template -double LaplaceKernelIntegration::nS(const Point &R) { - if(dim == 2) - return (-std::log(R.norm()) / numbers::PI); - else if(dim == 3) - return (1./(R.norm()*numbers::PI) ); - else { - Assert(false, ExcInternalError()); - } +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.; + } return 0.; } template -Point LaplaceKernelIntegration::nD(const Point &R) { +Point LaplaceKernel::double_layer(const Point &R, + bool factor_out_singularity) { Point D(R); - if(dim == 2) - D /= -numbers::PI * R.square(); - else if(dim == 3) - D /= -2*numbers::PI * R.square() * R.norm(); - else { + switch(dim) { + case 2: + factor_out_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; + default: Assert(false, ExcInternalError()); + break; } return D; } - - -template <> -LaplaceKernelIntegration<3>::LaplaceKernelIntegration(const FiniteElement<2,3> &fe) : - fe(&fe) -{ - // 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); - qps[2] = Point<2>(1,0); - qps[3] = Point<2>(1,1); - qps[4] = Point<2>(.5,.5); - vector ws(5,1.); - static Quadrature<2> quadrature(qps, ws); - fe_values = new FEValues<2,3>(fe,quadrature, - update_values | - update_jacobians | - update_cell_normal_vectors | - update_quadrature_points ); -} - - -// The one dimensional singular integration can be calculated -// exploiting QGaussLogR quadrature formula. The quadrature formula -// is constructed in each step, so the constructor is empty. -template <> -LaplaceKernelIntegration<2>::LaplaceKernelIntegration(const FiniteElement<1,2> &fe) : - fe(&fe) -{} - -template -LaplaceKernelIntegration::~LaplaceKernelIntegration() { - // 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. - if(fe_values) { - FEValues *fp = fe_values; - fe_values = 0; - delete fp; - } -} - - -template <> -double -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; - - cross_product(ra1,r,a1); - cross_product(ra2,r,a2); - cross_product(a12,a1,a2); - - double integral = - -1./2./numbers::PI - *( - - ra1*n/a1.norm() * asinh( r*a1/ra1.norm() ) - + ra2*n/a2.norm() * asinh( r*a2/ra2.norm() ) - + rn_c * atan( ra1*ra2 / (r.norm()* (r*(a12)))) - ); - - return integral; - -} - -template <> -double -LaplaceKernelIntegration<3>::term_D (const Point<3> &r, - const Point<3> &a1, - const Point<3> &a2) -{ - Point<3> ra1, ra2, a12; - - cross_product(ra1,r,a1); - cross_product(ra2,r,a2); - cross_product(a12,a1,a2); - - double integral = -1./2./numbers::PI - *atan( ra1*ra2 / (r.norm()* (r*(a12)))); - - return integral; - -} - -template <> -void -LaplaceKernelIntegration<3>::compute_SD_integral_on_cell(vector > > &dstvv, - DoFHandler<2,3>::active_cell_iterator &cell, - const vector > &q_points, - const Function<3> &rhs) -{ - fe_values->reinit(cell); - const vector > &jacobians = fe_values->get_jacobians(); - const vector > &quad_points = fe_values->get_quadrature_points(); - const vector > &normals = fe_values->get_cell_normal_vectors(); - - static vector > cell_wind - ( (*fe_values).n_quadrature_points, Vector(3) ); - static vector normal_wind(quad_points.size()); - - rhs.vector_value_list(quad_points, cell_wind); - - for(unsigned int q=0; q 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) * normal_wind[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() { @@ -583,14 +384,18 @@ void BEMProblem::refine_and_resize() { deallog << "Levels: " << tria.n_levels() << ", potential dofs: " << ndofs << endl; - if(system_matrix) { - LAPACKFullMatrix * p = system_matrix; - system_matrix = 0; - delete p; - } + // The matrix is a full matrix. Notwithstanding this fact, the + // SparseMatrix class coupled with the SparseDirectUMFPACK solver + // are still faster than Lapack solvers. The drawback is that we + // need to assemble a full SparsityPattern. + system_matrix.clear(); + sparsity.reinit(ndofs, ndofs, ndofs); + for(unsigned int i=0; i(ndofs, ndofs); - system_rhs.reinit(ndofs); phi.reinit(ndofs); } @@ -603,55 +408,63 @@ void BEMProblem::assemble_system() { cellj = dh.begin_active(), endc = dh.end(); - // 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 &outer_quadrature = *outer_quadrature_pointer; - Quadrature &inner_quadrature = *inner_quadrature_pointer; - - FEValues fe_outer(fe, outer_quadrature, - update_values | - update_cell_normal_vectors | - update_quadrature_points | - update_JxW_values); - - FEValues fe_inner(fe, inner_quadrature, - update_values | - update_cell_normal_vectors | - update_quadrature_points | - update_JxW_values); + // Quadrature formula for the integration of the kernel in non + // singular cells. This quadrature is selected with the parameter + // file, and should be quite precise, since the functions we are + // integrating are not polynomial functions. + Quadrature &quadrature = *quadrature_pointer; + + // We create initially the singular quadratures for the + // threedimensional problem, since in this case it is 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. + vector > sing_quadratures_3d; + for(unsigned int i=0; i<4; ++i) + sing_quadratures_3d.push_back(QGaussOneOverR<2>(quadrature.size(), i)); + + + FEValues fe_v(fe, quadrature, + update_values | + update_cell_normal_vectors | + update_quadrature_points | + update_JxW_values); - const unsigned int n_q_points_outer = fe_outer.n_quadrature_points; - const unsigned int n_q_points_inner = fe_inner.n_quadrature_points; + 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 > inner_cell_wind(n_q_points_inner, Vector(dim) ); - double inner_normal_wind; + 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); // The kernel. - LaplaceKernelIntegration kernel(fe); - - vector > > single_double_layer_potentials - (fe.dofs_per_cell, vector > - (n_q_points_outer, vector (2, 0.) ) ); + LaplaceKernel kernel; Point R; + // The index i runs on the collocation points, which are the + // 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 = + // GeometryInfo::vertices_per_cell. + AssertThrow(fe.dofs_per_cell == GeometryInfo::vertices_per_cell, + ExcDimensionMismatch(fe.dofs_per_cell, GeometryInfo::vertices_per_cell)); - // The index i runs on outer integration, while j runs on inner integration. for(; celli != endc; ++celli) { - fe_outer.reinit(celli); - const vector > &q_points_outer = fe_outer.get_quadrature_points(); - const vector > &normals_outer = fe_outer.get_cell_normal_vectors(); - + // 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 @@ -665,41 +478,36 @@ void BEMProblem::assemble_system() { local_rhs = 0; local_matrix = 0; - fe_inner.reinit(cellj); + fe_v.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, inner_cell_wind); + 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_inner=0; q_innervertex(i)-q_points[q]; + + local_rhs(i) += ( kernel.single_layer(R) * + normal_wind * + fe_v.JxW(q) ); - local_matrix(i,j) += ( fe_outer.shape_value(i,q_outer) * - fe_outer.JxW(q_outer) * - // - ( kernel.nD(R) * - normals_inner[q_inner] ) * - fe_inner.shape_value(j,q_inner) * - fe_inner.JxW(q_inner) ); - } + for(unsigned int j=0; j::assemble_system() { // same. In this case both the single and the double // layer potential are singular, and they require a // special treatment, as explained in the - // introduction. - if(dim == 3) { - kernel.compute_SD_integral_on_cell(single_double_layer_potentials, - cellj, q_points_outer, 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. + 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]); } - } else { - // In the two dimensional case we only need a - // QGaussLogR quadrature formula to correctly - // integrate the single layer potential. - for(unsigned int q_outer=0; q_outer singular_quad(inner_quadrature.size(), - outer_quadrature.point(q_outer), - 1./cellj->measure()); - FEValues<1,2> fe_v_singular(fe, singular_quad, - update_jacobians | - update_cell_normal_vectors | - update_quadrature_points ); - fe_v_singular.reinit(cellj); - - static vector > singular_cell_wind(singular_quad.size(), - Vector(dim) ); + + 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(); + 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); + 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 ones(dh.n_dofs()), alpha(dh.n_dofs()); + for(unsigned int i=0; i void BEMProblem::solve_system() { - phi.swap(system_rhs); - system_matrix->compute_lu_factorization(); - system_matrix->apply_lu_factorization(phi, false); + SparseDirectUMFPACK LU; + LU.initialize(system_matrix); + 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); } @@ -817,7 +661,7 @@ void BEMProblem::interpolate() { endc = dh.end(); - Quadrature &quadrature = *inner_quadrature_pointer; + Quadrature &quadrature = *quadrature_pointer; FEValues fe_v(fe, quadrature, update_values | @@ -833,7 +677,7 @@ void BEMProblem::interpolate() { vector > local_wind(n_q_points, Vector(dim) ); double normal_wind; - LaplaceKernelIntegration kernel(fe); + LaplaceKernel kernel; Point R; @@ -873,12 +717,12 @@ void BEMProblem::interpolate() { R = external_cell->vertex(i) - q_points[q]; - external_phi(external_dofs[i]) += ( ( - kernel.nS(R) * + external_phi(external_dofs[i]) += ( ( kernel.single_layer(R) * normal_wind - // - ( kernel.nD(R) * - normals[q] ) * - local_phi[q] ) * + ( kernel.double_layer(R) * + normals[q] ) * + local_phi[q] ) * fe_v.JxW(q) ); } } @@ -926,8 +770,8 @@ void BEMProblem::run() { solve_system(); output_results(cycle); } - - interpolate(); + if(extend_solution == true) + interpolate(); } @@ -937,10 +781,10 @@ int main () { deallog.depth_console (3); BEMProblem<2> laplace_problem_2d; - // BEMProblem<3> laplace_problem_3d; - laplace_problem_2d.run(); - // laplace_problem_3d.run(); + + BEMProblem<3> laplace_problem_3d; + laplace_problem_3d.run(); } catch (std::exception &exc) {