From a58d359d3fa280553a451ba422ed68297ef9b51a Mon Sep 17 00:00:00 2001 From: Nicola Giuliani Date: Wed, 25 Oct 2017 15:15:01 +0200 Subject: [PATCH] added test for singular integration --- doc/news/changes/minor/20171025NicolaGiuliani | 4 + .../integrate_one_over_r_telles_middle.cc | 160 ++++++++++++++++++ .../integrate_one_over_r_telles_middle.output | 146 ++++++++++++++++ 3 files changed, 310 insertions(+) create mode 100644 doc/news/changes/minor/20171025NicolaGiuliani create mode 100644 tests/codim_one/integrate_one_over_r_telles_middle.cc create mode 100644 tests/codim_one/integrate_one_over_r_telles_middle.output 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. + +
+(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 + +// all include files needed for the program +#include +#include +#include + +#include + +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 -- 2.39.5