1 /* ---------------------------------------------------------------------- 2 This is the 3 4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗ 5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝ 6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗ 7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║ 8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║ 9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝® 10 11 DEM simulation engine, released by 12 DCS Computing Gmbh, Linz, Austria 13 http://www.dcs-computing.com, office@dcs-computing.com 14 15 LIGGGHTS® is part of CFDEM®project: 16 http://www.liggghts.com | http://www.cfdem.com 17 18 Core developer and main author: 19 Christoph Kloss, christoph.kloss@dcs-computing.com 20 21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public 22 License, version 2 or later. It is distributed in the hope that it will 23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have 25 received a copy of the GNU General Public License along with LIGGGHTS®. 26 If not, see http://www.gnu.org/licenses . See also top-level README 27 and LICENSE files. 28 29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH, 30 the producer of the LIGGGHTS® software and the CFDEM®coupling software 31 See http://www.cfdem.com/terms-trademark-policy for details. 32 33 ------------------------------------------------------------------------- 34 Contributing author and copyright for this file: 35 36 Arno Mayrhofer (CFDEMresearch GmbH, Linz) 37 Christoph Kloss (DCS Computing GmbH, Linz) 38 Christoph Kloss (JKU Linz) 39 Richard Berger (JKU Linz) 40 41 Copyright 2016- CFDEMresearch GmbH, Linz 42 Copyright 2012- DCS Computing GmbH, Linz 43 Copyright 2009-2012 JKU Linz 44 ------------------------------------------------------------------------- */ 45 46 #include "global_properties.h" 47 #include <cmath> 48 #include <cstring> 49 #include <algorithm> 50 #include "update.h" 51 #include "atom.h" 52 #include "error.h" 53 #include "neighbor.h" 54 #include "modify.h" 55 #include "force.h" 56 57 using namespace std; 58 59 namespace MODEL_PARAMS 60 { 61 static const char * COALESCENCE_MODEL_SWITCHES = "coalescenceModelSwitches"; 62 static const char * COALESCENCE_MODEL_SETTINGS = "coalescenceModelSettings"; 63 static const char * COHESION_DISTANCE_SETTINGS = "cohesionDistanceSettings"; 64 static const char * COHESION_MODEL_SWITCHES = "cohesionModelSwitches"; 65 static const char * COHESION_ENERGY_DENSITY = "cohesionEnergyDensity"; 66 static const char * CHARACTERISTIC_VELOCITY = "characteristicVelocity"; 67 static const char * YOUNGS_MODULUS = "youngsModulus"; 68 static const char * POISSONS_RATIO = "poissonsRatio"; 69 static const char * COEFFICIENT_RESTITUTION = "coefficientRestitution"; 70 static const char * COEFFICIENT_RESTITUTION_LOG = "coefficientRestitutionLog"; 71 static const char * COEFFICIENT_FRICTION = "coefficientFriction"; 72 static const char * COEFFICIENT_ROLL_FRICTION = "coefficientRollingFriction"; 73 static const char * COEFFICIENT_ROLL_VISCOUS_DAMPING = "coefficientRollingViscousDamping"; 74 static const char * FLUID_VISCOSITY = "FluidViscosity"; 75 static const char * MAXIMUM_RESTITUTION = "MaximumRestitution"; 76 static const char * CRITITCAL_STOKES = "CriticalStokes"; 77 static const char * LIQUID_VOLUME = "liquidVolume"; 78 static const char * LIQUID_DENSITY = "liquidDensity"; 79 static const char * HISTORY_INDEX = "historyIndex"; 80 static const char * SURFACE_TENSION = "surfaceTension"; 81 static const char * SWITCH_MODEL = "switchModel"; 82 static const char * CONTACT_ANGLE = "contactAngle"; 83 static const char * COEFFICIENT_MAX_ELASTIC_STIFFNESS = "coefficientMaxElasticStiffness"; 84 static const char * COEFFICIENT_ADHESION_STIFFNESS = "coefficientAdhesionStiffness"; 85 static const char * COEFFICIENT_PLASTICITY_DEPTH = "coefficientPlasticityDepth"; 86 static const char * ROUGHNESS_ABSOLUTE = "roughnessAbsolute"; 87 static const char * ROUGHNESS_RELATIVE = "roughnessRelative"; 88 static const char * ROLLING_STIFFNESS = "rollingStiffness"; 89 static const char * NORMAL_DAMPING_COEFFICIENT = "normalDampingCoefficient"; 90 static const char * TANGENTIAL_STIFFNESS = "tangentialStiffness"; 91 static const char * INITIAL_COHESIVE_STRESS = "initialCohesiveStress"; 92 static const char * MAX_COHESIVE_STRESS = "maxCohesiveStress"; 93 static const char * COHESION_STRENGTH = "cohesionStrength"; 94 static const char * PULL_OFF_FORCE = "pullOffForce"; 95 static const char * OVERLAP_EXPONENT = "overlapExponent"; 96 static const char * ADHESION_EXPONENT = "adhesionExponent"; 97 static const char * SURFACE_ENERGY = "surfaceEnergy"; 98 static const char * TANGENTIAL_MULTIPLIER = "tangentialMultiplier"; 99 static const char * YIELD_RATIO = "coefficientYieldRatio"; 100 static const char * FRICTION_VISCOSITY = "FrictionViscosity"; 101 static const char * COEFFICIENT_FRICTION_STIFFNESS = "coeffFrictionStiffness"; 102 static const char * COEFFICIENT_ROLLING_FRICTION_STIFFNESS = "coeffRollingStiffness"; 103 static const char * UNLOADING_STIFFNESS = "UnloadingStiffness"; 104 static const char * LOADING_STIFFNESS = "LoadingStiffness"; 105 /* ----------------------------------------------------------------------- 106 * Utility functions 107 * ----------------------------------------------------------------------- */ 108 createScalarProperty(PropertyRegistry & registry,const char * name,const char * caller)109 ScalarProperty* createScalarProperty(PropertyRegistry & registry, const char* name, const char * caller) 110 { 111 ScalarProperty * scalar = new ScalarProperty(); 112 FixPropertyGlobal * property = registry.getGlobalProperty(name,"property/global","scalar",0,0,caller); 113 scalar->data = property->compute_scalar(); 114 return scalar; 115 } 116 createScalarProperty(PropertyRegistry & registry,const char * name,const char * caller,bool sanity_checks,const double lo,const double hi)117 ScalarProperty* createScalarProperty(PropertyRegistry & registry, const char* name, const char * caller, bool sanity_checks, const double lo, const double hi) 118 { 119 LAMMPS * lmp = registry.getLAMMPS(); 120 121 ScalarProperty * scalar = new ScalarProperty(); 122 FixPropertyGlobal * property = registry.getGlobalProperty(name,"property/global","scalar",0,0,caller); 123 124 const double value = property->compute_scalar(); 125 if (sanity_checks) 126 { 127 if (value < lo || value > hi) 128 { 129 char buf[200]; 130 sprintf(buf,"%s requires values between %g and %g \n",name,lo,hi); 131 lmp->error->all(FLERR,buf); 132 } 133 } 134 135 scalar->data = value; 136 return scalar; 137 } 138 createPerTypeProperty(PropertyRegistry & registry,const char * name,const char * caller)139 VectorProperty* createPerTypeProperty(PropertyRegistry & registry, const char* name, const char * caller) 140 { 141 const int max_type = registry.max_type(); 142 143 VectorProperty * vector = new VectorProperty(max_type+1); 144 FixPropertyGlobal * property = registry.getGlobalProperty(name,"property/global","vector",max_type,0,caller); 145 for(int i = 1; i < max_type+1; i++) 146 vector->data[i] = property->compute_vector(i-1); 147 148 return vector; 149 } 150 createPerTypeProperty(PropertyRegistry & registry,const char * name,const char * caller,bool sanity_checks,const double lo,const double hi)151 VectorProperty* createPerTypeProperty(PropertyRegistry & registry, const char* name, const char * caller, bool sanity_checks, const double lo, const double hi) 152 { 153 LAMMPS * lmp = registry.getLAMMPS(); 154 const int max_type = registry.max_type(); 155 156 VectorProperty * vector = new VectorProperty(max_type+1); 157 FixPropertyGlobal * property = registry.getGlobalProperty(name,"property/global","vector",max_type,0,caller); 158 for(int i = 1; i < max_type+1; i++) 159 { 160 const double value = property->compute_vector(i-1); 161 if (sanity_checks) 162 { 163 if (value < lo || value > hi) 164 { 165 char buf[200]; 166 sprintf(buf,"%s requires values between %g and %g \n",name,lo,hi); 167 lmp->error->all(FLERR,buf); 168 } 169 } 170 vector->data[i] = value; 171 } 172 173 return vector; 174 } 175 createPerTypePairProperty(PropertyRegistry & registry,const char * name,const char * caller)176 MatrixProperty* createPerTypePairProperty(PropertyRegistry & registry, const char * name, const char * caller) 177 { 178 const int max_type = registry.max_type(); 179 180 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 181 FixPropertyGlobal * property = registry.getGlobalProperty(name,"property/global","peratomtypepair",max_type,max_type,caller); 182 183 for(int i = 1; i < max_type+1; i++) 184 { 185 for(int j = 1; j < max_type+1; j++) 186 { 187 matrix->data[i][j] = property->compute_array(i-1,j-1); 188 } 189 } 190 191 return matrix; 192 } 193 createPerTypePairProperty(PropertyRegistry & registry,const char * name,const char * caller,bool sanity_checks,const double lo,const double hi)194 MatrixProperty* createPerTypePairProperty(PropertyRegistry & registry, const char * name, const char * caller, bool sanity_checks, const double lo, const double hi) 195 { 196 LAMMPS * lmp = registry.getLAMMPS(); 197 const int max_type = registry.max_type(); 198 199 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 200 FixPropertyGlobal * property = registry.getGlobalProperty(name,"property/global","peratomtypepair",max_type,max_type,caller); 201 202 for(int i = 1; i < max_type+1; i++) 203 { 204 for(int j = 1; j < max_type+1; j++) 205 { 206 const double value = property->compute_array(i-1,j-1); 207 if (sanity_checks) 208 { 209 if (value < lo || value > hi) 210 { 211 char buf[200]; 212 sprintf(buf,"%s requires values between %g and %g \n",name,lo,hi); 213 lmp->error->all(FLERR,buf); 214 } 215 } 216 matrix->data[i][j] = value; 217 } 218 } 219 220 return matrix; 221 } 222 223 /* ----------------------------------------------------------------------- 224 * Property Creators 225 * ----------------------------------------------------------------------- */ 226 createCharacteristicVelocity(PropertyRegistry & registry,const char * caller,bool sanity_checks)227 ScalarProperty* createCharacteristicVelocity(PropertyRegistry & registry, const char * caller, bool sanity_checks) 228 { 229 LAMMPS * lmp = registry.getLAMMPS(); 230 ScalarProperty* charVelScalar = createScalarProperty(registry, CHARACTERISTIC_VELOCITY, caller); 231 double charVel = charVelScalar->data; 232 233 if(sanity_checks) 234 { 235 if(strcmp(lmp->update->unit_style,"si") == 0 && charVel < 1e-2) 236 lmp->error->all(FLERR,"characteristicVelocity >= 1e-2 required for SI units"); 237 } 238 239 return charVelScalar; 240 } 241 242 /* ---------------------------------------------------------------------- */ createCoalescenceModelSwitches(PropertyRegistry & registry,const char * caller,bool sanity_checks)243 VectorProperty* createCoalescenceModelSwitches(PropertyRegistry & registry, const char * caller, bool sanity_checks) 244 { 245 LAMMPS * lmp = registry.getLAMMPS(); 246 int numSettings = 3; //use 2 settings, index = 0: for bubble-bubble, index = 1: for bubble-wall, 2: decideMethod 247 VectorProperty * vec = new VectorProperty(numSettings); 248 FixPropertyGlobal * ca = registry.getGlobalProperty(COALESCENCE_MODEL_SWITCHES,"property/global","vector",numSettings,0,caller); 249 250 for(int i=0; i <numSettings; i++) 251 { 252 const double aSetting = ca->compute_vector(i); 253 254 // error checks on v 255 if(sanity_checks) 256 { 257 if(aSetting < 0.0) 258 lmp->error->all(FLERR,"model switches for coalescence model must be all positive"); 259 } 260 vec->data[i] = aSetting; 261 } 262 263 return vec; 264 } 265 /* ---------------------------------------------------------------------- */ createCoalescenceModelSettings(PropertyRegistry & registry,const char * caller,bool sanity_checks)266 VectorProperty* createCoalescenceModelSettings(PropertyRegistry & registry, const char * caller, bool sanity_checks) 267 { 268 LAMMPS * lmp = registry.getLAMMPS(); 269 int numSettings = 6; //use 6 settings, starting index = 0 270 VectorProperty * vec = new VectorProperty(numSettings); 271 FixPropertyGlobal * ca = registry.getGlobalProperty(COALESCENCE_MODEL_SETTINGS,"property/global","vector",numSettings,0,caller); 272 273 for(int i=0; i <numSettings; i++) 274 { 275 const double aSetting = ca->compute_vector(i); 276 277 // error checks on v 278 if(sanity_checks) 279 { 280 if(aSetting < 0.0) 281 lmp->error->all(FLERR,"model settings for coalescence model must be all positive"); 282 } 283 284 vec->data[i] = aSetting; 285 } 286 287 return vec; 288 } 289 /* ---------------------------------------------------------------------- */ 290 createCohesionEnergyDensity(PropertyRegistry & registry,const char * caller,bool sanity_checks)291 MatrixProperty* createCohesionEnergyDensity(PropertyRegistry & registry, const char * caller, bool sanity_checks) 292 { 293 return createPerTypePairProperty(registry, COHESION_ENERGY_DENSITY, caller); 294 } 295 296 /* ---------------------------------------------------------------------- */ 297 createCohesionDistanceSettings(PropertyRegistry & registry,const char * caller,bool sanity_checks)298 VectorProperty * createCohesionDistanceSettings(PropertyRegistry & registry, const char * caller, bool sanity_checks) 299 { 300 LAMMPS * lmp = registry.getLAMMPS(); 301 int numSettings = 4; //use 4 settings, starting index = 0 302 VectorProperty * vec = new VectorProperty(numSettings); 303 FixPropertyGlobal * ca = registry.getGlobalProperty(COHESION_DISTANCE_SETTINGS,"property/global","vector",numSettings,0,caller); 304 305 for(int i=0; i <numSettings; i++) 306 { 307 const double aSetting = ca->compute_vector(i); 308 309 // error checks on v 310 if(sanity_checks) 311 { 312 if(aSetting < 0.0) 313 lmp->error->all(FLERR,"distance settings for cohesion model must be all positive"); 314 } 315 316 vec->data[i] = aSetting; 317 } 318 319 return vec; 320 } 321 322 /* ---------------------------------------------------------------------- */ 323 createCohesionModelSwitches(PropertyRegistry & registry,const char * caller,bool sanity_checks)324 VectorProperty * createCohesionModelSwitches(PropertyRegistry & registry, const char * caller, bool sanity_checks) 325 { 326 LAMMPS * lmp = registry.getLAMMPS(); 327 int numSettings = 2; //use 2 settings, index = 0: for particle-particle, index = 1: for particle-wall 328 VectorProperty * vec = new VectorProperty(numSettings); 329 FixPropertyGlobal * ca = registry.getGlobalProperty(COHESION_MODEL_SWITCHES,"property/global","vector",numSettings,0,caller); 330 331 for(int i=0; i <numSettings; i++) 332 { 333 const double aSetting = ca->compute_vector(i); 334 335 // error checks on v 336 if(sanity_checks) 337 { 338 if(aSetting < 0.0) 339 lmp->error->all(FLERR,"model switches for cohesion model must be all positive"); 340 } 341 vec->data[i] = aSetting; 342 } 343 344 return vec; 345 } 346 /* ---------------------------------------------------------------------- */ createYoungsModulus(PropertyRegistry & registry,const char * caller,bool sanity_checks)347 VectorProperty * createYoungsModulus(PropertyRegistry & registry, const char * caller, bool sanity_checks) 348 { 349 LAMMPS * lmp = registry.getLAMMPS(); 350 const int max_type = registry.max_type(); 351 352 VectorProperty * vec = new VectorProperty(max_type+1); 353 FixPropertyGlobal * Y = registry.getGlobalProperty(YOUNGS_MODULUS,"property/global","peratomtype",max_type,0,caller); 354 355 for(int i=1; i < max_type+1; i++) 356 { 357 const double Yi = Y->compute_vector(i-1); 358 359 // error checks on Y 360 361 if(sanity_checks && (0 == lmp->modify->n_fixes_style("bubble")) && (!registry.getLAMMPS()->atom->get_properties()->allow_soft_particles())) 362 { 363 if(strcmp(lmp->update->unit_style,"si") == 0 && Yi < 5e6) 364 lmp->error->all(FLERR,"youngsModulus >= 5e6 required for SI units (use command 'soft_particles yes' to override)"); 365 if(strcmp(lmp->update->unit_style,"cgs") == 0 && Yi < 5e7) 366 lmp->error->all(FLERR,"youngsModulus >= 5e7 required for CGS units (use command 'soft_particles yes' to override)"); 367 } 368 if(sanity_checks && !registry.getLAMMPS()->atom->get_properties()->allow_hard_particles()) 369 { 370 if(strcmp(lmp->update->unit_style,"si") == 0 && Yi > 1e9) 371 lmp->error->all(FLERR,"youngsModulus <= 1e9 required for SI units (use command 'hard_particles yes' to override)"); 372 if(strcmp(lmp->update->unit_style,"cgs") == 0 && Yi > 1e10) 373 lmp->error->all(FLERR,"youngsModulus <= 1e10 required for CGS units (use command 'hard_particles yes' to override)"); 374 } 375 376 vec->data[i] = Yi; 377 } 378 379 return vec; 380 } 381 382 /* ---------------------------------------------------------------------- */ 383 createPoissonsRatio(PropertyRegistry & registry,const char * caller,bool sanity_checks)384 VectorProperty * createPoissonsRatio(PropertyRegistry & registry, const char * caller, bool sanity_checks) 385 { 386 LAMMPS * lmp = registry.getLAMMPS(); 387 const int max_type = registry.max_type(); 388 389 VectorProperty * vec = new VectorProperty(max_type+1); 390 FixPropertyGlobal * v = registry.getGlobalProperty(POISSONS_RATIO,"property/global","peratomtype",max_type,0,caller); 391 392 for(int i=1; i < max_type+1; i++) 393 { 394 const double vi = v->compute_vector(i-1); 395 396 // error checks on v 397 if(sanity_checks) 398 { 399 if(vi < 0. || vi > 0.5) 400 lmp->error->all(FLERR,"0 <= poissonsRatio <= 0.5 required"); 401 } 402 403 vec->data[i] = vi; 404 } 405 406 return vec; 407 } 408 409 /* ---------------------------------------------------------------------- */ 410 createYeff(PropertyRegistry & registry,const char * caller,bool)411 MatrixProperty * createYeff(PropertyRegistry & registry, const char * caller, bool) 412 { 413 const int max_type = registry.max_type(); 414 415 registry.registerProperty(YOUNGS_MODULUS, &createYoungsModulus); 416 registry.registerProperty(POISSONS_RATIO, &createPoissonsRatio); 417 418 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 419 VectorProperty * youngsModulus = registry.getVectorProperty(YOUNGS_MODULUS,caller); 420 VectorProperty * poissonRatio = registry.getVectorProperty(POISSONS_RATIO,caller); 421 double * Y = youngsModulus->data; 422 double * v = poissonRatio->data; 423 424 for(int i=1;i< max_type+1; i++) 425 { 426 for(int j=1;j<max_type+1;j++) 427 { 428 const double Yi=Y[i]; 429 const double Yj=Y[j]; 430 const double vi=v[i]; 431 const double vj=v[j]; 432 matrix->data[i][j] = 1./((1.-pow(vi,2.))/Yi+(1.-pow(vj,2.))/Yj); 433 } 434 } 435 436 return matrix; 437 } 438 439 /* ---------------------------------------------------------------------- */ 440 createGeff(PropertyRegistry & registry,const char * caller,bool)441 MatrixProperty * createGeff(PropertyRegistry & registry, const char * caller, bool) 442 { 443 const int max_type = registry.max_type(); 444 445 registry.registerProperty(YOUNGS_MODULUS, &createYoungsModulus); 446 registry.registerProperty(POISSONS_RATIO, &createPoissonsRatio); 447 448 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 449 VectorProperty * youngsModulus = registry.getVectorProperty(YOUNGS_MODULUS,caller); 450 VectorProperty * poissonRatio = registry.getVectorProperty(POISSONS_RATIO,caller); 451 double * Y = youngsModulus->data; 452 double * v = poissonRatio->data; 453 454 for(int i=1;i< max_type+1; i++) 455 { 456 for(int j=1;j<max_type+1;j++) 457 { 458 const double Yi=Y[i]; 459 const double Yj=Y[j]; 460 const double vi=v[i]; 461 const double vj=v[j]; 462 463 matrix->data[i][j] = 1./(2.*(2.-vi)*(1.+vi)/Yi+2.*(2.-vj)*(1.+vj)/Yj); 464 } 465 } 466 467 return matrix; 468 } 469 470 /* ---------------------------------------------------------------------- */ 471 createCoeffRest(PropertyRegistry & registry,const char * caller,bool sanity_checks)472 MatrixProperty * createCoeffRest(PropertyRegistry & registry, const char * caller, bool sanity_checks) 473 { 474 LAMMPS * lmp = registry.getLAMMPS(); 475 const int max_type = registry.max_type(); 476 477 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 478 FixPropertyGlobal * coeffRest = registry.getGlobalProperty(COEFFICIENT_RESTITUTION,"property/global","peratomtypepair",max_type,max_type,caller); 479 480 for(int i=1;i< max_type+1; i++) 481 { 482 for(int j=1;j<max_type+1;j++) 483 { 484 const double cor = coeffRest->compute_array(i-1,j-1); 485 486 if(sanity_checks) 487 { 488 if(cor <= 0.05 || cor > 1) { 489 lmp->error->all(FLERR,"0.05 < coefficientRestitution <= 1 required"); 490 } 491 } 492 493 matrix->data[i][j] = cor; 494 } 495 } 496 497 return matrix; 498 } 499 500 /* ---------------------------------------------------------------------- */ 501 createCoeffRestLog(PropertyRegistry & registry,const char * caller,bool)502 MatrixProperty * createCoeffRestLog(PropertyRegistry & registry, const char * caller, bool) 503 { 504 const int max_type = registry.max_type(); 505 506 registry.registerProperty(COEFFICIENT_RESTITUTION, &createCoeffRest); 507 508 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 509 MatrixProperty * coeffRestProp = registry.getMatrixProperty(COEFFICIENT_RESTITUTION,caller); 510 double ** coeffRest = coeffRestProp->data; 511 512 for(int i=1;i< max_type+1; i++) 513 { 514 for(int j=1;j<max_type+1;j++) 515 { 516 matrix->data[i][j] = log(coeffRest[i][j]); 517 } 518 } 519 520 return matrix; 521 } 522 523 /* ---------------------------------------------------------------------- */ 524 createBetaEff(PropertyRegistry & registry,const char * caller,bool)525 MatrixProperty * createBetaEff(PropertyRegistry & registry, const char * caller, bool) 526 { 527 const int max_type = registry.max_type(); 528 529 registry.registerProperty(COEFFICIENT_RESTITUTION_LOG, &createCoeffRestLog); 530 531 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 532 MatrixProperty * coeffRestLogProp = registry.getMatrixProperty(COEFFICIENT_RESTITUTION_LOG,caller); 533 double ** coeffRestLog = coeffRestLogProp->data; 534 535 for(int i=1;i< max_type+1; i++) 536 { 537 for(int j=1;j<max_type+1;j++) 538 { 539 matrix->data[i][j] = coeffRestLog[i][j] / sqrt(pow(coeffRestLog[i][j],2.)+pow(M_PI,2.)); 540 } 541 } 542 543 return matrix; 544 } 545 546 /* ---------------------------------------------------------------------- */ 547 createCoeffFrict(PropertyRegistry & registry,const char * caller,bool)548 MatrixProperty* createCoeffFrict(PropertyRegistry & registry, const char * caller, bool) 549 { 550 return createPerTypePairProperty(registry, COEFFICIENT_FRICTION, caller); 551 } 552 553 /* ---------------------------------------------------------------------- */ 554 createCoeffRollFrict(PropertyRegistry & registry,const char * caller,bool)555 MatrixProperty* createCoeffRollFrict(PropertyRegistry & registry, const char * caller, bool) 556 { 557 return createPerTypePairProperty(registry, COEFFICIENT_ROLL_FRICTION, caller); 558 } 559 560 /* ---------------------------------------------------------------------- */ 561 createCoeffRollVisc(PropertyRegistry & registry,const char * caller,bool)562 MatrixProperty* createCoeffRollVisc(PropertyRegistry & registry, const char * caller, bool) 563 { 564 return createPerTypePairProperty(registry, COEFFICIENT_ROLL_VISCOUS_DAMPING, caller); 565 } 566 567 /* ---------------------------------------------------------------------- */ 568 createCoeffMu(PropertyRegistry & registry,const char * caller,bool sanity_checks)569 MatrixProperty * createCoeffMu(PropertyRegistry & registry, const char * caller, bool sanity_checks) 570 { 571 LAMMPS * lmp = registry.getLAMMPS(); 572 const int max_type = registry.max_type(); 573 574 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 575 FixPropertyGlobal * coeffMu = registry.getGlobalProperty(FLUID_VISCOSITY,"property/global","peratomtypepair",max_type,max_type,caller); 576 577 for(int i=1;i< max_type+1; i++) 578 { 579 for(int j=1;j<max_type+1;j++) 580 { 581 const double mu = coeffMu->compute_array(i-1,j-1); 582 583 if(sanity_checks) 584 { 585 if(mu <= 0.) 586 lmp->error->all(FLERR,"coeffMu > 0 required"); 587 } 588 589 matrix->data[i][j] = mu; 590 } 591 } 592 593 return matrix; 594 } 595 596 /* ---------------------------------------------------------------------- */ 597 createCoeffRestMax(PropertyRegistry & registry,const char * caller,bool sanity_checks)598 MatrixProperty * createCoeffRestMax(PropertyRegistry & registry, const char * caller, bool sanity_checks) 599 { 600 LAMMPS * lmp = registry.getLAMMPS(); 601 const int max_type = registry.max_type(); 602 603 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 604 FixPropertyGlobal * coeffRestMax = registry.getGlobalProperty(MAXIMUM_RESTITUTION,"property/global","peratomtypepair",max_type,max_type,caller); 605 606 for(int i=1;i< max_type+1; i++) 607 { 608 for(int j=1;j<max_type+1;j++) 609 { 610 const double restMax = coeffRestMax->compute_array(i-1,j-1); 611 612 if(sanity_checks) 613 { 614 if(restMax <= 0. || restMax > 1.0) 615 lmp->error->all(FLERR,"0 < MaximumRestitution <= 1 required"); 616 } 617 618 matrix->data[i][j] = restMax; 619 } 620 } 621 622 return matrix; 623 } 624 625 /* ---------------------------------------------------------------------- */ 626 createCoeffStc(PropertyRegistry & registry,const char * caller,bool sanity_checks)627 MatrixProperty * createCoeffStc(PropertyRegistry & registry, const char * caller, bool sanity_checks) 628 { 629 LAMMPS * lmp = registry.getLAMMPS(); 630 const int max_type = registry.max_type(); 631 632 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 633 FixPropertyGlobal * coeffStc = registry.getGlobalProperty(CRITITCAL_STOKES,"property/global","peratomtypepair",max_type,max_type,caller); 634 635 for(int i=1;i< max_type+1; i++) 636 { 637 for(int j=1;j<max_type+1;j++) 638 { 639 const double stc= coeffStc->compute_array(i-1,j-1); 640 641 if(sanity_checks) 642 { 643 if(stc <= 0.) 644 lmp->error->all(FLERR,"CriticalStokes > 0 required"); 645 } 646 647 matrix->data[i][j] = stc; 648 } 649 } 650 651 return matrix; 652 } 653 654 /* ---------------------------------------------------------------------- */ 655 createRollingStiffness(PropertyRegistry & registry,const char * caller,bool sanity_checks)656 ScalarProperty* createRollingStiffness(PropertyRegistry & registry, const char * caller, bool sanity_checks) 657 { 658 LAMMPS * lmp = registry.getLAMMPS(); 659 ScalarProperty* rollingStiffness = createScalarProperty(registry, ROLLING_STIFFNESS, caller); 660 double rollStiffness = rollingStiffness->data; 661 662 if(sanity_checks) 663 { 664 if(strcmp(lmp->update->unit_style,"si") == 0 && rollStiffness < 1.0 ) 665 lmp->error->all(FLERR,"rollingStiffness >= 1 required for SI units"); 666 } 667 668 return rollingStiffness; 669 } 670 671 /* ---------------------------------------------------------------------- */ 672 createLiquidVolume(PropertyRegistry & registry,const char * caller,bool)673 ScalarProperty* createLiquidVolume(PropertyRegistry & registry, const char * caller, bool) 674 { 675 return createScalarProperty(registry, LIQUID_VOLUME, caller); 676 } 677 678 /* ---------------------------------------------------------------------- */ 679 createLiquidDensity(PropertyRegistry & registry,const char * caller,bool)680 ScalarProperty* createLiquidDensity(PropertyRegistry & registry, const char * caller, bool) 681 { 682 return createScalarProperty(registry, LIQUID_DENSITY, caller); 683 } 684 685 /* ---------------------------------------------------------------------- */ 686 createSurfaceTension(PropertyRegistry & registry,const char * caller,bool)687 ScalarProperty* createSurfaceTension(PropertyRegistry & registry, const char * caller, bool) 688 { 689 return createScalarProperty(registry, SURFACE_TENSION, caller); 690 } 691 692 /* ---------------------------------------------------------------------- */ 693 createSwitchModel(PropertyRegistry & registry,const char * caller,bool)694 ScalarProperty* createSwitchModel(PropertyRegistry & registry, const char * caller, bool) 695 { 696 return createScalarProperty(registry, SWITCH_MODEL, caller); 697 } 698 699 /* ---------------------------------------------------------------------- */ 700 createHistoryIndex(PropertyRegistry & registry,const char * caller,bool)701 ScalarProperty* createHistoryIndex(PropertyRegistry & registry, const char * caller, bool) 702 { 703 return createScalarProperty(registry, HISTORY_INDEX, caller); 704 } 705 706 /* ---------------------------------------------------------------------- */ 707 createContactAngle(PropertyRegistry & registry,const char * caller,bool sanity_checks)708 VectorProperty * createContactAngle(PropertyRegistry & registry, const char * caller, bool sanity_checks) 709 { 710 LAMMPS * lmp = registry.getLAMMPS(); 711 const int max_type = registry.max_type(); 712 713 VectorProperty * vec = new VectorProperty(max_type+1); 714 FixPropertyGlobal * ca = registry.getGlobalProperty(CONTACT_ANGLE,"property/global","peratomtype",max_type,0,caller); 715 716 for(int i=1; i < max_type+1; i++) 717 { 718 const double vi = ca->compute_vector(i-1); 719 720 // error checks on v 721 if(sanity_checks) 722 { 723 if(vi < 0. || vi > 180.) 724 lmp->error->all(FLERR,"0 <= contactAngle <= 180° required"); 725 } 726 727 vec->data[i] = vi * M_PI / 180.; // from grad to rad 728 } 729 730 return vec; 731 } 732 733 /* ---------------------------------------------------------------------- */ 734 createKn(PropertyRegistry & registry,const char * caller,bool)735 MatrixProperty* createKn(PropertyRegistry & registry, const char * caller, bool) 736 { 737 LAMMPS * lmp = registry.getLAMMPS(); 738 const int max_type = registry.max_type(); 739 740 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 741 FixPropertyGlobal * k_n1 = registry.getGlobalProperty("kn","property/global","peratomtypepair",max_type,max_type,caller); 742 743 for(int i=1;i< max_type+1; i++) 744 { 745 const double cg_i = lmp->force->cg(i); 746 for(int j=1;j<max_type+1;j++) 747 { 748 const double cg_j = lmp->force->cg(j); 749 if(cg_i != cg_j) 750 lmp->error->all(FLERR,"per-type coarse-graining factors must be equal when property 'kn' is used"); 751 const double k_n = cg_i*k_n1->compute_array(i-1,j-1); 752 matrix->data[i][j] = k_n; 753 } 754 } 755 756 return matrix; 757 } 758 759 /* ---------------------------------------------------------------------- */ 760 createKt(PropertyRegistry & registry,const char * caller,bool)761 MatrixProperty* createKt(PropertyRegistry & registry, const char * caller, bool) 762 { 763 return createPerTypePairProperty(registry, "kt", caller); 764 } 765 766 /* ---------------------------------------------------------------------- */ 767 createGamman(PropertyRegistry & registry,const char * caller,bool)768 MatrixProperty* createGamman(PropertyRegistry & registry, const char * caller, bool) 769 { 770 LAMMPS * lmp = registry.getLAMMPS(); 771 const int max_type = registry.max_type(); 772 773 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 774 FixPropertyGlobal * gamma_n1 = registry.getGlobalProperty("gamman","property/global","peratomtypepair",max_type,max_type,caller); 775 776 for(int i=1;i< max_type+1; i++) 777 { 778 const double cg_i = lmp->force->cg(i); 779 for(int j=1;j<max_type+1;j++) 780 { 781 const double cg_j = lmp->force->cg(j); 782 if(cg_i != cg_j) 783 lmp->error->all(FLERR,"per-type coarse-graining factors must be equal when property 'gamman' is used"); 784 const double k_n = (1./cg_i)*gamma_n1->compute_array(i-1,j-1); 785 matrix->data[i][j] = k_n; 786 } 787 } 788 789 return matrix; 790 } 791 792 /* ---------------------------------------------------------------------- */ 793 createGammat(PropertyRegistry & registry,const char * caller,bool)794 MatrixProperty* createGammat(PropertyRegistry & registry, const char * caller, bool) 795 { 796 return createPerTypePairProperty(registry, "gammat", caller); 797 } 798 799 /* ---------------------------------------------------------------------- */ 800 createGammanAbs(PropertyRegistry & registry,const char * caller,bool)801 MatrixProperty* createGammanAbs(PropertyRegistry & registry, const char * caller, bool) 802 { 803 LAMMPS * lmp = registry.getLAMMPS(); 804 const int max_type = registry.max_type(); 805 806 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 807 FixPropertyGlobal * gamma_n1 = registry.getGlobalProperty("gamman_abs","property/global","peratomtypepair",max_type,max_type,caller); 808 809 for(int i=1;i< max_type+1; i++) 810 { 811 const double cg_i = lmp->force->cg(i); 812 for(int j=1;j<max_type+1;j++) 813 { 814 const double cg_j = lmp->force->cg(j); 815 if(cg_i != cg_j) 816 lmp->error->all(FLERR,"per-type coarse-graining factors must be equal when property 'gamman_abs' is used"); 817 const double gamma = cg_i*cg_i*gamma_n1->compute_array(i-1,j-1); 818 matrix->data[i][j] = gamma; 819 } 820 } 821 822 return matrix; 823 } 824 825 /* ---------------------------------------------------------------------- */ 826 createGammatAbs(PropertyRegistry & registry,const char * caller,bool)827 MatrixProperty* createGammatAbs(PropertyRegistry & registry, const char * caller, bool) 828 { 829 return createPerTypePairProperty(registry, "gammat_abs", caller); 830 } 831 832 /* ---------------------------------------------------------------------- */ 833 createCoeffMaxElasticStiffness(PropertyRegistry & registry,const char * caller,bool)834 MatrixProperty* createCoeffMaxElasticStiffness(PropertyRegistry & registry, const char * caller, bool) 835 { 836 return createPerTypePairProperty(registry, COEFFICIENT_MAX_ELASTIC_STIFFNESS, caller); 837 } 838 839 /* ---------------------------------------------------------------------- */ 840 createCoeffAdhesionStiffness(PropertyRegistry & registry,const char * caller,bool)841 MatrixProperty* createCoeffAdhesionStiffness(PropertyRegistry & registry, const char * caller, bool) 842 { 843 return createPerTypePairProperty(registry, COEFFICIENT_ADHESION_STIFFNESS, caller); 844 } 845 846 /* ---------------------------------------------------------------------- */ 847 createCoeffPlasticityDepth(PropertyRegistry & registry,const char * caller,bool)848 MatrixProperty* createCoeffPlasticityDepth(PropertyRegistry & registry, const char * caller, bool) 849 { 850 return createPerTypePairProperty(registry, COEFFICIENT_PLASTICITY_DEPTH, caller); 851 } 852 853 /* ---------------------------------------------------------------------- */ 854 createRoughnessAbsolute(PropertyRegistry & registry,const char * caller,bool)855 ScalarProperty* createRoughnessAbsolute(PropertyRegistry & registry, const char * caller, bool) 856 { 857 return createScalarProperty(registry, ROUGHNESS_ABSOLUTE, caller); 858 } 859 860 /* ---------------------------------------------------------------------- */ createPullOffForce(PropertyRegistry & registry,const char * caller,bool)861 MatrixProperty* createPullOffForce(PropertyRegistry & registry, const char * caller, bool) 862 { 863 return createPerTypePairProperty(registry, PULL_OFF_FORCE, caller); 864 } 865 createRoughnessRelative(PropertyRegistry & registry,const char * caller,bool)866 ScalarProperty* createRoughnessRelative(PropertyRegistry & registry, const char * caller, bool) 867 { 868 return createScalarProperty(registry, ROUGHNESS_RELATIVE, caller); 869 } 870 /* ---------------------------------------------------------------------- */ 871 createNormalDampingCoefficient(PropertyRegistry & registry,const char * caller,bool)872 MatrixProperty* createNormalDampingCoefficient(PropertyRegistry & registry, const char * caller, bool) 873 { 874 return createPerTypePairProperty(registry, NORMAL_DAMPING_COEFFICIENT, caller); 875 } 876 /* ---------------------------------------------------------------------- */ 877 createTangentialDampingCoefficient(PropertyRegistry & registry,const char * caller,bool)878 MatrixProperty* createTangentialDampingCoefficient(PropertyRegistry & registry, const char * caller, bool) 879 { 880 return createPerTypePairProperty(registry, NORMAL_DAMPING_COEFFICIENT, caller); 881 } 882 883 /* ---------------------------------------------------------------------- */ 884 createTangentialStiffness(PropertyRegistry & registry,const char * caller,bool)885 MatrixProperty* createTangentialStiffness(PropertyRegistry & registry, const char * caller, bool) 886 { 887 return createPerTypePairProperty(registry, TANGENTIAL_STIFFNESS, caller); 888 } 889 /* ---------------------------------------------------------------------- */ 890 createInitialCohesiveStress(PropertyRegistry & registry,const char * caller,bool)891 MatrixProperty* createInitialCohesiveStress(PropertyRegistry & registry, const char * caller, bool) 892 { 893 return createPerTypePairProperty(registry, INITIAL_COHESIVE_STRESS, caller); 894 } 895 /* ---------------------------------------------------------------------- */ createOverlapExponent(PropertyRegistry & registry,const char * caller,bool)896 ScalarProperty* createOverlapExponent(PropertyRegistry & registry, const char * caller, bool) 897 { 898 return createScalarProperty(registry, OVERLAP_EXPONENT, caller); 899 } 900 createAdhesionExponent(PropertyRegistry & registry,const char * caller,bool)901 ScalarProperty* createAdhesionExponent(PropertyRegistry & registry, const char * caller, bool) 902 { 903 return createScalarProperty(registry, ADHESION_EXPONENT, caller); 904 } 905 906 /* ---------------------------------------------------------------------- */ 907 createSurfaceEnergy(PropertyRegistry & registry,const char * caller,bool)908 MatrixProperty* createSurfaceEnergy(PropertyRegistry & registry, const char * caller, bool) 909 { 910 return createPerTypePairProperty(registry, SURFACE_ENERGY, caller); 911 } 912 913 /* ---------------------------------------------------------------------- */ 914 createTangentialMultiplier(PropertyRegistry & registry,const char * caller,bool)915 MatrixProperty* createTangentialMultiplier(PropertyRegistry & registry, const char * caller, bool) 916 { 917 return createPerTypePairProperty(registry, TANGENTIAL_MULTIPLIER, caller); 918 } 919 920 /* ---------------------------------------------------------------------- */ 921 922 // Thornton & Ning 923 createYieldRatio(PropertyRegistry & registry,const char * caller,bool sanity_checks)924 VectorProperty * createYieldRatio(PropertyRegistry & registry, const char * caller, bool sanity_checks) 925 { 926 LAMMPS * lmp = registry.getLAMMPS(); 927 const int max_type = registry.max_type(); 928 929 VectorProperty * vec = new VectorProperty(max_type+1); 930 FixPropertyGlobal * v = registry.getGlobalProperty(YIELD_RATIO,"property/global","peratomtype",max_type,0,caller); 931 932 for(int i=1; i < max_type+1; i++) 933 { 934 const double vi = v->compute_vector(i-1); 935 936 // error checks on v 937 if(sanity_checks) 938 { 939 if(vi < 0. || vi > 1) 940 lmp->error->all(FLERR,"0 <= poissonsRatio <= 1 required"); 941 } 942 943 vec->data[i] = vi; 944 } 945 946 return vec; 947 } 948 949 /* ----------------------Luding's tangential parameter-frictional viscous dampening------------------------------------------------ */ 950 createCoeffFricVisc(PropertyRegistry & registry,const char * caller,bool sanity_checks)951 MatrixProperty * createCoeffFricVisc(PropertyRegistry & registry, const char * caller, bool sanity_checks) 952 { 953 LAMMPS * lmp = registry.getLAMMPS(); 954 const int max_type = registry.max_type(); 955 956 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 957 FixPropertyGlobal * coeffFricVisc = registry.getGlobalProperty(FRICTION_VISCOSITY,"property/global","peratomtypepair",max_type,max_type,caller); 958 959 for(int i=1;i< max_type+1; i++) 960 { 961 for(int j=1;j<max_type+1;j++) 962 { 963 const double mu = coeffFricVisc->compute_array(i-1,j-1); 964 965 if(sanity_checks) 966 { 967 if(mu <= 0.) 968 lmp->error->all(FLERR,"coeffFricVisc > 0 required"); 969 } 970 971 matrix->data[i][j] = mu; 972 } 973 } 974 975 return matrix; 976 } 977 978 /* ---------------------------------------------------------------------- */ 979 createCoeffFrictionStiffness(PropertyRegistry & registry,const char * caller,bool)980 MatrixProperty* createCoeffFrictionStiffness(PropertyRegistry & registry, const char * caller, bool) 981 { 982 return createPerTypePairProperty(registry, COEFFICIENT_FRICTION_STIFFNESS, caller); 983 } 984 createCoeffRollingStiffness(PropertyRegistry & registry,const char * caller,bool)985 MatrixProperty* createCoeffRollingStiffness(PropertyRegistry & registry, const char * caller, bool) 986 { 987 return createPerTypePairProperty(registry, COEFFICIENT_ROLLING_FRICTION_STIFFNESS, caller); 988 } 989 990 /* --------------------- UnLoading Slope for Luding and Edinburgh Linear -----------------*/ createUnloadingStiffness(PropertyRegistry & registry,const char * caller,bool)991 MatrixProperty* createUnloadingStiffness(PropertyRegistry & registry, const char * caller, bool) 992 { 993 return createPerTypePairProperty(registry, UNLOADING_STIFFNESS, caller); 994 } 995 996 /* -------------------------- Loading Slope for Luding and Edinburgh Linear -------------------- */ createLoadingStiffness(PropertyRegistry & registry,const char * caller,bool)997 MatrixProperty* createLoadingStiffness(PropertyRegistry & registry, const char * caller, bool) 998 { 999 return createPerTypePairProperty(registry, LOADING_STIFFNESS, caller); 1000 } 1001 createMaxCohesiveStress(PropertyRegistry & registry,const char * caller,bool)1002 MatrixProperty* createMaxCohesiveStress(PropertyRegistry & registry, const char * caller, bool) 1003 { 1004 return createPerTypePairProperty(registry, MAX_COHESIVE_STRESS, caller); 1005 } 1006 1007 /* ---------------------------------------------------------------------- */ 1008 createCohesionStrength(PropertyRegistry & registry,const char * caller,bool)1009 MatrixProperty* createCohesionStrength(PropertyRegistry & registry, const char * caller, bool) 1010 { 1011 return createPerTypePairProperty(registry, COHESION_STRENGTH, caller); 1012 } 1013 1014 /* ---------------------------------------------------------------------- */ 1015 1016 // Dissipation matrix creation (fiber & bond model) (1/timeScale) 1017 createDissipationMatrix(PropertyRegistry & registry,const char * name,const char * caller,bool sanity_checks)1018 MatrixProperty* createDissipationMatrix(PropertyRegistry & registry, const char * name, const char * caller, bool sanity_checks) 1019 { 1020 LAMMPS * lmp = registry.getLAMMPS(); 1021 const int max_type = registry.max_type(); 1022 1023 MatrixProperty * matrix = new MatrixProperty(max_type+1, max_type+1); 1024 FixPropertyGlobal * property = registry.getGlobalProperty(name,"property/global","peratomtypepair",max_type,max_type,caller); 1025 1026 for(int i=1;i< max_type+1; i++) 1027 { 1028 for(int j=1;j<max_type+1;j++) 1029 { 1030 const double gamma_one = property->compute_array(i-1,j-1); 1031 1032 if(sanity_checks) 1033 { 1034 if(gamma_one < lmp->update->dt) 1035 { 1036 char buffer [100]; 1037 sprintf(buffer,"%s > time-step size required",name); 1038 lmp->error->all(FLERR,buffer); 1039 } 1040 } 1041 1042 matrix->data[i][j] = 1./gamma_one; //save already inverted parameter 1043 } 1044 } 1045 1046 return matrix; 1047 } 1048 1049 /* ---------------------------------------------------------------------- */ 1050 1051 // Attrition properties createAbrasiveWearSeverity(PropertyRegistry & registry,const char * caller,bool sanity_checks)1052 MatrixProperty* createAbrasiveWearSeverity(PropertyRegistry & registry, const char * caller, bool sanity_checks) 1053 { 1054 MatrixProperty* wearSeverityMatrix = MODEL_PARAMS::createPerTypePairProperty(registry, "abrasiveWearSeverity", caller); 1055 return wearSeverityMatrix; 1056 } 1057 createErosiveWearSeverity(PropertyRegistry & registry,const char * caller,bool sanity_checks)1058 MatrixProperty* createErosiveWearSeverity(PropertyRegistry & registry, const char * caller, bool sanity_checks) 1059 { 1060 MatrixProperty* wearSeverityMatrix = MODEL_PARAMS::createPerTypePairProperty(registry, "erosiveWearSeverity", caller); 1061 return wearSeverityMatrix; 1062 } 1063 createVelocityLimit(PropertyRegistry & registry,const char * caller,bool sanity_checks)1064 MatrixProperty* createVelocityLimit(PropertyRegistry & registry, const char * caller, bool sanity_checks) 1065 { 1066 MatrixProperty* velocityLimitMatrix = MODEL_PARAMS::createPerTypePairProperty(registry, "velocityLimit", caller); 1067 return velocityLimitMatrix; 1068 } 1069 createForceLimit(PropertyRegistry & registry,const char * caller,bool sanity_checks)1070 MatrixProperty* createForceLimit(PropertyRegistry & registry, const char * caller, bool sanity_checks) 1071 { 1072 MatrixProperty* forceLimitMatrix = MODEL_PARAMS::createPerTypePairProperty(registry, "forceLimit", caller); 1073 return forceLimitMatrix; 1074 } 1075 createAttritionLimit(PropertyRegistry & registry,const char * caller,bool sanity_checks)1076 VectorProperty* createAttritionLimit(PropertyRegistry & registry, const char * caller, bool sanity_checks) 1077 { 1078 VectorProperty* attritionLimit = MODEL_PARAMS::createPerTypeProperty(registry, "attritionLimit", caller); 1079 LAMMPS *lmp = registry.getLAMMPS(); 1080 for (int i = 1; i < registry.max_type()+1; i++) 1081 { 1082 const double limit = attritionLimit->data[i]; 1083 if (limit < 0 || limit > 1) 1084 lmp->error->all(FLERR, "Attrition limit needs to be between 0 and 1"); 1085 } 1086 return attritionLimit; 1087 } 1088 1089 /* ---------------------------------------------------------------------- */ 1090 } 1091