diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp index 3217117efb2f..7dab37ae7ba6 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp @@ -13,6 +13,8 @@ // Application includes #include "custom_elements/U_Pw_base_element.hpp" +#include "plane_strain_stress_state.h" +#include "three_dimensional_stress_state.h" namespace Kratos { @@ -700,7 +702,15 @@ double UPwBaseElement::CalculateIntegrationCoefficient( const GeometryType::IntegrationPointsArrayType& IntegrationPoints, unsigned int PointNumber, double detJ) { - return IntegrationPoints[PointNumber].Weight() * detJ; + std::unique_ptr p_stress_state_policy; + if constexpr (TDim == 2) { + p_stress_state_policy = std::make_unique(); + } else { + p_stress_state_policy = std::make_unique(); + } + + return p_stress_state_policy->CalculateIntegrationCoefficient(IntegrationPoints[PointNumber], + detJ, GetGeometry()); } //---------------------------------------------------------------------------------------- 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 515a0874f4aa..d1b4791cd8d0 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -13,6 +13,8 @@ // Application includes #include "custom_elements/U_Pw_small_strain_element.hpp" +#include "plane_strain_stress_state.h" +#include "three_dimensional_stress_state.h" namespace Kratos { @@ -1209,33 +1211,15 @@ void UPwSmallStrainElement::CalculateBMatrix(Matrix& rB, const { KRATOS_TRY - unsigned int index; - if constexpr (TDim > 2) { - for (unsigned int i = 0; i < TNumNodes; ++i) { - index = TDim * i; - - rB(INDEX_3D_XX, index + INDEX_X) = GradNpT(i, INDEX_X); - rB(INDEX_3D_YY, index + INDEX_Y) = GradNpT(i, INDEX_Y); - rB(INDEX_3D_ZZ, index + INDEX_Z) = GradNpT(i, INDEX_Z); - rB(INDEX_3D_XY, index + INDEX_X) = GradNpT(i, INDEX_Y); - rB(INDEX_3D_XY, index + INDEX_Y) = GradNpT(i, INDEX_X); - rB(INDEX_3D_YZ, index + INDEX_Y) = GradNpT(i, INDEX_Z); - rB(INDEX_3D_YZ, index + INDEX_Z) = GradNpT(i, INDEX_Y); - rB(INDEX_3D_XZ, index + INDEX_X) = GradNpT(i, INDEX_Z); - rB(INDEX_3D_XZ, index + INDEX_Z) = GradNpT(i, INDEX_X); - } + std::unique_ptr p_stress_state_policy; + if constexpr (TDim == 2) { + p_stress_state_policy = std::make_unique(); } else { - // 2D plane strain - for (unsigned int i = 0; i < TNumNodes; ++i) { - index = TDim * i; - - rB(INDEX_2D_PLANE_STRAIN_XX, index + INDEX_X) = GradNpT(i, INDEX_X); - rB(INDEX_2D_PLANE_STRAIN_YY, index + INDEX_Y) = GradNpT(i, INDEX_Y); - rB(INDEX_2D_PLANE_STRAIN_XY, index + INDEX_X) = GradNpT(i, INDEX_Y); - rB(INDEX_2D_PLANE_STRAIN_XY, index + INDEX_Y) = GradNpT(i, INDEX_X); - } + p_stress_state_policy = std::make_unique(); } + rB = p_stress_state_policy->CalculateBMatrix(GradNpT, Np, this->GetGeometry()); + KRATOS_CATCH("") } @@ -1581,28 +1565,14 @@ Vector UPwSmallStrainElement::CalculateGreenLagrangeStrain(cons { KRATOS_TRY - Vector result = ZeroVector(VoigtSize); - - //-Compute total deformation gradient - Matrix ETensor; - ETensor = prod(trans(rDeformationGradient), rDeformationGradient); - - for (unsigned int i = 0; i < TDim; ++i) - ETensor(i, i) -= 1.0; - ETensor *= 0.5; - + std::unique_ptr p_stress_state_policy; if constexpr (TDim == 2) { - Vector StrainVector; - StrainVector = MathUtils::StrainTensorToVector(ETensor); - result[INDEX_2D_PLANE_STRAIN_XX] = StrainVector[0]; - result[INDEX_2D_PLANE_STRAIN_YY] = StrainVector[1]; - result[INDEX_2D_PLANE_STRAIN_ZZ] = 0.0; - result[INDEX_2D_PLANE_STRAIN_XY] = StrainVector[2]; + p_stress_state_policy = std::make_unique(); } else { - result = MathUtils::StrainTensorToVector(ETensor); + p_stress_state_policy = std::make_unique(); } - return result; + return p_stress_state_policy->CalculateGreenLagrangeStrain(rDeformationGradient); 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 d79e0bfa4529..e9b92f7f1d5f 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -278,9 +278,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa double CalculateBiotCoefficient(const ElementVariables& rVariables, bool hasBiotCoefficient) const; 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); + virtual void CalculateCauchyStrain(ElementVariables& rVariables); + virtual void CalculateStrain(ElementVariables& rVariables, unsigned int GPoint); + virtual void CalculateDeformationGradient(ElementVariables& rVariables, unsigned int GPoint); void InitializeNodalDisplacementVariables(ElementVariables& rVariables); void InitializeNodalPorePressureVariables(ElementVariables& rVariables); diff --git a/applications/GeoMechanicsApplication/custom_elements/axisymmetric_stress_state.cpp b/applications/GeoMechanicsApplication/custom_elements/axisymmetric_stress_state.cpp index 0d974df9efb7..867a13e6bf3e 100644 --- a/applications/GeoMechanicsApplication/custom_elements/axisymmetric_stress_state.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/axisymmetric_stress_state.cpp @@ -54,12 +54,12 @@ double AxisymmetricStressState::CalculateIntegrationCoefficient(const Geometry AxisymmetricStressState::Clone() const +std::unique_ptr AxisymmetricStressState::Clone() const { return std::make_unique(); } -Vector AxisymmetricStressState::CalculateGreenLagrangeStrain(const Matrix& rTotalDeformationGradient) const +Vector AxisymmetricStressState::CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const { KRATOS_ERROR << "The calculation of Green Lagrange strain is not implemented for axisymmetric " "configurations.\n"; diff --git a/applications/GeoMechanicsApplication/custom_elements/axisymmetric_stress_state.h b/applications/GeoMechanicsApplication/custom_elements/axisymmetric_stress_state.h index 15f65cc99257..6754b9076bd6 100644 --- a/applications/GeoMechanicsApplication/custom_elements/axisymmetric_stress_state.h +++ b/applications/GeoMechanicsApplication/custom_elements/axisymmetric_stress_state.h @@ -27,8 +27,8 @@ class AxisymmetricStressState : public StressStatePolicy [[nodiscard]] Matrix CalculateBMatrix(const Matrix& rGradNpT, const Vector& rNp, const Geometry& rGeometry) const override; - [[nodiscard]] Vector CalculateGreenLagrangeStrain(const Matrix& rTotalDeformationGradient) const override; - [[nodiscard]] unique_ptr Clone() const override; + [[nodiscard]] Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const override; + [[nodiscard]] std::unique_ptr Clone() const override; }; } // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_elements/plane_strain_stress_state.cpp b/applications/GeoMechanicsApplication/custom_elements/plane_strain_stress_state.cpp new file mode 100644 index 000000000000..20479c13b70d --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_elements/plane_strain_stress_state.cpp @@ -0,0 +1,64 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Richard Faasse +// +#include "plane_strain_stress_state.h" +#include "custom_utilities/stress_strain_utilities.hpp" + +namespace Kratos +{ + +Matrix PlaneStrainStressState::CalculateBMatrix(const Matrix& rGradNpT, const Vector&, const Geometry& rGeometry) const +{ + const auto dimension = rGeometry.WorkingSpaceDimension(); + const auto number_of_nodes = rGeometry.size(); + Matrix result = ZeroMatrix(VOIGT_SIZE_2D_AXISYMMETRIC, dimension * number_of_nodes); + + for (unsigned int i = 0; i < number_of_nodes; ++i) { + const auto offset = dimension * i; + + result(INDEX_2D_PLANE_STRAIN_XX, offset + INDEX_X) = rGradNpT(i, INDEX_X); + result(INDEX_2D_PLANE_STRAIN_YY, offset + INDEX_Y) = rGradNpT(i, INDEX_Y); + result(INDEX_2D_PLANE_STRAIN_XY, offset + INDEX_X) = rGradNpT(i, INDEX_Y); + result(INDEX_2D_PLANE_STRAIN_XY, offset + INDEX_Y) = rGradNpT(i, INDEX_X); + } + + return result; +} + +double PlaneStrainStressState::CalculateIntegrationCoefficient(const Geometry::IntegrationPointType& rIntegrationPoint, + double DetJ, + const Geometry&) const +{ + return rIntegrationPoint.Weight() * DetJ; +} + +Vector PlaneStrainStressState::CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const +{ + return ConvertStrainTensorToVector(StressStrainUtilities::CalculateGreenLagrangeStrainTensor(rDeformationGradient)); +} + +std::unique_ptr PlaneStrainStressState::Clone() const +{ + return std::make_unique(); +} + +Vector PlaneStrainStressState::ConvertStrainTensorToVector(const Matrix& rStrainTensor) +{ + const auto strain_vector = MathUtils::StrainTensorToVector(rStrainTensor); + Vector result = ZeroVector(VOIGT_SIZE_2D_PLANE_STRAIN); + result[INDEX_2D_PLANE_STRAIN_XX] = strain_vector[0]; + result[INDEX_2D_PLANE_STRAIN_YY] = strain_vector[1]; + result[INDEX_2D_PLANE_STRAIN_ZZ] = 0.0; + result[INDEX_2D_PLANE_STRAIN_XY] = strain_vector[2]; + return result; +} + +} // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_elements/plane_strain_stress_state.h b/applications/GeoMechanicsApplication/custom_elements/plane_strain_stress_state.h new file mode 100644 index 000000000000..228e514a0283 --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_elements/plane_strain_stress_state.h @@ -0,0 +1,36 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Richard Faasse +// + +#pragma once + +#include "stress_state_policy.h" + +namespace Kratos +{ + +class PlaneStrainStressState : public StressStatePolicy +{ +public: + [[nodiscard]] Matrix CalculateBMatrix(const Matrix& rGradNpT, + const Vector& rNp, + const Geometry& rGeometry) const override; + [[nodiscard]] double CalculateIntegrationCoefficient(const Geometry::IntegrationPointType& rIntegrationPoint, + double DetJ, + const Geometry& rGeometry) const override; + [[nodiscard]] Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const override; + [[nodiscard]] std::unique_ptr Clone() const override; + +private: + [[nodiscard]] static Vector ConvertStrainTensorToVector(const Matrix& rStrainTensor); +}; + +} // namespace Kratos 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 60af48cadd28..f7c8bdd07c1e 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 @@ -21,6 +21,8 @@ // Application includes #include "custom_elements/small_strain_U_Pw_diff_order_element.hpp" #include "custom_utilities/element_utilities.hpp" +#include "plane_strain_stress_state.h" +#include "three_dimensional_stress_state.h" namespace Kratos { @@ -1797,38 +1799,15 @@ void SmallStrainUPwDiffOrderElement::CalculateBMatrix(Matrix& rB, const Matrix& { KRATOS_TRY - const GeometryType& rGeom = GetGeometry(); - const SizeType Dim = rGeom.WorkingSpaceDimension(); - const SizeType NumUNodes = rGeom.PointsNumber(); - - unsigned int index; - - if (Dim > 2) { - for (unsigned int i = 0; i < NumUNodes; ++i) { - index = Dim * i; - - rB(INDEX_3D_XX, index + INDEX_X) = DNp_DX(i, INDEX_X); - rB(INDEX_3D_YY, index + INDEX_Y) = DNp_DX(i, INDEX_Y); - rB(INDEX_3D_ZZ, index + INDEX_Z) = DNp_DX(i, INDEX_Z); - rB(INDEX_3D_XY, index + INDEX_X) = DNp_DX(i, INDEX_Y); - rB(INDEX_3D_XY, index + INDEX_Y) = DNp_DX(i, INDEX_X); - rB(INDEX_3D_YZ, index + INDEX_Y) = DNp_DX(i, INDEX_Z); - rB(INDEX_3D_YZ, index + INDEX_Z) = DNp_DX(i, INDEX_Y); - rB(INDEX_3D_XZ, index + INDEX_X) = DNp_DX(i, INDEX_Z); - rB(INDEX_3D_XZ, index + INDEX_Z) = DNp_DX(i, INDEX_X); - } + std::unique_ptr p_stress_state_policy; + if (GetGeometry().WorkingSpaceDimension() == 2) { + p_stress_state_policy = std::make_unique(); } else { - // 2D plane strain - for (unsigned int i = 0; i < NumUNodes; ++i) { - index = Dim * i; - - rB(INDEX_2D_PLANE_STRAIN_XX, index + INDEX_X) = DNp_DX(i, INDEX_X); - rB(INDEX_2D_PLANE_STRAIN_YY, index + INDEX_Y) = DNp_DX(i, INDEX_Y); - rB(INDEX_2D_PLANE_STRAIN_XY, index + INDEX_X) = DNp_DX(i, INDEX_Y); - rB(INDEX_2D_PLANE_STRAIN_XY, index + INDEX_Y) = DNp_DX(i, INDEX_X); - } + p_stress_state_policy = std::make_unique(); } + rB = p_stress_state_policy->CalculateBMatrix(DNp_DX, Np, this->GetGeometry()); + KRATOS_CATCH("") } @@ -1853,7 +1832,15 @@ void SmallStrainUPwDiffOrderElement::SetConstitutiveParameters(ElementVariables& double SmallStrainUPwDiffOrderElement::CalculateIntegrationCoefficient( const GeometryType::IntegrationPointsArrayType& IntegrationPoints, unsigned int GPoint, double detJ) { - return IntegrationPoints[GPoint].Weight() * detJ; + std::unique_ptr p_stress_state_policy; + if (GetGeometry().WorkingSpaceDimension() == 2) { + p_stress_state_policy = std::make_unique(); + } else { + p_stress_state_policy = std::make_unique(); + } + + return p_stress_state_policy->CalculateIntegrationCoefficient(IntegrationPoints[GPoint], detJ, + GetGeometry()); } void SmallStrainUPwDiffOrderElement::CalculateAndAddLHS(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables) @@ -2278,32 +2265,14 @@ Vector SmallStrainUPwDiffOrderElement::CalculateGreenLagrangeStrain(const Matrix { KRATOS_TRY - const GeometryType& rGeom = GetGeometry(); - const SizeType Dim = rGeom.WorkingSpaceDimension(); - - const SizeType VoigtSize = (Dim == 3 ? VOIGT_SIZE_3D : VOIGT_SIZE_2D_PLANE_STRAIN); - Vector result = ZeroVector(VoigtSize); - - //-Compute total deformation gradient - Matrix ETensor; - ETensor = prod(trans(rDeformationGradient), rDeformationGradient); - - for (unsigned int i = 0; i < Dim; ++i) - ETensor(i, i) -= 1.0; - ETensor *= 0.5; - - if (Dim == 2) { - Vector StrainVector; - StrainVector = MathUtils::StrainTensorToVector(ETensor); - result[INDEX_2D_PLANE_STRAIN_XX] = StrainVector[0]; - result[INDEX_2D_PLANE_STRAIN_YY] = StrainVector[1]; - result[INDEX_2D_PLANE_STRAIN_ZZ] = 0.0; - result[INDEX_2D_PLANE_STRAIN_XY] = StrainVector[2]; + std::unique_ptr p_stress_state_policy; + if (GetGeometry().WorkingSpaceDimension() == 2) { + p_stress_state_policy = std::make_unique(); } else { - result = MathUtils::StrainTensorToVector(ETensor); + p_stress_state_policy = std::make_unique(); } - return result; + return p_stress_state_policy->CalculateGreenLagrangeStrain(rDeformationGradient); 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 8645531fce71..9da39ae5c5e4 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 @@ -286,8 +286,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub void AssembleUBlockMatrix(Matrix& rLeftHandSideMatrix, const Matrix& StiffnessMatrix) const; virtual Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient); - virtual void CalculateCauchyStrain(ElementVariables& rVariables); - virtual void CalculateStrain(ElementVariables& rVariables, unsigned int GPoint); + virtual void CalculateCauchyStrain(ElementVariables& rVariables); + virtual void CalculateStrain(ElementVariables& rVariables, unsigned int GPoint); virtual void CalculateDeformationGradient(ElementVariables& rVariables, unsigned int GPoint); diff --git a/applications/GeoMechanicsApplication/custom_elements/stress_state_policy.h b/applications/GeoMechanicsApplication/custom_elements/stress_state_policy.h index 72ea67fdc7cf..4f98dfe31340 100644 --- a/applications/GeoMechanicsApplication/custom_elements/stress_state_policy.h +++ b/applications/GeoMechanicsApplication/custom_elements/stress_state_policy.h @@ -29,9 +29,9 @@ class StressStatePolicy const Vector& rNp, const Geometry& rGeometry) const = 0; [[nodiscard]] virtual double CalculateIntegrationCoefficient(const Geometry::IntegrationPointType& rIntegrationPoint, - double detJ, + double DetJ, const Geometry& rGeometry) const = 0; - [[nodiscard]] virtual Vector CalculateGreenLagrangeStrain(const Matrix& rTotalDeformationGradient) const = 0; + [[nodiscard]] virtual Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const = 0; [[nodiscard]] virtual std::unique_ptr Clone() const = 0; }; diff --git a/applications/GeoMechanicsApplication/custom_elements/three_dimensional_stress_state.cpp b/applications/GeoMechanicsApplication/custom_elements/three_dimensional_stress_state.cpp new file mode 100644 index 000000000000..387a798f7e31 --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_elements/three_dimensional_stress_state.cpp @@ -0,0 +1,61 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Marjan Fathian +// Richard Faasse +// +#include "three_dimensional_stress_state.h" +#include "custom_utilities/stress_strain_utilities.hpp" + +namespace Kratos +{ + +Matrix ThreeDimensionalStressState::CalculateBMatrix(const Matrix& rGradNpT, + const Vector&, + const Geometry& rGeometry) const +{ + const auto dimension = rGeometry.WorkingSpaceDimension(); + const auto number_of_nodes = rGeometry.size(); + Matrix result = ZeroMatrix(VOIGT_SIZE_3D, dimension * number_of_nodes); + + for (unsigned int i = 0; i < number_of_nodes; ++i) { + const auto offset = dimension * i; + + result(INDEX_3D_XX, offset + INDEX_X) = rGradNpT(i, INDEX_X); + result(INDEX_3D_YY, offset + INDEX_Y) = rGradNpT(i, INDEX_Y); + result(INDEX_3D_ZZ, offset + INDEX_Z) = rGradNpT(i, INDEX_Z); + result(INDEX_3D_XY, offset + INDEX_X) = rGradNpT(i, INDEX_Y); + result(INDEX_3D_XY, offset + INDEX_Y) = rGradNpT(i, INDEX_X); + result(INDEX_3D_YZ, offset + INDEX_Y) = rGradNpT(i, INDEX_Z); + result(INDEX_3D_YZ, offset + INDEX_Z) = rGradNpT(i, INDEX_Y); + result(INDEX_3D_XZ, offset + INDEX_X) = rGradNpT(i, INDEX_Z); + result(INDEX_3D_XZ, offset + INDEX_Z) = rGradNpT(i, INDEX_X); + } + + return result; +} + +double ThreeDimensionalStressState::CalculateIntegrationCoefficient( + const Geometry::IntegrationPointType& rIntegrationPoint, double DetJ, const Geometry&) const +{ + return rIntegrationPoint.Weight() * DetJ; +} + +Vector ThreeDimensionalStressState::CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const +{ + return MathUtils<>::StrainTensorToVector( + StressStrainUtilities::CalculateGreenLagrangeStrainTensor(rDeformationGradient)); +} + +std::unique_ptr ThreeDimensionalStressState::Clone() const +{ + return std::make_unique(); +} + +} // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_elements/three_dimensional_stress_state.h b/applications/GeoMechanicsApplication/custom_elements/three_dimensional_stress_state.h new file mode 100644 index 000000000000..9dd44bc6a5aa --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_elements/three_dimensional_stress_state.h @@ -0,0 +1,34 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Marjan Fathian +// Richard Faasse +// + +#pragma once + +#include "stress_state_policy.h" + +namespace Kratos +{ + +class ThreeDimensionalStressState : public StressStatePolicy +{ +public: + [[nodiscard]] Matrix CalculateBMatrix(const Matrix& rGradNpT, + const Vector&, + const Geometry& rGeometry) const override; + [[nodiscard]] double CalculateIntegrationCoefficient(const Geometry::IntegrationPointType& rIntegrationPoint, + double DetJ, + const Geometry& rGeometry) const override; + [[nodiscard]] Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const override; + [[nodiscard]] std::unique_ptr Clone() const override; +}; + +} // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.cpp b/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.cpp new file mode 100644 index 000000000000..438c085b5eb4 --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.cpp @@ -0,0 +1,24 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Richard Faasse +// + +#include "stress_strain_utilities.hpp" + +namespace Kratos +{ + +Matrix StressStrainUtilities::CalculateGreenLagrangeStrainTensor(const Matrix& rDeformationGradient) +{ + return 0.5 * (prod(trans(rDeformationGradient), rDeformationGradient) - + IdentityMatrix(rDeformationGradient.size1())); +} + +} diff --git a/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.hpp b/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.hpp index 3421d76f8528..ec00836e4904 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/stress_strain_utilities.hpp @@ -132,6 +132,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) StressStrainUtilities KRATOS_CATCH("") } + static Matrix CalculateGreenLagrangeStrainTensor(const Matrix& rDeformationGradient); }; } \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_axisymmetric_stress_state.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_axisymmetric_stress_state.cpp index 5f7e786dee3f..8e6bb7ccf8c4 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_axisymmetric_stress_state.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_axisymmetric_stress_state.cpp @@ -17,6 +17,7 @@ #include "geometries/geometry.h" #include "includes/checks.h" #include "testing/testing.h" +#include "tests/cpp_tests/test_utilities/model_setup_utilities.h" #include using namespace Kratos; @@ -24,27 +25,13 @@ using namespace Kratos; namespace Kratos::Testing { -ModelPart& CreateModelPart(Model& rModel) -{ - ModelPart& result = rModel.CreateModelPart("Main"); - - result.CreateNewNode(1, 0.0, 0.0, 0.0); - result.CreateNewNode(2, 1.0, 0.0, 0.0); - result.CreateNewNode(3, 1.0, 1.0, 0.0); - - std::vector node_ids{1, 2, 3}; - result.CreateNewElement("UPwSmallStrainElement2D3N", 1, node_ids, result.CreateNewProperties(0)); - - return result; -} - KRATOS_TEST_CASE_IN_SUITE(CalculateBMatrixWithValidGeometryReturnsCorrectResults, KratosGeoMechanicsFastSuite) { const std::unique_ptr p_stress_state_policy = std::make_unique(); Model model; - auto& model_part = CreateModelPart(model); + auto& r_model_part = ModelSetupUtilities::CreateModelPartWithASingle2D3NElement(model); Vector Np(3); Np <<= 1.0, 2.0, 3.0; @@ -57,7 +44,7 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateBMatrixWithValidGeometryReturnsCorrectResults // clang-format on const Matrix calculated_matrix = - p_stress_state_policy->CalculateBMatrix(GradNpT, Np, model_part.GetElement(1).GetGeometry()); + p_stress_state_policy->CalculateBMatrix(GradNpT, Np, r_model_part.GetElement(1).GetGeometry()); // clang-format off Matrix expected_matrix(4, 6); @@ -76,14 +63,14 @@ KRATOS_TEST_CASE_IN_SUITE(ReturnCorrectIntegrationCoefficient, KratosGeoMechanic std::make_unique(); Model model; - auto& model_part = CreateModelPart(model); + auto& r_model_part = ModelSetupUtilities::CreateModelPartWithASingle2D3NElement(model); // The shape function values for this integration point are 0.2, 0.5 and 0.3 for nodes 1, 2 and 3 respectively Geometry::IntegrationPointType integration_point(0.5, 0.3, 0.0, 0.5); const double detJ = 2.0; const double calculated_coefficient = p_stress_state_policy->CalculateIntegrationCoefficient( - integration_point, detJ, model_part.GetElement(1).GetGeometry()); + integration_point, detJ, r_model_part.GetElement(1).GetGeometry()); // The expected number is calculated as follows: // 2.0 * pi * 0.8 (radius) * 2.0 (detJ) * 0.5 (weight) = 5.02655 diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_plane_strain_stress_state.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_plane_strain_stress_state.cpp new file mode 100644 index 000000000000..646abcce33ce --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_plane_strain_stress_state.cpp @@ -0,0 +1,104 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Richard Faasse +// + +#include "containers/model.h" +#include "custom_elements/plane_strain_stress_state.h" +#include "includes/checks.h" +#include "testing/testing.h" +#include "tests/cpp_tests/test_utilities/model_setup_utilities.h" +#include + +using namespace Kratos; + +namespace Kratos::Testing +{ + +KRATOS_TEST_CASE_IN_SUITE(CalculateBMatrixGivesCorrectResults, KratosGeoMechanicsFastSuite) +{ + auto p_stress_state_policy = std::make_unique(); + + Model model; + auto& r_model_part = ModelSetupUtilities::CreateModelPartWithASingle2D3NElement(model); + + Vector Np(3); + Np <<= 1.0, 2.0, 3.0; + + // clang-format off + Matrix GradNpT(3, 2); + GradNpT <<= 1.0, 2.0, + 3.0, 4.0, + 5.0, 6.0; + // clang-format on + + const auto calculated_matrix = + p_stress_state_policy->CalculateBMatrix(GradNpT, Np, r_model_part.GetElement(1).GetGeometry()); + + // clang-format off + Matrix expected_matrix(4, 6); + expected_matrix <<= 1 ,0 ,3 ,0 ,5 ,0, // This row contains the first column of GradNpT on columns 1, 3 and 5 + 0 ,2 ,0 ,4 ,0 ,6, // This row contains the second column of GradNpT on columns 2, 4 and 6 + 0 ,0 ,0 ,0 ,0 ,0, // This row is not set, so remains 0. + 2 ,1 ,4 ,3 ,6 ,5; // This row contains the first and second columns of GradNpT, swapping x and y + // clang-format on + + KRATOS_CHECK_MATRIX_NEAR(calculated_matrix, expected_matrix, 1e-12) +} + +KRATOS_TEST_CASE_IN_SUITE(PlaneStrainStressState_ReturnsCorrectIntegrationCoefficient, KratosGeoMechanicsFastSuite) +{ + const auto p_stress_state_policy = std::make_unique(); + + Model model; + auto& r_model_part = ModelSetupUtilities::CreateModelPartWithASingle2D3NElement(model); + + // The shape function values for this integration point are 0.2, 0.5 and 0.3 for nodes 1, 2 and 3 respectively + Geometry::IntegrationPointType integration_point(0.5, 0.3, 0.0, 0.5); + + const auto detJ = 2.0; + const auto calculated_coefficient = p_stress_state_policy->CalculateIntegrationCoefficient( + integration_point, detJ, r_model_part.GetElement(1).GetGeometry()); + + // The expected number is calculated as follows: + // 2.0 (detJ) * 0.5 (weight) = 1.0 + KRATOS_EXPECT_NEAR(calculated_coefficient, 1.0, 1e-5); +} + +KRATOS_TEST_CASE_IN_SUITE(PlaneStrainStressState_ReturnsCorrectGreenLagrangeStrain, KratosGeoMechanicsFastSuite) +{ + const auto p_stress_state_policy = std::make_unique(); + + // clang-format off + Matrix deformation_gradient = ZeroMatrix(2, 2); + deformation_gradient <<= 1.0, 2.0, + 3.0, 4.0; + // clang-format on + + const auto calculated_strain = p_stress_state_policy->CalculateGreenLagrangeStrain(deformation_gradient); + + Vector expected_vector = ZeroVector(4); + expected_vector <<= 4.5, 9.5, 0, 14; + + // The expected strain is calculated as follows: + // 0.5 * (F^T * F - I) and then converted to a vector + KRATOS_CHECK_VECTOR_NEAR(calculated_strain, expected_vector, 1e-12) +} + +KRATOS_TEST_CASE_IN_SUITE(PlaneStrainStressState_GivesCorrectClone, KratosGeoMechanicsFastSuite) +{ + const std::unique_ptr p_stress_state_policy = + std::make_unique(); + const auto p_clone = p_stress_state_policy->Clone(); + + KRATOS_EXPECT_NE(dynamic_cast(p_clone.get()), nullptr); +} + +} // namespace Kratos::Testing \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_three_dimensional_stress_state.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_three_dimensional_stress_state.cpp new file mode 100644 index 000000000000..ddc425d8c2e9 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_three_dimensional_stress_state.cpp @@ -0,0 +1,115 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Marjan Fathian +// + +#include "containers/model.h" +#include "custom_elements/three_dimensional_stress_state.h" +#include "includes/checks.h" +#include "testing/testing.h" +#include "tests/cpp_tests/test_utilities/model_setup_utilities.h" +#include + +using namespace Kratos; + +namespace Kratos::Testing +{ + +KRATOS_TEST_CASE_IN_SUITE(ThreeDimensionalStressState_CalculateBMatrixReturnsCorrectResults, KratosGeoMechanicsFastSuite) +{ + auto p_stress_state_policy = std::make_unique(); + + Model model; + auto& r_model_part = ModelSetupUtilities::CreateModelPartWithASingle3D4NElement(model); + + Vector Np(4); + Np <<= 1.0, 2.0, 3.0, 4.0; + + // clang-format off + Matrix GradNpT(4, 3); + GradNpT <<= 1.0, 2.0, 3.0, + 4.0, 5.0, 6.0, + 7.0, 8.0, 9.0, + 10.0, 11.0, 12.0; + // clang-format on + + const auto calculated_matrix = + p_stress_state_policy->CalculateBMatrix(GradNpT, Np, r_model_part.GetElement(1).GetGeometry()); + + // clang-format off + Matrix expected_matrix(6, 12); + expected_matrix <<= // This row (INDEX_3D_XX) contains the first column of GradNpT as INDEX_X + 1, 0, 0, 4, 0, 0, 7, 0, 0, 10, 0, 0, + // This row (INDEX_3D_YY) contains the second column of GradNpT as INDEX_Y + 0, 2, 0, 0, 5, 0, 0, 8, 0, 0, 11, 0, + // This row (INDEX_3D_ZZ) contains the third column of GradNpT as INDEX_Z + 0, 0, 3, 0, 0, 6, 0, 0, 9, 0, 0, 12, + // This row (INDEX_3D_XY) contains the first and second columns of GradNpT as INDEX_Y and INDEX_X respectively + 2, 1, 0, 5, 4, 0, 8, 7, 0, 11, 10, 0, + // This row (INDEX_3D_YZ) contains the second and third column of GradNpT as INDEX_Z and INDEX_Y respectively + 0, 3, 2, 0, 6, 5, 0, 9, 8, 0, 12, 11, + // This row (INDEX_3D_XZ) contains the first and third column of GradNpT as INDEX_Z and INDEX_X respectively + 3, 0, 1, 6, 0, 4, 9, 0, 7, 12, 0, 10; + // clang-format on + + KRATOS_CHECK_MATRIX_NEAR(calculated_matrix, expected_matrix, 1e-12) +} + +KRATOS_TEST_CASE_IN_SUITE(ThreeDimensionalStressState_ReturnsCorrectIntegrationCoefficient, KratosGeoMechanicsFastSuite) +{ + const auto p_stress_state_policy = std::make_unique(); + + Model model; + auto& r_model_part = ModelSetupUtilities::CreateModelPartWithASingle3D4NElement(model); + + // The shape function values for this integration point are 0.2, 0.5 and 0.3 for nodes 1, 2 and 3 respectively + Geometry::IntegrationPointType integration_point(0.5, 0.3, 0.0, 0.5); + + const auto detJ = 2.0; + const auto calculated_coefficient = p_stress_state_policy->CalculateIntegrationCoefficient( + integration_point, detJ, r_model_part.GetElement(1).GetGeometry()); + + // The expected number is calculated as follows: + // 2.0 (detJ) * 0.5 (weight) = 1.0 + KRATOS_EXPECT_NEAR(calculated_coefficient, 1.0, 1e-5); +} + +KRATOS_TEST_CASE_IN_SUITE(ThreeDimensionalStressState_GivesCorrectClone, KratosGeoMechanicsFastSuite) +{ + const std::unique_ptr p_stress_state_policy = + std::make_unique(); + const auto p_clone = p_stress_state_policy->Clone(); + + KRATOS_EXPECT_NE(dynamic_cast(p_clone.get()), nullptr); +} + +KRATOS_TEST_CASE_IN_SUITE(ThreeDimensionalStressState_CalculateGreenLagrangeStrainReturnsCorrectResults, + KratosGeoMechanicsFastSuite) +{ + const auto p_stress_state_policy = std::make_unique(); + + // clang-format off + Matrix deformation_gradient = ZeroMatrix(3, 3); + deformation_gradient <<= 1.0, 2.0, 3.0, + 4.0, 5.0, 6.0, + 7.0, 8.0, 9.0; + // clang-format on + + // The expected strain is calculated as follows: + // 0.5 * (F^T * F - I) and then converted to a vector + Vector expected_vector = ZeroVector(6); + expected_vector <<= 32.5, 46, 62.5, 78, 108, 90; + + const auto calculated_vector = p_stress_state_policy->CalculateGreenLagrangeStrain(deformation_gradient); + + KRATOS_CHECK_VECTOR_NEAR(expected_vector, calculated_vector, 1e-12) +} + +} // namespace Kratos::Testing \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_utilities/model_setup_utilities.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_utilities/model_setup_utilities.cpp new file mode 100644 index 000000000000..a73eb6b46a61 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_utilities/model_setup_utilities.cpp @@ -0,0 +1,48 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Richard Faasse +// +#include "model_setup_utilities.h" +#include "includes/model_part.h" +#include "containers/model.h" + +namespace Kratos::Testing::ModelSetupUtilities +{ + +ModelPart& CreateModelPartWithASingle2D3NElement(Model& rModel) +{ + ModelPart& result = rModel.CreateModelPart("Main"); + + result.CreateNewNode(1, 0.0, 0.0, 0.0); + result.CreateNewNode(2, 1.0, 0.0, 0.0); + result.CreateNewNode(3, 1.0, 1.0, 0.0); + + const std::vector node_ids{1, 2, 3}; + result.CreateNewElement("UPwSmallStrainElement2D3N", 1, node_ids, result.CreateNewProperties(0)); + + return result; +} + +ModelPart& CreateModelPartWithASingle3D4NElement(Model& rModel) +{ + ModelPart& result = rModel.CreateModelPart("Main"); + + result.CreateNewNode(1, 0.0, 0.0, 0.0); + result.CreateNewNode(2, 1.0, 0.0, 0.0); + result.CreateNewNode(3, 0.0, 1.0, 0.0); + result.CreateNewNode(4, 0.0, 0.0, 1.0); + + const std::vector node_ids{1, 2, 3, 4}; + result.CreateNewElement("UPwSmallStrainElement3D4N", 1, node_ids, result.CreateNewProperties(0)); + + return result; +} + +} diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_utilities/model_setup_utilities.h b/applications/GeoMechanicsApplication/tests/cpp_tests/test_utilities/model_setup_utilities.h new file mode 100644 index 000000000000..efbe104a3d8a --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_utilities/model_setup_utilities.h @@ -0,0 +1,27 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Richard Faasse +// + +#pragma once + +namespace Kratos +{ +class Model; +class ModelPart; +} // namespace Kratos + +namespace Kratos::Testing::ModelSetupUtilities +{ + +ModelPart& CreateModelPartWithASingle2D3NElement(Model& rModel); +ModelPart& CreateModelPartWithASingle3D4NElement(Model& rModel); + +} // namespace Kratos::ModelSetupUtilities