add tensorProduct quadrature rule and [0,1] to [-1,1] transformation

ikarus/utils/CMakeLists.txt

@@ -23,6 +23,7 @@ install(
         polyfit.hh
         pythonautodiffdefinitions.hh
         tensorutils.hh
+        tensorproductquadrule.hh
         traits.hh
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/utils
 )

ikarus/utils/tensorproductquadrule.hh

@@ -0,0 +1,35 @@
+// SPDX-FileCopyrightText: 2023 The Ikarus Developers [email protected]
+// SPDX-License-Identifier: LGPL-2.1-or-later
+#pragma once
+
+#include <dune/geometry/quadraturerules.hh>
+namespace Ikarus {
+  template <class BaseQuadrature, class Quadrature>
+  auto tensorProductQuadrature(const BaseQuadrature& baseQuad, const Quadrature& onedQuad) {
+    constexpr int baseQuadDim       = BaseQuadrature::d;
+    auto rule                       = Dune::QuadratureRule<double, baseQuadDim + 1>();
+    const unsigned int baseQuadSize = baseQuad.size();
+    for (unsigned int bqi = 0; bqi < baseQuadSize; ++bqi) {
+      const typename Dune::QuadraturePoint<double, baseQuadDim>::Vector& basePoint = baseQuad[bqi].position();
+      const double& baseWeight                                                     = baseQuad[bqi].weight();
+
+      typename Dune::QuadraturePoint<double, baseQuadDim + 1>::Vector point;
+      for (unsigned int i = 0; i < baseQuadDim; ++i)
+        point[i] = basePoint[i];
+
+      const unsigned int onedQuadSize = onedQuad.size();
+      for (unsigned int oqi = 0; oqi < onedQuadSize; ++oqi) {
+        point[baseQuadDim] = onedQuad[oqi].position()[0];
+        rule.emplace_back(Dune::QuadraturePoint(point, baseWeight * onedQuad[oqi].weight()));
+      }
+    }
+    return rule;
+  }
+
+  auto transformFromZeroOneToMinusOneToOne(const Dune::QuadraturePoint<double, 1>& gp) {
+    typename Dune::QuadraturePoint<double, 1>::Vector mappedCoordinate = gp.position();
+    mappedCoordinate[0]                                                = 2 * mappedCoordinate[0] - 1.0;
+    double mappedWeight                                                = gp.weight() * 2.0;
+    return Dune::QuadraturePoint<double, 1>(mappedCoordinate, mappedWeight);
+  }
+}  // namespace Ikarus

tests/src/testshellthicknessquadrule.cpp

@@ -0,0 +1,29 @@
+// SPDX-FileCopyrightText: 2023 The Ikarus Developers [email protected]
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+#include <config.h>
+
+#include <dune/common/test/testsuite.hh>
+using Dune::TestSuite;
+#include <dune/common/float_cmp.hh>
+
+#include <ikarus/utils/init.hh>
+#include <ikarus/utils/tensorproductquadrule.hh>
+
+int main(int argc, char** argv) {
+  Ikarus::init(argc, argv);
+  TestSuite t;
+
+  for (int i = 0; i < 5; ++i) {
+    const auto& twoDRule = Dune::QuadratureRules<double, 2>::rule(Dune::GeometryTypes::quadrilateral, i);
+    const auto& oneDRule = Dune::QuadratureRules<double, 1>::rule(Dune::GeometryTypes::line, 1);
+
+    const auto rule = Ikarus::tensorProductQuadrature(twoDRule, oneDRule);
+    double vol      = 0;
+    for (int gpIndex = 0; auto& gp : rule) {
+      vol += gp.weight();
+      gpIndex++;
+    }
+    t.check(Dune::FloatCmp::eq(vol, 1.0));
+  }
+}