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 726ffbe69c9c..0307ea59ef4a 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 = + GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(Variables.FluidPressure); Variables.RelativePermeability = @@ -1736,16 +1737,6 @@ void UPwSmallStrainElement::SetConstitutiveParameters(ElementVa KRATOS_CATCH("") } -template -double UPwSmallStrainElement::CalculateFluidPressure(const ElementVariables& rVariables) const -{ - KRATOS_TRY - - return inner_prod(rVariables.Np, rVariables.PressureVector); - - KRATOS_CATCH("") -} - template void UPwSmallStrainElement::CalculateRetentionResponse(ElementVariables& rVariables, RetentionLaw::Parameters& rRetentionParameters, @@ -1753,7 +1744,8 @@ void UPwSmallStrainElement::CalculateRetentionResponse(ElementV { KRATOS_TRY - rVariables.FluidPressure = CalculateFluidPressure(rVariables); + rVariables.FluidPressure = + 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_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp index 139f01fa0a42..98581f4aba03 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -291,8 +291,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa void InitializeNodalPorePressureVariables(ElementVariables& rVariables); void InitializeNodalVolumeAccelerationVariables(ElementVariables& rVariables); - void InitializeProperties(ElementVariables& rVariables); - [[nodiscard]] double CalculateFluidPressure(const ElementVariables& rVariables) const; + 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 dab2b6049f85..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,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 = GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(Variables.FluidPressure); if (rVariable == DEGREE_OF_SATURATION) @@ -2092,23 +2093,14 @@ void UPwSmallStrainInterfaceElement::CalculateAndAddFluidBodyFl KRATOS_CATCH("") } -template -double UPwSmallStrainInterfaceElement::CalculateFluidPressure(const InterfaceElementVariables& rVariables) const -{ - 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 = + 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.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.hpp index b1f22f6a27c9..e963a8e7f2ec 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 @@ -281,8 +281,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainInterfaceElement template void InterpolateOutputValues(std::vector& rOutput, const std::vector& GPValues); - [[nodiscard]] double CalculateFluidPressure(const InterfaceElementVariables& rVariables) const; - double CalculateBulkModulus(const Matrix& ConstitutiveMatrix); void InitializeBiotCoefficients(InterfaceElementVariables& rVariables, const bool& hasBiotCoefficient = false); @@ -292,7 +290,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 d6d321b44d52..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,7 +917,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); - Variables.FluidPressure = CalculateFluidPressure(Variables); + Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(Variables.FluidPressure); if (rVariable == DEGREE_OF_SATURATION) @@ -1022,7 +1023,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable BodyAcceleration[idim] += Variables.Nu[i] * Variables.BodyAcceleration[Index++]; } - Variables.FluidPressure = CalculateFluidPressure(Variables); + Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(Variables.FluidPressure); const double RelativePermeability = @@ -2113,22 +2115,14 @@ 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::CalculateRetentionResponse(ElementVariables& rVariables, RetentionLaw::Parameters& rRetentionParameters, unsigned int GPoint) { KRATOS_TRY - rVariables.FluidPressure = CalculateFluidPressure(rVariables); + rVariables.FluidPressure = + 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.hpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp index dc6490f901cd..377e0739075a 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,7 +307,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub Matrix CalculateDeformationGradient(unsigned int GPoint) const; std::vector CalculateDeformationGradients() const; - [[nodiscard]] double CalculateFluidPressure(const ElementVariables& rVariables) const; void CalculateRetentionResponse(ElementVariables& rVariables, RetentionLaw::Parameters& rRetentionParameters, diff --git a/applications/GeoMechanicsApplication/custom_retention/retention_law.h b/applications/GeoMechanicsApplication/custom_retention/retention_law.h index 8cad6da9cb43..d239aadcaa66 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,14 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw { return mFluidPressure.value(); } - const ProcessInfo& GetProcessInfo() const - { - return mrCurrentProcessInfo; - } + const ProcessInfo& GetProcessInfo() const { return mrCurrentProcessInfo; } - const Properties& GetMaterialProperties() const - { - return mrMaterialProperties; - } + 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 +113,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 +133,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 +165,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 +178,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 +201,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..69f98888f409 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp @@ -106,12 +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 = inner_prod(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.*/ 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..3aa8daa4c5cf --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp @@ -0,0 +1,36 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Gennady Markelov +// + +#include "boost/numeric/ublas/assignment.hpp" +#include "custom_utilities/transport_equation_utilities.hpp" +#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 <<= 1.0, 2.0, 3.0, 4.0, 5.0; + + Vector pressure_vector(5); + pressure_vector <<= 0.5, 0.7, 0.8, 0.9, 0.4; + + auto fluid_pressure = GeoTransportEquationUtilities::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