Spring Network Model Implementation
See this page for the documentation of this contact model.
contactmodelrbsn.h
1#pragma once
2// contactmodelrbsn.h
3
4#include "contactmodel/src/contactmodelmechanical.h"
5
6#ifdef RBSN_LIB
7# define RBSN_EXPORT EXPORT_TAG
8#elif defined(NO_MODEL_IMPORT)
9# define RBSN_EXPORT
10#else
11# define RBSN_EXPORT IMPORT_TAG
12#endif
13
14namespace cmodelsxd {
15 using namespace itasca;
16
17 class ContactModelRBSN : public ContactModelMechanical {
18 public:
19 // Constructor: Set default values for contact model properties.
20 RBSN_EXPORT ContactModelRBSN();
21 RBSN_EXPORT ContactModelRBSN(const ContactModelRBSN &) noexcept;
22 RBSN_EXPORT const ContactModelRBSN & operator=(const ContactModelRBSN &);
23 RBSN_EXPORT void addToStorage(poly::vector<ContactModelMechanical> *s,int ww = -1) override;
24 // Destructor, called when contact is deleted: free allocated memory, etc.
25 RBSN_EXPORT virtual ~ContactModelRBSN();
26 // Contact model name (used as keyword for commands and FISH).
27 QString getName() const override { return "springnetwork"; }
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 , kwKrot
48 , kwFric
49 , kwBMul
50 , kwTMul
51 , kwSNMode
52 , kwSNF
53 , kwSNM
54 , kwSNS
55 , kwSNBS
56 , kwSNTS
57 , kwSNRMul
58 , kwSNRadius
59 , kwEmod
60 , kwKRatio
61 , kwDpNRatio
62 , kwDpSRatio
63 , kwDpMode
64 , kwDpF
65 , kwSNState
66 , kwSNTStr
67 , kwSNSStr
68 , kwSNCoh
69 , kwSNFa
70 , kwSNMCF
71 , kwSNSig
72 , kwSNTau
73 , kwSNArea
74 , kwUserArea
75 , kwRGap
76 , kwPForce
77 , kwPois
78 , kwPoisDiag
79 , kwSnCohRes
80 , kwSnDil
81 , kwSnDilZ
82 , kwSnNormDisp
83 , kwSnShearDisp
84 , kwSnCohDc
85 , kwSnHeal
86 , kwSnTenDc
87 , kwTenTable
88 , kwCohTable
89 , kwTablePos
90 , kwPorP
91 , kwStressNorm
92 };
93 // Contact model property names in a comma separated list. The order corresponds with
94 // the order of the PropertyKeys enumerator above. One can visualize any of these
95 // properties in PFC automatically.
96 QString getProperties() const override {
97 return "kn"
98 ",ks"
99 ",krot"
100 ",fric"
101 ",sn_bmul"
102 ",sn_tmul"
103 ",sn_mode"
104 ",sn_force"
105 ",sn_moment"
106 ",sn_slip"
107 ",sn_slipb"
108 ",sn_slipt"
109 ",sn_rmul"
110 ",sn_radius"
111 ",emod"
112 ",kratio"
113 ",dp_nratio"
114 ",dp_sratio"
115 ",dp_mode"
116 ",dp_force"
117 ",sn_state"
118 ",sn_ten"
119 ",sn_shear"
120 ",sn_coh"
121 ",sn_fa"
122 ",sn_mcf"
123 ",sn_sigma"
124 ",sn_tau"
125 ",sn_area"
126 ",user_area"
127 ",rgap"
128 ",sn_pois_force"
129 ",sn_pois"
130 ",sn_non_diag"
131 ",sn_cohres"
132 ",sn_dil"
133 ",sn_dilzero"
134 ",sn_ndisp"
135 ",sn_sdisp"
136 ",sn_cohdc"
137 ",sn_heal"
138 ",sn_tendc"
139 ",sn_tentab"
140 ",sn_cohtab"
141 ",sn_tabpos"
142 ",sn_porp"
143 ",sn_esigma"
144 ;
145 }
146 // Enumerator for the energies.
147 enum EnergyKeys {
148 kwEStrain=1
149 , kwESlip
150 , kwEDashpot
151 };
152 // Contact model energy names in a comma separated list. The order corresponds with
153 // the order of the EnergyKeys enumerator above.
154 QString getEnergies() const override {
155 return "energy-strain"
156 ",energy-slip"
157 ",energy-dashpot";
158 }
159 // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
160 double getEnergy(uint32 i) const override;
161 // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1)
162 // returns wther or not the estrain energy is accumulated which is false).
163 bool getEnergyAccumulate(uint32 i) const override;
164 // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
165 void setEnergy(uint32 i,const double &d) override; // Base 1
166 // Activate the energy. This is only called if the energy tracking is enabled.
167 void activateEnergy() override { if (energies_) return; energies_ = NEWC(Energies());}
168 // Returns whether or not the energy tracking has been enabled for this contact.
169 bool getEnergyActivated() const override {return (energies_ != 0);}
170
171 // Enumerator for contact model related FISH callback events.
172 enum FishCallEvents {
173 fActivated=0
174 ,fSlipChange
175 ,fBondBreak
176 };
177 // Contact model FISH callback event names in a comma separated list. The order corresponds with
178 // the order of the FishCallEvents enumerator above.
179 QString getFishCallEvents() const override {
180 return
181 "contact_activated"
182 ",slip_change"
183 ",bond_break";
184 }
185
186 // Return the specified contact model property.
187 QVariant getProperty(uint32 i,const IContact *) const override;
188 // The return value denotes whether or not the property corresponds to the global
189 // or local coordinate system (TRUE: global system, FALSE: local system). The
190 // local system is the contact-plane system (nst) defined as follows.
191 // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
192 // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
193 // and tc are unit vectors in directions of the nst axes.
194 // This is used when rendering contact model properties that are vectors.
195 bool getPropertyGlobal(uint32 i) const override;
196 // Set the specified contact model property, ensuring that it is of the correct type
197 // and within the correct range --- if not, then throw an exception.
198 // The return value denotes whether or not the update has affected the timestep
199 // computation (by having modified the translational or rotational tangent stiffnesses).
200 // If true is returned, then the timestep will be recomputed.
201 bool setProperty(uint32 i,const QVariant &v,IContact *) override;
202 // The return value denotes whether or not the property is read-only
203 // (TRUE: read-only, FALSE: read-write).
204 bool getPropertyReadOnly(uint32 i) const override;
205
206 // The return value denotes whether or not the property is inheritable
207 // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
208 // the endPropertyUpdated method.
209 bool supportsInheritance(uint32 i) const override;
210 // Return whether or not inheritance is enabled for the specified property.
211 bool getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i); return (inheritanceField_ & mask) ? true : false; }
212 // Set the inheritance flag for the specified property.
213 void setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i); if (b) inheritanceField_ |= mask; else inheritanceField_ &= ~mask; }
214
215 // Enumerator for contact model methods.
216 enum MethodKeys { kwAssignStiffness=1, kwStiffness, kwBond, kwUnbond, kwArea, kwResetDisp};
217 // Contact model methoid names in a comma separated list. The order corresponds with
218 // the order of the MethodKeys enumerator above.
219 QString getMethods() const override { return "assign-stiffness,compute-stiffness,bond,unbond,area,reset-disp";}
220 // Return a comma seprated list of the contact model method arguments (base 1).
221 QString getMethodArguments(uint32 i) const override;
222 // Set contact model method arguments (base 1).
223 // The return value denotes whether or not the update has affected the timestep
224 // computation (by having modified the translational or rotational tangent stiffnesses).
225 // If true is returned, then the timestep will be recomputed.
226 bool setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override;
227
228 // Prepare for entry into ForceDispLaw. The validate function is called when:
229 // (1) the contact is created, (2) a property of the contact that returns a true via
230 // the setProperty method has been modified and (3) when a set of cycles is executed
231 // via the {cycle N} command.
232 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
233 bool validate(ContactModelMechanicalState *state,const double ×tep) override;
234 // The endPropertyUpdated method is called whenever a surface property (with a name
235 // that matches an inheritable contact model property name) of one of the contacting
236 // pieces is modified. This allows the contact model to update its associated
237 // properties. The return value denotes whether or not the update has affected
238 // the time step computation (by having modified the translational or rotational
239 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
240 bool endPropertyUpdated(const QString &name,const IContactMechanical *c) override;
241 // The forceDisplacementLaw function is called during each cycle. Given the relative
242 // motion of the two contacting pieces (via
243 // state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
244 // state->relativeAngularIncrement_ (Dtt, Dtbs, Dtbt)
245 // Ddn : relative normal-displacement increment, Ddn > 0 is opening
246 // Ddss : relative shear-displacement increment (s-axis component)
247 // Ddst : relative shear-displacement increment (t-axis component)
248 // Dtt : relative twist-rotation increment
249 // Dtbs : relative bend-rotation increment (s-axis component)
250 // Dtbt : relative bend-rotation increment (t-axis component)
251 // The relative displacement and rotation increments:
252 // Dd = Ddn*nc + Ddss*sc + Ddst*tc
253 // Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
254 // where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
255 // [see {Table 1: Contact State Variables} in PFC Model Components:
256 // Contacts and Contact Models: Contact Resolution]
257 // ) and the contact properties, this function must update the contact force and
258 // moment.
259 // The force_ is acting on piece 2, and is expressed in the local coordinate system
260 // (defined in getPropertyGlobal) such that the first component positive denotes
261 // compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
262 // in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to
263 // get the total moment.
264 // The return value indicates the contact activity status (TRUE: active, FALSE:
265 // inactive) during the next cycle.
266 // Additional information:
267 // * If state->activated() is true, then the contact has just become active (it was
268 // inactive during the previous time step).
269 // * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
270 // the forceDispLaw handle the case of {state->canFail_ == true}.
271 bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) override;
272 bool thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
273
274 // The getEffectiveXStiffness functions return the translational and rotational
275 // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
276 // the translational tangent shear stiffness is zero (but this stiffness reduction
277 // is typically ignored when computing a stable time step). If the contact model
278 // includes a dashpot, then the translational stiffnesses must be increased (see
279 // Potyondy (2009)).
280 // [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
281 // Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
282 // December 7, 2009.]
283 DVect2 getEffectiveTranslationalStiffness() const override { return effectiveTranslationalStiffness_; }
284 DAVect getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness_;}
285
286 // Return a new instance of the contact model. This is used in the CMAT
287 // when a new contact is created.
288 ContactModelRBSN *clone() const override { return NEWC(ContactModelRBSN()); }
289 // The getActivityDistance function is called by the contact-resolution logic when
290 // the CMAT is modified. Return value is the activity distance used by the
291 // checkActivity function.
292 double getActivityDistance() const override {return rgap_;}
293 // The isOKToDelete function is called by the contact-resolution logic when...
294 // Return value indicates whether or not the contact may be deleted.
295 // If TRUE, then the contact may be deleted when it is inactive.
296 // If FALSE, then the contact may not be deleted (under any condition).
297 bool isOKToDelete() const override { return !isBonded(); }
298 // Zero the forces and moments stored in the contact model. This function is called
299 // when the contact becomes inactive.
300 void resetForcesAndMoments() override {
301 sn_F_ = DVect(0.0);
302 fictForce_ = DVect(0.0);
303 sn_M_ = DAVect(0.0);
304 dp_F(DVect(0.0));
305 if (energies_) {
306 energies_->estrain_ = 0.0;
307 }
308 }
309 void setForce(const DVect &v,IContact *c) override;
310 void setArea(const double &d) override { userArea_ = d; }
311 double getArea() const override { return userArea_; }
312
313 // The checkActivity function is called by the contact-resolution logic when...
314 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
315 bool checkActivity(const double &gap) override { return gap <= rgap_ || isBonded();}
316
317 // Returns the sliding state (FALSE is returned if not implemented).
318 bool isSliding() const override { return sn_S_; }
319 // Returns the bonding state (FALSE is returned if not implemented).
320 bool isBonded() const override { return sn_state_ >= 3; }
321 void unbond() override { sn_state_ = 0; }
322
323 // Both of these methods are called only for contacts with facets where the wall
324 // resolution scheme is set the full. In such cases one might wish to propagate
325 // contact state information (e.g., shear force) from one active contact to another.
326 // See the Faceted Wall section in the documentation.
327 void propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
328 void setNonForcePropsFrom(IContactModel *oldCM) override;
329
330 /// Return the total force that the contact model holds.
331 DVect getForce() const override;
332
333 /// Return the total moment on 1 that the contact model holds
334 DAVect getMomentOn1(const IContactMechanical *) const override;
335
336 /// Return the total moment on 1 that the contact model holds
337 DAVect getMomentOn2(const IContactMechanical *) const override;
338
339 DAVect getModelMomentOn1() const override;
340 DAVect getModelMomentOn2() const override;
341
342 // Used to efficiently get properties from the contact model for the object pane.
343 // List of properties for the object pane, comma separated.
344 // All properties will be cast to doubles for comparison. No other comparisons
345 // are supported. This may not be the same as the entire property list.
346 // Return property name and type for plotting.
347 void objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *) const override;
348 // All properties cast to doubles - this is what can be compared.
349 void objectPropValues(std::vector<double> *,const IContact *) const override;
350
351 // Methods to get and set properties.
352 double sn_Ten() const { return tenTable_[0].x(); }
353 void sn_Ten(const double &d) { tenTable_[0].rx() = d; }
354 double sn_Coh() const { return cohTable_[0].x(); }
355 void sn_Coh(const double &d) { cohTable_[0].rx() = d; }
356 void sn_MCF(const double &d) { sn_mcf_=d;}
357 double sn_cohdc() const {return cohTable_.back().y(); }
358 double sn_tendc() const {return tenTable_.back().y(); }
359
360 bool hasDamping() const {return dpProps_ ? true : false;}
361 double dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
362 void dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
363 double dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
364 void dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
365 int dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
366 void dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
367 DVect dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
368 void dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
369
370 bool hasEnergies() const {return energies_ ? true:false;}
371 double estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
372 void estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
373 double eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
374 void eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
375 double edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
376 void edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
377
378 uint32 inheritanceField() const {return inheritanceField_;}
379 void inheritanceField(uint32 i) {inheritanceField_ = i;}
380
381 const DVect2 & effectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
382 void effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
383 const DAVect & effectiveRotationalStiffness() const {return effectiveRotationalStiffness_;}
384 void effectiveRotationalStiffness(const DAVect &v ) {effectiveRotationalStiffness_=v;}
385
386 private:
387 // Index - used internally by PFC. Should be set to -1 in the cpp file.
388 static int index_;
389
390 bool FDLawBonded(ContactModelMechanicalState *state, const double ×tep);
391 bool FDLawUnBonded(ContactModelMechanicalState *state, const double ×tep);
392
393 // Structure to compute stiffness
394 struct StiffData {
395 DVect2 trans_ = DVect2(0.0);
396 DAVect ang_ = DAVect(0.0);
397 double reff_ = 0.0;
398 };
399
400 // Structure to store the energies.
401 struct Energies {
402 double estrain_ = 0.0; // elastic energy
403 double eslip_ = 0.0; // work dissipated by friction
404 double edashpot_ = 0.0; // work dissipated by dashpots
405 };
406
407 // Structure to store dashpot quantities.
408 struct dpProps {
409 double dp_nratio_ = 0.0; // normal viscous critical damping ratio
410 double dp_sratio_ = 0.0; // shear viscous critical damping ratio
411 int dp_mode_ = 0; // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
412 DVect dp_F_ = DVect(0.0); // Force in the dashpots
413 };
414
415 bool updateKn(const IContactMechanical *con);
416 bool updateKs(const IContactMechanical *con);
417 bool updateFric(const IContactMechanical *con);
418
419 StiffData computeStiffData(ContactModelMechanicalState *state) const;
420 DVect3 computeGeomData(const DVect2 &e1c,const DVect2 &e2c) const;
421 DVect2 SMax(const DVect2 &e1c,const DVect2 &e2c) const; // Maximum stress (tensile,shear) at bond periphery
422 double shearStrength(const double &pbArea) const; // Bond shear strength
423 double strainEnergy(double kn, double ks, double kb, double kt) const;
424
425 void updateStiffness(ContactModelMechanicalState *state);
426
427 // Contact model inheritance fields.
428 uint32 inheritanceField_;
429
430 // Effective translational stiffness.
431 DVect2 effectiveTranslationalStiffness_ = DVect2(0.0);
432 DAVect effectiveRotationalStiffness_ = DAVect(0.0); // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)
433
434 // linear model properties
435 DVect fictForce_ = DVect(0.0);// Ficticous force to be added
436 DVect sn_F_ = DVect(0.0); // Force carried in the model
437 DVect2 sn_sdisp_ = DVect2(0.0); // Accumulated total shear displacement
438 // The x component holds the current slip
439 DAVect sn_M_ = DAVect(0.0); // moment (bending + twisting in 3D)
440 DAVect kRot_ = DAVect(0.0); // Translational degrees of freedom
441 double kTran_ = 0.0; // Translational degrees of freedom
442 double kRatio_ = 1.0; // Ratio of normal to shear stiffness
443 double E_ = 0.0; // Young's modulus
444 double poisson_ = 0.0; // Poisson ratio
445 double fric_ = 0.0; // Coulomb friction coefficient
446 double sn_bmul_ = 1.0; // Bending friction multiplier
447 double sn_tmul_ = 1.0; // Twisting friction multiplier
448 double sn_rmul_ = 1.0; // Radius multiplier
449 double userArea_ = 0.0; // Area as specified by the user
450 double rgap_ = 0.0; // Reference gap
451 double sn_fa_ = 0.0; // friction angle (stored as tan(dDegrad*fa))
452 double sn_mcf_ = 1.0; // moment contribution factor
453 double sn_dil_ = 0.0; // Dilation (stored as tan(dDegrad*dil))
454 double sn_dilzero_ = 0.0; // Dilation zero
455 double transTen_ = 0.0; // Force for transition from tensile to compression
456 double sn_elong_ = 0.0; // Elongation (or normal displacement since softening)
457 double sn_ndisp_ = 0.0; // Accumulated normal displacement
458 double sn_cohres_ = 0.0; // Residual cohesion
459 double sn_por_ = 0.0; // Pore Pressure
460 uint sn_mode_ = 0; // specifies absolute (0) or incremental (1) behavior for the the normal force
461 uint sn_state_ = 0; // bond mode - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (B)
462 int sn_tabPos_ = 0; // Position in the table for query
463 bool poisOffDiag_ = false; // Add the off diagonal terms
464 bool sn_S_ = false; // The current slip state
465 bool sn_BS_ = false; // The bending slip state
466 bool sn_TS_ = false; // The twisting slip state
467 bool forceSet_ = false; // About setting the force
468 bool sn_heal_ = false; // Healing behavior
469
470 std::vector<DVect2> tenTable_ = { DVect2(0) };
471 std::vector<DVect2> cohTable_ = { DVect2(0) };
472
473 dpProps * dpProps_ = nullptr; // The viscous properties
474
475 Energies * energies_ = nullptr; // The energies
476
477 };
478} // namespace cmodelsxd
479// EoF
contactmodelrbsn.cpp
1// contactmodelrbsn.cpp
2#include "contactmodelrbsn.h"
3
4#include "module/interface/icontactmechanical.h"
5#include "module/interface/icontact.h"
6#include "module/interface/ipiece.h"
7#include "module/interface/ibody.h"
8#include "module/interface/ifishcalllist.h"
9
10#include "utility/src/tptr.h"
11#include "shared/src/mathutil.h"
12
13#include "kernel/interface/iprogram.h"
14#include "module/interface/icontactthermal.h"
15#include "contactmodel/src/contactmodelthermal.h"
16#include "../version.txt"
17#include "fish/src/parameter.h"
18
19#ifdef RBSN_LIB
20 int __stdcall DllMain(void *,unsigned, void *) {
21 return 1;
22 }
23
24 extern "C" EXPORT_TAG const char *getName() {
25#if DIM==3
26 return "contactmodelmechanical3drbsn";
27#else
28 return "contactmodelmechanical2drbsn";
29#endif
30 }
31
32 extern "C" EXPORT_TAG unsigned getMajorVersion() {
33 return MAJOR_VERSION;
34 }
35
36 extern "C" EXPORT_TAG unsigned getMinorVersion() {
37 return MINOR_VERSION;
38 }
39
40 extern "C" EXPORT_TAG void *createInstance() {
41 cmodelsxd::ContactModelRBSN *m = NEWC(cmodelsxd::ContactModelRBSN());
42 return (void *)m;
43 }
44#endif
45
46namespace cmodelsxd {
47 static const uint32 KnMask = 0x00000002; // Base 1!
48 static const uint32 KsMask = 0x00000004;
49 static const uint32 FricMask = 0x00000008;
50
51 using namespace itasca;
52
53 int ContactModelRBSN::index_ = -1;
54 uint32 ContactModelRBSN::getMinorVersion() const { return MINOR_VERSION; }
55
56 ContactModelRBSN::ContactModelRBSN() : inheritanceField_(KnMask|KsMask|FricMask) {
57 }
58
59 ContactModelRBSN::ContactModelRBSN(const ContactModelRBSN &m) noexcept {
60 inheritanceField(KnMask|KsMask|FricMask);
61 this->copy(&m);
62 }
63
64 const ContactModelRBSN & ContactModelRBSN::operator=(const ContactModelRBSN &m) {
65 inheritanceField(KnMask|KsMask|FricMask);
66 this->copy(&m);
67 return *this;
68 }
69
70 void ContactModelRBSN::addToStorage(poly::vector<ContactModelMechanical> *s,int ww) {
71 s->addToStorage<ContactModelRBSN>(*this,ww);
72 }
73
74 ContactModelRBSN::~ContactModelRBSN() {
75 // Make sure to clean up after yourself!
76 if (dpProps_)
77 delete dpProps_;
78 if (energies_)
79 delete energies_;
80 }
81
82 void ContactModelRBSN::archive(ArchiveStream &stream) {
83 if (stream.getRestoreVersion() > 4) {
84 // New version
85 stream & fictForce_;
86 stream & sn_F_;
87 stream & sn_sdisp_;
88 stream & sn_M_;
89 stream & kRot_;
90 stream & kTran_;
91 stream & kRatio_;
92 stream & E_;
93 stream & poisson_;
94 stream & fric_;
95 stream & sn_bmul_;
96 stream & sn_tmul_;
97 stream & sn_rmul_;
98 stream & userArea_;
99 stream & rgap_;
100 stream & sn_fa_;
101 stream & sn_mcf_;
102 stream & sn_dil_;
103 stream & sn_dilzero_;
104 stream & transTen_;
105 stream & sn_elong_;
106 stream & sn_ndisp_;
107 stream & sn_mode_;
108 stream & sn_state_;
109 stream & poisOffDiag_;
110 stream & sn_S_;
111 stream & sn_BS_;
112 stream & sn_TS_;
113 stream & forceSet_;
114 stream & sn_heal_;
115 stream & tenTable_;
116 stream & cohTable_;
117 stream & inheritanceField_;
118 stream & effectiveTranslationalStiffness_;
119 stream & effectiveRotationalStiffness_;
120 if (stream.getArchiveState()==ArchiveStream::Save) {
121 bool b = false;
122 if (dpProps_) {
123 b = true;
124 stream & b;
125 stream & dpProps_->dp_nratio_;
126 stream & dpProps_->dp_sratio_;
127 stream & dpProps_->dp_mode_;
128 stream & dpProps_->dp_F_;
129 }
130 else
131 stream & b;
132 b = false;
133 if (energies_) {
134 b = true;
135 stream & b;
136 stream & energies_->estrain_;
137 stream & energies_->eslip_;
138 stream & energies_->edashpot_;
139 }
140 else
141 stream & b;
142 } else {
143 bool b(false);
144 stream & b;
145 if (b) {
146 if (!dpProps_)
147 dpProps_ = NEWC(dpProps());
148 stream & dpProps_->dp_nratio_;
149 stream & dpProps_->dp_sratio_;
150 stream & dpProps_->dp_mode_;
151 stream & dpProps_->dp_F_;
152 }
153 stream & b;
154 if (b) {
155 if (!energies_)
156 energies_ = NEWC(Energies());
157 stream & energies_->estrain_;
158 stream & energies_->eslip_;
159 stream & energies_->edashpot_;
160 }
161 }
162
163 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 5) {
164 stream & sn_tabPos_;
165 stream & sn_cohres_;
166 }
167
168 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 6)
169 stream & sn_por_;
170
171 } else {
172 // Backward compatibility
173 stream & kTran_;
174 stream & E_;
175 stream & kRot_;
176 stream & fictForce_;
177 stream & poisson_;
178 stream & fric_;
179 stream & sn_mode_;
180 stream & sn_F_;
181 stream & sn_M_;
182 stream & sn_S_;
183 stream & sn_BS_;
184 stream & sn_TS_;
185 stream & sn_rmul_;
186
187 bool b(false);
188 stream & b;
189 if (b) {
190 if (!dpProps_)
191 dpProps_ = NEWC(dpProps());
192 stream & dpProps_->dp_nratio_;
193 stream & dpProps_->dp_sratio_;
194 stream & dpProps_->dp_mode_;
195 stream & dpProps_->dp_F_;
196 }
197 stream & b;
198 if (b) {
199 if (!energies_)
200 energies_ = NEWC(Energies());
201 stream & energies_->estrain_;
202 stream & energies_->eslip_;
203 stream & energies_->edashpot_;
204 }
205 stream & b;
206 if (b) {
207 int vi = 0;
208 stream & vi;
209 sn_state_ = abs(vi);
210 double val = 0.0;
211 stream & val;
212 tenTable_[0].rx() = val;
213 val = 0.0;
214 stream & val;
215 cohTable_[0].rx() = val;
216 stream & sn_fa_;
217 stream & sn_mcf_;
218 stream & val;
219 stream & val;
220 stream & val;
221 stream & val;
222 Quat q;
223 stream & q;
224 stream & val;
225 stream & val;
226 }
227
228 stream & inheritanceField_;
229 stream & effectiveTranslationalStiffness_;
230 stream & effectiveRotationalStiffness_;
231
232 stream & sn_bmul_;
233 stream & sn_tmul_;
234 stream & userArea_;
235 stream & rgap_;
236
237 if (stream.getRestoreVersion() > 1)
238 stream & kRatio_;
239
240 if (stream.getRestoreVersion() > 2) {
241 uint32 v;
242 stream & v;
243 poisOffDiag_ = v == 0 ? false : true;
244 }
245
246 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 3)
247 stream & forceSet_;
248
249 }
250 }
251
252 void ContactModelRBSN::copy(const ContactModel *cm) {
253 // Copy all of the contact model properties. Used in the CMAT
254 // when a new contact is created.
255 ContactModelMechanical::copy(cm);
256 const ContactModelRBSN *in = dynamic_cast<const ContactModelRBSN*>(cm);
257 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
258 fictForce_ = in->fictForce_;
259 sn_F_ = in->sn_F_;
260 sn_sdisp_ = in->sn_sdisp_;
261 sn_M_ = in->sn_M_;
262 kRot_ = in->kRot_;
263 kTran_ = in->kTran_;
264 kRatio_ = in->kRatio_;
265 E_ = in->E_;
266 poisson_ = in->poisson_;
267 fric_ = in->fric_;
268 sn_bmul_ = in->sn_bmul_;
269 sn_tmul_ = in->sn_tmul_;
270 sn_rmul_ = in->sn_rmul_;
271 userArea_ = in->userArea_;
272 rgap_ = in->rgap_;
273 sn_fa_ = in->sn_fa_;
274 sn_mcf_ = in->sn_mcf_;
275 sn_dil_ = in->sn_dil_;
276 sn_dilzero_ = in->sn_dilzero_;
277 transTen_ = in->transTen_;
278 sn_elong_ = in->sn_elong_;
279 sn_ndisp_ = in->sn_ndisp_;
280 sn_mode_ = in->sn_mode_;
281 sn_state_ = in->sn_state_;
282 poisOffDiag_ = in->poisOffDiag_;
283 sn_S_ = in->sn_S_;
284 sn_BS_ = in->sn_BS_;
285 sn_TS_ = in->sn_TS_;
286 forceSet_ = in->forceSet_;
287 sn_heal_ = in->sn_heal_;
288 tenTable_ = in->tenTable_;
289 cohTable_ = in->cohTable_;
290 sn_por_ = in->sn_por_;
291 if (in->hasDamping()) {
292 if (!dpProps_)
293 dpProps_ = NEWC(dpProps());
294 dp_nratio(in->dp_nratio());
295 dp_sratio(in->dp_sratio());
296 dp_mode(in->dp_mode());
297 dp_F(in->dp_F());
298 }
299 if (in->hasEnergies()) {
300 if (!energies_)
301 energies_ = NEWC(Energies());
302 estrain(in->estrain());
303 eslip(in->eslip());
304 edashpot(in->edashpot());
305 }
306 inheritanceField(in->inheritanceField());
307 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
308 effectiveRotationalStiffness(in->effectiveRotationalStiffness());
309 }
310
311
312 QVariant ContactModelRBSN::getProperty(uint32 i,const IContact *con) const {
313 // Return the property. The IContact pointer is provided so that
314 // more complicated properties, depending on contact characteristics,
315 // can be calcualted.
316 // The IContact pointer may be a nullptr!
317 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
318 QVariant var;
319 switch (i) {
320 case kwKn: return kTran_;
321 case kwKs: return kTran_ / kRatio_;
322 case kwKrot: var.setValue(kRot_); return var;
323 case kwFric: return fric_;
324 case kwBMul: return sn_bmul_;
325 case kwTMul: return sn_tmul_;
326 case kwSNMode: return sn_mode_;
327 case kwSNF: var.setValue(sn_F_); return var;
328 case kwSNM: var.setValue(sn_M_); return var;
329 case kwSNS: return sn_S_;
330 case kwSNBS: return sn_BS_;
331 case kwSNTS: return sn_TS_;
332 case kwPoisDiag: return poisOffDiag_ == false ? 0 : 1;
333 case kwSNRMul: return sn_rmul_;
334 case kwSNRadius: {
335 if (!c) return 0.0;
336 double Cmax1 = c->getEnd1Curvature().y();
337 double Cmax2 = c->getEnd2Curvature().y();
338 if (!userArea_)
339 return sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
340 else {
341#ifdef THREED
342 double rad = std::sqrt(userArea_ / dPi);
343#else
344 double rad = userArea_ / 2.0;
345#endif
346 return rad;
347 }
348
349 }
350 case kwEmod: return E_;
351 case kwKRatio: return kRatio_;
352 case kwDpNRatio: return dpProps_ ? dpProps_->dp_nratio_ : 0;
353 case kwDpSRatio: return dpProps_ ? dpProps_->dp_sratio_ : 0;
354 case kwDpMode: return dpProps_ ? dpProps_->dp_mode_ : 0;
355 case kwDpF: {
356 dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
357 return var;
358 }
359 case kwSNState: return sn_state_;
360 case kwSNTStr: return sn_Ten();
361 case kwSNSStr: {
362 if (!c) return 0.0;
363 double area = computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x();
364 return shearStrength(area);
365 }
366 case kwSNCoh: return cohTable_[0].x();
367 case kwSNFa: return std::atan(sn_fa_)/dDegrad;
368 case kwSNMCF: return sn_mcf_;
369 case kwSNSig: {
370 if (!c) return 0.0;
371 return SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x();
372 }
373 case kwSNTau: {
374 if (!c) return 0.0;
375 return SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).y();
376 }
377 case kwSNArea: {
378 if (userArea_) return userArea_;
379 if (!c)
380 return 0.0;
381 return computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x();
382 }
383 case kwUserArea:
384 return userArea_;
385 case kwRGap:
386 return rgap_;
387 case kwPForce: var.setValue(fictForce_); return var;
388 case kwPois: return poisson_;
389 case kwSnCohRes : return sn_cohres_;
390 //case kwSnTenRes : return sn_tenres();
391 case kwSnDil : return std::atan(sn_dil_)/dDegrad;
392 case kwSnDilZ : return sn_dilzero_;
393 case kwSnNormDisp: return sn_ndisp_;
394 case kwSnShearDisp: var.setValue(sn_sdisp_); return var;
395 case kwSnCohDc : return sn_cohdc();
396 case kwSnTenDc : return sn_tendc();
397 case kwSnHeal : return sn_heal_;
398 case kwTenTable:
399 if (sn_tabPos_ < tenTable_.size())
400 if (sn_tabPos_ == 0)
401 var.setValue(DVect2(1,0));
402 else
403 var.setValue(tenTable_[sn_tabPos_]);
404 else
405 var.setValue(DVect2(0,0));
406 return var;
407 case kwCohTable:
408 if (sn_tabPos_ < cohTable_.size())
409 if (sn_tabPos_ == 0)
410 var.setValue(DVect2(1,0));
411 else
412 var.setValue(cohTable_[sn_tabPos_]);
413 else
414 var.setValue(DVect2(0,0));
415 return var;
416 case kwTablePos: return sn_tabPos_+1;
417 case kwPorP: return sn_por_;
418 case kwStressNorm: {
419 if (!c)
420 return 0.0;
421 return SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x() + sn_por_;
422 }
423 }
424 assert(0);
425 return QVariant();
426 }
427
428 bool ContactModelRBSN::getPropertyGlobal(uint32 i) const {
429 // Returns whether or not a property is held in the global axis system (TRUE)
430 // or the local system (FALSE). Used by the plotting logic.
431 switch (i) {
432 case kwSNF:
433 case kwSNM:
434 case kwDpF:
435 return false;
436 }
437 return true;
438 }
439
440 bool ContactModelRBSN::setProperty(uint32 i,const QVariant &v,IContact *) {
441 // Set a contact model property. Return value indicates that the timestep
442 // should be recalculated.
443 dpProps dp;
444 switch (i) {
445 case kwKn: {
446 if (!v.canConvert<double>())
447 throw Exception("kn must be a double.");
448 double val(v.toDouble());
449 if (val<0.0)
450 throw Exception("Negative kn not allowed.");
451 kTran_ = val;
452 return true;
453 }
454 case kwKs: {
455 if (!v.canConvert<double>())
456 throw Exception("ks must be a double.");
457 double val(v.toDouble());
458 if (val<0.0)
459 throw Exception("Negative ks not allowed.");
460 if (!kTran_)
461 kTran_ = val;
462 else
463 kRatio_ = kTran_ / val;
464 return true;
465 }
466 case kwKrot: {
467 DAVect val(0.0);
468#ifdef TWOD
469 if (!v.canConvert<DAVect>() && !v.canConvert<double>())
470 throw Exception("krot must be an angular vector.");
471 if (v.canConvert<DAVect>())
472 val = DAVect(v.value<DAVect>());
473 else
474 val = DAVect(v.toDouble());
475#else
476 if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
477 throw Exception("krot must be an angular vector.");
478 if (v.canConvert<DAVect>())
479 val = DAVect(v.value<DAVect>());
480 else
481 val = DAVect(v.value<DVect>());
482#endif
483 kRot_ = val;
484 return false;
485 }
486 case kwFric: {
487 if (!v.canConvert<double>())
488 throw Exception("fric must be a double.");
489 double val(v.toDouble());
490 if (val<0.0)
491 throw Exception("Negative fric not allowed.");
492 fric_ = val;
493 //if (!sn_fa_)
494 // sn_fa_ = fric_;
495 return false;
496 }
497 case kwBMul: {
498 if (!v.canConvert<double>())
499 throw Exception("sn_bmul must be a double.");
500 double val(v.toDouble());
501 if (val<0.0)
502 throw Exception("Negative sn_bmul not allowed.");
503 sn_bmul_ = val;
504 return false;
505 }
506 case kwTMul: {
507 if (!v.canConvert<double>())
508 throw Exception("sn_tmul must be a double.");
509 double val(v.toDouble());
510 if (val<0.0)
511 throw Exception("Negative st_bmul not allowed.");
512 sn_tmul_ = val;
513 return false;
514 }
515 case kwSNMode: {
516 if (!v.canConvert<uint32>())
517 throw Exception("sn_mode must be 0 (absolute) or 1 (incremental).");
518 double val(v.toUInt());
519 if (val>1)
520 throw Exception("sn_mode must be 0 (absolute) or 1 (incremental).");
521 sn_mode_ = val;
522 return false;
523 }
524 case kwSNRMul: {
525 if (!v.canConvert<double>())
526 throw Exception("rmul must be a double.");
527 double val(v.toDouble());
528 if (val<0.0)
529 throw Exception("Negative rmul not allowed.");
530 sn_rmul_ = val;
531 return false;
532 }
533 case kwSNF: {
534 if (!v.canConvert<DVect>())
535 throw Exception("sn_force must be a vector.");
536 DVect val(v.value<DVect>());
537 sn_F_ = val;
538 return false;
539 }
540 case kwSNM: {
541 DAVect val(0.0);
542#ifdef TWOD
543 if (!v.canConvert<DAVect>() && !v.canConvert<double>())
544 throw Exception("res_moment must be an angular vector.");
545 if (v.canConvert<DAVect>())
546 val = DAVect(v.value<DAVect>());
547 else
548 val = DAVect(v.toDouble());
549#else
550 if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
551 throw Exception("res_moment must be an angular vector.");
552 if (v.canConvert<DAVect>())
553 val = DAVect(v.value<DAVect>());
554 else
555 val = DAVect(v.value<DVect>());
556#endif
557 sn_M_ = val;
558 return false;
559 }
560 case kwDpNRatio: {
561 if (!v.canConvert<double>())
562 throw Exception("dp_nratio must be a double.");
563 double val(v.toDouble());
564 if (val<0.0)
565 throw Exception("Negative dp_nratio not allowed.");
566 if (val == 0.0 && !dpProps_)
567 return false;
568 if (!dpProps_)
569 dpProps_ = NEWC(dpProps());
570 dpProps_->dp_nratio_ = val;
571 return true;
572 }
573 case kwDpSRatio: {
574 if (!v.canConvert<double>())
575 throw Exception("dp_sratio must be a double.");
576 double val(v.toDouble());
577 if (val<0.0)
578 throw Exception("Negative dp_sratio not allowed.");
579 if (val == 0.0 && !dpProps_)
580 return false;
581 if (!dpProps_)
582 dpProps_ = NEWC(dpProps());
583 dpProps_->dp_sratio_ = val;
584 return true;
585 }
586 case kwDpMode: {
587 if (!v.canConvert<int>())
588 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
589 int val(v.toInt());
590 if (val == 0 && !dpProps_)
591 return false;
592 if (val < 0 || val > 3)
593 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
594 if (!dpProps_)
595 dpProps_ = NEWC(dpProps());
596 dpProps_->dp_mode_ = val;
597 return false;
598 }
599 case kwDpF: {
600 if (!v.canConvert<DVect>())
601 throw Exception("dp_force must be a vector.");
602 DVect val(v.value<DVect>());
603 if (val.fsum() == 0.0 && !dpProps_)
604 return false;
605 if (!dpProps_)
606 dpProps_ = NEWC(dpProps());
607 dpProps_->dp_F_ = val;
608 return false;
609 }
610 case kwSNTStr: {
611 if (!v.canConvert<double>())
612 throw Exception("sn_ten must be a double.");
613 double val(v.toDouble());
614 if (val < 0.0)
615 throw Exception("Negative sn_ten not allowed.");
616 tenTable_[0].rx() = val;
617 return false;
618 }
619 case kwSNCoh: {
620 if (!v.canConvert<double>())
621 throw Exception("sn_coh must be a double.");
622 double val(v.toDouble());
623 if (val<0.0)
624 throw Exception("Negative pb_coh not allowed.");
625 cohTable_[0].rx() = val;
626 return false;
627 }
628 case kwSNFa: {
629 if (!v.canConvert<double>())
630 throw Exception("sn_fa must be a double.");
631 double val(v.toDouble());
632 if (val<0.0)
633 throw Exception("Negative sn_fa not allowed.");
634 if (val >= 90.0)
635 throw Exception("sn_fa must be lower than 90.0 degrees.");
636 sn_fa_ = std::tan(val*dDegrad);
637 if (!fric_)
638 fric_ = sn_fa_;
639 return false;
640 }
641 case kwSNMCF: {
642 if (!v.canConvert<double>())
643 throw Exception("sn_mcf must be a double.");
644 double val(v.toDouble());
645 if (val<0.0)
646 throw Exception("Negative sn_mcf not allowed.");
647 if (val > 1.0)
648 throw Exception("sn_mcf must be lower or equal to 1.0.");
649 sn_mcf_ = val;
650 return false;
651 }
652 case kwSNArea:
653 case kwUserArea: {
654 if (!v.canConvert<double>())
655 throw Exception("area must be a double.");
656 double val(v.toDouble());
657 if (val < 0.0)
658 throw Exception("Negative area not allowed.");
659 if (userArea_ && val) {
660 double rat = userArea_ / val;
661 kTran_ *= rat;
662 kRot_ *= rat;
663 }
664 userArea_ = val;
665 return true;
666 }
667 case kwRGap: {
668 if (!v.canConvert<double>())
669 throw Exception("Reference gap must be a double.");
670 double val(v.toDouble());
671 rgap_ = val;
672 return false;
673 }
674 case kwPois: {
675 if (!v.canConvert<double>())
676 throw Exception("Reference poisson must be a double.");
677 double val(v.toDouble());
678 poisson_ = val;
679 return false;
680 }
681 case kwPoisDiag: {
682 if (!v.canConvert<uint32>())
683 throw Exception("Reference diagonal must be an integer.");
684 uint32 b(v.toUInt());
685 if (b > 1)
686 throw Exception("diagonal must be 0 (diagonal terms only) or 1 (all terms).");
687 poisOffDiag_ = b == 0 ? false : true;
688 return false;
689 }
690 case kwSnCohRes: {
691 bool ok;
692 double val(v.toDouble(&ok));
693 if (!ok || val<0.0)
694 throw Exception("sn_cohres must be a positive double.");
695 sn_cohres_ = val;
696 return false;
697 }
698 case kwSnDil: {
699 bool ok;
700 double val(v.toDouble(&ok));
701 if (!ok || val<0.0)
702 throw Exception("sn_dil must be a positive double.");
703 sn_dil_ = std::tan(val*dDegrad);
704 return false;
705 }
706 case kwSnDilZ: {
707 bool ok;
708 double val(v.toDouble(&ok));
709 if (!ok || val<0.0)
710 throw Exception("sn_dil_zero must be a positive double.");
711 sn_dilzero_ = val;
712 return false;
713 }
714 case kwSnNormDisp: {
715 bool ok;
716 double val(v.toDouble(&ok));
717 if (!ok)
718 throw Exception("sn_ndisp must be a positive double.");
719 sn_ndisp_ = val;
720 return false;
721 }
722 case kwSnShearDisp: {
723 if (!v.canConvert<DVect2>())
724 throw Exception("sn_sdisp must be a vector.");
725 DVect2 val(v.value<DVect2>());
726 sn_sdisp_ = val;
727 return false;
728 }
729 case kwSnCohDc: {
730 bool ok;
731 double val(v.toDouble(&ok));
732 if (!ok || val<0.0)
733 throw Exception("sn_cohdc must be a positive double.");
734 if (cohTable_.size() == 1)
735 cohTable_.push_back(DVect2(0,val));
736 else {
737 cohTable_.back().ry() = val;
738 std::sort(cohTable_.begin(),cohTable_.end(), [](const DVect2& a, const DVect2& b) { return a.y() < b.y(); } );
739 while (cohTable_.back().y() > val)
740 cohTable_.pop_back();
741 }
742 if (sn_state_ == 0)
743 sn_state_ = 6;
744 return false;
745 }
746 case kwSnTenDc: {
747 bool ok;
748 double val(v.toDouble(&ok));
749 if (!ok || val<0.0)
750 throw Exception("sn_tendc must be a positive double.");
751 if (tenTable_.size() == 1)
752 tenTable_.push_back(DVect2(0,val));
753 else {
754 tenTable_.back().ry() = val;
755 std::sort(tenTable_.begin(),tenTable_.end(), [](const DVect2& a, const DVect2& b) { return a.y() < b.y(); } );
756 while (tenTable_.back().y() > val)
757 tenTable_.pop_back();
758 }
759 return false;
760 }
761 case kwSnHeal: {
762 bool ok;
763 int val(v.toInt(&ok));
764 if (!ok || (val != 0 && val != 1))
765 throw Exception("sn_heal must be 0 or 1.");
766 sn_heal_ = val == 0 ? false : true;
767 return false;
768 }
769 case kwTenTable: {
770 if (!v.canConvert<DVect2>())
771 throw Exception("sn_tentab entry must be a strength and distance.");
772 DVect2 vl(v.value<DVect2>());
773 if (vl.x() < 0 || vl.y() < 0)
774 throw Exception("The sn_tentab entries must be positive.");
775 if (vl.y() == 0)
776 throw Exception("Use sn_ten to set the tensile strength at 0 elongation.");
777 tenTable_.push_back(vl);
778 std::sort(tenTable_.begin(),tenTable_.end(), [](const DVect2& a, const DVect2& b) { return a.y() < b.y(); } );
779 }
780 return false;
781 case kwCohTable: {
782 if (!v.canConvert<DVect2>())
783 throw Exception("sn_cohtab entry must be a strength and distance.");
784 DVect2 vl(v.value<DVect2>());
785 if (vl.x() < 0 || vl.y() < 0)
786 throw Exception("The sn_cohtab entries must be positive.");
787 if (vl.y() == 0)
788 throw Exception("Use sn_coh to set the cohesive strength at 0 elongation.");
789 cohTable_.push_back(vl);
790 std::sort(cohTable_.begin(),cohTable_.end(), [](const DVect2& a, const DVect2& b) { return a.y() < b.y(); } );
791 }
792 return false;
793 case kwTablePos: {
794 bool ok;
795 int val(v.toInt(&ok));
796 if (!ok || val < 1)
797 throw Exception("sn_tabpos must be 1 or greater.");
798 sn_tabPos_ = val - 1;
799 return false;
800 }
801 case kwPorP: {
802 if (!v.canConvert<double>())
803 throw Exception("sn_porp must be a double.");
804 double val = v.toDouble();
805 sn_por_ = val;
806 return true;
807 }
808 }
809 return false;
810 }
811
812 bool ContactModelRBSN::getPropertyReadOnly(uint32 i) const {
813 // Returns TRUE if a property is read only or FALSE otherwise.
814 switch (i) {
815 case kwDpF:
816 case kwSNS:
817 case kwSNBS:
818 case kwSNTS:
819 case kwEmod:
820 case kwKRatio:
821 case kwSNState:
822 case kwSNRadius:
823 case kwSNSStr:
824 case kwSNSig:
825 case kwSNTau:
826 case kwPForce:
827 case kwStressNorm:
828 return true;
829 default:
830 break;
831 }
832 return false;
833 }
834
835 bool ContactModelRBSN::supportsInheritance(uint32 i) const {
836 // Returns TRUE if a property supports inheritance or FALSE otherwise.
837 switch (i) {
838 case kwKn:
839 case kwKs:
840 case kwFric:
841 return true;
842 default:
843 break;
844 }
845 return false;
846 }
847
848 QString ContactModelRBSN::getMethodArguments(uint32 i) const {
849 // Return a list of contact model method argument names.
850 switch (i) {
851 case kwAssignStiffness:
852 return "kn,kratio";
853 case kwStiffness:
854 return "emod,poisson";
855 case kwBond:
856 return "gap";
857 case kwUnbond:
858 return "gap";
859 case kwArea:
860 case kwResetDisp:
861 return QString();
862 }
863 assert(0);
864 return QString();
865 }
866
867 bool ContactModelRBSN::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
868 // Apply the specified method.
869 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
870 switch (i) {
871 case kwAssignStiffness: {
872 poisson_ = 0.0;
873 if (vl.at(0).isNull())
874 throw Exception("Argument kn must be specified with method assign-stiffness in contact model %1.",getName());
875 double knpa = vl.at(0).toDouble();
876 if (knpa<=0.0)
877 throw Exception("Negative or zero kn not allowed in contact model %1.",getName());
878 if (vl.at(1).isNull())
879 throw Exception("Argument kratio must be specified with method assign-stiffness in contact model %1.",getName());
880 kRatio_ = vl.at(1).toDouble();
881 if (kRatio_<0.0) {
882 kRatio_ = 0.0;
883 throw Exception("Negative kratio not allowed in contact model %1.",getName());
884 }
885 std::vector<DVect> pts;
886 c->getJointGeometry(&pts);
887 double area = 0.0;
888#ifdef THREED
889 // Step 1: get centroid and area
890 for (int i=1; i<pts.size()-1; ++i) {
891 double a = (pts[0] - pts[i]).mag();
892 double b = (pts[0] - pts[i+1]).mag();
893 double c = (pts[i] - pts[i+1]).mag();
894 double la = 0.0;
895 if (b > a)
896 std::swap(a,b);
897 if (c > a)
898 std::swap(a,c);
899 if (c > b)
900 std::swap(b,c);
901 if (c - (a - b) >= 0)
902 la = 0.25 * sqrt((a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c)));
903 area += la;
904 }
905#else
906 // Assume unit thickness in the out of plane direction
907 area = (pts[1] - pts[0]).mag();
908#endif
909 userArea_ = area;
910 kTran_ = knpa * area;
911 E_ = kTran_ / area;
912 kRot_ = DAVect(0.0);
913 setInheritance(1,false);
914 setInheritance(2,false);
915 sn_mode_ = 1.0;
916 return true;
917 }
918 case kwStiffness: {
919 FP_S;
920 poisson_ = 0.0;
921 if (vl.at(0).isNull())
922 throw Exception("Argument emod must be specified with method compute-stiffness in contact model %1.",getName());
923 E_ = vl.at(0).toDouble();
924 if (E_<=0.0)
925 throw Exception("Negative or zero emod not allowed in contact model %1.",getName());
926 if (vl.at(1).isNull())
927 throw Exception("Argument poisson must be specified with method compute-stiffness in contact model %1.",getName());
928 poisson_ = vl.at(1).toDouble();
929 if (poisson_ < 0.0) {
930 poisson_ = 0.0;
931 throw Exception("Negative poisson not allowed in contact model %1.",getName());
932 }
933 const IBody *b1 = con->getEnd1()->getIBody();
934 const IBody *b2 = con->getEnd2()->getIBody();
935 //double vol1 = b1->getVolume();
936 //double vol2 = b2->getVolume();
937 //if (std::max(vol1,vol2) > 10.*std::min(vol1,vol2))
938 // poisson_ = 0.0;
939 DVect pos1 = toDVect(b1->getIThing()->getLocation());
940 DVect pos2 = toDVect(b2->getIThing()->getLocation()) + con->getOffSet();
941 double dist = (pos1-pos2).mag();
942 if (con->withWall())
943 dist = (pos1 - con->getPosition()).mag();
944 double tol = std::max(pos1.abs().maxComp(),pos2.abs().maxComp())*limits<double>::epsilon()*1000;
945 if (dist < tol) {
946 poisson_ = 0;
947 userArea_ = 0;
948 kTran_ = 0;
949 kRot_.fill(0);
950 return true;
951 }
952 std::vector<DVect> pts;
953 FP_S;
954 c->getJointGeometry(&pts);
955 FP_S;
956 double area = 0.0;
957 DAVect inertia(0.0);
958#ifdef THREED
959 // Step 1: get centroid and area
960 DVect cm(0.0);
961 for (int i=1; i<pts.size()-1; ++i) {
962 DVect lcm = (pts[0] + pts[i] + pts[i+1])/3.0;
963 double a = (pts[0] - pts[i]).mag();
964 double b = (pts[0] - pts[i+1]).mag();
965 double c = (pts[i] - pts[i+1]).mag();
966 double la = 0.0;
967 if (b > a)
968 std::swap(a,b);
969 if (c > a)
970 std::swap(a,c);
971 if (c > b)
972 std::swap(b,c);
973 if (c - (a - b) >= 0)
974 la = 0.25 * sqrt((a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c)));
975 cm += lcm * la;
976 area += la;
977 }
978 FP_S;
979 if (area == 0.0) {
980 poisson_ = 0;
981 userArea_ = 0;
982 kTran_ = 0;
983 kRot_.fill(0);
984 return true;
985 }
986 cm /= area;
987 FP_S;
988 // Step 2 - center it and put in the local system
989 for (int i=0; i<pts.size(); ++i) {
990 pts[i] -= cm;
991 pts[i] = con->getLocalSystem().toLocal(pts[i]);
992 }
993 // Step 3: compute the polar inertia
994 for (int i=0; i<pts.size(); ++i) {
995 int j = i < pts.size()-1 ? i+1 : 0;
996 double xi = pts[i].y();
997 double xip1 = pts[j].y();
998 double yi = pts[i].z();
999 double yip1 = pts[j].z();
1000 double frnt = (xi*yip1-xip1*yi);
1001 inertia.ry() += frnt*(xi*xi+xi*xip1+xip1*xip1);
1002 inertia.rz() += frnt*(yi*yi+yi*yip1+yip1*yip1);
1003 }
1004 inertia.ry() = std::abs(inertia.y() / 12.);
1005 inertia.rz() = std::abs(inertia.z() / 12.);
1006 inertia.rx() = inertia.y() + inertia.z();
1007
1008#else
1009 // Assume unit thickness in the out of plane direction
1010 area = (pts[1] - pts[0]).mag();
1011 inertia.rz() = area*area*area/12.;
1012#endif
1013 userArea_ = area;
1014 kTran_ = E_ * area / dist;
1015 kRot_ = inertia *E_ / dist;
1016 setInheritance(1,false);
1017 setInheritance(2,false);
1018 sn_mode_ = 1.0;
1019 return true;
1020 }
1021 case kwBond: {
1022 if (sn_state_ == 3) return false;
1023 double mingap = -1.0 * limits<double>::max();
1024 double maxgap = 0;
1025 if (vl.at(0).canConvert<double>())
1026 maxgap = vl.at(0).toDouble();
1027 else if (vl.at(0).canConvert<DVect2>()) {
1028 DVect2 value = vl.at(0).value<DVect2>();
1029 mingap = value.minComp();
1030 maxgap = value.maxComp();
1031 }
1032 else if (!vl.at(0).isNull())
1033 throw Exception("gap value %1 not recognized in method bond in contact model %2.", vl.at(1), getName());
1034 double gap = c->getGap();
1035 if (gap >= mingap && gap <= maxgap) {
1036 sn_state_ = 3;
1037 sn_mode_ = 1;
1038 return true;
1039 }
1040 return false;
1041 }
1042 case kwUnbond: {
1043 if (sn_state_ == 0) return false;
1044 double mingap = -1.0 * limits<double>::max();
1045 double maxgap = 0;
1046 if (vl.at(0).canConvert<double>())
1047 maxgap = vl.at(0).toDouble();
1048 else if (vl.at(0).canConvert<DVect2>()) {
1049 DVect2 value = vl.at(0).value<DVect2>();
1050 mingap = value.minComp();
1051 maxgap = value.maxComp();
1052 }
1053 else if (!vl.at(0).isNull())
1054 throw Exception("gap value %1 not recognized in method unbond in contact model %2.", vl.at(0), getName());
1055 double gap = c->getGap();
1056 if (gap >= mingap && gap <= maxgap) {
1057 sn_state_ = 0;
1058 return true;
1059 }
1060 return false;
1061 }
1062 case kwArea: {
1063 if (!userArea_) {
1064 double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
1065#ifdef THREED
1066 userArea_ = rsq * rsq * dPi;
1067#else
1068 userArea_ = rsq * 2.0;
1069#endif
1070 }
1071 return true;
1072 }
1073 case kwResetDisp:
1074 sn_ndisp_ = 0.0;
1075 for (int i=1; i<dim; ++i)
1076 sn_sdisp_.rdof(i) = 0;
1077 break;
1078 }
1079 return false;
1080 }
1081
1082 double ContactModelRBSN::getEnergy(uint32 i) const {
1083 // Return an energy value.
1084 double ret(0.0);
1085 if (!energies_)
1086 return ret;
1087 switch (i) {
1088 case kwEStrain: return energies_->estrain_;
1089 case kwESlip: return energies_->eslip_;
1090 case kwEDashpot: return energies_->edashpot_;
1091 }
1092 assert(0);
1093 return ret;
1094 }
1095
1096 bool ContactModelRBSN::getEnergyAccumulate(uint32 i) const {
1097 // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
1098 switch (i) {
1099 case kwEStrain: return false;
1100 case kwESlip: return true;
1101 case kwEDashpot: return true;
1102 }
1103 assert(0);
1104 return false;
1105 }
1106
1107 void ContactModelRBSN::setEnergy(uint32 i,const double &d) {
1108 // Set an energy value.
1109 if (!energies_) return;
1110 switch (i) {
1111 case kwEStrain: energies_->estrain_ = d; return;
1112 case kwESlip: energies_->eslip_ = d; return;
1113 case kwEDashpot: energies_->edashpot_= d; return;
1114 }
1115 assert(0);
1116 return;
1117 }
1118
1119 bool ContactModelRBSN::validate(ContactModelMechanicalState *state,const double &) {
1120 // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
1121 // (1) the contact is created, (2) a property of the contact that returns a true via
1122 // the setProperty method has been modified and (3) when a set of cycles is executed
1123 // via the {cycle N} command.
1124 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
1125 assert(state);
1126 const IContactMechanical *c = state->getMechanicalContact();
1127 assert(c);
1128
1129 if (state->trackEnergy_)
1130 activateEnergy();
1131
1132 if (inheritanceField_ & KnMask)
1133 updateKn(c);
1134 if (inheritanceField_ & KsMask)
1135 updateKs(c);
1136 if (inheritanceField_ & FricMask)
1137 updateFric(c);
1138
1139 updateStiffness(state);
1140 return checkActivity(state->gap_);
1141 }
1142
1143 static const QString knstr("kn");
1144 bool ContactModelRBSN::updateKn(const IContactMechanical *con) {
1145 assert(con);
1146 QVariant v1 = con->getEnd1()->getProperty(knstr);
1147 QVariant v2 = con->getEnd2()->getProperty(knstr);
1148 if (!v1.isValid() || !v2.isValid())
1149 return false;
1150 double kn1 = v1.toDouble();
1151 double kn2 = v2.toDouble();
1152 double val = kTran_;
1153 if (kn1 && kn2)
1154 kTran_ = kn1*kn2/(kn1+kn2);
1155 else if (kn1)
1156 kTran_ = kn1;
1157 else if (kn2)
1158 kTran_ = kn2;
1159 return ( (kTran_ != val) );
1160 }
1161
1162 static const QString ksstr("ks");
1163 bool ContactModelRBSN::updateKs(const IContactMechanical *con) {
1164 assert(con);
1165 QVariant v1 = con->getEnd1()->getProperty(ksstr);
1166 QVariant v2 = con->getEnd2()->getProperty(ksstr);
1167 if (!v1.isValid() || !v2.isValid())
1168 return false;
1169 double ks1 = v1.toDouble();
1170 double ks2 = v2.toDouble();
1171 double val = kTran_;
1172 if (ks1 && ks2)
1173 kTran_ = ks1*ks2/(ks1+ks2);
1174 else if (ks1)
1175 kTran_ = ks1;
1176 else if (ks2)
1177 kTran_ = ks2;
1178 return ( (kTran_ != val) );
1179 }
1180
1181 static const QString fricstr("fric");
1182 bool ContactModelRBSN::updateFric(const IContactMechanical *con) {
1183 assert(con);
1184 QVariant v1 = con->getEnd1()->getProperty(fricstr);
1185 QVariant v2 = con->getEnd2()->getProperty(fricstr);
1186 if (!v1.isValid() || !v2.isValid())
1187 return false;
1188 double fric1 = std::max(0.0,v1.toDouble());
1189 double fric2 = std::max(0.0,v2.toDouble());
1190 double val = fric_;
1191 fric_ = std::min(fric1,fric2);
1192 return ( (fric_ != val) );
1193 }
1194
1195 bool ContactModelRBSN::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
1196 // The endPropertyUpdated method is called whenever a surface property (with a name
1197 // that matches an inheritable contact model property name) of one of the contacting
1198 // pieces is modified. This allows the contact model to update its associated
1199 // properties. The return value denotes whether or not the update has affected
1200 // the time step computation (by having modified the translational or rotational
1201 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
1202 assert(c);
1203 QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
1204 QRegExp rx(name,Qt::CaseInsensitive);
1205 int idx = availableProperties.indexOf(rx)+1;
1206 bool ret=false;
1207
1208 if (idx<=0)
1209 return ret;
1210
1211 switch(idx) {
1212 case kwKn: { //kn
1213 if (inheritanceField_ & KnMask)
1214 ret = updateKn(c);
1215 break;
1216 }
1217 case kwKs: { //ks
1218 if (inheritanceField_ & KsMask)
1219 ret =updateKs(c);
1220 break;
1221 }
1222 case kwFric: { //fric
1223 if (inheritanceField_ & FricMask)
1224 updateFric(c);
1225 break;
1226 }
1227 }
1228 return ret;
1229 }
1230
1231 ContactModelRBSN::StiffData ContactModelRBSN::computeStiffData(ContactModelMechanicalState *state) const {
1232 // Update contact data
1233 //double Cmin1 = state->end1Curvature_.x();
1234 double Cmax1 = state->end1Curvature_.y();
1235 double Cmax2 = state->end2Curvature_.y();
1236 double br = sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
1237 if (userArea_)
1238#ifdef THREED
1239 br = std::sqrt(userArea_ / dPi);
1240#else
1241 br = userArea_ / 2.0;
1242#endif
1243 StiffData ret;
1244 ret.reff_ = br;
1245 ret.trans_ = DVect2(kTran_,kTran_/kRatio_);
1246 ret.ang_ = kRot_;
1247 return ret;
1248 }
1249
1250 void ContactModelRBSN::updateStiffness(ContactModelMechanicalState *state) {
1251 // first compute stiffness data
1252 StiffData stiff = computeStiffData(state);
1253 // Now calculate effective stiffness
1254 DVect2 retT = stiff.trans_;
1255 // correction if viscous damping active
1256 if (dpProps_) {
1257 DVect2 correct(1.0);
1258 if (dpProps_->dp_nratio_)
1259 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
1260 if (dpProps_->dp_sratio_)
1261 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
1262 retT /= (correct*correct);
1263 }
1264 effectiveTranslationalStiffness_ = retT;
1265 // Effective rotational stiffness (bending and twisting)
1266 effectiveRotationalStiffness_ = stiff.ang_;
1267 }
1268
1269 bool ContactModelRBSN::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
1270 assert(state);
1271
1272 if (state->activated()) {
1273 // The contact was just activated from an inactive state
1274 // Trigger the FISH callback if one is hooked up to the
1275 // contact_activated event.
1276 if (cmEvents_[fActivated] >= 0) {
1277 auto c = state->getContact();
1278 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
1279 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1280 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
1281 }
1282 }
1283 updateStiffness(state);
1284 // accumulate shear displacement for dilation
1285 sn_ndisp_ += state->relativeTranslationalIncrement_.x();
1286 DVect shearInc = state->relativeTranslationalIncrement_;
1287 shearInc.rx() = 0;
1288 sn_sdisp_.ry() += shearInc.mag();
1289
1290 if (isBonded()) return FDLawBonded(state, timestep);
1291 else return FDLawUnBonded(state, timestep);
1292
1293 }
1294
1295 bool ContactModelRBSN::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
1296 // Account for thermal expansion in incremental mode
1297 if (sn_mode_ == 0 || ts->gapInc_ == 0.0) return false;
1298 DVect finc(0.0);
1299 finc.rx() = kTran_ * ts->gapInc_;
1300 sn_F_ -= finc;
1301 return true;
1302 }
1303
1304 bool ContactModelRBSN::FDLawBonded(ContactModelMechanicalState *state, const double ×tep) {
1305 // initialize ... get the geometry information
1306 DVect3 geom = computeGeomData(state->end1Curvature_,state->end2Curvature_);
1307 double area = geom.x();
1308 double bi = geom.y();
1309 double br = geom.z();
1310 double kn = kTran_;
1311 double ks = kn / kRatio_;
1312 double kb = kRot_.z();
1313#if DIM==3
1314 kb = std::sqrt(kb*kb + kRot_.y()*kRot_.y());
1315 double kt = kRot_.x();
1316#else
1317 double kt = 0.0;
1318#endif
1319
1320 DVect totalForce = sn_F_ + fictForce_;
1321
1322 //double nsmax0 = -(totalForce.x() / area) + sn_mcf_* sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z()) * br / bi;
1323
1324 // Relative translational/rotational displacement increments
1325 DVect trans = state->relativeTranslationalIncrement_;
1326 DAVect ang = state->relativeAngularIncrement_;
1327
1328 // Store previous force and moment
1329 DVect sn_F_old = totalForce;
1330 sn_F_old.rx() -= sn_por_ * area;
1331 DAVect sn_M_old = sn_M_;
1332
1333 DVect theStiff(ks);
1334 theStiff.rx() = kn;
1335 sn_F_ -= trans * theStiff;
1336 sn_M_ -= ang * kRot_;
1337 // Fictitious force captures Poisson's effect in RBSN models
1338 // (see common\contactmodel\rbsn\doc\manual\cmspringnetwork.rst)
1339 // Fictitious force between bonded bodies only.
1340 // No fictitious force between body and wall (zero volume contact end)
1341 if (poisson_ != 0.0 and (state->end1Volume_ and state->end2Volume_)) {
1342 auto stress11 = state->end1Stress_;
1343 auto stress22 = state->end2Stress_;
1344 double vol1 = state->end1Volume_;
1345 double vol2 = state->end2Volume_;
1346 double ms = 0.0;
1347 for (int i=0; i<stress11.size(); ++i) {
1348 stress11[i] = (stress11[i]*vol1 + stress22[i]*vol2)/(vol1 + vol2);
1349 ms = std::max(ms,abs(stress11[i]));
1350 }
1351 DMatrix<dim,dim> stresst(0.0);
1352#ifdef THREED
1353 stresst.get(0,0) = -poisson_*stress11[1] - poisson_*stress11[2];
1354 stresst.get(1,1) = -poisson_*stress11[0] - poisson_*stress11[2];
1355#else
1356 stresst.get(0,0) = -poisson_*stress11[1];
1357 stresst.get(1,1) = -poisson_*stress11[0];
1358#endif
1359#ifdef THREED
1360 stresst.get(2,2) = -poisson_*stress11[0] - poisson_*stress11[1];
1361#endif
1362 if (poisOffDiag_) {
1363#ifdef THREED
1364 double sxy = stress11[3];
1365 double szx = stress11[4];
1366 double syz = stress11[5];
1367 stresst.get(0,1) = poisson_ * sxy;
1368 stresst.get(1,0) = stresst.get(0,1);
1369 stresst.get(0,2) = poisson_ * szx;
1370 stresst.get(2,0) = stresst.get(0,2);
1371 stresst.get(1,2) = poisson_ * syz;
1372 stresst.get(2,1) = stresst.get(1,2);
1373#else
1374 double sxy = stress11[2];
1375 stresst.get(0,1) = poisson_ * sxy;
1376 stresst.get(1,0) = stresst.get(0,1);
1377#endif
1378 }
1379
1380
1381 DVect norm = state->axes_.e1();
1382 DVect traction = stresst * norm * userArea_;
1383 DVect shear(0.0);
1384 shear.ry() = 1.0;
1385 shear = state->axes_.toGlobal(shear);
1386#ifdef THREED
1387 DVect ns = state->axes_.toGlobal(DVect(0.,0.,1.));
1388 fictForce_ = DVect(norm|traction,shear|traction,ns|traction);
1389#else
1390 fictForce_ = DVect(norm|traction,shear|traction);
1391#endif
1392 if (forceSet_ && ms) {
1393 forceSet_ = false;
1394 sn_F_ -= fictForce_;
1395 }
1396 }
1397 FP_S;
1398 double porForce = sn_por_ * area;
1399 sn_F_.rx() -= porForce;
1400 totalForce = sn_F_ + fictForce_;
1401
1402 if (state->canFail_) {
1403 double dbend = sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z());
1404 double nsmax = -(totalForce.x() / area) + sn_mcf_*dbend * br / bi;
1405 bool failed = false;
1406 if (sn_state_ == 3 || sn_state_ == 5) {
1407 double compVal = sn_state_ == 3 ? tenTable_[0].x() : transTen_;
1408 if (nsmax >= compVal ) {
1409 if (tenTable_.back().y() < limits<double>::epsilon()*100)
1410 failed = true;
1411 else {
1412 if (sn_state_ == 3)
1413 sn_elong_ = 0;
1414 transTen_ = compVal;
1415 sn_state_ = 4;
1416 }
1417 }
1418 }
1419 FP_S;
1420
1421 if (sn_state_ == 4) {
1422 sn_elong_ += trans.x();
1423 sn_elong_ = std::max(0.0,sn_elong_);
1424 int ww = -1;
1425 if (sn_elong_ <= tenTable_.back().y()) {
1426 for (int i=0; i<tenTable_.size(); ++i) {
1427 if (tenTable_[i].y() >= sn_elong_) {
1428 ww = i;
1429 break;
1430 }
1431 }
1432 } else
1433 ww = tenTable_.size() - 1;
1434 if (ww > 0) {
1435 //double factor = std::min(1.0, sn_elong_ / tenTable_[ww].y());
1436 double prevVal = ww == 1 ? 1 : tenTable_[ww-1].x();
1437 double curVal = tenTable_[ww].x();
1438 double prevElong = tenTable_[ww-1].y();
1439 double curElong = tenTable_[ww].y();
1440 double slope = (curVal - prevVal)/(curElong - prevElong);
1441 FP_S;
1442 //y-y0 = m(x-x0)
1443 double nstren = slope * (sn_elong_ - prevElong) + prevVal;
1444 if (nstren <= 0)
1445 failed = true;
1446 else {
1447 nstren *= tenTable_[0].x();
1448 if (nsmax >= nstren || slope > 0) {
1449 double fac = nstren / nsmax;
1450 sn_F_.rx() *= fac;
1451#if DIM==3
1452 sn_M_.ry() *= fac;
1453#endif
1454 sn_M_.rz() *= fac;
1455 fictForce_.rx() *= fac;
1456 } else {
1457 sn_state_ = 5;
1458 transTen_ = -(sn_F_old.x() / area) + sn_mcf_* sqrt(sn_M_old.y()*sn_M_old.y() + sn_M_old.z()*sn_M_old.z()) * br / bi;
1459 }
1460 }
1461 }
1462 }
1463 if (sn_state_ == 6 && nsmax >= 0)
1464 failed = true;
1465 FP_S;
1466 if (failed) {
1467 // Failed in tension
1468 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1469 sn_state_ = 1;
1470 sn_F_.fill(0.0);
1471 sn_M_.fill(0.0);
1472 failed = true;
1473 fictForce_ = DVect(0.0);
1474 //sn_F_.rx() = -sn_tenres_ * area;
1475 if (cmEvents_[fBondBreak] >= 0) {
1476 auto c = state->getContact();
1477 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1478 fish::Parameter((qint64)sn_state_),
1479 fish::Parameter(nsmax),
1480 fish::Parameter(se)
1481 };
1482 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1483 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1484 }
1485 }
1486 FP_S;
1487
1488 if (!failed) {
1489 /* check for shear failure */
1490 double dtwist = sn_M_.x();
1491 DVect bfs(totalForce);
1492 bfs.rx() = 0.0;
1493 double dbfs = bfs.mag();
1494 double ssmax = dbfs / area + sn_mcf_*std::abs(dtwist) * 0.5* br / bi;
1495 double ss = shearStrength(area);
1496 FP_S;
1497 if (ss < 0)
1498 ss = 0;
1499 if (ss <= ssmax) {
1500 // strength when it breaks for
1501 // Failed in shear
1502 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1503 sn_state_ = 2;
1504 fictForce_ = DVect(0.0);
1505 FP_S;
1506 sn_F_ = totalForce;
1507 if (cmEvents_[fBondBreak] >= 0) {
1508 auto c = state->getContact();
1509 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1510 fish::Parameter((qint64)sn_state_),
1511 fish::Parameter(ss),
1512 fish::Parameter(se)
1513 };
1514 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1515 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1516 }
1517 double mm = 0.0;
1518 for (int i=1; i<dim; ++i)
1519 mm += trans[i]*trans[i];
1520 sn_sdisp_.rx() = std::sqrt(mm);
1521 double crit = sn_F_.x() * fric_ + sn_cohres_ * area;
1522 if (sn_cohdc() > std::numeric_limits<double>::epsilon()) {
1523 int ww = -1;
1524 if (sn_sdisp_.x() <= sn_cohdc()) {
1525 for (int i=0; i<cohTable_.size(); ++i) {
1526 if (cohTable_[i].y() >= sn_sdisp_.x()) {
1527 ww = i;
1528 break;
1529 }
1530 }
1531 } else
1532 ww = cohTable_.size() - 1;
1533 if (ww > 0) {
1534 double prevVal = ww == 1 ? 1: cohTable_[ww-1].x() ;//critBreak_ : cohTable_[ww-1].x() + sn_F_.x() * fric_;
1535 double curVal = cohTable_[ww].x();//cohTable_[ww].x() + sn_F_.x() * fric_;
1536 double prevShear = cohTable_[ww-1].y();
1537 double curShear = cohTable_[ww].y();
1538 double slope = (curVal - prevVal)/(curShear - prevShear);
1539 // should go from 0 to 1
1540 double mval = slope * (sn_sdisp_.x() - prevShear) + prevVal;
1541 double fact = std::max(0.,std::min(1.0,1.0 - mval));
1542 assert(fact >= 0 && fact <= 1.0);
1543 double theFric = fric_;
1544 if (sn_fa_)
1545 theFric = sn_fa_ + (fric_ - sn_fa_)*fact;
1546 double theCoh = sn_Coh();
1547 if (theCoh)
1548 theCoh = sn_Coh() + (sn_cohres_ - sn_Coh())*fact;
1549 crit = sn_F_.x() * theFric + theCoh * geom.x();
1550 }
1551 }
1552
1553 // Resolve sliding.
1554 FP_S;
1555 if (crit < 0)
1556 crit = 0;
1557 DVect sforce = sn_F_; sforce.rx() = 0.0;
1558 // The is the magnitude of the shear force.
1559 double sfmag = sforce.mag();
1560 // Sliding occurs when the magnitude of the shear force is greater than the
1561 // critical value.
1562 if (sfmag > crit) {
1563 // Lower the shear force to the critical value for sliding.
1564 double rat = crit / sfmag;
1565 sforce *= rat;
1566 sforce.rx() = sn_F_.x();
1567 sn_F_ = sforce;
1568 sn_S_ = true;
1569 }
1570 if (sn_S_) {
1571 if (sn_dil_ > 0) {
1572 double zdd = sn_dilzero_ != 0 ? sn_dilzero_ : limits<double>::max();
1573 double usm = sn_sdisp_.y();
1574 if (usm < zdd) {
1575 double sInc = 0.0;
1576 for (int i=1; i<dim; ++i)
1577 sInc += trans.dof(i)*trans.dof(i);
1578 sInc = std::sqrt(sInc);
1579 sn_F_.rx() += kTran_ * sn_dil_ * sInc;
1580 }
1581 }
1582 }
1583
1584 // Resolve bending
1585 crit = sn_bmul_*2.1*0.25*br*std::abs(sn_F_.x()); // Jiang 2015
1586 FP_S;
1587 DAVect test = sn_M_;
1588#if DIM==3
1589 test.rx() = 0.0;
1590#endif
1591 double tmag = test.mag();
1592 if (tmag > crit) {
1593 // Lower the bending moment to the critical value for sliding.
1594 double rat = crit / tmag;
1595 test *= rat;
1596 sn_BS_ = true;
1597 }
1598 sn_M_.rz() = test.z();
1599#if DIM==3
1600 sn_M_.ry() = test.y();
1601 // Resolve twisting
1602 crit = sn_tmul_ * 0.65*fric_* br*std::abs(sn_F_.x()) ; // Jiang 2015
1603 tmag = std::abs(sn_M_.x());
1604 if (tmag > crit) {
1605 // Lower the shear force to the critical value for sliding.
1606 double rat = crit / tmag;
1607 tmag = sn_M_.x() * rat;
1608 sn_TS_ = true;
1609 }
1610 sn_M_.rx() = tmag;
1611 FP_S;
1612#endif
1613 }
1614 }
1615 }
1616 sn_F_old.rx() += porForce;
1617 sn_F_.rx() += porForce;
1618 totalForce = sn_F_ + fictForce_;
1619 FP_S;
1620
1621 // Account for dashpot forces if the dashpot structure has been defined.
1622 if (dpProps_) {
1623 dpProps_->dp_F_.fill(0.0);
1624 double vcn(0.0), vcs(0.0);
1625 // Calculate the damping coefficients.
1626 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1627 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1628 // First damp the shear components
1629 for (int i = 1; i<dim; ++i)
1630 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1631 // Damp the normal component
1632 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1633 // Need to change behavior based on the dp_mode.
1634 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1635 // Limit in tension if not bonded.
1636 if (sn_state_ < 3 && (dpProps_->dp_F_.x() + totalForce.x() < 0))
1637 dpProps_->dp_F_.rx() = -totalForce.rx();
1638 }
1639 if (sn_S_ && dpProps_->dp_mode_ > 1) {
1640 // Limit in shear if sliding.
1641 double dfn = dpProps_->dp_F_.rx();
1642 dpProps_->dp_F_.fill(0.0);
1643 dpProps_->dp_F_.rx() = dfn;
1644 }
1645 }
1646 FP_S;
1647
1648 //Compute energies if energy tracking has been enabled.
1649 if (state->trackEnergy_) {
1650 assert(energies_);
1651 energies_->estrain_ = 0.0;
1652 if (kn)
1653 // Calculate the strain energy.
1654 energies_->estrain_ = 0.5*totalForce.x()*totalForce.x() / kn;
1655 if (ks) {
1656 DVect s = totalForce;
1657 s.rx() = 0.0;
1658 double smag2 = s.mag2();
1659 // Add the shear component of the strain energy.
1660 energies_->estrain_ += 0.5*smag2 / ks;
1661
1662 if (sn_S_) {
1663 // If sliding calculate the slip energy and accumulate it.
1664 sn_F_old.rx() = 0.0;
1665 DVect avg_F_s = (s + sn_F_old)*0.5;
1666 DVect u_s_el = (s - sn_F_old) / ks;
1667 DVect u_s(0.0);
1668 for (int i = 1; i < dim; ++i)
1669 u_s.rdof(i) = trans.dof(i);
1670 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1671 }
1672 }
1673 // Add the bending/twisting resistance energy contributions.
1674 if (kb) {
1675 DAVect tmp = sn_M_;
1676#ifdef THREED
1677 tmp.rx() = 0.0;
1678#endif
1679 energies_->estrain_ += 0.5*tmp.mag2() / kb;
1680 if (sn_BS_) {
1681 // accumulate bending slip energy.
1682 DAVect tmp_old = sn_M_old;
1683#ifdef THREED
1684 tmp_old.rx() = 0.0;
1685#endif
1686 DAVect avg_M = (tmp + tmp_old)*0.5;
1687 DAVect t_s_el = (tmp - tmp_old) / kb;
1688 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1689 }
1690 }
1691#ifdef THREED
1692 if (kt) {
1693 double mt = std::abs(sn_M_.x());
1694 energies_->estrain_ += 0.5*mt*mt / kt;
1695 if (sn_TS_) {
1696 // accumulate twisting slip energy.
1697 DAVect tmp(0.0);
1698 DAVect tmp_old(0.0);
1699 tmp.rx() = sn_M_.x();
1700 tmp_old.rx() = sn_M_old.x();
1701 DAVect avg_M = (tmp + tmp_old)*0.5;
1702 DAVect t_s_el = (tmp - tmp_old) / kt;
1703 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1704 }
1705 }
1706#endif
1707
1708 if (dpProps_) {
1709 // Calculate damping energy (accumulated) if the dashpots are active.
1710 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1711 }
1712 }
1713
1714 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1715 assert(sn_F_ == sn_F_);
1716 assert(sn_M_ == sn_M_);
1717 assert(fictForce_ == fictForce_);
1718 FP_S;
1719 return true;
1720 }
1721
1722 bool ContactModelRBSN::FDLawUnBonded(ContactModelMechanicalState *state, const double ×tep) {
1723 DVect3 geom = computeGeomData(state->end1Curvature_,state->end2Curvature_);
1724 double br = geom.z();
1725 // Relative translational/rotational displacement increments
1726 DVect trans = state->relativeTranslationalIncrement_;
1727 DAVect ang = state->relativeAngularIncrement_;
1728 double overlap = rgap_ - state->gap_;
1729 double correction = 1.0;
1730 if (state->activated() && sn_mode_ == 0 && trans.x()) {
1731 correction = -1.0*overlap / trans.x();
1732 if (correction < 0)
1733 correction = 1.0;
1734 }
1735
1736 // Store previous force and moment
1737 DVect sn_F_old = sn_F_;
1738 double porForce = sn_por_ * geom.x();
1739 sn_F_old.rx() -= porForce;
1740 DAVect sn_M_old = sn_M_;
1741
1742 double kb = kRot_.z();
1743#if DIM==3
1744 double kt = kRot_.x();
1745 //kb = std::sqrt(kb * kb + kRot_.y() * kRot_.y());
1746#endif
1747 // absolute/incremental normal force calculation
1748 DVect theStiff(kTran_/kRatio_);
1749 theStiff.rx() = kTran_;
1750
1751 if (sn_mode_==0)
1752 sn_F_.rx() = overlap * theStiff.x();
1753 else
1754 sn_F_.rx() -= trans.x() * theStiff.x();
1755 // Normal force can only be positive if unbonded
1756 sn_F_.rx() = std::max(0.0, sn_F_.x()) - porForce;
1757
1758 // Calculate the trial shear force.
1759 DVect sforce(0.0);
1760 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1761 // Loop over the shear components (note: the 0 component is the normal component)
1762 // and calculate the shear force.
1763 for (int i = 1; i<dim; ++i)
1764 sforce.rdof(i) = sn_F_.dof(i) - trans.dof(i) * theStiff.dof(i);
1765
1766 // Calculate the trial moment.
1767 DAVect mom = sn_M_ - ang*kRot_;
1768
1769 // If the SOLVE ELASTIC command is given then the
1770 // canFail state is set to FALSE. Otherwise it is always TRUE.
1771 if (state->canFail_) {
1772 bool changed = false;
1773 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
1774 bool slip_changed = false;
1775
1776 double crit = sn_F_.x() * fric_ + sn_cohres_ * geom.x();
1777 if (sn_state_ != 0) {
1778 bool doComp = true;
1779 if (!sn_S_) {
1780 if (sn_heal_) {
1781 sn_sdisp_.rx() = 0;
1782 crit = sn_F_.x() * sn_fa_ + cohTable_[0].x() * geom.x();
1783 doComp = false;
1784 }
1785 }
1786 if (doComp) {
1787 double mm = 0.0;
1788 for (int i=1; i<dim; ++i)
1789 mm += trans[i]*trans[i];
1790 sn_sdisp_.rx() += std::sqrt(mm);
1791 if (sn_cohdc() > std::numeric_limits<double>::epsilon()) {
1792 int ww = -1;
1793 if (sn_sdisp_.x() <= sn_cohdc()) {
1794 for (int i=0; i<cohTable_.size(); ++i) {
1795 if (cohTable_[i].y() >= sn_sdisp_.x()) {
1796 ww = i;
1797 break;
1798 }
1799 }
1800 } else
1801 ww = cohTable_.size() - 1;
1802 if (ww > 0) {
1803 double prevVal = ww == 1 ? 1: cohTable_[ww-1].x() ;//critBreak_ : cohTable_[ww-1].x() + sn_F_.x() * fric_;
1804 double curVal = cohTable_[ww].x();//cohTable_[ww].x() + sn_F_.x() * fric_;
1805 double prevShear = cohTable_[ww-1].y();
1806 double curShear = cohTable_[ww].y();
1807 double slope = (curVal - prevVal)/(curShear - prevShear);
1808 // should go from 0 to 1
1809 double mval = slope * (sn_sdisp_.x() - prevShear) + prevVal;
1810 double fact = std::max(0.,std::min(1.0,1.0 - mval));
1811 assert(fact >= 0 && fact <= 1.0);
1812 double theFric = fric_;
1813 if (sn_fa_)
1814 theFric = sn_fa_ + (fric_ - sn_fa_)*fact;
1815 double theCoh = sn_Coh();
1816 if (theCoh)
1817 theCoh = sn_Coh() + (sn_cohres_ - sn_Coh())*fact;
1818 crit = sn_F_.x() * theFric + theCoh * geom.x();
1819 }
1820 }
1821 }
1822 }
1823
1824 // Resolve sliding.
1825 if (crit < 0)
1826 crit = 0.0;
1827 // The is the magnitude of the shear force.
1828 double sfmag = sforce.mag();
1829 if (sfmag > crit) {
1830 // Lower the shear force to the critical value for sliding.
1831 double rat = crit / sfmag;
1832 sforce *= rat;
1833 if (!sn_S_) {
1834 slip_changed = true;
1835 changed = true;
1836 }
1837 sn_S_ = true;
1838 } else {
1839 if (sn_S_) {
1840 slip_changed = true;
1841 changed = true;
1842 }
1843 sn_S_ = false;
1844 }
1845 if (sn_S_) {
1846 if (sn_dil_ > 0) {
1847 double zdd = sn_dilzero_ != 0 ? sn_dilzero_ : limits<double>::max();
1848 double usm = sn_sdisp_.y();
1849 if (usm < zdd) {
1850 double sInc = 0.0;
1851 for (int i=1; i<dim; ++i)
1852 sInc += trans.dof(i)*trans.dof(i);
1853 sInc = std::sqrt(sInc);
1854 sn_F_.rx() += kTran_ * sn_dil_ * sInc;
1855 }
1856 }
1857 } else {
1858 if (sn_heal_) {
1859 if (shearStrength(geom.x()))
1860 sn_state_ = 6;
1861 }
1862 }
1863 // Resolve bending
1864 bool bslip_changed = false;
1865 crit = sn_bmul_ * 2.1*0.25*sn_F_.x() * br; // Jiang 2015
1866 DAVect test = mom;
1867#if DIM==3
1868 test.rx() = 0.0;
1869#endif
1870 double tmag = test.mag();
1871 if (tmag > crit) {
1872 // Lower the bending moment to the critical value for sliding.
1873 double rat = crit / tmag;
1874 test *= rat;
1875 if (!sn_BS_) {
1876 bslip_changed = true;
1877 changed = true;
1878 }
1879 sn_BS_ = true;
1880 }
1881 else {
1882 if (sn_BS_) {
1883 bslip_changed = true;
1884 changed = true;
1885 }
1886 sn_BS_ = false;
1887 }
1888 mom.rz() = test.z();
1889#if DIM==3
1890 mom.ry() = test.y();
1891 // Resolve twisting
1892 bool tslip_changed = false;
1893 crit = sn_tmul_ * 0.65*fric_*sn_F_.x() * br; // Jiang 2015
1894 tmag = std::abs(mom.x());
1895 if (tmag > crit) {
1896 // Lower the twisting moment to the critical value for sliding.
1897 double rat = crit / tmag;
1898 mom.rx() *= rat;
1899 if (!sn_TS_) {
1900 tslip_changed = true;
1901 changed = true;
1902 }
1903 sn_TS_ = true;
1904 }
1905 else {
1906 if (sn_TS_) {
1907 tslip_changed = true;
1908 changed = true;
1909 }
1910 sn_TS_ = false;
1911 }
1912#endif
1913 if (changed && cmEvents_[fSlipChange] >= 0) {
1914 qint64 code = 0;
1915 if (slip_changed) {
1916 code = 1;
1917 if (bslip_changed) {
1918 code = 4;
1919#if DIM==3
1920 if (tslip_changed)
1921 code = 7;
1922#endif
1923 }
1924 }
1925 else if (bslip_changed) {
1926 code = 2;
1927#if DIM==3
1928 if (tslip_changed)
1929 code = 6;
1930#endif
1931 }
1932#if DIM==3
1933 else if (tslip_changed) {
1934 code = 3;
1935 if (slip_changed)
1936 code = 5;
1937 }
1938#endif
1939 auto c = state->getContact();
1940 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1941 fish::Parameter(code),
1942 fish::Parameter(sn_S_),
1943 fish::Parameter(sn_BS_)
1944#ifdef THREED
1945 ,fish::Parameter(sn_TS_)
1946#endif
1947 };
1948 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1949 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1950 }
1951 }
1952
1953 sn_F_.rx() += porForce;
1954 sn_F_old.rx() += porForce;
1955
1956 // Set the shear components of the total force.
1957 for (int i = 1; i<dim; ++i)
1958 sn_F_.rdof(i) = sforce.dof(i);
1959
1960 // Set the moment.
1961 sn_M_ = mom;
1962
1963 // Account for dashpot forces if the dashpot structure has been defined.
1964 if (dpProps_) {
1965 dpProps_->dp_F_.fill(0.0);
1966 double vcn(0.0), vcs(0.0);
1967 // Calculate the damping coefficients.
1968 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kTran_)));
1969 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(kTran_/kRatio_)));
1970 // First damp the shear components
1971 for (int i = 1; i<dim; ++i)
1972 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1973 // Damp the normal component
1974 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1975 // Need to change behavior based on the dp_mode.
1976 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1977 // Limit in tension if not bonded.
1978 if (dpProps_->dp_F_.x() + sn_F_.x() < 0)
1979 dpProps_->dp_F_.rx() = -sn_F_.rx();
1980 }
1981 if (sn_S_ && dpProps_->dp_mode_ > 1) {
1982 // Limit in shear if not sliding.
1983 double dfn = dpProps_->dp_F_.rx();
1984 dpProps_->dp_F_.fill(0.0);
1985 dpProps_->dp_F_.rx() = dfn;
1986 }
1987 }
1988
1989 //Compute energies if energy tracking has been enabled.
1990 if (state->trackEnergy_) {
1991 assert(energies_);
1992 energies_->estrain_ = 0.0;
1993 if (kTran_)
1994 // Calculate the strain energy.
1995 energies_->estrain_ = 0.5*sn_F_.x()*sn_F_.x() / kTran_;
1996 if (kTran_) {
1997 DVect s = sn_F_;
1998 s.rx() = 0.0;
1999 double smag2 = s.mag2();
2000 // Add the shear component of the strain energy.
2001 energies_->estrain_ += 0.5*smag2 / (kTran_/kRatio_);
2002
2003 if (sn_S_) {
2004 // If sliding calculate the slip energy and accumulate it.
2005 sn_F_old.rx() = 0.0;
2006 DVect avg_F_s = (s + sn_F_old)*0.5;
2007 DVect u_s_el = (s - sn_F_old) / (kTran_/kRatio_);
2008 DVect u_s(0.0);
2009 for (int i = 1; i < dim; ++i)
2010 u_s.rdof(i) = trans.dof(i);
2011 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
2012 }
2013 }
2014 // Add the bending/twisting resistance energy contributions.
2015 if (kb) {
2016 DAVect tmp = sn_M_;
2017#ifdef THREED
2018 tmp.rx() = 0.0;
2019#endif
2020 energies_->estrain_ += 0.5*tmp.mag2() / kb;
2021 if (sn_BS_) {
2022 // accumulate bending slip energy.
2023 DAVect tmp_old = sn_M_old;
2024#ifdef THREED
2025 tmp_old.rx() = 0.0;
2026#endif
2027 DAVect avg_M = (tmp + tmp_old)*0.5;
2028 DAVect t_s_el = (tmp - tmp_old) / kb;
2029 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
2030 }
2031 }
2032#ifdef THREED
2033 if (kt) {
2034 double mt = std::abs(sn_M_.x());
2035 energies_->estrain_ += 0.5*mt*mt / kt;
2036 if (sn_TS_) {
2037 // accumulate twisting slip energy.
2038 DAVect tmp(0.0);
2039 DAVect tmp_old(0.0);
2040 tmp.rx() = sn_M_.x();
2041 tmp_old.rx() = sn_M_old.x();
2042 DAVect avg_M = (tmp + tmp_old)*0.5;
2043 DAVect t_s_el = (tmp - tmp_old) / kt;
2044 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
2045 }
2046 }
2047#endif
2048
2049 if (dpProps_) {
2050 // Calculate damping energy (accumulated) if the dashpots are active.
2051 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
2052 }
2053 }
2054
2055 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
2056 assert(sn_F_ == sn_F_);
2057 assert(fictForce_ == fictForce_);
2058 assert(sn_M_ == sn_M_);
2059 return true;
2060 }
2061
2062 void ContactModelRBSN::setForce(const DVect &v,IContact *c) {
2063 sn_F_ = v;
2064 fictForce_ = DVect(0.0);
2065 forceSet_ = true;
2066 if (v.x() > 0)
2067 rgap_ = c->getGap() + v.x() / kTran_;
2068 }
2069
2070 void ContactModelRBSN::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
2071 // Only called for contacts with wall facets when the wall resolution scheme
2072 // is set to full!
2073 // Only do something if the contact model is of the same type
2074 if (old->getContactModel()->getName().compare("springnetwork",Qt::CaseInsensitive) == 0 && !isBonded()) {
2075 ContactModelRBSN *oldCm = (ContactModelRBSN *)old;
2076#ifdef THREED
2077 // Need to rotate just the shear component from oldSystem to newSystem
2078
2079 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
2080 DVect axis = oldSystem.e1() & newSystem.e1();
2081 double c, ang, s;
2082 DVect re2;
2083 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
2084 axis = axis.unit();
2085 c = oldSystem.e1()|newSystem.e1();
2086 if (c > 0)
2087 c = std::min(c,1.0);
2088 else
2089 c = std::max(c,-1.0);
2090 ang = acos(c);
2091 s = sin(ang);
2092 double t = 1. - c;
2093 DMatrix<3,3> rm;
2094 rm.get(0,0) = t*axis.x()*axis.x() + c;
2095 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
2096 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
2097 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
2098 rm.get(1,1) = t*axis.y()*axis.y() + c;
2099 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
2100 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
2101 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
2102 rm.get(2,2) = t*axis.z()*axis.z() + c;
2103 re2 = rm*oldSystem.e2();
2104 }
2105 else
2106 re2 = oldSystem.e2();
2107 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
2108 axis = re2 & newSystem.e2();
2109 DVect2 tpf;
2110 DVect2 tpm;
2111 DMatrix<2,2> m;
2112 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
2113 axis = axis.unit();
2114 c = re2|newSystem.e2();
2115 if (c > 0)
2116 c = std::min(c,1.0);
2117 else
2118 c = std::max(c,-1.0);
2119 ang = acos(c);
2120 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
2121 ang *= -1;
2122 s = sin(ang);
2123 m.get(0,0) = c;
2124 m.get(1,0) = s;
2125 m.get(0,1) = -m.get(1,0);
2126 m.get(1,1) = m.get(0,0);
2127 tpf = m*DVect2(oldCm->sn_F_.y(),oldCm->sn_F_.z());
2128 tpm = m*DVect2(oldCm->sn_M_.y(),oldCm->sn_M_.z());
2129 } else {
2130 m.get(0,0) = 1.;
2131 m.get(0,1) = 0.;
2132 m.get(1,0) = 0.;
2133 m.get(1,1) = 1.;
2134 tpf = DVect2(oldCm->sn_F_.y(),oldCm->sn_F_.z());
2135 tpm = DVect2(oldCm->sn_M_.y(),oldCm->sn_M_.z());
2136 }
2137 DVect pforce = DVect(0,tpf.x(),tpf.y());
2138 //DVect pm = DVect(0,tpm.x(),tpm.y());
2139#else
2140 oldSystem;
2141 newSystem;
2142 DVect pforce = DVect(0,oldCm->sn_F_.y());
2143 //DVect pm = DVect(0,oldCm->sn_M_.y());
2144#endif
2145 for (int i=1; i<dim; ++i)
2146 sn_F_.rdof(i) += pforce.dof(i);
2147 if (sn_mode_ && oldCm->sn_mode_)
2148 sn_F_.rx() = oldCm->sn_F_.x();
2149 oldCm->sn_F_ = DVect(0.0);
2150 oldCm->sn_M_ = DAVect(0.0);
2151 if (dpProps_ && oldCm->dpProps_) {
2152#ifdef THREED
2153 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
2154 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
2155#else
2156 pforce = oldCm->dpProps_->dp_F_;
2157#endif
2158 dpProps_->dp_F_ += pforce;
2159 oldCm->dpProps_->dp_F_ = DVect(0.0);
2160 }
2161 if(oldCm->getEnergyActivated()) {
2162 activateEnergy();
2163 energies_->estrain_ = oldCm->energies_->estrain_;
2164 energies_->edashpot_ = oldCm->energies_->edashpot_;
2165 energies_->eslip_ = oldCm->energies_->eslip_;
2166 oldCm->energies_->estrain_ = 0.0;
2167 oldCm->energies_->edashpot_ = 0.0;
2168 oldCm->energies_->eslip_ = 0.0;
2169 }
2170 rgap_ = oldCm->rgap_;
2171 }
2172 assert(sn_F_ == sn_F_);
2173 }
2174
2175 void ContactModelRBSN::setNonForcePropsFrom(IContactModel *old) {
2176 // Only called for contacts with wall facets when the wall resolution scheme
2177 // is set to full!
2178 // Only do something if the contact model is of the same type
2179 if (old->getName().compare("springnetwork",Qt::CaseInsensitive) == 0 && !isBonded()) {
2180 ContactModelRBSN *oldCm = (ContactModelRBSN *)old;
2181
2182 fictForce_ = oldCm->fictForce_;
2183 sn_F_ = oldCm->sn_F_;
2184 sn_sdisp_ = oldCm->sn_sdisp_;
2185 sn_M_ = oldCm->sn_M_;
2186 kRot_ = oldCm->kRot_;
2187 kTran_ = oldCm->kTran_;
2188 kRatio_ = oldCm->kRatio_;
2189 E_ = oldCm->E_;
2190 poisson_ = oldCm->poisson_;
2191 fric_ = oldCm->fric_;
2192 sn_bmul_ = oldCm->sn_bmul_;
2193 sn_tmul_ = oldCm->sn_tmul_;
2194 sn_rmul_ = oldCm->sn_rmul_;
2195 userArea_ = oldCm->userArea_;
2196 rgap_ = oldCm->rgap_;
2197 sn_fa_ = oldCm->sn_fa_;
2198 sn_mcf_ = oldCm->sn_mcf_;
2199 sn_dil_ = oldCm->sn_dil_;
2200 sn_dilzero_ = oldCm->sn_dilzero_;
2201 transTen_ = oldCm->transTen_;
2202 sn_elong_ = oldCm->sn_elong_;
2203 sn_ndisp_ = oldCm->sn_ndisp_;
2204 sn_mode_ = oldCm->sn_mode_;
2205 sn_state_ = oldCm->sn_state_;
2206 poisOffDiag_ = oldCm->poisOffDiag_;
2207 sn_S_ = oldCm->sn_S_;
2208 sn_BS_ = oldCm->sn_BS_;
2209 sn_TS_ = oldCm->sn_TS_;
2210 forceSet_ = oldCm->forceSet_;
2211 sn_heal_ = oldCm->sn_heal_;
2212 tenTable_ = oldCm->tenTable_;
2213 cohTable_ = oldCm->cohTable_;
2214 if (oldCm->hasDamping()) {
2215 if (!dpProps_)
2216 dpProps_ = NEWC(dpProps());
2217 dp_nratio(oldCm->dp_nratio());
2218 dp_sratio(oldCm->dp_sratio());
2219 dp_mode(oldCm->dp_mode());
2220 dp_F(oldCm->dp_F());
2221 }
2222 }
2223 }
2224
2225 DVect ContactModelRBSN::getForce() const {
2226 DVect ret(sn_F_);
2227 ret += fictForce_;
2228 if (dpProps_)
2229 ret += dpProps_->dp_F_;
2230 return ret;
2231 }
2232
2233 DAVect ContactModelRBSN::getMomentOn1(const IContactMechanical *c) const {
2234 DAVect ret(sn_M_);
2235 if (c) {
2236 DVect force = getForce();
2237 c->updateResultingTorqueOn1Local(force,&ret);
2238 }
2239 return ret;
2240 }
2241
2242 DAVect ContactModelRBSN::getMomentOn2(const IContactMechanical *c) const {
2243 DAVect ret(sn_M_);
2244 if (c) {
2245 DVect force = getForce();
2246 c->updateResultingTorqueOn2Local(force,&ret);
2247 }
2248 return ret;
2249 }
2250
2251 DAVect ContactModelRBSN::getModelMomentOn1() const {
2252 DAVect ret(sn_M_);
2253 return ret;
2254 }
2255
2256 DAVect ContactModelRBSN::getModelMomentOn2() const {
2257 DAVect ret(sn_M_);
2258 return ret;
2259 }
2260
2261 void ContactModelRBSN::objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *ret) const {
2262 ret->clear();
2263 ret->push_back({"kn",ScalarInfo});
2264 ret->push_back({"ks",ScalarInfo});
2265 ret->push_back({"krot",VectorInfo});
2266 ret->push_back({"fric",ScalarInfo});
2267 ret->push_back({"sn_bmul",ScalarInfo});
2268 ret->push_back({"sn_tmul",ScalarInfo});
2269 ret->push_back({"sn_mode",ScalarInfo});
2270 ret->push_back({"sn_force",VectorInfo});
2271 ret->push_back({"sn_moment",VectorInfo});
2272 ret->push_back({"sn_slip",ScalarInfo});
2273 ret->push_back({"sn_slipb",ScalarInfo});
2274 ret->push_back({"sn_slipt",ScalarInfo});
2275 ret->push_back({"sn_rmul",ScalarInfo});
2276 ret->push_back({"sn_radius",ScalarInfo});
2277 ret->push_back({"emod",ScalarInfo});
2278 ret->push_back({"kratio",ScalarInfo});
2279 ret->push_back({"dp_nratio",ScalarInfo});
2280 ret->push_back({"dp_sratio",ScalarInfo});
2281 ret->push_back({"dp_mode",ScalarInfo});
2282 ret->push_back({"dp_force",VectorInfo});
2283 ret->push_back({"sn_state",ScalarInfo});
2284 ret->push_back({"sn_ten",ScalarInfo});
2285 ret->push_back({"sn_shear",ScalarInfo});
2286 ret->push_back({"sn_coh",ScalarInfo});
2287 ret->push_back({"sn_fa",ScalarInfo});
2288 ret->push_back({"sn_mcf",ScalarInfo});
2289 ret->push_back({"sn_sigma",ScalarInfo});
2290 ret->push_back({"sn_tau",ScalarInfo});
2291 ret->push_back({"sn_area",ScalarInfo});
2292 ret->push_back({"user_area",ScalarInfo});
2293 ret->push_back({"rgap",ScalarInfo});
2294 ret->push_back({"sn_pois_force",VectorInfo});
2295 ret->push_back({"sn_pois",ScalarInfo});
2296 ret->push_back({"sn_non_diag",ScalarInfo});
2297 ret->push_back({"sn_cohres",ScalarInfo});
2298 ret->push_back({"sn_dil",ScalarInfo});
2299 ret->push_back({"sn_dilzero",ScalarInfo});
2300 ret->push_back({"sn_ndisp",ScalarInfo});
2301 ret->push_back({"sn_sdisp",VectorInfo});
2302 ret->push_back({"sn_cohdc",ScalarInfo});
2303 ret->push_back({"sn_heal",ScalarInfo});
2304 ret->push_back({"sn_tendc",ScalarInfo});
2305 ret->push_back({"sn_tabpos",ScalarInfo});
2306 ret->push_back({"sn_porp",ScalarInfo});
2307 ret->push_back({"sn_esigma",ScalarInfo});
2308
2309 }
2310
2311 void ContactModelRBSN::objectPropValues(std::vector<double> *ret,const IContact *con) const {
2312 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
2313 ret->clear();
2314 ret->push_back(kTran_);
2315 ret->push_back(kTran_/kRatio_);
2316 ret->push_back(kRot_.mag());
2317 ret->push_back(fric_);
2318 ret->push_back(sn_bmul_);
2319 ret->push_back(sn_tmul_);
2320 ret->push_back(sn_mode_);
2321 ret->push_back(sn_F_.mag());
2322 ret->push_back(sn_M_.mag());
2323 ret->push_back(sn_S_);
2324 ret->push_back(sn_BS_);
2325 ret->push_back(sn_TS_);
2326 ret->push_back(sn_rmul_);
2327 double Cmax1 = c->getEnd1Curvature().y();
2328 double Cmax2 = c->getEnd2Curvature().y();
2329 double d;
2330 if (!userArea_)
2331 d= sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
2332 else {
2333#ifdef THREED
2334 d= std::sqrt(userArea_ / dPi);
2335#else
2336 d= userArea_ / 2.0;
2337#endif
2338 }
2339 ret->push_back(d);
2340 ret->push_back(E_);
2341 ret->push_back(kRatio_);
2342 ret->push_back(dpProps_ ? dpProps_->dp_nratio_ : 0);
2343 ret->push_back(dpProps_ ? dpProps_->dp_sratio_ : 0);
2344 ret->push_back(dpProps_ ? dpProps_->dp_mode_ : 0);
2345 ret->push_back(dpProps_ ? dpProps_->dp_F_.mag() : 0);
2346 ret->push_back(sn_state_);
2347 ret->push_back(sn_Ten());
2348 ret->push_back(shearStrength(computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x()));
2349 ret->push_back(cohTable_[0].x());
2350 ret->push_back(std::atan(sn_fa_)/dDegrad);
2351 ret->push_back(sn_mcf_);
2352 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x());
2353 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).y());
2354 ret->push_back(userArea_ ? userArea_ : computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x());
2355 ret->push_back(userArea_);
2356 ret->push_back(rgap_);
2357 ret->push_back(fictForce_.mag());
2358 ret->push_back(poisson_);
2359 ret->push_back(poisOffDiag_ == false ? 0 : 1);
2360 ret->push_back(sn_cohres_);
2361 ret->push_back(std::atan(sn_dil_)/dDegrad);
2362 ret->push_back(sn_dilzero_);
2363 ret->push_back(sn_ndisp_);
2364 ret->push_back(sn_sdisp_.mag());
2365 ret->push_back(sn_cohdc());
2366 ret->push_back(sn_heal_);
2367 ret->push_back(sn_tendc());
2368 ret->push_back(sn_tabPos_+1);
2369 ret->push_back(sn_por_);
2370 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x() + sn_por_);
2371 }
2372
2373 DVect3 ContactModelRBSN::computeGeomData(const DVect2 &end1Curvature,const DVect2 &end2Curvature) const {
2374 double Cmax1 = end1Curvature.y();
2375 double Cmax2 = end2Curvature.y();
2376 double br = sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
2377 if (userArea_)
2378#ifdef THREED
2379 br = std::sqrt(userArea_ / dPi);
2380#else
2381 br = userArea_ / 2.0;
2382#endif
2383 double br2 = br * br;
2384#ifdef THREED
2385 double area = dPi * br2;
2386 double bi = 0.25*area*br2;
2387#else
2388 double area = 2.0*br;
2389 double bi = 2.0*br*br2 / 3.0;
2390#endif
2391 return DVect3(area, bi, br);
2392 }
2393
2394 DVect2 ContactModelRBSN::SMax(const DVect2 &e1c,const DVect2 &e2c) const {
2395 DVect3 data = computeGeomData(e1c,e2c);
2396 double area = data.x();
2397 double bi = data.y();
2398 double br = data.z();
2399 /* maximum stresses */
2400 //double nsmax0 = -(totalForce.x() / area) + sn_mcf_* sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z()) * br / bi;
2401 double dbend = sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z());
2402 dbend *= sn_mcf_;
2403 double dtwist = sn_M_.x();
2404 DVect bfs(sn_F_);
2405 bfs.rx() = 0.0;
2406 double dbfs = bfs.mag();
2407 double nsmax = -((sn_F_.x()+fictForce_.x()) / area) + dbend * br / bi;
2408 double ssmax = dbfs / area + std::abs(dtwist) * 0.5* br / bi;
2409 return DVect2(nsmax, ssmax);
2410 }
2411
2412 double ContactModelRBSN::shearStrength(const double &area) const {
2413 double sig = -1.0*(sn_F_.x() + fictForce_.x()) / area;
2414 double nstr = (sn_state_ > 2 && sn_state_ != 6) ? tenTable_[0].x() : 0.0;
2415 return sig <= nstr ? cohTable_[0].x() - sn_fa_*sig
2416 : cohTable_[0].x() - sn_fa_*nstr;
2417 }
2418
2419
2420 double ContactModelRBSN::strainEnergy(double kn,double ks,double kb,double kt) const {
2421 double ret(0.0);
2422 if (kn)
2423 ret = 0.5 * (sn_F_.x()+fictForce_.x()) * (sn_F_.x()+fictForce_.x()) / kn;
2424 if (ks) {
2425 DVect tmp = sn_F_ + fictForce_;
2426 tmp.rx() = 0.0;
2427 double smag2 = tmp.mag2();
2428 ret += 0.5 * smag2 / ks;
2429 }
2430
2431 if (kt)
2432 ret += 0.5 * sn_M_.x() * sn_M_.x() / kt;
2433 if (kb) {
2434 DAVect tmp = sn_M_;
2435#ifdef THREED
2436 tmp.rx() = 0.0;
2437 double smag2 = tmp.mag2();
2438#else
2439 double smag2 = tmp.z() * tmp.z();
2440#endif
2441 ret += 0.5 * smag2 / kb;
2442 }
2443 return ret;
2444 }
2445
2446} // namespace cmodelsxd
2447// EoF
| Was this helpful? ... | Itasca Software © 2026, Itasca | Updated: Oct 23, 2026 |
