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