Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[GeoMechanicsApplication] Extract a static utility function for the calculation of the Damping Matrix (D) #12368

Merged
merged 7 commits into from
May 13, 2024
Merged

[GeoMechanicsApplication] Extract a static utility function for the calculation of the Damping Matrix (D) #12368

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
8e045e7
extracted static utility function CalculateDampingMatrix
May 10, 2024
500b6dd
changes to answer the review comments
May 10, 2024
528c509
removed resize and noalias
May 10, 2024
4b87c12
put tests for mass and damping matrices into one file
May 10, 2024
5edec92
removed size/noalias in geo_structural_base_element.cpp
May 10, 2024
28420e3
removed obvious comments
May 13, 2024
eb7d6e0
checked damping_matrix only with expected_damping_matrix
May 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 12 additions & 23 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down Expand Up @@ -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("")
}
Expand Down Expand Up @@ -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());
Expand Down
24 changes: 7 additions & 17 deletions applications/GeoMechanicsApplication/custom_elements/geo_structural_base_element.cpp
View file
Open in desktop
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for eliminating the duplication; this looks much better. I have the same suggestions about code comments as for the previous file.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

removed the comments

Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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);
Comment on lines +275 to +276
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The thing that struck me here was that apparently for "geo structural" elements the element properties need not be checked for the Rayleigh coefficients (contrary to the U-Pw elements). They are always taken from the "current process info".


KRATOS_CATCH("")
}
Expand Down
26 changes: 7 additions & 19 deletions ...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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("")
}
Expand Down Expand Up @@ -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);
}

Expand Down
8 changes: 8 additions & 0 deletions applications/GeoMechanicsApplication/custom_utilities/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 8 additions & 0 deletions applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,4 +58,12 @@ Vector GeoEquationOfMotionUtilities::CalculateDetJsInitialConfiguration(const Ge
return det_Js_initial_configuration;
}

Matrix GeoEquationOfMotionUtilities::CalculateDampingMatrix(double RayleighAlpha,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we ever want to check the inputs for the damping matrix (e.g. in this case I'd expect alpha and beta to be >0). We can also just assume the numbers to be correct, but since they are direct user inputs in this case, it might be a good idea to check them?

This might also be a broader discussion, since we don't do input validation for the jsons now, we might want to do that in a central place and avoid expensive checks for each time we call this function (which is done for each time step and each element of course)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Checks should be as early as possible and once, so not here but at element initialization and once per material property/constitutive law.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fully agreed with Wijtze-Pieter.

double RayleighBeta,
const Matrix& rMassMatrix,
const Matrix& rStiffnessMatrix)
{
return RayleighAlpha * rMassMatrix + RayleighBeta * rStiffnessMatrix;
}

} /* namespace Kratos.*/
5 changes: 5 additions & 0 deletions applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.*/
40 changes: 39 additions & 1 deletion ...tion/tests/cpp_tests/test_mass_matrix.cpp → ...sts/cpp_tests/test_equation_of_motion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Excellent unit test! The only thing I'm a bit hesitant about is repeating the implementation of the damping matrix computation here. Instead, we could also express the expected matrix like you did in the unit test above this one (see expected_mass_matrix).

Copy link
Contributor

@rfaasse rfaasse May 13, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since the expected matrix only has a single value, I would not express the matrix in the same way as expected_mass_matrix. We might want to give this double just the value it should have (which I guess is 15.0) and use expected_damping_matrix = scalar_matrix(n, n, expected_matrix_value) to fill the matrix as we do now. I'm not sure this approach is needed for such a simple calculation though. Not sure what you both think?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed with Richard. Now the calculated damping matrix is compared with expected damping matrix only. Each time the expected matrix is filled with a specific single 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
Loading