From 247784f1943b619e95d7e53f753187c80c30f453 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 21 Dec 2017 03:05:45 +0100 Subject: [PATCH] QTrianglePolar --- include/deal.II/base/quadrature_lib.h | 41 +++++ source/base/quadrature_lib.cc | 39 +++++ tests/base/quadrature_simplex_02.cc | 65 ++++++++ tests/base/quadrature_simplex_02.output | 207 ++++++++++++++++++++++++ 4 files changed, 352 insertions(+) create mode 100644 tests/base/quadrature_simplex_02.cc create mode 100644 tests/base/quadrature_simplex_02.output diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 9e8bf23a43..f61fab5418 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -684,6 +684,47 @@ public: }; +/** + * A two dimensional simplex quadrature with zero Jacobian on vertex zero. The + * quadrature is obtained through the following polar transformation: + * + * \f[ + * \begin{pmatrix} + * x \\ + * y + * \end{pmatrix} + * = + * \begin{pmatrix} + * \frac{\hat x}{\sin(\theta)+\cos(\theta)} cos(\theta) \\ + * \frac{\hat x}{\sin(\theta)+\cos(\theta)} sin(\theta) + * \end{pmatrix} + * \qquad \theta := \frac\pi 2 \hat y + * \f] + * + * @author Luca Heltai, 2017 + */ +class QTrianglePolar: public QSimplex<2> +{ +public: + /** + * Construct a QTrianglePolar quadrature, with different formulas in the + * radial and angular directions. + * + * @param radial_quadrature Radial quadrature + * @param angular_quadrature Angular quadrature + */ + QTrianglePolar(const Quadrature<1> &radial_quadrature, + const Quadrature<1> &angular_quadrature); + + /** + * Call the other constructor, with QGauss<1>(n) for both radial and + * angular quadrature. + * + * @param n Order of QGauss quadrature + */ + QTrianglePolar(const unsigned int &n); +}; + /*@}*/ /* -------------- declaration of explicit specializations ------------- */ diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index 9c0393b435..a430d1f4e7 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -1383,6 +1383,9 @@ QSimplex::affine_transformation(const std::array, dim+1>& vertic B = transpose(B); auto J = determinant(B); + AssertThrow(J>0, ExcMessage("The provided vertices generate an inverted Simplex." + " Make sure the vertices are oriented lexicographically.")) + for (unsigned int i=0; isize(); ++i) { qp[i] = Point(vertices[0]+B*this->point(i)); @@ -1393,6 +1396,42 @@ QSimplex::affine_transformation(const std::array, dim+1>& vertic +QTrianglePolar::QTrianglePolar(const Quadrature<1> &radial_quadrature, + const Quadrature<1> &angular_quadrature) : + 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>(r*ct, r*st); + this->weights[i] = w*J; + } +} + + + +QTrianglePolar::QTrianglePolar(const unsigned int &n) + :QTrianglePolar(QGauss<1>(n),QGauss<1>(n)) +{} + + + // explicit specialization // note that 1d formulae are specialized by implementation above template class QGauss<2>; diff --git a/tests/base/quadrature_simplex_02.cc b/tests/base/quadrature_simplex_02.cc new file mode 100644 index 0000000000..0db7e99461 --- /dev/null +++ b/tests/base/quadrature_simplex_02.cc @@ -0,0 +1,65 @@ +// --------------------------------------------------------------------- +// +// 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 a simplex simplex quadrature, and check that we can +// get an affine transformation out of it. + + +#include "../tests.h" + +#include + +#include + +using namespace dealii; + +void test(int n) +{ + const unsigned int dim = 2; + + QTrianglePolar quad(n); + + 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; + + std::array, dim+1> v = {{ + Point<2>(2,0), + Point<2>(5,4), + Point<2>(0,2) + } + }; + + auto quad2 = quad.affine_transformation(v); + + for (auto p:quad2.get_points()) + deallog << p << std::endl; + + + deallog << std::endl << "# Area 2: " << std::accumulate(quad2.get_weights().begin(), + quad2.get_weights().end(), + 0.0) << std::endl << std::endl; + +} + + +int main() +{ + initlog(); + test(10); +} diff --git a/tests/base/quadrature_simplex_02.output b/tests/base/quadrature_simplex_02.output new file mode 100644 index 0000000000..bd22a01c53 --- /dev/null +++ b/tests/base/quadrature_simplex_02.output @@ -0,0 +1,207 @@ + +DEAL::0.0127847 0.000262043 +DEAL::0.0661132 0.00135510 +DEAL::0.157076 0.00321952 +DEAL::0.277612 0.00569012 +DEAL::0.417015 0.00854741 +DEAL::0.562900 0.0115375 +DEAL::0.702303 0.0143948 +DEAL::0.822839 0.0168654 +DEAL::0.913802 0.0187299 +DEAL::0.967130 0.0198229 +DEAL::0.0117923 0.00125444 +DEAL::0.0609813 0.00648704 +DEAL::0.144883 0.0154123 +DEAL::0.256063 0.0272393 +DEAL::0.384645 0.0409176 +DEAL::0.519205 0.0552318 +DEAL::0.647788 0.0689101 +DEAL::0.758968 0.0807371 +DEAL::0.842869 0.0896624 +DEAL::0.892058 0.0948950 +DEAL::0.0103772 0.00266954 +DEAL::0.0536634 0.0138049 +DEAL::0.127497 0.0327986 +DEAL::0.225335 0.0579675 +DEAL::0.338487 0.0870760 +DEAL::0.456899 0.117538 +DEAL::0.570052 0.146646 +DEAL::0.667890 0.171815 +DEAL::0.741723 0.190809 +DEAL::0.785009 0.201944 +DEAL::0.00883377 0.00421296 +DEAL::0.0456819 0.0217864 +DEAL::0.108534 0.0517615 +DEAL::0.191820 0.0914821 +DEAL::0.288143 0.137420 +DEAL::0.388944 0.185493 +DEAL::0.485266 0.231431 +DEAL::0.568553 0.271152 +DEAL::0.631405 0.301127 +DEAL::0.668253 0.318700 +DEAL::0.00728961 0.00575712 +DEAL::0.0376966 0.0297717 +DEAL::0.0895619 0.0707334 +DEAL::0.158290 0.125013 +DEAL::0.237775 0.187788 +DEAL::0.320956 0.253481 +DEAL::0.400441 0.316257 +DEAL::0.469169 0.370536 +DEAL::0.521034 0.411498 +DEAL::0.551441 0.435512 +DEAL::0.00575712 0.00728961 +DEAL::0.0297717 0.0376966 +DEAL::0.0707334 0.0895619 +DEAL::0.125013 0.158290 +DEAL::0.187788 0.237775 +DEAL::0.253481 0.320956 +DEAL::0.316257 0.400441 +DEAL::0.370536 0.469169 +DEAL::0.411498 0.521034 +DEAL::0.435512 0.551441 +DEAL::0.00421296 0.00883377 +DEAL::0.0217864 0.0456819 +DEAL::0.0517615 0.108534 +DEAL::0.0914821 0.191820 +DEAL::0.137420 0.288143 +DEAL::0.185493 0.388944 +DEAL::0.231431 0.485266 +DEAL::0.271152 0.568553 +DEAL::0.301127 0.631405 +DEAL::0.318700 0.668253 +DEAL::0.00266954 0.0103772 +DEAL::0.0138049 0.0536634 +DEAL::0.0327986 0.127497 +DEAL::0.0579675 0.225335 +DEAL::0.0870760 0.338487 +DEAL::0.117538 0.456899 +DEAL::0.146646 0.570052 +DEAL::0.171815 0.667890 +DEAL::0.190809 0.741723 +DEAL::0.201944 0.785009 +DEAL::0.00125444 0.0117923 +DEAL::0.00648704 0.0609813 +DEAL::0.0154123 0.144883 +DEAL::0.0272393 0.256063 +DEAL::0.0409176 0.384645 +DEAL::0.0552318 0.519205 +DEAL::0.0689101 0.647788 +DEAL::0.0807371 0.758968 +DEAL::0.0896624 0.842869 +DEAL::0.0948950 0.892058 +DEAL::0.000262043 0.0127847 +DEAL::0.00135510 0.0661132 +DEAL::0.00321952 0.157076 +DEAL::0.00569012 0.277612 +DEAL::0.00854741 0.417015 +DEAL::0.0115375 0.562900 +DEAL::0.0143948 0.702303 +DEAL::0.0168654 0.822839 +DEAL::0.0187299 0.913802 +DEAL::0.0198229 0.967130 +DEAL:: +DEAL::# Area: 0.500000 +DEAL:: +DEAL::2.03783 0.0516629 +DEAL::2.19563 0.267163 +DEAL::2.46479 0.634742 +DEAL::2.82146 1.12183 +DEAL::3.23395 1.68516 +DEAL::3.66562 2.27467 +DEAL::4.07812 2.83800 +DEAL::4.43479 3.32509 +DEAL::4.70395 3.69267 +DEAL::4.86175 3.90817 +DEAL::2.03287 0.0496781 +DEAL::2.16997 0.256899 +DEAL::2.40382 0.610356 +DEAL::2.71371 1.07873 +DEAL::3.07210 1.62042 +DEAL::3.44715 2.18729 +DEAL::3.80554 2.72897 +DEAL::4.11543 3.19734 +DEAL::4.34928 3.55080 +DEAL::4.48638 3.75802 +DEAL::2.02579 0.0468479 +DEAL::2.13338 0.242263 +DEAL::2.31689 0.575584 +DEAL::2.56007 1.01727 +DEAL::2.84131 1.52810 +DEAL::3.13562 2.06267 +DEAL::3.41686 2.57350 +DEAL::3.66004 3.01519 +DEAL::3.84355 3.34851 +DEAL::3.95114 3.54392 +DEAL::2.01808 0.0437610 +DEAL::2.09347 0.226300 +DEAL::2.22208 0.537658 +DEAL::2.39250 0.950245 +DEAL::2.58959 1.42741 +DEAL::2.79584 1.92676 +DEAL::2.99294 2.40393 +DEAL::3.16336 2.81652 +DEAL::3.29196 3.12787 +DEAL::3.36736 3.31041 +DEAL::2.01035 0.0406727 +DEAL::2.05355 0.210330 +DEAL::2.12722 0.499714 +DEAL::2.22484 0.883184 +DEAL::2.33775 1.32668 +DEAL::2.45590 1.79079 +DEAL::2.56881 2.23428 +DEAL::2.66643 2.61775 +DEAL::2.74011 2.90713 +DEAL::2.78330 3.07679 +DEAL::2.00269 0.0376077 +DEAL::2.01392 0.194480 +DEAL::2.03308 0.462057 +DEAL::2.05846 0.816630 +DEAL::2.08781 1.22670 +DEAL::2.11853 1.65584 +DEAL::2.14789 2.06591 +DEAL::2.17327 2.42048 +DEAL::2.19242 2.68806 +DEAL::2.20365 2.84493 +DEAL::1.99497 0.0345194 +DEAL::1.97400 0.178509 +DEAL::1.93822 0.424113 +DEAL::1.89081 0.749569 +DEAL::1.83597 1.12597 +DEAL::1.77859 1.51986 +DEAL::1.72376 1.89626 +DEAL::1.67635 2.22171 +DEAL::1.64057 2.46732 +DEAL::1.61960 2.61131 +DEAL::1.98725 0.0314326 +DEAL::1.93409 0.162547 +DEAL::1.84340 0.386188 +DEAL::1.72323 0.682540 +DEAL::1.58425 1.02528 +DEAL::1.43881 1.38395 +DEAL::1.29984 1.72669 +DEAL::1.17967 2.02304 +DEAL::1.08898 2.24668 +DEAL::1.03581 2.37779 +DEAL::1.98018 0.0286023 +DEAL::1.89750 0.147911 +DEAL::1.75647 0.351415 +DEAL::1.56959 0.621083 +DEAL::1.35346 0.932961 +DEAL::1.12728 1.25934 +DEAL::0.911155 1.57122 +DEAL::0.724276 1.84088 +DEAL::0.583248 2.04439 +DEAL::0.500568 2.16370 +DEAL::1.97522 0.0266176 +DEAL::1.87184 0.137647 +DEAL::1.69551 0.327029 +DEAL::1.46185 0.577985 +DEAL::1.19161 0.868220 +DEAL::0.908813 1.17195 +DEAL::0.638579 1.46219 +DEAL::0.404918 1.71314 +DEAL::0.228586 1.90252 +DEAL::0.125208 2.01355 +DEAL:: +DEAL::# Area 2: 7.00000 +DEAL:: -- 2.39.5