diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 3ed82f8411bce855d2be453f60bdf025a8f08156..9e2804d38682166957787dcd2ba73c2629fa96e3 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -91,7 +91,8 @@ public:
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {
-        auto const& N = _ip_data[integration_point].N;
+        auto const& N = _shape_matrix_cache.NsHigherOrder<
+            typename ShapeFunction::MeshElement>()[integration_point];
 
         // assumes N is stored contiguously in memory
         return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
@@ -273,11 +274,15 @@ protected:
             *_process_data.media_map.getMedium(_element.getID());
         auto const& liquid_phase = medium.phase("AqueousLiquid");
 
+        auto const& Ns =
+            _shape_matrix_cache
+                .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
         for (unsigned ip = 0; ip < n_integration_points; ++ip)
         {
             auto const& ip_data = _ip_data[ip];
-            auto const& N = ip_data.N;
             auto const& dNdx = ip_data.dNdx;
+            auto const& N = Ns[ip];
 
             pos.setIntegrationPoint(ip);
 
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index 73474c5908ea2c114ca99986a564bdc43556f79d..ff7381ad074c6a110c48fcea1bb3ecfac2aa0210 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -123,13 +123,17 @@ public:
         double average_velocity_norm = 0.0;
         ip_flux_vector.reserve(n_integration_points);
 
+        auto const& Ns =
+            this->_shape_matrix_cache
+                .template NsHigherOrder<typename ShapeFunction::MeshElement>();
+
         for (unsigned ip(0); ip < n_integration_points; ip++)
         {
             pos.setIntegrationPoint(ip);
 
             auto const& ip_data = this->_ip_data[ip];
-            auto const& N = ip_data.N;
             auto const& dNdx = ip_data.dNdx;
+            auto const& N = Ns[ip];
             auto const& w = ip_data.integration_weight;
 
             double T_int_pt = 0.0;
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index 65c67a68b1ebc0dbfcf86acbc8aa760b743bfdc0..7a828cdddacd21a18fefe43b111026bf5304e10c 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -79,13 +79,17 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHydraulicEquation(
     unsigned const n_integration_points =
         this->_integration_method.getNumberOfPoints();
 
+    auto const& Ns =
+        this->_shape_matrix_cache
+            .template NsHigherOrder<typename ShapeFunction::MeshElement>();
+
     for (unsigned ip(0); ip < n_integration_points; ip++)
     {
         pos.setIntegrationPoint(ip);
 
         auto const& ip_data = this->_ip_data[ip];
-        auto const& N = ip_data.N;
         auto const& dNdx = ip_data.dNdx;
+        auto const& N = Ns[ip];
         auto const& w = ip_data.integration_weight;
 
         double p_int_pt = 0.0;
@@ -208,13 +212,17 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHeatTransportEquation(
     double average_velocity_norm = 0.0;
     ip_flux_vector.reserve(n_integration_points);
 
+    auto const& Ns =
+        this->_shape_matrix_cache
+            .template NsHigherOrder<typename ShapeFunction::MeshElement>();
+
     for (unsigned ip(0); ip < n_integration_points; ip++)
     {
         pos.setIntegrationPoint(ip);
 
         auto const& ip_data = this->_ip_data[ip];
-        auto const& N = ip_data.N;
         auto const& dNdx = ip_data.dNdx;
+        auto const& N = Ns[ip];
         auto const& w = ip_data.integration_weight;
 
         double p_at_xi = 0.;