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 QString 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 QString 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 QString 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_ = NEWC(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 QString getFishCallEvents() const override {
150 return
151 "contact_activated"
152 ",slip_change"
153 ",bond_break";
154 }
155
156 // Return the specified contact model property.
157 QVariant 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 QVariant &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 QString getMethods() const override { return "deformability,bond,unbond,area";}
190 // Return a comma seprated list of the contact model method arguments (base 1).
191 QString 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 QVector<QVariant> &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 QString &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 NEWC(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<QString,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 = NEWC(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_ = NEWC(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_ = NEWC(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_ = NEWC(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_ = NEWC(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_ = NEWC(Energies());
230 estrain(in->estrain());
231 eslip(in->eslip());
232 edashpot(in->edashpot());
233 }
234 if (in->hasBond()) {
235 if (!bProps_)
236 bProps_ = NEWC(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 QVariant 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 QVariant 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(0.0);
297 if (c->getEnd1Curvature().y())
298 rsum += 1.0/c->getEnd1Curvature().y();
299 if (c->getEnd2Curvature().y())
300 rsum += 1.0/c->getEnd2Curvature().y();
301 if (userArea_)
302#ifdef THREED
303 rsum = 2.0 * std::sqrt(userArea_ / dPi);
304#else
305 rsum = userArea_;
306#endif
307 return kn_ * rsum;
308 }
309 case kwKRatio: return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
310 case kwDpNRatio: return dpProps_ ? dpProps_->dp_nratio_ : 0;
311 case kwDpSRatio: return dpProps_ ? dpProps_->dp_sratio_ : 0;
312 case kwDpMode: return dpProps_ ? dpProps_->dp_mode_ : 0;
313 case kwDpF: {
314 dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
315 return var;
316 }
317 case kwSBState: return bProps_ ? bProps_->sb_state_ : 0;
318 case kwSBTStr: return bProps_ ? bProps_->sb_ten_ : 0.0;
319 case kwSBSStr: {
320 if (!bProps_ or not con) return 0.0;
321 double area = computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x();
322 return shearStrength(area);
323 }
324 case kwSBCoh: return bProps_ ? bProps_->sb_coh_ : 0;
325 case kwSBFa: return bProps_ ? bProps_->sb_fa_ : 0;
326 case kwSBMCF: return bProps_ ? bProps_->sb_mcf_ : 0;
327 case kwSBSig: {
328 if (!bProps_ or bProps_->sb_state_ < 3 or not con) return 0.0;
329 return SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x();
330 }
331 case kwSBTau: {
332 if (!bProps_ or bProps_->sb_state_ < 3 or not con) return 0.0;
333 return SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).y();
334 }
335 case kwSBSoft:
336 if (!bProps_) return 0.0;
337 return bProps_->sb_soft_;
338 case kwSBCut:
339 if (!bProps_) return 0.0;
340 return bProps_->sb_cut_;
341 case kwSBArea: {
342 if (userArea_) return userArea_;
343 if (!con)
344 return 0.0;
345 return computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x();
346 }
347 case kwUserArea:
348 return userArea_;
349 case kwRGap:
350 return rgap_;
351 }
352 assert(0);
353 return QVariant();
354 }
355
356 bool ContactModelSoftBond::getPropertyGlobal(uint32 i) const {
357 // Returns whether or not a property is held in the global axis system (TRUE)
358 // or the local system (FALSE). Used by the plotting logic.
359 switch (i) {
360 case kwSBF:
361 case kwSBM:
362 case kwDpF:
363 return false;
364 }
365 return true;
366 }
367
368 bool ContactModelSoftBond::setProperty(uint32 i,const QVariant &v,IContact *) {
369 // Set a contact model property. Return value indicates that the timestep
370 // should be recalculated.
371 dpProps dp;
372 switch (i) {
373 case kwKn: {
374 if (!v.canConvert<double>())
375 throw Exception("kn must be a double.");
376 double val(v.toDouble());
377 if (val<0.0)
378 throw Exception("Negative kn not allowed.");
379 kn_ = val;
380 return true;
381 }
382 case kwKs: {
383 if (!v.canConvert<double>())
384 throw Exception("ks must be a double.");
385 double val(v.toDouble());
386 if (val<0.0)
387 throw Exception("Negative ks not allowed.");
388 ks_ = val;
389 return true;
390 }
391 case kwFric: {
392 if (!v.canConvert<double>())
393 throw Exception("fric must be a double.");
394 double val(v.toDouble());
395 if (val<0.0)
396 throw Exception("Negative fric not allowed.");
397 fric_ = val;
398 return false;
399 }
400 case kwBMul: {
401 if (!v.canConvert<double>())
402 throw Exception("sb_bmul must be a double.");
403 double val(v.toDouble());
404 if (val<0.0)
405 throw Exception("Negative sb_bmul not allowed.");
406 sb_bmul_ = val;
407 return false;
408 }
409 case kwTMul: {
410 if (!v.canConvert<double>())
411 throw Exception("sb_tmul must be a double.");
412 double val(v.toDouble());
413 if (val<0.0)
414 throw Exception("Negative st_bmul not allowed.");
415 sb_tmul_ = val;
416 return false;
417 }
418 case kwSBMode: {
419 if (!v.canConvert<uint32>())
420 throw Exception("sb_mode must be 0 (absolute) or 1 (incremental).");
421 double val(v.toUInt());
422 if (val>1)
423 throw Exception("sb_mode must be 0 (absolute) or 1 (incremental).");
424 sb_mode_ = val;
425 return false;
426 }
427 case kwSBRMul: {
428 if (!v.canConvert<double>())
429 throw Exception("rmul must be a double.");
430 double val(v.toDouble());
431 if (val<0.0)
432 throw Exception("Negative rmul not allowed.");
433 sb_rmul_ = val;
434 return false;
435 }
436 case kwSBF: {
437 if (!v.canConvert<DVect>())
438 throw Exception("sb_force must be a vector.");
439 DVect val(v.value<DVect>());
440 sb_F_ = val;
441 return false;
442 }
443 case kwSBM: {
444 DAVect val(0.0);
445#ifdef TWOD
446 if (!v.canConvert<DAVect>() && !v.canConvert<double>())
447 throw Exception("res_moment must be an angular vector.");
448 if (v.canConvert<DAVect>())
449 val = DAVect(v.value<DAVect>());
450 else
451 val = DAVect(v.toDouble());
452#else
453 if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
454 throw Exception("res_moment must be an angular vector.");
455 if (v.canConvert<DAVect>())
456 val = DAVect(v.value<DAVect>());
457 else
458 val = DAVect(v.value<DVect>());
459#endif
460 sb_M_ = val;
461 return false;
462 }
463 case kwDpNRatio: {
464 if (!v.canConvert<double>())
465 throw Exception("dp_nratio must be a double.");
466 double val(v.toDouble());
467 if (val<0.0)
468 throw Exception("Negative dp_nratio not allowed.");
469 if (val == 0.0 && !dpProps_)
470 return false;
471 if (!dpProps_)
472 dpProps_ = NEWC(dpProps());
473 dpProps_->dp_nratio_ = val;
474 return true;
475 }
476 case kwDpSRatio: {
477 if (!v.canConvert<double>())
478 throw Exception("dp_sratio must be a double.");
479 double val(v.toDouble());
480 if (val<0.0)
481 throw Exception("Negative dp_sratio not allowed.");
482 if (val == 0.0 && !dpProps_)
483 return false;
484 if (!dpProps_)
485 dpProps_ = NEWC(dpProps());
486 dpProps_->dp_sratio_ = val;
487 return true;
488 }
489 case kwDpMode: {
490 if (!v.canConvert<int>())
491 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
492 int val(v.toInt());
493 if (val == 0 && !dpProps_)
494 return false;
495 if (val < 0 || val > 3)
496 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
497 if (!dpProps_)
498 dpProps_ = NEWC(dpProps());
499 dpProps_->dp_mode_ = val;
500 return false;
501 }
502 case kwDpF: {
503 if (!v.canConvert<DVect>())
504 throw Exception("dp_force must be a vector.");
505 DVect val(v.value<DVect>());
506 if (val.fsum() == 0.0 && !dpProps_)
507 return false;
508 if (!dpProps_)
509 dpProps_ = NEWC(dpProps());
510 dpProps_->dp_F_ = val;
511 return false;
512 }
513 case kwSBTStr: {
514 if (!v.canConvert<double>())
515 throw Exception("sb_ten must be a double.");
516 double val(v.toDouble());
517 if (val < 0.0)
518 throw Exception("Negative sb_ten not allowed.");
519 if (val == 0.0 && !bProps_)
520 return false;
521 if (!bProps_)
522 bProps_ = NEWC(bProps());
523 bProps_->sb_ten_ = val;
524 return false;
525 }
526 case kwSBCoh: {
527 if (!v.canConvert<double>())
528 throw Exception("sb_coh must be a double.");
529 double val(v.toDouble());
530 if (val<0.0)
531 throw Exception("Negative pb_coh not allowed.");
532 if (val == 0.0 && !bProps_)
533 return false;
534 if (!bProps_)
535 bProps_ = NEWC(bProps());
536 bProps_->sb_coh_ = val;
537 return false;
538 }
539 case kwSBFa: {
540 if (!v.canConvert<double>())
541 throw Exception("sb_fa must be a double.");
542 double val(v.toDouble());
543 if (val<0.0)
544 throw Exception("Negative sb_fa not allowed.");
545 if (val >= 90.0)
546 throw Exception("sb_fa must be lower than 90.0 degrees.");
547 if (val == 0.0 && !bProps_)
548 return false;
549 if (!bProps_)
550 bProps_ = NEWC(bProps());
551 bProps_->sb_fa_ = val;
552 return false;
553 }
554 case kwSBMCF: {
555 if (!v.canConvert<double>())
556 throw Exception("sb_mcf must be a double.");
557 double val(v.toDouble());
558 if (val<0.0)
559 throw Exception("Negative sb_mcf not allowed.");
560 if (val > 1.0)
561 throw Exception("sb_mcf must be lower or equal to 1.0.");
562 if (val == 1.0 && !bProps_)
563 return false;
564 if (!bProps_)
565 bProps_ = NEWC(bProps());
566 bProps_->sb_mcf_ = val;
567 return false;
568 }
569 case kwSBSoft: {
570 if (!v.canConvert<double>())
571 throw Exception("sb_soft must be a double.");
572 double val(v.toDouble());
573 if (val < 0.0)
574 throw Exception("Negative pb_soft not allowed.");
575 if (!bProps_)
576 bProps_ = NEWC(bProps());
577 bProps_->sb_soft_ = val;
578 return false;
579 }
580 case kwSBCut: {
581 if (!v.canConvert<double>())
582 throw Exception("sb_cut must be a double.");
583 double val(v.toDouble());
584 if (val < 0.0)
585 throw Exception("Negative sb_cut not allowed.");
586 if (!bProps_)
587 bProps_ = NEWC(bProps());
588 bProps_->sb_cut_ = val;
589 return false;
590 }
591 case kwSBArea:
592 case kwUserArea: {
593 if (!v.canConvert<double>())
594 throw Exception("area must be a double.");
595 double val(v.toDouble());
596 if (val < 0.0)
597 throw Exception("Negative area not allowed.");
598 userArea_ = val;
599 return true;
600 }
601 case kwRGap: {
602 if (!v.canConvert<double>())
603 throw Exception("Reference gap must be a double.");
604 double val(v.toDouble());
605 rgap_ = val;
606 return false;
607 }
608 }
609 return false;
610 }
611
612 bool ContactModelSoftBond::getPropertyReadOnly(uint32 i) const {
613 // Returns TRUE if a property is read only or FALSE otherwise.
614 switch (i) {
615 case kwDpF:
616 case kwSBS:
617 case kwSBBS:
618 case kwSBTS:
619 case kwEmod:
620 case kwKRatio:
621 case kwSBState:
622 case kwSBRadius:
623 case kwSBSStr:
624 case kwSBSig:
625 case kwSBTau:
626 return true;
627 default:
628 break;
629 }
630 return false;
631 }
632
633 bool ContactModelSoftBond::supportsInheritance(uint32 i) const {
634 // Returns TRUE if a property supports inheritance or FALSE otherwise.
635 switch (i) {
636 case kwKn:
637 case kwKs:
638 case kwFric:
639 return true;
640 default:
641 break;
642 }
643 return false;
644 }
645
646 QString ContactModelSoftBond::getMethodArguments(uint32 i) const {
647 // Return a list of contact model method argument names.
648 switch (i) {
649 case kwDeformability:
650 return "emod,kratio";
651 case kwBond:
652 return "gap,soft,cut";
653 case kwUnbond:
654 return "gap";
655 case kwArea:
656 return QString();
657 }
658 assert(0);
659 return QString();
660 }
661
662 bool ContactModelSoftBond::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
663 // Apply the specified method.
664 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
665 switch (i) {
666 case kwDeformability: {
667 double emod;
668 double krat;
669 if (vl.at(0).isNull())
670 throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
671 emod = vl.at(0).toDouble();
672 if (emod<0.0)
673 throw Exception("Negative emod not allowed in contact model %1.",getName());
674 if (vl.at(1).isNull())
675 throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
676 krat = vl.at(1).toDouble();
677 if (krat<0.0)
678 throw Exception("Negative stiffness ratio not allowed in contact model %1.",getName());
679 double rsum(0.0);
680 if (c->getEnd1Curvature().y())
681 rsum += 1.0 / c->getEnd1Curvature().y();
682 if (c->getEnd2Curvature().y())
683 rsum += 1.0 / c->getEnd2Curvature().y();
684 if (userArea_)
685#ifdef THREED
686 rsum = 2.0 * std::sqrt(userArea_ / dPi);
687#else
688 rsum = userArea_;
689#endif
690 kn_ = emod / rsum;
691 ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
692 setInheritance(1,false);
693 setInheritance(2,false);
694 return true;
695 }
696 case kwBond: {
697 if (bProps_ && bProps_->sb_state_ >= 3) return false;
698 double mingap = -1.0 * limits<double>::max();
699 double maxgap = 0;
700 if (vl.at(0).canConvert<double>())
701 maxgap = vl.at(0).toDouble();
702 else if (vl.at(0).canConvert<DVect2>()) {
703 DVect2 value = vl.at(0).value<DVect2>();
704 mingap = value.minComp();
705 maxgap = value.maxComp();
706 }
707 else if (!vl.at(0).isNull())
708 throw Exception("gap value %1 not recognized in method bond in contact model %2.", vl.at(1), getName());
709 double soft = bProps_ ? bProps_->sb_soft_ : 0.0;
710 if (!vl.at(1).isNull()) {
711 soft = vl.at(1).toDouble();
712 if (soft < 0.0)
713 throw Exception("Negative soft not allowed in contact model %1.", getName());
714 }
715 double cut = bProps_ ? bProps_->sb_cut_ : 1.0;
716 if (!vl.at(2).isNull()) {
717 if (vl.at(2).canConvert<double>())
718 cut = vl.at(2).toDouble();
719 if (cut < 0.0)
720 throw Exception("cut value %1 is negative, or not recognized in method bond in contact model %2.", vl.at(2), getName());
721 if (cut > 1.0)
722 throw Exception("cut value %1 must be in range [0,1] in method bond in contact model %2.", vl.at(2), getName());
723 }
724 double gap = c->getGap();
725 if (gap >= mingap && gap <= maxgap) {
726 if (!bProps_)
727 bProps_ = NEWC(bProps());
728 bProps_->sb_state_ = 3;
729 bProps_->sb_soft_ = soft;
730 // Update the critical distance
731 if (cut != -1)
732 bProps_->sb_cut_ = cut;
733 // seet to incremental normal force calculation
734 sb_mode_ = 1;
735 return true;
736 }
737 return false;
738 }
739 case kwUnbond: {
740 if (!bProps_ || bProps_->sb_state_ == 0) return false;
741 double mingap = -1.0 * limits<double>::max();
742 double maxgap = 0;
743 if (vl.at(0).canConvert<double>())
744 maxgap = vl.at(0).toDouble();
745 else if (vl.at(0).canConvert<DVect2>()) {
746 DVect2 value = vl.at(0).value<DVect2>();
747 mingap = value.minComp();
748 maxgap = value.maxComp();
749 }
750 else if (!vl.at(0).isNull())
751 throw Exception("gap value %1 not recognized in method unbond in contact model %2.", vl.at(0), getName());
752 double gap = c->getGap();
753 if (gap >= mingap && gap <= maxgap) {
754 bProps_->sb_state_ = 0;
755 return true;
756 }
757 return false;
758 }
759 case kwArea: {
760 if (!userArea_) {
761 double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
762#ifdef THREED
763 userArea_ = rsq * rsq * dPi;
764#else
765 userArea_ = rsq * 2.0;
766#endif
767 }
768 return true;
769 }
770
771 }
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 QString knstr("kn");
837 bool ContactModelSoftBond::updateKn(const IContactMechanical *con) {
838 assert(con);
839 QVariant v1 = con->getEnd1()->getProperty(knstr);
840 QVariant 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 QString ksstr("ks");
856 bool ContactModelSoftBond::updateKs(const IContactMechanical *con) {
857 assert(con);
858 QVariant v1 = con->getEnd1()->getProperty(ksstr);
859 QVariant 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 QString fricstr("fric");
875 bool ContactModelSoftBond::updateFric(const IContactMechanical *con) {
876 assert(con);
877 QVariant v1 = con->getEnd1()->getProperty(fricstr);
878 QVariant 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 QString &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 QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
897 QRegExp rx(name,Qt::CaseInsensitive);
898 int idx = availableProperties.indexOf(rx)+1;
899 bool ret=false;
900
901 if (idx<=0)
902 return ret;
903
904 switch(idx) {
905 case kwKn: { //kn
906 if (inheritanceField_ & KnMask)
907 ret = updateKn(c);
908 break;
909 }
910 case kwKs: { //ks
911 if (inheritanceField_ & KsMask)
912 ret =updateKs(c);
913 break;
914 }
915 case kwFric: { //fric
916 if (inheritanceField_ & FricMask)
917 updateFric(c);
918 break;
919 }
920 }
921 return ret;
922 }
923
924 ContactModelSoftBond::StiffData ContactModelSoftBond::computeStiffData(ContactModelMechanicalState *state) const {
925 // Update contact data
926 double Cmin1 = state->end1Curvature_.x();
927 double Cmax1 = state->end1Curvature_.y();
928 double Cmax2 = state->end2Curvature_.y();
929 double dthick = (Cmin1 == 0.0) ? 1.0 : 0.0;
930 double br = sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
931 if (userArea_)
932#ifdef THREED
933 br = std::sqrt(userArea_ / dPi);
934#else
935 br = userArea_ / 2.0;
936#endif
937 double br2 = br * br;
938 double area = dthick <= 0.0 ? dPi * br2 : 2.0*br*dthick;
939 double bi = dthick <= 0.0 ? 0.25*area*br2 : 2.0*br*br2*dthick / 3.0;
940 StiffData ret;
941 ret.reff_ = br;
942 ret.trans_ = DVect2(kn_ * area , ks_ * area);
943 ret.ang_ = DAVect(kn_ * bi);
944#if DIM==3
945 ret.ang_.rx() = ks_ * 2.0*bi;
946#endif
947 return ret;
948 }
949
950 void ContactModelSoftBond::updateStiffness(ContactModelMechanicalState *state) {
951 // first compute stiffness data
952 StiffData stiff = computeStiffData(state);
953 // Now calculate effective stiffness
954 DVect2 retT = stiff.trans_;
955 // correction if viscous damping active
956 if (dpProps_) {
957 DVect2 correct(1.0);
958 if (dpProps_->dp_nratio_)
959 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
960 if (dpProps_->dp_sratio_)
961 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
962 retT /= (correct*correct);
963 }
964 effectiveTranslationalStiffness_ = retT;
965 // Effective rotational stiffness (bending and twisting)
966 effectiveRotationalStiffness_ = stiff.ang_;
967 }
968
969 bool ContactModelSoftBond::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
970 assert(state);
971
972 if (state->activated()) {
973 // The contact was just activated from an inactive state
974 // Trigger the FISH callback if one is hooked up to the
975 // contact_activated event.
976 if (cmEvents_[fActivated] >= 0) {
977 auto c = state->getContact();
978 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
979 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
980 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
981 }
982 }
983 updateStiffness(state);
984 if (isBonded()) return FDLawBonded(state, timestep);
985 else return FDLawUnBonded(state, timestep);
986
987 }
988
989 bool ContactModelSoftBond::FDLawBonded(ContactModelMechanicalState *state, const double ×tep) {
990 // Relative translational/rotational displacement increments
991 DVect trans = state->relativeTranslationalIncrement_;
992 DAVect ang = state->relativeAngularIncrement_;
993
994 // Store previous force and moment
995 DVect sb_F_old = sb_F_;
996 DAVect sb_M_old = sb_M_;
997
998 // Update stiffness data
999 StiffData stiff = computeStiffData(state);
1000 DVect3 geom = computeGeomData(state->end1Curvature_,state->end2Curvature_);
1001 double area = geom.x();
1002 double bi = geom.y();
1003 double br = geom.z();
1004 double kn = stiff.trans_.x();
1005 double ks = stiff.trans_.y();
1006 double kb = stiff.ang_.z();
1007#if DIM==3
1008 double kt = stiff.ang_.x();
1009#else
1010 double kt = 0.0;
1011#endif
1012
1013 double nsmax0 = -(sb_F_.x() / area) + bProps_->sb_mcf_* sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z()) * br / bi;
1014
1015 // incremental normal force calculation
1016 sb_F_.rx() -= trans.x() * kn;
1017
1018 // shear force calculation
1019 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1020 // Loop over the shear components (note: the 0 component is the normal component)
1021 // and calculate the shear force.
1022 for (int i = 1; i<dim; ++i)
1023 sb_F_.rdof(i) -= trans.dof(i) * ks;
1024
1025 // moment calculation
1026 sb_M_ -= ang * stiff.ang_;
1027 double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1028
1029 // maximum tensile stress at bond periphery
1030 double nsmax = -(sb_F_.x() / area) + bProps_->sb_mcf_* dbend * br / bi;
1031
1032 bool softened = false;
1033 // Mode check
1034 if (state->canFail_) {
1035 if (bProps_->sb_state_ == 3 || bProps_->sb_state_ == 5) {
1036 double compVal = bProps_->sb_state_ == 3 ? bProps_->sb_ten_ : bProps_->sb_maxTen_;
1037 if (nsmax >= compVal ) {
1038 // enter softening regime
1039 // current bond elongation when softening starts
1040 // This is the elongation at the bond periphery
1041 double ls = - sb_F_.x() / kn + bProps_->sb_mcf_*dbend* br / kb;
1042 bProps_->sb_maxTen_ = compVal;
1043 if (bProps_->sb_state_ == 3)
1044 bProps_->sb_critu_ = ls /**(1.0+bProps_->sb_soft_)*/;
1045 bProps_->sb_delu_ = 0.0;
1046 bProps_->sb_delo_ = Quat::identity();
1047 if (bProps_->sb_state_ == 5 && nsmax < bProps_->sb_maxTen_)
1048 softened = true;
1049 bProps_->sb_state_ = 4;
1050 }
1051 }
1052 }
1053
1054 if (bProps_->sb_state_ == 4 && !softened && !checktol(bProps_->sb_soft_,0.0,1.0,100.0)) {
1055 double ls = bProps_->sb_critu_;
1056 double lc = ls * (1.0+bProps_->sb_soft_);
1057 DVect normal(0.0);
1058 normal.rx() = 1.0;
1059 DVect backNormal = (bProps_->sb_delo_.getConj().rotate(normal)).unit();
1060 double bend = acos(qBound(-1.0,normal|backNormal,1.0));
1061 double l0 = ls + bProps_->sb_maxu_ + bProps_->sb_delu_ + br*abs(bend);
1062 bProps_->sb_delu_ += trans.x();
1063 bProps_->sb_delo_.increment(ang);
1064 // Take the current contact normal and rotate it in the opposite direction of
1065 // the orientation - get the angle of bend from there
1066 backNormal = (bProps_->sb_delo_.getConj().rotate(normal)).unit();
1067 bend = acos(qBound(-1.0,normal|backNormal,1.0));
1068 double l = ls + bProps_->sb_maxu_ + bProps_->sb_delu_ + br*abs(bend);
1069 // target tensile stress
1070 double ns = bProps_->sb_ten_*(lc-l) / (bProps_->sb_soft_*ls);
1071 if (ns > 0) {
1072 if (nsmax >= ns) {
1073 double fac = ns / nsmax;
1074 sb_F_.rx() = fac*sb_F_.x();
1075#if DIM==3
1076 sb_M_.ry() = fac*sb_M_.y();
1077#endif
1078 sb_M_.rz() = fac*sb_M_.z();
1079 } else {
1080 bProps_->sb_state_ = 5;
1081 bProps_->sb_maxTen_ = nsmax0;
1082 bProps_->sb_maxu_ = (l0-ls);
1083 }
1084 } else {
1085 sb_F_.rx() = 0.0;
1086#if DIM==3
1087 sb_M_.ry() = 0.0;
1088#endif
1089 sb_M_.rz() = 0.0;
1090 }
1091 }
1092
1093 if (state->canFail_) {
1094 /* check for normal failure */
1095 bool failed = false;
1096 if (bProps_->sb_state_ == 4) {
1097 double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1098 double nsmax = -(sb_F_.x() / area) + bProps_->sb_mcf_*dbend * br / bi;
1099 if (nsmax <= bProps_->sb_ten_*bProps_->sb_cut_ || checktol(bProps_->sb_soft_,0.0,1.0,100.0)) {
1100 // Failed in tension
1101 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1102 bProps_->sb_state_ = 1;
1103 sb_F_.fill(0.0);
1104 sb_M_.fill(0.0);
1105 failed = true;
1106 if (cmEvents_[fBondBreak] >= 0) {
1107 auto c = state->getContact();
1108 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1109 fish::Parameter((qint64)bProps_->sb_state_),
1110 fish::Parameter(nsmax),
1111 fish::Parameter(se)
1112 };
1113 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1114 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1115 }
1116 }
1117 }
1118
1119 if (!failed) {
1120 /* check for shear failure */
1121 double dtwist = sb_M_.x();
1122 DVect bfs(sb_F_);
1123 bfs.rx() = 0.0;
1124 double dbfs = bfs.mag();
1125 double ssmax = dbfs / area + bProps_->sb_mcf_*std::abs(dtwist) * 0.5* br / bi;
1126 double ss = shearStrength(area);
1127 if (ss < 0)
1128 ss = 0;
1129 if (ss <= ssmax) {
1130 // Failed in shear
1131 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1132 bProps_->sb_state_ = 2;
1133 if (cmEvents_[fBondBreak] >= 0) {
1134 auto c = state->getContact();
1135 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1136 fish::Parameter((qint64)bProps_->sb_state_),
1137 fish::Parameter(ss),
1138 fish::Parameter(se)
1139 };
1140 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1141 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1142 }
1143 // Resolve sliding.
1144 double crit = sb_F_.x() * fric_;
1145 if (crit < 0)
1146 crit = 0;
1147 DVect sforce = sb_F_; sforce.rx() = 0.0;
1148 // The is the magnitude of the shear force.
1149 double sfmag = sforce.mag();
1150 // Sliding occurs when the magnitude of the shear force is greater than the
1151 // critical value.
1152 if (sfmag > crit) {
1153 // Lower the shear force to the critical value for sliding.
1154 double rat = crit / sfmag;
1155 sforce *= rat;
1156 sforce.rx() = sb_F_.x();
1157 sb_F_ = sforce;
1158 sb_S_ = true;
1159 }
1160
1161 // Resolve bending
1162 crit = sb_bmul_*2.1*0.25*stiff.reff_*std::abs(sb_F_.x()); // Jiang 2015
1163 DAVect test = sb_M_;
1164#if DIM==3
1165 test.rx() = 0.0;
1166#endif
1167 double tmag = test.mag();
1168 if (tmag > crit) {
1169 // Lower the bending moment to the critical value for sliding.
1170 double rat = crit / tmag;
1171 test *= rat;
1172 sb_BS_ = true;
1173 }
1174 sb_M_.rz() = test.z();
1175#if DIM==3
1176 sb_M_.ry() = test.y();
1177 // Resolve twisting
1178 crit = sb_tmul_ * 0.65*fric_* stiff.reff_*std::abs(sb_F_.x()) ; // Jiang 2015
1179 tmag = std::abs(sb_M_.x());
1180 if (tmag > crit) {
1181 // Lower the shear force to the critical value for sliding.
1182 double rat = crit / tmag;
1183 tmag = sb_M_.x() * rat;
1184 sb_TS_ = true;
1185 } else
1186 tmag = sb_M_.x();
1187 sb_M_.rx() = tmag;
1188#endif
1189 }
1190 }
1191 }
1192
1193 // Account for dashpot forces if the dashpot structure has been defined.
1194 if (dpProps_) {
1195 dpProps_->dp_F_.fill(0.0);
1196 double vcn(0.0), vcs(0.0);
1197 // Calculate the damping coefficients.
1198 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1199 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1200 // First damp the shear components
1201 for (int i = 1; i<dim; ++i)
1202 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1203 // Damp the normal component
1204 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1205 // Need to change behavior based on the dp_mode.
1206 if (bProps_->sb_state_ < 3 && (dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1207 // Limit in tension if not bonded.
1208 if (dpProps_->dp_F_.x() + sb_F_.x() < 0)
1209 dpProps_->dp_F_.rx() = -sb_F_.rx();
1210 }
1211 if (sb_S_ && dpProps_->dp_mode_ > 1) {
1212 // Limit in shear if sliding.
1213 double dfn = dpProps_->dp_F_.rx();
1214 dpProps_->dp_F_.fill(0.0);
1215 dpProps_->dp_F_.rx() = dfn;
1216 }
1217 }
1218
1219 //Compute energies if energy tracking has been enabled.
1220 if (state->trackEnergy_) {
1221 assert(energies_);
1222 energies_->estrain_ = 0.0;
1223 if (kn)
1224 // Calculate the strain energy.
1225 energies_->estrain_ = 0.5*sb_F_.x()*sb_F_.x() / kn;
1226 if (ks) {
1227 DVect s = sb_F_;
1228 s.rx() = 0.0;
1229 double smag2 = s.mag2();
1230 // Add the shear component of the strain energy.
1231 energies_->estrain_ += 0.5*smag2 / ks;
1232
1233 if (sb_S_) {
1234 // If sliding calculate the slip energy and accumulate it.
1235 sb_F_old.rx() = 0.0;
1236 DVect avg_F_s = (s + sb_F_old)*0.5;
1237 DVect u_s_el = (s - sb_F_old) / ks;
1238 DVect u_s(0.0);
1239 for (int i = 1; i < dim; ++i)
1240 u_s.rdof(i) = trans.dof(i);
1241 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1242 }
1243 }
1244 // Add the bending/twisting resistance energy contributions.
1245 if (kb) {
1246 DAVect tmp = sb_M_;
1247#ifdef THREED
1248 tmp.rx() = 0.0;
1249#endif
1250 energies_->estrain_ += 0.5*tmp.mag2() / kb;
1251 if (sb_BS_) {
1252 // accumulate bending slip energy.
1253 DAVect tmp_old = sb_M_old;
1254#ifdef THREED
1255 tmp_old.rx() = 0.0;
1256#endif
1257 DAVect avg_M = (tmp + tmp_old)*0.5;
1258 DAVect t_s_el = (tmp - tmp_old) / kb;
1259 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1260 }
1261 }
1262#ifdef THREED
1263 if (kt) {
1264 double mt = std::abs(sb_M_.x());
1265 energies_->estrain_ += 0.5*mt*mt / kt;
1266 if (sb_TS_) {
1267 // accumulate twisting slip energy.
1268 DAVect tmp(0.0);
1269 DAVect tmp_old(0.0);
1270 tmp.rx() = sb_M_.x();
1271 tmp_old.rx() = sb_M_old.x();
1272 DAVect avg_M = (tmp + tmp_old)*0.5;
1273 DAVect t_s_el = (tmp - tmp_old) / kt;
1274 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1275 }
1276 }
1277#endif
1278
1279 if (dpProps_) {
1280 // Calculate damping energy (accumulated) if the dashpots are active.
1281 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1282 }
1283 }
1284
1285 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1286 assert(sb_F_ == sb_F_);
1287 assert(sb_M_ == sb_M_);
1288 return true;
1289 }
1290
1291 bool ContactModelSoftBond::FDLawUnBonded(ContactModelMechanicalState *state, const double ×tep) {
1292
1293 // Relative translational/rotational displacement increments
1294 DVect trans = state->relativeTranslationalIncrement_;
1295 DAVect ang = state->relativeAngularIncrement_;
1296 double overlap = rgap_ - state->gap_;
1297 double correction = 1.0;
1298 if (state->activated() && sb_mode_ == 0 && trans.x()) {
1299 correction = -1.0*overlap / trans.x();
1300 if (correction < 0)
1301 correction = 1.0;
1302 }
1303
1304 // Store previous force and moment
1305 DVect sb_F_old = sb_F_;
1306 DAVect sb_M_old = sb_M_;
1307
1308 // Update stiffness data
1309 StiffData stiff = computeStiffData(state);
1310 double kn = stiff.trans_.x();
1311 double ks = stiff.trans_.y();
1312 double kb = stiff.ang_.z();
1313#if DIM==3
1314 double kt = stiff.ang_.x();
1315#endif
1316 // absolute/incremental normal force calculation
1317 if (sb_mode_==0)
1318 sb_F_.rx() = overlap * kn;
1319 else
1320 sb_F_.rx() -= trans.x() * kn;
1321 // Normal force can only be positive if unbonded
1322 sb_F_.rx() = std::max(0.0, sb_F_.x());
1323
1324 // Calculate the trial shear force.
1325 DVect sforce(0.0);
1326 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1327 // Loop over the shear components (note: the 0 component is the normal component)
1328 // and calculate the shear force.
1329 for (int i = 1; i<dim; ++i)
1330 sforce.rdof(i) = sb_F_.dof(i) - trans.dof(i) * ks;
1331
1332 // Calculate the trial moment.
1333 DAVect mom = sb_M_ - ang*stiff.ang_;
1334
1335 // If the SOLVE ELASTIC command is given then the
1336 // canFail state is set to FALSE. Otherwise it is always TRUE.
1337 if (state->canFail_) {
1338 bool changed = false;
1339 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
1340 bool slip_changed = false;
1341 double crit = sb_F_.x() * fric_;
1342 // The is the magnitude of the shear force.
1343 double sfmag = sforce.mag();
1344 // Sliding occurs when the magnitude of the shear force is greater than the
1345 // critical value.
1346 if (sfmag > crit) {
1347 // Lower the shear force to the critical value for sliding.
1348 double rat = crit / sfmag;
1349 sforce *= rat;
1350 if (!sb_S_) {
1351 slip_changed = true;
1352 changed = true;
1353 }
1354 sb_S_ = true;
1355 }
1356 else {
1357 if (sb_S_) {
1358 slip_changed = true;
1359 changed = true;
1360 }
1361 sb_S_ = false;
1362 }
1363
1364 // Resolve bending
1365 bool bslip_changed = false;
1366 crit = sb_bmul_ * 2.1*0.25*sb_F_.x() * stiff.reff_; // Jiang 2015
1367 DAVect test = mom;
1368#if DIM==3
1369 test.rx() = 0.0;
1370#endif
1371 double tmag = test.mag();
1372 if (tmag > crit) {
1373 // Lower the bending moment to the critical value for sliding.
1374 double rat = crit / tmag;
1375 test *= rat;
1376 if (!sb_BS_) {
1377 bslip_changed = true;
1378 changed = true;
1379 }
1380 sb_BS_ = true;
1381 }
1382 else {
1383 if (sb_BS_) {
1384 bslip_changed = true;
1385 changed = true;
1386 }
1387 sb_BS_ = false;
1388 }
1389 mom.rz() = test.z();
1390#if DIM==3
1391 mom.ry() = test.y();
1392 // Resolve twisting
1393 bool tslip_changed = false;
1394 crit = sb_tmul_ * 0.65*fric_*sb_F_.x() * stiff.reff_; // Jiang 2015
1395 tmag = std::abs(mom.x());
1396 if (tmag > crit) {
1397 // Lower the twisting moment to the critical value for sliding.
1398 double rat = crit / tmag;
1399 mom.rx() *= rat;
1400 if (!sb_TS_) {
1401 tslip_changed = true;
1402 changed = true;
1403 }
1404 sb_TS_ = true;
1405 } else {
1406 if (sb_TS_) {
1407 tslip_changed = true;
1408 changed = true;
1409 }
1410 sb_TS_ = false;
1411 }
1412#endif
1413 if (changed && cmEvents_[fSlipChange] >= 0) {
1414 qint64 code = 0;
1415 if (slip_changed) {
1416 code = 1;
1417 if (bslip_changed) {
1418 code = 4;
1419#if DIM==3
1420 if (tslip_changed)
1421 code = 7;
1422#endif
1423 }
1424 }
1425 else if (bslip_changed) {
1426 code = 2;
1427#if DIM==3
1428 if (tslip_changed)
1429 code = 6;
1430#endif
1431 }
1432#if DIM==3
1433 else if (tslip_changed) {
1434 code = 3;
1435 if (slip_changed)
1436 code = 5;
1437 }
1438#endif
1439 auto c = state->getContact();
1440 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1441 fish::Parameter(code),
1442 fish::Parameter(sb_S_),
1443 fish::Parameter(sb_BS_)
1444#ifdef THREED
1445 ,fish::Parameter(sb_TS_)
1446#endif
1447 };
1448 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1449 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1450 }
1451 }
1452
1453 // Set the shear components of the total force.
1454 for (int i = 1; i<dim; ++i)
1455 sb_F_.rdof(i) = sforce.dof(i);
1456
1457 // Set the moment.
1458 sb_M_ = mom;
1459
1460 // Account for dashpot forces if the dashpot structure has been defined.
1461 if (dpProps_) {
1462 dpProps_->dp_F_.fill(0.0);
1463 double vcn(0.0), vcs(0.0);
1464 // Calculate the damping coefficients.
1465 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1466 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1467 // First damp the shear components
1468 for (int i = 1; i<dim; ++i)
1469 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1470 // Damp the normal component
1471 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1472 // Need to change behavior based on the dp_mode.
1473 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1474 // Limit in tension if not bonded.
1475 if (dpProps_->dp_F_.x() + sb_F_.x() < 0)
1476 dpProps_->dp_F_.rx() = -sb_F_.rx();
1477 }
1478 if (sb_S_ && dpProps_->dp_mode_ > 1) {
1479 // Limit in shear if not sliding.
1480 double dfn = dpProps_->dp_F_.rx();
1481 dpProps_->dp_F_.fill(0.0);
1482 dpProps_->dp_F_.rx() = dfn;
1483 }
1484 }
1485
1486 //Compute energies if energy tracking has been enabled.
1487 if (state->trackEnergy_) {
1488 assert(energies_);
1489 energies_->estrain_ = 0.0;
1490 if (kn_)
1491 // Calculate the strain energy.
1492 energies_->estrain_ = 0.5*sb_F_.x()*sb_F_.x() / kn;
1493 if (ks_) {
1494 DVect s = sb_F_;
1495 s.rx() = 0.0;
1496 double smag2 = s.mag2();
1497 // Add the shear component of the strain energy.
1498 energies_->estrain_ += 0.5*smag2 / ks;
1499
1500 if (sb_S_) {
1501 // If sliding calculate the slip energy and accumulate it.
1502 sb_F_old.rx() = 0.0;
1503 DVect avg_F_s = (s + sb_F_old)*0.5;
1504 DVect u_s_el = (s - sb_F_old) / ks;
1505 DVect u_s(0.0);
1506 for (int i = 1; i < dim; ++i)
1507 u_s.rdof(i) = trans.dof(i);
1508 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1509 }
1510 }
1511 // Add the bending/twisting resistance energy contributions.
1512 if (kb) {
1513 DAVect tmp = sb_M_;
1514#ifdef THREED
1515 tmp.rx() = 0.0;
1516#endif
1517 energies_->estrain_ += 0.5*tmp.mag2() / kb;
1518 if (sb_BS_) {
1519 // accumulate bending slip energy.
1520 DAVect tmp_old = sb_M_old;
1521#ifdef THREED
1522 tmp_old.rx() = 0.0;
1523#endif
1524 DAVect avg_M = (tmp + tmp_old)*0.5;
1525 DAVect t_s_el = (tmp - tmp_old) / kb;
1526 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1527 }
1528 }
1529#ifdef THREED
1530 if (kt) {
1531 double mt = std::abs(sb_M_.x());
1532 energies_->estrain_ += 0.5*mt*mt / kt;
1533 if (sb_TS_) {
1534 // accumulate twisting slip energy.
1535 DAVect tmp(0.0);
1536 DAVect tmp_old(0.0);
1537 tmp.rx() = sb_M_.x();
1538 tmp_old.rx() = sb_M_old.x();
1539 DAVect avg_M = (tmp + tmp_old)*0.5;
1540 DAVect t_s_el = (tmp - tmp_old) / kt;
1541 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1542 }
1543 }
1544#endif
1545
1546 if (dpProps_) {
1547 // Calculate damping energy (accumulated) if the dashpots are active.
1548 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1549 }
1550 }
1551
1552 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1553 assert(sb_F_ == sb_F_);
1554 assert(sb_M_ == sb_M_);
1555 return true;
1556 }
1557
1558 bool ContactModelSoftBond::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
1559 // Account for thermal expansion in incremental mode
1560 if (sb_mode_ == 0 || ts->gapInc_ == 0.0) return false;
1561 DVect finc(0.0);
1562 finc.rx() = kn_ * ts->gapInc_;
1563 sb_F_ -= finc;
1564 return true;
1565 }
1566
1567 void ContactModelSoftBond::setForce(const DVect &v,IContact *c) {
1568 sb_F(v);
1569 if (v.x() > 0) {
1570 auto con = convert_getcast<IContactMechanical>(c);
1571 rgap_ = c->getGap() + v.x() / (kn_ * computeGeomData(con->getEnd1Curvature(),con->getEnd2Curvature()).x());
1572 }
1573 }
1574
1575 void ContactModelSoftBond::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
1576 // Only called for contacts with wall facets when the wall resolution scheme
1577 // is set to full!
1578 // Only do something if the contact model is of the same type
1579 if (old->getContactModel()->getName().compare("softbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
1580 ContactModelSoftBond *oldCm = (ContactModelSoftBond *)old;
1581#ifdef THREED
1582 // Need to rotate just the shear component from oldSystem to newSystem
1583
1584 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
1585 DVect axis = oldSystem.e1() & newSystem.e1();
1586 double c, ang, s;
1587 DVect re2;
1588 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1589 axis = axis.unit();
1590 c = oldSystem.e1()|newSystem.e1();
1591 if (c > 0)
1592 c = std::min(c,1.0);
1593 else
1594 c = std::max(c,-1.0);
1595 ang = acos(c);
1596 s = sin(ang);
1597 double t = 1. - c;
1598 DMatrix<3,3> rm;
1599 rm.get(0,0) = t*axis.x()*axis.x() + c;
1600 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
1601 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
1602 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
1603 rm.get(1,1) = t*axis.y()*axis.y() + c;
1604 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
1605 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
1606 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
1607 rm.get(2,2) = t*axis.z()*axis.z() + c;
1608 re2 = rm*oldSystem.e2();
1609 }
1610 else
1611 re2 = oldSystem.e2();
1612 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
1613 axis = re2 & newSystem.e2();
1614 DVect2 tpf;
1615 DVect2 tpm;
1616 DMatrix<2,2> m;
1617 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1618 axis = axis.unit();
1619 c = re2|newSystem.e2();
1620 if (c > 0)
1621 c = std::min(c,1.0);
1622 else
1623 c = std::max(c,-1.0);
1624 ang = acos(c);
1625 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
1626 ang *= -1;
1627 s = sin(ang);
1628 m.get(0,0) = c;
1629 m.get(1,0) = s;
1630 m.get(0,1) = -m.get(1,0);
1631 m.get(1,1) = m.get(0,0);
1632 tpf = m*DVect2(oldCm->sb_F_.y(),oldCm->sb_F_.z());
1633 tpm = m*DVect2(oldCm->sb_M_.y(),oldCm->sb_M_.z());
1634 } else {
1635 m.get(0,0) = 1.;
1636 m.get(0,1) = 0.;
1637 m.get(1,0) = 0.;
1638 m.get(1,1) = 1.;
1639 tpf = DVect2(oldCm->sb_F_.y(),oldCm->sb_F_.z());
1640 tpm = DVect2(oldCm->sb_M_.y(),oldCm->sb_M_.z());
1641 }
1642 DVect pforce = DVect(0,tpf.x(),tpf.y());
1643 //DVect pm = DVect(0,tpm.x(),tpm.y());
1644#else
1645 oldSystem;
1646 newSystem;
1647 DVect pforce = DVect(0,oldCm->sb_F_.y());
1648 //DVect pm = DVect(0,oldCm->sb_M_.y());
1649#endif
1650 for (int i=1; i<dim; ++i)
1651 sb_F_.rdof(i) += pforce.dof(i);
1652 if (sb_mode_ && oldCm->sb_mode_)
1653 sb_F_.rx() = oldCm->sb_F_.x();
1654 oldCm->sb_F_ = DVect(0.0);
1655 oldCm->sb_M_ = DAVect(0.0);
1656 if (dpProps_ && oldCm->dpProps_) {
1657#ifdef THREED
1658 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
1659 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
1660#else
1661 pforce = oldCm->dpProps_->dp_F_;
1662#endif
1663 dpProps_->dp_F_ += pforce;
1664 oldCm->dpProps_->dp_F_ = DVect(0.0);
1665 }
1666 if(oldCm->getEnergyActivated()) {
1667 activateEnergy();
1668 energies_->estrain_ = oldCm->energies_->estrain_;
1669 energies_->edashpot_ = oldCm->energies_->edashpot_;
1670 energies_->eslip_ = oldCm->energies_->eslip_;
1671 oldCm->energies_->estrain_ = 0.0;
1672 oldCm->energies_->edashpot_ = 0.0;
1673 oldCm->energies_->eslip_ = 0.0;
1674 }
1675 rgap_ = oldCm->rgap_;
1676 }
1677 assert(sb_F_ == sb_F_);
1678 }
1679
1680 void ContactModelSoftBond::setNonForcePropsFrom(IContactModel *old) {
1681 // Only called for contacts with wall facets when the wall resolution scheme
1682 // is set to full!
1683 // Only do something if the contact model is of the same type
1684 if (old->getName().compare("softbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
1685 ContactModelSoftBond *oldCm = (ContactModelSoftBond *)old;
1686 kn_ = oldCm->kn_;
1687 ks_ = oldCm->ks_;
1688 fric_ = oldCm->fric_;
1689 sb_bmul_ = oldCm->sb_bmul_;
1690 sb_tmul_ = oldCm->sb_tmul_;
1691 sb_mode_ = oldCm->sb_mode_;
1692 sb_rmul_ = oldCm->sb_rmul_;
1693 sb_S_ = oldCm->sb_S_;
1694 sb_BS_ = oldCm->sb_BS_;
1695 sb_TS_ = oldCm->sb_TS_;
1696 rgap_ = oldCm->rgap_;
1697 userArea_ = oldCm->userArea_;
1698
1699 if (oldCm->dpProps_) {
1700 if (!dpProps_)
1701 dpProps_ = NEWC(dpProps());
1702 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1703 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1704 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1705 }
1706 if (oldCm->bProps_) {
1707 if (!bProps_)
1708 bProps_ = NEWC(bProps());
1709 bProps_->sb_mcf_ = oldCm->bProps_->sb_mcf_;
1710 bProps_->sb_fa_ = oldCm->bProps_->sb_fa_;
1711 bProps_->sb_state_ = oldCm->bProps_->sb_state_;
1712 bProps_->sb_coh_ = oldCm->bProps_->sb_coh_;
1713 bProps_->sb_ten_ = oldCm->bProps_->sb_ten_;
1714 bProps_->sb_maxTen_ = oldCm->bProps_->sb_maxTen_;
1715 bProps_->sb_cut_ = oldCm->bProps_->sb_cut_;
1716 bProps_->sb_delu_ = oldCm->bProps_->sb_delu_;
1717 bProps_->sb_delo_ = oldCm->bProps_->sb_delo_;
1718 bProps_->sb_maxu_ = oldCm->bProps_->sb_maxu_;
1719 bProps_->sb_critu_ = oldCm->bProps_->sb_critu_;
1720 }
1721
1722 }
1723 }
1724
1725 DVect ContactModelSoftBond::getForce() const {
1726 DVect ret(sb_F_);
1727 if (dpProps_)
1728 ret += dpProps_->dp_F_;
1729 return ret;
1730 }
1731
1732 DAVect ContactModelSoftBond::getMomentOn1(const IContactMechanical *c) const {
1733 DAVect ret(sb_M_);
1734 if (c) {
1735 DVect force = getForce();
1736 c->updateResultingTorqueOn1Local(force,&ret);
1737 }
1738 return ret;
1739 }
1740
1741 DAVect ContactModelSoftBond::getMomentOn2(const IContactMechanical *c) const {
1742 DAVect ret(sb_M_);
1743 if (c) {
1744 DVect force = getForce();
1745 c->updateResultingTorqueOn2Local(force,&ret);
1746 }
1747 return ret;
1748 }
1749
1750 DAVect ContactModelSoftBond::getModelMomentOn1() const {
1751 DAVect ret(sb_M_);
1752 return ret;
1753 }
1754
1755 DAVect ContactModelSoftBond::getModelMomentOn2() const {
1756 DAVect ret(sb_M_);
1757 return ret;
1758 }
1759
1760 void ContactModelSoftBond::objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *ret) const {
1761 ret->clear();
1762 ret->push_back({"kn",ScalarInfo});
1763 ret->push_back({"ks",ScalarInfo});
1764 ret->push_back({"fric",ScalarInfo});
1765 ret->push_back({"sb_bmul",ScalarInfo});
1766 ret->push_back({"sb_tmul",ScalarInfo});
1767 ret->push_back({"sb_mode",ScalarInfo});
1768 ret->push_back({"sb_force",VectorInfo});
1769 ret->push_back({"sb_moment",VectorInfo});
1770 ret->push_back({"sb_slip",ScalarInfo});
1771 ret->push_back({"sb_slipb",ScalarInfo});
1772 ret->push_back({"sb_slipt",ScalarInfo});
1773 ret->push_back({"sb_rmul",ScalarInfo});
1774 ret->push_back({"sb_radius",ScalarInfo});
1775 ret->push_back({"emod",ScalarInfo});
1776 ret->push_back({"kratio",ScalarInfo});
1777 ret->push_back({"dp_nratio",ScalarInfo});
1778 ret->push_back({"dp_sratio",ScalarInfo});
1779 ret->push_back({"dp_mode",ScalarInfo});
1780 ret->push_back({"dp_force",VectorInfo});
1781 ret->push_back({"sb_state",ScalarInfo});
1782 ret->push_back({"sb_ten",ScalarInfo});
1783 ret->push_back({"sb_shear",ScalarInfo});
1784 ret->push_back({"sb_coh",ScalarInfo});
1785 ret->push_back({"sb_fa",ScalarInfo});
1786 ret->push_back({"sb_mcf",ScalarInfo});
1787 ret->push_back({"sb_sigma",ScalarInfo});
1788 ret->push_back({"sb_tau",ScalarInfo});
1789 ret->push_back({"sb_soft",ScalarInfo});
1790 ret->push_back({"sb_cut",ScalarInfo});
1791 ret->push_back({"sb_area",ScalarInfo});
1792 ret->push_back({"user_area",ScalarInfo});
1793 ret->push_back({"rgap",ScalarInfo});
1794
1795 }
1796
1797 void ContactModelSoftBond::objectPropValues(std::vector<double> *ret,const IContact *con) const {
1798 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
1799 ret->clear();
1800 ret->push_back(kn_);
1801 ret->push_back(ks_);
1802 ret->push_back(fric_);
1803 ret->push_back(sb_bmul_);
1804 ret->push_back(sb_tmul_);
1805 ret->push_back(sb_mode_);
1806 ret->push_back(sb_F_.mag());
1807 ret->push_back(sb_M_.mag());
1808 ret->push_back(sb_S_);
1809 ret->push_back(sb_BS_);
1810 ret->push_back(sb_TS_);
1811 ret->push_back(sb_rmul_);
1812 double d;
1813 double Cmax1 = c->getEnd1Curvature().y();
1814 double Cmax2 = c->getEnd2Curvature().y();
1815 if (!userArea_)
1816 d= sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
1817 else {
1818#ifdef THREED
1819 d = std::sqrt(userArea_ / dPi);
1820#else
1821 d = userArea_ / 2.0;
1822#endif
1823 }
1824 ret->push_back(d);
1825 double rsum(0.0);
1826 if (c->getEnd1Curvature().y())
1827 rsum += 1.0/c->getEnd1Curvature().y();
1828 if (c->getEnd2Curvature().y())
1829 rsum += 1.0/c->getEnd2Curvature().y();
1830 if (userArea_)
1831#ifdef THREED
1832 rsum = 2.0 * std::sqrt(userArea_ / dPi);
1833#else
1834 rsum = userArea_;
1835#endif
1836 d= kn_ * rsum;
1837 ret->push_back(d);
1838 ret->push_back((ks_ == 0.0) ? 0.0 : (kn_/ks_));
1839 ret->push_back(dpProps_ ? dpProps_->dp_nratio_ : 0);
1840 ret->push_back(dpProps_ ? dpProps_->dp_sratio_ : 0);
1841 ret->push_back(dpProps_ ? dpProps_->dp_mode_ : 0);
1842 ret->push_back(dpProps_ ? dpProps_->dp_F_.mag() : 0);
1843 ret->push_back(bProps_ ? bProps_->sb_state_ : 0);
1844 ret->push_back(bProps_ ? bProps_->sb_ten_ : 0);
1845 ret->push_back(shearStrength(computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x()));
1846 ret->push_back(bProps_ ? bProps_->sb_coh_ : 0);
1847 ret->push_back(bProps_ ? bProps_->sb_fa_ : 0);
1848 ret->push_back(bProps_ ? bProps_->sb_mcf_ : 0);
1849 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x());
1850 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).y());
1851 ret->push_back(bProps_ ? bProps_->sb_soft_ : 0);
1852 ret->push_back(bProps_ ? bProps_->sb_cut_ : 0);
1853 ret->push_back(userArea_ ? userArea_ : computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x());
1854 ret->push_back(userArea_);
1855 ret->push_back(rgap_);
1856 }
1857
1858 DVect3 ContactModelSoftBond::computeGeomData(const DVect2 &end1Curvature,const DVect2 &end2Curvature) const {
1859 double Cmax1 = end1Curvature.y();
1860 double Cmax2 = end2Curvature.y();
1861 double br = sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
1862 if (userArea_)
1863#ifdef THREED
1864 br = std::sqrt(userArea_ / dPi);
1865#else
1866 br = userArea_ / 2.0;
1867#endif
1868 double br2 = br * br;
1869#ifdef TWOD
1870 double area = 2.0*br;
1871 double bi = 2.0*br*br2 / 3.0;
1872#else
1873 double area = dPi * br2;
1874 double bi = 0.25*area*br2;
1875#endif
1876 return DVect3(area, bi, br);
1877 }
1878
1879 DVect2 ContactModelSoftBond::SMax(const DVect2 &end1Curvature,const DVect2 &end2Curvature) const {
1880 DVect3 data = computeGeomData(end1Curvature,end2Curvature);
1881 double area = data.x();
1882 double bi = data.y();
1883 double br = data.z();
1884 /* maximum stresses */
1885 double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1886 double dtwist = sb_M_.x();
1887 DVect bfs(sb_F_);
1888 bfs.rx() = 0.0;
1889 double dbfs = bfs.mag();
1890 double nsmax = -(sb_F_.x() / area) + dbend * br / bi;
1891 double ssmax = dbfs / area + std::abs(dtwist) * 0.5* br / bi;
1892 return DVect2(nsmax, ssmax);
1893 }
1894
1895 double ContactModelSoftBond::shearStrength(const double &area) const {
1896 if (!bProps_) return 0.0;
1897 double sig = -1.0*sb_F_.x() / area;
1898 double nstr = bProps_->sb_state_ > 2 ? bProps_->sb_ten_ : 0.0;
1899 return sig <= nstr ? bProps_->sb_coh_ - std::tan(dDegrad*bProps_->sb_fa_)*sig
1900 : bProps_->sb_coh_ - std::tan(dDegrad*bProps_->sb_fa_)*nstr;
1901 }
1902
1903
1904 double ContactModelSoftBond::strainEnergy(double kn,double ks,double kb,double kt) const {
1905 double ret(0.0);
1906 if (kn)
1907 ret = 0.5 * sb_F_.x() * sb_F_.x() / kn;
1908 if (ks) {
1909 DVect tmp = sb_F_;
1910 tmp.rx() = 0.0;
1911 double smag2 = tmp.mag2();
1912 ret += 0.5 * smag2 / ks;
1913 }
1914
1915 if (kt)
1916 ret += 0.5 * sb_M_.x() * sb_M_.x() / kt;
1917 if (kb) {
1918 DAVect tmp = sb_M_;
1919#ifdef THREED
1920 tmp.rx() = 0.0;
1921 double smag2 = tmp.mag2();
1922#else
1923 double smag2 = tmp.z() * tmp.z();
1924#endif
1925 ret += 0.5 * smag2 / kb;
1926 }
1927 return ret;
1928 }
1929
1930} // namespace cmodelsxd
1931// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Dec 14, 2024 |
