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 string 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 string 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 string 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_ = NEW 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 string getFishCallEvents() const override {
180 return
181 "contact_activated"
182 ",slip_change"
183 ",bond_break";
184 }
185
186 // Return the specified contact model property.
187 base::Property 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 base::Property &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 string 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 string 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 std::vector<base::Property> &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 string &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 NEW 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<string,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 uint32 sn_mode_ = 0; // specifies absolute (0) or incremental (1) behavior for the the normal force
461 uint32 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 = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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 base::Property 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 base::Property 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 base::Property();
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 base::Property &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<int64>())
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_ = NEW 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_ = NEW dpProps();
583 dpProps_->dp_sratio_ = val;
584 return true;
585 }
586 case kwDpMode: {
587 if (!v.canConvert<int64>())
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_ = NEW 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_ = NEW 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<int64>())
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 string 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 string();
862 }
863 assert(0);
864 return string();
865 }
866
867 bool ContactModelRBSN::setMethod(uint32 i,const std::vector<base::Property> &vl,IContact *con) {
868 FP_S;
869 // Apply the specified method.
870 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
871 switch (i) {
872 case kwAssignStiffness: {
873 poisson_ = 0.0;
874 if (vl.at(0).isNull())
875 throw Exception("Argument kn must be specified with method assign-stiffness in contact model {0}.",getName());
876 double knpa = vl.at(0).toDouble();
877 if (knpa<=0.0)
878 throw Exception("Negative or zero kn not allowed in contact model {0}.",getName());
879 if (vl.at(1).isNull())
880 throw Exception("Argument kratio must be specified with method assign-stiffness in contact model {0}.",getName());
881 kRatio_ = vl.at(1).toDouble();
882 if (kRatio_<0.0) {
883 kRatio_ = 0.0;
884 throw Exception("Negative kratio not allowed in contact model {0}.",getName());
885 }
886 std::vector<DVect> pts;
887 c->getJointGeometry(&pts);
888 double area = 0.0;
889#ifdef THREED
890 // Step 1: get centroid and area
891 for (int i=1; i<pts.size()-1; ++i) {
892 double a = (pts[0] - pts[i]).mag();
893 double b = (pts[0] - pts[i+1]).mag();
894 double c = (pts[i] - pts[i+1]).mag();
895 double la = 0.0;
896 if (b > a)
897 std::swap(a,b);
898 if (c > a)
899 std::swap(a,c);
900 if (c > b)
901 std::swap(b,c);
902 if (c - (a - b) >= 0)
903 la = 0.25 * sqrt((a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c)));
904 area += la;
905 }
906#else
907 // Assume unit thickness in the out of plane direction
908 area = (pts[1] - pts[0]).mag();
909#endif
910 userArea_ = area;
911 kTran_ = knpa * area;
912 E_ = kTran_ / area;
913 kRot_ = DAVect(0.0);
914 setInheritance(1,false);
915 setInheritance(2,false);
916 sn_mode_ = 1.0;
917 FP_S;
918 return true;
919 }
920 case kwStiffness: {
921 FP_S;
922 poisson_ = 0.0;
923 if (vl.at(0).isNull())
924 throw Exception("Argument emod must be specified with method compute-stiffness in contact model {0}.",getName());
925 E_ = vl.at(0).toDouble();
926 if (E_<=0.0)
927 throw Exception("Negative or zero emod not allowed in contact model {0}.",getName());
928 if (vl.at(1).isNull())
929 throw Exception("Argument poisson must be specified with method compute-stiffness in contact model {0}.",getName());
930 poisson_ = vl.at(1).toDouble();
931 if (poisson_ < 0.0) {
932 poisson_ = 0.0;
933 throw Exception("Negative poisson not allowed in contact model {0}.",getName());
934 }
935 const IBody *b1 = con->getEnd1()->getIBody();
936 const IBody *b2 = con->getEnd2()->getIBody();
937 //double vol1 = b1->getVolume();
938 //double vol2 = b2->getVolume();
939 //if (std::max(vol1,vol2) > 10.*std::min(vol1,vol2))
940 // poisson_ = 0.0;
941 DVect pos1 = toDVect(b1->getIThing()->getLocation());
942 DVect pos2 = toDVect(b2->getIThing()->getLocation()) + con->getOffSet();
943 double dist = (pos1-pos2).mag();
944 if (con->withWall())
945 dist = (pos1 - con->getPosition()).mag();
946 double tol = std::max(pos1.abs().maxComp(),pos2.abs().maxComp())*limits<double>::epsilon()*1000;
947 if (dist < tol) {
948 poisson_ = 0;
949 userArea_ = 0;
950 kTran_ = 0;
951 kRot_.fill(0);
952 FP_S;
953 return true;
954 }
955 std::vector<DVect> pts;
956 FP_S;
957 c->getJointGeometry(&pts);
958 FP_S;
959 double area = 0.0;
960 DAVect inertia(0.0);
961#ifdef THREED
962 // Step 1: get centroid and area
963 DVect cm(0.0);
964 for (int i=1; i<pts.size()-1; ++i) {
965 DVect lcm = (pts[0] + pts[i] + pts[i+1])/3.0;
966 double a = (pts[0] - pts[i]).mag();
967 double b = (pts[0] - pts[i+1]).mag();
968 double c = (pts[i] - pts[i+1]).mag();
969 double la = 0.0;
970 if (b > a)
971 std::swap(a,b);
972 if (c > a)
973 std::swap(a,c);
974 if (c > b)
975 std::swap(b,c);
976 if (c - (a - b) >= 0)
977 la = 0.25 * sqrt((a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c)));
978 cm += lcm * la;
979 area += la;
980 }
981 FP_S;
982 if (area == 0.0) {
983 poisson_ = 0;
984 userArea_ = 0;
985 kTran_ = 0;
986 kRot_.fill(0);
987 return true;
988 }
989 cm /= area;
990 FP_S;
991 // Step 2 - center it and put in the local system
992 for (int i=0; i<pts.size(); ++i) {
993 pts[i] -= cm;
994 pts[i] = con->getLocalSystem().toLocal(pts[i]);
995 }
996 // Step 3: compute the polar inertia
997 for (int i=0; i<pts.size(); ++i) {
998 int j = i < pts.size()-1 ? i+1 : 0;
999 double xi = pts[i].y();
1000 double xip1 = pts[j].y();
1001 double yi = pts[i].z();
1002 double yip1 = pts[j].z();
1003 double frnt = (xi*yip1-xip1*yi);
1004 inertia.ry() += frnt*(xi*xi+xi*xip1+xip1*xip1);
1005 inertia.rz() += frnt*(yi*yi+yi*yip1+yip1*yip1);
1006 }
1007 inertia.ry() = std::abs(inertia.y() / 12.);
1008 inertia.rz() = std::abs(inertia.z() / 12.);
1009 inertia.rx() = inertia.y() + inertia.z();
1010
1011#else
1012 // Assume unit thickness in the out of plane direction
1013 area = (pts[1] - pts[0]).mag();
1014 inertia.rz() = area*area*area/12.;
1015#endif
1016 userArea_ = area;
1017 kTran_ = E_ * area / dist;
1018 kRot_ = inertia *E_ / dist;
1019 setInheritance(1,false);
1020 setInheritance(2,false);
1021 sn_mode_ = 1.0;
1022 FP_S;
1023 return true;
1024 }
1025 case kwBond: {
1026 if (sn_state_ == 3) return false;
1027 double mingap = -1.0 * limits<double>::max();
1028 double maxgap = 0;
1029 if (vl.at(0).canConvert<double>())
1030 maxgap = vl.at(0).toDouble();
1031 else if (vl.at(0).canConvert<DVect2>()) {
1032 DVect2 value = vl.at(0).value<DVect2>();
1033 mingap = value.minComp();
1034 maxgap = value.maxComp();
1035 }
1036 else if (!vl.at(0).isNull())
1037 throw Exception("gap value {0} not recognized in method bond in contact model {1}.", vl.at(1), getName());
1038 double gap = c->getGap();
1039 if (gap >= mingap && gap <= maxgap) {
1040 sn_state_ = 3;
1041 sn_mode_ = 1;
1042 FP_S;
1043 return true;
1044 }
1045 FP_S;
1046 return false;
1047 }
1048 case kwUnbond: {
1049 if (sn_state_ == 0) return false;
1050 double mingap = -1.0 * limits<double>::max();
1051 double maxgap = 0;
1052 if (vl.at(0).canConvert<double>())
1053 maxgap = vl.at(0).toDouble();
1054 else if (vl.at(0).canConvert<DVect2>()) {
1055 DVect2 value = vl.at(0).value<DVect2>();
1056 mingap = value.minComp();
1057 maxgap = value.maxComp();
1058 }
1059 else if (!vl.at(0).isNull())
1060 throw Exception("gap value {0} not recognized in method unbond in contact model {1}.", vl.at(0), getName());
1061 double gap = c->getGap();
1062 if (gap >= mingap && gap <= maxgap) {
1063 sn_state_ = 0;
1064 FP_S;
1065 return true;
1066 }
1067 FP_S;
1068 return false;
1069 }
1070 case kwArea: {
1071 if (!userArea_) {
1072 double rsq = calcRSQ(c);
1073#ifdef THREED
1074 userArea_ = rsq * rsq * dPi;
1075#else
1076 userArea_ = rsq * 2.0;
1077#endif
1078 }
1079 FP_S;
1080 return true;
1081 }
1082 case kwResetDisp:
1083 sn_ndisp_ = 0.0;
1084 for (int i=1; i<dim; ++i)
1085 sn_sdisp_.rdof(i) = 0;
1086 break;
1087 }
1088 FP_S;
1089 return false;
1090 }
1091
1092 double ContactModelRBSN::getEnergy(uint32 i) const {
1093 // Return an energy value.
1094 double ret(0.0);
1095 if (!energies_)
1096 return ret;
1097 switch (i) {
1098 case kwEStrain: return energies_->estrain_;
1099 case kwESlip: return energies_->eslip_;
1100 case kwEDashpot: return energies_->edashpot_;
1101 }
1102 assert(0);
1103 return ret;
1104 }
1105
1106 bool ContactModelRBSN::getEnergyAccumulate(uint32 i) const {
1107 // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
1108 switch (i) {
1109 case kwEStrain: return false;
1110 case kwESlip: return true;
1111 case kwEDashpot: return true;
1112 }
1113 assert(0);
1114 return false;
1115 }
1116
1117 void ContactModelRBSN::setEnergy(uint32 i,const double &d) {
1118 // Set an energy value.
1119 if (!energies_) return;
1120 switch (i) {
1121 case kwEStrain: energies_->estrain_ = d; return;
1122 case kwESlip: energies_->eslip_ = d; return;
1123 case kwEDashpot: energies_->edashpot_= d; return;
1124 }
1125 assert(0);
1126 return;
1127 }
1128
1129 bool ContactModelRBSN::validate(ContactModelMechanicalState *state,const double &) {
1130 // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
1131 // (1) the contact is created, (2) a property of the contact that returns a true via
1132 // the setProperty method has been modified and (3) when a set of cycles is executed
1133 // via the {cycle N} command.
1134 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
1135 assert(state);
1136 const IContactMechanical *c = state->getMechanicalContact();
1137 assert(c);
1138
1139 if (state->trackEnergy_)
1140 activateEnergy();
1141
1142 if (inheritanceField_ & KnMask)
1143 updateKn(c);
1144 if (inheritanceField_ & KsMask)
1145 updateKs(c);
1146 if (inheritanceField_ & FricMask)
1147 updateFric(c);
1148
1149 updateStiffness(state);
1150 return checkActivity(state->gap_);
1151 }
1152
1153 static const string knstr("kn");
1154 bool ContactModelRBSN::updateKn(const IContactMechanical *con) {
1155 assert(con);
1156 base::Property v1 = con->getEnd1()->getProperty(knstr);
1157 base::Property v2 = con->getEnd2()->getProperty(knstr);
1158 if (!v1.isValid() || !v2.isValid())
1159 return false;
1160 double kn1 = v1.toDouble();
1161 double kn2 = v2.toDouble();
1162 double val = kTran_;
1163 if (kn1 && kn2)
1164 kTran_ = kn1*kn2/(kn1+kn2);
1165 else if (kn1)
1166 kTran_ = kn1;
1167 else if (kn2)
1168 kTran_ = kn2;
1169 return ( (kTran_ != val) );
1170 }
1171
1172 static const string ksstr("ks");
1173 bool ContactModelRBSN::updateKs(const IContactMechanical *con) {
1174 assert(con);
1175 base::Property v1 = con->getEnd1()->getProperty(ksstr);
1176 base::Property v2 = con->getEnd2()->getProperty(ksstr);
1177 if (!v1.isValid() || !v2.isValid())
1178 return false;
1179 double ks1 = v1.toDouble();
1180 double ks2 = v2.toDouble();
1181 double val = kTran_;
1182 if (ks1 && ks2)
1183 kTran_ = ks1*ks2/(ks1+ks2);
1184 else if (ks1)
1185 kTran_ = ks1;
1186 else if (ks2)
1187 kTran_ = ks2;
1188 return ( (kTran_ != val) );
1189 }
1190
1191 static const string fricstr("fric");
1192 bool ContactModelRBSN::updateFric(const IContactMechanical *con) {
1193 assert(con);
1194 base::Property v1 = con->getEnd1()->getProperty(fricstr);
1195 base::Property v2 = con->getEnd2()->getProperty(fricstr);
1196 if (!v1.isValid() || !v2.isValid())
1197 return false;
1198 double fric1 = std::max(0.0,v1.toDouble());
1199 double fric2 = std::max(0.0,v2.toDouble());
1200 double val = fric_;
1201 fric_ = std::min(fric1,fric2);
1202 return ( (fric_ != val) );
1203 }
1204
1205 bool ContactModelRBSN::endPropertyUpdated(const string &name,const IContactMechanical *c) {
1206 // The endPropertyUpdated method is called whenever a surface property (with a name
1207 // that matches an inheritable contact model property name) of one of the contacting
1208 // pieces is modified. This allows the contact model to update its associated
1209 // properties. The return value denotes whether or not the update has affected
1210 // the time step computation (by having modified the translational or rotational
1211 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
1212 assert(c);
1213 StringList availableProperties = split(replace(simplified(getProperties())," ",""),",");
1214 auto idx = findRegex(availableProperties,name);
1215 if (idx==string::npos) return false;
1216 idx += 1;
1217 bool ret=false;
1218 switch(idx) {
1219 case kwKn: { //kn
1220 if (inheritanceField_ & KnMask)
1221 ret = updateKn(c);
1222 break;
1223 }
1224 case kwKs: { //ks
1225 if (inheritanceField_ & KsMask)
1226 ret =updateKs(c);
1227 break;
1228 }
1229 case kwFric: { //fric
1230 if (inheritanceField_ & FricMask)
1231 updateFric(c);
1232 break;
1233 }
1234 }
1235 return ret;
1236 }
1237
1238 ContactModelRBSN::StiffData ContactModelRBSN::computeStiffData(ContactModelMechanicalState *state) const {
1239 // Update contact data
1240 //double Cmin1 = state->end1Curvature_.x();
1241 double Cmax1 = state->end1Curvature_.y();
1242 double Cmax2 = state->end2Curvature_.y();
1243 double br = sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
1244 if (userArea_)
1245#ifdef THREED
1246 br = std::sqrt(userArea_ / dPi);
1247#else
1248 br = userArea_ / 2.0;
1249#endif
1250 StiffData ret;
1251 ret.reff_ = br;
1252 ret.trans_ = DVect2(kTran_,kTran_/kRatio_);
1253 ret.ang_ = kRot_;
1254 return ret;
1255 }
1256
1257 void ContactModelRBSN::updateStiffness(ContactModelMechanicalState *state) {
1258 // first compute stiffness data
1259 StiffData stiff = computeStiffData(state);
1260 // Now calculate effective stiffness
1261 DVect2 retT = stiff.trans_;
1262 // correction if viscous damping active
1263 if (dpProps_) {
1264 DVect2 correct(1.0);
1265 if (dpProps_->dp_nratio_)
1266 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
1267 if (dpProps_->dp_sratio_)
1268 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
1269 retT /= (correct*correct);
1270 }
1271 effectiveTranslationalStiffness_ = retT;
1272 // Effective rotational stiffness (bending and twisting)
1273 effectiveRotationalStiffness_ = stiff.ang_;
1274 }
1275
1276 bool ContactModelRBSN::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
1277 assert(state);
1278
1279 if (state->activated()) {
1280 // The contact was just activated from an inactive state
1281 // Trigger the FISH callback if one is hooked up to the
1282 // contact_activated event.
1283 if (cmEvents_[fActivated] >= 0) {
1284 auto c = state->getContact();
1285 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
1286 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1287 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
1288 }
1289 }
1290 updateStiffness(state);
1291 // accumulate shear displacement for dilation
1292 sn_ndisp_ += state->relativeTranslationalIncrement_.x();
1293 DVect shearInc = state->relativeTranslationalIncrement_;
1294 shearInc.rx() = 0;
1295 sn_sdisp_.ry() += shearInc.mag();
1296
1297 if (isBonded()) return FDLawBonded(state, timestep);
1298 else return FDLawUnBonded(state, timestep);
1299
1300 }
1301
1302 bool ContactModelRBSN::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
1303 // Account for thermal expansion in incremental mode
1304 if (sn_mode_ == 0 || ts->gapInc_ == 0.0) return false;
1305 DVect finc(0.0);
1306 finc.rx() = kTran_ * ts->gapInc_;
1307 sn_F_ -= finc;
1308 return true;
1309 }
1310
1311 bool ContactModelRBSN::FDLawBonded(ContactModelMechanicalState *state, const double ×tep) {
1312 // initialize ... get the geometry information
1313 DVect3 geom = computeGeomData(state->end1Curvature_,state->end2Curvature_);
1314 double area = geom.x();
1315 double bi = geom.y();
1316 double br = geom.z();
1317 double kn = kTran_;
1318 double ks = kn / kRatio_;
1319 double kb = kRot_.z();
1320#if DIM==3
1321 kb = std::sqrt(kb*kb + kRot_.y()*kRot_.y());
1322 double kt = kRot_.x();
1323#else
1324 double kt = 0.0;
1325#endif
1326
1327 DVect totalForce = sn_F_ + fictForce_;
1328
1329 //double nsmax0 = -(totalForce.x() / area) + sn_mcf_* sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z()) * br / bi;
1330
1331 // Relative translational/rotational displacement increments
1332 DVect trans = state->relativeTranslationalIncrement_;
1333 DAVect ang = state->relativeAngularIncrement_;
1334
1335 // Store previous force and moment
1336 DVect sn_F_old = totalForce;
1337 sn_F_old.rx() -= sn_por_ * area;
1338 DAVect sn_M_old = sn_M_;
1339
1340 DVect theStiff(ks);
1341 theStiff.rx() = kn;
1342 sn_F_ -= trans * theStiff;
1343 sn_M_ -= ang * kRot_;
1344 // Fictitious force captures Poisson's effect in RBSN models
1345 // (see common\contactmodel\rbsn\doc\manual\cmspringnetwork.rst)
1346 // Fictitious force between bonded bodies only.
1347 // No fictitious force between body and wall (zero volume contact end)
1348 if (poisson_ != 0.0 and (state->end1Volume_ and state->end2Volume_)) {
1349 auto stress11 = state->end1Stress_;
1350 auto stress22 = state->end2Stress_;
1351 double vol1 = state->end1Volume_;
1352 double vol2 = state->end2Volume_;
1353 double ms = 0.0;
1354 for (int i=0; i<stress11.size(); ++i) {
1355 stress11[i] = (stress11[i]*vol1 + stress22[i]*vol2)/(vol1 + vol2);
1356 ms = std::max(ms,abs(stress11[i]));
1357 }
1358 DMatrix<dim,dim> stresst(0.0);
1359#ifdef THREED
1360 stresst.get(0,0) = -poisson_*stress11[1] - poisson_*stress11[2];
1361 stresst.get(1,1) = -poisson_*stress11[0] - poisson_*stress11[2];
1362#else
1363 stresst.get(0,0) = -poisson_*stress11[1];
1364 stresst.get(1,1) = -poisson_*stress11[0];
1365#endif
1366#ifdef THREED
1367 stresst.get(2,2) = -poisson_*stress11[0] - poisson_*stress11[1];
1368#endif
1369 if (poisOffDiag_) {
1370#ifdef THREED
1371 double sxy = stress11[3];
1372 double szx = stress11[4];
1373 double syz = stress11[5];
1374 stresst.get(0,1) = poisson_ * sxy;
1375 stresst.get(1,0) = stresst.get(0,1);
1376 stresst.get(0,2) = poisson_ * szx;
1377 stresst.get(2,0) = stresst.get(0,2);
1378 stresst.get(1,2) = poisson_ * syz;
1379 stresst.get(2,1) = stresst.get(1,2);
1380#else
1381 double sxy = stress11[2];
1382 stresst.get(0,1) = poisson_ * sxy;
1383 stresst.get(1,0) = stresst.get(0,1);
1384#endif
1385 }
1386
1387
1388 DVect norm = state->axes_.e1();
1389 DVect traction = stresst * norm * userArea_;
1390 DVect shear(0.0);
1391 shear.ry() = 1.0;
1392 shear = state->axes_.toGlobal(shear);
1393#ifdef THREED
1394 DVect ns = state->axes_.toGlobal(DVect(0.,0.,1.));
1395 fictForce_ = DVect(norm|traction,shear|traction,ns|traction);
1396#else
1397 fictForce_ = DVect(norm|traction,shear|traction);
1398#endif
1399 if (forceSet_ && ms) {
1400 forceSet_ = false;
1401 sn_F_ -= fictForce_;
1402 }
1403 }
1404 FP_S;
1405 double porForce = sn_por_ * area;
1406 sn_F_.rx() -= porForce;
1407 totalForce = sn_F_ + fictForce_;
1408
1409 if (state->canFail_) {
1410 double dbend = sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z());
1411 double nsmax = -(totalForce.x() / area) + sn_mcf_*dbend * br / bi;
1412 bool failed = false;
1413 if (sn_state_ == 3 || sn_state_ == 5) {
1414 double compVal = sn_state_ == 3 ? tenTable_[0].x() : transTen_;
1415 if (nsmax >= compVal ) {
1416 if (tenTable_.back().y() < limits<double>::epsilon()*100)
1417 failed = true;
1418 else {
1419 if (sn_state_ == 3)
1420 sn_elong_ = 0;
1421 transTen_ = compVal;
1422 sn_state_ = 4;
1423 }
1424 }
1425 }
1426 FP_S;
1427
1428 if (sn_state_ == 4) {
1429 sn_elong_ += trans.x();
1430 sn_elong_ = std::max(0.0,sn_elong_);
1431 int ww = -1;
1432 if (sn_elong_ <= tenTable_.back().y()) {
1433 for (int i=0; i<tenTable_.size(); ++i) {
1434 if (tenTable_[i].y() >= sn_elong_) {
1435 ww = i;
1436 break;
1437 }
1438 }
1439 } else
1440 ww = tenTable_.size() - 1;
1441 if (ww > 0) {
1442 //double factor = std::min(1.0, sn_elong_ / tenTable_[ww].y());
1443 double prevVal = ww == 1 ? 1 : tenTable_[ww-1].x();
1444 double curVal = tenTable_[ww].x();
1445 double prevElong = tenTable_[ww-1].y();
1446 double curElong = tenTable_[ww].y();
1447 double slope = (curVal - prevVal)/(curElong - prevElong);
1448 FP_S;
1449 //y-y0 = m(x-x0)
1450 double nstren = slope * (sn_elong_ - prevElong) + prevVal;
1451 if (nstren <= 0)
1452 failed = true;
1453 else {
1454 nstren *= tenTable_[0].x();
1455 if (nsmax >= nstren || slope > 0) {
1456 double fac = nstren / nsmax;
1457 sn_F_.rx() *= fac;
1458#if DIM==3
1459 sn_M_.ry() *= fac;
1460#endif
1461 sn_M_.rz() *= fac;
1462 fictForce_.rx() *= fac;
1463 } else {
1464 sn_state_ = 5;
1465 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;
1466 }
1467 }
1468 }
1469 }
1470 if (sn_state_ == 6 && nsmax >= 0)
1471 failed = true;
1472 FP_S;
1473 if (failed) {
1474 // Failed in tension
1475 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1476 sn_state_ = 1;
1477 sn_F_.fill(0.0);
1478 sn_M_.fill(0.0);
1479 failed = true;
1480 fictForce_ = DVect(0.0);
1481 //sn_F_.rx() = -sn_tenres_ * area;
1482 if (cmEvents_[fBondBreak] >= 0) {
1483 auto c = state->getContact();
1484 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1485 fish::Parameter((int64)sn_state_),
1486 fish::Parameter(nsmax),
1487 fish::Parameter(se)
1488 };
1489 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1490 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1491 }
1492 }
1493 FP_S;
1494
1495 if (!failed) {
1496 /* check for shear failure */
1497 double dtwist = sn_M_.x();
1498 DVect bfs(totalForce);
1499 bfs.rx() = 0.0;
1500 double dbfs = bfs.mag();
1501 double ssmax = dbfs / area + sn_mcf_*std::abs(dtwist) * 0.5* br / bi;
1502 double ss = shearStrength(area);
1503 FP_S;
1504 if (ss < 0)
1505 ss = 0;
1506 if (ss <= ssmax) {
1507 // strength when it breaks for
1508 // Failed in shear
1509 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1510 sn_state_ = 2;
1511 fictForce_ = DVect(0.0);
1512 FP_S;
1513 sn_F_ = totalForce;
1514 if (cmEvents_[fBondBreak] >= 0) {
1515 auto c = state->getContact();
1516 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1517 fish::Parameter((int64)sn_state_),
1518 fish::Parameter(ss),
1519 fish::Parameter(se)
1520 };
1521 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1522 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1523 }
1524 double mm = 0.0;
1525 for (int i=1; i<dim; ++i)
1526 mm += trans[i]*trans[i];
1527 sn_sdisp_.rx() = std::sqrt(mm);
1528 double crit = sn_F_.x() * fric_ + sn_cohres_ * area;
1529 if (sn_cohdc() > std::numeric_limits<double>::epsilon()) {
1530 int ww = -1;
1531 if (sn_sdisp_.x() <= sn_cohdc()) {
1532 for (int i=0; i<cohTable_.size(); ++i) {
1533 if (cohTable_[i].y() >= sn_sdisp_.x()) {
1534 ww = i;
1535 break;
1536 }
1537 }
1538 } else
1539 ww = cohTable_.size() - 1;
1540 if (ww > 0) {
1541 double prevVal = ww == 1 ? 1: cohTable_[ww-1].x() ;//critBreak_ : cohTable_[ww-1].x() + sn_F_.x() * fric_;
1542 double curVal = cohTable_[ww].x();//cohTable_[ww].x() + sn_F_.x() * fric_;
1543 double prevShear = cohTable_[ww-1].y();
1544 double curShear = cohTable_[ww].y();
1545 double slope = (curVal - prevVal)/(curShear - prevShear);
1546 // should go from 0 to 1
1547 double mval = slope * (sn_sdisp_.x() - prevShear) + prevVal;
1548 double fact = std::max(0.,std::min(1.0,1.0 - mval));
1549 assert(fact >= 0 && fact <= 1.0);
1550 double theFric = fric_;
1551 if (sn_fa_)
1552 theFric = sn_fa_ + (fric_ - sn_fa_)*fact;
1553 double theCoh = sn_Coh();
1554 if (theCoh)
1555 theCoh = sn_Coh() + (sn_cohres_ - sn_Coh())*fact;
1556 crit = sn_F_.x() * theFric + theCoh * geom.x();
1557 }
1558 }
1559
1560 // Resolve sliding.
1561 FP_S;
1562 if (crit < 0)
1563 crit = 0;
1564 DVect sforce = sn_F_; sforce.rx() = 0.0;
1565 // The is the magnitude of the shear force.
1566 double sfmag = sforce.mag();
1567 // Sliding occurs when the magnitude of the shear force is greater than the
1568 // critical value.
1569 if (sfmag > crit) {
1570 // Lower the shear force to the critical value for sliding.
1571 double rat = crit / sfmag;
1572 sforce *= rat;
1573 sforce.rx() = sn_F_.x();
1574 sn_F_ = sforce;
1575 sn_S_ = true;
1576 }
1577 if (sn_S_) {
1578 if (sn_dil_ > 0) {
1579 double zdd = sn_dilzero_ != 0 ? sn_dilzero_ : limits<double>::max();
1580 double usm = sn_sdisp_.y();
1581 if (usm < zdd) {
1582 double sInc = 0.0;
1583 for (int i=1; i<dim; ++i)
1584 sInc += trans.dof(i)*trans.dof(i);
1585 sInc = std::sqrt(sInc);
1586 sn_F_.rx() += kTran_ * sn_dil_ * sInc;
1587 }
1588 }
1589 }
1590
1591 // Resolve bending
1592 crit = sn_bmul_*2.1*0.25*br*std::abs(sn_F_.x()); // Jiang 2015
1593 FP_S;
1594 DAVect test = sn_M_;
1595#if DIM==3
1596 test.rx() = 0.0;
1597#endif
1598 double tmag = test.mag();
1599 if (tmag > crit) {
1600 // Lower the bending moment to the critical value for sliding.
1601 double rat = crit / tmag;
1602 test *= rat;
1603 sn_BS_ = true;
1604 }
1605 sn_M_.rz() = test.z();
1606#if DIM==3
1607 sn_M_.ry() = test.y();
1608 // Resolve twisting
1609 crit = sn_tmul_ * 0.65*fric_* br*std::abs(sn_F_.x()) ; // Jiang 2015
1610 tmag = std::abs(sn_M_.x());
1611 if (tmag > crit) {
1612 // Lower the shear force to the critical value for sliding.
1613 double rat = crit / tmag;
1614 tmag = sn_M_.x() * rat;
1615 sn_TS_ = true;
1616 }
1617 sn_M_.rx() = tmag;
1618 FP_S;
1619#endif
1620 }
1621 }
1622 }
1623 sn_F_old.rx() += porForce;
1624 sn_F_.rx() += porForce;
1625 totalForce = sn_F_ + fictForce_;
1626 FP_S;
1627
1628 // Account for dashpot forces if the dashpot structure has been defined.
1629 if (dpProps_) {
1630 dpProps_->dp_F_.fill(0.0);
1631 double vcn(0.0), vcs(0.0);
1632 // Calculate the damping coefficients.
1633 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1634 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1635 // First damp the shear components
1636 for (int i = 1; i<dim; ++i)
1637 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1638 // Damp the normal component
1639 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1640 // Need to change behavior based on the dp_mode.
1641 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1642 // Limit in tension if not bonded.
1643 if (sn_state_ < 3 && (dpProps_->dp_F_.x() + totalForce.x() < 0))
1644 dpProps_->dp_F_.rx() = -totalForce.rx();
1645 }
1646 if (sn_S_ && dpProps_->dp_mode_ > 1) {
1647 // Limit in shear if sliding.
1648 double dfn = dpProps_->dp_F_.rx();
1649 dpProps_->dp_F_.fill(0.0);
1650 dpProps_->dp_F_.rx() = dfn;
1651 }
1652 }
1653 FP_S;
1654
1655 //Compute energies if energy tracking has been enabled.
1656 if (state->trackEnergy_) {
1657 assert(energies_);
1658 energies_->estrain_ = 0.0;
1659 if (kn)
1660 // Calculate the strain energy.
1661 energies_->estrain_ = 0.5*totalForce.x()*totalForce.x() / kn;
1662 if (ks) {
1663 DVect s = totalForce;
1664 s.rx() = 0.0;
1665 double smag2 = s.mag2();
1666 // Add the shear component of the strain energy.
1667 energies_->estrain_ += 0.5*smag2 / ks;
1668
1669 if (sn_S_) {
1670 // If sliding calculate the slip energy and accumulate it.
1671 sn_F_old.rx() = 0.0;
1672 DVect avg_F_s = (s + sn_F_old)*0.5;
1673 DVect u_s_el = (s - sn_F_old) / ks;
1674 DVect u_s(0.0);
1675 for (int i = 1; i < dim; ++i)
1676 u_s.rdof(i) = trans.dof(i);
1677 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1678 }
1679 }
1680 // Add the bending/twisting resistance energy contributions.
1681 if (kb) {
1682 DAVect tmp = sn_M_;
1683#ifdef THREED
1684 tmp.rx() = 0.0;
1685#endif
1686 energies_->estrain_ += 0.5*tmp.mag2() / kb;
1687 if (sn_BS_) {
1688 // accumulate bending slip energy.
1689 DAVect tmp_old = sn_M_old;
1690#ifdef THREED
1691 tmp_old.rx() = 0.0;
1692#endif
1693 DAVect avg_M = (tmp + tmp_old)*0.5;
1694 DAVect t_s_el = (tmp - tmp_old) / kb;
1695 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1696 }
1697 }
1698#ifdef THREED
1699 if (kt) {
1700 double mt = std::abs(sn_M_.x());
1701 energies_->estrain_ += 0.5*mt*mt / kt;
1702 if (sn_TS_) {
1703 // accumulate twisting slip energy.
1704 DAVect tmp(0.0);
1705 DAVect tmp_old(0.0);
1706 tmp.rx() = sn_M_.x();
1707 tmp_old.rx() = sn_M_old.x();
1708 DAVect avg_M = (tmp + tmp_old)*0.5;
1709 DAVect t_s_el = (tmp - tmp_old) / kt;
1710 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1711 }
1712 }
1713#endif
1714
1715 if (dpProps_) {
1716 // Calculate damping energy (accumulated) if the dashpots are active.
1717 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1718 }
1719 }
1720
1721 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1722 assert(sn_F_ == sn_F_);
1723 assert(sn_M_ == sn_M_);
1724 assert(fictForce_ == fictForce_);
1725 FP_S;
1726 return true;
1727 }
1728
1729 bool ContactModelRBSN::FDLawUnBonded(ContactModelMechanicalState *state, const double ×tep) {
1730 DVect3 geom = computeGeomData(state->end1Curvature_,state->end2Curvature_);
1731 double br = geom.z();
1732 // Relative translational/rotational displacement increments
1733 DVect trans = state->relativeTranslationalIncrement_;
1734 DAVect ang = state->relativeAngularIncrement_;
1735 double overlap = rgap_ - state->gap_;
1736 double correction = 1.0;
1737 if (state->activated() && sn_mode_ == 0 && trans.x()) {
1738 correction = -1.0*overlap / trans.x();
1739 if (correction < 0)
1740 correction = 1.0;
1741 }
1742
1743 // Store previous force and moment
1744 DVect sn_F_old = sn_F_;
1745 double porForce = sn_por_ * geom.x();
1746 sn_F_old.rx() -= porForce;
1747 DAVect sn_M_old = sn_M_;
1748
1749 double kb = kRot_.z();
1750#if DIM==3
1751 double kt = kRot_.x();
1752 //kb = std::sqrt(kb * kb + kRot_.y() * kRot_.y());
1753#endif
1754 // absolute/incremental normal force calculation
1755 DVect theStiff(kTran_/kRatio_);
1756 theStiff.rx() = kTran_;
1757
1758 if (sn_mode_==0)
1759 sn_F_.rx() = overlap * theStiff.x();
1760 else
1761 sn_F_.rx() -= trans.x() * theStiff.x();
1762 // Normal force can only be positive if unbonded
1763 sn_F_.rx() = std::max(0.0, sn_F_.x()) - porForce;
1764
1765 // Calculate the trial shear force.
1766 DVect sforce(0.0);
1767 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1768 // Loop over the shear components (note: the 0 component is the normal component)
1769 // and calculate the shear force.
1770 for (int i = 1; i<dim; ++i)
1771 sforce.rdof(i) = sn_F_.dof(i) - trans.dof(i) * theStiff.dof(i);
1772
1773 // Calculate the trial moment.
1774 DAVect mom = sn_M_ - ang*kRot_;
1775
1776 // If the SOLVE ELASTIC command is given then the
1777 // canFail state is set to FALSE. Otherwise it is always TRUE.
1778 if (state->canFail_) {
1779 bool changed = false;
1780 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
1781 bool slip_changed = false;
1782
1783 double crit = sn_F_.x() * fric_ + sn_cohres_ * geom.x();
1784 if (sn_state_ != 0) {
1785 bool doComp = true;
1786 if (!sn_S_) {
1787 if (sn_heal_) {
1788 sn_sdisp_.rx() = 0;
1789 crit = sn_F_.x() * sn_fa_ + cohTable_[0].x() * geom.x();
1790 doComp = false;
1791 }
1792 }
1793 if (doComp) {
1794 double mm = 0.0;
1795 for (int i=1; i<dim; ++i)
1796 mm += trans[i]*trans[i];
1797 sn_sdisp_.rx() += std::sqrt(mm);
1798 if (sn_cohdc() > std::numeric_limits<double>::epsilon()) {
1799 int ww = -1;
1800 if (sn_sdisp_.x() <= sn_cohdc()) {
1801 for (int i=0; i<cohTable_.size(); ++i) {
1802 if (cohTable_[i].y() >= sn_sdisp_.x()) {
1803 ww = i;
1804 break;
1805 }
1806 }
1807 } else
1808 ww = cohTable_.size() - 1;
1809 if (ww > 0) {
1810 double prevVal = ww == 1 ? 1: cohTable_[ww-1].x() ;//critBreak_ : cohTable_[ww-1].x() + sn_F_.x() * fric_;
1811 double curVal = cohTable_[ww].x();//cohTable_[ww].x() + sn_F_.x() * fric_;
1812 double prevShear = cohTable_[ww-1].y();
1813 double curShear = cohTable_[ww].y();
1814 double slope = (curVal - prevVal)/(curShear - prevShear);
1815 // should go from 0 to 1
1816 double mval = slope * (sn_sdisp_.x() - prevShear) + prevVal;
1817 double fact = std::max(0.,std::min(1.0,1.0 - mval));
1818 assert(fact >= 0 && fact <= 1.0);
1819 double theFric = fric_;
1820 if (sn_fa_)
1821 theFric = sn_fa_ + (fric_ - sn_fa_)*fact;
1822 double theCoh = sn_Coh();
1823 if (theCoh)
1824 theCoh = sn_Coh() + (sn_cohres_ - sn_Coh())*fact;
1825 crit = sn_F_.x() * theFric + theCoh * geom.x();
1826 }
1827 }
1828 }
1829 }
1830
1831 // Resolve sliding.
1832 if (crit < 0)
1833 crit = 0.0;
1834 // The is the magnitude of the shear force.
1835 double sfmag = sforce.mag();
1836 if (sfmag > crit) {
1837 // Lower the shear force to the critical value for sliding.
1838 double rat = crit / sfmag;
1839 sforce *= rat;
1840 if (!sn_S_) {
1841 slip_changed = true;
1842 changed = true;
1843 }
1844 sn_S_ = true;
1845 } else {
1846 if (sn_S_) {
1847 slip_changed = true;
1848 changed = true;
1849 }
1850 sn_S_ = false;
1851 }
1852 if (sn_S_) {
1853 if (sn_dil_ > 0) {
1854 double zdd = sn_dilzero_ != 0 ? sn_dilzero_ : limits<double>::max();
1855 double usm = sn_sdisp_.y();
1856 if (usm < zdd) {
1857 double sInc = 0.0;
1858 for (int i=1; i<dim; ++i)
1859 sInc += trans.dof(i)*trans.dof(i);
1860 sInc = std::sqrt(sInc);
1861 sn_F_.rx() += kTran_ * sn_dil_ * sInc;
1862 }
1863 }
1864 } else {
1865 if (sn_heal_) {
1866 if (shearStrength(geom.x()))
1867 sn_state_ = 6;
1868 }
1869 }
1870 // Resolve bending
1871 bool bslip_changed = false;
1872 crit = sn_bmul_ * 2.1*0.25*sn_F_.x() * br; // Jiang 2015
1873 DAVect test = mom;
1874#if DIM==3
1875 test.rx() = 0.0;
1876#endif
1877 double tmag = test.mag();
1878 if (tmag > crit) {
1879 // Lower the bending moment to the critical value for sliding.
1880 double rat = crit / tmag;
1881 test *= rat;
1882 if (!sn_BS_) {
1883 bslip_changed = true;
1884 changed = true;
1885 }
1886 sn_BS_ = true;
1887 }
1888 else {
1889 if (sn_BS_) {
1890 bslip_changed = true;
1891 changed = true;
1892 }
1893 sn_BS_ = false;
1894 }
1895 mom.rz() = test.z();
1896#if DIM==3
1897 mom.ry() = test.y();
1898 // Resolve twisting
1899 bool tslip_changed = false;
1900 crit = sn_tmul_ * 0.65*fric_*sn_F_.x() * br; // Jiang 2015
1901 tmag = std::abs(mom.x());
1902 if (tmag > crit) {
1903 // Lower the twisting moment to the critical value for sliding.
1904 double rat = crit / tmag;
1905 mom.rx() *= rat;
1906 if (!sn_TS_) {
1907 tslip_changed = true;
1908 changed = true;
1909 }
1910 sn_TS_ = true;
1911 }
1912 else {
1913 if (sn_TS_) {
1914 tslip_changed = true;
1915 changed = true;
1916 }
1917 sn_TS_ = false;
1918 }
1919#endif
1920 if (changed && cmEvents_[fSlipChange] >= 0) {
1921 int64 code = 0;
1922 if (slip_changed) {
1923 code = 1;
1924 if (bslip_changed) {
1925 code = 4;
1926#if DIM==3
1927 if (tslip_changed)
1928 code = 7;
1929#endif
1930 }
1931 }
1932 else if (bslip_changed) {
1933 code = 2;
1934#if DIM==3
1935 if (tslip_changed)
1936 code = 6;
1937#endif
1938 }
1939#if DIM==3
1940 else if (tslip_changed) {
1941 code = 3;
1942 if (slip_changed)
1943 code = 5;
1944 }
1945#endif
1946 auto c = state->getContact();
1947 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1948 fish::Parameter(code),
1949 fish::Parameter(sn_S_),
1950 fish::Parameter(sn_BS_)
1951#ifdef THREED
1952 ,fish::Parameter(sn_TS_)
1953#endif
1954 };
1955 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1956 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1957 }
1958 }
1959
1960 sn_F_.rx() += porForce;
1961 sn_F_old.rx() += porForce;
1962
1963 // Set the shear components of the total force.
1964 for (int i = 1; i<dim; ++i)
1965 sn_F_.rdof(i) = sforce.dof(i);
1966
1967 // Set the moment.
1968 sn_M_ = mom;
1969
1970 // Account for dashpot forces if the dashpot structure has been defined.
1971 if (dpProps_) {
1972 dpProps_->dp_F_.fill(0.0);
1973 double vcn(0.0), vcs(0.0);
1974 // Calculate the damping coefficients.
1975 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kTran_)));
1976 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(kTran_/kRatio_)));
1977 // First damp the shear components
1978 for (int i = 1; i<dim; ++i)
1979 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1980 // Damp the normal component
1981 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1982 // Need to change behavior based on the dp_mode.
1983 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1984 // Limit in tension if not bonded.
1985 if (dpProps_->dp_F_.x() + sn_F_.x() < 0)
1986 dpProps_->dp_F_.rx() = -sn_F_.rx();
1987 }
1988 if (sn_S_ && dpProps_->dp_mode_ > 1) {
1989 // Limit in shear if not sliding.
1990 double dfn = dpProps_->dp_F_.rx();
1991 dpProps_->dp_F_.fill(0.0);
1992 dpProps_->dp_F_.rx() = dfn;
1993 }
1994 }
1995
1996 //Compute energies if energy tracking has been enabled.
1997 if (state->trackEnergy_) {
1998 assert(energies_);
1999 energies_->estrain_ = 0.0;
2000 if (kTran_)
2001 // Calculate the strain energy.
2002 energies_->estrain_ = 0.5*sn_F_.x()*sn_F_.x() / kTran_;
2003 if (kTran_) {
2004 DVect s = sn_F_;
2005 s.rx() = 0.0;
2006 double smag2 = s.mag2();
2007 // Add the shear component of the strain energy.
2008 energies_->estrain_ += 0.5*smag2 / (kTran_/kRatio_);
2009
2010 if (sn_S_) {
2011 // If sliding calculate the slip energy and accumulate it.
2012 sn_F_old.rx() = 0.0;
2013 DVect avg_F_s = (s + sn_F_old)*0.5;
2014 DVect u_s_el = (s - sn_F_old) / (kTran_/kRatio_);
2015 DVect u_s(0.0);
2016 for (int i = 1; i < dim; ++i)
2017 u_s.rdof(i) = trans.dof(i);
2018 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
2019 }
2020 }
2021 // Add the bending/twisting resistance energy contributions.
2022 if (kb) {
2023 DAVect tmp = sn_M_;
2024#ifdef THREED
2025 tmp.rx() = 0.0;
2026#endif
2027 energies_->estrain_ += 0.5*tmp.mag2() / kb;
2028 if (sn_BS_) {
2029 // accumulate bending slip energy.
2030 DAVect tmp_old = sn_M_old;
2031#ifdef THREED
2032 tmp_old.rx() = 0.0;
2033#endif
2034 DAVect avg_M = (tmp + tmp_old)*0.5;
2035 DAVect t_s_el = (tmp - tmp_old) / kb;
2036 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
2037 }
2038 }
2039#ifdef THREED
2040 if (kt) {
2041 double mt = std::abs(sn_M_.x());
2042 energies_->estrain_ += 0.5*mt*mt / kt;
2043 if (sn_TS_) {
2044 // accumulate twisting slip energy.
2045 DAVect tmp(0.0);
2046 DAVect tmp_old(0.0);
2047 tmp.rx() = sn_M_.x();
2048 tmp_old.rx() = sn_M_old.x();
2049 DAVect avg_M = (tmp + tmp_old)*0.5;
2050 DAVect t_s_el = (tmp - tmp_old) / kt;
2051 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
2052 }
2053 }
2054#endif
2055
2056 if (dpProps_) {
2057 // Calculate damping energy (accumulated) if the dashpots are active.
2058 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
2059 }
2060 }
2061
2062 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
2063 assert(sn_F_ == sn_F_);
2064 assert(fictForce_ == fictForce_);
2065 assert(sn_M_ == sn_M_);
2066 return true;
2067 }
2068
2069 void ContactModelRBSN::setForce(const DVect &v,IContact *c) {
2070 sn_F_ = v;
2071 fictForce_ = DVect(0.0);
2072 forceSet_ = true;
2073 if (v.x() > 0)
2074 rgap_ = c->getGap() + v.x() / kTran_;
2075 }
2076
2077 void ContactModelRBSN::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
2078 // Only called for contacts with wall facets when the wall resolution scheme
2079 // is set to full!
2080 // Only do something if the contact model is of the same type
2081 if (equal(old->getContactModel()->getName(),"springnetwork") && !isBonded()) {
2082 ContactModelRBSN *oldCm = (ContactModelRBSN *)old;
2083#ifdef THREED
2084 // Need to rotate just the shear component from oldSystem to newSystem
2085
2086 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
2087 DVect axis = oldSystem.e1() & newSystem.e1();
2088 double c, ang, s;
2089 DVect re2;
2090 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
2091 axis = axis.unit();
2092 c = oldSystem.e1()|newSystem.e1();
2093 if (c > 0)
2094 c = std::min(c,1.0);
2095 else
2096 c = std::max(c,-1.0);
2097 ang = acos(c);
2098 s = sin(ang);
2099 double t = 1. - c;
2100 DMatrix<3,3> rm;
2101 rm.get(0,0) = t*axis.x()*axis.x() + c;
2102 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
2103 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
2104 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
2105 rm.get(1,1) = t*axis.y()*axis.y() + c;
2106 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
2107 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
2108 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
2109 rm.get(2,2) = t*axis.z()*axis.z() + c;
2110 re2 = rm*oldSystem.e2();
2111 }
2112 else
2113 re2 = oldSystem.e2();
2114 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
2115 axis = re2 & newSystem.e2();
2116 DVect2 tpf;
2117 DVect2 tpm;
2118 DMatrix<2,2> m;
2119 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
2120 axis = axis.unit();
2121 c = re2|newSystem.e2();
2122 if (c > 0)
2123 c = std::min(c,1.0);
2124 else
2125 c = std::max(c,-1.0);
2126 ang = acos(c);
2127 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
2128 ang *= -1;
2129 s = sin(ang);
2130 m.get(0,0) = c;
2131 m.get(1,0) = s;
2132 m.get(0,1) = -m.get(1,0);
2133 m.get(1,1) = m.get(0,0);
2134 tpf = m*DVect2(oldCm->sn_F_.y(),oldCm->sn_F_.z());
2135 tpm = m*DVect2(oldCm->sn_M_.y(),oldCm->sn_M_.z());
2136 } else {
2137 m.get(0,0) = 1.;
2138 m.get(0,1) = 0.;
2139 m.get(1,0) = 0.;
2140 m.get(1,1) = 1.;
2141 tpf = DVect2(oldCm->sn_F_.y(),oldCm->sn_F_.z());
2142 tpm = DVect2(oldCm->sn_M_.y(),oldCm->sn_M_.z());
2143 }
2144 DVect pforce = DVect(0,tpf.x(),tpf.y());
2145 //DVect pm = DVect(0,tpm.x(),tpm.y());
2146#else
2147 oldSystem;
2148 newSystem;
2149 DVect pforce = DVect(0,oldCm->sn_F_.y());
2150 //DVect pm = DVect(0,oldCm->sn_M_.y());
2151#endif
2152 for (int i=1; i<dim; ++i)
2153 sn_F_.rdof(i) += pforce.dof(i);
2154 if (sn_mode_ && oldCm->sn_mode_)
2155 sn_F_.rx() = oldCm->sn_F_.x();
2156 oldCm->sn_F_ = DVect(0.0);
2157 oldCm->sn_M_ = DAVect(0.0);
2158 if (dpProps_ && oldCm->dpProps_) {
2159#ifdef THREED
2160 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
2161 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
2162#else
2163 pforce = oldCm->dpProps_->dp_F_;
2164#endif
2165 dpProps_->dp_F_ += pforce;
2166 oldCm->dpProps_->dp_F_ = DVect(0.0);
2167 }
2168 if(oldCm->getEnergyActivated()) {
2169 activateEnergy();
2170 energies_->estrain_ = oldCm->energies_->estrain_;
2171 energies_->edashpot_ = oldCm->energies_->edashpot_;
2172 energies_->eslip_ = oldCm->energies_->eslip_;
2173 oldCm->energies_->estrain_ = 0.0;
2174 oldCm->energies_->edashpot_ = 0.0;
2175 oldCm->energies_->eslip_ = 0.0;
2176 }
2177 rgap_ = oldCm->rgap_;
2178 }
2179 assert(sn_F_ == sn_F_);
2180 }
2181
2182 void ContactModelRBSN::setNonForcePropsFrom(IContactModel *old) {
2183 // Only called for contacts with wall facets when the wall resolution scheme
2184 // is set to full!
2185 // Only do something if the contact model is of the same type
2186 if (equal(old->getName(),"springnetwork") && !isBonded()) {
2187 ContactModelRBSN *oldCm = (ContactModelRBSN *)old;
2188
2189 fictForce_ = oldCm->fictForce_;
2190 sn_F_ = oldCm->sn_F_;
2191 sn_sdisp_ = oldCm->sn_sdisp_;
2192 sn_M_ = oldCm->sn_M_;
2193 kRot_ = oldCm->kRot_;
2194 kTran_ = oldCm->kTran_;
2195 kRatio_ = oldCm->kRatio_;
2196 E_ = oldCm->E_;
2197 poisson_ = oldCm->poisson_;
2198 fric_ = oldCm->fric_;
2199 sn_bmul_ = oldCm->sn_bmul_;
2200 sn_tmul_ = oldCm->sn_tmul_;
2201 sn_rmul_ = oldCm->sn_rmul_;
2202 userArea_ = oldCm->userArea_;
2203 rgap_ = oldCm->rgap_;
2204 sn_fa_ = oldCm->sn_fa_;
2205 sn_mcf_ = oldCm->sn_mcf_;
2206 sn_dil_ = oldCm->sn_dil_;
2207 sn_dilzero_ = oldCm->sn_dilzero_;
2208 transTen_ = oldCm->transTen_;
2209 sn_elong_ = oldCm->sn_elong_;
2210 sn_ndisp_ = oldCm->sn_ndisp_;
2211 sn_mode_ = oldCm->sn_mode_;
2212 sn_state_ = oldCm->sn_state_;
2213 poisOffDiag_ = oldCm->poisOffDiag_;
2214 sn_S_ = oldCm->sn_S_;
2215 sn_BS_ = oldCm->sn_BS_;
2216 sn_TS_ = oldCm->sn_TS_;
2217 forceSet_ = oldCm->forceSet_;
2218 sn_heal_ = oldCm->sn_heal_;
2219 tenTable_ = oldCm->tenTable_;
2220 cohTable_ = oldCm->cohTable_;
2221 if (oldCm->hasDamping()) {
2222 if (!dpProps_)
2223 dpProps_ = NEW dpProps();
2224 dp_nratio(oldCm->dp_nratio());
2225 dp_sratio(oldCm->dp_sratio());
2226 dp_mode(oldCm->dp_mode());
2227 dp_F(oldCm->dp_F());
2228 }
2229 }
2230 }
2231
2232 DVect ContactModelRBSN::getForce() const {
2233 DVect ret(sn_F_);
2234 ret += fictForce_;
2235 if (dpProps_)
2236 ret += dpProps_->dp_F_;
2237 return ret;
2238 }
2239
2240 DAVect ContactModelRBSN::getMomentOn1(const IContactMechanical *c) const {
2241 DAVect ret(sn_M_);
2242 if (c) {
2243 DVect force = getForce();
2244 c->updateResultingTorqueOn1Local(force,&ret);
2245 }
2246 return ret;
2247 }
2248
2249 DAVect ContactModelRBSN::getMomentOn2(const IContactMechanical *c) const {
2250 DAVect ret(sn_M_);
2251 if (c) {
2252 DVect force = getForce();
2253 c->updateResultingTorqueOn2Local(force,&ret);
2254 }
2255 return ret;
2256 }
2257
2258 DAVect ContactModelRBSN::getModelMomentOn1() const {
2259 DAVect ret(sn_M_);
2260 return ret;
2261 }
2262
2263 DAVect ContactModelRBSN::getModelMomentOn2() const {
2264 DAVect ret(sn_M_);
2265 return ret;
2266 }
2267
2268 void ContactModelRBSN::objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *ret) const {
2269 ret->clear();
2270 ret->push_back({"kn",ScalarInfo});
2271 ret->push_back({"ks",ScalarInfo});
2272 ret->push_back({"krot",VectorInfo});
2273 ret->push_back({"fric",ScalarInfo});
2274 ret->push_back({"sn_bmul",ScalarInfo});
2275 ret->push_back({"sn_tmul",ScalarInfo});
2276 ret->push_back({"sn_mode",ScalarInfo});
2277 ret->push_back({"sn_force",VectorInfo});
2278 ret->push_back({"sn_moment",VectorInfo});
2279 ret->push_back({"sn_slip",ScalarInfo});
2280 ret->push_back({"sn_slipb",ScalarInfo});
2281 ret->push_back({"sn_slipt",ScalarInfo});
2282 ret->push_back({"sn_rmul",ScalarInfo});
2283 ret->push_back({"sn_radius",ScalarInfo});
2284 ret->push_back({"emod",ScalarInfo});
2285 ret->push_back({"kratio",ScalarInfo});
2286 ret->push_back({"dp_nratio",ScalarInfo});
2287 ret->push_back({"dp_sratio",ScalarInfo});
2288 ret->push_back({"dp_mode",ScalarInfo});
2289 ret->push_back({"dp_force",VectorInfo});
2290 ret->push_back({"sn_state",ScalarInfo});
2291 ret->push_back({"sn_ten",ScalarInfo});
2292 ret->push_back({"sn_shear",ScalarInfo});
2293 ret->push_back({"sn_coh",ScalarInfo});
2294 ret->push_back({"sn_fa",ScalarInfo});
2295 ret->push_back({"sn_mcf",ScalarInfo});
2296 ret->push_back({"sn_sigma",ScalarInfo});
2297 ret->push_back({"sn_tau",ScalarInfo});
2298 ret->push_back({"sn_area",ScalarInfo});
2299 ret->push_back({"user_area",ScalarInfo});
2300 ret->push_back({"rgap",ScalarInfo});
2301 ret->push_back({"sn_pois_force",VectorInfo});
2302 ret->push_back({"sn_pois",ScalarInfo});
2303 ret->push_back({"sn_non_diag",ScalarInfo});
2304 ret->push_back({"sn_cohres",ScalarInfo});
2305 ret->push_back({"sn_dil",ScalarInfo});
2306 ret->push_back({"sn_dilzero",ScalarInfo});
2307 ret->push_back({"sn_ndisp",ScalarInfo});
2308 ret->push_back({"sn_sdisp",VectorInfo});
2309 ret->push_back({"sn_cohdc",ScalarInfo});
2310 ret->push_back({"sn_heal",ScalarInfo});
2311 ret->push_back({"sn_tendc",ScalarInfo});
2312 ret->push_back({"sn_tabpos",ScalarInfo});
2313 ret->push_back({"sn_porp",ScalarInfo});
2314 ret->push_back({"sn_esigma",ScalarInfo});
2315
2316 }
2317
2318 void ContactModelRBSN::objectPropValues(std::vector<double> *ret,const IContact *con) const {
2319 FP_S;
2320 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
2321 ret->clear();
2322 ret->push_back(kTran_);
2323 ret->push_back(kTran_/kRatio_);
2324 ret->push_back(kRot_.mag());
2325 ret->push_back(fric_);
2326 ret->push_back(sn_bmul_);
2327 ret->push_back(sn_tmul_);
2328 ret->push_back(sn_mode_);
2329 ret->push_back(sn_F_.mag());
2330 ret->push_back(sn_M_.mag());
2331 ret->push_back(sn_S_);
2332 ret->push_back(sn_BS_);
2333 ret->push_back(sn_TS_);
2334 ret->push_back(sn_rmul_);
2335 double Cmax1 = c->getEnd1Curvature().y();
2336 double Cmax2 = c->getEnd2Curvature().y();
2337 double d;
2338 if (!userArea_)
2339 d= sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
2340 else {
2341#ifdef THREED
2342 d= std::sqrt(userArea_ / dPi);
2343#else
2344 d= userArea_ / 2.0;
2345#endif
2346 }
2347 ret->push_back(d);
2348 ret->push_back(E_);
2349 ret->push_back(kRatio_);
2350 ret->push_back(dpProps_ ? dpProps_->dp_nratio_ : 0);
2351 ret->push_back(dpProps_ ? dpProps_->dp_sratio_ : 0);
2352 ret->push_back(dpProps_ ? dpProps_->dp_mode_ : 0);
2353 ret->push_back(dpProps_ ? dpProps_->dp_F_.mag() : 0);
2354 ret->push_back(sn_state_);
2355 ret->push_back(sn_Ten());
2356 ret->push_back(shearStrength(computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x()));
2357 ret->push_back(cohTable_[0].x());
2358 ret->push_back(std::atan(sn_fa_)/dDegrad);
2359 ret->push_back(sn_mcf_);
2360 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x());
2361 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).y());
2362 ret->push_back(userArea_ ? userArea_ : computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x());
2363 ret->push_back(userArea_);
2364 ret->push_back(rgap_);
2365 ret->push_back(fictForce_.mag());
2366 ret->push_back(poisson_);
2367 ret->push_back(poisOffDiag_ == false ? 0 : 1);
2368 ret->push_back(sn_cohres_);
2369 ret->push_back(std::atan(sn_dil_)/dDegrad);
2370 ret->push_back(sn_dilzero_);
2371 ret->push_back(sn_ndisp_);
2372 ret->push_back(sn_sdisp_.mag());
2373 ret->push_back(sn_cohdc());
2374 ret->push_back(sn_heal_);
2375 ret->push_back(sn_tendc());
2376 ret->push_back(sn_tabPos_+1);
2377 ret->push_back(sn_por_);
2378 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x() + sn_por_);
2379 FP_S;
2380 }
2381
2382 DVect3 ContactModelRBSN::computeGeomData(const DVect2 &end1Curvature,const DVect2 &end2Curvature) const {
2383 double Cmax1 = end1Curvature.y();
2384 double Cmax2 = end2Curvature.y();
2385 double br = sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
2386 if (userArea_)
2387#ifdef THREED
2388 br = std::sqrt(userArea_ / dPi);
2389#else
2390 br = userArea_ / 2.0;
2391#endif
2392 double br2 = br * br;
2393#ifdef THREED
2394 double area = dPi * br2;
2395 double bi = 0.25*area*br2;
2396#else
2397 double area = 2.0*br;
2398 double bi = 2.0*br*br2 / 3.0;
2399#endif
2400 return DVect3(area, bi, br);
2401 }
2402
2403 DVect2 ContactModelRBSN::SMax(const DVect2 &e1c,const DVect2 &e2c) const {
2404 DVect3 data = computeGeomData(e1c,e2c);
2405 double area = data.x();
2406 double bi = data.y();
2407 double br = data.z();
2408 /* maximum stresses */
2409 //double nsmax0 = -(totalForce.x() / area) + sn_mcf_* sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z()) * br / bi;
2410 double dbend = sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z());
2411 dbend *= sn_mcf_;
2412 double dtwist = sn_M_.x();
2413 DVect bfs(sn_F_);
2414 bfs.rx() = 0.0;
2415 double dbfs = bfs.mag();
2416 double nsmax = -((sn_F_.x()+fictForce_.x()) / area) + dbend * br / bi;
2417 double ssmax = dbfs / area + std::abs(dtwist) * 0.5* br / bi;
2418 return DVect2(nsmax, ssmax);
2419 }
2420
2421 double ContactModelRBSN::shearStrength(const double &area) const {
2422 double sig = -1.0*(sn_F_.x() + fictForce_.x()) / area;
2423 double nstr = (sn_state_ > 2 && sn_state_ != 6) ? tenTable_[0].x() : 0.0;
2424 return sig <= nstr ? cohTable_[0].x() - sn_fa_*sig
2425 : cohTable_[0].x() - sn_fa_*nstr;
2426 }
2427
2428
2429 double ContactModelRBSN::strainEnergy(double kn,double ks,double kb,double kt) const {
2430 double ret(0.0);
2431 if (kn)
2432 ret = 0.5 * (sn_F_.x()+fictForce_.x()) * (sn_F_.x()+fictForce_.x()) / kn;
2433 if (ks) {
2434 DVect tmp = sn_F_ + fictForce_;
2435 tmp.rx() = 0.0;
2436 double smag2 = tmp.mag2();
2437 ret += 0.5 * smag2 / ks;
2438 }
2439
2440 if (kt)
2441 ret += 0.5 * sn_M_.x() * sn_M_.x() / kt;
2442 if (kb) {
2443 DAVect tmp = sn_M_;
2444#ifdef THREED
2445 tmp.rx() = 0.0;
2446 double smag2 = tmp.mag2();
2447#else
2448 double smag2 = tmp.z() * tmp.z();
2449#endif
2450 ret += 0.5 * smag2 / kb;
2451 }
2452 return ret;
2453 }
2454
2455} // namespace cmodelsxd
2456// EoF
| Was this helpful? ... | Itasca Software © 2026, Itasca | Updated: May 13, 2026 |
