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