diff --git a/src/core/fem/src/general/utils/4C_fem_general_utils_gausspoints.cpp b/src/core/fem/src/general/utils/4C_fem_general_utils_gausspoints.cpp
index 3f2112e02d5..e7c453d80d6 100644
--- a/src/core/fem/src/general/utils/4C_fem_general_utils_gausspoints.cpp
+++ b/src/core/fem/src/general/utils/4C_fem_general_utils_gausspoints.cpp
@@ -7,14 +7,9 @@
 
 #include "4C_fem_general_utils_gausspoints.hpp"
 
-#include <Intrepid_Cubature.hpp>
-#include <Intrepid_DefaultCubatureFactory.hpp>
-#include <Intrepid_FieldContainer.hpp>
+#include <Intrepid2_DefaultCubatureFactory.hpp>
 #include <Shards_BasicTopologies.hpp>
 #include <Shards_CellTopology.hpp>
-#include <Shards_CellTopologyTraits.hpp>
-
-#include <stdexcept>
 
 FOUR_C_NAMESPACE_OPEN
 
@@ -22,39 +17,82 @@ namespace Core::FE
 {
   namespace
   {
+    void fill_pyramid5(Kokkos::DynRankView<double, Kokkos::HostSpace>& cub_points,
+        Kokkos::DynRankView<double, Kokkos::HostSpace>& cub_weights)
+    {
+      Kokkos::resize(cub_points, 8, 3);
+      Kokkos::resize(cub_weights, 8);
+      cub_points(0, 0) = -0.26318405556971;
+      cub_points(1, 0) = -0.50661630334979;
+      cub_points(2, 0) = -0.26318405556971;
+      cub_points(3, 0) = -0.50661630334979;
+      cub_points(4, 0) = 0.26318405556971;
+      cub_points(5, 0) = 0.50661630334979;
+      cub_points(6, 0) = 0.26318405556971;
+      cub_points(7, 0) = 0.50661630334979;
+      cub_points(0, 1) = -0.26318405556971;
+      cub_points(1, 1) = -0.50661630334979;
+      cub_points(2, 1) = 0.26318405556971;
+      cub_points(3, 1) = 0.50661630334979;
+      cub_points(4, 1) = -0.26318405556971;
+      cub_points(5, 1) = -0.50661630334979;
+      cub_points(6, 1) = 0.26318405556971;
+      cub_points(7, 1) = 0.50661630334979;
+      cub_points(0, 2) = 0.54415184401122;
+      cub_points(1, 2) = 0.12251482265544;
+      cub_points(2, 2) = 0.54415184401122;
+      cub_points(3, 2) = 0.12251482265544;
+      cub_points(4, 2) = 0.54415184401122;
+      cub_points(5, 2) = 0.12251482265544;
+      cub_points(6, 2) = 0.54415184401122;
+      cub_points(7, 2) = 0.12251482265544;
+
+      cub_weights(0) = 0.10078588207983;
+      cub_weights(1) = 0.23254745125351;
+      cub_weights(2) = 0.10078588207983;
+      cub_weights(3) = 0.23254745125351;
+      cub_weights(4) = 0.10078588207983;
+      cub_weights(5) = 0.23254745125351;
+      cub_weights(6) = 0.10078588207983;
+      cub_weights(7) = 0.23254745125351;
+    }
+
     /// wrapper to intrepid gauss point implementation
-    template <class Topology>
     class IntrepidGaussPoints : public GaussPoints
     {
      public:
-      explicit IntrepidGaussPoints(int cubDegree)
+      explicit IntrepidGaussPoints(const shards::CellTopology& cell_topology, int cubDegree)
+          : cub_points_("cubature_points", 1, 1),
+            cub_weights_("cubature_weights", 1),
+            cell_topology_(cell_topology)
       {
-        // cell type: tetrahedron
-        shards::CellTopology cellType = shards::getCellTopologyData<Topology>();
+        // Special case that Intrepid does not support
+        if (cell_topology ==
+            shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5>>()))
+        {
+          fill_pyramid5(cub_points_, cub_weights_);
+          return;
+        }
 
         // retrieve spatial dimension
-        int spaceDim = cellType.getDimension();
-        // int numNodes = cellType.getNodeCount();
+        const int spaceDim = cell_topology_.getDimension();
 
-        // create cubature factory
-        Intrepid::DefaultCubatureFactory<double> cubFactory;
-
-        // create default cubature
-        Teuchos::RCP<Intrepid::Cubature<double>> myCub = cubFactory.create(cellType, cubDegree);
+        Teuchos::RCP<Intrepid2::Cubature<Kokkos::HostSpace, double, double>> myCub =
+            Intrepid2::DefaultCubatureFactory::create<Kokkos::HostSpace>(cell_topology_, cubDegree);
 
         // retrieve number of cubature points
         int numCubPoints = myCub->getNumPoints();
 
-        cub_points_.resize(numCubPoints, spaceDim);
-        cub_weights_.resize(numCubPoints);
+        Kokkos::resize(cub_points_, numCubPoints, spaceDim);
+        Kokkos::resize(cub_weights_, numCubPoints);
 
         // retrieve cubature points and weights
         myCub->getCubature(cub_points_, cub_weights_);
       }
 
-      int num_points() const override { return cub_points_.dimension(0); }
+      int num_points() const override { return cub_points_.extent(0); }
 
-      int num_dimension() const override { return cub_points_.dimension(1); }
+      int num_dimension() const override { return cub_points_.extent(1); }
 
       const double* point(int point) const override { return &cub_points_(point, 0); }
 
@@ -62,253 +100,132 @@ namespace Core::FE
 
       void print() const override
       {
-        // cell type: tetrahedron
-        shards::CellTopology cellType = shards::getCellTopologyData<Topology>();
-
-        std::cout << cellType.getName() << " gauss points:\n";
+        std::cout << cell_topology_.getName() << " gauss points:\n";
         for (int i = 0; i < num_points(); ++i)
         {
-          std::cout << "    ";
+          std::cout << "coordinate: ";
           for (int j = 0; j < num_dimension(); ++j) std::cout << cub_points_(i, j) << " ";
-          std::cout << cub_weights_(i) << "\n";
+          std::cout << "weight: " << cub_weights_(i) << "\n";
         }
       }
 
      private:
-      Intrepid::FieldContainer<double> cub_points_;
-      Intrepid::FieldContainer<double> cub_weights_;
-    };
+      Kokkos::DynRankView<double, Kokkos::HostSpace> cub_points_;
+      Kokkos::DynRankView<double, Kokkos::HostSpace> cub_weights_;
 
+      shards::CellTopology cell_topology_;
+    };
 
-    /// special case that is not (yet) supported by intrepid
-    template <>
-    class IntrepidGaussPoints<shards::Pyramid<5>> : public GaussPoints
+    template <typename Topology>
+    std::shared_ptr<IntrepidGaussPoints> make_intrepid_gauss_points(int cub_degree)
     {
-     public:
-      explicit IntrepidGaussPoints(int cubDegree)
-      {
-        cub_points_.resize(8, 3);
-        cub_weights_.resize(8);
-
-        cub_points_(0, 0) = -0.26318405556971;
-        cub_points_(1, 0) = -0.50661630334979;
-        cub_points_(2, 0) = -0.26318405556971;
-        cub_points_(3, 0) = -0.50661630334979;
-        cub_points_(4, 0) = 0.26318405556971;
-        cub_points_(5, 0) = 0.50661630334979;
-        cub_points_(6, 0) = 0.26318405556971;
-        cub_points_(7, 0) = 0.50661630334979;
-        cub_points_(0, 1) = -0.26318405556971;
-        cub_points_(1, 1) = -0.50661630334979;
-        cub_points_(2, 1) = 0.26318405556971;
-        cub_points_(3, 1) = 0.50661630334979;
-        cub_points_(4, 1) = -0.26318405556971;
-        cub_points_(5, 1) = -0.50661630334979;
-        cub_points_(6, 1) = 0.26318405556971;
-        cub_points_(7, 1) = 0.50661630334979;
-        cub_points_(0, 2) = 0.54415184401122;
-        cub_points_(1, 2) = 0.12251482265544;
-        cub_points_(2, 2) = 0.54415184401122;
-        cub_points_(3, 2) = 0.12251482265544;
-        cub_points_(4, 2) = 0.54415184401122;
-        cub_points_(5, 2) = 0.12251482265544;
-        cub_points_(6, 2) = 0.54415184401122;
-        cub_points_(7, 2) = 0.12251482265544;
-
-        cub_weights_(0) = 0.10078588207983;
-        cub_weights_(1) = 0.23254745125351;
-        cub_weights_(2) = 0.10078588207983;
-        cub_weights_(3) = 0.23254745125351;
-        cub_weights_(4) = 0.10078588207983;
-        cub_weights_(5) = 0.23254745125351;
-        cub_weights_(6) = 0.10078588207983;
-        cub_weights_(7) = 0.23254745125351;
-      }
-
-      int num_points() const override { return cub_points_.dimension(0); }
-
-      int num_dimension() const override { return cub_points_.dimension(1); }
-
-      const double* point(int point) const override { return &cub_points_(point, 0); }
-
-      double weight(int point) const override { return cub_weights_(point); }
-
-      void print() const override
-      {
-        // cell type: tetrahedron
-        shards::CellTopology cellType = shards::getCellTopologyData<shards::Pyramid<5>>();
-
-        std::cout << cellType.getName() << " gauss points:\n";
-        for (int i = 0; i < num_points(); ++i)
-        {
-          std::cout << "    ";
-          for (int j = 0; j < num_dimension(); ++j) std::cout << cub_points_(i, j) << " ";
-          std::cout << cub_weights_(i) << "\n";
-        }
-      }
-
-     private:
-      Intrepid::FieldContainer<double> cub_points_;
-      Intrepid::FieldContainer<double> cub_weights_;
-    };
+      const shards::CellTopology cell_topology = shards::getCellTopologyData<Topology>();
+      return std::make_shared<IntrepidGaussPoints>(cell_topology, cub_degree);
+    }
 
   }  // namespace
 }  // namespace Core::FE
 
 Core::FE::GaussIntegration::GaussIntegration(Core::FE::CellType distype)
+    : gp_(create_gauss_points_default(distype))
+{
+}
+
+Core::FE::GaussIntegration::GaussIntegration(Core::FE::CellType distype, int degree)
+    : gp_(create_gauss_points(distype, degree))
+{
+}
+
+std::shared_ptr<Core::FE::GaussPoints> Core::FE::create_gauss_points(
+    Core::FE::CellType distype, int degree)
 {
   switch (distype)
   {
     case Core::FE::CellType::quad4:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::quad4, 3);
-      break;
+      return make_intrepid_gauss_points<shards::Quadrilateral<4>>(degree);
     case Core::FE::CellType::quad8:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::quad8, 4);
-      break;
+      return make_intrepid_gauss_points<shards::Quadrilateral<8>>(degree);
     case Core::FE::CellType::quad9:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::quad9, 4);
-      break;
+      return make_intrepid_gauss_points<shards::Quadrilateral<9>>(degree);
     case Core::FE::CellType::tri3:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::tri3, 3);
-      break;
+      return make_intrepid_gauss_points<shards::Triangle<3>>(degree);
     case Core::FE::CellType::tri6:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::tri6, 4);
-      break;
+      return make_intrepid_gauss_points<shards::Triangle<6>>(degree);
     case Core::FE::CellType::hex8:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::hex8, 3);
-      break;
+      return make_intrepid_gauss_points<shards::Hexahedron<8>>(degree);
     case Core::FE::CellType::hex20:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::hex20, 4);
-      break;
+      return make_intrepid_gauss_points<shards::Hexahedron<20>>(degree);
     case Core::FE::CellType::hex27:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::hex27, 4);
-      break;
+      return make_intrepid_gauss_points<shards::Hexahedron<27>>(degree);
     case Core::FE::CellType::tet4:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::tet4, 3);
-      break;
+      return make_intrepid_gauss_points<shards::Tetrahedron<4>>(degree);
     case Core::FE::CellType::tet10:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::tet10, 4);
-      break;
+      return make_intrepid_gauss_points<shards::Tetrahedron<10>>(degree);
     case Core::FE::CellType::wedge6:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::wedge6, 3);
-      break;
+      return make_intrepid_gauss_points<shards::Wedge<6>>(degree);
     case Core::FE::CellType::wedge15:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::wedge15, 4);
-      break;
+      return make_intrepid_gauss_points<shards::Wedge<15>>(degree);
     case Core::FE::CellType::pyramid5:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::pyramid5, 3);
-      break;
+      return make_intrepid_gauss_points<shards::Pyramid<5>>(degree);
     case Core::FE::CellType::line2:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::line2, 3);
-      break;
+      return make_intrepid_gauss_points<shards::Line<2>>(degree);
     case Core::FE::CellType::line3:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::line3, 4);
-      break;
-    case Core::FE::CellType::nurbs2:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::line2, 3);
-      break;
-    case Core::FE::CellType::nurbs3:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::line3, 4);
-      break;
-    case Core::FE::CellType::nurbs4:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::quad4, 3);
-      break;
-    case Core::FE::CellType::nurbs8:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::hex8, 3);
-      break;
-    case Core::FE::CellType::nurbs9:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::quad9, 4);
-      break;
-    case Core::FE::CellType::nurbs27:
-      gp_ = GaussPointCache::instance().create(Core::FE::CellType::hex27, 4);
-      break;
+      return make_intrepid_gauss_points<shards::Line<3>>(degree);
     default:
       FOUR_C_THROW("unsupported element shape");
   }
 }
 
-Core::FE::GaussIntegration::GaussIntegration(Core::FE::CellType distype, int degree)
-{
-  gp_ = GaussPointCache::instance().create(distype, degree);
-}
 
-Core::FE::GaussPointCache& Core::FE::GaussPointCache::instance()
+std::shared_ptr<Core::FE::GaussPoints> Core::FE::create_gauss_points_default(Core::FE::CellType distype)
 {
-  static std::unique_ptr<GaussPointCache> instance;
-  if (instance == nullptr)
-  {
-    instance = std::unique_ptr<GaussPointCache>(new GaussPointCache);
-  }
-  return *instance;
-}
-
-
-std::shared_ptr<Core::FE::GaussPoints> Core::FE::GaussPointCache::create(
-    Core::FE::CellType distype, int degree)
-{
-  std::map<std::pair<Core::FE::CellType, int>, std::shared_ptr<GaussPoints>>::iterator i =
-      gp_cache_.find(std::make_pair(distype, degree));
-  if (i != gp_cache_.end())
-  {
-    return i->second;
-  }
-
-  // this is expensive and should not be done too often
-  std::shared_ptr<GaussPoints> gp;
-
   switch (distype)
   {
     case Core::FE::CellType::quad4:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Quadrilateral<4>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::quad4, 3);
     case Core::FE::CellType::quad8:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Quadrilateral<8>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::quad8, 4);
     case Core::FE::CellType::quad9:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Quadrilateral<9>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::quad9, 4);
     case Core::FE::CellType::tri3:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Triangle<3>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::tri3, 3);
     case Core::FE::CellType::tri6:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Triangle<6>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::tri6, 4);
     case Core::FE::CellType::hex8:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Hexahedron<8>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::hex8, 3);
     case Core::FE::CellType::hex20:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Hexahedron<20>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::hex20, 4);
     case Core::FE::CellType::hex27:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Hexahedron<27>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::hex27, 4);
     case Core::FE::CellType::tet4:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Tetrahedron<4>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::tet4, 3);
     case Core::FE::CellType::tet10:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Tetrahedron<10>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::tet10, 4);
     case Core::FE::CellType::wedge6:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Wedge<6>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::wedge6, 3);
     case Core::FE::CellType::wedge15:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Wedge<15>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::wedge15, 4);
     case Core::FE::CellType::pyramid5:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Pyramid<5>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::pyramid5, 3);
     case Core::FE::CellType::line2:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Line<2>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::line2, 3);
     case Core::FE::CellType::line3:
-      gp = std::make_shared<IntrepidGaussPoints<shards::Line<3>>>(degree);
-      break;
+      return create_gauss_points(Core::FE::CellType::line3, 4);
+    case Core::FE::CellType::nurbs2:
+      return create_gauss_points(Core::FE::CellType::line2, 3);
+    case Core::FE::CellType::nurbs3:
+      return create_gauss_points(Core::FE::CellType::line3, 4);
+    case Core::FE::CellType::nurbs4:
+      return create_gauss_points(Core::FE::CellType::quad4, 3);
+    case Core::FE::CellType::nurbs8:
+      return create_gauss_points(Core::FE::CellType::hex8, 3);
+    case Core::FE::CellType::nurbs9:
+      return create_gauss_points(Core::FE::CellType::quad9, 4);
+    case Core::FE::CellType::nurbs27:
+      return create_gauss_points(Core::FE::CellType::hex27, 4);
     default:
       FOUR_C_THROW("unsupported element shape");
   }
-
-  gp_cache_[std::make_pair(distype, degree)] = gp;
-  return gp;
 }
 
 FOUR_C_NAMESPACE_CLOSE
diff --git a/src/core/fem/src/general/utils/4C_fem_general_utils_gausspoints.hpp b/src/core/fem/src/general/utils/4C_fem_general_utils_gausspoints.hpp
index 6214452fd8a..3cf458f53b6 100644
--- a/src/core/fem/src/general/utils/4C_fem_general_utils_gausspoints.hpp
+++ b/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);
 
diff --git a/src/core/fem/src/general/utils/4C_fem_general_utils_polynomial.cpp b/src/core/fem/src/general/utils/4C_fem_general_utils_polynomial.cpp
index 8633f6b7208..38fb7aba003 100644
--- a/src/core/fem/src/general/utils/4C_fem_general_utils_polynomial.cpp
+++ b/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);
diff --git a/src/core/fem/src/general/utils/4C_fem_general_utils_shapevalues_hdg.cpp b/src/core/fem/src/general/utils/4C_fem_general_utils_shapevalues_hdg.cpp
index 484effbd3c5..301d4852982 100644
--- a/src/core/fem/src/general/utils/4C_fem_general_utils_shapevalues_hdg.cpp
+++ b/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 =
diff --git a/src/core/fem/tests/geometry/4C_geometry_element_volume_test.cpp b/src/core/fem/tests/geometry/4C_geometry_element_volume_test.cpp
index 6d88de742ce..16becf589ed 100644
--- a/src/core/fem/tests/geometry/4C_geometry_element_volume_test.cpp
+++ b/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 {}
 
diff --git a/src/cut/4C_cut_elementhandle.cpp b/src/cut/4C_cut_elementhandle.cpp
index 42f3c09d716..e924692c449 100644
--- a/src/cut/4C_cut_elementhandle.cpp
+++ b/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();
 
diff --git a/src/elemag/4C_elemag_ele_calc.cpp b/src/elemag/4C_elemag_ele_calc.cpp
index 1b576c75071..a36a1fe0b72 100644
--- a/src/elemag/4C_elemag_ele_calc.cpp
+++ b/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;
diff --git a/src/fluid_ele/4C_fluid_ele_calc_intfaces_stab.cpp b/src/fluid_ele/4C_fluid_ele_calc_intfaces_stab.cpp
index 66ed902773f..1f79e0a5134 100644
--- a/src/fluid_ele/4C_fluid_ele_calc_intfaces_stab.cpp
+++ b/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();
 
diff --git a/src/scatra_ele/4C_scatra_ele.cpp b/src/scatra_ele/4C_scatra_ele.cpp
index 2077b95ad6c..2569d835c3e 100644
--- a/src/scatra_ele/4C_scatra_ele.cpp
+++ b/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
diff --git a/src/scatra_ele/4C_scatra_ele_calc_cardiac_monodomain_hdg.cpp b/src/scatra_ele/4C_scatra_ele_calc_cardiac_monodomain_hdg.cpp
index 1f0b26f5975..5a8cc78688f 100644
--- a/src/scatra_ele/4C_scatra_ele_calc_cardiac_monodomain_hdg.cpp
+++ b/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())
diff --git a/src/scatra_ele/4C_scatra_ele_hdg.cpp b/src/scatra_ele/4C_scatra_ele_hdg.cpp
index 5cea7a666f9..aa4ddcc6db9 100644
--- a/src/scatra_ele/4C_scatra_ele_hdg.cpp
+++ b/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
diff --git a/tests/cut_test/cut_test_main.cpp b/tests/cut_test/cut_test_main.cpp
index 699aadda4a4..b412074ffc4 100644
--- a/tests/cut_test/cut_test_main.cpp
+++ b/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;