From 7ddb7b5fd3daaf2008cd31598f7287f540e4171c Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Fri, 17 May 2024 15:17:42 +0200 Subject: [PATCH 01/10] Replaced two calls to `CalculateRetentionResponse` by more specific ones `CalculateRetentionResponse` calculates several parameters, even when not all of them are needed. In these two instances, it was clear how to replace these function calls, which makes the code easier to understand. --- .../U_Pw_small_strain_element.cpp | 26 +++++++++++++------ 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp index 1c16f05bfb9a..dc3d0f116b4b 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -520,15 +520,22 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); - - if (rVariable == DEGREE_OF_SATURATION) rOutput[GPoint] = Variables.DegreeOfSaturation; - if (rVariable == EFFECTIVE_SATURATION) rOutput[GPoint] = Variables.EffectiveSaturation; - if (rVariable == BISHOP_COEFFICIENT) rOutput[GPoint] = Variables.BishopCoefficient; + Variables.FluidPressure = + GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); + RetentionParameters.SetFluidPressure(Variables.FluidPressure); + + if (rVariable == DEGREE_OF_SATURATION) + rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateSaturation(RetentionParameters); + if (rVariable == EFFECTIVE_SATURATION) + rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(RetentionParameters); + if (rVariable == BISHOP_COEFFICIENT) + rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); if (rVariable == DERIVATIVE_OF_SATURATION) - rOutput[GPoint] = Variables.DerivativeOfSaturation; + rOutput[GPoint] = + mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(RetentionParameters); if (rVariable == RELATIVE_PERMEABILITY) - rOutput[GPoint] = Variables.RelativePermeability; + rOutput[GPoint] = + mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); } } else if (rVariable == HYDRAULIC_HEAD) { const PropertiesType& rProp = this->GetProperties(); @@ -714,7 +721,10 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const Variables.BiotCoefficient = CalculateBiotCoefficient(Variables, hasBiotCoefficient); - this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + Variables.FluidPressure = + GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); + RetentionParameters.SetFluidPressure(Variables.FluidPressure); + Variables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); noalias(TotalStressVector) = mStressVector[GPoint]; noalias(TotalStressVector) += PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * From 565f47eefeede3ac24c8ac1d51d69fe68d0971f5 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Fri, 17 May 2024 15:27:37 +0200 Subject: [PATCH 02/10] Removed two redundant `using` statements --- .../custom_elements/steady_state_Pw_element.hpp | 1 - .../custom_elements/transient_Pw_element.hpp | 1 - 2 files changed, 2 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.hpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.hpp index 4bdfab8a5172..d05eaa60b071 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.hpp @@ -45,7 +45,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SteadyStatePwElement : public Transi /// The definition of the sizetype using SizeType = std::size_t; - using BaseType::CalculateRetentionResponse; using BaseType::mConstitutiveLawVector; using BaseType::mIsInitialised; using BaseType::mRetentionLawVector; diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp index 0297b0917c98..bd2c46485e79 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp @@ -45,7 +45,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwElement : public UPwSmall /// The definition of the sizetype using SizeType = std::size_t; - using BaseType::CalculateRetentionResponse; using BaseType::mConstitutiveLawVector; using BaseType::mIsInitialised; using BaseType::mRetentionLawVector; From fdf5ca25c6afb7f7ca9b83fe8783fa10080e3921 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Fri, 17 May 2024 15:33:17 +0200 Subject: [PATCH 03/10] Use a local variable for the Bishop coefficient when possible This eliminates two usages of member `BishopCoefficient` of struct `ElementVariables`. --- .../custom_elements/U_Pw_small_strain_element.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp index dc3d0f116b4b..2af7a8f42cd7 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -724,11 +724,11 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(Variables.FluidPressure); - Variables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); + const auto bishop_coefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); noalias(TotalStressVector) = mStressVector[GPoint]; noalias(TotalStressVector) += PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * - Variables.BishopCoefficient * Variables.FluidPressure * VoigtVector; + bishop_coefficient * Variables.FluidPressure * VoigtVector; if (rOutput[GPoint].size() != TotalStressVector.size()) rOutput[GPoint].resize(TotalStressVector.size(), false); From 2eb434e1130130528576d3a868e553008b951705 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Fri, 17 May 2024 15:43:43 +0200 Subject: [PATCH 04/10] Eliminated member `FluidPressure` of struct `ElementVariables` Its usages could be replaced by using local variables. --- .../U_Pw_small_strain_element.cpp | 19 +++++++++---------- .../U_Pw_small_strain_element.hpp | 1 - .../custom_elements/transient_Pw_element.cpp | 1 - 3 files changed, 9 insertions(+), 12 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp index 2af7a8f42cd7..4b12692e52ed 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -520,9 +520,9 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - Variables.FluidPressure = + const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); - RetentionParameters.SetFluidPressure(Variables.FluidPressure); + RetentionParameters.SetFluidPressure(fluid_pressure); if (rVariable == DEGREE_OF_SATURATION) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateSaturation(RetentionParameters); @@ -721,14 +721,14 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const Variables.BiotCoefficient = CalculateBiotCoefficient(Variables, hasBiotCoefficient); - Variables.FluidPressure = + const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); - RetentionParameters.SetFluidPressure(Variables.FluidPressure); + RetentionParameters.SetFluidPressure(fluid_pressure); const auto bishop_coefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); noalias(TotalStressVector) = mStressVector[GPoint]; noalias(TotalStressVector) += PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * - bishop_coefficient * Variables.FluidPressure * VoigtVector; + bishop_coefficient * fluid_pressure * VoigtVector; if (rOutput[GPoint].size() != TotalStressVector.size()) rOutput[GPoint].resize(TotalStressVector.size(), false); @@ -1098,9 +1098,9 @@ std::vector> UPwSmallStrainElement::Calc GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - Variables.FluidPressure = + const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); - RetentionParameters.SetFluidPressure(Variables.FluidPressure); + RetentionParameters.SetFluidPressure(fluid_pressure); Variables.RelativePermeability = mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); @@ -1199,7 +1199,6 @@ void UPwSmallStrainElement::InitializeElementVariables(ElementV rVariables.UVoigtMatrix.resize(TNumNodes * TDim, VoigtSize, false); // Retention law - rVariables.FluidPressure = 0.0; rVariables.DegreeOfSaturation = 1.0; rVariables.DerivativeOfSaturation = 0.0; rVariables.RelativePermeability = 1.0; @@ -1638,9 +1637,9 @@ void UPwSmallStrainElement::CalculateRetentionResponse(ElementV { KRATOS_TRY - rVariables.FluidPressure = + const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); - rRetentionParameters.SetFluidPressure(rVariables.FluidPressure); + rRetentionParameters.SetFluidPressure(fluid_pressure); rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); rVariables.DerivativeOfSaturation = diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp index 11674f49bed8..c0eccffbb2f4 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -179,7 +179,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa GeometryType::ShapeFunctionsGradientsType DN_DXContainer; /// Retention Law parameters - double FluidPressure; double DegreeOfSaturation; double DerivativeOfSaturation; double RelativePermeability; diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp index 1a009644c290..1c4ca70a288f 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp @@ -522,7 +522,6 @@ void TransientPwElement::InitializeElementVariables(ElementVari rVariables.DN_DXContainer, rVariables.detJContainer, this->GetIntegrationMethod()); // Retention law - rVariables.FluidPressure = 0.0; rVariables.DegreeOfSaturation = 1.0; rVariables.DerivativeOfSaturation = 0.0; rVariables.RelativePermeability = 1.0; From 63a6b30fc763d2d4ce5f816d92585a65b137f3b1 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Fri, 17 May 2024 15:48:42 +0200 Subject: [PATCH 05/10] Removed member `EffectiveSaturation` of struct `ElementVariables` It was assigned to only once, but never read from (due to earlier changes). --- .../custom_elements/U_Pw_small_strain_element.cpp | 2 -- .../custom_elements/U_Pw_small_strain_element.hpp | 1 - 2 files changed, 3 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp index 4b12692e52ed..b0847e1719c7 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -1647,8 +1647,6 @@ void UPwSmallStrainElement::CalculateRetentionResponse(ElementV rVariables.RelativePermeability = mRetentionLawVector[GPoint]->CalculateRelativePermeability(rRetentionParameters); rVariables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(rRetentionParameters); - rVariables.EffectiveSaturation = - mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(rRetentionParameters); KRATOS_CATCH("") } diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp index c0eccffbb2f4..e9397115b85d 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -183,7 +183,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa double DerivativeOfSaturation; double RelativePermeability; double BishopCoefficient; - double EffectiveSaturation; // needed for updated Lagrangian: double detJ; From f650827e1917ccb7e6bda6a25e70297b49dfee6d Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Fri, 17 May 2024 16:00:09 +0200 Subject: [PATCH 06/10] Replaced a call to `CalculateRetentionResponse` by a more specific one In this case, only the fluid pressure and the Bishop coefficient were needed. --- .../custom_elements/small_strain_U_Pw_diff_order_element.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp index 73e4147146c7..8e6ac847dda7 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp @@ -1164,7 +1164,10 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable Variables.ConstitutiveMatrix = constitutive_matrices[GPoint]; Variables.BiotCoefficient = CalculateBiotCoefficient(Variables, hasBiotCoefficient); - this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + Variables.FluidPressure = + GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); + RetentionParameters.SetFluidPressure(Variables.FluidPressure); + Variables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); noalias(TotalStressVector) = mStressVector[GPoint]; noalias(TotalStressVector) += PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * From 17afc1b3b1e01246a7b6a18ada1cf5b119edb615 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Fri, 17 May 2024 16:13:45 +0200 Subject: [PATCH 07/10] Removed member `FluidPressure` from struct `ElementVariables` Replaced all of its usages by local variables. Also replaced two usages of member `BishopCoefficient` by a local variable. --- .../small_strain_U_Pw_diff_order_element.cpp | 21 +++++++++---------- .../small_strain_U_Pw_diff_order_element.hpp | 1 - 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp index 8e6ac847dda7..eaf9c2dff151 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp @@ -923,9 +923,9 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure( + const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure( Variables.Np, Variables.PressureVector); - RetentionParameters.SetFluidPressure(Variables.FluidPressure); + RetentionParameters.SetFluidPressure(fluid_pressure); if (rVariable == DEGREE_OF_SATURATION) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateSaturation(RetentionParameters); @@ -1030,9 +1030,9 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable BodyAcceleration[idim] += Variables.Nu[i] * Variables.BodyAcceleration[Index++]; } - Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure( + const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure( Variables.Np, Variables.PressureVector); - RetentionParameters.SetFluidPressure(Variables.FluidPressure); + RetentionParameters.SetFluidPressure(fluid_pressure); const double RelativePermeability = mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); @@ -1164,14 +1164,14 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable Variables.ConstitutiveMatrix = constitutive_matrices[GPoint]; Variables.BiotCoefficient = CalculateBiotCoefficient(Variables, hasBiotCoefficient); - Variables.FluidPressure = + const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); - RetentionParameters.SetFluidPressure(Variables.FluidPressure); - Variables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); + RetentionParameters.SetFluidPressure(fluid_pressure); + const auto bishop_coefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); noalias(TotalStressVector) = mStressVector[GPoint]; noalias(TotalStressVector) += PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * - Variables.BishopCoefficient * Variables.FluidPressure * VoigtVector; + bishop_coefficient * fluid_pressure * VoigtVector; if (rOutput[GPoint].size() != TotalStressVector.size()) { rOutput[GPoint].resize(TotalStressVector.size(), false); @@ -1470,7 +1470,6 @@ void SmallStrainUPwDiffOrderElement::InitializeElementVariables(ElementVariables rVariables.DtPressureCoefficient = rCurrentProcessInfo[DT_PRESSURE_COEFFICIENT]; // Retention law - rVariables.FluidPressure = 0.0; rVariables.DegreeOfSaturation = 1.0; rVariables.DerivativeOfSaturation = 0.0; rVariables.RelativePermeability = 1.0; @@ -2055,9 +2054,9 @@ void SmallStrainUPwDiffOrderElement::CalculateRetentionResponse(ElementVariables { KRATOS_TRY - rVariables.FluidPressure = + const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); - rRetentionParameters.SetFluidPressure(rVariables.FluidPressure); + rRetentionParameters.SetFluidPressure(fluid_pressure); rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); rVariables.DerivativeOfSaturation = diff --git a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp index 4893b9ada2d4..8424d282589e 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp @@ -193,7 +193,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub Vector PressureDtVector; /// Retention Law parameters - double FluidPressure; double DegreeOfSaturation; double DerivativeOfSaturation; double RelativePermeability; From 191efc5eea6b662a1b017dddcebcf2e852ea1436 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Fri, 17 May 2024 17:28:09 +0200 Subject: [PATCH 08/10] Removed member `Density` from struct `ElementVariables` Its single usage has been replaced by a local variable. --- .../custom_elements/small_strain_U_Pw_diff_order_element.cpp | 4 ++-- .../custom_elements/small_strain_U_Pw_diff_order_element.hpp | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp index eaf9c2dff151..c015bc17ac4b 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp @@ -1832,7 +1832,7 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddMixBodyForce(VectorType& rRi const SizeType Dim = rGeom.WorkingSpaceDimension(); const SizeType NumUNodes = rGeom.PointsNumber(); - rVariables.Density = GeoTransportEquationUtilities::CalculateSoilDensity( + const auto soil_density = GeoTransportEquationUtilities::CalculateSoilDensity( rVariables.DegreeOfSaturation, this->GetProperties()); Vector BodyAcceleration = ZeroVector(Dim); @@ -1847,7 +1847,7 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddMixBodyForce(VectorType& rRi Index = i * Dim; for (SizeType idim = 0; idim < Dim; ++idim) { rRightHandSideVector[Index + idim] += - rVariables.Nu[i] * rVariables.Density * BodyAcceleration[idim] * + rVariables.Nu[i] * soil_density * BodyAcceleration[idim] * rVariables.IntegrationCoefficientInitialConfiguration; } } diff --git a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp index 8424d282589e..c5f22acbf3ae 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp @@ -197,7 +197,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub double DerivativeOfSaturation; double RelativePermeability; double BishopCoefficient; - double Density; // Properties and processinfo variables bool IgnoreUndrained; From 9cf87eb578e65b17dd3a334d400b7a4fa1d45be2 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Tue, 21 May 2024 14:05:49 +0200 Subject: [PATCH 09/10] Fixed GCC builds by adding `this->` to `CalculateRetentionResponse` --- .../custom_elements/steady_state_Pw_element.cpp | 2 +- .../custom_elements/transient_Pw_element.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.cpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.cpp index f406c08bc01e..ef3f96f25813 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.cpp @@ -168,7 +168,7 @@ void SteadyStatePwElement::CalculateAll(MatrixType& rLef GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp index 1c4ca70a288f..250af15e7405 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp @@ -468,7 +468,7 @@ void TransientPwElement::CalculateAll(MatrixType& rLeftH GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); From 85da166b2fe2c20c886cec99e1bcdc871bd827a0 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Wed, 22 May 2024 09:26:35 +0200 Subject: [PATCH 10/10] Processed Wijtze Pieter's review comments - Made some selections mutually exclusive. - Eliminated a few intermediate variables. - Made the calculation of the total stress vector more straightforward. - Removed a few unnecessary resize operations. --- .../U_Pw_small_strain_element.cpp | 36 +++++++------------ .../small_strain_U_Pw_diff_order_element.cpp | 35 +++++++----------- 2 files changed, 25 insertions(+), 46 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp index b0847e1719c7..edcd25595a76 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -520,20 +520,19 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - const auto fluid_pressure = - GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); - RetentionParameters.SetFluidPressure(fluid_pressure); + RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector)); if (rVariable == DEGREE_OF_SATURATION) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateSaturation(RetentionParameters); - if (rVariable == EFFECTIVE_SATURATION) + else if (rVariable == EFFECTIVE_SATURATION) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(RetentionParameters); - if (rVariable == BISHOP_COEFFICIENT) + else if (rVariable == BISHOP_COEFFICIENT) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); - if (rVariable == DERIVATIVE_OF_SATURATION) + else if (rVariable == DERIVATIVE_OF_SATURATION) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(RetentionParameters); - if (rVariable == RELATIVE_PERMEABILITY) + else if (rVariable == RELATIVE_PERMEABILITY) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); } @@ -700,8 +699,6 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT); - Vector TotalStressVector(mStressVector[0].size()); - const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); auto strain_vectors = StressStrainUtilities::CalculateStrains( @@ -726,14 +723,8 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const RetentionParameters.SetFluidPressure(fluid_pressure); const auto bishop_coefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); - noalias(TotalStressVector) = mStressVector[GPoint]; - noalias(TotalStressVector) += PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * - bishop_coefficient * fluid_pressure * VoigtVector; - - if (rOutput[GPoint].size() != TotalStressVector.size()) - rOutput[GPoint].resize(TotalStressVector.size(), false); - - rOutput[GPoint] = TotalStressVector; + rOutput[GPoint] = mStressVector[GPoint] + PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * + bishop_coefficient * fluid_pressure * VoigtVector; } } else if (rVariable == ENGINEERING_STRAIN_VECTOR) { // Definition of variables @@ -1098,9 +1089,9 @@ std::vector> UPwSmallStrainElement::Calc GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - const auto fluid_pressure = - GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); - RetentionParameters.SetFluidPressure(fluid_pressure); + + RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector)); Variables.RelativePermeability = mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); @@ -1637,9 +1628,8 @@ void UPwSmallStrainElement::CalculateRetentionResponse(ElementV { KRATOS_TRY - const auto fluid_pressure = - GeoTransportEquationUtilities::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); - rRetentionParameters.SetFluidPressure(fluid_pressure); + rRetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure( + rVariables.Np, rVariables.PressureVector)); rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); rVariables.DerivativeOfSaturation = diff --git a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp index c015bc17ac4b..4e1447819a95 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp @@ -923,20 +923,19 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure( - Variables.Np, Variables.PressureVector); - RetentionParameters.SetFluidPressure(fluid_pressure); + RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector)); if (rVariable == DEGREE_OF_SATURATION) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateSaturation(RetentionParameters); - if (rVariable == EFFECTIVE_SATURATION) + else if (rVariable == EFFECTIVE_SATURATION) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(RetentionParameters); - if (rVariable == BISHOP_COEFFICIENT) + else if (rVariable == BISHOP_COEFFICIENT) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); - if (rVariable == DERIVATIVE_OF_SATURATION) + else if (rVariable == DERIVATIVE_OF_SATURATION) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(RetentionParameters); - if (rVariable == RELATIVE_PERMEABILITY) + else if (rVariable == RELATIVE_PERMEABILITY) rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); } @@ -1030,9 +1029,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable BodyAcceleration[idim] += Variables.Nu[i] * Variables.BodyAcceleration[Index++]; } - const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure( - Variables.Np, Variables.PressureVector); - RetentionParameters.SetFluidPressure(fluid_pressure); + RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector)); const double RelativePermeability = mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); @@ -1143,8 +1141,6 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT); - Vector TotalStressVector(mStressVector[0].size()); - const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); auto strain_vectors = StressStrainUtilities::CalculateStrains( @@ -1169,14 +1165,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable RetentionParameters.SetFluidPressure(fluid_pressure); const auto bishop_coefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); - noalias(TotalStressVector) = mStressVector[GPoint]; - noalias(TotalStressVector) += PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * - bishop_coefficient * fluid_pressure * VoigtVector; - - if (rOutput[GPoint].size() != TotalStressVector.size()) { - rOutput[GPoint].resize(TotalStressVector.size(), false); - } - rOutput[GPoint] = TotalStressVector; + rOutput[GPoint] = mStressVector[GPoint] + PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * + bishop_coefficient * fluid_pressure * VoigtVector; } } else { for (unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i) @@ -2054,9 +2044,8 @@ void SmallStrainUPwDiffOrderElement::CalculateRetentionResponse(ElementVariables { KRATOS_TRY - const auto fluid_pressure = - GeoTransportEquationUtilities::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); - rRetentionParameters.SetFluidPressure(fluid_pressure); + rRetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure( + rVariables.Np, rVariables.PressureVector)); rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); rVariables.DerivativeOfSaturation =