From 5be8639cd1bdd8fbdab701aa516a33eb393033ef Mon Sep 17 00:00:00 2001 From: Gennady Markelov Date: Tue, 14 May 2024 15:16:40 +0200 Subject: [PATCH 1/4] extracted CalculateFluidPressure function --- .../U_Pw_small_strain_element.cpp | 16 ++--- .../U_Pw_small_strain_element.hpp | 3 +- .../U_Pw_small_strain_interface_element.cpp | 16 ++--- .../U_Pw_small_strain_interface_element.hpp | 4 +- .../small_strain_U_Pw_diff_order_element.cpp | 18 ++--- .../small_strain_U_Pw_diff_order_element.hpp | 2 - .../custom_retention/retention_law.h | 69 ++++++++----------- .../transport_equation_utilities.hpp | 3 +- .../tests/cpp_tests/test_fluid_pressure.cpp | 43 ++++++++++++ 9 files changed, 88 insertions(+), 86 deletions(-) create mode 100644 applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp 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 915720c57a96..958e19eb449e 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -1153,7 +1153,8 @@ std::vector> UPwSmallStrainElement::Calc GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - Variables.FluidPressure = CalculateFluidPressure(Variables); + Variables.FluidPressure = + RetentionLaw::Parameters::CalculateFluidPressure(Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(Variables.FluidPressure); Variables.RelativePermeability = @@ -1747,16 +1748,6 @@ void UPwSmallStrainElement::SetRetentionParameters(const Elemen KRATOS_CATCH("") } -template -double UPwSmallStrainElement::CalculateFluidPressure(const ElementVariables& rVariables) -{ - KRATOS_TRY - - return inner_prod(rVariables.Np, rVariables.PressureVector); - - KRATOS_CATCH("") -} - template void UPwSmallStrainElement::CalculateRetentionResponse(ElementVariables& rVariables, RetentionLaw::Parameters& rRetentionParameters, @@ -1764,7 +1755,8 @@ void UPwSmallStrainElement::CalculateRetentionResponse(ElementV { KRATOS_TRY - rVariables.FluidPressure = CalculateFluidPressure(rVariables); + rVariables.FluidPressure = + RetentionLaw::Parameters::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); SetRetentionParameters(rVariables, rRetentionParameters); rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); 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 79d39cbfeec2..46178bf6d27b 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -293,8 +293,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa void InitializeNodalPorePressureVariables(ElementVariables& rVariables); void InitializeNodalVolumeAccelerationVariables(ElementVariables& rVariables); - void InitializeProperties(ElementVariables& rVariables); - double CalculateFluidPressure(const ElementVariables& rVariables); + void InitializeProperties(ElementVariables& rVariables); std::vector> CalculateFluidFluxes(const std::vector& rPermeabilityUpdateFactors, const ProcessInfo& rCurrentProcessInfo); diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp index 218b6875c4a9..801f850f2f9b 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp @@ -600,7 +600,8 @@ void UPwSmallStrainInterfaceElement::CalculateOnIntegrationPoin RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - Variables.FluidPressure = CalculateFluidPressure(Variables); + Variables.FluidPressure = + RetentionLaw::Parameters::CalculateFluidPressure(Variables.Np, Variables.PressureVector); SetRetentionParameters(Variables, RetentionParameters); if (rVariable == DEGREE_OF_SATURATION) @@ -2103,23 +2104,14 @@ void UPwSmallStrainInterfaceElement::SetRetentionParameters( KRATOS_CATCH("") } -template -double UPwSmallStrainInterfaceElement::CalculateFluidPressure(const InterfaceElementVariables& rVariables) -{ - KRATOS_TRY - - return inner_prod(rVariables.Np, rVariables.PressureVector); - - KRATOS_CATCH("") -} - template void UPwSmallStrainInterfaceElement::CalculateRetentionResponse( InterfaceElementVariables& rVariables, RetentionLaw::Parameters& rRetentionParameters, unsigned int GPoint) { KRATOS_TRY - rVariables.FluidPressure = CalculateFluidPressure(rVariables); + rVariables.FluidPressure = + RetentionLaw::Parameters::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); SetRetentionParameters(rVariables, rRetentionParameters); rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.hpp index ea526388e7ba..d21ff6c4ba87 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.hpp @@ -284,8 +284,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainInterfaceElement void SetRetentionParameters(const InterfaceElementVariables& rVariables, RetentionLaw::Parameters& rRetentionParameters); - double CalculateFluidPressure(const InterfaceElementVariables& rVariables); - double CalculateBulkModulus(const Matrix& ConstitutiveMatrix); void InitializeBiotCoefficients(InterfaceElementVariables& rVariables, const bool& hasBiotCoefficient = false); @@ -295,7 +293,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainInterfaceElement unsigned int GPoint); void CalculateSoilGamma(InterfaceElementVariables& rVariables); - + void SetConstitutiveParameters(InterfaceElementVariables& rVariables, ConstitutiveLaw::Parameters& rConstitutiveParameters); 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 390a1d5b14c0..3d40252dcc65 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 @@ -917,7 +917,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - Variables.FluidPressure = CalculateFluidPressure(Variables); + Variables.FluidPressure = + RetentionLaw::Parameters::CalculateFluidPressure(Variables.Np, Variables.PressureVector); SetRetentionParameters(Variables, RetentionParameters); if (rVariable == DEGREE_OF_SATURATION) @@ -1022,7 +1023,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable BodyAcceleration[idim] += Variables.Nu[i] * Variables.BodyAcceleration[Index++]; } - CalculateFluidPressure(Variables); + Variables.FluidPressure = + RetentionLaw::Parameters::CalculateFluidPressure(Variables.Np, Variables.PressureVector); SetRetentionParameters(Variables, RetentionParameters); const double RelativePermeability = @@ -2113,15 +2115,6 @@ void SmallStrainUPwDiffOrderElement::CalculateJacobianOnCurrentConfiguration(dou KRATOS_CATCH("") } -double SmallStrainUPwDiffOrderElement::CalculateFluidPressure(const ElementVariables& rVariables) const -{ - KRATOS_TRY - - return inner_prod(rVariables.Np, rVariables.PressureVector); - - KRATOS_CATCH("") -} - void SmallStrainUPwDiffOrderElement::SetRetentionParameters(const ElementVariables& rVariables, RetentionLaw::Parameters& rRetentionParameters) const { @@ -2138,7 +2131,8 @@ void SmallStrainUPwDiffOrderElement::CalculateRetentionResponse(ElementVariables { KRATOS_TRY - rVariables.FluidPressure = CalculateFluidPressure(rVariables); + rVariables.FluidPressure = + RetentionLaw::Parameters::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); SetRetentionParameters(rVariables, rRetentionParameters); rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); 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 d048274567b8..8d47a766e929 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 @@ -307,8 +307,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub Matrix CalculateDeformationGradient(unsigned int GPoint) const; std::vector CalculateDeformationGradients() const; - double CalculateFluidPressure(const ElementVariables& rVariables) const; - void SetRetentionParameters(const ElementVariables& rVariables, RetentionLaw::Parameters& rRetentionParameters) const; diff --git a/applications/GeoMechanicsApplication/custom_retention/retention_law.h b/applications/GeoMechanicsApplication/custom_retention/retention_law.h index 8cad6da9cb43..2f8d3557b968 100644 --- a/applications/GeoMechanicsApplication/custom_retention/retention_law.h +++ b/applications/GeoMechanicsApplication/custom_retention/retention_law.h @@ -24,27 +24,30 @@ #include "includes/serializer.h" #include -namespace Kratos { +namespace Kratos +{ /** * Base class of retention laws. */ -class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw { +class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw +{ public: /** * Type definitions * NOTE: geometries are assumed to be of type Node for all problems */ using ProcessInfoType = ProcessInfo; - using SizeType = std::size_t; - using GeometryType = Geometry; + using SizeType = std::size_t; + using GeometryType = Geometry; /** * Counted pointer of RetentionLaw */ KRATOS_CLASS_POINTER_DEFINITION(RetentionLaw); - class Parameters { + class Parameters + { KRATOS_CLASS_POINTER_DEFINITION(Parameters); /** @@ -64,17 +67,12 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw { */ public: - Parameters(const Properties& rMaterialProperties, - const ProcessInfo& rCurrentProcessInfo) - : mrCurrentProcessInfo(rCurrentProcessInfo), - mrMaterialProperties(rMaterialProperties){}; + Parameters(const Properties& rMaterialProperties, const ProcessInfo& rCurrentProcessInfo) + : mrCurrentProcessInfo(rCurrentProcessInfo), mrMaterialProperties(rMaterialProperties){}; ~Parameters() = default; - void SetFluidPressure(double FluidPressure) - { - mFluidPressure = FluidPressure; - }; + void SetFluidPressure(double FluidPressure) { mFluidPressure = FluidPressure; }; double GetFluidPressure() const { @@ -84,20 +82,19 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw { return mFluidPressure.value(); } - const ProcessInfo& GetProcessInfo() const + [[nodiscard]] static double CalculateFluidPressure(const Vector& rN, const Vector& rPressureVector) { - return mrCurrentProcessInfo; + return inner_prod(rN, rPressureVector); } - const Properties& GetMaterialProperties() const - { - return mrMaterialProperties; - } + const ProcessInfo& GetProcessInfo() const { return mrCurrentProcessInfo; } + + const Properties& GetMaterialProperties() const { return mrMaterialProperties; } private: std::optional mFluidPressure; - const ProcessInfo& mrCurrentProcessInfo; - const Properties& mrMaterialProperties; + const ProcessInfo& mrCurrentProcessInfo; + const Properties& mrMaterialProperties; }; // class Parameters end @@ -121,9 +118,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw { * @param rValue a reference to the returned value * @param rValue output: the value of the specified variable */ - virtual double& CalculateValue(Parameters& rParameters, - const Variable& rThisVariable, - double& rValue) = 0; + virtual double& CalculateValue(Parameters& rParameters, const Variable& rThisVariable, double& rValue) = 0; virtual double CalculateSaturation(Parameters& rParameters) = 0; @@ -143,9 +138,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw { * @param rElementGeometry the geometry of the current element * @param rCurrentProcessInfo process info */ - virtual void InitializeMaterial(const Properties& rMaterialProperties, + virtual void InitializeMaterial(const Properties& rMaterialProperties, const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues); + const Vector& rShapeFunctionsValues); virtual void Initialize(Parameters& rParameters); @@ -175,9 +170,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw { * @param rShapeFunctionsValues the shape functions values in the current integration point * @param the current ProcessInfo instance */ - virtual void ResetMaterial(const Properties& rMaterialProperties, + virtual void ResetMaterial(const Properties& rMaterialProperties, const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues); + const Vector& rShapeFunctionsValues); /** * This function is designed to be called once to perform all the checks @@ -188,8 +183,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw { * @param rCurrentProcessInfo * @return */ - virtual int Check(const Properties& rMaterialProperties, - const ProcessInfo& rCurrentProcessInfo) = 0; + virtual int Check(const Properties& rMaterialProperties, const ProcessInfo& rCurrentProcessInfo) = 0; /** * @brief This method is used to check that two Retention Laws are the same type (references) @@ -212,22 +206,13 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw { } /// Turn back information as a string. - virtual std::string Info() const - { - return "RetentionLaw"; - } + virtual std::string Info() const { return "RetentionLaw"; } /// Print information about this object. - virtual void PrintInfo(std::ostream& rOStream) const - { - rOStream << Info(); - } + virtual void PrintInfo(std::ostream& rOStream) const { rOStream << Info(); } /// Print object's data. - virtual void PrintData(std::ostream& rOStream) const - { - rOStream << "RetentionLaw has no data"; - } + virtual void PrintData(std::ostream& rOStream) const { rOStream << "RetentionLaw has no data"; } private: friend class Serializer; diff --git a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp index cf85fefbdb35..589391a61c61 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp @@ -106,7 +106,8 @@ class GeoTransportEquationUtilities RetentionLaw::Parameters retention_parameters(rProp, rCurrentProcessInfo); Vector result(rRetentionLawVector.size()); for (unsigned int g_point = 0; g_point < rRetentionLawVector.size(); ++g_point) { - const double fluid_pressure = inner_prod(row(rNContainer, g_point), rPressureSolution); + const double fluid_pressure = RetentionLaw::Parameters::CalculateFluidPressure( + row(rNContainer, g_point), rPressureSolution); retention_parameters.SetFluidPressure(fluid_pressure); result(g_point) = rRetentionLawVector[g_point]->CalculateSaturation(retention_parameters); } diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp new file mode 100644 index 000000000000..02813eb84450 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp @@ -0,0 +1,43 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Gennady Markelov +// + +#include "custom_retention/retention_law.h" +#include "includes/checks.h" +#include "testing/testing.h" + +using namespace Kratos; + +namespace Kratos::Testing +{ + +KRATOS_TEST_CASE_IN_SUITE(CalculateFluidPressureGivesCorrectResults, KratosGeoMechanicsFastSuite) +{ + Vector N(5); + N(0) = 1.0; + N(1) = 2.0; + N(2) = 3.0; + N(3) = 4.0; + N(4) = 5.0; + + Vector pressure_vector(5); + pressure_vector(0) = 0.5; + pressure_vector(1) = 0.7; + pressure_vector(2) = 0.8; + pressure_vector(3) = 0.9; + pressure_vector(4) = 0.4; + + auto fluid_pressure = RetentionLaw::Parameters::CalculateFluidPressure(N, pressure_vector); + double expected_value = 9.9; + KRATOS_CHECK_NEAR(fluid_pressure, expected_value, 1e-12); +} + +} // namespace Kratos::Testing \ No newline at end of file From dd7040b0a8367354836d2f31fe1aecb27c783bbb Mon Sep 17 00:00:00 2001 From: Gennady Markelov Date: Wed, 15 May 2024 14:10:07 +0200 Subject: [PATCH 2/4] used ublas assigment --- .../tests/cpp_tests/test_fluid_pressure.cpp | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp index 02813eb84450..75e44eeb7787 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp @@ -10,6 +10,7 @@ // Main authors: Gennady Markelov // +#include "boost/numeric/ublas/assignment.hpp" #include "custom_retention/retention_law.h" #include "includes/checks.h" #include "testing/testing.h" @@ -22,18 +23,10 @@ namespace Kratos::Testing KRATOS_TEST_CASE_IN_SUITE(CalculateFluidPressureGivesCorrectResults, KratosGeoMechanicsFastSuite) { Vector N(5); - N(0) = 1.0; - N(1) = 2.0; - N(2) = 3.0; - N(3) = 4.0; - N(4) = 5.0; + N <<= 1.0, 2.0, 3.0, 4.0, 5.0; Vector pressure_vector(5); - pressure_vector(0) = 0.5; - pressure_vector(1) = 0.7; - pressure_vector(2) = 0.8; - pressure_vector(3) = 0.9; - pressure_vector(4) = 0.4; + pressure_vector <<= 0.5, 0.7, 0.8, 0.9, 0.4; auto fluid_pressure = RetentionLaw::Parameters::CalculateFluidPressure(N, pressure_vector); double expected_value = 9.9; From 1f6eea2e2b00408ea1b83c75664d0cb08676f13a Mon Sep 17 00:00:00 2001 From: Gennady Markelov Date: Wed, 15 May 2024 15:19:54 +0200 Subject: [PATCH 3/4] moved CalculateFluidPressure function to transport_equation_utilities --- .../custom_elements/U_Pw_small_strain_element.cpp | 4 ++-- .../U_Pw_small_strain_interface_element.cpp | 6 +++--- .../small_strain_U_Pw_diff_order_element.cpp | 10 +++++----- .../custom_retention/retention_law.h | 5 ----- .../custom_utilities/transport_equation_utilities.hpp | 8 ++++++-- 5 files changed, 16 insertions(+), 17 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 d2a72f23e708..0307ea59ef4a 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -1154,7 +1154,7 @@ std::vector> UPwSmallStrainElement::Calc GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); Variables.FluidPressure = - RetentionLaw::Parameters::CalculateFluidPressure(Variables.Np, Variables.PressureVector); + GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(Variables.FluidPressure); Variables.RelativePermeability = @@ -1745,7 +1745,7 @@ void UPwSmallStrainElement::CalculateRetentionResponse(ElementV KRATOS_TRY rVariables.FluidPressure = - RetentionLaw::Parameters::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); + GeoTransportEquationUtilities::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); rRetentionParameters.SetFluidPressure(rVariables.FluidPressure); rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp index f5035b2ef2cc..9c19d36e4888 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp @@ -600,8 +600,8 @@ void UPwSmallStrainInterfaceElement::CalculateOnIntegrationPoin RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - Variables.FluidPressure = - RetentionLaw::Parameters::CalculateFluidPressure(Variables.Np, Variables.PressureVector); + Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(Variables.FluidPressure); if (rVariable == DEGREE_OF_SATURATION) @@ -2100,7 +2100,7 @@ void UPwSmallStrainInterfaceElement::CalculateRetentionResponse KRATOS_TRY rVariables.FluidPressure = - RetentionLaw::Parameters::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); + GeoTransportEquationUtilities::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); rRetentionParameters.SetFluidPressure(rVariables.FluidPressure); rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); 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 973d604ffd6a..f58d46de765c 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 @@ -917,8 +917,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - Variables.FluidPressure = - RetentionLaw::Parameters::CalculateFluidPressure(Variables.Np, Variables.PressureVector); + Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(Variables.FluidPressure); if (rVariable == DEGREE_OF_SATURATION) @@ -1023,8 +1023,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable BodyAcceleration[idim] += Variables.Nu[i] * Variables.BodyAcceleration[Index++]; } - Variables.FluidPressure = - RetentionLaw::Parameters::CalculateFluidPressure(Variables.Np, Variables.PressureVector); + Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(Variables.FluidPressure); const double RelativePermeability = @@ -2122,7 +2122,7 @@ void SmallStrainUPwDiffOrderElement::CalculateRetentionResponse(ElementVariables KRATOS_TRY rVariables.FluidPressure = - RetentionLaw::Parameters::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); + GeoTransportEquationUtilities::CalculateFluidPressure(rVariables.Np, rVariables.PressureVector); rRetentionParameters.SetFluidPressure(rVariables.FluidPressure); rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); diff --git a/applications/GeoMechanicsApplication/custom_retention/retention_law.h b/applications/GeoMechanicsApplication/custom_retention/retention_law.h index 2f8d3557b968..d239aadcaa66 100644 --- a/applications/GeoMechanicsApplication/custom_retention/retention_law.h +++ b/applications/GeoMechanicsApplication/custom_retention/retention_law.h @@ -82,11 +82,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw return mFluidPressure.value(); } - [[nodiscard]] static double CalculateFluidPressure(const Vector& rN, const Vector& rPressureVector) - { - return inner_prod(rN, rPressureVector); - } - const ProcessInfo& GetProcessInfo() const { return mrCurrentProcessInfo; } const Properties& GetMaterialProperties() const { return mrMaterialProperties; } diff --git a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp index 589391a61c61..69f98888f409 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp @@ -106,13 +106,17 @@ class GeoTransportEquationUtilities RetentionLaw::Parameters retention_parameters(rProp, rCurrentProcessInfo); Vector result(rRetentionLawVector.size()); for (unsigned int g_point = 0; g_point < rRetentionLawVector.size(); ++g_point) { - const double fluid_pressure = RetentionLaw::Parameters::CalculateFluidPressure( - row(rNContainer, g_point), rPressureSolution); + const double fluid_pressure = CalculateFluidPressure(row(rNContainer, g_point), rPressureSolution); retention_parameters.SetFluidPressure(fluid_pressure); result(g_point) = rRetentionLawVector[g_point]->CalculateSaturation(retention_parameters); } return result; } + [[nodiscard]] static double CalculateFluidPressure(const Vector& rN, const Vector& rPressureVector) + { + return inner_prod(rN, rPressureVector); + } + }; /* Class GeoTransportEquationUtilities*/ } /* namespace Kratos.*/ From 1bbad677a36c362f61691c4ab2e9e060b2188d00 Mon Sep 17 00:00:00 2001 From: Gennady Markelov Date: Wed, 15 May 2024 15:36:56 +0200 Subject: [PATCH 4/4] fixed a call of CalculateFixedPressure --- .../tests/cpp_tests/test_fluid_pressure.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp index 75e44eeb7787..3aa8daa4c5cf 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp @@ -11,7 +11,7 @@ // #include "boost/numeric/ublas/assignment.hpp" -#include "custom_retention/retention_law.h" +#include "custom_utilities/transport_equation_utilities.hpp" #include "includes/checks.h" #include "testing/testing.h" @@ -28,7 +28,7 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateFluidPressureGivesCorrectResults, KratosGeoMe Vector pressure_vector(5); pressure_vector <<= 0.5, 0.7, 0.8, 0.9, 0.4; - auto fluid_pressure = RetentionLaw::Parameters::CalculateFluidPressure(N, pressure_vector); + auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure(N, pressure_vector); double expected_value = 9.9; KRATOS_CHECK_NEAR(fluid_pressure, expected_value, 1e-12); }