From 1916ed2fbd756f10dd5452c8fc493ee11ead192d Mon Sep 17 00:00:00 2001 From: Nicola Giuliani Date: Thu, 21 Dec 2017 23:05:42 +0100 Subject: [PATCH] Added QDuffy --- doc/news/changes/minor/20171221LucaHeltai | 2 +- include/deal.II/base/quadrature_lib.h | 59 ++++++ source/base/quadrature_lib.cc | 34 ++++ tests/base/quadrature_simplex_08.cc | 98 ++++++++++ tests/base/quadrature_simplex_08.output | 222 ++++++++++++++++++++++ 5 files changed, 414 insertions(+), 1 deletion(-) create mode 100644 tests/base/quadrature_simplex_08.cc create mode 100644 tests/base/quadrature_simplex_08.output diff --git a/doc/news/changes/minor/20171221LucaHeltai b/doc/news/changes/minor/20171221LucaHeltai index 4415a09c93..b30f24b0c7 100644 --- a/doc/news/changes/minor/20171221LucaHeltai +++ b/doc/news/changes/minor/20171221LucaHeltai @@ -1,4 +1,4 @@ -New: Added QSimplex, QTrianglePolar, and QSplit classes to perform quadratures +New: Added QSimplex, QTrianglePolar, QDuffy, and QSplit classes to perform quadratures on reference simplices, on their affine transformations, and on hyper cubes split by a given point.
diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 017239139b..6ac9a23e3a 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -746,6 +746,65 @@ public: QTrianglePolar(const unsigned int &n); }; +/** + * A quadrature that implements the Duffy transformation from a square to a + * triangle to integrate singularities in the origin of the reference + * simplex. + * + * The Duffy transformation is defined as + * \f[ + * \begin{pmatrix} + * x\\ + * y + * \end{pmatrix} + * = + * \begin{pmatrix} + * \hat x^\beta (1-\hat y)\\ + * \hat x^\beta \hat y + * end{pmatrix} + * \f] + * with determinant of the Jacobian equal to $J= \beta \hat \x^{2\beta-1}$. + * Such transformation maps the reference square \$[0,1]\times[0,1]$ to the + * reference simplex, by collapsing the left \side of the square and + * squeezing quadrature points towards the orgin, and then shearing the + * resulting triangle to the reference one. This transformation, allows + * one to integrate singularities of order $1/R$ in the origin when $\beta = + * 1$, and higher when $1 < \beta \leq 2$. + * + * When $\beta = 1$, this transformation is also known as the Lachat-Watson + * transformation. + * + * @author Luca Heltai, Nicola Giuliani, 2017. + */ +class QDuffy: public QSimplex<2> +{ +public: + /** + * Constructor that allows the specificatino of different quadrature rules + * along the "radial" and "angular" directions. + * + * Since this quadrature is not based on a Polar change of coordinates, it + * is not fully proper to talk about radial and angular directions. However, + * the effect of the Duffy transformation is similar to a polar change + * of coordinates, since the resulting quadrature points are aligned radially + * with respect to the singularity. + * + * @param radial_quadrature Base quadrature to use in the radial direction + * @param angular_quadrature Base quadrature to use in the angular direction + */ + QDuffy(const Quadrature<1> &radial_quadrature, + const Quadrature<1> &angular_quadrature, + const double &beta = 1.0); + + /** + * Calls the above constructor with QGauss<1>(n) quadrature formulas for + * both the radial and angular quadratures. + * + * @param n + */ + QDuffy(const unsigned int &n, const double &beta); +}; + /** * A quadrature to use when the cell should be split in subregions to integrate * using one or more base quadratures. diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index c5c0202768..e53b2990ab 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -1434,6 +1434,40 @@ QTrianglePolar::QTrianglePolar(const unsigned int &n) +QDuffy::QDuffy(const Quadrature<1> &radial_quadrature, + const Quadrature<1> &angular_quadrature, + const double &beta) : + QSimplex<2>(Quadrature<2>()) +{ + QAnisotropic<2> base(radial_quadrature, angular_quadrature); + this->quadrature_points.resize(base.size()); + this->weights.resize(base.size()); + for (unsigned int i=0; iquadrature_points[i] = Point<2>(x,y); + this->weights[i] = w*J; + } +} + + + +QDuffy::QDuffy(const unsigned int &n, const double &beta) + :QDuffy(QGauss<1>(n), QGauss<1>(n), beta) +{} + + + template QSplit::QSplit(const QSimplex &base, const Point &split_point) diff --git a/tests/base/quadrature_simplex_08.cc b/tests/base/quadrature_simplex_08.cc new file mode 100644 index 0000000000..ce2d6eb45e --- /dev/null +++ b/tests/base/quadrature_simplex_08.cc @@ -0,0 +1,98 @@ +// --------------------------------------------------------------------- +// +// 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. +// Compare QTrianglePolar and QLachatWatson + +#include "../tests.h" +#include + +// all include files needed for the program +#include +#include +#include "simplex.h" + +#include + +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; + + const unsigned int max_order = 5; + + // m i j quadtype + double error[max_order][6][6][2] = {{{{0}}}}; + + for (unsigned int m=0; m(.5, .5); + + QSplit<2> quad(QTrianglePolar(m+1), split_point); + QSplit<2> quad_de(QDuffy(m+1, 1.0), 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; + double approx_integral_de = 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()); + } + + for (unsigned int q=0; q< quad_de.size(); ++q) + { + double x = quad_de.point(q)[0]; + double y = quad_de.point(q)[1]; + approx_integral_de += ( pow(x, (double)i) * + pow(y, (double)j) * + quad_de.weight(q) / + (quad_de.point(q)-split_point).norm()); + } + + error[m][i][j][0] = approx_integral - exact_integral; + error[m][i][j][1] = approx_integral_de - 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