[CLApp] Generalizing the Serial Parallel Rule of Mixtures to perform 2d simulations

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);

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;

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,10 +47,13 @@ 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
+ * The SP-RoM is able to work in 2D and in 3D.
  * @details
  * @author Alejandro Cornejo
  */
+template<unsigned int TDim>
 class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw
     : public ConstitutiveLaw
 {
@@ -59,10 +62,16 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw
     ///@{
 
     /// The node definition
-    typedef Node NodeType;
+    using NodeType = Node;
 
     /// The geometry definition
-    typedef Geometry<NodeType> GeometryType;
+    using GeometryType = Geometry<NodeType>;
+
+    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<double>::epsilon();
@@ -99,7 +108,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) 
     {
@@ -132,15 +141,15 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw
      */
     SizeType WorkingSpaceDimension() override
     {
-        return 3;
+        return Dimension;
     };
 
     /**
      * @brief Voigt tensor size:
      */
     SizeType GetStrainSize() const override
     {
-        return 6;
+        return VoigtSize;
     };
 
     /**
@@ -260,7 +269,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<double>& rThisVariable,
         const double& rValue,
         const ProcessInfo& rCurrentProcessInfo
@@ -539,23 +548,30 @@ 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);
 
     /**
      * @brief If the CL requires to initialize the material response, called by the element in InitializeSolutionStep.
      */
     bool RequiresInitializeMaterialResponse() override
     {
-        return true;
+        if (mpMatrixConstitutiveLaw->RequiresInitializeMaterialResponse() || mpFiberConstitutiveLaw->RequiresInitializeMaterialResponse()) {
+            return true;
+        }
+        return false;
     }
 
     /**
      * @brief If the CL requires to initialize the material response, called by the element in InitializeSolutionStep.
      */
     bool RequiresFinalizeMaterialResponse() override
     {
-        return true;
+        if (mpMatrixConstitutiveLaw->RequiresFinalizeMaterialResponse() || mpFiberConstitutiveLaw->RequiresFinalizeMaterialResponse()) {
+            return true;
+        }
+        return false;
     }
 
     /**
@@ -641,7 +657,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;
     }
 
     ///@}
@@ -672,8 +688,8 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) SerialParallelRuleOfMixturesLaw
     ConstitutiveLaw::Pointer mpMatrixConstitutiveLaw;
     ConstitutiveLaw::Pointer mpFiberConstitutiveLaw;
     double mFiberVolumetricParticipation;
-    array_1d<double, 6> mParallelDirections = ZeroVector(6);
-    array_1d<double, 6> mPreviousStrainVector = ZeroVector(6);
+    Vector mParallelDirections = ZeroVector(VoigtSize);
+    Vector mPreviousStrainVector = ZeroVector(VoigtSize);
     Vector mPreviousSerialStrainMatrix = ZeroVector(GetNumberOfSerialComponents());
     bool mIsPrestressed = false;
 
@@ -703,7 +719,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);
@@ -715,7 +731,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);

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  >

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<TVoigtSize>::CalculateOrthotropicElasticMa
 /***********************************************************************************/
 /***********************************************************************************/
 
+template <SizeType TVoigtSize>
+void AdvancedConstitutiveLawUtilities<TVoigtSize>::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<TangentOperatorEstimation>(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>;

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

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]
             },