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

[GeoMechanicsApplication] Extract biot modulus inverse list #12402

Merged
merged 20 commits into from
May 29, 2024
Merged

[GeoMechanicsApplication] Extract biot modulus inverse list #12402

Show file tree
Hide file tree
Changes from 10 commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
4d36bb7
Extracted a couple of lists for saturation and biot moduli (first imp…
rfaasse May 24, 2024
8f9a230
Extract some diff order functions for saturation lists
rfaasse May 24, 2024
e98bb12
Removed some calculations from the retention response function
rfaasse May 24, 2024
9cfb3e7
Removed detF from Variables
rfaasse May 24, 2024
688279b
Some minor rearrangement
rfaasse May 24, 2024
b578e71
Same exercise as for diff order elements
rfaasse May 24, 2024
c0ddbe1
Changed unit tests to use the list version instead of the single Calc…
rfaasse May 24, 2024
e6449a2
Merge remote-tracking branch 'origin/master' into geo/12319-extract-b…
rfaasse May 24, 2024
8596047
Add missing include in steady_state_Pw_element.cpp
rfaasse May 24, 2024
e41be53
Merge remote-tracking branch 'origin/master' into geo/12319-extract-b…
rfaasse May 27, 2024
d94aab2
Using auto to remove clang tidy warnings
rfaasse May 28, 2024
f723c2b
Using auto to remove clang tidy warnings and changed some arguments t…
rfaasse May 28, 2024
b0b6ba0
Added [[nodiscard]]
rfaasse May 28, 2024
56c4d94
Merge remote-tracking branch 'origin/master' into geo/12319-extract-b…
rfaasse May 28, 2024
82e1dc5
Clean up some includes in the small_strain_U_Pw_diff_order_element.cp…
rfaasse May 28, 2024
c04d7fa
Add some includes explicitly
rfaasse May 28, 2024
e772c29
Changed signature of new function
rfaasse May 28, 2024
ff2db9c
used calculate degrees of saturation of the element and removed the u…
rfaasse May 28, 2024
0f30211
Removed detF for normal element
rfaasse May 28, 2024
5f1e343
Added error if vector lengths are incompatible
rfaasse May 29, 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
11 changes: 9 additions & 2 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -456,6 +456,13 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
strain_vectors, mStressVector, constitutive_matrices);
const auto biot_coefficients =
GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, Prop);
const auto fluid_pressures = GeoTransportEquationUtilities::CalculateFluidPressures(
Copy link
Contributor

Choose a reason for hiding this comment

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

A very minor thing, which we can also do later: if you would move the calculation of the fluid pressures one line up, you could also pass it to CalculateRelativePermeabilityValues.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Will do in the next PR

Variables.NContainer, Variables.PressureVector);
const auto degrees_of_saturation = this->CalculateDegreesOfSaturation(fluid_pressures, RetentionParameters);
const auto derivatives_of_saturation =
this->CalculateDerivativesOfSaturation(fluid_pressures, RetentionParameters);
const auto biot_moduli_inverse = GeoTransportEquationUtilities::CalculateInverseBiotModuli(
biot_coefficients, degrees_of_saturation, derivatives_of_saturation, Prop);

Copy link
Contributor

Choose a reason for hiding this comment

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

This part is duplicated several times. It has to be shifted to a utilility. I think it can be done in a seperate PR.

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, but as discussed it's better to do it in a separate PR (same for the other comments)

for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
this->CalculateKinematics(Variables, GPoint);
Expand All @@ -474,8 +481,8 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
FICVariables.ShearModulus = CalculateShearModulus(Variables.ConstitutiveMatrix);

Variables.BiotCoefficient = biot_coefficients[GPoint];
Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse(
Variables.BiotCoefficient, Variables.DegreeOfSaturation, Variables.DerivativeOfSaturation, Prop);
Variables.BiotModulusInverse = biot_moduli_inverse[GPoint];
Variables.DegreeOfSaturation = degrees_of_saturation[GPoint];
rfaasse marked this conversation as resolved.
Show resolved Hide resolved

Variables.IntegrationCoefficient = integration_coefficients[GPoint];

Expand Down
52 changes: 42 additions & 10 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -978,6 +978,13 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLe
strain_vectors, mStressVector, constitutive_matrices);
const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients(
constitutive_matrices, this->GetProperties());
const auto fluid_pressures = GeoTransportEquationUtilities::CalculateFluidPressures(
Copy link
Contributor

Choose a reason for hiding this comment

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

Also here the calculated fluid_pressures could be passed to CalculateRelativePermeabilityValues.

Variables.NContainer, Variables.PressureVector);
const auto degrees_of_saturation = CalculateDegreesOfSaturation(fluid_pressures, RetentionParameters);
const auto derivatives_of_saturation =
CalculateDerivativesOfSaturation(fluid_pressures, RetentionParameters);
const auto biot_moduli_inverse = GeoTransportEquationUtilities::CalculateInverseBiotModuli(
biot_coefficients, degrees_of_saturation, derivatives_of_saturation, rProp);

Copy link
Contributor

Choose a reason for hiding this comment

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

Same as above comment. Can be shifted to utility.

for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
this->CalculateKinematics(Variables, GPoint);
Expand All @@ -994,9 +1001,8 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLe
CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

Variables.BiotCoefficient = biot_coefficients[GPoint];
Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse(
Variables.BiotCoefficient, Variables.DegreeOfSaturation,
Variables.DerivativeOfSaturation, this->GetProperties());
Variables.BiotModulusInverse = biot_moduli_inverse[GPoint];
Variables.DegreeOfSaturation = degrees_of_saturation[GPoint];
Variables.PermeabilityUpdateFactor = this->CalculatePermeabilityUpdateFactor(Variables.StrainVector);

Variables.IntegrationCoefficient = integration_coefficients[GPoint];
Expand All @@ -1015,6 +1021,36 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLe
KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDerivativesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
Copy link
Contributor

Choose a reason for hiding this comment

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

I would suggest to make RetentionParameters a local variable of this function, to keep its scope as small as possible. To create such an object, with my recent changes, all you need is a Properties instance.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done

{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Clang-tidy complains about this because it thinks the retention law is copied (it doesn't recognize this is a pointer). I could make it const auto&, but not sure if that makes it more clear.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Used auto to remove the clang-tidy issue

RetentionParameters.SetFluidPressure(fluid_pressure);
return pRetentionLaw->CalculateDerivativeOfSaturation(RetentionParameters);
});

return result;
}

template <unsigned int TDim, unsigned int TNumNodes>
std::vector<double> UPwSmallStrainElement<TDim, TNumNodes>::CalculateDegreesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
Copy link
Contributor

Choose a reason for hiding this comment

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

Also here, I would suggest to change RetentionParameters to a local variable rather than it being a function parameter (to reduce its scope as much as possible).

Also, please be aware that GeoTransportEquationUtilities::CalculateDegreesSaturation implements the same functionality, although it necessarily has more function parameters than your solution. I would suggest to replace that one by your new member function, since it seems more natural to have it in the element itself.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done, as per your suggestion I removed the GeoTransportEquationUtilities::CalculateDegreesSaturation and used the new function instead. Although there are then 2 versions of the function, we hope this duplication will be removed soon by merging the diff-order/non-diff order branch.

FYI @markelov208 since you separated out the utility. Please let me know if you have any reservations with this approach

{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
RetentionParameters.SetFluidPressure(fluid_pressure);
return pRetentionLaw->CalculateSaturation(RetentionParameters);
});

return result;
}
Copy link
Contributor Author

Choose a reason for hiding this comment

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

It's a shame this is a 100% duplication with diff order. We could extract a utility, but we'd need to supply the vector of retention laws. It's still open for discussion though, I'm doubting between the two possibilities

Copy link
Contributor

Choose a reason for hiding this comment

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

I've had the exact same struggle for the calculation of the relative permeability values. I have also opted for keeping the calculation in the element itself, since you need some information that is maintained by the element. Since our intention is to ultimately remove the different branches in the element hierarchy, I'm willing to accept the technical debt that arises from the code duplication.

Copy link
Contributor

Choose a reason for hiding this comment

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

As you mentioned, this part is also duplication, and it has to be shifted to a utility.


template <unsigned int TDim, unsigned int TNumNodes>
std::vector<array_1d<double, TDim>> UPwSmallStrainElement<TDim, TNumNodes>::CalculateFluidFluxes(
const std::vector<double>& rPermeabilityUpdateFactors, const ProcessInfo& rCurrentProcessInfo)
Expand Down Expand Up @@ -1129,10 +1165,9 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeElementVariables(ElementV
rVariables.UVoigtMatrix.resize(TNumNodes * TDim, VoigtSize, false);

// Retention law
rVariables.DegreeOfSaturation = 1.0;
rVariables.DerivativeOfSaturation = 0.0;
rVariables.RelativePermeability = 1.0;
rVariables.BishopCoefficient = 1.0;
rVariables.DegreeOfSaturation = 1.0;
rVariables.RelativePermeability = 1.0;
rVariables.BishopCoefficient = 1.0;

KRATOS_CATCH("")
}
Expand Down Expand Up @@ -1570,9 +1605,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateRetentionResponse(ElementV
rRetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure(
rVariables.Np, rVariables.PressureVector));

rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters);
rVariables.DerivativeOfSaturation =
mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(rRetentionParameters);
rVariables.RelativePermeability =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(rRetentionParameters);
rVariables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(rRetentionParameters);
Expand Down
6 changes: 5 additions & 1 deletion applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -180,7 +180,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa

/// Retention Law parameters
double DegreeOfSaturation;
double DerivativeOfSaturation;
rfaasse marked this conversation as resolved.
Show resolved Hide resolved
double RelativePermeability;
double BishopCoefficient;

Expand Down Expand Up @@ -281,6 +280,11 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa
RetentionLaw::Parameters& rRetentionParameters,
unsigned int GPoint);

std::vector<double> CalculateDegreesOfSaturation(const std::vector<double>& fluid_pressures,
RetentionLaw::Parameters& RetentionParameters);
std::vector<double> CalculateDerivativesOfSaturation(const std::vector<double>& fluid_pressures,
RetentionLaw::Parameters& RetentionParameters);

///
/// \brief This function calculates the constitutive matrices, stresses and strains depending on the
/// constitutive parameters. Note that depending on the settings in the rConstitutiveParameters
Expand Down
12 changes: 9 additions & 3 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,13 @@ void UPwUpdatedLagrangianFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& r
strain_vectors, mStressVector, constitutive_matrices);
const auto biot_coefficients =
GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, Prop);
const auto fluid_pressures = GeoTransportEquationUtilities::CalculateFluidPressures(
Variables.NContainer, Variables.PressureVector);
const auto degrees_of_saturation = this->CalculateDegreesOfSaturation(fluid_pressures, RetentionParameters);
const auto derivatives_of_saturation =
this->CalculateDerivativesOfSaturation(fluid_pressures, RetentionParameters);
const auto biot_moduli_inverse = GeoTransportEquationUtilities::CalculateInverseBiotModuli(
biot_coefficients, degrees_of_saturation, derivatives_of_saturation, Prop);

Copy link
Contributor

Choose a reason for hiding this comment

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

Duplication, as mentioned above

// Computing in all integrations points
for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) {
Expand Down Expand Up @@ -113,9 +120,8 @@ void UPwUpdatedLagrangianFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& r
FICVariables.ShearModulus = CalculateShearModulus(Variables.ConstitutiveMatrix);

Variables.BiotCoefficient = biot_coefficients[GPoint];
Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse(
Variables.BiotCoefficient, Variables.DegreeOfSaturation,
Variables.DerivativeOfSaturation, this->GetProperties());
Variables.BiotModulusInverse = biot_moduli_inverse[GPoint];
Variables.DegreeOfSaturation = degrees_of_saturation[GPoint];

Variables.IntegrationCoefficient = integration_coefficients[GPoint];

Expand Down
12 changes: 9 additions & 3 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,13 @@ void UPwUpdatedLagrangianElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLef
strain_vectors, mStressVector, constitutive_matrices);
const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients(
constitutive_matrices, this->GetProperties());
const auto fluid_pressures = GeoTransportEquationUtilities::CalculateFluidPressures(
Variables.NContainer, Variables.PressureVector);
const auto degrees_of_saturation = this->CalculateDegreesOfSaturation(fluid_pressures, RetentionParameters);
const auto derivatives_of_saturation =
this->CalculateDerivativesOfSaturation(fluid_pressures, RetentionParameters);
const auto biot_moduli_inverse = GeoTransportEquationUtilities::CalculateInverseBiotModuli(
biot_coefficients, degrees_of_saturation, derivatives_of_saturation, this->GetProperties());

for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) {
this->CalculateKinematics(Variables, GPoint);
Expand All @@ -101,9 +108,8 @@ void UPwUpdatedLagrangianElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLef
this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

Variables.BiotCoefficient = biot_coefficients[GPoint];
Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse(
Variables.BiotCoefficient, Variables.DegreeOfSaturation,
Variables.DerivativeOfSaturation, this->GetProperties());
Variables.BiotModulusInverse = biot_moduli_inverse[GPoint];
Variables.DegreeOfSaturation = degrees_of_saturation[GPoint];

Variables.IntegrationCoefficient = integration_coefficients[GPoint];

Expand Down
70 changes: 47 additions & 23 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 @@ -327,18 +327,14 @@ void SmallStrainUPwDiffOrderElement::InitializeSolutionStep(const ProcessInfo& r

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) {
// compute element kinematics (Np, gradNpT, |J|, B, strains)
this->CalculateKinematics(Variables, GPoint);
Variables.B = b_matrices[GPoint];
Variables.F = deformation_gradients[GPoint];

// Compute infinitesimal strain
Variables.detF = determinants_of_deformation_gradients[GPoint];
Variables.B = b_matrices[GPoint];
Variables.F = deformation_gradients[GPoint];
Variables.StrainVector = strain_vectors[GPoint];

ConstitutiveLawUtilities::SetConstitutiveParameters(
ConstitutiveParameters, Variables.StrainVector, Variables.ConstitutiveMatrix,
Variables.Nu, Variables.DNu_DX, Variables.F, Variables.detF);
ConstitutiveParameters, Variables.StrainVector, Variables.ConstitutiveMatrix, Variables.Nu,
Variables.DNu_DX, Variables.F, determinants_of_deformation_gradients[GPoint]);

// compute constitutive tensor and/or stresses
noalias(Variables.StressVector) = mStressVector[GPoint];
Expand Down Expand Up @@ -550,12 +546,11 @@ void SmallStrainUPwDiffOrderElement::FinalizeSolutionStep(const ProcessInfo& rCu

// Compute infinitesimal strain
Variables.F = deformation_gradients[GPoint];
Variables.detF = determinants_of_deformation_gradients[GPoint];
Variables.StrainVector = strain_vectors[GPoint];

ConstitutiveLawUtilities::SetConstitutiveParameters(
ConstitutiveParameters, Variables.StrainVector, Variables.ConstitutiveMatrix,
Variables.Nu, Variables.DNu_DX, Variables.F, Variables.detF);
ConstitutiveParameters, Variables.StrainVector, Variables.ConstitutiveMatrix, Variables.Nu,
Variables.DNu_DX, Variables.F, determinants_of_deformation_gradients[GPoint]);

// compute constitutive tensor and/or stresses
noalias(Variables.StressVector) = mStressVector[GPoint];
Expand Down Expand Up @@ -1278,8 +1273,16 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi
this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters,
Variables.NuContainer, Variables.DNu_DXContainer,
strain_vectors, mStressVector, constitutive_matrices);

const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients(
constitutive_matrices, this->GetProperties());
const auto fluid_pressures = GeoTransportEquationUtilities::CalculateFluidPressures(
Variables.NpContainer, Variables.PressureVector);
const auto degrees_of_saturation = CalculateDegreesOfSaturation(fluid_pressures, RetentionParameters);
const auto derivatives_of_saturation =
CalculateDerivativesOfSaturation(fluid_pressures, RetentionParameters);
const auto biot_moduli_inverse = GeoTransportEquationUtilities::CalculateInverseBiotModuli(
biot_coefficients, degrees_of_saturation, derivatives_of_saturation, rProp);

for (unsigned int GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) {
this->CalculateKinematics(Variables, GPoint);
Expand All @@ -1291,11 +1294,9 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi
CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

Variables.BiotCoefficient = biot_coefficients[GPoint];
Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse(
Variables.BiotCoefficient, Variables.DegreeOfSaturation,
Variables.DerivativeOfSaturation, this->GetProperties());
Variables.BiotModulusInverse = biot_moduli_inverse[GPoint];
Variables.PermeabilityUpdateFactor = this->CalculatePermeabilityUpdateFactor(Variables.StrainVector);

Variables.DegreeOfSaturation = degrees_of_saturation[GPoint];
Variables.IntegrationCoefficient = integration_coefficients[GPoint];

Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient(
Expand All @@ -1312,6 +1313,34 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi
KRATOS_CATCH("")
}

std::vector<double> SmallStrainUPwDiffOrderElement::CalculateDerivativesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
Copy link
Contributor

Choose a reason for hiding this comment

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

To comply with the Kratos Style Guide, we need to rename the first parameter:

Suggested change
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
const std::vector<double>& rFluidPressures, RetentionLaw::Parameters& RetentionParameters)

And as mentioned before, I would make RetentionParameters a local variable, to limit its scope as much as possible.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done

{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
Copy link
Contributor

Choose a reason for hiding this comment

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

You might consider to raise an error when the number of fluid pressure values does not match the number of constitutive law pointers.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done

std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
RetentionParameters.SetFluidPressure(fluid_pressure);
return pRetentionLaw->CalculateDerivativeOfSaturation(RetentionParameters);
});

return result;
}

std::vector<double> SmallStrainUPwDiffOrderElement::CalculateDegreesOfSaturation(
const std::vector<double>& fluid_pressures, RetentionLaw::Parameters& RetentionParameters)
Copy link
Contributor

Choose a reason for hiding this comment

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

The above comments also apply to this member function.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done

{
std::vector<double> result;
std::transform(fluid_pressures.begin(), fluid_pressures.end(), mRetentionLawVector.begin(),
std::back_inserter(result),
[&RetentionParameters](double fluid_pressure, RetentionLaw::Pointer pRetentionLaw) {
RetentionParameters.SetFluidPressure(fluid_pressure);
return pRetentionLaw->CalculateSaturation(RetentionParameters);
});

return result;
}

void SmallStrainUPwDiffOrderElement::CalculateMaterialStiffnessMatrix(MatrixType& rStiffnessMatrix,
Copy link
Contributor

Choose a reason for hiding this comment

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

Duplication, as mentioned

const ProcessInfo& rCurrentProcessInfo)
{
Expand Down Expand Up @@ -1428,7 +1457,6 @@ void SmallStrainUPwDiffOrderElement::InitializeElementVariables(ElementVariables
rVariables.StressVector.resize(VoigtSize, false);

// Needed parameters for consistency with the general constitutive law
rVariables.detF = 1.0;
rVariables.F.resize(Dim, Dim, false);
noalias(rVariables.F) = identity_matrix<double>(Dim);

Expand All @@ -1443,10 +1471,9 @@ void SmallStrainUPwDiffOrderElement::InitializeElementVariables(ElementVariables
rVariables.DtPressureCoefficient = rCurrentProcessInfo[DT_PRESSURE_COEFFICIENT];

// Retention law
rVariables.DegreeOfSaturation = 1.0;
rVariables.DerivativeOfSaturation = 0.0;
rVariables.RelativePermeability = 1.0;
rVariables.BishopCoefficient = 1.0;
rVariables.DegreeOfSaturation = 1.0;
rVariables.RelativePermeability = 1.0;
rVariables.BishopCoefficient = 1.0;

// permeability change
rVariables.PermeabilityUpdateFactor = 1.0;
Expand Down Expand Up @@ -1987,9 +2014,6 @@ void SmallStrainUPwDiffOrderElement::CalculateRetentionResponse(ElementVariables
rRetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure(
rVariables.Np, rVariables.PressureVector));

rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters);
rVariables.DerivativeOfSaturation =
mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(rRetentionParameters);
rVariables.RelativePermeability =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(rRetentionParameters);
rVariables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(rRetentionParameters);
Expand Down
6 changes: 5 additions & 1 deletion ...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp
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.

I believe you may now also be able to remove data member DerivativeOfSaturation of nested class ElementVariables (just as you did for class UPwSmallStrainElement::ElementVariables).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done

Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub
Matrix ConstitutiveMatrix;

// Variables needed for consistency with the general constitutive law
double detF;
Matrix F;

// needed for updated Lagrangian:
Expand Down Expand Up @@ -314,6 +313,11 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub

Vector GetPressureSolutionVector();

std::vector<double> CalculateDegreesOfSaturation(const std::vector<double>& fluid_pressures,
RetentionLaw::Parameters& RetentionParameters);
std::vector<double> CalculateDerivativesOfSaturation(const std::vector<double>& fluid_pressures,
RetentionLaw::Parameters& RetentionParameters);
Copy link
Contributor

Choose a reason for hiding this comment

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

Is there any reason why these member functions can't be const? And perhaps you may want to add attribute [[nodiscard]] to them?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, the function calls in retention law that we call in these are not const. Since that's a member, I assumed we couldn't make this one const. Turns out these members are ::Pointers, meaning we can do it.

However, I think it's quite confusing to have a const function which calls a non-const function on an object that this element has ownership over. What do you think in this case?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

But indeed we can add [[nodiscard]] regardless


//----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

private:
Expand Down
Loading
Loading