Smooth-Joint Model Implementation
See this page for the documentation of this contact model.
contactmodelsmoothjoint.h
1#pragma once
2// contactmodelsmoothjoint.h
3
4#include "contactmodel/src/contactmodelmechanical.h"
5
6#ifdef SMOOTHJOINT_LIB
7# define SMOOTHJOINT_EXPORT EXPORT_TAG
8#elif defined(NO_MODEL_IMPORT)
9# define SMOOTHJOINT_EXPORT
10#else
11# define SMOOTHJOINT_EXPORT IMPORT_TAG
12#endif
13
14namespace cmodelsxd {
15 using namespace itasca;
16
17 class ContactModelSmoothJoint : public ContactModelMechanical {
18 public:
19 enum PropertyKeys {
20 kwKn=1
21 , kwKs
22 , kwFric
23 , kwDA
24 , kwState
25 , kwTen
26 , kwBCoh
27 , kwBFA
28 , kwShear
29 , kwRMul
30 , kwLarge
31 , kwFn
32 , kwFs
33 , kwGap
34 , kwNorm
35 , kwArea
36 , kwRad
37 , kwUn
38 , kwUs
39 , kwSlip
40 , kwDip
41 , kwDD
42 };
43
44 enum FishCallEvents { fBondBreak=0, fSlipChange };
45
46 SMOOTHJOINT_EXPORT ContactModelSmoothJoint();
47 SMOOTHJOINT_EXPORT virtual ~ContactModelSmoothJoint();
48 virtual void copy(const ContactModel *c) override;
49 virtual void archive(ArchiveStream &);
50 virtual QString getName() const { return "smoothjoint"; }
51 virtual void setIndex(int i) { index_=i;}
52 virtual int getIndex() const {return index_;}
53
54 virtual QString getProperties() const {
55 return "sj_kn"
56 ",sj_ks"
57 ",sj_fric"
58 ",sj_da"
59 ",sj_state"
60 ",sj_ten"
61 ",sj_coh"
62 ",sj_fa"
63 ",sj_shear"
64 ",sj_rmul"
65 ",sj_large"
66 ",sj_fn"
67 ",sj_fs"
68 ",sj_gap"
69 ",sj_unorm"
70 ",sj_area"
71 ",sj_radius"
72 ",sj_un"
73 ",sj_us"
74 ",sj_slip"
75 ",dip"
76#ifdef THREED
77 ",ddir"
78#endif
79 ;
80 }
81
82 virtual QString getFishCallEvents() const { return "bond_break,slip_change"; }
83 virtual QVariant getProperty(uint32 i,const IContact *con=0) const;
84 virtual bool getPropertyGlobal(uint32 i) const;
85 virtual bool setProperty(uint32 i,const QVariant &v,IContact *con=0);
86 virtual bool getPropertyReadOnly(uint32 i) const;
87 virtual bool supportsInheritance(uint32 i) const;
88 virtual bool getInheritance(uint32 i) const { assert(i<32); uint32 mask = to<uint32>(1 << i); return (inheritanceField_ & mask) ? true : false; }
89 virtual void setInheritance(uint32 i,bool b) { assert(i<32); uint32 mask = to<uint32>(1 << i); if (b) inheritanceField_ |= mask; else inheritanceField_ &= ~mask; }
90
91 enum MethodKeys {
92 kwSetForce=1,kwBond,kwUnbond
93 };
94
95 virtual QString getMethods() const {
96 return "sj_setforce,bond,unbond";
97 }
98
99 virtual QString getMethodArguments(uint32 i) const;
100 virtual bool setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0); // Base 1 - returns true if timestep contributions need to be updated
101
102 virtual uint32 getMinorVersion() const;
103
104 enum EnergyKeys { kwEStrain=1,kwESlip};
105 virtual QString getEnergies() const { return "energy-strain,energy-slip";}
106 virtual double getEnergy(uint32 i) const; // Base 1
107 virtual bool getEnergyAccumulate(uint32 i) const; // Base 1
108 virtual void setEnergy(uint32 i,const double &d); // Base 1
109 virtual void activateEnergy() { if (energies_) return; energies_ = NEW Energies();}
110 virtual bool getEnergyActivated() const {return (energies_ !=0);}
111
112 virtual bool validate(ContactModelMechanicalState *state,const double ×tep);
113 virtual bool endPropertyUpdated(const QString &name,const IContactMechanical *c);
114 virtual bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep);
115 virtual DVect2 getEffectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
116 virtual DAVect getEffectiveRotationalStiffness() const {return DAVect(0.0);}
117
118 virtual ContactModelSmoothJoint *clone() const override { return NEW ContactModelSmoothJoint(); }
119 virtual double getActivityDistance() const {return 0.0;}
120 virtual void resetForcesAndMoments() { sj_Fn(0.0); sj_Fs(DVect(0.0)); }
121 virtual void setForce(const DVect &v,IContact *);
122 virtual void setArea(const double &) { throw Exception("The setArea method cannot be used with this contact model."); }
123 virtual double getArea() const { return 0.0; }
124
125 virtual bool isOKToDelete() const { return (!isBonded() && sj_large_); }
126
127 virtual bool checkActivity(const double &gap) { return ( !sj_large_ || gap <= 0.0 || isBonded()); }
128 virtual bool isSliding() const { return isjSliding_; }
129 virtual bool isBonded() const { return sj_state_ == 3; }
130 virtual void unbond() { sj_state_ = 0; }
131 virtual bool hasNormal() const { return true; }
132 virtual DVect3 getNormal() const { return toVect3(sj_unorm_); }
133
134 void sj_kn(const double &d) {sj_kn_ = d;}
135 void sj_ks(const double &d) {sj_ks_ = d;}
136 void sj_fric(const double &d) {sj_fric_ = d;}
137 void sj_da(const double &d) {sj_da_ = d;}
138 void sj_state(const double &d) {sj_state_ = d;}
139 void sj_bns(const double &d) {sj_bns_ = d;}
140 void sj_bcoh(const double &d) {sj_bcoh_ = d;}
141 void sj_bfa(const double &d) {sj_bfa_ = d;}
142 void dip(const double &d) {dip_ = d;}
143 void sj_rmul(const double &d) {sj_rmul_ = d;}
144 void sj_large(bool b) {sj_large_ = b;}
145
146 const double & sj_kn() const {return sj_kn_;}
147 const double & sj_ks() const {return sj_ks_;}
148 const double & sj_fric() const {return sj_fric_;}
149 const double & sj_da() const {return sj_da_;}
150 int sj_state() const {return sj_state_;}
151 const double & sj_bns() const {return sj_bns_;}
152 const double & sj_bcoh() const {return sj_bcoh_;}
153 const double & sj_bfa() const {return sj_bfa_;}
154 const double & dip() const {return dip_;}
155 const double & sj_rmul() const {return sj_rmul_;}
156 bool sj_large() const {return sj_large_;}
157
158#ifdef THREED
159 const double & dd() const {return dd_;}
160 void dd(const double &d) {dd_ =d;}
161#endif
162
163 void sj_unorm(const DVect &v) {sj_unorm_ = v;}
164 void sj_A(const double &d) {sj_A_ = d;}
165 void sj_rad(const double &d) {sj_rad_ = d;}
166 void sj_gap(const double &d) {sj_gap_ = d;}
167 void sj_Un(const double &d) {sj_Un_ = d;}
168 void sj_Us(const DVect &v) {sj_Us_ = v;}
169 void sj_Fn(const double &d) {sj_Fn_ = d;}
170 void sj_Fs(const DVect &v) {sj_Fs_ = v;}
171 void isjFlip(bool b) {isjFlip_ = b;}
172 void isjSliding(bool b) {isjSliding_ = b;}
173
174 const DVect & sj_unorm() const {return sj_unorm_;}
175 const double & sj_A() const {return sj_A_;}
176 const double & sj_rad() const {return sj_rad_;}
177 const double & sj_gap() const {return sj_gap_;}
178 const double & sj_Un() const {return sj_Un_;}
179 const DVect & sj_Us() const {return sj_Us_;}
180 const double & sj_Fn() const {return sj_Fn_;}
181 const DVect & sj_Fs() const {return sj_Fs_;}
182 bool isjFlip() const {return isjFlip_;}
183 bool isjSliding() const {return isjSliding_;}
184
185 bool hasEnergies() const {return energies_ ? true:false;}
186 double estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
187 void estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
188 double eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
189 void eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
190
191 uint32 inheritanceField() const {return inheritanceField_;}
192 void inheritanceField(uint32 i) {inheritanceField_ = i;}
193
194 const DVect2 & effectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
195 void effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
196
197 bool geomRecomp() const {return geomRecomp_ ;}
198 void geomRecomp(bool b) {geomRecomp_ = b;}
199
200 /// Return the total force that the contact model holds.
201 virtual DVect getForce(const IContactMechanical *) const;
202
203 /// Return the total moment on 1 that the contact model holds
204 virtual DAVect getMomentOn1(const IContactMechanical *) const;
205
206 /// Return the total moment on 1 that the contact model holds
207 virtual DAVect getMomentOn2(const IContactMechanical *) const;
208
209 private:
210 static int index_;
211
212 struct Energies {
213 Energies() : estrain_(0.0), eslip_(0.0) {}
214 double estrain_; // elastic energy stored in contact
215 double eslip_; // work dissipated by friction
216 };
217
218 bool updateKn(const IContactMechanical *con);
219 bool updateKs(const IContactMechanical *con);
220 bool updateFric(const IContactMechanical *con);
221 void updateAreaAndNormal(ContactModelMechanicalState *state);
222 void updateAreaAndNormalContact(const IContact *c);
223 double calcBSS() const;
224
225 void updateEffectiveTranslationalStiffness();
226
227 // property set fields
228 uint32 inheritanceField_;
229
230 // smooth joint model - contact model properties
231 double sj_kn_; // normal stiffness
232 double sj_ks_; // shear stiffness
233 double sj_fric_; // Coulomb friction coefficient
234 double sj_da_; // Dilation angle (stored in radians, returned/set in degrees)
235 int sj_state_; // Bond state - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (bonded)
236 double sj_bns_; // Bond normal (tensile) strength
237 double sj_bcoh_; // Bonded system cohesion
238 double sj_bfa_; // Bonded system friction angle (stored in radians, returned/set in degrees)
239 double dip_; // Dip angle (stored in radians, returned/set in degrees)
240 double sj_rmul_; // Radius multiplier
241 bool sj_large_; // Large strain indicator
242
243 // Internal state variables
244 DVect sj_unorm_; // Normal to the plane
245 double sj_A_; // Cross-sectional area
246 double sj_rad_; // Radius
247 double sj_gap_; // Gap - this can be user modified
248 double sj_Un_; // Normal displacement
249 DVect sj_Us_; // Shear displacement
250 double sj_Fn_; // Normal force
251 DVect sj_Fs_; // Shear force
252 bool isjFlip_; // In order to flip the order
253 bool isjSliding_;// True if this is sliding
254
255 DVect2 effectiveTranslationalStiffness_;
256 bool geomRecomp_; // If true then there must be a validate before FD!
257#ifdef THREED
258 double dd_; // Dip direction (stored in radians, returned/set in degrees)
259#endif
260 //DAVect effectiveRotationalStiffness_;
261
262 Energies *energies_;
263 };
264} // namespace itascaxd
265// EoF
contactmodelsmoothjoint.cpp
1// contactmodelsmoothjoint.cpp
2#include "contactmodelsmoothjoint.h"
3
4#include "../version.txt"
5#include "fish/src/parameter.h"
6#include "utility/src/tptr.h"
7#include "base/src/basetoqt.h"
8
9#include "kernel/interface/iprogram.h"
10#include "module/interface/icontact.h"
11#include "module/interface/icontactmechanical.h"
12#include "module/interface/ifishcalllist.h"
13#include "module/interface/ipiece.h"
14#include "module/interface/ipiecemechanical.h"
15
16#ifdef SMOOTHJOINT_LIB
17#ifdef _WIN32
18 int __stdcall DllMain(void *,unsigned, void *)
19 {
20 return 1;
21 }
22#endif
23
24 extern "C" EXPORT_TAG const char *getName()
25 {
26#if DIM==3
27 return "contactmodelmechanical3dsmoothjoint";
28#else
29 return "contactmodelmechanical2dsmoothjoint";
30#endif
31 }
32
33 extern "C" EXPORT_TAG unsigned getMajorVersion()
34 {
35 return MAJOR_VERSION;
36 }
37
38 extern "C" EXPORT_TAG unsigned getMinorVersion()
39 {
40 return MINOR_VERSION;
41 }
42
43 extern "C" EXPORT_TAG void *createInstance()
44 {
45 cmodelsxd::ContactModelSmoothJoint *m = NEW cmodelsxd::ContactModelSmoothJoint();
46 return (void *)m;
47 }
48#endif // SMOOTHJOINT_EXPORTS
49
50namespace cmodelsxd {
51 static const uint32 sjKnMask = 0x00002; // Base 1!
52 static const uint32 sjKsMask = 0x00004;
53 static const uint32 sjFricMask = 0x00008;
54
55 using namespace itasca;
56
57 int ContactModelSmoothJoint::index_ = -1;
58 uint32 ContactModelSmoothJoint::getMinorVersion() const { return MINOR_VERSION; }
59
60 ContactModelSmoothJoint::ContactModelSmoothJoint() : inheritanceField_(sjKnMask|sjKsMask|sjFricMask)
61 , sj_kn_(0.0), sj_ks_(0.0), sj_fric_(0.0)
62 , sj_da_(0.0), sj_state_(0), sj_bns_(0.0)
63 , sj_bcoh_(0.0), sj_bfa_(0.0), dip_(0.0)
64 , sj_rmul_(1.0), sj_large_(false)
65 , sj_unorm_(0.0), sj_A_(0.0), sj_rad_(0.0)
66 , sj_gap_(0.0), sj_Un_(0.0), sj_Us_(0.0)
67 , sj_Fn_(0.0), sj_Fs_(0.0), isjFlip_(false)
68 , isjSliding_(false)
69 , effectiveTranslationalStiffness_(0.0)
70 , geomRecomp_(true)
71#ifdef THREED
72 , dd_(0.0)
73#endif
74 , energies_(0) {
75// setFromParent(ContactModelMechanicalList::instance()->find(getName()));
76 }
77
78 ContactModelSmoothJoint::~ContactModelSmoothJoint() {
79 if (energies_)
80 delete energies_;
81 }
82
83 void ContactModelSmoothJoint::archive(ArchiveStream &stream) {
84 stream & sj_kn_;
85 stream & sj_ks_;
86 stream & sj_fric_;
87 stream & sj_da_;
88 stream & sj_state_;
89 stream & sj_bns_;
90 stream & sj_bcoh_;
91 stream & sj_bfa_;
92 stream & dip_;
93 stream & sj_rmul_;
94 stream & sj_large_;
95 stream & sj_unorm_;
96 stream & sj_A_;
97 stream & sj_rad_;
98 stream & sj_gap_;
99 stream & sj_Un_;
100 stream & sj_Us_;
101 stream & sj_Fn_;
102 stream & sj_Fs_;
103 stream & isjFlip_;
104 stream & isjSliding_;
105#ifdef THREED
106 stream & dd_;
107#endif
108 if (stream.getArchiveState()==ArchiveStream::Save) {
109 bool b = false;
110 if (energies_) {
111 b = true;
112 stream & b;
113 stream & energies_->estrain_;
114 stream & energies_->eslip_;
115 }
116 else
117 stream & b;
118 } else {
119 bool b(false);
120 stream & b;
121 if (b) {
122 if (!energies_)
123 energies_ = NEW Energies();
124 stream & energies_->estrain_;
125 stream & energies_->eslip_;
126 }
127 }
128
129 stream & inheritanceField_;
130 stream & geomRecomp_;
131 stream & effectiveTranslationalStiffness_;
132
133 }
134
135 void ContactModelSmoothJoint::copy(const ContactModel *cm) {
136 ContactModelMechanical::copy(cm);
137 const ContactModelSmoothJoint *in = dynamic_cast<const ContactModelSmoothJoint*>(cm);
138 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
139 sj_kn(in->sj_kn());
140 sj_ks(in->sj_ks());
141 sj_fric(in->sj_fric());
142 sj_da(in->sj_da());
143 sj_state(in->sj_state());
144 sj_bns(in->sj_bns());
145 sj_bcoh(in->sj_bcoh());
146 sj_bfa(in->sj_bfa());
147 dip(in->dip());
148 sj_rmul(in->sj_rmul());
149 sj_large(in->sj_large());
150#ifdef THREED
151 dd(in->dd());
152#endif
153 sj_unorm(in->sj_unorm());
154 sj_A(in->sj_A());
155 sj_rad(in->sj_rad());
156 sj_gap(in->sj_gap());
157 sj_Un(in->sj_Un());
158 sj_Us(in->sj_Us());
159 sj_Fn(in->sj_Fn());
160 sj_Fs(in->sj_Fs());
161 isjFlip(in->isjFlip());
162 isjSliding(in->isjSliding());
163
164 if (in->hasEnergies()) {
165 if (!energies_)
166 energies_ = NEW Energies();
167 estrain(in->estrain());
168 eslip(in->eslip());
169 }
170 inheritanceField(in->inheritanceField());
171 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
172 geomRecomp(in->geomRecomp());
173 }
174
175 QVariant ContactModelSmoothJoint::getProperty(uint32 i,const IContact *) const {
176 QVariant var;
177 switch (i) {
178 case kwKn: return sj_kn_;
179 case kwKs: return sj_ks_;
180 case kwFric: return sj_fric_;
181 case kwDA: return sj_da_/dDegrad; // Returned in degrees
182 case kwState: return sj_state_;
183 case kwTen: return sj_bns_;
184 case kwBCoh: return sj_bcoh_;
185 case kwBFA: return sj_bfa_/dDegrad; // Returned in degrees
186 case kwShear: return calcBSS();
187 case kwDip: return dip_/dDegrad; // Returned in degrees
188#ifdef THREED
189 case kwDD: return dd_/dDegrad; // Returned in degrees
190#endif
191 case kwRMul: return sj_rmul_;
192 case kwLarge: return sj_large_;
193 case kwFn: return sj_Fn_;
194 case kwFs: {
195 var.setValue(sj_Fs_);
196 return var;
197 }
198 case kwGap: return sj_gap_;
199 case kwNorm: {
200 var.setValue(sj_unorm_);
201 return var;
202 }
203 case kwArea: return sj_A_;
204 case kwRad: return sj_rad_;
205 case kwUn: return sj_Un_;
206 case kwUs: {
207 var.setValue(sj_Us_);
208 return var;
209 }
210 case kwSlip: return isjSliding_;
211 }
212 assert(0);
213 return QVariant();
214 }
215
216 bool ContactModelSmoothJoint::getPropertyGlobal(uint32 ) const {
217 return true;
218 }
219
220 bool ContactModelSmoothJoint::setProperty(uint32 i,const QVariant &v,IContact *c) {
221 switch (i) {
222 case kwKn: {
223 if (!v.canConvert<double>())
224 throw Exception("sj_kn must be a double.");
225 double val = v.toDouble();
226 if (val<0.0)
227 throw Exception("Negative sj_kn not allowed.");
228 sj_kn_ = val;
229 updateEffectiveTranslationalStiffness();
230 return true;
231 }
232 case kwKs: {
233 if (!v.canConvert<double>())
234 throw Exception("sj_ks must be a double.");
235 double val = v.toDouble();
236 if (val<0.0)
237 throw Exception("Negative sj_ks not allowed.");
238 sj_ks_ = val;
239 updateEffectiveTranslationalStiffness();
240 return true;
241 }
242 case kwFric: {
243 if (!v.canConvert<double>())
244 throw Exception("sj_fric must be a double.");
245 double val = v.toDouble();
246 if (val<0.0)
247 throw Exception("Negative friction not allowed.");
248 sj_fric_ = val;
249 return false;
250 }
251 case kwDA: {
252 if (!v.canConvert<double>())
253 throw Exception("sj_da must be a double.");
254 sj_da_ = v.toDouble()*dDegrad; // Given in degrees
255 return false;
256 }
257 case kwState: {
258 if (!v.canConvert<uint32>())
259 throw Exception("sj_state must be must be an integer between 0 and 3.");
260 int val = v.toInt();
261 if (val<0 || val>3)
262 throw Exception("The bond state must be an integer between 0 and 3.");
263 sj_state_ = val;
264 return false;
265 }
266 case kwTen: {
267 if (!v.canConvert<double>())
268 throw Exception("sj_ten must be a double.");
269 double val = v.toDouble();
270 if (val<0.0)
271 throw Exception("Negative bond normal (tensile) strength not allowed.");
272 sj_bns_ = val;
273 return false;
274 }
275 case kwBCoh: {
276 if (!v.canConvert<double>())
277 throw Exception("sj_coh must be a double.");
278 double val = v.toDouble();
279 if (val<0.0)
280 throw Exception("Negative bond system cohesion not allowed.");
281 sj_bcoh_ = val;
282 return false;
283 }
284 case kwBFA: {
285 if (!v.canConvert<double>())
286 throw Exception("sj_bfa must be a double.");
287 sj_bfa_ = v.toDouble()*dDegrad; // Given in degrees
288 return false;
289 }
290 case kwDip: {
291 if (!v.canConvert<double>())
292 throw Exception("dip must be a double.");
293 dip_ = v.toDouble()*dDegrad; // Given in degrees
294 if (c)
295 updateAreaAndNormalContact(c);
296 else
297 geomRecomp_ = true;
298 return true;
299 }
300 case kwDD: {
301#ifdef THREED
302 if (!v.canConvert<double>())
303 throw Exception("ddir must be a double.");
304 dd_ = v.toDouble()*dDegrad; // Given in degrees
305 if (c)
306 updateAreaAndNormalContact(c);
307 else
308 geomRecomp_ = true;
309#endif
310 return true;
311 }
312 case kwRMul: {
313 if (!v.canConvert<double>())
314 throw Exception("sj_rmul must be a double.");
315 double val = v.toDouble();
316 if (val<0.0)
317 throw Exception("Negative radius multiplier not allowed.");
318 sj_rmul_ = val;
319 if (!geomRecomp_) geomRecomp_ = true;
320 return false;
321 }
322 case kwLarge: {
323 if (!v.canConvert<bool>())
324 throw Exception("Large-strain flag must be a boolean.");
325 sj_large_ = v.toBool();
326 return false;
327 }
328 case kwFn: {
329 if (!v.canConvert<double>())
330 throw Exception("sj_fn must be a double.");
331 sj_Fn_ = v.toDouble();
332 return false;
333 }
334 case kwFs: {
335 if (!v.canConvert<DVect>())
336 throw Exception("sj_fs must be a vector.");
337 sj_Fs_ = v.value<DVect>();
338 return false;
339 }
340 case kwGap: {
341 if (!v.canConvert<double>())
342 throw Exception("sj_gap must be a double.");
343 sj_gap_ = v.toDouble();
344 return false;
345 }
346
347
348 // Following are read-only !
349 case kwNorm:
350 case kwRad:
351 case kwShear:
352 case kwUn:
353 case kwUs:
354 case kwSlip:
355 return false;
356 }
357 assert(0);
358 return false;
359 }
360
361 bool ContactModelSmoothJoint::getPropertyReadOnly(uint32 i) const {
362 switch (i) {
363 case kwNorm:
364 case kwRad:
365 case kwUn:
366 case kwUs:
367 case kwSlip:
368 case kwShear:
369 case kwArea:
370 return true;
371 default:
372 break;
373 }
374 return false;
375 }
376
377 bool ContactModelSmoothJoint::supportsInheritance(uint32 i) const {
378 switch (i) {
379 case kwKn:
380 case kwKs:
381 case kwFric:
382 return true;
383 default:
384 break;
385 }
386 return false;
387 }
388
389 QString ContactModelSmoothJoint::getMethodArguments(uint32 i) const {
390 QString def = QString();
391 switch (i) {
392 case kwBond:
393 return "gap";
394 case kwUnbond:
395 return "gap";
396 }
397 return def;
398 }
399
400 bool ContactModelSmoothJoint::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
401 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
402 switch (i) {
403 case kwSetForce: {
404 DVect gf = c->getGlobalForce();
405 if (isjFlip_)
406 gf *= -1.0;
407 sj_Fn_ = gf | sj_unorm_;
408 return false;
409 }
410 case kwBond: {
411 if (sj_state_ == 3) return false;
412 double mingap = -1.0 * limits<double>::max();
413 double maxgap = 0;
414 if (vl.at(0).canConvert<double>())
415 maxgap = vl.at(0).toDouble();
416 else if (vl.at(0).canConvert<DVect2>()) {
417 DVect2 value = vl.at(0).value<DVect2>();
418 mingap = value.minComp();
419 maxgap = value.maxComp();
420 } else if (!vl.at(0).isNull())
421 throw Exception("gap value %1 not recognized in method bond in contact model %2.",vl.at(0),getName());
422 if ( sj_gap_ >= mingap && sj_gap_ <= maxgap) {
423 sj_state_ = 3;
424 return true;
425 }
426 return false;
427 }
428 case kwUnbond: {
429 if (sj_state_ == 0) return false;
430 double mingap = -1.0 * limits<double>::max();
431 double maxgap = 0;
432 if (vl.at(0).canConvert<double>())
433 maxgap = vl.at(0).toDouble();
434 else if (vl.at(0).canConvert<DVect2>()) {
435 DVect2 value = vl.at(0).value<DVect2>();
436 mingap = value.minComp();
437 maxgap = value.maxComp();
438 }
439 else if (!vl.at(0).isNull())
440 throw Exception("gap value %1 not recognized in method unbond in contact model %2.",vl.at(0),getName());
441 if ( sj_gap_ >= mingap && sj_gap_ <= maxgap) {
442 sj_state_ = 0;
443 return true;
444 }
445 return false;
446 }
447
448 }
449 return false;
450 }
451
452 double ContactModelSmoothJoint::getEnergy(uint32 i) const {
453 double ret(0.0);
454 if (!energies_)
455 return ret;
456 switch (i) {
457 case kwEStrain: return energies_->estrain_;
458 case kwESlip: return energies_->eslip_;
459 }
460 assert(0);
461 return ret;
462 }
463
464 bool ContactModelSmoothJoint::getEnergyAccumulate(uint32 i) const {
465 bool ret(false);
466 if (!energies_) return ret;
467 switch (i) {
468 case kwEStrain: return false;
469 case kwESlip: return true;
470 }
471 assert(0);
472 return ret;
473 }
474
475 void ContactModelSmoothJoint::setEnergy(uint32 i,const double &d) {
476 if (!energies_) return;
477 switch (i) {
478 case kwEStrain: energies_->estrain_ = d; return;
479 case kwESlip: energies_->eslip_ = d; return;
480 }
481 assert(0);
482 return;
483 }
484
485 bool ContactModelSmoothJoint::validate(ContactModelMechanicalState *state,const double &) {
486 assert(state);
487 const IContactMechanical *c = state->getMechanicalContact();
488 assert(c);
489
490 if (state->trackEnergy_)
491 activateEnergy();
492
493 // The radius multiplier has been set previously so calculate the sj_A_
494 // Need to do this regardless of whether or not the radius multiplier has been updated as this is the
495 // first place with the contact state to get the ball radii!
496 updateAreaAndNormal(state);
497
498 if (inheritanceField_ & sjKnMask)
499 updateKn(c);
500 if (inheritanceField_ & sjKsMask)
501 updateKs(c);
502 if (inheritanceField_ & sjFricMask)
503 updateFric(c);
504
505 updateEffectiveTranslationalStiffness();
506
507 return checkActivity(state->gap_);
508 }
509
510 void ContactModelSmoothJoint::updateAreaAndNormal(ContactModelMechanicalState *state) {
511 updateAreaAndNormalContact(state->getContact());
512 }
513
514 void ContactModelSmoothJoint::updateAreaAndNormalContact(const IContact *c) {
515 assert(c);
516 if (geomRecomp_) geomRecomp_ = false;
517 // Use the maximum value of the curvature here - not the min!
518 auto *mc = convert_getcast<IContactMechanical>(c);
519 sj_rad_ = 1.0 / std::max(mc->getEnd1Curvature().y(),mc->getEnd2Curvature().y());
520 sj_rad_ = sj_rmul_ * sj_rad_;
521#ifdef THREED
522 sj_A_ = dPi * sj_rad_ * sj_rad_;
523 sj_unorm_.rx() = sin(dip_) * sin(dd_);
524 sj_unorm_.ry() = sin(dip_) * cos(dd_);
525 sj_unorm_.rz() = cos(dip_);
526#else
527 sj_A_ = 2.0 * sj_rad_; // Assumes thickness of 1 in 2D
528 sj_unorm_.rx() = sin(dip_);
529 sj_unorm_.ry() = cos(dip_);
530#endif
531 DVect nc = toVect(c->getNormal());
532 // Set the flip boolean so that the ordering is correct from side1 to side2
533 isjFlip_ = ( ((nc | sj_unorm_) >= 0.0 ) ? false : true );
534 }
535
536 static const QString kn("sj_kn");
537 bool ContactModelSmoothJoint::updateKn(const IContactMechanical *con) {
538 assert(con);
539 QVariant v1 = con->getEnd1()->getProperty(kn);
540 QVariant v2 = con->getEnd2()->getProperty(kn);
541 if (!v1.isValid() || !v2.isValid())
542 return false;
543 double kn1 = v1.toDouble();
544 double kn2 = v2.toDouble();
545 double val = sj_kn_;
546 if (kn1 && kn2)
547 sj_kn_ = kn1*kn2/(kn1+kn2);
548 else if (kn1)
549 sj_kn_ = kn1;
550 else if (kn2)
551 sj_kn_ = kn2;
552 return ( (sj_kn_ != val) );
553 }
554
555 static const QString ks("sj_ks");
556 bool ContactModelSmoothJoint::updateKs(const IContactMechanical *con) {
557 assert(con);
558 QVariant v1 = con->getEnd1()->getProperty(ks);
559 QVariant v2 = con->getEnd2()->getProperty(ks);
560 if (!v1.isValid() || !v2.isValid())
561 return false;
562 double kn1 = v1.toDouble();
563 double kn2 = v2.toDouble();
564 double val = sj_ks_;
565 if (kn1 && kn2)
566 sj_ks_ = kn1*kn2/(kn1+kn2);
567 else if (kn1)
568 sj_ks_ = kn1;
569 else if (kn2)
570 sj_ks_ = kn2;
571 return ( (sj_ks_ != val) );
572 }
573
574 static const QString fric("sj_fric");
575 bool ContactModelSmoothJoint::updateFric(const IContactMechanical *con) {
576 assert(con);
577 QVariant v1 = con->getEnd1()->getProperty(fric);
578 QVariant v2 = con->getEnd2()->getProperty(fric);
579 if (!v1.isValid() || !v2.isValid())
580 return false;
581 double fric1 = std::max(0.0,v1.toDouble());
582 double fric2 = std::max(0.0,v2.toDouble());
583 double val = sj_fric_;
584 sj_fric_ = std::min(fric1,fric2);
585 return ( (sj_fric_ != val) );
586 }
587
588 bool ContactModelSmoothJoint::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
589 assert(c);
590 QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
591 QRegularExpression rx(name, QRegularExpression::CaseInsensitiveOption);
592 int idx = availableProperties.indexOf(rx) + 1;
593 bool ret=false;
594
595 if (idx<=0)
596 return ret;
597
598 switch(idx) {
599 case kwKn: { //sj_kn_
600 if (inheritanceField_ & sjKnMask)
601 ret = updateKn(c);
602 break;
603 }
604 case kwKs: { //sj_ks_
605 if (inheritanceField_ & sjKsMask)
606 ret =updateKs(c);
607 break;
608 }
609 case kwFric: { //sj_fric_
610 if (inheritanceField_ & sjFricMask)
611 updateFric(c);
612 break;
613 }
614 }
615
616 if (ret)
617 updateEffectiveTranslationalStiffness();
618 return ret;
619 }
620
621 void ContactModelSmoothJoint::updateEffectiveTranslationalStiffness() {
622 effectiveTranslationalStiffness_ = DVect2(sj_kn_,sj_ks_)*sj_A_;
623 }
624
625 bool ContactModelSmoothJoint::forceDisplacementLaw(ContactModelMechanicalState *state,const double &) {
626 assert(state);
627
628 if (geomRecomp_)
629 updateAreaAndNormal(state);
630
631 // Skip out if this should not be active
632 if (sj_large_ && state->gap_ > 0.0 && sj_state_ != 3) {
633 sj_Fn_ = 0.0;
634 sj_Fs_ = DVect(0.0);
635 sj_gap_ = 0.0;
636 sj_Un_ = 0.0;
637 sj_Us_ = DVect(0.0);
638 return false;
639 }
640
641 // Get the translational increment in the global system
642 CAxes localSys = state->getMechanicalContact()->getContact()->getLocalSystem();
643 DVect tInc = localSys.toGlobal(state->relativeTranslationalIncrement_);
644
645 // Now have the translational increment and the normal in the global coordinate system
646 // Compute the normal/shear displacement increments along the sj
647 if (isjFlip_)
648 tInc *= -1.0;
649 double nInc = tInc | sj_unorm_;
650 DVect sInc = tInc - sj_unorm_ * nInc;
651 sj_Un_ -= nInc;
652 sj_Us_ += sInc;
653
654 double nElInc(0.0);
655 DVect sElInc(0.0);
656
657 double g0 = sj_gap_, g1 = sj_gap_ + nInc;
658 sj_gap_ = g1;
659
660 if (!state->canFail_ || sj_state_ == 3 ) { // Bonded
661 nElInc = nInc;
662 sElInc = sInc;
663 } else {
664 if ( g0 <= 0.0 ) {
665 if ( g1 <= 0.0 ) {
666 nElInc = nInc;
667 sElInc = sInc;
668 } else { // g1 > 0.0
669 double xi = -g0 / (g1 - g0);
670 nElInc = nInc * xi;
671 sElInc = sInc * xi;
672 }
673 } else { // g0 > 0.0
674 if ( g1 >= 0.0 ) {
675 nElInc = 0.0;
676 sElInc.fill(0.0);
677 } else { // g1 < 0.0
678 double xi = -g0 / (g1 - g0);
679 nElInc = nInc * (1.0 - xi);
680 sElInc = sInc * (1.0 - xi);
681 }
682 }
683 }
684
685 double del_Fn = -sj_kn_ * sj_A_ * nElInc ;
686 sj_Fn_ += del_Fn ;
687 int slideChange(-1);
688 DVect sj_Fs_old = sj_Fs_;
689
690 if ( state->canFail_ && sj_state_ < 3 ) { // Coulomb sliding with dilation
691 if ( sj_Fn_ <= 0.0 ) {
692 sj_Fn_ = 0.0;
693 sj_Fs_.fill(0.0);
694 } else {
695 DVect del_Fs = sElInc * (-sj_ks_ * sj_A_);
696 sj_Fs_ += del_Fs;
697 double max_Fs = sj_Fn_ * sj_fric_;
698 double magFs = sj_Fs_.mag();
699 if ( magFs > max_Fs ) { // sliding
700 if (!isjSliding_) {
701 isjSliding_ = true;
702 slideChange = 0;
703 }
704 if ( sj_ks_ > 0.0 ) // dilation
705 sj_Fn_ += ( (magFs - max_Fs) / sj_ks_ ) * sj_kn_ * tan(sj_da_);
706 double rat = max_Fs / magFs ;
707 sj_Fs_ *= rat ;
708 } else {
709 if (isjSliding_) {
710 isjSliding_ = false;
711 slideChange = 1;
712 }
713 }
714 }
715 } else { // bonded behavior
716 if ( state->canFail_ && sj_Fn_ <= -sj_bns_ * sj_A_) {
717 sj_state_ = 1;
718 if (cmEvents_[fBondBreak] >= 0) {
719 auto c = state->getContact();
720 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
721 fish::Parameter((qint64)sj_state_),
722 fish::Parameter(sj_bns_) };
723 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
724 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
725 }
726 sj_Fn_ = 0.0;
727 sj_Fs_.fill(0.0);
728 } else {
729 DVect del_Fs = sElInc * (-sj_ks_ * sj_A_);
730 sj_Fs_ += del_Fs;
731 double magFs = sj_Fs_.mag();
732 double ss = calcBSS();
733 if ( state->canFail_ && magFs >= ss * sj_A_) { // break in shear
734 sj_state_ = 2;
735 if (cmEvents_[fBondBreak] >= 0) {
736 auto c = state->getContact();
737 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
738 fish::Parameter((qint64)sj_state_),
739 fish::Parameter(ss) };
740 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
741 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
742 }
743
744 if ( sj_Fn_ < 0.0 ) {
745 sj_Fn_ = 0.0;
746 sj_Fs_.fill(0.0);
747 } else { // was in compression // was in tension
748 double max_Fs = sj_Fn_ * sj_fric_ ;
749 if ( magFs > max_Fs) { // sliding, but no dilation
750 if (!isjSliding_) {
751 isjSliding_ = true;
752 slideChange = 0;
753 }
754 double rat = max_Fs / magFs ;
755 sj_Fs_ *= rat ;
756 } else {
757 if (isjSliding_ == true) {
758 isjSliding_ = false;
759 slideChange = 1;
760 }
761 }
762 }
763 }
764 }
765 }
766 if (slideChange >= 0 && cmEvents_[fSlipChange] >= 0) {
767 auto c = state->getContact();
768 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
769 fish::Parameter((qint64)slideChange) };
770 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
771 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
772 }
773
774 // Have updated the normal and shear forces so need to put them into the contact local
775 // coordinate system and update the forces
776 DVect Fj = sj_unorm_ * sj_Fn_ + sj_Fs_;
777 if (isjFlip_)
778 Fj *= -1.0;
779
780 // Return the correct activity status
781 bool isactive = true;
782 if (sj_large_ && state->gap_ > 0.0 && sj_state_ != 3) {
783 sj_Fn_ = 0.0;
784 sj_Fs_ = DVect(0.0);
785 sj_gap_ = 0.0;
786 sj_Un_ = 0.0;
787 sj_Us_ = DVect(0.0);
788 isactive = false;
789 }
790
791 // update energies
792 if (state->trackEnergy_) {
793 assert(energies_);
794 energies_->estrain_ = 0.0;
795 if (sj_kn_)
796 energies_->estrain_ = 0.5*sj_Fn_*sj_Fn_/sj_kn_;
797 if (sj_ks_) {
798 double smag2 = sj_Fs_.mag2();
799 energies_->estrain_ += 0.5*smag2 / sj_ks_;
800 if (isjSliding_) {
801 DVect avg_F_s = (sj_Fs_old + sj_Fs_)*0.5;
802 DVect u_s_el = (sj_Fs_ - sj_Fs_old) / sj_ks_;
803 energies_->eslip_ -= std::min(0.0,(avg_F_s | (sInc + u_s_el)));
804 }
805 }
806 }
807
808 return isactive ;
809 }
810
811 void ContactModelSmoothJoint::setForce(const DVect &v,IContact *c) {
812 // this is in the local coordinate system
813 CAxes localSys = c->getLocalSystem();
814 DVect globForce = localSys.toGlobal(v);
815 if (isjFlip_)
816 globForce *= -1.0;
817 sj_Fn_ = (sj_unorm_ | globForce);
818 sj_Fs_ = globForce - sj_unorm_ * sj_Fn_;
819 }
820
821 DVect ContactModelSmoothJoint::getForce(const IContactMechanical *c) const {
822 CAxes localSys = c->getContact()->getLocalSystem();
823 DVect Fj = sj_unorm_ * sj_Fn_ + sj_Fs_;
824 if (isjFlip_)
825 Fj *= -1.0;
826 DVect ret(localSys.toLocal(Fj));
827 return ret;
828 }
829
830 DAVect ContactModelSmoothJoint::getMomentOn1(const IContactMechanical *c) const {
831 DVect force = getForce(c);
832 DAVect ret(0.0);
833 c->updateResultingTorqueOn1Local(force,&ret);
834 return ret;
835 }
836
837 DAVect ContactModelSmoothJoint::getMomentOn2(const IContactMechanical *c) const {
838 DVect force = getForce(c);
839 DAVect ret(0.0);
840 c->updateResultingTorqueOn2Local(force,&ret);
841 return ret;
842 }
843
844 double ContactModelSmoothJoint::calcBSS() const {
845 if (sj_A_ > 0) {
846 double dSigma = sj_Fn_ / sj_A_;
847 return dSigma >= -sj_bns_ ? sj_bcoh_ + (tan(sj_bfa_) * dSigma)
848 : sj_bcoh_ + (tan(sj_bfa_) * (-sj_bns_));
849 }
850 else
851 return 0.0;
852 }
853
854
855} // namespace itascaxd
856// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Sep 26, 2024 |
