Merge pull request #12368 from KratosMultiphysics/geo/12366-extract-u…

…tility-damping-matrix [GeoMechanicsApplication] Extract a static utility function for the calculation of the Damping Matrix (D)
KratosMultiphysics · May 13, 2024 · 8c20a67 · 8c20a67
2 parents 0435380 + eb7d6e0
commit 8c20a67
Show file tree

Hide file tree

Showing 7 changed files with 86 additions and 60 deletions.
diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp
@@ -14,6 +14,7 @@
 // Application includes
 #include "custom_elements/U_Pw_base_element.hpp"
 #include "custom_utilities/dof_utilities.h"
+#include "custom_utilities/equation_of_motion_utilities.h"
 
 namespace Kratos
 {
@@ -322,31 +323,19 @@ void UPwBaseElement<TDim, TNumNodes>::CalculateDampingMatrix(MatrixType&
 {
     KRATOS_TRY
 
-    // Rayleigh Method: Damping Matrix = alpha*M + beta*K
-
     const unsigned int N_DOF = this->GetNumberOfDOF();
 
-    // Compute Mass Matrix
-    MatrixType MassMatrix(N_DOF, N_DOF);
-
-    this->CalculateMassMatrix(MassMatrix, rCurrentProcessInfo);
-
-    // Compute Stiffness matrix
-    MatrixType StiffnessMatrix(N_DOF, N_DOF);
-
-    this->CalculateMaterialStiffnessMatrix(StiffnessMatrix, rCurrentProcessInfo);
-
-    // Compute Damping Matrix
-    if (rDampingMatrix.size1() != N_DOF) rDampingMatrix.resize(N_DOF, N_DOF, false);
-    noalias(rDampingMatrix) = ZeroMatrix(N_DOF, N_DOF);
-
-    const PropertiesType& rProp = this->GetProperties();
+    MatrixType mass_matrix(N_DOF, N_DOF);
+    this->CalculateMassMatrix(mass_matrix, rCurrentProcessInfo);
 
-    if (rProp.Has(RAYLEIGH_ALPHA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_ALPHA] * MassMatrix;
-    else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_ALPHA] * MassMatrix;
+    MatrixType stiffness_matrix(N_DOF, N_DOF);
+    this->CalculateMaterialStiffnessMatrix(stiffness_matrix, rCurrentProcessInfo);
 
-    if (rProp.Has(RAYLEIGH_BETA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_BETA] * StiffnessMatrix;
-    else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_BETA] * StiffnessMatrix;
+    const PropertiesType& r_prop = this->GetProperties();
+    rDampingMatrix               = GeoEquationOfMotionUtilities::CalculateDampingMatrix(
+        r_prop.Has(RAYLEIGH_ALPHA) ? r_prop[RAYLEIGH_ALPHA] : rCurrentProcessInfo[RAYLEIGH_ALPHA],
+        r_prop.Has(RAYLEIGH_BETA) ? r_prop[RAYLEIGH_BETA] : rCurrentProcessInfo[RAYLEIGH_BETA],
+        mass_matrix, stiffness_matrix);
 
     KRATOS_CATCH("")
 }
@@ -527,8 +516,8 @@ void UPwBaseElement<TDim, TNumNodes>::CalculateAll(MatrixType&        rLeftHandS
 
 //----------------------------------------------------------------------------------------
 template <unsigned int TDim, unsigned int TNumNodes>
-double UPwBaseElement<TDim, TNumNodes>::CalculateIntegrationCoefficient(
-    const GeometryType::IntegrationPointType& rIntegrationPoint, double detJ) const
+double UPwBaseElement<TDim, TNumNodes>::CalculateIntegrationCoefficient(const GeometryType::IntegrationPointType& rIntegrationPoint,
+                                                                        double detJ) const
 
 {
     return mpStressStatePolicy->CalculateIntegrationCoefficient(rIntegrationPoint, detJ, GetGeometry());

diff --git a/applications/GeoMechanicsApplication/custom_elements/geo_structural_base_element.cpp b/applications/GeoMechanicsApplication/custom_elements/geo_structural_base_element.cpp
@@ -13,6 +13,7 @@
 // Application includes
 #include "custom_elements/geo_structural_base_element.hpp"
 #include "custom_utilities/dof_utilities.h"
+#include "custom_utilities/equation_of_motion_utilities.h"
 #include "geo_mechanics_application_variables.h"
 
 namespace Kratos
@@ -265,25 +266,14 @@ void GeoStructuralBaseElement<TDim, TNumNodes>::CalculateDampingMatrix(MatrixTyp
 {
     KRATOS_TRY
 
-    // Rayleigh Method: Damping Matrix = alpha*M + beta*K
+    MatrixType mass_matrix(N_DOF_ELEMENT, N_DOF_ELEMENT);
+    this->CalculateMassMatrix(mass_matrix, rCurrentProcessInfo);
 
-    // Compute Mass Matrix
-    MatrixType MassMatrix(N_DOF_ELEMENT, N_DOF_ELEMENT);
+    MatrixType stiffness_matrix(N_DOF_ELEMENT, N_DOF_ELEMENT);
+    this->CalculateStiffnessMatrix(stiffness_matrix, rCurrentProcessInfo);
 
-    this->CalculateMassMatrix(MassMatrix, rCurrentProcessInfo);
-
-    // Compute Stiffness matrix
-    MatrixType StiffnessMatrix(N_DOF_ELEMENT, N_DOF_ELEMENT);
-
-    this->CalculateStiffnessMatrix(StiffnessMatrix, rCurrentProcessInfo);
-
-    // Compute Damping Matrix
-    if (rDampingMatrix.size1() != N_DOF_ELEMENT)
-        rDampingMatrix.resize(N_DOF_ELEMENT, N_DOF_ELEMENT, false);
-    noalias(rDampingMatrix) = ZeroMatrix(N_DOF_ELEMENT, N_DOF_ELEMENT);
-
-    noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_ALPHA] * MassMatrix;
-    noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_BETA] * StiffnessMatrix;
+    rDampingMatrix = GeoEquationOfMotionUtilities::CalculateDampingMatrix(
+        rCurrentProcessInfo[RAYLEIGH_ALPHA], rCurrentProcessInfo[RAYLEIGH_BETA], mass_matrix, stiffness_matrix);
 
     KRATOS_CATCH("")
 }

diff --git a/...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp b/...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp
@@ -478,33 +478,21 @@ void SmallStrainUPwDiffOrderElement::CalculateDampingMatrix(MatrixType&        r
 {
     KRATOS_TRY
 
-    // Rayleigh Method: Damping Matrix = alpha*M + beta*K
-
     const GeometryType&   r_geom = GetGeometry();
     const PropertiesType& r_prop = this->GetProperties();
     const SizeType        element_size =
         r_geom.PointsNumber() * r_geom.WorkingSpaceDimension() + mpPressureGeometry->PointsNumber();
 
-    // Compute Mass Matrix
     MatrixType mass_matrix = ZeroMatrix(element_size, element_size);
     this->CalculateMassMatrix(mass_matrix, rCurrentProcessInfo);
 
-    // Compute Stiffness matrix
-    MatrixType StiffnessMatrix(element_size, element_size);
-
-    this->CalculateMaterialStiffnessMatrix(StiffnessMatrix, rCurrentProcessInfo);
-
-    // Compute Damping Matrix
-    if (rDampingMatrix.size1() != element_size)
-        rDampingMatrix.resize(element_size, element_size, false);
-    noalias(rDampingMatrix) = ZeroMatrix(element_size, element_size);
+    MatrixType stiffness_matrix(element_size, element_size);
+    this->CalculateMaterialStiffnessMatrix(stiffness_matrix, rCurrentProcessInfo);
 
-    if (r_prop.Has(RAYLEIGH_ALPHA)) noalias(rDampingMatrix) += r_prop[RAYLEIGH_ALPHA] * mass_matrix;
-    else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_ALPHA] * mass_matrix;
-
-    if (r_prop.Has(RAYLEIGH_BETA))
-        noalias(rDampingMatrix) += r_prop[RAYLEIGH_BETA] * StiffnessMatrix;
-    else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_BETA] * StiffnessMatrix;
+    rDampingMatrix = GeoEquationOfMotionUtilities::CalculateDampingMatrix(
+        r_prop.Has(RAYLEIGH_ALPHA) ? r_prop[RAYLEIGH_ALPHA] : rCurrentProcessInfo[RAYLEIGH_ALPHA],
+        r_prop.Has(RAYLEIGH_BETA) ? r_prop[RAYLEIGH_BETA] : rCurrentProcessInfo[RAYLEIGH_BETA],
+        mass_matrix, stiffness_matrix);
 
     KRATOS_CATCH("")
 }
@@ -2036,7 +2024,7 @@ Vector SmallStrainUPwDiffOrderElement::CalculateStrain(const Matrix& rDeformatio
         const SizeType VoigtSize = (Dim == N_DIM_3D ? VOIGT_SIZE_3D : VOIGT_SIZE_2D_PLANE_STRAIN);
         return StressStrainUtilities::CalculateHenckyStrain(rDeformationGradient, VoigtSize);
     }
-    
+
     return this->CalculateCauchyStrain(rB, rDisplacements);
 }
 

diff --git a/applications/GeoMechanicsApplication/custom_utilities/README.md b/applications/GeoMechanicsApplication/custom_utilities/README.md
@@ -58,8 +58,16 @@ $$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.
 
+### 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 
 -  CalculateMassMatrix function
+-  CalculateDampingMatrix function
 -  CalculateIntegrationCoefficientsInitialConfiguration function that calculates integration coefficient for all integration points
 
 ## Stress strain utilities

diff --git a/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.cpp b/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.cpp
@@ -58,4 +58,12 @@ Vector GeoEquationOfMotionUtilities::CalculateDetJsInitialConfiguration(const Ge
     return det_Js_initial_configuration;
 }
 
+Matrix GeoEquationOfMotionUtilities::CalculateDampingMatrix(double        RayleighAlpha,
+                                                            double        RayleighBeta,
+                                                            const Matrix& rMassMatrix,
+                                                            const Matrix& rStiffnessMatrix)
+{
+    return RayleighAlpha * rMassMatrix + RayleighBeta * rStiffnessMatrix;
+}
+
 } /* namespace Kratos.*/
diff --git a/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.h b/applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.h
@@ -36,5 +36,10 @@ class GeoEquationOfMotionUtilities
     static Vector CalculateDetJsInitialConfiguration(const Geometry<Node>& rGeom,
                                                      const GeometryData::IntegrationMethod IntegrationMethod);
 
+    static Matrix CalculateDampingMatrix(double        RayleighAlpha,
+                                         double        RayleighBeta,
+                                         const Matrix& rMassMatrix,
+                                         const Matrix& rStiffnessMatrix);
+
 }; /* Class GeoTransportEquationUtilities*/
 } /* namespace Kratos.*/
diff --git a/...tion/tests/cpp_tests/test_mass_matrix.cpp → ...sts/cpp_tests/test_equation_of_motion.cpp b/...tion/tests/cpp_tests/test_mass_matrix.cpp → ...sts/cpp_tests/test_equation_of_motion.cpp
@@ -10,7 +10,6 @@
 //  Main authors:    Gennady Markelov
 //
 
-#include "custom_elements/plane_strain_stress_state.h"
 #include "custom_utilities/element_utilities.hpp"
 #include "custom_utilities/equation_of_motion_utilities.h"
 #include "testing/testing.h"
@@ -120,4 +119,43 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateMassMatrix3D4NGivesCorrectResults, KratosGeoM
     KRATOS_CHECK_MATRIX_NEAR(mass_matrix, expected_mass_matrix, 1e-4)
 }
 
+KRATOS_TEST_CASE_IN_SUITE(CalculateDampingMatrixGivesCorrectResults, KratosGeoMechanicsFastSuite)
+{
+    constexpr std::size_t n = 10;
+
+    constexpr double mass_matrix_value = 10;
+    const auto       mass_matrix       = scalar_matrix(n, n, mass_matrix_value);
+
+    constexpr double stiffness_matrix_value = 20;
+    const auto       stiffness_matrix       = scalar_matrix(n, n, stiffness_matrix_value);
+
+    double rayleigh_alpha = 0.0;
+    double rayleigh_beta  = 1.0;
+    auto   damping_matrix = GeoEquationOfMotionUtilities::CalculateDampingMatrix(
+        rayleigh_alpha, rayleigh_beta, mass_matrix, stiffness_matrix);
+
+    auto expected_damping_matrix = scalar_matrix(n, n, stiffness_matrix_value);
+
+    KRATOS_CHECK_MATRIX_NEAR(damping_matrix, expected_damping_matrix, 1e-4)
+
+    rayleigh_alpha = 1.0;
+    rayleigh_beta  = 0.0;
+    damping_matrix = GeoEquationOfMotionUtilities::CalculateDampingMatrix(
+        rayleigh_alpha, rayleigh_beta, mass_matrix, stiffness_matrix);
+
+    expected_damping_matrix = scalar_matrix(n, n, mass_matrix_value);
+
+    KRATOS_CHECK_MATRIX_NEAR(damping_matrix, expected_damping_matrix, 1e-4)
+
+    rayleigh_alpha = 0.5;
+    rayleigh_beta  = 0.5;
+    damping_matrix = GeoEquationOfMotionUtilities::CalculateDampingMatrix(
+        rayleigh_alpha, rayleigh_beta, mass_matrix, stiffness_matrix);
+
+    const double expected_matrix_value = rayleigh_alpha * mass_matrix_value + rayleigh_beta * stiffness_matrix_value;
+    expected_damping_matrix = scalar_matrix(n, n, expected_matrix_value);
+
+    KRATOS_CHECK_MATRIX_NEAR(damping_matrix, expected_damping_matrix, 1e-4)
+}
+
 } // namespace Kratos::Testing