diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 16d2116f6f4..4b3e41fffcc 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -143,12 +143,15 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
         _process_data.element_rotation_matrices[_element.getID()].transpose() *
         _process_data.specific_body_force;
 
+    auto const& N = _shape_matrix_cache
+                        .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
         auto const& ip_data = _ip_data[ip];
 
         double p = 0.;
-        NumLib::shapeFunctionInterpolate(local_x, ip_data.N, p);
+        NumLib::shapeFunctionInterpolate(local_x, N[ip], p);
         vars.liquid_phase_pressure = p;
 
         // Compute density:
@@ -173,7 +176,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
         // Assemble mass matrix, M
         local_M.noalias() +=
             (porosity * ddensity_dpressure / fluid_density + storage) *
-            ip_data.N.transpose() * ip_data.N * ip_data.integration_weight;
+            N[ip].transpose() * N[ip] * ip_data.integration_weight;
 
         // Compute viscosity:
         auto const viscosity =
@@ -292,11 +295,14 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
         _process_data.element_rotation_matrices[_element.getID()].transpose() *
         _process_data.specific_body_force;
 
+    auto const& N = _shape_matrix_cache
+                        .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
         auto const& ip_data = _ip_data[ip];
         double p = 0.;
-        NumLib::shapeFunctionInterpolate(local_x, ip_data.N, p);
+        NumLib::shapeFunctionInterpolate(local_x, N[ip], p);
         vars.liquid_phase_pressure = p;
 
         // Compute density:
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 87650cd538b..86ae1838af3 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -37,14 +37,12 @@ namespace LiquidFlow
 template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
 struct IntegrationPointData final
 {
-    explicit IntegrationPointData(NodalRowVectorType const& N_,
-                                  GlobalDimNodalMatrixType const& dNdx_,
+    explicit IntegrationPointData(GlobalDimNodalMatrixType const& dNdx_,
                                   double const& integration_weight_)
-        : N(N_), dNdx(dNdx_), integration_weight(integration_weight_)
+        : dNdx(dNdx_), integration_weight(integration_weight_)
     {
     }
 
-    NodalRowVectorType const N;
     GlobalDimNodalMatrixType const dNdx;
     double const integration_weight;
 
@@ -115,7 +113,7 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             _ip_data.emplace_back(
-                shape_matrices[ip].N, shape_matrices[ip].dNdx,
+                shape_matrices[ip].dNdx,
                 _integration_method.getWeightedPoint(ip).getWeight() *
                     shape_matrices[ip].integralMeasure *
                     shape_matrices[ip].detJ * aperture_size);
@@ -138,10 +136,13 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface
     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>();
 
         // assumes N is stored contiguously in memory
-        return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+        return Eigen::Map<const Eigen::RowVectorXd>(
+            N[integration_point].data(), N[integration_point].size());
     }
 
     std::vector<double> const& getIntPtDarcyVelocity(