Skip to content

Commit a030ba3

Browse files
WPK4FEMWPK4FEM
authored
Geo/introduce geo drainage type (#14378)
* Introduce drainage type enum variable * Use Enum drainage type variable for the fully coupled cases. Made the material input variable a std::string. Cleaned some redundant things from material .json files. * IGNORE_UNDRAINED is converted to GEO_DRAINAGE_TYPEfor the constant Pw field computation * Replaced IGNORE_UNDRAINED true in material .json files. * Switched unit tests from IGNORE_UNDRAINED to GEO_DRAINAGE_TYPE input. * Adapted test README.md to GEO_DRAINAGE_TYPE input. * Renamed unit test to get rid of ignore undrained in names * Renamed IgnoreUndrained variables. * Sonar Cloud & clang tidy remarks in files I touched. * Use replacing utility for IGNORE_UNDRAINED. * Do not hide function in base class InitializeProperties now virtual in U_Pw_small_strain_element with an override in transient_Pw_element. * Conversion from IGNORE_UNDRAINED to GEO_DRAINAGE_TYPE from constructor to Initialize
1 parent f77e6f2 commit a030ba3

202 files changed

Lines changed: 693 additions & 913 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎applications/GeoMechanicsApplication/benchmarks/upw_diff_order_element_benchmark.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@ std::shared_ptr<Properties> CreatePropertiesForUPwDiffOrderElementBenchmark()
4242
p_properties->SetValue(BIOT_COEFFICIENT, 1.000000e+00);
4343
p_properties->SetValue(RETENTION_LAW, "SaturatedLaw");
4444
p_properties->SetValue(SATURATED_SATURATION, 1.000000e+00);
45-
p_properties->SetValue(IGNORE_UNDRAINED, false);
45+
p_properties->SetValue(GEO_DRAINAGE_TYPE, "FULLY_COUPLED");
4646

4747
return p_properties;
4848
}

‎applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
#include "custom_elements/U_Pw_base_element.h"
1616
#include "custom_retention/retention_law_factory.h"
1717
#include "custom_utilities/check_utilities.hpp"
18+
#include "custom_utilities/constitutive_law_utilities.h"
1819
#include "custom_utilities/dof_utilities.hpp"
1920
#include "custom_utilities/element_utilities.hpp"
2021
#include "custom_utilities/equation_of_motion_utilities.hpp"
@@ -94,6 +95,9 @@ void UPwBaseElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
9495
{
9596
KRATOS_TRY
9697

98+
// IGNORE_UNDRAINED is deprecated.
99+
ConstitutiveLawUtilities::ReplaceIgnoreUndrainedByDrainageType(GetProperties());
100+
97101
const auto& r_properties = this->GetProperties();
98102
const auto& r_geometry = this->GetGeometry();
99103
const auto number_of_integration_points = r_geometry.IntegrationPointsNumber(mThisIntegrationMethod);

‎applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.cpp‎

Lines changed: 15 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -88,11 +88,6 @@ Geo::ProcessInfoGetter CreateProcessInfoGetter(const ProcessInfo& rProcessInfo)
8888
return [&rProcessInfo]() -> const ProcessInfo& { return rProcessInfo; };
8989
}
9090

91-
bool GetIgnoreUndrained(const Properties& rProperties)
92-
{
93-
return rProperties.Has(IGNORE_UNDRAINED) ? rProperties[IGNORE_UNDRAINED] : false;
94-
}
95-
9691
} // namespace
9792

9893
namespace Kratos
@@ -158,9 +153,9 @@ void UPwInterfaceElement::EquationIdVector(EquationIdVectorType& rResult, const
158153

159154
void UPwInterfaceElement::CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, const ProcessInfo& rProcessInfo)
160155
{
161-
const auto number_of_dofs = GetDofs().size();
162-
const auto ignore_undrained = GetIgnoreUndrained(GetProperties());
163-
rLeftHandSideMatrix = ZeroMatrix{number_of_dofs, number_of_dofs};
156+
const auto number_of_dofs = GetDofs().size();
157+
const auto is_constant_pw_field = ConstitutiveLawUtilities::IsConstantWaterPressure(GetProperties());
158+
rLeftHandSideMatrix = ZeroMatrix{number_of_dofs, number_of_dofs};
164159

165160
for (auto contribution : mContributions) {
166161
switch (contribution) {
@@ -172,10 +167,10 @@ void UPwInterfaceElement::CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix,
172167
CalculateAndAssignUPCouplingMatrix(rLeftHandSideMatrix);
173168
break;
174169
case PUCoupling:
175-
if (!ignore_undrained) CalculateAndAssignPUCouplingMatrix(rLeftHandSideMatrix);
170+
if (!is_constant_pw_field) CalculateAndAssignPUCouplingMatrix(rLeftHandSideMatrix);
176171
break;
177172
case Permeability:
178-
if (!ignore_undrained) CalculateAndAssignPermeabilityMatrix(rLeftHandSideMatrix);
173+
if (!is_constant_pw_field) CalculateAndAssignPermeabilityMatrix(rLeftHandSideMatrix);
179174
break;
180175
case FluidBodyFlow:
181176
break;
@@ -267,8 +262,8 @@ void UPwInterfaceElement::CalculateAndAssignPermeabilityMatrix(Element::MatrixTy
267262
void UPwInterfaceElement::CalculateRightHandSide(Element::VectorType& rRightHandSideVector,
268263
const ProcessInfo& rProcessInfo)
269264
{
270-
const auto ignore_undrained = GetIgnoreUndrained(GetProperties());
271-
rRightHandSideVector = ZeroVector{GetDofs().size()};
265+
const auto is_constant_pw_field = ConstitutiveLawUtilities::IsConstantWaterPressure(GetProperties());
266+
rRightHandSideVector = ZeroVector{GetDofs().size()};
272267

273268
for (auto contribution : mContributions) {
274269
switch (contribution) {
@@ -280,13 +275,16 @@ void UPwInterfaceElement::CalculateRightHandSide(Element::VectorType& rRightHand
280275
CalculateAndAssembleUPCouplingForceVector(rRightHandSideVector);
281276
break;
282277
case PUCoupling:
283-
if (!ignore_undrained) CalculateAndAssemblePUCouplingForceVector(rRightHandSideVector);
278+
if (!is_constant_pw_field)
279+
CalculateAndAssemblePUCouplingForceVector(rRightHandSideVector);
284280
break;
285281
case Permeability:
286-
if (!ignore_undrained) CalculateAndAssemblePermeabilityFlowVector(rRightHandSideVector);
282+
if (!is_constant_pw_field)
283+
CalculateAndAssemblePermeabilityFlowVector(rRightHandSideVector);
287284
break;
288285
case FluidBodyFlow:
289-
if (!ignore_undrained) CalculateAndAssembleFluidBodyFlowVector(rRightHandSideVector);
286+
if (!is_constant_pw_field)
287+
CalculateAndAssembleFluidBodyFlowVector(rRightHandSideVector);
290288
break;
291289
default:
292290
KRATOS_ERROR << "This contribution is not supported \n";
@@ -464,6 +462,8 @@ void UPwInterfaceElement::GetDofList(DofsVectorType& rElementalDofList, const Pr
464462
void UPwInterfaceElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
465463
{
466464
Element::Initialize(rCurrentProcessInfo);
465+
// IGNORE_UNDRAINED is deprecated
466+
ConstitutiveLawUtilities::ReplaceIgnoreUndrainedByDrainageType(GetProperties());
467467

468468
mConstitutiveLaws.clear();
469469
mRetentionLaws.clear();

‎applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp‎

Lines changed: 6 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313

1414
// Application includes
1515
#include "custom_elements/U_Pw_small_strain_FIC_element.h"
16+
#include "custom_utilities/constitutive_law_utilities.h"
1617
#include "custom_utilities/equation_of_motion_utilities.hpp"
1718
#include "custom_utilities/extrapolation_utilities.h"
1819

@@ -121,13 +122,10 @@ int UPwSmallStrainFICElement<TDim, TNumNodes>::Check(const ProcessInfo& rCurrent
121122
int ierr = UPwSmallStrainElement<TDim, TNumNodes>::Check(rCurrentProcessInfo);
122123
if (ierr != 0) return ierr;
123124

124-
const PropertiesType& Prop = this->GetProperties();
125-
126-
// Verify specific properties
127-
if (Prop[IGNORE_UNDRAINED])
128-
KRATOS_ERROR << "IGNORE_UNDRAINED cannot be used in FIC elements. Use "
129-
"Non FIC elements instead"
130-
<< this->Id() << std::endl;
125+
KRATOS_ERROR_IF(ConstitutiveLawUtilities::IsConstantWaterPressure(this->GetProperties()))
126+
<< "Constant water pressure fields cannot be used in FIC elements. "
127+
"Use Non FIC elements instead"
128+
<< this->Id() << std::endl;
131129

132130
return ierr;
133131

@@ -522,7 +520,7 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
522520
template <unsigned int TDim, unsigned int TNumNodes>
523521
double UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateShearModulus(const Matrix& ConstitutiveMatrix) const
524522
{
525-
const int IndexG = ConstitutiveMatrix.size1() - 1;
523+
const auto IndexG = ConstitutiveMatrix.size1() - 1;
526524
return ConstitutiveMatrix(IndexG, IndexG);
527525
}
528526

‎applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp‎

Lines changed: 7 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@
2424
#include "custom_utilities/stress_strain_utilities.h"
2525
#include "custom_utilities/transport_equation_utilities.hpp"
2626
#include "custom_utilities/variables_utilities.hpp"
27+
#include "geo_mechanics_application_constants.h"
2728
#include "geo_mechanics_application_variables.h"
2829
#include "includes/cfd_variables.h"
2930

@@ -96,8 +97,7 @@ int UPwSmallStrainElement<TDim, TNumNodes>::Check(const ProcessInfo& rCurrentPro
9697

9798
const CheckProperties check_properties(r_properties, "property", this->Id(),
9899
CheckProperties::Bounds::AllInclusive);
99-
check_properties.CheckAvailability(IGNORE_UNDRAINED);
100-
if (!r_properties[IGNORE_UNDRAINED]) {
100+
if (!ConstitutiveLawUtilities::IsConstantWaterPressure(r_properties)) {
101101
check_properties.SingleUseBounds(CheckProperties::Bounds::AllExclusive).Check(BULK_MODULUS_FLUID);
102102
check_properties.SingleUseBounds(CheckProperties::Bounds::AllExclusive).Check(DYNAMIC_VISCOSITY);
103103
check_properties.CheckPermeabilityProperties(r_geometry.WorkingSpaceDimension());
@@ -765,8 +765,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLe
765765
ElementVariables Variables;
766766
this->InitializeElementVariables(Variables, rCurrentProcessInfo);
767767

768-
RetentionLaw::Parameters RetentionParameters(r_properties);
769-
770768
const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
771769
const auto integration_coefficients =
772770
this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer);
@@ -962,7 +960,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddLHS(MatrixType& rLef
962960

963961
this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables);
964962

965-
if (!rVariables.IgnoreUndrained) {
963+
if (!rVariables.IsConstantWaterPressure) {
966964
const auto permeability_matrix =
967965
GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
968966
rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix,
@@ -1000,7 +998,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddCouplingMatrix(Matri
1000998
rVariables.BiotCoefficient, rVariables.BishopCoefficient, rVariables.IntegrationCoefficient);
1001999
GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, coupling_matrix);
10021000

1003-
if (!rVariables.IgnoreUndrained) {
1001+
if (!rVariables.IsConstantWaterPressure) {
10041002
const auto p_coupling_matrix = GeoTransportEquationUtilities::CalculateCouplingMatrix<TDim, TNumNodes>(
10051003
rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
10061004
rVariables.BiotCoefficient, rVariables.DegreeOfSaturation, rVariables.IntegrationCoefficient);
@@ -1041,7 +1039,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddRHS(VectorType& rRig
10411039

10421040
this->CalculateAndAddCouplingTerms(rRightHandSideVector, rVariables);
10431041

1044-
if (!rVariables.IgnoreUndrained) {
1042+
if (!rVariables.IsConstantWaterPressure) {
10451043
this->CalculateAndAddCompressibilityFlow(rRightHandSideVector, rVariables);
10461044

10471045
this->CalculateAndAddPermeabilityFlow(rRightHandSideVector, rVariables);
@@ -1107,7 +1105,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddCouplingTerms(Vector
11071105
const array_1d<double, TNumNodes * TDim> coupling_force = prod(coupling_matrix, rVariables.PressureVector);
11081106
GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, (-1.0) * coupling_force);
11091107

1110-
if (!rVariables.IgnoreUndrained) {
1108+
if (!rVariables.IsConstantWaterPressure) {
11111109
const auto p_coupling_matrix = GeoTransportEquationUtilities::CalculateCouplingMatrix<TDim, TNumNodes>(
11121110
rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
11131111
rVariables.BiotCoefficient, rVariables.DegreeOfSaturation, rVariables.IntegrationCoefficient);
@@ -1311,7 +1309,7 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeProperties(ElementVariabl
13111309

13121310
const auto& r_properties = this->GetProperties();
13131311

1314-
rVariables.IgnoreUndrained = r_properties[IGNORE_UNDRAINED];
1312+
rVariables.IsConstantWaterPressure = ConstitutiveLawUtilities::IsConstantWaterPressure(r_properties);
13151313
rVariables.UseHenckyStrain = r_properties.Has(USE_HENCKY_STRAIN) ? r_properties[USE_HENCKY_STRAIN] : false;
13161314

13171315
rVariables.ConsiderGeometricStiffness =

‎applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -119,7 +119,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa
119119
protected:
120120
struct ElementVariables {
121121
/// Properties variables
122-
bool IgnoreUndrained;
122+
bool IsConstantWaterPressure;
123123
bool UseHenckyStrain;
124124
bool ConsiderGeometricStiffness;
125125
double DynamicViscosityInverse;
@@ -234,7 +234,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa
234234
void InitializeNodalPorePressureVariables(ElementVariables& rVariables);
235235
void InitializeNodalVolumeAccelerationVariables(ElementVariables& rVariables);
236236

237-
void InitializeProperties(ElementVariables& rVariables);
237+
virtual void InitializeProperties(ElementVariables& rVariables);
238238

239239
[[nodiscard]] std::vector<double> CalculateDegreesOfSaturation(const std::vector<double>& rFluidPressures) const;
240240
[[nodiscard]] std::vector<double> CalculateDerivativesOfSaturation(const std::vector<double>& rFluidPressures) const;

‎applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp‎

Lines changed: 10 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -111,8 +111,7 @@ int UPwSmallStrainInterfaceElement<TDim, TNumNodes>::Check(const ProcessInfo& rC
111111
const CheckProperties check_properties(r_properties, "property", element_Id,
112112
CheckProperties::Bounds::AllInclusive);
113113
check_properties.SingleUseBounds(CheckProperties::Bounds::AllExclusive).Check(MINIMUM_JOINT_WIDTH);
114-
check_properties.CheckAvailability(IGNORE_UNDRAINED);
115-
if (!r_properties[IGNORE_UNDRAINED]) {
114+
if (!ConstitutiveLawUtilities::IsConstantWaterPressure(r_properties)) {
116115
check_properties.Check(TRANSVERSAL_PERMEABILITY);
117116
check_properties.SingleUseBounds(CheckProperties::Bounds::AllExclusive).Check(BULK_MODULUS_FLUID);
118117
check_properties.SingleUseBounds(CheckProperties::Bounds::AllExclusive).Check(DYNAMIC_VISCOSITY);
@@ -1183,8 +1182,7 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::InitializeElementVariables
11831182
KRATOS_TRY
11841183

11851184
// Properties variables
1186-
rVariables.IgnoreUndrained = rProperties[IGNORE_UNDRAINED];
1187-
1185+
rVariables.IsConstantWaterPressure = ConstitutiveLawUtilities::IsConstantWaterPressure(rProperties);
11881186
rVariables.DynamicViscosityInverse = 1.0 / rProperties[DYNAMIC_VISCOSITY];
11891187

11901188
// ProcessInfo variables
@@ -1613,7 +1611,7 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAndAddLHS(MatrixT
16131611

16141612
this->CalculateAndAddStiffnessMatrix(rLeftHandSideMatrix, rVariables);
16151613

1616-
if (!rVariables.IgnoreUndrained) {
1614+
if (!rVariables.IsConstantWaterPressure) {
16171615
this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables);
16181616

16191617
this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix, rVariables);
@@ -1657,7 +1655,7 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAndAddCouplingMat
16571655

16581656
GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, coupling_matrix);
16591657

1660-
if (!rVariables.IgnoreUndrained) {
1658+
if (!rVariables.IsConstantWaterPressure) {
16611659
const double SaturationCoefficient = rVariables.DegreeOfSaturation / rVariables.BishopCoefficient;
16621660
const BoundedMatrix<double, TNumNodes, TNumNodes * TDim> transposed_coupling_matrix =
16631661
PORE_PRESSURE_SIGN_FACTOR * SaturationCoefficient * rVariables.VelocityCoefficient *
@@ -1712,7 +1710,7 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAndAddRHS(VectorT
17121710

17131711
this->CalculateAndAddCouplingTerms(rRightHandSideVector, rVariables);
17141712

1715-
if (!rVariables.IgnoreUndrained) {
1713+
if (!rVariables.IsConstantWaterPressure) {
17161714
this->CalculateAndAddCompressibilityFlow(rRightHandSideVector, rVariables);
17171715

17181716
this->CalculateAndAddPermeabilityFlow(rRightHandSideVector, rVariables);
@@ -1787,7 +1785,7 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAndAddCouplingTer
17871785

17881786
GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, u_vector);
17891787

1790-
if (!rVariables.IgnoreUndrained) {
1788+
if (!rVariables.IsConstantWaterPressure) {
17911789
const double SaturationCoefficient = rVariables.DegreeOfSaturation / rVariables.BishopCoefficient;
17921790
const array_1d<double, TNumNodes> coupling_flow_vector =
17931791
PORE_PRESSURE_SIGN_FACTOR * SaturationCoefficient *
@@ -1883,11 +1881,11 @@ double UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateBulkModulus(con
18831881
{
18841882
KRATOS_TRY
18851883

1886-
const int IndexM = ConstitutiveMatrix.size1() - 1;
1887-
const double M = ConstitutiveMatrix(IndexM, IndexM);
1888-
const double G = ConstitutiveMatrix(0, 0);
1884+
const auto IndexM = ConstitutiveMatrix.size1() - 1;
1885+
const auto normal_stiffnes = ConstitutiveMatrix(IndexM, IndexM);
1886+
const auto shear_stiffnes = ConstitutiveMatrix(0, 0);
18891887

1890-
return M - (4.0 / 3.0) * G;
1888+
return normal_stiffnes - (4.0 / 3.0) * shear_stiffnes;
18911889

18921890
KRATOS_CATCH("")
18931891
}

‎applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -108,7 +108,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainInterfaceElement : pub
108108

109109
struct InterfaceElementVariables {
110110
/// Properties variables
111-
bool IgnoreUndrained;
111+
bool IsConstantWaterPressure;
112112
double DynamicViscosityInverse;
113113
double BiotCoefficient;
114114
double BiotModulusInverse;

‎applications/GeoMechanicsApplication/custom_elements/contribution_calculators/compressibility_calculator.hpp‎

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414

1515
#include "contribution_calculator.h"
1616
#include "custom_retention/retention_law.h"
17+
#include "custom_utilities/constitutive_law_utilities.h"
1718
#include "custom_utilities/transport_equation_utilities.hpp"
1819
#include "geo_aliases.h"
1920
#include "geo_mechanics_application_variables.h"
@@ -103,15 +104,15 @@ class CompressibilityCalculator : public ContributionCalculator<TNumNodes>
103104
[[nodiscard]] double CalculateBiotModulusInverse(const RetentionLaw::Pointer& rRetentionLaw,
104105
double FluidPresssure) const
105106
{
106-
const auto& r_properties = mInputProvider.GetElementProperties();
107-
const double biot_coefficient = r_properties[BIOT_COEFFICIENT];
107+
const auto& r_properties = mInputProvider.GetElementProperties();
108+
const auto biot_coefficient = r_properties[BIOT_COEFFICIENT];
108109

109-
double bulk_fluid = TINY;
110-
if (!r_properties[IGNORE_UNDRAINED]) {
111-
bulk_fluid = r_properties[BULK_MODULUS_FLUID];
112-
}
113-
double result = (biot_coefficient - r_properties[POROSITY]) / r_properties[BULK_MODULUS_SOLID] +
114-
r_properties[POROSITY] / bulk_fluid;
110+
const auto bulk_fluid = ConstitutiveLawUtilities::IsConstantWaterPressure(r_properties)
111+
? TINY
112+
: r_properties[BULK_MODULUS_FLUID];
113+
114+
auto result = (biot_coefficient - r_properties[POROSITY]) / r_properties[BULK_MODULUS_SOLID] +
115+
r_properties[POROSITY] / bulk_fluid;
115116

116117
RetentionLaw::Parameters retention_parameters(r_properties);
117118
retention_parameters.SetFluidPressure(FluidPresssure);

0 commit comments

Comments
 (0)