From: Luca Heltai Date: Thu, 21 Dec 2017 12:12:07 +0000 (+0100) Subject: Added QSplit X-Git-Tag: v9.0.0-rc1~583^2~9 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f5c8c0f5b802f7b5a3f411801df3deb4251212f0;p=dealii.git Added QSplit --- diff --git a/doc/doxygen/images/split_quadrature.png b/doc/doxygen/images/split_quadrature.png new file mode 100644 index 0000000000..f222d8ebc6 Binary files /dev/null and b/doc/doxygen/images/split_quadrature.png differ diff --git a/doc/news/changes/minor/20171221LucaHeltai b/doc/news/changes/minor/20171221LucaHeltai index a785488185..4415a09c93 100644 --- a/doc/news/changes/minor/20171221LucaHeltai +++ b/doc/news/changes/minor/20171221LucaHeltai @@ -1,3 +1,5 @@ -New: Added QSimplex class to perform quadratures on reference simplices, and on their affine transformations. +New: Added QSimplex, QTrianglePolar, and QSplit classes to perform quadratures +on reference simplices, on their affine transformations, and on hyper cubes +split by a given point.
(Luca Heltai, 2017/12/21) diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 5d213d901a..fe0a51bdbb 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -746,6 +746,43 @@ public: QTrianglePolar(const unsigned int &n); }; +/** + * A quadrature to use when the cell should be split in subregions to integrate + * using one or more base quadratures. + * + * @author Luca Heltai, 2017. + */ +template +class QSplit : public Quadrature +{ +public: + /** + * Construct a quadrature formula by splitting the reference hyper cube into + * the minimum number of simplices that have vertex zero coinciding with + * `split_point`, and patch together affine transformations of the `base` + * quadrature. + * + * The resulting quadrature can be used, for example, to integrate functions + * with integrable singularities at the split point, provided that you select + * as base quadrature one that can integrate singular points on vertex zero + * of the reference simplex. + * + * An example usage in dimension two is given by: + * @code + * const unsigned int order = 5; + * QSplit<2> quad(QTrianglePolar(order), Point<2>(.3,.4)); + * @endcode + * + * The resulting quadrature will look like the following: + * @image html split_quadrature.png "" + * + * @param base Base QSimplex quadrature to use + * @param split_point Where to split the hyper cube + */ + QSplit(const QSimplex &base, + const Point &split_point); +}; + /*@}*/ /* -------------- declaration of explicit specializations ------------- */ diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index 9a3ae2c92a..967150fa27 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -1434,6 +1434,46 @@ QTrianglePolar::QTrianglePolar(const unsigned int &n) +template +QSplit::QSplit(const QSimplex &base, + const Point &split_point) +{ + + std::array, dim+1> vertices; + vertices[0] = split_point; + switch (dim) + { + case 3: + Assert(false, ExcNotImplemented()); + break; + default: + { + for (unsigned int f=0; f< GeometryInfo::faces_per_cell; ++f) + { + for (unsigned int i=0; i::vertices_per_face; ++i) + vertices[i+1] = + GeometryInfo::unit_cell_vertex( + GeometryInfo::face_to_cell_vertices(f, i) + ); + auto quad = base.compute_affine_transformation(vertices); + if (quad.size()) + { + this->quadrature_points.insert( + this->quadrature_points.end(), + quad.get_points().begin(), + quad.get_points().end()); + this->weights.insert( + this->weights.end(), + quad.get_weights().begin(), + quad.get_weights().end()); + } + } + } + } +} + + + // explicit specialization // note that 1d formulae are specialized by implementation above template class QGauss<2>; @@ -1476,4 +1516,8 @@ template class QSimplex<1>; template class QSimplex<2>; template class QSimplex<3>; +template class QSplit<1>; +template class QSplit<2>; +template class QSplit<3>; + DEAL_II_NAMESPACE_CLOSE diff --git a/tests/base/quadrature_simplex_04.cc b/tests/base/quadrature_simplex_04.cc new file mode 100644 index 0000000000..de9e52fd9f --- /dev/null +++ b/tests/base/quadrature_simplex_04.cc @@ -0,0 +1,81 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 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. +// +// --------------------------------------------------------------------- + +// construct an anisotropic simplex quadrature, and check that we can +// get an affine transformation out of it. + + +#include "../tests.h" +#include +#include +#include "simplex.h" + +void test(int n) +{ + const unsigned int dim = 2; + + { + QSplit<2> quad(QTrianglePolar(n), Point<2>(.3,.4)); + + for (auto p:quad.get_points()) + deallog << p << std::endl; + + deallog << std::endl << "# Area: " << std::accumulate(quad.get_weights().begin(), + quad.get_weights().end(), + 0.0) << std::endl << std::endl; + + if (quad.size() == n*n*4) + deallog << "# Size OK" << std::endl; + else + deallog << "# Size NOT OK" << std::endl; + + } + { + QSplit<2> quad(QTrianglePolar(n), Point<2>(0,.2)); + + for (auto p:quad.get_points()) + deallog << p << std::endl; + + deallog << std::endl << "# Area: " << std::accumulate(quad.get_weights().begin(), + quad.get_weights().end(), + 0.0) << std::endl << std::endl; + if (quad.size() == n*n*3) + deallog << "# Size OK" << std::endl; + else + deallog << "# Size NOT OK" << std::endl; + + } + { + QSplit<2> quad(QTrianglePolar(n), Point<2>(1,0)); + + for (auto p:quad.get_points()) + deallog << p << std::endl; + + deallog << std::endl << "# Area: " << std::accumulate(quad.get_weights().begin(), + quad.get_weights().end(), + 0.0) << std::endl << std::endl; + if (quad.size() == n*n*2) + deallog << "# Size OK" << std::endl; + else + deallog << "# Size NOT OK" << std::endl; + } +} + + +int main() +{ + initlog(); + test(5); +} diff --git a/tests/base/quadrature_simplex_04.output b/tests/base/quadrature_simplex_04.output new file mode 100644 index 0000000000..2d45ed5616 --- /dev/null +++ b/tests/base/quadrature_simplex_04.output @@ -0,0 +1,238 @@ + +DEAL::0.285927 0.384461 +DEAL::0.230770 0.323558 +DEAL::0.150000 0.234373 +DEAL::0.0692296 0.145187 +DEAL::0.0140730 0.0842843 +DEAL::0.285927 0.394135 +DEAL::0.230770 0.371146 +DEAL::0.150000 0.337482 +DEAL::0.0692296 0.303818 +DEAL::0.0140730 0.280830 +DEAL::0.285927 0.404691 +DEAL::0.230770 0.423077 +DEAL::0.150000 0.450000 +DEAL::0.0692296 0.476923 +DEAL::0.0140730 0.495309 +DEAL::0.285927 0.415247 +DEAL::0.230770 0.475007 +DEAL::0.150000 0.562518 +DEAL::0.0692296 0.650028 +DEAL::0.0140730 0.709788 +DEAL::0.285927 0.424921 +DEAL::0.230770 0.522595 +DEAL::0.150000 0.665627 +DEAL::0.0692296 0.808660 +DEAL::0.0140730 0.906334 +DEAL::0.332837 0.384461 +DEAL::0.461536 0.323558 +DEAL::0.650000 0.234373 +DEAL::0.838464 0.145187 +DEAL::0.967163 0.0842843 +DEAL::0.332837 0.394135 +DEAL::0.461536 0.371146 +DEAL::0.650000 0.337482 +DEAL::0.838464 0.303818 +DEAL::0.967163 0.280830 +DEAL::0.332837 0.404691 +DEAL::0.461536 0.423077 +DEAL::0.650000 0.450000 +DEAL::0.838464 0.476923 +DEAL::0.967163 0.495309 +DEAL::0.332837 0.415247 +DEAL::0.461536 0.475007 +DEAL::0.650000 0.562518 +DEAL::0.838464 0.650028 +DEAL::0.967163 0.709788 +DEAL::0.332837 0.424921 +DEAL::0.461536 0.522595 +DEAL::0.650000 0.665627 +DEAL::0.838464 0.808660 +DEAL::0.967163 0.906334 +DEAL::0.289152 0.381236 +DEAL::0.246634 0.307694 +DEAL::0.184373 0.200000 +DEAL::0.122111 0.0923061 +DEAL::0.0795933 0.0187640 +DEAL::0.298826 0.381236 +DEAL::0.294223 0.307694 +DEAL::0.287482 0.200000 +DEAL::0.280742 0.0923061 +DEAL::0.276139 0.0187640 +DEAL::0.309382 0.381236 +DEAL::0.346153 0.307694 +DEAL::0.400000 0.200000 +DEAL::0.453847 0.0923061 +DEAL::0.490618 0.0187640 +DEAL::0.319938 0.381236 +DEAL::0.398083 0.307694 +DEAL::0.512518 0.200000 +DEAL::0.626952 0.0923061 +DEAL::0.705097 0.0187640 +DEAL::0.329612 0.381236 +DEAL::0.445672 0.307694 +DEAL::0.615627 0.200000 +DEAL::0.785583 0.0923061 +DEAL::0.901643 0.0187640 +DEAL::0.289152 0.428146 +DEAL::0.246634 0.538459 +DEAL::0.184373 0.700000 +DEAL::0.122111 0.861541 +DEAL::0.0795933 0.971854 +DEAL::0.298826 0.428146 +DEAL::0.294223 0.538459 +DEAL::0.287482 0.700000 +DEAL::0.280742 0.861541 +DEAL::0.276139 0.971854 +DEAL::0.309382 0.428146 +DEAL::0.346153 0.538459 +DEAL::0.400000 0.700000 +DEAL::0.453847 0.861541 +DEAL::0.490618 0.971854 +DEAL::0.319938 0.428146 +DEAL::0.398083 0.538459 +DEAL::0.512518 0.700000 +DEAL::0.626952 0.861541 +DEAL::0.705097 0.971854 +DEAL::0.329612 0.428146 +DEAL::0.445672 0.538459 +DEAL::0.615627 0.700000 +DEAL::0.785583 0.861541 +DEAL::0.901643 0.971854 +DEAL:: +DEAL::# Area: 0.999975 +DEAL:: +DEAL::# Size OK +DEAL::0.0469101 0.193843 +DEAL::0.230765 0.169711 +DEAL::0.500000 0.134373 +DEAL::0.769235 0.0990342 +DEAL::0.953090 0.0749023 +DEAL::0.0469101 0.203517 +DEAL::0.230765 0.217299 +DEAL::0.500000 0.237482 +DEAL::0.769235 0.257665 +DEAL::0.953090 0.271448 +DEAL::0.0469101 0.214073 +DEAL::0.230765 0.269230 +DEAL::0.500000 0.350000 +DEAL::0.769235 0.430770 +DEAL::0.953090 0.485927 +DEAL::0.0469101 0.224629 +DEAL::0.230765 0.321160 +DEAL::0.500000 0.462518 +DEAL::0.769235 0.603875 +DEAL::0.953090 0.700406 +DEAL::0.0469101 0.234303 +DEAL::0.230765 0.368748 +DEAL::0.500000 0.565627 +DEAL::0.769235 0.762507 +DEAL::0.953090 0.896952 +DEAL::0.00322484 0.190618 +DEAL::0.0158640 0.153847 +DEAL::0.0343725 0.100000 +DEAL::0.0528811 0.0461531 +DEAL::0.0655202 0.00938202 +DEAL::0.0128986 0.190618 +DEAL::0.0634523 0.153847 +DEAL::0.137482 0.100000 +DEAL::0.211512 0.0461531 +DEAL::0.262066 0.00938202 +DEAL::0.0234550 0.190618 +DEAL::0.115383 0.153847 +DEAL::0.250000 0.100000 +DEAL::0.384617 0.0461531 +DEAL::0.476545 0.00938202 +DEAL::0.0340115 0.190618 +DEAL::0.167313 0.153847 +DEAL::0.362518 0.100000 +DEAL::0.557722 0.0461531 +DEAL::0.691024 0.00938202 +DEAL::0.0436852 0.190618 +DEAL::0.214901 0.153847 +DEAL::0.465627 0.100000 +DEAL::0.716354 0.0461531 +DEAL::0.887570 0.00938202 +DEAL::0.00322484 0.237528 +DEAL::0.0158640 0.384612 +DEAL::0.0343725 0.600000 +DEAL::0.0528811 0.815388 +DEAL::0.0655202 0.962472 +DEAL::0.0128986 0.237528 +DEAL::0.0634523 0.384612 +DEAL::0.137482 0.600000 +DEAL::0.211512 0.815388 +DEAL::0.262066 0.962472 +DEAL::0.0234550 0.237528 +DEAL::0.115383 0.384612 +DEAL::0.250000 0.600000 +DEAL::0.384617 0.815388 +DEAL::0.476545 0.962472 +DEAL::0.0340115 0.237528 +DEAL::0.167313 0.384612 +DEAL::0.362518 0.600000 +DEAL::0.557722 0.815388 +DEAL::0.691024 0.962472 +DEAL::0.0436852 0.237528 +DEAL::0.214901 0.384612 +DEAL::0.465627 0.600000 +DEAL::0.716354 0.815388 +DEAL::0.887570 0.962472 +DEAL:: +DEAL::# Area: 0.999975 +DEAL:: +DEAL::# Size OK +DEAL::0.953090 0.00322484 +DEAL::0.769235 0.0158640 +DEAL::0.500000 0.0343725 +DEAL::0.230765 0.0528811 +DEAL::0.0469101 0.0655202 +DEAL::0.953090 0.0128986 +DEAL::0.769235 0.0634523 +DEAL::0.500000 0.137482 +DEAL::0.230765 0.211512 +DEAL::0.0469101 0.262066 +DEAL::0.953090 0.0234550 +DEAL::0.769235 0.115383 +DEAL::0.500000 0.250000 +DEAL::0.230765 0.384617 +DEAL::0.0469101 0.476545 +DEAL::0.953090 0.0340115 +DEAL::0.769235 0.167313 +DEAL::0.500000 0.362518 +DEAL::0.230765 0.557722 +DEAL::0.0469101 0.691024 +DEAL::0.953090 0.0436852 +DEAL::0.769235 0.214901 +DEAL::0.500000 0.465627 +DEAL::0.230765 0.716354 +DEAL::0.0469101 0.887570 +DEAL::0.956315 0.0469101 +DEAL::0.785099 0.230765 +DEAL::0.534373 0.500000 +DEAL::0.283646 0.769235 +DEAL::0.112430 0.953090 +DEAL::0.965989 0.0469101 +DEAL::0.832687 0.230765 +DEAL::0.637482 0.500000 +DEAL::0.442278 0.769235 +DEAL::0.308976 0.953090 +DEAL::0.976545 0.0469101 +DEAL::0.884617 0.230765 +DEAL::0.750000 0.500000 +DEAL::0.615383 0.769235 +DEAL::0.523455 0.953090 +DEAL::0.987101 0.0469101 +DEAL::0.936548 0.230765 +DEAL::0.862518 0.500000 +DEAL::0.788488 0.769235 +DEAL::0.737934 0.953090 +DEAL::0.996775 0.0469101 +DEAL::0.984136 0.230765 +DEAL::0.965627 0.500000 +DEAL::0.947119 0.769235 +DEAL::0.934480 0.953090 +DEAL:: +DEAL::# Area: 0.999975 +DEAL:: +DEAL::# Size OK