Draft: Refactor Gauss points and move to Intrepid2

src/core/fem/src/general/utils/4C_fem_general_utils_gausspoints.hpp

@@ -180,18 +180,15 @@ namespace Core::FE
     std::vector<std::shared_ptr<GaussPoints>> gp_;
   };
 
-  /// remember calculated gauss points so we do not need to calculate again
-  class GaussPointCache
-  {
-   public:
-    static GaussPointCache& instance();
-
-    std::shared_ptr<GaussPoints> create(Core::FE::CellType distype, int degree);
+  /**
+   * Create GaussPoints object for a given element type and degree.
+   */
+  std::shared_ptr<GaussPoints> create_gauss_points(Core::FE::CellType distype, int degree);
 
-   private:
-    /// cache of already created gauss rules
-    std::map<std::pair<Core::FE::CellType, int>, std::shared_ptr<GaussPoints>> gp_cache_;
-  };
+  /**
+   * Create GaussPoints object for a given element type with a default degree.
+   */
+  std::shared_ptr<GaussPoints> create_gauss_points_default(Core::FE::CellType distype);
 
   /// gauss integration interface
   class GaussIntegration
@@ -243,7 +240,7 @@ namespace Core::FE
     typedef GaussPointIterator const_iterator;
 
     /// construct the optimal (normal) rule for a given element shape
-    GaussIntegration(Core::FE::CellType distype);
+    explicit GaussIntegration(Core::FE::CellType distype);
 
     /// construct rule for a given element shape
     GaussIntegration(Core::FE::CellType distype, int degree);
@@ -286,9 +283,8 @@ namespace Core::FE
         const Core::LinAlg::Matrix<Core::FE::dim<distype>, Core::FE::num_nodes<distype>>& xie,
         int degree)
     {
-      std::shared_ptr<GaussPoints> gp = GaussPointCache::instance().create(distype, degree);
-      std::shared_ptr<CollectedGaussPoints> cgp =
-          std::make_shared<CollectedGaussPoints>(gp->num_points());
+      std::shared_ptr<GaussPoints> gp = create_gauss_points(distype, degree);
+      std::shared_ptr<CollectedGaussPoints> cgp = std::make_shared<CollectedGaussPoints>(gp->num_points());
 
       GaussIntegration intpoints(gp);
 

src/core/fem/src/general/utils/4C_fem_general_utils_polynomial.cpp

@@ -10,8 +10,8 @@
 #include "4C_fem_general_utils_integration.hpp"
 #include "4C_utils_singleton_owner.hpp"
 
-#include <Intrepid_FieldContainer.hpp>
-#include <Intrepid_PointTools.hpp>
+#include <Intrepid2_PointTools.hpp>
+#include <Kokkos_DynRankView.hpp>
 #include <Shards_CellTopology.hpp>
 
 FOUR_C_NAMESPACE_OPEN
@@ -591,12 +591,12 @@ namespace Core::FE
     }
     else
     {
-      Intrepid::FieldContainer<double> wb_points(size(degree), 2);
+      Kokkos::DynRankView<double, Kokkos::HostSpace> wb_points("fekete_points", size(degree), 2);
 
       //  CellTopologyData
       shards::CellTopology myTri(shards::getCellTopologyData<shards::Triangle<3>>());
-      Intrepid::PointTools::getLattice<double, Intrepid::FieldContainer<double>>(
-          wb_points, myTri, degree, 0, Intrepid::POINTTYPE_WARPBLEND);
+      Intrepid2::PointTools::getLattice(
+          wb_points, myTri, degree, 0, Intrepid2::POINTTYPE_WARPBLEND);
 
       for (unsigned int i = 0; i < size(degree); ++i)
         for (int j = 0; j < 2; ++j) fekete_points_(j, i) = wb_points(i, j);
@@ -629,12 +629,12 @@ namespace Core::FE
     }
     else
     {
-      Intrepid::FieldContainer<double> wb_points(size(degree), 3);
+      Kokkos::DynRankView<double, Kokkos::HostSpace> wb_points("fekete_points", size(degree), 3);
 
       //  CellTopologyData
       shards::CellTopology myTet(shards::getCellTopologyData<shards::Tetrahedron<4>>());
-      Intrepid::PointTools::getLattice<double, Intrepid::FieldContainer<double>>(
-          wb_points, myTet, degree, 0, Intrepid::POINTTYPE_WARPBLEND);
+      Intrepid2::PointTools::getLattice(
+          wb_points, myTet, degree, 0, Intrepid2::POINTTYPE_WARPBLEND);
 
       for (unsigned int i = 0; i < size(degree); ++i)
         for (int j = 0; j < 3; ++j) fekete_points_(j, i) = wb_points(i, j);

src/core/fem/src/general/utils/4C_fem_general_utils_shapevalues_hdg.cpp

@@ -21,7 +21,7 @@ template <Core::FE::CellType distype>
 Core::FE::ShapeValues<distype>::ShapeValues(
     const unsigned int degree, const bool completepoly, const unsigned int quadratureDegree)
     : degree_(degree),
-      quadrature_(Core::FE::GaussPointCache::instance().create(distype, quadratureDegree)),
+      quadrature_(Core::FE::create_gauss_points(distype, quadratureDegree)),
       usescompletepoly_(completepoly),
       nqpoints_(quadrature_->num_points())
 {
@@ -160,7 +160,7 @@ Core::FE::ShapeValuesFace<distype>::ShapeValuesFace(ShapeValuesFaceParams params
   polySpace_ = Core::FE::PolynomialSpaceCache<nsd_ - 1>::instance().create(polyparams);
 
   nfdofs_ = polySpace_->size();
-  quadrature_ = Core::FE::GaussPointCache::instance().create(
+  quadrature_ = Core::FE::create_gauss_points(
       Core::FE::DisTypeToFaceShapeType<distype>::shape, params.quadraturedegree_);
   nqpoints_ = quadrature_->num_points();
 
@@ -600,7 +600,7 @@ Core::FE::ShapeValuesInteriorOnFaceCache<distype>::create(ShapeValuesFaceParams
   std::shared_ptr<PolynomialSpace<nsd>> polySpace =
       Core::FE::PolynomialSpaceCache<nsd>::instance().create(polyparams);
   Core::LinAlg::SerialDenseVector(polySpace->size());
-  std::shared_ptr<Core::FE::GaussPoints> quadrature = Core::FE::GaussPointCache::instance().create(
+  std::shared_ptr<Core::FE::GaussPoints> quadrature = Core::FE::create_gauss_points(
       Core::FE::get_ele_face_shape_type(distype, params.face_), params.quadraturedegree_);
 
   std::shared_ptr<ShapeValuesInteriorOnFace> container =

src/core/fem/tests/geometry/4C_geometry_element_volume_test.cpp

@@ -9,12 +9,19 @@
 
 #include "4C_fem_geometry_element_volume.hpp"
 
+#include <Kokkos_Core.hpp>
+
 namespace
 {
   using namespace FourC;
 
   class ElementVolumeTest : public ::testing::Test
   {
+   public:
+    static void SetUpTestSuite() { Kokkos::initialize(); }
+
+    static void TearDownTestSuite() { Kokkos::finalize(); }
+
    protected:
     void SetUp() override {}
 

src/cut/4C_cut_elementhandle.cpp

@@ -139,7 +139,7 @@ void Cut::ElementHandle::append_volume_cell_gauss_points_tessellation(
     Cut::IntegrationCell* ic = *i;
 
     std::shared_ptr<Core::FE::GaussPoints> gp_ic =
-        Core::FE::GaussPointCache::instance().create(ic->shape(), ic->cubature_degree(ic->shape()));
+        Core::FE::create_gauss_points(ic->shape(), ic->cubature_degree(ic->shape()));
     const std::vector<Cut::Point*>& cpoints = ic->points();
 
     switch (ic->shape())
@@ -285,7 +285,7 @@ std::shared_ptr<Core::FE::GaussPointsComposite> Cut::ElementHandle::gauss_points
       {
         Cut::IntegrationCell* ic = *i;
 
-        std::shared_ptr<Core::FE::GaussPoints> gp_ic = Core::FE::GaussPointCache::instance().create(
+        std::shared_ptr<Core::FE::GaussPoints> gp_ic = Core::FE::create_gauss_points(
             ic->shape(), ic->cubature_degree(ic->shape()));
         const std::vector<Cut::Point*>& cpoints = ic->points();
 

src/elemag/4C_elemag_ele_calc.cpp

@@ -599,7 +599,7 @@ void Discret::Elements::ElemagEleCalc<distype>::LocalSolver::compute_error(
   const double time = params.get<double>("time");
   // for the calculation of the error, we use a higher integration rule
   std::shared_ptr<Core::FE::GaussPoints> highquad =
-      Core::FE::GaussPointCache::instance().create(distype, (ele->degree() + 2) * 2);
+      Core::FE::create_gauss_points(distype, (ele->degree() + 2) * 2);
   Core::LinAlg::Matrix<nsd_, 1> xsi;
   Core::LinAlg::SerialDenseVector values(shapes_.ndofs_);
   Core::LinAlg::Matrix<nsd_, nen_> deriv;

src/fluid_ele/4C_fluid_ele_calc_intfaces_stab.cpp

@@ -324,7 +324,7 @@ Discret::Elements::FluidInternalSurfaceStab<distype, pdistype, ndistype>::FluidI
 
   // creating a singleton instance ensures that the object will be deleted at the end
   // create intpoints with computed degree
-  intpoints_ = Core::FE::GaussPointCache::instance().create(distype, patch_degree);
+  intpoints_ = Core::FE::create_gauss_points(distype, patch_degree);
 
   numgp_ = intpoints_->num_points();
 

src/scatra_ele/4C_scatra_ele.cpp

@@ -689,7 +689,7 @@ int Discret::Elements::Transport::initialize()
     else
       deg = 3 * this->degree();
     std::shared_ptr<Core::FE::GaussPoints> quadrature(
-        Core::FE::GaussPointCache::instance().create(this->shape(), deg));
+        Core::FE::create_gauss_points(this->shape(), deg));
     int gp = quadrature->num_points();
     if (actmat->parameter() != nullptr and
         !actmat->myocard_mat())  // in case we are not in post-process mode

src/scatra_ele/4C_scatra_ele_calc_cardiac_monodomain_hdg.cpp

@@ -330,7 +330,7 @@ void Discret::Elements::ScaTraEleCalcHDGCardiacMonodomain<distype, probdim>::mat
       deg = 3 * this->shapes_->degree_;
 
     std::shared_ptr<Core::FE::GaussPoints> quadrature_(
-        Core::FE::GaussPointCache::instance().create(distype, deg));
+        Core::FE::create_gauss_points(distype, deg));
     nqpoints = quadrature_->num_points();
 
     if (nqpoints != actmat->get_number_of_gp())

src/scatra_ele/4C_scatra_ele_hdg.cpp

@@ -321,7 +321,7 @@ int Discret::Elements::ScaTraHDG::initialize()
     else
     {
       std::shared_ptr<Core::FE::GaussPoints> quadrature_(
-          Core::FE::GaussPointCache::instance().create(this->shape(), deg));
+          Core::FE::create_gauss_points(this->shape(), deg));
       gp = quadrature_->num_points();
     }
     if (actmat->parameter() != nullptr and

tests/cut_test/cut_test_main.cpp

@@ -14,6 +14,7 @@
 #include "4C_utils_singleton_owner.hpp"
 
 #include <fenv.h>
+#include <Kokkos_Core.hpp>
 #include <Teuchos_CommandLineProcessor.hpp>
 
 #include <map>
@@ -440,6 +441,7 @@ int main(int argc, char** argv)
 {
   Core::Utils::SingletonOwnerRegistry::ScopeGuard guard;
   MPI_Init(&argc, &argv);
+  Kokkos::ScopeGuard kokkos_guard{};
 
   std::map<std::string, testfunct> functable;