Merge branch 'ShapeMatrixCacheForLF' into 'master'

Use shape matrix cache in LiquidFlow process to reduce memory consumption

See merge request ogs/ogs!4962

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:

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h

@@ -24,6 +24,7 @@
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/Integration/GenericIntegrationMethod.h"
+#include "NumLib/Fem/ShapeMatrixCache.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
@@ -36,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;
 
@@ -87,10 +86,12 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface
         std::size_t const /*local_matrix_size*/,
         NumLib::GenericIntegrationMethod const& integration_method,
         bool const is_axially_symmetric,
-        LiquidFlowData const& process_data)
+        LiquidFlowData const& process_data,
+        NumLib::ShapeMatrixCache const& shape_matrix_cache)
         : _element(element),
           _integration_method(integration_method),
-          _process_data(process_data)
+          _process_data(process_data),
+          _shape_matrix_cache(shape_matrix_cache)
     {
         unsigned const n_integration_points =
             _integration_method.getNumberOfPoints();
@@ -112,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);
@@ -135,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(
@@ -227,6 +231,7 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface
                               VelocityCacheType& darcy_velocity_at_ips) const;
 
     const LiquidFlowData& _process_data;
+    NumLib::ShapeMatrixCache const& _shape_matrix_cache;
 };
 
 }  // namespace LiquidFlow

ProcessLib/LiquidFlow/LiquidFlowProcess.cpp

@@ -41,6 +41,7 @@ LiquidFlowProcess::LiquidFlowProcess(
               integration_order, std::move(process_variables),
               std::move(secondary_variables)),
       _process_data(std::move(process_data)),
+      _shape_matrix_cache{integration_order},
       _surfaceflux(std::move(surfaceflux)),
       _is_linear(is_linear)
 {
@@ -59,7 +60,7 @@ void LiquidFlowProcess::initializeConcreteProcess(
     ProcessLib::createLocalAssemblers<LiquidFlowLocalAssembler>(
         mesh_space_dimension, mesh.getElements(), dof_table, _local_assemblers,
         NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
-        _process_data);
+        _process_data, _shape_matrix_cache);
 
     _secondary_variables.addSecondaryVariable(
         "darcy_velocity",

ProcessLib/LiquidFlow/LiquidFlowProcess.h

@@ -14,9 +14,10 @@
 
 #include <memory>
 
-#include "LiquidFlowLocalAssembler.h"
 #include "LiquidFlowData.h"
+#include "LiquidFlowLocalAssembler.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/ShapeMatrixCache.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/SurfaceFlux/SurfaceFluxData.h"
 
@@ -108,6 +109,8 @@ class LiquidFlowProcess final : public Process
     std::vector<std::unique_ptr<LiquidFlowLocalAssemblerInterface>>
         _local_assemblers;
 
+    NumLib::ShapeMatrixCache _shape_matrix_cache;
+
     std::unique_ptr<ProcessLib::SurfaceFluxData> _surfaceflux;
     MeshLib::PropertyVector<double>* _hydraulic_flow = nullptr;
     bool _is_linear = false;