]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implemented the Telles quadrature formulas and their test 708/head
authorNicola Giuliani <ngiuliani@sissa.it>
Mon, 30 Mar 2015 16:22:05 +0000 (18:22 +0200)
committerNicola Giuliani <ngiuliani@sissa.it>
Thu, 2 Apr 2015 07:42:16 +0000 (09:42 +0200)
implemented some suggestions from the authors

applied the comments on the pull request

bug removed

class made more readable

class made more readable

testing over m

changes.h ok

corrected the telles test

adjust the tolerance

corrected the formula

corrected the change.h

added a description for the dim > 1 cases

added the reference for Telles

doc/news/changes.h
include/deal.II/base/quadrature_lib.h
source/base/quadrature_lib.cc
tests/codim_one/integrate_one_over_r_telles.cc [new file with mode: 0644]
tests/codim_one/integrate_one_over_r_telles.output [new file with mode: 0644]

index 1dabada2bb3a299157693cdb09a539f479dddb0c..50c4e72e2cf5d840e598b6fb0f2530b9b52faa4d 100644 (file)
@@ -402,9 +402,16 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
-  <li> New: A function to get a map with all vertices at boundaries
-  has been added at GridTools::get_all_vertices_at_boundary.
-  This function will return a map which can be used in functions
+  <li> New: There is now a new quadrature formula in quadrature_lib. It is
+  now possible to use Telles' quadrature rules through the function QTelles
+  to integrate singular integrals
+  <br>
+  (Nicola Giuliani, 2015/04/01)
+  <li>
+
+  <li> New: A new function to get a map with all vertices at boundaries
+  has been added in GridTools::get_all_vertices_at_boundary.
+  This function will return a map wich can be used in functions
   like GridTools::Laplace_transform.
   <br>
   (Fernando Posada, 2015/03/31)
index b484edb34cac26118478decc5c9ea045b460b753..0400b13eacd9eb3a679407b1ec0fb528c6c7895a 100644 (file)
@@ -439,6 +439,75 @@ public:
                   const std::pair<double, Point<dim> > &b);
 };
 
+/**
+ * Telles quadrature of arbitrary order.
+ *
+ * The coefficients of these quadrature rules are computed using
+ * a non linear change of variables starting from a Gauss-Legendre
+ * quadrature formula.
+ * This is done using a cubic polynomial,
+ * $n = a x^3 + b x^2 + c x + d$
+ * in order to integrate
+ * a singular integral, with singularity at a given point x_0.
+ *
+ * We start from a Gauss Quadrature Formula with arbitrary
+ * function. Then we apply the cubic variable change.
+ * In the paper, J.C.F.Telles:A Self-Adaptive Co-ordinate Transformation
+ * For Efficient Numerical Evaluation of General Boundary Element Integrals.
+ * International Journal for Numerical Methods in Engineering, vol 24,
+ * pages 959–973. year 1987, the author applies the transformation on the
+ * reference cell $[-1, 1]$ getting
+ * @f{align*} n(1) &= 1, \\ n(-1) &= -1, \\ dn/dx $= 0 at x = x_0, \\
+ d2n/dx2 = 0 at x &= x_0 @f}
+ * We get
+ * @f{align*} a &= 1/q, \\ b &= -3gamma_bar/q, \\ c &= 3gamma_bar/q, \\
+ d &= -b, @f}
+ * with
+ * @f{align*} eta_star &= eta_bar^2 - 1, \\ gamma_bar 6 &= nthroot( eta_bar
+ * eta_star + abs(eta_star) ,3) + nthroot(eta_bar*eta_star - abs(eta_star),3)
+ + eta_bar, \\ q &= ((gamma-gamma_bar).^3 + gamma_bar*(gamma_bar^2+3))
+ /(1+3*gamma_bar^2) @f}.
+ * Since the library assumes [0,1] as reference interval, we will map
+ * these values on the proper reference interval in the implementation.
+ *
+ * This variable change can be used to integrate singular integrals.
+ * One example is $\f(x)/abs(x-x_0)$ on the reference interval $[0,1]$,
+ * where $x_0$ is given at construction time, and is the location of the
+ * singularity $x_0$, and $f(x)$ is a smooth non singular function.
+ *
+ * Singular quadrature formula are rather expensive, nevertheless Telles'
+ * quadrature formula are much easier to compute with respect to other singular
+ * integration techniques as Lachat-Watson.
+ *
+ * We have implemented the case for $dim = 1$. When we deal the case $dim >1$
+ * we have computed the quadrature formula has a tensorial product of one
+ * dimensional Telles' quadrature formulas considering the different components
+ * of the singularity.
+ *
+ * The weights and functions for Gauss Legendre formula have been tabulated up
+ * to order 12.
+ *
+ * @author Nicola Giuliani, Luca Heltai 2015
+ */
+template <int dim>
+class QTelles: public Quadrature<dim>
+{
+public:
+  /**
+  * A constructor that takes a quadrature formula and a singular point as
+  * argument. The quadrature formula will be mapped using Telles' rule. Make
+  * sure that the order of the quadrature rule is appropriate for the
+  * singularity in question.
+  **/
+  QTelles (const Quadrature<1> &base_quad, const Point<dim> &singularity);
+  /**
+  * A variant of above constructor that takes as parameters the order @p n
+  * and location of a singularity. A Gauss Legendre quadrature of order n
+  * will be used
+  **/
+  QTelles (const unsigned int n, const Point<dim> &singularity);
+
+};
 
 /*@}*/
 
@@ -470,6 +539,7 @@ template <> QWeddle<1>::QWeddle ();
 template <> QGaussLog<1>::QGaussLog (const unsigned int n, const bool revert);
 template <> QGaussLogR<1>::QGaussLogR (const unsigned int n, const Point<1> x0, const double alpha, const bool flag);
 template <> QGaussOneOverR<2>::QGaussOneOverR (const unsigned int n, const unsigned int index, const bool flag);
+template <> QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> &singularity);
 
 
 
index 98304d1e656ebb951c6d1c777008ea59fde024cd..829bb52fb3e58330fe61677013e8f07fafb5dfbc 100644 (file)
@@ -1014,6 +1014,118 @@ QWeddle<dim>::QWeddle ()
   Quadrature<dim> (QWeddle<dim-1>(), QWeddle<1>())
 {}
 
+template <int dim>
+QTelles<dim>::QTelles (
+  const Quadrature<1> &base_quad, const Point<dim> &singularity)
+  :
+/**
+* We need the explicit implementation if dim == 1. If dim > 1 we use the
+* former implementation and apply a tensorial product to obtain the higher
+* dimensions.
+**/
+  Quadrature<dim>(
+    dim == 2 ?
+    QAnisotropic<dim>(
+      QTelles<1>(base_quad, Point<1>(singularity[0])),
+      QTelles<1>(base_quad, Point<1>(singularity[1]))) :
+    dim == 3 ?
+    QAnisotropic<dim>(
+      QTelles<1>(base_quad, Point<1>(singularity[0])),
+      QTelles<1>(base_quad, Point<1>(singularity[1])),
+      QTelles<1>(base_quad, Point<1>(singularity[2]))) :
+    Quadrature<dim>())
+{
+}
+
+template <int dim>
+QTelles<dim>::QTelles (
+  const unsigned int n, const Point<dim> &singularity)
+  :
+/**
+* In this case we map the standard Gauss Legendre formula using the given
+* singularity point coordinates.
+**/
+  Quadrature<dim>(QTelles<dim>(QGauss<1>(n), singularity))
+{}
+
+
+
+template <>
+QTelles<1>::QTelles (
+  const Quadrature<1> &base_quad, const Point<1> &singularity)
+  :
+/**
+* We explicitly implement the Telles' variable change if dim == 1.
+**/
+  Quadrature<1>(base_quad)
+{
+  /**
+  * We define all the constants to be used in the implementation of
+  * Telles' rule
+  **/
+  const double eta_bar = singularity[0] * 2. - 1.;
+  const double eta_star = eta_bar * eta_bar - 1.;
+  double gamma_bar;
+
+  std::vector<Point<1>> quadrature_points_dummy(quadrature_points.size());
+  std::vector<double> weights_dummy(weights.size());
+  unsigned int cont = 0;
+  const double tol = 1e-10;
+  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
+    {
+      if (std::abs(quadrature_points[d][0] - singularity[0]) > tol)
+        {
+          quadrature_points_dummy[d-cont] = quadrature_points[d];
+          weights_dummy[d-cont] = weights[d];
+        }
+      else
+        {
+          // We need to remove the singularity point from the quadrature point
+          // list. To do so we use the variable cont.
+          cont = 1;
+        }
+
+    }
+  if (cont == 1)
+    {
+      quadrature_points.resize(quadrature_points_dummy.size()-1);
+      weights.resize(weights_dummy.size()-1);
+      for (unsigned int d = 0; d < quadrature_points.size()-1; ++d)
+        {
+          quadrature_points[d] = quadrature_points_dummy[d];
+          weights[d] = weights_dummy[d];
+        }
+    }
+  // We need to check if the singularity is at the boundary of the interval.
+  if (std::abs(eta_star) <= tol)
+    {
+      gamma_bar = std::pow((eta_bar * eta_star + std::abs(eta_star)),1.0 / 3.0)
+                  + std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0)
+                  + eta_bar;
+    }
+  else
+    {
+      gamma_bar = (eta_bar * eta_star + std::abs(eta_star))/std::abs(eta_bar * eta_star + std::abs(eta_star))*
+                  std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)),1.0 / 3.0)
+                  + (eta_bar * eta_star - std::abs(eta_star))/std::abs(eta_bar * eta_star - std::abs(eta_star))*
+                  std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0)
+                  + eta_bar;
+    }
+  for (unsigned int q = 0; q < quadrature_points.size(); ++q)
+    {
+      double gamma = quadrature_points[q][0] * 2 - 1;
+      double eta = (std::pow(gamma - gamma_bar, 3.0)
+                    + gamma_bar * (gamma_bar * gamma_bar + 3))
+                   / (1 + 3 * gamma_bar * gamma_bar);
+
+      double J = 3 * ((gamma - gamma_bar) *(gamma - gamma_bar))
+                 / (1 + 3 * gamma_bar * gamma_bar);
+
+      quadrature_points[q][0] = (eta + 1) / 2.0;
+      weights[q] = J * weights[q];
+
+    }
+}
 
 // explicit specialization
 // note that 1d formulae are specialized by implementation above
@@ -1037,4 +1149,8 @@ template class QSorted<1>;
 template class QSorted<2>;
 template class QSorted<3>;
 
+template class QTelles<1> ;
+template class QTelles<2> ;
+template class QTelles<3> ;
+
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/codim_one/integrate_one_over_r_telles.cc b/tests/codim_one/integrate_one_over_r_telles.cc
new file mode 100644 (file)
index 0000000..46d1613
--- /dev/null
@@ -0,0 +1,302 @@
+// ---------------------------------------------------------------------
+// $Id: integrate_one_over_r.cc 30338 2013-08-18 22:02:27Z heltai $
+//
+// Copyright (C) 2005 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// integrates the function *f(x,y)/R, where f(x,y) is a power of x and
+// y on the set [0,1]x[0,1]. dim = 2 only.
+
+#include <fstream>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+
+// all include files needed for the program
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <fstream>
+#include <string>
+
+#include <math.h>
+
+using namespace std;
+using namespace dealii;
+
+// We test the integration of singular kernels with a singularity of kind 1/R
+// We multiply this function with a polynomial up to degree 6.
+
+double
+exact_integral_one_over_r (
+    const unsigned int i, const unsigned int j,
+    const unsigned int vertex_index);
+
+ofstream logfile("output");
+
+int
+main ()
+{
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog << std::fixed;
+
+  deallog << endl
+          << "Calculation of the integral of f(x,y)*1/R on [0,1]x[0,1]" << endl
+          << "for f(x,y) = x^i y^j, with i,j ranging from 0 to 5, and R being" << endl
+          << "the distance from (x,y) to four vertices of the square." << endl
+          << endl;
+
+
+  std::vector<Point<2> > vertices = FE_Q<2>(1).get_unit_support_points();
+
+  for (unsigned int m=1; m<7; ++m)
+    {
+      deallog << " =========Quadrature Order: " << m
+              << " =============================== " << endl;
+      deallog << " ============================================================ " << endl;
+      unsigned int index=0;
+        {
+          deallog << " ===============Vertex: " << vertices[index]
+                  << " ============================= " << endl;
+          QTelles<2> quad(m, vertices[index]);
+          QGaussOneOverR<2> quad2(m, vertices[index]);
+
+
+          for (unsigned int i = 0; i < 6; ++i)
+            {
+              for (unsigned int j = 0; j < 6; ++j)
+                {
+
+                  double exact_integral  = exact_integral_one_over_r(index, i,j);
+                  double approx_integral = 0;
+                  double approx_integral_2 = 0;
+
+                  for (unsigned int q=0; q< quad.size(); ++q)
+                    {
+                      double x = quad.point(q)[0];
+                      double y = quad.point(q)[1];
+                      double R = sqrt(x*x+y*y);
+                      approx_integral += ( pow(x, (double)i) *
+                                           pow(y, (double)j) / R *
+                                           quad.weight(q) );
+                    }
+
+                  for (unsigned int q=0; q< quad2.size(); ++q)
+                    {
+                      double x = quad2.point(q)[0];
+                      double y = quad2.point(q)[1];
+                      double R = sqrt(x*x+y*y);
+                      approx_integral_2 += ( pow(x, (double)i) *
+                                             pow(y, (double)j) *
+                                             quad2.weight(q) );
+                    }
+
+                  deallog << "f(x,y) = x^" << i
+                          << " y^" << j
+                          << ", Errors = "
+                          << approx_integral - exact_integral
+                          << ", "
+                          << approx_integral_2 - exact_integral
+                          << std::endl;
+                }
+            }
+        }
+    }
+    }
+
+double exact_integral_one_over_r(const unsigned int vertex_index,
+    const unsigned int i,
+    const unsigned int j)
+  {
+    Assert(vertex_index < 4, ExcInternalError());
+    Assert(i<6, ExcNotImplemented());
+    Assert(j<6, ExcNotImplemented());
+
+// The integrals are computed using the following maple snippet of
+// code:
+//
+//      singint := proc(index, N, M)
+//         if index = 0 then
+//            return int(int(x^N *y^M/sqrt(x^2+y^2), x=0.0..1.0), y=0.0..1.0);
+//         elif index = 1 then
+//            return int(int(x^N *y^M/sqrt((x-1)^2+y^2), x=0.0..1.0), y=0.0..1.0);
+//         elif index = 2 then
+//            return int(int(x^N *y^M/sqrt(x^2+(y-1)^2), x=0.0..1.0), y=0.0..1.0);
+//         elif index = 3 then
+//            return int(int((1-x)^N *(1-y)^M/sqrt(x^2+y^2), x=0.0..1.0), y=0.0..1.0);
+//         end if;
+//      end proc;
+//      Digits := 20;
+//      for i from 3 to 3 do
+//         for n from 0 to 5 do
+//          for m from 0 to 5 do
+//               v[i+1][n+1][m+1] = sing_int(i, n, m);
+//            end do;
+//         end do;
+//      end do;
+//      C(v)
+
+    static double v[4][6][6] =
+      {
+          {
+              { 0}}};
+    if (v[0][0][0] == 0)
+      {
+        v[0][0][0] = 0.17627471740390860505e1;
+        v[0][0][1] = 0.64779357469631903702e0;
+        v[0][0][2] = 0.38259785823210634567e0;
+        v[0][0][3] = 0.26915893322379450224e0;
+        v[0][0][4] = 0.20702239737104695572e0;
+        v[0][0][5] = 0.16800109713227567467e0;
+        v[0][1][0] = 0.64779357469631903702e0;
+        v[0][1][1] = 0.27614237491539669920e0;
+        v[0][1][2] = 0.17015838751246776515e0;
+        v[0][1][3] = 0.12189514164974600651e0;
+        v[0][1][4] = 0.94658660368131133694e-1;
+        v[0][1][5] = 0.77263794021029438797e-1;
+        v[0][2][0] = 0.38259785823210634567e0;
+        v[0][2][1] = 0.17015838751246776515e0;
+        v[0][2][2] = 0.10656799507071040471e0;
+        v[0][2][3] = 0.76947022258735165920e-1;
+        v[0][2][4] = 0.60022626787495395021e-1;
+        v[0][2][5] = 0.49131622931360879320e-1;
+        v[0][3][0] = 0.26915893322379450224e0;
+        v[0][3][1] = 0.12189514164974600651e0;
+        v[0][3][2] = 0.76947022258735165919e-1;
+        v[0][3][3] = 0.55789184535895709637e-1;
+        v[0][3][4] = 0.43625068213915842136e-1;
+        v[0][3][5] = 0.35766126849971778500e-1;
+        v[0][4][0] = 0.20702239737104695572e0;
+        v[0][4][1] = 0.94658660368131133694e-1;
+        v[0][4][2] = 0.60022626787495395021e-1;
+        v[0][4][3] = 0.43625068213915842137e-1;
+        v[0][4][4] = 0.34164088852375945192e-1;
+        v[0][4][5] = 0.28037139560980277614e-1;
+        v[0][5][0] = 0.16800109713227567467e0;
+        v[0][5][1] = 0.77263794021029438797e-1;
+        v[0][5][2] = 0.49131622931360879320e-1;
+        v[0][5][3] = 0.35766126849971778501e-1;
+        v[0][5][4] = 0.28037139560980277614e-1;
+        v[0][5][5] = 0.23024181049838367777e-1;
+        v[1][0][0] = 0.17627471740390860505e1;
+        v[1][0][1] = 0.64779357469631903702e0;
+        v[1][0][2] = 0.38259785823210634567e0;
+        v[1][0][3] = 0.26915893322379450224e0;
+        v[1][0][4] = 0.20702239737104695572e0;
+        v[1][0][5] = 0.16800109713227567467e0;
+        v[1][1][0] = 0.11149535993427670134e1;
+        v[1][1][1] = 0.37165119978092233782e0;
+        v[1][1][2] = 0.21243947071963858053e0;
+        v[1][1][3] = 0.14726379157404849573e0;
+        v[1][1][4] = 0.11236373700291582202e0;
+        v[1][1][5] = 0.90737303111246235871e-1;
+        v[1][2][0] = 0.84975788287855432210e0;
+        v[1][2][1] = 0.26566721237799340376e0;
+        v[1][2][2] = 0.14884907827788122009e0;
+        v[1][2][3] = 0.10231567218303765515e0;
+        v[1][2][4] = 0.77727703422280083352e-1;
+        v[1][2][5] = 0.62605132021577676395e-1;
+        v[1][3][0] = 0.69800109142265347423e0;
+        v[1][3][1] = 0.20794647083778622837e0;
+        v[1][3][2] = 0.11487965864809909847e0;
+        v[1][3][3] = 0.78525390514866270852e-1;
+        v[1][3][4] = 0.59489228415223897572e-1;
+        v[1][3][5] = 0.47838457013298217744e-1;
+        v[1][4][0] = 0.59754668912231692323e0;
+        v[1][4][1] = 0.17125249387868593878e0;
+        v[1][4][2] = 0.93606816359052444729e-1;
+        v[1][4][3] = 0.63728830247554475330e-1;
+        v[1][4][4] = 0.48187332620207367724e-1;
+        v[1][4][5] = 0.38708290797416359020e-1;
+        v[1][5][0] = 0.52527944036356840363e0;
+        v[1][5][1] = 0.14574366656617935708e0;
+        v[1][5][2] = 0.78997159795636003667e-1;
+        v[1][5][3] = 0.53620816423066464705e-1;
+        v[1][5][4] = 0.40487985967086264433e-1;
+        v[1][5][5] = 0.32498604596082509165e-1;
+        v[2][0][0] = 0.17627471740390860505e1;
+        v[2][0][1] = 0.11149535993427670134e1;
+        v[2][0][2] = 0.84975788287855432210e0;
+        v[2][0][3] = 0.69800109142265347419e0;
+        v[2][0][4] = 0.59754668912231692318e0;
+        v[2][0][5] = 0.52527944036356840362e0;
+        v[2][1][0] = 0.64779357469631903702e0;
+        v[2][1][1] = 0.37165119978092233782e0;
+        v[2][1][2] = 0.26566721237799340376e0;
+        v[2][1][3] = 0.20794647083778622835e0;
+        v[2][1][4] = 0.17125249387868593876e0;
+        v[2][1][5] = 0.14574366656617935708e0;
+        v[2][2][0] = 0.38259785823210634567e0;
+        v[2][2][1] = 0.21243947071963858053e0;
+        v[2][2][2] = 0.14884907827788122009e0;
+        v[2][2][3] = 0.11487965864809909845e0;
+        v[2][2][4] = 0.93606816359052444712e-1;
+        v[2][2][5] = 0.78997159795636003667e-1;
+        v[2][3][0] = 0.26915893322379450223e0;
+        v[2][3][1] = 0.14726379157404849572e0;
+        v[2][3][2] = 0.10231567218303765514e0;
+        v[2][3][3] = 0.78525390514866270835e-1;
+        v[2][3][4] = 0.63728830247554475311e-1;
+        v[2][3][5] = 0.53620816423066464702e-1;
+        v[2][4][0] = 0.20702239737104695572e0;
+        v[2][4][1] = 0.11236373700291582202e0;
+        v[2][4][2] = 0.77727703422280083352e-1;
+        v[2][4][3] = 0.59489228415223897563e-1;
+        v[2][4][4] = 0.48187332620207367713e-1;
+        v[2][4][5] = 0.40487985967086264434e-1;
+        v[2][5][0] = 0.16800109713227567468e0;
+        v[2][5][1] = 0.90737303111246235879e-1;
+        v[2][5][2] = 0.62605132021577676399e-1;
+        v[2][5][3] = 0.47838457013298217740e-1;
+        v[2][5][4] = 0.38708290797416359014e-1;
+        v[2][5][5] = 0.32498604596082509169e-1;
+        v[3][0][0] = 0.17627471740390860505e1;
+        v[3][0][1] = 0.11149535993427670134e1;
+        v[3][0][2] = 0.84975788287855432210e0;
+        v[3][0][3] = 0.69800109142265347419e0;
+        v[3][0][4] = 0.59754668912231692318e0;
+        v[3][0][5] = 0.52527944036356840362e0;
+        v[3][1][0] = 0.11149535993427670134e1;
+        v[3][1][1] = 0.74330239956184467563e0;
+        v[3][1][2] = 0.58409067050056091834e0;
+        v[3][1][3] = 0.49005462058486724584e0;
+        v[3][1][4] = 0.42629419524363098443e0;
+        v[3][1][5] = 0.37953577379738904654e0;
+        v[3][2][0] = 0.84975788287855432210e0;
+        v[3][2][1] = 0.58409067050056091834e0;
+        v[3][2][2] = 0.46727253640044873467e0;
+        v[3][2][3] = 0.39698780839518011595e0;
+        v[3][2][4] = 0.34864851772399749038e0;
+        v[3][2][5] = 0.31278926702684569312e0;
+        v[3][3][0] = 0.69800109142265347423e0;
+        v[3][3][1] = 0.49005462058486724586e0;
+        v[3][3][2] = 0.39698780839518011599e0;
+        v[3][3][3] = 0.34027526433872581371e0;
+        v[3][3][4] = 0.30088082631586196583e0;
+        v[3][3][5] = 0.27141910362887187844e0;
+        v[3][4][0] = 0.59754668912231692323e0;
+        v[3][4][1] = 0.42629419524363098445e0;
+        v[3][4][2] = 0.34864851772399749044e0;
+        v[3][4][3] = 0.30088082631586196576e0;
+        v[3][4][4] = 0.26744962339187730308e0;
+        v[3][4][5] = 0.24229245314748740295e0;
+        v[3][5][0] = 0.52527944036356840363e0;
+        v[3][5][1] = 0.37953577379738904655e0;
+        v[3][5][2] = 0.31278926702684569301e0;
+        v[3][5][3] = 0.27141910362887187862e0;
+        v[3][5][4] = 0.24229245314748740263e0;
+        v[3][5][5] = 0.22026586649771582089e0;
+      }
+    return v[vertex_index][i][j];
+  }
diff --git a/tests/codim_one/integrate_one_over_r_telles.output b/tests/codim_one/integrate_one_over_r_telles.output
new file mode 100644 (file)
index 0000000..5afde93
--- /dev/null
@@ -0,0 +1,240 @@
+
+DEAL::
+DEAL::Calculation of the integral of f(x,y)*1/R on [0,1]x[0,1]
+DEAL::for f(x,y) = x^i y^j, with i,j ranging from 0 to 5, and R being
+DEAL::the distance from (x,y) to four vertices of the square.
+DEAL::
+DEAL:: =========Quadrature Order: 1 =============================== 
+DEAL:: ============================================================ 
+DEAL:: ===============Vertex: 0.000000 0.000000 ============================= 
+DEAL::f(x,y) = x^0 y^0, Errors = 1.419233, -0.062529
+DEAL::f(x,y) = x^0 y^1, Errors = -0.250046, -0.046676
+DEAL::f(x,y) = x^0 y^2, Errors = -0.332879, -0.133607
+DEAL::f(x,y) = x^0 y^3, Errors = -0.262944, -0.155343
+DEAL::f(x,y) = x^0 y^4, Errors = -0.206246, -0.152327
+DEAL::f(x,y) = x^0 y^5, Errors = -0.167904, -0.141111
+DEAL::f(x,y) = x^1 y^0, Errors = -0.250046, -0.046676
+DEAL::f(x,y) = x^1 y^1, Errors = -0.226424, -0.100079
+DEAL::f(x,y) = x^1 y^2, Errors = -0.163944, -0.107911
+DEAL::f(x,y) = x^1 y^3, Errors = -0.121118, -0.096111
+DEAL::f(x,y) = x^1 y^4, Errors = -0.094562, -0.082873
+DEAL::f(x,y) = x^1 y^5, Errors = -0.077252, -0.071600
+DEAL::f(x,y) = x^2 y^0, Errors = -0.332879, -0.133607
+DEAL::f(x,y) = x^2 y^1, Errors = -0.163944, -0.107911
+DEAL::f(x,y) = x^2 y^2, Errors = -0.105791, -0.088336
+DEAL::f(x,y) = x^2 y^3, Errors = -0.076850, -0.070501
+DEAL::f(x,y) = x^2 y^4, Errors = -0.060010, -0.057353
+DEAL::f(x,y) = x^2 y^5, Errors = -0.049130, -0.047911
+DEAL::f(x,y) = x^3 y^0, Errors = -0.262944, -0.155343
+DEAL::f(x,y) = x^3 y^1, Errors = -0.121118, -0.096111
+DEAL::f(x,y) = x^3 y^2, Errors = -0.076850, -0.070501
+DEAL::f(x,y) = x^3 y^3, Errors = -0.055777, -0.053901
+DEAL::f(x,y) = x^3 y^4, Errors = -0.043624, -0.042958
+DEAL::f(x,y) = x^3 y^5, Errors = -0.035766, -0.035490
+DEAL::f(x,y) = x^4 y^0, Errors = -0.206246, -0.152327
+DEAL::f(x,y) = x^4 y^1, Errors = -0.094562, -0.082873
+DEAL::f(x,y) = x^4 y^2, Errors = -0.060010, -0.057353
+DEAL::f(x,y) = x^4 y^3, Errors = -0.043624, -0.042958
+DEAL::f(x,y) = x^4 y^4, Errors = -0.034164, -0.033969
+DEAL::f(x,y) = x^4 y^5, Errors = -0.028037, -0.027968
+DEAL::f(x,y) = x^5 y^0, Errors = -0.167904, -0.141111
+DEAL::f(x,y) = x^5 y^1, Errors = -0.077252, -0.071600
+DEAL::f(x,y) = x^5 y^2, Errors = -0.049130, -0.047911
+DEAL::f(x,y) = x^5 y^3, Errors = -0.035766, -0.035490
+DEAL::f(x,y) = x^5 y^4, Errors = -0.028037, -0.027968
+DEAL::f(x,y) = x^5 y^5, Errors = -0.023024, -0.023004
+DEAL:: =========Quadrature Order: 2 =============================== 
+DEAL:: ============================================================ 
+DEAL:: ===============Vertex: 0.000000 0.000000 ============================= 
+DEAL::f(x,y) = x^0 y^0, Errors = 0.083005, -0.001795
+DEAL::f(x,y) = x^0 y^1, Errors = 0.034615, -0.002254
+DEAL::f(x,y) = x^0 y^2, Errors = -0.049939, -0.003643
+DEAL::f(x,y) = x^0 y^3, Errors = -0.105989, -0.004860
+DEAL::f(x,y) = x^0 y^4, Errors = -0.126977, -0.011516
+DEAL::f(x,y) = x^0 y^5, Errors = -0.128734, -0.019894
+DEAL::f(x,y) = x^1 y^0, Errors = 0.034615, -0.002254
+DEAL::f(x,y) = x^1 y^1, Errors = 0.027031, -0.002407
+DEAL::f(x,y) = x^1 y^2, Errors = -0.021731, -0.003410
+DEAL::f(x,y) = x^1 y^3, Errors = -0.049085, -0.007693
+DEAL::f(x,y) = x^1 y^4, Errors = -0.058941, -0.012832
+DEAL::f(x,y) = x^1 y^5, Errors = -0.059742, -0.017132
+DEAL::f(x,y) = x^2 y^0, Errors = -0.049939, -0.003643
+DEAL::f(x,y) = x^2 y^1, Errors = -0.021731, -0.003410
+DEAL::f(x,y) = x^2 y^2, Errors = -0.033894, -0.006861
+DEAL::f(x,y) = x^2 y^3, Errors = -0.041297, -0.010777
+DEAL::f(x,y) = x^2 y^4, Errors = -0.042534, -0.014116
+DEAL::f(x,y) = x^2 y^5, Errors = -0.040552, -0.016356
+DEAL::f(x,y) = x^3 y^0, Errors = -0.105989, -0.004860
+DEAL::f(x,y) = x^3 y^1, Errors = -0.049085, -0.007693
+DEAL::f(x,y) = x^3 y^2, Errors = -0.041297, -0.010777
+DEAL::f(x,y) = x^3 y^3, Errors = -0.038301, -0.013245
+DEAL::f(x,y) = x^3 y^4, Errors = -0.035046, -0.014991
+DEAL::f(x,y) = x^3 y^5, Errors = -0.031558, -0.015882
+DEAL::f(x,y) = x^4 y^0, Errors = -0.126977, -0.011516
+DEAL::f(x,y) = x^4 y^1, Errors = -0.058941, -0.012832
+DEAL::f(x,y) = x^4 y^2, Errors = -0.042534, -0.014116
+DEAL::f(x,y) = x^4 y^3, Errors = -0.035046, -0.014991
+DEAL::f(x,y) = x^4 y^4, Errors = -0.029956, -0.015454
+DEAL::f(x,y) = x^4 y^5, Errors = -0.025973, -0.015409
+DEAL::f(x,y) = x^5 y^0, Errors = -0.128734, -0.019894
+DEAL::f(x,y) = x^5 y^1, Errors = -0.059742, -0.017132
+DEAL::f(x,y) = x^5 y^2, Errors = -0.040552, -0.016356
+DEAL::f(x,y) = x^5 y^3, Errors = -0.031558, -0.015882
+DEAL::f(x,y) = x^5 y^4, Errors = -0.025973, -0.015409
+DEAL::f(x,y) = x^5 y^5, Errors = -0.022011, -0.014742
+DEAL:: =========Quadrature Order: 3 =============================== 
+DEAL:: ============================================================ 
+DEAL:: ===============Vertex: 0.000000 0.000000 ============================= 
+DEAL::f(x,y) = x^0 y^0, Errors = 0.049500, -0.000052
+DEAL::f(x,y) = x^0 y^1, Errors = -0.000455, -0.000092
+DEAL::f(x,y) = x^0 y^2, Errors = 0.000338, -0.000200
+DEAL::f(x,y) = x^0 y^3, Errors = -0.010300, -0.000369
+DEAL::f(x,y) = x^0 y^4, Errors = -0.027273, -0.000603
+DEAL::f(x,y) = x^0 y^5, Errors = -0.042569, -0.000866
+DEAL::f(x,y) = x^1 y^0, Errors = -0.000455, -0.000092
+DEAL::f(x,y) = x^1 y^1, Errors = 0.000150, -0.000105
+DEAL::f(x,y) = x^1 y^2, Errors = 0.001770, -0.000183
+DEAL::f(x,y) = x^1 y^3, Errors = -0.004425, -0.000322
+DEAL::f(x,y) = x^1 y^4, Errors = -0.012927, -0.000525
+DEAL::f(x,y) = x^1 y^5, Errors = -0.020210, -0.000953
+DEAL::f(x,y) = x^2 y^0, Errors = 0.000338, -0.000200
+DEAL::f(x,y) = x^2 y^1, Errors = 0.001770, -0.000183
+DEAL::f(x,y) = x^2 y^2, Errors = 0.002045, -0.000229
+DEAL::f(x,y) = x^2 y^3, Errors = -0.002509, -0.000337
+DEAL::f(x,y) = x^2 y^4, Errors = -0.008202, -0.000658
+DEAL::f(x,y) = x^2 y^5, Errors = -0.012954, -0.001202
+DEAL::f(x,y) = x^3 y^0, Errors = -0.010300, -0.000369
+DEAL::f(x,y) = x^3 y^1, Errors = -0.004425, -0.000322
+DEAL::f(x,y) = x^3 y^2, Errors = -0.002509, -0.000337
+DEAL::f(x,y) = x^3 y^3, Errors = -0.004742, -0.000553
+DEAL::f(x,y) = x^3 y^4, Errors = -0.008084, -0.000986
+DEAL::f(x,y) = x^3 y^5, Errors = -0.010953, -0.001572
+DEAL::f(x,y) = x^4 y^0, Errors = -0.027273, -0.000603
+DEAL::f(x,y) = x^4 y^1, Errors = -0.012927, -0.000525
+DEAL::f(x,y) = x^4 y^2, Errors = -0.008202, -0.000658
+DEAL::f(x,y) = x^4 y^3, Errors = -0.008084, -0.000986
+DEAL::f(x,y) = x^4 y^4, Errors = -0.009419, -0.001464
+DEAL::f(x,y) = x^4 y^5, Errors = -0.010761, -0.002024
+DEAL::f(x,y) = x^5 y^0, Errors = -0.042569, -0.000866
+DEAL::f(x,y) = x^5 y^1, Errors = -0.020210, -0.000953
+DEAL::f(x,y) = x^5 y^2, Errors = -0.012954, -0.001202
+DEAL::f(x,y) = x^5 y^3, Errors = -0.010953, -0.001572
+DEAL::f(x,y) = x^5 y^4, Errors = -0.010761, -0.002024
+DEAL::f(x,y) = x^5 y^5, Errors = -0.010963, -0.002508
+DEAL:: =========Quadrature Order: 4 =============================== 
+DEAL:: ============================================================ 
+DEAL:: ===============Vertex: 0.000000 0.000000 ============================= 
+DEAL::f(x,y) = x^0 y^0, Errors = 0.021072, -0.000002
+DEAL::f(x,y) = x^0 y^1, Errors = 0.000947, -0.000003
+DEAL::f(x,y) = x^0 y^2, Errors = -0.000061, -0.000009
+DEAL::f(x,y) = x^0 y^3, Errors = 0.000186, -0.000022
+DEAL::f(x,y) = x^0 y^4, Errors = -0.001562, -0.000045
+DEAL::f(x,y) = x^0 y^5, Errors = -0.005915, -0.000080
+DEAL::f(x,y) = x^1 y^0, Errors = 0.000947, -0.000003
+DEAL::f(x,y) = x^1 y^1, Errors = 0.000090, -0.000004
+DEAL::f(x,y) = x^1 y^2, Errors = -0.000191, -0.000008
+DEAL::f(x,y) = x^1 y^3, Errors = -0.000007, -0.000019
+DEAL::f(x,y) = x^1 y^4, Errors = -0.000944, -0.000038
+DEAL::f(x,y) = x^1 y^5, Errors = -0.003100, -0.000069
+DEAL::f(x,y) = x^2 y^0, Errors = -0.000061, -0.000009
+DEAL::f(x,y) = x^2 y^1, Errors = -0.000191, -0.000008
+DEAL::f(x,y) = x^2 y^2, Errors = -0.000083, -0.000011
+DEAL::f(x,y) = x^2 y^3, Errors = 0.000084, -0.000019
+DEAL::f(x,y) = x^2 y^4, Errors = -0.000575, -0.000036
+DEAL::f(x,y) = x^2 y^5, Errors = -0.002016, -0.000063
+DEAL::f(x,y) = x^3 y^0, Errors = 0.000186, -0.000022
+DEAL::f(x,y) = x^3 y^1, Errors = -0.000007, -0.000019
+DEAL::f(x,y) = x^3 y^2, Errors = 0.000084, -0.000019
+DEAL::f(x,y) = x^3 y^3, Errors = 0.000171, -0.000025
+DEAL::f(x,y) = x^3 y^4, Errors = -0.000360, -0.000039
+DEAL::f(x,y) = x^3 y^5, Errors = -0.001451, -0.000070
+DEAL::f(x,y) = x^4 y^0, Errors = -0.001562, -0.000045
+DEAL::f(x,y) = x^4 y^1, Errors = -0.000944, -0.000038
+DEAL::f(x,y) = x^4 y^2, Errors = -0.000575, -0.000036
+DEAL::f(x,y) = x^4 y^3, Errors = -0.000360, -0.000039
+DEAL::f(x,y) = x^4 y^4, Errors = -0.000688, -0.000057
+DEAL::f(x,y) = x^4 y^5, Errors = -0.001478, -0.000099
+DEAL::f(x,y) = x^5 y^0, Errors = -0.005915, -0.000080
+DEAL::f(x,y) = x^5 y^1, Errors = -0.003100, -0.000069
+DEAL::f(x,y) = x^5 y^2, Errors = -0.002016, -0.000063
+DEAL::f(x,y) = x^5 y^3, Errors = -0.001451, -0.000070
+DEAL::f(x,y) = x^5 y^4, Errors = -0.001478, -0.000099
+DEAL::f(x,y) = x^5 y^5, Errors = -0.001951, -0.000155
+DEAL:: =========Quadrature Order: 5 =============================== 
+DEAL:: ============================================================ 
+DEAL:: ===============Vertex: 0.000000 0.000000 ============================= 
+DEAL::f(x,y) = x^0 y^0, Errors = 0.009674, 0.000000
+DEAL::f(x,y) = x^0 y^1, Errors = 0.000228, 0.000000
+DEAL::f(x,y) = x^0 y^2, Errors = -0.000007, 0.000000
+DEAL::f(x,y) = x^0 y^3, Errors = -0.000038, -0.000001
+DEAL::f(x,y) = x^0 y^4, Errors = -0.000015, -0.000003
+DEAL::f(x,y) = x^0 y^5, Errors = -0.000326, -0.000006
+DEAL::f(x,y) = x^1 y^0, Errors = 0.000228, 0.000000
+DEAL::f(x,y) = x^1 y^1, Errors = 0.000081, 0.000000
+DEAL::f(x,y) = x^1 y^2, Errors = 0.000029, 0.000000
+DEAL::f(x,y) = x^1 y^3, Errors = -0.000001, -0.000001
+DEAL::f(x,y) = x^1 y^4, Errors = 0.000013, -0.000002
+DEAL::f(x,y) = x^1 y^5, Errors = -0.000154, -0.000005
+DEAL::f(x,y) = x^2 y^0, Errors = -0.000007, 0.000000
+DEAL::f(x,y) = x^2 y^1, Errors = 0.000029, 0.000000
+DEAL::f(x,y) = x^2 y^2, Errors = -0.000006, 0.000000
+DEAL::f(x,y) = x^2 y^3, Errors = -0.000019, -0.000001
+DEAL::f(x,y) = x^2 y^4, Errors = -0.000005, -0.000002
+DEAL::f(x,y) = x^2 y^5, Errors = -0.000118, -0.000005
+DEAL::f(x,y) = x^3 y^0, Errors = -0.000038, -0.000001
+DEAL::f(x,y) = x^3 y^1, Errors = -0.000001, -0.000001
+DEAL::f(x,y) = x^3 y^2, Errors = -0.000019, -0.000001
+DEAL::f(x,y) = x^3 y^3, Errors = -0.000017, -0.000001
+DEAL::f(x,y) = x^3 y^4, Errors = -0.000002, -0.000002
+DEAL::f(x,y) = x^3 y^5, Errors = -0.000087, -0.000004
+DEAL::f(x,y) = x^4 y^0, Errors = -0.000015, -0.000003
+DEAL::f(x,y) = x^4 y^1, Errors = 0.000013, -0.000002
+DEAL::f(x,y) = x^4 y^2, Errors = -0.000005, -0.000002
+DEAL::f(x,y) = x^4 y^3, Errors = -0.000002, -0.000002
+DEAL::f(x,y) = x^4 y^4, Errors = 0.000009, -0.000003
+DEAL::f(x,y) = x^4 y^5, Errors = -0.000062, -0.000005
+DEAL::f(x,y) = x^5 y^0, Errors = -0.000326, -0.000006
+DEAL::f(x,y) = x^5 y^1, Errors = -0.000154, -0.000005
+DEAL::f(x,y) = x^5 y^2, Errors = -0.000118, -0.000005
+DEAL::f(x,y) = x^5 y^3, Errors = -0.000087, -0.000004
+DEAL::f(x,y) = x^5 y^4, Errors = -0.000062, -0.000005
+DEAL::f(x,y) = x^5 y^5, Errors = -0.000110, -0.000007
+DEAL:: =========Quadrature Order: 6 =============================== 
+DEAL:: ============================================================ 
+DEAL:: ===============Vertex: 0.000000 0.000000 ============================= 
+DEAL::f(x,y) = x^0 y^0, Errors = 0.004879, 0.000000
+DEAL::f(x,y) = x^0 y^1, Errors = 0.000082, 0.000000
+DEAL::f(x,y) = x^0 y^2, Errors = -0.000001, 0.000000
+DEAL::f(x,y) = x^0 y^3, Errors = -0.000001, 0.000000
+DEAL::f(x,y) = x^0 y^4, Errors = -0.000004, 0.000000
+DEAL::f(x,y) = x^0 y^5, Errors = -0.000008, 0.000000
+DEAL::f(x,y) = x^1 y^0, Errors = 0.000082, 0.000000
+DEAL::f(x,y) = x^1 y^1, Errors = 0.000009, 0.000000
+DEAL::f(x,y) = x^1 y^2, Errors = 0.000001, 0.000000
+DEAL::f(x,y) = x^1 y^3, Errors = 0.000000, 0.000000
+DEAL::f(x,y) = x^1 y^4, Errors = -0.000003, 0.000000
+DEAL::f(x,y) = x^1 y^5, Errors = -0.000003, 0.000000
+DEAL::f(x,y) = x^2 y^0, Errors = -0.000001, 0.000000
+DEAL::f(x,y) = x^2 y^1, Errors = 0.000001, 0.000000
+DEAL::f(x,y) = x^2 y^2, Errors = 0.000004, 0.000000
+DEAL::f(x,y) = x^2 y^3, Errors = 0.000002, 0.000000
+DEAL::f(x,y) = x^2 y^4, Errors = 0.000000, 0.000000
+DEAL::f(x,y) = x^2 y^5, Errors = 0.000000, 0.000000
+DEAL::f(x,y) = x^3 y^0, Errors = -0.000001, 0.000000
+DEAL::f(x,y) = x^3 y^1, Errors = 0.000000, 0.000000
+DEAL::f(x,y) = x^3 y^2, Errors = 0.000002, 0.000000
+DEAL::f(x,y) = x^3 y^3, Errors = 0.000000, 0.000000
+DEAL::f(x,y) = x^3 y^4, Errors = -0.000002, 0.000000
+DEAL::f(x,y) = x^3 y^5, Errors = -0.000001, 0.000000
+DEAL::f(x,y) = x^4 y^0, Errors = -0.000004, 0.000000
+DEAL::f(x,y) = x^4 y^1, Errors = -0.000003, 0.000000
+DEAL::f(x,y) = x^4 y^2, Errors = 0.000000, 0.000000
+DEAL::f(x,y) = x^4 y^3, Errors = -0.000002, 0.000000
+DEAL::f(x,y) = x^4 y^4, Errors = -0.000002, 0.000000
+DEAL::f(x,y) = x^4 y^5, Errors = -0.000002, 0.000000
+DEAL::f(x,y) = x^5 y^0, Errors = -0.000008, 0.000000
+DEAL::f(x,y) = x^5 y^1, Errors = -0.000003, 0.000000
+DEAL::f(x,y) = x^5 y^2, Errors = 0.000000, 0.000000
+DEAL::f(x,y) = x^5 y^3, Errors = -0.000001, 0.000000
+DEAL::f(x,y) = x^5 y^4, Errors = -0.000002, 0.000000
+DEAL::f(x,y) = x^5 y^5, Errors = -0.000001, 0.000000

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.