Linear Contact Bond Model Implementation
See this file for the documentation of this contact model.
contactmodellinearcbond.h
1#pragma once
2// contactmodellinearcbond.h
3
4#include "contactmodel/src/contactmodelmechanical.h"
5
6#ifdef LINEARCBOND_LIB
7# define LINEARCBOND_EXPORT EXPORT_TAG
8#elif defined(NO_MODEL_IMPORT)
9# define LINEARCBOND_EXPORT
10#else
11# define LINEARCBOND_EXPORT IMPORT_TAG
12#endif
13
14namespace cmodelsxd {
15 using namespace itasca;
16
17 class ContactModelLinearCBond : public ContactModelMechanical {
18 public:
19 LINEARCBOND_EXPORT ContactModelLinearCBond();
20 LINEARCBOND_EXPORT virtual ~ContactModelLinearCBond();
21 void copy(const ContactModel *c) override;
22 void archive(ArchiveStream &) override;
23
24 QString getName() const override { return "linearcbond"; }
25 void setIndex(int i) override { index_=i;}
26 int getIndex() const override {return index_;}
27
28 enum PropertyKeys {
29 kwKn=1
30 , kwKs
31 , kwFric
32 , kwLinF
33 , kwLinS
34 , kwLinMode
35 , kwRGap
36 , kwEmod
37 , kwKRatio
38 , kwDpNRatio
39 , kwDpSRatio
40 , kwDpMode
41 , kwDpF
42 , kwCbState
43 , kwCbTenF
44 , kwCbShearF
45 , kwCbTStr
46 , kwCbSStr
47 , kwUserArea
48 };
49
50 QString getProperties() const override {
51 return "kn"
52 ",ks"
53 ",fric"
54 ",lin_force"
55 ",lin_slip"
56 ",lin_mode"
57 ",rgap"
58 ",emod"
59 ",kratio"
60 ",dp_nratio"
61 ",dp_sratio"
62 ",dp_mode"
63 ",dp_force"
64 ",cb_state"
65 ",cb_tenf"
66 ",cb_shearf"
67 ",cb_tens"
68 ",cb_shears"
69 ",user_area";
70 }
71
72 enum EnergyKeys { kwEStrain=1,kwESlip,kwEDashpot};
73 QString getEnergies() const override { return "energy-strain,energy-slip,energy-dashpot";}
74 double getEnergy(uint32 i) const override; // Base 1
75 bool getEnergyAccumulate(uint32 i) const override; // Base 1
76 void setEnergy(uint32 i,const double &d) override; // Base 1
77 void activateEnergy() override { if (energies_) return; energies_ = NEWC(Energies());}
78 bool getEnergyActivated() const override {return (energies_ !=0);}
79
80 enum FishCallEvents {fActivated=0,fBondBreak, fSlipChange };
81 QString getFishCallEvents() const override { return "contact_activated,bond_break,slip_change"; }
82 QVariant getProperty(uint32 i,const IContact *) const override;
83 bool getPropertyGlobal(uint32 i) const override;
84 bool setProperty(uint32 i,const QVariant &v,IContact *) override;
85 bool getPropertyReadOnly(uint32 i) const override;
86
87 bool supportsInheritance(uint32 i) const override;
88 bool getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i); return (inheritanceField_ & mask) ? true : false; }
89 void setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i); if (b) inheritanceField_ |= mask; else inheritanceField_ &= ~mask; }
90
91 enum MethodKeys {
92 kwDeformability=1
93 , kwCbBond
94 , kwCbStrength
95 , kwCbUnbond
96 , kwArea
97 };
98
99 QString getMethods() const override {
100 return "deformability"
101 ",bond"
102 ",cb_strength"
103 ",unbond"
104 ",area";
105 }
106
107 QString getMethodArguments(uint32 i) const override;
108
109 bool setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override; // Base 1 - returns true if timestep contributions need to be updated
110
111 uint32 getMinorVersion() const override;
112
113 bool validate(ContactModelMechanicalState *state,const double ×tep) override;
114 bool endPropertyUpdated(const QString &name,const IContactMechanical *c) override;
115 bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) override;
116 DVect2 getEffectiveTranslationalStiffness() const override { DVect2 ret = effectiveTranslationalStiffness_; return ret;}
117 bool thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
118
119 DAVect getEffectiveRotationalStiffness() const override { return DAVect(0.0);}
120
121 ContactModelLinearCBond *clone() const override { return NEWC(ContactModelLinearCBond()); }
122 double getActivityDistance() const override {return rgap_;}
123 bool isOKToDelete() const override { return !isBonded(); }
124 void resetForcesAndMoments() override { lin_F(DVect(0.0)); dp_F(DVect(0.0)); if (energies_) energies_->estrain_ = 0.0; }
125 void setForce(const DVect &v,IContact *c) override;
126 void setArea(const double &d) override { userArea_ = d; }
127 double getArea() const override { return userArea_; }
128
129 bool checkActivity(const double &gap) override { return (gap <= rgap_ || isBonded()); }
130
131 bool isSliding() const override { return lin_S_; }
132 bool isBonded() const override { return (cb_state_==3); }
133 void unbond() override { cb_state_ = 0; }
134 void propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
135 void setNonForcePropsFrom(IContactModel *oldCM) override;
136 /// Return the total force that the contact model holds.
137 DVect getForce(const IContactMechanical*) const override;
138 /// Return the total moment on 1 that the contact model holds
139 DAVect getMomentOn1(const IContactMechanical*) const override;
140 /// Return the total moment on 1 that the contact model holds
141 DAVect getMomentOn2(const IContactMechanical*) const override;
142
143 const double & kn() const {return kn_;}
144 void kn(const double &d) {kn_=d;}
145 const double & ks() const {return ks_;}
146 void ks(const double &d) {ks_=d;}
147 const double & fric() const {return fric_;}
148 void fric(const double &d) {fric_=d;}
149 const DVect & lin_F() const {return lin_F_;}
150 void lin_F(const DVect &f) { lin_F_=f;}
151 bool lin_S() const {return lin_S_;}
152 void lin_S(bool b) { lin_S_=b;}
153 uint32 lin_mode() const {return lin_mode_;}
154 void lin_mode(uint32 i) { lin_mode_=i;}
155 const double & rgap() const {return rgap_;}
156 void rgap(const double &d) {rgap_=d;}
157 uint32 cb_state() const {return cb_state_;}
158 void cb_state(uint32 b) { cb_state_=b;}
159 const double & cb_tenF() const {return cb_tenF_;}
160 void cb_tenF(const double &d) {cb_tenF_=d;}
161 const double & cb_shearF() const {return cb_shearF_;}
162 void cb_shearF(const double &d) {cb_shearF_=d;}
163
164 bool hasDamping() const {return dpProps_ ? true : false;}
165 double dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
166 void dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
167 double dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
168 void dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
169 int dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
170 void dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
171 DVect dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
172 void dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
173
174 bool hasEnergies() const {return energies_ ? true:false;}
175 double estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
176 void estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
177 double eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
178 void eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
179 double edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
180 void edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
181
182 uint32 inheritanceField() const {return inheritanceField_;}
183 void inheritanceField(uint32 i) {inheritanceField_ = i;}
184
185 const DVect2 & effectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
186 void effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
187
188 private:
189 static int index_;
190
191 struct Energies {
192 Energies() : estrain_(0.0), eslip_(0.0),edashpot_(0.0) {}
193 double estrain_; // elastic energy stored in contact
194 double eslip_; // work dissipated by friction
195 double edashpot_; // work dissipated by dashpots
196 };
197
198 struct dpProps {
199 dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)) {}
200 double dp_nratio_; // normal viscous critical damping ratio
201 double dp_sratio_; // shear viscous critical damping ratio
202 int dp_mode_; // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
203 DVect dp_F_; // Force in the dashpots
204 };
205
206 bool updateKn(const IContactMechanical *con);
207 bool updateKs(const IContactMechanical *con);
208 bool updateFric(const IContactMechanical *con);
209
210 void updateEffectiveStiffness(ContactModelMechanicalState *state);
211
212 void setDampCoefficients(const double &mass,double *vcn,double *vcs);
213
214 // inheritance fields
215 uint32 inheritanceField_;
216
217 // linear model
218 double kn_; // normal stiffness
219 double ks_; // shear stiffness
220 double fric_; // Coulomb friction coefficient
221 DVect lin_F_; // Force carried in the linear model
222 bool lin_S_; // current slip state
223 uint32 lin_mode_; // specifies incremental or absolute for the the linear part
224 double rgap_; // reference gap
225
226 uint32 cb_state_; // Bond state - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (B)
227 double cb_tenF_;
228 double cb_shearF_;
229
230 dpProps * dpProps_; // The viscous properties
231
232 double userArea_; // Area as specified by the user
233
234 Energies * energies_; // energies
235
236 DVect2 effectiveTranslationalStiffness_;
237
238 };
239} // namespace itascaxd
240// EoF
contactmodellinearcbond.cpp
1// contactmodellinearcbond.cpp
2#include "contactmodellinearcbond.h"
3
4#include "../version.txt"
5#include "contactmodel/src/contactmodelthermal.h"
6#include "fish/src/parameter.h"
7#include "utility/src/tptr.h"
8#include "shared/src/mathutil.h"
9#include "base/src/basetoqt.h"
10
11#include "kernel/interface/iprogram.h"
12#include "module/interface/icontact.h"
13#include "module/interface/icontactmechanical.h"
14#include "module/interface/icontactthermal.h"
15#include "module/interface/ifishcalllist.h"
16#include "module/interface/ipiece.h"
17#include "module/interface/ipiecemechanical.h"
18
19
20
21
22#ifdef LINEARCBOND_LIB
23#ifdef _WIN32
24 int __stdcall DllMain(void *,unsigned, void *)
25 {
26 return 1;
27 }
28#endif
29
30 extern "C" EXPORT_TAG const char *getName()
31 {
32#if DIM==3
33 return "contactmodelmechanical3dlinearcbond";
34#else
35 return "contactmodelmechanical2dlinearcbond";
36#endif
37 }
38
39 extern "C" EXPORT_TAG unsigned getMajorVersion()
40 {
41 return MAJOR_VRESION;
42 }
43
44 extern "C" EXPORT_TAG unsigned getMinorVersion()
45 {
46 return MINOR_VERSION;
47 }
48
49 extern "C" EXPORT_TAG void *createInstance()
50 {
51 cmodelsxd::ContactModelLinearCBond *m = NEWC(cmodelsxd::ContactModelLinearCBond());
52 return (void *)m;
53 }
54#endif // LINEARCBOND_EXPORTS
55
56namespace cmodelsxd {
57 static const uint32 linKnMask = 0x00002; // Base 1!
58 static const uint32 linKsMask = 0x00004;
59 static const uint32 linFricMask = 0x00008;
60
61 using namespace itasca;
62
63 int ContactModelLinearCBond::index_ = -1;
64 uint32 ContactModelLinearCBond::getMinorVersion() const { return MINOR_VERSION;}
65
66 ContactModelLinearCBond::ContactModelLinearCBond() : inheritanceField_(linKnMask|linKsMask|linFricMask)
67 , kn_(0.0)
68 , ks_(0.0)
69 , fric_(0.0)
70 , lin_F_(DVect(0.0))
71 , lin_S_(false)
72 , lin_mode_(0)
73 , rgap_(0.0)
74 , cb_state_(0)
75 , cb_tenF_(0.0)
76 , cb_shearF_(0.0)
77 , dpProps_(0)
78 , userArea_(0)
79 , energies_(0)
80 , effectiveTranslationalStiffness_(DVect2(0.0)) {
81// setFromParent(ContactModelMechanicalList::instance()->find(getName()));
82 }
83
84 ContactModelLinearCBond::~ContactModelLinearCBond() {
85 if (dpProps_)
86 delete dpProps_;
87 if (energies_)
88 delete energies_;
89 }
90
91 void ContactModelLinearCBond::archive(ArchiveStream &stream) {
92 stream & kn_;
93 stream & ks_;
94 stream & fric_;
95 stream & lin_F_;
96 stream & lin_S_;
97 stream & lin_mode_;
98 stream & cb_state_;
99 stream & cb_tenF_;
100 stream & cb_shearF_;
101
102 if (stream.getArchiveState()==ArchiveStream::Save) {
103 bool b = false;
104 if (dpProps_) {
105 b = true;
106 stream & b;
107 stream & dpProps_->dp_nratio_;
108 stream & dpProps_->dp_sratio_;
109 stream & dpProps_->dp_mode_;
110 stream & dpProps_->dp_F_;
111 }
112 else
113 stream & b;
114
115 b = false;
116 if (energies_) {
117 b = true;
118 stream & b;
119 stream & energies_->estrain_;
120 stream & energies_->eslip_;
121 stream & energies_->edashpot_;
122 }
123 else
124 stream & b;
125 } else {
126 bool b(false);
127 stream & b;
128 if (b) {
129 if (!dpProps_)
130 dpProps_ = NEWC(dpProps());
131 stream & dpProps_->dp_nratio_;
132 stream & dpProps_->dp_sratio_;
133 stream & dpProps_->dp_mode_;
134 stream & dpProps_->dp_F_;
135 }
136 stream & b;
137 if (b) {
138 if (!energies_)
139 energies_ = NEWC(Energies());
140 stream & energies_->estrain_;
141 stream & energies_->eslip_;
142 stream & energies_->edashpot_;
143 }
144 }
145
146 stream & inheritanceField_;
147 stream & effectiveTranslationalStiffness_;
148
149 if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() == getMinorVersion())
150 stream & rgap_;
151
152 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 2)
153 stream & userArea_;
154 }
155
156 void ContactModelLinearCBond::copy(const ContactModel *cm) {
157 ContactModelMechanical::copy(cm);
158 const ContactModelLinearCBond *in = dynamic_cast<const ContactModelLinearCBond*>(cm);
159 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
160 kn(in->kn());
161 ks(in->ks());
162 fric(in->fric());
163 lin_F(in->lin_F());
164 lin_S(in->lin_S());
165 lin_mode(in->lin_mode());
166 rgap(in->rgap());
167 cb_state(in->cb_state());
168 cb_tenF(in->cb_tenF());
169 cb_shearF(in->cb_shearF());
170 if (in->hasDamping()) {
171 if (!dpProps_)
172 dpProps_ = NEWC(dpProps());
173 dp_nratio(in->dp_nratio());
174 dp_sratio(in->dp_sratio());
175 dp_mode(in->dp_mode());
176 dp_F(in->dp_F());
177 }
178 if (in->hasEnergies()) {
179 if (!energies_)
180 energies_ = NEWC(Energies());
181 estrain(in->estrain());
182 eslip(in->eslip());
183 edashpot(in->edashpot());
184 }
185 userArea_ = in->userArea_;
186 inheritanceField(in->inheritanceField());
187 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
188 }
189
190 QVariant ContactModelLinearCBond::getProperty(uint32 i,const IContact *con) const {
191 QVariant var;
192 bool nstr = false;
193 switch (i) {
194 case kwKn: return kn_;
195 case kwKs: return ks_;
196 case kwFric: return fric_;
197 case kwLinF: var.setValue(lin_F_); return var;
198 case kwLinS: return lin_S_;
199 case kwLinMode: return lin_mode_;
200 case kwRGap: return rgap_;
201 case kwEmod: {
202 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
203 if (c ==nullptr) return 0.0;
204 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
205 double rsum(0.0);
206 if (c->getEnd1Curvature().y())
207 rsum += 1.0/c->getEnd1Curvature().y();
208 if (c->getEnd2Curvature().y())
209 rsum += 1.0/c->getEnd2Curvature().y();
210 if (userArea_) {
211#ifdef THREED
212 rsq = std::sqrt(userArea_ / dPi);
213#else
214 rsq = userArea_ / 2.0;
215#endif
216 rsum = rsq + rsq;
217 rsq = 1. / rsq;
218 }
219#ifdef TWOD
220 return (kn_ * rsum * rsq / 2.0);
221#else
222 return (kn_ * rsum * rsq * rsq) / dPi;
223#endif
224 }
225 case kwKRatio: return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
226 case kwDpNRatio: return dpProps_ ? dpProps_->dp_nratio_ : 0;
227 case kwDpSRatio: return dpProps_ ? dpProps_->dp_sratio_ : 0;
228 case kwDpMode: return dpProps_ ? dpProps_->dp_mode_ : 0;
229 case kwDpF: {
230 dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
231 return var;
232 }
233 case kwCbState: return cb_state_;
234 case kwCbTenF: return cb_tenF_;
235 case kwCbShearF: return cb_shearF_;
236 case kwCbTStr: nstr = true;
237 [[fallthrough]];
238 case kwCbSStr: {
239 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
240 if (c ==nullptr) return 0.0;
241 double tmp(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
242 if (userArea_) {
243#ifdef THREED
244 tmp = std::sqrt(userArea_ / dPi);
245#else
246 tmp = userArea_ / 2.0;
247#endif
248 tmp = 1. / tmp;
249 }
250 if (nstr) {
251#ifdef TWOD
252 return (cb_tenF_ * tmp / 2.0);
253#else
254 return (cb_tenF_ * tmp * tmp / dPi);
255#endif
256 } else {
257#ifdef TWOD
258 return (cb_shearF_ * tmp / 2.0);
259#else
260 return (cb_shearF_ * tmp * tmp / dPi);
261#endif
262 }
263 }
264 case kwUserArea: return userArea_;
265 }
266 assert(0);
267 return QVariant();
268 }
269
270 bool ContactModelLinearCBond::getPropertyGlobal(uint32 i) const {
271 switch (i) {
272 case kwLinF:
273 case kwDpF:
274 return false;
275 }
276 return true;
277 }
278
279 bool ContactModelLinearCBond::setProperty(uint32 i,const QVariant &v,IContact *) {
280 dpProps dp;
281 switch (i) {
282 case kwKn: {
283 if (!v.canConvert<double>())
284 throw Exception("kn must be a double.");
285 double val(v.toDouble());
286 if (val<0.0)
287 throw Exception("Negative kn not allowed.");
288 kn_ = val;
289 return true;
290 }
291 case kwKs: {
292 if (!v.canConvert<double>())
293 throw Exception("ks must be a double.");
294 double val(v.toDouble());
295 if (val<0.0)
296 throw Exception("Negative ks not allowed.");
297 ks_ = val;
298 return true;
299 }
300 case kwFric: {
301 if (!v.canConvert<double>())
302 throw Exception("fric must be a double.");
303 double val(v.toDouble());
304 if (val<0.0)
305 throw Exception("Negative fric not allowed.");
306 fric_ = val;
307 return false;
308 }
309 case kwCbTenF: {
310 if (!v.canConvert<double>())
311 throw Exception("cb_tenf must be a double.");
312 double val(v.toDouble());
313 if (val<0.0)
314 throw Exception("Negative cb_tenf not allowed.");
315 cb_tenF_ = val;
316 return false;
317 }
318 case kwCbShearF: {
319 if (!v.canConvert<double>())
320 throw Exception("cb_shearf must be a double.");
321 double val(v.toDouble());
322 if (val<0.0)
323 throw Exception("Negative cb_shearf not allowed.");
324 cb_shearF_ = val;
325 return false;
326 }
327 case kwLinF: {
328 if (!v.canConvert<DVect>())
329 throw Exception("lin_force must be a vector.");
330 DVect val(v.value<DVect>());
331 lin_F_ = val;
332 return false;
333 }
334 case kwLinMode: {
335 if (!v.canConvert<uint32>())
336 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
337 uint32 val(v.toUInt());
338 if (val>1)
339 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
340 lin_mode_ = val;
341 return false;
342 }
343 case kwRGap: {
344 if (!v.canConvert<double>())
345 throw Exception("Reference gap must be a double.");
346 double val(v.toDouble());
347 rgap_ = val;
348 return false;
349 }
350 case kwDpNRatio: {
351 if (!v.canConvert<double>())
352 throw Exception("dp_nratio must be a double.");
353 double val(v.toDouble());
354 if (val<0.0)
355 throw Exception("Negative dp_nratio not allowed.");
356 if (val == 0.0 && !dpProps_)
357 return false;
358 if (!dpProps_)
359 dpProps_ = NEWC(dpProps());
360 dpProps_->dp_nratio_ = val;
361 return true;
362 }
363 case kwDpSRatio: {
364 if (!v.canConvert<double>())
365 throw Exception("dp_sratio must be a double.");
366 double val(v.toDouble());
367 if (val<0.0)
368 throw Exception("Negative dp_sratio not allowed.");
369 if (val == 0.0 && !dpProps_)
370 return false;
371 if (!dpProps_)
372 dpProps_ = NEWC(dpProps());
373 dpProps_->dp_sratio_ = val;
374 return true;
375 }
376 case kwDpMode: {
377 if (!v.canConvert<int>())
378 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
379 int val(v.toInt());
380 if (val == 0 && !dpProps_)
381 return false;
382 if (val < 0 || val > 3)
383 throw Exception("The dashpot mode dp_mode must be 0, 1, 2, or 3.");
384 if (!dpProps_)
385 dpProps_ = NEWC(dpProps());
386 dpProps_->dp_mode_ = val;
387 return false;
388 }
389 case kwDpF: {
390 if (!v.canConvert<DVect>())
391 throw Exception("dp_force must be a vector.");
392 DVect val(v.value<DVect>());
393 if (val.fsum() == 0.0 && !dpProps_)
394 return false;
395 if (!dpProps_)
396 dpProps_ = NEWC(dpProps());
397 dpProps_->dp_F_ = val;
398 return false;
399 }
400 case kwUserArea: {
401 if (!v.canConvert<double>())
402 throw Exception("user_area must be a double.");
403 double val(v.toDouble());
404 if (val < 0.0)
405 throw Exception("Negative user_area not allowed.");
406 userArea_ = val;
407 return true;
408 }
409 }
410 return false;
411 }
412
413 bool ContactModelLinearCBond::getPropertyReadOnly(uint32 i) const {
414 switch (i) {
415 case kwDpF:
416 case kwLinS:
417 case kwEmod:
418 case kwKRatio:
419 case kwCbState:
420 case kwCbTStr:
421 case kwCbSStr:
422 return true;
423 default:
424 break;
425 }
426 return false;
427 }
428
429 bool ContactModelLinearCBond::supportsInheritance(uint32 i) const {
430 switch (i) {
431 case kwKn:
432 case kwKs:
433 case kwFric:
434 return true;
435 default:
436 break;
437 }
438 return false;
439 }
440
441 QString ContactModelLinearCBond::getMethodArguments(uint32 i) const {
442 switch (i) {
443 case kwCbBond:
444 return "gap";
445 case kwDeformability:
446 return "emod,kratio";
447 case kwCbStrength:
448 return "tensile,shear";
449 case kwCbUnbond:
450 return "gap";
451 case kwArea:
452 return QString();
453 }
454 assert(0);
455 return "";
456 }
457
458 bool ContactModelLinearCBond::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
459 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
460 switch (i) {
461 case kwCbBond: {
462 if (cb_state_ == 3) return false;
463 double mingap = -1.0 * limits<double>::max();
464 double maxgap = 0;
465 if (vl.at(0).canConvert<double>())
466 maxgap = vl.at(0).toDouble();
467 else if (vl.at(0).canConvert<DVect2>()) {
468 DVect2 value = vl.at(0).value<DVect2>();
469 mingap = value.minComp();
470 maxgap = value.maxComp();
471 } else if (!vl.at(0).isNull())
472 throw Exception("gap value %1 not recognized in method bond in contact model %2.",vl.at(0),getName());
473
474 double gap = c->getGap();
475 if ( gap >= mingap && gap <= maxgap)
476 cb_state_ = 3;
477 return false;
478 }
479 case kwCbUnbond: {
480 if (cb_state_ == 0) return false;
481 double mingap = -1.0 * limits<double>::max();
482 double maxgap = 0;
483 if (vl.at(0).canConvert<double>())
484 maxgap = vl.at(0).toDouble();
485 else if (vl.at(0).canConvert<DVect2>()) {
486 DVect2 value = vl.at(0).value<DVect2>();
487 mingap = value.minComp();
488 maxgap = value.maxComp();
489 }
490 else if (!vl.at(0).isNull())
491 throw Exception("gap value %1 not recognized in method unbond in contact model %2.",vl.at(0),getName());
492
493 double gap = c->getGap();
494 if ( gap >= mingap && gap <= maxgap)
495 cb_state_ = 0;
496 return false;
497 }
498 case kwDeformability: {
499 double emod(0.0);
500 double krat(0.0);
501 if (vl.at(0).isNull())
502 throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
503 emod = vl.at(0).toDouble();
504 if (emod<0.0)
505 throw Exception("Negative emod not allowed in contact model %1.",getName());
506 if (vl.at(1).isNull())
507 throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
508 krat = vl.at(1).toDouble();
509 if (krat<0.0)
510 throw Exception("Negative linear stiffness ratio not allowed in contact model %1.",getName());
511 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
512 double rsum(0.0);
513 if (c->getEnd1Curvature().y())
514 rsum += 1.0/c->getEnd1Curvature().y();
515 if (c->getEnd2Curvature().y())
516 rsum += 1.0/c->getEnd2Curvature().y();
517 if (userArea_) {
518#ifdef THREED
519 rsq = std::sqrt(userArea_ / dPi);
520#else
521 rsq = userArea_ / 2.0;
522#endif
523 rsum = rsq + rsq;
524 rsq = 1. / rsq;
525 }
526#ifdef TWOD
527 kn_ = 2.0 * emod / (rsq * rsum);
528#else
529 kn_ = dPi * emod / (rsq * rsq * rsum);
530#endif
531 ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
532 setInheritance(1,false);
533 setInheritance(2,false);
534 return true;
535 }
536 case kwCbStrength: {
537 if (cb_state_ != 3) return false;
538 double nval(0.0);
539 double sval(0.0);
540 if (vl.at(0).isNull())
541 throw Exception("tensile value must be specified with method cb_strength in contact model %1.",getName());
542 nval = vl.at(0).toDouble();
543 if (nval<0.0)
544 throw Exception("Negative tensile strength not allowed in contact model %1.",getName());
545 if (vl.at(1).isNull())
546 throw Exception("shear value must be specified with method cb_strength in contact model %1.",getName());
547 sval = vl.at(1).toDouble();
548 if (sval<0.0)
549 throw Exception("Negative shear strength not allowed in contact model %1.",getName());
550 double tmp(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
551 if (userArea_) {
552#ifdef THREED
553 tmp = std::sqrt(userArea_ / dPi);
554#else
555 tmp = userArea_ / 2.0;
556#endif
557 tmp = 1. / tmp;
558 }
559#ifdef TWOD
560 cb_tenF_ = nval * 2.0 / tmp;
561 cb_shearF_ = sval * 2.0 / tmp;
562#else
563 cb_tenF_ = nval * dPi / ( tmp * tmp );
564 cb_shearF_ = sval * dPi / (tmp * tmp);
565#endif
566 return false;
567 }
568 case kwArea: {
569 if (!userArea_) {
570 double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
571#ifdef THREED
572 userArea_ = rsq * rsq * dPi;
573#else
574 userArea_ = rsq * 2.0;
575#endif
576 }
577 return true;
578 }
579
580 }
581 return false;
582 }
583
584 double ContactModelLinearCBond::getEnergy(uint32 i) const {
585 double ret(0.0);
586 if (!energies_)
587 return ret;
588 switch (i) {
589 case kwEStrain: return energies_->estrain_;
590 case kwESlip: return energies_->eslip_;
591 case kwEDashpot: return energies_->edashpot_;
592 }
593 assert(0);
594 return ret;
595 }
596
597 bool ContactModelLinearCBond::getEnergyAccumulate(uint32 i) const {
598 switch (i) {
599 case kwEStrain: return false;
600 case kwESlip: return true;
601 case kwEDashpot: return true;
602 }
603 assert(0);
604 return false;
605 }
606
607 void ContactModelLinearCBond::setEnergy(uint32 i,const double &d) {
608 if (!energies_) return;
609 switch (i) {
610 case kwEStrain: energies_->estrain_ = d; return;
611 case kwESlip: energies_->eslip_ = d; return;
612 case kwEDashpot: energies_->edashpot_= d; return;
613 }
614 assert(0);
615 return;
616 }
617
618 bool ContactModelLinearCBond::validate(ContactModelMechanicalState *state,const double &) {
619 assert(state);
620 const IContactMechanical *c = state->getMechanicalContact();
621 assert(c);
622
623 if (state->trackEnergy_)
624 activateEnergy();
625
626 if (inheritanceField_ & linKnMask)
627 updateKn(c);
628 if (inheritanceField_ & linKsMask)
629 updateKs(c);
630 if (inheritanceField_ & linFricMask)
631 updateFric(c);
632
633 updateEffectiveStiffness(state);
634 return checkActivity(state->gap_);
635 }
636
637 static const QString knstr("kn");
638 bool ContactModelLinearCBond::updateKn(const IContactMechanical *con) {
639 assert(con);
640 QVariant v1 = con->getEnd1()->getProperty(knstr);
641 QVariant v2 = con->getEnd2()->getProperty(knstr);
642 if (!v1.isValid() || !v2.isValid())
643 return false;
644 double kn1 = v1.toDouble();
645 double kn2 = v2.toDouble();
646 double val = kn_;
647 if (kn1 && kn2)
648 kn_ = kn1*kn2/(kn1+kn2);
649 else if (kn1)
650 kn_ = kn1;
651 else if (kn2)
652 kn_ = kn2;
653 return ( (kn_ != val) );
654 }
655
656 static const QString ksstr("ks");
657 bool ContactModelLinearCBond::updateKs(const IContactMechanical *con) {
658 assert(con);
659 QVariant v1 = con->getEnd1()->getProperty(ksstr);
660 QVariant v2 = con->getEnd2()->getProperty(ksstr);
661 if (!v1.isValid() || !v2.isValid())
662 return false;
663 double ks1 = v1.toDouble();
664 double ks2 = v2.toDouble();
665 double val = ks_;
666 if (ks1 && ks2)
667 ks_ = ks1*ks2/(ks1+ks2);
668 else if (ks1)
669 ks_ = ks1;
670 else if (ks2)
671 ks_ = ks2;
672 return ( (ks_ != val) );
673 }
674
675 static const QString fricstr("fric");
676 bool ContactModelLinearCBond::updateFric(const IContactMechanical *con) {
677 assert(con);
678 QVariant v1 = con->getEnd1()->getProperty(fricstr);
679 QVariant v2 = con->getEnd2()->getProperty(fricstr);
680 if (!v1.isValid() || !v2.isValid())
681 return false;
682 double fric1 = std::max(0.0,v1.toDouble());
683 double fric2 = std::max(0.0,v2.toDouble());
684 double val = fric_;
685 fric_ = std::min(fric1,fric2);
686 return ( (fric_ != val) );
687 }
688
689 bool ContactModelLinearCBond::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
690 assert(c);
691 QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
692 QRegularExpression rx(name, QRegularExpression::CaseInsensitiveOption);
693 int idx = availableProperties.indexOf(rx)+1;
694 bool ret=false;
695
696 if (idx<=0)
697 return ret;
698
699 switch(idx) {
700 case kwKn: { //kn
701 if (inheritanceField_ & linKnMask)
702 ret = updateKn(c);
703 break;
704 }
705 case kwKs: { //ks
706 if (inheritanceField_ & linKsMask)
707 ret =updateKs(c);
708 break;
709 }
710 case kwFric: { //fric
711 if (inheritanceField_ & linFricMask)
712 updateFric(c);
713 break;
714 }
715 }
716 return ret;
717 }
718
719 void ContactModelLinearCBond::updateEffectiveStiffness(ContactModelMechanicalState *) {
720 DVect2 ret(kn_,ks_);
721 // correction if viscous damping active
722 if (dpProps_) {
723 DVect2 correct(1.0);
724 if (dpProps_->dp_nratio_)
725 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
726 if (dpProps_->dp_sratio_)
727 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
728 ret /= (correct*correct);
729 }
730 effectiveTranslationalStiffness_ = ret;
731 }
732
733 bool ContactModelLinearCBond::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
734 assert(state);
735
736 double overlap = rgap_ - state->gap_;
737 DVect trans = state->relativeTranslationalIncrement_;
738 double correction = 1.0;
739
740 if (state->activated()) {
741 if (cmEvents_[fActivated] >= 0) {
742 auto c = state->getContact();
743 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
744 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
745 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
746 }
747 if (lin_mode_ == 0 && trans.x()) {
748 correction = -1.0*overlap / trans.x();
749 if (correction < 0)
750 correction = 1.0;
751 }
752 }
753
754#ifdef THREED
755 DVect norm(trans.x(),0.0,0.0);
756#else
757 DVect norm(trans.x(),0.0);
758#endif
759 //DAVect ang = state->relativeAngularIncrement_;
760 DVect lin_F_old = lin_F_;
761
762 if (lin_mode_ == 0)
763 lin_F_.rx() = overlap * kn_;
764 else
765 lin_F_.rx() -= correction * norm.x() * kn_;
766
767 DVect u_s = trans;
768 u_s.rx() = 0.0;
769 DVect sforce = lin_F_ - u_s * ks_ * correction;
770 sforce.rx() = 0.0;
771
772 // Resolve failure (contact bonds and friction)
773 if (state->canFail_) {
774 // Resolve contact bond failure - done first so that this way, even if breaks, one can ensure a valid sliding state
775 if (cb_state_ == 3) { // bonded - Note: this means that isSliding is false!
776 if (lin_F_.x() <= -cb_tenF_) {
777 // Broke in tension
778 cb_state_ = 1;
779 if (cmEvents_[fBondBreak] >= 0) {
780 auto c = state->getContact();
781 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
782 fish::Parameter((qint64)cb_state_),
783 fish::Parameter(cb_tenF_) };
784 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
785 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
786 }
787 } else if (sforce.mag() >= cb_shearF_) {
788 // Broke in shear
789 cb_state_ = 2;
790 if (cmEvents_[fBondBreak] >= 0) {
791 auto c = state->getContact();
792 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
793 fish::Parameter((qint64)cb_state_),
794 fish::Parameter(cb_shearF_) };
795 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
796 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
797 }
798 }
799 }
800
801 // 2) Resolve sliding if no contact bond exists
802 if (cb_state_ < 3) {
803 // No contact bond - normal force is positive only
804 lin_F_.rx() = std::max(0.0,lin_F_.x());
805 // No contact bond - sliding can occur
806 double crit = lin_F_.x() * fric_;
807 double sfmag = sforce.mag();
808 if (sfmag > crit) {
809 double rat = crit / sfmag;
810 sforce *= rat;
811 if (!lin_S_ && cmEvents_[fSlipChange] >= 0) {
812 auto c = state->getContact();
813 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
814 fish::Parameter() };
815 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
816 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
817 }
818 lin_S_ = true;
819 } else {
820 if (lin_S_) {
821 if (cmEvents_[fSlipChange] >= 0) {
822 auto c = state->getContact();
823 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
824 fish::Parameter((qint64)1) };
825 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
826 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
827 }
828 lin_S_ = false;
829 }
830 }
831 }
832 }
833
834 sforce.rx() = lin_F_.x();
835 lin_F_ = sforce; // total force in linear contact model
836
837 // 3) Account for dashpot forces
838 if (dpProps_) {
839 dpProps_->dp_F_.fill(0.0);
840 double vcn(0.0), vcs(0.0);
841 setDampCoefficients(state->inertialMass_,&vcn,&vcs);
842 // First damp all components
843 dpProps_->dp_F_ = u_s * (-1.0* vcs) / timestep; // shear component
844 dpProps_->dp_F_ -= norm * vcn / timestep; // normal component
845 // Need to change behavior based on the dp_mode
846 if (cb_state_ !=3 && (dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) { // limit the tensile if not bonded
847 if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
848 dpProps_->dp_F_.rx() = - lin_F_.rx();
849 }
850 if (lin_S_ && dpProps_->dp_mode_ > 1) { // limit the shear if not sliding
851 double dfn = dpProps_->dp_F_.rx();
852 dpProps_->dp_F_.fill(0.0);
853 dpProps_->dp_F_.rx() = dfn;
854 }
855 }
856
857 // 5) Compute energies
858 if (state->trackEnergy_) {
859 assert(energies_);
860 energies_->estrain_ = 0.0;
861 if (kn_)
862 energies_->estrain_ = 0.5*lin_F_.x()*lin_F_.x()/kn_;
863 if (ks_) {
864 DVect s = lin_F_;
865 s.rx() = 0.0;
866 double smag2 = s.mag2();
867 energies_->estrain_ += 0.5*smag2 / ks_;
868
869 if (lin_S_) {
870 lin_F_old.rx() = 0.0;
871 DVect avg_F_s = (s + lin_F_old)*0.5;
872 DVect u_s_el = (s - lin_F_old) / ks_;
873 energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
874 }
875 }
876 if (dpProps_) {
877 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
878 }
879 }
880 assert(lin_F_ == lin_F_);
881 return checkActivity(state->gap_);
882 }
883
884 bool ContactModelLinearCBond::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
885 // Account for thermal expansion in incremental mode
886 if (lin_mode_ == 0 || ts->gapInc_ == 0.0) return false;
887 DVect finc(0.0);
888 finc.rx() = kn_ * ts->gapInc_;
889 lin_F_ -= finc;
890 return true;
891 }
892
893 void ContactModelLinearCBond::setForce(const DVect &v,IContact *c) {
894 lin_F(v);
895 if (v.x() > 0)
896 rgap_ = c->getGap() + v.x() / kn_;
897 }
898
899 void ContactModelLinearCBond::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
900 // Only do something if the contact model is of the same type
901 if (old->getContactModel()->getName().compare("linearcbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
902 ContactModelLinearCBond *oldCm = (ContactModelLinearCBond *)old;
903#ifdef THREED
904 // Need to rotate just the shear component from oldSystem to newSystem
905
906 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
907 DVect axis = oldSystem.e1() & newSystem.e1();
908 double c, ang, s;
909 DVect re2;
910 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
911 axis = axis.unit();
912 c = oldSystem.e1()|newSystem.e1();
913 if (c > 0)
914 c = std::min(c,1.0);
915 else
916 c = std::max(c,-1.0);
917 ang = acos(c);
918 s = sin(ang);
919 double t = 1. - c;
920 DMatrix<3,3> rm;
921 rm.get(0,0) = t*axis.x()*axis.x() + c;
922 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
923 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
924 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
925 rm.get(1,1) = t*axis.y()*axis.y() + c;
926 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
927 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
928 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
929 rm.get(2,2) = t*axis.z()*axis.z() + c;
930 re2 = rm*oldSystem.e2();
931 }
932 else
933 re2 = oldSystem.e2();
934 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
935 axis = re2 & newSystem.e2();
936 DVect2 tpf;
937 DMatrix<2,2> m;
938 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
939 axis = axis.unit();
940 c = re2|newSystem.e2();
941 if (c > 0)
942 c = std::min(c,1.0);
943 else
944 c = std::max(c,-1.0);
945 ang = acos(c);
946 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
947 ang *= -1;
948 s = sin(ang);
949 m.get(0,0) = c;
950 m.get(1,0) = s;
951 m.get(0,1) = -m.get(1,0);
952 m.get(1,1) = m.get(0,0);
953 tpf = m*DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
954 } else {
955 m.get(0,0) = 1.;
956 m.get(0,1) = 0.;
957 m.get(1,0) = 0.;
958 m.get(1,1) = 1.;
959 tpf = DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
960 }
961 DVect pforce = DVect(0,tpf.x(),tpf.y());
962#else
963 oldSystem;
964 newSystem;
965 DVect pforce = DVect(0,oldCm->lin_F_.y());
966#endif
967 for (int i=1; i<dim; ++i)
968 lin_F_.rdof(i) += pforce.dof(i);
969 if (lin_mode_ && oldCm->lin_mode_)
970 lin_F_.rx() = oldCm->lin_F_.x();
971 oldCm->lin_F_ = DVect(0.0);
972 if (dpProps_ && oldCm->dpProps_) {
973#ifdef THREED
974 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
975 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
976#else
977 pforce = oldCm->dpProps_->dp_F_;
978#endif
979 dpProps_->dp_F_ += pforce;
980 oldCm->dpProps_->dp_F_ = DVect(0.0);
981 }
982 if(oldCm->getEnergyActivated()) {
983 activateEnergy();
984 energies_->estrain_ = oldCm->energies_->estrain_;
985 energies_->eslip_ = oldCm->energies_->eslip_;
986 energies_->edashpot_ = oldCm->energies_->edashpot_;
987 oldCm->energies_->estrain_ = 0.0;
988 oldCm->energies_->edashpot_ = 0.0;
989 oldCm->energies_->eslip_ = 0.0;
990 }
991 rgap_ = oldCm->rgap_;
992 }
993 assert(lin_F_ == lin_F_);
994 }
995
996 void ContactModelLinearCBond::setNonForcePropsFrom(IContactModel *old) {
997 // Only do something if the contact model is of the same type
998 if (old->getName().compare("linearcbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
999 ContactModelLinearCBond *oldCm = (ContactModelLinearCBond *)old;
1000 kn_ = oldCm->kn_;
1001 ks_ = oldCm->ks_;
1002 fric_ = oldCm->fric_;
1003 lin_mode_ = oldCm->lin_mode_;
1004 rgap_ = oldCm->rgap_;
1005 userArea_ = oldCm->userArea_;
1006
1007 if (oldCm->dpProps_) {
1008 if (!dpProps_)
1009 dpProps_ = NEWC(dpProps());
1010 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1011 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1012 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1013 }
1014 }
1015 }
1016
1017 DVect ContactModelLinearCBond::getForce(const IContactMechanical *) const {
1018 DVect ret(lin_F_);
1019 if (dpProps_)
1020 ret += dpProps_->dp_F_;
1021 return ret;
1022 }
1023
1024 DAVect ContactModelLinearCBond::getMomentOn1(const IContactMechanical *c) const {
1025 DVect force = getForce(c);
1026 DAVect ret(0.0);
1027 c->updateResultingTorqueOn1Local(force,&ret);
1028 return ret;
1029 }
1030
1031 DAVect ContactModelLinearCBond::getMomentOn2(const IContactMechanical *c) const {
1032 DVect force = getForce(c);
1033 DAVect ret(0.0);
1034 c->updateResultingTorqueOn2Local(force,&ret);
1035 return ret;
1036 }
1037
1038 void ContactModelLinearCBond::setDampCoefficients(const double &mass,double *vcn,double *vcs) {
1039 *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(kn_));
1040 *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(ks_));
1041 }
1042
1043} // namespace itascaxd
1044// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Aug 13, 2024 |
