Merge pull request #12388 from KratosMultiphysics/geo/12371-extract-s…

…tiffness-matrix-utility

[GeoMechanicsApplication] extract stiffness matrix utility

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp

@@ -856,7 +856,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateMaterialStiffnessMatrix(Ma
     const GeometryType&                             rGeom = this->GetGeometry();
     const GeometryType::IntegrationPointsArrayType& IntegrationPoints =
         rGeom.IntegrationPoints(mThisIntegrationMethod);
-    const IndexType NumGPoints = IntegrationPoints.size();
 
     // Constitutive Law parameters
     ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom, rProp, rCurrentProcessInfo);
@@ -873,21 +872,13 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateMaterialStiffnessMatrix(Ma
     this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters,
                                          Variables.NContainer, Variables.DN_DXContainer,
                                          strain_vectors, mStressVector, constitutive_matrices);
+    const auto integration_coefficients =
+        this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer);
 
-    for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
-        this->CalculateKinematics(Variables, GPoint);
-        Variables.B                  = b_matrices[GPoint];
-        Variables.F                  = deformation_gradients[GPoint];
-        Variables.StrainVector       = strain_vectors[GPoint];
-        Variables.ConstitutiveMatrix = constitutive_matrices[GPoint];
-
-        // Compute weighting coefficient for integration
-        Variables.IntegrationCoefficient =
-            this->CalculateIntegrationCoefficient(IntegrationPoints[GPoint], Variables.detJ);
+    const auto stiffness_matrix = GeoEquationOfMotionUtilities::CalculateStiffnessMatrix(
+        b_matrices, constitutive_matrices, integration_coefficients);
 
-        // Compute stiffness matrix
-        this->CalculateAndAddStiffnessMatrix(rStiffnessMatrix, Variables);
-    }
+    GeoElementUtilities::AssembleUUBlockMatrix(rStiffnessMatrix, stiffness_matrix);
 
     KRATOS_CATCH("")
 }
@@ -1248,11 +1239,9 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddStiffnessMatrix(Matr
 {
     KRATOS_TRY
 
-    noalias(rVariables.UVoigtMatrix) = prod(trans(rVariables.B), rVariables.ConstitutiveMatrix);
-    noalias(rVariables.UUMatrix) = prod(rVariables.UVoigtMatrix, rVariables.B) * rVariables.IntegrationCoefficient;
-
-    // Distribute stiffness block matrix into the elemental matrix
-    GeoElementUtilities::AssembleUUBlockMatrix(rLeftHandSideMatrix, rVariables.UUMatrix);
+    const auto stiffness_matrix = GeoEquationOfMotionUtilities::CalculateStiffnessMatrixGPoint(
+        rVariables.B, rVariables.ConstitutiveMatrix, rVariables.IntegrationCoefficient);
+    GeoElementUtilities::AssembleUUBlockMatrix(rLeftHandSideMatrix, stiffness_matrix);
 
     KRATOS_CATCH("")
 }
@@ -1348,7 +1337,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddStiffnessForce(Vecto
     noalias(rVariables.UVector) =
         -1.0 * prod(trans(rVariables.B), mStressVector[GPoint]) * rVariables.IntegrationCoefficient;
 
-    // Distribute stiffness block vector into elemental vector
     GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, rVariables.UVector);
 
     KRATOS_CATCH("")

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp

@@ -1355,24 +1355,13 @@ void SmallStrainUPwDiffOrderElement::CalculateMaterialStiffnessMatrix(MatrixType
     this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters,
                                          Variables.NuContainer, Variables.DNu_DXContainer,
                                          strain_vectors, mStressVector, constitutive_matrices);
+    const auto integration_coefficients =
+        CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJuContainer);
 
-    for (unsigned int GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) {
-        // compute element kinematics (Np, gradNpT, |J|, B, strains)
-        this->CalculateKinematics(Variables, GPoint);
-        Variables.B = b_matrices[GPoint];
-
-        // Compute infinitesimal strain
-        Variables.F                  = deformation_gradients[GPoint];
-        Variables.StrainVector       = strain_vectors[GPoint];
-        Variables.ConstitutiveMatrix = constitutive_matrices[GPoint];
-
-        // calculating weighting coefficient for integration
-        Variables.IntegrationCoefficient =
-            this->CalculateIntegrationCoefficient(IntegrationPoints[GPoint], Variables.detJ);
+    const auto stiffness_matrix = GeoEquationOfMotionUtilities::CalculateStiffnessMatrix(
+        b_matrices, constitutive_matrices, integration_coefficients);
 
-        // Contributions of material stiffness to the left hand side
-        this->CalculateAndAddStiffnessMatrix(rStiffnessMatrix, Variables);
-    }
+    GeoElementUtilities::AssembleUUBlockMatrix(rStiffnessMatrix, stiffness_matrix);
 
     KRATOS_CATCH("")
 }
@@ -1718,16 +1707,14 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddLHS(MatrixType& rLeftHandSid
 }
 
 void SmallStrainUPwDiffOrderElement::CalculateAndAddStiffnessMatrix(MatrixType& rLeftHandSideMatrix,
-                                                                    ElementVariables& rVariables) const
+                                                                    const ElementVariables& rVariables) const
 {
     KRATOS_TRY
 
-    Matrix StiffnessMatrix =
-        prod(trans(rVariables.B), Matrix(prod(rVariables.ConstitutiveMatrix, rVariables.B))) *
-        rVariables.IntegrationCoefficient;
+    const auto stiffness_matrix = GeoEquationOfMotionUtilities::CalculateStiffnessMatrixGPoint(
+        rVariables.B, rVariables.ConstitutiveMatrix, rVariables.IntegrationCoefficient);
 
-    // Distribute stiffness block matrix into the elemental matrix
-    GeoElementUtilities::AssembleUUBlockMatrix(rLeftHandSideMatrix, StiffnessMatrix);
+    GeoElementUtilities::AssembleUUBlockMatrix(rLeftHandSideMatrix, stiffness_matrix);
 
     KRATOS_CATCH("")
 }

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp

@@ -256,7 +256,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub
 
     void CalculateAndAddLHS(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables);
 
-    void CalculateAndAddStiffnessMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables) const;
+    void CalculateAndAddStiffnessMatrix(MatrixType& rLeftHandSideMatrix, const ElementVariables& rVariables) const;
 
     void CalculateAndAddCouplingMatrix(MatrixType& rLeftHandSideMatrix, const ElementVariables& rVariables);
 

applications/GeoMechanicsApplication/custom_utilities/README.md

@@ -58,15 +58,25 @@ $$M = \int_\Omega N_{u}^T \rho N_u d\Omega$$
 
 Where $\Omega$ is the domain, $N_u$ is the displacement shape function and $\rho$ is the density matrix that holds density for all directions.
 
+### Stiffness Matrix (K)
+
+The mathematical definition is:
+$$K = \int_\Omega B^T C_{constitutive} B d\Omega$$
+
+Where $\Omega$ is the domain, $B$ is the B-matrix and $C_{constitutive}$ is the constitutive matrix.
+
+
 ### Damping Matrix (D)
 
 The mathematical definition is:
 $$D = \alpha_R M + \beta_R K$$
 
 Where $M$ and $K$ are the mass and stiffness  matrices respectively and $\alpha_R$ and $\beta_R$ are the coefficients from the Rayleigh Method.
 
-File equation_of_motion_utilities.hpp includes 
+File equation_of_motion_utilities.h includes 
 -  CalculateMassMatrix function
+-  CalculateStiffnessMatrixGPoint provides a stiffness matrix for a specific integration point
+-  CalculateStiffnessMatrix provides a stiffness matrix for an element
 -  CalculateDampingMatrix function
 -  CalculateIntegrationCoefficientsInitialConfiguration function that calculates integration coefficient for all integration points
 

applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.cpp

@@ -66,4 +66,23 @@ Matrix GeoEquationOfMotionUtilities::CalculateDampingMatrix(double        Raylei
     return RayleighAlpha * rMassMatrix + RayleighBeta * rStiffnessMatrix;
 }
 
+Matrix GeoEquationOfMotionUtilities::CalculateStiffnessMatrixGPoint(const Matrix& rB,
+                                                                    const Matrix& rConstitutiveMatrix,
+                                                                    double IntegrationCoefficient)
+{
+    return prod(trans(rB), Matrix(prod(rConstitutiveMatrix, rB))) * IntegrationCoefficient;
+}
+
+Matrix GeoEquationOfMotionUtilities::CalculateStiffnessMatrix(const std::vector<Matrix>& rBs,
+                                                              const std::vector<Matrix>& rConstitutiveMatrices,
+                                                              const std::vector<double>& rIntegrationCoefficients)
+{
+    Matrix result = ZeroMatrix(rBs[0].size2(), rBs[0].size2());
+    for (unsigned int GPoint = 0; GPoint < rBs.size(); ++GPoint) {
+        result += CalculateStiffnessMatrixGPoint(rBs[GPoint], rConstitutiveMatrices[GPoint],
+                                                 rIntegrationCoefficients[GPoint]);
+    }
+    return result;
+}
+
 } /* namespace Kratos.*/

applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.h

@@ -41,5 +41,13 @@ class GeoEquationOfMotionUtilities
                                          const Matrix& rMassMatrix,
                                          const Matrix& rStiffnessMatrix);
 
+    static Matrix CalculateStiffnessMatrixGPoint(const Matrix& rB,
+                                                 const Matrix& rConstitutiveMatrix,
+                                                 double        IntegrationCoefficient);
+
+    static Matrix CalculateStiffnessMatrix(const std::vector<Matrix>& rBs,
+                                           const std::vector<Matrix>& rConstitutiveMatrices,
+                                           const std::vector<double>& rIntegrationCoefficients);
+
 }; /* Class GeoTransportEquationUtilities*/
 } /* namespace Kratos.*/

applications/GeoMechanicsApplication/tests/cpp_tests/test_equation_of_motion.cpp

@@ -158,4 +158,29 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateDampingMatrixGivesCorrectResults, KratosGeoMe
     KRATOS_CHECK_MATRIX_NEAR(damping_matrix, expected_damping_matrix, 1e-4)
 }
 
+KRATOS_TEST_CASE_IN_SUITE(CalculateStiffnessMatrixGivesCorrectResults, KratosGeoMechanicsFastSuite)
+{
+    constexpr std::size_t voigt_size = 4;
+    constexpr std::size_t n          = 3;
+
+    const Matrix        b = ScalarMatrix(voigt_size, n, 0.5);
+    std::vector<Matrix> b_matrices;
+    b_matrices.push_back(b);
+    b_matrices.push_back(b);
+
+    const Matrix        constitutive = ScalarMatrix(voigt_size, voigt_size, 2.0);
+    std::vector<Matrix> constitutive_matrices;
+    constitutive_matrices.push_back(constitutive);
+    constitutive_matrices.push_back(constitutive);
+
+    std::vector<double> integration_coefficients{1.0, 2.0};
+
+    auto stiffness_matrix = GeoEquationOfMotionUtilities::CalculateStiffnessMatrix(
+        b_matrices, constitutive_matrices, integration_coefficients);
+
+    const auto expected_stiffness_matrix = ScalarMatrix(n, n, 24.0);
+
+    KRATOS_CHECK_MATRIX_NEAR(stiffness_matrix, expected_stiffness_matrix, 1e-4)
+}
+
 } // namespace Kratos::Testing