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 string 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 string 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 string 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_ = NEW Energies();}
81 bool getEnergyActivated() const override {return (energies_ !=0);}
82
83 enum FishCallEvents {fActivated=0,fBondBreak, fSlipChange };
84 string getFishCallEvents() const override { return "contact_activated,bond_break,slip_change"; }
85 base::Property getProperty(uint32 i,const IContact *) const override;
86 bool getPropertyGlobal(uint32 i) const override;
87 bool setProperty(uint32 i,const base::Property &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 string getMethods() const override {
103 return "deformability"
104 ",bond"
105 ",cb_strength"
106 ",unbond"
107 ",area";
108 }
109
110 string getMethodArguments(uint32 i) const override;
111
112 bool setMethod(uint32 i,const std::vector<base::Property> &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 string &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 NEW 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<string,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 = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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_ = NEW 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 base::Property 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 base::Property 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 = calcRSum(c);
203 if (userArea_) {
204#ifdef THREED
205 rsq = std::sqrt(userArea_ / dPi);
206#else
207 rsq = userArea_ / 2.0;
208#endif
209 rsum = rsq + rsq;
210 rsq = 1. / rsq;
211 }
212#ifdef TWOD
213 return (kn_ * rsum * rsq / 2.0);
214#else
215 return (kn_ * rsum * rsq * rsq) / dPi;
216#endif
217 }
218 case kwKRatio: return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
219 case kwDpNRatio: return dpProps_ ? dpProps_->dp_nratio_ : 0;
220 case kwDpSRatio: return dpProps_ ? dpProps_->dp_sratio_ : 0;
221 case kwDpMode: return dpProps_ ? dpProps_->dp_mode_ : 0;
222 case kwDpF: {
223 dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
224 return var;
225 }
226 case kwCbState: return cb_state_;
227 case kwCbTenF: return cb_tenF_;
228 case kwCbShearF: return cb_shearF_;
229 case kwCbTStr: nstr = true;
230 [[fallthrough]];
231 case kwCbSStr: {
232 if (c ==nullptr) return 0.0;
233 double tmp(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
234 if (userArea_) {
235#ifdef THREED
236 tmp = std::sqrt(userArea_ / dPi);
237#else
238 tmp = userArea_ / 2.0;
239#endif
240 tmp = 1. / tmp;
241 }
242 if (nstr) {
243#ifdef TWOD
244 return (cb_tenF_ * tmp / 2.0);
245#else
246 return (cb_tenF_ * tmp * tmp / dPi);
247#endif
248 } else {
249#ifdef TWOD
250 return (cb_shearF_ * tmp / 2.0);
251#else
252 return (cb_shearF_ * tmp * tmp / dPi);
253#endif
254 }
255 }
256 case kwUserArea: return userArea_;
257 }
258 assert(0);
259 return base::Property();
260 }
261
262 bool ContactModelLinearCBond::getPropertyGlobal(uint32 i) const {
263 switch (i) {
264 case kwLinF:
265 case kwDpF:
266 return false;
267 }
268 return true;
269 }
270
271 bool ContactModelLinearCBond::setProperty(uint32 i,const base::Property &v,IContact *) {
272 dpProps dp;
273 switch (i) {
274 case kwKn: {
275 if (!v.canConvert<double>())
276 throw Exception("kn must be a double.");
277 double val(v.toDouble());
278 if (val<0.0)
279 throw Exception("Negative kn not allowed.");
280 kn_ = val;
281 return true;
282 }
283 case kwKs: {
284 if (!v.canConvert<double>())
285 throw Exception("ks must be a double.");
286 double val(v.toDouble());
287 if (val<0.0)
288 throw Exception("Negative ks not allowed.");
289 ks_ = val;
290 return true;
291 }
292 case kwFric: {
293 if (!v.canConvert<double>())
294 throw Exception("fric must be a double.");
295 double val(v.toDouble());
296 if (val<0.0)
297 throw Exception("Negative fric not allowed.");
298 fric_ = val;
299 return false;
300 }
301 case kwCbTenF: {
302 if (!v.canConvert<double>())
303 throw Exception("cb_tenf must be a double.");
304 double val(v.toDouble());
305 if (val<0.0)
306 throw Exception("Negative cb_tenf not allowed.");
307 cb_tenF_ = val;
308 return false;
309 }
310 case kwCbShearF: {
311 if (!v.canConvert<double>())
312 throw Exception("cb_shearf must be a double.");
313 double val(v.toDouble());
314 if (val<0.0)
315 throw Exception("Negative cb_shearf not allowed.");
316 cb_shearF_ = val;
317 return false;
318 }
319 case kwLinF: {
320 if (!v.canConvert<DVect>())
321 throw Exception("lin_force must be a vector.");
322 DVect val(v.value<DVect>());
323 lin_F_ = val;
324 return false;
325 }
326 case kwLinMode: {
327 if (!v.canConvert<int64>())
328 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
329 uint32 val(v.toUInt());
330 if (val>1)
331 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
332 lin_mode_ = val;
333 return false;
334 }
335 case kwRGap: {
336 if (!v.canConvert<double>())
337 throw Exception("Reference gap must be a double.");
338 double val(v.toDouble());
339 rgap_ = val;
340 return false;
341 }
342 case kwDpNRatio: {
343 if (!v.canConvert<double>())
344 throw Exception("dp_nratio must be a double.");
345 double val(v.toDouble());
346 if (val<0.0)
347 throw Exception("Negative dp_nratio not allowed.");
348 if (val == 0.0 && !dpProps_)
349 return false;
350 if (!dpProps_)
351 dpProps_ = NEW dpProps();
352 dpProps_->dp_nratio_ = val;
353 return true;
354 }
355 case kwDpSRatio: {
356 if (!v.canConvert<double>())
357 throw Exception("dp_sratio must be a double.");
358 double val(v.toDouble());
359 if (val<0.0)
360 throw Exception("Negative dp_sratio not allowed.");
361 if (val == 0.0 && !dpProps_)
362 return false;
363 if (!dpProps_)
364 dpProps_ = NEW dpProps();
365 dpProps_->dp_sratio_ = val;
366 return true;
367 }
368 case kwDpMode: {
369 if (!v.canConvert<int64>())
370 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
371 int val(v.toInt());
372 if (val == 0 && !dpProps_)
373 return false;
374 if (val < 0 || val > 3)
375 throw Exception("The dashpot mode dp_mode must be 0, 1, 2, or 3.");
376 if (!dpProps_)
377 dpProps_ = NEW dpProps();
378 dpProps_->dp_mode_ = val;
379 return false;
380 }
381 case kwDpF: {
382 if (!v.canConvert<DVect>())
383 throw Exception("dp_force must be a vector.");
384 DVect val(v.value<DVect>());
385 if (val.fsum() == 0.0 && !dpProps_)
386 return false;
387 if (!dpProps_)
388 dpProps_ = NEW dpProps();
389 dpProps_->dp_F_ = val;
390 return false;
391 }
392 case kwUserArea: {
393 if (!v.canConvert<double>())
394 throw Exception("user_area must be a double.");
395 double val(v.toDouble());
396 if (val < 0.0)
397 throw Exception("Negative user_area not allowed.");
398 userArea_ = val;
399 return true;
400 }
401 }
402 return false;
403 }
404
405 bool ContactModelLinearCBond::getPropertyReadOnly(uint32 i) const {
406 switch (i) {
407 case kwDpF:
408 case kwLinS:
409 case kwEmod:
410 case kwKRatio:
411 case kwCbState:
412 case kwCbTStr:
413 case kwCbSStr:
414 return true;
415 default:
416 break;
417 }
418 return false;
419 }
420
421 bool ContactModelLinearCBond::supportsInheritance(uint32 i) const {
422 switch (i) {
423 case kwKn:
424 case kwKs:
425 case kwFric:
426 return true;
427 default:
428 break;
429 }
430 return false;
431 }
432
433 string ContactModelLinearCBond::getMethodArguments(uint32 i) const {
434 switch (i) {
435 case kwCbBond:
436 return "gap";
437 case kwDeformability:
438 return "emod,kratio";
439 case kwCbStrength:
440 return "tensile,shear";
441 case kwCbUnbond:
442 return "gap";
443 case kwArea:
444 return string();
445 }
446 assert(0);
447 return "";
448 }
449
450 bool ContactModelLinearCBond::setMethod(uint32 i,const std::vector<base::Property> &vl,IContact *con) {
451 FP_S;
452 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
453 switch (i) {
454 case kwCbBond: {
455 if (cb_state_ == 3) return false;
456 double mingap = -1.0 * limits<double>::max();
457 double maxgap = 0;
458 if (vl.at(0).canConvert<double>())
459 maxgap = vl.at(0).toDouble();
460 else if (vl.at(0).canConvert<DVect2>()) {
461 DVect2 value = vl.at(0).value<DVect2>();
462 mingap = value.minComp();
463 maxgap = value.maxComp();
464 } else if (!vl.at(0).isNull())
465 throw Exception("gap value {0} not recognized in method bond in contact model {1}.",vl.at(0),getName());
466
467 double gap = c->getGap();
468 if ( gap >= mingap && gap <= maxgap)
469 cb_state_ = 3;
470 FP_S;
471 return false;
472 }
473 case kwCbUnbond: {
474 if (cb_state_ == 0) return false;
475 double mingap = -1.0 * limits<double>::max();
476 double maxgap = 0;
477 if (vl.at(0).canConvert<double>())
478 maxgap = vl.at(0).toDouble();
479 else if (vl.at(0).canConvert<DVect2>()) {
480 DVect2 value = vl.at(0).value<DVect2>();
481 mingap = value.minComp();
482 maxgap = value.maxComp();
483 }
484 else if (!vl.at(0).isNull())
485 throw Exception("gap value {0} not recognized in method unbond in contact model {1}.",vl.at(0),getName());
486
487 double gap = c->getGap();
488 if ( gap >= mingap && gap <= maxgap)
489 cb_state_ = 0;
490 FP_S;
491 return false;
492 }
493 case kwDeformability: {
494 double emod(0.0);
495 double krat(0.0);
496 if (vl.at(0).isNull())
497 throw Exception("Argument emod must be specified with method deformability in contact model {0}.",getName());
498 emod = vl.at(0).toDouble();
499 if (emod<0.0)
500 throw Exception("Negative emod not allowed in contact model {0}.",getName());
501 if (vl.at(1).isNull())
502 throw Exception("Argument kratio must be specified with method deformability in contact model {0}.",getName());
503 krat = vl.at(1).toDouble();
504 if (krat<0.0)
505 throw Exception("Negative linear stiffness ratio not allowed in contact model {0}.",getName());
506 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
507 double rsum = calcRSum(c);
508 if (userArea_) {
509#ifdef THREED
510 rsq = std::sqrt(userArea_ / dPi);
511#else
512 rsq = userArea_ / 2.0;
513#endif
514 rsum = rsq + rsq;
515 rsq = 1. / rsq;
516 }
517#ifdef TWOD
518 kn_ = 2.0 * emod / (rsq * rsum);
519#else
520 kn_ = dPi * emod / (rsq * rsq * rsum);
521#endif
522 ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
523 setInheritance(1,false);
524 setInheritance(2,false);
525 FP_S;
526 return true;
527 }
528 case kwCbStrength: {
529 if (cb_state_ != 3) return false;
530 double nval(0.0);
531 double sval(0.0);
532 if (vl.at(0).isNull())
533 throw Exception("tensile value must be specified with method cb_strength in contact model {0}.",getName());
534 nval = vl.at(0).toDouble();
535 if (nval<0.0)
536 throw Exception("Negative tensile strength not allowed in contact model {0}.",getName());
537 if (vl.at(1).isNull())
538 throw Exception("shear value must be specified with method cb_strength in contact model {0}.",getName());
539 sval = vl.at(1).toDouble();
540 if (sval<0.0)
541 throw Exception("Negative shear strength not allowed in contact model {0}.",getName());
542 double tmp(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
543 if (userArea_) {
544#ifdef THREED
545 tmp = std::sqrt(userArea_ / dPi);
546#else
547 tmp = userArea_ / 2.0;
548#endif
549 tmp = 1. / tmp;
550 }
551#ifdef TWOD
552 cb_tenF_ = nval * 2.0 / tmp;
553 cb_shearF_ = sval * 2.0 / tmp;
554#else
555 cb_tenF_ = nval * dPi / ( tmp * tmp );
556 cb_shearF_ = sval * dPi / (tmp * tmp);
557#endif
558 FP_S;
559 return false;
560 }
561 case kwArea: {
562 if (!userArea_) {
563 double rsq = calcRSQ(c);
564#ifdef THREED
565 userArea_ = rsq * rsq * dPi;
566#else
567 userArea_ = rsq * 2.0;
568#endif
569 }
570 FP_S;
571 return true;
572 }
573
574 }
575 FP_S;
576 return false;
577 }
578
579 double ContactModelLinearCBond::getEnergy(uint32 i) const {
580 double ret(0.0);
581 if (!energies_)
582 return ret;
583 switch (i) {
584 case kwEStrain: return energies_->estrain_;
585 case kwESlip: return energies_->eslip_;
586 case kwEDashpot: return energies_->edashpot_;
587 }
588 assert(0);
589 return ret;
590 }
591
592 bool ContactModelLinearCBond::getEnergyAccumulate(uint32 i) const {
593 switch (i) {
594 case kwEStrain: return false;
595 case kwESlip: return true;
596 case kwEDashpot: return true;
597 }
598 assert(0);
599 return false;
600 }
601
602 void ContactModelLinearCBond::setEnergy(uint32 i,const double &d) {
603 if (!energies_) return;
604 switch (i) {
605 case kwEStrain: energies_->estrain_ = d; return;
606 case kwESlip: energies_->eslip_ = d; return;
607 case kwEDashpot: energies_->edashpot_= d; return;
608 }
609 assert(0);
610 return;
611 }
612
613 bool ContactModelLinearCBond::validate(ContactModelMechanicalState *state,const double &) {
614 assert(state);
615 const IContactMechanical *c = state->getMechanicalContact();
616 assert(c);
617
618 if (state->trackEnergy_)
619 activateEnergy();
620
621 if (inheritanceField_ & linKnMask)
622 updateKn(c);
623 if (inheritanceField_ & linKsMask)
624 updateKs(c);
625 if (inheritanceField_ & linFricMask)
626 updateFric(c);
627
628 updateEffectiveStiffness(state);
629 return checkActivity(state->gap_);
630 }
631
632 static const string knstr("kn");
633 bool ContactModelLinearCBond::updateKn(const IContactMechanical *con) {
634 assert(con);
635 base::Property v1 = con->getEnd1()->getProperty(knstr);
636 base::Property v2 = con->getEnd2()->getProperty(knstr);
637 if (!v1.isValid() || !v2.isValid())
638 return false;
639 double kn1 = v1.toDouble();
640 double kn2 = v2.toDouble();
641 double val = kn_;
642 if (kn1 && kn2)
643 kn_ = kn1*kn2/(kn1+kn2);
644 else if (kn1)
645 kn_ = kn1;
646 else if (kn2)
647 kn_ = kn2;
648 return ( (kn_ != val) );
649 }
650
651 static const string ksstr("ks");
652 bool ContactModelLinearCBond::updateKs(const IContactMechanical *con) {
653 assert(con);
654 base::Property v1 = con->getEnd1()->getProperty(ksstr);
655 base::Property v2 = con->getEnd2()->getProperty(ksstr);
656 if (!v1.isValid() || !v2.isValid())
657 return false;
658 double ks1 = v1.toDouble();
659 double ks2 = v2.toDouble();
660 double val = ks_;
661 if (ks1 && ks2)
662 ks_ = ks1*ks2/(ks1+ks2);
663 else if (ks1)
664 ks_ = ks1;
665 else if (ks2)
666 ks_ = ks2;
667 return ( (ks_ != val) );
668 }
669
670 static const string fricstr("fric");
671 bool ContactModelLinearCBond::updateFric(const IContactMechanical *con) {
672 assert(con);
673 base::Property v1 = con->getEnd1()->getProperty(fricstr);
674 base::Property v2 = con->getEnd2()->getProperty(fricstr);
675 if (!v1.isValid() || !v2.isValid())
676 return false;
677 double fric1 = std::max(0.0,v1.toDouble());
678 double fric2 = std::max(0.0,v2.toDouble());
679 double val = fric_;
680 fric_ = std::min(fric1,fric2);
681 return ( (fric_ != val) );
682 }
683
684 bool ContactModelLinearCBond::endPropertyUpdated(const string &name,const IContactMechanical *c) {
685 assert(c);
686 StringList availableProperties = split(replace(simplified(getProperties())," ",""),",");
687 auto idx = findRegex(availableProperties,name);
688 if (idx==string::npos) return false;
689 idx += 1;
690 bool ret=false;
691 switch(idx) {
692 case kwKn: { //kn
693 if (inheritanceField_ & linKnMask)
694 ret = updateKn(c);
695 break;
696 }
697 case kwKs: { //ks
698 if (inheritanceField_ & linKsMask)
699 ret =updateKs(c);
700 break;
701 }
702 case kwFric: { //fric
703 if (inheritanceField_ & linFricMask)
704 updateFric(c);
705 break;
706 }
707 }
708 return ret;
709 }
710
711 void ContactModelLinearCBond::updateEffectiveStiffness(ContactModelMechanicalState *) {
712 DVect2 ret(kn_,ks_);
713 // correction if viscous damping active
714 if (dpProps_) {
715 DVect2 correct(1.0);
716 if (dpProps_->dp_nratio_)
717 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
718 if (dpProps_->dp_sratio_)
719 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
720 ret /= (correct*correct);
721 }
722 effectiveTranslationalStiffness_ = ret;
723 }
724
725 bool ContactModelLinearCBond::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
726 assert(state);
727
728 double overlap = rgap_ - state->gap_;
729 DVect trans = state->relativeTranslationalIncrement_;
730 double correction = 1.0;
731
732 if (state->activated()) {
733 if (cmEvents_[fActivated] >= 0) {
734 auto c = state->getContact();
735 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
736 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
737 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
738 }
739 if (lin_mode_ == 0 && trans.x()) {
740 correction = -1.0*overlap / trans.x();
741 if (correction < 0)
742 correction = 1.0;
743 }
744 }
745
746#ifdef THREED
747 DVect norm(trans.x(),0.0,0.0);
748#else
749 DVect norm(trans.x(),0.0);
750#endif
751 //DAVect ang = state->relativeAngularIncrement_;
752 DVect lin_F_old = lin_F_;
753
754 if (lin_mode_ == 0)
755 lin_F_.rx() = overlap * kn_;
756 else
757 lin_F_.rx() -= correction * norm.x() * kn_;
758
759 DVect u_s = trans;
760 u_s.rx() = 0.0;
761 DVect sforce = lin_F_ - u_s * ks_ * correction;
762 sforce.rx() = 0.0;
763
764 // Resolve failure (contact bonds and friction)
765 if (state->canFail_) {
766 // Resolve contact bond failure - done first so that this way, even if breaks, one can ensure a valid sliding state
767 if (cb_state_ == 3) { // bonded - Note: this means that isSliding is false!
768 if (lin_F_.x() <= -cb_tenF_) {
769 // Broke in tension
770 cb_state_ = 1;
771 if (cmEvents_[fBondBreak] >= 0) {
772 auto c = state->getContact();
773 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
774 fish::Parameter((int64)cb_state_),
775 fish::Parameter(cb_tenF_) };
776 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
777 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
778 }
779 } else if (sforce.mag() >= cb_shearF_) {
780 // Broke in shear
781 cb_state_ = 2;
782 if (cmEvents_[fBondBreak] >= 0) {
783 auto c = state->getContact();
784 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
785 fish::Parameter((int64)cb_state_),
786 fish::Parameter(cb_shearF_) };
787 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
788 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
789 }
790 }
791 }
792
793 // 2) Resolve sliding if no contact bond exists
794 if (cb_state_ < 3) {
795 // No contact bond - normal force is positive only
796 lin_F_.rx() = std::max(0.0,lin_F_.x());
797 // No contact bond - sliding can occur
798 double crit = lin_F_.x() * fric_;
799 double sfmag = sforce.mag();
800 if (sfmag > crit) {
801 double rat = crit / sfmag;
802 sforce *= rat;
803 if (!lin_S_ && cmEvents_[fSlipChange] >= 0) {
804 auto c = state->getContact();
805 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
806 fish::Parameter() };
807 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
808 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
809 }
810 lin_S_ = true;
811 } else {
812 if (lin_S_) {
813 if (cmEvents_[fSlipChange] >= 0) {
814 auto c = state->getContact();
815 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
816 fish::Parameter((int64)1) };
817 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
818 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
819 }
820 lin_S_ = false;
821 }
822 }
823 }
824 }
825
826 sforce.rx() = lin_F_.x();
827 lin_F_ = sforce; // total force in linear contact model
828
829 // 3) Account for dashpot forces
830 if (dpProps_) {
831 dpProps_->dp_F_.fill(0.0);
832 double vcn(0.0), vcs(0.0);
833 setDampCoefficients(state->inertialMass_,&vcn,&vcs);
834 // First damp all components
835 dpProps_->dp_F_ = u_s * (-1.0* vcs) / timestep; // shear component
836 dpProps_->dp_F_ -= norm * vcn / timestep; // normal component
837 // Need to change behavior based on the dp_mode
838 if (cb_state_ !=3 && (dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) { // limit the tensile if not bonded
839 if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
840 dpProps_->dp_F_.rx() = - lin_F_.rx();
841 }
842 if (lin_S_ && dpProps_->dp_mode_ > 1) { // limit the shear if not sliding
843 double dfn = dpProps_->dp_F_.rx();
844 dpProps_->dp_F_.fill(0.0);
845 dpProps_->dp_F_.rx() = dfn;
846 }
847 }
848
849 // 5) Compute energies
850 if (state->trackEnergy_) {
851 assert(energies_);
852 energies_->estrain_ = 0.0;
853 if (kn_)
854 energies_->estrain_ = 0.5*lin_F_.x()*lin_F_.x()/kn_;
855 if (ks_) {
856 DVect s = lin_F_;
857 s.rx() = 0.0;
858 double smag2 = s.mag2();
859 energies_->estrain_ += 0.5*smag2 / ks_;
860
861 if (lin_S_) {
862 lin_F_old.rx() = 0.0;
863 DVect avg_F_s = (s + lin_F_old)*0.5;
864 DVect u_s_el = (s - lin_F_old) / ks_;
865 energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
866 }
867 }
868 if (dpProps_) {
869 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
870 }
871 }
872 assert(lin_F_ == lin_F_);
873 return checkActivity(state->gap_);
874 }
875
876 bool ContactModelLinearCBond::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
877 // Account for thermal expansion in incremental mode
878 if (lin_mode_ == 0 || ts->gapInc_ == 0.0) return false;
879 DVect finc(0.0);
880 finc.rx() = kn_ * ts->gapInc_;
881 lin_F_ -= finc;
882 return true;
883 }
884
885 void ContactModelLinearCBond::setForce(const DVect &v,IContact *c) {
886 lin_F(v);
887 if (v.x() > 0)
888 rgap_ = c->getGap() + v.x() / kn_;
889 }
890
891 void ContactModelLinearCBond::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
892 // Only do something if the contact model is of the same type
893 if (equal(old->getContactModel()->getName(),"linearcbond") && !isBonded()) {
894 ContactModelLinearCBond *oldCm = (ContactModelLinearCBond *)old;
895#ifdef THREED
896 // Need to rotate just the shear component from oldSystem to newSystem
897
898 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
899 DVect axis = oldSystem.e1() & newSystem.e1();
900 double c, ang, s;
901 DVect re2;
902 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
903 axis = axis.unit();
904 c = oldSystem.e1()|newSystem.e1();
905 if (c > 0)
906 c = std::min(c,1.0);
907 else
908 c = std::max(c,-1.0);
909 ang = acos(c);
910 s = sin(ang);
911 double t = 1. - c;
912 DMatrix<3,3> rm;
913 rm.get(0,0) = t*axis.x()*axis.x() + c;
914 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
915 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
916 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
917 rm.get(1,1) = t*axis.y()*axis.y() + c;
918 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
919 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
920 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
921 rm.get(2,2) = t*axis.z()*axis.z() + c;
922 re2 = rm*oldSystem.e2();
923 }
924 else
925 re2 = oldSystem.e2();
926 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
927 axis = re2 & newSystem.e2();
928 DVect2 tpf;
929 DMatrix<2,2> m;
930 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
931 axis = axis.unit();
932 c = re2|newSystem.e2();
933 if (c > 0)
934 c = std::min(c,1.0);
935 else
936 c = std::max(c,-1.0);
937 ang = acos(c);
938 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
939 ang *= -1;
940 s = sin(ang);
941 m.get(0,0) = c;
942 m.get(1,0) = s;
943 m.get(0,1) = -m.get(1,0);
944 m.get(1,1) = m.get(0,0);
945 tpf = m*DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
946 } else {
947 m.get(0,0) = 1.;
948 m.get(0,1) = 0.;
949 m.get(1,0) = 0.;
950 m.get(1,1) = 1.;
951 tpf = DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
952 }
953 DVect pforce = DVect(0,tpf.x(),tpf.y());
954#else
955 oldSystem;
956 newSystem;
957 DVect pforce = DVect(0,oldCm->lin_F_.y());
958#endif
959 for (int i=1; i<dim; ++i)
960 lin_F_.rdof(i) += pforce.dof(i);
961 if (lin_mode_ && oldCm->lin_mode_)
962 lin_F_.rx() = oldCm->lin_F_.x();
963 oldCm->lin_F_ = DVect(0.0);
964 if (dpProps_ && oldCm->dpProps_) {
965#ifdef THREED
966 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
967 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
968#else
969 pforce = oldCm->dpProps_->dp_F_;
970#endif
971 dpProps_->dp_F_ += pforce;
972 oldCm->dpProps_->dp_F_ = DVect(0.0);
973 }
974 if(oldCm->getEnergyActivated()) {
975 activateEnergy();
976 energies_->estrain_ = oldCm->energies_->estrain_;
977 energies_->eslip_ = oldCm->energies_->eslip_;
978 energies_->edashpot_ = oldCm->energies_->edashpot_;
979 oldCm->energies_->estrain_ = 0.0;
980 oldCm->energies_->edashpot_ = 0.0;
981 oldCm->energies_->eslip_ = 0.0;
982 }
983 rgap_ = oldCm->rgap_;
984 }
985 assert(lin_F_ == lin_F_);
986 }
987
988 void ContactModelLinearCBond::setNonForcePropsFrom(IContactModel *old) {
989 // Only do something if the contact model is of the same type
990 if (equal(old->getName(),"linearcbond") && !isBonded()) {
991 ContactModelLinearCBond *oldCm = (ContactModelLinearCBond *)old;
992 kn_ = oldCm->kn_;
993 ks_ = oldCm->ks_;
994 fric_ = oldCm->fric_;
995 lin_mode_ = oldCm->lin_mode_;
996 rgap_ = oldCm->rgap_;
997 userArea_ = oldCm->userArea_;
998
999 if (oldCm->dpProps_) {
1000 if (!dpProps_)
1001 dpProps_ = NEW dpProps();
1002 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1003 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1004 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1005 }
1006 }
1007 }
1008
1009 DVect ContactModelLinearCBond::getForce() const {
1010 DVect ret(lin_F_);
1011 if (dpProps_)
1012 ret += dpProps_->dp_F_;
1013 return ret;
1014 }
1015
1016 DAVect ContactModelLinearCBond::getMomentOn1(const IContactMechanical *c) const {
1017 DVect force = getForce();
1018 DAVect ret(0.0);
1019 c->updateResultingTorqueOn1Local(force,&ret);
1020 return ret;
1021 }
1022
1023 DAVect ContactModelLinearCBond::getMomentOn2(const IContactMechanical *c) const {
1024 DVect force = getForce();
1025 DAVect ret(0.0);
1026 c->updateResultingTorqueOn2Local(force,&ret);
1027 return ret;
1028 }
1029
1030 DAVect ContactModelLinearCBond::getModelMomentOn1() const {
1031 DAVect ret(0.0);
1032 return ret;
1033 }
1034
1035 DAVect ContactModelLinearCBond::getModelMomentOn2() const {
1036 DAVect ret(0.0);
1037 return ret;
1038 }
1039
1040 void ContactModelLinearCBond::objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *ret) const {
1041 ret->clear();
1042 ret->push_back({"kn",ScalarInfo});
1043 ret->push_back({"ks",ScalarInfo});
1044 ret->push_back({"fric",ScalarInfo});
1045 ret->push_back({"lin_force",VectorInfo});
1046 ret->push_back({"lin_slip",ScalarInfo});
1047 ret->push_back({"lin_mode",ScalarInfo});
1048 ret->push_back({"rgap",ScalarInfo});
1049 ret->push_back({"emod",ScalarInfo});
1050 ret->push_back({"kratio",ScalarInfo});
1051 ret->push_back({"dp_nratio",ScalarInfo});
1052 ret->push_back({"dp_sratio",ScalarInfo});
1053 ret->push_back({"dp_mode",ScalarInfo});
1054 ret->push_back({"dp_force",VectorInfo});
1055 ret->push_back({"cb_state",ScalarInfo});
1056 ret->push_back({"cb_tenf",ScalarInfo});
1057 ret->push_back({"cb_shearf",ScalarInfo});
1058 ret->push_back({"cb_tens",ScalarInfo});
1059 ret->push_back({"cb_shears",ScalarInfo});
1060 ret->push_back({"user_area",ScalarInfo});
1061 }
1062
1063 void ContactModelLinearCBond::objectPropValues(std::vector<double> *ret,const IContact *con) const {
1064 FP_S;
1065 ret->clear();
1066 ret->push_back(kn());
1067 ret->push_back(ks());
1068 ret->push_back(fric());
1069 ret->push_back(lin_F().mag());
1070 ret->push_back(lin_S());
1071 ret->push_back(lin_mode());
1072 ret->push_back(rgap());
1073 double d;
1074 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
1075 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
1076 double rsum = calcRSum(c);
1077 if (userArea_) {
1078#ifdef THREED
1079 rsq = std::sqrt(userArea_ / dPi);
1080#else
1081 rsq = userArea_ / 2.0;
1082#endif
1083 rsum = rsq + rsq;
1084 rsq = 1. / rsq;
1085 }
1086#ifdef TWOD
1087 d= (kn_ * rsum * rsq / 2.0);
1088#else
1089 d= (kn_ * rsum * rsq * rsq) / dPi;
1090#endif
1091 ret->push_back(d);
1092 ret->push_back((ks_ == 0.0) ? 0.0 : (kn_/ks_));
1093 ret->push_back(dp_nratio());
1094 ret->push_back(dp_sratio());
1095 ret->push_back(dp_mode());
1096 ret->push_back(dp_F().mag());
1097 ret->push_back(cb_state_);
1098 ret->push_back(cb_tenF_);
1099 ret->push_back(cb_shearF_);
1100 double tmp(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
1101 if (userArea_) {
1102#ifdef THREED
1103 tmp = std::sqrt(userArea_ / dPi);
1104#else
1105 tmp = userArea_ / 2.0;
1106#endif
1107 tmp = 1. / tmp;
1108 }
1109#ifdef TWOD
1110 ret->push_back(cb_tenF_ * tmp / 2.0);
1111 ret->push_back(cb_shearF_ * tmp / 2.0);
1112#else
1113 ret->push_back(cb_tenF_ * tmp * tmp / dPi);
1114 ret->push_back(cb_shearF_ * tmp * tmp / dPi);
1115#endif
1116 ret->push_back(getArea());
1117 FP_S;
1118 }
1119
1120 void ContactModelLinearCBond::setDampCoefficients(const double &mass,double *vcn,double *vcs) {
1121 *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(kn_));
1122 *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(ks_));
1123 }
1124
1125} // namespace itascaxd
1126// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Jun 07, 2025 |
