From 079392e4f57916afd0720cc082e983ddc29ca03a Mon Sep 17 00:00:00 2001 From: Alejandro Date: Sat, 4 Jan 2025 15:23:24 +0100 Subject: [PATCH 01/16] update DOI and refs --- .../serial_parallel_rule_of_mixtures_law.cpp | 6 +++--- .../composites/serial_parallel_rule_of_mixtures_law.h | 11 ++++++----- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index c1e38913e84e..5391182c84ff 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -9,9 +9,9 @@ // license: structural_mechanics_application/license.txt // // Main authors: Alejandro Cornejo -// Vicente Mataix Ferrandiz -// Fernando Rastellini -// Collaborator: Lucia Barbu +// Vicente Mataix +// +// // // System includes diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h index e3711859670e..abf75f0dcfe9 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h @@ -10,8 +10,8 @@ // // Main authors: Alejandro Cornejo // Vicente Mataix -// Fernando Rastellini -// Collaborator: Lucia Barbu +// +// // #pragma once @@ -47,7 +47,8 @@ namespace Kratos /** * @class SerialParallelRuleOfMixturesLaw * @ingroup StructuralMechanicsApplication - * @brief This CL implements the serial-parallel rule of mixtures developed by F.Rastellini + * @brief This CL implements the serial-parallel rule of mixtures detailed in Cornejo et al. "Methodology for the analysis of post-tensioned structures using a constitutive serial-parallel rule of mixtures" + * DOI: https://doi.org/10.1016/j.compstruct.2018.05.123 * @details * @author Alejandro Cornejo */ @@ -59,10 +60,10 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw ///@{ /// The node definition - typedef Node NodeType; + using NodeType = Node; /// The geometry definition - typedef Geometry GeometryType; + using GeometryType = Geometry; /// Definition of the machine precision tolerance static constexpr double machine_tolerance = std::numeric_limits::epsilon(); From 486dc72975e18256c6a2cb84f3fa1bd7735bfaaf Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 08:57:57 +0100 Subject: [PATCH 02/16] improvements --- .../serial_parallel_rule_of_mixtures_law.cpp | 3 +++ .../serial_parallel_rule_of_mixtures_law.h | 20 ++++++++++++++----- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 5391182c84ff..02711c012a57 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -1641,6 +1641,9 @@ int SerialParallelRuleOfMixturesLaw::Check( KRATOS_ERROR << "A wrong fiber volumetric participation has been set: Greater than 1 or lower than 0... " << std::to_string(mFiberVolumetricParticipation) << std::endl; aux_out += 1; } + KRATOS_ERROR_IF(mpMatrixConstitutiveLaw->WorkingSpaceDimension() != mpFiberConstitutiveLaw->WorkingSpaceDimension()) << "The WorkingSpaceDimension of the fiber and matrix mismatch..." << std::endl; + KRATOS_ERROR_IF(mpMatrixConstitutiveLaw->GetStrainSize() != mpFiberConstitutiveLaw->GetStrainSize()) << "The GetStrainSize of the fiber and matrix mismatch..." << std::endl; + return aux_out; } diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h index abf75f0dcfe9..29a7290c4412 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h @@ -49,6 +49,7 @@ namespace Kratos * @ingroup StructuralMechanicsApplication * @brief This CL implements the serial-parallel rule of mixtures detailed in Cornejo et al. "Methodology for the analysis of post-tensioned structures using a constitutive serial-parallel rule of mixtures" * DOI: https://doi.org/10.1016/j.compstruct.2018.05.123 + * The SP-RoM is able to work in 2D and in 3D. * @details * @author Alejandro Cornejo */ @@ -133,7 +134,8 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw */ SizeType WorkingSpaceDimension() override { - return 3; + return mpMatrixConstitutiveLaw->WorkingSpaceDimension(); + KRATOS_DEBUG_ERROR_IF(mpMatrixConstitutiveLaw->WorkingSpaceDimension() != mpFiberConstitutiveLaw->WorkingSpaceDimension()) << "Th WorkingSpaceDimension of the fiber and matrix mismatch..." << std::endl; }; /** @@ -141,7 +143,9 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw */ SizeType GetStrainSize() const override { - return 6; + return mpMatrixConstitutiveLaw->GetStrainSize(); + KRATOS_DEBUG_ERROR_IF(mpMatrixConstitutiveLaw->GetStrainSize() != mpFiberConstitutiveLaw->GetStrainSize()) << "Th GetStrainSize of the fiber and matrix mismatch..." << std::endl; + }; /** @@ -261,7 +265,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw * @param rValue new value of the specified variable * @param rCurrentProcessInfo the process info */ - void SetValue( + void SetValue( const Variable& rThisVariable, const double& rValue, const ProcessInfo& rCurrentProcessInfo @@ -548,7 +552,10 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw */ bool RequiresInitializeMaterialResponse() override { - return true; + if (mpMatrixConstitutiveLaw->RequiresInitializeMaterialResponse() || mpFiberConstitutiveLaw->RequiresInitializeMaterialResponse()) { + return true; + } + return false; } /** @@ -556,7 +563,10 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw */ bool RequiresFinalizeMaterialResponse() override { - return true; + if (mpMatrixConstitutiveLaw->RequiresFinalizeMaterialResponse() || mpFiberConstitutiveLaw->RequiresFinalizeMaterialResponse()) { + return true; + } + return false; } /** From 332cf8514969ecdba90aa90c07b9125b013f9ca9 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 09:33:24 +0100 Subject: [PATCH 03/16] unify the perturbation method in the utils --- .../advanced_constitutive_law_utilities.cpp | 25 +++++++++++++++++++ .../advanced_constitutive_law_utilities.h | 9 +++++++ 2 files changed, 34 insertions(+) diff --git a/applications/ConstitutiveLawsApplication/custom_utilities/advanced_constitutive_law_utilities.cpp b/applications/ConstitutiveLawsApplication/custom_utilities/advanced_constitutive_law_utilities.cpp index 7f13d7309d54..60755c038738 100644 --- a/applications/ConstitutiveLawsApplication/custom_utilities/advanced_constitutive_law_utilities.cpp +++ b/applications/ConstitutiveLawsApplication/custom_utilities/advanced_constitutive_law_utilities.cpp @@ -19,6 +19,7 @@ #include "includes/global_variables.h" #include "utilities/math_utils.h" #include "custom_utilities/advanced_constitutive_law_utilities.h" +#include "custom_utilities/tangent_operator_calculator_utility.h" #include "constitutive_laws_application_variables.h" namespace Kratos @@ -991,6 +992,30 @@ void AdvancedConstitutiveLawUtilities::CalculateOrthotropicElasticMa /***********************************************************************************/ /***********************************************************************************/ +template +void AdvancedConstitutiveLawUtilities::CalculateTangentTensorByPerturbation( + ConstitutiveLaw::Parameters& rValues, + ConstitutiveLaw* pConstitutiveLaw, + const ConstitutiveLaw::StressMeasure& rStressMeasure + ) +{ + const Properties& r_material_properties = rValues.GetMaterialProperties(); + + const bool consider_perturbation_threshold = r_material_properties.Has(CONSIDER_PERTURBATION_THRESHOLD) ? r_material_properties[CONSIDER_PERTURBATION_THRESHOLD] : true; + const TangentOperatorEstimation tangent_operator_estimation = r_material_properties.Has(TANGENT_OPERATOR_ESTIMATION) ? static_cast(r_material_properties[TANGENT_OPERATOR_ESTIMATION]) : TangentOperatorEstimation::SecondOrderPerturbation; + + if (tangent_operator_estimation == TangentOperatorEstimation::FirstOrderPerturbation) { + // Calculates the Tangent Constitutive Tensor by perturbation (first order) + TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, pConstitutiveLaw, rStressMeasure, consider_perturbation_threshold, 1); + } else if (tangent_operator_estimation == TangentOperatorEstimation::SecondOrderPerturbation) { + // Calculates the Tangent Constitutive Tensor by perturbation (second order) + TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, pConstitutiveLaw, rStressMeasure, consider_perturbation_threshold, 2); + } else if (tangent_operator_estimation == TangentOperatorEstimation::SecondOrderPerturbationV2) { + // Calculates the Tangent Constitutive Tensor by perturbation (second order) + TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, pConstitutiveLaw, rStressMeasure, consider_perturbation_threshold, 4); + } +} + template class AdvancedConstitutiveLawUtilities<3>; template class AdvancedConstitutiveLawUtilities<6>; diff --git a/applications/ConstitutiveLawsApplication/custom_utilities/advanced_constitutive_law_utilities.h b/applications/ConstitutiveLawsApplication/custom_utilities/advanced_constitutive_law_utilities.h index 265cbe7587fa..141e687ba20f 100644 --- a/applications/ConstitutiveLawsApplication/custom_utilities/advanced_constitutive_law_utilities.h +++ b/applications/ConstitutiveLawsApplication/custom_utilities/advanced_constitutive_law_utilities.h @@ -533,5 +533,14 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) AdvancedConstitutiveLawUtilities const double Temperature ); + /** + * @brief This method computes the tangent tensor + * @param rValues The constitutive law parameters and flags + */ + static void CalculateTangentTensorByPerturbation( + ConstitutiveLaw::Parameters& rValues, + ConstitutiveLaw* pConstitutiveLaw, + const ConstitutiveLaw::StressMeasure& rStressMeasure = ConstitutiveLaw::StressMeasure_Cauchy); + }; // class AdvancedConstitutiveLawUtilities } // namespace Kratos \ No newline at end of file From 98059c60993a45f060541633e0052131f2e43859 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 09:38:32 +0100 Subject: [PATCH 04/16] more improvements --- .../serial_parallel_rule_of_mixtures_law.cpp | 24 +++++++------------ .../serial_parallel_rule_of_mixtures_law.h | 4 ++-- 2 files changed, 11 insertions(+), 17 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 02711c012a57..2cf82e960af6 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -621,7 +621,13 @@ void SerialParallelRuleOfMixturesLaw::CalculateAlmansiStrain(ConstitutiveLaw::Pa B_tensor.resize(dimension, dimension, false); noalias(B_tensor) = prod(F, trans(F)); - AdvancedConstitutiveLawUtilities<6>::CalculateAlmansiStrain(B_tensor, r_strain_vector); + const SizeType strain_size = GetStrainSize(); + if (strain_size == 3) { + AdvancedConstitutiveLawUtilities<3>::CalculateAlmansiStrain(B_tensor, r_strain_vector); + } else { + AdvancedConstitutiveLawUtilities<6>::CalculateAlmansiStrain(B_tensor, r_strain_vector); + } + } /***********************************************************************************/ @@ -1606,20 +1612,8 @@ void SerialParallelRuleOfMixturesLaw::CalculateTangentTensor( ConstitutiveLaw::Parameters& rValues, const ConstitutiveLaw::StressMeasure& rStressMeasure) { - const Properties& r_material_properties = rValues.GetMaterialProperties(); - - const bool consider_perturbation_threshold = r_material_properties.Has(CONSIDER_PERTURBATION_THRESHOLD) ? r_material_properties[CONSIDER_PERTURBATION_THRESHOLD] : true; - const TangentOperatorEstimation tangent_operator_estimation = r_material_properties.Has(TANGENT_OPERATOR_ESTIMATION) ? static_cast(r_material_properties[TANGENT_OPERATOR_ESTIMATION]) : TangentOperatorEstimation::SecondOrderPerturbation; - - if (tangent_operator_estimation == TangentOperatorEstimation::Analytic) { - KRATOS_ERROR << "Analytic solution not available" << std::endl; - } else if (tangent_operator_estimation == TangentOperatorEstimation::FirstOrderPerturbation) { - // Calculates the Tangent Constitutive Tensor by perturbation (first order) - TangentOperatorCalculatorUtility::CalculateTangentTensor(rValues, this, rStressMeasure, consider_perturbation_threshold, 1); - } 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); - } + // Independent from the Voigt Size + AdvancedConstitutiveLawUtilities<6>::CalculateTangentTensorByPerturbation(rValues, this, rStressMeasure); } /***********************************************************************************/ /***********************************************************************************/ diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h index 29a7290c4412..076d52d424ad 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h @@ -135,7 +135,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw SizeType WorkingSpaceDimension() override { return mpMatrixConstitutiveLaw->WorkingSpaceDimension(); - KRATOS_DEBUG_ERROR_IF(mpMatrixConstitutiveLaw->WorkingSpaceDimension() != mpFiberConstitutiveLaw->WorkingSpaceDimension()) << "Th WorkingSpaceDimension of the fiber and matrix mismatch..." << std::endl; + KRATOS_DEBUG_ERROR_IF(mpMatrixConstitutiveLaw->WorkingSpaceDimension() != mpFiberConstitutiveLaw->WorkingSpaceDimension()) << "The WorkingSpaceDimension of the fiber and matrix mismatch..." << std::endl; }; /** @@ -144,7 +144,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw SizeType GetStrainSize() const override { return mpMatrixConstitutiveLaw->GetStrainSize(); - KRATOS_DEBUG_ERROR_IF(mpMatrixConstitutiveLaw->GetStrainSize() != mpFiberConstitutiveLaw->GetStrainSize()) << "Th GetStrainSize of the fiber and matrix mismatch..." << std::endl; + KRATOS_DEBUG_ERROR_IF(mpMatrixConstitutiveLaw->GetStrainSize() != mpFiberConstitutiveLaw->GetStrainSize()) << "The GetStrainSize of the fiber and matrix mismatch..." << std::endl; }; From ced3789af9945a6aa2b2060e3fe38d6f9b5d94cd Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 10:10:35 +0100 Subject: [PATCH 05/16] more... --- .../serial_parallel_rule_of_mixtures_law.cpp | 1 - .../serial_parallel_rule_of_mixtures_law.h | 21 +++++++++++-------- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 2cf82e960af6..6ea2d102f693 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -22,7 +22,6 @@ #include "utilities/math_utils.h" #include "constitutive_laws_application_variables.h" #include "serial_parallel_rule_of_mixtures_law.h" -#include "custom_utilities/tangent_operator_calculator_utility.h" #include "custom_utilities/advanced_constitutive_law_utilities.h" #include "custom_utilities/constitutive_law_utilities.h" diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h index 076d52d424ad..cae73bfb6a49 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h @@ -66,6 +66,8 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw /// The geometry definition using GeometryType = Geometry; + using BaseType = ConstitutiveLaw; + /// Definition of the machine precision tolerance static constexpr double machine_tolerance = std::numeric_limits::epsilon(); @@ -101,7 +103,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw // Copy constructor SerialParallelRuleOfMixturesLaw(SerialParallelRuleOfMixturesLaw const& rOther) - : ConstitutiveLaw(rOther), mpMatrixConstitutiveLaw(rOther.mpMatrixConstitutiveLaw), mpFiberConstitutiveLaw(rOther.mpFiberConstitutiveLaw), + : BaseType(rOther), mpMatrixConstitutiveLaw(rOther.mpMatrixConstitutiveLaw), mpFiberConstitutiveLaw(rOther.mpFiberConstitutiveLaw), mFiberVolumetricParticipation(rOther.mFiberVolumetricParticipation), mParallelDirections(rOther.mParallelDirections) , mPreviousStrainVector(rOther.mPreviousStrainVector) , mPreviousSerialStrainMatrix(rOther.mPreviousSerialStrainMatrix) , mIsPrestressed(rOther.mIsPrestressed) { @@ -134,8 +136,8 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw */ SizeType WorkingSpaceDimension() override { - return mpMatrixConstitutiveLaw->WorkingSpaceDimension(); KRATOS_DEBUG_ERROR_IF(mpMatrixConstitutiveLaw->WorkingSpaceDimension() != mpFiberConstitutiveLaw->WorkingSpaceDimension()) << "The WorkingSpaceDimension of the fiber and matrix mismatch..." << std::endl; + return mpMatrixConstitutiveLaw->WorkingSpaceDimension(); }; /** @@ -143,8 +145,8 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw */ SizeType GetStrainSize() const override { - return mpMatrixConstitutiveLaw->GetStrainSize(); KRATOS_DEBUG_ERROR_IF(mpMatrixConstitutiveLaw->GetStrainSize() != mpFiberConstitutiveLaw->GetStrainSize()) << "The GetStrainSize of the fiber and matrix mismatch..." << std::endl; + return mpMatrixConstitutiveLaw->GetStrainSize(); }; @@ -544,7 +546,8 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw * @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 ConstitutiveLaw::StressMeasure& rStressMeasure = ConstitutiveLaw::StressMeasure_Cauchy); /** @@ -652,7 +655,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw int GetNumberOfSerialComponents() { const int parallel_components = inner_prod(mParallelDirections, mParallelDirections); - return this->GetStrainSize() - parallel_components; + return GetStrainSize() - parallel_components; } ///@} @@ -683,8 +686,8 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw ConstitutiveLaw::Pointer mpMatrixConstitutiveLaw; ConstitutiveLaw::Pointer mpFiberConstitutiveLaw; double mFiberVolumetricParticipation; - array_1d mParallelDirections = ZeroVector(6); - array_1d mPreviousStrainVector = ZeroVector(6); + Vector mParallelDirections = ZeroVector(GetStrainSize()); + Vector mPreviousStrainVector = ZeroVector(GetStrainSize()); Vector mPreviousSerialStrainMatrix = ZeroVector(GetNumberOfSerialComponents()); bool mIsPrestressed = false; @@ -714,7 +717,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw void save(Serializer& rSerializer) const override { - KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, ConstitutiveLaw) + KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, BaseType) rSerializer.save("MatrixConstitutiveLaw", mpMatrixConstitutiveLaw); rSerializer.save("FiberConstitutiveLaw", mpFiberConstitutiveLaw); rSerializer.save("FiberVolumetricParticipation", mFiberVolumetricParticipation); @@ -726,7 +729,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw void load(Serializer& rSerializer) override { - KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, ConstitutiveLaw) + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, BaseType) rSerializer.load("MatrixConstitutiveLaw", mpMatrixConstitutiveLaw); rSerializer.load("FiberConstitutiveLaw", mpFiberConstitutiveLaw); rSerializer.load("FiberVolumetricParticipation", mFiberVolumetricParticipation); From 5321a0ec57554bac57996851b214ca1614e8104e Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 10:40:34 +0100 Subject: [PATCH 06/16] improved more --- .../serial_parallel_rule_of_mixtures_law.cpp | 96 ++++++++++--------- 1 file changed, 53 insertions(+), 43 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 6ea2d102f693..1ec65651a372 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -34,7 +34,7 @@ ConstitutiveLaw::Pointer SerialParallelRuleOfMixturesLaw::Create(Kratos::Paramet if (fiber_volumetric_participation < 0.0 || fiber_volumetric_participation > 1.0) { KRATOS_ERROR << "A wrong fiber volumetric participation has been set: Greater than 1 or lower than 0..." << std::endl; } - const int voigt_size = 6; + const int voigt_size = GetStrainSize(); Vector parallel_directions(voigt_size); for (IndexType i_comp = 0; i_comp < voigt_size; ++i_comp) { parallel_directions[i_comp] = NewParameters["parallel_behaviour_directions"][i_comp].GetInt(); @@ -253,18 +253,18 @@ void SerialParallelRuleOfMixturesLaw::IntegrateStrainSerialParallelBehaviour( const ConstitutiveLaw::StressMeasure& rStressMeasure ) { - const std::size_t voigt_size = this->GetStrainSize(); - const int num_parallel_components = inner_prod(mParallelDirections, mParallelDirections); - const int num_serial_components = voigt_size - num_parallel_components; + const SizeType voigt_size = GetStrainSize(); + const SizeType num_parallel_components = inner_prod(mParallelDirections, mParallelDirections); + const SizeType num_serial_components = voigt_size - num_parallel_components; Matrix parallel_projector(voigt_size, num_parallel_components), serial_projector(num_serial_components, voigt_size); - this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); + CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); bool is_converged = false; int iteration = 0, max_iterations = 150; - Vector parallel_strain_matrix(num_parallel_components), stress_residual(rSerialStrainMatrix.size()); + Vector parallel_strain_matrix(num_parallel_components), stress_residual(num_serial_components); Matrix constitutive_tensor_matrix_ss(num_serial_components, num_serial_components), constitutive_tensor_fiber_ss(num_serial_components, num_serial_components); @@ -286,7 +286,7 @@ void SerialParallelRuleOfMixturesLaw::IntegrateStrainSerialParallelBehaviour( break; } else { // We correct the independent var: serial_strain_matrix - this->CorrectSerialStrainMatrix(rValues, stress_residual, rSerialStrainMatrix, serial_projector, rStressMeasure); + CorrectSerialStrainMatrix(rValues, stress_residual, rSerialStrainMatrix, serial_projector, rStressMeasure); iteration++; } } @@ -303,7 +303,7 @@ void SerialParallelRuleOfMixturesLaw::CorrectSerialStrainMatrix( const ConstitutiveLaw::StressMeasure& rStressMeasure ) { - const std::size_t voigt_size = this->GetStrainSize(); + const SizeType voigt_size = this->GetStrainSize(); const int num_parallel_components = inner_prod(mParallelDirections, mParallelDirections); const int num_serial_components = voigt_size - num_parallel_components; @@ -423,8 +423,16 @@ void SerialParallelRuleOfMixturesLaw::IntegrateStressesOfFiberAndMatrix( const ConstitutiveLaw::StressMeasure& rStressMeasure ) { - rMatrixStressVector.resize(GetStrainSize(), false); - rFiberStressVector.resize(GetStrainSize(), false); + const SizeType voigt_size = GetStrainSize(); + + if (rMatrixStressVector.size() != voigt_size) { + rMatrixStressVector.resize(voigt_size, false); + } + + if (rFiberStressVector.size() != voigt_size) { + rFiberStressVector.resize(voigt_size, false); + } + auto& r_material_properties = rValues.GetMaterialProperties(); const auto it_cl_begin = r_material_properties.GetSubProperties().begin(); const auto& r_props_matrix_cl = *(it_cl_begin); @@ -462,7 +470,7 @@ void SerialParallelRuleOfMixturesLaw::CalculateInitialApproximationSerialStrainM const ConstitutiveLaw::StressMeasure& rStressMeasure ) { - const std::size_t voigt_size = this->GetStrainSize(); + const SizeType voigt_size = GetStrainSize(); const Vector& r_total_strain_vector_parallel = prod(trans(rParallelProjector), rStrainVector); const Vector& r_total_strain_vector_serial = prod(rSerialProjector, rStrainVector); @@ -561,21 +569,21 @@ void SerialParallelRuleOfMixturesLaw::CalculateSerialParallelProjectionMatrices( Matrix& rSerialProjector ) { - const std::size_t voigt_size = this->GetStrainSize(); - const int num_parallel_components = inner_prod(mParallelDirections, mParallelDirections); + const SizeType voigt_size = GetStrainSize(); + const SizeType num_parallel_components = inner_prod(mParallelDirections, mParallelDirections); KRATOS_ERROR_IF(num_parallel_components == 0) << "There is no parallel direction!" << std::endl; - const int num_serial_components = voigt_size - num_parallel_components; + const SizeType num_serial_components = voigt_size - num_parallel_components; - if (rParallelProjector.size1() != voigt_size) + if (rParallelProjector.size1() != voigt_size || rParallelProjector.size2() != num_parallel_components) rParallelProjector.resize(voigt_size, num_parallel_components, false); - if (rSerialProjector.size1() != voigt_size) + if (rSerialProjector.size1() != num_serial_components || rParallelProjector.size2() != voigt_size) rSerialProjector.resize(num_serial_components, voigt_size, false); noalias(rParallelProjector) = ZeroMatrix(voigt_size, num_parallel_components); noalias(rSerialProjector) = ZeroMatrix(num_serial_components, voigt_size); - IndexType parallel_counter = 0, serial_counter = 0; + SizeType parallel_counter = 0, serial_counter = 0; for (IndexType i_comp = 0; i_comp < voigt_size; ++i_comp) { if (mParallelDirections[i_comp] == 1) { rParallelProjector(i_comp, parallel_counter) = 1.0; @@ -598,11 +606,15 @@ void SerialParallelRuleOfMixturesLaw::CalculateGreenLagrangeStrain(ConstitutiveL Matrix F(dimension, dimension); noalias(F) = rValues.GetDeformationGradientF(); - Matrix C_tensor; - C_tensor.resize(dimension, dimension, false); + Matrix C_tensor(dimension, dimension); noalias(C_tensor) = prod(trans(F),F); - ConstitutiveLawUtilities<6>::CalculateGreenLagrangianStrain(C_tensor, r_strain_vector); + const SizeType strain_size = GetStrainSize(); + if (strain_size == 3) { + AdvancedConstitutiveLawUtilities<3>::CalculateGreenLagrangianStrain(C_tensor, r_strain_vector); + } else { + AdvancedConstitutiveLawUtilities<6>::CalculateGreenLagrangianStrain(C_tensor, r_strain_vector); + } } /***********************************************************************************/ @@ -616,8 +628,7 @@ void SerialParallelRuleOfMixturesLaw::CalculateAlmansiStrain(ConstitutiveLaw::Pa Matrix F(dimension, dimension); noalias(F) = rValues.GetDeformationGradientF(); - Matrix B_tensor; - B_tensor.resize(dimension, dimension, false); + Matrix B_tensor(dimension, dimension) noalias(B_tensor) = prod(F, trans(F)); const SizeType strain_size = GetStrainSize(); @@ -626,7 +637,6 @@ void SerialParallelRuleOfMixturesLaw::CalculateAlmansiStrain(ConstitutiveLaw::Pa } else { AdvancedConstitutiveLawUtilities<6>::CalculateAlmansiStrain(B_tensor, r_strain_vector); } - } /***********************************************************************************/ @@ -724,7 +734,7 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK2(ConstitutiveLa // Total strain vector Vector fiber_stress_vector, matrix_stress_vector; - this->IntegrateStrainSerialParallelBehaviour(r_strain_vector, fiber_stress_vector, matrix_stress_vector, r_material_properties, rValues, mPreviousSerialStrainMatrix, ConstitutiveLaw::StressMeasure_PK2); + IntegrateStrainSerialParallelBehaviour(r_strain_vector, fiber_stress_vector, matrix_stress_vector, r_material_properties, rValues, mPreviousSerialStrainMatrix, ConstitutiveLaw::StressMeasure_PK2); // We call the FinalizeMaterialResponse of the matrix and fiber CL auto& r_material_properties = rValues.GetMaterialProperties(); @@ -739,10 +749,10 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK2(ConstitutiveLa values_fiber.SetMaterialProperties(r_props_fiber_cl); Matrix parallel_projector, serial_projector; - this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); + CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); - this->CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); + CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); values_fiber.SetStrainVector(fiber_strain_vector); values_matrix.SetStrainVector(matrix_strain_vector); @@ -788,7 +798,7 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseKirchhoff(Constitu // Total strain vector Vector fiber_stress_vector, matrix_stress_vector; - this->IntegrateStrainSerialParallelBehaviour(r_strain_vector, fiber_stress_vector, matrix_stress_vector, r_material_properties, rValues, mPreviousSerialStrainMatrix, ConstitutiveLaw::StressMeasure_PK2); + IntegrateStrainSerialParallelBehaviour(r_strain_vector, fiber_stress_vector, matrix_stress_vector, r_material_properties, rValues, mPreviousSerialStrainMatrix, ConstitutiveLaw::StressMeasure_PK2); // We call the FinalizeMaterialResponse of the matrix and fiber CL auto& r_material_properties = rValues.GetMaterialProperties(); @@ -803,10 +813,10 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseKirchhoff(Constitu values_fiber.SetMaterialProperties(r_props_fiber_cl); Matrix parallel_projector, serial_projector; - this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); + CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); - this->CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); + CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); values_fiber.SetStrainVector(fiber_strain_vector); values_matrix.SetStrainVector(matrix_strain_vector); @@ -852,7 +862,7 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseCauchy(Constitutiv // Total strain vector Vector fiber_stress_vector, matrix_stress_vector; - this->IntegrateStrainSerialParallelBehaviour(r_strain_vector, fiber_stress_vector, matrix_stress_vector, r_material_properties, rValues, mPreviousSerialStrainMatrix, ConstitutiveLaw::StressMeasure_PK2); + IntegrateStrainSerialParallelBehaviour(r_strain_vector, fiber_stress_vector, matrix_stress_vector, r_material_properties, rValues, mPreviousSerialStrainMatrix, ConstitutiveLaw::StressMeasure_PK2); // We call the FinalizeMaterialResponse of the matrix and fiber CL auto& r_material_properties = rValues.GetMaterialProperties(); @@ -867,10 +877,10 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseCauchy(Constitutiv values_fiber.SetMaterialProperties(r_props_fiber_cl); Matrix parallel_projector, serial_projector; - this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); + CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); - this->CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); + CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); values_fiber.SetStrainVector(fiber_strain_vector); values_matrix.SetStrainVector(matrix_strain_vector); @@ -1431,11 +1441,11 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( // We compute the constitutive matrix if (rThisVariable == CONSTITUTIVE_MATRIX) { - this->CalculateMaterialResponse(rParameterValues, this->GetStressMeasure()); + CalculateMaterialResponse(rParameterValues, this->GetStressMeasure()); } else if (rThisVariable == CONSTITUTIVE_MATRIX_PK2) { - this->CalculateMaterialResponsePK2(rParameterValues); + CalculateMaterialResponsePK2(rParameterValues); } else if (rThisVariable == CONSTITUTIVE_MATRIX_KIRCHHOFF) { - this->CalculateMaterialResponseKirchhoff(rParameterValues); + CalculateMaterialResponseKirchhoff(rParameterValues); } if (rValue.size1() != voigt_size) @@ -1544,11 +1554,11 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( // We compute the stress if (rThisVariable == CAUCHY_STRESS_TENSOR) { - this->CalculateMaterialResponseCauchy(rParameterValues); + CalculateMaterialResponseCauchy(rParameterValues); } else if (rThisVariable == PK2_STRESS_TENSOR) { - this->CalculateMaterialResponsePK2(rParameterValues); + CalculateMaterialResponsePK2(rParameterValues); } else if (rThisVariable == KIRCHHOFF_STRESS_TENSOR) { - this->CalculateMaterialResponseKirchhoff(rParameterValues); + CalculateMaterialResponseKirchhoff(rParameterValues); } if (rValue.size1() != voigt_size) @@ -1561,23 +1571,23 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( return rValue; } else if (rThisVariable == GREEN_LAGRANGE_STRAIN_TENSOR_MATRIX) { Matrix parallel_projector, serial_projector; - this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); + CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector& r_strain_vector = rParameterValues.GetStrainVector(); Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); - this->CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); + CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); if (rValue.size1() != voigt_size) rValue.resize(voigt_size, voigt_size, false); noalias(rValue) = MathUtils::StrainVectorToTensor(matrix_strain_vector); return rValue; } else if (rThisVariable == GREEN_LAGRANGE_STRAIN_TENSOR_FIBER) { - const std::size_t voigt_size = this->GetStrainSize(); + const std::size_t voigt_size = GetStrainSize(); Matrix parallel_projector, serial_projector; - this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); + CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector& r_strain_vector = rParameterValues.GetStrainVector(); Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); - this->CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); + CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); if (rValue.size1() != voigt_size) rValue.resize(voigt_size, voigt_size, false); noalias(rValue) = MathUtils::StrainVectorToTensor(fiber_strain_vector); From 2073d8a96ed90a88fd3a1d2fc2f25cfdc052adf2 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 10:44:28 +0100 Subject: [PATCH 07/16] now compiles --- .../composites/serial_parallel_rule_of_mixtures_law.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 1ec65651a372..1acbd828c62a 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -611,9 +611,9 @@ void SerialParallelRuleOfMixturesLaw::CalculateGreenLagrangeStrain(ConstitutiveL const SizeType strain_size = GetStrainSize(); if (strain_size == 3) { - AdvancedConstitutiveLawUtilities<3>::CalculateGreenLagrangianStrain(C_tensor, r_strain_vector); + ConstitutiveLawUtilities<3>::CalculateGreenLagrangianStrain(C_tensor, r_strain_vector); } else { - AdvancedConstitutiveLawUtilities<6>::CalculateGreenLagrangianStrain(C_tensor, r_strain_vector); + ConstitutiveLawUtilities<6>::CalculateGreenLagrangianStrain(C_tensor, r_strain_vector); } } @@ -628,7 +628,7 @@ void SerialParallelRuleOfMixturesLaw::CalculateAlmansiStrain(ConstitutiveLaw::Pa Matrix F(dimension, dimension); noalias(F) = rValues.GetDeformationGradientF(); - Matrix B_tensor(dimension, dimension) + Matrix B_tensor(dimension, dimension); noalias(B_tensor) = prod(F, trans(F)); const SizeType strain_size = GetStrainSize(); From 872c020ceebf6398144bdde4c9dd6e770d3cd3b9 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 11:11:03 +0100 Subject: [PATCH 08/16] avoid seg fault --- .../composites/serial_parallel_rule_of_mixtures_law.cpp | 6 ++++++ .../composites/serial_parallel_rule_of_mixtures_law.h | 8 +++----- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 1acbd828c62a..504af000be13 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -1410,6 +1410,12 @@ void SerialParallelRuleOfMixturesLaw::InitializeMaterial( mpMatrixConstitutiveLaw->InitializeMaterial(r_props_matrix_cl, rElementGeometry, rShapeFunctionsValues); mpFiberConstitutiveLaw ->InitializeMaterial(r_props_fiber_cl, rElementGeometry, rShapeFunctionsValues); + const SizeType strain_size = GetStrainSize(); + mPreviousStrainVector.resize(strain_size, false); + mPreviousStrainVector.clear(); + mPreviousSerialStrainMatrix.resize(GetNumberOfSerialComponents(), false); + mPreviousSerialStrainMatrix.clear(); + if (rElementGeometry.Has(SERIAL_PARALLEL_IMPOSED_STRAIN)) mIsPrestressed = true; } diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h index cae73bfb6a49..a9f819a18f26 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h @@ -136,7 +136,6 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw */ SizeType WorkingSpaceDimension() override { - KRATOS_DEBUG_ERROR_IF(mpMatrixConstitutiveLaw->WorkingSpaceDimension() != mpFiberConstitutiveLaw->WorkingSpaceDimension()) << "The WorkingSpaceDimension of the fiber and matrix mismatch..." << std::endl; return mpMatrixConstitutiveLaw->WorkingSpaceDimension(); }; @@ -145,7 +144,6 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw */ SizeType GetStrainSize() const override { - KRATOS_DEBUG_ERROR_IF(mpMatrixConstitutiveLaw->GetStrainSize() != mpFiberConstitutiveLaw->GetStrainSize()) << "The GetStrainSize of the fiber and matrix mismatch..." << std::endl; return mpMatrixConstitutiveLaw->GetStrainSize(); }; @@ -686,9 +684,9 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw ConstitutiveLaw::Pointer mpMatrixConstitutiveLaw; ConstitutiveLaw::Pointer mpFiberConstitutiveLaw; double mFiberVolumetricParticipation; - Vector mParallelDirections = ZeroVector(GetStrainSize()); - Vector mPreviousStrainVector = ZeroVector(GetStrainSize()); - Vector mPreviousSerialStrainMatrix = ZeroVector(GetNumberOfSerialComponents()); + Vector mParallelDirections; + Vector mPreviousStrainVector; + Vector mPreviousSerialStrainMatrix; bool mIsPrestressed = false; ///@} From d660886f403de76cb84e8236d50654a0064f5e6e Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 11:38:07 +0100 Subject: [PATCH 09/16] minor --- .../composites/serial_parallel_rule_of_mixtures_law.cpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 504af000be13..53230c791680 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -34,11 +34,8 @@ ConstitutiveLaw::Pointer SerialParallelRuleOfMixturesLaw::Create(Kratos::Paramet if (fiber_volumetric_participation < 0.0 || fiber_volumetric_participation > 1.0) { KRATOS_ERROR << "A wrong fiber volumetric participation has been set: Greater than 1 or lower than 0..." << std::endl; } - const int voigt_size = GetStrainSize(); - Vector parallel_directions(voigt_size); - for (IndexType i_comp = 0; i_comp < voigt_size; ++i_comp) { - parallel_directions[i_comp] = NewParameters["parallel_behaviour_directions"][i_comp].GetInt(); - } + + Vector parallel_directions = NewParameters["parallel_behaviour_directions"].GetVector(); return Kratos::make_shared(fiber_volumetric_participation, parallel_directions); } @@ -979,7 +976,7 @@ Vector& SerialParallelRuleOfMixturesLaw::GetValue( const bool matrix_has = mpMatrixConstitutiveLaw->Has(rThisVariable); const bool fiber_has = mpFiberConstitutiveLaw->Has(rThisVariable); const SizeType voigt_size = GetStrainSize(); - rValue.resize(GetStrainSize(), false); + rValue.resize(voigt_size, false); rValue.clear(); if (matrix_has && fiber_has) { Vector r_vector_matrix(voigt_size), r_vector_fiber(voigt_size); From 92bcc18833ea4479e12061f50752afc59f6a20d9 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 11:42:21 +0100 Subject: [PATCH 10/16] space --- .../composites/serial_parallel_rule_of_mixtures_law.h | 1 - 1 file changed, 1 deletion(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h index a9f819a18f26..a876bf9c7b3b 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h @@ -145,7 +145,6 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw SizeType GetStrainSize() const override { return mpMatrixConstitutiveLaw->GetStrainSize(); - }; /** From 70d792859fc9eceef41d9182a493b0b7908c7b5d Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 13:19:01 +0100 Subject: [PATCH 11/16] adding template to CL --- .../constitutive_laws_application.cpp | 3 +- .../constitutive_laws_application.h | 3 +- .../serial_parallel_rule_of_mixtures_law.cpp | 154 ++++++++++++------ .../serial_parallel_rule_of_mixtures_law.h | 5 + ...add_custom_constitutive_laws_to_python.cpp | 8 +- 5 files changed, 120 insertions(+), 53 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/constitutive_laws_application.cpp b/applications/ConstitutiveLawsApplication/constitutive_laws_application.cpp index f872f72efa07..fb0bcefa4fad 100644 --- a/applications/ConstitutiveLawsApplication/constitutive_laws_application.cpp +++ b/applications/ConstitutiveLawsApplication/constitutive_laws_application.cpp @@ -69,7 +69,8 @@ void KratosConstitutiveLawsApplication::Register() // Custom Constitutive laws // Serial-Parallel Rule Of Mixtures - KRATOS_REGISTER_CONSTITUTIVE_LAW("SerialParallelRuleOfMixturesLaw", mSerialParallelRuleOfMixturesLaw); + KRATOS_REGISTER_CONSTITUTIVE_LAW("SerialParallelRuleOfMixturesLaw3D", mSerialParallelRuleOfMixturesLaw3D); + KRATOS_REGISTER_CONSTITUTIVE_LAW("SerialParallelRuleOfMixturesLaw2D", mSerialParallelRuleOfMixturesLaw2D); // Anisotropic law KRATOS_REGISTER_CONSTITUTIVE_LAW("GenericAnisotropicPlaneStrain2DLaw", mGenericAnisotropicPlaneStrain2DLaw); diff --git a/applications/ConstitutiveLawsApplication/constitutive_laws_application.h b/applications/ConstitutiveLawsApplication/constitutive_laws_application.h index 6c1947ea8412..d7d9cae639c1 100644 --- a/applications/ConstitutiveLawsApplication/constitutive_laws_application.h +++ b/applications/ConstitutiveLawsApplication/constitutive_laws_application.h @@ -301,7 +301,8 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) KratosConstitutiveLawsApplicatio const MultiLinearIsotropicPlaneStress2D mMultiLinearIsotropicPlaneStress2D; // Damage and plasticity laws - const SerialParallelRuleOfMixturesLaw mSerialParallelRuleOfMixturesLaw; + const SerialParallelRuleOfMixturesLaw<3> mSerialParallelRuleOfMixturesLaw3D; + const SerialParallelRuleOfMixturesLaw<2> mSerialParallelRuleOfMixturesLaw2D; const SmallStrainIsotropicPlasticityFactory mSmallStrainIsotropicPlasticityFactory; const SmallStrainKinematicPlasticityFactory mSmallStrainKinematicPlasticityFactory; const FiniteStrainIsotropicPlasticityFactory mFiniteStrainIsotropicPlasticityFactory; diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 53230c791680..63e6b4e692a1 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -28,7 +28,8 @@ namespace Kratos { -ConstitutiveLaw::Pointer SerialParallelRuleOfMixturesLaw::Create(Kratos::Parameters NewParameters) const +template +ConstitutiveLaw::Pointer SerialParallelRuleOfMixturesLaw::Create(Kratos::Parameters NewParameters) const { const double fiber_volumetric_participation = NewParameters["combination_factors"][1].GetDouble(); if (fiber_volumetric_participation < 0.0 || fiber_volumetric_participation > 1.0) { @@ -36,13 +37,14 @@ ConstitutiveLaw::Pointer SerialParallelRuleOfMixturesLaw::Create(Kratos::Paramet } Vector parallel_directions = NewParameters["parallel_behaviour_directions"].GetVector(); - return Kratos::make_shared(fiber_volumetric_participation, parallel_directions); + return Kratos::make_shared>(fiber_volumetric_participation, parallel_directions); } /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponsePK1(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponsePK1(ConstitutiveLaw::Parameters& rValues) { this->InitializeMaterialResponseCauchy(rValues); } @@ -50,7 +52,8 @@ void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponsePK1(Constitutive /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponsePK2(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponsePK2(ConstitutiveLaw::Parameters& rValues) { this->InitializeMaterialResponseCauchy(rValues); } @@ -58,7 +61,8 @@ void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponsePK2(Constitutive /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponseKirchhoff(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponseKirchhoff(ConstitutiveLaw::Parameters& rValues) { this->InitializeMaterialResponseCauchy(rValues); } @@ -66,7 +70,8 @@ void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponseKirchhoff(Consti /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) { auto& r_material_properties = rValues.GetMaterialProperties(); const auto it_cl_begin = r_material_properties.GetSubProperties().begin(); @@ -97,7 +102,8 @@ void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponseCauchy(Constitut /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponsePK1(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponsePK1(ConstitutiveLaw::Parameters& rValues) { this->CalculateMaterialResponsePK2(rValues); @@ -112,7 +118,8 @@ void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponsePK1(ConstitutiveL /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponsePK2(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponsePK2(ConstitutiveLaw::Parameters& rValues) { // Get Values to compute the constitutive law: Flags& r_flags = rValues.GetOptions(); @@ -162,7 +169,8 @@ void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponsePK2(ConstitutiveL /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponseKirchhoff(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponseKirchhoff(ConstitutiveLaw::Parameters& rValues) { // Get Values to compute the constitutive law: Flags& r_flags = rValues.GetOptions(); @@ -223,7 +231,8 @@ void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponseKirchhoff(Constit /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) { this->CalculateMaterialResponseKirchhoff(rValues); @@ -240,7 +249,9 @@ void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponseCauchy(Constituti /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::IntegrateStrainSerialParallelBehaviour( + +template +void SerialParallelRuleOfMixturesLaw::IntegrateStrainSerialParallelBehaviour( const Vector& rStrainVector, Vector& rFiberStressVector, Vector& rMatrixStressVector, @@ -269,16 +280,16 @@ void SerialParallelRuleOfMixturesLaw::IntegrateStrainSerialParallelBehaviour( while (!is_converged && iteration <= max_iterations) { if (iteration == 0) { // Computes an initial approximation of the independent var: rSerialStrainMatrix - this->CalculateInitialApproximationSerialStrainMatrix(rStrainVector, mPreviousStrainVector, rMaterialProperties, parallel_projector, serial_projector, constitutive_tensor_matrix_ss, constitutive_tensor_fiber_ss, rSerialStrainMatrix, rValues, rStressMeasure); + CalculateInitialApproximationSerialStrainMatrix(rStrainVector, mPreviousStrainVector, rMaterialProperties, parallel_projector, serial_projector, constitutive_tensor_matrix_ss, constitutive_tensor_fiber_ss, rSerialStrainMatrix, rValues, rStressMeasure); } // This method computes the strain vector for the matrix & fiber - this->CalculateStrainsOnEachComponent(rStrainVector, parallel_projector, serial_projector, rSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues, iteration); + CalculateStrainsOnEachComponent(rStrainVector, parallel_projector, serial_projector, rSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues, iteration); // This method integrates the stress according to each simple material CL - this->IntegrateStressesOfFiberAndMatrix(rValues, matrix_strain_vector, fiber_strain_vector, rMatrixStressVector, rFiberStressVector, rStressMeasure); + IntegrateStressesOfFiberAndMatrix(rValues, matrix_strain_vector, fiber_strain_vector, rMatrixStressVector, rFiberStressVector, rStressMeasure); // Here we check the convergence of the loop -> serial stresses equilibrium - this->CheckStressEquilibrium(rValues, rStrainVector, serial_projector, rMatrixStressVector, rFiberStressVector, stress_residual, is_converged, constitutive_tensor_matrix_ss, constitutive_tensor_fiber_ss); + CheckStressEquilibrium(rValues, rStrainVector, serial_projector, rMatrixStressVector, rFiberStressVector, stress_residual, is_converged, constitutive_tensor_matrix_ss, constitutive_tensor_fiber_ss); if (is_converged) { break; } else { @@ -292,7 +303,9 @@ void SerialParallelRuleOfMixturesLaw::IntegrateStrainSerialParallelBehaviour( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CorrectSerialStrainMatrix( + +template +void SerialParallelRuleOfMixturesLaw::CorrectSerialStrainMatrix( ConstitutiveLaw::Parameters& rValues, const Vector& rResidualStresses, Vector& rSerialStrainMatrix, @@ -359,7 +372,9 @@ void SerialParallelRuleOfMixturesLaw::CorrectSerialStrainMatrix( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CheckStressEquilibrium( + +template +void SerialParallelRuleOfMixturesLaw::CheckStressEquilibrium( ConstitutiveLaw::Parameters& rValues, const Vector& rStrainVector, const Matrix& rSerialProjector, @@ -411,7 +426,9 @@ void SerialParallelRuleOfMixturesLaw::CheckStressEquilibrium( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::IntegrateStressesOfFiberAndMatrix( + +template +void SerialParallelRuleOfMixturesLaw::IntegrateStressesOfFiberAndMatrix( ConstitutiveLaw::Parameters& rValues, Vector& rMatrixStrainVector, Vector& rFiberStrainVector, @@ -454,7 +471,9 @@ void SerialParallelRuleOfMixturesLaw::IntegrateStressesOfFiberAndMatrix( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CalculateInitialApproximationSerialStrainMatrix( + +template +void SerialParallelRuleOfMixturesLaw::CalculateInitialApproximationSerialStrainMatrix( const Vector& rStrainVector, const Vector& rPreviousStrainVector, const Properties& rMaterialProperties, @@ -528,7 +547,9 @@ void SerialParallelRuleOfMixturesLaw::CalculateInitialApproximationSerialStrainM /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CalculateStrainsOnEachComponent( + +template +void SerialParallelRuleOfMixturesLaw::CalculateStrainsOnEachComponent( const Vector& rStrainVector, const Matrix& rParallelProjector, const Matrix& rSerialProjector, @@ -561,7 +582,9 @@ void SerialParallelRuleOfMixturesLaw::CalculateStrainsOnEachComponent( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CalculateSerialParallelProjectionMatrices( + +template +void SerialParallelRuleOfMixturesLaw::CalculateSerialParallelProjectionMatrices( Matrix& rParallelProjector, Matrix& rSerialProjector ) @@ -595,7 +618,8 @@ void SerialParallelRuleOfMixturesLaw::CalculateSerialParallelProjectionMatrices( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CalculateGreenLagrangeStrain(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::CalculateGreenLagrangeStrain(ConstitutiveLaw::Parameters& rValues) { // Some auxiliary values const SizeType dimension = WorkingSpaceDimension(); @@ -617,7 +641,8 @@ void SerialParallelRuleOfMixturesLaw::CalculateGreenLagrangeStrain(ConstitutiveL /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CalculateAlmansiStrain(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::CalculateAlmansiStrain(ConstitutiveLaw::Parameters& rValues) { // Some auxiliary values const SizeType dimension = WorkingSpaceDimension(); @@ -639,7 +664,8 @@ void SerialParallelRuleOfMixturesLaw::CalculateAlmansiStrain(ConstitutiveLaw::Pa /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK1(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK1(ConstitutiveLaw::Parameters& rValues) { Flags& r_flags = rValues.GetOptions(); // Some auxiliary values @@ -703,7 +729,8 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK1(ConstitutiveLa /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK2(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK2(ConstitutiveLaw::Parameters& rValues) { Flags& r_flags = rValues.GetOptions(); // Some auxiliary values @@ -767,7 +794,8 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK2(ConstitutiveLa /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseKirchhoff(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseKirchhoff(ConstitutiveLaw::Parameters& rValues) { Flags& r_flags = rValues.GetOptions(); // Some auxiliary values @@ -831,7 +859,8 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseKirchhoff(Constitu /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) +template +void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) { Flags& r_flags = rValues.GetOptions(); // Some auxiliary values @@ -895,7 +924,8 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseCauchy(Constitutiv /***********************************************************************************/ /***********************************************************************************/ -bool& SerialParallelRuleOfMixturesLaw::GetValue( +template +bool& SerialParallelRuleOfMixturesLaw::GetValue( const Variable& rThisVariable, bool& rValue ) @@ -915,7 +945,8 @@ bool& SerialParallelRuleOfMixturesLaw::GetValue( /***********************************************************************************/ /***********************************************************************************/ -int& SerialParallelRuleOfMixturesLaw::GetValue( +template +int& SerialParallelRuleOfMixturesLaw::GetValue( const Variable& rThisVariable, int& rValue ) @@ -932,7 +963,8 @@ int& SerialParallelRuleOfMixturesLaw::GetValue( /***********************************************************************************/ /***********************************************************************************/ -double& SerialParallelRuleOfMixturesLaw::GetValue( +template +double& SerialParallelRuleOfMixturesLaw::GetValue( const Variable& rThisVariable, double& rValue ) @@ -968,7 +1000,8 @@ double& SerialParallelRuleOfMixturesLaw::GetValue( /***********************************************************************************/ /***********************************************************************************/ -Vector& SerialParallelRuleOfMixturesLaw::GetValue( +template +Vector& SerialParallelRuleOfMixturesLaw::GetValue( const Variable& rThisVariable, Vector& rValue ) @@ -994,7 +1027,8 @@ Vector& SerialParallelRuleOfMixturesLaw::GetValue( /***********************************************************************************/ /***********************************************************************************/ -Matrix& SerialParallelRuleOfMixturesLaw::GetValue( +template +Matrix& SerialParallelRuleOfMixturesLaw::GetValue( const Variable& rThisVariable, Matrix& rValue ) @@ -1011,7 +1045,8 @@ Matrix& SerialParallelRuleOfMixturesLaw::GetValue( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::SetValue( +template +void SerialParallelRuleOfMixturesLaw::SetValue( const Variable& rThisVariable, const bool& rValue, const ProcessInfo& rCurrentProcessInfo @@ -1032,7 +1067,8 @@ void SerialParallelRuleOfMixturesLaw::SetValue( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::SetValue( +template +void SerialParallelRuleOfMixturesLaw::SetValue( const Variable& rThisVariable, const int& rValue, const ProcessInfo& rCurrentProcessInfo @@ -1049,7 +1085,8 @@ void SerialParallelRuleOfMixturesLaw::SetValue( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::SetValue( +template +void SerialParallelRuleOfMixturesLaw::SetValue( const Variable& rThisVariable, const double& rValue, const ProcessInfo& rCurrentProcessInfo @@ -1069,7 +1106,8 @@ void SerialParallelRuleOfMixturesLaw::SetValue( /***********************************************************************************/ /***********************************************************************************/ -bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) +template +bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) { if (mpMatrixConstitutiveLaw->Has(rThisVariable)) { return true; @@ -1086,7 +1124,8 @@ bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) /***********************************************************************************/ /***********************************************************************************/ -bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) +template +bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) { if (mpMatrixConstitutiveLaw->Has(rThisVariable)) { return true; @@ -1100,7 +1139,8 @@ bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) /***********************************************************************************/ /***********************************************************************************/ -bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) +template +bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) { if (mpMatrixConstitutiveLaw->Has(rThisVariable)) { return true; @@ -1120,7 +1160,8 @@ bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) /***********************************************************************************/ /***********************************************************************************/ -bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) +template +bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) { if (mpMatrixConstitutiveLaw->Has(rThisVariable)) { return true; @@ -1134,7 +1175,8 @@ bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) /***********************************************************************************/ /***********************************************************************************/ -bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) +template +bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) { if (mpMatrixConstitutiveLaw->Has(rThisVariable)) { return true; @@ -1148,7 +1190,8 @@ bool SerialParallelRuleOfMixturesLaw::Has(const Variable& rThisVariable) /***********************************************************************************/ /***********************************************************************************/ -bool& SerialParallelRuleOfMixturesLaw::CalculateValue( +template +bool& SerialParallelRuleOfMixturesLaw::CalculateValue( Parameters& rParameterValues, const Variable& rThisVariable, bool& rValue) @@ -1162,7 +1205,8 @@ bool& SerialParallelRuleOfMixturesLaw::CalculateValue( /***********************************************************************************/ /***********************************************************************************/ -int& SerialParallelRuleOfMixturesLaw::CalculateValue( +template +int& SerialParallelRuleOfMixturesLaw::CalculateValue( Parameters& rParameterValues, const Variable& rThisVariable, int& rValue) @@ -1176,7 +1220,8 @@ int& SerialParallelRuleOfMixturesLaw::CalculateValue( /***********************************************************************************/ /***********************************************************************************/ -double& SerialParallelRuleOfMixturesLaw::CalculateValue( +template +double& SerialParallelRuleOfMixturesLaw::CalculateValue( Parameters& rParameterValues, const Variable& rThisVariable, double& rValue) @@ -1226,7 +1271,8 @@ double& SerialParallelRuleOfMixturesLaw::CalculateValue( /***********************************************************************************/ /***********************************************************************************/ -Vector& SerialParallelRuleOfMixturesLaw::CalculateValue( +template +Vector& SerialParallelRuleOfMixturesLaw::CalculateValue( Parameters& rParameterValues, const Variable& rThisVariable, Vector& rValue) @@ -1390,7 +1436,8 @@ Vector& SerialParallelRuleOfMixturesLaw::CalculateValue( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::InitializeMaterial( +template +void SerialParallelRuleOfMixturesLaw::InitializeMaterial( const Properties& rMaterialProperties, const GeometryType& rElementGeometry, const Vector& rShapeFunctionsValues) @@ -1420,7 +1467,8 @@ void SerialParallelRuleOfMixturesLaw::InitializeMaterial( /***********************************************************************************/ /***********************************************************************************/ -Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( +template +Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( ConstitutiveLaw::Parameters& rParameterValues, const Variable& rThisVariable, Matrix& rValue @@ -1620,18 +1668,20 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( /***********************************************************************************/ /***********************************************************************************/ -void SerialParallelRuleOfMixturesLaw::CalculateTangentTensor( +template +void SerialParallelRuleOfMixturesLaw::CalculateTangentTensor( ConstitutiveLaw::Parameters& rValues, const ConstitutiveLaw::StressMeasure& rStressMeasure) { // Independent from the Voigt Size AdvancedConstitutiveLawUtilities<6>::CalculateTangentTensorByPerturbation(rValues, this, rStressMeasure); } + /***********************************************************************************/ /***********************************************************************************/ - -int SerialParallelRuleOfMixturesLaw::Check( +template +int SerialParallelRuleOfMixturesLaw::Check( const Properties& rMaterialProperties, const GeometryType& rElementGeometry, const ProcessInfo& rCurrentProcessInfo @@ -1653,4 +1703,10 @@ int SerialParallelRuleOfMixturesLaw::Check( return aux_out; } +/***********************************************************************************/ +/***********************************************************************************/ + +template class SerialParallelRuleOfMixturesLaw<2>; +template class SerialParallelRuleOfMixturesLaw<3>; + } // namespace Kratos diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h index a876bf9c7b3b..5c9acb648401 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h @@ -53,6 +53,7 @@ namespace Kratos * @details * @author Alejandro Cornejo */ +template class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw : public ConstitutiveLaw { @@ -68,6 +69,10 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw using BaseType = ConstitutiveLaw; + static constexpr SizeType Dimension = TDim; + + static constexpr SizeType VoigtSize = (TDim == 3) ? 6 : 3; + /// Definition of the machine precision tolerance static constexpr double machine_tolerance = std::numeric_limits::epsilon(); diff --git a/applications/ConstitutiveLawsApplication/custom_python/add_custom_constitutive_laws_to_python.cpp b/applications/ConstitutiveLawsApplication/custom_python/add_custom_constitutive_laws_to_python.cpp index f0764e9d8133..bf94ead50627 100644 --- a/applications/ConstitutiveLawsApplication/custom_python/add_custom_constitutive_laws_to_python.cpp +++ b/applications/ConstitutiveLawsApplication/custom_python/add_custom_constitutive_laws_to_python.cpp @@ -1353,8 +1353,12 @@ void AddCustomConstitutiveLawsToPython(pybind11::module& m) (m,"ParallelRuleOfMixturesLaw2D").def(py::init<>()) ; - py::class_< SerialParallelRuleOfMixturesLaw, typename SerialParallelRuleOfMixturesLaw::Pointer, ConstitutiveLaw > - (m,"SerialParallelRuleOfMixturesLaw").def(py::init<>()) + py::class_< SerialParallelRuleOfMixturesLaw<3>, typename SerialParallelRuleOfMixturesLaw<3>::Pointer, ConstitutiveLaw > + (m,"SerialParallelRuleOfMixturesLaw3D").def(py::init<>()) + ; + + py::class_< SerialParallelRuleOfMixturesLaw<2>, typename SerialParallelRuleOfMixturesLaw<2>::Pointer, ConstitutiveLaw > + (m,"SerialParallelRuleOfMixturesLaw2D").def(py::init<>()) ; py::class_< GenericAnisotropicLaw<3>, typename GenericAnisotropicLaw<3>::Pointer, ConstitutiveLaw > From c132da59362519f438b61e34eb23c0de33edb337 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 13:46:55 +0100 Subject: [PATCH 12/16] now using the temaplte --- .../serial_parallel_rule_of_mixtures_law.cpp | 64 ++++++++----------- .../serial_parallel_rule_of_mixtures_law.h | 10 +-- ...serial_parallel_damage_test_materials.json | 2 +- 3 files changed, 33 insertions(+), 43 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 63e6b4e692a1..8ae7f8b6ae78 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -46,7 +46,7 @@ ConstitutiveLaw::Pointer SerialParallelRuleOfMixturesLaw::Create(Kratos::P template void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponsePK1(ConstitutiveLaw::Parameters& rValues) { - this->InitializeMaterialResponseCauchy(rValues); + InitializeMaterialResponseCauchy(rValues); } /***********************************************************************************/ @@ -55,7 +55,7 @@ void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponsePK1(Consti template void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponsePK2(ConstitutiveLaw::Parameters& rValues) { - this->InitializeMaterialResponseCauchy(rValues); + InitializeMaterialResponseCauchy(rValues); } /***********************************************************************************/ @@ -64,7 +64,7 @@ void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponsePK2(Consti template void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponseKirchhoff(ConstitutiveLaw::Parameters& rValues) { - this->InitializeMaterialResponseCauchy(rValues); + InitializeMaterialResponseCauchy(rValues); } /***********************************************************************************/ @@ -105,7 +105,7 @@ void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponseCauchy(Con template void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponsePK1(ConstitutiveLaw::Parameters& rValues) { - this->CalculateMaterialResponsePK2(rValues); + CalculateMaterialResponsePK2(rValues); if (rValues.IsSetDeterminantF()) { Vector& stress_vector = rValues.GetStressVector(); @@ -208,15 +208,15 @@ void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponseKirchhoff(C if (rValues.IsSetDeterminantF()) { // we push forward the stress - Matrix stress_matrix(3, 3); + Matrix stress_matrix(Dimension, Dimension); noalias(stress_matrix) = MathUtils::StressVectorToTensor(r_integrated_stress_vector); - ContraVariantPushForward (stress_matrix, rValues.GetDeformationGradientF()); //Kirchhoff - noalias(r_integrated_stress_vector) = MathUtils::StressTensorToVector( stress_matrix, r_integrated_stress_vector.size() ); + ContraVariantPushForward(stress_matrix, rValues.GetDeformationGradientF()); // Kirchhoff + noalias(r_integrated_stress_vector) = MathUtils::StressTensorToVector(stress_matrix, r_integrated_stress_vector.size()); } if (flag_const_tensor) { - this->CalculateTangentTensor(rValues, ConstitutiveLaw::StressMeasure_PK2); - // push forward Constitutive tangent tensor + CalculateTangentTensor(rValues, ConstitutiveLaw::StressMeasure_PK2); + // Push forward Constitutive tangent tensor if (rValues.IsSetDeterminantF()) PushForwardConstitutiveMatrix(rValues.GetConstitutiveMatrix(), rValues.GetDeformationGradientF()); } @@ -234,7 +234,7 @@ void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponseKirchhoff(C template void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) { - this->CalculateMaterialResponseKirchhoff(rValues); + CalculateMaterialResponseKirchhoff(rValues); if (rValues.IsSetDeterminantF()) { Vector& stress_vector = rValues.GetStressVector(); @@ -273,8 +273,7 @@ void SerialParallelRuleOfMixturesLaw::IntegrateStrainSerialParallelBehavio bool is_converged = false; int iteration = 0, max_iterations = 150; Vector parallel_strain_matrix(num_parallel_components), stress_residual(num_serial_components); - Matrix constitutive_tensor_matrix_ss(num_serial_components, num_serial_components), - constitutive_tensor_fiber_ss(num_serial_components, num_serial_components); + Matrix constitutive_tensor_matrix_ss(num_serial_components, num_serial_components), constitutive_tensor_fiber_ss(num_serial_components, num_serial_components); // Iterative procedure until the equilibrium is reached in the serial stresses while (!is_converged && iteration <= max_iterations) { @@ -313,7 +312,7 @@ void SerialParallelRuleOfMixturesLaw::CorrectSerialStrainMatrix( const ConstitutiveLaw::StressMeasure& rStressMeasure ) { - const SizeType voigt_size = this->GetStrainSize(); + const SizeType voigt_size = GetStrainSize(); const int num_parallel_components = inner_prod(mParallelDirections, mParallelDirections); const int num_serial_components = voigt_size - num_parallel_components; @@ -630,12 +629,7 @@ void SerialParallelRuleOfMixturesLaw::CalculateGreenLagrangeStrain(Constit Matrix C_tensor(dimension, dimension); noalias(C_tensor) = prod(trans(F),F); - const SizeType strain_size = GetStrainSize(); - if (strain_size == 3) { - ConstitutiveLawUtilities<3>::CalculateGreenLagrangianStrain(C_tensor, r_strain_vector); - } else { - ConstitutiveLawUtilities<6>::CalculateGreenLagrangianStrain(C_tensor, r_strain_vector); - } + ConstitutiveLawUtilities::CalculateGreenLagrangianStrain(C_tensor, r_strain_vector); } /***********************************************************************************/ @@ -653,12 +647,7 @@ void SerialParallelRuleOfMixturesLaw::CalculateAlmansiStrain(ConstitutiveL Matrix B_tensor(dimension, dimension); noalias(B_tensor) = prod(F, trans(F)); - const SizeType strain_size = GetStrainSize(); - if (strain_size == 3) { - AdvancedConstitutiveLawUtilities<3>::CalculateAlmansiStrain(B_tensor, r_strain_vector); - } else { - AdvancedConstitutiveLawUtilities<6>::CalculateAlmansiStrain(B_tensor, r_strain_vector); - } + AdvancedConstitutiveLawUtilities::CalculateAlmansiStrain(B_tensor, r_strain_vector); } /***********************************************************************************/ @@ -708,10 +697,10 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK1(Constitu values_fiber.SetMaterialProperties(r_props_fiber_cl); Matrix parallel_projector, serial_projector; - this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); + CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); - this->CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); + CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); values_fiber.SetStrainVector(fiber_strain_vector); values_matrix.SetStrainVector(matrix_strain_vector); @@ -1229,10 +1218,10 @@ double& SerialParallelRuleOfMixturesLaw::CalculateValue( if (rThisVariable == UNIAXIAL_STRESS_MATRIX) { const SizeType voigt_size = GetStrainSize(); Matrix parallel_projector, serial_projector; - this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); + CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector strain_vector = rParameterValues.GetStrainVector(); Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); - this->CalculateStrainsOnEachComponent(strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); + CalculateStrainsOnEachComponent(strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); const auto &props = rParameterValues.GetMaterialProperties(); const auto it_cl_begin = props.GetSubProperties().begin(); @@ -1248,10 +1237,10 @@ double& SerialParallelRuleOfMixturesLaw::CalculateValue( } else if (rThisVariable == UNIAXIAL_STRESS_FIBER) { const SizeType voigt_size = GetStrainSize(); Matrix parallel_projector, serial_projector; - this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); + CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector strain_vector = rParameterValues.GetStrainVector(); Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); - this->CalculateStrainsOnEachComponent(strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); + CalculateStrainsOnEachComponent(strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); const auto &props = rParameterValues.GetMaterialProperties(); const auto it_cl_begin = props.GetSubProperties().begin(); @@ -1400,13 +1389,13 @@ Vector& SerialParallelRuleOfMixturesLaw::CalculateValue( noalias(rValue) = matrix_strain_vector; return rValue; } else if (rThisVariable == GREEN_LAGRANGE_STRAIN_VECTOR_FIBER) { - const std::size_t voigt_size = this->GetStrainSize(); + const std::size_t voigt_size = GetStrainSize(); Matrix parallel_projector, serial_projector; - this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); + CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector& r_strain_vector = rParameterValues.GetStrainVector(); Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); - this->CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); + CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); if (rValue.size() != voigt_size) rValue.resize(voigt_size, false); noalias(rValue) = fiber_strain_vector; @@ -1582,7 +1571,7 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( Vector& r_strain_vector = rParameterValues.GetStrainVector(); Vector serial_strain_matrix_old = mPreviousSerialStrainMatrix; Vector fiber_stress_vector, matrix_stress_vector; - this->IntegrateStrainSerialParallelBehaviour(r_strain_vector, fiber_stress_vector, matrix_stress_vector, r_material_properties, rParameterValues, serial_strain_matrix_old); + IntegrateStrainSerialParallelBehaviour(r_strain_vector, fiber_stress_vector, matrix_stress_vector, r_material_properties, rParameterValues, serial_strain_matrix_old); // Previous flags restored r_flags.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN, flag_strain); r_flags.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, flag_const_tensor); @@ -1673,8 +1662,7 @@ void SerialParallelRuleOfMixturesLaw::CalculateTangentTensor( ConstitutiveLaw::Parameters& rValues, const ConstitutiveLaw::StressMeasure& rStressMeasure) { - // Independent from the Voigt Size - AdvancedConstitutiveLawUtilities<6>::CalculateTangentTensorByPerturbation(rValues, this, rStressMeasure); + AdvancedConstitutiveLawUtilities::CalculateTangentTensorByPerturbation(rValues, this, rStressMeasure); } /***********************************************************************************/ @@ -1697,6 +1685,8 @@ int SerialParallelRuleOfMixturesLaw::Check( KRATOS_ERROR << "A wrong fiber volumetric participation has been set: Greater than 1 or lower than 0... " << std::to_string(mFiberVolumetricParticipation) << std::endl; aux_out += 1; } + KRATOS_ERROR_IF(VoigtSize != mpMatrixConstitutiveLaw->GetStrainSize()) << "The Strain size of the mpMatrixConstitutiveLaw is not consistent with the Srial Parallel Rule of Mixtures..." << std::endl; + KRATOS_ERROR_IF(VoigtSize != mpFiberConstitutiveLaw->GetStrainSize()) << "The Strain size of the mpFiberConstitutiveLaw is not consistent with the Srial Parallel Rule of Mixtures..." << std::endl; KRATOS_ERROR_IF(mpMatrixConstitutiveLaw->WorkingSpaceDimension() != mpFiberConstitutiveLaw->WorkingSpaceDimension()) << "The WorkingSpaceDimension of the fiber and matrix mismatch..." << std::endl; KRATOS_ERROR_IF(mpMatrixConstitutiveLaw->GetStrainSize() != mpFiberConstitutiveLaw->GetStrainSize()) << "The GetStrainSize of the fiber and matrix mismatch..." << std::endl; diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h index 5c9acb648401..dcf695257d24 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h @@ -141,7 +141,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw */ SizeType WorkingSpaceDimension() override { - return mpMatrixConstitutiveLaw->WorkingSpaceDimension(); + return Dimension; }; /** @@ -149,7 +149,7 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw */ SizeType GetStrainSize() const override { - return mpMatrixConstitutiveLaw->GetStrainSize(); + return VoigtSize; }; /** @@ -688,9 +688,9 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw ConstitutiveLaw::Pointer mpMatrixConstitutiveLaw; ConstitutiveLaw::Pointer mpFiberConstitutiveLaw; double mFiberVolumetricParticipation; - Vector mParallelDirections; - Vector mPreviousStrainVector; - Vector mPreviousSerialStrainMatrix; + Vector mParallelDirections = ZeroVector(VoigtSize); + Vector mPreviousStrainVector = ZeroVector(VoigtSize); + Vector mPreviousSerialStrainMatrix = ZeroVector(GetNumberOfSerialComponents()); bool mIsPrestressed = false; ///@} diff --git a/applications/ConstitutiveLawsApplication/tests/SerialParallelRuleOfMixturesCube/serial_parallel_damage_test_materials.json b/applications/ConstitutiveLawsApplication/tests/SerialParallelRuleOfMixturesCube/serial_parallel_damage_test_materials.json index 4a38a5aa90c3..0f735b3b3721 100644 --- a/applications/ConstitutiveLawsApplication/tests/SerialParallelRuleOfMixturesCube/serial_parallel_damage_test_materials.json +++ b/applications/ConstitutiveLawsApplication/tests/SerialParallelRuleOfMixturesCube/serial_parallel_damage_test_materials.json @@ -4,7 +4,7 @@ "properties_id" : 1, "Material" : { "constitutive_law" : { - "name" : "SerialParallelRuleOfMixturesLaw", + "name" : "SerialParallelRuleOfMixturesLaw3D", "combination_factors" : [0.7, 0.3], "parallel_behaviour_directions" : [1,0,0,0,0,0] }, From 0af5eab6792a32e78e9965090603ad493f983413 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Tue, 7 Jan 2025 16:11:16 +0100 Subject: [PATCH 13/16] final touch :) --- .../composites/serial_parallel_rule_of_mixtures_law.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 8ae7f8b6ae78..f05529a9b07d 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -37,6 +37,7 @@ ConstitutiveLaw::Pointer SerialParallelRuleOfMixturesLaw::Create(Kratos::P } Vector parallel_directions = NewParameters["parallel_behaviour_directions"].GetVector(); + KRATOS_ERROR_IF_NOT(parallel_directions.size() == VoigtSize) << "The parallel_behaviour_directions vector is not consistent with the VoigtSize of the Serial Parallel Rule of Mixtures CL..." << std::endl; return Kratos::make_shared>(fiber_volumetric_participation, parallel_directions); } From 57d92f7b95c9f2012446b0891c5c1b216e6d2008 Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Thu, 9 Jan 2025 15:57:25 +0100 Subject: [PATCH 14/16] surname is BACK --- .../composites/serial_parallel_rule_of_mixtures_law.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index f05529a9b07d..a768a48ffcf6 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -9,7 +9,7 @@ // license: structural_mechanics_application/license.txt // // Main authors: Alejandro Cornejo -// Vicente Mataix +// Vicente Mataix Ferrandiz // // // From d7676154f0c79605f723f9f2eff6c93d26aef6c2 Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Thu, 9 Jan 2025 16:04:57 +0100 Subject: [PATCH 15/16] using the static variable of VoigtSize --- .../serial_parallel_rule_of_mixtures_law.cpp | 147 ++++++++---------- .../serial_parallel_rule_of_mixtures_law.h | 3 +- 2 files changed, 66 insertions(+), 84 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index a768a48ffcf6..0c5fe5214707 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -28,6 +28,10 @@ namespace Kratos { + +/***********************************************************************************/ +/***********************************************************************************/ + template ConstitutiveLaw::Pointer SerialParallelRuleOfMixturesLaw::Create(Kratos::Parameters NewParameters) const { @@ -37,7 +41,7 @@ ConstitutiveLaw::Pointer SerialParallelRuleOfMixturesLaw::Create(Kratos::P } Vector parallel_directions = NewParameters["parallel_behaviour_directions"].GetVector(); - KRATOS_ERROR_IF_NOT(parallel_directions.size() == VoigtSize) << "The parallel_behaviour_directions vector is not consistent with the VoigtSize of the Serial Parallel Rule of Mixtures CL..." << std::endl; + KRATOS_ERROR_IF_NOT(parallel_directions.size() == VoigtSize) << "The parallel_behaviour_directions vector is not consistent with the VoigtSize of the Serial Parallel Rule of Mixtures CL. parallel_directions.size(): " << parallel_directions.size() << " vs VoigtSize: " << VoigtSize << std::endl; return Kratos::make_shared>(fiber_volumetric_participation, parallel_directions); } @@ -79,7 +83,6 @@ void SerialParallelRuleOfMixturesLaw::InitializeMaterialResponseCauchy(Con const auto& r_props_matrix_cl = *(it_cl_begin); const auto& r_props_fiber_cl = *(it_cl_begin + 1); - // Get Values to compute the constitutive law: Flags& r_flags = rValues.GetOptions(); const bool flag_strain = r_flags.Is(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN ); @@ -238,13 +241,13 @@ void SerialParallelRuleOfMixturesLaw::CalculateMaterialResponseCauchy(Cons CalculateMaterialResponseKirchhoff(rValues); if (rValues.IsSetDeterminantF()) { - Vector& stress_vector = rValues.GetStressVector(); - Matrix& constitutive_matrix = rValues.GetConstitutiveMatrix(); + Vector& r_stress_vector = rValues.GetStressVector(); + Matrix& r_constitutive_matrix = rValues.GetConstitutiveMatrix(); const double determinant_f = rValues.GetDeterminantF(); // Set to Cauchy Stress: - stress_vector /= determinant_f; - constitutive_matrix /= determinant_f; + r_stress_vector /= determinant_f; + r_constitutive_matrix /= determinant_f; } } @@ -262,14 +265,13 @@ void SerialParallelRuleOfMixturesLaw::IntegrateStrainSerialParallelBehavio const ConstitutiveLaw::StressMeasure& rStressMeasure ) { - const SizeType voigt_size = GetStrainSize(); const SizeType num_parallel_components = inner_prod(mParallelDirections, mParallelDirections); - const SizeType num_serial_components = voigt_size - num_parallel_components; + const SizeType num_serial_components = VoigtSize - num_parallel_components; - Matrix parallel_projector(voigt_size, num_parallel_components), serial_projector(num_serial_components, voigt_size); + Matrix parallel_projector(VoigtSize, num_parallel_components), serial_projector(num_serial_components, VoigtSize); CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); bool is_converged = false; int iteration = 0, max_iterations = 150; @@ -313,9 +315,8 @@ void SerialParallelRuleOfMixturesLaw::CorrectSerialStrainMatrix( const ConstitutiveLaw::StressMeasure& rStressMeasure ) { - const SizeType voigt_size = GetStrainSize(); const int num_parallel_components = inner_prod(mParallelDirections, mParallelDirections); - const int num_serial_components = voigt_size - num_parallel_components; + const int num_serial_components = VoigtSize - num_parallel_components; auto& r_material_properties = rValues.GetMaterialProperties(); const auto it_cl_begin = r_material_properties.GetSubProperties().begin(); @@ -338,7 +339,7 @@ void SerialParallelRuleOfMixturesLaw::CorrectSerialStrainMatrix( ConstitutiveLaw::Parameters values_fiber = rValues; ConstitutiveLaw::Parameters values_matrix = rValues; - Matrix fiber_tangent_tensor(voigt_size, voigt_size), matrix_tangent_tensor(voigt_size, voigt_size); + Matrix fiber_tangent_tensor(VoigtSize, VoigtSize), matrix_tangent_tensor(VoigtSize, VoigtSize); Matrix fiber_tangent_tensor_ss(num_serial_components, num_serial_components), matrix_tangent_tensor_ss(num_serial_components, num_serial_components); // Compute the tangent tensor of the matrix @@ -437,14 +438,12 @@ void SerialParallelRuleOfMixturesLaw::IntegrateStressesOfFiberAndMatrix( const ConstitutiveLaw::StressMeasure& rStressMeasure ) { - const SizeType voigt_size = GetStrainSize(); - - if (rMatrixStressVector.size() != voigt_size) { - rMatrixStressVector.resize(voigt_size, false); + if (rMatrixStressVector.size() != VoigtSize) { + rMatrixStressVector.resize(VoigtSize, false); } - if (rFiberStressVector.size() != voigt_size) { - rFiberStressVector.resize(voigt_size, false); + if (rFiberStressVector.size() != VoigtSize) { + rFiberStressVector.resize(VoigtSize, false); } auto& r_material_properties = rValues.GetMaterialProperties(); @@ -486,7 +485,6 @@ void SerialParallelRuleOfMixturesLaw::CalculateInitialApproximationSerialS const ConstitutiveLaw::StressMeasure& rStressMeasure ) { - const SizeType voigt_size = GetStrainSize(); const Vector& r_total_strain_vector_parallel = prod(trans(rParallelProjector), rStrainVector); const Vector& r_total_strain_vector_serial = prod(rSerialProjector, rStrainVector); @@ -495,7 +493,7 @@ void SerialParallelRuleOfMixturesLaw::CalculateInitialApproximationSerialS const double k_f = mFiberVolumetricParticipation; const double k_m = 1.0 - mFiberVolumetricParticipation; - Matrix constitutive_tensor_matrix(voigt_size, voigt_size), constitutive_tensor_fiber(voigt_size, voigt_size); + Matrix constitutive_tensor_matrix(VoigtSize, VoigtSize), constitutive_tensor_fiber(VoigtSize, VoigtSize); const auto it_cl_begin = rMaterialProperties.GetSubProperties().begin(); const auto& r_props_matrix_cl = *(it_cl_begin); @@ -589,22 +587,21 @@ void SerialParallelRuleOfMixturesLaw::CalculateSerialParallelProjectionMat Matrix& rSerialProjector ) { - const SizeType voigt_size = GetStrainSize(); const SizeType num_parallel_components = inner_prod(mParallelDirections, mParallelDirections); KRATOS_ERROR_IF(num_parallel_components == 0) << "There is no parallel direction!" << std::endl; - const SizeType num_serial_components = voigt_size - num_parallel_components; + const SizeType num_serial_components = VoigtSize - num_parallel_components; - if (rParallelProjector.size1() != voigt_size || rParallelProjector.size2() != num_parallel_components) - rParallelProjector.resize(voigt_size, num_parallel_components, false); + if (rParallelProjector.size1() != VoigtSize || rParallelProjector.size2() != num_parallel_components) + rParallelProjector.resize(VoigtSize, num_parallel_components, false); - if (rSerialProjector.size1() != num_serial_components || rParallelProjector.size2() != voigt_size) - rSerialProjector.resize(num_serial_components, voigt_size, false); + if (rSerialProjector.size1() != num_serial_components || rParallelProjector.size2() != VoigtSize) + rSerialProjector.resize(num_serial_components, VoigtSize, false); - noalias(rParallelProjector) = ZeroMatrix(voigt_size, num_parallel_components); - noalias(rSerialProjector) = ZeroMatrix(num_serial_components, voigt_size); + noalias(rParallelProjector) = ZeroMatrix(VoigtSize, num_parallel_components); + noalias(rSerialProjector) = ZeroMatrix(num_serial_components, VoigtSize); SizeType parallel_counter = 0, serial_counter = 0; - for (IndexType i_comp = 0; i_comp < voigt_size; ++i_comp) { + for (IndexType i_comp = 0; i_comp < VoigtSize; ++i_comp) { if (mParallelDirections[i_comp] == 1) { rParallelProjector(i_comp, parallel_counter) = 1.0; parallel_counter++; @@ -658,8 +655,6 @@ template void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK1(ConstitutiveLaw::Parameters& rValues) { Flags& r_flags = rValues.GetOptions(); - // Some auxiliary values - const SizeType voigt_size = GetStrainSize(); // In case the element has not computed the Strain if (r_flags.IsNot(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) { @@ -699,7 +694,7 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK1(Constitu Matrix parallel_projector, serial_projector; CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); @@ -723,8 +718,6 @@ template void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK2(ConstitutiveLaw::Parameters& rValues) { Flags& r_flags = rValues.GetOptions(); - // Some auxiliary values - const SizeType voigt_size = GetStrainSize(); // In case the element has not computed the Strain if (r_flags.IsNot(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) { @@ -764,7 +757,7 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponsePK2(Constitu Matrix parallel_projector, serial_projector; CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); @@ -788,8 +781,6 @@ template void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseKirchhoff(ConstitutiveLaw::Parameters& rValues) { Flags& r_flags = rValues.GetOptions(); - // Some auxiliary values - const SizeType voigt_size = GetStrainSize(); // In case the element has not computed the Strain if (r_flags.IsNot(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) { @@ -829,7 +820,7 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseKirchhoff(Co Matrix parallel_projector, serial_projector; CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); @@ -853,8 +844,6 @@ template void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) { Flags& r_flags = rValues.GetOptions(); - // Some auxiliary values - const SizeType voigt_size = GetStrainSize(); // In case the element has not computed the Strain if (r_flags.IsNot(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) { @@ -894,7 +883,7 @@ void SerialParallelRuleOfMixturesLaw::FinalizeMaterialResponseCauchy(Const Matrix parallel_projector, serial_projector; CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rValues); @@ -998,11 +987,10 @@ Vector& SerialParallelRuleOfMixturesLaw::GetValue( { const bool matrix_has = mpMatrixConstitutiveLaw->Has(rThisVariable); const bool fiber_has = mpFiberConstitutiveLaw->Has(rThisVariable); - const SizeType voigt_size = GetStrainSize(); - rValue.resize(voigt_size, false); + rValue.resize(VoigtSize, false); rValue.clear(); if (matrix_has && fiber_has) { - Vector r_vector_matrix(voigt_size), r_vector_fiber(voigt_size); + Vector r_vector_matrix(VoigtSize), r_vector_fiber(VoigtSize); mpMatrixConstitutiveLaw->GetValue(rThisVariable, r_vector_matrix); mpFiberConstitutiveLaw->GetValue(rThisVariable, r_vector_fiber); noalias(rValue) = mFiberVolumetricParticipation * r_vector_fiber + (1.0 - mFiberVolumetricParticipation) * r_vector_matrix; @@ -1217,11 +1205,10 @@ double& SerialParallelRuleOfMixturesLaw::CalculateValue( double& rValue) { if (rThisVariable == UNIAXIAL_STRESS_MATRIX) { - const SizeType voigt_size = GetStrainSize(); Matrix parallel_projector, serial_projector; CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector strain_vector = rParameterValues.GetStrainVector(); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); CalculateStrainsOnEachComponent(strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); const auto &props = rParameterValues.GetMaterialProperties(); @@ -1236,11 +1223,10 @@ double& SerialParallelRuleOfMixturesLaw::CalculateValue( noalias(rParameterValues.GetStrainVector()) = strain_vector; return rValue; } else if (rThisVariable == UNIAXIAL_STRESS_FIBER) { - const SizeType voigt_size = GetStrainSize(); Matrix parallel_projector, serial_projector; CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector strain_vector = rParameterValues.GetStrainVector(); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); CalculateStrainsOnEachComponent(strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); const auto &props = rParameterValues.GetMaterialProperties(); @@ -1267,7 +1253,6 @@ Vector& SerialParallelRuleOfMixturesLaw::CalculateValue( const Variable& rThisVariable, Vector& rValue) { - const SizeType voigt_size = GetStrainSize(); // We do some special operations for constitutive matrices if (rThisVariable == CAUCHY_STRESS_VECTOR_FIBER) { @@ -1305,8 +1290,8 @@ Vector& SerialParallelRuleOfMixturesLaw::CalculateValue( r_flags.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, flag_const_tensor); r_flags.Set(ConstitutiveLaw::COMPUTE_STRESS, flag_stress); - if (rValue.size() != voigt_size) - rValue.resize(voigt_size, false); + if (rValue.size() != VoigtSize) + rValue.resize(VoigtSize, false); noalias(rValue) = fiber_stress_vector; return rValue; @@ -1345,8 +1330,8 @@ Vector& SerialParallelRuleOfMixturesLaw::CalculateValue( r_flags.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN, flag_strain); r_flags.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, flag_const_tensor); r_flags.Set(ConstitutiveLaw::COMPUTE_STRESS, flag_stress); - if (rValue.size() != voigt_size) - rValue.resize(voigt_size, false); + if (rValue.size() != VoigtSize) + rValue.resize(VoigtSize, false); noalias(rValue) = matrix_stress_vector; return rValue; @@ -1370,8 +1355,8 @@ Vector& SerialParallelRuleOfMixturesLaw::CalculateValue( this->CalculateMaterialResponseKirchhoff(rParameterValues); } - if (rValue.size() != voigt_size) - rValue.resize(voigt_size, false); + if (rValue.size() != VoigtSize) + rValue.resize(VoigtSize, false); noalias(rValue) = rParameterValues.GetStressVector(); // Previous flags restored @@ -1383,22 +1368,22 @@ Vector& SerialParallelRuleOfMixturesLaw::CalculateValue( this->CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector& r_strain_vector = rParameterValues.GetStrainVector(); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); this->CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); - if (rValue.size() != voigt_size) - rValue.resize(voigt_size, false); + if (rValue.size() != VoigtSize) + rValue.resize(VoigtSize, false); noalias(rValue) = matrix_strain_vector; return rValue; } else if (rThisVariable == GREEN_LAGRANGE_STRAIN_VECTOR_FIBER) { - const std::size_t voigt_size = GetStrainSize(); + const std::size_t VoigtSize = GetStrainSize(); Matrix parallel_projector, serial_projector; CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector& r_strain_vector = rParameterValues.GetStrainVector(); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); - if (rValue.size() != voigt_size) - rValue.resize(voigt_size, false); + if (rValue.size() != VoigtSize) + rValue.resize(VoigtSize, false); noalias(rValue) = fiber_strain_vector; return rValue; } else { @@ -1464,8 +1449,6 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( Matrix& rValue ) { - const SizeType dimension = WorkingSpaceDimension(); - const SizeType voigt_size = GetStrainSize(); // We do some special operations for constitutive matrices if (rThisVariable == CONSTITUTIVE_MATRIX || rThisVariable == CONSTITUTIVE_MATRIX_PK2 || rThisVariable == CONSTITUTIVE_MATRIX_KIRCHHOFF) { // Get Values to compute the constitutive law: @@ -1489,8 +1472,8 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( CalculateMaterialResponseKirchhoff(rParameterValues); } - if (rValue.size1() != voigt_size) - rValue.resize(voigt_size, voigt_size, false); + if (rValue.size1() != VoigtSize) + rValue.resize(VoigtSize, VoigtSize, false); noalias(rValue) = rParameterValues.GetConstitutiveMatrix(); // Previous flags restored @@ -1498,8 +1481,8 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( r_flags.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, flag_const_tensor); r_flags.Set(ConstitutiveLaw::COMPUTE_STRESS, flag_stress); } else if (rThisVariable == DEFORMATION_GRADIENT) { - if (rValue.size1() != dimension) - rValue.resize(dimension, dimension, false); + if (rValue.size1() != Dimension) + rValue.resize(Dimension, Dimension, false); noalias(rValue) = rParameterValues.GetDeformationGradientF(); } else if (rThisVariable == CAUCHY_STRESS_TENSOR_FIBER) { // TODO: Make in the future modifications for take into account different layers combinations @@ -1537,8 +1520,8 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( r_flags.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, flag_const_tensor); r_flags.Set(ConstitutiveLaw::COMPUTE_STRESS, flag_stress); - if (rValue.size1() != voigt_size) - rValue.resize(voigt_size, voigt_size, false); + if (rValue.size1() != VoigtSize) + rValue.resize(VoigtSize, VoigtSize, false); noalias(rValue) = MathUtils::StressVectorToTensor(fiber_stress_vector); return rValue; @@ -1577,8 +1560,8 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( r_flags.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN, flag_strain); r_flags.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, flag_const_tensor); r_flags.Set(ConstitutiveLaw::COMPUTE_STRESS, flag_stress); - if (rValue.size1() != voigt_size) - rValue.resize(voigt_size, voigt_size, false); + if (rValue.size1() != VoigtSize) + rValue.resize(VoigtSize, VoigtSize, false); noalias(rValue) = MathUtils::StressVectorToTensor(matrix_stress_vector); return rValue; @@ -1602,8 +1585,8 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( CalculateMaterialResponseKirchhoff(rParameterValues); } - if (rValue.size1() != voigt_size) - rValue.resize(voigt_size, voigt_size, false); + if (rValue.size1() != VoigtSize) + rValue.resize(VoigtSize, VoigtSize, false); noalias(rValue) = MathUtils::StressVectorToTensor(rParameterValues.GetStressVector()); // Previous flags restored @@ -1615,22 +1598,22 @@ Matrix& SerialParallelRuleOfMixturesLaw::CalculateValue( CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector& r_strain_vector = rParameterValues.GetStrainVector(); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); - if (rValue.size1() != voigt_size) - rValue.resize(voigt_size, voigt_size, false); + if (rValue.size1() != VoigtSize) + rValue.resize(VoigtSize, VoigtSize, false); noalias(rValue) = MathUtils::StrainVectorToTensor(matrix_strain_vector); return rValue; } else if (rThisVariable == GREEN_LAGRANGE_STRAIN_TENSOR_FIBER) { - const std::size_t voigt_size = GetStrainSize(); + const std::size_t VoigtSize = GetStrainSize(); Matrix parallel_projector, serial_projector; CalculateSerialParallelProjectionMatrices(parallel_projector, serial_projector); const Vector& r_strain_vector = rParameterValues.GetStrainVector(); - Vector matrix_strain_vector(voigt_size), fiber_strain_vector(voigt_size); + Vector matrix_strain_vector(VoigtSize), fiber_strain_vector(VoigtSize); CalculateStrainsOnEachComponent(r_strain_vector, parallel_projector, serial_projector, mPreviousSerialStrainMatrix, matrix_strain_vector, fiber_strain_vector, rParameterValues); - if (rValue.size1() != voigt_size) - rValue.resize(voigt_size, voigt_size, false); + if (rValue.size1() != VoigtSize) + rValue.resize(VoigtSize, VoigtSize, false); noalias(rValue) = MathUtils::StrainVectorToTensor(fiber_strain_vector); return rValue; } else { diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h index dcf695257d24..3b3ac8cab14c 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.h @@ -9,8 +9,7 @@ // license: structural_mechanics_application/license.txt // // Main authors: Alejandro Cornejo -// Vicente Mataix -// +// Vicente Mataix Ferrandiz // // From 666926c805b9d21fbefc607ce752db9c17879d69 Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Thu, 9 Jan 2025 16:13:18 +0100 Subject: [PATCH 16/16] addressing comments --- .../composites/serial_parallel_rule_of_mixtures_law.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp index 0c5fe5214707..c5bf811e3308 100644 --- a/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp +++ b/applications/ConstitutiveLawsApplication/custom_constitutive/composites/serial_parallel_rule_of_mixtures_law.cpp @@ -1669,10 +1669,10 @@ int SerialParallelRuleOfMixturesLaw::Check( KRATOS_ERROR << "A wrong fiber volumetric participation has been set: Greater than 1 or lower than 0... " << std::to_string(mFiberVolumetricParticipation) << std::endl; aux_out += 1; } - KRATOS_ERROR_IF(VoigtSize != mpMatrixConstitutiveLaw->GetStrainSize()) << "The Strain size of the mpMatrixConstitutiveLaw is not consistent with the Srial Parallel Rule of Mixtures..." << std::endl; - KRATOS_ERROR_IF(VoigtSize != mpFiberConstitutiveLaw->GetStrainSize()) << "The Strain size of the mpFiberConstitutiveLaw is not consistent with the Srial Parallel Rule of Mixtures..." << std::endl; - KRATOS_ERROR_IF(mpMatrixConstitutiveLaw->WorkingSpaceDimension() != mpFiberConstitutiveLaw->WorkingSpaceDimension()) << "The WorkingSpaceDimension of the fiber and matrix mismatch..." << std::endl; - KRATOS_ERROR_IF(mpMatrixConstitutiveLaw->GetStrainSize() != mpFiberConstitutiveLaw->GetStrainSize()) << "The GetStrainSize of the fiber and matrix mismatch..." << std::endl; + KRATOS_ERROR_IF(VoigtSize != mpMatrixConstitutiveLaw->GetStrainSize()) << "The Strain size of the mpMatrixConstitutiveLaw: " << mpMatrixConstitutiveLaw->GetStrainSize() << "is not consistent with the Serial Parallel Rule of Mixtures: " << VoigtSize << std::endl; + KRATOS_ERROR_IF(VoigtSize != mpFiberConstitutiveLaw->GetStrainSize()) << "The Strain size of the mpFiberConstitutiveLaw: " << mpFiberConstitutiveLaw->GetStrainSize() << "is not consistent with the Serial Parallel Rule of Mixtures: " << VoigtSize << std::endl; + KRATOS_ERROR_IF(mpMatrixConstitutiveLaw->WorkingSpaceDimension() != mpFiberConstitutiveLaw->WorkingSpaceDimension()) << "The WorkingSpaceDimension of the matrix (" << mpMatrixConstitutiveLaw->WorkingSpaceDimension() << ") and fiber (" << mpFiberConstitutiveLaw->WorkingSpaceDimension() << ") mismatch..." << std::endl; + KRATOS_ERROR_IF(mpMatrixConstitutiveLaw->GetStrainSize() != mpFiberConstitutiveLaw->GetStrainSize()) << "The GetStrainSize of the matrix (" << mpMatrixConstitutiveLaw->GetStrainSize() << ") and fiber (" << mpFiberConstitutiveLaw->GetStrainSize() << ") mismatch..." << std::endl; return aux_out; }