From 356492913ed6d05e0cd34ae640e5a2713fc9dc79 Mon Sep 17 00:00:00 2001 From: Richard Faasse Date: Thu, 2 May 2024 15:13:00 +0200 Subject: [PATCH 1/7] Copied fluid flux vector calculation to pw element --- .../custom_elements/transient_Pw_element.cpp | 44 ++++++++++++++++++- 1 file changed, 43 insertions(+), 1 deletion(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp index 7b7c60338049..432f4768e47d 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp @@ -383,8 +383,50 @@ void TransientPwElement::CalculateOnIntegrationPoints(const Var { KRATOS_TRY + const GeometryType& rGeom = this->GetGeometry(); + const IndexType NumGPoints = rGeom.IntegrationPointsNumber(this->GetIntegrationMethod()); + if (rOutput.size() != NumGPoints) rOutput.resize(NumGPoints); + if (rVariable == FLUID_FLUX_VECTOR) { - BaseType::CalculateOnIntegrationPoints(rVariable, rOutput, rCurrentProcessInfo); + ElementVariables Variables; + this->InitializeElementVariables(Variables, rCurrentProcessInfo); + + const PropertiesType& rProp = this->GetProperties(); + + array_1d GradPressureTerm; + array_1d FluidFlux; + + // Create general parameters of retention law + RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + + for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { + // Compute Np, GradNpT, B and StrainVector + this->CalculateKinematics(Variables, GPoint); + + // Compute strain, needed for update permeability + this->CalculateStrain(Variables, GPoint); + this->CalculatePermeabilityUpdateFactor(Variables); + + GeoElementUtilities::InterpolateVariableWithComponents( + Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); + Variables.FluidPressure = CalculateFluidPressure(Variables); + + RetentionParameters.SetFluidPressure(Variables.FluidPressure); + + Variables.RelativePermeability = + mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); + + noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector); + + noalias(GradPressureTerm) += + PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; + + noalias(FluidFlux) = PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * + Variables.RelativePermeability * Variables.PermeabilityUpdateFactor * + prod(Variables.PermeabilityMatrix, GradPressureTerm); + + GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], FluidFlux); + } } else { if (rOutput.size() != mRetentionLawVector.size()) rOutput.resize(mRetentionLawVector.size()); From 4c027647d8e65f505caff7d49c814976746f95e9 Mon Sep 17 00:00:00 2001 From: Richard Faasse Date: Thu, 2 May 2024 16:27:08 +0200 Subject: [PATCH 2/7] Extract `CalculateFluidFluxes` which takes permeability update factors as an input list --- .../U_Pw_small_strain_element.cpp | 42 +++++++++++++++++++ .../U_Pw_small_strain_element.hpp | 5 ++- .../custom_elements/transient_Pw_element.cpp | 39 ++--------------- 3 files changed, 49 insertions(+), 37 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 4e6aff89455a..e1afe5bc5720 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -1112,6 +1112,48 @@ void UPwSmallStrainElement::InitializeBiotCoefficients(ElementV KRATOS_CATCH("") } +template +std::vector> UPwSmallStrainElement::CalculateFluidFluxes( + const std::vector& permeability_update_factors, const ProcessInfo& rCurrentProcessInfo) +{ + const GeometryType& rGeom = this->GetGeometry(); + const IndexType NumGPoints = rGeom.IntegrationPointsNumber(this->GetIntegrationMethod()); + + std::vector> FluidFluxes; + ElementVariables Variables; + this->InitializeElementVariables(Variables, rCurrentProcessInfo); + + const PropertiesType& rProp = this->GetProperties(); + + array_1d GradPressureTerm; + + // Create general parameters of retention law + RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + + for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { + this->CalculateKinematics(Variables, GPoint); + Variables.PermeabilityUpdateFactor = permeability_update_factors[GPoint]; + + GeoElementUtilities::InterpolateVariableWithComponents( + Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); + Variables.FluidPressure = CalculateFluidPressure(Variables); + RetentionParameters.SetFluidPressure(Variables.FluidPressure); + + Variables.RelativePermeability = + mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); + + noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector); + + noalias(GradPressureTerm) += PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; + + FluidFluxes.push_back(PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * + Variables.RelativePermeability * Variables.PermeabilityUpdateFactor * + prod(Variables.PermeabilityMatrix, GradPressureTerm)); + } + + return FluidFluxes; +} + template void UPwSmallStrainElement::CalculatePermeabilityUpdateFactor(ElementVariables& rVariables) { 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 3217c3b4527e..c42752041445 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -240,7 +240,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa virtual void CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables); - virtual void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, const ElementVariables& rVariables); + virtual void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, + const ElementVariables& rVariables); virtual void CalculateAndAddRHS(VectorType& rRightHandSideVector, ElementVariables& rVariables, unsigned int GPoint); @@ -283,6 +284,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa void InitializeProperties(ElementVariables& rVariables); double CalculateFluidPressure(const ElementVariables& rVariables); + std::vector> CalculateFluidFluxes(const std::vector& permeability_update_factors, + const ProcessInfo& rCurrentProcessInfo); void CalculateRetentionResponse(ElementVariables& rVariables, RetentionLaw::Parameters& rRetentionParameters, diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp index 432f4768e47d..668324dacce2 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp @@ -388,44 +388,11 @@ void TransientPwElement::CalculateOnIntegrationPoints(const Var if (rOutput.size() != NumGPoints) rOutput.resize(NumGPoints); if (rVariable == FLUID_FLUX_VECTOR) { - ElementVariables Variables; - this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - const PropertiesType& rProp = this->GetProperties(); - - array_1d GradPressureTerm; - array_1d FluidFlux; - - // Create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + std::vector permeability_update_factors(NumGPoints, 1.0); + const auto FluidFluxes = this->CalculateFluidFluxes(permeability_update_factors, rCurrentProcessInfo); for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - // Compute Np, GradNpT, B and StrainVector - this->CalculateKinematics(Variables, GPoint); - - // Compute strain, needed for update permeability - this->CalculateStrain(Variables, GPoint); - this->CalculatePermeabilityUpdateFactor(Variables); - - GeoElementUtilities::InterpolateVariableWithComponents( - Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - Variables.FluidPressure = CalculateFluidPressure(Variables); - - RetentionParameters.SetFluidPressure(Variables.FluidPressure); - - Variables.RelativePermeability = - mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); - - noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector); - - noalias(GradPressureTerm) += - PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; - - noalias(FluidFlux) = PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * - Variables.RelativePermeability * Variables.PermeabilityUpdateFactor * - prod(Variables.PermeabilityMatrix, GradPressureTerm); - - GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], FluidFlux); + GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], FluidFluxes[GPoint]); } } else { if (rOutput.size() != mRetentionLawVector.size()) From 4a5056ee9514daea3887190c2a3255c40a28d039 Mon Sep 17 00:00:00 2001 From: Richard Faasse Date: Thu, 2 May 2024 17:16:10 +0200 Subject: [PATCH 3/7] Use `CalculateFluidFluxes` in U_Pw_small_strain_element.cpp --- .../U_Pw_small_strain_element.cpp | 39 +++++-------------- .../custom_elements/transient_Pw_element.cpp | 4 +- 2 files changed, 12 insertions(+), 31 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 e1afe5bc5720..d2d90863c606 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -620,15 +620,8 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints( if (rVariable == FLUID_FLUX_VECTOR) { ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - const PropertiesType& rProp = this->GetProperties(); - - array_1d GradPressureTerm; - array_1d FluidFlux; - - // Create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); - + + std::vector permeability_update_factors; for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); @@ -636,28 +629,16 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints( // Compute strain, needed for update permeability this->CalculateStrain(Variables, GPoint); this->CalculatePermeabilityUpdateFactor(Variables); + permeability_update_factors.push_back(Variables.PermeabilityUpdateFactor); + } - GeoElementUtilities::InterpolateVariableWithComponents( - Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - Variables.FluidPressure = CalculateFluidPressure(Variables); - - RetentionParameters.SetFluidPressure(Variables.FluidPressure); - - Variables.RelativePermeability = - mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); - - noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector); - - noalias(GradPressureTerm) += - PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; - - noalias(FluidFlux) = PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * - Variables.RelativePermeability * Variables.PermeabilityUpdateFactor * - prod(Variables.PermeabilityMatrix, GradPressureTerm); - - GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], FluidFlux); + const auto fluid_fluxes = CalculateFluidFluxes(permeability_update_factors, rCurrentProcessInfo); + for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { + GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], fluid_fluxes[GPoint]); } - } else { + } + + else { if (rOutput.size() != mConstitutiveLawVector.size()) rOutput.resize(mConstitutiveLawVector.size()); diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp index 668324dacce2..1a009644c290 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp @@ -389,10 +389,10 @@ void TransientPwElement::CalculateOnIntegrationPoints(const Var if (rVariable == FLUID_FLUX_VECTOR) { std::vector permeability_update_factors(NumGPoints, 1.0); - const auto FluidFluxes = this->CalculateFluidFluxes(permeability_update_factors, rCurrentProcessInfo); + const auto fluid_fluxes = this->CalculateFluidFluxes(permeability_update_factors, rCurrentProcessInfo); for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], FluidFluxes[GPoint]); + GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], fluid_fluxes[GPoint]); } } else { if (rOutput.size() != mRetentionLawVector.size()) From 14f82805482c083fcc8a76a78e9cfd783383e7d1 Mon Sep 17 00:00:00 2001 From: Richard Faasse Date: Thu, 2 May 2024 17:29:57 +0200 Subject: [PATCH 4/7] Cleaning up and simplifying CalculateFluidFluxes --- .../U_Pw_small_strain_element.cpp | 27 ++++++++----------- 1 file changed, 11 insertions(+), 16 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 d2d90863c606..daa3f4af2cc6 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -620,7 +620,7 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints( if (rVariable == FLUID_FLUX_VECTOR) { ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - + std::vector permeability_update_factors; for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { // Compute Np, GradNpT, B and StrainVector @@ -1100,39 +1100,34 @@ std::vector> UPwSmallStrainElement::Calc const GeometryType& rGeom = this->GetGeometry(); const IndexType NumGPoints = rGeom.IntegrationPointsNumber(this->GetIntegrationMethod()); - std::vector> FluidFluxes; - ElementVariables Variables; + ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); const PropertiesType& rProp = this->GetProperties(); - array_1d GradPressureTerm; - // Create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + std::vector> result; for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { this->CalculateKinematics(Variables, GPoint); - Variables.PermeabilityUpdateFactor = permeability_update_factors[GPoint]; GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - Variables.FluidPressure = CalculateFluidPressure(Variables); - RetentionParameters.SetFluidPressure(Variables.FluidPressure); + RetentionParameters.SetFluidPressure(CalculateFluidPressure(Variables)); - Variables.RelativePermeability = + const auto relative_permeability = mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); - noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector); - - noalias(GradPressureTerm) += PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; + const auto grad_pressure_term = prod(trans(Variables.GradNpT), Variables.PressureVector) + + PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; - FluidFluxes.push_back(PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * - Variables.RelativePermeability * Variables.PermeabilityUpdateFactor * - prod(Variables.PermeabilityMatrix, GradPressureTerm)); + result.push_back(PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * + relative_permeability * permeability_update_factors[GPoint] * + prod(Variables.PermeabilityMatrix, grad_pressure_term)); } - return FluidFluxes; + return result; } template From e3d770cf00986da5433d9d7ff7e9e73d58b94504 Mon Sep 17 00:00:00 2001 From: Richard Faasse Date: Thu, 2 May 2024 17:40:23 +0200 Subject: [PATCH 5/7] Minor formatting commit --- .../custom_elements/U_Pw_small_strain_element.cpp | 4 +--- 1 file changed, 1 insertion(+), 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 daa3f4af2cc6..1b0a96b2208b 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -636,9 +636,7 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints( for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], fluid_fluxes[GPoint]); } - } - - else { + } else { if (rOutput.size() != mConstitutiveLawVector.size()) rOutput.resize(mConstitutiveLawVector.size()); From f4e155afe2f1e00ab7d0456b7c79dd1135cf4f54 Mon Sep 17 00:00:00 2001 From: Richard Faasse Date: Fri, 3 May 2024 12:58:56 +0200 Subject: [PATCH 6/7] Revert "Cleaning up and simplifying CalculateFluidFluxes" This reverts commit 14f82805482c083fcc8a76a78e9cfd783383e7d1. --- .../U_Pw_small_strain_element.cpp | 27 +++++++++++-------- 1 file changed, 16 insertions(+), 11 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 1b0a96b2208b..a703074b7f41 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -620,7 +620,7 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints( if (rVariable == FLUID_FLUX_VECTOR) { ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - + std::vector permeability_update_factors; for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { // Compute Np, GradNpT, B and StrainVector @@ -1098,34 +1098,39 @@ std::vector> UPwSmallStrainElement::Calc const GeometryType& rGeom = this->GetGeometry(); const IndexType NumGPoints = rGeom.IntegrationPointsNumber(this->GetIntegrationMethod()); - ElementVariables Variables; + std::vector> FluidFluxes; + ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); const PropertiesType& rProp = this->GetProperties(); + array_1d GradPressureTerm; + // Create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); - std::vector> result; for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { this->CalculateKinematics(Variables, GPoint); + Variables.PermeabilityUpdateFactor = permeability_update_factors[GPoint]; GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - RetentionParameters.SetFluidPressure(CalculateFluidPressure(Variables)); + Variables.FluidPressure = CalculateFluidPressure(Variables); + RetentionParameters.SetFluidPressure(Variables.FluidPressure); - const auto relative_permeability = + Variables.RelativePermeability = mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); - const auto grad_pressure_term = prod(trans(Variables.GradNpT), Variables.PressureVector) + - PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; + noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector); + + noalias(GradPressureTerm) += PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; - result.push_back(PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * - relative_permeability * permeability_update_factors[GPoint] * - prod(Variables.PermeabilityMatrix, grad_pressure_term)); + FluidFluxes.push_back(PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * + Variables.RelativePermeability * Variables.PermeabilityUpdateFactor * + prod(Variables.PermeabilityMatrix, GradPressureTerm)); } - return result; + return FluidFluxes; } template From c25296eb2d505b5730b9389a9f9e5c2c68f5991a Mon Sep 17 00:00:00 2001 From: Richard Faasse Date: Fri, 3 May 2024 17:19:07 +0200 Subject: [PATCH 7/7] Processing review comments --- .../custom_elements/U_Pw_small_strain_element.cpp | 7 +++---- .../custom_elements/U_Pw_small_strain_element.hpp | 2 +- 2 files changed, 4 insertions(+), 5 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 dcc7b4849c92..57c8c11581dd 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -620,7 +620,7 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints( if (rVariable == FLUID_FLUX_VECTOR) { ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - + std::vector permeability_update_factors; for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { // Compute Np, GradNpT, B and StrainVector @@ -1093,7 +1093,7 @@ void UPwSmallStrainElement::InitializeBiotCoefficients(ElementV template std::vector> UPwSmallStrainElement::CalculateFluidFluxes( - const std::vector& permeability_update_factors, const ProcessInfo& rCurrentProcessInfo) + const std::vector& rPermeabilityUpdateFactors, const ProcessInfo& rCurrentProcessInfo) { const GeometryType& rGeom = this->GetGeometry(); const IndexType NumGPoints = rGeom.IntegrationPointsNumber(this->GetIntegrationMethod()); @@ -1111,7 +1111,7 @@ std::vector> UPwSmallStrainElement::Calc for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { this->CalculateKinematics(Variables, GPoint); - Variables.PermeabilityUpdateFactor = permeability_update_factors[GPoint]; + Variables.PermeabilityUpdateFactor = rPermeabilityUpdateFactors[GPoint]; GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); @@ -1122,7 +1122,6 @@ std::vector> UPwSmallStrainElement::Calc mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector); - noalias(GradPressureTerm) += PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; FluidFluxes.push_back(PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * 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 637a1f55c524..933c9d54cee8 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -284,7 +284,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa void InitializeProperties(ElementVariables& rVariables); double CalculateFluidPressure(const ElementVariables& rVariables); - std::vector> CalculateFluidFluxes(const std::vector& permeability_update_factors, + std::vector> CalculateFluidFluxes(const std::vector& rPermeabilityUpdateFactors, const ProcessInfo& rCurrentProcessInfo); void CalculateRetentionResponse(ElementVariables& rVariables,