From e4778442273df30ca15b66310d2792e59a1755e6 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Fri, 4 Mar 2022 22:48:54 +0100 Subject: [PATCH 1/9] computing C also when element not provides Strain --- .../hyper_elastic_isotropic_neo_hookean_3d.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp index 889f85fe8cb8..d9e4b04242a3 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp @@ -98,17 +98,19 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponsePK2(Constituti const double lame_lambda = (young_modulus * poisson_coefficient)/((1.0 + poisson_coefficient)*(1.0 - 2.0 * poisson_coefficient)); const double lame_mu = young_modulus/(2.0 * (1.0 + poisson_coefficient)); - // We compute the right Cauchy-Green tensor (C): - const Matrix C_tensor = prod(trans( deformation_gradient_f), deformation_gradient_f); - - // Inverse of the right Cauchy-Green tensor (C): double aux_det; - Matrix inverse_C_tensor(dimension, dimension); - MathUtils::InvertMatrix( C_tensor, inverse_C_tensor, aux_det); + Matrix C_tensor(dimension, dimension), inverse_C_tensor(dimension, dimension); if(r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN )) { this->CalculateGreenLagrangianStrain(rValues, strain_vector); + Matrix strain_tensor(dimension, dimension); + noalias(strain_tensor) = MathUtils StrainVectorToTensor(strain_vector); + noalias(C_tensor) = 2.0 * strain_tensor + IdentityMatrix(dimension); + } else { + noalias(C_tensor) = prod(trans(deformation_gradient_f), deformation_gradient_f); } + double aux_det; + MathUtils::InvertMatrix( C_tensor, inverse_C_tensor, aux_det); if( r_flags.Is( ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR ) ){ Matrix& constitutive_matrix = rValues.GetConstitutiveMatrix(); From b807f354878bd21d08cc1601105e50a3b77421c5 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Fri, 4 Mar 2022 22:52:38 +0100 Subject: [PATCH 2/9] minor --- .../hyper_elastic_isotropic_neo_hookean_3d.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp index d9e4b04242a3..4a65c1639b0d 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp @@ -98,13 +98,12 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponsePK2(Constituti const double lame_lambda = (young_modulus * poisson_coefficient)/((1.0 + poisson_coefficient)*(1.0 - 2.0 * poisson_coefficient)); const double lame_mu = young_modulus/(2.0 * (1.0 + poisson_coefficient)); - double aux_det; Matrix C_tensor(dimension, dimension), inverse_C_tensor(dimension, dimension); if(r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN )) { this->CalculateGreenLagrangianStrain(rValues, strain_vector); Matrix strain_tensor(dimension, dimension); - noalias(strain_tensor) = MathUtils StrainVectorToTensor(strain_vector); + noalias(strain_tensor) = MathUtils::StrainVectorToTensor(strain_vector); noalias(C_tensor) = 2.0 * strain_tensor + IdentityMatrix(dimension); } else { noalias(C_tensor) = prod(trans(deformation_gradient_f), deformation_gradient_f); From 2f44b26b5400d551c92289c9fa531ca45f0d7e1e Mon Sep 17 00:00:00 2001 From: Alejandro Date: Fri, 4 Mar 2022 23:49:14 +0100 Subject: [PATCH 3/9] minor optimization and adaption --- ...hyper_elastic_isotropic_neo_hookean_3d.cpp | 26 ++++++++++--------- 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp index 4a65c1639b0d..ebf325cb567e 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp @@ -91,7 +91,7 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponsePK2(Constituti // The deformation gradient const Matrix& deformation_gradient_f = rValues.GetDeformationGradientF(); - const double determinant_f = rValues.GetDeterminantF(); + double determinant_f = rValues.GetDeterminantF(); KRATOS_ERROR_IF(determinant_f < 0.0) << "Deformation gradient determinant (detF) < 0.0 : " << determinant_f << std::endl; // The LAME parameters @@ -102,14 +102,16 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponsePK2(Constituti if(r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN )) { this->CalculateGreenLagrangianStrain(rValues, strain_vector); + noalias(C_tensor) = prod(trans(deformation_gradient_f), deformation_gradient_f); + } else { Matrix strain_tensor(dimension, dimension); noalias(strain_tensor) = MathUtils::StrainVectorToTensor(strain_vector); noalias(C_tensor) = 2.0 * strain_tensor + IdentityMatrix(dimension); - } else { - noalias(C_tensor) = prod(trans(deformation_gradient_f), deformation_gradient_f); + determinant_f = std::sqrt(MathUtils::Det(C_tensor)); } + double aux_det; - MathUtils::InvertMatrix( C_tensor, inverse_C_tensor, aux_det); + MathUtils::InvertMatrix(C_tensor, inverse_C_tensor, aux_det); if( r_flags.Is( ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR ) ){ Matrix& constitutive_matrix = rValues.GetConstitutiveMatrix(); @@ -524,11 +526,10 @@ void HyperElasticIsotropicNeoHookean3D::CalculatePK2Stress( const double LameMu ) { - Matrix stress_matrix; - const SizeType dimension = WorkingSpaceDimension(); + Matrix stress_matrix(dimension, dimension); - stress_matrix = LameLambda * std::log(DeterminantF) * rInvCTensor + LameMu * ( IdentityMatrix(dimension) - rInvCTensor ); + noalias(stress_matrix) = LameLambda * std::log(DeterminantF) * rInvCTensor + LameMu * ( IdentityMatrix(dimension) - rInvCTensor ); rStressVector = MathUtils::StressTensorToVector( stress_matrix, GetStrainSize() ); } @@ -544,11 +545,10 @@ void HyperElasticIsotropicNeoHookean3D::CalculateKirchhoffStress( const double LameMu ) { - Matrix stress_matrix; - const SizeType dimension = WorkingSpaceDimension(); + Matrix stress_matrix(dimension, dimension); - stress_matrix = LameLambda * std::log(DeterminantF) * IdentityMatrix(dimension) + LameMu * ( rBTensor - IdentityMatrix(dimension) ); + noalias(stress_matrix) = LameLambda * std::log(DeterminantF) * IdentityMatrix(dimension) + LameMu * ( rBTensor - IdentityMatrix(dimension) ); rStressVector = MathUtils::StressTensorToVector( stress_matrix, rStressVector.size() ); } @@ -565,7 +565,8 @@ void HyperElasticIsotropicNeoHookean3D::CalculateGreenLagrangianStrain( const Matrix& F = rValues.GetDeformationGradientF(); // 2.-Compute e = 0.5*(inv(C) - I) - const Matrix C_tensor = prod(trans(F),F); + Matrix C_tensor(Dimension, Dimension); + noalias(C_tensor) = prod(trans(F),F); ConstitutiveLawUtilities::CalculateGreenLagrangianStrain(C_tensor, rStrainVector); } @@ -581,7 +582,8 @@ void HyperElasticIsotropicNeoHookean3D::CalculateAlmansiStrain( const Matrix& F = rValues.GetDeformationGradientF(); // 2.-COmpute e = 0.5*(1-inv(B)) - const Matrix B_tensor = prod(F,trans(F)); + Matrix B_tensor(Dimension, Dimension); + noalias(B_tensor) = prod(F,trans(F)); AdvancedConstitutiveLawUtilities::CalculateAlmansiStrain(B_tensor, rStrainVector); } From 7709a0a7148e1eeb646da9c61e720c1d78e43a67 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Sat, 5 Mar 2022 00:17:07 +0100 Subject: [PATCH 4/9] more opt --- ...hyper_elastic_isotropic_neo_hookean_3d.cpp | 20 ++++++++----------- 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp index ebf325cb567e..8fbc3bebbbcc 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp @@ -526,12 +526,10 @@ void HyperElasticIsotropicNeoHookean3D::CalculatePK2Stress( const double LameMu ) { - const SizeType dimension = WorkingSpaceDimension(); - Matrix stress_matrix(dimension, dimension); - - noalias(stress_matrix) = LameLambda * std::log(DeterminantF) * rInvCTensor + LameMu * ( IdentityMatrix(dimension) - rInvCTensor ); - - rStressVector = MathUtils::StressTensorToVector( stress_matrix, GetStrainSize() ); + Matrix stress_matrix(Dimension, Dimension); + const Matrix Id = IdentityMatrix(Dimension); + noalias(stress_matrix) = LameLambda * std::log(DeterminantF) * rInvCTensor + LameMu * ( Id - rInvCTensor ); + noalias(rStressVector) = MathUtils::StressTensorToVector( stress_matrix, GetStrainSize() ); } /***********************************************************************************/ @@ -545,12 +543,10 @@ void HyperElasticIsotropicNeoHookean3D::CalculateKirchhoffStress( const double LameMu ) { - const SizeType dimension = WorkingSpaceDimension(); - Matrix stress_matrix(dimension, dimension); - - noalias(stress_matrix) = LameLambda * std::log(DeterminantF) * IdentityMatrix(dimension) + LameMu * ( rBTensor - IdentityMatrix(dimension) ); - - rStressVector = MathUtils::StressTensorToVector( stress_matrix, rStressVector.size() ); + Matrix stress_matrix(Dimension, Dimension); + const Matrix Id = IdentityMatrix(Dimension); + noalias(stress_matrix) = LameLambda * std::log(DeterminantF) * Id + LameMu * (rBTensor - Id); + noalias(rStressVector) = MathUtils::StressTensorToVector(stress_matrix, rStressVector.size()); } /***********************************************************************************/ From 121697911dd2a1549ecb71658e5658bebb74cde6 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Sat, 5 Mar 2022 22:19:31 +0100 Subject: [PATCH 5/9] adapting kirchoff method to strain provided --- .../hyper_elastic_isotropic_neo_hookean_3d.cpp | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp index 8fbc3bebbbcc..71a520736b3d 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp @@ -136,7 +136,6 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponseKirchhoff (Cons const Properties& material_properties = rValues.GetMaterialProperties(); Vector& strain_vector = rValues.GetStrainVector(); - Vector& stress_vector = rValues.GetStressVector(); // The material properties const double young_modulus = material_properties[YOUNG_MODULUS]; @@ -144,15 +143,25 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponseKirchhoff (Cons // The deformation gradient const Matrix& deformation_gradient_f = rValues.GetDeformationGradientF(); - const double determinant_f = rValues.GetDeterminantF(); + double determinant_f = rValues.GetDeterminantF(); KRATOS_ERROR_IF(determinant_f < 0.0) << "Deformation gradient determinant (detF) < 0.0 : " << determinant_f << std::endl; // The LAME parameters const double lame_lambda = (young_modulus * poisson_coefficient)/((1.0 + poisson_coefficient)*(1.0 - 2.0 * poisson_coefficient)); const double lame_mu = young_modulus/(2.0 * (1.0 + poisson_coefficient)); + Matrix B_tensor(Dimension, Dimension), inverse_B_tensor(Dimension, Dimension); + if(r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN )) { CalculateAlmansiStrain(rValues, strain_vector); + noalias(B_tensor) = prod(deformation_gradient_f, trans( deformation_gradient_f)); + } else { + Matrix strain_tensor(Dimension, Dimension); + noalias(strain_tensor) = MathUtils::StrainVectorToTensor(strain_vector); + noalias(inverse_B_tensor) = IdentityMatrix(Dimension) - 2.0 * strain_tensor; + double aux_det; + MathUtils::InvertMatrix(inverse_B_tensor, B_tensor, aux_det); + determinant_f = std::sqrt(MathUtils::Det(B_tensor)); } if( r_flags.Is( ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR ) ) { @@ -162,7 +171,7 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponseKirchhoff (Cons if( r_flags.Is( ConstitutiveLaw::COMPUTE_STRESS ) ) { // We compute the left Cauchy-Green tensor (B): - const Matrix B_tensor = prod(deformation_gradient_f, trans( deformation_gradient_f)); + Vector& stress_vector = rValues.GetStressVector(); CalculateKirchhoffStress( B_tensor, stress_vector, determinant_f, lame_lambda, lame_mu ); } } From 1e609a680e74bb099a74a3b226140b33a9eaf9bc Mon Sep 17 00:00:00 2001 From: Alejandro Date: Sat, 5 Mar 2022 22:24:04 +0100 Subject: [PATCH 6/9] minor opt --- .../hyper_elastic_isotropic_neo_hookean_3d.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp index 71a520736b3d..e580c3dc3b87 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp @@ -150,7 +150,7 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponseKirchhoff (Cons const double lame_lambda = (young_modulus * poisson_coefficient)/((1.0 + poisson_coefficient)*(1.0 - 2.0 * poisson_coefficient)); const double lame_mu = young_modulus/(2.0 * (1.0 + poisson_coefficient)); - Matrix B_tensor(Dimension, Dimension), inverse_B_tensor(Dimension, Dimension); + Matrix B_tensor(Dimension, Dimension); if(r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN )) { CalculateAlmansiStrain(rValues, strain_vector); @@ -158,6 +158,7 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponseKirchhoff (Cons } else { Matrix strain_tensor(Dimension, Dimension); noalias(strain_tensor) = MathUtils::StrainVectorToTensor(strain_vector); + Matrix inverse_B_tensor(Dimension, Dimension); noalias(inverse_B_tensor) = IdentityMatrix(Dimension) - 2.0 * strain_tensor; double aux_det; MathUtils::InvertMatrix(inverse_B_tensor, B_tensor, aux_det); From 82cb86721f9b7495cf3d9c68bee4cd14cb0bbb9f Mon Sep 17 00:00:00 2001 From: Alejandro Date: Sun, 6 Mar 2022 10:16:21 +0100 Subject: [PATCH 7/9] using the constexpre --- .../hyper_elastic_isotropic_neo_hookean_3d.cpp | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp index e580c3dc3b87..67f3b9492b9a 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp @@ -80,8 +80,6 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponsePK2(Constituti // Get Values to compute the constitutive law: Flags& r_flags=rValues.GetOptions(); - const SizeType dimension = WorkingSpaceDimension(); - const Properties& material_properties = rValues.GetMaterialProperties(); Vector& strain_vector = rValues.GetStrainVector(); @@ -98,15 +96,15 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponsePK2(Constituti const double lame_lambda = (young_modulus * poisson_coefficient)/((1.0 + poisson_coefficient)*(1.0 - 2.0 * poisson_coefficient)); const double lame_mu = young_modulus/(2.0 * (1.0 + poisson_coefficient)); - Matrix C_tensor(dimension, dimension), inverse_C_tensor(dimension, dimension); + Matrix C_tensor(Dimension, Dimension), inverse_C_tensor(Dimension, Dimension); if(r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN )) { this->CalculateGreenLagrangianStrain(rValues, strain_vector); noalias(C_tensor) = prod(trans(deformation_gradient_f), deformation_gradient_f); } else { - Matrix strain_tensor(dimension, dimension); + Matrix strain_tensor(Dimension, Dimension); noalias(strain_tensor) = MathUtils::StrainVectorToTensor(strain_vector); - noalias(C_tensor) = 2.0 * strain_tensor + IdentityMatrix(dimension); + noalias(C_tensor) = 2.0 * strain_tensor + IdentityMatrix(Dimension); determinant_f = std::sqrt(MathUtils::Det(C_tensor)); } @@ -441,7 +439,7 @@ void HyperElasticIsotropicNeoHookean3D::GetLawFeatures(Features& rFeatures) rFeatures.mStrainMeasures.push_back(StrainMeasure_Deformation_Gradient); //Set the strain size - rFeatures.mStrainSize = GetStrainSize(); + rFeatures.mStrainSize = VoigtSize; //Set the spacedimension rFeatures.mSpaceDimension = WorkingSpaceDimension(); @@ -539,7 +537,7 @@ void HyperElasticIsotropicNeoHookean3D::CalculatePK2Stress( Matrix stress_matrix(Dimension, Dimension); const Matrix Id = IdentityMatrix(Dimension); noalias(stress_matrix) = LameLambda * std::log(DeterminantF) * rInvCTensor + LameMu * ( Id - rInvCTensor ); - noalias(rStressVector) = MathUtils::StressTensorToVector( stress_matrix, GetStrainSize() ); + noalias(rStressVector) = MathUtils::StressTensorToVector( stress_matrix, VoigtSize ); } /***********************************************************************************/ From 6ec81dd9dc9f031dbe6fdc347df62b9c213734a0 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Sun, 6 Mar 2022 11:02:21 +0100 Subject: [PATCH 8/9] opt plane strain neo hooke --- ...yper_elastic_isotropic_neo_hookean_plane_strain_2d.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_plane_strain_2d.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_plane_strain_2d.cpp index 27771ec369cc..e87b48ff98b3 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_plane_strain_2d.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_plane_strain_2d.cpp @@ -141,7 +141,8 @@ void HyperElasticIsotropicNeoHookeanPlaneStrain2D::CalculateGreenLagrangianStrai const Matrix& F = rValues.GetDeformationGradientF(); // e = 0.5*(inv(C) - I) - Matrix C_tensor = prod(trans(F),F); + Matrix C_tensor(Dimension, Dimension); + noalias(C_tensor) = prod(trans(F), F); rStrainVector[0] = 0.5 * ( C_tensor( 0, 0 ) - 1.00 ); rStrainVector[1] = 0.5 * ( C_tensor( 1, 1 ) - 1.00 ); @@ -160,10 +161,11 @@ void HyperElasticIsotropicNeoHookeanPlaneStrain2D::CalculateAlmansiStrain( const Matrix& F = rValues.GetDeformationGradientF(); // e = 0.5*(1-inv(B)) - Matrix B_tensor = prod(F,trans(F)); + Matrix B_tensor(Dimension, Dimension); + noalias(B_tensor) = prod(F,trans(F)); //Calculating the inverse of the jacobian - Matrix inverse_B_tensor ( 2, 2 ); + Matrix inverse_B_tensor ( Dimension, Dimension ); double aux_det_b = 0; MathUtils::InvertMatrix( B_tensor, inverse_B_tensor, aux_det_b); From f7156cd23b0fde39d19e2e4b2bd0243c08d2501d Mon Sep 17 00:00:00 2001 From: Alejandro Date: Mon, 7 Mar 2022 09:58:53 +0100 Subject: [PATCH 9/9] now corrected --- ...hyper_elastic_isotropic_neo_hookean_3d.cpp | 40 ++++++++++++------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp index 67f3b9492b9a..d5f9de72ff59 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/hyper_elastic_isotropic_neo_hookean_3d.cpp @@ -80,6 +80,8 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponsePK2(Constituti // Get Values to compute the constitutive law: Flags& r_flags=rValues.GetOptions(); + const SizeType dimension = WorkingSpaceDimension(); + const Properties& material_properties = rValues.GetMaterialProperties(); Vector& strain_vector = rValues.GetStrainVector(); @@ -96,15 +98,15 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponsePK2(Constituti const double lame_lambda = (young_modulus * poisson_coefficient)/((1.0 + poisson_coefficient)*(1.0 - 2.0 * poisson_coefficient)); const double lame_mu = young_modulus/(2.0 * (1.0 + poisson_coefficient)); - Matrix C_tensor(Dimension, Dimension), inverse_C_tensor(Dimension, Dimension); + Matrix C_tensor(dimension, dimension), inverse_C_tensor(dimension, dimension); if(r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN )) { this->CalculateGreenLagrangianStrain(rValues, strain_vector); noalias(C_tensor) = prod(trans(deformation_gradient_f), deformation_gradient_f); } else { - Matrix strain_tensor(Dimension, Dimension); + Matrix strain_tensor(dimension, dimension); noalias(strain_tensor) = MathUtils::StrainVectorToTensor(strain_vector); - noalias(C_tensor) = 2.0 * strain_tensor + IdentityMatrix(Dimension); + noalias(C_tensor) = 2.0 * strain_tensor + IdentityMatrix(dimension); determinant_f = std::sqrt(MathUtils::Det(C_tensor)); } @@ -132,6 +134,8 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponseKirchhoff (Cons // Get Values to compute the constitutive law: Flags& r_flags=rValues.GetOptions(); + const SizeType dimension = WorkingSpaceDimension(); + const Properties& material_properties = rValues.GetMaterialProperties(); Vector& strain_vector = rValues.GetStrainVector(); @@ -148,16 +152,16 @@ void HyperElasticIsotropicNeoHookean3D::CalculateMaterialResponseKirchhoff (Cons const double lame_lambda = (young_modulus * poisson_coefficient)/((1.0 + poisson_coefficient)*(1.0 - 2.0 * poisson_coefficient)); const double lame_mu = young_modulus/(2.0 * (1.0 + poisson_coefficient)); - Matrix B_tensor(Dimension, Dimension); + Matrix B_tensor(dimension, dimension); if(r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN )) { CalculateAlmansiStrain(rValues, strain_vector); noalias(B_tensor) = prod(deformation_gradient_f, trans( deformation_gradient_f)); } else { - Matrix strain_tensor(Dimension, Dimension); + Matrix strain_tensor(dimension, dimension); noalias(strain_tensor) = MathUtils::StrainVectorToTensor(strain_vector); - Matrix inverse_B_tensor(Dimension, Dimension); - noalias(inverse_B_tensor) = IdentityMatrix(Dimension) - 2.0 * strain_tensor; + Matrix inverse_B_tensor(dimension, dimension); + noalias(inverse_B_tensor) = IdentityMatrix(dimension) - 2.0 * strain_tensor; double aux_det; MathUtils::InvertMatrix(inverse_B_tensor, B_tensor, aux_det); determinant_f = std::sqrt(MathUtils::Det(B_tensor)); @@ -534,10 +538,11 @@ void HyperElasticIsotropicNeoHookean3D::CalculatePK2Stress( const double LameMu ) { - Matrix stress_matrix(Dimension, Dimension); - const Matrix Id = IdentityMatrix(Dimension); + const SizeType dimension = WorkingSpaceDimension(); + Matrix stress_matrix(dimension, dimension); + const Matrix Id = IdentityMatrix(dimension); noalias(stress_matrix) = LameLambda * std::log(DeterminantF) * rInvCTensor + LameMu * ( Id - rInvCTensor ); - noalias(rStressVector) = MathUtils::StressTensorToVector( stress_matrix, VoigtSize ); + noalias(rStressVector) = MathUtils::StressTensorToVector( stress_matrix, GetStrainSize()); } /***********************************************************************************/ @@ -551,10 +556,11 @@ void HyperElasticIsotropicNeoHookean3D::CalculateKirchhoffStress( const double LameMu ) { - Matrix stress_matrix(Dimension, Dimension); - const Matrix Id = IdentityMatrix(Dimension); + const SizeType dimension = WorkingSpaceDimension(); + Matrix stress_matrix(dimension, dimension); + const Matrix Id = IdentityMatrix(dimension); noalias(stress_matrix) = LameLambda * std::log(DeterminantF) * Id + LameMu * (rBTensor - Id); - noalias(rStressVector) = MathUtils::StressTensorToVector(stress_matrix, rStressVector.size()); + noalias(rStressVector) = MathUtils::StressTensorToVector(stress_matrix, GetStrainSize()); } /***********************************************************************************/ @@ -565,11 +571,13 @@ void HyperElasticIsotropicNeoHookean3D::CalculateGreenLagrangianStrain( Vector& rStrainVector ) { + const SizeType dimension = WorkingSpaceDimension(); + // 1.-Compute total deformation gradient const Matrix& F = rValues.GetDeformationGradientF(); // 2.-Compute e = 0.5*(inv(C) - I) - Matrix C_tensor(Dimension, Dimension); + Matrix C_tensor(dimension, dimension); noalias(C_tensor) = prod(trans(F),F); ConstitutiveLawUtilities::CalculateGreenLagrangianStrain(C_tensor, rStrainVector); } @@ -582,11 +590,13 @@ void HyperElasticIsotropicNeoHookean3D::CalculateAlmansiStrain( Vector& rStrainVector ) { + const SizeType dimension = WorkingSpaceDimension(); + // 1.-Compute total deformation gradient const Matrix& F = rValues.GetDeformationGradientF(); // 2.-COmpute e = 0.5*(1-inv(B)) - Matrix B_tensor(Dimension, Dimension); + Matrix B_tensor(dimension, dimension); noalias(B_tensor) = prod(F,trans(F)); AdvancedConstitutiveLawUtilities::CalculateAlmansiStrain(B_tensor, rStrainVector); }