Add additional quadrature rule for Simplex

diff --git a/doc/news/changes/minor/20200702Munch b/doc/news/changes/minor/20200702Munch
new file mode 100644 (file)
index 0000000..b9d4d55
--- /dev/null
+++ b/doc/news/changes/minor/20200702Munch
@@ -0,0 +1,3 @@
+New: A new quadrature rule for simplex geometric entities has been added.
+<br>
+(Peter Munch, 2020/07/02)

diff --git a/include/deal.II/simplex/quadrature_lib.h b/include/deal.II/simplex/quadrature_lib.h
new file mode 100644 (file)
index 0000000..fbe4f88
--- /dev/null
+++ b/include/deal.II/simplex/quadrature_lib.h
@@ -0,0 +1,48 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_tet_quadrature_lib_h
+#define dealii_tet_quadrature_lib_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/quadrature_lib.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+  /**
+   * Integration rule for simplex entities.
+   *
+   * Following number of quadrature points are currently supported:
+   *   - 2D: 1, 3, 7
+   *   - 3D: 1, 4, 10
+   */
+  template <int dim>
+  class PGauss : public QSimplex<dim>
+  {
+  public:
+    /**
+     * Constructor taking the number of quadrature points @p n_points.
+     */
+    explicit PGauss(const unsigned int n_points);
+  };
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt
index 7173d25b85b730b911fa2218b57e3ca4da96e513..a6d85d8b68025eb3a16929cf0e20b7eb5dab8c50 100644 (file)
--- a/source/CMakeLists.txt
+++ b/source/CMakeLists.txt
@@ -71,6 +71,7 @@ ADD_SUBDIRECTORY(differentiation)
 ADD_SUBDIRECTORY(physics)
 ADD_SUBDIRECTORY(optimization/rol)
 ADD_SUBDIRECTORY(non_matching)
+ADD_SUBDIRECTORY(simplex)
 ADD_SUBDIRECTORY(sundials)
 
 FOREACH(build ${DEAL_II_BUILD_TYPES})

diff --git a/source/simplex/CMakeLists.txt b/source/simplex/CMakeLists.txt
new file mode 100644 (file)
index 0000000..cfc0183
--- /dev/null
+++ b/source/simplex/CMakeLists.txt
@@ -0,0 +1,42 @@
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2020 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.md at
+## the top level directory of deal.II.
+##
+## ---------------------------------------------------------------------
+
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
+
+SET(_unity_include_src
+  quadrature_lib.cc
+  )
+
+SET(_separate_src
+  )
+
+# include all files in the unity file
+SET(_n_includes_per_unity_file 20)
+
+SETUP_SOURCE_LIST("${_unity_include_src}"
+  "${_separate_src}"
+  ${_n_includes_per_unity_file}
+  _src
+  )
+
+SET(_inst
+)
+
+FILE(GLOB _header
+  ${CMAKE_SOURCE_DIR}/include/deal.II/tet/*.h
+  )
+
+DEAL_II_ADD_LIBRARY(obj_simplex OBJECT ${_src} ${_header} ${_inst})
+EXPAND_INSTANTIATIONS(obj_simplex "${_inst}")

diff --git a/source/simplex/quadrature_lib.cc b/source/simplex/quadrature_lib.cc
new file mode 100644 (file)
index 0000000..94bd3ad
--- /dev/null
+++ b/source/simplex/quadrature_lib.cc
@@ -0,0 +1,146 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/polynomial.h>
+
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include <algorithm>
+#include <cmath>
+#include <functional>
+#include <limits>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+  template <int dim>
+  PGauss<dim>::PGauss(const unsigned int n_points)
+    : QSimplex<dim>(Quadrature<dim>())
+  {
+    // fill quadrature points and quadrature weights
+
+    if (dim == 2)
+      {
+        if (n_points == 1)
+          {
+            const double p = 1.0 / 3.0;
+            this->quadrature_points.emplace_back(p, p);
+            this->weights.emplace_back(1.0);
+          }
+        else if (n_points == 3)
+          {
+            const double Q23 = 2.0 / 3.0;
+            const double Q16 = 1.0 / 6.0;
+
+            this->quadrature_points.emplace_back(Q23, Q16);
+            this->quadrature_points.emplace_back(Q16, Q23);
+            this->quadrature_points.emplace_back(Q16, Q16);
+            this->weights.emplace_back(Q16);
+            this->weights.emplace_back(Q16);
+            this->weights.emplace_back(Q16);
+          }
+        else if (n_points == 7)
+          {
+            const double q12 = 0.5;
+
+            // clang-format off
+            this->quadrature_points.emplace_back(0.3333333333330, 0.3333333333330);
+            this->quadrature_points.emplace_back(0.7974269853530, 0.1012865073230);
+            this->quadrature_points.emplace_back(0.1012865073230, 0.7974269853530);
+            this->quadrature_points.emplace_back(0.1012865073230, 0.1012865073230);
+            this->quadrature_points.emplace_back(0.0597158717898, 0.4701420641050);
+            this->quadrature_points.emplace_back(0.4701420641050, 0.0597158717898);
+            this->quadrature_points.emplace_back(0.4701420641050, 0.4701420641050);
+            // clang-format on
+
+            this->weights.emplace_back(q12 * 0.225);
+            this->weights.emplace_back(q12 * 0.125939180545);
+            this->weights.emplace_back(q12 * 0.125939180545);
+            this->weights.emplace_back(q12 * 0.125939180545);
+            this->weights.emplace_back(q12 * 0.132394152789);
+            this->weights.emplace_back(q12 * 0.132394152789);
+            this->weights.emplace_back(q12 * 0.132394152789);
+          }
+      }
+    else if (dim == 3)
+      {
+        if (n_points == 1)
+          {
+            const double Q14 = 1.0 / 4.0;
+            const double Q16 = 1.0 / 6.0;
+
+            this->quadrature_points.emplace_back(Q14, Q14, Q14);
+            this->weights.emplace_back(Q16);
+          }
+        else if (n_points == 4)
+          {
+            const double Q124 = 1.0 / 6.0 / 4.0;
+
+            const double palpha = (5.0 + 3.0 * sqrt(5.0)) / 20.0;
+            const double pbeta  = (5.0 - sqrt(5.0)) / 20.0;
+            this->quadrature_points.emplace_back(pbeta, pbeta, pbeta);
+            this->quadrature_points.emplace_back(palpha, pbeta, pbeta);
+            this->quadrature_points.emplace_back(pbeta, palpha, pbeta);
+            this->quadrature_points.emplace_back(pbeta, pbeta, palpha);
+            this->weights.emplace_back(Q124);
+            this->weights.emplace_back(Q124);
+            this->weights.emplace_back(Q124);
+            this->weights.emplace_back(Q124);
+          }
+        else if (n_points == 10)
+          {
+            const double Q16 = 1.0 / 6.0;
+
+            // clang-format off
+            this->quadrature_points.emplace_back(0.5684305841968444, 0.1438564719343852, 0.1438564719343852);
+            this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.1438564719343852);
+            this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.5684305841968444);
+            this->quadrature_points.emplace_back(0.1438564719343852, 0.5684305841968444, 0.1438564719343852);
+            this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.5000000000000000);
+            this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.5000000000000000);
+            this->quadrature_points.emplace_back(0.5000000000000000, 0.5000000000000000, 0.0000000000000000);
+            this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.0000000000000000);
+            this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.0000000000000000);
+            this->quadrature_points.emplace_back(0.0000000000000000, 0.0000000000000000, 0.5000000000000000);
+            // clang-format on
+
+            this->weights.emplace_back(0.2177650698804054 * Q16);
+            this->weights.emplace_back(0.2177650698804054 * Q16);
+            this->weights.emplace_back(0.2177650698804054 * Q16);
+            this->weights.emplace_back(0.2177650698804054 * Q16);
+            this->weights.emplace_back(0.0214899534130631 * Q16);
+            this->weights.emplace_back(0.0214899534130631 * Q16);
+            this->weights.emplace_back(0.0214899534130631 * Q16);
+            this->weights.emplace_back(0.0214899534130631 * Q16);
+            this->weights.emplace_back(0.0214899534130631 * Q16);
+            this->weights.emplace_back(0.0214899534130631 * Q16);
+          }
+      }
+
+    AssertDimension(this->quadrature_points.size(), this->weights.size());
+    Assert(this->quadrature_points.size() > 0,
+           ExcMessage("No valid quadrature points!"));
+  }
+
+} // namespace Simplex
+
+
+template class Simplex::PGauss<2>;
+template class Simplex::PGauss<3>;
+
+DEAL_II_NAMESPACE_CLOSE

diff --git a/tests/simplex/quadrature_lib_01.cc b/tests/simplex/quadrature_lib_01.cc
new file mode 100644 (file)
index 0000000..a59a34f
--- /dev/null
+++ b/tests/simplex/quadrature_lib_01.cc
@@ -0,0 +1,76 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test Simplex::PGauss: output its quadrature points and weights.
+
+
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test(const unsigned int n_points)
+{
+  Simplex::PGauss<dim> quad(n_points);
+
+  for (unsigned int q = 0; q < quad.size(); ++q)
+    {
+      deallog << quad.point(q) << " ";
+      deallog << quad.weight(q) << " ";
+      deallog << std::endl;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  {
+    deallog.push("2d-1");
+    test<2>(1 /*n_points*/);
+    deallog.pop();
+  }
+  {
+    deallog.push("2d-3");
+    test<2>(3);
+    deallog.pop();
+  }
+  {
+    deallog.push("2d-7");
+    test<2>(7);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("3d-1");
+    test<3>(1);
+    deallog.pop();
+  }
+  {
+    deallog.push("3d-4");
+    test<3>(4);
+    deallog.pop();
+  }
+  {
+    deallog.push("3d-10");
+    test<3>(10);
+    deallog.pop();
+  }
+}

diff --git a/tests/simplex/quadrature_lib_01.output b/tests/simplex/quadrature_lib_01.output
new file mode 100644 (file)
index 0000000..0325881
--- /dev/null
+++ b/tests/simplex/quadrature_lib_01.output
@@ -0,0 +1,27 @@
+
+DEAL:2d-1::0.333333 0.333333 1.00000 
+DEAL:2d-3::0.666667 0.166667 0.166667 
+DEAL:2d-3::0.166667 0.666667 0.166667 
+DEAL:2d-3::0.166667 0.166667 0.166667 
+DEAL:2d-7::0.333333 0.333333 0.112500 
+DEAL:2d-7::0.797427 0.101287 0.0629696 
+DEAL:2d-7::0.101287 0.797427 0.0629696 
+DEAL:2d-7::0.101287 0.101287 0.0629696 
+DEAL:2d-7::0.0597159 0.470142 0.0661971 
+DEAL:2d-7::0.470142 0.0597159 0.0661971 
+DEAL:2d-7::0.470142 0.470142 0.0661971 
+DEAL:3d-1::0.250000 0.250000 0.250000 0.166667 
+DEAL:3d-4::0.138197 0.138197 0.138197 0.0416667 
+DEAL:3d-4::0.585410 0.138197 0.138197 0.0416667 
+DEAL:3d-4::0.138197 0.585410 0.138197 0.0416667 
+DEAL:3d-4::0.138197 0.138197 0.585410 0.0416667 
+DEAL:3d-10::0.568431 0.143856 0.143856 0.0362942 
+DEAL:3d-10::0.143856 0.143856 0.143856 0.0362942 
+DEAL:3d-10::0.143856 0.143856 0.568431 0.0362942 
+DEAL:3d-10::0.143856 0.568431 0.143856 0.0362942 
+DEAL:3d-10::0.00000 0.500000 0.500000 0.00358166 
+DEAL:3d-10::0.500000 0.00000 0.500000 0.00358166 
+DEAL:3d-10::0.500000 0.500000 0.00000 0.00358166 
+DEAL:3d-10::0.500000 0.00000 0.00000 0.00358166 
+DEAL:3d-10::0.00000 0.500000 0.00000 0.00358166 
+DEAL:3d-10::0.00000 0.00000 0.500000 0.00358166