Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[GeoMechanicsApplication] Isolated Hencky strain computation, such that it appears once. #11667

Merged
merged 3 commits into from
Oct 6, 2023
Merged

[GeoMechanicsApplication] Isolated Hencky strain computation, such that it appears once. #11667

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0615c6e
Isolated Hencky strain computation, such that it appears once. Made s…
WPK4FEM Oct 6, 2023
03dc2ff
Review comments from Richard processed
WPK4FEM Oct 6, 2023
44198b1
Review comments by Anne taken up.
WPK4FEM Oct 6, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 1 addition & 45 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1759,7 +1759,7 @@ void UPwSmallStrainElement<TDim,TNumNodes>::
{
if (rVariables.UseHenckyStrain) {
this->CalculateDeformationGradient(rVariables, GPoint);
this->CalculateHenckyStrain( rVariables );
noalias(rVariables.StrainVector) = StressStrainUtilities::CalculateHenckyStrain(rVariables.F, VoigtSize);
WPK4FEM marked this conversation as resolved.
Show resolved Hide resolved
} else {
this->CalculateCauchyStrain( rVariables );
}
Expand Down Expand Up @@ -1838,50 +1838,6 @@ void UPwSmallStrainElement<TDim,TNumNodes>::

}

//----------------------------------------------------------------------------------------
template< unsigned int TDim, unsigned int TNumNodes >
void UPwSmallStrainElement<TDim,TNumNodes>::
CalculateHenckyStrain( ElementVariables& rVariables )
{
KRATOS_TRY

//-Compute total deformation gradient
const Matrix& F = rVariables.F;

Matrix CMatrix;
CMatrix = prod(trans(F), F);

// Declare the different matrix
Matrix EigenValuesMatrix = ZeroMatrix(TDim, TDim);
Matrix EigenVectorsMatrix = ZeroMatrix(TDim, TDim);

// Decompose matrix
MathUtils<double>::GaussSeidelEigenSystem(CMatrix, EigenVectorsMatrix, EigenValuesMatrix, 1.0e-16, 20);

// Calculate the eigenvalues of the E matrix
for (IndexType i = 0; i < TDim; ++i) {
EigenValuesMatrix(i, i) = 0.5 * std::log(EigenValuesMatrix(i, i));
}

// Calculate E matrix
Matrix ETensor = ZeroMatrix(TDim, TDim);
MathUtils<double>::BDBtProductOperation(ETensor, EigenValuesMatrix, EigenVectorsMatrix);

// Hencky Strain Calculation
if constexpr (TDim==2) {
Vector StrainVector;
StrainVector = MathUtils<double>::StrainTensorToVector(ETensor);
rVariables.StrainVector[INDEX_2D_PLANE_STRAIN_XX] = StrainVector[0];
rVariables.StrainVector[INDEX_2D_PLANE_STRAIN_YY] = StrainVector[1];
rVariables.StrainVector[INDEX_2D_PLANE_STRAIN_ZZ] = 0.0;
rVariables.StrainVector[INDEX_2D_PLANE_STRAIN_XY] = StrainVector[2];
} else {
noalias(rVariables.StrainVector) = MathUtils<double>::StrainTensorToVector(ETensor);
}

KRATOS_CATCH( "" )
}

//----------------------------------------------------------------------------------------
template< unsigned int TDim, unsigned int TNumNodes >
void UPwSmallStrainElement<TDim,TNumNodes>::
Expand Down
1 change: 0 additions & 1 deletion applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -311,7 +311,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement :
virtual void CalculateCauchyAlmansiStrain( ElementVariables &rVariables );
virtual void CalculateCauchyGreenStrain( ElementVariables &rVariables );
virtual void CalculateCauchyStrain( ElementVariables &rVariables );
virtual void CalculateHenckyStrain( ElementVariables& rVariables );
virtual void CalculateStrain( ElementVariables &rVariables, unsigned int GPoint );
virtual void CalculateDeformationGradient( ElementVariables& rVariables,
unsigned int GPoint );
Expand Down
51 changes: 3 additions & 48 deletions ...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2451,7 +2451,6 @@ GeometryData::IntegrationMethod
break;
}

//return GeometryData::IntegrationMethod::GI_GAUSS_2;
return GI_GAUSS;
}

Expand All @@ -2461,7 +2460,9 @@ void SmallStrainUPwDiffOrderElement::
{
if (rVariables.UseHenckyStrain) {
this->CalculateDeformationGradient(rVariables, GPoint);
this->CalculateHenckyStrain( rVariables );
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);
} else {
this->CalculateCauchyStrain( rVariables );
}
Expand Down Expand Up @@ -2585,52 +2586,6 @@ void SmallStrainUPwDiffOrderElement::
KRATOS_CATCH( "" )
}

//----------------------------------------------------------------------------------------
void SmallStrainUPwDiffOrderElement::
CalculateHenckyStrain( ElementVariables& rVariables )
{
KRATOS_TRY

const GeometryType& rGeom = GetGeometry();
const SizeType Dim = rGeom.WorkingSpaceDimension();

//-Compute total deformation gradient
const Matrix& F = rVariables.F;

Matrix CMatrix;
CMatrix = prod(trans(F), F);

// Declare the different matrix
Matrix EigenValuesMatrix = ZeroMatrix(Dim, Dim);
Matrix EigenVectorsMatrix = ZeroMatrix(Dim, Dim);

// Decompose matrix
MathUtils<double>::GaussSeidelEigenSystem(CMatrix, EigenVectorsMatrix, EigenValuesMatrix, 1.0e-16, 20);

// Calculate the eigenvalues of the E matrix
for (IndexType i = 0; i < Dim; ++i) {
EigenValuesMatrix(i, i) = 0.5 * std::log(EigenValuesMatrix(i, i));
}

// Calculate E matrix
Matrix ETensor = ZeroMatrix(Dim, Dim);
MathUtils<double>::BDBtProductOperation(ETensor, EigenValuesMatrix, EigenVectorsMatrix);

// Hencky Strain Calculation
if (Dim==2) {
Vector StrainVector;
StrainVector = MathUtils<double>::StrainTensorToVector(ETensor);
rVariables.StrainVector[INDEX_2D_PLANE_STRAIN_XX] = StrainVector[0];
rVariables.StrainVector[INDEX_2D_PLANE_STRAIN_YY] = StrainVector[1];
rVariables.StrainVector[INDEX_2D_PLANE_STRAIN_ZZ] = 0.0;
rVariables.StrainVector[INDEX_2D_PLANE_STRAIN_XY] = StrainVector[2];
} else {
noalias(rVariables.StrainVector) = MathUtils<double>::StrainTensorToVector(ETensor);
}

KRATOS_CATCH( "" )
}

//----------------------------------------------------------------------------------------
void SmallStrainUPwDiffOrderElement::
CalculateJacobianOnCurrentConfiguration(double& detJ,
Expand Down
1 change: 0 additions & 1 deletion ...ications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -309,7 +309,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub
virtual void CalculateCauchyGreenStrain( ElementVariables& rVariables );
virtual void CalculateCauchyStrain( ElementVariables& rVariables );
virtual void CalculateStrain( ElementVariables& rVariables, unsigned int GPoint );
virtual void CalculateHenckyStrain( ElementVariables& rVariables );

virtual void CalculateDeformationGradient( ElementVariables& rVariables,
unsigned int GPoint );
Expand Down
189 changes: 40 additions & 149 deletions applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,51 +15,14 @@

/* Project includes */
#include "custom_utilities/math_utilities.hpp"
#include "geo_mechanics_application_constants.h"

namespace Kratos
{

/**@name Kratos Globals */
WPK4FEM marked this conversation as resolved.
Show resolved Hide resolved
/*@{ */


/*@} */
/**@name Type Definitions */
/*@{ */

/*@} */


/**@name Enum's */
/*@{ */


/*@} */
/**@name Functions */
/*@{ */



/*@} */
/**@name Kratos Classes */
/*@{ */


class KRATOS_API(GEO_MECHANICS_APPLICATION) StressStrainUtilities
{
public:
/**@name Type Definitions */
/*@{ */


/*@} */
/**@name Life Cycle
*/
/*@{ */

/** Operators.
*/

static double CalculateStressNorm(const Vector& StressVector)
{
KRATOS_TRY
Expand All @@ -74,7 +37,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) StressStrainUtilities
}
return std::sqrt(StressNorm);

KRATOS_CATCH( "" )
KRATOS_CATCH("")
}

static double CalculateVonMisesStress(const Vector& StressVector)
Expand All @@ -87,7 +50,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) StressStrainUtilities
noalias(StressTensor) = ZeroMatrix(3,3);
for (std::size_t i=0; i < LocalStressTensor.size1(); ++i) {
for (std::size_t j=0; j < LocalStressTensor.size2(); ++j) {
StressTensor(i,j) = LocalStressTensor(i,j);
StressTensor(i,j) = LocalStressTensor(i,j);
}
}

Expand All @@ -100,7 +63,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) StressStrainUtilities

return std::sqrt(std::max(SigmaEquivalent, 0.));

KRATOS_CATCH( "" )
KRATOS_CATCH("")
}

static double CalculateTrace(const Vector& StressVector)
Expand All @@ -116,131 +79,59 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) StressStrainUtilities

return trace;

KRATOS_CATCH( "" )
KRATOS_CATCH("")
}

static double CalculateMeanStress(const Vector& StressVector)
{
KRATOS_TRY

return CalculateTrace(StressVector) / (StressVector.size() == 3 ? 2.0 : 3.0);

KRATOS_CATCH( "" )
KRATOS_CATCH("")
}

static double CalculateVonMisesStrain(const Vector& StrainVector)
{
KRATOS_TRY

return (2.0/3.0) * CalculateVonMisesStress(StrainVector);

KRATOS_CATCH( "" )
KRATOS_CATCH("")
}

/*@} */
/**@name Operations */
/*@{ */


/*@} */
/**@name Access */
/*@{ */


/*@} */
/**@name Inquiry */
/*@{ */


/*@} */
/**@name Friends */
/*@{ */


/*@} */

protected:
/**@name Protected static Member Variables */
/*@{ */


/*@} */
/**@name Protected member Variables */
/*@{ */


/*@} */
/**@name Protected Operators*/
/*@{ */


/*@} */
/**@name Protected Operations*/
/*@{ */


/*@} */
/**@name Protected Access */
/*@{ */


/*@} */
/**@name Protected Inquiry */
/*@{ */


/*@} */
/**@name Protected LifeCycle */
/*@{ */



/*@} */

private:
/**@name Static Member Variables */
/*@{ */


/*@} */
/**@name Member Variables */
/*@{ */

/*@} */
/**@name Private Operators*/
/*@{ */


/*@} */
/**@name Private Operations*/
/*@{ */


/*@} */
/**@name Private Access */
/*@{ */


/*@} */
/**@name Private Inquiry */
/*@{ */


/*@} */
/**@name Un accessible methods */
/*@{ */


/*@} */

}; /* Class StressStrainUtilities */
static Vector CalculateHenckyStrain(const Matrix& DeformationGradient, size_t VoigtSize)
{
KRATOS_TRY

/*@} */
// right Cauchy Green deformation tensor C
Matrix C = prod(trans(DeformationGradient), DeformationGradient);
// Eigenvalues of C matrix, so principal right Cauchy Green deformation tensor C
Matrix EigenValuesMatrix, EigenVectorsMatrix;
MathUtils<double>::GaussSeidelEigenSystem(C, EigenVectorsMatrix, EigenValuesMatrix, 1.0e-16, 20);
// Compute natural strain == Logarithmic strain == Hencky strain from principal strains
for (std::size_t i = 0; i < DeformationGradient.size1(); ++i){
EigenValuesMatrix(i,i) = 0.5 * std::log(EigenValuesMatrix(i,i));
}

/**@name Type Definitions */
/*@{ */
// Rotate from principal strains back to the used coordinate system
Matrix ETensor;
MathUtils<double>::BDBtProductOperation(ETensor, EigenValuesMatrix, EigenVectorsMatrix);

// From tensor to vector
if (DeformationGradient.size1()==2 && VoigtSize == 4) {
// Plane strain
Vector StrainVector2D;
StrainVector2D = MathUtils<double>::StrainTensorToVector(ETensor, 3);
Vector StrainVector(4);
StrainVector[INDEX_2D_PLANE_STRAIN_XX] = StrainVector2D[0];
StrainVector[INDEX_2D_PLANE_STRAIN_YY] = StrainVector2D[1];
StrainVector[INDEX_2D_PLANE_STRAIN_ZZ] = 0.0;
StrainVector[INDEX_2D_PLANE_STRAIN_XY] = StrainVector2D[2];
return StrainVector;
} else {
return MathUtils<double>::StrainTensorToVector(ETensor, VoigtSize);
}

KRATOS_CATCH("")
}

/*@} */
};

}
Loading
Loading