Skip to content

Commit

Permalink
Merge branch 'master' into geo/extract-set-constitutive-parameters
Browse files Browse the repository at this point in the history 
# Conflicts:
#	applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
#	applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp
  • Loading branch information
Gennady Markelov committed May 15, 2024
2 parents 4fc57b8 + 085a2af commit a4cec35
Show file tree
Hide file tree
Showing 43 changed files with 886 additions and 1,165 deletions.
3 changes: 2 additions & 1 deletion applications/ConstitutiveLawsApplication/constitutive_laws_application_variables.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,8 @@ namespace Kratos
{
enum class SofteningType {Linear = 0, Exponential = 1, HardeningDamage = 2, CurveFittingDamage = 3};

enum class TangentOperatorEstimation {Analytic = 0, FirstOrderPerturbation = 1, SecondOrderPerturbation = 2, Secant = 3, SecondOrderPerturbationV2 = 4, InitialStiffness = 5};
enum class TangentOperatorEstimation {Analytic = 0, FirstOrderPerturbation = 1, SecondOrderPerturbation = 2,
Secant = 3, SecondOrderPerturbationV2 = 4, InitialStiffness = 5, OrthogonalSecant = 6};

///@name Functions
///@{
Expand Down
6 changes: 3 additions & 3 deletions ..._constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_kinematic_plasticity.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,7 +243,7 @@ class GenericConstitutiveLawIntegratorKinematicPlasticity
const double Denominator
)
{
rTangent = rElasticMatrix - outer_prod(Vector(prod(rElasticMatrix, rGFluxVector)), Vector(prod(rElasticMatrix, rFFluxVector))) * Denominator;
noalias(rTangent) = rElasticMatrix - outer_prod(Vector(prod(rElasticMatrix, rGFluxVector)), Vector(prod(rElasticMatrix, rFFluxVector))) * Denominator;
}


Expand Down Expand Up @@ -281,10 +281,10 @@ class GenericConstitutiveLawIntegratorKinematicPlasticity
{
BoundedArrayType deviator = ZeroVector(VoigtSize);
BoundedArrayType h_capa = ZeroVector(VoigtSize);
double J2, tensile_indicator_factor, compression_indicator_factor, slope, hardening_parameter, equivalent_plastic_strain;
double J2, tensile_indicator_factor, compression_indicator_factor, slope, hardening_parameter, equivalent_plastic_strain, I1;

YieldSurfaceType::CalculateEquivalentStress( rPredictiveStressVector, rStrainVector, rUniaxialStress, rValues);
const double I1 = rPredictiveStressVector[0] + rPredictiveStressVector[1] + rPredictiveStressVector[2];
AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateI1Invariant(rPredictiveStressVector, I1);
AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateJ2Invariant(rPredictiveStressVector, I1, deviator, J2);
CalculateDerivativeYieldSurface(rPredictiveStressVector, deviator, J2, rYieldSurfaceDerivative, rValues);
CalculateDerivativePlasticPotential(rPredictiveStressVector, deviator, J2, rDerivativePlasticPotential, rValues);
Expand Down
6 changes: 3 additions & 3 deletions ...ion/custom_constitutive/auxiliary_files/cl_integrators/generic_cl_integrator_plasticity.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,10 +245,10 @@ class GenericConstitutiveLawIntegratorPlasticity
{
array_1d<double, VoigtSize> deviator = ZeroVector(VoigtSize);
array_1d<double, VoigtSize> h_capa = ZeroVector(VoigtSize);
double J2, tensile_indicator_factor, compression_indicator_factor, slope, hardening_parameter, equivalent_plastic_strain;
double J2, tensile_indicator_factor, compression_indicator_factor, slope, hardening_parameter, equivalent_plastic_strain, I1;

YieldSurfaceType::CalculateEquivalentStress( rPredictiveStressVector, rStrainVector, rUniaxialStress, rValues);
const double I1 = rPredictiveStressVector[0] + rPredictiveStressVector[1] + rPredictiveStressVector[2];
AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateI1Invariant(rPredictiveStressVector, I1);
AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateJ2Invariant(rPredictiveStressVector, I1, deviator, J2);
CalculateFFluxVector(rPredictiveStressVector, deviator, J2, rFflux, rValues);
CalculateGFluxVector(rPredictiveStressVector, deviator, J2, rGflux, rValues);
Expand Down Expand Up @@ -278,7 +278,7 @@ class GenericConstitutiveLawIntegratorPlasticity
const double Denominator
)
{
rTangent = rElasticMatrix - outer_prod(Vector(prod(rElasticMatrix, rGFluxVector)), Vector(prod(rElasticMatrix, rFFluxVector))) * Denominator;
noalias(rTangent) = rElasticMatrix - outer_prod(Vector(prod(rElasticMatrix, rGFluxVector)), Vector(prod(rElasticMatrix, rFFluxVector))) * Denominator;
}


Expand Down
8 changes: 1 addition & 7 deletions ...tion/custom_constitutive/auxiliary_files/plastic_potentials/von_mises_plastic_potential.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,16 +116,10 @@ class VonMisesPlasticPotential
)
{
array_1d<double, VoigtSize> first_vector, second_vector, third_vector;

AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateFirstVector(first_vector);
AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateSecondVector(rDeviator, J2, second_vector);
AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateThirdVector(rDeviator, J2, third_vector);

const double c1 = 0.0;
const double c2 = std::sqrt(3.0);
const double c3 = 0.0;

noalias(rGFlux) = c1 * first_vector + c2 * second_vector + c3 * third_vector;
noalias(rGFlux) = c2 * second_vector;
}

/**
Expand Down
2 changes: 1 addition & 1 deletion ...ication/custom_constitutive/auxiliary_files/yield_surfaces/drucker_prager_yield_surface.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,7 @@ class DruckerPragerYieldSurface
AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateFirstVector(first_vector);
AdvancedConstitutiveLawUtilities<VoigtSize>::CalculateSecondVector(rDeviator, J2, second_vector);

const double friction_angle = r_material_properties[FRICTION_ANGLE] * Globals::Pi / 180.0;;
const double friction_angle = r_material_properties[FRICTION_ANGLE] * Globals::Pi / 180.0;
const double sin_phi = std::sin(friction_angle);
const double Root3 = std::sqrt(3.0);

Expand Down
11 changes: 7 additions & 4 deletions ...tom_constitutive/finite_strains/plasticity/generic_finite_strain_isotropic_plasticity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,9 +164,7 @@ void GenericFiniteStrainIsotropicPlasticity<TConstLawIntegratorType>::
noalias(r_integrated_stress_vector) = predictive_stress_vector;

if (r_constitutive_law_options.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
this->CalculateTangentTensor(rValues, ConstitutiveLaw::StressMeasure_Kirchhoff);
} else {
BaseType::CalculateElasticMatrix( r_constitutive_matrix, rValues);
this->CalculateTangentTensor(rValues, ConstitutiveLaw::StressMeasure_Kirchhoff, plastic_strain);
}
}
}
Expand Down Expand Up @@ -301,7 +299,8 @@ template<class TConstLawIntegratorType>
void GenericFiniteStrainIsotropicPlasticity<TConstLawIntegratorType>::
CalculateTangentTensor(
ConstitutiveLaw::Parameters& rValues,
const ConstitutiveLaw::StressMeasure& rStressMeasure
const ConstitutiveLaw::StressMeasure& rStressMeasure,
const Vector& rPlasticStrain
)
{
const Properties& r_material_properties = rValues.GetMaterialProperties();
Expand All @@ -317,6 +316,10 @@ void GenericFiniteStrainIsotropicPlasticity<TConstLawIntegratorType>::
} else if (tangent_operator_estimation == TangentOperatorEstimation::SecondOrderPerturbation) {
// Calculates the Tangent Constitutive Tensor by perturbation (second order)
TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, this, rStressMeasure, consider_perturbation_threshold, 2);
} else if (tangent_operator_estimation == TangentOperatorEstimation::Secant) {
const Vector num = prod(rValues.GetConstitutiveMatrix(), rPlasticStrain);
const double denom = inner_prod(rValues.GetStrainVector(), num);
noalias(rValues.GetConstitutiveMatrix()) -= outer_prod(num, num) / denom;
}
}

Expand Down
3 changes: 2 additions & 1 deletion ...ustom_constitutive/finite_strains/plasticity/generic_finite_strain_isotropic_plasticity.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -315,7 +315,8 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) GenericFiniteStrainIsotropicPlas
*/
void CalculateTangentTensor(
ConstitutiveLaw::Parameters &rValues,
const ConstitutiveLaw::StressMeasure& rStressMeasure = ConstitutiveLaw::StressMeasure_Cauchy
const ConstitutiveLaw::StressMeasure& rStressMeasure,
const Vector& rPlasticStrain
);

///@}
Expand Down
11 changes: 7 additions & 4 deletions ...tom_constitutive/finite_strains/plasticity/generic_finite_strain_kinematic_plasticity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,9 +170,7 @@ void GenericFiniteStrainKinematicPlasticity<TConstLawIntegratorType>::
noalias(r_integrated_stress_vector) = predictive_stress_vector;

if (r_constitutive_law_options.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
this->CalculateTangentTensor(rValues, ConstitutiveLaw::StressMeasure_Kirchhoff);
} else {
BaseType::CalculateElasticMatrix( r_constitutive_matrix, rValues);
this->CalculateTangentTensor(rValues, ConstitutiveLaw::StressMeasure_Kirchhoff, plastic_strain);
}
}
}
Expand Down Expand Up @@ -314,7 +312,8 @@ template<class TConstLawIntegratorType>
void GenericFiniteStrainKinematicPlasticity<TConstLawIntegratorType>::
CalculateTangentTensor(
ConstitutiveLaw::Parameters& rValues,
const ConstitutiveLaw::StressMeasure& rStressMeasure
const ConstitutiveLaw::StressMeasure& rStressMeasure,
const Vector& rPlasticStrain
)
{
const Properties& r_material_properties = rValues.GetMaterialProperties();
Expand All @@ -330,6 +329,10 @@ void GenericFiniteStrainKinematicPlasticity<TConstLawIntegratorType>::
} else if (tangent_operator_estimation == TangentOperatorEstimation::SecondOrderPerturbation) {
// Calculates the Tangent Constitutive Tensor by perturbation (second order)
TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, this, ConstitutiveLaw::StressMeasure_Cauchy, consider_perturbation_threshold, 2);
} else if (tangent_operator_estimation == TangentOperatorEstimation::Secant) {
const Vector num = prod(rValues.GetConstitutiveMatrix(), rPlasticStrain);
const double denom = inner_prod(rValues.GetStrainVector(), num);
noalias(rValues.GetConstitutiveMatrix()) -= outer_prod(num, num) / denom;
}
}

Expand Down
3 changes: 2 additions & 1 deletion ...ustom_constitutive/finite_strains/plasticity/generic_finite_strain_kinematic_plasticity.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -315,7 +315,8 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) GenericFiniteStrainKinematicPlas
*/
void CalculateTangentTensor(
ConstitutiveLaw::Parameters &rValues,
const ConstitutiveLaw::StressMeasure& rStressMeasure = ConstitutiveLaw::StressMeasure_Cauchy
const ConstitutiveLaw::StressMeasure& rStressMeasure,
const Vector& rPlasticStrain
);

///@}
Expand Down
3 changes: 2 additions & 1 deletion ...cation/custom_constitutive/small_strains/damage/generic_small_strain_isotropic_damage.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,6 +129,7 @@ void GenericSmallStrainIsotropicDamage<TConstLawIntegratorType>::CalculateMateri
noalias(r_integrated_stress_vector) = predictive_stress_vector;

if (r_constitutive_law_options.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
r_constitutive_matrix *= (1.0 - damage);
this->CalculateTangentTensor(rValues);
}
}
Expand Down Expand Up @@ -163,7 +164,7 @@ void GenericSmallStrainIsotropicDamage<TConstLawIntegratorType>::CalculateTangen
// Calculates the Tangent Constitutive Tensor by perturbation (second order)
TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, this, ConstitutiveLaw::StressMeasure_Cauchy, consider_perturbation_threshold, 2);
} else if (tangent_operator_estimation == TangentOperatorEstimation::Secant) {
rValues.GetConstitutiveMatrix() *= (1.0 - mDamage);
return;
} else if (tangent_operator_estimation == TangentOperatorEstimation::InitialStiffness) {
return;
}
Expand Down
23 changes: 14 additions & 9 deletions ...ustom_constitutive/small_strains/plasticity/generic_small_strain_isotropic_plasticity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,8 +114,8 @@ void GenericSmallStrainIsotropicPlasticity<TConstLawIntegratorType>::CalculateMa
this->template AddInitialStrainVectorContribution<Vector>(r_strain_vector);

// We compute the stress or the constitutive matrix
if (r_constitutive_law_options.Is( ConstitutiveLaw::COMPUTE_STRESS) ||
r_constitutive_law_options.Is( ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
if (r_constitutive_law_options.Is(ConstitutiveLaw::COMPUTE_STRESS) ||
r_constitutive_law_options.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {

// We get some variables
double threshold = this->GetThreshold();
Expand Down Expand Up @@ -167,9 +167,7 @@ void GenericSmallStrainIsotropicPlasticity<TConstLawIntegratorType>::CalculateMa
noalias(r_integrated_stress_vector) = predictive_stress_vector;

if (r_constitutive_law_options.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
this->CalculateTangentTensor(rValues); // this modifies the ConstitutiveMatrix
} else {
BaseType::CalculateElasticMatrix( r_constitutive_matrix, rValues);
this->CalculateTangentTensor(rValues, plastic_strain); // this modifies the ConstitutiveMatrix
}
}
}
Expand All @@ -180,7 +178,9 @@ void GenericSmallStrainIsotropicPlasticity<TConstLawIntegratorType>::CalculateMa
/***********************************************************************************/

template <class TConstLawIntegratorType>
void GenericSmallStrainIsotropicPlasticity<TConstLawIntegratorType>::CalculateTangentTensor(ConstitutiveLaw::Parameters& rValues)
void GenericSmallStrainIsotropicPlasticity<TConstLawIntegratorType>::CalculateTangentTensor(
ConstitutiveLaw::Parameters& rValues,
const Vector& rPlasticStrain)
{
const Properties& r_material_properties = rValues.GetMaterialProperties();

Expand All @@ -195,11 +195,17 @@ void GenericSmallStrainIsotropicPlasticity<TConstLawIntegratorType>::CalculateTa
} else if (tangent_operator_estimation == TangentOperatorEstimation::SecondOrderPerturbation) {
// Calculates the Tangent Constitutive Tensor by perturbation (second order)
TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, this, ConstitutiveLaw::StressMeasure_Cauchy, consider_perturbation_threshold, 2);
} else if (tangent_operator_estimation == TangentOperatorEstimation::Secant) {
const Vector num = prod(rValues.GetConstitutiveMatrix(), rPlasticStrain);
const double denom = inner_prod(rValues.GetStrainVector(), num);
noalias(rValues.GetConstitutiveMatrix()) -= outer_prod(num, num) / denom;
} else if (tangent_operator_estimation == TangentOperatorEstimation::SecondOrderPerturbationV2) {
// Calculates the Tangent Constitutive Tensor by perturbation (second order)
TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, this, ConstitutiveLaw::StressMeasure_Cauchy, consider_perturbation_threshold, 4);
} else if (tangent_operator_estimation == TangentOperatorEstimation::InitialStiffness) {
BaseType::CalculateElasticMatrix( rValues.GetConstitutiveMatrix(), rValues);
} else if (tangent_operator_estimation == TangentOperatorEstimation::InitialStiffness) {
BaseType::CalculateElasticMatrix(rValues.GetConstitutiveMatrix(), rValues);
} else if (tangent_operator_estimation == TangentOperatorEstimation::OrthogonalSecant) {
TangentOperatorCalculatorUtility::CalculateOrthogonalSecantTensor(rValues);
}
}

Expand Down Expand Up @@ -358,7 +364,6 @@ void GenericSmallStrainIsotropicPlasticity<TConstLawIntegratorType>::FinalizeMat
plastic_dissipation, plastic_strain_increment,
r_constitutive_matrix, plastic_strain, rValues,
characteristic_length);
BaseType::CalculateElasticMatrix(r_constitutive_matrix, rValues);
}

mPlasticDissipation = plastic_dissipation;
Expand Down
4 changes: 3 additions & 1 deletion .../custom_constitutive/small_strains/plasticity/generic_small_strain_isotropic_plasticity.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -441,7 +441,9 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) GenericSmallStrainIsotropicPlast
* @brief This method computes the tangent tensor
* @param rValues The constitutive law parameters and flags
*/
void CalculateTangentTensor(ConstitutiveLaw::Parameters &rValues);
void CalculateTangentTensor(
ConstitutiveLaw::Parameters &rValues,
const Vector& rPlasticStrain);

///@}
///@name Private Access
Expand Down
13 changes: 11 additions & 2 deletions ...ustom_constitutive/small_strains/plasticity/generic_small_strain_kinematic_plasticity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,7 +179,7 @@ void GenericSmallStrainKinematicPlasticity<TConstLawIntegratorType>::CalculateMa
noalias(r_integrated_stress_vector) = predictive_stress_vector;

if (r_constitutive_law_options.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
this->CalculateTangentTensor(rValues); // this modifies the ConstitutiveMatrix
this->CalculateTangentTensor(rValues, plastic_strain); // this modifies the ConstitutiveMatrix
}
}
}
Expand All @@ -191,7 +191,8 @@ void GenericSmallStrainKinematicPlasticity<TConstLawIntegratorType>::CalculateMa

template <class TConstLawIntegratorType>
void GenericSmallStrainKinematicPlasticity<TConstLawIntegratorType>::CalculateTangentTensor(
ConstitutiveLaw::Parameters& rValues
ConstitutiveLaw::Parameters& rValues,
const Vector& rPlasticStrain
)
{
const Properties& r_material_properties = rValues.GetMaterialProperties();
Expand All @@ -207,9 +208,17 @@ void GenericSmallStrainKinematicPlasticity<TConstLawIntegratorType>::CalculateTa
} else if (tangent_operator_estimation == TangentOperatorEstimation::SecondOrderPerturbation) {
// Calculates the Tangent Constitutive Tensor by perturbation (second order)
TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, this, ConstitutiveLaw::StressMeasure_Cauchy, consider_perturbation_threshold, 2);
} else if (tangent_operator_estimation == TangentOperatorEstimation::Secant) {
const Vector num = prod(rValues.GetConstitutiveMatrix(), rPlasticStrain);
const double denom = inner_prod(rValues.GetStrainVector(), num);
noalias(rValues.GetConstitutiveMatrix()) -= outer_prod(num, num) / denom;
} else if (tangent_operator_estimation == TangentOperatorEstimation::SecondOrderPerturbationV2) {
// Calculates the Tangent Constitutive Tensor by perturbation (second order)
TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, this, ConstitutiveLaw::StressMeasure_Cauchy, consider_perturbation_threshold, 4);
} else if (tangent_operator_estimation == TangentOperatorEstimation::InitialStiffness) {
BaseType::CalculateElasticMatrix(rValues.GetConstitutiveMatrix(), rValues);
} else if (tangent_operator_estimation == TangentOperatorEstimation::OrthogonalSecant) {
TangentOperatorCalculatorUtility::CalculateOrthogonalSecantTensor(rValues);
}
}

Expand Down
4 changes: 3 additions & 1 deletion .../custom_constitutive/small_strains/plasticity/generic_small_strain_kinematic_plasticity.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -430,7 +430,9 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) GenericSmallStrainKinematicPlast
* @brief This method computes the tangent tensor
* @param rValues The constitutive law parameters and flags
*/
void CalculateTangentTensor(ConstitutiveLaw::Parameters &rValues);
void CalculateTangentTensor(
ConstitutiveLaw::Parameters &rValues,
const Vector& rPlasticStrain);

///@}
///@name Private Access
Expand Down
7 changes: 3 additions & 4 deletions ...ions/ConstitutiveLawsApplication/custom_utilities/advanced_constitutive_law_utilities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,10 +137,9 @@ void AdvancedConstitutiveLawUtilities<TVoigtSize>::CalculateThirdVector(
const double J2thirds = J2 / 3.0;

if constexpr (Dimension == 2) {
rThirdVector[0] = rDeviator[1] * rDeviator[2] + J2thirds;
rThirdVector[1] = rDeviator[0] * rDeviator[2] + J2thirds;
rThirdVector[2] = rDeviator[0] * rDeviator[1] - std::pow(rDeviator[3], 2) + J2thirds;
rThirdVector[3] = -2.0 * rDeviator[3] * rDeviator[2];
rThirdVector[0] = J2thirds;
rThirdVector[1] = J2thirds;
rThirdVector[2] = 0.0; // The szz should be added when 4-size is ready in plane strain
} else {
rThirdVector[0] = rDeviator[1] * rDeviator[2] - rDeviator[4] * rDeviator[4] + J2thirds;
rThirdVector[1] = rDeviator[0] * rDeviator[2] - rDeviator[5] * rDeviator[5] + J2thirds;
Expand Down
Loading

0 comments on commit a4cec35

Please sign in to comment.