[PL] SD: Compute B-matrices on the fly.

This is slightly slower, but also uses less memory.

ProcessLib/SmallDeformation/SmallDeformationFEM.h

@@ -31,7 +31,8 @@ namespace ProcessLib
 {
 namespace SmallDeformation
 {
-template <typename BMatricesType, int DisplacementDim>
+template <typename BMatricesType, typename ShapeMatricesType,
+          int DisplacementDim>
 struct IntegrationPointData final
 {
     explicit IntegrationPointData(
@@ -46,7 +47,6 @@ struct IntegrationPointData final
     // The default generated move-ctor is correctly generated for other
     // compilers.
     explicit IntegrationPointData(IntegrationPointData&& other)
-        : b_matrices(std::move(other.b_matrices)),
           sigma(std::move(other.sigma)),
           sigma_prev(std::move(other.sigma_prev)),
           eps(std::move(other.eps)),
@@ -58,7 +58,6 @@ struct IntegrationPointData final
     }
 #endif  // _MSC_VER
 
-    typename BMatricesType::BMatrixType b_matrices;
     typename BMatricesType::KelvinVectorType sigma, sigma_prev;
     typename BMatricesType::KelvinVectorType eps, eps_prev;
 
@@ -68,6 +67,8 @@ struct IntegrationPointData final
         material_state_variables;
 
     double integration_weight;
+    typename ShapeMatricesType::NodalRowVectorType N;
+    typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
 
     void pushBackState()
     {
@@ -153,12 +154,13 @@ public:
     SmallDeformationLocalAssembler(
         MeshLib::Element const& e,
         std::size_t const /*local_matrix_size*/,
-        bool is_axially_symmetric,
+        bool const is_axially_symmetric,
         unsigned const integration_order,
         SmallDeformationProcessData<DisplacementDim>& process_data)
         : _process_data(process_data),
           _integration_method(integration_order),
-          _element(e)
+          _element(e),
+          _is_axially_symmetric(is_axially_symmetric)
     {
         unsigned const n_integration_points =
             _integration_method.getNumberOfPoints();
@@ -179,25 +181,15 @@ public:
             _ip_data[ip].integration_weight =
                 _integration_method.getWeightedPoint(ip).getWeight() *
                 sm.integralMeasure * sm.detJ;
-            ip_data.b_matrices.resize(
-                KelvinVectorDimensions<DisplacementDim>::value,
-                ShapeFunction::NPOINTS * DisplacementDim);
 
-            auto const x_coord =
-                interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(e,
-                                                                         sm.N);
-            LinearBMatrix::computeBMatrix<DisplacementDim,
-                                          ShapeFunction::NPOINTS>(
-                sm.dNdx, ip_data.b_matrices, is_axially_symmetric, sm.N,
-                x_coord);
+            ip_data.N = sm.N;
+            ip_data.dNdx = sm.dNdx;
 
-            ip_data.sigma.resize(
-                KelvinVectorDimensions<DisplacementDim>::value);
+            ip_data.sigma.resize(KelvinVectorDimensions<DisplacementDim>::value);
             ip_data.sigma_prev.resize(
                 KelvinVectorDimensions<DisplacementDim>::value);
             ip_data.eps.resize(KelvinVectorDimensions<DisplacementDim>::value);
-            ip_data.eps_prev.resize(
-                KelvinVectorDimensions<DisplacementDim>::value);
+            ip_data.eps_prev.resize(KelvinVectorDimensions<DisplacementDim>::value);
 
             _secondary_data.N[ip] = shape_matrices[ip].N;
         }
@@ -240,8 +232,19 @@ public:
         {
             x_position.setIntegrationPoint(ip);
             auto const& w = _ip_data[ip].integration_weight;
+            auto const& N = _ip_data[ip].N;
+            auto const& dNdx = _ip_data[ip].dNdx;
+
+            typename BMatricesType::BMatrixType B(
+                KelvinVectorDimensions<DisplacementDim>(),
+                ShapeFunction::NPOINTS * DisplacementDim);
+            auto const x_coord =
+                interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+                    _element, N);
+            LinearBMatrix::computeBMatrix<DisplacementDim,
+                                          ShapeFunction::NPOINTS>(
+                dNdx, B, _is_axially_symmetric, N, x_coord);
 
-            auto const& B = _ip_data[ip].b_matrices;
             auto const& eps_prev = _ip_data[ip].eps_prev;
             auto const& sigma_prev = _ip_data[ip].sigma_prev;
 
@@ -286,9 +289,10 @@ public:
         std::vector<double>& nodal_values) const override
     {
         return ProcessLib::SmallDeformation::getNodalForces<
-            DisplacementDim, ShapeFunction::NPOINTS,
-            NodalDisplacementVectorType>(nodal_values, _integration_method,
-                                         _ip_data, _element.getID());
+            DisplacementDim, ShapeFunction, ShapeMatricesType,
+            NodalDisplacementVectorType, typename BMatricesType::BMatrixType>(
+            nodal_values, _integration_method, _ip_data, _element,
+            _is_axially_symmetric);
     }
 
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
@@ -413,14 +417,16 @@ private:
 
     SmallDeformationProcessData<DisplacementDim>& _process_data;
 
-    std::vector<IntegrationPointData<BMatricesType, DisplacementDim>,
-                Eigen::aligned_allocator<
-                    IntegrationPointData<BMatricesType, DisplacementDim>>>
+    std::vector<
+        IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>,
+        Eigen::aligned_allocator<IntegrationPointData<
+            BMatricesType, ShapeMatricesType, DisplacementDim>>>
         _ip_data;
 
     IntegrationMethod _integration_method;
     MeshLib::Element const& _element;
     SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+    bool const _is_axially_symmetric;
 };
 
 template <typename ShapeFunction, typename IntegrationMethod,

ProcessLib/SmallDeformationCommon/Common.h

@@ -16,6 +16,9 @@
 
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
 
 namespace ProcessLib
 {
@@ -29,30 +32,37 @@ struct NodalForceCalculationInterface
     virtual ~NodalForceCalculationInterface() = default;
 };
 
-template <int DisplacementDim, int NPoints,
-          typename NodalDisplacementVectorType, typename IPData,
-          typename IntegrationMethod>
+template <int DisplacementDim, typename ShapeFunction,
+          typename ShapeMatricesType, typename NodalDisplacementVectorType,
+          typename BMatrixType, typename IPData, typename IntegrationMethod>
 std::vector<double> const& getNodalForces(
     std::vector<double>& nodal_values,
     IntegrationMethod const& _integration_method, IPData const& _ip_data,
-    int element_id)
+    MeshLib::Element const& element, bool const is_axially_symmetric)
 {
     nodal_values.clear();
     auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
-        nodal_values, NPoints * DisplacementDim);
+        nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
 
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
 
     SpatialPosition x_position;
-    x_position.setElementID(element_id);
+    x_position.setElementID(element.getID());
 
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
         x_position.setIntegrationPoint(ip);
         auto const& w = _ip_data[ip].integration_weight;
 
-        auto const& B = _ip_data[ip].b_matrices;
+        BMatrixType B(KelvinVectorDimensions<DisplacementDim>(),
+                      ShapeFunction::NPOINTS * DisplacementDim);
+        auto const x_coord =
+            interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+                element, _ip_data[ip].N);
+        LinearBMatrix::computeBMatrix<DisplacementDim, ShapeFunction::NPOINTS>(
+            _ip_data[ip].dNdx, B, is_axially_symmetric, _ip_data[ip].N,
+            x_coord);
         auto& sigma = _ip_data[ip].sigma;
 
         local_b.noalias() += B.transpose() * sigma * w;