Skip to content

Commit

Permalink
removed CalculateMassMatrix
Browse files Browse the repository at this point in the history
  • Loading branch information
Gennady Markelov committed Apr 19, 2024
1 parent 83f27a1 commit 122e145
Show file tree
Hide file tree
Showing 2 changed files with 3 additions and 90 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -429,87 +429,6 @@ void SmallStrainUPwDiffOrderElement::CalculateRightHandSide(VectorType& r
KRATOS_CATCH("")
}

void SmallStrainUPwDiffOrderElement::CalculateMassMatrix(MatrixType& rMassMatrix, const ProcessInfo& rCurrentProcessInfo)
{
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);

this->CalculateSoilDensity(Variables);

// 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;
}

// Distribute mass block matrix into the elemental matrix
const SizeType NumPNodes = mpPressureGeometry->PointsNumber();

const SizeType ElementSize = BlockElementSize + NumPNodes;

if (rMassMatrix.size1() != ElementSize || rMassMatrix.size2() != ElementSize)
rMassMatrix.resize(ElementSize, ElementSize, false);
noalias(rMassMatrix) = ZeroMatrix(ElementSize, ElementSize);

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);
}
}
}
}

KRATOS_CATCH("")
}

void SmallStrainUPwDiffOrderElement::CalculateDampingMatrix(MatrixType& rDampingMatrix,
const ProcessInfo& rCurrentProcessInfo)
{
Expand All @@ -522,12 +441,8 @@ void SmallStrainUPwDiffOrderElement::CalculateDampingMatrix(MatrixType& r
const SizeType NumPNodes = mpPressureGeometry->PointsNumber();

const SizeType ElementSize = NumUNodes * Dim + NumPNodes;

// Compute Mass Matrix
MatrixType MassMatrix(ElementSize, ElementSize);

// this->CalculateMassMatrix(MassMatrix, rCurrentProcessInfo);
MassMatrix = GeoTransportEquationUtilities::CalculateMassMatrix(

MatrixType MassMatrix = GeoTransportEquationUtilities::CalculateMassMatrix(
this->GetGeometry(), mpPressureGeometry, this->GetIntegrationMethod(), mpStressStatePolicy,
mRetentionLawVector, this->GetProperties(), rCurrentProcessInfo);

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -86,9 +86,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub
void CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, const ProcessInfo& rCurrentProcessInfo) override;

void CalculateRightHandSide(VectorType& rRightHandSideVector, const ProcessInfo& rCurrentProcessInfo) override;

void CalculateMassMatrix(MatrixType& rMassMatrix, const ProcessInfo& rCurrentProcessInfo) override;


void EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo&) const override;

void GetFirstDerivativesVector(Vector& rValues, int Step = 0) const override;
Expand Down

0 comments on commit 122e145

Please sign in to comment.