From 3ba631274a98d2d936c869e15842f690f2719c1b Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Tue, 7 Feb 2023 19:26:40 +0100
Subject: [PATCH] [PL] Use eigenBlockMatrixView for N_u_op

This saves significant (up to 1/2) memory of the integration point data.
---
 .../HydroMechanics/HydroMechanicsFEM-impl.h   | 19 ++------------
 ProcessLib/HydroMechanics/HydroMechanicsFEM.h |  9 ++++---
 .../RichardsMechanics/IntegrationPointData.h  |  3 ---
 .../RichardsMechanicsFEM-impl.h               | 22 +++-------------
 .../RichardsMechanics/RichardsMechanicsFEM.h  |  5 ++++
 .../SmallDeformation/SmallDeformationFEM.h    | 20 ++++-----------
 ProcessLib/TH2M/IntegrationPointData.h        |  3 ---
 ProcessLib/TH2M/TH2MFEM-impl.h                | 25 ++++++-------------
 ProcessLib/TH2M/TH2MFEM.h                     |  5 ++++
 .../IntegrationPointData.h                    |  4 ---
 .../ThermoRichardsMechanicsFEM-impl.h         | 16 ++----------
 .../ThermoRichardsMechanicsFEM.h              |  5 ++++
 12 files changed, 41 insertions(+), 95 deletions(-)

diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 750f733d9bb..35f1400e6da 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -85,17 +85,6 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         ip_data.eps_prev.resize(kelvin_vector_size);
         ip_data.sigma_eff_prev.resize(kelvin_vector_size);
 
-        ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
-            DisplacementDim, displacement_size>::Zero(DisplacementDim,
-                                                      displacement_size);
-        for (int i = 0; i < DisplacementDim; ++i)
-        {
-            ip_data.N_u_op
-                .template block<1, displacement_size / DisplacementDim>(
-                    i, i * displacement_size / DisplacementDim)
-                .noalias() = sm_u.N;
-        }
-
         ip_data.N_u = sm_u.N;
         ip_data.dNdx_u = sm_u.dNdx;
 
@@ -208,8 +197,6 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         x_position.setIntegrationPoint(ip);
         auto const& w = _ip_data[ip].integration_weight;
 
-        auto const& N_u_op = _ip_data[ip].N_u_op;
-
         auto const& N_u = _ip_data[ip].N_u;
         auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
@@ -314,7 +301,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
         local_rhs.template segment<displacement_size>(displacement_index)
             .noalias() -=
-            (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
+            (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
 
         //
         // displacement equation, pressure part
@@ -711,8 +698,6 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         x_position.setIntegrationPoint(ip);
         auto const& w = _ip_data[ip].integration_weight;
 
-        auto const& N_u_op = _ip_data[ip].N_u_op;
-
         auto const& N_u = _ip_data[ip].N_u;
         auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
@@ -778,7 +763,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
         local_rhs.noalias() -=
             (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
-             N_u_op.transpose() * rho * b) *
+             N_u_op(N_u).transpose() * rho * b) *
             w;
     }
 }
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index f3139b7e0db..af485aa1033 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -17,6 +17,7 @@
 #include "LocalAssemblerInterface.h"
 #include "MaterialLib/PhysicalConstant.h"
 #include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/EigenBlockMatrixView.h"
 #include "MathLib/KelvinVector.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
@@ -48,10 +49,6 @@ struct IntegrationPointData final
     {
     }
 
-    typename ShapeMatricesTypeDisplacement::template MatrixType<
-        DisplacementDim, NPoints * DisplacementDim>
-        N_u_op; /**< for interpolation of the displacement vector, whereas N_u
-                   interpolates one component (scalar) */
     typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
     typename BMatricesType::KelvinVectorType eps, eps_prev;
 
@@ -180,6 +177,10 @@ public:
 
     using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
 
+    static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
+        DisplacementDim,
+        typename ShapeMatricesTypeDisplacement::NodalRowVectorType>;
+
     HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler const&) = delete;
     HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler&&) = delete;
 
diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h
index 64971c5bfc5..86c64530ab6 100644
--- a/ProcessLib/RichardsMechanics/IntegrationPointData.h
+++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h
@@ -43,9 +43,6 @@ struct IntegrationPointData final
         sigma_eff_prev.resize(kelvin_vector_size);
     }
 
-    typename ShapeMatrixTypeDisplacement::template MatrixType<
-        DisplacementDim, NPoints * DisplacementDim>
-        N_u_op;
     typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
     typename BMatricesType::KelvinVectorType sigma_sw, sigma_sw_prev;
     typename BMatricesType::KelvinVectorType eps, eps_prev;
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index e69b8227e66..42500cc7132 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -19,6 +19,7 @@
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Utils/FormEigenTensor.h"
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "MathLib/EigenBlockMatrixView.h"
 #include "MathLib/KelvinVector.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
@@ -162,17 +163,6 @@ RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
             _integration_method.getWeightedPoint(ip).getWeight() *
             sm_u.integralMeasure * sm_u.detJ;
 
-        ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
-            DisplacementDim, displacement_size>::Zero(DisplacementDim,
-                                                      displacement_size);
-        for (int i = 0; i < DisplacementDim; ++i)
-        {
-            ip_data.N_u_op
-                .template block<1, displacement_size / DisplacementDim>(
-                    i, i * displacement_size / DisplacementDim)
-                .noalias() = sm_u.N;
-        }
-
         ip_data.N_u = sm_u.N;
         ip_data.dNdx_u = sm_u.dNdx;
 
@@ -432,8 +422,6 @@ void RichardsMechanicsLocalAssembler<
         x_position.setIntegrationPoint(ip);
         auto const& w = _ip_data[ip].integration_weight;
 
-        auto const& N_u_op = _ip_data[ip].N_u_op;
-
         auto const& N_u = _ip_data[ip].N_u;
         auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
@@ -645,7 +633,7 @@ void RichardsMechanicsLocalAssembler<
         //
         double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
         rhs.template segment<displacement_size>(displacement_index).noalias() -=
-            (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
+            (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
 
         //
         // pressure equation, pressure part.
@@ -789,8 +777,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         x_position.setIntegrationPoint(ip);
         auto const& w = _ip_data[ip].integration_weight;
 
-        auto const& N_u_op = _ip_data[ip].N_u_op;
-
         auto const& N_u = _ip_data[ip].N_u;
         auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
@@ -984,7 +970,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
         local_rhs.template segment<displacement_size>(displacement_index)
             .noalias() -=
-            (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
+            (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
 
         //
         // displacement equation, pressure part
@@ -1007,7 +993,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
             .template block<displacement_size, pressure_size>(
                 displacement_index, pressure_index)
             .noalias() +=
-            N_u_op.transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
+            N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
 
         // For the swelling stress with double structure model the corresponding
         // Jacobian u-p entry would be required, but it does not improve
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index c458948496d..2a28693ed49 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -17,6 +17,7 @@
 #include "LocalAssemblerInterface.h"
 #include "MaterialLib/MPL/VariableType.h"
 #include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/EigenBlockMatrixView.h"
 #include "MathLib/KelvinVector.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
@@ -72,6 +73,10 @@ public:
 
     using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
 
+    static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
+        DisplacementDim,
+        typename ShapeMatricesTypeDisplacement::NodalRowVectorType>;
+
     RichardsMechanicsLocalAssembler(RichardsMechanicsLocalAssembler const&) =
         delete;
     RichardsMechanicsLocalAssembler(RichardsMechanicsLocalAssembler&&) = delete;
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index b7ec52b90f6..779ed47fcf5 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -18,6 +18,7 @@
 #include "LocalAssemblerInterface.h"
 #include "MaterialLib/PhysicalConstant.h"
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "MathLib/EigenBlockMatrixView.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
@@ -107,6 +108,9 @@ public:
     using IpData =
         IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>;
 
+    static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
+        DisplacementDim, typename ShapeMatricesType::NodalRowVectorType>;
+
     SmallDeformationLocalAssembler(SmallDeformationLocalAssembler const&) =
         delete;
     SmallDeformationLocalAssembler(SmallDeformationLocalAssembler&&) = delete;
@@ -272,20 +276,6 @@ public:
             auto const& N = _ip_data[ip].N;
             auto const& dNdx = _ip_data[ip].dNdx;
 
-            typename ShapeMatricesType::template MatrixType<DisplacementDim,
-                                                            displacement_size>
-                N_u_op = ShapeMatricesType::template MatrixType<
-                    DisplacementDim,
-                    displacement_size>::Zero(DisplacementDim,
-                                             displacement_size);
-            for (int i = 0; i < DisplacementDim; ++i)
-            {
-                N_u_op
-                    .template block<1, displacement_size / DisplacementDim>(
-                        i, i * displacement_size / DisplacementDim)
-                    .noalias() = N;
-            }
-
             auto const x_coord =
                 NumLib::interpolateXCoordinate<ShapeFunction,
                                                ShapeMatricesType>(_element, N);
@@ -336,7 +326,7 @@ public:
 
             auto const rho = _process_data.solid_density(t, x_position)[0];
             local_b.noalias() -=
-                (B.transpose() * sigma - N_u_op.transpose() * rho * b) * w;
+                (B.transpose() * sigma - N_u_op(N).transpose() * rho * b) * w;
             local_Jac.noalias() += B.transpose() * C * B * w;
         }
     }
diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index fed694cbf1a..1ceeda5f4ae 100644
--- a/ProcessLib/TH2M/IntegrationPointData.h
+++ b/ProcessLib/TH2M/IntegrationPointData.h
@@ -50,9 +50,6 @@ struct IntegrationPointData final
         sigma_eff_prev.resize(kelvin_vector_size);
     }
 
-    typename ShapeMatrixTypeDisplacement::template MatrixType<
-        DisplacementDim, NPoints * DisplacementDim>
-        N_u_op;
     typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
     typename BMatricesType::KelvinVectorType sigma_sw, sigma_sw_prev;
     typename BMatricesType::KelvinVectorType eps, eps_prev;
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 58849432443..ec0601f539b 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -17,6 +17,7 @@
 #include "MaterialLib/MPL/Utils/FormEigenTensor.h"
 #include "MaterialLib/PhysicalConstant.h"
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "MathLib/EigenBlockMatrixView.h"
 #include "MathLib/KelvinVector.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
@@ -74,17 +75,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             _integration_method.getWeightedPoint(ip).getWeight() *
             sm_u.integralMeasure * sm_u.detJ;
 
-        ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
-            DisplacementDim, displacement_size>::Zero(DisplacementDim,
-                                                      displacement_size);
-        for (int i = 0; i < DisplacementDim; ++i)
-        {
-            ip_data.N_u_op
-                .template block<1, displacement_size / DisplacementDim>(
-                    i, i * displacement_size / DisplacementDim)
-                .noalias() = sm_u.N;
-        }
-
         ip_data.N_u = sm_u.N;
         ip_data.dNdx_u = sm_u.dNdx;
 
@@ -1089,7 +1079,6 @@ void TH2MLocalAssembler<
         auto const& gradNpT = gradNp.transpose().eval();
         auto const& gradNTT = gradNT.transpose().eval();
 
-        auto const& Nu_op = ip.N_u_op;
         auto const& w = ip.integration_weight;
 
         auto const& m = Invariants::identity2;
@@ -1349,7 +1338,8 @@ void TH2MLocalAssembler<
 
         KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_s_L * m * Np) * w;
 
-        fU.noalias() -= (BuT * ip.sigma_eff - Nu_op.transpose() * rho * b) * w;
+        fU.noalias() -=
+            (BuT * ip.sigma_eff - N_u_op(Nu).transpose() * rho * b) * w;
 
         if (_process_data.apply_mass_lumping)
         {
@@ -1515,7 +1505,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         auto const& gradNpT = gradNp.transpose().eval();
         auto const& gradNTT = gradNT.transpose().eval();
 
-        auto const& Nu_op = ip.N_u_op;
         auto const& w = ip.integration_weight;
 
         auto const& m = Invariants::identity2;
@@ -2047,7 +2036,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             .noalias() += BuT * ip_cv.C * Bu * w;
 
         // fU_1
-        fU.noalias() -= (BuT * ip.sigma_eff - Nu_op.transpose() * rho * b) * w;
+        fU.noalias() -=
+            (BuT * ip.sigma_eff - N_u_op(Nu).transpose() * rho * b) * w;
 
         // KuT
         local_Jac
@@ -2059,8 +2049,9 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         local_Jac
             .template block<displacement_size, temperature_size>(
                 displacement_index, temperature_index)
-            .noalias() += Nu_op.transpose() * ip_cv.drho_dT * b * Nu_op * w;
-            */
+            .noalias() += N_u_op(Nu).transpose() * ip_cv.drho_dT * b *
+                          N_u_op(Nu).transpose() * w;
+         */
 
         if (_process_data.apply_mass_lumping)
         {
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index 09fc3854c29..3eb21370d8d 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -18,6 +18,7 @@
 #include "LocalAssemblerInterface.h"
 #include "MaterialLib/PhysicalConstant.h"
 #include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/EigenBlockMatrixView.h"
 #include "MathLib/KelvinVector.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
@@ -73,6 +74,10 @@ public:
         MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
     using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
 
+    static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
+        DisplacementDim,
+        typename ShapeMatricesTypeDisplacement::NodalRowVectorType>;
+
     using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
 
     TH2MLocalAssembler(TH2MLocalAssembler const&) = delete;
diff --git a/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h b/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h
index e7a31d16bea..a4a910a70fa 100644
--- a/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h
+++ b/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h
@@ -26,10 +26,6 @@ template <typename ShapeMatrixTypeDisplacement,
           typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints>
 struct IntegrationPointData final
 {
-    typename ShapeMatrixTypeDisplacement::template MatrixType<
-        DisplacementDim, NPoints * DisplacementDim>
-        N_u_op;
-
     typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
     typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
 
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
index d925a4bd769..43c1b899258 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
@@ -69,17 +69,6 @@ ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
             integration_method.getWeightedPoint(ip).getWeight() *
             sm_u.integralMeasure * sm_u.detJ;
 
-        ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
-            DisplacementDim, displacement_size>::Zero(DisplacementDim,
-                                                      displacement_size);
-        for (int i = 0; i < DisplacementDim; ++i)
-        {
-            ip_data.N_u_op
-                .template block<1, displacement_size / DisplacementDim>(
-                    i, i * displacement_size / DisplacementDim)
-                .noalias() = sm_u.N;
-        }
-
         ip_data.N_u = sm_u.N;
         ip_data.dNdx_u = sm_u.dNdx;
 
@@ -316,7 +305,6 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         typename ConstitutiveTraits::OutputData& output_data) const
 {
     auto const& N_u = ip_data.N_u;
-    auto const& N_u_op = ip_data.N_u_op;
     auto const& dNdx_u = ip_data.dNdx_u;
 
     // N and dNdx are used for both p and T variables
@@ -403,7 +391,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
     block_u(out.res).noalias() =
         B.transpose() * sigma_total -
         static_cast<int>(this->process_data_.apply_body_force_for_deformation) *
-            N_u_op.transpose() * CD.grav_data.volumetric_body_force;
+            N_u_op(N_u).transpose() * CD.grav_data.volumetric_body_force;
 
     // Storage matrices
     out.storage_p_a_p.noalias() = CD.eq_p_data.storage_p_a_p_X_NTN * NTN;
@@ -447,7 +435,7 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         B.transpose() * CD.s_mech_data.J_uT_BT_K_N * N;
     block_up(out.Jac).noalias() =
         B.transpose() * CD.s_mech_data.J_up_BT_K_N * N +
-        N_u_op.transpose() * CD.grav_data.J_up_HT_V_N * N;
+        N_u_op(N_u).transpose() * CD.grav_data.J_up_HT_V_N * N;
     block_uu(out.Jac).noalias() =
         B.transpose() * CD.s_mech_data.stiffness_tensor * B;
 
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h
index 45b30aef7df..49f4768d7d2 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h
@@ -15,6 +15,7 @@
 
 #include "IntegrationPointData.h"
 #include "LocalAssemblerInterface.h"
+#include "MathLib/EigenBlockMatrixView.h"
 #include "MathLib/KelvinVector.h"
 #include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/Integration/GenericIntegrationMethod.h"
@@ -66,6 +67,10 @@ public:
 
     using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
 
+    static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
+        DisplacementDim,
+        typename ShapeMatricesTypeDisplacement::NodalRowVectorType>;
+
     ThermoRichardsMechanicsLocalAssembler(
         ThermoRichardsMechanicsLocalAssembler const&) = delete;
     ThermoRichardsMechanicsLocalAssembler(
-- 
GitLab