]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Modified step-34 to implement exclusively the collocation boundary element
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 23 Mar 2009 18:00:48 +0000 (18:00 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 23 Mar 2009 18:00:48 +0000 (18:00 +0000)
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

deal.II/examples/step-34/doc/intro.dox
deal.II/examples/step-34/parameters.prm
deal.II/examples/step-34/step-34.cc

index 98881f9775a229507d54fb8225ec2c5c07a10ff2..1d559028a036d7f688367ee5ff2f6cace1719273 100644 (file)
@@ -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 &mdash; 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.
 
+<h3> Collocation boundary element method </h3>
 
-<h3> Galerkin boundary element method </h3>
+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$.
 <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
+$\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
index 9bd40dc8374cbc611f8440cbfead36825b8a6be8..b61ad552754cacbb812a533f8b6efb9efffbd3a3 100644 (file)
@@ -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  = 
index 5fdf13e53f0efa53266f3546a597e6dc8b6071a4..0b7b0924b92e12e89bac05773264485664fe6c16 100644 (file)
@@ -45,6 +45,7 @@
 #include <lac/vector.h>
 #include <lac/sparse_direct.h>
 #include <lac/lapack_full_matrix.h>
+#include <lac/matrix_lib.h>
 
 #include <numerics/data_out.h>
 
@@ -59,15 +60,14 @@ using namespace dealii;
 
 
 template <int dim>
-class LaplaceKernelIntegration;
+class LaplaceKernel;
 
 
 template <int dim> 
 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<dim-1, dim>  tria;
-    FE_DGP<dim-1,dim>          fe;
+    FE_Q<dim-1,dim>            fe;
     DoFHandler<dim-1,dim>      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<LAPACKFullMatrix<double> >    system_matrix;    
-    Vector<double>                             system_rhs;
-    Vector<double>                             phi;
+    SparsityPattern            sparsity;
+    SparseMatrix<double>       system_matrix;    
+    Vector<double>             system_rhs;
+    Vector<double>             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<double>     external_phi;
     
     // The following variables are the one that we fill through a
-    // parameter file.
-    // The new objects that we use in this example 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<dim> wind;
-    SmartPointer<Quadrature<dim-1> > outer_quadrature_pointer;
-    SmartPointer<Quadrature<dim-1> > inner_quadrature_pointer;
+    SmartPointer<Quadrature<dim-1> > quadrature_pointer;
     unsigned int n_cycles;
     unsigned int external_refinement;
+    bool extend_solution;
 };
 
 
 
 template <int dim>
-class LaplaceKernelIntegration
+class LaplaceKernel
 {
 public:
-
-    LaplaceKernelIntegration(const FiniteElement<dim-1,dim> &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<vector<vector<double> > > &dst,
-                                    typename DoFHandler<dim-1,dim>::active_cell_iterator &cell,
-                                    const vector<Point<dim> > &q,
-                                    const Function<dim> &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<dim> &R);
-    Point<dim> nD(const Point<dim> &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<const FiniteElement<dim-1, dim> > fe;
-    SmartPointer<FEValues<dim-1,dim> > 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<dim> &R, 
+                       bool factor_out_singularity = false);
+    Point<dim> double_layer(const Point<dim> &R, 
+                           bool factor_out_singularity = false);
 };
 
 
 
 template <int dim>
-BEMProblem<dim>::BEMProblem(const unsigned int degree) :
-    fe(degree),
+BEMProblem<dim>::BEMProblem() :
+    fe(1),
     dh(tria),
     external_fe(1),
     external_dh(external_tria),
     wind(dim)
 {}
 
-template <int dim>
-BEMProblem<dim>::~BEMProblem() {
-    LAPACKFullMatrix<double> * p = system_matrix;
-    system_matrix = 0;
-    delete p;
-}
-
-
 template <int dim> 
 void BEMProblem<dim>::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<dim>::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<dim-1> 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<dim-1> inner_quadrature
+    prm.enter_subsection("Quadrature rule");
+    static QuadratureSelector<dim-1> quadrature
                      (prm.get("Quadrature type"),
                       prm.get_integer("Quadrature order"));
     prm.leave_subsection();
@@ -330,200 +283,48 @@ void BEMProblem<dim>::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 <int dim>
-double LaplaceKernelIntegration<dim>::nS(const Point<dim> &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<dim>::single_layer(const Point<dim> &R, 
+                                       bool factor_out_singularity) {
+    if(factor_out_singularity == true) 
+       return (dim == 2 ? -1. : 1.)/(2*(dim-1)*numbers::PI);
+    else
+       if(dim == 2)
+           return (-std::log(R.norm()) / (2*numbers::PI) );
+       else if(dim == 3)
+           return (1./( R.norm()*4*numbers::PI ) );
+       else {
+           Assert(false, ExcInternalError());
+           return 0.;
+       }
     return 0.;
 }
        
 
 
 template <int dim>
-Point<dim> LaplaceKernelIntegration<dim>::nD(const Point<dim> &R) {
+Point<dim> LaplaceKernel<dim>::double_layer(const Point<dim> &R,
+                                           bool factor_out_singularity) {
     Point<dim> 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<Point<2> > 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<double> 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 <int dim>
-LaplaceKernelIntegration<dim>::~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<dim-1,dim> *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<vector<vector<double> > > &dstvv,
-                                                        DoFHandler<2,3>::active_cell_iterator &cell,
-                                                        const vector<Point<3> > &q_points, 
-                                                        const Function<3> &rhs)
-{
-    fe_values->reinit(cell);
-    const vector<Tensor<2,3> > &jacobians = fe_values->get_jacobians();
-    const vector<Point<3> > &quad_points = fe_values->get_quadrature_points();
-    const vector<Point<3> > &normals = fe_values->get_cell_normal_vectors();
-    
-    static vector<Vector<double> > cell_wind
-       ( (*fe_values).n_quadrature_points, Vector<double>(3) );
-    static vector<double> normal_wind(quad_points.size());
-    
-    rhs.vector_value_list(quad_points, cell_wind);
-    
-    for(unsigned int q=0; q<quad_points.size(); ++q) {
-       normal_wind[q] = 0;
-       for(unsigned int d=0; d<3; ++d)
-           normal_wind[q] += normals[q][d] * cell_wind[q](d);
-    }
-    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) * 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 <int dim>
 void BEMProblem<dim>::read_domain() {
     
@@ -583,14 +384,18 @@ void BEMProblem<dim>::refine_and_resize() {
     deallog << "Levels: " << tria.n_levels()
            << ", potential dofs: " << ndofs <<  endl;
     
-    if(system_matrix) {
-       LAPACKFullMatrix<double> * 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;++i)
+       for(unsigned int j=0; j<ndofs; ++j)
+           sparsity.add(i,j);
+    sparsity.compress();
+    system_matrix.reinit(sparsity);
     
-    system_matrix = new LAPACKFullMatrix<double>(ndofs, ndofs);
-
     system_rhs.reinit(ndofs);
     phi.reinit(ndofs);
 }    
@@ -603,55 +408,63 @@ void BEMProblem<dim>::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<dim-1> &outer_quadrature = *outer_quadrature_pointer;
-    Quadrature<dim-1> &inner_quadrature = *inner_quadrature_pointer;
-    
-    FEValues<dim-1,dim> fe_outer(fe, outer_quadrature,
-                                update_values |
-                                update_cell_normal_vectors |
-                                update_quadrature_points |
-                                update_JxW_values);
-    
-    FEValues<dim-1,dim> 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<dim-1> &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<QGaussOneOverR<2> > sing_quadratures_3d; 
+    for(unsigned int i=0; i<4; ++i)
+       sing_quadratures_3d.push_back(QGaussOneOverR<2>(quadrature.size(), i));
+       
+
+    FEValues<dim-1,dim> 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<unsigned int> dofs_i(fe.dofs_per_cell);
     vector<unsigned int> dofs_j(fe.dofs_per_cell);
 
-    vector<Vector<double> > inner_cell_wind(n_q_points_inner, Vector<double>(dim) );
-    double inner_normal_wind;
+    vector<Vector<double> > cell_wind(n_q_points, Vector<double>(dim) );
+    double normal_wind;
     
     Vector<double>     local_rhs(fe.dofs_per_cell);
     FullMatrix<double>  local_matrix(fe.dofs_per_cell, fe.dofs_per_cell);
     
     // The kernel.
-    LaplaceKernelIntegration<dim> kernel(fe);
-    
-    vector<vector<vector<double> > > single_double_layer_potentials
-       (fe.dofs_per_cell, vector<vector<double> >
-        (n_q_points_outer, vector<double> (2, 0.) ) ); 
+    LaplaceKernel<dim> kernel;
     
     Point<dim> 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<dim>::vertices_per_cell.
+    AssertThrow(fe.dofs_per_cell == GeometryInfo<dim-1>::vertices_per_cell,
+               ExcDimensionMismatch(fe.dofs_per_cell, GeometryInfo<dim-1>::vertices_per_cell));
     
-    // The index i runs on outer integration, while j runs on inner integration.
     for(; celli != endc; ++celli) {
-       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();
-
+       // 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<dim>::assemble_system() {
            local_rhs = 0;
            local_matrix = 0;
            
-           fe_inner.reinit(cellj);
+           fe_v.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, inner_cell_wind);
+           const vector<Point<dim> > &q_points = fe_v.get_quadrature_points();
+           const vector<Point<dim> > &normals = fe_v.get_cell_normal_vectors();
+           wind.vector_value_list(q_pointscell_wind);
            
            if(is_singular == false) {
-               for(unsigned int q_inner=0; q_inner<n_q_points_inner; ++q_inner) {
-                   inner_normal_wind = 0;
+               for(unsigned int q=0; q<n_q_points; ++q) {
+                   normal_wind = 0;
                    for(unsigned int d=0; d<dim; ++d) 
-                       inner_normal_wind += normals_inner[q_inner][d]*inner_cell_wind[q_inner](d);
-                   
-                   for(unsigned int q_outer=0; q_outer<n_q_points_outer; ++q_outer) {
-                   
-                       R = q_points_outer[q_outer]-q_points_inner[q_inner];
+                       normal_wind += normals[q][d]*cell_wind[q](d);
                        
-                       for(unsigned int i=0; i<fe.dofs_per_cell; ++i) {
-                           local_rhs(i) += ( fe_outer.shape_value(i,q_outer)   *
-                                             fe_outer.JxW(q_outer)             *
-                                             //
-                                             kernel.nS(R)                      * 
-                                             inner_normal_wind                 *
-                                             fe_inner.JxW(q_inner) );
-                               
-                           for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
+                   for(unsigned int i=0; i<fe.dofs_per_cell; ++i) {
+                       
+                       // Distance between the external support point
+                       // and the quadrature point on the internal
+                       // cell.
+                       R = celli->vertex(i)-q_points[q];
+                           
+                       local_rhs(i) += ( kernel.single_layer(R)        * 
+                                         normal_wind                   *
+                                         fe_v.JxW(q) );
                                
-                               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<fe.dofs_per_cell; ++j) {
+                           
+                           local_matrix(i,j) += ( ( kernel.double_layer(R)     * 
+                                                    normals[q] )               *
+                                                  fe_v.shape_value(j,q)        *
+                                                  fe_v.JxW(q)  );
                        }
                    }
                }
@@ -710,94 +518,130 @@ void BEMProblem<dim>::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<fe.dofs_per_cell; ++i) {
-                       for(unsigned int q_outer=0; q_outer<n_q_points_outer; ++q_outer) {
-                           local_rhs(i) += ( - single_double_layer_potentials[0][q_outer][0]  * 
-                                             fe_outer.shape_value(i,q_outer)              *
-                                             fe_outer.JxW(q_outer) );
-                           
-                           for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
-                               
-                               // When the indices are the same, we
-                               // assemble also the mass matrix.
-                               local_matrix(i,j) += ( fe_outer.shape_value(i,q_outer)       * 
-                                                      fe_outer.shape_value(j,q_outer)       *
-                                                      fe_outer.JxW(q_outer) );
-                               
-                               local_matrix(i,j) += ( -single_double_layer_potentials[j][q_outer][1] * 
-                                                      fe_outer.shape_value(i,q_outer)               *
-                                                      fe_outer.JxW(q_outer) );
-                           }
-                       }
+               // introduction.
+               // 
+               // In the two dimensional case we perform the
+               // integration using a QGaussLogR quadrature formula,
+               // which is specifically designed to integrate
+               // logarithmic singularities on the unit interval,
+               // while in three dimensions we use the
+               // QGaussOneOverR, which allows us to integrate 1/R
+               // singularities on the vertices of the reference
+               // element. Since we don't want to rebuild the two
+               // dimensional quadrature formula at each singular
+               // integration, we built them outside the loop on the
+               // cells, and we only use a pointer to that quadrature
+               // here.
+               //
+               // Notice that in one dimensional integration this is
+               // not possible, since we need to know the scaling
+               // parameter for the quadrature, which is not known a
+               // priori.
+               //
+               // Dimension independent programming here is a little
+               // tricky, but can be achieved via dynamic casting. We
+               // check that everything went ok with an assertion at
+               // the end of this block. Notice that the dynamic cast
+               // will only work when the dimension is the correct
+               // one, in which case it is possible to cast a
+               // QGaussLogR and QGaussOneOverR to a Quadrature<1>
+               // and Quadrature<2> object.
+               //
+               // In the other cases this won't be called, and even
+               // if it was, the dynamic_cast function would just
+               // return a null pointer. We check that this is not
+               // the case with the Assert at the end.
+               //
+               // Notice that in two dimensions the singular
+               // quadrature rule depends also on the size of the
+               // current cell. For this reason, it is necessary to
+               // create a new quadrature for each singular
+               // integration. Since we create it using the new
+               // operator of C++, we also need to destroy it using
+               // the dual of new: delete. This is done at the end,
+               // and only if dim == 2.
+               Quadrature<dim-1> * singular_quadrature;
+               for(unsigned int i=0; i<fe.dofs_per_cell; ++i) {
+                   if(dim == 2) {
+                       singular_quadrature = dynamic_cast<Quadrature<dim-1> *>(
+                           new QGaussLogR<1>(quadrature.size(),
+                                             Point<1>((double)i),
+                                             1./cellj->measure()));
+                   } else {
+                       singular_quadrature = dynamic_cast<Quadrature<dim-1> *>(
+                           & sing_quadratures_3d[i]);
                    }
-               } 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<n_q_points_outer; ++q_outer) {
-                       QGaussLogR<1> 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<Vector<double> > singular_cell_wind(singular_quad.size(), 
-                                                                         Vector<double>(dim) );
+                   
+                   Assert(singular_quadrature, ExcInternalError());
+
+                   
+                   FEValues<dim-1,dim> fe_v_singular(fe, *singular_quadrature, 
+                                                     update_jacobians |
+                                                     update_values |
+                                                     update_cell_normal_vectors |
+                                                     update_quadrature_points );
+                   fe_v_singular.reinit(cellj);
+                   
+                   static vector<Vector<double> > singular_cell_wind( (*singular_quadrature).size(), 
+                                                                      Vector<double>(dim) );
                        
-                       const vector<Point<dim> > &singular_normals = fe_v_singular.get_cell_normal_vectors();
-                       const vector<Point<dim> > &singular_q_points = fe_v_singular.get_quadrature_points();
+                   const vector<Point<dim> > &singular_normals = fe_v_singular.get_cell_normal_vectors();
+                   const vector<Point<dim> > &singular_q_points = fe_v_singular.get_quadrature_points();
                        
-                       wind.vector_value_list(singular_q_points, singular_cell_wind);
+                   wind.vector_value_list(singular_q_points, singular_cell_wind);
+                   
+                   for(unsigned int q=0; q<singular_quadrature->size(); ++q) {
+                       R = celli->vertex(i)-singular_q_points[q];
+                       double normal_wind = 0;
+                       for(unsigned int d=0; d<dim; ++d)
+                           normal_wind += (singular_cell_wind[q](d)*
+                                           singular_normals[q][d]);
                        
-                       for(unsigned int i=0; i<fe.dofs_per_cell; ++i) {
-                           for(unsigned int q_inner=0; q_inner<singular_quad.size(); ++q_inner) {
-                               double normal_wind = 0;
-                               for(unsigned int d=0; d<dim; ++d)
-                                   normal_wind += (singular_cell_wind[q_inner](d)*
-                                                   singular_normals[q_inner][d]);
-                               
-                               local_rhs(i) -= ( normal_wind *
-                                                 fe_v_singular.JxW(q_inner)            /
-                                                 numbers::PI                           *
-                                                 fe_outer.shape_value(i,q_outer)       *
-                                                 fe_outer.JxW(q_outer) );
-                               
-                               for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
-                                   
-                                   // When the indices are the same, we
-                                   // assemble also the mass matrix.
-                                   local_matrix(i,j) += ( fe_outer.shape_value(i,q_outer)           * 
-                                                          fe_outer.shape_value(j,q_outer)           *
-                                                          fe_outer.JxW(q_outer) );
-                               }
-                           }
+                       local_rhs(i) += ( kernel.single_layer(R, is_singular) *
+                                         normal_wind                         *
+                                         fe_v_singular.JxW(q) );
+
+                       for(unsigned int j=0; j<fe.dofs_per_cell; ++j) {
+                           local_matrix(i,j) += (( kernel.double_layer(R, is_singular)  *
+                                                   singular_normals[q])                 *
+                                                 fe_v_singular.shape_value(j,q)         *
+                                                 fe_v_singular.JxW(q)  );
                        }
                    }
+                   if(dim==2) delete singular_quadrature;
                }
            }
            // Move the local matrix and rhs to the global one.
            for(unsigned int i=0; i<fe.dofs_per_cell; ++i) {
                system_rhs(dofs_i[i]) += local_rhs(i);
                for(unsigned int j=0; j<fe.dofs_per_cell; ++j) 
-                   (*system_matrix)(dofs_i[i],dofs_j[j]) += local_matrix(i,j);
+                   system_matrix.add(dofs_i[i],dofs_j[j], local_matrix(i,j));
            }
        }
     }
+    // One quick way to compute the matrix of the solid angles, is to
+    // use the Neumann matrix itself. It is enough to multiply the
+    // matrix with the vector of ones, to get the diagonal matrix of
+    // the alpha solid angles.
+    Vector<double> ones(dh.n_dofs()), alpha(dh.n_dofs());
+    for(unsigned int i=0; i<dh.n_dofs(); ++i) 
+       ones(i) = 1.;
+    system_matrix.vmult(alpha, ones);
+    for(unsigned int i=0; i<dh.n_dofs(); ++i) 
+       system_matrix.add(i,i,-alpha(i));
 }
 
 template <int dim>
 void BEMProblem<dim>::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<dim>::interpolate() {
        endc = dh.end();
 
     
-    Quadrature<dim-1> &quadrature = *inner_quadrature_pointer;
+    Quadrature<dim-1> &quadrature = *quadrature_pointer;
     
     FEValues<dim-1,dim> fe_v(fe, quadrature,
                             update_values |
@@ -833,7 +677,7 @@ void BEMProblem<dim>::interpolate() {
     vector<Vector<double> > local_wind(n_q_points, Vector<double>(dim) );
     double normal_wind;
     
-    LaplaceKernelIntegration<dim> kernel(fe);
+    LaplaceKernel<dim> kernel;
     Point<dim> R;
 
 
@@ -873,12 +717,12 @@ void BEMProblem<dim>::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<dim>::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)
     {

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.