Soft-Bond Model Implementation
See this page for the documentation of this contact model.
contactmodelsoftbond.h
1#pragma once
2// contactmodelsoftbond.h
3
4#include "contactmodel/src/contactmodelmechanical.h"
5
6#ifdef SOFTBOND_LIB
7# define SOFTBOND_EXPORT EXPORT_TAG
8#elif defined(NO_MODEL_IMPORT)
9# define SOFTBOND_EXPORT
10#else
11# define SOFTBOND_EXPORT IMPORT_TAG
12#endif
13
14namespace cmodelsxd {
15 using namespace itasca;
16
17 class ContactModelSoftBond : public ContactModelMechanical {
18 public:
19 // Constructor: Set default values for contact model properties.
20 SOFTBOND_EXPORT ContactModelSoftBond();
21 SOFTBOND_EXPORT ContactModelSoftBond(const ContactModelSoftBond &) noexcept;
22 SOFTBOND_EXPORT const ContactModelSoftBond & operator=(const ContactModelSoftBond &);
23 SOFTBOND_EXPORT void addToStorage(poly::vector<ContactModelMechanical> *s,int ww = -1) override;
24 // Destructor, called when contact is deleted: free allocated memory, etc.
25 SOFTBOND_EXPORT virtual ~ContactModelSoftBond();
26 // Contact model name (used as keyword for commands and FISH).
27 string getName() const override { return "softbond"; }
28 // The index provides a quick way to determine the type of contact model.
29 // Each type of contact model in PFC must have a unique index; this is assigned
30 // by PFC when the contact model is loaded. This index should be set to -1
31 void setIndex(int i) override { index_=i;}
32 int getIndex() const override {return index_;}
33 // Contact model version number (e.g., MyModel05_1). The version number can be
34 // accessed during the save-restore operation (within the archive method,
35 // testing {stream.getRestoreVersion() == getMinorVersion()} to allow for
36 // future modifications to the contact model data structure.
37 uint32 getMinorVersion() const override;
38 // Copy the state information to a newly created contact model.
39 // Provide access to state information, for use by copy method.
40 void copy(const ContactModel *c) override;
41 // Provide save-restore capability for the state information.
42 void archive(ArchiveStream &) override;
43 // Enumerator for the properties.
44 enum PropertyKeys {
45 kwKn=1
46 , kwKs
47 , kwFric
48 , kwBMul
49 , kwTMul
50 , kwSBMode
51 , kwSBF
52 , kwSBM
53 , kwSBS
54 , kwSBBS
55 , kwSBTS
56 , kwSBRMul
57 , kwSBRadius
58 , kwEmod
59 , kwKRatio
60 , kwDpNRatio
61 , kwDpSRatio
62 , kwDpMode
63 , kwDpF
64 , kwSBState
65 , kwSBTStr
66 , kwSBSStr
67 , kwSBCoh
68 , kwSBFa
69 , kwSBMCF
70 , kwSBSig
71 , kwSBTau
72 , kwSBSoft
73 , kwSBCut
74 , kwSBArea
75 , kwUserArea
76 , kwRGap
77 };
78 // Contact model property names in a comma separated list. The order corresponds with
79 // the order of the PropertyKeys enumerator above. One can visualize any of these
80 // properties in PFC automatically.
81 string getProperties() const override {
82 return "kn"
83 ",ks"
84 ",fric"
85 ",sb_bmul"
86 ",sb_tmul"
87 ",sb_mode"
88 ",sb_force"
89 ",sb_moment"
90 ",sb_slip"
91 ",sb_slipb"
92 ",sb_slipt"
93 ",sb_rmul"
94 ",sb_radius"
95 ",emod"
96 ",kratio"
97 ",dp_nratio"
98 ",dp_sratio"
99 ",dp_mode"
100 ",dp_force"
101 ",sb_state"
102 ",sb_ten"
103 ",sb_shear"
104 ",sb_coh"
105 ",sb_fa"
106 ",sb_mcf"
107 ",sb_sigma"
108 ",sb_tau"
109 ",sb_soft"
110 ",sb_cut"
111 ",sb_area"
112 ",user_area"
113 ",rgap"
114 ;
115 }
116 // Enumerator for the energies.
117 enum EnergyKeys {
118 kwEStrain=1
119 , kwESlip
120 , kwEDashpot
121 };
122 // Contact model energy names in a comma separated list. The order corresponds with
123 // the order of the EnergyKeys enumerator above.
124 string getEnergies() const override {
125 return "energy-strain"
126 ",energy-slip"
127 ",energy-dashpot";
128 }
129 // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
130 double getEnergy(uint32 i) const override;
131 // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1)
132 // returns wther or not the estrain energy is accumulated which is false).
133 bool getEnergyAccumulate(uint32 i) const override;
134 // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
135 void setEnergy(uint32 i,const double &d) override; // Base 1
136 // Activate the energy. This is only called if the energy tracking is enabled.
137 void activateEnergy() override { if (energies_) return; energies_ = NEW Energies();}
138 // Returns whether or not the energy tracking has been enabled for this contact.
139 bool getEnergyActivated() const override {return (energies_ != 0);}
140
141 // Enumerator for contact model related FISH callback events.
142 enum FishCallEvents {
143 fActivated=0
144 ,fSlipChange
145 ,fBondBreak
146 };
147 // Contact model FISH callback event names in a comma separated list. The order corresponds with
148 // the order of the FishCallEvents enumerator above.
149 string getFishCallEvents() const override {
150 return
151 "contact_activated"
152 ",slip_change"
153 ",bond_break";
154 }
155
156 // Return the specified contact model property.
157 base::Property getProperty(uint32 i,const IContact *) const override;
158 // The return value denotes whether or not the property corresponds to the global
159 // or local coordinate system (TRUE: global system, FALSE: local system). The
160 // local system is the contact-plane system (nst) defined as follows.
161 // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
162 // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
163 // and tc are unit vectors in directions of the nst axes.
164 // This is used when rendering contact model properties that are vectors.
165 bool getPropertyGlobal(uint32 i) const override;
166 // Set the specified contact model property, ensuring that it is of the correct type
167 // and within the correct range --- if not, then throw an exception.
168 // The return value denotes whether or not the update has affected the timestep
169 // computation (by having modified the translational or rotational tangent stiffnesses).
170 // If true is returned, then the timestep will be recomputed.
171 bool setProperty(uint32 i,const base::Property &v,IContact *) override;
172 // The return value denotes whether or not the property is read-only
173 // (TRUE: read-only, FALSE: read-write).
174 bool getPropertyReadOnly(uint32 i) const override;
175
176 // The return value denotes whether or not the property is inheritable
177 // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
178 // the endPropertyUpdated method.
179 bool supportsInheritance(uint32 i) const override;
180 // Return whether or not inheritance is enabled for the specified property.
181 bool getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i); return (inheritanceField_ & mask) ? true : false; }
182 // Set the inheritance flag for the specified property.
183 void setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i); if (b) inheritanceField_ |= mask; else inheritanceField_ &= ~mask; }
184
185 // Enumerator for contact model methods.
186 enum MethodKeys { kwDeformability=1, kwBond, kwUnbond, kwArea};
187 // Contact model methoid names in a comma separated list. The order corresponds with
188 // the order of the MethodKeys enumerator above.
189 string getMethods() const override { return "deformability,bond,unbond,area";}
190 // Return a comma seprated list of the contact model method arguments (base 1).
191 string getMethodArguments(uint32 i) const override;
192 // Set contact model method arguments (base 1).
193 // The return value denotes whether or not the update has affected the timestep
194 // computation (by having modified the translational or rotational tangent stiffnesses).
195 // If true is returned, then the timestep will be recomputed.
196 bool setMethod(uint32 i,const std::vector<base::Property> &vl,IContact *con=0) override;
197
198 // Prepare for entry into ForceDispLaw. The validate function is called when:
199 // (1) the contact is created, (2) a property of the contact that returns a true via
200 // the setProperty method has been modified and (3) when a set of cycles is executed
201 // via the {cycle N} command.
202 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
203 bool validate(ContactModelMechanicalState *state,const double ×tep) override;
204 // The endPropertyUpdated method is called whenever a surface property (with a name
205 // that matches an inheritable contact model property name) of one of the contacting
206 // pieces is modified. This allows the contact model to update its associated
207 // properties. The return value denotes whether or not the update has affected
208 // the time step computation (by having modified the translational or rotational
209 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
210 bool endPropertyUpdated(const string &name,const IContactMechanical *c) override;
211 // The forceDisplacementLaw function is called during each cycle. Given the relative
212 // motion of the two contacting pieces (via
213 // state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
214 // state->relativeAngularIncrement_ (Dtt, Dtbs, Dtbt)
215 // Ddn : relative normal-displacement increment, Ddn > 0 is opening
216 // Ddss : relative shear-displacement increment (s-axis component)
217 // Ddst : relative shear-displacement increment (t-axis component)
218 // Dtt : relative twist-rotation increment
219 // Dtbs : relative bend-rotation increment (s-axis component)
220 // Dtbt : relative bend-rotation increment (t-axis component)
221 // The relative displacement and rotation increments:
222 // Dd = Ddn*nc + Ddss*sc + Ddst*tc
223 // Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
224 // where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
225 // [see {Table 1: Contact State Variables} in PFC Model Components:
226 // Contacts and Contact Models: Contact Resolution]
227 // ) and the contact properties, this function must update the contact force and
228 // moment.
229 // The force_ is acting on piece 2, and is expressed in the local coordinate system
230 // (defined in getPropertyGlobal) such that the first component positive denotes
231 // compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
232 // in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to
233 // get the total moment.
234 // The return value indicates the contact activity status (TRUE: active, FALSE:
235 // inactive) during the next cycle.
236 // Additional information:
237 // * If state->activated() is true, then the contact has just become active (it was
238 // inactive during the previous time step).
239 // * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
240 // the forceDispLaw handle the case of {state->canFail_ == true}.
241 bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) override;
242 bool thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
243 // The getEffectiveXStiffness functions return the translational and rotational
244 // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
245 // the translational tangent shear stiffness is zero (but this stiffness reduction
246 // is typically ignored when computing a stable time step). If the contact model
247 // includes a dashpot, then the translational stiffnesses must be increased (see
248 // Potyondy (2009)).
249 // [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
250 // Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
251 // December 7, 2009.]
252 DVect2 getEffectiveTranslationalStiffness() const override { return effectiveTranslationalStiffness_; }
253 DAVect getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness_;}
254
255 // Return a new instance of the contact model. This is used in the CMAT
256 // when a new contact is created.
257 ContactModelSoftBond *clone() const override { return NEW ContactModelSoftBond(); }
258 // The getActivityDistance function is called by the contact-resolution logic when
259 // the CMAT is modified. Return value is the activity distance used by the
260 // checkActivity function.
261 double getActivityDistance() const override {return rgap_;}
262 // The isOKToDelete function is called by the contact-resolution logic when...
263 // Return value indicates whether or not the contact may be deleted.
264 // If TRUE, then the contact may be deleted when it is inactive.
265 // If FALSE, then the contact may not be deleted (under any condition).
266 bool isOKToDelete() const override { return !isBonded(); }
267 // Zero the forces and moments stored in the contact model. This function is called
268 // when the contact becomes inactive.
269 void resetForcesAndMoments() override {
270 sb_F(DVect(0.0));
271 dp_F(DVect(0.0));
272 sb_M(DAVect(0.0));
273 if (energies_) {
274 energies_->estrain_ = 0.0;
275 }
276 }
277 void setForce(const DVect &v,IContact *c) override;
278 void setArea(const double &d) override { userArea_ = d; }
279 double getArea() const override { return userArea_; }
280
281 // The checkActivity function is called by the contact-resolution logic when...
282 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
283 bool checkActivity(const double &gap) override { return gap <= rgap_ || isBonded();}
284
285 // Returns the sliding state (FALSE is returned if not implemented).
286 bool isSliding() const override { return sb_S_; }
287 // Returns the bonding state (FALSE is returned if not implemented).
288 bool isBonded() const override { return bProps_ ? (bProps_->sb_state_ >= 3) : false; }
289 void unbond() override { if (bProps_) bProps_->sb_state_ = 0; }
290
291 // Both of these methods are called only for contacts with facets where the wall
292 // resolution scheme is set the full. In such cases one might wish to propagate
293 // contact state information (e.g., shear force) from one active contact to another.
294 // See the Faceted Wall section in the documentation.
295 void propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
296 void setNonForcePropsFrom(IContactModel *oldCM) override;
297 /// Return the total force that the contact model holds.
298 DVect getForce() const override;
299 /// Return the total moment on 1 that the contact model holds
300 DAVect getMomentOn1(const IContactMechanical *) const override;
301 /// Return the total moment on 1 that the contact model holds
302 DAVect getMomentOn2(const IContactMechanical *) const override;
303
304 /// Return moments without torque
305 DAVect getModelMomentOn1() const override;
306 DAVect getModelMomentOn2() const override;
307
308 // Used to efficiently get properties from the contact model for the object pane.
309 // List of properties for the object pane, comma separated.
310 // All properties will be cast to doubles for comparison. No other comparisons
311 // are supported. This may not be the same as the entire property list.
312 // Return property name and type for plotting.
313 void objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *) const override;
314 // All properties cast to doubles - this is what can be compared.
315 void objectPropValues(std::vector<double> *,const IContact *) const override;
316
317 // Methods to get and set properties.
318 const double & kn() const {return kn_;}
319 void kn(const double &d) {kn_=d;}
320 const double & ks() const {return ks_;}
321 void ks(const double &d) {ks_=d;}
322 const double & fric() const {return fric_;}
323 void fric(const double &d) {fric_=d;}
324 const double & sb_bmul() const { return sb_bmul_; }
325 void sb_bmul(const double &d) { sb_bmul_ = d; }
326 const double & sb_tmul() const { return sb_tmul_; }
327 void sb_tmul(const double &d) { sb_tmul_ = d; }
328 const DVect & sb_F() const {return sb_F_;}
329 void sb_F(const DVect &f) { sb_F_=f;}
330 const DAVect & sb_M() const { return sb_M_; }
331 void sb_M(const DAVect &f) { sb_M_ = f; }
332 bool sb_S() const {return sb_S_;}
333 void sb_S(bool b) { sb_S_=b;}
334 bool sb_BS() const { return sb_BS_; }
335 void sb_BS(bool b) { sb_BS_ = b; }
336 bool sb_TS() const { return sb_TS_; }
337 void sb_TS(bool b) { sb_TS_ = b; }
338 const double & sb_rmul() const { return sb_rmul_; }
339 void sb_rmul(const double &d) { sb_rmul_ = d; }
340 uint32 sb_mode() const {return sb_mode_;}
341 void sb_mode(uint32 i) { sb_mode_=i;}
342
343 bool hasBond() const { return bProps_ ? true : false; }
344 int sb_state() const { return (hasBond() ? bProps_->sb_state_ : 0); }
345 void sb_state(int i) { if (!hasBond()) return; bProps_->sb_state_ = i; }
346 double sb_Ten() const { return (hasBond() ? (bProps_->sb_ten_) : 0.0); }
347 void sb_Ten(const double &d) { if (!hasBond()) return; bProps_->sb_ten_ = d; }
348 double sb_Coh() const { return (hasBond() ? (bProps_->sb_coh_) : 0.0); }
349 void sb_Coh(const double &d) { if (!hasBond()) return; bProps_->sb_coh_ = d; }
350 double sb_FA() const { return (hasBond() ? (bProps_->sb_fa_) : 0.0); }
351 void sb_FA(const double &d) { if (!hasBond()) return; bProps_->sb_fa_ = d; }
352 double sb_MCF() const {return (hasBond() ? (bProps_->sb_mcf_) : 0.0);}
353 void sb_MCF(const double &d) { if(!hasBond()) return; bProps_->sb_mcf_=d;}
354 double sb_soft() const { return (hasBond() ? (bProps_->sb_soft_) : 0.0); }
355 void sb_soft(const double &d) { if (!hasBond()) return; bProps_->sb_soft_ = d; }
356 double sb_cut() const { return (hasBond() ? (bProps_->sb_cut_) : 0.0); }
357 void sb_cut(const double &d) { if (!hasBond()) return; bProps_->sb_cut_ = d; }
358 double sb_maxTen() const { return (hasBond() ? (bProps_->sb_maxTen_) : 0.0); }
359 void sb_maxTen(const double &d) { if (!hasBond()) return; bProps_->sb_maxTen_ = d; }
360 double sb_delu() const { return (hasBond() ? (bProps_->sb_delu_) : 0.0); }
361 void sb_delu(const double &d) { if (!hasBond()) return; bProps_->sb_delu_ = d; }
362 Quat sb_delo() const { return (hasBond() ? (bProps_->sb_delo_) : Quat::identity()); }
363 void sb_delo(const Quat &d) { if (!hasBond()) return; bProps_->sb_delo_ = d; }
364 double sb_maxu() const { return (hasBond() ? (bProps_->sb_maxu_) : 0.0); }
365 void sb_maxu(const double &d) { if (!hasBond()) return; bProps_->sb_maxu_ = d; }
366 double sb_critu() const { return (hasBond() ? (bProps_->sb_critu_) : 0.0); }
367 void sb_critu(const double &d) { if (!hasBond()) return; bProps_->sb_critu_ = d; }
368
369 bool hasDamping() const {return dpProps_ ? true : false;}
370 double dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
371 void dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
372 double dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
373 void dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
374 int dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
375 void dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
376 DVect dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
377 void dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
378
379 bool hasEnergies() const {return energies_ ? true:false;}
380 double estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
381 void estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
382 double eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
383 void eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
384 double edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
385 void edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
386
387 uint32 inheritanceField() const {return inheritanceField_;}
388 void inheritanceField(uint32 i) {inheritanceField_ = i;}
389
390 const DVect2 & effectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
391 void effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
392 const DAVect & effectiveRotationalStiffness() const {return effectiveRotationalStiffness_;}
393 void effectiveRotationalStiffness(const DAVect &v ) {effectiveRotationalStiffness_=v;}
394
395 private:
396 // Index - used internally by PFC. Should be set to -1 in the cpp file.
397 static int index_;
398
399 bool FDLawBonded(ContactModelMechanicalState *state, const double ×tep);
400 bool FDLawUnBonded(ContactModelMechanicalState *state, const double ×tep);
401
402 // Structure to compute stiffness
403 struct StiffData {
404 DVect2 trans_ = DVect2(0.0);
405 DAVect ang_ = DAVect(0.0);
406 double reff_ = 0.0;
407 };
408
409 // Structure to store the energies.
410 struct Energies {
411 double estrain_ = 0.0; // elastic energy
412 double eslip_ = 0.0; // work dissipated by friction
413 double edashpot_ = 0.0; // work dissipated by dashpots
414 };
415
416 // Structure to store dashpot quantities.
417 struct dpProps {
418 double dp_nratio_ = 0.0; // normal viscous critical damping ratio
419 double dp_sratio_ = 0.0; // shear viscous critical damping ratio
420 int dp_mode_ = 0; // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
421 DVect dp_F_ = DVect(0.0); // Force in the dashpots
422 };
423
424 // Structure to store bond-related quantities.
425 struct bProps {
426 int sb_state_ = 0; // bond mode - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (B), 4 (B-Softening), 5 (B-Compression from Softening)
427 double sb_ten_ = 0.0; // normal strength
428 double sb_coh_ = 0.0; // cohesion
429 double sb_fa_ = 0.0; // friction angle
430 double sb_mcf_ = 1.0; // moment contribution factor
431 double sb_soft_ = 0.0; // softening factor
432 double sb_cut_ = 1.0; // critical bond length
433 double sb_maxTen_ = 0.0; // tensile strength one needs to reach for softening
434 double sb_delu_ = 0.0; // incremental elongation in softening
435 Quat sb_delo_ = Quat::identity(); // incremental orientation in softening
436 double sb_maxu_ = 0.0; // max elongation for softening
437 double sb_critu_ = 0.0; // critical elongation for softening
438 };
439
440
441 bool updateKn(const IContactMechanical *con);
442 bool updateKs(const IContactMechanical *con);
443 bool updateFric(const IContactMechanical *con);
444
445 StiffData computeStiffData(ContactModelMechanicalState *state) const;
446 DVect3 computeGeomData(const DVect2 &end1C,const DVect2 &end2C) const;
447 DVect2 SMax(const DVect2 &end1C,const DVect2 &end2C) const; // Maximum stress (tensile,shear) at bond periphery
448 double shearStrength(const double &pbArea) const; // Bond shear strength
449 double strainEnergy(double kn, double ks, double kb, double kt) const;
450
451 void updateStiffness(ContactModelMechanicalState *state);
452
453 // Contact model inheritance fields.
454 uint32 inheritanceField_;
455
456 // Effective translational stiffness.
457 DVect2 effectiveTranslationalStiffness_ = DVect2(0.0);
458 DAVect effectiveRotationalStiffness_ = DAVect(0.0); // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)
459
460 // linear model properties
461 double kn_ = 0.0; // Normal stiffness
462 double ks_ = 0.0; // Shear stiffness
463 double fric_ = 0.0; // Coulomb friction coefficient
464 double sb_bmul_ = 1.0; // Bending friction multiplier
465 double sb_tmul_ = 1.0; // Twisting friction multiplier
466 uint32 sb_mode_ = 0; // specifies absolute (0) or incremental (1) behavior for the the normal force
467 DVect sb_F_ = DVect(0); // Force carried in the model
468 DAVect sb_M_ = DAVect(0); // moment (bending + twisting in 3D)
469 bool sb_S_ = false; // The current slip state
470 bool sb_BS_ = false; // The bending slip state
471 bool sb_TS_ = false; // The twisting slip state
472 double sb_rmul_ = 1.0; // Radius multiplier
473 double userArea_ = 0.0; // Area as specified by the user
474 double rgap_ = 0.0; // Reference gap
475
476 dpProps * dpProps_ = nullptr; // The viscous properties
477 bProps * bProps_ = nullptr; // The bond properties
478
479 Energies * energies_ = nullptr; // The energies
480
481 };
482} // namespace cmodelsxd
483// EoF
contactmodelsoftbond.cpp
1// contactmodelsoftbond.cpp
2#include "contactmodelsoftbond.h"
3
4#include "module/interface/icontactmechanical.h"
5#include "module/interface/icontact.h"
6#include "module/interface/ipiece.h"
7#include "module/interface/ifishcalllist.h"
8
9#include "utility/src/tptr.h"
10#include "shared/src/mathutil.h"
11
12#include "kernel/interface/iprogram.h"
13#include "module/interface/icontactthermal.h"
14#include "contactmodel/src/contactmodelthermal.h"
15#include "contactmodel/src/contactmodelfluid.h"
16#include "../version.txt"
17#include "fish/src/parameter.h"
18
19#ifdef SOFTBOND_LIB
20#ifdef _WIN32
21 int __stdcall DllMain(void *,unsigned, void *) {
22 return 1;
23 }
24#endif
25 extern "C" EXPORT_TAG const char *getName() {
26#if DIM==3
27 return "contactmodelmechanical3dsoftbond";
28#else
29 return "contactmodelmechanical2dsoftbond";
30#endif
31 }
32
33 extern "C" EXPORT_TAG unsigned getMajorVersion() {
34 return MAJOR_VERSION;
35 }
36
37 extern "C" EXPORT_TAG unsigned getMinorVersion() {
38 return MINOR_VERSION;
39 }
40
41 extern "C" EXPORT_TAG void *createInstance() {
42 cmodelsxd::ContactModelSoftBond *m = NEW cmodelsxd::ContactModelSoftBond();
43 return (void *)m;
44 }
45#endif
46
47namespace cmodelsxd {
48 static const uint32 KnMask = 0x00000002; // Base 1!
49 static const uint32 KsMask = 0x00000004;
50 static const uint32 FricMask = 0x00000008;
51
52 using namespace itasca;
53
54 int ContactModelSoftBond::index_ = -1;
55 uint32 ContactModelSoftBond::getMinorVersion() const { return MINOR_VERSION; }
56
57 ContactModelSoftBond::ContactModelSoftBond() : inheritanceField_(KnMask|KsMask|FricMask) {
58 }
59
60 ContactModelSoftBond::ContactModelSoftBond(const ContactModelSoftBond &m) noexcept {
61 inheritanceField(KnMask|KsMask|FricMask);
62 this->copy(&m);
63 }
64
65 const ContactModelSoftBond & ContactModelSoftBond::operator=(const ContactModelSoftBond &m) {
66 inheritanceField(KnMask|KsMask|FricMask);
67 this->copy(&m);
68 return *this;
69 }
70
71 void ContactModelSoftBond::addToStorage(poly::vector<ContactModelMechanical> *s,int ww) {
72 s->addToStorage<ContactModelSoftBond>(*this,ww);
73 }
74
75 ContactModelSoftBond::~ContactModelSoftBond() {
76 // Make sure to clean up after yourself!
77 if (dpProps_)
78 delete dpProps_;
79 if (bProps_)
80 delete bProps_;
81 if (energies_)
82 delete energies_;
83 }
84
85 void ContactModelSoftBond::archive(ArchiveStream &stream) {
86 // The stream allows one to archive the values of the contact model
87 // so that it can be saved and restored. The minor version can be
88 // used here to allow for incremental changes to the contact model too.
89 stream & kn_;
90 stream & ks_;
91 stream & fric_;
92 stream & sb_mode_;
93 stream & sb_F_;
94 stream & sb_M_;
95 stream & sb_S_;
96 stream & sb_BS_;
97 stream & sb_TS_;
98 stream & sb_rmul_;
99
100 if (stream.getArchiveState()==ArchiveStream::Save) {
101 bool b = false;
102 if (dpProps_) {
103 b = true;
104 stream & b;
105 stream & dpProps_->dp_nratio_;
106 stream & dpProps_->dp_sratio_;
107 stream & dpProps_->dp_mode_;
108 stream & dpProps_->dp_F_;
109 }
110 else
111 stream & b;
112
113 b = false;
114 if (energies_) {
115 b = true;
116 stream & b;
117 stream & energies_->estrain_;
118 stream & energies_->eslip_;
119 stream & energies_->edashpot_;
120 }
121 else
122 stream & b;
123
124 b = false;
125 if (bProps_) {
126 b = true;
127 stream & b;
128 stream & bProps_->sb_state_;
129 stream & bProps_->sb_ten_;
130 stream & bProps_->sb_coh_;
131 stream & bProps_->sb_fa_;
132 stream & bProps_->sb_mcf_;
133 stream & bProps_->sb_soft_;
134 stream & bProps_->sb_cut_;
135 stream & bProps_->sb_maxTen_;
136 stream & bProps_->sb_delu_;
137 stream & bProps_->sb_delo_;
138 stream & bProps_->sb_maxu_;
139 stream & bProps_->sb_critu_;
140 }
141 else
142 stream & b;
143
144 } else {
145 bool b(false);
146 stream & b;
147 if (b) {
148 if (!dpProps_)
149 dpProps_ = NEW dpProps();
150 stream & dpProps_->dp_nratio_;
151 stream & dpProps_->dp_sratio_;
152 stream & dpProps_->dp_mode_;
153 stream & dpProps_->dp_F_;
154 }
155 stream & b;
156 if (b) {
157 if (!energies_)
158 energies_ = NEW Energies();
159 stream & energies_->estrain_;
160 stream & energies_->eslip_;
161 stream & energies_->edashpot_;
162 }
163 stream & b;
164 if (b) {
165 if (!bProps_)
166 bProps_ = NEW bProps();
167 stream & bProps_->sb_state_;
168 stream & bProps_->sb_ten_;
169 stream & bProps_->sb_coh_;
170 stream & bProps_->sb_fa_;
171 stream & bProps_->sb_mcf_;
172 stream & bProps_->sb_soft_;
173 stream & bProps_->sb_cut_;
174 stream & bProps_->sb_maxTen_;
175 stream & bProps_->sb_delu_;
176 stream & bProps_->sb_delo_;
177 stream & bProps_->sb_maxu_;
178 stream & bProps_->sb_critu_;
179 }
180
181 }
182
183 stream & inheritanceField_;
184 stream & effectiveTranslationalStiffness_;
185 stream & effectiveRotationalStiffness_;
186
187 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 1) {
188 stream & sb_bmul_;
189 stream & sb_tmul_;
190 }
191
192 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 2)
193 stream & userArea_;
194
195 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 3)
196 stream & rgap_;
197
198 }
199
200 void ContactModelSoftBond::copy(const ContactModel *cm) {
201 // Copy all of the contact model properties. Used in the CMAT
202 // when a new contact is created.
203 ContactModelMechanical::copy(cm);
204 const ContactModelSoftBond *in = dynamic_cast<const ContactModelSoftBond*>(cm);
205 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
206 kn(in->kn());
207 ks(in->ks());
208 fric(in->fric());
209 sb_bmul(in->sb_bmul());
210 sb_tmul(in->sb_tmul());
211 sb_mode(in->sb_mode());
212 sb_F(in->sb_F());
213 sb_S(in->sb_S());
214 sb_BS(in->sb_BS());
215 sb_TS(in->sb_TS());
216 sb_rmul(in->sb_rmul());
217 sb_M(in->sb_M());
218
219 if (in->hasDamping()) {
220 if (!dpProps_)
221 dpProps_ = NEW dpProps();
222 dp_nratio(in->dp_nratio());
223 dp_sratio(in->dp_sratio());
224 dp_mode(in->dp_mode());
225 dp_F(in->dp_F());
226 }
227 if (in->hasEnergies()) {
228 if (!energies_)
229 energies_ = NEW Energies();
230 estrain(in->estrain());
231 eslip(in->eslip());
232 edashpot(in->edashpot());
233 }
234 if (in->hasBond()) {
235 if (!bProps_)
236 bProps_ = NEW bProps();
237 sb_state(in->sb_state());
238 sb_Ten(in->sb_Ten());
239 sb_Coh(in->sb_Coh());
240 sb_FA(in->sb_FA());
241 sb_MCF(in->sb_MCF());
242 sb_soft(in->sb_soft());
243 sb_cut(in->sb_cut());
244 sb_maxTen(in->sb_maxTen());
245 sb_delu(in->sb_delu());
246 sb_delo(in->sb_delo());
247 sb_maxu(in->sb_maxu());
248 sb_critu(in->sb_critu());
249 }
250 userArea_ = in->userArea_;
251 rgap_ = in->rgap_;
252 inheritanceField(in->inheritanceField());
253 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
254 effectiveRotationalStiffness(in->effectiveRotationalStiffness());
255 }
256
257
258 base::Property ContactModelSoftBond::getProperty(uint32 i,const IContact *con) const {
259 // Return the property. The IContact pointer is provided so that
260 // more complicated properties, depending on contact characteristics,
261 // can be calcualted.
262 // The IContact pointer may be a nullptr!
263 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
264 base::Property var;
265 switch (i) {
266 case kwKn: return kn_;
267 case kwKs: return ks_;
268 case kwFric: return fric_;
269 case kwBMul: return sb_bmul_;
270 case kwTMul: return sb_tmul_;
271 case kwSBMode: return sb_mode_;
272 case kwSBF: var.setValue(sb_F_); return var;
273 case kwSBM: var.setValue(sb_M_); return var;
274 case kwSBS: return sb_S_;
275 case kwSBBS: return sb_BS_;
276 case kwSBTS: return sb_TS_;
277 case kwSBRMul: return sb_rmul_;
278 case kwSBRadius: {
279 if (!c) return 0.0;
280 double Cmax1 = c->getEnd1Curvature().y();
281 double Cmax2 = c->getEnd2Curvature().y();
282 if (!userArea_)
283 return sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
284 else {
285#ifdef THREED
286 double rad = std::sqrt(userArea_ / dPi);
287#else
288 double rad = userArea_ / 2.0;
289#endif
290 return rad;
291 }
292
293 }
294 case kwEmod: {
295 if (!c) return 0.0;
296 double rsum = calcRSum(c);
297 if (userArea_)
298#ifdef THREED
299 rsum = 2.0 * std::sqrt(userArea_ / dPi);
300#else
301 rsum = userArea_;
302#endif
303 return kn_ * rsum;
304 }
305 case kwKRatio: return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
306 case kwDpNRatio: return dpProps_ ? dpProps_->dp_nratio_ : 0;
307 case kwDpSRatio: return dpProps_ ? dpProps_->dp_sratio_ : 0;
308 case kwDpMode: return dpProps_ ? dpProps_->dp_mode_ : 0;
309 case kwDpF: {
310 dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
311 return var;
312 }
313 case kwSBState: return bProps_ ? bProps_->sb_state_ : 0;
314 case kwSBTStr: return bProps_ ? bProps_->sb_ten_ : 0.0;
315 case kwSBSStr: {
316 if (!bProps_ or not con) return 0.0;
317 double area = computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x();
318 return shearStrength(area);
319 }
320 case kwSBCoh: return bProps_ ? bProps_->sb_coh_ : 0;
321 case kwSBFa: return bProps_ ? bProps_->sb_fa_ : 0;
322 case kwSBMCF: return bProps_ ? bProps_->sb_mcf_ : 0;
323 case kwSBSig: {
324 if (!bProps_ or bProps_->sb_state_ < 3 or not con) return 0.0;
325 return SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x();
326 }
327 case kwSBTau: {
328 if (!bProps_ or bProps_->sb_state_ < 3 or not con) return 0.0;
329 return SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).y();
330 }
331 case kwSBSoft:
332 if (!bProps_) return 0.0;
333 return bProps_->sb_soft_;
334 case kwSBCut:
335 if (!bProps_) return 0.0;
336 return bProps_->sb_cut_;
337 case kwSBArea: {
338 if (userArea_) return userArea_;
339 if (!con)
340 return 0.0;
341 return computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x();
342 }
343 case kwUserArea:
344 return userArea_;
345 case kwRGap:
346 return rgap_;
347 }
348 assert(0);
349 return base::Property();
350 }
351
352 bool ContactModelSoftBond::getPropertyGlobal(uint32 i) const {
353 // Returns whether or not a property is held in the global axis system (TRUE)
354 // or the local system (FALSE). Used by the plotting logic.
355 switch (i) {
356 case kwSBF:
357 case kwSBM:
358 case kwDpF:
359 return false;
360 }
361 return true;
362 }
363
364 bool ContactModelSoftBond::setProperty(uint32 i,const base::Property &v,IContact *) {
365 // Set a contact model property. Return value indicates that the timestep
366 // should be recalculated.
367 dpProps dp;
368 switch (i) {
369 case kwKn: {
370 if (!v.canConvert<double>())
371 throw Exception("kn must be a double.");
372 double val(v.toDouble());
373 if (val<0.0)
374 throw Exception("Negative kn not allowed.");
375 kn_ = val;
376 return true;
377 }
378 case kwKs: {
379 if (!v.canConvert<double>())
380 throw Exception("ks must be a double.");
381 double val(v.toDouble());
382 if (val<0.0)
383 throw Exception("Negative ks not allowed.");
384 ks_ = val;
385 return true;
386 }
387 case kwFric: {
388 if (!v.canConvert<double>())
389 throw Exception("fric must be a double.");
390 double val(v.toDouble());
391 if (val<0.0)
392 throw Exception("Negative fric not allowed.");
393 fric_ = val;
394 return false;
395 }
396 case kwBMul: {
397 if (!v.canConvert<double>())
398 throw Exception("sb_bmul must be a double.");
399 double val(v.toDouble());
400 if (val<0.0)
401 throw Exception("Negative sb_bmul not allowed.");
402 sb_bmul_ = val;
403 return false;
404 }
405 case kwTMul: {
406 if (!v.canConvert<double>())
407 throw Exception("sb_tmul must be a double.");
408 double val(v.toDouble());
409 if (val<0.0)
410 throw Exception("Negative st_bmul not allowed.");
411 sb_tmul_ = val;
412 return false;
413 }
414 case kwSBMode: {
415 if (!v.canConvert<int64>())
416 throw Exception("sb_mode must be 0 (absolute) or 1 (incremental).");
417 double val(v.toUInt());
418 if (val>1)
419 throw Exception("sb_mode must be 0 (absolute) or 1 (incremental).");
420 sb_mode_ = val;
421 return false;
422 }
423 case kwSBRMul: {
424 if (!v.canConvert<double>())
425 throw Exception("rmul must be a double.");
426 double val(v.toDouble());
427 if (val<0.0)
428 throw Exception("Negative rmul not allowed.");
429 sb_rmul_ = val;
430 return false;
431 }
432 case kwSBF: {
433 if (!v.canConvert<DVect>())
434 throw Exception("sb_force must be a vector.");
435 DVect val(v.value<DVect>());
436 sb_F_ = val;
437 return false;
438 }
439 case kwSBM: {
440 DAVect val(0.0);
441#ifdef TWOD
442 if (!v.canConvert<DAVect>() && !v.canConvert<double>())
443 throw Exception("res_moment must be an angular vector.");
444 if (v.canConvert<DAVect>())
445 val = DAVect(v.value<DAVect>());
446 else
447 val = DAVect(v.toDouble());
448#else
449 if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
450 throw Exception("res_moment must be an angular vector.");
451 if (v.canConvert<DAVect>())
452 val = DAVect(v.value<DAVect>());
453 else
454 val = DAVect(v.value<DVect>());
455#endif
456 sb_M_ = val;
457 return false;
458 }
459 case kwDpNRatio: {
460 if (!v.canConvert<double>())
461 throw Exception("dp_nratio must be a double.");
462 double val(v.toDouble());
463 if (val<0.0)
464 throw Exception("Negative dp_nratio not allowed.");
465 if (val == 0.0 && !dpProps_)
466 return false;
467 if (!dpProps_)
468 dpProps_ = NEW dpProps();
469 dpProps_->dp_nratio_ = val;
470 return true;
471 }
472 case kwDpSRatio: {
473 if (!v.canConvert<double>())
474 throw Exception("dp_sratio must be a double.");
475 double val(v.toDouble());
476 if (val<0.0)
477 throw Exception("Negative dp_sratio not allowed.");
478 if (val == 0.0 && !dpProps_)
479 return false;
480 if (!dpProps_)
481 dpProps_ = NEW dpProps();
482 dpProps_->dp_sratio_ = val;
483 return true;
484 }
485 case kwDpMode: {
486 if (!v.canConvert<int64>())
487 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
488 int val(v.toInt());
489 if (val == 0 && !dpProps_)
490 return false;
491 if (val < 0 || val > 3)
492 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
493 if (!dpProps_)
494 dpProps_ = NEW dpProps();
495 dpProps_->dp_mode_ = val;
496 return false;
497 }
498 case kwDpF: {
499 if (!v.canConvert<DVect>())
500 throw Exception("dp_force must be a vector.");
501 DVect val(v.value<DVect>());
502 if (val.fsum() == 0.0 && !dpProps_)
503 return false;
504 if (!dpProps_)
505 dpProps_ = NEW dpProps();
506 dpProps_->dp_F_ = val;
507 return false;
508 }
509 case kwSBTStr: {
510 if (!v.canConvert<double>())
511 throw Exception("sb_ten must be a double.");
512 double val(v.toDouble());
513 if (val < 0.0)
514 throw Exception("Negative sb_ten not allowed.");
515 if (val == 0.0 && !bProps_)
516 return false;
517 if (!bProps_)
518 bProps_ = NEW bProps();
519 bProps_->sb_ten_ = val;
520 return false;
521 }
522 case kwSBCoh: {
523 if (!v.canConvert<double>())
524 throw Exception("sb_coh must be a double.");
525 double val(v.toDouble());
526 if (val<0.0)
527 throw Exception("Negative pb_coh not allowed.");
528 if (val == 0.0 && !bProps_)
529 return false;
530 if (!bProps_)
531 bProps_ = NEW bProps();
532 bProps_->sb_coh_ = val;
533 return false;
534 }
535 case kwSBFa: {
536 if (!v.canConvert<double>())
537 throw Exception("sb_fa must be a double.");
538 double val(v.toDouble());
539 if (val<0.0)
540 throw Exception("Negative sb_fa not allowed.");
541 if (val >= 90.0)
542 throw Exception("sb_fa must be lower than 90.0 degrees.");
543 if (val == 0.0 && !bProps_)
544 return false;
545 if (!bProps_)
546 bProps_ = NEW bProps();
547 bProps_->sb_fa_ = val;
548 return false;
549 }
550 case kwSBMCF: {
551 if (!v.canConvert<double>())
552 throw Exception("sb_mcf must be a double.");
553 double val(v.toDouble());
554 if (val<0.0)
555 throw Exception("Negative sb_mcf not allowed.");
556 if (val > 1.0)
557 throw Exception("sb_mcf must be lower or equal to 1.0.");
558 if (val == 1.0 && !bProps_)
559 return false;
560 if (!bProps_)
561 bProps_ = NEW bProps();
562 bProps_->sb_mcf_ = val;
563 return false;
564 }
565 case kwSBSoft: {
566 if (!v.canConvert<double>())
567 throw Exception("sb_soft must be a double.");
568 double val(v.toDouble());
569 if (val < 0.0)
570 throw Exception("Negative pb_soft not allowed.");
571 if (!bProps_)
572 bProps_ = NEW bProps();
573 bProps_->sb_soft_ = val;
574 return false;
575 }
576 case kwSBCut: {
577 if (!v.canConvert<double>())
578 throw Exception("sb_cut must be a double.");
579 double val(v.toDouble());
580 if (val < 0.0)
581 throw Exception("Negative sb_cut not allowed.");
582 if (!bProps_)
583 bProps_ = NEW bProps();
584 bProps_->sb_cut_ = val;
585 return false;
586 }
587 case kwSBArea:
588 case kwUserArea: {
589 if (!v.canConvert<double>())
590 throw Exception("area must be a double.");
591 double val(v.toDouble());
592 if (val < 0.0)
593 throw Exception("Negative area not allowed.");
594 userArea_ = val;
595 return true;
596 }
597 case kwRGap: {
598 if (!v.canConvert<double>())
599 throw Exception("Reference gap must be a double.");
600 double val(v.toDouble());
601 rgap_ = val;
602 return false;
603 }
604 }
605 return false;
606 }
607
608 bool ContactModelSoftBond::getPropertyReadOnly(uint32 i) const {
609 // Returns TRUE if a property is read only or FALSE otherwise.
610 switch (i) {
611 case kwDpF:
612 case kwSBS:
613 case kwSBBS:
614 case kwSBTS:
615 case kwEmod:
616 case kwKRatio:
617 case kwSBState:
618 case kwSBRadius:
619 case kwSBSStr:
620 case kwSBSig:
621 case kwSBTau:
622 return true;
623 default:
624 break;
625 }
626 return false;
627 }
628
629 bool ContactModelSoftBond::supportsInheritance(uint32 i) const {
630 // Returns TRUE if a property supports inheritance or FALSE otherwise.
631 switch (i) {
632 case kwKn:
633 case kwKs:
634 case kwFric:
635 return true;
636 default:
637 break;
638 }
639 return false;
640 }
641
642 string ContactModelSoftBond::getMethodArguments(uint32 i) const {
643 // Return a list of contact model method argument names.
644 switch (i) {
645 case kwDeformability:
646 return "emod,kratio";
647 case kwBond:
648 return "gap,soft,cut";
649 case kwUnbond:
650 return "gap";
651 case kwArea:
652 return string();
653 }
654 assert(0);
655 return string();
656 }
657
658 bool ContactModelSoftBond::setMethod(uint32 i,const std::vector<base::Property> &vl,IContact *con) {
659 FP_S;
660 // Apply the specified method.
661 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
662 switch (i) {
663 case kwDeformability: {
664 double emod;
665 double krat;
666 if (vl.at(0).isNull())
667 throw Exception("Argument emod must be specified with method deformability in contact model {0}.",getName());
668 emod = vl.at(0).toDouble();
669 if (emod<0.0)
670 throw Exception("Negative emod not allowed in contact model {0}.",getName());
671 if (vl.at(1).isNull())
672 throw Exception("Argument kratio must be specified with method deformability in contact model {0}.",getName());
673 krat = vl.at(1).toDouble();
674 if (krat<0.0)
675 throw Exception("Negative stiffness ratio not allowed in contact model {0}.",getName());
676 double rsum = calcRSum(c);
677 if (userArea_)
678#ifdef THREED
679 rsum = 2.0 * std::sqrt(userArea_ / dPi);
680#else
681 rsum = userArea_;
682#endif
683 kn_ = emod / rsum;
684 ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
685 setInheritance(1,false);
686 setInheritance(2,false);
687 FP_S;
688 return true;
689 }
690 case kwBond: {
691 if (bProps_ && bProps_->sb_state_ >= 3) return false;
692 double mingap = -1.0 * limits<double>::max();
693 double maxgap = 0;
694 if (vl.at(0).canConvert<double>())
695 maxgap = vl.at(0).toDouble();
696 else if (vl.at(0).canConvert<DVect2>()) {
697 DVect2 value = vl.at(0).value<DVect2>();
698 mingap = value.minComp();
699 maxgap = value.maxComp();
700 }
701 else if (!vl.at(0).isNull())
702 throw Exception("gap value {0} not recognized in method bond in contact model {1}.", vl.at(1), getName());
703 double soft = bProps_ ? bProps_->sb_soft_ : 0.0;
704 if (!vl.at(1).isNull()) {
705 soft = vl.at(1).toDouble();
706 if (soft < 0.0)
707 throw Exception("Negative soft not allowed in contact model {0}.", getName());
708 }
709 double cut = bProps_ ? bProps_->sb_cut_ : 1.0;
710 if (!vl.at(2).isNull()) {
711 if (vl.at(2).canConvert<double>())
712 cut = vl.at(2).toDouble();
713 if (cut < 0.0)
714 throw Exception("cut value {0} is negative, or not recognized in method bond in contact model {1}.", vl.at(2), getName());
715 if (cut > 1.0)
716 throw Exception("cut value {0} must be in range [0,1] in method bond in contact model {1}.", vl.at(2), getName());
717 }
718 double gap = c->getGap();
719 if (gap >= mingap && gap <= maxgap) {
720 if (!bProps_)
721 bProps_ = NEW bProps();
722 bProps_->sb_state_ = 3;
723 bProps_->sb_soft_ = soft;
724 // Update the critical distance
725 if (cut != -1)
726 bProps_->sb_cut_ = cut;
727 // seet to incremental normal force calculation
728 sb_mode_ = 1;
729 FP_S;
730 return true;
731 }
732 FP_S;
733 return false;
734 }
735 case kwUnbond: {
736 if (!bProps_ || bProps_->sb_state_ == 0) return false;
737 double mingap = -1.0 * limits<double>::max();
738 double maxgap = 0;
739 if (vl.at(0).canConvert<double>())
740 maxgap = vl.at(0).toDouble();
741 else if (vl.at(0).canConvert<DVect2>()) {
742 DVect2 value = vl.at(0).value<DVect2>();
743 mingap = value.minComp();
744 maxgap = value.maxComp();
745 }
746 else if (!vl.at(0).isNull())
747 throw Exception("gap value {0} not recognized in method unbond in contact model {1}.", vl.at(0), getName());
748 double gap = c->getGap();
749 if (gap >= mingap && gap <= maxgap) {
750 bProps_->sb_state_ = 0;
751 FP_S;
752 return true;
753 }
754 FP_S;
755 return false;
756 }
757 case kwArea: {
758 if (!userArea_) {
759 double rsq = calcRSQ(c);
760#ifdef THREED
761 userArea_ = rsq * rsq * dPi;
762#else
763 userArea_ = rsq * 2.0;
764#endif
765 }
766 FP_S;
767 return true;
768 }
769
770 }
771 FP_S;
772 return false;
773 }
774
775 double ContactModelSoftBond::getEnergy(uint32 i) const {
776 // Return an energy value.
777 double ret(0.0);
778 if (!energies_)
779 return ret;
780 switch (i) {
781 case kwEStrain: return energies_->estrain_;
782 case kwESlip: return energies_->eslip_;
783 case kwEDashpot: return energies_->edashpot_;
784 }
785 assert(0);
786 return ret;
787 }
788
789 bool ContactModelSoftBond::getEnergyAccumulate(uint32 i) const {
790 // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
791 switch (i) {
792 case kwEStrain: return false;
793 case kwESlip: return true;
794 case kwEDashpot: return true;
795 }
796 assert(0);
797 return false;
798 }
799
800 void ContactModelSoftBond::setEnergy(uint32 i,const double &d) {
801 // Set an energy value.
802 if (!energies_) return;
803 switch (i) {
804 case kwEStrain: energies_->estrain_ = d; return;
805 case kwESlip: energies_->eslip_ = d; return;
806 case kwEDashpot: energies_->edashpot_= d; return;
807 }
808 assert(0);
809 return;
810 }
811
812 bool ContactModelSoftBond::validate(ContactModelMechanicalState *state,const double &) {
813 // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
814 // (1) the contact is created, (2) a property of the contact that returns a true via
815 // the setProperty method has been modified and (3) when a set of cycles is executed
816 // via the {cycle N} command.
817 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
818 assert(state);
819 const IContactMechanical *c = state->getMechanicalContact();
820 assert(c);
821
822 if (state->trackEnergy_)
823 activateEnergy();
824
825 if (inheritanceField_ & KnMask)
826 updateKn(c);
827 if (inheritanceField_ & KsMask)
828 updateKs(c);
829 if (inheritanceField_ & FricMask)
830 updateFric(c);
831
832 updateStiffness(state);
833 return checkActivity(state->gap_);
834 }
835
836 static const string knstr("kn");
837 bool ContactModelSoftBond::updateKn(const IContactMechanical *con) {
838 assert(con);
839 base::Property v1 = con->getEnd1()->getProperty(knstr);
840 base::Property v2 = con->getEnd2()->getProperty(knstr);
841 if (!v1.isValid() || !v2.isValid())
842 return false;
843 double kn1 = v1.toDouble();
844 double kn2 = v2.toDouble();
845 double val = kn_;
846 if (kn1 && kn2)
847 kn_ = kn1*kn2/(kn1+kn2);
848 else if (kn1)
849 kn_ = kn1;
850 else if (kn2)
851 kn_ = kn2;
852 return ( (kn_ != val) );
853 }
854
855 static const string ksstr("ks");
856 bool ContactModelSoftBond::updateKs(const IContactMechanical *con) {
857 assert(con);
858 base::Property v1 = con->getEnd1()->getProperty(ksstr);
859 base::Property v2 = con->getEnd2()->getProperty(ksstr);
860 if (!v1.isValid() || !v2.isValid())
861 return false;
862 double ks1 = v1.toDouble();
863 double ks2 = v2.toDouble();
864 double val = ks_;
865 if (ks1 && ks2)
866 ks_ = ks1*ks2/(ks1+ks2);
867 else if (ks1)
868 ks_ = ks1;
869 else if (ks2)
870 ks_ = ks2;
871 return ( (ks_ != val) );
872 }
873
874 static const string fricstr("fric");
875 bool ContactModelSoftBond::updateFric(const IContactMechanical *con) {
876 assert(con);
877 base::Property v1 = con->getEnd1()->getProperty(fricstr);
878 base::Property v2 = con->getEnd2()->getProperty(fricstr);
879 if (!v1.isValid() || !v2.isValid())
880 return false;
881 double fric1 = std::max(0.0,v1.toDouble());
882 double fric2 = std::max(0.0,v2.toDouble());
883 double val = fric_;
884 fric_ = std::min(fric1,fric2);
885 return ( (fric_ != val) );
886 }
887
888 bool ContactModelSoftBond::endPropertyUpdated(const string &name,const IContactMechanical *c) {
889 // The endPropertyUpdated method is called whenever a surface property (with a name
890 // that matches an inheritable contact model property name) of one of the contacting
891 // pieces is modified. This allows the contact model to update its associated
892 // properties. The return value denotes whether or not the update has affected
893 // the time step computation (by having modified the translational or rotational
894 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
895 assert(c);
896 StringList availableProperties = split(replace(simplified(getProperties())," ",""),",");
897 auto idx = findRegex(availableProperties,name);
898 if (idx==string::npos) return false;
899 idx += 1;
900 bool ret=false;
901 switch(idx) {
902 case kwKn: { //kn
903 if (inheritanceField_ & KnMask)
904 ret = updateKn(c);
905 break;
906 }
907 case kwKs: { //ks
908 if (inheritanceField_ & KsMask)
909 ret =updateKs(c);
910 break;
911 }
912 case kwFric: { //fric
913 if (inheritanceField_ & FricMask)
914 updateFric(c);
915 break;
916 }
917 }
918 return ret;
919 }
920
921 ContactModelSoftBond::StiffData ContactModelSoftBond::computeStiffData(ContactModelMechanicalState *state) const {
922 // Update contact data
923 double Cmin1 = state->end1Curvature_.x();
924 double Cmax1 = state->end1Curvature_.y();
925 double Cmax2 = state->end2Curvature_.y();
926 double dthick = (Cmin1 == 0.0) ? 1.0 : 0.0;
927 double br = sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
928 if (userArea_)
929#ifdef THREED
930 br = std::sqrt(userArea_ / dPi);
931#else
932 br = userArea_ / 2.0;
933#endif
934 double br2 = br * br;
935 double area = dthick <= 0.0 ? dPi * br2 : 2.0*br*dthick;
936 double bi = dthick <= 0.0 ? 0.25*area*br2 : 2.0*br*br2*dthick / 3.0;
937 StiffData ret;
938 ret.reff_ = br;
939 ret.trans_ = DVect2(kn_ * area , ks_ * area);
940 ret.ang_ = DAVect(kn_ * bi);
941#if DIM==3
942 ret.ang_.rx() = ks_ * 2.0*bi;
943#endif
944 return ret;
945 }
946
947 void ContactModelSoftBond::updateStiffness(ContactModelMechanicalState *state) {
948 // first compute stiffness data
949 StiffData stiff = computeStiffData(state);
950 // Now calculate effective stiffness
951 DVect2 retT = stiff.trans_;
952 // correction if viscous damping active
953 if (dpProps_) {
954 DVect2 correct(1.0);
955 if (dpProps_->dp_nratio_)
956 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
957 if (dpProps_->dp_sratio_)
958 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
959 retT /= (correct*correct);
960 }
961 effectiveTranslationalStiffness_ = retT;
962 // Effective rotational stiffness (bending and twisting)
963 effectiveRotationalStiffness_ = stiff.ang_;
964 }
965
966 bool ContactModelSoftBond::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
967 assert(state);
968
969 if (state->activated()) {
970 // The contact was just activated from an inactive state
971 // Trigger the FISH callback if one is hooked up to the
972 // contact_activated event.
973 if (cmEvents_[fActivated] >= 0) {
974 auto c = state->getContact();
975 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
976 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
977 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
978 }
979 }
980 updateStiffness(state);
981 if (isBonded()) return FDLawBonded(state, timestep);
982 else return FDLawUnBonded(state, timestep);
983
984 }
985
986 bool ContactModelSoftBond::FDLawBonded(ContactModelMechanicalState *state, const double ×tep) {
987 // Relative translational/rotational displacement increments
988 DVect trans = state->relativeTranslationalIncrement_;
989 DAVect ang = state->relativeAngularIncrement_;
990
991 // Store previous force and moment
992 DVect sb_F_old = sb_F_;
993 DAVect sb_M_old = sb_M_;
994
995 // Update stiffness data
996 StiffData stiff = computeStiffData(state);
997 DVect3 geom = computeGeomData(state->end1Curvature_,state->end2Curvature_);
998 double area = geom.x();
999 double bi = geom.y();
1000 double br = geom.z();
1001 double kn = stiff.trans_.x();
1002 double ks = stiff.trans_.y();
1003 double kb = stiff.ang_.z();
1004#if DIM==3
1005 double kt = stiff.ang_.x();
1006#else
1007 double kt = 0.0;
1008#endif
1009
1010 double nsmax0 = -(sb_F_.x() / area) + bProps_->sb_mcf_* sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z()) * br / bi;
1011
1012 // incremental normal force calculation
1013 sb_F_.rx() -= trans.x() * kn;
1014
1015 // shear force calculation
1016 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1017 // Loop over the shear components (note: the 0 component is the normal component)
1018 // and calculate the shear force.
1019 for (int i = 1; i<dim; ++i)
1020 sb_F_.rdof(i) -= trans.dof(i) * ks;
1021
1022 // moment calculation
1023 sb_M_ -= ang * stiff.ang_;
1024 double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1025
1026 // maximum tensile stress at bond periphery
1027 double nsmax = -(sb_F_.x() / area) + bProps_->sb_mcf_* dbend * br / bi;
1028
1029 bool softened = false;
1030 // Mode check
1031 if (state->canFail_) {
1032 if (bProps_->sb_state_ == 3 || bProps_->sb_state_ == 5) {
1033 double compVal = bProps_->sb_state_ == 3 ? bProps_->sb_ten_ : bProps_->sb_maxTen_;
1034 if (nsmax >= compVal ) {
1035 // enter softening regime
1036 // current bond elongation when softening starts
1037 // This is the elongation at the bond periphery
1038 double ls = - sb_F_.x() / kn + bProps_->sb_mcf_*dbend* br / kb;
1039 bProps_->sb_maxTen_ = compVal;
1040 if (bProps_->sb_state_ == 3)
1041 bProps_->sb_critu_ = ls /**(1.0+bProps_->sb_soft_)*/;
1042 bProps_->sb_delu_ = 0.0;
1043 bProps_->sb_delo_ = Quat::identity();
1044 if (bProps_->sb_state_ == 5 && nsmax < bProps_->sb_maxTen_)
1045 softened = true;
1046 bProps_->sb_state_ = 4;
1047 }
1048 }
1049 }
1050
1051 if (bProps_->sb_state_ == 4 && !softened && !checktol(bProps_->sb_soft_,0.0,1.0,100.0)) {
1052 double ls = bProps_->sb_critu_;
1053 double lc = ls * (1.0+bProps_->sb_soft_);
1054 DVect normal(0.0);
1055 normal.rx() = 1.0;
1056 DVect backNormal = (bProps_->sb_delo_.getConj().rotate(normal)).unit();
1057 double bend = acos(std::clamp(normal|backNormal,-1.0,1.0));
1058 double l0 = ls + bProps_->sb_maxu_ + bProps_->sb_delu_ + br*abs(bend);
1059 bProps_->sb_delu_ += trans.x();
1060 bProps_->sb_delo_.increment(ang);
1061 // Take the current contact normal and rotate it in the opposite direction of
1062 // the orientation - get the angle of bend from there
1063 backNormal = (bProps_->sb_delo_.getConj().rotate(normal)).unit();
1064 bend = acos(std::clamp(normal|backNormal,-1.0,1.0));
1065 double l = ls + bProps_->sb_maxu_ + bProps_->sb_delu_ + br*abs(bend);
1066 // target tensile stress
1067 double ns = bProps_->sb_ten_*(lc-l) / (bProps_->sb_soft_*ls);
1068 if (ns > 0) {
1069 if (nsmax >= ns) {
1070 double fac = ns / nsmax;
1071 sb_F_.rx() = fac*sb_F_.x();
1072#if DIM==3
1073 sb_M_.ry() = fac*sb_M_.y();
1074#endif
1075 sb_M_.rz() = fac*sb_M_.z();
1076 } else {
1077 bProps_->sb_state_ = 5;
1078 bProps_->sb_maxTen_ = nsmax0;
1079 bProps_->sb_maxu_ = (l0-ls);
1080 }
1081 } else {
1082 sb_F_.rx() = 0.0;
1083#if DIM==3
1084 sb_M_.ry() = 0.0;
1085#endif
1086 sb_M_.rz() = 0.0;
1087 }
1088 }
1089
1090 if (state->canFail_) {
1091 /* check for normal failure */
1092 bool failed = false;
1093 if (bProps_->sb_state_ == 4) {
1094 double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1095 double nsmax = -(sb_F_.x() / area) + bProps_->sb_mcf_*dbend * br / bi;
1096 if (nsmax <= bProps_->sb_ten_*bProps_->sb_cut_ || checktol(bProps_->sb_soft_,0.0,1.0,100.0)) {
1097 // Failed in tension
1098 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1099 bProps_->sb_state_ = 1;
1100 sb_F_.fill(0.0);
1101 sb_M_.fill(0.0);
1102 failed = true;
1103 if (cmEvents_[fBondBreak] >= 0) {
1104 auto c = state->getContact();
1105 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1106 fish::Parameter((int64)bProps_->sb_state_),
1107 fish::Parameter(nsmax),
1108 fish::Parameter(se)
1109 };
1110 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1111 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1112 }
1113 }
1114 }
1115
1116 if (!failed) {
1117 /* check for shear failure */
1118 double dtwist = sb_M_.x();
1119 DVect bfs(sb_F_);
1120 bfs.rx() = 0.0;
1121 double dbfs = bfs.mag();
1122 double ssmax = dbfs / area + bProps_->sb_mcf_*std::abs(dtwist) * 0.5* br / bi;
1123 double ss = shearStrength(area);
1124 if (ss < 0)
1125 ss = 0;
1126 if (ss <= ssmax) {
1127 // Failed in shear
1128 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1129 bProps_->sb_state_ = 2;
1130 if (cmEvents_[fBondBreak] >= 0) {
1131 auto c = state->getContact();
1132 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1133 fish::Parameter((int64)bProps_->sb_state_),
1134 fish::Parameter(ss),
1135 fish::Parameter(se)
1136 };
1137 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1138 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1139 }
1140 // Resolve sliding.
1141 double crit = sb_F_.x() * fric_;
1142 if (crit < 0)
1143 crit = 0;
1144 DVect sforce = sb_F_; sforce.rx() = 0.0;
1145 // The is the magnitude of the shear force.
1146 double sfmag = sforce.mag();
1147 // Sliding occurs when the magnitude of the shear force is greater than the
1148 // critical value.
1149 if (sfmag > crit) {
1150 // Lower the shear force to the critical value for sliding.
1151 double rat = crit / sfmag;
1152 sforce *= rat;
1153 sforce.rx() = sb_F_.x();
1154 sb_F_ = sforce;
1155 sb_S_ = true;
1156 }
1157
1158 // Resolve bending
1159 crit = sb_bmul_*2.1*0.25*stiff.reff_*std::abs(sb_F_.x()); // Jiang 2015
1160 DAVect test = sb_M_;
1161#if DIM==3
1162 test.rx() = 0.0;
1163#endif
1164 double tmag = test.mag();
1165 if (tmag > crit) {
1166 // Lower the bending moment to the critical value for sliding.
1167 double rat = crit / tmag;
1168 test *= rat;
1169 sb_BS_ = true;
1170 }
1171 sb_M_.rz() = test.z();
1172#if DIM==3
1173 sb_M_.ry() = test.y();
1174 // Resolve twisting
1175 crit = sb_tmul_ * 0.65*fric_* stiff.reff_*std::abs(sb_F_.x()) ; // Jiang 2015
1176 tmag = std::abs(sb_M_.x());
1177 if (tmag > crit) {
1178 // Lower the shear force to the critical value for sliding.
1179 double rat = crit / tmag;
1180 tmag = sb_M_.x() * rat;
1181 sb_TS_ = true;
1182 } else
1183 tmag = sb_M_.x();
1184 sb_M_.rx() = tmag;
1185#endif
1186 }
1187 }
1188 }
1189
1190 // Account for dashpot forces if the dashpot structure has been defined.
1191 if (dpProps_) {
1192 dpProps_->dp_F_.fill(0.0);
1193 double vcn(0.0), vcs(0.0);
1194 // Calculate the damping coefficients.
1195 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1196 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1197 // First damp the shear components
1198 for (int i = 1; i<dim; ++i)
1199 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1200 // Damp the normal component
1201 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1202 // Need to change behavior based on the dp_mode.
1203 if (bProps_->sb_state_ < 3 && (dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1204 // Limit in tension if not bonded.
1205 if (dpProps_->dp_F_.x() + sb_F_.x() < 0)
1206 dpProps_->dp_F_.rx() = -sb_F_.rx();
1207 }
1208 if (sb_S_ && dpProps_->dp_mode_ > 1) {
1209 // Limit in shear if sliding.
1210 double dfn = dpProps_->dp_F_.rx();
1211 dpProps_->dp_F_.fill(0.0);
1212 dpProps_->dp_F_.rx() = dfn;
1213 }
1214 }
1215
1216 //Compute energies if energy tracking has been enabled.
1217 if (state->trackEnergy_) {
1218 assert(energies_);
1219 energies_->estrain_ = 0.0;
1220 if (kn)
1221 // Calculate the strain energy.
1222 energies_->estrain_ = 0.5*sb_F_.x()*sb_F_.x() / kn;
1223 if (ks) {
1224 DVect s = sb_F_;
1225 s.rx() = 0.0;
1226 double smag2 = s.mag2();
1227 // Add the shear component of the strain energy.
1228 energies_->estrain_ += 0.5*smag2 / ks;
1229
1230 if (sb_S_) {
1231 // If sliding calculate the slip energy and accumulate it.
1232 sb_F_old.rx() = 0.0;
1233 DVect avg_F_s = (s + sb_F_old)*0.5;
1234 DVect u_s_el = (s - sb_F_old) / ks;
1235 DVect u_s(0.0);
1236 for (int i = 1; i < dim; ++i)
1237 u_s.rdof(i) = trans.dof(i);
1238 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1239 }
1240 }
1241 // Add the bending/twisting resistance energy contributions.
1242 if (kb) {
1243 DAVect tmp = sb_M_;
1244#ifdef THREED
1245 tmp.rx() = 0.0;
1246#endif
1247 energies_->estrain_ += 0.5*tmp.mag2() / kb;
1248 if (sb_BS_) {
1249 // accumulate bending slip energy.
1250 DAVect tmp_old = sb_M_old;
1251#ifdef THREED
1252 tmp_old.rx() = 0.0;
1253#endif
1254 DAVect avg_M = (tmp + tmp_old)*0.5;
1255 DAVect t_s_el = (tmp - tmp_old) / kb;
1256 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1257 }
1258 }
1259#ifdef THREED
1260 if (kt) {
1261 double mt = std::abs(sb_M_.x());
1262 energies_->estrain_ += 0.5*mt*mt / kt;
1263 if (sb_TS_) {
1264 // accumulate twisting slip energy.
1265 DAVect tmp(0.0);
1266 DAVect tmp_old(0.0);
1267 tmp.rx() = sb_M_.x();
1268 tmp_old.rx() = sb_M_old.x();
1269 DAVect avg_M = (tmp + tmp_old)*0.5;
1270 DAVect t_s_el = (tmp - tmp_old) / kt;
1271 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1272 }
1273 }
1274#endif
1275
1276 if (dpProps_) {
1277 // Calculate damping energy (accumulated) if the dashpots are active.
1278 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1279 }
1280 }
1281
1282 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1283 assert(sb_F_ == sb_F_);
1284 assert(sb_M_ == sb_M_);
1285 return true;
1286 }
1287
1288 bool ContactModelSoftBond::FDLawUnBonded(ContactModelMechanicalState *state, const double ×tep) {
1289
1290 // Relative translational/rotational displacement increments
1291 DVect trans = state->relativeTranslationalIncrement_;
1292 DAVect ang = state->relativeAngularIncrement_;
1293 double overlap = rgap_ - state->gap_;
1294 double correction = 1.0;
1295 if (state->activated() && sb_mode_ == 0 && trans.x()) {
1296 correction = -1.0*overlap / trans.x();
1297 if (correction < 0)
1298 correction = 1.0;
1299 }
1300
1301 // Store previous force and moment
1302 DVect sb_F_old = sb_F_;
1303 DAVect sb_M_old = sb_M_;
1304
1305 // Update stiffness data
1306 StiffData stiff = computeStiffData(state);
1307 double kn = stiff.trans_.x();
1308 double ks = stiff.trans_.y();
1309 double kb = stiff.ang_.z();
1310#if DIM==3
1311 double kt = stiff.ang_.x();
1312#endif
1313 // absolute/incremental normal force calculation
1314 if (sb_mode_==0)
1315 sb_F_.rx() = overlap * kn;
1316 else
1317 sb_F_.rx() -= trans.x() * kn;
1318 // Normal force can only be positive if unbonded
1319 sb_F_.rx() = std::max(0.0, sb_F_.x());
1320
1321 // Calculate the trial shear force.
1322 DVect sforce(0.0);
1323 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1324 // Loop over the shear components (note: the 0 component is the normal component)
1325 // and calculate the shear force.
1326 for (int i = 1; i<dim; ++i)
1327 sforce.rdof(i) = sb_F_.dof(i) - trans.dof(i) * ks;
1328
1329 // Calculate the trial moment.
1330 DAVect mom = sb_M_ - ang*stiff.ang_;
1331
1332 // If the SOLVE ELASTIC command is given then the
1333 // canFail state is set to FALSE. Otherwise it is always TRUE.
1334 if (state->canFail_) {
1335 bool changed = false;
1336 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
1337 bool slip_changed = false;
1338 double crit = sb_F_.x() * fric_;
1339 // The is the magnitude of the shear force.
1340 double sfmag = sforce.mag();
1341 // Sliding occurs when the magnitude of the shear force is greater than the
1342 // critical value.
1343 if (sfmag > crit) {
1344 // Lower the shear force to the critical value for sliding.
1345 double rat = crit / sfmag;
1346 sforce *= rat;
1347 if (!sb_S_) {
1348 slip_changed = true;
1349 changed = true;
1350 }
1351 sb_S_ = true;
1352 }
1353 else {
1354 if (sb_S_) {
1355 slip_changed = true;
1356 changed = true;
1357 }
1358 sb_S_ = false;
1359 }
1360
1361 // Resolve bending
1362 bool bslip_changed = false;
1363 crit = sb_bmul_ * 2.1*0.25*sb_F_.x() * stiff.reff_; // Jiang 2015
1364 DAVect test = mom;
1365#if DIM==3
1366 test.rx() = 0.0;
1367#endif
1368 double tmag = test.mag();
1369 if (tmag > crit) {
1370 // Lower the bending moment to the critical value for sliding.
1371 double rat = crit / tmag;
1372 test *= rat;
1373 if (!sb_BS_) {
1374 bslip_changed = true;
1375 changed = true;
1376 }
1377 sb_BS_ = true;
1378 }
1379 else {
1380 if (sb_BS_) {
1381 bslip_changed = true;
1382 changed = true;
1383 }
1384 sb_BS_ = false;
1385 }
1386 mom.rz() = test.z();
1387#if DIM==3
1388 mom.ry() = test.y();
1389 // Resolve twisting
1390 bool tslip_changed = false;
1391 crit = sb_tmul_ * 0.65*fric_*sb_F_.x() * stiff.reff_; // Jiang 2015
1392 tmag = std::abs(mom.x());
1393 if (tmag > crit) {
1394 // Lower the twisting moment to the critical value for sliding.
1395 double rat = crit / tmag;
1396 mom.rx() *= rat;
1397 if (!sb_TS_) {
1398 tslip_changed = true;
1399 changed = true;
1400 }
1401 sb_TS_ = true;
1402 } else {
1403 if (sb_TS_) {
1404 tslip_changed = true;
1405 changed = true;
1406 }
1407 sb_TS_ = false;
1408 }
1409#endif
1410 if (changed && cmEvents_[fSlipChange] >= 0) {
1411 int64 code = 0;
1412 if (slip_changed) {
1413 code = 1;
1414 if (bslip_changed) {
1415 code = 4;
1416#if DIM==3
1417 if (tslip_changed)
1418 code = 7;
1419#endif
1420 }
1421 }
1422 else if (bslip_changed) {
1423 code = 2;
1424#if DIM==3
1425 if (tslip_changed)
1426 code = 6;
1427#endif
1428 }
1429#if DIM==3
1430 else if (tslip_changed) {
1431 code = 3;
1432 if (slip_changed)
1433 code = 5;
1434 }
1435#endif
1436 auto c = state->getContact();
1437 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1438 fish::Parameter(code),
1439 fish::Parameter(sb_S_),
1440 fish::Parameter(sb_BS_)
1441#ifdef THREED
1442 ,fish::Parameter(sb_TS_)
1443#endif
1444 };
1445 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1446 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1447 }
1448 }
1449
1450 // Set the shear components of the total force.
1451 for (int i = 1; i<dim; ++i)
1452 sb_F_.rdof(i) = sforce.dof(i);
1453
1454 // Set the moment.
1455 sb_M_ = mom;
1456
1457 // Account for dashpot forces if the dashpot structure has been defined.
1458 if (dpProps_) {
1459 dpProps_->dp_F_.fill(0.0);
1460 double vcn(0.0), vcs(0.0);
1461 // Calculate the damping coefficients.
1462 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1463 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1464 // First damp the shear components
1465 for (int i = 1; i<dim; ++i)
1466 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1467 // Damp the normal component
1468 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1469 // Need to change behavior based on the dp_mode.
1470 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1471 // Limit in tension if not bonded.
1472 if (dpProps_->dp_F_.x() + sb_F_.x() < 0)
1473 dpProps_->dp_F_.rx() = -sb_F_.rx();
1474 }
1475 if (sb_S_ && dpProps_->dp_mode_ > 1) {
1476 // Limit in shear if not sliding.
1477 double dfn = dpProps_->dp_F_.rx();
1478 dpProps_->dp_F_.fill(0.0);
1479 dpProps_->dp_F_.rx() = dfn;
1480 }
1481 }
1482
1483 //Compute energies if energy tracking has been enabled.
1484 if (state->trackEnergy_) {
1485 assert(energies_);
1486 energies_->estrain_ = 0.0;
1487 if (kn_)
1488 // Calculate the strain energy.
1489 energies_->estrain_ = 0.5*sb_F_.x()*sb_F_.x() / kn;
1490 if (ks_) {
1491 DVect s = sb_F_;
1492 s.rx() = 0.0;
1493 double smag2 = s.mag2();
1494 // Add the shear component of the strain energy.
1495 energies_->estrain_ += 0.5*smag2 / ks;
1496
1497 if (sb_S_) {
1498 // If sliding calculate the slip energy and accumulate it.
1499 sb_F_old.rx() = 0.0;
1500 DVect avg_F_s = (s + sb_F_old)*0.5;
1501 DVect u_s_el = (s - sb_F_old) / ks;
1502 DVect u_s(0.0);
1503 for (int i = 1; i < dim; ++i)
1504 u_s.rdof(i) = trans.dof(i);
1505 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1506 }
1507 }
1508 // Add the bending/twisting resistance energy contributions.
1509 if (kb) {
1510 DAVect tmp = sb_M_;
1511#ifdef THREED
1512 tmp.rx() = 0.0;
1513#endif
1514 energies_->estrain_ += 0.5*tmp.mag2() / kb;
1515 if (sb_BS_) {
1516 // accumulate bending slip energy.
1517 DAVect tmp_old = sb_M_old;
1518#ifdef THREED
1519 tmp_old.rx() = 0.0;
1520#endif
1521 DAVect avg_M = (tmp + tmp_old)*0.5;
1522 DAVect t_s_el = (tmp - tmp_old) / kb;
1523 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1524 }
1525 }
1526#ifdef THREED
1527 if (kt) {
1528 double mt = std::abs(sb_M_.x());
1529 energies_->estrain_ += 0.5*mt*mt / kt;
1530 if (sb_TS_) {
1531 // accumulate twisting slip energy.
1532 DAVect tmp(0.0);
1533 DAVect tmp_old(0.0);
1534 tmp.rx() = sb_M_.x();
1535 tmp_old.rx() = sb_M_old.x();
1536 DAVect avg_M = (tmp + tmp_old)*0.5;
1537 DAVect t_s_el = (tmp - tmp_old) / kt;
1538 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1539 }
1540 }
1541#endif
1542
1543 if (dpProps_) {
1544 // Calculate damping energy (accumulated) if the dashpots are active.
1545 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1546 }
1547 }
1548
1549 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1550 assert(sb_F_ == sb_F_);
1551 assert(sb_M_ == sb_M_);
1552 return true;
1553 }
1554
1555 bool ContactModelSoftBond::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
1556 // Account for thermal expansion in incremental mode
1557 if (sb_mode_ == 0 || ts->gapInc_ == 0.0) return false;
1558 DVect finc(0.0);
1559 finc.rx() = kn_ * ts->gapInc_;
1560 sb_F_ -= finc;
1561 return true;
1562 }
1563
1564 void ContactModelSoftBond::setForce(const DVect &v,IContact *c) {
1565 sb_F(v);
1566 if (v.x() > 0) {
1567 auto con = convert_getcast<IContactMechanical>(c);
1568 rgap_ = c->getGap() + v.x() / (kn_ * computeGeomData(con->getEnd1Curvature(),con->getEnd2Curvature()).x());
1569 }
1570 }
1571
1572 void ContactModelSoftBond::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
1573 // Only called for contacts with wall facets when the wall resolution scheme
1574 // is set to full!
1575 // Only do something if the contact model is of the same type
1576 if (equal(old->getContactModel()->getName(),"softbond") && !isBonded()) {
1577 ContactModelSoftBond *oldCm = (ContactModelSoftBond *)old;
1578#ifdef THREED
1579 // Need to rotate just the shear component from oldSystem to newSystem
1580
1581 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
1582 DVect axis = oldSystem.e1() & newSystem.e1();
1583 double c, ang, s;
1584 DVect re2;
1585 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1586 axis = axis.unit();
1587 c = oldSystem.e1()|newSystem.e1();
1588 if (c > 0)
1589 c = std::min(c,1.0);
1590 else
1591 c = std::max(c,-1.0);
1592 ang = acos(c);
1593 s = sin(ang);
1594 double t = 1. - c;
1595 DMatrix<3,3> rm;
1596 rm.get(0,0) = t*axis.x()*axis.x() + c;
1597 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
1598 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
1599 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
1600 rm.get(1,1) = t*axis.y()*axis.y() + c;
1601 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
1602 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
1603 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
1604 rm.get(2,2) = t*axis.z()*axis.z() + c;
1605 re2 = rm*oldSystem.e2();
1606 }
1607 else
1608 re2 = oldSystem.e2();
1609 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
1610 axis = re2 & newSystem.e2();
1611 DVect2 tpf;
1612 DVect2 tpm;
1613 DMatrix<2,2> m;
1614 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1615 axis = axis.unit();
1616 c = re2|newSystem.e2();
1617 if (c > 0)
1618 c = std::min(c,1.0);
1619 else
1620 c = std::max(c,-1.0);
1621 ang = acos(c);
1622 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
1623 ang *= -1;
1624 s = sin(ang);
1625 m.get(0,0) = c;
1626 m.get(1,0) = s;
1627 m.get(0,1) = -m.get(1,0);
1628 m.get(1,1) = m.get(0,0);
1629 tpf = m*DVect2(oldCm->sb_F_.y(),oldCm->sb_F_.z());
1630 tpm = m*DVect2(oldCm->sb_M_.y(),oldCm->sb_M_.z());
1631 } else {
1632 m.get(0,0) = 1.;
1633 m.get(0,1) = 0.;
1634 m.get(1,0) = 0.;
1635 m.get(1,1) = 1.;
1636 tpf = DVect2(oldCm->sb_F_.y(),oldCm->sb_F_.z());
1637 tpm = DVect2(oldCm->sb_M_.y(),oldCm->sb_M_.z());
1638 }
1639 DVect pforce = DVect(0,tpf.x(),tpf.y());
1640 //DVect pm = DVect(0,tpm.x(),tpm.y());
1641#else
1642 oldSystem;
1643 newSystem;
1644 DVect pforce = DVect(0,oldCm->sb_F_.y());
1645 //DVect pm = DVect(0,oldCm->sb_M_.y());
1646#endif
1647 for (int i=1; i<dim; ++i)
1648 sb_F_.rdof(i) += pforce.dof(i);
1649 if (sb_mode_ && oldCm->sb_mode_)
1650 sb_F_.rx() = oldCm->sb_F_.x();
1651 oldCm->sb_F_ = DVect(0.0);
1652 oldCm->sb_M_ = DAVect(0.0);
1653 if (dpProps_ && oldCm->dpProps_) {
1654#ifdef THREED
1655 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
1656 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
1657#else
1658 pforce = oldCm->dpProps_->dp_F_;
1659#endif
1660 dpProps_->dp_F_ += pforce;
1661 oldCm->dpProps_->dp_F_ = DVect(0.0);
1662 }
1663 if(oldCm->getEnergyActivated()) {
1664 activateEnergy();
1665 energies_->estrain_ = oldCm->energies_->estrain_;
1666 energies_->edashpot_ = oldCm->energies_->edashpot_;
1667 energies_->eslip_ = oldCm->energies_->eslip_;
1668 oldCm->energies_->estrain_ = 0.0;
1669 oldCm->energies_->edashpot_ = 0.0;
1670 oldCm->energies_->eslip_ = 0.0;
1671 }
1672 rgap_ = oldCm->rgap_;
1673 }
1674 assert(sb_F_ == sb_F_);
1675 }
1676
1677 void ContactModelSoftBond::setNonForcePropsFrom(IContactModel *old) {
1678 // Only called for contacts with wall facets when the wall resolution scheme
1679 // is set to full!
1680 // Only do something if the contact model is of the same type
1681 if (equal(old->getName(),"softbond") && !isBonded()) {
1682 ContactModelSoftBond *oldCm = (ContactModelSoftBond *)old;
1683 kn_ = oldCm->kn_;
1684 ks_ = oldCm->ks_;
1685 fric_ = oldCm->fric_;
1686 sb_bmul_ = oldCm->sb_bmul_;
1687 sb_tmul_ = oldCm->sb_tmul_;
1688 sb_mode_ = oldCm->sb_mode_;
1689 sb_rmul_ = oldCm->sb_rmul_;
1690 sb_S_ = oldCm->sb_S_;
1691 sb_BS_ = oldCm->sb_BS_;
1692 sb_TS_ = oldCm->sb_TS_;
1693 rgap_ = oldCm->rgap_;
1694 userArea_ = oldCm->userArea_;
1695
1696 if (oldCm->dpProps_) {
1697 if (!dpProps_)
1698 dpProps_ = NEW dpProps();
1699 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1700 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1701 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1702 }
1703 if (oldCm->bProps_) {
1704 if (!bProps_)
1705 bProps_ = NEW bProps();
1706 bProps_->sb_mcf_ = oldCm->bProps_->sb_mcf_;
1707 bProps_->sb_fa_ = oldCm->bProps_->sb_fa_;
1708 bProps_->sb_state_ = oldCm->bProps_->sb_state_;
1709 bProps_->sb_coh_ = oldCm->bProps_->sb_coh_;
1710 bProps_->sb_ten_ = oldCm->bProps_->sb_ten_;
1711 bProps_->sb_maxTen_ = oldCm->bProps_->sb_maxTen_;
1712 bProps_->sb_cut_ = oldCm->bProps_->sb_cut_;
1713 bProps_->sb_delu_ = oldCm->bProps_->sb_delu_;
1714 bProps_->sb_delo_ = oldCm->bProps_->sb_delo_;
1715 bProps_->sb_maxu_ = oldCm->bProps_->sb_maxu_;
1716 bProps_->sb_critu_ = oldCm->bProps_->sb_critu_;
1717 }
1718
1719 }
1720 }
1721
1722 DVect ContactModelSoftBond::getForce() const {
1723 DVect ret(sb_F_);
1724 if (dpProps_)
1725 ret += dpProps_->dp_F_;
1726 return ret;
1727 }
1728
1729 DAVect ContactModelSoftBond::getMomentOn1(const IContactMechanical *c) const {
1730 DAVect ret(sb_M_);
1731 if (c) {
1732 DVect force = getForce();
1733 c->updateResultingTorqueOn1Local(force,&ret);
1734 }
1735 return ret;
1736 }
1737
1738 DAVect ContactModelSoftBond::getMomentOn2(const IContactMechanical *c) const {
1739 DAVect ret(sb_M_);
1740 if (c) {
1741 DVect force = getForce();
1742 c->updateResultingTorqueOn2Local(force,&ret);
1743 }
1744 return ret;
1745 }
1746
1747 DAVect ContactModelSoftBond::getModelMomentOn1() const {
1748 DAVect ret(sb_M_);
1749 return ret;
1750 }
1751
1752 DAVect ContactModelSoftBond::getModelMomentOn2() const {
1753 DAVect ret(sb_M_);
1754 return ret;
1755 }
1756
1757 void ContactModelSoftBond::objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *ret) const {
1758 ret->clear();
1759 ret->push_back({"kn",ScalarInfo});
1760 ret->push_back({"ks",ScalarInfo});
1761 ret->push_back({"fric",ScalarInfo});
1762 ret->push_back({"sb_bmul",ScalarInfo});
1763 ret->push_back({"sb_tmul",ScalarInfo});
1764 ret->push_back({"sb_mode",ScalarInfo});
1765 ret->push_back({"sb_force",VectorInfo});
1766 ret->push_back({"sb_moment",VectorInfo});
1767 ret->push_back({"sb_slip",ScalarInfo});
1768 ret->push_back({"sb_slipb",ScalarInfo});
1769 ret->push_back({"sb_slipt",ScalarInfo});
1770 ret->push_back({"sb_rmul",ScalarInfo});
1771 ret->push_back({"sb_radius",ScalarInfo});
1772 ret->push_back({"emod",ScalarInfo});
1773 ret->push_back({"kratio",ScalarInfo});
1774 ret->push_back({"dp_nratio",ScalarInfo});
1775 ret->push_back({"dp_sratio",ScalarInfo});
1776 ret->push_back({"dp_mode",ScalarInfo});
1777 ret->push_back({"dp_force",VectorInfo});
1778 ret->push_back({"sb_state",ScalarInfo});
1779 ret->push_back({"sb_ten",ScalarInfo});
1780 ret->push_back({"sb_shear",ScalarInfo});
1781 ret->push_back({"sb_coh",ScalarInfo});
1782 ret->push_back({"sb_fa",ScalarInfo});
1783 ret->push_back({"sb_mcf",ScalarInfo});
1784 ret->push_back({"sb_sigma",ScalarInfo});
1785 ret->push_back({"sb_tau",ScalarInfo});
1786 ret->push_back({"sb_soft",ScalarInfo});
1787 ret->push_back({"sb_cut",ScalarInfo});
1788 ret->push_back({"sb_area",ScalarInfo});
1789 ret->push_back({"user_area",ScalarInfo});
1790 ret->push_back({"rgap",ScalarInfo});
1791
1792 }
1793
1794 void ContactModelSoftBond::objectPropValues(std::vector<double> *ret,const IContact *con) const {
1795 FP_S;
1796 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
1797 ret->clear();
1798 ret->push_back(kn_);
1799 ret->push_back(ks_);
1800 ret->push_back(fric_);
1801 ret->push_back(sb_bmul_);
1802 ret->push_back(sb_tmul_);
1803 ret->push_back(sb_mode_);
1804 ret->push_back(sb_F_.mag());
1805 ret->push_back(sb_M_.mag());
1806 ret->push_back(sb_S_);
1807 ret->push_back(sb_BS_);
1808 ret->push_back(sb_TS_);
1809 ret->push_back(sb_rmul_);
1810 double d;
1811 double Cmax1 = c->getEnd1Curvature().y();
1812 double Cmax2 = c->getEnd2Curvature().y();
1813 if (!userArea_)
1814 d= sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
1815 else {
1816#ifdef THREED
1817 d = std::sqrt(userArea_ / dPi);
1818#else
1819 d = userArea_ / 2.0;
1820#endif
1821 }
1822 ret->push_back(d);
1823 double rsum = calcRSum(c);
1824 if (userArea_)
1825#ifdef THREED
1826 rsum = 2.0 * std::sqrt(userArea_ / dPi);
1827#else
1828 rsum = userArea_;
1829#endif
1830 d= kn_ * rsum;
1831 ret->push_back(d);
1832 ret->push_back((ks_ == 0.0) ? 0.0 : (kn_/ks_));
1833 ret->push_back(dpProps_ ? dpProps_->dp_nratio_ : 0);
1834 ret->push_back(dpProps_ ? dpProps_->dp_sratio_ : 0);
1835 ret->push_back(dpProps_ ? dpProps_->dp_mode_ : 0);
1836 ret->push_back(dpProps_ ? dpProps_->dp_F_.mag() : 0);
1837 ret->push_back(bProps_ ? bProps_->sb_state_ : 0);
1838 ret->push_back(bProps_ ? bProps_->sb_ten_ : 0);
1839 ret->push_back(shearStrength(computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x()));
1840 ret->push_back(bProps_ ? bProps_->sb_coh_ : 0);
1841 ret->push_back(bProps_ ? bProps_->sb_fa_ : 0);
1842 ret->push_back(bProps_ ? bProps_->sb_mcf_ : 0);
1843 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x());
1844 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).y());
1845 ret->push_back(bProps_ ? bProps_->sb_soft_ : 0);
1846 ret->push_back(bProps_ ? bProps_->sb_cut_ : 0);
1847 ret->push_back(userArea_ ? userArea_ : computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x());
1848 ret->push_back(userArea_);
1849 ret->push_back(rgap_);
1850 FP_S;
1851 }
1852
1853 DVect3 ContactModelSoftBond::computeGeomData(const DVect2 &end1Curvature,const DVect2 &end2Curvature) const {
1854 double Cmax1 = end1Curvature.y();
1855 double Cmax2 = end2Curvature.y();
1856 double br = sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
1857 if (userArea_)
1858#ifdef THREED
1859 br = std::sqrt(userArea_ / dPi);
1860#else
1861 br = userArea_ / 2.0;
1862#endif
1863 double br2 = br * br;
1864#ifdef TWOD
1865 double area = 2.0*br;
1866 double bi = 2.0*br*br2 / 3.0;
1867#else
1868 double area = dPi * br2;
1869 double bi = 0.25*area*br2;
1870#endif
1871 return DVect3(area, bi, br);
1872 }
1873
1874 DVect2 ContactModelSoftBond::SMax(const DVect2 &end1Curvature,const DVect2 &end2Curvature) const {
1875 DVect3 data = computeGeomData(end1Curvature,end2Curvature);
1876 double area = data.x();
1877 double bi = data.y();
1878 double br = data.z();
1879 /* maximum stresses */
1880 double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1881 double dtwist = sb_M_.x();
1882 DVect bfs(sb_F_);
1883 bfs.rx() = 0.0;
1884 double dbfs = bfs.mag();
1885 double nsmax = -(sb_F_.x() / area) + dbend * br / bi;
1886 double ssmax = dbfs / area + std::abs(dtwist) * 0.5* br / bi;
1887 return DVect2(nsmax, ssmax);
1888 }
1889
1890 double ContactModelSoftBond::shearStrength(const double &area) const {
1891 if (!bProps_) return 0.0;
1892 double sig = -1.0*sb_F_.x() / area;
1893 double nstr = bProps_->sb_state_ > 2 ? bProps_->sb_ten_ : 0.0;
1894 return sig <= nstr ? bProps_->sb_coh_ - std::tan(dDegrad*bProps_->sb_fa_)*sig
1895 : bProps_->sb_coh_ - std::tan(dDegrad*bProps_->sb_fa_)*nstr;
1896 }
1897
1898
1899 double ContactModelSoftBond::strainEnergy(double kn,double ks,double kb,double kt) const {
1900 double ret(0.0);
1901 if (kn)
1902 ret = 0.5 * sb_F_.x() * sb_F_.x() / kn;
1903 if (ks) {
1904 DVect tmp = sb_F_;
1905 tmp.rx() = 0.0;
1906 double smag2 = tmp.mag2();
1907 ret += 0.5 * smag2 / ks;
1908 }
1909
1910 if (kt)
1911 ret += 0.5 * sb_M_.x() * sb_M_.x() / kt;
1912 if (kb) {
1913 DAVect tmp = sb_M_;
1914#ifdef THREED
1915 tmp.rx() = 0.0;
1916 double smag2 = tmp.mag2();
1917#else
1918 double smag2 = tmp.z() * tmp.z();
1919#endif
1920 ret += 0.5 * smag2 / kb;
1921 }
1922 return ret;
1923 }
1924
1925} // namespace cmodelsxd
1926// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Jun 07, 2025 |
