Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[geoMechanicsApplication] extract CalculateFluidPressure function #12381

Merged
merged 5 commits into from
May 16, 2024
Merged

[geoMechanicsApplication] extract CalculateFluidPressure function #12381

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5be8639
extracted CalculateFluidPressure function
May 14, 2024
d04250f
Merge branch 'master' into geo/extract-calculate-fluid-pressure
May 15, 2024
dd7040b
used ublas assigment
May 15, 2024
1f6eea2
moved CalculateFluidPressure function to transport_equation_utilities
May 15, 2024
1bbad67
fixed a call of CalculateFixedPressure
May 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 4 additions & 12 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1153,7 +1153,8 @@ std::vector<array_1d<double, TDim>> UPwSmallStrainElement<TDim, TNumNodes>::Calc

GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
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 =
Expand Down Expand Up @@ -1736,24 +1737,15 @@ void UPwSmallStrainElement<TDim, TNumNodes>::SetConstitutiveParameters(ElementVa
KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
double UPwSmallStrainElement<TDim, TNumNodes>::CalculateFluidPressure(const ElementVariables& rVariables) const
{
KRATOS_TRY

return inner_prod(rVariables.Np, rVariables.PressureVector);

KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
void UPwSmallStrainElement<TDim, TNumNodes>::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);
Expand Down
3 changes: 1 addition & 2 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<array_1d<double, TDim>> CalculateFluidFluxes(const std::vector<double>& rPermeabilityUpdateFactors,
const ProcessInfo& rCurrentProcessInfo);

Expand Down
16 changes: 4 additions & 12 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -600,7 +600,8 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::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)
Expand Down Expand Up @@ -2092,23 +2093,14 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateAndAddFluidBodyFl
KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
double UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateFluidPressure(const InterfaceElementVariables& rVariables) const
{
KRATOS_TRY

return inner_prod(rVariables.Np, rVariables.PressureVector);

KRATOS_CATCH("")
}

template <unsigned int TDim, unsigned int TNumNodes>
void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::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);
Expand Down
4 changes: 1 addition & 3 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -281,8 +281,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainInterfaceElement
template <class TValueType>
void InterpolateOutputValues(std::vector<TValueType>& rOutput, const std::vector<TValueType>& GPValues);

[[nodiscard]] double CalculateFluidPressure(const InterfaceElementVariables& rVariables) const;

double CalculateBulkModulus(const Matrix& ConstitutiveMatrix);

void InitializeBiotCoefficients(InterfaceElementVariables& rVariables, const bool& hasBiotCoefficient = false);
Expand All @@ -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);

Expand Down
18 changes: 6 additions & 12 deletions ...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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 =
Expand Down Expand Up @@ -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);
Expand Down
1 change: 0 additions & 1 deletion ...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -307,7 +307,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub
Matrix CalculateDeformationGradient(unsigned int GPoint) const;
std::vector<Matrix> CalculateDeformationGradients() const;

[[nodiscard]] double CalculateFluidPressure(const ElementVariables& rVariables) const;

void CalculateRetentionResponse(ElementVariables& rVariables,
RetentionLaw::Parameters& rRetentionParameters,
Expand Down
68 changes: 24 additions & 44 deletions applications/GeoMechanicsApplication/custom_retention/retention_law.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,27 +24,30 @@
#include "includes/serializer.h"
#include <optional>

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<Node>;
using SizeType = std::size_t;
using GeometryType = Geometry<Node>;

/**
* Counted pointer of RetentionLaw
*/
KRATOS_CLASS_POINTER_DEFINITION(RetentionLaw);

class Parameters {
class Parameters
{
KRATOS_CLASS_POINTER_DEFINITION(Parameters);

/**
Expand All @@ -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
{
Expand All @@ -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<double> mFluidPressure;
const ProcessInfo& mrCurrentProcessInfo;
const Properties& mrMaterialProperties;
const ProcessInfo& mrCurrentProcessInfo;
const Properties& mrMaterialProperties;

}; // class Parameters end

Expand All @@ -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<double>& rThisVariable,
double& rValue) = 0;
virtual double& CalculateValue(Parameters& rParameters, const Variable<double>& rThisVariable, double& rValue) = 0;

virtual double CalculateSaturation(Parameters& rParameters) = 0;

Expand All @@ -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);

Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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;
Expand Down
7 changes: 6 additions & 1 deletion applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.*/
36 changes: 36 additions & 0 deletions applications/GeoMechanicsApplication/tests/cpp_tests/test_fluid_pressure.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
Copy link
Contributor

Choose a reason for hiding this comment

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

For an real interpolation function based on 5 nodes, these values are strange. However, for testing they are fine.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It is just for testing. Put something into and to get this later back.

Copy link
Contributor

Choose a reason for hiding this comment

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

No need to change.


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
Loading