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 if (poisson_ != 0.0 and (state->end1Volume_ or state->end2Volume_)) {
1345 //const IContact *con = state->getContact();
1346 //const IBody *b1 = con->getEnd1()->getIBody();
1347 //const IBody *b2 = con->getEnd2()->getIBody();
1348 auto stress11 = state->end1Stress_;
1349 auto stress22 = state->end2Stress_;
1350 double vol1 = state->end1Volume_;
1351 double vol2 = state->end2Volume_;
1352 double ms = 0.0;
1353 for (int i=0; i<stress11.size(); ++i) {
1354 stress11[i] = (stress11[i]*vol1 + stress22[i]*vol2)/(vol1 + vol2);
1355 ms = std::max(ms,abs(stress11[i]));
1356 }
1357 DMatrix<dim,dim> stresst(0.0);
1358#ifdef THREED
1359 stresst.get(0,0) = -poisson_*stress11[1] - poisson_*stress11[2];
1360 stresst.get(1,1) = -poisson_*stress11[0] - poisson_*stress11[2];
1361#else
1362 stresst.get(0,0) = -poisson_*stress11[1];
1363 stresst.get(1,1) = -poisson_*stress11[0];
1364#endif
1365#ifdef THREED
1366 stresst.get(2,2) = -poisson_*stress11[0] - poisson_*stress11[1];
1367#endif
1368 if (poisOffDiag_) {
1369#ifdef THREED
1370 double sxy = stress11[3];
1371 double szx = stress11[4];
1372 double syz = stress11[5];
1373 stresst.get(0,1) = poisson_ * sxy;
1374 stresst.get(1,0) = stresst.get(0,1);
1375 stresst.get(0,2) = poisson_ * szx;
1376 stresst.get(2,0) = stresst.get(0,2);
1377 stresst.get(1,2) = poisson_ * syz;
1378 stresst.get(2,1) = stresst.get(1,2);
1379#else
1380 double sxy = stress11[2];
1381 stresst.get(0,1) = poisson_ * sxy;
1382 stresst.get(1,0) = stresst.get(0,1);
1383#endif
1384 }
1385
1386
1387 DVect norm = state->axes_.e1();
1388 DVect traction = stresst * norm * userArea_;
1389 DVect shear(0.0);
1390 shear.ry() = 1.0;
1391 shear = state->axes_.toGlobal(shear);
1392#ifdef THREED
1393 DVect ns = state->axes_.toGlobal(DVect(0.,0.,1.));
1394 fictForce_ = DVect(norm|traction,shear|traction,ns|traction);
1395#else
1396 fictForce_ = DVect(norm|traction,shear|traction);
1397#endif
1398 if (forceSet_ && ms) {
1399 forceSet_ = false;
1400 sn_F_ -= fictForce_;
1401 }
1402 }
1403 FP_S;
1404 double porForce = sn_por_ * area;
1405 sn_F_.rx() -= porForce;
1406 totalForce = sn_F_ + fictForce_;
1407
1408 if (state->canFail_) {
1409 double dbend = sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z());
1410 double nsmax = -(totalForce.x() / area) + sn_mcf_*dbend * br / bi;
1411 bool failed = false;
1412 if (sn_state_ == 3 || sn_state_ == 5) {
1413 double compVal = sn_state_ == 3 ? tenTable_[0].x() : transTen_;
1414 if (nsmax >= compVal ) {
1415 if (tenTable_.back().y() < limits<double>::epsilon()*100)
1416 failed = true;
1417 else {
1418 if (sn_state_ == 3)
1419 sn_elong_ = 0;
1420 transTen_ = compVal;
1421 sn_state_ = 4;
1422 }
1423 }
1424 }
1425 FP_S;
1426
1427 if (sn_state_ == 4) {
1428 sn_elong_ += trans.x();
1429 sn_elong_ = std::max(0.0,sn_elong_);
1430 int ww = -1;
1431 if (sn_elong_ <= tenTable_.back().y()) {
1432 for (int i=0; i<tenTable_.size(); ++i) {
1433 if (tenTable_[i].y() >= sn_elong_) {
1434 ww = i;
1435 break;
1436 }
1437 }
1438 } else
1439 ww = tenTable_.size() - 1;
1440 if (ww > 0) {
1441 //double factor = std::min(1.0, sn_elong_ / tenTable_[ww].y());
1442 double prevVal = ww == 1 ? 1 : tenTable_[ww-1].x();
1443 double curVal = tenTable_[ww].x();
1444 double prevElong = tenTable_[ww-1].y();
1445 double curElong = tenTable_[ww].y();
1446 double slope = (curVal - prevVal)/(curElong - prevElong);
1447 FP_S;
1448 //y-y0 = m(x-x0)
1449 double nstren = slope * (sn_elong_ - prevElong) + prevVal;
1450 if (nstren <= 0)
1451 failed = true;
1452 else {
1453 nstren *= tenTable_[0].x();
1454 if (nsmax >= nstren || slope > 0) {
1455 double fac = nstren / nsmax;
1456 sn_F_.rx() *= fac;
1457#if DIM==3
1458 sn_M_.ry() *= fac;
1459#endif
1460 sn_M_.rz() *= fac;
1461 fictForce_.rx() *= fac;
1462 } else {
1463 sn_state_ = 5;
1464 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;
1465 }
1466 }
1467 }
1468 }
1469 if (sn_state_ == 6 && nsmax >= 0)
1470 failed = true;
1471 FP_S;
1472 if (failed) {
1473 // Failed in tension
1474 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1475 sn_state_ = 1;
1476 sn_F_.fill(0.0);
1477 sn_M_.fill(0.0);
1478 failed = true;
1479 fictForce_ = DVect(0.0);
1480 //sn_F_.rx() = -sn_tenres_ * area;
1481 if (cmEvents_[fBondBreak] >= 0) {
1482 auto c = state->getContact();
1483 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1484 fish::Parameter((int64)sn_state_),
1485 fish::Parameter(nsmax),
1486 fish::Parameter(se)
1487 };
1488 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1489 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1490 }
1491 }
1492 FP_S;
1493
1494 if (!failed) {
1495 /* check for shear failure */
1496 double dtwist = sn_M_.x();
1497 DVect bfs(totalForce);
1498 bfs.rx() = 0.0;
1499 double dbfs = bfs.mag();
1500 double ssmax = dbfs / area + sn_mcf_*std::abs(dtwist) * 0.5* br / bi;
1501 double ss = shearStrength(area);
1502 FP_S;
1503 if (ss < 0)
1504 ss = 0;
1505 if (ss <= ssmax) {
1506 // strength when it breaks for
1507 // Failed in shear
1508 double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1509 sn_state_ = 2;
1510 fictForce_ = DVect(0.0);
1511 FP_S;
1512 sn_F_ = totalForce;
1513 if (cmEvents_[fBondBreak] >= 0) {
1514 auto c = state->getContact();
1515 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1516 fish::Parameter((int64)sn_state_),
1517 fish::Parameter(ss),
1518 fish::Parameter(se)
1519 };
1520 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1521 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1522 }
1523 double mm = 0.0;
1524 for (int i=1; i<dim; ++i)
1525 mm += trans[i]*trans[i];
1526 sn_sdisp_.rx() = std::sqrt(mm);
1527 double crit = sn_F_.x() * fric_ + sn_cohres_ * area;
1528 if (sn_cohdc() > std::numeric_limits<double>::epsilon()) {
1529 int ww = -1;
1530 if (sn_sdisp_.x() <= sn_cohdc()) {
1531 for (int i=0; i<cohTable_.size(); ++i) {
1532 if (cohTable_[i].y() >= sn_sdisp_.x()) {
1533 ww = i;
1534 break;
1535 }
1536 }
1537 } else
1538 ww = cohTable_.size() - 1;
1539 if (ww > 0) {
1540 double prevVal = ww == 1 ? 1: cohTable_[ww-1].x() ;//critBreak_ : cohTable_[ww-1].x() + sn_F_.x() * fric_;
1541 double curVal = cohTable_[ww].x();//cohTable_[ww].x() + sn_F_.x() * fric_;
1542 double prevShear = cohTable_[ww-1].y();
1543 double curShear = cohTable_[ww].y();
1544 double slope = (curVal - prevVal)/(curShear - prevShear);
1545 // should go from 0 to 1
1546 double mval = slope * (sn_sdisp_.x() - prevShear) + prevVal;
1547 double fact = std::max(0.,std::min(1.0,1.0 - mval));
1548 assert(fact >= 0 && fact <= 1.0);
1549 double theFric = fric_;
1550 if (sn_fa_)
1551 theFric = sn_fa_ + (fric_ - sn_fa_)*fact;
1552 double theCoh = sn_Coh();
1553 if (theCoh)
1554 theCoh = sn_Coh() + (sn_cohres_ - sn_Coh())*fact;
1555 crit = sn_F_.x() * theFric + theCoh * geom.x();
1556 }
1557 }
1558
1559 // Resolve sliding.
1560 FP_S;
1561 if (crit < 0)
1562 crit = 0;
1563 DVect sforce = sn_F_; sforce.rx() = 0.0;
1564 // The is the magnitude of the shear force.
1565 double sfmag = sforce.mag();
1566 // Sliding occurs when the magnitude of the shear force is greater than the
1567 // critical value.
1568 if (sfmag > crit) {
1569 // Lower the shear force to the critical value for sliding.
1570 double rat = crit / sfmag;
1571 sforce *= rat;
1572 sforce.rx() = sn_F_.x();
1573 sn_F_ = sforce;
1574 sn_S_ = true;
1575 }
1576 if (sn_S_) {
1577 if (sn_dil_ > 0) {
1578 double zdd = sn_dilzero_ != 0 ? sn_dilzero_ : limits<double>::max();
1579 double usm = sn_sdisp_.y();
1580 if (usm < zdd) {
1581 double sInc = 0.0;
1582 for (int i=1; i<dim; ++i)
1583 sInc += trans.dof(i)*trans.dof(i);
1584 sInc = std::sqrt(sInc);
1585 sn_F_.rx() += kTran_ * sn_dil_ * sInc;
1586 }
1587 }
1588 }
1589
1590 // Resolve bending
1591 crit = sn_bmul_*2.1*0.25*br*std::abs(sn_F_.x()); // Jiang 2015
1592 FP_S;
1593 DAVect test = sn_M_;
1594#if DIM==3
1595 test.rx() = 0.0;
1596#endif
1597 double tmag = test.mag();
1598 if (tmag > crit) {
1599 // Lower the bending moment to the critical value for sliding.
1600 double rat = crit / tmag;
1601 test *= rat;
1602 sn_BS_ = true;
1603 }
1604 sn_M_.rz() = test.z();
1605#if DIM==3
1606 sn_M_.ry() = test.y();
1607 // Resolve twisting
1608 crit = sn_tmul_ * 0.65*fric_* br*std::abs(sn_F_.x()) ; // Jiang 2015
1609 tmag = std::abs(sn_M_.x());
1610 if (tmag > crit) {
1611 // Lower the shear force to the critical value for sliding.
1612 double rat = crit / tmag;
1613 tmag = sn_M_.x() * rat;
1614 sn_TS_ = true;
1615 }
1616 sn_M_.rx() = tmag;
1617 FP_S;
1618#endif
1619 }
1620 }
1621 }
1622 sn_F_old.rx() += porForce;
1623 sn_F_.rx() += porForce;
1624 totalForce = sn_F_ + fictForce_;
1625 FP_S;
1626
1627 // Account for dashpot forces if the dashpot structure has been defined.
1628 if (dpProps_) {
1629 dpProps_->dp_F_.fill(0.0);
1630 double vcn(0.0), vcs(0.0);
1631 // Calculate the damping coefficients.
1632 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1633 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1634 // First damp the shear components
1635 for (int i = 1; i<dim; ++i)
1636 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1637 // Damp the normal component
1638 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1639 // Need to change behavior based on the dp_mode.
1640 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1641 // Limit in tension if not bonded.
1642 if (sn_state_ < 3 && (dpProps_->dp_F_.x() + totalForce.x() < 0))
1643 dpProps_->dp_F_.rx() = -totalForce.rx();
1644 }
1645 if (sn_S_ && dpProps_->dp_mode_ > 1) {
1646 // Limit in shear if sliding.
1647 double dfn = dpProps_->dp_F_.rx();
1648 dpProps_->dp_F_.fill(0.0);
1649 dpProps_->dp_F_.rx() = dfn;
1650 }
1651 }
1652 FP_S;
1653
1654 //Compute energies if energy tracking has been enabled.
1655 if (state->trackEnergy_) {
1656 assert(energies_);
1657 energies_->estrain_ = 0.0;
1658 if (kn)
1659 // Calculate the strain energy.
1660 energies_->estrain_ = 0.5*totalForce.x()*totalForce.x() / kn;
1661 if (ks) {
1662 DVect s = totalForce;
1663 s.rx() = 0.0;
1664 double smag2 = s.mag2();
1665 // Add the shear component of the strain energy.
1666 energies_->estrain_ += 0.5*smag2 / ks;
1667
1668 if (sn_S_) {
1669 // If sliding calculate the slip energy and accumulate it.
1670 sn_F_old.rx() = 0.0;
1671 DVect avg_F_s = (s + sn_F_old)*0.5;
1672 DVect u_s_el = (s - sn_F_old) / ks;
1673 DVect u_s(0.0);
1674 for (int i = 1; i < dim; ++i)
1675 u_s.rdof(i) = trans.dof(i);
1676 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1677 }
1678 }
1679 // Add the bending/twisting resistance energy contributions.
1680 if (kb) {
1681 DAVect tmp = sn_M_;
1682#ifdef THREED
1683 tmp.rx() = 0.0;
1684#endif
1685 energies_->estrain_ += 0.5*tmp.mag2() / kb;
1686 if (sn_BS_) {
1687 // accumulate bending slip energy.
1688 DAVect tmp_old = sn_M_old;
1689#ifdef THREED
1690 tmp_old.rx() = 0.0;
1691#endif
1692 DAVect avg_M = (tmp + tmp_old)*0.5;
1693 DAVect t_s_el = (tmp - tmp_old) / kb;
1694 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1695 }
1696 }
1697#ifdef THREED
1698 if (kt) {
1699 double mt = std::abs(sn_M_.x());
1700 energies_->estrain_ += 0.5*mt*mt / kt;
1701 if (sn_TS_) {
1702 // accumulate twisting slip energy.
1703 DAVect tmp(0.0);
1704 DAVect tmp_old(0.0);
1705 tmp.rx() = sn_M_.x();
1706 tmp_old.rx() = sn_M_old.x();
1707 DAVect avg_M = (tmp + tmp_old)*0.5;
1708 DAVect t_s_el = (tmp - tmp_old) / kt;
1709 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1710 }
1711 }
1712#endif
1713
1714 if (dpProps_) {
1715 // Calculate damping energy (accumulated) if the dashpots are active.
1716 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1717 }
1718 }
1719
1720 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1721 assert(sn_F_ == sn_F_);
1722 assert(sn_M_ == sn_M_);
1723 assert(fictForce_ == fictForce_);
1724 FP_S;
1725 return true;
1726 }
1727
1728 bool ContactModelRBSN::FDLawUnBonded(ContactModelMechanicalState *state, const double ×tep) {
1729 DVect3 geom = computeGeomData(state->end1Curvature_,state->end2Curvature_);
1730 double br = geom.z();
1731 // Relative translational/rotational displacement increments
1732 DVect trans = state->relativeTranslationalIncrement_;
1733 DAVect ang = state->relativeAngularIncrement_;
1734 double overlap = rgap_ - state->gap_;
1735 double correction = 1.0;
1736 if (state->activated() && sn_mode_ == 0 && trans.x()) {
1737 correction = -1.0*overlap / trans.x();
1738 if (correction < 0)
1739 correction = 1.0;
1740 }
1741
1742 // Store previous force and moment
1743 DVect sn_F_old = sn_F_;
1744 double porForce = sn_por_ * geom.x();
1745 sn_F_old.rx() -= porForce;
1746 DAVect sn_M_old = sn_M_;
1747
1748 double kb = kRot_.z();
1749#if DIM==3
1750 double kt = kRot_.x();
1751 //kb = std::sqrt(kb * kb + kRot_.y() * kRot_.y());
1752#endif
1753 // absolute/incremental normal force calculation
1754 DVect theStiff(kTran_/kRatio_);
1755 theStiff.rx() = kTran_;
1756
1757 if (sn_mode_==0)
1758 sn_F_.rx() = overlap * theStiff.x();
1759 else
1760 sn_F_.rx() -= trans.x() * theStiff.x();
1761 // Normal force can only be positive if unbonded
1762 sn_F_.rx() = std::max(0.0, sn_F_.x()) - porForce;
1763
1764 // Calculate the trial shear force.
1765 DVect sforce(0.0);
1766 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1767 // Loop over the shear components (note: the 0 component is the normal component)
1768 // and calculate the shear force.
1769 for (int i = 1; i<dim; ++i)
1770 sforce.rdof(i) = sn_F_.dof(i) - trans.dof(i) * theStiff.dof(i);
1771
1772 // Calculate the trial moment.
1773 DAVect mom = sn_M_ - ang*kRot_;
1774
1775 // If the SOLVE ELASTIC command is given then the
1776 // canFail state is set to FALSE. Otherwise it is always TRUE.
1777 if (state->canFail_) {
1778 bool changed = false;
1779 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
1780 bool slip_changed = false;
1781
1782 double crit = sn_F_.x() * fric_ + sn_cohres_ * geom.x();
1783 if (sn_state_ != 0) {
1784 bool doComp = true;
1785 if (!sn_S_) {
1786 if (sn_heal_) {
1787 sn_sdisp_.rx() = 0;
1788 crit = sn_F_.x() * sn_fa_ + cohTable_[0].x() * geom.x();
1789 doComp = false;
1790 }
1791 }
1792 if (doComp) {
1793 double mm = 0.0;
1794 for (int i=1; i<dim; ++i)
1795 mm += trans[i]*trans[i];
1796 sn_sdisp_.rx() += std::sqrt(mm);
1797 if (sn_cohdc() > std::numeric_limits<double>::epsilon()) {
1798 int ww = -1;
1799 if (sn_sdisp_.x() <= sn_cohdc()) {
1800 for (int i=0; i<cohTable_.size(); ++i) {
1801 if (cohTable_[i].y() >= sn_sdisp_.x()) {
1802 ww = i;
1803 break;
1804 }
1805 }
1806 } else
1807 ww = cohTable_.size() - 1;
1808 if (ww > 0) {
1809 double prevVal = ww == 1 ? 1: cohTable_[ww-1].x() ;//critBreak_ : cohTable_[ww-1].x() + sn_F_.x() * fric_;
1810 double curVal = cohTable_[ww].x();//cohTable_[ww].x() + sn_F_.x() * fric_;
1811 double prevShear = cohTable_[ww-1].y();
1812 double curShear = cohTable_[ww].y();
1813 double slope = (curVal - prevVal)/(curShear - prevShear);
1814 // should go from 0 to 1
1815 double mval = slope * (sn_sdisp_.x() - prevShear) + prevVal;
1816 double fact = std::max(0.,std::min(1.0,1.0 - mval));
1817 assert(fact >= 0 && fact <= 1.0);
1818 double theFric = fric_;
1819 if (sn_fa_)
1820 theFric = sn_fa_ + (fric_ - sn_fa_)*fact;
1821 double theCoh = sn_Coh();
1822 if (theCoh)
1823 theCoh = sn_Coh() + (sn_cohres_ - sn_Coh())*fact;
1824 crit = sn_F_.x() * theFric + theCoh * geom.x();
1825 }
1826 }
1827 }
1828 }
1829
1830 // Resolve sliding.
1831 if (crit < 0)
1832 crit = 0.0;
1833 // The is the magnitude of the shear force.
1834 double sfmag = sforce.mag();
1835 if (sfmag > crit) {
1836 // Lower the shear force to the critical value for sliding.
1837 double rat = crit / sfmag;
1838 sforce *= rat;
1839 if (!sn_S_) {
1840 slip_changed = true;
1841 changed = true;
1842 }
1843 sn_S_ = true;
1844 } else {
1845 if (sn_S_) {
1846 slip_changed = true;
1847 changed = true;
1848 }
1849 sn_S_ = false;
1850 }
1851 if (sn_S_) {
1852 if (sn_dil_ > 0) {
1853 double zdd = sn_dilzero_ != 0 ? sn_dilzero_ : limits<double>::max();
1854 double usm = sn_sdisp_.y();
1855 if (usm < zdd) {
1856 double sInc = 0.0;
1857 for (int i=1; i<dim; ++i)
1858 sInc += trans.dof(i)*trans.dof(i);
1859 sInc = std::sqrt(sInc);
1860 sn_F_.rx() += kTran_ * sn_dil_ * sInc;
1861 }
1862 }
1863 } else {
1864 if (sn_heal_) {
1865 if (shearStrength(geom.x()))
1866 sn_state_ = 6;
1867 }
1868 }
1869 // Resolve bending
1870 bool bslip_changed = false;
1871 crit = sn_bmul_ * 2.1*0.25*sn_F_.x() * br; // Jiang 2015
1872 DAVect test = mom;
1873#if DIM==3
1874 test.rx() = 0.0;
1875#endif
1876 double tmag = test.mag();
1877 if (tmag > crit) {
1878 // Lower the bending moment to the critical value for sliding.
1879 double rat = crit / tmag;
1880 test *= rat;
1881 if (!sn_BS_) {
1882 bslip_changed = true;
1883 changed = true;
1884 }
1885 sn_BS_ = true;
1886 }
1887 else {
1888 if (sn_BS_) {
1889 bslip_changed = true;
1890 changed = true;
1891 }
1892 sn_BS_ = false;
1893 }
1894 mom.rz() = test.z();
1895#if DIM==3
1896 mom.ry() = test.y();
1897 // Resolve twisting
1898 bool tslip_changed = false;
1899 crit = sn_tmul_ * 0.65*fric_*sn_F_.x() * br; // Jiang 2015
1900 tmag = std::abs(mom.x());
1901 if (tmag > crit) {
1902 // Lower the twisting moment to the critical value for sliding.
1903 double rat = crit / tmag;
1904 mom.rx() *= rat;
1905 if (!sn_TS_) {
1906 tslip_changed = true;
1907 changed = true;
1908 }
1909 sn_TS_ = true;
1910 }
1911 else {
1912 if (sn_TS_) {
1913 tslip_changed = true;
1914 changed = true;
1915 }
1916 sn_TS_ = false;
1917 }
1918#endif
1919 if (changed && cmEvents_[fSlipChange] >= 0) {
1920 int64 code = 0;
1921 if (slip_changed) {
1922 code = 1;
1923 if (bslip_changed) {
1924 code = 4;
1925#if DIM==3
1926 if (tslip_changed)
1927 code = 7;
1928#endif
1929 }
1930 }
1931 else if (bslip_changed) {
1932 code = 2;
1933#if DIM==3
1934 if (tslip_changed)
1935 code = 6;
1936#endif
1937 }
1938#if DIM==3
1939 else if (tslip_changed) {
1940 code = 3;
1941 if (slip_changed)
1942 code = 5;
1943 }
1944#endif
1945 auto c = state->getContact();
1946 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1947 fish::Parameter(code),
1948 fish::Parameter(sn_S_),
1949 fish::Parameter(sn_BS_)
1950#ifdef THREED
1951 ,fish::Parameter(sn_TS_)
1952#endif
1953 };
1954 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1955 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1956 }
1957 }
1958
1959 sn_F_.rx() += porForce;
1960 sn_F_old.rx() += porForce;
1961
1962 // Set the shear components of the total force.
1963 for (int i = 1; i<dim; ++i)
1964 sn_F_.rdof(i) = sforce.dof(i);
1965
1966 // Set the moment.
1967 sn_M_ = mom;
1968
1969 // Account for dashpot forces if the dashpot structure has been defined.
1970 if (dpProps_) {
1971 dpProps_->dp_F_.fill(0.0);
1972 double vcn(0.0), vcs(0.0);
1973 // Calculate the damping coefficients.
1974 vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kTran_)));
1975 vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(kTran_/kRatio_)));
1976 // First damp the shear components
1977 for (int i = 1; i<dim; ++i)
1978 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1979 // Damp the normal component
1980 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1981 // Need to change behavior based on the dp_mode.
1982 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1983 // Limit in tension if not bonded.
1984 if (dpProps_->dp_F_.x() + sn_F_.x() < 0)
1985 dpProps_->dp_F_.rx() = -sn_F_.rx();
1986 }
1987 if (sn_S_ && dpProps_->dp_mode_ > 1) {
1988 // Limit in shear if not sliding.
1989 double dfn = dpProps_->dp_F_.rx();
1990 dpProps_->dp_F_.fill(0.0);
1991 dpProps_->dp_F_.rx() = dfn;
1992 }
1993 }
1994
1995 //Compute energies if energy tracking has been enabled.
1996 if (state->trackEnergy_) {
1997 assert(energies_);
1998 energies_->estrain_ = 0.0;
1999 if (kTran_)
2000 // Calculate the strain energy.
2001 energies_->estrain_ = 0.5*sn_F_.x()*sn_F_.x() / kTran_;
2002 if (kTran_) {
2003 DVect s = sn_F_;
2004 s.rx() = 0.0;
2005 double smag2 = s.mag2();
2006 // Add the shear component of the strain energy.
2007 energies_->estrain_ += 0.5*smag2 / (kTran_/kRatio_);
2008
2009 if (sn_S_) {
2010 // If sliding calculate the slip energy and accumulate it.
2011 sn_F_old.rx() = 0.0;
2012 DVect avg_F_s = (s + sn_F_old)*0.5;
2013 DVect u_s_el = (s - sn_F_old) / (kTran_/kRatio_);
2014 DVect u_s(0.0);
2015 for (int i = 1; i < dim; ++i)
2016 u_s.rdof(i) = trans.dof(i);
2017 energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
2018 }
2019 }
2020 // Add the bending/twisting resistance energy contributions.
2021 if (kb) {
2022 DAVect tmp = sn_M_;
2023#ifdef THREED
2024 tmp.rx() = 0.0;
2025#endif
2026 energies_->estrain_ += 0.5*tmp.mag2() / kb;
2027 if (sn_BS_) {
2028 // accumulate bending slip energy.
2029 DAVect tmp_old = sn_M_old;
2030#ifdef THREED
2031 tmp_old.rx() = 0.0;
2032#endif
2033 DAVect avg_M = (tmp + tmp_old)*0.5;
2034 DAVect t_s_el = (tmp - tmp_old) / kb;
2035 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
2036 }
2037 }
2038#ifdef THREED
2039 if (kt) {
2040 double mt = std::abs(sn_M_.x());
2041 energies_->estrain_ += 0.5*mt*mt / kt;
2042 if (sn_TS_) {
2043 // accumulate twisting slip energy.
2044 DAVect tmp(0.0);
2045 DAVect tmp_old(0.0);
2046 tmp.rx() = sn_M_.x();
2047 tmp_old.rx() = sn_M_old.x();
2048 DAVect avg_M = (tmp + tmp_old)*0.5;
2049 DAVect t_s_el = (tmp - tmp_old) / kt;
2050 energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
2051 }
2052 }
2053#endif
2054
2055 if (dpProps_) {
2056 // Calculate damping energy (accumulated) if the dashpots are active.
2057 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
2058 }
2059 }
2060
2061 // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
2062 assert(sn_F_ == sn_F_);
2063 assert(fictForce_ == fictForce_);
2064 assert(sn_M_ == sn_M_);
2065 return true;
2066 }
2067
2068 void ContactModelRBSN::setForce(const DVect &v,IContact *c) {
2069 sn_F_ = v;
2070 fictForce_ = DVect(0.0);
2071 forceSet_ = true;
2072 if (v.x() > 0)
2073 rgap_ = c->getGap() + v.x() / kTran_;
2074 }
2075
2076 void ContactModelRBSN::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
2077 // Only called for contacts with wall facets when the wall resolution scheme
2078 // is set to full!
2079 // Only do something if the contact model is of the same type
2080 if (equal(old->getContactModel()->getName(),"springnetwork") && !isBonded()) {
2081 ContactModelRBSN *oldCm = (ContactModelRBSN *)old;
2082#ifdef THREED
2083 // Need to rotate just the shear component from oldSystem to newSystem
2084
2085 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
2086 DVect axis = oldSystem.e1() & newSystem.e1();
2087 double c, ang, s;
2088 DVect re2;
2089 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
2090 axis = axis.unit();
2091 c = oldSystem.e1()|newSystem.e1();
2092 if (c > 0)
2093 c = std::min(c,1.0);
2094 else
2095 c = std::max(c,-1.0);
2096 ang = acos(c);
2097 s = sin(ang);
2098 double t = 1. - c;
2099 DMatrix<3,3> rm;
2100 rm.get(0,0) = t*axis.x()*axis.x() + c;
2101 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
2102 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
2103 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
2104 rm.get(1,1) = t*axis.y()*axis.y() + c;
2105 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
2106 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
2107 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
2108 rm.get(2,2) = t*axis.z()*axis.z() + c;
2109 re2 = rm*oldSystem.e2();
2110 }
2111 else
2112 re2 = oldSystem.e2();
2113 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
2114 axis = re2 & newSystem.e2();
2115 DVect2 tpf;
2116 DVect2 tpm;
2117 DMatrix<2,2> m;
2118 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
2119 axis = axis.unit();
2120 c = re2|newSystem.e2();
2121 if (c > 0)
2122 c = std::min(c,1.0);
2123 else
2124 c = std::max(c,-1.0);
2125 ang = acos(c);
2126 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
2127 ang *= -1;
2128 s = sin(ang);
2129 m.get(0,0) = c;
2130 m.get(1,0) = s;
2131 m.get(0,1) = -m.get(1,0);
2132 m.get(1,1) = m.get(0,0);
2133 tpf = m*DVect2(oldCm->sn_F_.y(),oldCm->sn_F_.z());
2134 tpm = m*DVect2(oldCm->sn_M_.y(),oldCm->sn_M_.z());
2135 } else {
2136 m.get(0,0) = 1.;
2137 m.get(0,1) = 0.;
2138 m.get(1,0) = 0.;
2139 m.get(1,1) = 1.;
2140 tpf = DVect2(oldCm->sn_F_.y(),oldCm->sn_F_.z());
2141 tpm = DVect2(oldCm->sn_M_.y(),oldCm->sn_M_.z());
2142 }
2143 DVect pforce = DVect(0,tpf.x(),tpf.y());
2144 //DVect pm = DVect(0,tpm.x(),tpm.y());
2145#else
2146 oldSystem;
2147 newSystem;
2148 DVect pforce = DVect(0,oldCm->sn_F_.y());
2149 //DVect pm = DVect(0,oldCm->sn_M_.y());
2150#endif
2151 for (int i=1; i<dim; ++i)
2152 sn_F_.rdof(i) += pforce.dof(i);
2153 if (sn_mode_ && oldCm->sn_mode_)
2154 sn_F_.rx() = oldCm->sn_F_.x();
2155 oldCm->sn_F_ = DVect(0.0);
2156 oldCm->sn_M_ = DAVect(0.0);
2157 if (dpProps_ && oldCm->dpProps_) {
2158#ifdef THREED
2159 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
2160 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
2161#else
2162 pforce = oldCm->dpProps_->dp_F_;
2163#endif
2164 dpProps_->dp_F_ += pforce;
2165 oldCm->dpProps_->dp_F_ = DVect(0.0);
2166 }
2167 if(oldCm->getEnergyActivated()) {
2168 activateEnergy();
2169 energies_->estrain_ = oldCm->energies_->estrain_;
2170 energies_->edashpot_ = oldCm->energies_->edashpot_;
2171 energies_->eslip_ = oldCm->energies_->eslip_;
2172 oldCm->energies_->estrain_ = 0.0;
2173 oldCm->energies_->edashpot_ = 0.0;
2174 oldCm->energies_->eslip_ = 0.0;
2175 }
2176 rgap_ = oldCm->rgap_;
2177 }
2178 assert(sn_F_ == sn_F_);
2179 }
2180
2181 void ContactModelRBSN::setNonForcePropsFrom(IContactModel *old) {
2182 // Only called for contacts with wall facets when the wall resolution scheme
2183 // is set to full!
2184 // Only do something if the contact model is of the same type
2185 if (equal(old->getName(),"springnetwork") && !isBonded()) {
2186 ContactModelRBSN *oldCm = (ContactModelRBSN *)old;
2187
2188 fictForce_ = oldCm->fictForce_;
2189 sn_F_ = oldCm->sn_F_;
2190 sn_sdisp_ = oldCm->sn_sdisp_;
2191 sn_M_ = oldCm->sn_M_;
2192 kRot_ = oldCm->kRot_;
2193 kTran_ = oldCm->kTran_;
2194 kRatio_ = oldCm->kRatio_;
2195 E_ = oldCm->E_;
2196 poisson_ = oldCm->poisson_;
2197 fric_ = oldCm->fric_;
2198 sn_bmul_ = oldCm->sn_bmul_;
2199 sn_tmul_ = oldCm->sn_tmul_;
2200 sn_rmul_ = oldCm->sn_rmul_;
2201 userArea_ = oldCm->userArea_;
2202 rgap_ = oldCm->rgap_;
2203 sn_fa_ = oldCm->sn_fa_;
2204 sn_mcf_ = oldCm->sn_mcf_;
2205 sn_dil_ = oldCm->sn_dil_;
2206 sn_dilzero_ = oldCm->sn_dilzero_;
2207 transTen_ = oldCm->transTen_;
2208 sn_elong_ = oldCm->sn_elong_;
2209 sn_ndisp_ = oldCm->sn_ndisp_;
2210 sn_mode_ = oldCm->sn_mode_;
2211 sn_state_ = oldCm->sn_state_;
2212 poisOffDiag_ = oldCm->poisOffDiag_;
2213 sn_S_ = oldCm->sn_S_;
2214 sn_BS_ = oldCm->sn_BS_;
2215 sn_TS_ = oldCm->sn_TS_;
2216 forceSet_ = oldCm->forceSet_;
2217 sn_heal_ = oldCm->sn_heal_;
2218 tenTable_ = oldCm->tenTable_;
2219 cohTable_ = oldCm->cohTable_;
2220 if (oldCm->hasDamping()) {
2221 if (!dpProps_)
2222 dpProps_ = NEW dpProps();
2223 dp_nratio(oldCm->dp_nratio());
2224 dp_sratio(oldCm->dp_sratio());
2225 dp_mode(oldCm->dp_mode());
2226 dp_F(oldCm->dp_F());
2227 }
2228 }
2229 }
2230
2231 DVect ContactModelRBSN::getForce() const {
2232 DVect ret(sn_F_);
2233 ret += fictForce_;
2234 if (dpProps_)
2235 ret += dpProps_->dp_F_;
2236 return ret;
2237 }
2238
2239 DAVect ContactModelRBSN::getMomentOn1(const IContactMechanical *c) const {
2240 DAVect ret(sn_M_);
2241 if (c) {
2242 DVect force = getForce();
2243 c->updateResultingTorqueOn1Local(force,&ret);
2244 }
2245 return ret;
2246 }
2247
2248 DAVect ContactModelRBSN::getMomentOn2(const IContactMechanical *c) const {
2249 DAVect ret(sn_M_);
2250 if (c) {
2251 DVect force = getForce();
2252 c->updateResultingTorqueOn2Local(force,&ret);
2253 }
2254 return ret;
2255 }
2256
2257 DAVect ContactModelRBSN::getModelMomentOn1() const {
2258 DAVect ret(sn_M_);
2259 return ret;
2260 }
2261
2262 DAVect ContactModelRBSN::getModelMomentOn2() const {
2263 DAVect ret(sn_M_);
2264 return ret;
2265 }
2266
2267 void ContactModelRBSN::objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *ret) const {
2268 ret->clear();
2269 ret->push_back({"kn",ScalarInfo});
2270 ret->push_back({"ks",ScalarInfo});
2271 ret->push_back({"krot",VectorInfo});
2272 ret->push_back({"fric",ScalarInfo});
2273 ret->push_back({"sn_bmul",ScalarInfo});
2274 ret->push_back({"sn_tmul",ScalarInfo});
2275 ret->push_back({"sn_mode",ScalarInfo});
2276 ret->push_back({"sn_force",VectorInfo});
2277 ret->push_back({"sn_moment",VectorInfo});
2278 ret->push_back({"sn_slip",ScalarInfo});
2279 ret->push_back({"sn_slipb",ScalarInfo});
2280 ret->push_back({"sn_slipt",ScalarInfo});
2281 ret->push_back({"sn_rmul",ScalarInfo});
2282 ret->push_back({"sn_radius",ScalarInfo});
2283 ret->push_back({"emod",ScalarInfo});
2284 ret->push_back({"kratio",ScalarInfo});
2285 ret->push_back({"dp_nratio",ScalarInfo});
2286 ret->push_back({"dp_sratio",ScalarInfo});
2287 ret->push_back({"dp_mode",ScalarInfo});
2288 ret->push_back({"dp_force",VectorInfo});
2289 ret->push_back({"sn_state",ScalarInfo});
2290 ret->push_back({"sn_ten",ScalarInfo});
2291 ret->push_back({"sn_shear",ScalarInfo});
2292 ret->push_back({"sn_coh",ScalarInfo});
2293 ret->push_back({"sn_fa",ScalarInfo});
2294 ret->push_back({"sn_mcf",ScalarInfo});
2295 ret->push_back({"sn_sigma",ScalarInfo});
2296 ret->push_back({"sn_tau",ScalarInfo});
2297 ret->push_back({"sn_area",ScalarInfo});
2298 ret->push_back({"user_area",ScalarInfo});
2299 ret->push_back({"rgap",ScalarInfo});
2300 ret->push_back({"sn_pois_force",VectorInfo});
2301 ret->push_back({"sn_pois",ScalarInfo});
2302 ret->push_back({"sn_non_diag",ScalarInfo});
2303 ret->push_back({"sn_cohres",ScalarInfo});
2304 ret->push_back({"sn_dil",ScalarInfo});
2305 ret->push_back({"sn_dilzero",ScalarInfo});
2306 ret->push_back({"sn_ndisp",ScalarInfo});
2307 ret->push_back({"sn_sdisp",VectorInfo});
2308 ret->push_back({"sn_cohdc",ScalarInfo});
2309 ret->push_back({"sn_heal",ScalarInfo});
2310 ret->push_back({"sn_tendc",ScalarInfo});
2311 ret->push_back({"sn_tabpos",ScalarInfo});
2312 ret->push_back({"sn_porp",ScalarInfo});
2313 ret->push_back({"sn_esigma",ScalarInfo});
2314
2315 }
2316
2317 void ContactModelRBSN::objectPropValues(std::vector<double> *ret,const IContact *con) const {
2318 FP_S;
2319 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
2320 ret->clear();
2321 ret->push_back(kTran_);
2322 ret->push_back(kTran_/kRatio_);
2323 ret->push_back(kRot_.mag());
2324 ret->push_back(fric_);
2325 ret->push_back(sn_bmul_);
2326 ret->push_back(sn_tmul_);
2327 ret->push_back(sn_mode_);
2328 ret->push_back(sn_F_.mag());
2329 ret->push_back(sn_M_.mag());
2330 ret->push_back(sn_S_);
2331 ret->push_back(sn_BS_);
2332 ret->push_back(sn_TS_);
2333 ret->push_back(sn_rmul_);
2334 double Cmax1 = c->getEnd1Curvature().y();
2335 double Cmax2 = c->getEnd2Curvature().y();
2336 double d;
2337 if (!userArea_)
2338 d= sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
2339 else {
2340#ifdef THREED
2341 d= std::sqrt(userArea_ / dPi);
2342#else
2343 d= userArea_ / 2.0;
2344#endif
2345 }
2346 ret->push_back(d);
2347 ret->push_back(E_);
2348 ret->push_back(kRatio_);
2349 ret->push_back(dpProps_ ? dpProps_->dp_nratio_ : 0);
2350 ret->push_back(dpProps_ ? dpProps_->dp_sratio_ : 0);
2351 ret->push_back(dpProps_ ? dpProps_->dp_mode_ : 0);
2352 ret->push_back(dpProps_ ? dpProps_->dp_F_.mag() : 0);
2353 ret->push_back(sn_state_);
2354 ret->push_back(sn_Ten());
2355 ret->push_back(shearStrength(computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x()));
2356 ret->push_back(cohTable_[0].x());
2357 ret->push_back(std::atan(sn_fa_)/dDegrad);
2358 ret->push_back(sn_mcf_);
2359 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x());
2360 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).y());
2361 ret->push_back(userArea_ ? userArea_ : computeGeomData(c->getEnd1Curvature(),c->getEnd2Curvature()).x());
2362 ret->push_back(userArea_);
2363 ret->push_back(rgap_);
2364 ret->push_back(fictForce_.mag());
2365 ret->push_back(poisson_);
2366 ret->push_back(poisOffDiag_ == false ? 0 : 1);
2367 ret->push_back(sn_cohres_);
2368 ret->push_back(std::atan(sn_dil_)/dDegrad);
2369 ret->push_back(sn_dilzero_);
2370 ret->push_back(sn_ndisp_);
2371 ret->push_back(sn_sdisp_.mag());
2372 ret->push_back(sn_cohdc());
2373 ret->push_back(sn_heal_);
2374 ret->push_back(sn_tendc());
2375 ret->push_back(sn_tabPos_+1);
2376 ret->push_back(sn_por_);
2377 ret->push_back(SMax(c->getEnd1Curvature(),c->getEnd2Curvature()).x() + sn_por_);
2378 FP_S;
2379 }
2380
2381 DVect3 ContactModelRBSN::computeGeomData(const DVect2 &end1Curvature,const DVect2 &end2Curvature) const {
2382 double Cmax1 = end1Curvature.y();
2383 double Cmax2 = end2Curvature.y();
2384 double br = sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
2385 if (userArea_)
2386#ifdef THREED
2387 br = std::sqrt(userArea_ / dPi);
2388#else
2389 br = userArea_ / 2.0;
2390#endif
2391 double br2 = br * br;
2392#ifdef THREED
2393 double area = dPi * br2;
2394 double bi = 0.25*area*br2;
2395#else
2396 double area = 2.0*br;
2397 double bi = 2.0*br*br2 / 3.0;
2398#endif
2399 return DVect3(area, bi, br);
2400 }
2401
2402 DVect2 ContactModelRBSN::SMax(const DVect2 &e1c,const DVect2 &e2c) const {
2403 DVect3 data = computeGeomData(e1c,e2c);
2404 double area = data.x();
2405 double bi = data.y();
2406 double br = data.z();
2407 /* maximum stresses */
2408 //double nsmax0 = -(totalForce.x() / area) + sn_mcf_* sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z()) * br / bi;
2409 double dbend = sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z());
2410 dbend *= sn_mcf_;
2411 double dtwist = sn_M_.x();
2412 DVect bfs(sn_F_);
2413 bfs.rx() = 0.0;
2414 double dbfs = bfs.mag();
2415 double nsmax = -((sn_F_.x()+fictForce_.x()) / area) + dbend * br / bi;
2416 double ssmax = dbfs / area + std::abs(dtwist) * 0.5* br / bi;
2417 return DVect2(nsmax, ssmax);
2418 }
2419
2420 double ContactModelRBSN::shearStrength(const double &area) const {
2421 double sig = -1.0*(sn_F_.x() + fictForce_.x()) / area;
2422 double nstr = (sn_state_ > 2 && sn_state_ != 6) ? tenTable_[0].x() : 0.0;
2423 return sig <= nstr ? cohTable_[0].x() - sn_fa_*sig
2424 : cohTable_[0].x() - sn_fa_*nstr;
2425 }
2426
2427
2428 double ContactModelRBSN::strainEnergy(double kn,double ks,double kb,double kt) const {
2429 double ret(0.0);
2430 if (kn)
2431 ret = 0.5 * (sn_F_.x()+fictForce_.x()) * (sn_F_.x()+fictForce_.x()) / kn;
2432 if (ks) {
2433 DVect tmp = sn_F_ + fictForce_;
2434 tmp.rx() = 0.0;
2435 double smag2 = tmp.mag2();
2436 ret += 0.5 * smag2 / ks;
2437 }
2438
2439 if (kt)
2440 ret += 0.5 * sn_M_.x() * sn_M_.x() / kt;
2441 if (kb) {
2442 DAVect tmp = sn_M_;
2443#ifdef THREED
2444 tmp.rx() = 0.0;
2445 double smag2 = tmp.mag2();
2446#else
2447 double smag2 = tmp.z() * tmp.z();
2448#endif
2449 ret += 0.5 * smag2 / kb;
2450 }
2451 return ret;
2452 }
2453
2454} // namespace cmodelsxd
2455// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Jun 07, 2025 |
