From: Nicola Giuliani <ngiuliani@sissa.it>
Date: Wed, 25 Oct 2017 13:15:01 +0000 (+0200)
Subject: added test for singular integration
X-Git-Tag: v9.0.0-rc1~875^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5310%2Fhead;p=dealii.git

added test for singular integration
---

diff --git a/doc/news/changes/minor/20171025NicolaGiuliani b/doc/news/changes/minor/20171025NicolaGiuliani
new file mode 100644
index 0000000000..0654b9eb27
--- /dev/null
+++ b/doc/news/changes/minor/20171025NicolaGiuliani
@@ -0,0 +1,4 @@
+New: Removed a typo in QTelles and added a new singular integration test.
+
+<br>
+(Nicola Giuliani, 25/10/2017
diff --git a/tests/codim_one/integrate_one_over_r_telles_middle.cc b/tests/codim_one/integrate_one_over_r_telles_middle.cc
new file mode 100644
index 0000000000..d236633411
--- /dev/null
+++ b/tests/codim_one/integrate_one_over_r_telles_middle.cc
@@ -0,0 +1,160 @@
+// ---------------------------------------------------------------------
+// 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];
+}
diff --git a/tests/codim_one/integrate_one_over_r_telles_middle.output b/tests/codim_one/integrate_one_over_r_telles_middle.output
new file mode 100644
index 0000000000..971d0b885a
--- /dev/null
+++ b/tests/codim_one/integrate_one_over_r_telles_middle.output
@@ -0,0 +1,146 @@
+
+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