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.

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

[GeoMechanicsApplication] Extract computation of integration coefficients #12334

Merged
merged 10 commits into from
May 1, 2024
Original file line number Diff line number Diff line change
Expand Up @@ -528,11 +528,10 @@ void UPwBaseElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHandS
//----------------------------------------------------------------------------------------
template <unsigned int TDim, unsigned int TNumNodes>
double UPwBaseElement<TDim, TNumNodes>::CalculateIntegrationCoefficient(
const GeometryType::IntegrationPointsArrayType& IntegrationPoints, unsigned int PointNumber, double detJ)
const GeometryType::IntegrationPointType& rIntegrationPoint, double detJ)

{
return mpStressStatePolicy->CalculateIntegrationCoefficient(IntegrationPoints[PointNumber],
detJ, GetGeometry());
return mpStressStatePolicy->CalculateIntegrationCoefficient(rIntegrationPoint, detJ, GetGeometry());
}

//----------------------------------------------------------------------------------------
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -164,9 +164,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwBaseElement : public Element
bool CalculateStiffnessMatrixFlag,
bool CalculateResidualVectorFlag);

virtual double CalculateIntegrationCoefficient(const GeometryType::IntegrationPointsArrayType& IntegrationPoints,
unsigned int PointNumber,
double detJ);
virtual double CalculateIntegrationCoefficient(const GeometryType::IntegrationPointType& rIntegrationPoint,
double detJ);

void CalculateDerivativesOnInitialConfiguration(
double& detJ, Matrix& J0, Matrix& InvJ0, Matrix& DN_DX, unsigned int PointNumber) const;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -438,6 +438,13 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa

const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT);

auto integration_coefficients = std::vector<double>{};
std::transform(IntegrationPoints.begin(), IntegrationPoints.end(),
Variables.detJContainer.begin(), std::back_inserter(integration_coefficients),
[this](const auto& rIntegrationPoint, const auto& rDetJ) {
return this->CalculateIntegrationCoefficient(rIntegrationPoint, rDetJ);
});

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
// Compute Np, GradNpT, B and StrainVector
Expand Down Expand Up @@ -468,12 +475,10 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
// calculate Bulk modulus from stiffness matrix
this->InitializeBiotCoefficients(Variables, hasBiotCoefficient);

// Compute weighting coefficient for integration
Variables.IntegrationCoefficient =
this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, Variables.detJ);
Variables.IntegrationCoefficient = integration_coefficients[GPoint];

Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(
IntegrationPoints, GPoint, Variables.detJInitialConfiguration);
Comment on lines -471 to -476
Copy link
Contributor

Choose a reason for hiding this comment

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

This piece of code seems to use the same function with the same input. Therefore I think that detJ == detJInitialConfiguration. In the new implementation storing the list of detJ for the integrationpoints would then avoid the second call to CalculateIntegrationCoefficient.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I'm not sure. ElementVariables::detJ is taken from ElementVariables::detJContainer, which is computed by calling member InitializeElementVariables. ElementVariables::detJInitialConfiguration on the other hand is computed by calling member CalculateKinematics (which in turn calls member CalculateDerivativesOnInitialConfiguration). So they are probably two different things. However, it is unclear to me where ElementVariables::detJInitialConfiguration is being used and whether that is correct. To be investigated.

IntegrationPoints[GPoint], Variables.detJInitialConfiguration);

if (CalculateStiffnessMatrixFlag) {
// Contributions to the left hand side
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -217,7 +217,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateHydraulicDischarge(const P

// Compute weighting coefficient for integration
Variables.IntegrationCoefficient =
this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, Variables.detJ);
this->CalculateIntegrationCoefficient(IntegrationPoints[GPoint], Variables.detJ);

for (unsigned int node = 0; node < TNumNodes; ++node) {
double HydraulicDischarge = 0;
Expand Down Expand Up @@ -906,7 +906,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateMaterialStiffnessMatrix(Ma

// Compute weighting coefficient for integration
Variables.IntegrationCoefficient =
this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, Variables.detJ);
this->CalculateIntegrationCoefficient(IntegrationPoints[GPoint], Variables.detJ);

// Compute stiffness matrix
this->CalculateAndAddStiffnessMatrix(rStiffnessMatrix, Variables);
Expand Down Expand Up @@ -973,7 +973,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateMassMatrix(MatrixType& rMa

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

CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

Expand Down Expand Up @@ -1023,6 +1023,13 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLe

const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT);

auto integration_coefficients = std::vector<double>{};
std::transform(IntegrationPoints.begin(), IntegrationPoints.end(),
Variables.detJContainer.begin(), std::back_inserter(integration_coefficients),
[this](const auto& rIntegrationPoint, const auto& rDetJ) {
return this->CalculateIntegrationCoefficient(rIntegrationPoint, rDetJ);
});

for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
// Compute GradNpT, B and StrainVector
this->CalculateKinematics(Variables, GPoint);
Expand All @@ -1047,12 +1054,10 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLe
this->InitializeBiotCoefficients(Variables, hasBiotCoefficient);
this->CalculatePermeabilityUpdateFactor(Variables);

// Compute weighting coefficient for integration
Variables.IntegrationCoefficient =
this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, Variables.detJ);
Variables.IntegrationCoefficient = integration_coefficients[GPoint];

Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(
IntegrationPoints, GPoint, Variables.detJInitialConfiguration);
IntegrationPoints[GPoint], Variables.detJInitialConfiguration);

// Contributions to the left hand side
if (CalculateStiffnessMatrixFlag) this->CalculateAndAddLHS(rLeftHandSideMatrix, Variables);
Expand Down Expand Up @@ -1287,16 +1292,16 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddCompressibilityMatri

template <unsigned int TDim, unsigned int TNumNodes>
void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix,
ElementVariables& rVariables)
const ElementVariables& rVariables)
{
KRATOS_TRY

rVariables.PPMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
const auto permeability_matrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
rfaasse marked this conversation as resolved.
Show resolved Hide resolved
rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix,
rVariables.RelativePermeability, rVariables.PermeabilityUpdateFactor, rVariables.IntegrationCoefficient);

// Distribute permeability block matrix into the elemental matrix
GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, rVariables.PPMatrix);
GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, permeability_matrix);

KRATOS_CATCH("")
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -240,7 +240,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa

virtual void CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables);

virtual void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables);
virtual void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, const ElementVariables& rVariables);

virtual void CalculateAndAddRHS(VectorType& rRightHandSideVector, ElementVariables& rVariables, unsigned int GPoint);

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -195,7 +195,7 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateMassMatrix(Matrix

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

CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

Expand Down Expand Up @@ -1278,7 +1278,7 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateMaterialStiffness

// Compute weighting coefficient for integration
Variables.IntegrationCoefficient =
this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, detJContainer[GPoint]);
this->CalculateIntegrationCoefficient(IntegrationPoints[GPoint], detJContainer[GPoint]);

// Compute stiffness matrix
this->CalculateAndAddStiffnessMatrix(rStiffnessMatrix, Variables);
Expand Down Expand Up @@ -1336,6 +1336,13 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAll(MatrixType& r

const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT);

auto integration_coefficients = std::vector<double>{};
std::transform(IntegrationPoints.begin(), IntegrationPoints.end(),
detJContainer.begin(), std::back_inserter(integration_coefficients),
[this](const auto& rIntegrationPoint, const auto& rDetJ) {
return this->CalculateIntegrationCoefficient(rIntegrationPoint, rDetJ);
});

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
// Compute Np, StrainVector, JointWidth, GradNpT
Expand Down Expand Up @@ -1371,9 +1378,7 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAll(MatrixType& r

this->InitializeBiotCoefficients(Variables, hasBiotCoefficient);

// Compute weighting coefficient for integration
Variables.IntegrationCoefficient =
this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, detJContainer[GPoint]);
Variables.IntegrationCoefficient = integration_coefficients[GPoint];

// Contributions to the left hand side
if (CalculateStiffnessMatrixFlag) this->CalculateAndAddLHS(rLeftHandSideMatrix, Variables);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -333,6 +333,13 @@ void UPwSmallStrainLinkInterfaceElement<TDim, TNumNodes>::CalculateAll(MatrixTyp

const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT);

auto integration_coefficients = std::vector<double>{};
std::transform(IntegrationPoints.begin(), IntegrationPoints.end(),
detJContainer.begin(), std::back_inserter(integration_coefficients),
[this](const auto& rIntegrationPoint, const auto& rDetJ) {
return this->CalculateIntegrationCoefficient(rIntegrationPoint, rDetJ);
});

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
// Compute Np, StrainVector, JointWidth, GradNpT
Expand Down Expand Up @@ -362,9 +369,7 @@ void UPwSmallStrainLinkInterfaceElement<TDim, TNumNodes>::CalculateAll(MatrixTyp

this->InitializeBiotCoefficients(Variables, hasBiotCoefficient);

// Compute weighting coefficient for integration
Variables.IntegrationCoefficient =
this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, detJContainer[GPoint]);
Variables.IntegrationCoefficient = integration_coefficients[GPoint];

// Contributions to the left hand side
if (CalculateStiffnessMatrixFlag) this->CalculateAndAddLHS(rLeftHandSideMatrix, Variables);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,13 @@ void UPwUpdatedLagrangianFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& r

const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT);

auto integration_coefficients = std::vector<double>{};
std::transform(IntegrationPoints.begin(), IntegrationPoints.end(),
Variables.detJContainer.begin(), std::back_inserter(integration_coefficients),
[this](const auto& rIntegrationPoint, const auto& rDetJ) {
return this->CalculateIntegrationCoefficient(rIntegrationPoint, rDetJ);
});

// Computing in all integrations points
for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) {
// Compute element kinematics B, F, GradNpT ...
Expand Down Expand Up @@ -105,12 +112,10 @@ void UPwUpdatedLagrangianFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& r
// calculate Bulk modulus from stiffness matrix
this->InitializeBiotCoefficients(Variables, hasBiotCoefficient);

// Compute weighting coefficient for integration
Variables.IntegrationCoefficient =
this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, Variables.detJ);
Variables.IntegrationCoefficient = integration_coefficients[GPoint];

Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(
IntegrationPoints, GPoint, Variables.detJInitialConfiguration);
IntegrationPoints[GPoint], Variables.detJInitialConfiguration);

if (CalculateStiffnessMatrixFlag) {
// Contributions to stiffness matrix calculated on the reference config
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,13 @@ void UPwUpdatedLagrangianElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLef

const bool hasBiotCoefficient = this->GetProperties().Has(BIOT_COEFFICIENT);

auto integration_coefficients = std::vector<double>{};
std::transform(IntegrationPoints.begin(), IntegrationPoints.end(),
Variables.detJContainer.begin(), std::back_inserter(integration_coefficients),
[this](const auto& rIntegrationPoint, const auto& rDetJ) {
return this->CalculateIntegrationCoefficient(rIntegrationPoint, rDetJ);
});

// Computing in all integrations points
for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) {
// Compute element kinematics B, F, GradNpT ...
Expand All @@ -110,12 +117,10 @@ void UPwUpdatedLagrangianElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLef
// calculate Bulk modulus from stiffness matrix
this->InitializeBiotCoefficients(Variables, hasBiotCoefficient);

// Calculating weights for integration on the reference configuration
Variables.IntegrationCoefficient =
this->CalculateIntegrationCoefficient(IntegrationPoints, GPoint, Variables.detJ);
Variables.IntegrationCoefficient = integration_coefficients[GPoint];

Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(
IntegrationPoints, GPoint, Variables.detJInitialConfiguration);
IntegrationPoints[GPoint], Variables.detJInitialConfiguration);

if (CalculateStiffnessMatrixFlag) {
// Contributions to stiffness matrix calculated on the reference config
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -461,7 +461,7 @@ void SmallStrainUPwDiffOrderElement::CalculateMassMatrix(MatrixType& rMassMatrix

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

CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

Expand Down Expand Up @@ -1315,6 +1315,13 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi

const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT);

auto integration_coefficients = std::vector<double>{};
std::transform(IntegrationPoints.begin(), IntegrationPoints.end(),
Variables.detJuContainer.begin(), std::back_inserter(integration_coefficients),
[this](const auto& rIntegrationPoint, const auto& rDetJ) {
return this->CalculateIntegrationCoefficient(rIntegrationPoint, rDetJ);
});

for (unsigned int GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) {
// compute element kinematics (Np, gradNpT, |J|, B, strains)
this->CalculateKinematics(Variables, GPoint);
Expand All @@ -1334,12 +1341,10 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi
this->InitializeBiotCoefficients(Variables, hasBiotCoefficient);
this->CalculatePermeabilityUpdateFactor(Variables);

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

Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(
IntegrationPoints, GPoint, Variables.detJInitialConfiguration);
IntegrationPoints[GPoint], Variables.detJInitialConfiguration);

// Contributions to the left hand side
if (CalculateStiffnessMatrixFlag) this->CalculateAndAddLHS(rLeftHandSideMatrix, Variables);
Expand Down Expand Up @@ -1389,7 +1394,7 @@ void SmallStrainUPwDiffOrderElement::CalculateMaterialStiffnessMatrix(MatrixType

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

// Contributions of material stiffness to the left hand side
this->CalculateAndAddStiffnessMatrix(rStiffnessMatrix, Variables);
Expand Down Expand Up @@ -1718,10 +1723,10 @@ void SmallStrainUPwDiffOrderElement::SetConstitutiveParameters(ElementVariables&
KRATOS_CATCH("")
}

double SmallStrainUPwDiffOrderElement::CalculateIntegrationCoefficient(
const GeometryType::IntegrationPointsArrayType& IntegrationPoints, unsigned int GPoint, double detJ)
double SmallStrainUPwDiffOrderElement::CalculateIntegrationCoefficient(const GeometryType::IntegrationPointType& rIntegrationPoint,
double detJ)
{
return mpStressStatePolicy->CalculateIntegrationCoefficient(IntegrationPoints[GPoint], detJ, GetGeometry());
return mpStressStatePolicy->CalculateIntegrationCoefficient(rIntegrationPoint, detJ, GetGeometry());
}

void SmallStrainUPwDiffOrderElement::CalculateAndAddLHS(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -252,9 +252,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub
void SetConstitutiveParameters(ElementVariables& rVariables,
ConstitutiveLaw::Parameters& rConstitutiveParameters) const;

virtual double CalculateIntegrationCoefficient(const GeometryType::IntegrationPointsArrayType& IntegrationPoints,
unsigned int PointNumber,
double detJ);
virtual double CalculateIntegrationCoefficient(const GeometryType::IntegrationPointType& rIntegrationPoint,
double detJ);

void CalculateAndAddLHS(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables);

Expand Down
Loading
Loading