Added tests with singular point in the middle.

diff --git a/tests/base/quadrature_simplex_07.cc b/tests/base/quadrature_simplex_07.cc
new file mode 100644 (file)
index 0000000..16dca16
--- /dev/null
+++ b/tests/base/quadrature_simplex_07.cc
@@ -0,0 +1,78 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 "../tests.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 "simplex.h"
+
+
+int main()
+{
+  initlog();
+
+  deallog << std::endl
+          << "Calculation of the integral of f(x,y)*1/R on [0,1]x[0,1]" << std::endl
+          << "for f(x,y) = x^i y^j, with i,j ranging from 0 to 5, and R being" << std::endl
+          << "the distance from (x,y) to [0.5,0.5]." << std::endl
+          << std::endl;
+
+  double eps = 1e-10;
+
+  //           m  i  j
+  double error[5][6][6] = {{{0}}};
+
+  for (unsigned int m=0; m<5; ++m)
+    {
+      auto split_point = Point<2>(.5, .5);
+
+      QSplit<2> quad(QTrianglePolar(m+1), split_point);
+
+      for (unsigned int i=0; i<6; ++i)
+        for (unsigned int j=0; j<6; ++j)
+          {
+            double exact_integral  = exact_integral_one_over_r_middle(i,j);
+            double approx_integral = 0;
+
+            for (unsigned int q=0; q< quad.size(); ++q)
+              {
+                double x = quad.point(q)[0];
+                double y = quad.point(q)[1];
+                approx_integral += ( pow(x, (double)i) *
+                                     pow(y, (double)j) *
+                                     quad.weight(q) /
+                                     (quad.point(q)-split_point).norm());
+              }
+            error[m][i][j] = approx_integral - exact_integral;
+          }
+    }
+
+  for (unsigned int i=0; i<6; ++i)
+    for (unsigned int j=0; j<6; ++j)
+      {
+        deallog << "======= f(x,y) = x^" << i
+                << " y^" << j << std::endl;
+
+        for (unsigned int m=0; m<5; ++m)
+          deallog << "Order[" << m + 1 << "], error = "
+                  << error[m][i][j]  << std::endl;
+      }
+}

diff --git a/tests/base/quadrature_simplex_07.output b/tests/base/quadrature_simplex_07.output
new file mode 100644 (file)
index 0000000..f86a6ba
--- /dev/null
+++ b/tests/base/quadrature_simplex_07.output
@@ -0,0 +1,222 @@
+
+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::======= f(x,y) = x^0 y^0
+DEAL::Order[1], error = -0.383902
+DEAL::Order[2], error = -0.0307244
+DEAL::Order[3], error = -0.00229163
+DEAL::Order[4], error = -0.000167737
+DEAL::Order[5], error = -1.21893e-05
+DEAL::======= f(x,y) = x^0 y^1
+DEAL::Order[1], error = -0.191951
+DEAL::Order[2], error = -0.0153622
+DEAL::Order[3], error = -0.00114581
+DEAL::Order[4], error = -8.38636e-05
+DEAL::Order[5], error = -6.08965e-06
+DEAL::======= f(x,y) = x^0 y^2
+DEAL::Order[1], error = -0.189100
+DEAL::Order[2], error = -0.0187842
+DEAL::Order[3], error = -0.00208103
+DEAL::Order[4], error = -0.000217194
+DEAL::Order[5], error = -2.15931e-05
+DEAL::======= f(x,y) = x^0 y^3
+DEAL::Order[1], error = -0.187674
+DEAL::Order[2], error = -0.0204952
+DEAL::Order[3], error = -0.00254864
+DEAL::Order[4], error = -0.000283851
+DEAL::Order[5], error = -2.93368e-05
+DEAL::======= f(x,y) = x^0 y^4
+DEAL::Order[1], error = -0.183422
+DEAL::Order[2], error = -0.0220195
+DEAL::Order[3], error = -0.00313498
+DEAL::Order[4], error = -0.000409087
+DEAL::Order[5], error = -4.93578e-05
+DEAL::======= f(x,y) = x^0 y^5
+DEAL::Order[1], error = 0.236233
+DEAL::Order[2], error = 0.390539
+DEAL::Order[3], error = 0.410209
+DEAL::Order[4], error = 0.413426
+DEAL::Order[5], error = 0.413914
+DEAL::======= f(x,y) = x^1 y^0
+DEAL::Order[1], error = -0.191951
+DEAL::Order[2], error = -0.0153622
+DEAL::Order[3], error = -0.00114581
+DEAL::Order[4], error = -8.38636e-05
+DEAL::Order[5], error = -6.08965e-06
+DEAL::======= f(x,y) = x^1 y^1
+DEAL::Order[1], error = -0.0959754
+DEAL::Order[2], error = -0.00768110
+DEAL::Order[3], error = -0.000572906
+DEAL::Order[4], error = -4.19338e-05
+DEAL::Order[5], error = -3.04683e-06
+DEAL::======= f(x,y) = x^1 y^2
+DEAL::Order[1], error = -0.0945498
+DEAL::Order[2], error = -0.00939208
+DEAL::Order[3], error = -0.00104052
+DEAL::Order[4], error = -0.000108595
+DEAL::Order[5], error = -1.07945e-05
+DEAL::======= f(x,y) = x^1 y^3
+DEAL::Order[1], error = -0.0938370
+DEAL::Order[2], error = -0.0102476
+DEAL::Order[3], error = -0.00127432
+DEAL::Order[4], error = -0.000141926
+DEAL::Order[5], error = -1.46689e-05
+DEAL::======= f(x,y) = x^1 y^4
+DEAL::Order[1], error = -0.0917110
+DEAL::Order[2], error = -0.0110098
+DEAL::Order[3], error = -0.00156749
+DEAL::Order[4], error = -0.000204544
+DEAL::Order[5], error = -2.46789e-05
+DEAL::======= f(x,y) = x^1 y^5
+DEAL::Order[1], error = 0.118117
+DEAL::Order[2], error = 0.195270
+DEAL::Order[3], error = 0.205105
+DEAL::Order[4], error = 0.206713
+DEAL::Order[5], error = 0.206957
+DEAL::======= f(x,y) = x^2 y^0
+DEAL::Order[1], error = -0.189100
+DEAL::Order[2], error = -0.0187842
+DEAL::Order[3], error = -0.00208103
+DEAL::Order[4], error = -0.000217194
+DEAL::Order[5], error = -2.15931e-05
+DEAL::======= f(x,y) = x^2 y^1
+DEAL::Order[1], error = -0.0945498
+DEAL::Order[2], error = -0.00939208
+DEAL::Order[3], error = -0.00104052
+DEAL::Order[4], error = -0.000108595
+DEAL::Order[5], error = -1.07945e-05
+DEAL::======= f(x,y) = x^2 y^2
+DEAL::Order[1], error = -0.0838769
+DEAL::Order[2], error = -0.0107068
+DEAL::Order[3], error = -0.00132108
+DEAL::Order[4], error = -0.000148592
+DEAL::Order[5], error = -1.54436e-05
+DEAL::======= f(x,y) = x^2 y^3
+DEAL::Order[1], error = -0.0785405
+DEAL::Order[2], error = -0.0113642
+DEAL::Order[3], error = -0.00146136
+DEAL::Order[4], error = -0.000168590
+DEAL::Order[5], error = -1.77676e-05
+DEAL::======= f(x,y) = x^2 y^4
+DEAL::Order[1], error = -0.0735330
+DEAL::Order[2], error = -0.0119613
+DEAL::Order[3], error = -0.00168353
+DEAL::Order[4], error = -0.000217502
+DEAL::Order[5], error = -2.57948e-05
+DEAL::======= f(x,y) = x^2 y^5
+DEAL::Order[1], error = 0.0621262
+DEAL::Order[2], error = 0.118288
+DEAL::Order[3], error = 0.128869
+DEAL::Order[4], error = 0.130535
+DEAL::Order[5], error = 0.130779
+DEAL::======= f(x,y) = x^3 y^0
+DEAL::Order[1], error = -0.187674
+DEAL::Order[2], error = -0.0204952
+DEAL::Order[3], error = -0.00254864
+DEAL::Order[4], error = -0.000283851
+DEAL::Order[5], error = -2.93368e-05
+DEAL::======= f(x,y) = x^3 y^1
+DEAL::Order[1], error = -0.0938370
+DEAL::Order[2], error = -0.0102476
+DEAL::Order[3], error = -0.00127432
+DEAL::Order[4], error = -0.000141926
+DEAL::Order[5], error = -1.46689e-05
+DEAL::======= f(x,y) = x^3 y^2
+DEAL::Order[1], error = -0.0785405
+DEAL::Order[2], error = -0.0113642
+DEAL::Order[3], error = -0.00146136
+DEAL::Order[4], error = -0.000168590
+DEAL::Order[5], error = -1.77676e-05
+DEAL::======= f(x,y) = x^3 y^3
+DEAL::Order[1], error = -0.0708923
+DEAL::Order[2], error = -0.0119226
+DEAL::Order[3], error = -0.00155489
+DEAL::Order[4], error = -0.000181922
+DEAL::Order[5], error = -1.93169e-05
+DEAL::======= f(x,y) = x^3 y^4
+DEAL::Order[1], error = -0.0644439
+DEAL::Order[2], error = -0.0124370
+DEAL::Order[3], error = -0.00174155
+DEAL::Order[4], error = -0.000223980
+DEAL::Order[5], error = -2.63528e-05
+DEAL::======= f(x,y) = x^3 y^5
+DEAL::Order[1], error = 0.0341311
+DEAL::Order[2], error = 0.0797970
+DEAL::Order[3], error = 0.0907519
+DEAL::Order[4], error = 0.0924462
+DEAL::Order[5], error = 0.0926905
+DEAL::======= f(x,y) = x^4 y^0
+DEAL::Order[1], error = -0.183422
+DEAL::Order[2], error = -0.0220195
+DEAL::Order[3], error = -0.00313498
+DEAL::Order[4], error = -0.000409087
+DEAL::Order[5], error = -4.93578e-05
+DEAL::======= f(x,y) = x^4 y^1
+DEAL::Order[1], error = -0.0917110
+DEAL::Order[2], error = -0.0110098
+DEAL::Order[3], error = -0.00156749
+DEAL::Order[4], error = -0.000204544
+DEAL::Order[5], error = -2.46789e-05
+DEAL::======= f(x,y) = x^4 y^2
+DEAL::Order[1], error = -0.0735330
+DEAL::Order[2], error = -0.0119613
+DEAL::Order[3], error = -0.00168353
+DEAL::Order[4], error = -0.000217502
+DEAL::Order[5], error = -2.57948e-05
+DEAL::======= f(x,y) = x^4 y^3
+DEAL::Order[1], error = -0.0644439
+DEAL::Order[2], error = -0.0124370
+DEAL::Order[3], error = -0.00174155
+DEAL::Order[4], error = -0.000223980
+DEAL::Order[5], error = -2.63528e-05
+DEAL::======= f(x,y) = x^4 y^4
+DEAL::Order[1], error = -0.0572944
+DEAL::Order[2], error = -0.0128503
+DEAL::Order[3], error = -0.00190082
+DEAL::Order[4], error = -0.000259133
+DEAL::Order[5], error = -3.22754e-05
+DEAL::======= f(x,y) = x^4 y^5
+DEAL::Order[1], error = 0.0195583
+DEAL::Order[2], error = 0.0574407
+DEAL::Order[3], error = 0.0685621
+DEAL::Order[4], error = 0.0703642
+DEAL::Order[5], error = 0.0706320
+DEAL::======= f(x,y) = x^5 y^0
+DEAL::Order[1], error = 0.236233
+DEAL::Order[2], error = 0.390539
+DEAL::Order[3], error = 0.410209
+DEAL::Order[4], error = 0.413426
+DEAL::Order[5], error = 0.413914
+DEAL::======= f(x,y) = x^5 y^1
+DEAL::Order[1], error = 0.118117
+DEAL::Order[2], error = 0.195270
+DEAL::Order[3], error = 0.205105
+DEAL::Order[4], error = 0.206713
+DEAL::Order[5], error = 0.206957
+DEAL::======= f(x,y) = x^5 y^2
+DEAL::Order[1], error = 0.0621262
+DEAL::Order[2], error = 0.118288
+DEAL::Order[3], error = 0.128869
+DEAL::Order[4], error = 0.130535
+DEAL::Order[5], error = 0.130779
+DEAL::======= f(x,y) = x^5 y^3
+DEAL::Order[1], error = 0.0341311
+DEAL::Order[2], error = 0.0797970
+DEAL::Order[3], error = 0.0907519
+DEAL::Order[4], error = 0.0924462
+DEAL::Order[5], error = 0.0926905
+DEAL::======= f(x,y) = x^5 y^4
+DEAL::Order[1], error = 0.0195583
+DEAL::Order[2], error = 0.0574407
+DEAL::Order[3], error = 0.0685621
+DEAL::Order[4], error = 0.0703642
+DEAL::Order[5], error = 0.0706320
+DEAL::======= f(x,y) = x^5 y^5
+DEAL::Order[1], error = 0.0116966
+DEAL::Order[2], error = 0.0431516
+DEAL::Order[3], error = 0.0543363
+DEAL::Order[4], error = 0.0562857
+DEAL::Order[5], error = 0.0565886

diff --git a/tests/base/simplex.h b/tests/base/simplex.h
index a6e92c1d710b9c9fc7823a7c5a72c99fee554f5a..c972949410a9e4a8e11abde814b2e4eded9fa80a 100644 (file)
--- a/tests/base/simplex.h
+++ b/tests/base/simplex.h
@@ -236,4 +236,57 @@ double exact_integral_one_over_r(const unsigned int vertex_index,
 }
 
 
+
+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[", 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[6][6] = {{0}};
+
+  if (v[0][0] == 0)
+    {
+      v[0][0] = 3.52549435;;
+      v[1][0] = 1.76274717;
+      v[2][0]=1.07267252;
+      v[3][0]=0.727635187;
+      v[4][0]=0.53316959;
+      v[0][1]=1.76274717;
+      v[1][1]=0.881373587;
+      v[2][1]=0.536336258;
+      v[3][1]=0.363817594;
+      v[4][1]=0.266584795;
+      v[0][2]=1.07267252;
+      v[1][2]=0.536336258;
+      v[2][2]=0.329313861;
+      v[3][2]=0.225802662;
+      v[4][2]=0.167105787;
+      v[0][3]=0.727635187;
+      v[1][3]=0.363817594;
+      v[2][3]=0.225802662;
+      v[3][3]=0.156795196;
+      v[4][3]=0.117366283;
+      v[0][4]=0.53316959;
+      v[1][4]=0.266584795;
+      v[2][4]=0.167105787;
+      v[3][4]=0.117366283;
+      v[4][4]=0.0887410133;
+    }
+  return v[i][j];
+}
+
 #endif