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 string getName() const { return "smoothjoint"; }
54 virtual void setIndex(int i) { index_=i;}
55 virtual int getIndex() const {return index_;}
56
57 virtual string 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 string getFishCallEvents() const { return "bond_break,slip_change"; }
86 virtual base::Property getProperty(uint32 i,const IContact *con=0) const;
87 virtual bool getPropertyGlobal(uint32 i) const;
88 virtual bool setProperty(uint32 i,const base::Property &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 string getMethods() const {
99 return "sj_setforce,bond,unbond";
100 }
101
102 virtual string getMethodArguments(uint32 i) const;
103 virtual bool setMethod(uint32 i,const std::vector<base::Property> &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 string 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_ = NEW Energies();}
113 virtual bool getEnergyActivated() const {return (energies_ !=0);}
114
115 virtual bool validate(ContactModelMechanicalState *state,const double ×tep);
116 virtual bool endPropertyUpdated(const string &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 NEW 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<string,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 = NEW 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_ = NEW 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_ = NEW 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 base::Property ContactModelSmoothJoint::getProperty(uint32 i,const IContact *) const {
177 // The IContact pointer may be a nullptr!
178 base::Property 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 base::Property();
216 }
217
218 bool ContactModelSmoothJoint::getPropertyGlobal(uint32 ) const {
219 return true;
220 }
221
222 bool ContactModelSmoothJoint::setProperty(uint32 i,const base::Property &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<int64>())
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.to<bool>();
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 string ContactModelSmoothJoint::getMethodArguments(uint32 i) const {
394 string def = string();
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 std::vector<base::Property> &vl,IContact *con) {
405 FP_S;
406 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
407 switch (i) {
408 case kwSetForce: {
409 DVect gf = c->getGlobalForce();
410 if (isjFlip_)
411 gf *= -1.0;
412 sj_Fn_ = gf | sj_unorm_;
413 FP_S;
414 return false;
415 }
416 case kwBond: {
417 if (sj_state_ == 3) return false;
418 double mingap = -1.0 * limits<double>::max();
419 double maxgap = 0;
420 if (vl.at(0).canConvert<double>())
421 maxgap = vl.at(0).toDouble();
422 else if (vl.at(0).canConvert<DVect2>()) {
423 DVect2 value = vl.at(0).value<DVect2>();
424 mingap = value.minComp();
425 maxgap = value.maxComp();
426 } else if (!vl.at(0).isNull())
427 throw Exception("gap value {0} not recognized in method bond in contact model {1}.",vl.at(0),getName());
428 if ( sj_gap_ >= mingap && sj_gap_ <= maxgap) {
429 sj_state_ = 3;
430 FP_S;
431 return true;
432 }
433 FP_S;
434 return false;
435 }
436 case kwUnbond: {
437 if (sj_state_ == 0) return false;
438 double mingap = -1.0 * limits<double>::max();
439 double maxgap = 0;
440 if (vl.at(0).canConvert<double>())
441 maxgap = vl.at(0).toDouble();
442 else if (vl.at(0).canConvert<DVect2>()) {
443 DVect2 value = vl.at(0).value<DVect2>();
444 mingap = value.minComp();
445 maxgap = value.maxComp();
446 }
447 else if (!vl.at(0).isNull())
448 throw Exception("gap value {0} not recognized in method unbond in contact model {1}.",vl.at(0),getName());
449 if ( sj_gap_ >= mingap && sj_gap_ <= maxgap) {
450 sj_state_ = 0;
451 FP_S;
452 return true;
453 }
454 FP_S;
455 return false;
456 }
457
458 }
459 FP_S;
460 return false;
461 }
462
463 double ContactModelSmoothJoint::getEnergy(uint32 i) const {
464 double ret(0.0);
465 if (!energies_)
466 return ret;
467 switch (i) {
468 case kwEStrain: return energies_->estrain_;
469 case kwESlip: return energies_->eslip_;
470 }
471 assert(0);
472 return ret;
473 }
474
475 bool ContactModelSmoothJoint::getEnergyAccumulate(uint32 i) const {
476 bool ret(false);
477 if (!energies_) return ret;
478 switch (i) {
479 case kwEStrain: return false;
480 case kwESlip: return true;
481 }
482 assert(0);
483 return ret;
484 }
485
486 void ContactModelSmoothJoint::setEnergy(uint32 i,const double &d) {
487 if (!energies_) return;
488 switch (i) {
489 case kwEStrain: energies_->estrain_ = d; return;
490 case kwESlip: energies_->eslip_ = d; return;
491 }
492 assert(0);
493 return;
494 }
495
496 bool ContactModelSmoothJoint::validate(ContactModelMechanicalState *state,const double &) {
497 assert(state);
498 const IContactMechanical *c = state->getMechanicalContact();
499 assert(c);
500 localSystem_ = c->getContact()->getLocalSystem();
501
502 if (state->trackEnergy_)
503 activateEnergy();
504
505 // The radius multiplier has been set previously so calculate the sj_A_
506 // Need to do this regardless of whether or not the radius multiplier has been updated as this is the
507 // first place with the contact state to get the ball radii!
508 updateAreaAndNormal(state);
509
510 if (inheritanceField_ & sjKnMask)
511 updateKn(c);
512 if (inheritanceField_ & sjKsMask)
513 updateKs(c);
514 if (inheritanceField_ & sjFricMask)
515 updateFric(c);
516
517 updateEffectiveTranslationalStiffness();
518
519 return checkActivity(state->gap_);
520 }
521
522 void ContactModelSmoothJoint::updateAreaAndNormal(ContactModelMechanicalState *state) {
523 updateAreaAndNormalContact(state->end1Curvature_,state->end2Curvature_,state->axes_.e1());
524 }
525
526 void ContactModelSmoothJoint::updateAreaAndNormalContact(const DVect2 &end1Curvature,const DVect2 &end2Curvature,const DVect &norm) {
527 if (geomRecomp_) geomRecomp_ = false;
528 // Use the maximum value of the curvature here - not the min!
529 sj_rad_ = 1.0 / std::max(end1Curvature.y(),end2Curvature.y());
530 sj_rad_ = sj_rmul_ * sj_rad_;
531#ifdef THREED
532 sj_A_ = dPi * sj_rad_ * sj_rad_;
533 sj_unorm_.rx() = sin(dip_) * sin(dd_);
534 sj_unorm_.ry() = sin(dip_) * cos(dd_);
535 sj_unorm_.rz() = cos(dip_);
536#else
537 sj_A_ = 2.0 * sj_rad_; // Assumes thickness of 1 in 2D
538 sj_unorm_.rx() = sin(dip_);
539 sj_unorm_.ry() = cos(dip_);
540#endif
541 DVect nc = norm;
542 // Set the flip boolean so that the ordering is correct from side1 to side2
543 isjFlip_ = ( ((nc | sj_unorm_) >= 0.0 ) ? false : true );
544 }
545
546 static const string kn("sj_kn");
547 bool ContactModelSmoothJoint::updateKn(const IContactMechanical *con) {
548 assert(con);
549 base::Property v1 = con->getEnd1()->getProperty(kn);
550 base::Property v2 = con->getEnd2()->getProperty(kn);
551 if (!v1.isValid() || !v2.isValid())
552 return false;
553 double kn1 = v1.toDouble();
554 double kn2 = v2.toDouble();
555 double val = sj_kn_;
556 if (kn1 && kn2)
557 sj_kn_ = kn1*kn2/(kn1+kn2);
558 else if (kn1)
559 sj_kn_ = kn1;
560 else if (kn2)
561 sj_kn_ = kn2;
562 return ( (sj_kn_ != val) );
563 }
564
565 static const string ks("sj_ks");
566 bool ContactModelSmoothJoint::updateKs(const IContactMechanical *con) {
567 assert(con);
568 base::Property v1 = con->getEnd1()->getProperty(ks);
569 base::Property v2 = con->getEnd2()->getProperty(ks);
570 if (!v1.isValid() || !v2.isValid())
571 return false;
572 double kn1 = v1.toDouble();
573 double kn2 = v2.toDouble();
574 double val = sj_ks_;
575 if (kn1 && kn2)
576 sj_ks_ = kn1*kn2/(kn1+kn2);
577 else if (kn1)
578 sj_ks_ = kn1;
579 else if (kn2)
580 sj_ks_ = kn2;
581 return ( (sj_ks_ != val) );
582 }
583
584 static const string fric("sj_fric");
585 bool ContactModelSmoothJoint::updateFric(const IContactMechanical *con) {
586 assert(con);
587 base::Property v1 = con->getEnd1()->getProperty(fric);
588 base::Property v2 = con->getEnd2()->getProperty(fric);
589 if (!v1.isValid() || !v2.isValid())
590 return false;
591 double fric1 = std::max(0.0,v1.toDouble());
592 double fric2 = std::max(0.0,v2.toDouble());
593 double val = sj_fric_;
594 sj_fric_ = std::min(fric1,fric2);
595 return ( (sj_fric_ != val) );
596 }
597
598 bool ContactModelSmoothJoint::endPropertyUpdated(const string &name,const IContactMechanical *c) {
599 assert(c);
600 StringList availableProperties = split(replace(simplified(getProperties())," ",""),",");
601 auto idx = findRegex(availableProperties,name);
602 if (idx==string::npos) return false;
603 idx += 1;
604 bool ret=false;
605 switch(idx) {
606 case kwKn: { //sj_kn_
607 if (inheritanceField_ & sjKnMask)
608 ret = updateKn(c);
609 break;
610 }
611 case kwKs: { //sj_ks_
612 if (inheritanceField_ & sjKsMask)
613 ret =updateKs(c);
614 break;
615 }
616 case kwFric: { //sj_fric_
617 if (inheritanceField_ & sjFricMask)
618 updateFric(c);
619 break;
620 }
621 }
622
623 if (ret)
624 updateEffectiveTranslationalStiffness();
625 return ret;
626 }
627
628 void ContactModelSmoothJoint::updateEffectiveTranslationalStiffness() {
629 effectiveTranslationalStiffness_ = DVect2(sj_kn_,sj_ks_)*sj_A_;
630 }
631
632 bool ContactModelSmoothJoint::forceDisplacementLaw(ContactModelMechanicalState *state,const double &) {
633 assert(state);
634
635 if (geomRecomp_)
636 updateAreaAndNormal(state);
637
638 // Skip out if this should not be active
639 if (sj_large_ && state->gap_ > 0.0 && sj_state_ != 3) {
640 sj_Fn_ = 0.0;
641 sj_Fs_ = DVect(0.0);
642 sj_gap_ = 0.0;
643 sj_Un_ = 0.0;
644 sj_Us_ = DVect(0.0);
645 return false;
646 }
647
648 // Get the translational increment in the global system
649 localSystem_ = state->axes_;
650 DVect tInc = localSystem_.toGlobal(state->relativeTranslationalIncrement_);
651
652 // Now have the translational increment and the normal in the global coordinate system
653 // Compute the normal/shear displacement increments along the sj
654 if (isjFlip_)
655 tInc *= -1.0;
656 double nInc = tInc | sj_unorm_;
657 DVect sInc = tInc - sj_unorm_ * nInc;
658 sj_Un_ -= nInc;
659 sj_Us_ += sInc;
660
661 double nElInc(0.0);
662 DVect sElInc(0.0);
663
664 double g0 = sj_gap_, g1 = sj_gap_ + nInc;
665 sj_gap_ = g1;
666
667 if (!state->canFail_ || sj_state_ == 3 ) { // Bonded
668 nElInc = nInc;
669 sElInc = sInc;
670 } else {
671 if ( g0 <= 0.0 ) {
672 if ( g1 <= 0.0 ) {
673 nElInc = nInc;
674 sElInc = sInc;
675 } else { // g1 > 0.0
676 double xi = -g0 / (g1 - g0);
677 nElInc = nInc * xi;
678 sElInc = sInc * xi;
679 }
680 } else { // g0 > 0.0
681 if ( g1 >= 0.0 ) {
682 nElInc = 0.0;
683 sElInc.fill(0.0);
684 } else { // g1 < 0.0
685 double xi = -g0 / (g1 - g0);
686 nElInc = nInc * (1.0 - xi);
687 sElInc = sInc * (1.0 - xi);
688 }
689 }
690 }
691
692 double del_Fn = -sj_kn_ * sj_A_ * nElInc ;
693 sj_Fn_ += del_Fn ;
694 int slideChange(-1);
695 DVect sj_Fs_old = sj_Fs_;
696
697 if ( state->canFail_ && sj_state_ < 3 ) { // Coulomb sliding with dilation
698 if ( sj_Fn_ <= 0.0 ) {
699 sj_Fn_ = 0.0;
700 sj_Fs_.fill(0.0);
701 } else {
702 DVect del_Fs = sElInc * (-sj_ks_ * sj_A_);
703 sj_Fs_ += del_Fs;
704 double max_Fs = sj_Fn_ * sj_fric_;
705 double magFs = sj_Fs_.mag();
706 if ( magFs > max_Fs ) { // sliding
707 if (!isjSliding_) {
708 isjSliding_ = true;
709 slideChange = 0;
710 }
711 if ( sj_ks_ > 0.0 ) // dilation
712 sj_Fn_ += ( (magFs - max_Fs) / sj_ks_ ) * sj_kn_ * tan(sj_da_);
713 double rat = max_Fs / magFs ;
714 sj_Fs_ *= rat ;
715 } else {
716 if (isjSliding_) {
717 isjSliding_ = false;
718 slideChange = 1;
719 }
720 }
721 }
722 } else { // bonded behavior
723 if ( state->canFail_ && sj_Fn_ <= -sj_bns_ * sj_A_) {
724 sj_state_ = 1;
725 if (cmEvents_[fBondBreak] >= 0) {
726 auto c = state->getContact();
727 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
728 fish::Parameter((int64)sj_state_),
729 fish::Parameter(sj_bns_) };
730 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
731 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
732 }
733 sj_Fn_ = 0.0;
734 sj_Fs_.fill(0.0);
735 } else {
736 DVect del_Fs = sElInc * (-sj_ks_ * sj_A_);
737 sj_Fs_ += del_Fs;
738 double magFs = sj_Fs_.mag();
739 double ss = calcBSS();
740 if ( state->canFail_ && magFs >= ss * sj_A_) { // break in shear
741 sj_state_ = 2;
742 if (cmEvents_[fBondBreak] >= 0) {
743 auto c = state->getContact();
744 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
745 fish::Parameter((int64)sj_state_),
746 fish::Parameter(ss) };
747 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
748 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
749 }
750
751 if ( sj_Fn_ < 0.0 ) {
752 sj_Fn_ = 0.0;
753 sj_Fs_.fill(0.0);
754 } else { // was in compression // was in tension
755 double max_Fs = sj_Fn_ * sj_fric_ ;
756 if ( magFs > max_Fs) { // sliding, but no dilation
757 if (!isjSliding_) {
758 isjSliding_ = true;
759 slideChange = 0;
760 }
761 double rat = max_Fs / magFs ;
762 sj_Fs_ *= rat ;
763 } else {
764 if (isjSliding_ == true) {
765 isjSliding_ = false;
766 slideChange = 1;
767 }
768 }
769 }
770 }
771 }
772 }
773 if (slideChange >= 0 && cmEvents_[fSlipChange] >= 0) {
774 auto c = state->getContact();
775 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
776 fish::Parameter((int64)slideChange) };
777 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
778 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
779 }
780
781 // Have updated the normal and shear forces so need to put them into the contact local
782 // coordinate system and update the forces
783 DVect Fj = sj_unorm_ * sj_Fn_ + sj_Fs_;
784 if (isjFlip_)
785 Fj *= -1.0;
786
787 // Return the correct activity status
788 bool isactive = true;
789 if (sj_large_ && state->gap_ > 0.0 && sj_state_ != 3) {
790 sj_Fn_ = 0.0;
791 sj_Fs_ = DVect(0.0);
792 sj_gap_ = 0.0;
793 sj_Un_ = 0.0;
794 sj_Us_ = DVect(0.0);
795 isactive = false;
796 }
797
798 // update energies
799 if (state->trackEnergy_) {
800 assert(energies_);
801 energies_->estrain_ = 0.0;
802 if (sj_kn_)
803 energies_->estrain_ = 0.5*sj_Fn_*sj_Fn_/sj_kn_;
804 if (sj_ks_) {
805 double smag2 = sj_Fs_.mag2();
806 energies_->estrain_ += 0.5*smag2 / sj_ks_;
807 if (isjSliding_) {
808 DVect avg_F_s = (sj_Fs_old + sj_Fs_)*0.5;
809 DVect u_s_el = (sj_Fs_ - sj_Fs_old) / sj_ks_;
810 energies_->eslip_ -= std::min(0.0,(avg_F_s | (sInc + u_s_el)));
811 }
812 }
813 }
814
815 return isactive ;
816 }
817
818 void ContactModelSmoothJoint::setForce(const DVect &v,IContact *c) {
819 // this is in the local coordinate system
820 localSystem_ = c->getLocalSystem();
821 DVect globForce = localSystem_.toGlobal(v);
822 if (isjFlip_)
823 globForce *= -1.0;
824 sj_Fn_ = (sj_unorm_ | globForce);
825 sj_Fs_ = globForce - sj_unorm_ * sj_Fn_;
826 }
827
828 DVect ContactModelSmoothJoint::getForce() const {
829 DVect Fj = sj_unorm_ * sj_Fn_ + sj_Fs_;
830 if (isjFlip_)
831 Fj *= -1.0;
832 DVect ret(localSystem_.toLocal(Fj));
833 return ret;
834 }
835
836 DAVect ContactModelSmoothJoint::getMomentOn1(const IContactMechanical *c) const {
837 DVect force = getForce();
838 DAVect ret(0.0);
839 c->updateResultingTorqueOn1Local(force,&ret);
840 return ret;
841 }
842
843 DAVect ContactModelSmoothJoint::getMomentOn2(const IContactMechanical *c) const {
844 DVect force = getForce();
845 DAVect ret(0.0);
846 c->updateResultingTorqueOn2Local(force,&ret);
847 return ret;
848 }
849
850 DAVect ContactModelSmoothJoint::getModelMomentOn1() const {
851 DAVect ret(0.0);
852 return ret;
853 }
854
855 DAVect ContactModelSmoothJoint::getModelMomentOn2() const {
856 DAVect ret(0.0);
857 return ret;
858 }
859
860 void ContactModelSmoothJoint::objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *ret) const {
861 ret->clear();
862 ret->push_back({"sj_kn",ScalarInfo});
863 ret->push_back({"sj_ks",ScalarInfo});
864 ret->push_back({"sj_fric",ScalarInfo});
865 ret->push_back({"sj_da",ScalarInfo});
866 ret->push_back({"sj_state",ScalarInfo});
867 ret->push_back({"sj_ten",ScalarInfo});
868 ret->push_back({"sj_coh",ScalarInfo});
869 ret->push_back({"sj_fa",ScalarInfo});
870 ret->push_back({"sj_shear",ScalarInfo});
871 ret->push_back({"sj_rmul",ScalarInfo});
872 ret->push_back({"sj_large",ScalarInfo});
873 ret->push_back({"sj_fn",ScalarInfo});
874 ret->push_back({"sj_fs",VectorInfo});
875 ret->push_back({"sj_gap",ScalarInfo});
876 ret->push_back({"sj_unorm",VectorInfo});
877 ret->push_back({"sj_area",ScalarInfo});
878 ret->push_back({"sj_radius",ScalarInfo});
879 ret->push_back({"sj_un",ScalarInfo});
880 ret->push_back({"sj_us",VectorInfo});
881 ret->push_back({"sj_slip",ScalarInfo});
882 ret->push_back({"dip",ScalarInfo});
883#ifdef THREED
884 ret->push_back({"ddir",ScalarInfo});
885#endif
886
887 }
888
889 void ContactModelSmoothJoint::objectPropValues(std::vector<double> *ret,const IContact *) const {
890 FP_S;
891 ret->clear();
892 ret->push_back(sj_kn_);
893 ret->push_back(sj_ks_);
894 ret->push_back(sj_fric_);
895 ret->push_back(sj_da_/dDegrad);
896 ret->push_back(sj_state_);
897 ret->push_back(sj_bns_);
898 ret->push_back(sj_bcoh_);
899 ret->push_back(sj_bfa_/dDegrad);
900 ret->push_back(calcBSS());
901 ret->push_back(sj_rmul_);
902 ret->push_back(sj_large_);
903 ret->push_back(sj_Fn_);
904 ret->push_back(sj_Fs_.mag());
905 ret->push_back(sj_gap_);
906 ret->push_back(sj_unorm_.mag());
907 ret->push_back(sj_A_);
908 ret->push_back(sj_rad_);
909 ret->push_back(sj_Un_);
910 ret->push_back(sj_Us_.mag());
911 ret->push_back(isjSliding_);
912 ret->push_back(dip_/dDegrad);
913#ifdef THREED
914 ret->push_back(dd_/dDegrad);
915#endif
916 FP_S;
917 }
918
919 double ContactModelSmoothJoint::calcBSS() const {
920 if (sj_A_ > 0) {
921 double dSigma = sj_Fn_ / sj_A_;
922 return dSigma >= -sj_bns_ ? sj_bcoh_ + (tan(sj_bfa_) * dSigma)
923 : sj_bcoh_ + (tan(sj_bfa_) * (-sj_bns_));
924 }
925 else
926 return 0.0;
927 }
928
929
930} // namespace itascaxd
931// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Jun 07, 2025 |
