diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp index d4f1e98fbbba..dd51f42ae3e9 100644 --- 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::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::CalculateAll(MatrixType& rLeftHandS //---------------------------------------------------------------------------------------- template -double UPwBaseElement::CalculateIntegrationCoefficient( - const GeometryType::IntegrationPointType& rIntegrationPoint, double detJ) const +double UPwBaseElement::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 index 11632006d5ed..9661c761c267 100644 --- 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::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/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp index ad53ebc86fda..cf6596c0952f 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp +++ b/applications/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 index 80ce94a8611c..ef3ce4d2fc81 100644 --- 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 index 3eda9360e4b1..b46c7f6194cc 100644 --- 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 index 9a3b5b6ac264..847c78f4687e 100644 --- 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& 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/applications/GeoMechanicsApplication/tests/cpp_tests/test_mass_matrix.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_equation_of_motion.cpp similarity index 80% rename from applications/GeoMechanicsApplication/tests/cpp_tests/test_mass_matrix.cpp rename to applications/GeoMechanicsApplication/tests/cpp_tests/test_equation_of_motion.cpp index 3326fe2dbd3c..8f92fdba1370 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_mass_matrix.cpp +++ b/applications/GeoMechanicsApplication/tests/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 \ No newline at end of file