Skip to content

Commit

Permalink
Merge pull request #12299 from KratosMultiphysics/geo/12279-extract-f…
Browse files Browse the repository at this point in the history 
…unction-mass-matrix

[GeoMechanicsApplication] Extract a static utility function for the calculation of the Mass Matrix (M)
  • Loading branch information
@markelov208
markelov208 authored May 10, 2024
2 parents 71872f4 + cf9f0d9 commit 9c4a7b0
Show file tree
Hide file tree
Showing 15 changed files with 728 additions and 398 deletions.
77 changes: 30 additions & 47 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@

// Application includes
#include "custom_elements/U_Pw_small_strain_element.hpp"
#include "custom_utilities/equation_of_motion_utilities.h"
#include "custom_utilities/transport_equation_utilities.hpp"

namespace Kratos
Expand Down Expand Up @@ -945,50 +946,33 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateMassMatrix(MatrixType& rMa
KRATOS_TRY

const IndexType N_DOF = this->GetNumberOfDOF();

if (rMassMatrix.size1() != N_DOF) rMassMatrix.resize(N_DOF, N_DOF, false);
noalias(rMassMatrix) = ZeroMatrix(N_DOF, N_DOF);

const GeometryType& rGeom = this->GetGeometry();
const GeometryType::IntegrationPointsArrayType& IntegrationPoints =
rGeom.IntegrationPoints(this->mThisIntegrationMethod);
const IndexType NumGPoints = IntegrationPoints.size();
const GeometryType& r_geom = this->GetGeometry();
const auto integration_method = this->GetIntegrationMethod();
const GeometryType::IntegrationPointsArrayType& integration_points =
r_geom.IntegrationPoints(integration_method);
const auto N_container = r_geom.ShapeFunctionsValues(integration_method);

ElementVariables Variables;
this->InitializeElementVariables(Variables, rCurrentProcessInfo);

// Create general parameters of retention law
RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);

// Defining shape functions at all integration points
// Defining necessary variables
BoundedMatrix<double, TDim, TNumNodes * TDim> Nut = ZeroMatrix(TDim, TNumNodes * TDim);
BoundedMatrix<double, TDim, TNumNodes * TDim> AuxDensityMatrix = ZeroMatrix(TDim, TNumNodes * TDim);
BoundedMatrix<double, TDim, TDim> DensityMatrix = ZeroMatrix(TDim, TDim);

for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
GeoElementUtilities::CalculateNuMatrix<TDim, TNumNodes>(Nut, Variables.NContainer, GPoint);

Matrix J0, InvJ0, DNu_DX0;
this->CalculateDerivativesOnInitialConfiguration(Variables.detJInitialConfiguration, J0,
InvJ0, DNu_DX0, GPoint);

// Calculating weighting coefficient for integration
Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(
IntegrationPoints[GPoint], Variables.detJInitialConfiguration);
const auto degrees_saturation = GeoTransportEquationUtilities::CalculateDegreesSaturation(
this->GetPressureSolutionVector(), N_container, mRetentionLawVector, this->GetProperties(),
rCurrentProcessInfo);

CalculateRetentionResponse(Variables, RetentionParameters, GPoint);
const auto solid_densities =
GeoTransportEquationUtilities::CalculateSoilDensities(degrees_saturation, this->GetProperties());

this->CalculateSoilDensity(Variables);
const auto det_Js_initial_configuration =
GeoEquationOfMotionUtilities::CalculateDetJsInitialConfiguration(r_geom, integration_method);

GeoElementUtilities::AssembleDensityMatrix(DensityMatrix, Variables.Density);
const auto integration_coefficients =
this->CalculateIntegrationCoefficients(integration_points, det_Js_initial_configuration);

noalias(AuxDensityMatrix) = prod(DensityMatrix, Nut);
const auto mass_matrix_u = GeoEquationOfMotionUtilities::CalculateMassMatrix(
r_geom.WorkingSpaceDimension(), r_geom.PointsNumber(), integration_points.size(),
r_geom.ShapeFunctionsValues(integration_method), solid_densities, integration_coefficients);

// Adding contribution to Mass matrix
GeoElementUtilities::AssembleUUBlockMatrix(
rMassMatrix, prod(trans(Nut), AuxDensityMatrix) * Variables.IntegrationCoefficientInitialConfiguration);
}
GeoElementUtilities::AssembleUUBlockMatrix(rMassMatrix, mass_matrix_u);

KRATOS_CATCH("")
}
Expand Down Expand Up @@ -1417,23 +1401,13 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddMixBodyForce(VectorT
KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
void UPwSmallStrainElement<TDim, TNumNodes>::CalculateSoilDensity(ElementVariables& rVariables)
{
KRATOS_TRY

rVariables.Density = rVariables.DegreeOfSaturation * rVariables.Porosity * rVariables.FluidDensity +
(1.0 - rVariables.Porosity) * rVariables.SolidDensity;

KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
void UPwSmallStrainElement<TDim, TNumNodes>::CalculateSoilGamma(ElementVariables& rVariables)
{
KRATOS_TRY

this->CalculateSoilDensity(rVariables);
rVariables.Density = GeoTransportEquationUtilities::CalculateSoilDensity(
rVariables.DegreeOfSaturation, this->GetProperties());

noalias(rVariables.SoilGamma) = rVariables.Density * rVariables.BodyAcceleration;

Expand Down Expand Up @@ -1909,6 +1883,15 @@ void UPwSmallStrainElement<3, 8>::CalculateExtrapolationMatrix(BoundedMatrix<dou
rExtrapolationMatrix(7, 7) = 2.549038105676658;
}

template <unsigned int TDim, unsigned int TNumNodes>
Vector UPwSmallStrainElement<TDim, TNumNodes>::GetPressureSolutionVector()
{
Vector result(TNumNodes);
std::transform(this->GetGeometry().begin(), this->GetGeometry().end(), result.begin(),
[](const auto& node) { return node.FastGetSolutionStepValue(WATER_PRESSURE); });
return result;
}

template class UPwSmallStrainElement<2, 3>;
template class UPwSmallStrainElement<2, 4>;
template class UPwSmallStrainElement<3, 4>;
Expand Down
9 changes: 5 additions & 4 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,15 +295,16 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa

void CalculateExtrapolationMatrix(BoundedMatrix<double, TNumNodes, TNumNodes>& rExtrapolationMatrix);

void ResetHydraulicDischarge();
void CalculateHydraulicDischarge(const ProcessInfo& rCurrentProcessInfo);
void CalculateSoilGamma(ElementVariables& rVariables);
virtual void CalculateSoilDensity(ElementVariables& rVariables);
void ResetHydraulicDischarge();
void CalculateHydraulicDischarge(const ProcessInfo& rCurrentProcessInfo);
void CalculateSoilGamma(ElementVariables& rVariables);

virtual void CalculateAndAddGeometricStiffnessMatrix(MatrixType& rLeftHandSideMatrix,
ElementVariables& rVariables,
unsigned int GPoint);

VectorType GetPressureSolutionVector();

private:
friend class Serializer;

Expand Down
17 changes: 4 additions & 13 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,8 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateMassMatrix(Matrix

CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

this->CalculateSoilDensity(Variables);
Variables.Density = GeoTransportEquationUtilities::CalculateSoilDensity(
Variables.DegreeOfSaturation, this->GetProperties());

GeoElementUtilities::AssembleDensityMatrix(DensityMatrix, Variables.Density);

Expand Down Expand Up @@ -1988,23 +1989,13 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAndAddMixBodyForc
KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateSoilDensity(InterfaceElementVariables& rVariables)
{
KRATOS_TRY

rVariables.Density = rVariables.DegreeOfSaturation * rVariables.Porosity * rVariables.FluidDensity +
(1.0 - rVariables.Porosity) * rVariables.SolidDensity;

KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateSoilGamma(InterfaceElementVariables& rVariables)
{
KRATOS_TRY

this->CalculateSoilDensity(rVariables);
rVariables.Density = GeoTransportEquationUtilities::CalculateSoilDensity(
rVariables.DegreeOfSaturation, this->GetProperties());

noalias(rVariables.SoilGamma) = rVariables.Density * rVariables.BodyAcceleration;

Expand Down
4 changes: 1 addition & 3 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,9 +295,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainInterfaceElement
unsigned int GPoint);

void CalculateSoilGamma(InterfaceElementVariables& rVariables);

void CalculateSoilDensity(InterfaceElementVariables& rVariables);


void SetConstitutiveParameters(InterfaceElementVariables& rVariables,
ConstitutiveLaw::Parameters& rConstitutiveParameters);

Expand Down
146 changes: 49 additions & 97 deletions ...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
#include "custom_elements/small_strain_U_Pw_diff_order_element.hpp"
#include "custom_utilities/dof_utilities.h"
#include "custom_utilities/element_utilities.hpp"
#include "custom_utilities/equation_of_motion_utilities.h"
#include "custom_utilities/transport_equation_utilities.hpp"

namespace Kratos
Expand Down Expand Up @@ -436,79 +437,35 @@ void SmallStrainUPwDiffOrderElement::CalculateMassMatrix(MatrixType& rMassMatrix
{
KRATOS_TRY

const GeometryType& rGeom = GetGeometry();
const SizeType Dim = rGeom.WorkingSpaceDimension();
const SizeType NumUNodes = rGeom.PointsNumber();
const SizeType BlockElementSize = NumUNodes * Dim;
const GeometryType::IntegrationPointsArrayType& IntegrationPoints =
rGeom.IntegrationPoints(this->GetIntegrationMethod());

ElementVariables Variables;
this->InitializeElementVariables(Variables, rCurrentProcessInfo);

// create general parameters of retention law
RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);

Matrix MassMatrixContribution = ZeroMatrix(BlockElementSize, BlockElementSize);

// Defining shape functions and the determinant of the jacobian at all integration points

// Loop over integration points
Matrix Nu = ZeroMatrix(Dim, NumUNodes * Dim);
Matrix AuxDensityMatrix = ZeroMatrix(Dim, NumUNodes * Dim);
Matrix DensityMatrix = ZeroMatrix(Dim, Dim);

for (unsigned int GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) {
// compute element kinematics (Np, gradNpT, |J|, B)
this->CalculateKinematics(Variables, GPoint);

// calculating weighting coefficient for integration
Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(
IntegrationPoints[GPoint], Variables.detJInitialConfiguration);

CalculateRetentionResponse(Variables, RetentionParameters, GPoint);
const GeometryType& r_geom = GetGeometry();
const auto integration_method = this->GetIntegrationMethod();
const GeometryType::IntegrationPointsArrayType& integration_points =
r_geom.IntegrationPoints(integration_method);
const MatrixType Np_container = mpPressureGeometry->ShapeFunctionsValues(integration_method);
const PropertiesType& r_prop = this->GetProperties();

this->CalculateSoilDensity(Variables);
const auto degrees_saturation = GeoTransportEquationUtilities::CalculateDegreesSaturation(
this->GetPressureSolutionVector(), Np_container, mRetentionLawVector, r_prop, rCurrentProcessInfo);

// Setting the shape function matrix
SizeType Index = 0;
for (SizeType i = 0; i < NumUNodes; ++i) {
for (SizeType iDim = 0; iDim < Dim; ++iDim) {
Nu(iDim, Index++) = Variables.Nu(i);
}
}

GeoElementUtilities::AssembleDensityMatrix(DensityMatrix, Variables.Density);

noalias(AuxDensityMatrix) = prod(DensityMatrix, Nu);

// Adding contribution to Mass matrix
noalias(MassMatrixContribution) +=
prod(trans(Nu), AuxDensityMatrix) * Variables.IntegrationCoefficientInitialConfiguration;
}
const auto solid_densities =
GeoTransportEquationUtilities::CalculateSoilDensities(degrees_saturation, r_prop);

// Distribute mass block matrix into the elemental matrix
const SizeType NumPNodes = mpPressureGeometry->PointsNumber();
const auto det_Js_initial_configuration =
GeoEquationOfMotionUtilities::CalculateDetJsInitialConfiguration(r_geom, integration_method);

const SizeType ElementSize = BlockElementSize + NumPNodes;
const auto integration_coefficients =
CalculateIntegrationCoefficients(integration_points, det_Js_initial_configuration);

if (rMassMatrix.size1() != ElementSize || rMassMatrix.size2() != ElementSize)
rMassMatrix.resize(ElementSize, ElementSize, false);
noalias(rMassMatrix) = ZeroMatrix(ElementSize, ElementSize);
const auto mass_matrix_u = GeoEquationOfMotionUtilities::CalculateMassMatrix(
r_geom.WorkingSpaceDimension(), r_geom.PointsNumber(), integration_points.size(),
r_geom.ShapeFunctionsValues(integration_method), solid_densities, integration_coefficients);

for (SizeType i = 0; i < NumUNodes; ++i) {
SizeType Index_i = i * Dim;

for (SizeType j = 0; j < NumUNodes; ++j) {
SizeType Index_j = j * Dim;
for (SizeType idim = 0; idim < Dim; ++idim) {
for (SizeType jdim = 0; jdim < Dim; ++jdim) {
rMassMatrix(Index_i + idim, Index_j + jdim) +=
MassMatrixContribution(Index_i + idim, Index_j + jdim);
}
}
}
}
const SizeType element_size =
r_geom.PointsNumber() * r_geom.WorkingSpaceDimension() + mpPressureGeometry->PointsNumber();
if (rMassMatrix.size1() != element_size || rMassMatrix.size2() != element_size)
rMassMatrix.resize(element_size, element_size, false);
noalias(rMassMatrix) = ZeroMatrix(element_size, element_size);
GeoElementUtilities::AssembleUUBlockMatrix(rMassMatrix, mass_matrix_u);

KRATOS_CATCH("")
}
Expand All @@ -519,34 +476,31 @@ void SmallStrainUPwDiffOrderElement::CalculateDampingMatrix(MatrixType& r
KRATOS_TRY

// Rayleigh Method: Damping Matrix = alpha*M + beta*K
const GeometryType& rGeom = GetGeometry();
const SizeType Dim = rGeom.WorkingSpaceDimension();
const SizeType NumUNodes = rGeom.PointsNumber();
const SizeType NumPNodes = mpPressureGeometry->PointsNumber();

const SizeType ElementSize = NumUNodes * Dim + NumPNodes;
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 MassMatrix(ElementSize, ElementSize);

this->CalculateMassMatrix(MassMatrix, rCurrentProcessInfo);
MatrixType mass_matrix = ZeroMatrix(element_size, element_size);
this->CalculateMassMatrix(mass_matrix, rCurrentProcessInfo);

// Compute Stiffness matrix
MatrixType StiffnessMatrix(ElementSize, ElementSize);
MatrixType StiffnessMatrix(element_size, element_size);

this->CalculateMaterialStiffnessMatrix(StiffnessMatrix, rCurrentProcessInfo);

// Compute Damping Matrix
if (rDampingMatrix.size1() != ElementSize)
rDampingMatrix.resize(ElementSize, ElementSize, false);
noalias(rDampingMatrix) = ZeroMatrix(ElementSize, ElementSize);

const PropertiesType& rProp = this->GetProperties();
if (rDampingMatrix.size1() != element_size)
rDampingMatrix.resize(element_size, element_size, false);
noalias(rDampingMatrix) = ZeroMatrix(element_size, element_size);

if (rProp.Has(RAYLEIGH_ALPHA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_ALPHA] * MassMatrix;
else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_ALPHA] * MassMatrix;
if (r_prop.Has(RAYLEIGH_ALPHA)) noalias(rDampingMatrix) += r_prop[RAYLEIGH_ALPHA] * mass_matrix;
else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_ALPHA] * mass_matrix;

if (rProp.Has(RAYLEIGH_BETA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_BETA] * StiffnessMatrix;
if (r_prop.Has(RAYLEIGH_BETA))
noalias(rDampingMatrix) += r_prop[RAYLEIGH_BETA] * StiffnessMatrix;
else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_BETA] * StiffnessMatrix;

KRATOS_CATCH("")
Expand Down Expand Up @@ -1902,7 +1856,8 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddMixBodyForce(VectorType& rRi
const SizeType Dim = rGeom.WorkingSpaceDimension();
const SizeType NumUNodes = rGeom.PointsNumber();

this->CalculateSoilDensity(rVariables);
rVariables.Density = GeoTransportEquationUtilities::CalculateSoilDensity(
rVariables.DegreeOfSaturation, this->GetProperties());

Vector BodyAcceleration = ZeroVector(Dim);
SizeType Index = 0;
Expand All @@ -1924,18 +1879,6 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddMixBodyForce(VectorType& rRi
KRATOS_CATCH("")
}

void SmallStrainUPwDiffOrderElement::CalculateSoilDensity(ElementVariables& rVariables)
{
KRATOS_TRY

const PropertiesType& rProp = this->GetProperties();

rVariables.Density = (rVariables.DegreeOfSaturation * rProp[POROSITY] * rProp[DENSITY_WATER]) +
(1.0 - rProp[POROSITY]) * rProp[DENSITY_SOLID];

KRATOS_CATCH("")
}

void SmallStrainUPwDiffOrderElement::CalculateAndAddCouplingTerms(VectorType& rRightHandSideVector,
ElementVariables& rVariables)
{
Expand Down Expand Up @@ -2180,4 +2123,13 @@ const StressStatePolicy& SmallStrainUPwDiffOrderElement::GetStressStatePolicy()
return *mpStressStatePolicy;
}

Vector SmallStrainUPwDiffOrderElement::GetPressureSolutionVector()
{
Vector result(mpPressureGeometry->PointsNumber());
std::transform(this->GetGeometry().begin(),
this->GetGeometry().begin() + mpPressureGeometry->PointsNumber(), result.begin(),
[](const auto& node) { return node.FastGetSolutionStepValue(WATER_PRESSURE); });
return result;
}

} // Namespace Kratos
Loading

0 comments on commit 9c4a7b0

Please sign in to comment.