--- /dev/null
+// ---------------------------------------------------------------------
+// Copyright (C) 2005 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+#include "../tests.h"
+
+// 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 <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 <string>
+
+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_middle (
+ const unsigned int i, const unsigned int j);
+
+int
+main ()
+{
+ initlog();
+ 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 [0.5,0.5]." << endl
+ << endl;
+
+
+ Point <2> center(0.5,0.5);
+ for (unsigned int m=1; m<6; ++m)
+ {
+ deallog << " =========Quadrature Order: " << m
+ << " =============================== " << endl;
+ deallog << " ============================================================ " << endl;
+ deallog << " ===============Vertex: " << center
+ << " ============================= " << endl;
+ QTelles<2> quad(4*m, center);
+ QGaussOneOverR<2> quad_2(m, center, true);
+
+
+ for (unsigned int i = 0; i < 5; ++i)
+ {
+ for (unsigned int j = 0; j < 5; ++j)
+ {
+
+ double exact_integral = exact_integral_one_over_r_middle(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-center[0])*(x-center[0])+(y-center[1])*(y-center[1]));
+ approx_integral += ( pow(x, (double)i) *
+ pow(y, (double)j) / R *
+ quad.weight(q) );
+ }
+
+ for (unsigned int q=0; q< quad_2.size(); ++q)
+ {
+ double x = quad_2.point(q)[0];
+ double y = quad_2.point(q)[1];
+ double R = sqrt((x-center[0])*(x-center[0])+(y-center[1])*(y-center[1]));
+ approx_integral_2 += ( pow(x, (double)i) *
+ pow(y, (double)j) / R *
+ quad_2.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_middle(const unsigned int i,
+ const unsigned int j)
+{
+ Assert(i<6, ExcNotImplemented());
+ Assert(j<6, ExcNotImplemented());
+
+// The integrals are computed using the following Mathematica snippet of
+// code:
+//
+// x0 = 0.5
+// y0 = 0.5
+// Do[Do[Print["v[0][", n, "][", m, "]=",
+// NumberForm[
+// NIntegrate[
+// x^n*y^m/Sqrt[(x - x0)^2 + (y - y0)^2], {x, 0, 1}, {y, 0, 1},
+// MaxRecursion -> 10000, PrecisionGoal -> 9], 9], ";"], {n, 0,
+// 4}], {m, 0, 4}]
+
+
+ static double v[1][6][6] =
+ {
+ {
+ { 0}
+ }
+ };
+ if (v[0][0][0] == 0)
+ {
+ v[0][0][0] = 3.52549435;;
+ v[0][1][0] = 1.76274717;
+ v[0][2][0]=1.07267252;
+ v[0][3][0]=0.727635187;
+ v[0][4][0]=0.53316959;
+ v[0][0][1]=1.76274717;
+ v[0][1][1]=0.881373587;
+ v[0][2][1]=0.536336258;
+ v[0][3][1]=0.363817594;
+ v[0][4][1]=0.266584795;
+ v[0][0][2]=1.07267252;
+ v[0][1][2]=0.536336258;
+ v[0][2][2]=0.329313861;
+ v[0][3][2]=0.225802662;
+ v[0][4][2]=0.167105787;
+ v[0][0][3]=0.727635187;
+ v[0][1][3]=0.363817594;
+ v[0][2][3]=0.225802662;
+ v[0][3][3]=0.156795196;
+ v[0][4][3]=0.117366283;
+ v[0][0][4]=0.53316959;
+ v[0][1][4]=0.266584795;
+ v[0][2][4]=0.167105787;
+ v[0][3][4]=0.117366283;
+ v[0][4][4]=0.0887410133;
+ }
+ return v[0][i][j];
+}
--- /dev/null
+
+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 [0.5,0.5].
+DEAL::
+DEAL:: =========Quadrature Order: 1 ===============================
+DEAL:: ============================================================
+DEAL:: ===============Vertex: 0.500000 0.500000 =============================
+DEAL::f(x,y) = x^0 y^0, Errors = 0.735219, -0.125059
+DEAL::f(x,y) = x^0 y^1, Errors = 0.367609, -0.062529
+DEAL::f(x,y) = x^0 y^2, Errors = 0.184404, -0.098068
+DEAL::f(x,y) = x^0 y^3, Errors = 0.092801, -0.115837
+DEAL::f(x,y) = x^0 y^4, Errors = 0.040442, -0.127062
+DEAL::f(x,y) = x^1 y^0, Errors = 0.367609, -0.062529
+DEAL::f(x,y) = x^1 y^1, Errors = 0.183805, -0.031265
+DEAL::f(x,y) = x^1 y^2, Errors = 0.092202, -0.049034
+DEAL::f(x,y) = x^1 y^3, Errors = 0.046401, -0.057919
+DEAL::f(x,y) = x^1 y^4, Errors = 0.020221, -0.063531
+DEAL::f(x,y) = x^2 y^0, Errors = 0.184404, -0.098068
+DEAL::f(x,y) = x^2 y^1, Errors = 0.092202, -0.049034
+DEAL::f(x,y) = x^2 y^2, Errors = 0.046757, -0.052260
+DEAL::f(x,y) = x^2 y^3, Errors = 0.024034, -0.053873
+DEAL::f(x,y) = x^2 y^4, Errors = 0.010439, -0.054296
+DEAL::f(x,y) = x^3 y^0, Errors = 0.092801, -0.115837
+DEAL::f(x,y) = x^3 y^1, Errors = 0.046401, -0.057919
+DEAL::f(x,y) = x^3 y^2, Errors = 0.024034, -0.053873
+DEAL::f(x,y) = x^3 y^3, Errors = 0.012851, -0.051850
+DEAL::f(x,y) = x^3 y^4, Errors = 0.005548, -0.049678
+DEAL::f(x,y) = x^4 y^0, Errors = 0.040442, -0.127062
+DEAL::f(x,y) = x^4 y^1, Errors = 0.020221, -0.063531
+DEAL::f(x,y) = x^4 y^2, Errors = 0.010439, -0.054296
+DEAL::f(x,y) = x^4 y^3, Errors = 0.005548, -0.049678
+DEAL::f(x,y) = x^4 y^4, Errors = 0.001794, -0.045881
+DEAL:: =========Quadrature Order: 2 ===============================
+DEAL:: ============================================================
+DEAL:: ===============Vertex: 0.500000 0.500000 =============================
+DEAL::f(x,y) = x^0 y^0, Errors = 0.123997, -0.003591
+DEAL::f(x,y) = x^0 y^1, Errors = 0.061998, -0.001795
+DEAL::f(x,y) = x^0 y^2, Errors = 0.030989, -0.002719
+DEAL::f(x,y) = x^0 y^3, Errors = 0.015484, -0.003181
+DEAL::f(x,y) = x^0 y^4, Errors = 0.007730, -0.004396
+DEAL::f(x,y) = x^1 y^0, Errors = 0.061998, -0.001795
+DEAL::f(x,y) = x^1 y^1, Errors = 0.030999, -0.000898
+DEAL::f(x,y) = x^1 y^2, Errors = 0.015494, -0.001360
+DEAL::f(x,y) = x^1 y^3, Errors = 0.007742, -0.001590
+DEAL::f(x,y) = x^1 y^4, Errors = 0.003865, -0.002198
+DEAL::f(x,y) = x^2 y^0, Errors = 0.030989, -0.002719
+DEAL::f(x,y) = x^2 y^1, Errors = 0.015494, -0.001360
+DEAL::f(x,y) = x^2 y^2, Errors = 0.007749, -0.001993
+DEAL::f(x,y) = x^2 y^3, Errors = 0.003876, -0.002309
+DEAL::f(x,y) = x^2 y^4, Errors = 0.001939, -0.002940
+DEAL::f(x,y) = x^3 y^0, Errors = 0.015484, -0.003181
+DEAL::f(x,y) = x^3 y^1, Errors = 0.007742, -0.001590
+DEAL::f(x,y) = x^3 y^2, Errors = 0.003876, -0.002309
+DEAL::f(x,y) = x^3 y^3, Errors = 0.001944, -0.002669
+DEAL::f(x,y) = x^3 y^4, Errors = 0.000975, -0.003312
+DEAL::f(x,y) = x^4 y^0, Errors = 0.007730, -0.004396
+DEAL::f(x,y) = x^4 y^1, Errors = 0.003865, -0.002198
+DEAL::f(x,y) = x^4 y^2, Errors = 0.001939, -0.002940
+DEAL::f(x,y) = x^4 y^3, Errors = 0.000975, -0.003312
+DEAL::f(x,y) = x^4 y^4, Errors = 0.000492, -0.003909
+DEAL:: =========Quadrature Order: 3 ===============================
+DEAL:: ============================================================
+DEAL:: ===============Vertex: 0.500000 0.500000 =============================
+DEAL::f(x,y) = x^0 y^0, Errors = 0.040837, -0.000104
+DEAL::f(x,y) = x^0 y^1, Errors = 0.020419, -0.000052
+DEAL::f(x,y) = x^0 y^2, Errors = 0.010209, -0.000126
+DEAL::f(x,y) = x^0 y^3, Errors = 0.005104, -0.000163
+DEAL::f(x,y) = x^0 y^4, Errors = 0.002552, -0.000232
+DEAL::f(x,y) = x^1 y^0, Errors = 0.020419, -0.000052
+DEAL::f(x,y) = x^1 y^1, Errors = 0.010209, -0.000026
+DEAL::f(x,y) = x^1 y^2, Errors = 0.005104, -0.000063
+DEAL::f(x,y) = x^1 y^3, Errors = 0.002552, -0.000081
+DEAL::f(x,y) = x^1 y^4, Errors = 0.001276, -0.000116
+DEAL::f(x,y) = x^2 y^0, Errors = 0.010209, -0.000126
+DEAL::f(x,y) = x^2 y^1, Errors = 0.005104, -0.000063
+DEAL::f(x,y) = x^2 y^2, Errors = 0.002552, -0.000085
+DEAL::f(x,y) = x^2 y^3, Errors = 0.001276, -0.000096
+DEAL::f(x,y) = x^2 y^4, Errors = 0.000638, -0.000128
+DEAL::f(x,y) = x^3 y^0, Errors = 0.005104, -0.000163
+DEAL::f(x,y) = x^3 y^1, Errors = 0.002552, -0.000081
+DEAL::f(x,y) = x^3 y^2, Errors = 0.001276, -0.000096
+DEAL::f(x,y) = x^3 y^3, Errors = 0.000638, -0.000103
+DEAL::f(x,y) = x^3 y^4, Errors = 0.000319, -0.000134
+DEAL::f(x,y) = x^4 y^0, Errors = 0.002552, -0.000232
+DEAL::f(x,y) = x^4 y^1, Errors = 0.001276, -0.000116
+DEAL::f(x,y) = x^4 y^2, Errors = 0.000638, -0.000128
+DEAL::f(x,y) = x^4 y^3, Errors = 0.000319, -0.000134
+DEAL::f(x,y) = x^4 y^4, Errors = 0.000159, -0.000166
+DEAL:: =========Quadrature Order: 4 ===============================
+DEAL:: ============================================================
+DEAL:: ===============Vertex: 0.500000 0.500000 =============================
+DEAL::f(x,y) = x^0 y^0, Errors = 0.018129, -0.000003
+DEAL::f(x,y) = x^0 y^1, Errors = 0.009064, -0.000002
+DEAL::f(x,y) = x^0 y^2, Errors = 0.004532, -0.000006
+DEAL::f(x,y) = x^0 y^3, Errors = 0.002266, -0.000007
+DEAL::f(x,y) = x^0 y^4, Errors = 0.001133, -0.000013
+DEAL::f(x,y) = x^1 y^0, Errors = 0.009064, -0.000002
+DEAL::f(x,y) = x^1 y^1, Errors = 0.004532, -0.000001
+DEAL::f(x,y) = x^1 y^2, Errors = 0.002266, -0.000003
+DEAL::f(x,y) = x^1 y^3, Errors = 0.001133, -0.000004
+DEAL::f(x,y) = x^1 y^4, Errors = 0.000566, -0.000006
+DEAL::f(x,y) = x^2 y^0, Errors = 0.004532, -0.000006
+DEAL::f(x,y) = x^2 y^1, Errors = 0.002266, -0.000003
+DEAL::f(x,y) = x^2 y^2, Errors = 0.001133, -0.000004
+DEAL::f(x,y) = x^2 y^3, Errors = 0.000566, -0.000005
+DEAL::f(x,y) = x^2 y^4, Errors = 0.000283, -0.000007
+DEAL::f(x,y) = x^3 y^0, Errors = 0.002266, -0.000007
+DEAL::f(x,y) = x^3 y^1, Errors = 0.001133, -0.000004
+DEAL::f(x,y) = x^3 y^2, Errors = 0.000566, -0.000005
+DEAL::f(x,y) = x^3 y^3, Errors = 0.000283, -0.000005
+DEAL::f(x,y) = x^3 y^4, Errors = 0.000142, -0.000007
+DEAL::f(x,y) = x^4 y^0, Errors = 0.001133, -0.000013
+DEAL::f(x,y) = x^4 y^1, Errors = 0.000566, -0.000006
+DEAL::f(x,y) = x^4 y^2, Errors = 0.000283, -0.000007
+DEAL::f(x,y) = x^4 y^3, Errors = 0.000142, -0.000007
+DEAL::f(x,y) = x^4 y^4, Errors = 0.000071, -0.000009
+DEAL:: =========Quadrature Order: 5 ===============================
+DEAL:: ============================================================
+DEAL:: ===============Vertex: 0.500000 0.500000 =============================
+DEAL::f(x,y) = x^0 y^0, Errors = 0.009556, 0.000000
+DEAL::f(x,y) = x^0 y^1, Errors = 0.004778, 0.000000
+DEAL::f(x,y) = x^0 y^2, Errors = 0.002389, 0.000000
+DEAL::f(x,y) = x^0 y^3, Errors = 0.001194, 0.000000
+DEAL::f(x,y) = x^0 y^4, Errors = 0.000597, -0.000001
+DEAL::f(x,y) = x^1 y^0, Errors = 0.004778, 0.000000
+DEAL::f(x,y) = x^1 y^1, Errors = 0.002389, 0.000000
+DEAL::f(x,y) = x^1 y^2, Errors = 0.001194, 0.000000
+DEAL::f(x,y) = x^1 y^3, Errors = 0.000597, 0.000000
+DEAL::f(x,y) = x^1 y^4, Errors = 0.000299, 0.000000
+DEAL::f(x,y) = x^2 y^0, Errors = 0.002389, 0.000000
+DEAL::f(x,y) = x^2 y^1, Errors = 0.001194, 0.000000
+DEAL::f(x,y) = x^2 y^2, Errors = 0.000597, 0.000000
+DEAL::f(x,y) = x^2 y^3, Errors = 0.000299, 0.000000
+DEAL::f(x,y) = x^2 y^4, Errors = 0.000149, 0.000000
+DEAL::f(x,y) = x^3 y^0, Errors = 0.001194, 0.000000
+DEAL::f(x,y) = x^3 y^1, Errors = 0.000597, 0.000000
+DEAL::f(x,y) = x^3 y^2, Errors = 0.000299, 0.000000
+DEAL::f(x,y) = x^3 y^3, Errors = 0.000149, 0.000000
+DEAL::f(x,y) = x^3 y^4, Errors = 0.000075, 0.000000
+DEAL::f(x,y) = x^4 y^0, Errors = 0.000597, -0.000001
+DEAL::f(x,y) = x^4 y^1, Errors = 0.000299, 0.000000
+DEAL::f(x,y) = x^4 y^2, Errors = 0.000149, 0.000000
+DEAL::f(x,y) = x^4 y^3, Errors = 0.000075, 0.000000
+DEAL::f(x,y) = x^4 y^4, Errors = 0.000037, 0.000000