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