From e1d781a38dfdbb3d518dd47954c8a46cb6196334 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Tue, 7 Feb 2023 18:00:56 +0100
Subject: [PATCH] [PL/THM] Use eigen block matrix view in N_u_op expr.

Further reduces memory allocations and slightly faster too.
---
 .../ThermoHydroMechanicsFEM-impl.h                 | 14 ++------------
 .../ThermoHydroMechanics/ThermoHydroMechanicsFEM.h |  5 +++++
 2 files changed, 7 insertions(+), 12 deletions(-)

diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
index a18efd608a3..a04d34336f0 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -18,6 +18,7 @@
 #include "MaterialLib/MPL/Utils/FormKelvinVector.h"
 #include "MaterialLib/MPL/Utils/GetLiquidThermalExpansivity.h"
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "MathLib/EigenBlockMatrixView.h"
 #include "MathLib/KelvinVector.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
@@ -399,20 +400,9 @@ void ThermoHydroMechanicsLocalAssembler<
         auto const rho =
             solid_density * (1 - porosity) + porosity * fluid_density;
 
-        typename ShapeMatricesTypeDisplacement::template MatrixType<
-            DisplacementDim, displacement_size>
-            N_u_op = ShapeMatricesTypeDisplacement::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_u;
-
         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 (K_up)
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h
index 01487f32c7e..fe6b242e113 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.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"
@@ -69,6 +70,10 @@ public:
 
     using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
 
+    static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
+        DisplacementDim,
+        typename ShapeMatricesTypeDisplacement::NodalRowVectorType>;
+
     ThermoHydroMechanicsLocalAssembler(
         ThermoHydroMechanicsLocalAssembler const&) = delete;
     ThermoHydroMechanicsLocalAssembler(ThermoHydroMechanicsLocalAssembler&&) =
-- 
GitLab