Merge pull request #12545 from KratosMultiphysics/dem/improvement_part2

[DEMApplication] Fixed bugs related to the bonded particle models and moment calculation [improvement part 2]

applications/DEMApplication/CMakeLists.txt

@@ -91,6 +91,8 @@ set(KRATOS_DEM_APPLICATION_CORE
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_parallel_bond_CL.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_smooth_joint_CL.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_parallel_bond_for_membrane_CL.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/inlet.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/fast_filling_creator.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/force_based_inlet.cpp

applications/DEMApplication/DEM_application.cpp

@@ -49,6 +49,8 @@
 #include "custom_constitutive/DEM_parallel_bond_CL.h"
 #include "custom_constitutive/DEM_smooth_joint_CL.h"
 #include "custom_constitutive/DEM_parallel_bond_for_membrane_CL.h"
+#include "custom_constitutive/DEM_parallel_bond_bilinear_damage_CL.h"
+#include "custom_constitutive/DEM_parallel_bond_bilinear_damage_mixed_CL.h"
 #include "custom_constitutive/DEM_rolling_friction_model.h"
 #include "custom_constitutive/DEM_rolling_friction_model_constant_torque.h"
 #include "custom_constitutive/DEM_rolling_friction_model_viscous_torque.h"
@@ -107,6 +109,7 @@ KRATOS_CREATE_VARIABLE(bool, GLOBAL_COORDINATION_NUMBER_OPTION)
 KRATOS_CREATE_VARIABLE(bool, AUTOMATIC_SKIN_COMPUTATION)
 KRATOS_CREATE_VARIABLE(double, SKIN_FACTOR_RADIUS)
 KRATOS_CREATE_VARIABLE(int, CLEAN_INDENT_OPTION)
+KRATOS_CREATE_VARIABLE(int, CLEAN_INDENT_V2_OPTION)
 KRATOS_CREATE_VARIABLE(int, TRIHEDRON_OPTION)
 KRATOS_CREATE_VARIABLE(int, ROLLING_FRICTION_OPTION)
 KRATOS_CREATE_VARIABLE(int, POISSON_EFFECT_OPTION)
@@ -134,6 +137,7 @@ KRATOS_CREATE_VARIABLE(double, DEM_ENGINE_PERFORMANCE)
 KRATOS_CREATE_VARIABLE(double, DEM_DRAG_CONSTANT_X)
 KRATOS_CREATE_VARIABLE(double, DEM_DRAG_CONSTANT_Y)
 KRATOS_CREATE_VARIABLE(double, DEM_DRAG_CONSTANT_Z)
+KRATOS_CREATE_VARIABLE(bool, ENERGY_CALCULATION_OPTION)
 
 KRATOS_CREATE_VARIABLE(double, INITIAL_VELOCITY_X_VALUE)
 KRATOS_CREATE_VARIABLE(double, INITIAL_VELOCITY_Y_VALUE)
@@ -253,6 +257,9 @@ KRATOS_CREATE_VARIABLE(bool, DEBUG_PRINTING_OPTION)
 KRATOS_CREATE_VARIABLE(int, DEBUG_PRINTING_ID_1)
 KRATOS_CREATE_VARIABLE(int, DEBUG_PRINTING_ID_2)
 KRATOS_CREATE_VARIABLE(double, FRACTURE_ENERGY)
+KRATOS_CREATE_VARIABLE(double, FRACTURE_ENERGY_NORMAL)
+KRATOS_CREATE_VARIABLE(double, FRACTURE_ENERGY_TANGENTIAL)
+KRATOS_CREATE_VARIABLE(double, FRACTURE_ENERGY_EXPONENT)
 KRATOS_CREATE_VARIABLE(double, SIGMA_SLOPE_CHANGE_THRESHOLD)
 KRATOS_CREATE_VARIABLE(double, INTERNAL_FRICTION_AFTER_THRESHOLD)
 KRATOS_CREATE_VARIABLE(double, SEARCH_RADIUS_INCREMENT_FOR_BONDS_CREATION)
@@ -325,6 +332,12 @@ KRATOS_CREATE_VARIABLE(Quaternion<double>, AUX_ORIENTATION)
 KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(LOCAL_AUX_ANGULAR_VELOCITY)
 // ******************* Quaternion Integration END *******************
 
+// ****************Radius expansion method BEGIN*******************
+KRATOS_CREATE_VARIABLE(bool, IS_RADIUS_EXPANSION)
+KRATOS_CREATE_VARIABLE(double, RADIUS_EXPANSION_RATE)
+KRATOS_CREATE_VARIABLE(double, RADIUS_MULTIPLIER_MAX)
+// *****************Radius expansion method END********************
+
 // FORCE AND MOMENTUM
 KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(CONTACT_IMPULSE)
 KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(PARTICLE_MOMENT)
@@ -347,6 +360,7 @@ KRATOS_CREATE_VARIABLE(double, PARTICLE_GRAVITATIONAL_ENERGY)
 KRATOS_CREATE_VARIABLE(double, PARTICLE_INELASTIC_VISCODAMPING_ENERGY)
 KRATOS_CREATE_VARIABLE(double, PARTICLE_INELASTIC_FRICTIONAL_ENERGY)
 KRATOS_CREATE_VARIABLE(double, PARTICLE_INELASTIC_ROLLING_RESISTANCE_ENERGY)
+KRATOS_CREATE_VARIABLE(double, PARTICLE_MAX_NORMAL_BALL_TO_BALL_FORCE_TIMES_RADIUS)
 KRATOS_CREATE_VARIABLE(int, COMPUTE_ENERGY_OPTION)
 KRATOS_CREATE_VARIABLE(double, GLOBAL_DAMPING)
 KRATOS_CREATE_VARIABLE(double, GLOBAL_VISCOUS_DAMPING)
@@ -539,6 +553,14 @@ template class DEM_compound_constitutive_law_for_PBM<DEM_parallel_bond, DEM_D_He
 template class DEM_compound_constitutive_law_for_PBM<DEM_parallel_bond, DEM_D_Linear_classic>;
 template class DEM_compound_constitutive_law_for_PBM<DEM_parallel_bond, DEM_D_Quadratic>;
 
+template class DEM_compound_constitutive_law_for_PBM<DEM_parallel_bond_bilinear_damage, DEM_D_Hertz_viscous_Coulomb>;
+template class DEM_compound_constitutive_law_for_PBM<DEM_parallel_bond_bilinear_damage, DEM_D_Linear_classic>;
+template class DEM_compound_constitutive_law_for_PBM<DEM_parallel_bond_bilinear_damage, DEM_D_Quadratic>;
+
+template class DEM_compound_constitutive_law_for_PBM<DEM_parallel_bond_bilinear_damage_mixed, DEM_D_Hertz_viscous_Coulomb>;
+template class DEM_compound_constitutive_law_for_PBM<DEM_parallel_bond_bilinear_damage_mixed, DEM_D_Linear_classic>;
+template class DEM_compound_constitutive_law_for_PBM<DEM_parallel_bond_bilinear_damage_mixed, DEM_D_Quadratic>;
+
 void KratosDEMApplication::Register() {
     KRATOS_INFO("") << "\n"
                     << "     KRATOS |  _ \\| ____|  \\/  |  _ \\ __ _  ___| | __      \n"
@@ -593,6 +615,7 @@ void KratosDEMApplication::Register() {
     KRATOS_REGISTER_VARIABLE(AUTOMATIC_SKIN_COMPUTATION)
     KRATOS_REGISTER_VARIABLE(SKIN_FACTOR_RADIUS)
     KRATOS_REGISTER_VARIABLE(CLEAN_INDENT_OPTION)
+    KRATOS_REGISTER_VARIABLE(CLEAN_INDENT_V2_OPTION)
     KRATOS_REGISTER_VARIABLE(TRIHEDRON_OPTION)
     KRATOS_REGISTER_VARIABLE(NEIGH_INITIALIZED)
     KRATOS_REGISTER_VARIABLE(TRIAXIAL_TEST_OPTION)
@@ -616,6 +639,7 @@ void KratosDEMApplication::Register() {
     KRATOS_REGISTER_VARIABLE(DEM_DRAG_CONSTANT_X)
     KRATOS_REGISTER_VARIABLE(DEM_DRAG_CONSTANT_Y)
     KRATOS_REGISTER_VARIABLE(DEM_DRAG_CONSTANT_Z)
+    KRATOS_REGISTER_VARIABLE(ENERGY_CALCULATION_OPTION)
 
     KRATOS_REGISTER_VARIABLE(INITIAL_VELOCITY_X_VALUE)
     KRATOS_REGISTER_VARIABLE(INITIAL_VELOCITY_Y_VALUE)
@@ -729,6 +753,9 @@ void KratosDEMApplication::Register() {
     KRATOS_REGISTER_VARIABLE(DEBUG_PRINTING_ID_1)
     KRATOS_REGISTER_VARIABLE(DEBUG_PRINTING_ID_2)
     KRATOS_REGISTER_VARIABLE(FRACTURE_ENERGY)
+    KRATOS_REGISTER_VARIABLE(FRACTURE_ENERGY_NORMAL)
+    KRATOS_REGISTER_VARIABLE(FRACTURE_ENERGY_TANGENTIAL)
+    KRATOS_REGISTER_VARIABLE(FRACTURE_ENERGY_EXPONENT)
     KRATOS_REGISTER_VARIABLE(SIGMA_SLOPE_CHANGE_THRESHOLD)
     KRATOS_REGISTER_VARIABLE(INTERNAL_FRICTION_AFTER_THRESHOLD)
     KRATOS_REGISTER_VARIABLE(SEARCH_RADIUS_INCREMENT_FOR_BONDS_CREATION)
@@ -808,6 +835,12 @@ void KratosDEMApplication::Register() {
     KRATOS_REGISTER_VARIABLE(LOCAL_AUX_ANGULAR_VELOCITY)
     // ******************* Quaternion Integration END *******************
 
+    // ****************Radius expansion method BEGIN*******************
+    KRATOS_REGISTER_VARIABLE(IS_RADIUS_EXPANSION)
+    KRATOS_REGISTER_VARIABLE(RADIUS_EXPANSION_RATE)
+    KRATOS_REGISTER_VARIABLE(RADIUS_MULTIPLIER_MAX)
+    // *****************Radius expansion method END********************
+
     // FORCE AND MOMENTUM
     KRATOS_REGISTER_3D_VARIABLE_WITH_COMPONENTS(CONTACT_IMPULSE)
     KRATOS_REGISTER_3D_VARIABLE_WITH_COMPONENTS(PARTICLE_MOMENT)
@@ -830,6 +863,7 @@ void KratosDEMApplication::Register() {
     KRATOS_REGISTER_VARIABLE(PARTICLE_INELASTIC_VISCODAMPING_ENERGY)
     KRATOS_REGISTER_VARIABLE(PARTICLE_INELASTIC_FRICTIONAL_ENERGY)
     KRATOS_REGISTER_VARIABLE(PARTICLE_INELASTIC_ROLLING_RESISTANCE_ENERGY)
+    KRATOS_REGISTER_VARIABLE(PARTICLE_MAX_NORMAL_BALL_TO_BALL_FORCE_TIMES_RADIUS)
     KRATOS_REGISTER_VARIABLE(COMPUTE_ENERGY_OPTION)
     KRATOS_REGISTER_VARIABLE(GLOBAL_DAMPING)
     KRATOS_REGISTER_VARIABLE(GLOBAL_VISCOUS_DAMPING)
@@ -1026,6 +1060,8 @@ void KratosDEMApplication::Register() {
     Serializer::Register("DEM_parallel_bond", DEM_parallel_bond());
     Serializer::Register("DEM_smooth_joint", DEM_smooth_joint());
     Serializer::Register("DEM_parallel_bond_for_membrane", DEM_parallel_bond_for_membrane());
+    Serializer::Register("DEM_parallel_bond_bilinear_damage", DEM_parallel_bond_bilinear_damage());
+    Serializer::Register("DEM_parallel_bond_bilinear_damage_mixed", DEM_parallel_bond_bilinear_damage_mixed());
     Serializer::Register("DEMRollingFrictionModelConstantTorque", DEMRollingFrictionModelConstantTorque());
     Serializer::Register("DEMRollingFrictionModelViscousTorque", DEMRollingFrictionModelViscousTorque());
     Serializer::Register("DEMRollingFrictionModelBounded", DEMRollingFrictionModelBounded());

applications/DEMApplication/DEM_application_variables.h

@@ -76,6 +76,7 @@ namespace Kratos
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, bool, AUTOMATIC_SKIN_COMPUTATION)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, SKIN_FACTOR_RADIUS)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, CLEAN_INDENT_OPTION)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, CLEAN_INDENT_V2_OPTION)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, TRIHEDRON_OPTION)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, ROLLING_FRICTION_OPTION)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, POISSON_EFFECT_OPTION)
@@ -102,6 +103,7 @@ namespace Kratos
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, DEM_DRAG_CONSTANT_X)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, DEM_DRAG_CONSTANT_Y)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, DEM_DRAG_CONSTANT_Z)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, bool, ENERGY_CALCULATION_OPTION)
 
 
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, INITIAL_VELOCITY_X_VALUE)
@@ -217,6 +219,9 @@ namespace Kratos
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, DEBUG_PRINTING_ID_1)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, DEBUG_PRINTING_ID_2)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, FRACTURE_ENERGY)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, FRACTURE_ENERGY_NORMAL)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, FRACTURE_ENERGY_TANGENTIAL)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, FRACTURE_ENERGY_EXPONENT)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, SIGMA_SLOPE_CHANGE_THRESHOLD)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, INTERNAL_FRICTION_AFTER_THRESHOLD)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, SEARCH_RADIUS_INCREMENT_FOR_BONDS_CREATION)
@@ -301,6 +306,12 @@ namespace Kratos
   KRATOS_DEFINE_3D_APPLICATION_VARIABLE_WITH_COMPONENTS(DEM_APPLICATION, LOCAL_AUX_ANGULAR_VELOCITY)
   // ******************* Quaternion Integration END *******************
 
+  // ****************Radius expansion method BEGIN*******************
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, bool, IS_RADIUS_EXPANSION)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, RADIUS_EXPANSION_RATE)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, RADIUS_MULTIPLIER_MAX)
+  // *****************Radius expansion method END********************
+
   // FORCE AND MOMENTUM
 
   KRATOS_DEFINE_3D_APPLICATION_VARIABLE_WITH_COMPONENTS(DEM_APPLICATION, PARTICLE_MOMENT)
@@ -324,6 +335,7 @@ namespace Kratos
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, PARTICLE_INELASTIC_VISCODAMPING_ENERGY)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, PARTICLE_INELASTIC_FRICTIONAL_ENERGY)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, PARTICLE_INELASTIC_ROLLING_RESISTANCE_ENERGY)
+  KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, PARTICLE_MAX_NORMAL_BALL_TO_BALL_FORCE_TIMES_RADIUS)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, int, COMPUTE_ENERGY_OPTION)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, GLOBAL_DAMPING)
   KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, GLOBAL_VISCOUS_DAMPING)

applications/DEMApplication/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.cpp

@@ -161,6 +161,7 @@ namespace Kratos {
         ViscoDampingLocalContactForce[0] = - equiv_visco_damp_coeff_tangential * LocalRelVel[0];
         ViscoDampingLocalContactForce[1] = - equiv_visco_damp_coeff_tangential * LocalRelVel[1];
         ViscoDampingLocalContactForce[2] = - equiv_visco_damp_coeff_normal     * LocalRelVel[2];
+
     }
 
     /////////////////////////
@@ -392,4 +393,8 @@ namespace Kratos {
         inelastic_viscodamping_energy += viscodamping_energy;
     }
 
+    double DEM_D_Hertz_viscous_Coulomb::GetTangentialStiffness() {
+        return mKt;
+    }
+
 } // namespace Kratos

applications/DEMApplication/custom_constitutive/DEM_D_Hertz_viscous_Coulomb_CL.h

@@ -129,6 +129,7 @@ namespace Kratos {
         void CalculateInelasticViscodampingEnergyFEM(double& inelastic_viscodamping_energy,
                                                      double ViscoDampingLocalContactForce[3],
                                                      double LocalDeltDisp[3]);
+        double GetTangentialStiffness() override;
 
     private:
 

applications/DEMApplication/custom_constitutive/DEM_D_Linear_classic_CL.cpp

@@ -19,13 +19,15 @@ namespace Kratos{
         return Kratos::make_unique<DEM_D_Linear_classic>();
     }
 
-    void DEM_D_Linear_classic::InitializeContact(SphericParticle* const element1, SphericParticle* const element2, const double indentation) {
+    void DEM_D_Linear_classic::InitializeContact(SphericParticle* const element1, SphericParticle* const element2, const double indentation_particle) {
 
         //Get equivalent Radius
         const double my_radius       = element1->GetRadius();
         const double other_radius    = element2->GetRadius();
         const double radius_sum      = my_radius + other_radius;
         const double min_radius      = std::min(my_radius, other_radius);
+        //const double radius_sum_inv  = 1.0 / radius_sum;
+        //const double equiv_radius    = my_radius * other_radius * radius_sum_inv;
 
         //Get equivalent Young's Modulus
         const double my_young        = element1->GetYoung();
@@ -40,9 +42,10 @@ namespace Kratos{
         //const double other_shear_modulus = 0.5 * other_young / (1.0 + other_poisson);
         //const double equiv_shear = 1.0 / ((2.0 - my_poisson)/my_shear_modulus + (2.0 - other_poisson)/other_shear_modulus);
 
-        //Literature [Cundall, 2004, "A bonded particle model for rock"] [PFC 7.0 manual]
+        //Literature [Cundall, 2004, "A bonded particle model for rock"] [PFC 7.0 manual][mKt Mindlin model]
         mKn = equiv_young * Globals::Pi * min_radius * min_radius / radius_sum;
         mKt = mKn / (2.0 * (equiv_poisson + 1.0));
+        //mKt = 8.0 * equiv_shear * sqrt(equiv_radius * indentation_particle);
     }
 
     void DEM_D_Linear_classic::InitializeContactWithFEM(SphericParticle* const element, Condition* const wall, const double indentation, const double ini_delta) {
@@ -68,4 +71,8 @@ namespace Kratos{
         mKn = equiv_young * Globals::Pi * effective_radius;
         mKt = mKn / (2.0 * (equiv_poisson + 1.0));
     }
+
+    double DEM_D_Linear_classic::GetTangentialStiffness() {
+        return mKt;
+    }
 } // namespace Kratos

applications/DEMApplication/custom_constitutive/DEM_D_Linear_classic_CL.h

@@ -33,6 +33,8 @@ namespace Kratos {
 
         void InitializeContactWithFEM(SphericParticle* const element, Condition* const wall, const double indentation, const double ini_delta = 0.0) override;
 
+        double GetTangentialStiffness() override;
+
     protected:
 
     };

applications/DEMApplication/custom_constitutive/DEM_D_Linear_viscous_Coulomb_CL.cpp

@@ -170,6 +170,7 @@ namespace Kratos {
         ViscoDampingLocalContactForce[0] = - equiv_visco_damp_coeff_tangential * LocalRelVel[0];
         ViscoDampingLocalContactForce[1] = - equiv_visco_damp_coeff_tangential * LocalRelVel[1];
         ViscoDampingLocalContactForce[2] = - equiv_visco_damp_coeff_normal     * LocalRelVel[2];
+
     }
 
     /////////////////////////

applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.cpp

@@ -98,4 +98,8 @@ namespace Kratos {
         mKt = mKn / (2.0 * (equiv_poisson + 1.0));
     }
 
+    double DEM_D_Quadratic::GetTangentialStiffness() {
+        return mKt;
+    }
+
 } //namespace Kratos

applications/DEMApplication/custom_constitutive/DEM_D_Quadratic_CL.h

@@ -45,6 +45,8 @@ namespace Kratos{
 
         void InitializeContactWithFEM(SphericParticle* const element, Condition* const wall, const double indentation, const double ini_delta = 0.0) override;
 
+        double GetTangentialStiffness() override;
+
     private:
 
     }; //CLASS DEM_D_QUADRATIC

applications/DEMApplication/custom_constitutive/DEM_Dempack_CL.cpp

@@ -146,6 +146,7 @@ namespace Kratos {
                                     double equiv_young,
                                     double equiv_shear,
                                     double indentation,
+                                    double indentation_particle,
                                     double calculation_area,
                                     double& acumulated_damage,
                                     SphericContinuumParticle* element1,
@@ -164,6 +165,7 @@ namespace Kratos {
                 kn_el,
                 equiv_young,
                 indentation,
+                indentation_particle,
                 calculation_area,
                 acumulated_damage,
                 element1,
@@ -214,6 +216,7 @@ namespace Kratos {
             const double kn_el,
             double equiv_young,
             double indentation,
+            double indentation_particle,
             double calculation_area,
             double& acumulated_damage,
             SphericContinuumParticle* element1,
@@ -505,7 +508,8 @@ namespace Kratos {
                     double ElasticLocalRotationalMoment[3], 
                     double ViscoLocalRotationalMoment[3], 
                     double equiv_poisson, 
-                    double indentation, 
+                    double indentation,
+                    double indentation_particle,
                     double LocalElasticContactForce[3],
                     double normalLocalContactForce,
                     double GlobalElasticContactForces[3],

applications/DEMApplication/custom_constitutive/DEM_Dempack_CL.h

@@ -80,6 +80,7 @@ namespace Kratos {
                                     double equiv_young,
                                     double equiv_shear,
                                     double indentation,
+                                    double indentation_particle,
                                     double calculation_area,
                                     double& acumulated_damage,
                                     SphericContinuumParticle* element1,
@@ -96,6 +97,7 @@ namespace Kratos {
                                             const double kn_el,
                                             double equiv_young,
                                             double indentation,
+                                            double indentation_particle,
                                             double calculation_area,
                                             double& acumulated_damage,
                                             SphericContinuumParticle* element1,
@@ -134,7 +136,8 @@ namespace Kratos {
                                       double ElasticLocalRotationalMoment[3], 
                                       double ViscoLocalRotationalMoment[3], 
                                       double equiv_poisson, 
-                                      double indentation, 
+                                      double indentation,
+                                      double indentation_particle, 
                                       double LocalElasticContactForce[3],
                                       double normalLocalContactForce,
                                       double GlobalElasticContactForces[3],

applications/DEMApplication/custom_constitutive/DEM_Dempack_dev_CL.cpp

@@ -77,6 +77,7 @@ namespace Kratos {
             const double kn_el,
             double equiv_young,
             double indentation,
+            double indentation_particle,
             double calculation_area,
             double& acumulated_damage,
             SphericContinuumParticle* element1,

applications/DEMApplication/custom_constitutive/DEM_Dempack_dev_CL.h

@@ -69,6 +69,7 @@ namespace Kratos {
                 const double kn_el,
                 double equiv_young,
                 double indentation,
+                double indentation_particle,
                 double calculation_area,
                 double& acumulated_damage,
                 SphericContinuumParticle* element1,

applications/DEMApplication/custom_constitutive/DEM_ExponentialHC_CL.cpp

@@ -35,6 +35,7 @@ namespace Kratos {
             const double kn_el,
             double equiv_young,
             double indentation,
+            double indentation_particle,
             double calculation_area,
             double& acumulated_damage,
             SphericContinuumParticle* element1,