QTrianglePolar

diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h
index 9e8bf23a43a1f40e9b31ad8a2e43521c511535d9..f61fab54180833949c6c988bd98a4aa6f8413709 100644 (file)
--- 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 9c0393b435e8183d16edaa6055faebe85e66ffff..a430d1f4e7818356c33f28e3384b2dfa0b9bdbd8 100644 (file)
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -1383,6 +1383,9 @@ QSimplex<dim>::affine_transformation(const std::array<Point<dim>, 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; i<this->size(); ++i)
     {
       qp[i] = Point<dim>(vertices[0]+B*this->point(i));
@@ -1393,6 +1396,42 @@ QSimplex<dim>::affine_transformation(const std::array<Point<dim>, 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; i<base.size(); ++i)
+    {
+      const auto &q = base.point(i);
+      const auto &w = base.weight(i);
+
+      const auto &xhat = q[0];
+      const auto &yhat = q[1];
+
+      const double t = numbers::PI_2*yhat;
+      const double &pi = numbers::PI;
+      const double st = std::sin(t);
+      const double ct = std::cos(t);
+      const double r = xhat/(st+ct);
+
+      const double J = pi*xhat/(2*(std::sin(pi*yhat) + 1));
+
+      this->quadrature_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 (file)
index 0000000..0db7e99
--- /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 <deal.II/base/quadrature_lib.h>
+
+#include <numeric>
+
+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<Point<dim>, 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 (file)
index 0000000..bd22a01
--- /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::