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] Unit tests for TransientPwElement #13012

Merged
merged 24 commits into from
Jan 20, 2025
Merged

[GeoMechanicsApplication] Unit tests for TransientPwElement #13012

Show file tree
Hide file tree
Changes from 22 commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
bc860ea
added unit test
markelov208 Jan 16, 2025
63d6c30
added forgotten file
markelov208 Jan 16, 2025
2be611a
removed lines added for debugging
markelov208 Jan 16, 2025
e5d6494
kp format
markelov208 Jan 16, 2025
989ebb6
kp format and Z permeabilities
markelov208 Jan 16, 2025
e89e3b9
removed commented lines
markelov208 Jan 16, 2025
7987f55
applied Kratos style for variables
markelov208 Jan 16, 2025
464163f
cleaning and typeid
markelov208 Jan 16, 2025
2a91f96
replaced typeid with dynamic_cast
markelov208 Jan 16, 2025
98ea7d7
created SetBasicPropertiesAndVariables
markelov208 Jan 16, 2025
21b0e2c
typo in RemoveThreeNOdes
markelov208 Jan 16, 2025
5cf570a
renamed GPoint
markelov208 Jan 16, 2025
25a62a9
used ThreeDimensionalStressState for 3D element
markelov208 Jan 16, 2025
af842de
used for loop
markelov208 Jan 16, 2025
be4cd89
used for each
markelov208 Jan 16, 2025
da14762
used template functions
markelov208 Jan 16, 2025
8a4fe5c
clarified name CreateCoincidentNodes
markelov208 Jan 16, 2025
619b882
minor changes
markelov208 Jan 16, 2025
7b5c979
changed node to r_node and enable a call of RetentionLaw Check
markelov208 Jan 17, 2025
55098bc
add nit test for zero return and empty functions
markelov208 Jan 17, 2025
b53a76b
used all_of and none_of
markelov208 Jan 17, 2025
73a2139
replaced mRetentionLawVector.size() with number_of_integration_points
markelov208 Jan 17, 2025
f836652
removed ; to avoid code smells
markelov208 Jan 17, 2025
2990678
removed unit tests for empty functions
markelov208 Jan 20, 2025
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
200 changes: 100 additions & 100 deletions applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,11 +43,11 @@ void TransientPwElement<TDim, TNumNodes>::CalculateMassMatrix(MatrixType& rMassM
{
KRATOS_TRY

const unsigned int N_DOF = this->GetNumberOfDOF();
const unsigned int n_DoF = this->GetNumberOfDOF();

// Resizing mass matrix
if (rMassMatrix.size1() != N_DOF) rMassMatrix.resize(N_DOF, N_DOF, false);
noalias(rMassMatrix) = ZeroMatrix(N_DOF, N_DOF);
if (rMassMatrix.size1() != n_DoF) rMassMatrix.resize(n_DoF, n_DoF, false);
noalias(rMassMatrix) = ZeroMatrix(n_DoF, n_DoF);

KRATOS_CATCH("")
}
Expand All @@ -58,11 +58,11 @@ void TransientPwElement<TDim, TNumNodes>::CalculateDampingMatrix(MatrixType& rDa
{
KRATOS_TRY

const unsigned int N_DOF = this->GetNumberOfDOF();
const unsigned int n_DoF = this->GetNumberOfDOF();

// Compute Damping Matrix
if (rDampingMatrix.size1() != N_DOF) rDampingMatrix.resize(N_DOF, N_DOF, false);
noalias(rDampingMatrix) = ZeroMatrix(N_DOF, N_DOF);
if (rDampingMatrix.size1() != n_DoF) rDampingMatrix.resize(n_DoF, n_DoF, false);
noalias(rDampingMatrix) = ZeroMatrix(n_DoF, n_DoF);

KRATOS_CATCH("")
}
Expand All @@ -72,9 +72,9 @@ void TransientPwElement<TDim, TNumNodes>::GetValuesVector(Vector& rValues, int S
{
KRATOS_TRY

const unsigned int N_DOF = this->GetNumberOfDOF();
const unsigned int n_DoF = this->GetNumberOfDOF();

if (rValues.size() != N_DOF) rValues.resize(N_DOF, false);
if (rValues.size() != n_DoF) rValues.resize(n_DoF, false);

// Why are we constructing a zero vector here?
for (unsigned int i = 0; i < TNumNodes; ++i) {
Expand All @@ -89,9 +89,9 @@ void TransientPwElement<TDim, TNumNodes>::GetFirstDerivativesVector(Vector& rVal
{
KRATOS_TRY

const unsigned int N_DOF = this->GetNumberOfDOF();
const unsigned int n_DoF = this->GetNumberOfDOF();

if (rValues.size() != N_DOF) rValues.resize(N_DOF, false);
if (rValues.size() != n_DoF) rValues.resize(n_DoF, false);

// Why are we constructing a zero vector here?
for (unsigned int i = 0; i < TNumNodes; ++i) {
Expand All @@ -106,9 +106,9 @@ void TransientPwElement<TDim, TNumNodes>::GetSecondDerivativesVector(Vector& rVa
{
KRATOS_TRY

const unsigned int N_DOF = this->GetNumberOfDOF();
const unsigned int n_DoF = this->GetNumberOfDOF();

if (rValues.size() != N_DOF) rValues.resize(N_DOF, false);
if (rValues.size() != n_DoF) rValues.resize(n_DoF, false);

// Why are we constructing a zero vector here?
for (unsigned int i = 0; i < TNumNodes; ++i) {
Expand All @@ -123,20 +123,21 @@ void TransientPwElement<TDim, TNumNodes>::Initialize(const ProcessInfo& rCurrent
{
KRATOS_TRY

const PropertiesType& Prop = this->GetProperties();
const GeometryType& Geom = this->GetGeometry();
const unsigned int NumGPoints = Geom.IntegrationPointsNumber(this->GetIntegrationMethod());
const PropertiesType& r_properties = this->GetProperties();
const GeometryType& r_geom = this->GetGeometry();
const unsigned int number_of_integration_points =
r_geom.IntegrationPointsNumber(this->GetIntegrationMethod());

// pointer to constitutive laws
if (mConstitutiveLawVector.size() != NumGPoints) mConstitutiveLawVector.resize(NumGPoints);

for (unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i) {
mConstitutiveLawVector[i] = nullptr;
if (mConstitutiveLawVector.size() != number_of_integration_points)
mConstitutiveLawVector.resize(number_of_integration_points);
for (auto& constitutive_law : mConstitutiveLawVector) {
constitutive_law = nullptr;
}

if (mRetentionLawVector.size() != NumGPoints) mRetentionLawVector.resize(NumGPoints);
for (unsigned int i = 0; i < mRetentionLawVector.size(); ++i) {
mRetentionLawVector[i] = RetentionLawFactory::Clone(Prop);
if (mRetentionLawVector.size() != number_of_integration_points)
mRetentionLawVector.resize(number_of_integration_points);
for (auto& r_retention_law : mRetentionLawVector) {
r_retention_law = RetentionLawFactory::Clone(r_properties);
}

mIsInitialised = true;
Expand All @@ -149,105 +150,102 @@ int TransientPwElement<TDim, TNumNodes>::Check(const ProcessInfo& rCurrentProces
{
KRATOS_TRY

const PropertiesType& Prop = this->GetProperties();
const GeometryType& Geom = this->GetGeometry();
const PropertiesType& r_properties = this->GetProperties();
const GeometryType& r_geom = this->GetGeometry();

if (Geom.DomainSize() < 1.0e-15)
if (r_geom.DomainSize() < 1.0e-15)
KRATOS_ERROR << "DomainSize < 1.0e-15 for the element " << this->Id() << std::endl;

for (unsigned int i = 0; i < TNumNodes; ++i) {
if (Geom[i].SolutionStepsDataHas(WATER_PRESSURE) == false)
KRATOS_ERROR << "missing variable WATER_PRESSURE on node " << Geom[i].Id() << std::endl;
if (r_geom[i].SolutionStepsDataHas(WATER_PRESSURE) == false)
KRATOS_ERROR << "Missing variable WATER_PRESSURE on node " << r_geom[i].Id() << std::endl;

if (Geom[i].SolutionStepsDataHas(DT_WATER_PRESSURE) == false)
KRATOS_ERROR << "missing variable DT_WATER_PRESSURE on node " << Geom[i].Id() << std::endl;
if (r_geom[i].SolutionStepsDataHas(DT_WATER_PRESSURE) == false)
KRATOS_ERROR << "Missing variable DT_WATER_PRESSURE on node " << r_geom[i].Id() << std::endl;

if (Geom[i].SolutionStepsDataHas(VOLUME_ACCELERATION) == false)
KRATOS_ERROR << "missing variable VOLUME_ACCELERATION on node " << Geom[i].Id() << std::endl;
if (r_geom[i].SolutionStepsDataHas(VOLUME_ACCELERATION) == false)
KRATOS_ERROR << "Missing variable VOLUME_ACCELERATION on node " << r_geom[i].Id() << std::endl;

if (Geom[i].HasDofFor(WATER_PRESSURE) == false)
KRATOS_ERROR << "missing variable WATER_PRESSURE on node " << Geom[i].Id() << std::endl;
if (r_geom[i].HasDofFor(WATER_PRESSURE) == false)
KRATOS_ERROR << "Missing variable WATER_PRESSURE on node " << r_geom[i].Id() << std::endl;
}

// Verify ProcessInfo variables

// Verify properties
if (Prop.Has(DENSITY_WATER) == false || Prop[DENSITY_WATER] < 0.0)
if (r_properties.Has(DENSITY_WATER) == false || r_properties[DENSITY_WATER] < 0.0)
KRATOS_ERROR << "DENSITY_WATER does not exist in the material "
"properties or has an invalid value at element"
"properties or has an invalid value at element "
<< this->Id() << std::endl;

if (Prop.Has(BULK_MODULUS_SOLID) == false || Prop[BULK_MODULUS_SOLID] < 0.0)
if (r_properties.Has(BULK_MODULUS_SOLID) == false || r_properties[BULK_MODULUS_SOLID] < 0.0)
KRATOS_ERROR << "BULK_MODULUS_SOLID does not exist in the material "
"properties or has an invalid value at element"
"properties or has an invalid value at element "
<< this->Id() << std::endl;

if (Prop.Has(POROSITY) == false || Prop[POROSITY] < 0.0 || Prop[POROSITY] > 1.0)
if (r_properties.Has(POROSITY) == false || r_properties[POROSITY] < 0.0 || r_properties[POROSITY] > 1.0)
KRATOS_ERROR << "POROSITY does not exist in the material properties or "
"has an invalid value at element"
"has an invalid value at element "
<< this->Id() << std::endl;

if (TDim == 2) {
// If this is a 2D problem, nodes must be in XY plane
for (unsigned int i = 0; i < TNumNodes; ++i) {
if (Geom[i].Z() != 0.0)
KRATOS_ERROR << " Node with non-zero Z coordinate found. Id: " << Geom[i].Id() << std::endl;
if (r_geom[i].Z() != 0.0)
KRATOS_ERROR << " Node with non-zero Z coordinate found. Id: " << r_geom[i].Id() << std::endl;
}
}

// Verify specific properties
if (Prop.Has(BULK_MODULUS_FLUID) == false || Prop[BULK_MODULUS_FLUID] < 0.0)
if (r_properties.Has(BULK_MODULUS_FLUID) == false || r_properties[BULK_MODULUS_FLUID] < 0.0)
KRATOS_ERROR << "BULK_MODULUS_FLUID does not exist in the material "
"properties or has an invalid value at element"
"properties or has an invalid value at element "
<< this->Id() << std::endl;

if (Prop.Has(DYNAMIC_VISCOSITY) == false || Prop[DYNAMIC_VISCOSITY] < 0.0)
if (r_properties.Has(DYNAMIC_VISCOSITY) == false || r_properties[DYNAMIC_VISCOSITY] < 0.0)
KRATOS_ERROR << "DYNAMIC_VISCOSITY does not exist in the material "
"properties or has an invalid value at element"
"properties or has an invalid value at element "
<< this->Id() << std::endl;

if (Prop.Has(PERMEABILITY_XX) == false || Prop[PERMEABILITY_XX] < 0.0)
if (r_properties.Has(PERMEABILITY_XX) == false || r_properties[PERMEABILITY_XX] < 0.0)
KRATOS_ERROR << "PERMEABILITY_XX does not exist in the material "
"properties or has an invalid value at element"
"properties or has an invalid value at element "
<< this->Id() << std::endl;

if (Prop.Has(PERMEABILITY_YY) == false || Prop[PERMEABILITY_YY] < 0.0)
if (r_properties.Has(PERMEABILITY_YY) == false || r_properties[PERMEABILITY_YY] < 0.0)
KRATOS_ERROR << "PERMEABILITY_YY does not exist in the material "
"properties or has an invalid value at element"
"properties or has an invalid value at element "
<< this->Id() << std::endl;

if (Prop.Has(PERMEABILITY_XY) == false || Prop[PERMEABILITY_XY] < 0.0)
if (r_properties.Has(PERMEABILITY_XY) == false || r_properties[PERMEABILITY_XY] < 0.0)
KRATOS_ERROR << "PERMEABILITY_XY does not exist in the material "
"properties or has an invalid value at element"
"properties or has an invalid value at element "
<< this->Id() << std::endl;

if (!Prop.Has(BIOT_COEFFICIENT))
if (!r_properties.Has(BIOT_COEFFICIENT))
KRATOS_ERROR << "BIOT_COEFFICIENT does not exist in the material "
"properties in element"
"properties in element "
<< this->Id() << std::endl;

if constexpr (TDim > 2) {
if (Prop.Has(PERMEABILITY_ZZ) == false || Prop[PERMEABILITY_ZZ] < 0.0)
if (r_properties.Has(PERMEABILITY_ZZ) == false || r_properties[PERMEABILITY_ZZ] < 0.0)
KRATOS_ERROR << "PERMEABILITY_ZZ does not exist in the material "
"properties or has an invalid value at element"
"properties or has an invalid value at element "
<< this->Id() << std::endl;

if (Prop.Has(PERMEABILITY_YZ) == false || Prop[PERMEABILITY_YZ] < 0.0)
if (r_properties.Has(PERMEABILITY_YZ) == false || r_properties[PERMEABILITY_YZ] < 0.0)
KRATOS_ERROR << "PERMEABILITY_YZ does not exist in the material "
"properties or has an invalid value at element"
"properties or has an invalid value at element "
<< this->Id() << std::endl;

if (Prop.Has(PERMEABILITY_ZX) == false || Prop[PERMEABILITY_ZX] < 0.0)
if (r_properties.Has(PERMEABILITY_ZX) == false || r_properties[PERMEABILITY_ZX] < 0.0)
KRATOS_ERROR << "PERMEABILITY_ZX does not exist in the material "
"properties or has an invalid value at element"
"properties or has an invalid value at element "
<< this->Id() << std::endl;
}

// Verify that the constitutive law has the correct dimension

// Check constitutive law
if (mRetentionLawVector.size() > 0) {
return mRetentionLawVector[0]->Check(Prop, rCurrentProcessInfo);
return mRetentionLawVector[0]->Check(r_properties, rCurrentProcessInfo);
}

return 0;
Expand Down Expand Up @@ -318,23 +316,24 @@ void TransientPwElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(const Var
{
KRATOS_TRY

const GeometryType& rGeom = this->GetGeometry();
const IndexType NumGPoints = rGeom.IntegrationPointsNumber(this->GetIntegrationMethod());
if (rOutput.size() != NumGPoints) rOutput.resize(NumGPoints);
const GeometryType& r_geom = this->GetGeometry();
const IndexType number_of_integration_points =
r_geom.IntegrationPointsNumber(this->GetIntegrationMethod());
if (rOutput.size() != number_of_integration_points)
rOutput.resize(number_of_integration_points);

if (rVariable == FLUID_FLUX_VECTOR) {
std::vector<double> permeability_update_factors(NumGPoints, 1.0);
std::vector<double> permeability_update_factors(number_of_integration_points, 1.0);
const auto fluid_fluxes = this->CalculateFluidFluxes(permeability_update_factors, rCurrentProcessInfo);

for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], fluid_fluxes[GPoint]);
for (unsigned int integration_point = 0; integration_point < number_of_integration_points;
++integration_point) {
GeoElementUtilities::FillArray1dOutput(rOutput[integration_point], fluid_fluxes[integration_point]);
}
} else {
if (rOutput.size() != mRetentionLawVector.size())
rOutput.resize(mRetentionLawVector.size());

for (unsigned int i = 0; i < mRetentionLawVector.size(); ++i) {
noalias(rOutput[i]) = ZeroVector(3);
for (unsigned int integration_point = 0; integration_point < number_of_integration_points;
++integration_point) {
noalias(rOutput[integration_point]) = ZeroVector(3);
}
}

Expand Down Expand Up @@ -373,11 +372,11 @@ void TransientPwElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftH
KRATOS_TRY

// Previous definitions
const PropertiesType& Prop = this->GetProperties();
const GeometryType& Geom = this->GetGeometry();
const PropertiesType& r_properties = this->GetProperties();
const GeometryType& r_geom = this->GetGeometry();
const GeometryType::IntegrationPointsArrayType& IntegrationPoints =
Geom.IntegrationPoints(this->GetIntegrationMethod());
const unsigned int NumGPoints = IntegrationPoints.size();
r_geom.IntegrationPoints(this->GetIntegrationMethod());
const unsigned int number_of_integration_points = IntegrationPoints.size();

// Element variables
ElementVariables Variables;
Expand All @@ -391,37 +390,37 @@ void TransientPwElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftH
const auto bishop_coefficients = this->CalculateBishopCoefficients(fluid_pressures);
const auto integration_coefficients =
this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer);
std::vector<double> biot_coefficients(NumGPoints, Prop[BIOT_COEFFICIENT]);
const auto degrees_of_saturation = this->CalculateDegreesOfSaturation(fluid_pressures);
std::vector<double> biot_coefficients(number_of_integration_points, r_properties[BIOT_COEFFICIENT]);
const auto degrees_of_saturation = this->CalculateDegreesOfSaturation(fluid_pressures);
const auto derivatives_of_saturation = this->CalculateDerivativesOfSaturation(fluid_pressures);
const auto biot_moduli_inverse = GeoTransportEquationUtilities::CalculateInverseBiotModuli(
biot_coefficients, degrees_of_saturation, derivatives_of_saturation, Prop);
biot_coefficients, degrees_of_saturation, derivatives_of_saturation, r_properties);

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
for (unsigned int integration_point = 0; integration_point < number_of_integration_points; ++integration_point) {
// Compute GradNpT, B and StrainVector
this->CalculateKinematics(Variables, GPoint);
this->CalculateKinematics(Variables, integration_point);

// Compute Nu and BodyAcceleration
GeoElementUtilities::CalculateNuMatrix<TDim, TNumNodes>(Variables.Nu, Variables.NContainer, GPoint);
GeoElementUtilities::CalculateNuMatrix<TDim, TNumNodes>(Variables.Nu, Variables.NContainer, integration_point);
GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint);
Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, integration_point);

Variables.RelativePermeability = relative_permeability_values[GPoint];
Variables.BishopCoefficient = bishop_coefficients[GPoint];
Variables.RelativePermeability = relative_permeability_values[integration_point];
Variables.BishopCoefficient = bishop_coefficients[integration_point];

Variables.BiotCoefficient = biot_coefficients[GPoint];
Variables.BiotModulusInverse = biot_moduli_inverse[GPoint];
Variables.DegreeOfSaturation = degrees_of_saturation[GPoint];
Variables.BiotCoefficient = biot_coefficients[integration_point];
Variables.BiotModulusInverse = biot_moduli_inverse[integration_point];
Variables.DegreeOfSaturation = degrees_of_saturation[integration_point];

Variables.IntegrationCoefficient = integration_coefficients[GPoint];
Variables.IntegrationCoefficient = integration_coefficients[integration_point];

// Contributions to the left hand side
if (CalculateStiffnessMatrixFlag) this->CalculateAndAddLHS(rLeftHandSideMatrix, Variables);

// Contributions to the right hand side
if (CalculateResidualVectorFlag)
this->CalculateAndAddRHS(rRightHandSideVector, Variables, GPoint);
this->CalculateAndAddRHS(rRightHandSideVector, Variables, integration_point);
}

KRATOS_CATCH("")
Expand Down Expand Up @@ -449,17 +448,18 @@ void TransientPwElement<TDim, TNumNodes>::InitializeElementVariables(ElementVari
rVariables.GradNpT.resize(TNumNodes, TDim, false);

// General Variables
const GeometryType& Geom = this->GetGeometry();
const unsigned int NumGPoints = Geom.IntegrationPointsNumber(this->GetIntegrationMethod());
const GeometryType& r_geom = this->GetGeometry();
const unsigned int number_of_integration_points =
r_geom.IntegrationPointsNumber(this->GetIntegrationMethod());

// shape functions
(rVariables.NContainer).resize(NumGPoints, TNumNodes, false);
rVariables.NContainer = Geom.ShapeFunctionsValues(this->GetIntegrationMethod());
(rVariables.NContainer).resize(number_of_integration_points, TNumNodes, false);
rVariables.NContainer = r_geom.ShapeFunctionsValues(this->GetIntegrationMethod());

// gradient of shape functions and determinant of Jacobian
(rVariables.detJContainer).resize(NumGPoints, false);
(rVariables.detJContainer).resize(number_of_integration_points, false);

Geom.ShapeFunctionsIntegrationPointsGradients(
r_geom.ShapeFunctionsIntegrationPointsGradients(
rVariables.DN_DXContainer, rVariables.detJContainer, this->GetIntegrationMethod());

// Retention law
Expand Down Expand Up @@ -503,7 +503,7 @@ void TransientPwElement<TDim, TNumNodes>::CalculateAndAddCompressibilityMatrix(M
template <unsigned int TDim, unsigned int TNumNodes>
void TransientPwElement<TDim, TNumNodes>::CalculateAndAddRHS(VectorType& rRightHandSideVector,
ElementVariables& rVariables,
unsigned int GPoint)
unsigned int integration_point)
{
KRATOS_TRY

Expand Down
Loading
Loading