From: Luca Heltai Date: Thu, 21 Dec 2017 16:49:09 +0000 (+0100) Subject: Added tests with singular point in the middle. X-Git-Tag: v9.0.0-rc1~583^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fa366c2e88d54a6dbc8016c61be38356da756813;p=dealii.git 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 index 0000000000..16dca16f26 --- /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 + +// all include files needed for the program +#include +#include +#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 index 0000000000..f86a6bacbc --- /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 a6e92c1d71..c972949410 100644 --- 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