Added QSplit

diff --git a/doc/doxygen/images/split_quadrature.png b/doc/doxygen/images/split_quadrature.png
new file mode 100644 (file)
index 0000000..f222d8e
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 a7854881859a335ed221284b8daac21a346e3600..4415a09c937139cfb67c4f23eaf9fec8dafe2023 100644 (file)
--- 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.
 <br> 
 (Luca Heltai, 2017/12/21)

diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h
index 5d213d901a06a4bab2fdf8c6bfeb3f45fde26a26..fe0a51bdbba568af799ab2a691649435bbdf2570 100644 (file)
--- 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<int dim>
+class QSplit : public Quadrature<dim>
+{
+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<dim> &base,
+         const Point<dim> &split_point);
+};
+
 /*@}*/
 
 /* -------------- declaration of explicit specializations ------------- */

diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc
index 9a3ae2c92a7776349c5f8316a23b1cb75e94b31b..967150fa27c61a922f7f25cdcb8977cdc8e24bb9 100644 (file)
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -1434,6 +1434,46 @@ QTrianglePolar::QTrianglePolar(const unsigned int &n)
 
 
 
+template<int dim>
+QSplit<dim>::QSplit(const QSimplex<dim> &base,
+                    const Point<dim> &split_point)
+{
+
+  std::array<Point<dim>, dim+1> vertices;
+  vertices[0] = split_point;
+  switch (dim)
+    {
+    case 3:
+      Assert(false, ExcNotImplemented());
+      break;
+    default:
+    {
+      for (unsigned int f=0; f< GeometryInfo<dim>::faces_per_cell; ++f)
+        {
+          for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i)
+            vertices[i+1] =
+              GeometryInfo<dim>::unit_cell_vertex(
+                GeometryInfo<dim>::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 (file)
index 0000000..de9e52f
--- /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 <deal.II/base/quadrature_lib.h>
+#include <numeric>
+#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 (file)
index 0000000..2d45ed5
--- /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