From 2150e5caad229a7fb7488b44bcffbe9f0f1bfdcf Mon Sep 17 00:00:00 2001 From: Nicola Giuliani Date: Mon, 30 Mar 2015 18:22:05 +0200 Subject: [PATCH] Implemented the Telles quadrature formulas and their test implemented some suggestions from the authors applied the comments on the pull request bug removed class made more readable class made more readable testing over m changes.h ok corrected the telles test adjust the tolerance corrected the formula corrected the change.h added a description for the dim > 1 cases added the reference for Telles --- doc/news/changes.h | 13 +- include/deal.II/base/quadrature_lib.h | 70 ++++ source/base/quadrature_lib.cc | 116 +++++++ .../codim_one/integrate_one_over_r_telles.cc | 302 ++++++++++++++++++ .../integrate_one_over_r_telles.output | 240 ++++++++++++++ 5 files changed, 738 insertions(+), 3 deletions(-) create mode 100644 tests/codim_one/integrate_one_over_r_telles.cc create mode 100644 tests/codim_one/integrate_one_over_r_telles.output diff --git a/doc/news/changes.h b/doc/news/changes.h index 1dabada2bb..50c4e72e2c 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -402,9 +402,16 @@ inconvenience this causes.

Specific improvements

    -
  1. New: A function to get a map with all vertices at boundaries - has been added at GridTools::get_all_vertices_at_boundary. - This function will return a map which can be used in functions +
  2. New: There is now a new quadrature formula in quadrature_lib. It is + now possible to use Telles' quadrature rules through the function QTelles + to integrate singular integrals +
    + (Nicola Giuliani, 2015/04/01) +
  3. + +
  4. New: A new function to get a map with all vertices at boundaries + has been added in GridTools::get_all_vertices_at_boundary. + This function will return a map wich can be used in functions like GridTools::Laplace_transform.
    (Fernando Posada, 2015/03/31) diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index b484edb34c..0400b13eac 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -439,6 +439,75 @@ public: const std::pair > &b); }; +/** + * Telles quadrature of arbitrary order. + * + * The coefficients of these quadrature rules are computed using + * a non linear change of variables starting from a Gauss-Legendre + * quadrature formula. + * This is done using a cubic polynomial, + * $n = a x^3 + b x^2 + c x + d$ + * in order to integrate + * a singular integral, with singularity at a given point x_0. + * + * We start from a Gauss Quadrature Formula with arbitrary + * function. Then we apply the cubic variable change. + * In the paper, J.C.F.Telles:A Self-Adaptive Co-ordinate Transformation + * For Efficient Numerical Evaluation of General Boundary Element Integrals. + * International Journal for Numerical Methods in Engineering, vol 24, + * pages 959–973. year 1987, the author applies the transformation on the + * reference cell $[-1, 1]$ getting + * @f{align*} n(1) &= 1, \\ n(-1) &= -1, \\ dn/dx $= 0 at x = x_0, \\ + d2n/dx2 = 0 at x &= x_0 @f} + * We get + * @f{align*} a &= 1/q, \\ b &= -3gamma_bar/q, \\ c &= 3gamma_bar/q, \\ + d &= -b, @f} + * with + * @f{align*} eta_star &= eta_bar^2 - 1, \\ gamma_bar 6 &= nthroot( eta_bar + * eta_star + abs(eta_star) ,3) + nthroot(eta_bar*eta_star - abs(eta_star),3) + + eta_bar, \\ q &= ((gamma-gamma_bar).^3 + gamma_bar*(gamma_bar^2+3)) + /(1+3*gamma_bar^2) @f}. + * Since the library assumes [0,1] as reference interval, we will map + * these values on the proper reference interval in the implementation. + * + * This variable change can be used to integrate singular integrals. + * One example is $\f(x)/abs(x-x_0)$ on the reference interval $[0,1]$, + * where $x_0$ is given at construction time, and is the location of the + * singularity $x_0$, and $f(x)$ is a smooth non singular function. + * + * Singular quadrature formula are rather expensive, nevertheless Telles' + * quadrature formula are much easier to compute with respect to other singular + * integration techniques as Lachat-Watson. + * + * We have implemented the case for $dim = 1$. When we deal the case $dim >1$ + * we have computed the quadrature formula has a tensorial product of one + * dimensional Telles' quadrature formulas considering the different components + * of the singularity. + * + * The weights and functions for Gauss Legendre formula have been tabulated up + * to order 12. + * + * @author Nicola Giuliani, Luca Heltai 2015 + */ +template +class QTelles: public Quadrature +{ +public: + /** + * A constructor that takes a quadrature formula and a singular point as + * argument. The quadrature formula will be mapped using Telles' rule. Make + * sure that the order of the quadrature rule is appropriate for the + * singularity in question. + **/ + QTelles (const Quadrature<1> &base_quad, const Point &singularity); + /** + * A variant of above constructor that takes as parameters the order @p n + * and location of a singularity. A Gauss Legendre quadrature of order n + * will be used + **/ + QTelles (const unsigned int n, const Point &singularity); + +}; /*@}*/ @@ -470,6 +539,7 @@ template <> QWeddle<1>::QWeddle (); template <> QGaussLog<1>::QGaussLog (const unsigned int n, const bool revert); template <> QGaussLogR<1>::QGaussLogR (const unsigned int n, const Point<1> x0, const double alpha, const bool flag); template <> QGaussOneOverR<2>::QGaussOneOverR (const unsigned int n, const unsigned int index, const bool flag); +template <> QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> &singularity); diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index 98304d1e65..829bb52fb3 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -1014,6 +1014,118 @@ QWeddle::QWeddle () Quadrature (QWeddle(), QWeddle<1>()) {} +template +QTelles::QTelles ( + const Quadrature<1> &base_quad, const Point &singularity) + : +/** +* We need the explicit implementation if dim == 1. If dim > 1 we use the +* former implementation and apply a tensorial product to obtain the higher +* dimensions. +**/ + Quadrature( + dim == 2 ? + QAnisotropic( + QTelles<1>(base_quad, Point<1>(singularity[0])), + QTelles<1>(base_quad, Point<1>(singularity[1]))) : + dim == 3 ? + QAnisotropic( + QTelles<1>(base_quad, Point<1>(singularity[0])), + QTelles<1>(base_quad, Point<1>(singularity[1])), + QTelles<1>(base_quad, Point<1>(singularity[2]))) : + Quadrature()) +{ +} + +template +QTelles::QTelles ( + const unsigned int n, const Point &singularity) + : +/** +* In this case we map the standard Gauss Legendre formula using the given +* singularity point coordinates. +**/ + Quadrature(QTelles(QGauss<1>(n), singularity)) +{} + + + +template <> +QTelles<1>::QTelles ( + const Quadrature<1> &base_quad, const Point<1> &singularity) + : +/** +* We explicitly implement the Telles' variable change if dim == 1. +**/ + Quadrature<1>(base_quad) +{ + /** + * We define all the constants to be used in the implementation of + * Telles' rule + **/ + const double eta_bar = singularity[0] * 2. - 1.; + const double eta_star = eta_bar * eta_bar - 1.; + double gamma_bar; + + std::vector> quadrature_points_dummy(quadrature_points.size()); + std::vector weights_dummy(weights.size()); + unsigned int cont = 0; + const double tol = 1e-10; + for (unsigned int d = 0; d < quadrature_points.size(); ++d) + { + if (std::abs(quadrature_points[d][0] - singularity[0]) > tol) + { + quadrature_points_dummy[d-cont] = quadrature_points[d]; + weights_dummy[d-cont] = weights[d]; + } + else + { + // We need to remove the singularity point from the quadrature point + // list. To do so we use the variable cont. + cont = 1; + } + + } + if (cont == 1) + { + quadrature_points.resize(quadrature_points_dummy.size()-1); + weights.resize(weights_dummy.size()-1); + for (unsigned int d = 0; d < quadrature_points.size()-1; ++d) + { + quadrature_points[d] = quadrature_points_dummy[d]; + weights[d] = weights_dummy[d]; + } + } + // We need to check if the singularity is at the boundary of the interval. + if (std::abs(eta_star) <= tol) + { + gamma_bar = std::pow((eta_bar * eta_star + std::abs(eta_star)),1.0 / 3.0) + + std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0) + + eta_bar; + } + else + { + gamma_bar = (eta_bar * eta_star + std::abs(eta_star))/std::abs(eta_bar * eta_star + std::abs(eta_star))* + std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)),1.0 / 3.0) + + (eta_bar * eta_star - std::abs(eta_star))/std::abs(eta_bar * eta_star - std::abs(eta_star))* + std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0) + + eta_bar; + } + for (unsigned int q = 0; q < quadrature_points.size(); ++q) + { + double gamma = quadrature_points[q][0] * 2 - 1; + double eta = (std::pow(gamma - gamma_bar, 3.0) + + gamma_bar * (gamma_bar * gamma_bar + 3)) + / (1 + 3 * gamma_bar * gamma_bar); + + double J = 3 * ((gamma - gamma_bar) *(gamma - gamma_bar)) + / (1 + 3 * gamma_bar * gamma_bar); + + quadrature_points[q][0] = (eta + 1) / 2.0; + weights[q] = J * weights[q]; + + } +} // explicit specialization // note that 1d formulae are specialized by implementation above @@ -1037,4 +1149,8 @@ template class QSorted<1>; template class QSorted<2>; template class QSorted<3>; +template class QTelles<1> ; +template class QTelles<2> ; +template class QTelles<3> ; + DEAL_II_NAMESPACE_CLOSE diff --git a/tests/codim_one/integrate_one_over_r_telles.cc b/tests/codim_one/integrate_one_over_r_telles.cc new file mode 100644 index 0000000000..46d1613530 --- /dev/null +++ b/tests/codim_one/integrate_one_over_r_telles.cc @@ -0,0 +1,302 @@ +// --------------------------------------------------------------------- +// $Id: integrate_one_over_r.cc 30338 2013-08-18 22:02:27Z heltai $ +// +// Copyright (C) 2005 - 2013 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 +#include +#include + +// all include files needed for the program +#include +#include +#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 ( + const unsigned int i, const unsigned int j, + const unsigned int vertex_index); + +ofstream logfile("output"); + +int +main () +{ + deallog.attach(logfile); + deallog.depth_console(0); + 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 four vertices of the square." << endl + << endl; + + + std::vector > vertices = FE_Q<2>(1).get_unit_support_points(); + + for (unsigned int m=1; m<7; ++m) + { + deallog << " =========Quadrature Order: " << m + << " =============================== " << endl; + deallog << " ============================================================ " << endl; + unsigned int index=0; + { + deallog << " ===============Vertex: " << vertices[index] + << " ============================= " << endl; + QTelles<2> quad(m, vertices[index]); + QGaussOneOverR<2> quad2(m, vertices[index]); + + + for (unsigned int i = 0; i < 6; ++i) + { + for (unsigned int j = 0; j < 6; ++j) + { + + double exact_integral = exact_integral_one_over_r(index, 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*x+y*y); + approx_integral += ( pow(x, (double)i) * + pow(y, (double)j) / R * + quad.weight(q) ); + } + + for (unsigned int q=0; q< quad2.size(); ++q) + { + double x = quad2.point(q)[0]; + double y = quad2.point(q)[1]; + double R = sqrt(x*x+y*y); + approx_integral_2 += ( pow(x, (double)i) * + pow(y, (double)j) * + quad2.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(const unsigned int vertex_index, + const unsigned int i, + const unsigned int j) + { + Assert(vertex_index < 4, ExcInternalError()); + Assert(i<6, ExcNotImplemented()); + Assert(j<6, ExcNotImplemented()); + +// The integrals are computed using the following maple snippet of +// code: +// +// singint := proc(index, N, M) +// if index = 0 then +// return int(int(x^N *y^M/sqrt(x^2+y^2), x=0.0..1.0), y=0.0..1.0); +// elif index = 1 then +// return int(int(x^N *y^M/sqrt((x-1)^2+y^2), x=0.0..1.0), y=0.0..1.0); +// elif index = 2 then +// return int(int(x^N *y^M/sqrt(x^2+(y-1)^2), x=0.0..1.0), y=0.0..1.0); +// elif index = 3 then +// return int(int((1-x)^N *(1-y)^M/sqrt(x^2+y^2), x=0.0..1.0), y=0.0..1.0); +// end if; +// end proc; +// Digits := 20; +// for i from 3 to 3 do +// for n from 0 to 5 do +// for m from 0 to 5 do +// v[i+1][n+1][m+1] = sing_int(i, n, m); +// end do; +// end do; +// end do; +// C(v) + + static double v[4][6][6] = + { + { + { 0}}}; + if (v[0][0][0] == 0) + { + v[0][0][0] = 0.17627471740390860505e1; + v[0][0][1] = 0.64779357469631903702e0; + v[0][0][2] = 0.38259785823210634567e0; + v[0][0][3] = 0.26915893322379450224e0; + v[0][0][4] = 0.20702239737104695572e0; + v[0][0][5] = 0.16800109713227567467e0; + v[0][1][0] = 0.64779357469631903702e0; + v[0][1][1] = 0.27614237491539669920e0; + v[0][1][2] = 0.17015838751246776515e0; + v[0][1][3] = 0.12189514164974600651e0; + v[0][1][4] = 0.94658660368131133694e-1; + v[0][1][5] = 0.77263794021029438797e-1; + v[0][2][0] = 0.38259785823210634567e0; + v[0][2][1] = 0.17015838751246776515e0; + v[0][2][2] = 0.10656799507071040471e0; + v[0][2][3] = 0.76947022258735165920e-1; + v[0][2][4] = 0.60022626787495395021e-1; + v[0][2][5] = 0.49131622931360879320e-1; + v[0][3][0] = 0.26915893322379450224e0; + v[0][3][1] = 0.12189514164974600651e0; + v[0][3][2] = 0.76947022258735165919e-1; + v[0][3][3] = 0.55789184535895709637e-1; + v[0][3][4] = 0.43625068213915842136e-1; + v[0][3][5] = 0.35766126849971778500e-1; + v[0][4][0] = 0.20702239737104695572e0; + v[0][4][1] = 0.94658660368131133694e-1; + v[0][4][2] = 0.60022626787495395021e-1; + v[0][4][3] = 0.43625068213915842137e-1; + v[0][4][4] = 0.34164088852375945192e-1; + v[0][4][5] = 0.28037139560980277614e-1; + v[0][5][0] = 0.16800109713227567467e0; + v[0][5][1] = 0.77263794021029438797e-1; + v[0][5][2] = 0.49131622931360879320e-1; + v[0][5][3] = 0.35766126849971778501e-1; + v[0][5][4] = 0.28037139560980277614e-1; + v[0][5][5] = 0.23024181049838367777e-1; + v[1][0][0] = 0.17627471740390860505e1; + v[1][0][1] = 0.64779357469631903702e0; + v[1][0][2] = 0.38259785823210634567e0; + v[1][0][3] = 0.26915893322379450224e0; + v[1][0][4] = 0.20702239737104695572e0; + v[1][0][5] = 0.16800109713227567467e0; + v[1][1][0] = 0.11149535993427670134e1; + v[1][1][1] = 0.37165119978092233782e0; + v[1][1][2] = 0.21243947071963858053e0; + v[1][1][3] = 0.14726379157404849573e0; + v[1][1][4] = 0.11236373700291582202e0; + v[1][1][5] = 0.90737303111246235871e-1; + v[1][2][0] = 0.84975788287855432210e0; + v[1][2][1] = 0.26566721237799340376e0; + v[1][2][2] = 0.14884907827788122009e0; + v[1][2][3] = 0.10231567218303765515e0; + v[1][2][4] = 0.77727703422280083352e-1; + v[1][2][5] = 0.62605132021577676395e-1; + v[1][3][0] = 0.69800109142265347423e0; + v[1][3][1] = 0.20794647083778622837e0; + v[1][3][2] = 0.11487965864809909847e0; + v[1][3][3] = 0.78525390514866270852e-1; + v[1][3][4] = 0.59489228415223897572e-1; + v[1][3][5] = 0.47838457013298217744e-1; + v[1][4][0] = 0.59754668912231692323e0; + v[1][4][1] = 0.17125249387868593878e0; + v[1][4][2] = 0.93606816359052444729e-1; + v[1][4][3] = 0.63728830247554475330e-1; + v[1][4][4] = 0.48187332620207367724e-1; + v[1][4][5] = 0.38708290797416359020e-1; + v[1][5][0] = 0.52527944036356840363e0; + v[1][5][1] = 0.14574366656617935708e0; + v[1][5][2] = 0.78997159795636003667e-1; + v[1][5][3] = 0.53620816423066464705e-1; + v[1][5][4] = 0.40487985967086264433e-1; + v[1][5][5] = 0.32498604596082509165e-1; + v[2][0][0] = 0.17627471740390860505e1; + v[2][0][1] = 0.11149535993427670134e1; + v[2][0][2] = 0.84975788287855432210e0; + v[2][0][3] = 0.69800109142265347419e0; + v[2][0][4] = 0.59754668912231692318e0; + v[2][0][5] = 0.52527944036356840362e0; + v[2][1][0] = 0.64779357469631903702e0; + v[2][1][1] = 0.37165119978092233782e0; + v[2][1][2] = 0.26566721237799340376e0; + v[2][1][3] = 0.20794647083778622835e0; + v[2][1][4] = 0.17125249387868593876e0; + v[2][1][5] = 0.14574366656617935708e0; + v[2][2][0] = 0.38259785823210634567e0; + v[2][2][1] = 0.21243947071963858053e0; + v[2][2][2] = 0.14884907827788122009e0; + v[2][2][3] = 0.11487965864809909845e0; + v[2][2][4] = 0.93606816359052444712e-1; + v[2][2][5] = 0.78997159795636003667e-1; + v[2][3][0] = 0.26915893322379450223e0; + v[2][3][1] = 0.14726379157404849572e0; + v[2][3][2] = 0.10231567218303765514e0; + v[2][3][3] = 0.78525390514866270835e-1; + v[2][3][4] = 0.63728830247554475311e-1; + v[2][3][5] = 0.53620816423066464702e-1; + v[2][4][0] = 0.20702239737104695572e0; + v[2][4][1] = 0.11236373700291582202e0; + v[2][4][2] = 0.77727703422280083352e-1; + v[2][4][3] = 0.59489228415223897563e-1; + v[2][4][4] = 0.48187332620207367713e-1; + v[2][4][5] = 0.40487985967086264434e-1; + v[2][5][0] = 0.16800109713227567468e0; + v[2][5][1] = 0.90737303111246235879e-1; + v[2][5][2] = 0.62605132021577676399e-1; + v[2][5][3] = 0.47838457013298217740e-1; + v[2][5][4] = 0.38708290797416359014e-1; + v[2][5][5] = 0.32498604596082509169e-1; + v[3][0][0] = 0.17627471740390860505e1; + v[3][0][1] = 0.11149535993427670134e1; + v[3][0][2] = 0.84975788287855432210e0; + v[3][0][3] = 0.69800109142265347419e0; + v[3][0][4] = 0.59754668912231692318e0; + v[3][0][5] = 0.52527944036356840362e0; + v[3][1][0] = 0.11149535993427670134e1; + v[3][1][1] = 0.74330239956184467563e0; + v[3][1][2] = 0.58409067050056091834e0; + v[3][1][3] = 0.49005462058486724584e0; + v[3][1][4] = 0.42629419524363098443e0; + v[3][1][5] = 0.37953577379738904654e0; + v[3][2][0] = 0.84975788287855432210e0; + v[3][2][1] = 0.58409067050056091834e0; + v[3][2][2] = 0.46727253640044873467e0; + v[3][2][3] = 0.39698780839518011595e0; + v[3][2][4] = 0.34864851772399749038e0; + v[3][2][5] = 0.31278926702684569312e0; + v[3][3][0] = 0.69800109142265347423e0; + v[3][3][1] = 0.49005462058486724586e0; + v[3][3][2] = 0.39698780839518011599e0; + v[3][3][3] = 0.34027526433872581371e0; + v[3][3][4] = 0.30088082631586196583e0; + v[3][3][5] = 0.27141910362887187844e0; + v[3][4][0] = 0.59754668912231692323e0; + v[3][4][1] = 0.42629419524363098445e0; + v[3][4][2] = 0.34864851772399749044e0; + v[3][4][3] = 0.30088082631586196576e0; + v[3][4][4] = 0.26744962339187730308e0; + v[3][4][5] = 0.24229245314748740295e0; + v[3][5][0] = 0.52527944036356840363e0; + v[3][5][1] = 0.37953577379738904655e0; + v[3][5][2] = 0.31278926702684569301e0; + v[3][5][3] = 0.27141910362887187862e0; + v[3][5][4] = 0.24229245314748740263e0; + v[3][5][5] = 0.22026586649771582089e0; + } + return v[vertex_index][i][j]; + } diff --git a/tests/codim_one/integrate_one_over_r_telles.output b/tests/codim_one/integrate_one_over_r_telles.output new file mode 100644 index 0000000000..5afde93dd1 --- /dev/null +++ b/tests/codim_one/integrate_one_over_r_telles.output @@ -0,0 +1,240 @@ + +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 four vertices of the square. +DEAL:: +DEAL:: =========Quadrature Order: 1 =============================== +DEAL:: ============================================================ +DEAL:: ===============Vertex: 0.000000 0.000000 ============================= +DEAL::f(x,y) = x^0 y^0, Errors = 1.419233, -0.062529 +DEAL::f(x,y) = x^0 y^1, Errors = -0.250046, -0.046676 +DEAL::f(x,y) = x^0 y^2, Errors = -0.332879, -0.133607 +DEAL::f(x,y) = x^0 y^3, Errors = -0.262944, -0.155343 +DEAL::f(x,y) = x^0 y^4, Errors = -0.206246, -0.152327 +DEAL::f(x,y) = x^0 y^5, Errors = -0.167904, -0.141111 +DEAL::f(x,y) = x^1 y^0, Errors = -0.250046, -0.046676 +DEAL::f(x,y) = x^1 y^1, Errors = -0.226424, -0.100079 +DEAL::f(x,y) = x^1 y^2, Errors = -0.163944, -0.107911 +DEAL::f(x,y) = x^1 y^3, Errors = -0.121118, -0.096111 +DEAL::f(x,y) = x^1 y^4, Errors = -0.094562, -0.082873 +DEAL::f(x,y) = x^1 y^5, Errors = -0.077252, -0.071600 +DEAL::f(x,y) = x^2 y^0, Errors = -0.332879, -0.133607 +DEAL::f(x,y) = x^2 y^1, Errors = -0.163944, -0.107911 +DEAL::f(x,y) = x^2 y^2, Errors = -0.105791, -0.088336 +DEAL::f(x,y) = x^2 y^3, Errors = -0.076850, -0.070501 +DEAL::f(x,y) = x^2 y^4, Errors = -0.060010, -0.057353 +DEAL::f(x,y) = x^2 y^5, Errors = -0.049130, -0.047911 +DEAL::f(x,y) = x^3 y^0, Errors = -0.262944, -0.155343 +DEAL::f(x,y) = x^3 y^1, Errors = -0.121118, -0.096111 +DEAL::f(x,y) = x^3 y^2, Errors = -0.076850, -0.070501 +DEAL::f(x,y) = x^3 y^3, Errors = -0.055777, -0.053901 +DEAL::f(x,y) = x^3 y^4, Errors = -0.043624, -0.042958 +DEAL::f(x,y) = x^3 y^5, Errors = -0.035766, -0.035490 +DEAL::f(x,y) = x^4 y^0, Errors = -0.206246, -0.152327 +DEAL::f(x,y) = x^4 y^1, Errors = -0.094562, -0.082873 +DEAL::f(x,y) = x^4 y^2, Errors = -0.060010, -0.057353 +DEAL::f(x,y) = x^4 y^3, Errors = -0.043624, -0.042958 +DEAL::f(x,y) = x^4 y^4, Errors = -0.034164, -0.033969 +DEAL::f(x,y) = x^4 y^5, Errors = -0.028037, -0.027968 +DEAL::f(x,y) = x^5 y^0, Errors = -0.167904, -0.141111 +DEAL::f(x,y) = x^5 y^1, Errors = -0.077252, -0.071600 +DEAL::f(x,y) = x^5 y^2, Errors = -0.049130, -0.047911 +DEAL::f(x,y) = x^5 y^3, Errors = -0.035766, -0.035490 +DEAL::f(x,y) = x^5 y^4, Errors = -0.028037, -0.027968 +DEAL::f(x,y) = x^5 y^5, Errors = -0.023024, -0.023004 +DEAL:: =========Quadrature Order: 2 =============================== +DEAL:: ============================================================ +DEAL:: ===============Vertex: 0.000000 0.000000 ============================= +DEAL::f(x,y) = x^0 y^0, Errors = 0.083005, -0.001795 +DEAL::f(x,y) = x^0 y^1, Errors = 0.034615, -0.002254 +DEAL::f(x,y) = x^0 y^2, Errors = -0.049939, -0.003643 +DEAL::f(x,y) = x^0 y^3, Errors = -0.105989, -0.004860 +DEAL::f(x,y) = x^0 y^4, Errors = -0.126977, -0.011516 +DEAL::f(x,y) = x^0 y^5, Errors = -0.128734, -0.019894 +DEAL::f(x,y) = x^1 y^0, Errors = 0.034615, -0.002254 +DEAL::f(x,y) = x^1 y^1, Errors = 0.027031, -0.002407 +DEAL::f(x,y) = x^1 y^2, Errors = -0.021731, -0.003410 +DEAL::f(x,y) = x^1 y^3, Errors = -0.049085, -0.007693 +DEAL::f(x,y) = x^1 y^4, Errors = -0.058941, -0.012832 +DEAL::f(x,y) = x^1 y^5, Errors = -0.059742, -0.017132 +DEAL::f(x,y) = x^2 y^0, Errors = -0.049939, -0.003643 +DEAL::f(x,y) = x^2 y^1, Errors = -0.021731, -0.003410 +DEAL::f(x,y) = x^2 y^2, Errors = -0.033894, -0.006861 +DEAL::f(x,y) = x^2 y^3, Errors = -0.041297, -0.010777 +DEAL::f(x,y) = x^2 y^4, Errors = -0.042534, -0.014116 +DEAL::f(x,y) = x^2 y^5, Errors = -0.040552, -0.016356 +DEAL::f(x,y) = x^3 y^0, Errors = -0.105989, -0.004860 +DEAL::f(x,y) = x^3 y^1, Errors = -0.049085, -0.007693 +DEAL::f(x,y) = x^3 y^2, Errors = -0.041297, -0.010777 +DEAL::f(x,y) = x^3 y^3, Errors = -0.038301, -0.013245 +DEAL::f(x,y) = x^3 y^4, Errors = -0.035046, -0.014991 +DEAL::f(x,y) = x^3 y^5, Errors = -0.031558, -0.015882 +DEAL::f(x,y) = x^4 y^0, Errors = -0.126977, -0.011516 +DEAL::f(x,y) = x^4 y^1, Errors = -0.058941, -0.012832 +DEAL::f(x,y) = x^4 y^2, Errors = -0.042534, -0.014116 +DEAL::f(x,y) = x^4 y^3, Errors = -0.035046, -0.014991 +DEAL::f(x,y) = x^4 y^4, Errors = -0.029956, -0.015454 +DEAL::f(x,y) = x^4 y^5, Errors = -0.025973, -0.015409 +DEAL::f(x,y) = x^5 y^0, Errors = -0.128734, -0.019894 +DEAL::f(x,y) = x^5 y^1, Errors = -0.059742, -0.017132 +DEAL::f(x,y) = x^5 y^2, Errors = -0.040552, -0.016356 +DEAL::f(x,y) = x^5 y^3, Errors = -0.031558, -0.015882 +DEAL::f(x,y) = x^5 y^4, Errors = -0.025973, -0.015409 +DEAL::f(x,y) = x^5 y^5, Errors = -0.022011, -0.014742 +DEAL:: =========Quadrature Order: 3 =============================== +DEAL:: ============================================================ +DEAL:: ===============Vertex: 0.000000 0.000000 ============================= +DEAL::f(x,y) = x^0 y^0, Errors = 0.049500, -0.000052 +DEAL::f(x,y) = x^0 y^1, Errors = -0.000455, -0.000092 +DEAL::f(x,y) = x^0 y^2, Errors = 0.000338, -0.000200 +DEAL::f(x,y) = x^0 y^3, Errors = -0.010300, -0.000369 +DEAL::f(x,y) = x^0 y^4, Errors = -0.027273, -0.000603 +DEAL::f(x,y) = x^0 y^5, Errors = -0.042569, -0.000866 +DEAL::f(x,y) = x^1 y^0, Errors = -0.000455, -0.000092 +DEAL::f(x,y) = x^1 y^1, Errors = 0.000150, -0.000105 +DEAL::f(x,y) = x^1 y^2, Errors = 0.001770, -0.000183 +DEAL::f(x,y) = x^1 y^3, Errors = -0.004425, -0.000322 +DEAL::f(x,y) = x^1 y^4, Errors = -0.012927, -0.000525 +DEAL::f(x,y) = x^1 y^5, Errors = -0.020210, -0.000953 +DEAL::f(x,y) = x^2 y^0, Errors = 0.000338, -0.000200 +DEAL::f(x,y) = x^2 y^1, Errors = 0.001770, -0.000183 +DEAL::f(x,y) = x^2 y^2, Errors = 0.002045, -0.000229 +DEAL::f(x,y) = x^2 y^3, Errors = -0.002509, -0.000337 +DEAL::f(x,y) = x^2 y^4, Errors = -0.008202, -0.000658 +DEAL::f(x,y) = x^2 y^5, Errors = -0.012954, -0.001202 +DEAL::f(x,y) = x^3 y^0, Errors = -0.010300, -0.000369 +DEAL::f(x,y) = x^3 y^1, Errors = -0.004425, -0.000322 +DEAL::f(x,y) = x^3 y^2, Errors = -0.002509, -0.000337 +DEAL::f(x,y) = x^3 y^3, Errors = -0.004742, -0.000553 +DEAL::f(x,y) = x^3 y^4, Errors = -0.008084, -0.000986 +DEAL::f(x,y) = x^3 y^5, Errors = -0.010953, -0.001572 +DEAL::f(x,y) = x^4 y^0, Errors = -0.027273, -0.000603 +DEAL::f(x,y) = x^4 y^1, Errors = -0.012927, -0.000525 +DEAL::f(x,y) = x^4 y^2, Errors = -0.008202, -0.000658 +DEAL::f(x,y) = x^4 y^3, Errors = -0.008084, -0.000986 +DEAL::f(x,y) = x^4 y^4, Errors = -0.009419, -0.001464 +DEAL::f(x,y) = x^4 y^5, Errors = -0.010761, -0.002024 +DEAL::f(x,y) = x^5 y^0, Errors = -0.042569, -0.000866 +DEAL::f(x,y) = x^5 y^1, Errors = -0.020210, -0.000953 +DEAL::f(x,y) = x^5 y^2, Errors = -0.012954, -0.001202 +DEAL::f(x,y) = x^5 y^3, Errors = -0.010953, -0.001572 +DEAL::f(x,y) = x^5 y^4, Errors = -0.010761, -0.002024 +DEAL::f(x,y) = x^5 y^5, Errors = -0.010963, -0.002508 +DEAL:: =========Quadrature Order: 4 =============================== +DEAL:: ============================================================ +DEAL:: ===============Vertex: 0.000000 0.000000 ============================= +DEAL::f(x,y) = x^0 y^0, Errors = 0.021072, -0.000002 +DEAL::f(x,y) = x^0 y^1, Errors = 0.000947, -0.000003 +DEAL::f(x,y) = x^0 y^2, Errors = -0.000061, -0.000009 +DEAL::f(x,y) = x^0 y^3, Errors = 0.000186, -0.000022 +DEAL::f(x,y) = x^0 y^4, Errors = -0.001562, -0.000045 +DEAL::f(x,y) = x^0 y^5, Errors = -0.005915, -0.000080 +DEAL::f(x,y) = x^1 y^0, Errors = 0.000947, -0.000003 +DEAL::f(x,y) = x^1 y^1, Errors = 0.000090, -0.000004 +DEAL::f(x,y) = x^1 y^2, Errors = -0.000191, -0.000008 +DEAL::f(x,y) = x^1 y^3, Errors = -0.000007, -0.000019 +DEAL::f(x,y) = x^1 y^4, Errors = -0.000944, -0.000038 +DEAL::f(x,y) = x^1 y^5, Errors = -0.003100, -0.000069 +DEAL::f(x,y) = x^2 y^0, Errors = -0.000061, -0.000009 +DEAL::f(x,y) = x^2 y^1, Errors = -0.000191, -0.000008 +DEAL::f(x,y) = x^2 y^2, Errors = -0.000083, -0.000011 +DEAL::f(x,y) = x^2 y^3, Errors = 0.000084, -0.000019 +DEAL::f(x,y) = x^2 y^4, Errors = -0.000575, -0.000036 +DEAL::f(x,y) = x^2 y^5, Errors = -0.002016, -0.000063 +DEAL::f(x,y) = x^3 y^0, Errors = 0.000186, -0.000022 +DEAL::f(x,y) = x^3 y^1, Errors = -0.000007, -0.000019 +DEAL::f(x,y) = x^3 y^2, Errors = 0.000084, -0.000019 +DEAL::f(x,y) = x^3 y^3, Errors = 0.000171, -0.000025 +DEAL::f(x,y) = x^3 y^4, Errors = -0.000360, -0.000039 +DEAL::f(x,y) = x^3 y^5, Errors = -0.001451, -0.000070 +DEAL::f(x,y) = x^4 y^0, Errors = -0.001562, -0.000045 +DEAL::f(x,y) = x^4 y^1, Errors = -0.000944, -0.000038 +DEAL::f(x,y) = x^4 y^2, Errors = -0.000575, -0.000036 +DEAL::f(x,y) = x^4 y^3, Errors = -0.000360, -0.000039 +DEAL::f(x,y) = x^4 y^4, Errors = -0.000688, -0.000057 +DEAL::f(x,y) = x^4 y^5, Errors = -0.001478, -0.000099 +DEAL::f(x,y) = x^5 y^0, Errors = -0.005915, -0.000080 +DEAL::f(x,y) = x^5 y^1, Errors = -0.003100, -0.000069 +DEAL::f(x,y) = x^5 y^2, Errors = -0.002016, -0.000063 +DEAL::f(x,y) = x^5 y^3, Errors = -0.001451, -0.000070 +DEAL::f(x,y) = x^5 y^4, Errors = -0.001478, -0.000099 +DEAL::f(x,y) = x^5 y^5, Errors = -0.001951, -0.000155 +DEAL:: =========Quadrature Order: 5 =============================== +DEAL:: ============================================================ +DEAL:: ===============Vertex: 0.000000 0.000000 ============================= +DEAL::f(x,y) = x^0 y^0, Errors = 0.009674, 0.000000 +DEAL::f(x,y) = x^0 y^1, Errors = 0.000228, 0.000000 +DEAL::f(x,y) = x^0 y^2, Errors = -0.000007, 0.000000 +DEAL::f(x,y) = x^0 y^3, Errors = -0.000038, -0.000001 +DEAL::f(x,y) = x^0 y^4, Errors = -0.000015, -0.000003 +DEAL::f(x,y) = x^0 y^5, Errors = -0.000326, -0.000006 +DEAL::f(x,y) = x^1 y^0, Errors = 0.000228, 0.000000 +DEAL::f(x,y) = x^1 y^1, Errors = 0.000081, 0.000000 +DEAL::f(x,y) = x^1 y^2, Errors = 0.000029, 0.000000 +DEAL::f(x,y) = x^1 y^3, Errors = -0.000001, -0.000001 +DEAL::f(x,y) = x^1 y^4, Errors = 0.000013, -0.000002 +DEAL::f(x,y) = x^1 y^5, Errors = -0.000154, -0.000005 +DEAL::f(x,y) = x^2 y^0, Errors = -0.000007, 0.000000 +DEAL::f(x,y) = x^2 y^1, Errors = 0.000029, 0.000000 +DEAL::f(x,y) = x^2 y^2, Errors = -0.000006, 0.000000 +DEAL::f(x,y) = x^2 y^3, Errors = -0.000019, -0.000001 +DEAL::f(x,y) = x^2 y^4, Errors = -0.000005, -0.000002 +DEAL::f(x,y) = x^2 y^5, Errors = -0.000118, -0.000005 +DEAL::f(x,y) = x^3 y^0, Errors = -0.000038, -0.000001 +DEAL::f(x,y) = x^3 y^1, Errors = -0.000001, -0.000001 +DEAL::f(x,y) = x^3 y^2, Errors = -0.000019, -0.000001 +DEAL::f(x,y) = x^3 y^3, Errors = -0.000017, -0.000001 +DEAL::f(x,y) = x^3 y^4, Errors = -0.000002, -0.000002 +DEAL::f(x,y) = x^3 y^5, Errors = -0.000087, -0.000004 +DEAL::f(x,y) = x^4 y^0, Errors = -0.000015, -0.000003 +DEAL::f(x,y) = x^4 y^1, Errors = 0.000013, -0.000002 +DEAL::f(x,y) = x^4 y^2, Errors = -0.000005, -0.000002 +DEAL::f(x,y) = x^4 y^3, Errors = -0.000002, -0.000002 +DEAL::f(x,y) = x^4 y^4, Errors = 0.000009, -0.000003 +DEAL::f(x,y) = x^4 y^5, Errors = -0.000062, -0.000005 +DEAL::f(x,y) = x^5 y^0, Errors = -0.000326, -0.000006 +DEAL::f(x,y) = x^5 y^1, Errors = -0.000154, -0.000005 +DEAL::f(x,y) = x^5 y^2, Errors = -0.000118, -0.000005 +DEAL::f(x,y) = x^5 y^3, Errors = -0.000087, -0.000004 +DEAL::f(x,y) = x^5 y^4, Errors = -0.000062, -0.000005 +DEAL::f(x,y) = x^5 y^5, Errors = -0.000110, -0.000007 +DEAL:: =========Quadrature Order: 6 =============================== +DEAL:: ============================================================ +DEAL:: ===============Vertex: 0.000000 0.000000 ============================= +DEAL::f(x,y) = x^0 y^0, Errors = 0.004879, 0.000000 +DEAL::f(x,y) = x^0 y^1, Errors = 0.000082, 0.000000 +DEAL::f(x,y) = x^0 y^2, Errors = -0.000001, 0.000000 +DEAL::f(x,y) = x^0 y^3, Errors = -0.000001, 0.000000 +DEAL::f(x,y) = x^0 y^4, Errors = -0.000004, 0.000000 +DEAL::f(x,y) = x^0 y^5, Errors = -0.000008, 0.000000 +DEAL::f(x,y) = x^1 y^0, Errors = 0.000082, 0.000000 +DEAL::f(x,y) = x^1 y^1, Errors = 0.000009, 0.000000 +DEAL::f(x,y) = x^1 y^2, Errors = 0.000001, 0.000000 +DEAL::f(x,y) = x^1 y^3, Errors = 0.000000, 0.000000 +DEAL::f(x,y) = x^1 y^4, Errors = -0.000003, 0.000000 +DEAL::f(x,y) = x^1 y^5, Errors = -0.000003, 0.000000 +DEAL::f(x,y) = x^2 y^0, Errors = -0.000001, 0.000000 +DEAL::f(x,y) = x^2 y^1, Errors = 0.000001, 0.000000 +DEAL::f(x,y) = x^2 y^2, Errors = 0.000004, 0.000000 +DEAL::f(x,y) = x^2 y^3, Errors = 0.000002, 0.000000 +DEAL::f(x,y) = x^2 y^4, Errors = 0.000000, 0.000000 +DEAL::f(x,y) = x^2 y^5, Errors = 0.000000, 0.000000 +DEAL::f(x,y) = x^3 y^0, Errors = -0.000001, 0.000000 +DEAL::f(x,y) = x^3 y^1, Errors = 0.000000, 0.000000 +DEAL::f(x,y) = x^3 y^2, Errors = 0.000002, 0.000000 +DEAL::f(x,y) = x^3 y^3, Errors = 0.000000, 0.000000 +DEAL::f(x,y) = x^3 y^4, Errors = -0.000002, 0.000000 +DEAL::f(x,y) = x^3 y^5, Errors = -0.000001, 0.000000 +DEAL::f(x,y) = x^4 y^0, Errors = -0.000004, 0.000000 +DEAL::f(x,y) = x^4 y^1, Errors = -0.000003, 0.000000 +DEAL::f(x,y) = x^4 y^2, Errors = 0.000000, 0.000000 +DEAL::f(x,y) = x^4 y^3, Errors = -0.000002, 0.000000 +DEAL::f(x,y) = x^4 y^4, Errors = -0.000002, 0.000000 +DEAL::f(x,y) = x^4 y^5, Errors = -0.000002, 0.000000 +DEAL::f(x,y) = x^5 y^0, Errors = -0.000008, 0.000000 +DEAL::f(x,y) = x^5 y^1, Errors = -0.000003, 0.000000 +DEAL::f(x,y) = x^5 y^2, Errors = 0.000000, 0.000000 +DEAL::f(x,y) = x^5 y^3, Errors = -0.000001, 0.000000 +DEAL::f(x,y) = x^5 y^4, Errors = -0.000002, 0.000000 +DEAL::f(x,y) = x^5 y^5, Errors = -0.000001, 0.000000 -- 2.39.5