From ac6f95d64572e2648172f1c23b19ccd271d61864 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Fri, 3 May 2024 16:45:19 +0200 Subject: [PATCH] Clean up calculation of deformation gradients (#12341) Reduced the use of `ElementVariables` in the context of calculating the deformation gradients. This improves local reasoning. Furthermore, the code was simplified significantly, which eases understanding. Other changes include: - Changed the signatures of members `CalculateDeformationGradient`: * They are no longer `virtual`, since no overrides were defined. * The return type has been changed to `Matrix` to return the result. * They no longer need a reference to data structure `ElementVariables`, since the result is now returned using the return type rather than a member of the output argument. * They have been made `const`, since no other members need to be modified. - Members `CalculateDeformationGradient` now only compute the deformation gradient itself, and no longer its determinant as well. To be on the safe side, the determinant is computed explicitly in a few places where it was not obvious whether the determinant is being used or not. - In several places, the use of `ElementVariables` had become redundant as a result of simplifying the code (i.e. after eliminating usages of members of `ElementVariables`). In most cases, we could use temporaries. - There is no need to explicitly check sizes before resizing a container: the `resize` operation of the container already takes care of that. - Removed two redundant calls to member `CalculateKinematics` (which only became apparent after carrying out the above changes). - Removed several comments that had no additional value. --- .../U_Pw_small_strain_element.cpp | 10 ++-- .../U_Pw_small_strain_element.hpp | 2 +- .../U_Pw_updated_lagrangian_FIC_element.cpp | 43 +++-------------- .../U_Pw_updated_lagrangian_element.cpp | 44 +++--------------- .../small_strain_U_Pw_diff_order_element.cpp | 9 ++-- .../small_strain_U_Pw_diff_order_element.hpp | 2 +- ...ted_lagrangian_U_Pw_diff_order_element.cpp | 46 +++---------------- 7 files changed, 28 insertions(+), 128 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp index 4e6aff89455a..dfd1bb4acdc5 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -1505,7 +1505,8 @@ template void UPwSmallStrainElement::CalculateStrain(ElementVariables& rVariables, unsigned int GPoint) { if (rVariables.UseHenckyStrain) { - this->CalculateDeformationGradient(rVariables, GPoint); + rVariables.F = this->CalculateDeformationGradient(GPoint); + rVariables.detF = MathUtils<>::Det(rVariables.F); noalias(rVariables.StrainVector) = StressStrainUtilities::CalculateHenckyStrain(rVariables.F, VoigtSize); } else { this->CalculateCauchyStrain(rVariables); @@ -1525,8 +1526,7 @@ Vector UPwSmallStrainElement::CalculateGreenLagrangeStrain(cons } template -void UPwSmallStrainElement::CalculateDeformationGradient(ElementVariables& rVariables, - unsigned int GPoint) +Matrix UPwSmallStrainElement::CalculateDeformationGradient(unsigned int GPoint) const { KRATOS_TRY @@ -1544,9 +1544,7 @@ void UPwSmallStrainElement::CalculateDeformationGradient(Elemen KRATOS_ERROR_IF(detJ < 0.0) << "ERROR:: ELEMENT ID: " << this->Id() << " INVERTED. DETJ: " << detJ << " nodes:" << this->GetGeometry() << std::endl; - // Deformation gradient - noalias(rVariables.F) = prod(J, InvJ0); - rVariables.detF = MathUtils::Det(rVariables.F); + return prod(J, InvJ0); KRATOS_CATCH("") } diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp index 3217c3b4527e..587ec3af2afc 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -275,7 +275,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa virtual Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient); virtual void CalculateCauchyStrain(ElementVariables& rVariables); virtual void CalculateStrain(ElementVariables& rVariables, unsigned int GPoint); - virtual void CalculateDeformationGradient(ElementVariables& rVariables, unsigned int GPoint); + Matrix CalculateDeformationGradient(unsigned int GPoint) const; void InitializeNodalDisplacementVariables(ElementVariables& rVariables); void InitializeNodalPorePressureVariables(ElementVariables& rVariables); 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 19ce40174713..f54b3c335868 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 @@ -140,16 +140,9 @@ void UPwUpdatedLagrangianFICElement::CalculateOnIntegrationPoin const Variable& rVariable, std::vector& rOutput, const ProcessInfo& rCurrentProcessInfo) { if (rVariable == REFERENCE_DEFORMATION_GRADIENT_DETERMINANT) { - if (rOutput.size() != mConstitutiveLawVector.size()) - rOutput.resize(mConstitutiveLawVector.size()); - - ElementVariables Variables; - this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - // Loop over integration points + rOutput.clear(); for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { - this->CalculateDeformationGradient(Variables, GPoint); - rOutput[GPoint] = Variables.detF; + rOutput.emplace_back(MathUtils<>::Det(this->CalculateDeformationGradient(GPoint))); } } else { UPwSmallStrainFICElement::CalculateOnIntegrationPoints( @@ -162,40 +155,16 @@ template void UPwUpdatedLagrangianFICElement::CalculateOnIntegrationPoints( const Variable& rVariable, std::vector& rOutput, const ProcessInfo& rCurrentProcessInfo) { - // Defining necessary variables - const GeometryType& rGeom = this->GetGeometry(); - const unsigned int NumGPoints = rGeom.IntegrationPointsNumber(mThisIntegrationMethod); - - if (rOutput.size() != NumGPoints) rOutput.resize(NumGPoints); + rOutput.resize(this->GetGeometry().IntegrationPointsNumber(mThisIntegrationMethod)); if (rVariable == REFERENCE_DEFORMATION_GRADIENT) { - ElementVariables Variables; - this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { - this->CalculateDeformationGradient(Variables, GPoint); - - if (rOutput[GPoint].size2() != TDim) rOutput[GPoint].resize(TDim, TDim, false); - rOutput[GPoint] = Variables.F; + rOutput[GPoint] = this->CalculateDeformationGradient(GPoint); } } else if (rVariable == GREEN_LAGRANGE_STRAIN_TENSOR) { - // Definition of variables - ElementVariables Variables; - this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { - // compute element kinematics (Np, gradNpT, |J|, B, strains) - this->CalculateKinematics(Variables, GPoint); - - // Compute strain - this->CalculateDeformationGradient(Variables, GPoint); - Variables.StrainVector = this->CalculateGreenLagrangeStrain(Variables.F); - - if (rOutput[GPoint].size2() != TDim) rOutput[GPoint].resize(TDim, TDim, false); - - rOutput[GPoint] = MathUtils::StrainVectorToTensor(Variables.StrainVector); + rOutput[GPoint] = MathUtils<>::StrainVectorToTensor( + this->CalculateGreenLagrangeStrain(this->CalculateDeformationGradient(GPoint))); } } else { UPwSmallStrainFICElement::CalculateOnIntegrationPoints( 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 a8b0a3114d55..d6b45cb3a3f9 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp @@ -143,16 +143,9 @@ void UPwUpdatedLagrangianElement::CalculateOnIntegrationPoints( const Variable& rVariable, std::vector& rOutput, const ProcessInfo& rCurrentProcessInfo) { if (rVariable == REFERENCE_DEFORMATION_GRADIENT_DETERMINANT) { - if (rOutput.size() != mConstitutiveLawVector.size()) - rOutput.resize(mConstitutiveLawVector.size()); - - ElementVariables Variables; - this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - // Loop over integration points + rOutput.clear(); for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { - this->CalculateDeformationGradient(Variables, GPoint); - rOutput[GPoint] = Variables.detF; + rOutput.emplace_back(MathUtils<>::Det(this->CalculateDeformationGradient(GPoint))); } } else { UPwSmallStrainElement::CalculateOnIntegrationPoints(rVariable, rOutput, rCurrentProcessInfo); @@ -164,41 +157,16 @@ template void UPwUpdatedLagrangianElement::CalculateOnIntegrationPoints( const Variable& rVariable, std::vector& rOutput, const ProcessInfo& rCurrentProcessInfo) { - // Defining necessary variables - const GeometryType& rGeom = this->GetGeometry(); - const unsigned int NumGPoints = rGeom.IntegrationPointsNumber(mThisIntegrationMethod); - - if (rOutput.size() != NumGPoints) rOutput.resize(NumGPoints); + rOutput.resize(this->GetGeometry().IntegrationPointsNumber(mThisIntegrationMethod)); if (rVariable == REFERENCE_DEFORMATION_GRADIENT) { - ElementVariables Variables; - this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { - this->CalculateDeformationGradient(Variables, GPoint); - - if (rOutput[GPoint].size2() != TDim) rOutput[GPoint].resize(TDim, TDim, false); - - rOutput[GPoint] = Variables.F; + rOutput[GPoint] = this->CalculateDeformationGradient(GPoint); } } else if (rVariable == GREEN_LAGRANGE_STRAIN_TENSOR) { - // Definition of variables - ElementVariables Variables; - this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { - // compute element kinematics (Np, gradNpT, |J|, B, strains) - this->CalculateKinematics(Variables, GPoint); - - // Compute strain - this->CalculateDeformationGradient(Variables, GPoint); - Variables.StrainVector = this->CalculateGreenLagrangeStrain(Variables.F); - - if (rOutput[GPoint].size2() != TDim) rOutput[GPoint].resize(TDim, TDim, false); - - rOutput[GPoint] = MathUtils::StrainVectorToTensor(Variables.StrainVector); + rOutput[GPoint] = MathUtils<>::StrainVectorToTensor( + this->CalculateGreenLagrangeStrain(this->CalculateDeformationGradient(GPoint))); } } else { UPwSmallStrainElement::CalculateOnIntegrationPoints(rVariable, rOutput, rCurrentProcessInfo); 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 e54c98bb91c3..eaa6031801ed 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 @@ -2042,7 +2042,8 @@ GeometryData::IntegrationMethod SmallStrainUPwDiffOrderElement::GetIntegrationMe void SmallStrainUPwDiffOrderElement::CalculateStrain(ElementVariables& rVariables, unsigned int GPoint) { if (rVariables.UseHenckyStrain) { - this->CalculateDeformationGradient(rVariables, GPoint); + rVariables.F = this->CalculateDeformationGradient(GPoint); + rVariables.detF = MathUtils<>::Det(rVariables.F); const SizeType Dim = GetGeometry().WorkingSpaceDimension(); const SizeType VoigtSize = (Dim == N_DIM_3D ? VOIGT_SIZE_3D : VOIGT_SIZE_2D_PLANE_STRAIN); noalias(rVariables.StrainVector) = StressStrainUtilities::CalculateHenckyStrain(rVariables.F, VoigtSize); @@ -2061,7 +2062,7 @@ Vector SmallStrainUPwDiffOrderElement::CalculateGreenLagrangeStrain(const Matrix return mpStressStatePolicy->CalculateGreenLagrangeStrain(rDeformationGradient); } -void SmallStrainUPwDiffOrderElement::CalculateDeformationGradient(ElementVariables& rVariables, unsigned int GPoint) +Matrix SmallStrainUPwDiffOrderElement::CalculateDeformationGradient(unsigned int GPoint) const { KRATOS_TRY @@ -2085,9 +2086,7 @@ void SmallStrainUPwDiffOrderElement::CalculateDeformationGradient(ElementVariabl "mesh size." << std::endl; - // Deformation gradient - noalias(rVariables.F) = prod(J, InvJ0); - rVariables.detF = MathUtils::Det(rVariables.F); + return prod(J, InvJ0); 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 3b6e4f43d2fa..a767eaa622d8 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 @@ -295,7 +295,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub virtual void CalculateCauchyStrain(ElementVariables& rVariables); virtual void CalculateStrain(ElementVariables& rVariables, unsigned int GPoint); - virtual void CalculateDeformationGradient(ElementVariables& rVariables, unsigned int GPoint); + Matrix CalculateDeformationGradient(unsigned int GPoint) const; double CalculateFluidPressure(const ElementVariables& rVariables) 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 8a890993e8b9..37b5fbd59d86 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 @@ -151,16 +151,9 @@ void UpdatedLagrangianUPwDiffOrderElement::CalculateOnIntegrationPoints(const Va const ProcessInfo& rCurrentProcessInfo) { if (rVariable == REFERENCE_DEFORMATION_GRADIENT_DETERMINANT) { - if (rOutput.size() != mConstitutiveLawVector.size()) - rOutput.resize(mConstitutiveLawVector.size()); - - ElementVariables Variables; - this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - // Loop over integration points + rOutput.clear(); for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { - this->CalculateDeformationGradient(Variables, GPoint); - rOutput[GPoint] = Variables.detF; + rOutput.emplace_back(MathUtils<>::Det(this->CalculateDeformationGradient(GPoint))); } } else { SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(rVariable, rOutput, rCurrentProcessInfo); @@ -174,26 +167,11 @@ void UpdatedLagrangianUPwDiffOrderElement::CalculateOnIntegrationPoints(const Va { KRATOS_TRY - const GeometryType& rGeom = GetGeometry(); - const unsigned int& IntegrationPointsNumber = - rGeom.IntegrationPointsNumber(this->GetIntegrationMethod()); - - if (rOutput.size() != IntegrationPointsNumber) rOutput.resize(IntegrationPointsNumber); + rOutput.resize(this->GetGeometry().IntegrationPointsNumber(this->GetIntegrationMethod())); if (rVariable == GREEN_LAGRANGE_STRAIN_VECTOR) { - // Definition of variables - ElementVariables Variables; - this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { - this->CalculateDeformationGradient(Variables, GPoint); - Variables.StrainVector = this->CalculateGreenLagrangeStrain(Variables.F); - - if (rOutput[GPoint].size() != Variables.StrainVector.size()) - rOutput[GPoint].resize(Variables.StrainVector.size(), false); - - rOutput[GPoint] = Variables.StrainVector; + rOutput[GPoint] = this->CalculateGreenLagrangeStrain(this->CalculateDeformationGradient(GPoint)); } } else { SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(rVariable, rOutput, rCurrentProcessInfo); @@ -207,23 +185,11 @@ void UpdatedLagrangianUPwDiffOrderElement::CalculateOnIntegrationPoints(const Va std::vector& rOutput, const ProcessInfo& rCurrentProcessInfo) { - const GeometryType& rGeom = GetGeometry(); - const SizeType Dim = rGeom.WorkingSpaceDimension(); - - if (rOutput.size() != mConstitutiveLawVector.size()) - rOutput.resize(mConstitutiveLawVector.size()); + rOutput.resize(mConstitutiveLawVector.size()); if (rVariable == REFERENCE_DEFORMATION_GRADIENT) { - ElementVariables Variables; - this->InitializeElementVariables(Variables, rCurrentProcessInfo); - - // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { - this->CalculateDeformationGradient(Variables, GPoint); - - if (rOutput[GPoint].size2() != Dim) rOutput[GPoint].resize(Dim, Dim, false); - - rOutput[GPoint] = Variables.F; + rOutput[GPoint] = this->CalculateDeformationGradient(GPoint); } } else { SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(rVariable, rOutput, rCurrentProcessInfo);