diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 82404610c262315987be0e364bb7207d784e9c79..99bc458dfdc0408c641f56203e39b390f2a9b4c5 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -143,15 +143,16 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
         _process_data.element_rotation_matrices[_element.getID()].transpose() *
         _process_data.specific_body_force;
 
-    auto const& N = _process_data.shape_matrix_cache
-                        .NsHigherOrder<typename ShapeFunction::MeshElement>();
+    auto const& Ns = _process_data.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 = Ns[ip];
 
         double p = 0.;
-        NumLib::shapeFunctionInterpolate(local_x, N[ip], p);
+        NumLib::shapeFunctionInterpolate(local_x, N, p);
         vars.liquid_phase_pressure = p;
 
         // Compute density:
@@ -176,7 +177,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
         // Assemble mass matrix, M
         local_M.noalias() +=
             (porosity * ddensity_dpressure / fluid_density + storage) *
-            N[ip].transpose() * N[ip] * ip_data.integration_weight;
+            N.transpose() * N * ip_data.integration_weight;
 
         // Compute viscosity:
         auto const viscosity =
@@ -295,14 +296,15 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
         _process_data.element_rotation_matrices[_element.getID()].transpose() *
         _process_data.specific_body_force;
 
-    auto const& N = _process_data.shape_matrix_cache
-                        .NsHigherOrder<typename ShapeFunction::MeshElement>();
+    auto const& Ns = _process_data.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 = Ns[ip];
         double p = 0.;
-        NumLib::shapeFunctionInterpolate(local_x, N[ip], p);
+        NumLib::shapeFunctionInterpolate(local_x, N, p);
         vars.liquid_phase_pressure = p;
 
         // Compute density: