-
Notifications
You must be signed in to change notification settings - Fork 242
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.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We鈥檒l occasionally send you account related emails.
Already on GitHub? Sign in to your account
[GeoMechanicsApplication] Extract a static utility function for the calculation of the Damping Matrix (D) #12368
Changes from 1 commit
8e045e7
500b6dd
528c509
4b87c12
5edec92
28420e3
eb7d6e0
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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 | ||
{ | ||
|
@@ -327,26 +328,22 @@ void UPwBaseElement<TDim, TNumNodes>::CalculateDampingMatrix(MatrixType& | |
const unsigned int N_DOF = this->GetNumberOfDOF(); | ||
|
||
// Compute Mass Matrix | ||
MatrixType MassMatrix(N_DOF, N_DOF); | ||
|
||
this->CalculateMassMatrix(MassMatrix, rCurrentProcessInfo); | ||
MatrixType mass_matrix(N_DOF, N_DOF); | ||
this->CalculateMassMatrix(mass_matrix, rCurrentProcessInfo); | ||
|
||
// Compute Stiffness matrix | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Feel free to also remove this comment. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. removed |
||
MatrixType StiffnessMatrix(N_DOF, N_DOF); | ||
|
||
this->CalculateMaterialStiffnessMatrix(StiffnessMatrix, rCurrentProcessInfo); | ||
MatrixType stiffness_matrix(N_DOF, N_DOF); | ||
this->CalculateMaterialStiffnessMatrix(stiffness_matrix, rCurrentProcessInfo); | ||
|
||
// Compute Damping Matrix | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. And the same thing applies to this comment. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. removed |
||
if (rDampingMatrix.size1() != N_DOF) rDampingMatrix.resize(N_DOF, N_DOF, false); | ||
noalias(rDampingMatrix) = ZeroMatrix(N_DOF, N_DOF); | ||
|
||
const PropertiesType& rProp = this->GetProperties(); | ||
|
||
if (rProp.Has(RAYLEIGH_ALPHA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_ALPHA] * MassMatrix; | ||
else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_ALPHA] * MassMatrix; | ||
|
||
if (rProp.Has(RAYLEIGH_BETA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_BETA] * StiffnessMatrix; | ||
else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_BETA] * StiffnessMatrix; | ||
const PropertiesType& r_prop = this->GetProperties(); | ||
noalias(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); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Here, I think it might be nice to simplify by just assigning rDampingMatrix to GeoEquationOfMotionUtilities::CalculateDampingMatrix without There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Let's start removing resize and noalias ;) Indeed it is not clear a benefit of noalias nowadays. |
||
|
||
KRATOS_CATCH("") | ||
} | ||
|
@@ -527,8 +524,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()); | ||
|
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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. There was a problem hiding this comment. Choose a reason for hiding this commentThe 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 |
---|---|---|
|
@@ -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 | ||
|
@@ -268,22 +269,18 @@ void GeoStructuralBaseElement<TDim, TNumNodes>::CalculateDampingMatrix(MatrixTyp | |
// Rayleigh Method: Damping Matrix = alpha*M + beta*K | ||
|
||
// Compute Mass Matrix | ||
MatrixType MassMatrix(N_DOF_ELEMENT, N_DOF_ELEMENT); | ||
|
||
this->CalculateMassMatrix(MassMatrix, rCurrentProcessInfo); | ||
MatrixType mass_matrix(N_DOF_ELEMENT, N_DOF_ELEMENT); | ||
this->CalculateMassMatrix(mass_matrix, rCurrentProcessInfo); | ||
|
||
// Compute Stiffness matrix | ||
MatrixType StiffnessMatrix(N_DOF_ELEMENT, N_DOF_ELEMENT); | ||
|
||
this->CalculateStiffnessMatrix(StiffnessMatrix, rCurrentProcessInfo); | ||
MatrixType stiffness_matrix(N_DOF_ELEMENT, N_DOF_ELEMENT); | ||
this->CalculateStiffnessMatrix(stiffness_matrix, 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; | ||
noalias(rDampingMatrix) = GeoEquationOfMotionUtilities::CalculateDampingMatrix( | ||
rCurrentProcessInfo[RAYLEIGH_ALPHA], rCurrentProcessInfo[RAYLEIGH_BETA], mass_matrix, stiffness_matrix); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Apparently, the damping for the "structural" elements within geo have an implementation that differs from the continuum elements. I wonder whether we should make the possibility for a material based alpha and beta possible or to avoid that effort as we hope to let the "structural" elements within geo slowly disappear. |
||
|
||
KRATOS_CATCH("") | ||
} | ||
|
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Same suggestions about the code comments. There was a problem hiding this comment. Choose a reason for hiding this commentThe 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 |
---|---|---|
|
@@ -58,4 +58,12 @@ Vector GeoEquationOfMotionUtilities::CalculateDetJsInitialConfiguration(const Ge | |
return det_Js_initial_configuration; | ||
} | ||
|
||
Matrix GeoEquationOfMotionUtilities::CalculateDampingMatrix(double RayleighAlpha, | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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 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) There was a problem hiding this comment. Choose a reason for hiding this commentThe 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. There was a problem hiding this comment. Choose a reason for hiding this commentThe 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 rMassMatrix * RayleighAlpha + rStiffnessMatrix * RayleighBeta; | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Correct, but the formula is usually written as: alphaM + betaK Is there a need for checking the sizes of M and K before doing this operation, or do we accept that will be o.k. as those are generated for the same element and therefore have the same and correct size. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Thank you. Indeed ublas provides both options. Changed it. Yes the same size. In all three cases the matrices are generated just before CalculateDampingMatrix function. |
||
} | ||
|
||
} /* namespace Kratos.*/ |
rfaasse marked this conversation as resolved.
Show resolved
Hide resolved
|
Original file line number | Diff line number | Diff line change | ||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
@@ -0,0 +1,58 @@ | ||||||||||||||||||
// KRATOS___ | ||||||||||||||||||
// // ) ) | ||||||||||||||||||
// // ___ ___ | ||||||||||||||||||
// // ____ //___) ) // ) ) | ||||||||||||||||||
// // / / // // / / | ||||||||||||||||||
// ((____/ / ((____ ((___/ / MECHANICS | ||||||||||||||||||
// | ||||||||||||||||||
// License: geo_mechanics_application/license.txt | ||||||||||||||||||
// | ||||||||||||||||||
// Main authors: Gennady Markelov | ||||||||||||||||||
// | ||||||||||||||||||
|
||||||||||||||||||
#include "custom_utilities/equation_of_motion_utilities.h" | ||||||||||||||||||
#include "includes/checks.h" | ||||||||||||||||||
#include "testing/testing.h" | ||||||||||||||||||
#include <boost/numeric/ublas/assignment.hpp> | ||||||||||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I don't see assignment being used, could we remove this include? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Thank you. It is removed. |
||||||||||||||||||
|
||||||||||||||||||
using namespace Kratos; | ||||||||||||||||||
|
||||||||||||||||||
namespace Kratos::Testing | ||||||||||||||||||
{ | ||||||||||||||||||
|
||||||||||||||||||
KRATOS_TEST_CASE_IN_SUITE(CalculateDampingMatrixGivesCorrectResults, KratosGeoMechanicsFastSuite) | ||||||||||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Very clear test! |
||||||||||||||||||
{ | ||||||||||||||||||
constexpr std::size_t n = 10; | ||||||||||||||||||
|
||||||||||||||||||
constexpr double mass_matrix_value = 10; | ||||||||||||||||||
auto mass_matrix = scalar_matrix(n, n, mass_matrix_value); | ||||||||||||||||||
|
||||||||||||||||||
constexpr double stiffness_matrix_value = 20; | ||||||||||||||||||
auto stiffness_matrix = scalar_matrix(n, n, stiffness_matrix_value); | ||||||||||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. These can be const
Suggested change
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. fixed. |
||||||||||||||||||
|
||||||||||||||||||
double rayleigh_alpha = 0.0; | ||||||||||||||||||
double rayleigh_beta = 1.0; | ||||||||||||||||||
auto damping_matrix = GeoEquationOfMotionUtilities::CalculateDampingMatrix( | ||||||||||||||||||
rayleigh_alpha, rayleigh_beta, mass_matrix, stiffness_matrix); | ||||||||||||||||||
|
||||||||||||||||||
KRATOS_CHECK_MATRIX_NEAR(damping_matrix, stiffness_matrix, 1e-4); | ||||||||||||||||||
|
||||||||||||||||||
rayleigh_alpha = 1.0; | ||||||||||||||||||
rayleigh_beta = 0.0; | ||||||||||||||||||
damping_matrix = GeoEquationOfMotionUtilities::CalculateDampingMatrix( | ||||||||||||||||||
rayleigh_alpha, rayleigh_beta, mass_matrix, stiffness_matrix); | ||||||||||||||||||
|
||||||||||||||||||
KRATOS_CHECK_MATRIX_NEAR(damping_matrix, mass_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; | ||||||||||||||||||
auto expected_damping_matrix = scalar_matrix(n, n, expected_matrix_value); | ||||||||||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Can be const
Suggested change
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. fixed |
||||||||||||||||||
|
||||||||||||||||||
KRATOS_CHECK_MATRIX_NEAR(damping_matrix, expected_damping_matrix, 1e-4) | ||||||||||||||||||
} | ||||||||||||||||||
|
||||||||||||||||||
} // namespace Kratos::Testing |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Feel free to remove this comment, since it doesn't add anything that the code already tells.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you. Indeed the comments just repeat the function names. Removed