diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp index a5e3ffbad047..f27964482d6d 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp @@ -135,10 +135,8 @@ void UPwBaseElement::Initialize(const ProcessInfo& rCurrentProcessInfo) } mRetentionLawVector.resize(number_of_integration_points); - for (unsigned int i = 0; i < mRetentionLawVector.size(); ++i) { - mRetentionLawVector[i] = RetentionLawFactory::Clone(r_properties); - mRetentionLawVector[i]->InitializeMaterial( - r_properties, r_geometry, row(r_geometry.ShapeFunctionsValues(mThisIntegrationMethod), i)); + for (auto& r_retention_law : mRetentionLawVector) { + r_retention_law = RetentionLawFactory::Clone(r_properties); } if (mStressVector.size() != number_of_integration_points) { @@ -151,9 +149,8 @@ void UPwBaseElement::Initialize(const ProcessInfo& rCurrentProcessInfo) mStateVariablesFinalized.resize(number_of_integration_points); for (unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i) { - int nStateVariables = 0; - nStateVariables = mConstitutiveLawVector[i]->GetValue(NUMBER_OF_UMAT_STATE_VARIABLES, nStateVariables); - if (nStateVariables > 0) { + if (auto nStateVariables = 0; + mConstitutiveLawVector[i]->GetValue(NUMBER_OF_UMAT_STATE_VARIABLES, nStateVariables) > 0) { mConstitutiveLawVector[i]->SetValue(STATE_VARIABLES, mStateVariablesFinalized[i], rCurrentProcessInfo); } } 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 0e890df34ecd..95279b35bc45 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -258,8 +258,6 @@ void UPwSmallStrainElement::InitializeSolutionStep(const Proces ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - RetentionLaw::Parameters RetentionParameters(this->GetProperties()); - const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); const auto determinants_of_deformation_gradients = @@ -283,8 +281,6 @@ void UPwSmallStrainElement::InitializeSolutionStep(const Proces noalias(Variables.StressVector) = mStressVector[GPoint]; ConstitutiveParameters.SetStressVector(Variables.StressVector); mConstitutiveLawVector[GPoint]->InitializeMaterialResponseCauchy(ConstitutiveParameters); - - mRetentionLawVector[GPoint]->InitializeSolutionStep(RetentionParameters); } // Reset hydraulic discharge @@ -299,9 +295,8 @@ void UPwSmallStrainElement::ResetHydraulicDischarge() KRATOS_TRY // Reset hydraulic discharge - GeometryType& rGeom = this->GetGeometry(); - for (unsigned int i = 0; i < TNumNodes; ++i) { - ThreadSafeNodeWrite(rGeom[i], HYDRAULIC_DISCHARGE, 0.0); + for (auto& r_node : this->GetGeometry()) { + ThreadSafeNodeWrite(r_node, HYDRAULIC_DISCHARGE, 0.0); } KRATOS_CATCH("") @@ -357,8 +352,8 @@ void UPwSmallStrainElement::InitializeNonLinearIteration(const { KRATOS_TRY - const GeometryType& rGeom = this->GetGeometry(); - ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom, this->GetProperties(), rCurrentProcessInfo); + ConstitutiveLaw::Parameters ConstitutiveParameters(this->GetGeometry(), this->GetProperties(), + rCurrentProcessInfo); ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_STRESS); ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); ConstitutiveParameters.Set(ConstitutiveLaw::INITIALIZE_MATERIAL_RESPONSE); // Note: this is for nonlocal damage @@ -403,8 +398,6 @@ void UPwSmallStrainElement::FinalizeSolutionStep(const ProcessI ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - RetentionLaw::Parameters RetentionParameters(this->GetProperties()); - const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); const auto determinants_of_deformation_gradients = @@ -434,8 +427,6 @@ void UPwSmallStrainElement::FinalizeSolutionStep(const ProcessI mStateVariablesFinalized[GPoint] = mConstitutiveLawVector[GPoint]->GetValue(STATE_VARIABLES, mStateVariablesFinalized[GPoint]); - mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters); - if (rCurrentProcessInfo[NODAL_SMOOTHING]) this->SaveGPStress(StressContainer, mStressVector[GPoint], GPoint); } @@ -622,19 +613,8 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure( Variables.Np, Variables.PressureVector)); - - if (rVariable == DEGREE_OF_SATURATION) - rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateSaturation(RetentionParameters); - else if (rVariable == EFFECTIVE_SATURATION) - rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(RetentionParameters); - else if (rVariable == BISHOP_COEFFICIENT) - rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); - else if (rVariable == DERIVATIVE_OF_SATURATION) - rOutput[GPoint] = - mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(RetentionParameters); - else if (rVariable == RELATIVE_PERMEABILITY) - rOutput[GPoint] = - mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); + rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateValue( + RetentionParameters, rVariable, rOutput[GPoint]); } } else if (rVariable == HYDRAULIC_HEAD) { const PropertiesType& rProp = this->GetProperties(); @@ -1246,16 +1226,17 @@ void UPwSmallStrainElement::CalculateAndAddLHS(MatrixType& rLef KRATOS_TRY this->CalculateAndAddStiffnessMatrix(rLeftHandSideMatrix, rVariables); - this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix, rVariables); - if (!rVariables.IgnoreUndrained) { - this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables); + this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables); + if (!rVariables.IgnoreUndrained) { const auto permeability_matrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix, rVariables.RelativePermeability, rVariables.IntegrationCoefficient); GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, permeability_matrix); + + this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix, rVariables); } KRATOS_CATCH("") @@ -1280,17 +1261,18 @@ void UPwSmallStrainElement::CalculateAndAddCouplingMatrix(Matri { KRATOS_TRY - const auto coupling_matrix = GeoTransportEquationUtilities::CalculateCouplingMatrix( - rVariables.B, this->GetStressStatePolicy().GetVoigtVector(), rVariables.Np, + const Matrix coupling_matrix_up = GeoTransportEquationUtilities::CalculateCouplingMatrix( + rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np, rVariables.BiotCoefficient, rVariables.BishopCoefficient, rVariables.IntegrationCoefficient); - GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, coupling_matrix); + GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, coupling_matrix_up); if (!rVariables.IgnoreUndrained) { - const double SaturationCoefficient = rVariables.DegreeOfSaturation / rVariables.BishopCoefficient; - const BoundedMatrix transposed_coupling_matrix = - PORE_PRESSURE_SIGN_FACTOR * SaturationCoefficient * rVariables.VelocityCoefficient * - trans(coupling_matrix); - GeoElementUtilities::AssemblePUBlockMatrix(rLeftHandSideMatrix, transposed_coupling_matrix); + const Matrix coupling_matrix_pu = GeoTransportEquationUtilities::CalculateCouplingMatrix( + rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np, + rVariables.BiotCoefficient, rVariables.DegreeOfSaturation, rVariables.IntegrationCoefficient); + GeoElementUtilities::AssemblePUBlockMatrix( + rLeftHandSideMatrix, + PORE_PRESSURE_SIGN_FACTOR * rVariables.VelocityCoefficient * trans(coupling_matrix_pu)); } KRATOS_CATCH("") @@ -1389,15 +1371,18 @@ void UPwSmallStrainElement::CalculateAndAddCouplingTerms(Vector rVariables.B, this->GetStressStatePolicy().GetVoigtVector(), rVariables.Np, rVariables.BiotCoefficient, rVariables.BishopCoefficient, rVariables.IntegrationCoefficient); - const array_1d coupling_terms = prod(coupling_matrix, rVariables.PressureVector); - GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, coupling_terms); + const array_1d coupling_force = prod(coupling_matrix, rVariables.PressureVector); + GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, coupling_force); if (!rVariables.IgnoreUndrained) { - const double SaturationCoefficient = rVariables.DegreeOfSaturation / rVariables.BishopCoefficient; - const array_1d transposed_coupling_matrix = - PORE_PRESSURE_SIGN_FACTOR * SaturationCoefficient * - prod(trans(coupling_matrix), rVariables.VelocityVector); - GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, transposed_coupling_matrix); + const Matrix p_coupling_matrix = + (-1.0) * GeoTransportEquationUtilities::CalculateCouplingMatrix( + rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np, + rVariables.BiotCoefficient, rVariables.DegreeOfSaturation, + rVariables.IntegrationCoefficient); + const array_1d coupling_flow = + PORE_PRESSURE_SIGN_FACTOR * prod(trans(p_coupling_matrix), rVariables.VelocityVector); + GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, coupling_flow); } KRATOS_CATCH("") @@ -1493,7 +1478,7 @@ array_1d UPwSmallStrainElement::CalculateFlu const BoundedMatrix temp_matrix = prod(rVariables.GradNpT, rVariables.PermeabilityMatrix) * rVariables.IntegrationCoefficient; - return rVariables.DynamicViscosityInverse * this->GetProperties()[DENSITY_WATER] * + return rVariables.DynamicViscosityInverse * this->GetProperties()[DENSITY_WATER] * rVariables.BishopCoefficient * rVariables.RelativePermeability * prod(temp_matrix, rVariables.BodyAcceleration); KRATOS_CATCH("") 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 5a815b97e38e..da678158304a 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 @@ -192,8 +192,6 @@ void UPwSmallStrainInterfaceElement::InitializeSolutionStep(con unsigned int NumGPoints = mConstitutiveLawVector.size(); std::vector JointWidthContainer(NumGPoints); - RetentionLaw::Parameters RetentionParameters(this->GetProperties()); - // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { InterfaceElementUtilities::CalculateNuMatrix(Nu, NContainer, GPoint); @@ -206,9 +204,6 @@ void UPwSmallStrainInterfaceElement::InitializeSolutionStep(con noalias(StressVector) = mStressVector[GPoint]; ConstitutiveParameters.SetStressVector(StressVector); mConstitutiveLawVector[GPoint]->InitializeMaterialResponseCauchy(ConstitutiveParameters); - - // Initialize retention law - mRetentionLawVector[GPoint]->InitializeSolutionStep(RetentionParameters); } KRATOS_CATCH("") @@ -253,8 +248,6 @@ void UPwSmallStrainInterfaceElement::FinalizeSolutionStep(const unsigned int NumGPoints = mConstitutiveLawVector.size(); std::vector JointWidthContainer(NumGPoints); - RetentionLaw::Parameters RetentionParameters(this->GetProperties()); - // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { InterfaceElementUtilities::CalculateNuMatrix(Nu, NContainer, GPoint); @@ -279,9 +272,6 @@ void UPwSmallStrainInterfaceElement::FinalizeSolutionStep(const mStateVariablesFinalized[GPoint] = mConstitutiveLawVector[GPoint]->GetValue(STATE_VARIABLES, mStateVariablesFinalized[GPoint]); - - // retention law - mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters); } if (rCurrentProcessInfo[NODAL_SMOOTHING]) this->ExtrapolateGPValues(JointWidthContainer); 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 8c8d588186cd..0f893eb53018 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 @@ -24,7 +24,6 @@ // Application includes #include "custom_elements/small_strain_U_Pw_diff_order_element.hpp" -#include "custom_retention/retention_law_factory.h" #include "custom_utilities/constitutive_law_utilities.hpp" #include "custom_utilities/dof_utilities.h" #include "custom_utilities/element_utilities.hpp" @@ -40,25 +39,24 @@ Element::Pointer SmallStrainUPwDiffOrderElement::Create(IndexType NodesArrayType const& ThisNodes, PropertiesType::Pointer pProperties) const { - return Element::Pointer(new SmallStrainUPwDiffOrderElement( - NewId, GetGeometry().Create(ThisNodes), pProperties, this->GetStressStatePolicy().Clone())); + return Create(NewId, GetGeometry().Create(ThisNodes), pProperties); } Element::Pointer SmallStrainUPwDiffOrderElement::Create(IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const { - return Element::Pointer(new SmallStrainUPwDiffOrderElement( - NewId, pGeom, pProperties, this->GetStressStatePolicy().Clone())); + return make_intrusive(NewId, pGeom, pProperties, + this->GetStressStatePolicy().Clone()); } int SmallStrainUPwDiffOrderElement::Check(const ProcessInfo& rCurrentProcessInfo) const { KRATOS_TRY - const GeometryType& rGeom = GetGeometry(); + const auto& r_geom = GetGeometry(); - if (rGeom.DomainSize() < 1.0e-15) + if (r_geom.DomainSize() < 1.0e-15) KRATOS_ERROR << "DomainSize < 1.0e-15 for the element " << this->Id() << std::endl; // check pressure geometry pointer @@ -66,39 +64,39 @@ int SmallStrainUPwDiffOrderElement::Check(const ProcessInfo& rCurrentProcessInfo // verify that the variables are correctly initialized // Verify specific properties - const PropertiesType& rProp = this->GetProperties(); + const auto& r_prop = this->GetProperties(); - if (!rProp.Has(IGNORE_UNDRAINED)) + if (!r_prop.Has(IGNORE_UNDRAINED)) KRATOS_ERROR << "IGNORE_UNDRAINED does not exist in the parameter list" << this->Id() << std::endl; - if (!rProp[IGNORE_UNDRAINED]) { - if (!rProp.Has(PERMEABILITY_XX) || rProp[PERMEABILITY_XX] < 0.0) + if (!r_prop[IGNORE_UNDRAINED]) { + if (!r_prop.Has(PERMEABILITY_XX) || r_prop[PERMEABILITY_XX] < 0.0) KRATOS_ERROR << "PERMEABILITY_XX has Key zero, is not defined or " "has an invalid value at element" << this->Id() << std::endl; - if (!rProp.Has(PERMEABILITY_YY) || rProp[PERMEABILITY_YY] < 0.0) + if (!r_prop.Has(PERMEABILITY_YY) || r_prop[PERMEABILITY_YY] < 0.0) KRATOS_ERROR << "PERMEABILITY_YY has Key zero, is not defined or " "has an invalid value at element" << this->Id() << std::endl; - if (!rProp.Has(PERMEABILITY_XY) || rProp[PERMEABILITY_XY] < 0.0) + if (!r_prop.Has(PERMEABILITY_XY) || r_prop[PERMEABILITY_XY] < 0.0) KRATOS_ERROR << "PERMEABILITY_XY has Key zero, is not defined or " "has an invalid value at element" << this->Id() << std::endl; - if (rGeom.WorkingSpaceDimension() > 2) { - if (!rProp.Has(PERMEABILITY_ZZ) || rProp[PERMEABILITY_ZZ] < 0.0) + if (r_geom.WorkingSpaceDimension() > 2) { + if (!r_prop.Has(PERMEABILITY_ZZ) || r_prop[PERMEABILITY_ZZ] < 0.0) KRATOS_ERROR << "PERMEABILITY_ZZ has Key zero, is not defined " "or has an invalid value at element" << this->Id() << std::endl; - if (!rProp.Has(PERMEABILITY_YZ) || rProp[PERMEABILITY_YZ] < 0.0) + if (!r_prop.Has(PERMEABILITY_YZ) || r_prop[PERMEABILITY_YZ] < 0.0) KRATOS_ERROR << "PERMEABILITY_YZ has Key zero, is not defined " "or has an invalid value at element" << this->Id() << std::endl; - if (!rProp.Has(PERMEABILITY_ZX) || rProp[PERMEABILITY_ZX] < 0.0) + if (!r_prop.Has(PERMEABILITY_ZX) || r_prop[PERMEABILITY_ZX] < 0.0) KRATOS_ERROR << "PERMEABILITY_ZX has Key zero, is not defined " "or has an invalid value at element" << this->Id() << std::endl; @@ -106,62 +104,59 @@ int SmallStrainUPwDiffOrderElement::Check(const ProcessInfo& rCurrentProcessInfo } // verify that the dofs exist - for (unsigned int i = 0; i < rGeom.size(); ++i) { - if (!rGeom[i].SolutionStepsDataHas(DISPLACEMENT)) - KRATOS_ERROR << "missing variable DISPLACEMENT on node " << rGeom[i].Id() << std::endl; - - if (!rGeom[i].HasDofFor(DISPLACEMENT_X) || !rGeom[i].HasDofFor(DISPLACEMENT_Y) || - !rGeom[i].HasDofFor(DISPLACEMENT_Z)) - KRATOS_ERROR << "missing one of the dofs for the variable " - "DISPLACEMENT on node " - << rGeom[i].Id() << std::endl; - - if (!rGeom[i].SolutionStepsDataHas(WATER_PRESSURE)) - KRATOS_ERROR << "missing variable WATER_PRESSURE on node " << rGeom[i].Id() << std::endl; - - if (!rGeom[i].HasDofFor(WATER_PRESSURE)) - KRATOS_ERROR << "missing the dof for the variable WATER_PRESSURE " - "on node " - << rGeom[i].Id() << std::endl; + for (const auto& r_node : r_geom) { + if (!r_node.SolutionStepsDataHas(DISPLACEMENT)) + KRATOS_ERROR << "missing variable DISPLACEMENT on node " << r_node.Id() << std::endl; + + if (!r_node.HasDofFor(DISPLACEMENT_X) || !r_node.HasDofFor(DISPLACEMENT_Y) || + !r_node.HasDofFor(DISPLACEMENT_Z)) + KRATOS_ERROR << "missing one of the dofs for the variable DISPLACEMENT on node " + << r_node.Id() << std::endl; + + if (!r_node.SolutionStepsDataHas(WATER_PRESSURE)) + KRATOS_ERROR << "missing variable WATER_PRESSURE on node " << r_node.Id() << std::endl; + + if (!r_node.HasDofFor(WATER_PRESSURE)) + KRATOS_ERROR << "missing the dof for the variable WATER_PRESSURE on node " + << r_node.Id() << std::endl; } // Verify that the constitutive law exists - KRATOS_ERROR_IF_NOT(rProp.Has(CONSTITUTIVE_LAW)) - << "Constitutive law not provided for property " << rProp.Id() << std::endl; + KRATOS_ERROR_IF_NOT(r_prop.Has(CONSTITUTIVE_LAW)) + << "Constitutive law not provided for property " << r_prop.Id() << std::endl; // verify compatibility with the constitutive law ConstitutiveLaw::Features LawFeatures; - rProp.GetValue(CONSTITUTIVE_LAW)->GetLawFeatures(LawFeatures); + r_prop.GetValue(CONSTITUTIVE_LAW)->GetLawFeatures(LawFeatures); - bool correct_strain_measure = false; - for (unsigned int i = 0; i < LawFeatures.mStrainMeasures.size(); ++i) { - if (LawFeatures.mStrainMeasures[i] == ConstitutiveLaw::StrainMeasure_Infinitesimal) - correct_strain_measure = true; - } + KRATOS_ERROR_IF(std::find(LawFeatures.mStrainMeasures.cbegin(), LawFeatures.mStrainMeasures.cend(), + ConstitutiveLaw::StrainMeasure_Infinitesimal) == + LawFeatures.mStrainMeasures.cend()) + << "In element " << this->Id() + << " the constitutive law is not compatible with the element " + "type StrainMeasure_Infinitesimal." + << std::endl; - if (!correct_strain_measure) - KRATOS_ERROR << "constitutive law is not compatible with the element " - "type StrainMeasure_Infinitesimal " - << this->Id() << std::endl; - - rProp.GetValue(CONSTITUTIVE_LAW)->Check(rProp, rGeom, rCurrentProcessInfo); + r_prop.GetValue(CONSTITUTIVE_LAW)->Check(r_prop, r_geom, rCurrentProcessInfo); // Verify that the constitutive law has the correct dimension const SizeType strainSize = this->GetProperties().GetValue(CONSTITUTIVE_LAW)->GetStrainSize(); - if (rGeom.WorkingSpaceDimension() > 2) { + if (r_geom.WorkingSpaceDimension() > 2) { KRATOS_ERROR_IF_NOT(strainSize == VOIGT_SIZE_3D) - << "Wrong constitutive law used. This is a 3D element! expected " - "strain size is " - << VOIGT_SIZE_3D << " But received: " << strainSize << " in element id: " << this->Id() - << std::endl; + << "Wrong constitutive law used. This is a 3D element! expected strain size is " << VOIGT_SIZE_3D + << " But received: " << strainSize << " in element id: " << this->Id() << std::endl; } else { KRATOS_ERROR_IF_NOT(strainSize == VOIGT_SIZE_2D_PLANE_STRAIN) - << "Wrong constitutive law used. This is a 2D element! expected " - "strain size is " + << "Wrong constitutive law used. This is a 2D element! expected strain size is " << VOIGT_SIZE_2D_PLANE_STRAIN << " But received: " << strainSize << " in element id: " << this->Id() << std::endl; } + // Check retention law + if (!mRetentionLawVector.empty()) { + return mRetentionLawVector[0]->Check(r_prop, rCurrentProcessInfo); + } + return 0; KRATOS_CATCH("") @@ -180,8 +175,6 @@ void SmallStrainUPwDiffOrderElement::InitializeSolutionStep(const ProcessInfo& r ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - RetentionLaw::Parameters RetentionParameters(GetProperties()); - const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); const auto determinants_of_deformation_gradients = @@ -204,8 +197,6 @@ void SmallStrainUPwDiffOrderElement::InitializeSolutionStep(const ProcessInfo& r noalias(Variables.StressVector) = mStressVector[GPoint]; ConstitutiveParameters.SetStressVector(Variables.StressVector); mConstitutiveLawVector[GPoint]->InitializeMaterialResponseCauchy(ConstitutiveParameters); - - mRetentionLawVector[GPoint]->InitializeSolutionStep(RetentionParameters); } KRATOS_CATCH("") @@ -254,8 +245,6 @@ void SmallStrainUPwDiffOrderElement::FinalizeSolutionStep(const ProcessInfo& rCu ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - RetentionLaw::Parameters RetentionParameters(GetProperties()); - const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); const auto determinants_of_deformation_gradients = @@ -281,8 +270,6 @@ void SmallStrainUPwDiffOrderElement::FinalizeSolutionStep(const ProcessInfo& rCu mConstitutiveLawVector[GPoint]->FinalizeMaterialResponseCauchy(ConstitutiveParameters); mStateVariablesFinalized[GPoint] = mConstitutiveLawVector[GPoint]->GetValue(STATE_VARIABLES, mStateVariablesFinalized[GPoint]); - - mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters); } // Assign pressure values to the intermediate nodes for post-processing @@ -291,182 +278,151 @@ void SmallStrainUPwDiffOrderElement::FinalizeSolutionStep(const ProcessInfo& rCu KRATOS_CATCH("") } +Vector SmallStrainUPwDiffOrderElement::GetPressures(const size_t n_nodes) const +{ + const auto& r_geom = GetGeometry(); + Vector pressure(n_nodes); + std::transform(r_geom.begin(), r_geom.begin() + n_nodes, pressure.begin(), + [](const auto& node) { return node.FastGetSolutionStepValue(WATER_PRESSURE); }); + return pressure; +} + void SmallStrainUPwDiffOrderElement::AssignPressureToIntermediateNodes() { // Assign pressure values to the intermediate nodes for post-processing KRATOS_TRY - GeometryType& rGeom = GetGeometry(); - const SizeType NumUNodes = rGeom.PointsNumber(); - const SizeType NumDim = rGeom.WorkingSpaceDimension(); + GeometryType& r_geom = GetGeometry(); + const SizeType num_u_nodes = r_geom.PointsNumber(); + const SizeType n_dim = r_geom.WorkingSpaceDimension(); - switch (NumUNodes) { + switch (num_u_nodes) { case 6: // 2D T6P3 { - const double p0 = rGeom[0].FastGetSolutionStepValue(WATER_PRESSURE); - const double p1 = rGeom[1].FastGetSolutionStepValue(WATER_PRESSURE); - const double p2 = rGeom[2].FastGetSolutionStepValue(WATER_PRESSURE); - ThreadSafeNodeWrite(rGeom[3], WATER_PRESSURE, 0.5 * (p0 + p1)); - ThreadSafeNodeWrite(rGeom[4], WATER_PRESSURE, 0.5 * (p1 + p2)); - ThreadSafeNodeWrite(rGeom[5], WATER_PRESSURE, 0.5 * (p2 + p0)); + const Vector p = GetPressures(3); + ThreadSafeNodeWrite(r_geom[3], WATER_PRESSURE, 0.5 * (p[0] + p[1])); + ThreadSafeNodeWrite(r_geom[4], WATER_PRESSURE, 0.5 * (p[1] + p[2])); + ThreadSafeNodeWrite(r_geom[5], WATER_PRESSURE, 0.5 * (p[2] + p[0])); break; } case 8: // 2D Q8P4 { - const double p0 = rGeom[0].FastGetSolutionStepValue(WATER_PRESSURE); - const double p1 = rGeom[1].FastGetSolutionStepValue(WATER_PRESSURE); - const double p2 = rGeom[2].FastGetSolutionStepValue(WATER_PRESSURE); - const double p3 = rGeom[3].FastGetSolutionStepValue(WATER_PRESSURE); - ThreadSafeNodeWrite(rGeom[4], WATER_PRESSURE, 0.5 * (p0 + p1)); - ThreadSafeNodeWrite(rGeom[5], WATER_PRESSURE, 0.5 * (p1 + p2)); - ThreadSafeNodeWrite(rGeom[6], WATER_PRESSURE, 0.5 * (p2 + p3)); - ThreadSafeNodeWrite(rGeom[7], WATER_PRESSURE, 0.5 * (p3 + p0)); + const Vector p = GetPressures(4); + ThreadSafeNodeWrite(r_geom[4], WATER_PRESSURE, 0.5 * (p[0] + p[1])); + ThreadSafeNodeWrite(r_geom[5], WATER_PRESSURE, 0.5 * (p[1] + p[2])); + ThreadSafeNodeWrite(r_geom[6], WATER_PRESSURE, 0.5 * (p[2] + p[3])); + ThreadSafeNodeWrite(r_geom[7], WATER_PRESSURE, 0.5 * (p[3] + p[0])); break; } case 9: // 2D Q9P4 { - const double p0 = rGeom[0].FastGetSolutionStepValue(WATER_PRESSURE); - const double p1 = rGeom[1].FastGetSolutionStepValue(WATER_PRESSURE); - const double p2 = rGeom[2].FastGetSolutionStepValue(WATER_PRESSURE); - const double p3 = rGeom[3].FastGetSolutionStepValue(WATER_PRESSURE); - ThreadSafeNodeWrite(rGeom[4], WATER_PRESSURE, 0.5 * (p0 + p1)); - ThreadSafeNodeWrite(rGeom[5], WATER_PRESSURE, 0.5 * (p1 + p2)); - ThreadSafeNodeWrite(rGeom[6], WATER_PRESSURE, 0.5 * (p2 + p3)); - ThreadSafeNodeWrite(rGeom[7], WATER_PRESSURE, 0.5 * (p3 + p0)); - ThreadSafeNodeWrite(rGeom[8], WATER_PRESSURE, 0.25 * (p0 + p1 + p2 + p3)); + const Vector p = GetPressures(4); + ThreadSafeNodeWrite(r_geom[4], WATER_PRESSURE, 0.5 * (p[0] + p[1])); + ThreadSafeNodeWrite(r_geom[5], WATER_PRESSURE, 0.5 * (p[1] + p[2])); + ThreadSafeNodeWrite(r_geom[6], WATER_PRESSURE, 0.5 * (p[2] + p[3])); + ThreadSafeNodeWrite(r_geom[7], WATER_PRESSURE, 0.5 * (p[3] + p[0])); + ThreadSafeNodeWrite(r_geom[8], WATER_PRESSURE, 0.25 * (p[0] + p[1] + p[2] + p[3])); break; } case 10: // 3D T10P4 //2D T10P6 { - if (NumDim == 3) { - const double p0 = rGeom[0].FastGetSolutionStepValue(WATER_PRESSURE); - const double p1 = rGeom[1].FastGetSolutionStepValue(WATER_PRESSURE); - const double p2 = rGeom[2].FastGetSolutionStepValue(WATER_PRESSURE); - const double p3 = rGeom[3].FastGetSolutionStepValue(WATER_PRESSURE); - ThreadSafeNodeWrite(rGeom[4], WATER_PRESSURE, 0.5 * (p0 + p1)); - ThreadSafeNodeWrite(rGeom[5], WATER_PRESSURE, 0.5 * (p1 + p2)); - ThreadSafeNodeWrite(rGeom[6], WATER_PRESSURE, 0.5 * (p2 + p0)); - ThreadSafeNodeWrite(rGeom[7], WATER_PRESSURE, 0.5 * (p0 + p3)); - ThreadSafeNodeWrite(rGeom[8], WATER_PRESSURE, 0.5 * (p1 + p3)); - ThreadSafeNodeWrite(rGeom[9], WATER_PRESSURE, 0.5 * (p2 + p3)); - } else if (NumDim == 2) { + if (n_dim == 3) { + const Vector p = GetPressures(4); + ThreadSafeNodeWrite(r_geom[4], WATER_PRESSURE, 0.5 * (p[0] + p[1])); + ThreadSafeNodeWrite(r_geom[5], WATER_PRESSURE, 0.5 * (p[1] + p[2])); + ThreadSafeNodeWrite(r_geom[6], WATER_PRESSURE, 0.5 * (p[2] + p[0])); + ThreadSafeNodeWrite(r_geom[7], WATER_PRESSURE, 0.5 * (p[0] + p[3])); + ThreadSafeNodeWrite(r_geom[8], WATER_PRESSURE, 0.5 * (p[1] + p[3])); + ThreadSafeNodeWrite(r_geom[9], WATER_PRESSURE, 0.5 * (p[2] + p[3])); + } else if (n_dim == 2) { constexpr double c1 = 1.0 / 9.0; - const double p0 = rGeom[0].FastGetSolutionStepValue(WATER_PRESSURE); - const double p1 = rGeom[1].FastGetSolutionStepValue(WATER_PRESSURE); - const double p2 = rGeom[2].FastGetSolutionStepValue(WATER_PRESSURE); - const double p3 = rGeom[3].FastGetSolutionStepValue(WATER_PRESSURE); - const double p4 = rGeom[4].FastGetSolutionStepValue(WATER_PRESSURE); - const double p5 = rGeom[5].FastGetSolutionStepValue(WATER_PRESSURE); - ThreadSafeNodeWrite(rGeom[0], WATER_PRESSURE, p0); - ThreadSafeNodeWrite(rGeom[1], WATER_PRESSURE, p1); - ThreadSafeNodeWrite(rGeom[2], WATER_PRESSURE, p2); - ThreadSafeNodeWrite(rGeom[3], WATER_PRESSURE, (2.0 * p0 - p1 + 8.0 * p3) * c1); - ThreadSafeNodeWrite(rGeom[4], WATER_PRESSURE, (2.0 * p1 - p0 + 8.0 * p3) * c1); - ThreadSafeNodeWrite(rGeom[5], WATER_PRESSURE, (2.0 * p1 - p2 + 8.0 * p4) * c1); - ThreadSafeNodeWrite(rGeom[6], WATER_PRESSURE, (2.0 * p2 - p1 + 8.0 * p4) * c1); - ThreadSafeNodeWrite(rGeom[7], WATER_PRESSURE, (2.0 * p2 - p0 + 8.0 * p5) * c1); - ThreadSafeNodeWrite(rGeom[8], WATER_PRESSURE, (2.0 * p0 - p2 + 8.0 * p5) * c1); - ThreadSafeNodeWrite(rGeom[9], WATER_PRESSURE, (4.0 * (p3 + p4 + p5) - (p0 + p1 + p2)) * c1); + const Vector p = GetPressures(6); + ThreadSafeNodeWrite(r_geom[3], WATER_PRESSURE, (2.0 * p[0] - p[1] + 8.0 * p[3]) * c1); + ThreadSafeNodeWrite(r_geom[4], WATER_PRESSURE, (2.0 * p[1] - p[0] + 8.0 * p[3]) * c1); + ThreadSafeNodeWrite(r_geom[5], WATER_PRESSURE, (2.0 * p[1] - p[2] + 8.0 * p[4]) * c1); + ThreadSafeNodeWrite(r_geom[6], WATER_PRESSURE, (2.0 * p[2] - p[1] + 8.0 * p[4]) * c1); + ThreadSafeNodeWrite(r_geom[7], WATER_PRESSURE, (2.0 * p[2] - p[0] + 8.0 * p[5]) * c1); + ThreadSafeNodeWrite(r_geom[8], WATER_PRESSURE, (2.0 * p[0] - p[2] + 8.0 * p[5]) * c1); + ThreadSafeNodeWrite(r_geom[9], WATER_PRESSURE, + (4.0 * (p[3] + p[4] + p[5]) - (p[0] + p[1] + p[2])) * c1); } break; } case 15: // 2D T15P10 { constexpr double c1 = 0.0390625; - const double p0 = rGeom[0].FastGetSolutionStepValue(WATER_PRESSURE); - const double p1 = rGeom[1].FastGetSolutionStepValue(WATER_PRESSURE); - const double p2 = rGeom[2].FastGetSolutionStepValue(WATER_PRESSURE); - const double p3 = rGeom[3].FastGetSolutionStepValue(WATER_PRESSURE); - const double p4 = rGeom[4].FastGetSolutionStepValue(WATER_PRESSURE); - const double p5 = rGeom[5].FastGetSolutionStepValue(WATER_PRESSURE); - const double p6 = rGeom[6].FastGetSolutionStepValue(WATER_PRESSURE); - const double p7 = rGeom[7].FastGetSolutionStepValue(WATER_PRESSURE); - const double p8 = rGeom[8].FastGetSolutionStepValue(WATER_PRESSURE); - const double p9 = rGeom[9].FastGetSolutionStepValue(WATER_PRESSURE); - ThreadSafeNodeWrite(rGeom[0], WATER_PRESSURE, p0); - ThreadSafeNodeWrite(rGeom[1], WATER_PRESSURE, p1); - ThreadSafeNodeWrite(rGeom[2], WATER_PRESSURE, p2); - ThreadSafeNodeWrite(rGeom[3], WATER_PRESSURE, (3.0 * p0 + p1 + 27.0 * p3 - 5.4 * p4) * c1); - ThreadSafeNodeWrite(rGeom[4], WATER_PRESSURE, (14.4 * (p3 + p4) - 1.6 * (p0 + p1)) * c1); - ThreadSafeNodeWrite(rGeom[5], WATER_PRESSURE, (3.0 * p1 + p0 + 27.0 * p4 - 5.4 * p3) * c1); - ThreadSafeNodeWrite(rGeom[6], WATER_PRESSURE, (3.0 * p1 + p2 + 27.0 * p5 - 5.4 * p6) * c1); - ThreadSafeNodeWrite(rGeom[7], WATER_PRESSURE, (14.4 * (p5 + p6) - 1.6 * (p1 + p2)) * c1); - ThreadSafeNodeWrite(rGeom[8], WATER_PRESSURE, (3.0 * p2 + p1 + 27.0 * p6 - 5.4 * p5) * c1); - ThreadSafeNodeWrite(rGeom[9], WATER_PRESSURE, (3.0 * p2 + p0 + 27.0 * p7 - 5.4 * p8) * c1); - ThreadSafeNodeWrite(rGeom[10], WATER_PRESSURE, (14.4 * (p7 + p8) - 1.6 * (p0 + p2)) * c1); - ThreadSafeNodeWrite(rGeom[11], WATER_PRESSURE, (3.0 * p0 + p2 + 27.0 * p8 - 5.4 * p7) * c1); - ThreadSafeNodeWrite( - rGeom[12], WATER_PRESSURE, - (p1 + p2 + 7.2 * (p3 + p8) - 3.6 * (p4 + p7) - 1.8 * (p5 + p6) + 21.6 * p9 - 1.6 * p0) * c1); - ThreadSafeNodeWrite( - rGeom[13], WATER_PRESSURE, - (p0 + p2 + 7.2 * (p4 + p5) - 3.6 * (p3 + p6) - 1.8 * (p7 + p8) + 21.6 * p9 - 1.6 * p1) * c1); - ThreadSafeNodeWrite( - rGeom[14], WATER_PRESSURE, - (p0 + p1 + 7.2 * (p6 + p7) - 3.6 * (p5 + p8) - 1.8 * (p3 + p4) + 21.6 * p9 - 1.6 * p2) * c1); + const Vector p = GetPressures(10); + ThreadSafeNodeWrite(r_geom[3], WATER_PRESSURE, (3.0 * p[0] + p[1] + 27.0 * p[3] - 5.4 * p[4]) * c1); + ThreadSafeNodeWrite(r_geom[4], WATER_PRESSURE, (14.4 * (p[3] + p[4]) - 1.6 * (p[0] + p[1])) * c1); + ThreadSafeNodeWrite(r_geom[5], WATER_PRESSURE, (3.0 * p[1] + p[0] + 27.0 * p[4] - 5.4 * p[3]) * c1); + ThreadSafeNodeWrite(r_geom[6], WATER_PRESSURE, (3.0 * p[1] + p[2] + 27.0 * p[5] - 5.4 * p[6]) * c1); + ThreadSafeNodeWrite(r_geom[7], WATER_PRESSURE, (14.4 * (p[5] + p[6]) - 1.6 * (p[1] + p[2])) * c1); + ThreadSafeNodeWrite(r_geom[8], WATER_PRESSURE, (3.0 * p[2] + p[1] + 27.0 * p[6] - 5.4 * p[5]) * c1); + ThreadSafeNodeWrite(r_geom[9], WATER_PRESSURE, (3.0 * p[2] + p[0] + 27.0 * p[7] - 5.4 * p[8]) * c1); + ThreadSafeNodeWrite(r_geom[10], WATER_PRESSURE, (14.4 * (p[7] + p[8]) - 1.6 * (p[0] + p[2])) * c1); + ThreadSafeNodeWrite(r_geom[11], WATER_PRESSURE, (3.0 * p[0] + p[2] + 27.0 * p[8] - 5.4 * p[7]) * c1); + ThreadSafeNodeWrite(r_geom[12], WATER_PRESSURE, + (p[1] + p[2] + 7.2 * (p[3] + p[8]) - 3.6 * (p[4] + p[7]) - + 1.8 * (p[5] + p[6]) + 21.6 * p[9] - 1.6 * p[0]) * + c1); + ThreadSafeNodeWrite(r_geom[13], WATER_PRESSURE, + (p[0] + p[2] + 7.2 * (p[4] + p[5]) - 3.6 * (p[3] + p[6]) - + 1.8 * (p[7] + p[8]) + 21.6 * p[9] - 1.6 * p[1]) * + c1); + ThreadSafeNodeWrite(r_geom[14], WATER_PRESSURE, + (p[0] + p[1] + 7.2 * (p[6] + p[7]) - 3.6 * (p[5] + p[8]) - + 1.8 * (p[3] + p[4]) + 21.6 * p[9] - 1.6 * p[2]) * + c1); break; } case 20: // 3D H20P8 { - const double p0 = rGeom[0].FastGetSolutionStepValue(WATER_PRESSURE); - const double p1 = rGeom[1].FastGetSolutionStepValue(WATER_PRESSURE); - const double p2 = rGeom[2].FastGetSolutionStepValue(WATER_PRESSURE); - const double p3 = rGeom[3].FastGetSolutionStepValue(WATER_PRESSURE); - const double p4 = rGeom[4].FastGetSolutionStepValue(WATER_PRESSURE); - const double p5 = rGeom[5].FastGetSolutionStepValue(WATER_PRESSURE); - const double p6 = rGeom[6].FastGetSolutionStepValue(WATER_PRESSURE); - const double p7 = rGeom[7].FastGetSolutionStepValue(WATER_PRESSURE); + const Vector p = GetPressures(8); // edges -- bottom - ThreadSafeNodeWrite(rGeom[8], WATER_PRESSURE, 0.5 * (p0 + p1)); - ThreadSafeNodeWrite(rGeom[9], WATER_PRESSURE, 0.5 * (p1 + p2)); - ThreadSafeNodeWrite(rGeom[10], WATER_PRESSURE, 0.5 * (p2 + p3)); - ThreadSafeNodeWrite(rGeom[11], WATER_PRESSURE, 0.5 * (p3 + p0)); + ThreadSafeNodeWrite(r_geom[8], WATER_PRESSURE, 0.5 * (p[0] + p[1])); + ThreadSafeNodeWrite(r_geom[9], WATER_PRESSURE, 0.5 * (p[1] + p[2])); + ThreadSafeNodeWrite(r_geom[10], WATER_PRESSURE, 0.5 * (p[2] + p[3])); + ThreadSafeNodeWrite(r_geom[11], WATER_PRESSURE, 0.5 * (p[3] + p[0])); // edges -- middle - ThreadSafeNodeWrite(rGeom[12], WATER_PRESSURE, 0.5 * (p4 + p0)); - ThreadSafeNodeWrite(rGeom[13], WATER_PRESSURE, 0.5 * (p5 + p1)); - ThreadSafeNodeWrite(rGeom[14], WATER_PRESSURE, 0.5 * (p6 + p2)); - ThreadSafeNodeWrite(rGeom[15], WATER_PRESSURE, 0.5 * (p7 + p3)); + ThreadSafeNodeWrite(r_geom[12], WATER_PRESSURE, 0.5 * (p[4] + p[0])); + ThreadSafeNodeWrite(r_geom[13], WATER_PRESSURE, 0.5 * (p[5] + p[1])); + ThreadSafeNodeWrite(r_geom[14], WATER_PRESSURE, 0.5 * (p[6] + p[2])); + ThreadSafeNodeWrite(r_geom[15], WATER_PRESSURE, 0.5 * (p[7] + p[3])); // edges -- top - ThreadSafeNodeWrite(rGeom[16], WATER_PRESSURE, 0.5 * (p4 + p5)); - ThreadSafeNodeWrite(rGeom[17], WATER_PRESSURE, 0.5 * (p5 + p6)); - ThreadSafeNodeWrite(rGeom[18], WATER_PRESSURE, 0.5 * (p6 + p7)); - ThreadSafeNodeWrite(rGeom[19], WATER_PRESSURE, 0.5 * (p7 + p4)); + ThreadSafeNodeWrite(r_geom[16], WATER_PRESSURE, 0.5 * (p[4] + p[5])); + ThreadSafeNodeWrite(r_geom[17], WATER_PRESSURE, 0.5 * (p[5] + p[6])); + ThreadSafeNodeWrite(r_geom[18], WATER_PRESSURE, 0.5 * (p[6] + p[7])); + ThreadSafeNodeWrite(r_geom[19], WATER_PRESSURE, 0.5 * (p[7] + p[4])); break; } case 27: // 3D H27P8 { - const double p0 = rGeom[0].FastGetSolutionStepValue(WATER_PRESSURE); - const double p1 = rGeom[1].FastGetSolutionStepValue(WATER_PRESSURE); - const double p2 = rGeom[2].FastGetSolutionStepValue(WATER_PRESSURE); - const double p3 = rGeom[3].FastGetSolutionStepValue(WATER_PRESSURE); - const double p4 = rGeom[4].FastGetSolutionStepValue(WATER_PRESSURE); - const double p5 = rGeom[5].FastGetSolutionStepValue(WATER_PRESSURE); - const double p6 = rGeom[6].FastGetSolutionStepValue(WATER_PRESSURE); - const double p7 = rGeom[7].FastGetSolutionStepValue(WATER_PRESSURE); + const Vector p = GetPressures(8); // edges -- bottom - ThreadSafeNodeWrite(rGeom[8], WATER_PRESSURE, 0.5 * (p0 + p1)); - ThreadSafeNodeWrite(rGeom[9], WATER_PRESSURE, 0.5 * (p1 + p2)); - ThreadSafeNodeWrite(rGeom[10], WATER_PRESSURE, 0.5 * (p2 + p3)); - ThreadSafeNodeWrite(rGeom[11], WATER_PRESSURE, 0.5 * (p3 + p0)); + ThreadSafeNodeWrite(r_geom[8], WATER_PRESSURE, 0.5 * (p[0] + p[1])); + ThreadSafeNodeWrite(r_geom[9], WATER_PRESSURE, 0.5 * (p[1] + p[2])); + ThreadSafeNodeWrite(r_geom[10], WATER_PRESSURE, 0.5 * (p[2] + p[3])); + ThreadSafeNodeWrite(r_geom[11], WATER_PRESSURE, 0.5 * (p[3] + p[0])); // edges -- middle - ThreadSafeNodeWrite(rGeom[12], WATER_PRESSURE, 0.5 * (p4 + p0)); - ThreadSafeNodeWrite(rGeom[13], WATER_PRESSURE, 0.5 * (p5 + p1)); - ThreadSafeNodeWrite(rGeom[14], WATER_PRESSURE, 0.5 * (p6 + p2)); - ThreadSafeNodeWrite(rGeom[15], WATER_PRESSURE, 0.5 * (p7 + p3)); + ThreadSafeNodeWrite(r_geom[12], WATER_PRESSURE, 0.5 * (p[4] + p[0])); + ThreadSafeNodeWrite(r_geom[13], WATER_PRESSURE, 0.5 * (p[5] + p[1])); + ThreadSafeNodeWrite(r_geom[14], WATER_PRESSURE, 0.5 * (p[6] + p[2])); + ThreadSafeNodeWrite(r_geom[15], WATER_PRESSURE, 0.5 * (p[7] + p[3])); // edges -- top - ThreadSafeNodeWrite(rGeom[16], WATER_PRESSURE, 0.5 * (p4 + p5)); - ThreadSafeNodeWrite(rGeom[17], WATER_PRESSURE, 0.5 * (p5 + p6)); - ThreadSafeNodeWrite(rGeom[18], WATER_PRESSURE, 0.5 * (p6 + p7)); - ThreadSafeNodeWrite(rGeom[19], WATER_PRESSURE, 0.5 * (p7 + p0)); + ThreadSafeNodeWrite(r_geom[16], WATER_PRESSURE, 0.5 * (p[4] + p[5])); + ThreadSafeNodeWrite(r_geom[17], WATER_PRESSURE, 0.5 * (p[5] + p[6])); + ThreadSafeNodeWrite(r_geom[18], WATER_PRESSURE, 0.5 * (p[6] + p[7])); + ThreadSafeNodeWrite(r_geom[19], WATER_PRESSURE, 0.5 * (p[7] + p[0])); // face centers - ThreadSafeNodeWrite(rGeom[20], WATER_PRESSURE, 0.25 * (p0 + p1 + p2 + p3)); - ThreadSafeNodeWrite(rGeom[21], WATER_PRESSURE, 0.25 * (p0 + p1 + p4 + p5)); - ThreadSafeNodeWrite(rGeom[22], WATER_PRESSURE, 0.25 * (p1 + p2 + p5 + p6)); - ThreadSafeNodeWrite(rGeom[23], WATER_PRESSURE, 0.25 * (p2 + p3 + p6 + p7)); - ThreadSafeNodeWrite(rGeom[24], WATER_PRESSURE, 0.25 * (p3 + p0 + p7 + p4)); - ThreadSafeNodeWrite(rGeom[25], WATER_PRESSURE, 0.25 * (p4 + p5 + p6 + p7)); + ThreadSafeNodeWrite(r_geom[20], WATER_PRESSURE, 0.25 * (p[0] + p[1] + p[2] + p[3])); + ThreadSafeNodeWrite(r_geom[21], WATER_PRESSURE, 0.25 * (p[0] + p[1] + p[4] + p[5])); + ThreadSafeNodeWrite(r_geom[22], WATER_PRESSURE, 0.25 * (p[1] + p[2] + p[5] + p[6])); + ThreadSafeNodeWrite(r_geom[23], WATER_PRESSURE, 0.25 * (p[2] + p[3] + p[6] + p[7])); + ThreadSafeNodeWrite(r_geom[24], WATER_PRESSURE, 0.25 * (p[3] + p[0] + p[7] + p[4])); + ThreadSafeNodeWrite(r_geom[25], WATER_PRESSURE, 0.25 * (p[4] + p[5] + p[6] + p[7])); // element center - ThreadSafeNodeWrite(rGeom[26], WATER_PRESSURE, 0.125 * (p0 + p1 + p2 + p3 + p4 + p5 + p6 + p7)); + ThreadSafeNodeWrite(r_geom[26], WATER_PRESSURE, + 0.125 * (p[0] + p[1] + p[2] + p[3] + p[4] + p[5] + p[6] + p[7])); break; } default: @@ -526,8 +482,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable { KRATOS_TRY - const GeometryType& rGeom = GetGeometry(); - const auto number_of_integration_points = rGeom.IntegrationPointsNumber(this->GetIntegrationMethod()); + const GeometryType& r_geom = GetGeometry(); + const auto number_of_integration_points = r_geom.IntegrationPointsNumber(this->GetIntegrationMethod()); rOutput.resize(number_of_integration_points); @@ -589,55 +545,46 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure( Variables.Np, Variables.PressureVector)); - if (rVariable == DEGREE_OF_SATURATION) - rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateSaturation(RetentionParameters); - else if (rVariable == EFFECTIVE_SATURATION) - rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(RetentionParameters); - else if (rVariable == BISHOP_COEFFICIENT) - rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); - else if (rVariable == DERIVATIVE_OF_SATURATION) - rOutput[GPoint] = - mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(RetentionParameters); - else if (rVariable == RELATIVE_PERMEABILITY) - rOutput[GPoint] = - mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); + rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateValue( + RetentionParameters, rVariable, rOutput[GPoint]); } } else if (rVariable == HYDRAULIC_HEAD) { - const double NumericalLimit = std::numeric_limits::epsilon(); - const PropertiesType& rProp = this->GetProperties(); + constexpr auto numerical_limit = std::numeric_limits::epsilon(); + const PropertiesType& r_prop = this->GetProperties(); // Defining the shape functions, the Jacobian and the shape functions local gradients Containers - const Matrix& NContainer = rGeom.ShapeFunctionsValues(this->GetIntegrationMethod()); - const SizeType NumUNodes = rGeom.PointsNumber(); + const Matrix& n_container = r_geom.ShapeFunctionsValues(this->GetIntegrationMethod()); + const SizeType num_u_nodes = r_geom.PointsNumber(); // Defining necessary variables - Vector NodalHydraulicHead = ZeroVector(NumUNodes); - for (unsigned int node = 0; node < NumUNodes; ++node) { + Vector nodal_hydraulic_head = ZeroVector(num_u_nodes); + for (unsigned int node = 0; node < num_u_nodes; ++node) { Vector NodeVolumeAcceleration(3); - noalias(NodeVolumeAcceleration) = rGeom[node].FastGetSolutionStepValue(VOLUME_ACCELERATION, 0); + noalias(NodeVolumeAcceleration) = r_geom[node].FastGetSolutionStepValue(VOLUME_ACCELERATION, 0); const double g = norm_2(NodeVolumeAcceleration); - if (g > NumericalLimit) { - const double FluidWeight = g * rProp[DENSITY_WATER]; - - Vector NodeCoordinates(3); - noalias(NodeCoordinates) = rGeom[node].Coordinates(); - Vector NodeVolumeAccelerationUnitVector(3); - noalias(NodeVolumeAccelerationUnitVector) = NodeVolumeAcceleration / g; - - const double WaterPressure = rGeom[node].FastGetSolutionStepValue(WATER_PRESSURE); - NodalHydraulicHead[node] = -inner_prod(NodeCoordinates, NodeVolumeAccelerationUnitVector) - - PORE_PRESSURE_SIGN_FACTOR * WaterPressure / FluidWeight; + if (g > numerical_limit) { + const auto fluid_weight = g * r_prop[DENSITY_WATER]; + + Vector node_coordinates(3); + noalias(node_coordinates) = r_geom[node].Coordinates(); + Vector node_volume_acceleration_unit_vector(3); + noalias(node_volume_acceleration_unit_vector) = NodeVolumeAcceleration / g; + + const auto water_pressure = r_geom[node].FastGetSolutionStepValue(WATER_PRESSURE); + nodal_hydraulic_head[node] = + -inner_prod(node_coordinates, node_volume_acceleration_unit_vector) - + PORE_PRESSURE_SIGN_FACTOR * water_pressure / fluid_weight; } else { - NodalHydraulicHead[node] = 0.0; + nodal_hydraulic_head[node] = 0.0; } } - for (unsigned int GPoint = 0; GPoint < number_of_integration_points; GPoint++) { - double HydraulicHead = 0.0; - for (unsigned int node = 0; node < NumUNodes; ++node) - HydraulicHead += NContainer(GPoint, node) * NodalHydraulicHead[node]; + for (unsigned int g_point = 0; g_point < number_of_integration_points; g_point++) { + double hydraulic_head = 0.0; + for (unsigned int node = 0; node < num_u_nodes; ++node) + hydraulic_head += n_container(g_point, node) * nodal_hydraulic_head[node]; - rOutput[GPoint] = HydraulicHead; + rOutput[g_point] = hydraulic_head; } } else { for (unsigned int i = 0; i < number_of_integration_points; ++i) { @@ -677,10 +624,10 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable std::multiplies<>{}); // Loop over integration points - for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { + for (unsigned int g_point = 0; g_point < mConstitutiveLawVector.size(); ++g_point) { // compute element kinematics (Np, gradNpT, |J|, B, strains) - this->CalculateKinematics(Variables, GPoint); - Variables.B = b_matrices[GPoint]; + this->CalculateKinematics(Variables, g_point); + Variables.B = b_matrices[g_point]; // Compute FluidFlux vector q [m/s] const SizeType Dim = r_geometry.WorkingSpaceDimension(); @@ -693,28 +640,27 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable BodyAcceleration[idim] += Variables.Nu[i] * Variables.BodyAcceleration[Index++]; } - const auto relative_permeability = relative_permeability_values[GPoint]; + const auto relative_permeability = relative_permeability_values[g_point]; // Compute strain, need to update porosity - Variables.F = deformation_gradients[GPoint]; - Variables.StrainVector = strain_vectors[GPoint]; + Variables.F = deformation_gradients[g_point]; + Variables.StrainVector = strain_vectors[g_point]; Vector GradPressureTerm(Dim); noalias(GradPressureTerm) = prod(trans(Variables.DNp_DX), Variables.PressureVector); noalias(GradPressureTerm) += PORE_PRESSURE_SIGN_FACTOR * GetProperties()[DENSITY_WATER] * BodyAcceleration; - Vector AuxFluidFlux = ZeroVector(Dim); - AuxFluidFlux = PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * - relative_permeability * prod(Variables.IntrinsicPermeability, GradPressureTerm); + Vector aux_fluid_flux = ZeroVector(Dim); + aux_fluid_flux = PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * + relative_permeability * prod(Variables.IntrinsicPermeability, GradPressureTerm); - Vector FluidFlux = ZeroVector(3); + Vector fluid_flux = ZeroVector(3); for (unsigned int idim = 0; idim < Dim; ++idim) - FluidFlux[idim] = AuxFluidFlux[idim]; + fluid_flux[idim] = aux_fluid_flux[idim]; - if (rOutput[GPoint].size() != 3) rOutput[GPoint].resize(3, false); - - rOutput[GPoint] = FluidFlux; + rOutput[g_point].resize(3, false); + rOutput[g_point] = fluid_flux; } } @@ -727,9 +673,9 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable { KRATOS_TRY - const GeometryType& rGeom = GetGeometry(); + const GeometryType& r_geom = GetGeometry(); - rOutput.resize(rGeom.IntegrationPointsNumber(this->GetIntegrationMethod())); + rOutput.resize(r_geom.IntegrationPointsNumber(this->GetIntegrationMethod())); if (rVariable == CAUCHY_STRESS_VECTOR) { for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { @@ -742,7 +688,7 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom, GetProperties(), rCurrentProcessInfo); + ConstitutiveLaw::Parameters ConstitutiveParameters(r_geom, GetProperties(), rCurrentProcessInfo); ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR); ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); @@ -861,12 +807,12 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi { KRATOS_TRY - const PropertiesType& rProp = this->GetProperties(); - const GeometryType& rGeom = GetGeometry(); - const GeometryType::IntegrationPointsArrayType& IntegrationPoints = - rGeom.IntegrationPoints(this->GetIntegrationMethod()); + const PropertiesType& r_prop = this->GetProperties(); + const GeometryType& r_geom = GetGeometry(); + const GeometryType::IntegrationPointsArrayType& r_integration_points = + r_geom.IntegrationPoints(this->GetIntegrationMethod()); - ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom, rProp, rCurrentProcessInfo); + ConstitutiveLaw::Parameters ConstitutiveParameters(r_geom, r_prop, rCurrentProcessInfo); // Stiffness matrix is needed to calculate Biot coefficient ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR); @@ -876,17 +822,15 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - RetentionLaw::Parameters RetentionParameters(rProp); - const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto integration_coefficients = - this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJuContainer); + this->CalculateIntegrationCoefficients(r_integration_points, Variables.detJuContainer); - const auto det_Js_initial_configuration = - GeoEquationOfMotionUtilities::CalculateDetJsInitialConfiguration(rGeom, this->GetIntegrationMethod()); + const auto det_Js_initial_configuration = GeoEquationOfMotionUtilities::CalculateDetJsInitialConfiguration( + r_geom, this->GetIntegrationMethod()); const auto integration_coefficients_on_initial_configuration = - this->CalculateIntegrationCoefficients(IntegrationPoints, det_Js_initial_configuration); + this->CalculateIntegrationCoefficients(r_integration_points, det_Js_initial_configuration); const auto deformation_gradients = CalculateDeformationGradients(); auto strain_vectors = StressStrainUtilities::CalculateStrains( @@ -903,7 +847,7 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi const auto degrees_of_saturation = CalculateDegreesOfSaturation(fluid_pressures); const auto derivatives_of_saturation = CalculateDerivativesOfSaturation(fluid_pressures); const auto biot_moduli_inverse = GeoTransportEquationUtilities::CalculateInverseBiotModuli( - biot_coefficients, degrees_of_saturation, derivatives_of_saturation, rProp); + biot_coefficients, degrees_of_saturation, derivatives_of_saturation, r_prop); auto relative_permeability_values = CalculateRelativePermeabilityValues(fluid_pressures); const auto permeability_update_factors = GetOptionalPermeabilityUpdateFactors(strain_vectors); std::transform(permeability_update_factors.cbegin(), permeability_update_factors.cend(), @@ -912,7 +856,7 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi const auto bishop_coefficients = CalculateBishopCoefficients(fluid_pressures); - for (unsigned int GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) { + for (unsigned int GPoint = 0; GPoint < r_integration_points.size(); ++GPoint) { this->CalculateKinematics(Variables, GPoint); Variables.B = b_matrices[GPoint]; Variables.F = deformation_gradients[GPoint]; @@ -981,20 +925,20 @@ void SmallStrainUPwDiffOrderElement::CalculateMaterialStiffnessMatrix(MatrixType { KRATOS_TRY - const GeometryType& rGeom = GetGeometry(); + const GeometryType& r_geom = GetGeometry(); // Definition of variables ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); // Create constitutive law parameters: - ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom, GetProperties(), rCurrentProcessInfo); + ConstitutiveLaw::Parameters ConstitutiveParameters(r_geom, GetProperties(), rCurrentProcessInfo); ConstitutiveParameters.GetOptions().Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); ConstitutiveParameters.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR); // Loop over integration points - const GeometryType::IntegrationPointsArrayType& IntegrationPoints = - rGeom.IntegrationPoints(this->GetIntegrationMethod()); + const GeometryType::IntegrationPointsArrayType& r_integration_points = + r_geom.IntegrationPoints(this->GetIntegrationMethod()); const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); @@ -1006,7 +950,7 @@ void SmallStrainUPwDiffOrderElement::CalculateMaterialStiffnessMatrix(MatrixType Variables.NuContainer, Variables.DNu_DXContainer, strain_vectors, mStressVector, constitutive_matrices); const auto integration_coefficients = - this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJuContainer); + this->CalculateIntegrationCoefficients(r_integration_points, Variables.detJuContainer); const auto stiffness_matrix = GeoEquationOfMotionUtilities::CalculateStiffnessMatrix( b_matrices, constitutive_matrices, integration_coefficients); @@ -1021,31 +965,31 @@ void SmallStrainUPwDiffOrderElement::InitializeElementVariables(ElementVariables { KRATOS_TRY - const GeometryType& rGeom = GetGeometry(); - const SizeType NumUNodes = rGeom.PointsNumber(); - const SizeType NumPNodes = mpPressureGeometry->PointsNumber(); - const SizeType NumGPoints = rGeom.IntegrationPointsNumber(this->GetIntegrationMethod()); - const SizeType Dim = rGeom.WorkingSpaceDimension(); + const GeometryType& r_geom = GetGeometry(); + const SizeType num_u_nodes = r_geom.PointsNumber(); + const SizeType num_p_nodes = mpPressureGeometry->PointsNumber(); + const SizeType num_g_points = r_geom.IntegrationPointsNumber(this->GetIntegrationMethod()); + const SizeType n_dim = r_geom.WorkingSpaceDimension(); // Variables at all integration points - rVariables.NuContainer.resize(NumGPoints, NumUNodes, false); - rVariables.NuContainer = rGeom.ShapeFunctionsValues(this->GetIntegrationMethod()); + rVariables.NuContainer.resize(num_g_points, num_u_nodes, false); + rVariables.NuContainer = r_geom.ShapeFunctionsValues(this->GetIntegrationMethod()); - rVariables.NpContainer.resize(NumGPoints, NumPNodes, false); + rVariables.NpContainer.resize(num_g_points, num_p_nodes, false); rVariables.NpContainer = mpPressureGeometry->ShapeFunctionsValues(this->GetIntegrationMethod()); - rVariables.Nu.resize(NumUNodes, false); - rVariables.Np.resize(NumPNodes, false); + rVariables.Nu.resize(num_u_nodes, false); + rVariables.Np.resize(num_p_nodes, false); - rVariables.DNu_DXContainer.resize(NumGPoints, false); - for (SizeType i = 0; i < NumGPoints; ++i) - ((rVariables.DNu_DXContainer)[i]).resize(NumUNodes, Dim, false); - rVariables.DNu_DX.resize(NumUNodes, Dim, false); - rVariables.DNu_DXInitialConfiguration.resize(NumUNodes, Dim, false); - rVariables.detJuContainer.resize(NumGPoints, false); + rVariables.DNu_DXContainer.resize(num_g_points, false); + for (SizeType i = 0; i < num_g_points; ++i) + ((rVariables.DNu_DXContainer)[i]).resize(num_u_nodes, n_dim, false); + rVariables.DNu_DX.resize(num_u_nodes, n_dim, false); + rVariables.DNu_DXInitialConfiguration.resize(num_u_nodes, n_dim, false); + rVariables.detJuContainer.resize(num_g_points, false); try { - rGeom.ShapeFunctionsIntegrationPointsGradients( + r_geom.ShapeFunctionsIntegrationPointsGradients( rVariables.DNu_DXContainer, rVariables.detJuContainer, this->GetIntegrationMethod()); } catch (Kratos::Exception& e) { KRATOS_INFO("Original error message") << e.what() << std::endl; @@ -1059,11 +1003,11 @@ void SmallStrainUPwDiffOrderElement::InitializeElementVariables(ElementVariables << this->Id() << std::endl; } - (rVariables.DNp_DXContainer).resize(NumGPoints, false); - for (SizeType i = 0; i < NumGPoints; ++i) - ((rVariables.DNp_DXContainer)[i]).resize(NumPNodes, Dim, false); - (rVariables.DNp_DX).resize(NumPNodes, Dim, false); - Vector detJpContainer = ZeroVector(NumGPoints); + (rVariables.DNp_DXContainer).resize(num_g_points, false); + for (SizeType i = 0; i < num_g_points; ++i) + ((rVariables.DNp_DXContainer)[i]).resize(num_p_nodes, n_dim, false); + (rVariables.DNp_DX).resize(num_p_nodes, n_dim, false); + Vector detJpContainer = ZeroVector(num_g_points); try { mpPressureGeometry->ShapeFunctionsIntegrationPointsGradients( @@ -1083,8 +1027,8 @@ void SmallStrainUPwDiffOrderElement::InitializeElementVariables(ElementVariables // Variables computed at each integration point const SizeType VoigtSize = this->GetStressStatePolicy().GetVoigtSize(); - rVariables.B.resize(VoigtSize, NumUNodes * Dim, false); - noalias(rVariables.B) = ZeroMatrix(VoigtSize, NumUNodes * Dim); + rVariables.B.resize(VoigtSize, num_u_nodes * n_dim, false); + noalias(rVariables.B) = ZeroMatrix(VoigtSize, num_u_nodes * n_dim); rVariables.StrainVector.resize(VoigtSize, false); rVariables.ConstitutiveMatrix.resize(VoigtSize, VoigtSize, false); @@ -1092,8 +1036,8 @@ void SmallStrainUPwDiffOrderElement::InitializeElementVariables(ElementVariables rVariables.StressVector.resize(VoigtSize, false); // Needed parameters for consistency with the general constitutive law - rVariables.F.resize(Dim, Dim, false); - noalias(rVariables.F) = identity_matrix(Dim); + rVariables.F.resize(n_dim, n_dim, false); + noalias(rVariables.F) = identity_matrix(n_dim); // Nodal variables this->InitializeNodalVariables(rVariables); @@ -1117,43 +1061,43 @@ void SmallStrainUPwDiffOrderElement::InitializeNodalVariables(ElementVariables& { KRATOS_TRY - const GeometryType& rGeom = GetGeometry(); - const SizeType Dim = rGeom.WorkingSpaceDimension(); - const SizeType NumUNodes = rGeom.PointsNumber(); - const SizeType NumPNodes = mpPressureGeometry->PointsNumber(); + const GeometryType& r_geom = GetGeometry(); + const SizeType n_dim = r_geom.WorkingSpaceDimension(); + const SizeType num_u_nodes = r_geom.PointsNumber(); + const SizeType num_p_nodes = mpPressureGeometry->PointsNumber(); Vector BodyAccelerationAux = ZeroVector(3); - rVariables.BodyAcceleration.resize(NumUNodes * Dim, false); - rVariables.DisplacementVector.resize(NumUNodes * Dim, false); - rVariables.VelocityVector.resize(NumUNodes * Dim, false); + rVariables.BodyAcceleration.resize(num_u_nodes * n_dim, false); + rVariables.DisplacementVector.resize(num_u_nodes * n_dim, false); + rVariables.VelocityVector.resize(num_u_nodes * n_dim, false); - for (SizeType i = 0; i < NumUNodes; ++i) { - SizeType Local_i = i * Dim; - BodyAccelerationAux = rGeom[i].FastGetSolutionStepValue(VOLUME_ACCELERATION); + for (SizeType i = 0; i < num_u_nodes; ++i) { + SizeType Local_i = i * n_dim; + BodyAccelerationAux = r_geom[i].FastGetSolutionStepValue(VOLUME_ACCELERATION); rVariables.BodyAcceleration[Local_i] = BodyAccelerationAux[0]; - rVariables.DisplacementVector[Local_i] = rGeom[i].FastGetSolutionStepValue(DISPLACEMENT_X); - rVariables.VelocityVector[Local_i] = rGeom[i].FastGetSolutionStepValue(VELOCITY_X); + rVariables.DisplacementVector[Local_i] = r_geom[i].FastGetSolutionStepValue(DISPLACEMENT_X); + rVariables.VelocityVector[Local_i] = r_geom[i].FastGetSolutionStepValue(VELOCITY_X); rVariables.BodyAcceleration[Local_i + 1] = BodyAccelerationAux[1]; - rVariables.DisplacementVector[Local_i + 1] = rGeom[i].FastGetSolutionStepValue(DISPLACEMENT_Y); - rVariables.VelocityVector[Local_i + 1] = rGeom[i].FastGetSolutionStepValue(VELOCITY_Y); + rVariables.DisplacementVector[Local_i + 1] = r_geom[i].FastGetSolutionStepValue(DISPLACEMENT_Y); + rVariables.VelocityVector[Local_i + 1] = r_geom[i].FastGetSolutionStepValue(VELOCITY_Y); - if (Dim > 2) { + if (n_dim > 2) { rVariables.BodyAcceleration[Local_i + 2] = BodyAccelerationAux[2]; - rVariables.DisplacementVector[Local_i + 2] = rGeom[i].FastGetSolutionStepValue(DISPLACEMENT_Z); - rVariables.VelocityVector[Local_i + 2] = rGeom[i].FastGetSolutionStepValue(VELOCITY_Z); + rVariables.DisplacementVector[Local_i + 2] = r_geom[i].FastGetSolutionStepValue(DISPLACEMENT_Z); + rVariables.VelocityVector[Local_i + 2] = r_geom[i].FastGetSolutionStepValue(VELOCITY_Z); } } - rVariables.PressureVector.resize(NumPNodes, false); - rVariables.PressureDtVector.resize(NumPNodes, false); - rVariables.DeltaPressureVector.resize(NumPNodes, false); - for (SizeType i = 0; i < NumPNodes; ++i) { - rVariables.PressureVector[i] = rGeom[i].FastGetSolutionStepValue(WATER_PRESSURE); - rVariables.PressureDtVector[i] = rGeom[i].FastGetSolutionStepValue(DT_WATER_PRESSURE); - rVariables.DeltaPressureVector[i] = rGeom[i].FastGetSolutionStepValue(WATER_PRESSURE) - - rGeom[i].FastGetSolutionStepValue(WATER_PRESSURE, 1); + rVariables.PressureVector.resize(num_p_nodes, false); + rVariables.PressureDtVector.resize(num_p_nodes, false); + rVariables.DeltaPressureVector.resize(num_p_nodes, false); + for (SizeType i = 0; i < num_p_nodes; ++i) { + rVariables.PressureVector[i] = r_geom[i].FastGetSolutionStepValue(WATER_PRESSURE); + rVariables.PressureDtVector[i] = r_geom[i].FastGetSolutionStepValue(DT_WATER_PRESSURE); + rVariables.DeltaPressureVector[i] = r_geom[i].FastGetSolutionStepValue(WATER_PRESSURE) - + r_geom[i].FastGetSolutionStepValue(WATER_PRESSURE, 1); } KRATOS_CATCH("") @@ -1194,7 +1138,6 @@ void SmallStrainUPwDiffOrderElement::InitializeProperties(ElementVariables& rVar } void SmallStrainUPwDiffOrderElement::CalculateKinematics(ElementVariables& rVariables, unsigned int GPoint) - { KRATOS_TRY @@ -1237,15 +1180,16 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddLHS(MatrixType& rLeftH KRATOS_TRY this->CalculateAndAddStiffnessMatrix(rLeftHandSideMatrix, rVariables); - this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix, rVariables); - if (!rVariables.IgnoreUndrained) { - this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables); + this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables); + if (!rVariables.IgnoreUndrained) { const auto permeability_matrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( rVariables.DNp_DX, rVariables.DynamicViscosityInverse, rVariables.IntrinsicPermeability, rVariables.RelativePermeability, rVariables.IntegrationCoefficient); GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, permeability_matrix); + + this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix, rVariables); } KRATOS_CATCH("") @@ -1269,16 +1213,18 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddCouplingMatrix(MatrixType& r { KRATOS_TRY - Matrix CouplingMatrix = GeoTransportEquationUtilities::CalculateCouplingMatrix( + const Matrix coupling_matrix_up = GeoTransportEquationUtilities::CalculateCouplingMatrix( rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np, rVariables.BiotCoefficient, rVariables.BishopCoefficient, rVariables.IntegrationCoefficient); - GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, CouplingMatrix); + GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, coupling_matrix_up); if (!rVariables.IgnoreUndrained) { - const double SaturationCoefficient = rVariables.DegreeOfSaturation / rVariables.BishopCoefficient; - Matrix CouplingMatrixT = PORE_PRESSURE_SIGN_FACTOR * SaturationCoefficient * - rVariables.VelocityCoefficient * trans(CouplingMatrix); - GeoElementUtilities::AssemblePUBlockMatrix(rLeftHandSideMatrix, CouplingMatrixT); + const Matrix coupling_matrix_pu = GeoTransportEquationUtilities::CalculateCouplingMatrix( + rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np, + rVariables.BiotCoefficient, rVariables.DegreeOfSaturation, rVariables.IntegrationCoefficient); + GeoElementUtilities::AssemblePUBlockMatrix( + rLeftHandSideMatrix, + PORE_PRESSURE_SIGN_FACTOR * rVariables.VelocityCoefficient * trans(coupling_matrix_pu)); } KRATOS_CATCH("") @@ -1370,17 +1316,21 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddCouplingTerms(VectorType& rR { KRATOS_TRY - const Matrix coupling_matrix = + const Matrix u_coupling_matrix = (-1.0) * GeoTransportEquationUtilities::CalculateCouplingMatrix( rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np, rVariables.BiotCoefficient, rVariables.BishopCoefficient, rVariables.IntegrationCoefficient); - const Vector coupling_force = prod(coupling_matrix, rVariables.PressureVector); + const Vector coupling_force = prod(u_coupling_matrix, rVariables.PressureVector); GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, coupling_force); if (!rVariables.IgnoreUndrained) { - const double SaturationCoefficient = rVariables.DegreeOfSaturation / rVariables.BishopCoefficient; - const Vector coupling_flow = PORE_PRESSURE_SIGN_FACTOR * SaturationCoefficient * - prod(trans(coupling_matrix), rVariables.VelocityVector); + const Matrix p_coupling_matrix = + (-1.0) * GeoTransportEquationUtilities::CalculateCouplingMatrix( + rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np, + rVariables.BiotCoefficient, rVariables.DegreeOfSaturation, + rVariables.IntegrationCoefficient); + const Vector coupling_flow = + PORE_PRESSURE_SIGN_FACTOR * prod(trans(p_coupling_matrix), rVariables.VelocityVector); GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, coupling_flow); } @@ -1392,9 +1342,9 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddCompressibilityFlow(VectorTy { KRATOS_TRY - Matrix CompressibilityMatrix = GeoTransportEquationUtilities::CalculateCompressibilityMatrix( + Matrix compressibility_matrix = GeoTransportEquationUtilities::CalculateCompressibilityMatrix( rVariables.Np, rVariables.BiotModulusInverse, rVariables.IntegrationCoefficient); - Vector CompressibilityFlow = -prod(CompressibilityMatrix, rVariables.PressureDtVector); + Vector CompressibilityFlow = -prod(compressibility_matrix, rVariables.PressureDtVector); GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, CompressibilityFlow); KRATOS_CATCH("") @@ -1439,9 +1389,7 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddPermeabilityFlow(VectorType& -PORE_PRESSURE_SIGN_FACTOR * rVariables.DynamicViscosityInverse * rVariables.RelativePermeability * prod(rVariables.DNp_DX, Matrix(prod(rVariables.IntrinsicPermeability, trans(rVariables.DNp_DX)))) * rVariables.IntegrationCoefficient; - const Vector permeability_flow = -prod(permeability_matrix, rVariables.PressureVector); - GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, permeability_flow); KRATOS_CATCH("") @@ -1452,14 +1400,14 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddFluidBodyFlow(VectorType& rR { KRATOS_TRY - const Matrix grad_Np_T_perm = - rVariables.DynamicViscosityInverse * GetProperties()[DENSITY_WATER] * rVariables.RelativePermeability * - prod(rVariables.DNp_DX, rVariables.IntrinsicPermeability) * rVariables.IntegrationCoefficient; + const Matrix grad_Np_T_perm = rVariables.DynamicViscosityInverse * rVariables.BishopCoefficient * + GetProperties()[DENSITY_WATER] * rVariables.RelativePermeability * + prod(rVariables.DNp_DX, rVariables.IntrinsicPermeability) * + rVariables.IntegrationCoefficient; const GeometryType& r_geom = GetGeometry(); const SizeType dimension = r_geom.WorkingSpaceDimension(); const SizeType num_U_nodes = r_geom.PointsNumber(); - const SizeType num_P_nodes = mpPressureGeometry->PointsNumber(); Vector body_acceleration = ZeroVector(dimension); @@ -1471,10 +1419,8 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddFluidBodyFlow(VectorType& rR } } - for (SizeType i = 0; i < num_P_nodes; ++i) { - rRightHandSideVector[num_U_nodes * dimension + i] += - inner_prod(row(grad_Np_T_perm, i), body_acceleration); - } + const Vector fluid_body_flow = prod(grad_Np_T_perm, body_acceleration); + GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, fluid_body_flow); KRATOS_CATCH("") } 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 0aaf1fc9febf..88d8533762cc 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 @@ -224,7 +224,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub std::vector CalculateBMatrices(const GeometryType::ShapeFunctionsGradientsType& rDN_DXContainer, const Matrix& rNContainer) const; - void AssignPressureToIntermediateNodes(); + Vector GetPressures(size_t n_nodes) const; + void AssignPressureToIntermediateNodes(); virtual Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const; diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp index ffc103bbe2db..07f60f231d6b 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp @@ -137,8 +137,6 @@ void TransientPwElement::Initialize(const ProcessInfo& rCurrent if (mRetentionLawVector.size() != NumGPoints) mRetentionLawVector.resize(NumGPoints); for (unsigned int i = 0; i < mRetentionLawVector.size(); ++i) { mRetentionLawVector[i] = RetentionLawFactory::Clone(Prop); - mRetentionLawVector[i]->InitializeMaterial( - Prop, Geom, row(Geom.ShapeFunctionsValues(this->GetIntegrationMethod()), i)); } mIsInitialised = true; @@ -264,18 +262,6 @@ void TransientPwElement::InitializeSolutionStep(const ProcessIn if (!mIsInitialised) this->Initialize(rCurrentProcessInfo); - // Defining necessary variables - const GeometryType& Geom = this->GetGeometry(); - const unsigned int NumGPoints = Geom.IntegrationPointsNumber(this->GetIntegrationMethod()); - - RetentionLaw::Parameters RetentionParameters(this->GetProperties()); - - // Loop over integration points - for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - // Initialize retention law - mRetentionLawVector[GPoint]->InitializeSolutionStep(RetentionParameters); - } - // reset hydraulic discharge this->ResetHydraulicDischarge(); @@ -301,18 +287,6 @@ void TransientPwElement::FinalizeSolutionStep(const ProcessInfo this->CalculateHydraulicDischarge(rCurrentProcessInfo); - // Defining necessary variables - const GeometryType& Geom = this->GetGeometry(); - const unsigned int NumGPoints = Geom.IntegrationPointsNumber(this->GetIntegrationMethod()); - - RetentionLaw::Parameters RetentionParameters(this->GetProperties()); - - // Loop over integration points - for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - // retention law - mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters); - } - KRATOS_CATCH("") } diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.cpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.cpp index 88e984d28542..93976fcf0614 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.cpp @@ -146,33 +146,13 @@ void TransientPwInterfaceElement::CalculateMassMatrix(MatrixTyp template void TransientPwInterfaceElement::InitializeSolutionStep(const ProcessInfo&) { - KRATOS_TRY - - RetentionLaw::Parameters RetentionParameters(this->GetProperties()); - - // Loop over integration points - for (unsigned int GPoint = 0; GPoint < mRetentionLawVector.size(); ++GPoint) { - // Initialize retention law - mRetentionLawVector[GPoint]->InitializeSolutionStep(RetentionParameters); - } - - KRATOS_CATCH("") + // Intentionally empty } template void TransientPwInterfaceElement::FinalizeSolutionStep(const ProcessInfo&) { - KRATOS_TRY - - RetentionLaw::Parameters RetentionParameters(this->GetProperties()); - - // Loop over integration points - for (unsigned int GPoint = 0; GPoint < mRetentionLawVector.size(); ++GPoint) { - // retention law - mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters); - } - - KRATOS_CATCH("") + // Intentionally empty } template diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h index 168ac01c9dce..0b62c7713622 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h @@ -75,27 +75,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem { mRetentionLawVector.resize(GetGeometry().IntegrationPointsNumber(GetIntegrationMethod())); - for (unsigned int i = 0; i < mRetentionLawVector.size(); ++i) { - mRetentionLawVector[i] = RetentionLawFactory::Clone(GetProperties()); - mRetentionLawVector[i]->InitializeMaterial( - GetProperties(), GetGeometry(), - row(GetGeometry().ShapeFunctionsValues(GetIntegrationMethod()), i)); - } - } - - void InitializeSolutionStep(const ProcessInfo&) override - { - RetentionLaw::Parameters RetentionParameters(this->GetProperties()); - for (const auto& retention_law : mRetentionLawVector) { - retention_law->InitializeSolutionStep(RetentionParameters); - } - } - - void FinalizeSolutionStep(const ProcessInfo&) override - { - RetentionLaw::Parameters RetentionParameters(this->GetProperties()); - for (const auto& retention_law : mRetentionLawVector) { - retention_law->FinalizeSolutionStep(RetentionParameters); + for (auto& r_retention_law : mRetentionLawVector) { + r_retention_law = RetentionLawFactory::Clone(GetProperties()); } } diff --git a/applications/GeoMechanicsApplication/custom_retention/retention_law.cpp b/applications/GeoMechanicsApplication/custom_retention/retention_law.cpp index 0f58109062b1..06bd9ef2aeb4 100644 --- a/applications/GeoMechanicsApplication/custom_retention/retention_law.cpp +++ b/applications/GeoMechanicsApplication/custom_retention/retention_law.cpp @@ -11,45 +11,11 @@ // /* Project includes */ -#include "custom_retention/saturated_law.h" +#include "custom_retention/retention_law.h" namespace Kratos { -void RetentionLaw::InitializeMaterial(const Properties& rMaterialProperties, - const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues) -{ - // nothing -} - -void RetentionLaw::Initialize(Parameters& rParameters) -{ - // nothing -} - -void RetentionLaw::InitializeSolutionStep(Parameters& rParameters) -{ - // nothing -} - -void RetentionLaw::FinalizeSolutionStep(Parameters& rParameters) -{ - // nothing -} - -void RetentionLaw::Finalize(Parameters& rParameters) -{ - // nothing -} - -void RetentionLaw::ResetMaterial(const Properties& rMaterialProperties, - const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues) -{ - // nothing -} - void RetentionLaw::save(Serializer& rSerializer) const { // there is no member variables to be saved diff --git a/applications/GeoMechanicsApplication/custom_retention/retention_law.h b/applications/GeoMechanicsApplication/custom_retention/retention_law.h index 679cdcd31c98..73a6bd1f67f9 100644 --- a/applications/GeoMechanicsApplication/custom_retention/retention_law.h +++ b/applications/GeoMechanicsApplication/custom_retention/retention_law.h @@ -11,12 +11,6 @@ // #pragma once - -/* System includes */ - -/* External includes */ - -/* Project includes */ #include "geometries/geometry.h" #include "includes/define.h" #include "includes/process_info.h" @@ -33,17 +27,9 @@ namespace Kratos 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 GeometryType = Geometry; - /** - * Counted pointer of RetentionLaw - */ + // Counted pointer of RetentionLaw KRATOS_CLASS_POINTER_DEFINITION(RetentionLaw); class Parameters @@ -52,19 +38,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw /** * Structure "Parameters" to be used by the element to pass the parameters into the retention law * - - * KINEMATIC PARAMETERS: - - *** NOTE: Pointers are used only to point to a certain variable, - * no "new" or "malloc" can be used for this Parameters *** - - * MATERIAL PROPERTIES: - * @param mrMaterialProperties reference to the material's Properties object (input data) - - * PROCESS PROPERTIES: - * @param mrCurrentProcessInfo reference to current ProcessInfo instance (input data) - - */ + */ public: explicit Parameters(const Properties& rMaterialProperties) @@ -76,7 +50,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw void SetFluidPressure(double FluidPressure) { mFluidPressure = FluidPressure; }; - double GetFluidPressure() const + [[nodiscard]] double GetFluidPressure() const { KRATOS_ERROR_IF_NOT(mFluidPressure.has_value()) << "Fluid pressure is not yet set in the retention " @@ -84,7 +58,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw return mFluidPressure.value(); } - const Properties& GetMaterialProperties() const { return mrMaterialProperties; } + [[nodiscard]] const Properties& GetMaterialProperties() const + { + return mrMaterialProperties; + } private: std::optional mFluidPressure; @@ -92,8 +69,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw }; // class Parameters end - RetentionLaw() = default; - virtual ~RetentionLaw() = default; /** @@ -103,11 +78,11 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw * RetentionLaw::Pointer p_clone(new RetentionLaw()); * return p_clone; */ - virtual RetentionLaw::Pointer Clone() const = 0; + [[nodiscard]] virtual Pointer Clone() const = 0; /** * @brief Calculates the value of a specified variable (double) - * @param rParameterValues the needed parameters for the CL calculation + * @param rParameters the needed parameters for the CL calculation * @param rThisVariable the variable to be returned * @param rValue a reference to the returned value * @param rValue output: the value of the specified variable @@ -124,56 +99,11 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw virtual double CalculateBishopCoefficient(Parameters& rParameters) const = 0; - /** - * This is to be called at the very beginning of the calculation - * (e.g. from InitializeElement) in order to initialize all relevant - * attributes of the retention law - * @param rMaterialProperties the Properties instance of the current element - * @param rElementGeometry the geometry of the current element - * @param rCurrentProcessInfo process info - */ - virtual void InitializeMaterial(const Properties& rMaterialProperties, - const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues); - - virtual void Initialize(Parameters& rParameters); - - /** - * to be called at the beginning of each solution step - * (e.g. from Element::InitializeSolutionStep) - */ - virtual void InitializeSolutionStep(Parameters& rParameters); - - /** - * to be called at the end of each solution step - * (e.g. from Element::FinalizeSolutionStep) - */ - virtual void FinalizeSolutionStep(Parameters& rParameters); - - /** - * Finalize the material response in terms of Cauchy stresses - * @see Parameters - */ - virtual void Finalize(Parameters& rParameters); - - /** - * This can be used in order to reset all internal variables of the - * retention law (e.g. if a model should be reset to its reference state) - * @param rMaterialProperties the Properties instance of the current element - * @param rElementGeometry the geometry of the current element - * @param rShapeFunctionsValues the shape functions values in the current integration point - * @param the current ProcessInfo instance - */ - virtual void ResetMaterial(const Properties& rMaterialProperties, - const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues); - /** * This function is designed to be called once to perform all the checks * needed on the input provided. Checks can be "expensive" as the function * is designed to catch user's errors. * @param rMaterialProperties - * @param rElementGeometry * @param rCurrentProcessInfo * @return */ @@ -196,16 +126,13 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw */ inline static bool HasSameType(const RetentionLaw* rLHS, const RetentionLaw* rRHS) { - return RetentionLaw::HasSameType(*rLHS, *rRHS); + return HasSameType(*rLHS, *rRHS); } - /// Turn back information as a string. - virtual std::string Info() const { return "RetentionLaw"; } + [[nodiscard]] virtual std::string Info() const { return "RetentionLaw"; } - /// Print information about this object. 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"; } private: @@ -217,10 +144,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw }; /* Class RetentionLaw */ -/// input stream function -inline std::istream& operator>>(std::istream& rIStream, RetentionLaw& rThis); - -/// output stream function +// output stream function inline std::ostream& operator<<(std::ostream& rOStream, const RetentionLaw& rThis) { rThis.PrintInfo(rOStream); @@ -230,4 +154,4 @@ inline std::ostream& operator<<(std::ostream& rOStream, const RetentionLaw& rThi return rOStream; } -} /* namespace Kratos.*/ +} /* namespace Kratos.*/ \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/custom_retention/saturated_below_phreatic_level_law.cpp b/applications/GeoMechanicsApplication/custom_retention/saturated_below_phreatic_level_law.cpp index 4bb447ec0c5a..e659f732bd98 100644 --- a/applications/GeoMechanicsApplication/custom_retention/saturated_below_phreatic_level_law.cpp +++ b/applications/GeoMechanicsApplication/custom_retention/saturated_below_phreatic_level_law.cpp @@ -11,155 +11,92 @@ // Main authors: Vahid Galavi // -// System includes -#include - -// External includes - -// Project includes #include "custom_retention/saturated_below_phreatic_level_law.h" +#include "geo_mechanics_application_variables.h" namespace Kratos { -SaturatedBelowPhreaticLevelLaw::SaturatedBelowPhreaticLevelLaw() : RetentionLaw() {} - -SaturatedBelowPhreaticLevelLaw::SaturatedBelowPhreaticLevelLaw(const SaturatedBelowPhreaticLevelLaw& rOther) - : RetentionLaw(rOther) -{ -} - RetentionLaw::Pointer SaturatedBelowPhreaticLevelLaw::Clone() const { return Kratos::make_shared(*this); } -SaturatedBelowPhreaticLevelLaw::~SaturatedBelowPhreaticLevelLaw() {} - double SaturatedBelowPhreaticLevelLaw::CalculateSaturation(Parameters& rParameters) const { - KRATOS_TRY - - const double& p = rParameters.GetFluidPressure(); - - if (p < 0.0) { - return rParameters.GetMaterialProperties()[SATURATED_SATURATION]; - } else { - return rParameters.GetMaterialProperties()[RESIDUAL_SATURATION]; - } - - KRATOS_CATCH("") + return (rParameters.GetFluidPressure() < 0.0) + ? rParameters.GetMaterialProperties()[SATURATED_SATURATION] + : rParameters.GetMaterialProperties()[RESIDUAL_SATURATION]; } double SaturatedBelowPhreaticLevelLaw::CalculateEffectiveSaturation(Parameters& rParameters) const { - KRATOS_TRY - - return CalculateSaturation(rParameters); + const auto& r_material_properties = rParameters.GetMaterialProperties(); + const auto sat_max = r_material_properties[SATURATED_SATURATION]; + const auto sat_min = r_material_properties[RESIDUAL_SATURATION]; - KRATOS_CATCH("") + return (CalculateSaturation(rParameters) - sat_min) / (sat_max - sat_min); } -double SaturatedBelowPhreaticLevelLaw::CalculateDerivativeOfSaturation(Parameters& rParameters) const +double SaturatedBelowPhreaticLevelLaw::CalculateDerivativeOfSaturation(Parameters&) const { return 0.0; } double SaturatedBelowPhreaticLevelLaw::CalculateRelativePermeability(Parameters& rParameters) const { - KRATOS_TRY - - const double& p = rParameters.GetFluidPressure(); - - if (p < 0.0) { - return 1.0; - } else { - return rParameters.GetMaterialProperties()[MINIMUM_RELATIVE_PERMEABILITY]; - } - - KRATOS_CATCH("") + return rParameters.GetFluidPressure() < 0.0 + ? 1.0 + : rParameters.GetMaterialProperties()[MINIMUM_RELATIVE_PERMEABILITY]; } double SaturatedBelowPhreaticLevelLaw::CalculateBishopCoefficient(Parameters& rParameters) const { - KRATOS_TRY - return CalculateEffectiveSaturation(rParameters); - - KRATOS_CATCH("") } -double& SaturatedBelowPhreaticLevelLaw::CalculateValue(RetentionLaw::Parameters& rParameterValues, - const Variable& rThisVariable, - double& rValue) +double& SaturatedBelowPhreaticLevelLaw::CalculateValue(Parameters& rParameterValues, + const Variable& rThisVariable, + double& rValue) { if (rThisVariable == DEGREE_OF_SATURATION) { rValue = this->CalculateSaturation(rParameterValues); - return rValue; } else if (rThisVariable == EFFECTIVE_SATURATION) { rValue = this->CalculateEffectiveSaturation(rParameterValues); - return rValue; } else if (rThisVariable == BISHOP_COEFFICIENT) { rValue = this->CalculateBishopCoefficient(rParameterValues); - return rValue; } else if (rThisVariable == DERIVATIVE_OF_SATURATION) { rValue = this->CalculateDerivativeOfSaturation(rParameterValues); - return rValue; } else if (rThisVariable == RELATIVE_PERMEABILITY) { rValue = this->CalculateRelativePermeability(rParameterValues); - return rValue; } - return rValue; } -void SaturatedBelowPhreaticLevelLaw::InitializeMaterial(const Properties& rMaterialProperties, - const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues) -{ - // nothing is needed -} - -void SaturatedBelowPhreaticLevelLaw::Initialize(Parameters& rParameters) -{ - // nothing is needed -} - -void SaturatedBelowPhreaticLevelLaw::InitializeSolutionStep(Parameters& rParameters) -{ - // nothing is needed -} - -void SaturatedBelowPhreaticLevelLaw::Finalize(Parameters& rParameters) -{ - // nothing is needed -} - -void SaturatedBelowPhreaticLevelLaw::FinalizeSolutionStep(Parameters& rParameters) -{ - // nothing is needed -} - -int SaturatedBelowPhreaticLevelLaw::Check(const Properties& rMaterialProperties, const ProcessInfo& rCurrentProcessInfo) +int SaturatedBelowPhreaticLevelLaw::Check(const Properties& rMaterialProperties, const ProcessInfo&) { KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(SATURATED_SATURATION)) - << "SATURATED_SATURATION is not available in material parameters" << std::endl; - KRATOS_ERROR_IF(rMaterialProperties[SATURATED_SATURATION] < 0.0) - << "SATURATED_SATURATION cannot be less than 0 " << std::endl; + << "SATURATED_SATURATION is not available in the parameters of material " + << rMaterialProperties.Id() << "." << std::endl; + KRATOS_ERROR_IF(rMaterialProperties[SATURATED_SATURATION] < 0.0 || rMaterialProperties[SATURATED_SATURATION] > 1.0) + << "SATURATED_SATURATION (" << rMaterialProperties[SATURATED_SATURATION] + << ") must be in the range [0.0, 1.0] for material " << rMaterialProperties.Id() << "." + << std::endl; KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(RESIDUAL_SATURATION)) - << "RESIDUAL_SATURATION is not available in material parameters" << std::endl; - KRATOS_DEBUG_ERROR_IF_NOT(rMaterialProperties[RESIDUAL_SATURATION] > 0.0) - << "RESIDUAL_SATURATION must be greater than 0 " << std::endl; - KRATOS_ERROR_IF(rMaterialProperties[RESIDUAL_SATURATION] > 1.0) - << "RESIDUAL_SATURATION cannot be greater than 1.0 " << std::endl; - - KRATOS_ERROR_IF(rMaterialProperties[SATURATED_SATURATION] < rMaterialProperties[RESIDUAL_SATURATION]) - << "RESIDUAL_SATURATION cannot be greater than SATURATED_SATURATION " << std::endl; + << "RESIDUAL_SATURATION is not available in the parameters of material " + << rMaterialProperties.Id() << "." << std::endl; + KRATOS_ERROR_IF(rMaterialProperties[RESIDUAL_SATURATION] < 0.0 || rMaterialProperties[RESIDUAL_SATURATION] >= rMaterialProperties[SATURATED_SATURATION]) + << "RESIDUAL_SATURATION (" << rMaterialProperties[RESIDUAL_SATURATION] + << ") must be in the range [0.0, " << rMaterialProperties[SATURATED_SATURATION] + << "> for material " << rMaterialProperties.Id() << "." << std::endl; KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(MINIMUM_RELATIVE_PERMEABILITY)) - << "MINIMUM_RELATIVE_PERMEABILITY is not available in material parameters" << std::endl; - KRATOS_ERROR_IF_NOT((rMaterialProperties[MINIMUM_RELATIVE_PERMEABILITY] > 0.0)) - << "MINIMUM_RELATIVE_PERMEABILITY must be greater than 0 " << std::endl; + << "MINIMUM_RELATIVE_PERMEABILITY is not available in the parameters of material " + << rMaterialProperties.Id() << "." << std::endl; + KRATOS_ERROR_IF(rMaterialProperties[MINIMUM_RELATIVE_PERMEABILITY] < 0.0 || rMaterialProperties[MINIMUM_RELATIVE_PERMEABILITY] > 1.0) + << "MINIMUM_RELATIVE_PERMEABILITY (" << rMaterialProperties[MINIMUM_RELATIVE_PERMEABILITY] + << ") must be in the range [0.0, 1.0] for material " << rMaterialProperties.Id() << "." + << std::endl; return 0; } diff --git a/applications/GeoMechanicsApplication/custom_retention/saturated_below_phreatic_level_law.h b/applications/GeoMechanicsApplication/custom_retention/saturated_below_phreatic_level_law.h index ea7b6378933a..ef9bb4349229 100644 --- a/applications/GeoMechanicsApplication/custom_retention/saturated_below_phreatic_level_law.h +++ b/applications/GeoMechanicsApplication/custom_retention/saturated_below_phreatic_level_law.h @@ -12,20 +12,10 @@ #pragma once -// System includes -#include "includes/define.h" -#include -#include - -// External includes - // Project includes #include "custom_retention/retention_law.h" #include "includes/serializer.h" -// Application includes -#include "geo_mechanics_application_variables.h" - namespace Kratos { /** @@ -38,32 +28,12 @@ namespace Kratos class KRATOS_API(GEO_MECHANICS_APPLICATION) SaturatedBelowPhreaticLevelLaw : public RetentionLaw { public: - /// The base class RetentionLaw type definition - using BaseType = RetentionLaw; - using GeometryType = Geometry; - /// The size type definition - using SizeType = std::size_t; - - /// Counted pointer of SaturatedBelowPhreaticLevelLaw + // Counted pointer of SaturatedBelowPhreaticLevelLaw KRATOS_CLASS_POINTER_DEFINITION(SaturatedBelowPhreaticLevelLaw); - SaturatedBelowPhreaticLevelLaw(); - - RetentionLaw::Pointer Clone() const override; - - SaturatedBelowPhreaticLevelLaw(const SaturatedBelowPhreaticLevelLaw& rOther); - - ~SaturatedBelowPhreaticLevelLaw() override; - - void InitializeMaterial(const Properties& rMaterialProperties, - const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues) override; - - void Initialize(Parameters& rParameters) override; - - void InitializeSolutionStep(Parameters& rParameters) override; + [[nodiscard]] RetentionLaw::Pointer Clone() const override; double CalculateSaturation(Parameters& rParameters) const override; @@ -75,10 +45,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SaturatedBelowPhreaticLevelLaw : pub double CalculateBishopCoefficient(Parameters& rParameters) const override; - void Finalize(Parameters& rParameters) override; - - void FinalizeSolutionStep(Parameters& rParameters) override; - /** * @brief It calculates the value of a specified variable (double case) * @param rParameterValues the needed parameters for the CL calculation @@ -94,7 +60,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SaturatedBelowPhreaticLevelLaw : pub * @brief This function provides the place to perform checks on the completeness of the input. * @details It is designed to be called only once (or anyway, not often) typically at the beginning of the calculations, so to verify that nothing is missing from the input or that no common error is found. * @param rMaterialProperties The properties of the material - * @param rElementGeometry The geometry of the element * @param rCurrentProcessInfo The current process info instance * @return 0 if OK, 1 otherwise */ diff --git a/applications/GeoMechanicsApplication/custom_retention/saturated_law.cpp b/applications/GeoMechanicsApplication/custom_retention/saturated_law.cpp index ec6ab2bfd3db..fc2a01b3670b 100644 --- a/applications/GeoMechanicsApplication/custom_retention/saturated_law.cpp +++ b/applications/GeoMechanicsApplication/custom_retention/saturated_law.cpp @@ -11,13 +11,8 @@ // Main authors: Vahid Galavi // -// System includes -#include - -// External includes - -// Project includes #include "custom_retention/saturated_law.h" +#include "geo_mechanics_application_variables.h" namespace Kratos { @@ -29,17 +24,8 @@ RetentionLaw::Pointer SaturatedLaw::Clone() const double SaturatedLaw::CalculateSaturation(Parameters& rParameters) const { - KRATOS_TRY - const Properties& rMaterialProperties = rParameters.GetMaterialProperties(); - - if (rMaterialProperties.Has(SATURATED_SATURATION)) { - return rMaterialProperties[SATURATED_SATURATION]; - } else { - return 1.0; - } - - KRATOS_CATCH("") + return rMaterialProperties.Has(SATURATED_SATURATION) ? rMaterialProperties[SATURATED_SATURATION] : 1.0; } double SaturatedLaw::CalculateEffectiveSaturation(Parameters& rParameters) const { return 1.0; } @@ -50,11 +36,7 @@ double SaturatedLaw::CalculateRelativePermeability(Parameters& rParameters) cons double SaturatedLaw::CalculateBishopCoefficient(Parameters& rParameters) const { - KRATOS_TRY - return CalculateEffectiveSaturation(rParameters); - - KRATOS_CATCH("") } double& SaturatedLaw::CalculateValue(RetentionLaw::Parameters& rParameterValues, @@ -63,59 +45,25 @@ double& SaturatedLaw::CalculateValue(RetentionLaw::Parameters& rParameterValues, { if (rThisVariable == DEGREE_OF_SATURATION) { rValue = this->CalculateSaturation(rParameterValues); - return rValue; } else if (rThisVariable == EFFECTIVE_SATURATION) { rValue = this->CalculateEffectiveSaturation(rParameterValues); - return rValue; } else if (rThisVariable == BISHOP_COEFFICIENT) { rValue = this->CalculateBishopCoefficient(rParameterValues); - return rValue; } else if (rThisVariable == DERIVATIVE_OF_SATURATION) { rValue = this->CalculateDerivativeOfSaturation(rParameterValues); - return rValue; } else if (rThisVariable == RELATIVE_PERMEABILITY) { rValue = this->CalculateRelativePermeability(rParameterValues); - return rValue; } - return rValue; } -void SaturatedLaw::InitializeMaterial(const Properties& rMaterialProperties, - const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues) -{ - // nothing is needed -} - -void SaturatedLaw::Initialize(Parameters& rParameters) -{ - // nothing is needed -} - -void SaturatedLaw::InitializeSolutionStep(Parameters& rParameters) -{ - // nothing is needed -} - -void SaturatedLaw::Finalize(Parameters& rParameters) -{ - // nothing is needed -} - -void SaturatedLaw::FinalizeSolutionStep(Parameters& rParameters) -{ - // nothing is needed -} - int SaturatedLaw::Check(const Properties& rMaterialProperties, const ProcessInfo& rCurrentProcessInfo) { if (rMaterialProperties.Has(SATURATED_SATURATION)) { - KRATOS_ERROR_IF(rMaterialProperties[SATURATED_SATURATION] < 0.0) - << "SATURATED_SATURATION cannot be less than 0 " << std::endl; - - KRATOS_ERROR_IF(rMaterialProperties[SATURATED_SATURATION] > 1.0) - << "SATURATED_SATURATION cannot be greater than 1.0 " << std::endl; + KRATOS_ERROR_IF(rMaterialProperties[SATURATED_SATURATION] < 0.0 || rMaterialProperties[SATURATED_SATURATION] > 1.0) + << "SATURATED_SATURATION (" << rMaterialProperties[SATURATED_SATURATION] + << ") must be in the range [0.0, 1.0] for material " << rMaterialProperties.Id() << "." + << std::endl; } return 0; diff --git a/applications/GeoMechanicsApplication/custom_retention/saturated_law.h b/applications/GeoMechanicsApplication/custom_retention/saturated_law.h index c37244e55890..611e926341db 100644 --- a/applications/GeoMechanicsApplication/custom_retention/saturated_law.h +++ b/applications/GeoMechanicsApplication/custom_retention/saturated_law.h @@ -12,20 +12,10 @@ #pragma once -// System includes -#include "includes/define.h" -#include -#include - -// External includes - -// Project includes #include "custom_retention/retention_law.h" +#include "includes/define.h" #include "includes/serializer.h" -// Application includes -#include "geo_mechanics_application_variables.h" - namespace Kratos { /** @@ -38,26 +28,12 @@ namespace Kratos class KRATOS_API(GEO_MECHANICS_APPLICATION) SaturatedLaw : public RetentionLaw { public: - /// The base class RetentionLaw type definition - using BaseType = RetentionLaw; - using GeometryType = Geometry; - /// The size type definition - using SizeType = std::size_t; - - /// Counted pointer of SaturatedLaw + // Counted pointer of SaturatedLaw KRATOS_CLASS_POINTER_DEFINITION(SaturatedLaw); - RetentionLaw::Pointer Clone() const override; - - void InitializeMaterial(const Properties& rMaterialProperties, - const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues) override; - - void Initialize(Parameters& rParameters) override; - - void InitializeSolutionStep(Parameters& rParameters) override; + [[nodiscard]] RetentionLaw::Pointer Clone() const override; double CalculateSaturation(Parameters& rParameters) const override; @@ -69,10 +45,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SaturatedLaw : public RetentionLaw double CalculateBishopCoefficient(Parameters& rParameters) const override; - void Finalize(Parameters& rParameters) override; - - void FinalizeSolutionStep(Parameters& rParameters) override; - /** * @brief It calculates the value of a specified variable (double case) * @param rParameterValues the needed parameters for the CL calculation @@ -80,15 +52,12 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SaturatedLaw : public RetentionLaw * @param rValue a reference to the returned value * @return rValue output: the value of the specified variable */ - double& CalculateValue(RetentionLaw::Parameters& rParameterValues, - const Variable& rThisVariable, - double& rValue) override; + double& CalculateValue(Parameters& rParameterValues, const Variable& rThisVariable, double& rValue) override; /** * @brief This function provides the place to perform checks on the completeness of the input. * @details It is designed to be called only once (or anyway, not often) typically at the beginning of the calculations, so to verify that nothing is missing from the input or that no common error is found. * @param rMaterialProperties The properties of the material - * @param rElementGeometry The geometry of the element * @param rCurrentProcessInfo The current process info instance * @return 0 if OK, 1 otherwise */ diff --git a/applications/GeoMechanicsApplication/custom_retention/van_genuchten_law.cpp b/applications/GeoMechanicsApplication/custom_retention/van_genuchten_law.cpp index f71ee4af5cab..95a910c70de3 100644 --- a/applications/GeoMechanicsApplication/custom_retention/van_genuchten_law.cpp +++ b/applications/GeoMechanicsApplication/custom_retention/van_genuchten_law.cpp @@ -11,11 +11,8 @@ // Main authors: Vahid Galavi // -// System includes -#include - -// Project includes #include "custom_retention/van_genuchten_law.h" +#include "geo_mechanics_application_variables.h" namespace Kratos { @@ -29,19 +26,19 @@ double VanGenuchtenLaw::CalculateSaturation(Parameters& rParameters) const { KRATOS_TRY - const double& p = rParameters.GetFluidPressure(); - const Properties& rMaterialProperties = rParameters.GetMaterialProperties(); + const auto p = rParameters.GetFluidPressure(); + const auto& r_material_properties = rParameters.GetMaterialProperties(); if (p > 0.0) { - const double& satMax = rMaterialProperties[SATURATED_SATURATION]; - const double& satMin = rMaterialProperties[RESIDUAL_SATURATION]; - const double& pb = rMaterialProperties[VAN_GENUCHTEN_AIR_ENTRY_PRESSURE]; - const double& gn = rMaterialProperties[VAN_GENUCHTEN_GN]; - const double gc = (1.0 - gn) / gn; + const auto sat_max = r_material_properties[SATURATED_SATURATION]; + const auto sat_min = r_material_properties[RESIDUAL_SATURATION]; + const auto pb = r_material_properties[VAN_GENUCHTEN_AIR_ENTRY_PRESSURE]; + const auto gn = r_material_properties[VAN_GENUCHTEN_GN]; + const auto gc = (1.0 - gn) / gn; - return satMin + (satMax - satMin) * pow(1.0 + pow(p / pb, gn), gc); + return sat_min + (sat_max - sat_min) * std::pow(1.0 + std::pow(p / pb, gn), gc); } else { - return rMaterialProperties[SATURATED_SATURATION]; + return r_material_properties[SATURATED_SATURATION]; } KRATOS_CATCH("") @@ -51,11 +48,11 @@ double VanGenuchtenLaw::CalculateEffectiveSaturation(Parameters& rParameters) co { KRATOS_TRY - const auto& rMaterialProperties = rParameters.GetMaterialProperties(); - const double& satMax = rMaterialProperties[SATURATED_SATURATION]; - const double& satMin = rMaterialProperties[RESIDUAL_SATURATION]; + const auto& r_material_properties = rParameters.GetMaterialProperties(); + const auto sat_max = r_material_properties[SATURATED_SATURATION]; + const auto sat_min = r_material_properties[RESIDUAL_SATURATION]; - return (CalculateSaturation(rParameters) - satMin) / (satMax - satMin); + return (CalculateSaturation(rParameters) - sat_min) / (sat_max - sat_min); KRATOS_CATCH("") } @@ -63,18 +60,18 @@ double VanGenuchtenLaw::CalculateEffectiveSaturation(Parameters& rParameters) co double VanGenuchtenLaw::CalculateDerivativeOfSaturation(Parameters& rParameters) const { KRATOS_TRY - const double& p = rParameters.GetFluidPressure(); + const auto p = rParameters.GetFluidPressure(); if (p > 0.0) { - const auto& rMaterialProperties = rParameters.GetMaterialProperties(); - const double& satMax = rMaterialProperties[SATURATED_SATURATION]; - const double& satMin = rMaterialProperties[RESIDUAL_SATURATION]; - const double& pb = rMaterialProperties[VAN_GENUCHTEN_AIR_ENTRY_PRESSURE]; - const double& gn = rMaterialProperties[VAN_GENUCHTEN_GN]; - const double gc = (1.0 - gn) / gn; - - return (satMax - satMin) * gc * pow((1.0 + pow(p / pb, gn)), gc - 1.0) * gn * pow(pb, -gn) * - pow(p, gn - 1.0); + const auto& r_material_properties = rParameters.GetMaterialProperties(); + const auto sat_max = r_material_properties[SATURATED_SATURATION]; + const auto sat_min = r_material_properties[RESIDUAL_SATURATION]; + const auto pb = r_material_properties[VAN_GENUCHTEN_AIR_ENTRY_PRESSURE]; + const auto gn = r_material_properties[VAN_GENUCHTEN_GN]; + const auto gc = (1.0 - gn) / gn; + + return (sat_max - sat_min) * gc * std::pow((1.0 + std::pow(p / pb, gn)), gc - 1.0) * gn * + std::pow(pb, -gn) * std::pow(p, gn - 1.0); } else { return 0.0; } @@ -86,108 +83,91 @@ double VanGenuchtenLaw::CalculateRelativePermeability(Parameters& rParameters) c { KRATOS_TRY - const double effSat = CalculateEffectiveSaturation(rParameters); + const auto eff_sat = CalculateEffectiveSaturation(rParameters); - const auto& rMaterialProperties = rParameters.GetMaterialProperties(); - const double& gl = rMaterialProperties[VAN_GENUCHTEN_GL]; - const double& gn = rMaterialProperties[VAN_GENUCHTEN_GN]; + const auto& r_material_properties = rParameters.GetMaterialProperties(); + const auto gl = r_material_properties[VAN_GENUCHTEN_GL]; + const auto gn = r_material_properties[VAN_GENUCHTEN_GN]; - double relPerm = - pow(effSat, gl) * pow(1.0 - pow(1.0 - pow(effSat, gn / (gn - 1.0)), (gn - 1.0) / gn), 2); + const auto rel_perm = + std::pow(eff_sat, gl) * + std::pow(1.0 - std::pow(1.0 - std::pow(eff_sat, gn / (gn - 1.0)), (gn - 1.0) / gn), 2); - return std::max(relPerm, rMaterialProperties[MINIMUM_RELATIVE_PERMEABILITY]); + return std::max(rel_perm, r_material_properties[MINIMUM_RELATIVE_PERMEABILITY]); KRATOS_CATCH("") } double VanGenuchtenLaw::CalculateBishopCoefficient(Parameters& rParameters) const { - KRATOS_TRY - return CalculateEffectiveSaturation(rParameters); - - KRATOS_CATCH("") } -double& VanGenuchtenLaw::CalculateValue(RetentionLaw::Parameters& rParameterValues, - const Variable& rThisVariable, - double& rValue) +double& VanGenuchtenLaw::CalculateValue(Parameters& rParameterValues, const Variable& rThisVariable, double& rValue) { if (rThisVariable == DEGREE_OF_SATURATION) { rValue = this->CalculateSaturation(rParameterValues); - return rValue; } else if (rThisVariable == EFFECTIVE_SATURATION) { rValue = this->CalculateEffectiveSaturation(rParameterValues); - return rValue; } else if (rThisVariable == BISHOP_COEFFICIENT) { rValue = this->CalculateBishopCoefficient(rParameterValues); - return rValue; } else if (rThisVariable == DERIVATIVE_OF_SATURATION) { rValue = this->CalculateDerivativeOfSaturation(rParameterValues); - return rValue; } else if (rThisVariable == RELATIVE_PERMEABILITY) { rValue = this->CalculateRelativePermeability(rParameterValues); - return rValue; } return rValue; } -void VanGenuchtenLaw::InitializeMaterial(const Properties& rMaterialProperties, - const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues) -{ - // nothing is needed -} - -void VanGenuchtenLaw::Initialize(Parameters& rParameters) -{ - // nothing is needed -} - -void VanGenuchtenLaw::InitializeSolutionStep(Parameters& rParameters) -{ - // nothing is needed -} - -void VanGenuchtenLaw::Finalize(Parameters& rParameters) -{ - // nothing is needed -} - -void VanGenuchtenLaw::FinalizeSolutionStep(Parameters& rParameters) -{ - // nothing is needed -} - int VanGenuchtenLaw::Check(const Properties& rMaterialProperties, const ProcessInfo& rCurrentProcessInfo) { KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(SATURATED_SATURATION)) - << "SATURATED_SATURATION is not available in material parameters" << std::endl; - KRATOS_ERROR_IF(rMaterialProperties[SATURATED_SATURATION] < 0.0) - << "SATURATED_SATURATION cannot be less than 0 " << std::endl; - KRATOS_ERROR_IF(rMaterialProperties[SATURATED_SATURATION] > 1.0) - << "SATURATED_SATURATION cannot be greater than 1.0 " << std::endl; + << "SATURATED_SATURATION is not available in the parameters of material " + << rMaterialProperties.Id() << "." << std::endl; + KRATOS_ERROR_IF(rMaterialProperties[SATURATED_SATURATION] < 0.0 || rMaterialProperties[SATURATED_SATURATION] > 1.0) + << "SATURATED_SATURATION (" << rMaterialProperties[SATURATED_SATURATION] + << ") must be in the range [0.0, 1.0] for material " << rMaterialProperties.Id() << "." + << std::endl; KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(RESIDUAL_SATURATION)) - << "RESIDUAL_SATURATION is not available in material parameters" << std::endl; - KRATOS_DEBUG_ERROR_IF_NOT(rMaterialProperties[RESIDUAL_SATURATION] > 0.0) - << "RESIDUAL_SATURATION must be greater than 0 " << std::endl; - KRATOS_ERROR_IF(rMaterialProperties[RESIDUAL_SATURATION] > 1.0) - << "RESIDUAL_SATURATION cannot be greater than 1.0 " << std::endl; + << "RESIDUAL_SATURATION is not available in the parameters of material " + << rMaterialProperties.Id() << "." << std::endl; + KRATOS_ERROR_IF(rMaterialProperties[RESIDUAL_SATURATION] < 0.0 || rMaterialProperties[RESIDUAL_SATURATION] >= rMaterialProperties[SATURATED_SATURATION]) + << "RESIDUAL_SATURATION (" << rMaterialProperties[RESIDUAL_SATURATION] + << ") must be in the range [0.0, " << rMaterialProperties[SATURATED_SATURATION] + << "> for material " << rMaterialProperties.Id() << "." << std::endl; - KRATOS_ERROR_IF(rMaterialProperties[SATURATED_SATURATION] < rMaterialProperties[RESIDUAL_SATURATION]) - << "RESIDUAL_SATURATION cannot be greater than SATURATED_SATURATION " << std::endl; + KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(MINIMUM_RELATIVE_PERMEABILITY)) + << "MINIMUM_RELATIVE_PERMEABILITY is not available in the parameters of material " + << rMaterialProperties.Id() << "." << std::endl; + KRATOS_ERROR_IF(rMaterialProperties[MINIMUM_RELATIVE_PERMEABILITY] < 0.0 || rMaterialProperties[MINIMUM_RELATIVE_PERMEABILITY] > 1.0) + << "MINIMUM_RELATIVE_PERMEABILITY (" << rMaterialProperties[MINIMUM_RELATIVE_PERMEABILITY] + << ") must be in the range [0.0, 1.0] for material " << rMaterialProperties.Id() << "." + << std::endl; KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(VAN_GENUCHTEN_AIR_ENTRY_PRESSURE)) - << "VAN_GENUCHTEN_AIR_ENTRY_PRESSURE is not available in material parameters" << std::endl; + << "VAN_GENUCHTEN_AIR_ENTRY_PRESSURE is not available in the parameters of material " + << rMaterialProperties.Id() << "." << std::endl; KRATOS_ERROR_IF_NOT((rMaterialProperties[VAN_GENUCHTEN_AIR_ENTRY_PRESSURE] > 0.0)) - << "VAN_GENUCHTEN_AIR_ENTRY_PRESSURE must be greater than 0 " << std::endl; - - KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(MINIMUM_RELATIVE_PERMEABILITY)) - << "MINIMUM_RELATIVE_PERMEABILITY is not available in material parameters" << std::endl; - KRATOS_ERROR_IF_NOT((rMaterialProperties[MINIMUM_RELATIVE_PERMEABILITY] > 0.0)) - << "MINIMUM_RELATIVE_PERMEABILITY must be greater than 0 " << std::endl; + << "VAN_GENUCHTEN_AIR_ENTRY_PRESSURE (" + << rMaterialProperties[VAN_GENUCHTEN_AIR_ENTRY_PRESSURE] << ") must be greater than 0 " + << "for material " << rMaterialProperties.Id() << "." << std::endl; + + KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(VAN_GENUCHTEN_GN)) + << "VAN_GENUCHTEN_GN is not available in the parameters of material " + << rMaterialProperties.Id() << "." << std::endl; + KRATOS_ERROR_IF_NOT((rMaterialProperties[VAN_GENUCHTEN_GN] > 0.0)) + << "VAN_GENUCHTEN_GN (" << rMaterialProperties[VAN_GENUCHTEN_GN] << ") must be greater than 0 " + << "for material " << rMaterialProperties.Id() << "." << std::endl; + + KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(VAN_GENUCHTEN_GL)) + << "VAN_GENUCHTEN_GL is not available in the parameters of material " + << rMaterialProperties.Id() << "." << std::endl; + KRATOS_ERROR_IF_NOT((rMaterialProperties[VAN_GENUCHTEN_GL] >= 0.0)) + << "VAN_GENUCHTEN_GL (" << rMaterialProperties[VAN_GENUCHTEN_GL] + << ") must be greater than or equal to 0 for material " << rMaterialProperties.Id() << "." + << std::endl; return 0; } diff --git a/applications/GeoMechanicsApplication/custom_retention/van_genuchten_law.h b/applications/GeoMechanicsApplication/custom_retention/van_genuchten_law.h index 9174370beedf..ca008a74c8e9 100644 --- a/applications/GeoMechanicsApplication/custom_retention/van_genuchten_law.h +++ b/applications/GeoMechanicsApplication/custom_retention/van_genuchten_law.h @@ -12,20 +12,9 @@ #pragma once -// System includes -#include "includes/define.h" -#include -#include - -// External includes - -// Project includes #include "custom_retention/retention_law.h" #include "includes/serializer.h" -// Application includes -#include "geo_mechanics_application_variables.h" - namespace Kratos { /** @@ -38,26 +27,12 @@ namespace Kratos class KRATOS_API(GEO_MECHANICS_APPLICATION) VanGenuchtenLaw : public RetentionLaw { public: - /// The base class RetentionLaw type definition - using BaseType = RetentionLaw; - using GeometryType = Geometry; - /// The size type definition - using SizeType = std::size_t; - - /// Counted pointer of VanGenuchtenLaw + // Counted pointer of VanGenuchtenLaw KRATOS_CLASS_POINTER_DEFINITION(VanGenuchtenLaw); - RetentionLaw::Pointer Clone() const override; - - void InitializeMaterial(const Properties& rMaterialProperties, - const GeometryType& rElementGeometry, - const Vector& rShapeFunctionsValues) override; - - void Initialize(Parameters& rParameters) override; - - void InitializeSolutionStep(Parameters& rParameters) override; + [[nodiscard]] RetentionLaw::Pointer Clone() const override; double CalculateSaturation(Parameters& rParameters) const override; @@ -69,10 +44,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) VanGenuchtenLaw : public RetentionLa double CalculateBishopCoefficient(Parameters& rParameters) const override; - void Finalize(Parameters& rParameters) override; - - void FinalizeSolutionStep(Parameters& rParameters) override; - /** * @brief It calculates the value of a specified variable (double case) * @param rParameterValues the needed parameters for the CL calculation @@ -80,15 +51,12 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) VanGenuchtenLaw : public RetentionLa * @param rValue a reference to the returned value * @return rValue output: the value of the specified variable */ - double& CalculateValue(RetentionLaw::Parameters& rParameterValues, - const Variable& rThisVariable, - double& rValue) override; + double& CalculateValue(Parameters& rParameterValues, const Variable& rThisVariable, double& rValue) override; /** * @brief This function provides the place to perform checks on the completeness of the input. * @details It is designed to be called only once (or anyway, not often) typically at the beginning of the calculations, so to verify that nothing is missing from the input or that no common error is found. * @param rMaterialProperties The properties of the material - * @param rElementGeometry The geometry of the element * @param rCurrentProcessInfo The current process info instance * @return 0 if OK, 1 otherwise */ diff --git a/applications/GeoMechanicsApplication/tests/Interface/interface_on_beam/MaterialParameters.json b/applications/GeoMechanicsApplication/tests/Interface/interface_on_beam/MaterialParameters.json index 7e6612f8d59e..a495c0291ba9 100644 --- a/applications/GeoMechanicsApplication/tests/Interface/interface_on_beam/MaterialParameters.json +++ b/applications/GeoMechanicsApplication/tests/Interface/interface_on_beam/MaterialParameters.json @@ -47,7 +47,7 @@ "name" : "SmallStrainUDSM2DInterfaceLaw" }, "Variables": { - "IGNORE_UNDRAINED" : false, + "IGNORE_UNDRAINED" : true, "DENSITY_SOLID" : 2.0e3, "DENSITY_WATER" : 1.0e3, "POROSITY" : 0.3, @@ -78,7 +78,7 @@ "name" : "SmallStrainUDSM2DInterfaceLaw" }, "Variables": { - "IGNORE_UNDRAINED" : false, + "IGNORE_UNDRAINED" : true, "DENSITY_SOLID" : 2.0e3, "DENSITY_WATER" : 1.0e3, "POROSITY" : 0.3, diff --git a/applications/GeoMechanicsApplication/tests/Interface/weak_interface_on_beam/MaterialParameters.json b/applications/GeoMechanicsApplication/tests/Interface/weak_interface_on_beam/MaterialParameters.json index 1d09e969c7f0..0b7f1e22f148 100644 --- a/applications/GeoMechanicsApplication/tests/Interface/weak_interface_on_beam/MaterialParameters.json +++ b/applications/GeoMechanicsApplication/tests/Interface/weak_interface_on_beam/MaterialParameters.json @@ -47,7 +47,7 @@ "name" : "SmallStrainUDSM2DInterfaceLaw" }, "Variables": { - "IGNORE_UNDRAINED" : false, + "IGNORE_UNDRAINED" : true, "DENSITY_SOLID" : 2.0e3, "DENSITY_WATER" : 1.0e3, "POROSITY" : 0.3, @@ -78,7 +78,7 @@ "name" : "SmallStrainUDSM2DInterfaceLaw" }, "Variables": { - "IGNORE_UNDRAINED" : false, + "IGNORE_UNDRAINED" : true, "DENSITY_SOLID" : 2.0e3, "DENSITY_WATER" : 1.0e3, "POROSITY" : 0.3, diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_retention/test_saturated_below_phreatic_level_law.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_retention/test_saturated_below_phreatic_level_law.cpp new file mode 100644 index 000000000000..5e44adc01c13 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_retention/test_saturated_below_phreatic_level_law.cpp @@ -0,0 +1,101 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Wijtze Pieter Kikstra +// + +#include "custom_retention/saturated_below_phreatic_level_law.h" +#include "geo_mechanics_application_variables.h" +#include "tests/cpp_tests/geo_mechanics_fast_suite.h" + +namespace Kratos::Testing +{ + +KRATOS_TEST_CASE_IN_SUITE(SaturatedBelowPhreaticLevelLawReturnsCloneOfCorrectType, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto law = SaturatedBelowPhreaticLevelLaw(); + const auto p_law_clone = law.Clone(); + KRATOS_EXPECT_NE(&law, p_law_clone.get()); + KRATOS_EXPECT_NE(dynamic_cast(p_law_clone.get()), nullptr); +} + +KRATOS_TEST_CASE_IN_SUITE(SaturatedBelowPhreaticLevelLawReturnsCalculatedValues, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law = SaturatedBelowPhreaticLevelLaw(); + Properties properties; + properties.SetValue(SATURATED_SATURATION, 0.9); + properties.SetValue(RESIDUAL_SATURATION, 0.1); + properties.SetValue(MINIMUM_RELATIVE_PERMEABILITY, 0.05); + + auto retention_law_parameters = RetentionLaw::Parameters{properties}; + + retention_law_parameters.SetFluidPressure(-10.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateSaturation(retention_law_parameters), 0.9); + double value = 0.0; + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DEGREE_OF_SATURATION, value), 0.9); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateEffectiveSaturation(retention_law_parameters), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, EFFECTIVE_SATURATION, value), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateBishopCoefficient(retention_law_parameters), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, BISHOP_COEFFICIENT, value), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateDerivativeOfSaturation(retention_law_parameters), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DERIVATIVE_OF_SATURATION, value), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateRelativePermeability(retention_law_parameters), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, RELATIVE_PERMEABILITY, value), 1.0); + + retention_law_parameters.SetFluidPressure(10.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateSaturation(retention_law_parameters), 0.1); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DEGREE_OF_SATURATION, value), 0.1); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateEffectiveSaturation(retention_law_parameters), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, EFFECTIVE_SATURATION, value), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateBishopCoefficient(retention_law_parameters), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, BISHOP_COEFFICIENT, value), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateDerivativeOfSaturation(retention_law_parameters), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DERIVATIVE_OF_SATURATION, value), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateRelativePermeability(retention_law_parameters), 0.05); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, RELATIVE_PERMEABILITY, value), 0.05); +} + +KRATOS_TEST_CASE_IN_SUITE(SaturatedBelowPhreaticLevelLawChecksInputParameters, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + Properties properties; + properties.SetId(1); + const auto process_info = ProcessInfo{}; + auto law = SaturatedBelowPhreaticLevelLaw(); + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "SATURATED_SATURATION is not available in the parameters of material 1."); + properties.SetValue(SATURATED_SATURATION, 1.1); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "SATURATED_SATURATION (1.1) must be in the range [0.0, 1.0] for material 1."); + properties.SetValue(SATURATED_SATURATION, 0.9); + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "RESIDUAL_SATURATION is not available in the parameters of material 1."); + properties.SetValue(RESIDUAL_SATURATION, 1.1); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "RESIDUAL_SATURATION (1.1) must be in the range [0.0, 0.9> for material 1."); + properties.SetValue(RESIDUAL_SATURATION, 0.1); + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "MINIMUM_RELATIVE_PERMEABILITY is not available in the parameters of material 1."); + properties.SetValue(MINIMUM_RELATIVE_PERMEABILITY, 1.1); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "MINIMUM_RELATIVE_PERMEABILITY (1.1) must be in the range [0.0, 1.0] for material 1."); + properties.SetValue(MINIMUM_RELATIVE_PERMEABILITY, 0.05); + + KRATOS_EXPECT_EQ(law.Check(properties, process_info), 0); +} + +} // namespace Kratos::Testing \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_retention/test_saturated_law.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_retention/test_saturated_law.cpp new file mode 100644 index 000000000000..d5b901cd0f0a --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_retention/test_saturated_law.cpp @@ -0,0 +1,69 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Wijtze Pieter Kikstra +// + +#include "custom_retention/saturated_law.h" +#include "geo_mechanics_application_variables.h" +#include "tests/cpp_tests/geo_mechanics_fast_suite.h" + +namespace Kratos::Testing +{ + +KRATOS_TEST_CASE_IN_SUITE(SaturatedLawReturnsCloneOfCorrectType, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto law = SaturatedLaw(); + const auto p_law_clone = law.Clone(); + KRATOS_EXPECT_NE(&law, p_law_clone.get()); + KRATOS_EXPECT_NE(dynamic_cast(p_law_clone.get()), nullptr); +} + +KRATOS_TEST_CASE_IN_SUITE(SaturatedLawReturnsCalculatedValues, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law = SaturatedLaw(); + Properties properties; + auto retention_law_parameters = RetentionLaw::Parameters{properties}; + + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateSaturation(retention_law_parameters), 1.0); + double value = 0.0; + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DEGREE_OF_SATURATION, value), 1.0); + + properties.SetValue(SATURATED_SATURATION, 0.9); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateSaturation(retention_law_parameters), 0.9); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DEGREE_OF_SATURATION, value), 0.9); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateEffectiveSaturation(retention_law_parameters), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, EFFECTIVE_SATURATION, value), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateBishopCoefficient(retention_law_parameters), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, BISHOP_COEFFICIENT, value), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateDerivativeOfSaturation(retention_law_parameters), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DERIVATIVE_OF_SATURATION, value), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateRelativePermeability(retention_law_parameters), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, RELATIVE_PERMEABILITY, value), 1.0); +} + +KRATOS_TEST_CASE_IN_SUITE(SaturatedLawChecksInputParameters, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + Properties properties; + properties.SetId(1); + const auto process_info = ProcessInfo{}; + auto law = SaturatedLaw(); + + KRATOS_EXPECT_EQ(law.Check(properties, process_info), 0); + + properties.SetValue(SATURATED_SATURATION, 1.1); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "SATURATED_SATURATION (1.1) must be in the range [0.0, 1.0] for material 1."); + + properties.SetValue(SATURATED_SATURATION, 0.9); + KRATOS_EXPECT_EQ(law.Check(properties, process_info), 0); +} + +} // namespace Kratos::Testing \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_retention/test_van_genuchten_law.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_retention/test_van_genuchten_law.cpp new file mode 100644 index 000000000000..39081bc17a86 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_retention/test_van_genuchten_law.cpp @@ -0,0 +1,137 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Wijtze Pieter Kikstra +// + +#include "custom_retention/van_genuchten_law.h" +#include "geo_mechanics_application_variables.h" +#include "tests/cpp_tests/geo_mechanics_fast_suite.h" + +namespace Kratos::Testing +{ + +KRATOS_TEST_CASE_IN_SUITE(VanGenuchtenLawReturnsCloneOfCorrectType, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto law = VanGenuchtenLaw(); + const auto p_law_clone = law.Clone(); + KRATOS_EXPECT_NE(&law, p_law_clone.get()); + KRATOS_EXPECT_NE(dynamic_cast(p_law_clone.get()), nullptr); +} + +KRATOS_TEST_CASE_IN_SUITE(VanGenuchtenLawReturnsCalculatedValues, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law = VanGenuchtenLaw(); + Properties properties; + properties.SetValue(SATURATED_SATURATION, 0.9); + properties.SetValue(RESIDUAL_SATURATION, 0.1); + properties.SetValue(MINIMUM_RELATIVE_PERMEABILITY, 0.05); + properties.SetValue(VAN_GENUCHTEN_AIR_ENTRY_PRESSURE, 2.5); + properties.SetValue(VAN_GENUCHTEN_GN, 2.5); + properties.SetValue(VAN_GENUCHTEN_GL, 1.5); + + auto retention_law_parameters = RetentionLaw::Parameters{properties}; + + retention_law_parameters.SetFluidPressure(-10.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateSaturation(retention_law_parameters), 0.9); + double value = 0.0; + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DEGREE_OF_SATURATION, value), 0.9); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateEffectiveSaturation(retention_law_parameters), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, EFFECTIVE_SATURATION, value), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateBishopCoefficient(retention_law_parameters), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, BISHOP_COEFFICIENT, value), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateDerivativeOfSaturation(retention_law_parameters), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DERIVATIVE_OF_SATURATION, value), 0.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateRelativePermeability(retention_law_parameters), 1.0); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, RELATIVE_PERMEABILITY, value), 1.0); + + // Values below are regression values, avoiding reimplementation here of the Van Genuchten Law + retention_law_parameters.SetFluidPressure(1.5); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateSaturation(retention_law_parameters), 0.79023542376288392); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DEGREE_OF_SATURATION, value), + 0.79023542376288392); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateEffectiveSaturation(retention_law_parameters), 0.86279427970360489); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, EFFECTIVE_SATURATION, value), + 0.86279427970360489); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateBishopCoefficient(retention_law_parameters), 0.86279427970360489); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, BISHOP_COEFFICIENT, value), + 0.86279427970360489); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateDerivativeOfSaturation(retention_law_parameters), -0.15050611026881838); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, DERIVATIVE_OF_SATURATION, value), + -0.15050611026881838); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateRelativePermeability(retention_law_parameters), 0.28755984470352691); + KRATOS_EXPECT_DOUBLE_EQ(law.CalculateValue(retention_law_parameters, RELATIVE_PERMEABILITY, value), + 0.28755984470352691); +} + +KRATOS_TEST_CASE_IN_SUITE(VanGenuchtenLawChecksInputParameters, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + Properties properties; + properties.SetId(1); + const auto process_info = ProcessInfo{}; + auto law = VanGenuchtenLaw(); + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "SATURATED_SATURATION is not available in the parameters of material 1."); + properties.SetValue(SATURATED_SATURATION, 1.1); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "SATURATED_SATURATION (1.1) must be in the range [0.0, 1.0] for material 1."); + properties.SetValue(SATURATED_SATURATION, 0.9); + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "RESIDUAL_SATURATION is not available in the parameters of material 1."); + properties.SetValue(RESIDUAL_SATURATION, 1.1); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "RESIDUAL_SATURATION (1.1) must be in the range [0.0, 0.9> for material 1."); + properties.SetValue(RESIDUAL_SATURATION, 0.1); + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "MINIMUM_RELATIVE_PERMEABILITY is not available in the parameters of material 1."); + properties.SetValue(MINIMUM_RELATIVE_PERMEABILITY, 1.1); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "MINIMUM_RELATIVE_PERMEABILITY (1.1) must be in the range [0.0, 1.0] for material 1."); + properties.SetValue(MINIMUM_RELATIVE_PERMEABILITY, 0.05); + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "VAN_GENUCHTEN_AIR_ENTRY_PRESSURE is not available in the parameters of material 1."); + properties.SetValue(VAN_GENUCHTEN_AIR_ENTRY_PRESSURE, -4.0); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "VAN_GENUCHTEN_AIR_ENTRY_PRESSURE (-4) must be greater than 0 for material 1."); + properties.SetValue(VAN_GENUCHTEN_AIR_ENTRY_PRESSURE, 4.0); + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "VAN_GENUCHTEN_GN is not available in the parameters of material 1."); + properties.SetValue(VAN_GENUCHTEN_GN, -2.5); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "VAN_GENUCHTEN_GN (-2.5) must be greater than 0 for material 1."); + properties.SetValue(VAN_GENUCHTEN_GN, 2.5); + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "VAN_GENUCHTEN_GL is not available in the parameters of material 1."); + properties.SetValue(VAN_GENUCHTEN_GL, -1.5); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, process_info), + "VAN_GENUCHTEN_GL (-1.5) must be greater than or equal to 0 for material 1."); + properties.SetValue(VAN_GENUCHTEN_GL, 1.5); + + KRATOS_EXPECT_EQ(law.Check(properties, process_info), 0); +} + +} // namespace Kratos::Testing \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/tests/test_GeoMechanicsApplication.py b/applications/GeoMechanicsApplication/tests/test_GeoMechanicsApplication.py index bd80f661a99c..6920cee8d98f 100644 --- a/applications/GeoMechanicsApplication/tests/test_GeoMechanicsApplication.py +++ b/applications/GeoMechanicsApplication/tests/test_GeoMechanicsApplication.py @@ -42,6 +42,7 @@ from test_rotation_with_moving_load import KratosGeoMechanicsRotationWithMovingLoadTests from test_time_integration import KratosGeoMechanicsTimeIntegrationTests from c_phi_reduction_process import KratosGeoMechanicsCPhiReductionProcess +from test_partial_saturation import KratosGeoMechanicsPartialSaturation from test_conditions import KratosGeoMechanicsConditionTests from test_prescribed_derivatives import KratosGeoMechanicsPrescribedDerivatives @@ -124,7 +125,8 @@ def AssembleTestSuites(): KratosGeoMechanicsTransientThermalTests, KratosGeoMechanicsTimeIntegrationTests, KratosGeoMechanicsTransientPressureLineElementTests, - KratosGeoMechanicsTransientPressurePointFluxTests + KratosGeoMechanicsTransientPressurePointFluxTests, + KratosGeoMechanicsPartialSaturation ] night_test_cases.extend(small_test_cases) diff --git a/applications/GeoMechanicsApplication/tests/test_Transient_Case_A1_2D3N.gid/MaterialParameters_stage1.json b/applications/GeoMechanicsApplication/tests/test_Transient_Case_A1_2D3N.gid/MaterialParameters_stage1.json index 6c034d079b48..32295e74fc15 100644 --- a/applications/GeoMechanicsApplication/tests/test_Transient_Case_A1_2D3N.gid/MaterialParameters_stage1.json +++ b/applications/GeoMechanicsApplication/tests/test_Transient_Case_A1_2D3N.gid/MaterialParameters_stage1.json @@ -21,13 +21,8 @@ "DYNAMIC_VISCOSITY" : 9.81, "THICKNESS" : 1.0, "BIOT_COEFFICIENT" : 1.0, - "RETENTION_LAW" : "SaturatedLaw", - "SATURATED_SATURATION" : 1.0, - "RESIDUAL_SATURATION" : 0.06203, - "VAN_GENUCHTEN_AIR_ENTRY_PRESSURE" : 4.379464286, - "VAN_GENUCHTEN_GN" : 2.286, - "VAN_GENUCHTEN_GL" : 0, - "MINIMUM_RELATIVE_PERMEABILITY" : 0.0001 + "RETENTION_LAW" : "SaturatedLaw", + "SATURATED_SATURATION" : 1.0 }, "Tables": {} } diff --git a/applications/GeoMechanicsApplication/tests/test_Transient_Case_A1_2D6N.gid/MaterialParameters_stage1.json b/applications/GeoMechanicsApplication/tests/test_Transient_Case_A1_2D6N.gid/MaterialParameters_stage1.json index 6c034d079b48..32295e74fc15 100644 --- a/applications/GeoMechanicsApplication/tests/test_Transient_Case_A1_2D6N.gid/MaterialParameters_stage1.json +++ b/applications/GeoMechanicsApplication/tests/test_Transient_Case_A1_2D6N.gid/MaterialParameters_stage1.json @@ -21,13 +21,8 @@ "DYNAMIC_VISCOSITY" : 9.81, "THICKNESS" : 1.0, "BIOT_COEFFICIENT" : 1.0, - "RETENTION_LAW" : "SaturatedLaw", - "SATURATED_SATURATION" : 1.0, - "RESIDUAL_SATURATION" : 0.06203, - "VAN_GENUCHTEN_AIR_ENTRY_PRESSURE" : 4.379464286, - "VAN_GENUCHTEN_GN" : 2.286, - "VAN_GENUCHTEN_GL" : 0, - "MINIMUM_RELATIVE_PERMEABILITY" : 0.0001 + "RETENTION_LAW" : "SaturatedLaw", + "SATURATED_SATURATION" : 1.0 }, "Tables": {} } diff --git a/applications/GeoMechanicsApplication/tests/test_partial_saturation.py b/applications/GeoMechanicsApplication/tests/test_partial_saturation.py new file mode 100644 index 000000000000..a8b736f4f22c --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partial_saturation.py @@ -0,0 +1,79 @@ +import os + +import KratosMultiphysics as Kratos +import KratosMultiphysics.GeoMechanicsApplication.geomechanics_analysis as analysis + +import KratosMultiphysics.KratosUnittest as KratosUnittest + +import test_helper + +class KratosGeoMechanicsPartialSaturation(KratosUnittest.TestCase): + """ + This class contains benchmark tests which are checked with the analytical solution + """ + + def __test_saturated_below_phreatic_level_pw(self, test_name): + n_stages = 2 + + # get the parameter file names for all stages + file_path = test_helper.get_file_path(os.path.join('test_partially_saturated', test_name)) + parameter_file_names = [os.path.join(file_path, f'ProjectParameters_stage{i+1}.json') for i in + range(n_stages)] + + # set stage parameters + parameters_stages = [] + os.chdir(file_path) + for parameter_file_name in parameter_file_names: + with open(parameter_file_name, 'r') as parameter_file: + parameters_stages.append(Kratos.Parameters(parameter_file.read())) + + model = Kratos.Model() + + # run stages and get water pressure/displacement results per stage + stage_water_pressure = [] + coords = [] + for stage_parameters in parameters_stages: + stage = analysis.GeoMechanicsAnalysis(model, stage_parameters) + stage.Run() + stage_water_pressure.append(test_helper.get_water_pressure(stage)) + coords.append(test_helper.get_nodal_coordinates(stage)) + + # get y coords of all the nodes + y_coords = [coord[1] for coord in coords[0]] + + # calculate water pressure analytical solution for all stages and calculate the error + rel_p_stage = [self.__compute_hydrostatic_water_pressure(y_coord) for y_coord in y_coords] + + errors_stage = [actual_pressure - expected_pressure for actual_pressure, expected_pressure in + zip(stage_water_pressure[1], rel_p_stage)] + rmse_stages = (sum([error ** 2 for error in errors_stage]) / len(errors_stage)) ** 0.5 + + # assert if average error in all stages is below accuracy + accuracy = 1.0e-7 + self.assertLess(rmse_stages, accuracy) + + def test_saturated_below_phreatic_level_pw_triangle3N(self): + self.__test_saturated_below_phreatic_level_pw('test_saturated_below_phreatic_level_pw_triangle3N') + + def test_saturated_below_phreatic_level_pw_triangle6N(self): + self.__test_saturated_below_phreatic_level_pw('test_saturated_below_phreatic_level_pw_triangle6N') + + def test_saturated_below_phreatic_level_upw_difforder_triangle6n(self): + self.__test_saturated_below_phreatic_level_pw('test_saturated_below_phreatic_level_upw_difforder_triangle6n') + + def test_saturated_below_phreatic_level_upw_smallstrain_triangle3n(self): + self.__test_saturated_below_phreatic_level_pw('test_saturated_below_phreatic_level_upw_smallstrain_triangle3n') + + def test_saturated_below_phreatic_level_upw_smallstrain_triangle6n(self): + self.__test_saturated_below_phreatic_level_pw('test_saturated_below_phreatic_level_upw_smallstrain_triangle6n') + + def __compute_hydrostatic_water_pressure(self, y_coord): + water_density = 1019.367991845056 + gravity = -9.81 + phreatic_level = -2.0 + result = gravity * water_density * (phreatic_level - y_coord) + return min([result, 0.0]) + + +if __name__ == '__main__': + KratosUnittest.main() diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/Common/MaterialParameters_stage1.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/Common/MaterialParameters_stage1.json new file mode 100644 index 000000000000..fc2e750c8e3b --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/Common/MaterialParameters_stage1.json @@ -0,0 +1,32 @@ +{ + "properties": [ + { + "properties_id": 0, + "model_part_name": "PorousDomain.Parts_Solid_Grond", + "Material": { + "constitutive_law": { + "name": "GeoLinearElasticPlaneStrain2DLaw" + }, + "Variables": { + "IGNORE_UNDRAINED": true, + "BIOT_COEFFICIENT": 1.0, + "DENSITY_SOLID": 2038.735983690112, + "DENSITY_WATER": 1019.367991845056, + "POROSITY": 0.3, + "BULK_MODULUS_SOLID": 20000000000.0, + "BULK_MODULUS_FLUID": 2200000000.0, + "PERMEABILITY_XX": 1.157E-05, + "PERMEABILITY_YY": 1.157E-05, + "PERMEABILITY_XY": 0.0, + "DYNAMIC_VISCOSITY": 0.001, + "RETENTION_LAW": "SaturatedBelowPhreaticLevelLaw", + "RESIDUAL_SATURATION": 0.001, + "SATURATED_SATURATION": 1.0, + "MINIMUM_RELATIVE_PERMEABILITY": 0.001, + "YOUNG_MODULUS": 1000000000.0, + "POISSON_RATIO": 0.2 + } + } + } + ] +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/Common/MaterialParameters_stage2.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/Common/MaterialParameters_stage2.json new file mode 100644 index 000000000000..a59003318063 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/Common/MaterialParameters_stage2.json @@ -0,0 +1,32 @@ +{ + "properties": [ + { + "properties_id": 0, + "model_part_name": "PorousDomain.Parts_Solid_Grond", + "Material": { + "constitutive_law": { + "name": "GeoLinearElasticPlaneStrain2DLaw" + }, + "Variables": { + "IGNORE_UNDRAINED": false, + "BIOT_COEFFICIENT": 1.0, + "DENSITY_SOLID": 2038.735983690112, + "DENSITY_WATER": 1019.367991845056, + "POROSITY": 0.3, + "BULK_MODULUS_SOLID": 20000000000.0, + "BULK_MODULUS_FLUID": 2200000000.0, + "PERMEABILITY_XX": 1.157E-05, + "PERMEABILITY_YY": 1.157E-05, + "PERMEABILITY_XY": 0.0, + "DYNAMIC_VISCOSITY": 0.001, + "RETENTION_LAW": "SaturatedBelowPhreaticLevelLaw", + "RESIDUAL_SATURATION": 0.001, + "SATURATED_SATURATION": 1.0, + "MINIMUM_RELATIVE_PERMEABILITY": 0.001, + "YOUNG_MODULUS": 1000000000.0, + "POISSON_RATIO": 0.2 + } + } + } + ] +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/README.md b/applications/GeoMechanicsApplication/tests/test_partially_saturated/README.md new file mode 100644 index 000000000000..3ad042bfef70 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/README.md @@ -0,0 +1,21 @@ +# Test Cases for partially saturated flow +## Saturated below phreatic level + +**Author:** [Mohamed Nabi](https://github.com/mnabideltares) + +**Source files:** [Partially saturated flow](https://github.com/KratosMultiphysics/Kratos/tree/master/applications/GeoMechanicsApplication/tests/test_partially_saturated) + +## Case Specification +In this test case, a column of 1 x 5 m soil is considered (set between $`y = -5 \mathrm{m}`$ and $`y = 0 \mathrm{m}`$). A phreatic line is set at the level of $y = -2$ m. The water pore pressure field is then calculated. +The simulation is done on a double stage process, and steady state Pw and U-Pw elements are considered. This test is conducted for various configurations, including + +- 2D3N: for Pw and U-Pw small strain elements +- 2D6N: for Pw and U-Pw small strain and U-Pw different order elements. + +The pressure distribution along the column is then compared with its analytical result. + +## Results + +The picture below illustrates the pressure contours resulting from the simulation (as an example the 2D6N for Pw is shown below). + +Pressure field for case of saturation below phreatic level at stage 2 \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/documentation_data/test_saturated_below_phreatic_level_pw_triangle6n_results.png b/applications/GeoMechanicsApplication/tests/test_partially_saturated/documentation_data/test_saturated_below_phreatic_level_pw_triangle6n_results.png new file mode 100644 index 000000000000..00bac40d18d9 Binary files /dev/null and b/applications/GeoMechanicsApplication/tests/test_partially_saturated/documentation_data/test_saturated_below_phreatic_level_pw_triangle6n_results.png differ diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle3N/ProjectParameters_stage1.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle3N/ProjectParameters_stage1.json new file mode 100644 index 000000000000..077a19ee822d --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle3N/ProjectParameters_stage1.json @@ -0,0 +1,179 @@ +{ + "problem_data": { + "problem_name" : "test_saturated_below_phreatic_level_pw_triangle3N_stage1", + "start_time" : -0.1, + "end_time" : 0.0, + "echo_level" : 1, + "parallel_type" : "OpenMP", + "number_of_threads" : 1 + }, + "solver_settings": { + "solver_type" : "Pw", + "model_part_name" : "PorousDomain", + "domain_size" : 2, + "model_import_settings": { + "input_type" : "mdpa", + "input_filename" : "column" + }, + "material_import_settings": { + "materials_filename" : "../Common/MaterialParameters_stage1.json" + }, + "time_stepping": { + "time_step" : 0.1, + "max_delta_time_factor" : 1.0 + }, + "buffer_size" : 2, + "echo_level" : 1, + "clear_storage" : false, + "compute_reactions" : true, + "reform_dofs_at_each_step" : false, + "nodal_smoothing" : false, + "block_builder" : true, + "solution_type" : "Steady_State_Groundwater_Flow", + "scheme_type" : "Backward_Euler", + "reset_displacements" : true, + "strategy_type" : "newton_raphson", + "convergence_criterion" : "water_pressure_criterion", + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1E-09, + "displacement_relative_tolerance" : 0.001, + "displacement_absolute_tolerance" : 1E-09, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1E-09, + "min_iterations" : 6, + "max_iterations" : 30, + "number_cycles" : 1, + "reduction_factor" : 1.0, + "increase_factor" : 1.0, + "desired_iterations" : 4, + "max_radius_factor" : 10.0, + "min_radius_factor" : 0.1, + "calculate_reactions" : true, + "max_line_search_iterations" : 5, + "first_alpha_value" : 0.5, + "second_alpha_value" : 1.0, + "min_alpha" : 0.1, + "max_alpha" : 2.0, + "line_search_tolerance" : 0.5, + "rotation_dofs" : false, + "linear_solver_settings": { + "solver_type" : "LinearSolversApplication.sparse_lu", + "scaling" : true + }, + "problem_domain_sub_model_part_list": ["Parts_Solid_Grond"], + "processes_sub_model_part_list" : ["gravity", "WATER_PRESSURE_Initial", "WATER_PRESSURE_Bottom", "WATER_PRESSURE_Top"], + "body_domain_sub_model_part_list" : ["Parts_Solid_Grond"], + "newmark_beta" : 0.25, + "newmark_gamma" : 0.5, + "newmark_theta" : 0.5, + "rayleigh_m" : 0.0, + "rayleigh_k" : 0.0, + "move_mesh_flag" : false + }, + "output_processes": { + "gid_output": [ + { + "python_module" : "gid_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "GiDOutputProcess", + "Parameters": { + "model_part_name" : "PorousDomain.porous_computational_model_part", + "output_name" : "stage1", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag" : "WriteElementsOnly", + "GiDPostMode" : "GiD_PostAscii", + "MultiFileFlag" : "SingleFile" + }, + "file_label" : "step", + "output_control_type" : "step", + "output_interval" : 1, + "body_output" : true, + "node_output" : false, + "skin_output" : false, + "plane_output" : [], + "nodal_results" : ["WATER_PRESSURE"], + "gauss_point_results" : [] + }, + "point_data_configuration" : [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Initial", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : false, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Bottom", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Top", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + } + ], + "loads_process_list": [ + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "VOLUME_ACCELERATION", + "active" : [false, true, false], + "value" : [0.0, -9.81, 0.0], + "table" : [0, 0, 0] + } + } + ], + "auxiliar_process_list": [] + } +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle3N/ProjectParameters_stage2.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle3N/ProjectParameters_stage2.json new file mode 100644 index 000000000000..676c1c307769 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle3N/ProjectParameters_stage2.json @@ -0,0 +1,161 @@ +{ + "problem_data": { + "problem_name" : "test_saturated_below_phreatic_level_pw_triangle3N_stage2", + "start_time" : 0.0, + "end_time" : 1.0, + "echo_level" : 1, + "parallel_type" : "OpenMP", + "number_of_threads" : 1 + }, + "solver_settings": { + "solver_type" : "Pw", + "model_part_name" : "PorousDomain", + "domain_size" : 2, + "model_import_settings": { + "input_type" : "mdpa", + "input_filename" : "column" + }, + "material_import_settings": { + "materials_filename" : "../Common/MaterialParameters_stage2.json" + }, + "time_stepping": { + "time_step" : 0.1, + "max_delta_time_factor" : 1.0 + }, + "buffer_size" : 2, + "echo_level" : 1, + "clear_storage" : false, + "compute_reactions" : true, + "reform_dofs_at_each_step" : false, + "nodal_smoothing" : false, + "block_builder" : true, + "solution_type" : "Steady_State_Groundwater_Flow", + "scheme_type" : "Backward_Euler", + "reset_displacements" : false, + "strategy_type" : "newton_raphson", + "convergence_criterion" : "water_pressure_criterion", + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1E-09, + "displacement_relative_tolerance" : 0.0001, + "displacement_absolute_tolerance" : 1E-09, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1E-09, + "min_iterations" : 6, + "max_iterations" : 30, + "number_cycles" : 1, + "reduction_factor" : 1.0, + "increase_factor" : 1.0, + "desired_iterations" : 4, + "max_radius_factor" : 10.0, + "min_radius_factor" : 0.1, + "calculate_reactions" : true, + "max_line_search_iterations" : 5, + "first_alpha_value" : 0.5, + "second_alpha_value" : 1.0, + "min_alpha" : 0.1, + "max_alpha" : 2.0, + "line_search_tolerance" : 0.5, + "rotation_dofs" : false, + "linear_solver_settings": { + "solver_type" : "LinearSolversApplication.sparse_lu", + "scaling" : true + }, + "problem_domain_sub_model_part_list": ["Parts_Solid_Grond"], + "processes_sub_model_part_list" : ["gravity", "WATER_PRESSURE_Bottom", "WATER_PRESSURE_Top"], + "body_domain_sub_model_part_list" : ["Parts_Solid_Grond"], + "newmark_beta" : 0.25, + "newmark_gamma" : 0.5, + "newmark_theta" : 0.5, + "rayleigh_m" : 0.0, + "rayleigh_k" : 0.0, + "move_mesh_flag" : false + }, + "output_processes": { + "gid_output": [ + { + "python_module" : "gid_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "GiDOutputProcess", + "Parameters": { + "model_part_name" : "PorousDomain.porous_computational_model_part", + "output_name" : "stage2", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag" : "WriteElementsOnly", + "GiDPostMode" : "GiD_PostAscii", + "MultiFileFlag" : "SingleFile" + }, + "file_label" : "step", + "output_control_type" : "step", + "output_interval" : 1, + "body_output" : true, + "node_output" : false, + "skin_output" : false, + "plane_output" : [], + "nodal_results" : ["WATER_PRESSURE"], + "gauss_point_results" : [] + }, + "point_data_configuration" : [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Bottom", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Top", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + } + ], + "loads_process_list": [ + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "VOLUME_ACCELERATION", + "active" : [false, true, false], + "value" : [0.0, -9.81, 0.0], + "table" : [0, 0, 0] + } + } + ], + "auxiliar_process_list": [] + } +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle3N/column.mdpa b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle3N/column.mdpa new file mode 100644 index 000000000000..a2d7c31ff6c1 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle3N/column.mdpa @@ -0,0 +1,172 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties + +Begin Nodes + 1 0.0000000000 -5.0000000000 0.0000000000 + 2 1.0000000000 -5.0000000000 0.0000000000 + 3 0.0000000000 -4.5000000000 0.0000000000 + 4 1.0000000000 -4.5000000000 0.0000000000 + 5 0.0000000000 -4.0000000000 0.0000000000 + 6 1.0000000000 -4.0000000000 0.0000000000 + 7 0.0000000000 -3.5000000000 0.0000000000 + 8 1.0000000000 -3.5000000000 0.0000000000 + 9 0.0000000000 -3.0000000000 0.0000000000 + 10 1.0000000000 -3.0000000000 0.0000000000 + 11 0.0000000000 -2.5000000000 0.0000000000 + 12 1.0000000000 -2.5000000000 0.0000000000 + 13 0.0000000000 -2.0000000000 0.0000000000 + 14 1.0000000000 -2.0000000000 0.0000000000 + 15 0.0000000000 -1.5000000000 0.0000000000 + 16 1.0000000000 -1.5000000000 0.0000000000 + 17 0.0000000000 -1.0000000000 0.0000000000 + 18 1.0000000000 -1.0000000000 0.0000000000 + 19 0.0000000000 -0.5000000000 0.0000000000 + 20 1.0000000000 -0.5000000000 0.0000000000 + 21 0.0000000000 0.0000000000 0.0000000000 + 22 1.0000000000 0.0000000000 0.0000000000 +End Nodes + + +Begin Elements SteadyStatePwElement2D3N // GUI group identifier: Grond + 1 0 1 2 4 + 2 0 4 3 1 + 3 0 3 4 6 + 4 0 6 5 3 + 5 0 5 6 8 + 6 0 8 7 5 + 7 0 7 8 10 + 8 0 10 9 7 + 9 0 9 10 12 + 10 0 12 11 9 + 11 0 11 12 14 + 12 0 14 13 11 + 13 0 13 14 16 + 14 0 16 15 13 + 15 0 15 16 18 + 16 0 18 17 15 + 17 0 17 18 20 + 18 0 20 19 17 + 19 0 19 20 22 + 20 0 22 21 19 +End Elements + +Begin SubModelPart Parts_Solid_Grond // Group Grond // Subtree Parts_Solid + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + End SubModelPartElements +End SubModelPart + +Begin SubModelPart gravity + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Initial + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Top + Begin SubModelPartNodes + 21 + 22 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Bottom + Begin SubModelPartNodes + 1 + 2 + End SubModelPartNodes +End SubModelPart diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle6N/ProjectParameters_stage1.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle6N/ProjectParameters_stage1.json new file mode 100644 index 000000000000..077a19ee822d --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle6N/ProjectParameters_stage1.json @@ -0,0 +1,179 @@ +{ + "problem_data": { + "problem_name" : "test_saturated_below_phreatic_level_pw_triangle3N_stage1", + "start_time" : -0.1, + "end_time" : 0.0, + "echo_level" : 1, + "parallel_type" : "OpenMP", + "number_of_threads" : 1 + }, + "solver_settings": { + "solver_type" : "Pw", + "model_part_name" : "PorousDomain", + "domain_size" : 2, + "model_import_settings": { + "input_type" : "mdpa", + "input_filename" : "column" + }, + "material_import_settings": { + "materials_filename" : "../Common/MaterialParameters_stage1.json" + }, + "time_stepping": { + "time_step" : 0.1, + "max_delta_time_factor" : 1.0 + }, + "buffer_size" : 2, + "echo_level" : 1, + "clear_storage" : false, + "compute_reactions" : true, + "reform_dofs_at_each_step" : false, + "nodal_smoothing" : false, + "block_builder" : true, + "solution_type" : "Steady_State_Groundwater_Flow", + "scheme_type" : "Backward_Euler", + "reset_displacements" : true, + "strategy_type" : "newton_raphson", + "convergence_criterion" : "water_pressure_criterion", + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1E-09, + "displacement_relative_tolerance" : 0.001, + "displacement_absolute_tolerance" : 1E-09, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1E-09, + "min_iterations" : 6, + "max_iterations" : 30, + "number_cycles" : 1, + "reduction_factor" : 1.0, + "increase_factor" : 1.0, + "desired_iterations" : 4, + "max_radius_factor" : 10.0, + "min_radius_factor" : 0.1, + "calculate_reactions" : true, + "max_line_search_iterations" : 5, + "first_alpha_value" : 0.5, + "second_alpha_value" : 1.0, + "min_alpha" : 0.1, + "max_alpha" : 2.0, + "line_search_tolerance" : 0.5, + "rotation_dofs" : false, + "linear_solver_settings": { + "solver_type" : "LinearSolversApplication.sparse_lu", + "scaling" : true + }, + "problem_domain_sub_model_part_list": ["Parts_Solid_Grond"], + "processes_sub_model_part_list" : ["gravity", "WATER_PRESSURE_Initial", "WATER_PRESSURE_Bottom", "WATER_PRESSURE_Top"], + "body_domain_sub_model_part_list" : ["Parts_Solid_Grond"], + "newmark_beta" : 0.25, + "newmark_gamma" : 0.5, + "newmark_theta" : 0.5, + "rayleigh_m" : 0.0, + "rayleigh_k" : 0.0, + "move_mesh_flag" : false + }, + "output_processes": { + "gid_output": [ + { + "python_module" : "gid_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "GiDOutputProcess", + "Parameters": { + "model_part_name" : "PorousDomain.porous_computational_model_part", + "output_name" : "stage1", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag" : "WriteElementsOnly", + "GiDPostMode" : "GiD_PostAscii", + "MultiFileFlag" : "SingleFile" + }, + "file_label" : "step", + "output_control_type" : "step", + "output_interval" : 1, + "body_output" : true, + "node_output" : false, + "skin_output" : false, + "plane_output" : [], + "nodal_results" : ["WATER_PRESSURE"], + "gauss_point_results" : [] + }, + "point_data_configuration" : [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Initial", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : false, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Bottom", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Top", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + } + ], + "loads_process_list": [ + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "VOLUME_ACCELERATION", + "active" : [false, true, false], + "value" : [0.0, -9.81, 0.0], + "table" : [0, 0, 0] + } + } + ], + "auxiliar_process_list": [] + } +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle6N/ProjectParameters_stage2.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle6N/ProjectParameters_stage2.json new file mode 100644 index 000000000000..676c1c307769 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle6N/ProjectParameters_stage2.json @@ -0,0 +1,161 @@ +{ + "problem_data": { + "problem_name" : "test_saturated_below_phreatic_level_pw_triangle3N_stage2", + "start_time" : 0.0, + "end_time" : 1.0, + "echo_level" : 1, + "parallel_type" : "OpenMP", + "number_of_threads" : 1 + }, + "solver_settings": { + "solver_type" : "Pw", + "model_part_name" : "PorousDomain", + "domain_size" : 2, + "model_import_settings": { + "input_type" : "mdpa", + "input_filename" : "column" + }, + "material_import_settings": { + "materials_filename" : "../Common/MaterialParameters_stage2.json" + }, + "time_stepping": { + "time_step" : 0.1, + "max_delta_time_factor" : 1.0 + }, + "buffer_size" : 2, + "echo_level" : 1, + "clear_storage" : false, + "compute_reactions" : true, + "reform_dofs_at_each_step" : false, + "nodal_smoothing" : false, + "block_builder" : true, + "solution_type" : "Steady_State_Groundwater_Flow", + "scheme_type" : "Backward_Euler", + "reset_displacements" : false, + "strategy_type" : "newton_raphson", + "convergence_criterion" : "water_pressure_criterion", + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1E-09, + "displacement_relative_tolerance" : 0.0001, + "displacement_absolute_tolerance" : 1E-09, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1E-09, + "min_iterations" : 6, + "max_iterations" : 30, + "number_cycles" : 1, + "reduction_factor" : 1.0, + "increase_factor" : 1.0, + "desired_iterations" : 4, + "max_radius_factor" : 10.0, + "min_radius_factor" : 0.1, + "calculate_reactions" : true, + "max_line_search_iterations" : 5, + "first_alpha_value" : 0.5, + "second_alpha_value" : 1.0, + "min_alpha" : 0.1, + "max_alpha" : 2.0, + "line_search_tolerance" : 0.5, + "rotation_dofs" : false, + "linear_solver_settings": { + "solver_type" : "LinearSolversApplication.sparse_lu", + "scaling" : true + }, + "problem_domain_sub_model_part_list": ["Parts_Solid_Grond"], + "processes_sub_model_part_list" : ["gravity", "WATER_PRESSURE_Bottom", "WATER_PRESSURE_Top"], + "body_domain_sub_model_part_list" : ["Parts_Solid_Grond"], + "newmark_beta" : 0.25, + "newmark_gamma" : 0.5, + "newmark_theta" : 0.5, + "rayleigh_m" : 0.0, + "rayleigh_k" : 0.0, + "move_mesh_flag" : false + }, + "output_processes": { + "gid_output": [ + { + "python_module" : "gid_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "GiDOutputProcess", + "Parameters": { + "model_part_name" : "PorousDomain.porous_computational_model_part", + "output_name" : "stage2", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag" : "WriteElementsOnly", + "GiDPostMode" : "GiD_PostAscii", + "MultiFileFlag" : "SingleFile" + }, + "file_label" : "step", + "output_control_type" : "step", + "output_interval" : 1, + "body_output" : true, + "node_output" : false, + "skin_output" : false, + "plane_output" : [], + "nodal_results" : ["WATER_PRESSURE"], + "gauss_point_results" : [] + }, + "point_data_configuration" : [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Bottom", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Top", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + } + ], + "loads_process_list": [ + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "VOLUME_ACCELERATION", + "active" : [false, true, false], + "value" : [0.0, -9.81, 0.0], + "table" : [0, 0, 0] + } + } + ], + "auxiliar_process_list": [] + } +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle6N/column.mdpa b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle6N/column.mdpa new file mode 100644 index 000000000000..f5420a424359 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_pw_triangle6N/column.mdpa @@ -0,0 +1,338 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties + +Begin Nodes + 1 0.0000000000 -5.0000000000 0.0000000000 + 2 1.0000000000 -5.0000000000 0.0000000000 + 3 0.0000000000 -4.5000000000 0.0000000000 + 4 1.0000000000 -4.5000000000 0.0000000000 + 5 0.0000000000 -4.0000000000 0.0000000000 + 6 1.0000000000 -4.0000000000 0.0000000000 + 7 0.0000000000 -3.5000000000 0.0000000000 + 8 1.0000000000 -3.5000000000 0.0000000000 + 9 0.0000000000 -3.0000000000 0.0000000000 + 10 1.0000000000 -3.0000000000 0.0000000000 + 11 0.0000000000 -2.5000000000 0.0000000000 + 12 1.0000000000 -2.5000000000 0.0000000000 + 13 0.0000000000 -2.0000000000 0.0000000000 + 14 1.0000000000 -2.0000000000 0.0000000000 + 15 0.0000000000 -1.5000000000 0.0000000000 + 16 1.0000000000 -1.5000000000 0.0000000000 + 17 0.0000000000 -1.0000000000 0.0000000000 + 18 1.0000000000 -1.0000000000 0.0000000000 + 19 0.0000000000 -0.5000000000 0.0000000000 + 20 1.0000000000 -0.5000000000 0.0000000000 + 21 0.0000000000 0.0000000000 0.0000000000 + 22 1.0000000000 0.0000000000 0.0000000000 + 23 0.0000000000 -4.7500000000 0.0000000000 + 24 0.0000000000 -4.2500000000 0.0000000000 + 25 0.0000000000 -3.7500000000 0.0000000000 + 26 0.0000000000 -3.2500000000 0.0000000000 + 27 0.0000000000 -2.7500000000 0.0000000000 + 28 0.0000000000 -2.2500000000 0.0000000000 + 29 0.0000000000 -1.7500000000 0.0000000000 + 30 0.0000000000 -1.2500000000 0.0000000000 + 31 0.0000000000 -0.7500000000 0.0000000000 + 32 0.0000000000 -0.2500000000 0.0000000000 + 33 1.0000000000 -4.7500000000 0.0000000000 + 34 1.0000000000 -4.2500000000 0.0000000000 + 35 1.0000000000 -3.7500000000 0.0000000000 + 36 1.0000000000 -3.2500000000 0.0000000000 + 37 1.0000000000 -2.7500000000 0.0000000000 + 38 1.0000000000 -2.2500000000 0.0000000000 + 39 1.0000000000 -1.7500000000 0.0000000000 + 40 1.0000000000 -1.2500000000 0.0000000000 + 41 1.0000000000 -0.7500000000 0.0000000000 + 42 1.0000000000 -0.2500000000 0.0000000000 + 43 0.5000000000 -5.0000000000 0.0000000000 + 44 0.5000000000 -4.7500000000 0.0000000000 + 45 0.5000000000 -4.5000000000 0.0000000000 + 46 0.5000000000 -4.2500000000 0.0000000000 + 47 0.5000000000 -4.0000000000 0.0000000000 + 48 0.5000000000 -3.7500000000 0.0000000000 + 49 0.5000000000 -3.5000000000 0.0000000000 + 50 0.5000000000 -3.2500000000 0.0000000000 + 51 0.5000000000 -3.0000000000 0.0000000000 + 52 0.5000000000 -2.7500000000 0.0000000000 + 53 0.5000000000 -2.5000000000 0.0000000000 + 54 0.5000000000 -2.2500000000 0.0000000000 + 55 0.5000000000 -2.0000000000 0.0000000000 + 56 0.5000000000 -1.7500000000 0.0000000000 + 57 0.5000000000 -1.5000000000 0.0000000000 + 58 0.5000000000 -1.2500000000 0.0000000000 + 59 0.5000000000 -1.0000000000 0.0000000000 + 60 0.5000000000 -0.7500000000 0.0000000000 + 61 0.5000000000 -0.5000000000 0.0000000000 + 62 0.5000000000 -0.2500000000 0.0000000000 + 63 0.5000000000 0.0000000000 0.0000000000 +End Nodes + + +Begin Elements SteadyStatePwElement2D6N // GUI group identifier: Grond + 1 0 1 2 4 43 33 44 + 2 0 4 3 1 45 23 44 + 3 0 3 4 6 45 34 46 + 4 0 6 5 3 47 24 46 + 5 0 5 6 8 47 35 48 + 6 0 8 7 5 49 25 48 + 7 0 7 8 10 49 36 50 + 8 0 10 9 7 51 26 50 + 9 0 9 10 12 51 37 52 + 10 0 12 11 9 53 27 52 + 11 0 11 12 14 53 38 54 + 12 0 14 13 11 55 28 54 + 13 0 13 14 16 55 39 56 + 14 0 16 15 13 57 29 56 + 15 0 15 16 18 57 40 58 + 16 0 18 17 15 59 30 58 + 17 0 17 18 20 59 41 60 + 18 0 20 19 17 61 31 60 + 19 0 19 20 22 61 42 62 + 20 0 22 21 19 63 32 62 +End Elements + +Begin SubModelPart Parts_Solid_Grond // Group Grond // Subtree Parts_Solid + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + 54 + 55 + 56 + 57 + 58 + 59 + 60 + 61 + 62 + 63 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + End SubModelPartElements +End SubModelPart + +Begin SubModelPart gravity + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + 54 + 55 + 56 + 57 + 58 + 59 + 60 + 61 + 62 + 63 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Initial + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + 54 + 55 + 56 + 57 + 58 + 59 + 60 + 61 + 62 + 63 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Top + Begin SubModelPartNodes + 21 + 22 + 63 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Bottom + Begin SubModelPartNodes + 1 + 2 + 43 + End SubModelPartNodes +End SubModelPart diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_difforder_triangle6n/ProjectParameters_stage1.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_difforder_triangle6n/ProjectParameters_stage1.json new file mode 100644 index 000000000000..00f7bcaf6c87 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_difforder_triangle6n/ProjectParameters_stage1.json @@ -0,0 +1,192 @@ +{ + "problem_data": { + "problem_name" : "test_saturated_below_phreatic_level_upw_difforder_triangle6n_stage1", + "start_time" : -0.1, + "end_time" : 0.0, + "echo_level" : 1, + "parallel_type" : "OpenMP", + "number_of_threads" : 1 + }, + "solver_settings": { + "solver_type" : "U_Pw", + "model_part_name" : "PorousDomain", + "domain_size" : 2, + "model_import_settings": { + "input_type" : "mdpa", + "input_filename" : "column" + }, + "material_import_settings": { + "materials_filename" : "../Common/MaterialParameters_stage1.json" + }, + "time_stepping": { + "time_step" : 0.1, + "max_delta_time_factor" : 1.0 + }, + "buffer_size" : 2, + "echo_level" : 1, + "clear_storage" : false, + "compute_reactions" : true, + "reform_dofs_at_each_step" : false, + "nodal_smoothing" : false, + "block_builder" : true, + "solution_type" : "Quasi-Static", + "scheme_type" : "Newmark", + "reset_displacements" : true, + "strategy_type" : "newton_raphson", + "convergence_criterion" : "water_pressure_criterion", + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1E-09, + "displacement_relative_tolerance" : 0.001, + "displacement_absolute_tolerance" : 1E-09, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1E-09, + "min_iterations" : 6, + "max_iterations" : 30, + "number_cycles" : 1, + "reduction_factor" : 1.0, + "increase_factor" : 1.0, + "desired_iterations" : 4, + "max_radius_factor" : 10.0, + "min_radius_factor" : 0.1, + "calculate_reactions" : true, + "max_line_search_iterations" : 5, + "first_alpha_value" : 0.5, + "second_alpha_value" : 1.0, + "min_alpha" : 0.1, + "max_alpha" : 2.0, + "line_search_tolerance" : 0.5, + "rotation_dofs" : false, + "linear_solver_settings": { + "solver_type" : "LinearSolversApplication.sparse_lu", + "scaling" : true + }, + "problem_domain_sub_model_part_list": ["Parts_Solid_Grond"], + "processes_sub_model_part_list" : ["gravity", "WATER_PRESSURE_Initial", "WATER_PRESSURE_Bottom", "WATER_PRESSURE_Top"], + "body_domain_sub_model_part_list" : ["Parts_Solid_Grond"], + "newmark_beta" : 0.25, + "newmark_gamma" : 0.5, + "newmark_theta" : 0.5, + "rayleigh_m" : 0.0, + "rayleigh_k" : 0.0, + "move_mesh_flag" : false + }, + "output_processes": { + "gid_output": [ + { + "python_module" : "gid_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "GiDOutputProcess", + "Parameters": { + "model_part_name" : "PorousDomain.porous_computational_model_part", + "output_name" : "stage1", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag" : "WriteElementsOnly", + "GiDPostMode" : "GiD_PostAscii", + "MultiFileFlag" : "SingleFile" + }, + "file_label" : "step", + "output_control_type" : "step", + "output_interval" : 1, + "body_output" : true, + "node_output" : false, + "skin_output" : false, + "plane_output" : [], + "nodal_results" : ["WATER_PRESSURE"], + "gauss_point_results" : [] + }, + "point_data_configuration" : [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Initial", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : false, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Bottom", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Top", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "DISPLACEMENT", + "active" : [true, true, true], + "is_fixed" : [true, true, true], + "value" : [0.0, 0.0, 0.0], + "table" : [0, 0, 0] + } + } + ], + "loads_process_list": [ + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "VOLUME_ACCELERATION", + "active" : [false, true, false], + "value" : [0.0, -9.81, 0.0], + "table" : [0, 0, 0] + } + } + ], + "auxiliar_process_list": [] + } +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_difforder_triangle6n/ProjectParameters_stage2.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_difforder_triangle6n/ProjectParameters_stage2.json new file mode 100644 index 000000000000..89eac7fb285f --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_difforder_triangle6n/ProjectParameters_stage2.json @@ -0,0 +1,174 @@ +{ + "problem_data": { + "problem_name" : "test_saturated_below_phreatic_level_upw_difforder_triangle6n_stage2", + "start_time" : 0.0, + "end_time" : 1.0, + "echo_level" : 1, + "parallel_type" : "OpenMP", + "number_of_threads" : 1 + }, + "solver_settings": { + "solver_type" : "U_Pw", + "model_part_name" : "PorousDomain", + "domain_size" : 2, + "model_import_settings": { + "input_type" : "mdpa", + "input_filename" : "column" + }, + "material_import_settings": { + "materials_filename" : "../Common/MaterialParameters_stage2.json" + }, + "time_stepping": { + "time_step" : 0.1, + "max_delta_time_factor" : 1.0 + }, + "buffer_size" : 2, + "echo_level" : 1, + "clear_storage" : false, + "compute_reactions" : true, + "reform_dofs_at_each_step" : false, + "nodal_smoothing" : false, + "block_builder" : true, + "solution_type" : "Quasi-Static", + "scheme_type" : "Newmark", + "reset_displacements" : false, + "strategy_type" : "newton_raphson", + "convergence_criterion" : "water_pressure_criterion", + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1E-09, + "displacement_relative_tolerance" : 0.0001, + "displacement_absolute_tolerance" : 1E-09, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1E-09, + "min_iterations" : 6, + "max_iterations" : 30, + "number_cycles" : 1, + "reduction_factor" : 1.0, + "increase_factor" : 1.0, + "desired_iterations" : 4, + "max_radius_factor" : 10.0, + "min_radius_factor" : 0.1, + "calculate_reactions" : true, + "max_line_search_iterations" : 5, + "first_alpha_value" : 0.5, + "second_alpha_value" : 1.0, + "min_alpha" : 0.1, + "max_alpha" : 2.0, + "line_search_tolerance" : 0.5, + "rotation_dofs" : false, + "linear_solver_settings": { + "solver_type" : "LinearSolversApplication.sparse_lu", + "scaling" : true + }, + "problem_domain_sub_model_part_list": ["Parts_Solid_Grond" ], + "processes_sub_model_part_list" : ["gravity", "WATER_PRESSURE_Bottom", "WATER_PRESSURE_Top"], + "body_domain_sub_model_part_list" : ["Parts_Solid_Grond" ], + "newmark_beta" : 0.25, + "newmark_gamma" : 0.5, + "newmark_theta" : 0.5, + "rayleigh_m" : 0.0, + "rayleigh_k" : 0.0, + "move_mesh_flag" : false + }, + "output_processes": { + "gid_output": [ + { + "python_module" : "gid_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "GiDOutputProcess", + "Parameters": { + "model_part_name" : "PorousDomain.porous_computational_model_part", + "output_name" : "stage2", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag" : "WriteElementsOnly", + "GiDPostMode" : "GiD_PostAscii", + "MultiFileFlag" : "SingleFile" + }, + "file_label" : "step", + "output_control_type" : "step", + "output_interval" : 1, + "body_output" : true, + "node_output" : false, + "skin_output" : false, + "plane_output" : [], + "nodal_results" : ["WATER_PRESSURE"], + "gauss_point_results" : [] + }, + "point_data_configuration" : [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Bottom", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Top", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "DISPLACEMENT", + "active" : [true, true, true], + "is_fixed" : [true, true, true], + "value" : [0.0, 0.0, 0.0], + "table" : [0, 0, 0] + } + } + ], + "loads_process_list": [ + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "VOLUME_ACCELERATION", + "active" : [false, true, false], + "value" : [0.0, -9.81, 0.0], + "table" : [0, 0, 0] + } + } + ], + "auxiliar_process_list": [] + } +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_difforder_triangle6n/column.mdpa b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_difforder_triangle6n/column.mdpa new file mode 100644 index 000000000000..fc3ff8c95ae3 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_difforder_triangle6n/column.mdpa @@ -0,0 +1,393 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties + +Begin Nodes + 1 0.0000000000 -5.0000000000 0.0000000000 + 2 1.0000000000 -5.0000000000 0.0000000000 + 3 0.0000000000 -4.5000000000 0.0000000000 + 4 1.0000000000 -4.5000000000 0.0000000000 + 5 0.0000000000 -4.0000000000 0.0000000000 + 6 1.0000000000 -4.0000000000 0.0000000000 + 7 0.0000000000 -3.5000000000 0.0000000000 + 8 1.0000000000 -3.5000000000 0.0000000000 + 9 0.0000000000 -3.0000000000 0.0000000000 + 10 1.0000000000 -3.0000000000 0.0000000000 + 11 0.0000000000 -2.5000000000 0.0000000000 + 12 1.0000000000 -2.5000000000 0.0000000000 + 13 0.0000000000 -2.0000000000 0.0000000000 + 14 1.0000000000 -2.0000000000 0.0000000000 + 15 0.0000000000 -1.5000000000 0.0000000000 + 16 1.0000000000 -1.5000000000 0.0000000000 + 17 0.0000000000 -1.0000000000 0.0000000000 + 18 1.0000000000 -1.0000000000 0.0000000000 + 19 0.0000000000 -0.5000000000 0.0000000000 + 20 1.0000000000 -0.5000000000 0.0000000000 + 21 0.0000000000 0.0000000000 0.0000000000 + 22 1.0000000000 0.0000000000 0.0000000000 + 23 0.0000000000 -4.7500000000 0.0000000000 + 24 0.0000000000 -4.2500000000 0.0000000000 + 25 0.0000000000 -3.7500000000 0.0000000000 + 26 0.0000000000 -3.2500000000 0.0000000000 + 27 0.0000000000 -2.7500000000 0.0000000000 + 28 0.0000000000 -2.2500000000 0.0000000000 + 29 0.0000000000 -1.7500000000 0.0000000000 + 30 0.0000000000 -1.2500000000 0.0000000000 + 31 0.0000000000 -0.7500000000 0.0000000000 + 32 0.0000000000 -0.2500000000 0.0000000000 + 33 1.0000000000 -4.7500000000 0.0000000000 + 34 1.0000000000 -4.2500000000 0.0000000000 + 35 1.0000000000 -3.7500000000 0.0000000000 + 36 1.0000000000 -3.2500000000 0.0000000000 + 37 1.0000000000 -2.7500000000 0.0000000000 + 38 1.0000000000 -2.2500000000 0.0000000000 + 39 1.0000000000 -1.7500000000 0.0000000000 + 40 1.0000000000 -1.2500000000 0.0000000000 + 41 1.0000000000 -0.7500000000 0.0000000000 + 42 1.0000000000 -0.2500000000 0.0000000000 + 43 0.5000000000 -5.0000000000 0.0000000000 + 44 0.5000000000 -4.7500000000 0.0000000000 + 45 0.5000000000 -4.5000000000 0.0000000000 + 46 0.5000000000 -4.2500000000 0.0000000000 + 47 0.5000000000 -4.0000000000 0.0000000000 + 48 0.5000000000 -3.7500000000 0.0000000000 + 49 0.5000000000 -3.5000000000 0.0000000000 + 50 0.5000000000 -3.2500000000 0.0000000000 + 51 0.5000000000 -3.0000000000 0.0000000000 + 52 0.5000000000 -2.7500000000 0.0000000000 + 53 0.5000000000 -2.5000000000 0.0000000000 + 54 0.5000000000 -2.2500000000 0.0000000000 + 55 0.5000000000 -2.0000000000 0.0000000000 + 56 0.5000000000 -1.7500000000 0.0000000000 + 57 0.5000000000 -1.5000000000 0.0000000000 + 58 0.5000000000 -1.2500000000 0.0000000000 + 59 0.5000000000 -1.0000000000 0.0000000000 + 60 0.5000000000 -0.7500000000 0.0000000000 + 61 0.5000000000 -0.5000000000 0.0000000000 + 62 0.5000000000 -0.2500000000 0.0000000000 + 63 0.5000000000 0.0000000000 0.0000000000 +End Nodes + + +Begin Elements SteadyStatePwElement2D6N // GUI group identifier: Grond + 1 0 1 2 4 43 33 44 + 2 0 4 3 1 45 23 44 + 3 0 3 4 6 45 34 46 + 4 0 6 5 3 47 24 46 + 5 0 5 6 8 47 35 48 + 6 0 8 7 5 49 25 48 + 7 0 7 8 10 49 36 50 + 8 0 10 9 7 51 26 50 + 9 0 9 10 12 51 37 52 + 10 0 12 11 9 53 27 52 + 11 0 11 12 14 53 38 54 + 12 0 14 13 11 55 28 54 + 13 0 13 14 16 55 39 56 + 14 0 16 15 13 57 29 56 + 15 0 15 16 18 57 40 58 + 16 0 18 17 15 59 30 58 + 17 0 17 18 20 59 41 60 + 18 0 20 19 17 61 31 60 + 19 0 19 20 22 61 42 62 + 20 0 22 21 19 63 32 62 +End Elements + +Begin SubModelPart Parts_Solid_Grond // Group Grond // Subtree Parts_Solid + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + 54 + 55 + 56 + 57 + 58 + 59 + 60 + 61 + 62 + 63 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + End SubModelPartElements +End SubModelPart + +Begin SubModelPart gravity + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + 54 + 55 + 56 + 57 + 58 + 59 + 60 + 61 + 62 + 63 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart DISPLACEMENT_Bottom // Group Bottom // Subtree DISPLACEMENT + Begin SubModelPartNodes + 1 + 2 + 43 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart DISPLACEMENT_Sides // Group Sides // Subtree DISPLACEMENT + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Initial + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + 54 + 55 + 56 + 57 + 58 + 59 + 60 + 61 + 62 + 63 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Top + Begin SubModelPartNodes + 21 + 22 + 63 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Bottom + Begin SubModelPartNodes + 1 + 2 + 43 + End SubModelPartNodes +End SubModelPart diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle3n/ProjectParameters_stage1.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle3n/ProjectParameters_stage1.json new file mode 100644 index 000000000000..27dcb79d2498 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle3n/ProjectParameters_stage1.json @@ -0,0 +1,192 @@ +{ + "problem_data": { + "problem_name" : "test_saturated_below_phreatic_level_upw_smallstrain_triangle3n_stage1", + "start_time" : -0.1, + "end_time" : 0.0, + "echo_level" : 1, + "parallel_type" : "OpenMP", + "number_of_threads" : 1 + }, + "solver_settings": { + "solver_type" : "U_Pw", + "model_part_name" : "PorousDomain", + "domain_size" : 2, + "model_import_settings": { + "input_type" : "mdpa", + "input_filename" : "column" + }, + "material_import_settings": { + "materials_filename" : "../Common/MaterialParameters_stage1.json" + }, + "time_stepping": { + "time_step" : 0.1, + "max_delta_time_factor" : 1.0 + }, + "buffer_size" : 2, + "echo_level" : 1, + "clear_storage" : false, + "compute_reactions" : true, + "reform_dofs_at_each_step" : false, + "nodal_smoothing" : false, + "block_builder" : true, + "solution_type" : "Quasi-Static", + "scheme_type" : "Newmark", + "reset_displacements" : true, + "strategy_type" : "newton_raphson", + "convergence_criterion" : "water_pressure_criterion", + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1E-09, + "displacement_relative_tolerance" : 0.001, + "displacement_absolute_tolerance" : 1E-09, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1E-09, + "min_iterations" : 6, + "max_iterations" : 30, + "number_cycles" : 1, + "reduction_factor" : 1.0, + "increase_factor" : 1.0, + "desired_iterations" : 4, + "max_radius_factor" : 10.0, + "min_radius_factor" : 0.1, + "calculate_reactions" : true, + "max_line_search_iterations" : 5, + "first_alpha_value" : 0.5, + "second_alpha_value" : 1.0, + "min_alpha" : 0.1, + "max_alpha" : 2.0, + "line_search_tolerance" : 0.5, + "rotation_dofs" : false, + "linear_solver_settings": { + "solver_type" : "LinearSolversApplication.sparse_lu", + "scaling" : true + }, + "problem_domain_sub_model_part_list": ["Parts_Solid_Grond"], + "processes_sub_model_part_list" : ["gravity", "WATER_PRESSURE_Initial", "WATER_PRESSURE_Bottom", "WATER_PRESSURE_Top"], + "body_domain_sub_model_part_list" : ["Parts_Solid_Grond"], + "newmark_beta" : 0.25, + "newmark_gamma" : 0.5, + "newmark_theta" : 0.5, + "rayleigh_m" : 0.0, + "rayleigh_k" : 0.0, + "move_mesh_flag" : false + }, + "output_processes": { + "gid_output": [ + { + "python_module" : "gid_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "GiDOutputProcess", + "Parameters": { + "model_part_name" : "PorousDomain.porous_computational_model_part", + "output_name" : "stage1", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag" : "WriteElementsOnly", + "GiDPostMode" : "GiD_PostAscii", + "MultiFileFlag" : "SingleFile" + }, + "file_label" : "step", + "output_control_type" : "step", + "output_interval" : 1, + "body_output" : true, + "node_output" : false, + "skin_output" : false, + "plane_output" : [], + "nodal_results" : ["WATER_PRESSURE"], + "gauss_point_results" : [] + }, + "point_data_configuration" : [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Initial", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : false, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Bottom", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Top", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "DISPLACEMENT", + "active" : [true, true, true], + "is_fixed" : [true, true, true], + "value" : [0.0, 0.0, 0.0], + "table" : [0, 0, 0] + } + } + ], + "loads_process_list": [ + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "VOLUME_ACCELERATION", + "active" : [false, true, false], + "value" : [0.0, -9.81, 0.0], + "table" : [0, 0, 0] + } + } + ], + "auxiliar_process_list": [] + } +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle3n/ProjectParameters_stage2.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle3n/ProjectParameters_stage2.json new file mode 100644 index 000000000000..2097c94444a5 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle3n/ProjectParameters_stage2.json @@ -0,0 +1,174 @@ +{ + "problem_data": { + "problem_name" : "test_saturated_below_phreatic_level_upw_smallstrain_triangle3n_stage2", + "start_time" : 0.0, + "end_time" : 1.0, + "echo_level" : 1, + "parallel_type" : "OpenMP", + "number_of_threads" : 1 + }, + "solver_settings": { + "solver_type" : "U_Pw", + "model_part_name" : "PorousDomain", + "domain_size" : 2, + "model_import_settings": { + "input_type" : "mdpa", + "input_filename" : "column" + }, + "material_import_settings": { + "materials_filename" : "../Common/MaterialParameters_stage2.json" + }, + "time_stepping": { + "time_step" : 0.1, + "max_delta_time_factor" : 1.0 + }, + "buffer_size" : 2, + "echo_level" : 1, + "clear_storage" : false, + "compute_reactions" : true, + "reform_dofs_at_each_step" : false, + "nodal_smoothing" : false, + "block_builder" : true, + "solution_type" : "Quasi-Static", + "scheme_type" : "Newmark", + "reset_displacements" : false, + "strategy_type" : "newton_raphson", + "convergence_criterion" : "water_pressure_criterion", + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1E-09, + "displacement_relative_tolerance" : 0.0001, + "displacement_absolute_tolerance" : 1E-09, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1E-09, + "min_iterations" : 6, + "max_iterations" : 30, + "number_cycles" : 1, + "reduction_factor" : 1.0, + "increase_factor" : 1.0, + "desired_iterations" : 4, + "max_radius_factor" : 10.0, + "min_radius_factor" : 0.1, + "calculate_reactions" : true, + "max_line_search_iterations" : 5, + "first_alpha_value" : 0.5, + "second_alpha_value" : 1.0, + "min_alpha" : 0.1, + "max_alpha" : 2.0, + "line_search_tolerance" : 0.5, + "rotation_dofs" : false, + "linear_solver_settings": { + "solver_type" : "LinearSolversApplication.sparse_lu", + "scaling" : true + }, + "problem_domain_sub_model_part_list": ["Parts_Solid_Grond" ], + "processes_sub_model_part_list" : ["gravity", "WATER_PRESSURE_Bottom", "WATER_PRESSURE_Top"], + "body_domain_sub_model_part_list" : ["Parts_Solid_Grond" ], + "newmark_beta" : 0.25, + "newmark_gamma" : 0.5, + "newmark_theta" : 0.5, + "rayleigh_m" : 0.0, + "rayleigh_k" : 0.0, + "move_mesh_flag" : false + }, + "output_processes": { + "gid_output": [ + { + "python_module" : "gid_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "GiDOutputProcess", + "Parameters": { + "model_part_name" : "PorousDomain.porous_computational_model_part", + "output_name" : "stage2", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag" : "WriteElementsOnly", + "GiDPostMode" : "GiD_PostAscii", + "MultiFileFlag" : "SingleFile" + }, + "file_label" : "step", + "output_control_type" : "step", + "output_interval" : 1, + "body_output" : true, + "node_output" : false, + "skin_output" : false, + "plane_output" : [], + "nodal_results" : ["WATER_PRESSURE"], + "gauss_point_results" : [] + }, + "point_data_configuration" : [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Bottom", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Top", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "DISPLACEMENT", + "active" : [true, true, true], + "is_fixed" : [true, true, true], + "value" : [0.0, 0.0, 0.0], + "table" : [0, 0, 0] + } + } + ], + "loads_process_list": [ + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "VOLUME_ACCELERATION", + "active" : [false, true, false], + "value" : [0.0, -9.81, 0.0], + "table" : [0, 0, 0] + } + } + ], + "auxiliar_process_list": [] + } +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle3n/column.mdpa b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle3n/column.mdpa new file mode 100644 index 000000000000..f6e16007f1df --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle3n/column.mdpa @@ -0,0 +1,206 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties + +Begin Nodes + 1 0.0000000000 -5.0000000000 0.0000000000 + 2 1.0000000000 -5.0000000000 0.0000000000 + 3 0.0000000000 -4.5000000000 0.0000000000 + 4 1.0000000000 -4.5000000000 0.0000000000 + 5 0.0000000000 -4.0000000000 0.0000000000 + 6 1.0000000000 -4.0000000000 0.0000000000 + 7 0.0000000000 -3.5000000000 0.0000000000 + 8 1.0000000000 -3.5000000000 0.0000000000 + 9 0.0000000000 -3.0000000000 0.0000000000 + 10 1.0000000000 -3.0000000000 0.0000000000 + 11 0.0000000000 -2.5000000000 0.0000000000 + 12 1.0000000000 -2.5000000000 0.0000000000 + 13 0.0000000000 -2.0000000000 0.0000000000 + 14 1.0000000000 -2.0000000000 0.0000000000 + 15 0.0000000000 -1.5000000000 0.0000000000 + 16 1.0000000000 -1.5000000000 0.0000000000 + 17 0.0000000000 -1.0000000000 0.0000000000 + 18 1.0000000000 -1.0000000000 0.0000000000 + 19 0.0000000000 -0.5000000000 0.0000000000 + 20 1.0000000000 -0.5000000000 0.0000000000 + 21 0.0000000000 0.0000000000 0.0000000000 + 22 1.0000000000 0.0000000000 0.0000000000 +End Nodes + + +Begin Elements UPwSmallStrainElement2D3N // GUI group identifier: Grond + 1 0 1 2 4 + 2 0 4 3 1 + 3 0 3 4 6 + 4 0 6 5 3 + 5 0 5 6 8 + 6 0 8 7 5 + 7 0 7 8 10 + 8 0 10 9 7 + 9 0 9 10 12 + 10 0 12 11 9 + 11 0 11 12 14 + 12 0 14 13 11 + 13 0 13 14 16 + 14 0 16 15 13 + 15 0 15 16 18 + 16 0 18 17 15 + 17 0 17 18 20 + 18 0 20 19 17 + 19 0 19 20 22 + 20 0 22 21 19 +End Elements + +Begin SubModelPart Parts_Solid_Grond // Group Grond // Subtree Parts_Solid + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + End SubModelPartElements +End SubModelPart + +Begin SubModelPart gravity + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart DISPLACEMENT_Bottom // Group Bottom // Subtree DISPLACEMENT + Begin SubModelPartNodes + 1 + 2 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart DISPLACEMENT_Sides // Group Sides // Subtree DISPLACEMENT + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Initial + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Top + Begin SubModelPartNodes + 21 + 22 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Bottom + Begin SubModelPartNodes + 1 + 2 + End SubModelPartNodes +End SubModelPart diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle6n/ProjectParameters_stage1.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle6n/ProjectParameters_stage1.json new file mode 100644 index 000000000000..27dcb79d2498 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle6n/ProjectParameters_stage1.json @@ -0,0 +1,192 @@ +{ + "problem_data": { + "problem_name" : "test_saturated_below_phreatic_level_upw_smallstrain_triangle3n_stage1", + "start_time" : -0.1, + "end_time" : 0.0, + "echo_level" : 1, + "parallel_type" : "OpenMP", + "number_of_threads" : 1 + }, + "solver_settings": { + "solver_type" : "U_Pw", + "model_part_name" : "PorousDomain", + "domain_size" : 2, + "model_import_settings": { + "input_type" : "mdpa", + "input_filename" : "column" + }, + "material_import_settings": { + "materials_filename" : "../Common/MaterialParameters_stage1.json" + }, + "time_stepping": { + "time_step" : 0.1, + "max_delta_time_factor" : 1.0 + }, + "buffer_size" : 2, + "echo_level" : 1, + "clear_storage" : false, + "compute_reactions" : true, + "reform_dofs_at_each_step" : false, + "nodal_smoothing" : false, + "block_builder" : true, + "solution_type" : "Quasi-Static", + "scheme_type" : "Newmark", + "reset_displacements" : true, + "strategy_type" : "newton_raphson", + "convergence_criterion" : "water_pressure_criterion", + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1E-09, + "displacement_relative_tolerance" : 0.001, + "displacement_absolute_tolerance" : 1E-09, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1E-09, + "min_iterations" : 6, + "max_iterations" : 30, + "number_cycles" : 1, + "reduction_factor" : 1.0, + "increase_factor" : 1.0, + "desired_iterations" : 4, + "max_radius_factor" : 10.0, + "min_radius_factor" : 0.1, + "calculate_reactions" : true, + "max_line_search_iterations" : 5, + "first_alpha_value" : 0.5, + "second_alpha_value" : 1.0, + "min_alpha" : 0.1, + "max_alpha" : 2.0, + "line_search_tolerance" : 0.5, + "rotation_dofs" : false, + "linear_solver_settings": { + "solver_type" : "LinearSolversApplication.sparse_lu", + "scaling" : true + }, + "problem_domain_sub_model_part_list": ["Parts_Solid_Grond"], + "processes_sub_model_part_list" : ["gravity", "WATER_PRESSURE_Initial", "WATER_PRESSURE_Bottom", "WATER_PRESSURE_Top"], + "body_domain_sub_model_part_list" : ["Parts_Solid_Grond"], + "newmark_beta" : 0.25, + "newmark_gamma" : 0.5, + "newmark_theta" : 0.5, + "rayleigh_m" : 0.0, + "rayleigh_k" : 0.0, + "move_mesh_flag" : false + }, + "output_processes": { + "gid_output": [ + { + "python_module" : "gid_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "GiDOutputProcess", + "Parameters": { + "model_part_name" : "PorousDomain.porous_computational_model_part", + "output_name" : "stage1", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag" : "WriteElementsOnly", + "GiDPostMode" : "GiD_PostAscii", + "MultiFileFlag" : "SingleFile" + }, + "file_label" : "step", + "output_control_type" : "step", + "output_interval" : 1, + "body_output" : true, + "node_output" : false, + "skin_output" : false, + "plane_output" : [], + "nodal_results" : ["WATER_PRESSURE"], + "gauss_point_results" : [] + }, + "point_data_configuration" : [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Initial", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : false, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Bottom", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Top", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "DISPLACEMENT", + "active" : [true, true, true], + "is_fixed" : [true, true, true], + "value" : [0.0, 0.0, 0.0], + "table" : [0, 0, 0] + } + } + ], + "loads_process_list": [ + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "VOLUME_ACCELERATION", + "active" : [false, true, false], + "value" : [0.0, -9.81, 0.0], + "table" : [0, 0, 0] + } + } + ], + "auxiliar_process_list": [] + } +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle6n/ProjectParameters_stage2.json b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle6n/ProjectParameters_stage2.json new file mode 100644 index 000000000000..5c101f84528a --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle6n/ProjectParameters_stage2.json @@ -0,0 +1,174 @@ +{ + "problem_data": { + "problem_name" : "test_saturated_below_phreatic_level_upw_smallstrain_triangle6n_stage2", + "start_time" : 0.0, + "end_time" : 1.0, + "echo_level" : 1, + "parallel_type" : "OpenMP", + "number_of_threads" : 1 + }, + "solver_settings": { + "solver_type" : "U_Pw", + "model_part_name" : "PorousDomain", + "domain_size" : 2, + "model_import_settings": { + "input_type" : "mdpa", + "input_filename" : "column" + }, + "material_import_settings": { + "materials_filename" : "../Common/MaterialParameters_stage2.json" + }, + "time_stepping": { + "time_step" : 0.1, + "max_delta_time_factor" : 1.0 + }, + "buffer_size" : 2, + "echo_level" : 1, + "clear_storage" : false, + "compute_reactions" : true, + "reform_dofs_at_each_step" : false, + "nodal_smoothing" : false, + "block_builder" : true, + "solution_type" : "Quasi-Static", + "scheme_type" : "Newmark", + "reset_displacements" : false, + "strategy_type" : "newton_raphson", + "convergence_criterion" : "water_pressure_criterion", + "water_pressure_relative_tolerance": 0.0001, + "water_pressure_absolute_tolerance": 1E-09, + "displacement_relative_tolerance" : 0.0001, + "displacement_absolute_tolerance" : 1E-09, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1E-09, + "min_iterations" : 6, + "max_iterations" : 30, + "number_cycles" : 1, + "reduction_factor" : 1.0, + "increase_factor" : 1.0, + "desired_iterations" : 4, + "max_radius_factor" : 10.0, + "min_radius_factor" : 0.1, + "calculate_reactions" : true, + "max_line_search_iterations" : 5, + "first_alpha_value" : 0.5, + "second_alpha_value" : 1.0, + "min_alpha" : 0.1, + "max_alpha" : 2.0, + "line_search_tolerance" : 0.5, + "rotation_dofs" : false, + "linear_solver_settings": { + "solver_type" : "LinearSolversApplication.sparse_lu", + "scaling" : true + }, + "problem_domain_sub_model_part_list": ["Parts_Solid_Grond" ], + "processes_sub_model_part_list" : ["gravity", "WATER_PRESSURE_Bottom", "WATER_PRESSURE_Top"], + "body_domain_sub_model_part_list" : ["Parts_Solid_Grond" ], + "newmark_beta" : 0.25, + "newmark_gamma" : 0.5, + "newmark_theta" : 0.5, + "rayleigh_m" : 0.0, + "rayleigh_k" : 0.0, + "move_mesh_flag" : false + }, + "output_processes": { + "gid_output": [ + { + "python_module" : "gid_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "GiDOutputProcess", + "Parameters": { + "model_part_name" : "PorousDomain.porous_computational_model_part", + "output_name" : "stage2", + "postprocess_parameters": { + "result_file_configuration": { + "gidpost_flags": { + "WriteDeformedMeshFlag": "WriteUndeformed", + "WriteConditionsFlag" : "WriteElementsOnly", + "GiDPostMode" : "GiD_PostAscii", + "MultiFileFlag" : "SingleFile" + }, + "file_label" : "step", + "output_control_type" : "step", + "output_interval" : 1, + "body_output" : true, + "node_output" : false, + "skin_output" : false, + "plane_output" : [], + "nodal_results" : ["WATER_PRESSURE"], + "gauss_point_results" : [] + }, + "point_data_configuration" : [] + } + } + } + ] + }, + "processes": { + "constraints_process_list": [ + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Bottom", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_scalar_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.WATER_PRESSURE_Top", + "variable_name" : "WATER_PRESSURE", + "is_fixed" : true, + "table" : 0, + "fluid_pressure_type" : "Phreatic_Multi_Line", + "gravity_direction" : 1, + "out_of_plane_direction" : 2, + "x_coordinates" : [ 0.0, 1.0], + "y_coordinates" : [-2.0, -2.0], + "z_coordinates" : [ 0.0, 0.0], + "specific_weight" : 10000.0 + } + }, + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "DISPLACEMENT", + "active" : [true, true, true], + "is_fixed" : [true, true, true], + "value" : [0.0, 0.0, 0.0], + "table" : [0, 0, 0] + } + } + ], + "loads_process_list": [ + { + "python_module" : "apply_vector_constraint_table_process", + "kratos_module" : "KratosMultiphysics.GeoMechanicsApplication", + "process_name" : "ApplyVectorConstraintTableProcess", + "Parameters": { + "model_part_name" : "PorousDomain.gravity", + "variable_name" : "VOLUME_ACCELERATION", + "active" : [false, true, false], + "value" : [0.0, -9.81, 0.0], + "table" : [0, 0, 0] + } + } + ], + "auxiliar_process_list": [] + } +} diff --git a/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle6n/column.mdpa b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle6n/column.mdpa new file mode 100644 index 000000000000..0a051af8658a --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/test_partially_saturated/test_saturated_below_phreatic_level_upw_smallstrain_triangle6n/column.mdpa @@ -0,0 +1,393 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties + +Begin Nodes + 1 0.0000000000 -5.0000000000 0.0000000000 + 2 1.0000000000 -5.0000000000 0.0000000000 + 3 0.0000000000 -4.5000000000 0.0000000000 + 4 1.0000000000 -4.5000000000 0.0000000000 + 5 0.0000000000 -4.0000000000 0.0000000000 + 6 1.0000000000 -4.0000000000 0.0000000000 + 7 0.0000000000 -3.5000000000 0.0000000000 + 8 1.0000000000 -3.5000000000 0.0000000000 + 9 0.0000000000 -3.0000000000 0.0000000000 + 10 1.0000000000 -3.0000000000 0.0000000000 + 11 0.0000000000 -2.5000000000 0.0000000000 + 12 1.0000000000 -2.5000000000 0.0000000000 + 13 0.0000000000 -2.0000000000 0.0000000000 + 14 1.0000000000 -2.0000000000 0.0000000000 + 15 0.0000000000 -1.5000000000 0.0000000000 + 16 1.0000000000 -1.5000000000 0.0000000000 + 17 0.0000000000 -1.0000000000 0.0000000000 + 18 1.0000000000 -1.0000000000 0.0000000000 + 19 0.0000000000 -0.5000000000 0.0000000000 + 20 1.0000000000 -0.5000000000 0.0000000000 + 21 0.0000000000 0.0000000000 0.0000000000 + 22 1.0000000000 0.0000000000 0.0000000000 + 23 0.0000000000 -4.7500000000 0.0000000000 + 24 0.0000000000 -4.2500000000 0.0000000000 + 25 0.0000000000 -3.7500000000 0.0000000000 + 26 0.0000000000 -3.2500000000 0.0000000000 + 27 0.0000000000 -2.7500000000 0.0000000000 + 28 0.0000000000 -2.2500000000 0.0000000000 + 29 0.0000000000 -1.7500000000 0.0000000000 + 30 0.0000000000 -1.2500000000 0.0000000000 + 31 0.0000000000 -0.7500000000 0.0000000000 + 32 0.0000000000 -0.2500000000 0.0000000000 + 33 1.0000000000 -4.7500000000 0.0000000000 + 34 1.0000000000 -4.2500000000 0.0000000000 + 35 1.0000000000 -3.7500000000 0.0000000000 + 36 1.0000000000 -3.2500000000 0.0000000000 + 37 1.0000000000 -2.7500000000 0.0000000000 + 38 1.0000000000 -2.2500000000 0.0000000000 + 39 1.0000000000 -1.7500000000 0.0000000000 + 40 1.0000000000 -1.2500000000 0.0000000000 + 41 1.0000000000 -0.7500000000 0.0000000000 + 42 1.0000000000 -0.2500000000 0.0000000000 + 43 0.5000000000 -5.0000000000 0.0000000000 + 44 0.5000000000 -4.7500000000 0.0000000000 + 45 0.5000000000 -4.5000000000 0.0000000000 + 46 0.5000000000 -4.2500000000 0.0000000000 + 47 0.5000000000 -4.0000000000 0.0000000000 + 48 0.5000000000 -3.7500000000 0.0000000000 + 49 0.5000000000 -3.5000000000 0.0000000000 + 50 0.5000000000 -3.2500000000 0.0000000000 + 51 0.5000000000 -3.0000000000 0.0000000000 + 52 0.5000000000 -2.7500000000 0.0000000000 + 53 0.5000000000 -2.5000000000 0.0000000000 + 54 0.5000000000 -2.2500000000 0.0000000000 + 55 0.5000000000 -2.0000000000 0.0000000000 + 56 0.5000000000 -1.7500000000 0.0000000000 + 57 0.5000000000 -1.5000000000 0.0000000000 + 58 0.5000000000 -1.2500000000 0.0000000000 + 59 0.5000000000 -1.0000000000 0.0000000000 + 60 0.5000000000 -0.7500000000 0.0000000000 + 61 0.5000000000 -0.5000000000 0.0000000000 + 62 0.5000000000 -0.2500000000 0.0000000000 + 63 0.5000000000 0.0000000000 0.0000000000 +End Nodes + + +Begin Elements UPwSmallStrainElement2D6N // GUI group identifier: Grond + 1 0 1 2 4 43 33 44 + 2 0 4 3 1 45 23 44 + 3 0 3 4 6 45 34 46 + 4 0 6 5 3 47 24 46 + 5 0 5 6 8 47 35 48 + 6 0 8 7 5 49 25 48 + 7 0 7 8 10 49 36 50 + 8 0 10 9 7 51 26 50 + 9 0 9 10 12 51 37 52 + 10 0 12 11 9 53 27 52 + 11 0 11 12 14 53 38 54 + 12 0 14 13 11 55 28 54 + 13 0 13 14 16 55 39 56 + 14 0 16 15 13 57 29 56 + 15 0 15 16 18 57 40 58 + 16 0 18 17 15 59 30 58 + 17 0 17 18 20 59 41 60 + 18 0 20 19 17 61 31 60 + 19 0 19 20 22 61 42 62 + 20 0 22 21 19 63 32 62 +End Elements + +Begin SubModelPart Parts_Solid_Grond // Group Grond // Subtree Parts_Solid + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + 54 + 55 + 56 + 57 + 58 + 59 + 60 + 61 + 62 + 63 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + End SubModelPartElements +End SubModelPart + +Begin SubModelPart gravity + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + 54 + 55 + 56 + 57 + 58 + 59 + 60 + 61 + 62 + 63 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart DISPLACEMENT_Bottom // Group Bottom // Subtree DISPLACEMENT + Begin SubModelPartNodes + 1 + 2 + 43 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart DISPLACEMENT_Sides // Group Sides // Subtree DISPLACEMENT + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Initial + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + 15 + 16 + 17 + 18 + 19 + 20 + 21 + 22 + 23 + 24 + 25 + 26 + 27 + 28 + 29 + 30 + 31 + 32 + 33 + 34 + 35 + 36 + 37 + 38 + 39 + 40 + 41 + 42 + 43 + 44 + 45 + 46 + 47 + 48 + 49 + 50 + 51 + 52 + 53 + 54 + 55 + 56 + 57 + 58 + 59 + 60 + 61 + 62 + 63 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Top + Begin SubModelPartNodes + 21 + 22 + 63 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart WATER_PRESSURE_Bottom + Begin SubModelPartNodes + 1 + 2 + 43 + End SubModelPartNodes +End SubModelPart diff --git a/applications/GeoMechanicsApplication/tests/test_transient_groundwater_flow.py b/applications/GeoMechanicsApplication/tests/test_transient_groundwater_flow.py index 8ed90589da34..897b69ea1876 100644 --- a/applications/GeoMechanicsApplication/tests/test_transient_groundwater_flow.py +++ b/applications/GeoMechanicsApplication/tests/test_transient_groundwater_flow.py @@ -26,7 +26,7 @@ def test_Transient_Case_B1_2D3N(self): p_height_1_72_m = water_pressure_stage_2[32] - self.assertAlmostEqual(1.9762325593327588, p_height_1_72_m, 3) + self.assertAlmostEqual(1.8902712667818171, p_height_1_72_m, 3) def test_Transient_Case_A1_2D3N(self): test_name = 'test_Transient_Case_A1_2D3N' @@ -37,7 +37,7 @@ def test_Transient_Case_A1_2D3N(self): p_height_1_72_m = water_pressure_stage_2[32] - self.assertAlmostEqual(3.7224838666844984, p_height_1_72_m) + self.assertAlmostEqual(3.730389835521539, p_height_1_72_m) def test_Transient_Case_A1_2D6N(self): test_name = 'test_Transient_Case_A1_2D6N' @@ -48,7 +48,7 @@ def test_Transient_Case_A1_2D6N(self): p_height_1_72_m = water_pressure[87] - self.assertAlmostEqual(2.6423560731957005, p_height_1_72_m) + self.assertAlmostEqual(2.617220832597392, p_height_1_72_m) if __name__ == '__main__':