diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp index 85fded5aef17..d1006505f1b8 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp @@ -139,8 +139,8 @@ void UPwSmallStrainFICElement::InitializeNonLinearIteration(con const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = this->CalculateDeformationGradients(); - auto strain_vectors = this->CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, @@ -194,8 +194,8 @@ void UPwSmallStrainFICElement::FinalizeNonLinearIteration(const const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = this->CalculateDeformationGradients(); - auto strain_vectors = this->CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, @@ -449,8 +449,8 @@ void UPwSmallStrainFICElement::CalculateAll(MatrixType& rLeftHa const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); const auto deformation_gradients = this->CalculateDeformationGradients(); - auto strain_vectors = this->CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, 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 dbebb87a5bd0..823c63ae7751 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -155,8 +155,8 @@ void UPwSmallStrainElement::InitializeSolutionStep(const Proces const auto deformation_gradients = CalculateDeformationGradients(); const auto determinants_of_deformation_gradients = GeoMechanicsMathUtilities::CalculateDeterminants(deformation_gradients); - const auto strain_vectors = CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + const auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { @@ -262,8 +262,8 @@ void UPwSmallStrainElement::InitializeNonLinearIteration(const const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); - auto strain_vectors = CalculateStrains(deformation_gradients, b_matrices, - Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, @@ -307,8 +307,8 @@ void UPwSmallStrainElement::FinalizeSolutionStep(const ProcessI const auto deformation_gradients = CalculateDeformationGradients(); const auto determinants_of_deformation_gradients = GeoMechanicsMathUtilities::CalculateDeterminants(deformation_gradients); - const auto strain_vectors = CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + const auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); Matrix StressContainer(NumGPoints, mStressVector[0].size()); // Loop over integration points @@ -573,8 +573,8 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); - auto strain_vectors = CalculateStrains(deformation_gradients, b_matrices, - Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); const PropertiesType& rProp = this->GetProperties(); @@ -623,8 +623,8 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints( const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); - const auto strain_vectors = CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + const auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); std::vector permeability_update_factors; for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { @@ -694,8 +694,8 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); - auto strain_vectors = CalculateStrains(deformation_gradients, b_matrices, - Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, @@ -738,7 +738,8 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const Variables.B = this->CalculateBMatrix(Variables.GradNpTInitialConfiguration, Variables.Np); // Compute infinitesimal strain - Variables.StrainVector = this->CalculateCauchyStrain(Variables.B, Variables.DisplacementVector); + Variables.StrainVector = + StressStrainUtilities::CalculateCauchyStrain(Variables.B, Variables.DisplacementVector); if (rOutput[GPoint].size() != Variables.StrainVector.size()) rOutput[GPoint].resize(Variables.StrainVector.size(), false); @@ -751,8 +752,9 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); - rOutput = CalculateStrains(deformation_gradients, b_matrices, Variables.DisplacementVector, - Variables.UseHenckyStrain); + rOutput = StressStrainUtilities::CalculateStrains(deformation_gradients, b_matrices, + Variables.DisplacementVector, + Variables.UseHenckyStrain, VoigtSize); } else if (rProp.Has(rVariable)) { // Map initial material property to Gauss points, as required for the output rOutput.clear(); @@ -862,8 +864,8 @@ void UPwSmallStrainElement::CalculateMaterialStiffnessMatrix(Ma this->InitializeElementVariables(Variables, rCurrentProcessInfo); const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); - auto strain_vectors = CalculateStrains(deformation_gradients, b_matrices, - Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, @@ -981,9 +983,9 @@ void UPwSmallStrainElement::CalculateAll(MatrixType& rLe const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); - const auto deformation_gradients = CalculateDeformationGradients(); - auto strain_vectors = CalculateStrains(deformation_gradients, b_matrices, - Variables.DisplacementVector, Variables.UseHenckyStrain); + const auto deformation_gradients = CalculateDeformationGradients(); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, VoigtSize); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, @@ -1498,38 +1500,6 @@ void UPwSmallStrainElement::CalculateAndAddFluidBodyFlow(Vector KRATOS_CATCH("") } -template -Vector UPwSmallStrainElement::CalculateStrain(const Matrix& rDeformationGradient, - const Matrix& rB, - const Vector& rDisplacements, - bool UseHenckyStrain) const -{ - return UseHenckyStrain ? StressStrainUtilities::CalculateHenckyStrain(rDeformationGradient, VoigtSize) - : this->CalculateCauchyStrain(rB, rDisplacements); -} - -template -std::vector UPwSmallStrainElement::CalculateStrains(const std::vector& rDeformationGradients, - const std::vector& rBs, - const Vector& rDisplacements, - bool UseHenckyStrain) const -{ - std::vector result; - std::transform( - rDeformationGradients.begin(), rDeformationGradients.end(), rBs.begin(), std::back_inserter(result), - [this, &rDisplacements, UseHenckyStrain](const auto& rDeformationGradient, const auto& rB) { - return CalculateStrain(rDeformationGradient, rB, rDisplacements, UseHenckyStrain); - }); - - return result; -} - -template -Vector UPwSmallStrainElement::CalculateCauchyStrain(const Matrix& rB, const Vector& rDisplacements) const -{ - return prod(rB, rDisplacements); -} - template Vector UPwSmallStrainElement::CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const { diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp index fcdbd8a33af1..860b0d0cf30d 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -273,16 +273,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa double CalculateBulkModulus(const Matrix& ConstitutiveMatrix) const; double CalculateBiotCoefficient(const ElementVariables& rVariables, bool hasBiotCoefficient) const; - virtual Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const; - virtual Vector CalculateCauchyStrain(const Matrix& rB, const Vector& rDisplacements) const; - virtual Vector CalculateStrain(const Matrix& rDeformationGradient, - const Matrix& rB, - const Vector& rDisplacements, - bool UseHenckyStrain) const; - std::vector CalculateStrains(const std::vector& rDeformationGradients, - const std::vector& rBs, - const Vector& rDisplacements, - bool UseHenckyStrain) const; + virtual Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const; Matrix CalculateDeformationGradient(unsigned int GPoint) const; std::vector CalculateDeformationGradients() const; diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp index 5378da3187df..24335ab04e1c 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp @@ -77,8 +77,9 @@ void UPwUpdatedLagrangianFICElement::CalculateAll(MatrixType& r const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); const auto deformation_gradients = this->CalculateDeformationGradients(); - auto strain_vectors = this->CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, + this->VoigtSize); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp index a97166f67ad7..bc02c59201ec 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp @@ -73,8 +73,9 @@ void UPwUpdatedLagrangianElement::CalculateAll(MatrixType& rLef const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); const auto deformation_gradients = this->CalculateDeformationGradients(); - auto strain_vectors = this->CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, + this->VoigtSize); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, 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 a59e95a17721..a225b95cf37a 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 @@ -320,8 +320,9 @@ void SmallStrainUPwDiffOrderElement::InitializeSolutionStep(const ProcessInfo& r const auto deformation_gradients = CalculateDeformationGradients(); const auto determinants_of_deformation_gradients = GeoMechanicsMathUtilities::CalculateDeterminants(deformation_gradients); - const auto strain_vectors = CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + const auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, + GetVoigtSize()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { @@ -535,8 +536,9 @@ void SmallStrainUPwDiffOrderElement::FinalizeSolutionStep(const ProcessInfo& rCu const auto deformation_gradients = CalculateDeformationGradients(); const auto determinants_of_deformation_gradients = GeoMechanicsMathUtilities::CalculateDeterminants(deformation_gradients); - const auto strain_vectors = CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + const auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, + GetVoigtSize()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { @@ -1003,8 +1005,9 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable RetentionLaw::Parameters RetentionParameters(GetProperties(), rCurrentProcessInfo); const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); - const auto strain_vectors = CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + const auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, + Variables.UseHenckyStrain, GetVoigtSize()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { @@ -1096,7 +1099,8 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable Variables.B = this->CalculateBMatrix(Variables.DNu_DXInitialConfiguration, Variables.Nu); // Compute infinitesimal strain - Variables.StrainVector = this->CalculateCauchyStrain(Variables.B, Variables.DisplacementVector); + Variables.StrainVector = + StressStrainUtilities::CalculateCauchyStrain(Variables.B, Variables.DisplacementVector); if (rOutput[GPoint].size() != Variables.StrainVector.size()) rOutput[GPoint].resize(Variables.StrainVector.size(), false); @@ -1109,8 +1113,9 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable this->InitializeElementVariables(Variables, rCurrentProcessInfo); const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); - rOutput = CalculateStrains(deformation_gradients, b_matrices, Variables.DisplacementVector, - Variables.UseHenckyStrain); + rOutput = StressStrainUtilities::CalculateStrains(deformation_gradients, b_matrices, + Variables.DisplacementVector, + Variables.UseHenckyStrain, GetVoigtSize()); } else if (rVariable == TOTAL_STRESS_VECTOR) { // Definition of variables ElementVariables Variables; @@ -1137,9 +1142,10 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable Vector TotalStressVector(mStressVector[0].size()); const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); - const auto deformation_gradients = CalculateDeformationGradients(); - auto strain_vectors = CalculateStrains(deformation_gradients, b_matrices, - Variables.DisplacementVector, Variables.UseHenckyStrain); + const auto deformation_gradients = CalculateDeformationGradients(); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, + Variables.UseHenckyStrain, GetVoigtSize()); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NuContainer, Variables.DNu_DXContainer, @@ -1277,8 +1283,9 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi const auto deformation_gradients = CalculateDeformationGradients(); const auto integration_coefficients = CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJuContainer); - auto strain_vectors = CalculateStrains(deformation_gradients, b_matrices, - Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, + GetVoigtSize()); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NuContainer, Variables.DNu_DXContainer, @@ -1337,8 +1344,9 @@ void SmallStrainUPwDiffOrderElement::CalculateMaterialStiffnessMatrix(MatrixType const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); - auto strain_vectors = CalculateStrains(deformation_gradients, b_matrices, - Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, + GetVoigtSize()); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NuContainer, Variables.DNu_DXContainer, @@ -2000,38 +2008,9 @@ GeometryData::IntegrationMethod SmallStrainUPwDiffOrderElement::GetIntegrationMe return GI_GAUSS; } -std::vector SmallStrainUPwDiffOrderElement::CalculateStrains(const std::vector& rDeformationGradients, - const std::vector& rBs, - const Vector& rDisplacements, - bool UseHenckyStrain) const -{ - std::vector result; - std::transform( - rDeformationGradients.begin(), rDeformationGradients.end(), rBs.begin(), std::back_inserter(result), - [this, &rDisplacements, UseHenckyStrain](const auto& rDeformationGradient, const auto& rB) { - return CalculateStrain(rDeformationGradient, rB, rDisplacements, UseHenckyStrain); - }); - - return result; -} - -Vector SmallStrainUPwDiffOrderElement::CalculateStrain(const Matrix& rDeformationGradient, - const Matrix& rB, - const Vector& rDisplacements, - bool UseHenckyStrain) const -{ - if (UseHenckyStrain) { - const SizeType Dim = GetGeometry().WorkingSpaceDimension(); - const SizeType VoigtSize = (Dim == N_DIM_3D ? VOIGT_SIZE_3D : VOIGT_SIZE_2D_PLANE_STRAIN); - return StressStrainUtilities::CalculateHenckyStrain(rDeformationGradient, VoigtSize); - } - - return this->CalculateCauchyStrain(rB, rDisplacements); -} - -Vector SmallStrainUPwDiffOrderElement::CalculateCauchyStrain(const Matrix& rB, const Vector& rDisplacements) const +SizeType SmallStrainUPwDiffOrderElement::GetVoigtSize() const { - return prod(rB, rDisplacements); + return (GetGeometry().WorkingSpaceDimension() == N_DIM_3D ? VOIGT_SIZE_3D : VOIGT_SIZE_2D_PLANE_STRAIN); } Vector SmallStrainUPwDiffOrderElement::CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const 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 61d0e28b9634..7f2642c141f5 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 @@ -293,16 +293,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub void AssignPressureToIntermediateNodes(); - virtual Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const; - virtual Vector CalculateCauchyStrain(const Matrix& rB, const Vector& rDisplacements) const; - virtual Vector CalculateStrain(const Matrix& rDeformationGradient, - const Matrix& rB, - const Vector& rDisplacements, - bool UseHenckyStrain) const; - std::vector CalculateStrains(const std::vector& rDeformationGradients, - const std::vector& rBs, - const Vector& rDisplacements, - bool UseHenckyStrain) const; + virtual Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const; + + SizeType GetVoigtSize() const; Matrix CalculateDeformationGradient(unsigned int GPoint) const; std::vector CalculateDeformationGradients() const; diff --git a/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp b/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp index 30449a8679ab..d2de2fa02705 100644 --- a/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp @@ -75,8 +75,9 @@ void UpdatedLagrangianUPwDiffOrderElement::CalculateAll(MatrixType& rLeft const auto deformation_gradients = this->CalculateDeformationGradients(); const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJuContainer); - auto strain_vectors = this->CalculateStrains( - deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain); + auto strain_vectors = StressStrainUtilities::CalculateStrains( + deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, + this->GetVoigtSize()); std::vector constitutive_matrices; this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NuContainer, Variables.DNu_DXContainer, diff --git a/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.hpp b/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.hpp index 6c992886e133..1fbdec4e411f 100644 --- a/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.hpp @@ -67,7 +67,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UpdatedLagrangianUPwDiffOrderElement /// The definition of the sizetype using SizeType = std::size_t; - using SmallStrainUPwDiffOrderElement::CalculateCauchyStrain; using SmallStrainUPwDiffOrderElement::CalculateDerivativesOnInitialConfiguration; using SmallStrainUPwDiffOrderElement::CalculateGreenLagrangeStrain; using SmallStrainUPwDiffOrderElement::mConstitutiveLawVector; diff --git a/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.cpp b/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.cpp index ffb6dc5c55e3..4f3d7ce3dc24 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.cpp +++ b/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.cpp @@ -152,4 +152,26 @@ Matrix StressStrainUtilities::CalculateGreenLagrangeStrainTensor(const Matrix& r IdentityMatrix(rDeformationGradient.size1())); } +Vector StressStrainUtilities::CalculateCauchyStrain(const Matrix& rB, const Vector& rDisplacements) +{ + return prod(rB, rDisplacements); +} + +std::vector StressStrainUtilities::CalculateStrains(const std::vector& rDeformationGradients, + const std::vector& rBs, + const Vector& rDisplacements, + bool UseHenckyStrain, + std::size_t VoigtSize) +{ + std::vector result; + std::transform( + rDeformationGradients.begin(), rDeformationGradients.end(), rBs.begin(), std::back_inserter(result), + [&rDisplacements, UseHenckyStrain, VoigtSize](const auto& rDeformationGradient, const auto& rB) { + return UseHenckyStrain ? CalculateHenckyStrain(rDeformationGradient, VoigtSize) + : CalculateCauchyStrain(rB, rDisplacements); + }); + + return result; +} + } // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.h b/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.h index 1e3e1c5d2aa6..37240ea31f73 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.h +++ b/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.h @@ -32,6 +32,12 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) StressStrainUtilities static double CalculateVonMisesStrain(const Vector& rStrainVector); static Vector CalculateHenckyStrain(const Matrix& rDeformationGradient, size_t VoigtSize); static Matrix CalculateGreenLagrangeStrainTensor(const Matrix& rDeformationGradient); + static Vector CalculateCauchyStrain(const Matrix& rB, const Vector& rDisplacements); + static std::vector CalculateStrains(const std::vector& rDeformationGradients, + const std::vector& rBs, + const Vector& rDisplacements, + bool UseHenckyStrain, + std::size_t VoigtSize); private: static double CalculateQMohrCoulomb(const Vector& rStressVector, double C, double Phi); diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_stress_strain_utitlities.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_stress_strain_utitlities.cpp index 40745b0cfd7f..fe6b32a2486e 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_stress_strain_utitlities.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_stress_strain_utitlities.cpp @@ -8,6 +8,7 @@ // License: geo_mechanics_application/license.txt // // Main authors: Wijtze Pieter Kikstra +// Gennady Markelov // #include "custom_utilities/math_utilities.h" @@ -118,4 +119,75 @@ KRATOS_TEST_CASE_IN_SUITE(CheckCalculateMohrCoulombPressureCapacityShearOnly, Kr 1.E-10); } +KRATOS_TEST_CASE_IN_SUITE(CheckCalculateCauchyStrain, KratosGeoMechanicsFastSuite) +{ + Matrix B(5, 5); + // clang-format off + B <<=1.0, 2.0, 3.0, 4.0, 5.0, + 1.1, 1.2, 1.3, 1.4, 1.5, + 1.0, 0.9, 0.8, 0.7, 0.6, + 0.0, 1.0, 2.0, 3.0, 4.0, + 5.0, 4.0, 3.0, 2.0, 1.0; + // clang-format on + + Vector displacements(5); + displacements <<= 0.01, 0.02, 0.03, 0.04, 0.05; + + const auto strain = StressStrainUtilities::CalculateCauchyStrain(B, displacements); + + Vector expected_strain(5); + expected_strain <<= 0.55, 0.205, 0.11, 0.4, 0.35; + + KRATOS_EXPECT_VECTOR_NEAR(strain, expected_strain, 1.E-10); +} + +KRATOS_TEST_CASE_IN_SUITE(CheckCalculateStrains, KratosGeoMechanicsFastSuite) +{ + Matrix B(4, 4); + // clang-format off + B <<=1.0, 2.0, 3.0, 4.0, + 1.1, 1.2, 1.3, 1.4, + 1.0, 0.9, 0.8, 0.7, + 0.0, 1.0, 2.0, 3.0; + // clang-format on + + Vector displacements(4); + displacements <<= 0.01, 0.02, 0.03, 0.04; + + bool use_hencky_strain = false; + std::size_t voigt_size = 4; + + std::vector Bs; + Bs.push_back(B); + Bs.push_back(B); + + std::vector deformation_gradients; + Matrix deformation_gradient = std::sqrt(3.) * IdentityMatrix(3); + deformation_gradients.push_back(deformation_gradient); + deformation_gradients.push_back(deformation_gradient); + + auto strains = StressStrainUtilities::CalculateStrains(deformation_gradients, Bs, displacements, + use_hencky_strain, voigt_size); + Vector expected_strain(4); + expected_strain <<= 0.3, 0.13, 0.08, 0.2; + std::vector expected_strains; + expected_strains.push_back(expected_strain); + expected_strains.push_back(expected_strain); + + for (size_t i = 0; i < strains.size(); ++i) + KRATOS_EXPECT_VECTOR_NEAR(strains[i], expected_strains[i], 1.E-10); + + use_hencky_strain = true; + strains = StressStrainUtilities::CalculateStrains(deformation_gradients, Bs, displacements, + use_hencky_strain, voigt_size); + + expected_strain <<= 0.549306, 0.549306, 0.549306, 0; + expected_strains.clear(); + expected_strains.push_back(expected_strain); + expected_strains.push_back(expected_strain); + + for (size_t i = 0; i < strains.size(); ++i) + KRATOS_EXPECT_VECTOR_NEAR(strains[i], expected_strains[i], 1.E-6); +} + } // namespace Kratos::Testing \ No newline at end of file