Linear Parallel Bond Model Implementation
See this page for the documentation of this contact model.
contactmodellinearpbond.h
1#pragma once
2// contactmodellinearpbond.h
3
4#include "contactmodel/src/contactmodelmechanical.h"
5
6#ifdef LINEARPBOND_LIB
7# define LINEARPBOND_EXPORT EXPORT_TAG
8#elif defined(NO_MODEL_IMPORT)
9# define LINEARPBOND_EXPORT
10#else
11# define LINEARPBOND_EXPORT IMPORT_TAG
12#endif
13
14namespace cmodelsxd {
15 using namespace itasca;
16
17 class ContactModelLinearPBond : public ContactModelMechanical {
18 public:
19 LINEARPBOND_EXPORT ContactModelLinearPBond();
20 LINEARPBOND_EXPORT ContactModelLinearPBond(const ContactModelLinearPBond &) noexcept;
21 LINEARPBOND_EXPORT const ContactModelLinearPBond & operator=(const ContactModelLinearPBond &);
22 LINEARPBOND_EXPORT void addToStorage(poly::vector<ContactModelMechanical> *s,int ww = -1) override;
23 LINEARPBOND_EXPORT virtual ~ContactModelLinearPBond();
24 void copy(const ContactModel *c) override;
25 void archive(ArchiveStream &) override;
26 QString getName() const override { return "linearpbond"; }
27 void setIndex(int i) override { index_=i;}
28 int getIndex() const override {return index_;}
29
30 enum PropertyKeys {
31 kwLinKn=1
32 , kwLinKs
33 , kwLinFric
34 , kwLinF
35 , kwLinS
36 , kwLinMode
37 , kwRGap
38 , kwEmod
39 , kwKRatio
40 , kwDpNRatio
41 , kwDpSRatio
42 , kwDpMode
43 , kwDpF
44 , kwPbState
45 , kwPbRMul
46 , kwPbKn
47 , kwPbKs
48 , kwPbMcf
49 , kwPbTStrength
50 , kwPbSStrength
51 , kwPbCoh
52 , kwPbFa
53 , kwPbSig
54 , kwPbTau
55 , kwPbF
56 , kwPbM
57 , kwPbRadius
58 , kwPbEmod
59 , kwPbKRatio
60 , kwUserArea
61 };
62
63 QString getProperties() const override {
64 return "kn"
65 ",ks"
66 ",fric"
67 ",lin_force"
68 ",lin_slip"
69 ",lin_mode"
70 ",rgap"
71 ",emod"
72 ",kratio"
73 ",dp_nratio"
74 ",dp_sratio"
75 ",dp_mode"
76 ",dp_force"
77 ",pb_state"
78 ",pb_rmul"
79 ",pb_kn"
80 ",pb_ks"
81 ",pb_mcf"
82 ",pb_ten"
83 ",pb_shear"
84 ",pb_coh"
85 ",pb_fa"
86 ",pb_sigma"
87 ",pb_tau"
88 ",pb_force"
89 ",pb_moment"
90 ",pb_radius"
91 ",pb_emod"
92 ",pb_kratio"
93 ",user_area";
94 }
95
96 enum EnergyKeys { kwEStrain=1,kwESlip,kwEDashpot,kwEPbStrain};
97 QString getEnergies() const override { return "energy-strain,energy-slip,energy-dashpot,energy-pbstrain";}
98 double getEnergy(uint32 i) const override; // Base 1
99 bool getEnergyAccumulate(uint32 i) const override; // Base 1
100 void setEnergy(uint32 i,const double &d) override; // Base 1
101 void activateEnergy() override { if (energies_) return; energies_ = NEWC(Energies());}
102 bool getEnergyActivated() const override {return (energies_ !=0);}
103
104 enum FishCallEvents {fActivated=0,fBondBreak,fSlipChange};
105 QString getFishCallEvents() const override { return "contact_activated,bond_break,slip_change"; }
106 QVariant getProperty(uint32 i,const IContact *con=0) const override;
107 bool getPropertyGlobal(uint32 i) const override;
108 bool setProperty(uint32 i,const QVariant &v,IContact *con=0);
109 bool getPropertyReadOnly(uint32 i) const override;
110 bool supportsInheritance(uint32 i) const override;
111 bool getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i); return (inheritanceField_ & mask) ? true : false; }
112 void setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i); if (b) inheritanceField_ |= mask; else inheritanceField_ &= ~mask; }
113
114 enum MethodKeys { kwDeformability=1
115 , kwPbDeformability
116 , kwPbBond
117 , kwPbUnbond
118 , kwArea
119 };
120
121 QString getMethods() const override {
122 return "deformability"
123 ",pb_deformability"
124 ",bond"
125 ",unbond"
126 ",area";
127 }
128
129 QString getMethodArguments(uint32 i) const override;
130
131 bool setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override; // Base 1 - returns true if timestep contributions need to be updated
132
133 uint32 getMinorVersion() const override;
134
135 bool validate(ContactModelMechanicalState *state,const double ×tep) override;
136 bool endPropertyUpdated(const QString &name,const IContactMechanical *c) override;
137 bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) override;
138 DVect2 getEffectiveTranslationalStiffness() const override { DVect2 ret = effectiveTranslationalStiffness_; if(pbProps_) ret+= pbProps_->pbTransStiff_ ;return ret;}
139 DAVect getEffectiveRotationalStiffness() const override {if (!pbProps_) return DAVect(0.0); return pbProps_->pbAngStiff_;}
140
141 bool thermalCoupling(ContactModelMechanicalState *, ContactModelThermalState * , IContactThermal *,const double &) override;
142
143 ContactModelLinearPBond *clone() const override { return NEWC(ContactModelLinearPBond()); }
144 double getActivityDistance() const override {return rgap_;}
145 bool isOKToDelete() const override { return !isBonded(); }
146 void resetForcesAndMoments() override { lin_F(DVect(0.0)); dp_F(DVect(0.0)); pbF(DVect(0.0)); pbM(DAVect(0.0)); if (energies_) { energies_->estrain_ = 0.0; if (energies_) energies_->epbstrain_ = 0.0;}}
147 void setForce(const DVect &v,IContact *c) override;
148 void setArea(const double &d) override { userArea_ = d; }
149 double getArea() const override { return userArea_; }
150 bool checkActivity(const double &gap) override { return (gap <= rgap_ || isBonded()); }
151 bool isSliding() const override { return lin_S_; }
152 bool isBonded() const override { return pbProps_ ? (pbProps_->pb_state_==3) : false; }
153 void unbond() override { if (pbProps_) pbProps_->pb_state_= 0; }
154 void propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
155 void setNonForcePropsFrom(IContactModel *oldCM) override;
156 /// Return the total force that the contact model holds.
157 DVect getForce() const override;
158 /// Return the total moment on 1 that the contact model holds
159 DAVect getMomentOn1(const IContactMechanical*) const override;
160 /// Return the total moment on 1 that the contact model holds
161 DAVect getMomentOn2(const IContactMechanical*) const override;
162
163 DAVect getModelMomentOn1() const override;
164 DAVect getModelMomentOn2() const override;
165
166 // Used to efficiently get properties from the contact model for the object pane.
167 // List of properties for the object pane, comma separated.
168 // All properties will be cast to doubles for comparison. No other comparisons
169 // are supported. This may not be the same as the entire property list.
170 // Return property name and type for plotting.
171 void objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *) const override;
172 // All properties cast to doubles - this is what can be compared.
173 void objectPropValues(std::vector<double> *,const IContact *) const override;
174
175 const double & kn() const {return kn_;}
176 void kn(const double &d) {kn_=d;}
177 const double & ks() const {return ks_;}
178 void ks(const double &d) {ks_=d;}
179 const double & fric() const {return fric_;}
180 void fric(const double &d) {fric_=d;}
181 const DVect & lin_F() const {return lin_F_;}
182 void lin_F(const DVect &f) { lin_F_=f;}
183 bool lin_S() const {return lin_S_;}
184 void lin_S(bool b) { lin_S_=b;}
185 uint32 lin_mode() const {return lin_mode_;}
186 void lin_mode(uint32 i) { lin_mode_=i;}
187 const double & rgap() const {return rgap_;}
188 void rgap(const double &d) {rgap_=d;}
189
190 bool hasDamping() const {return dpProps_ ? true : false;}
191 double dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
192 void dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
193 double dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
194 void dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
195 int dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
196 void dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
197 DVect dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
198 void dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
199
200 bool hasEnergies() const {return energies_ ? true:false;}
201 double estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
202 void estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
203 double eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
204 void eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
205 double edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
206 void edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
207 double epbstrain() const {return hasEnergies() ? energies_->epbstrain_: 0.0;}
208 void epbstrain(const double &d) { if(!hasEnergies()) return; energies_->epbstrain_=d;}
209
210 bool hasPBond() const {return pbProps_ ? true:false;}
211 int pbState() const {return hasPBond() ? pbProps_->pb_state_: 0;}
212 void pbState(int i) { if(!hasPBond()) return; pbProps_->pb_state_=i;}
213 double pbRmul() const {return (hasPBond() ? (pbProps_->pb_rmul_) : 0.0);}
214 void pbRmul(const double &d) { if(!hasPBond()) return; pbProps_->pb_rmul_=d;}
215 double pbKn() const {return (hasPBond() ? (pbProps_->pb_kn_) : 0.0);}
216 void pbKn(const double &d) { if(!hasPBond()) return; pbProps_->pb_kn_=d;}
217 double pbKs() const {return (hasPBond() ? (pbProps_->pb_ks_) : 0.0);}
218 void pbKs(const double &d) { if(!hasPBond()) return; pbProps_->pb_ks_=d;}
219 double pbMCF() const {return (hasPBond() ? (pbProps_->pb_mcf_) : 0.0);}
220 void pbMCF(const double &d) { if(!hasPBond()) return; pbProps_->pb_mcf_=d;}
221 double pbTen() const {return (hasPBond() ? (pbProps_->pb_ten_) : 0.0);}
222 void pbTen(const double &d) { if(!hasPBond()) return; pbProps_->pb_ten_=d;}
223 double pbCoh() const {return (hasPBond() ? (pbProps_->pb_coh_) : 0.0);}
224 void pbCoh(const double &d) { if(!hasPBond()) return; pbProps_->pb_coh_=d;}
225 double pbFA() const {return (hasPBond() ? (pbProps_->pb_fa_) : 0.0);}
226 void pbFA(const double &d) { if(!hasPBond()) return; pbProps_->pb_fa_=d;}
227 DVect pbF() const {return hasPBond() ? pbProps_->pb_F_: DVect(0.0);}
228 void pbF(const DVect &f) { if(!hasPBond()) return; pbProps_->pb_F_=f;}
229 DAVect pbM() const {return hasPBond() ? pbProps_->pb_M_: DAVect(0.0);}
230 void pbM(const DAVect &m) { if(!hasPBond()) return; pbProps_->pb_M_=m;}
231 DVect2 pbTransStiff() const {return hasPBond() ? pbProps_->pbTransStiff_: DVect2(0.0);}
232 void pbTransStiff(const DVect2 &f) { if(!hasPBond()) return; pbProps_->pbTransStiff_=f;}
233 DAVect pbAngStiff() const {return hasPBond() ? pbProps_->pbAngStiff_: DAVect(0.0);}
234 void pbAngStiff(const DAVect &m) { if(!hasPBond()) return; pbProps_->pbAngStiff_=m;}
235
236 uint32 inheritanceField() const {return inheritanceField_;}
237 void inheritanceField(uint32 i) {inheritanceField_ = i;}
238
239 const DVect2 & effectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
240 void effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
241
242 private:
243 static int index_;
244
245 struct Energies {
246 Energies() : estrain_(0.0), eslip_(0.0), edashpot_(0.0), epbstrain_(0.0) {}
247 double estrain_; // elastic energy stored in contact
248 double eslip_; // work dissipated by friction
249 double edashpot_; // work dissipated by dashpots
250 double epbstrain_; // parallel bond strain energy
251 };
252
253 struct dpProps {
254 dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)) {}
255 double dp_nratio_; // normal viscous critical damping ratio
256 double dp_sratio_; // shear viscous critical damping ratio
257 int dp_mode_; // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
258 DVect dp_F_; // Force in the dashpots
259 };
260
261 struct pbProps {
262 pbProps() : pb_state_(0), pb_rmul_(1.0), pb_kn_(0.0), pb_ks_(0.0),
263 pb_mcf_(1.0), pb_ten_(0.0), pb_coh_(0.0), pb_fa_(0.0), pb_F_(DVect(0.0)), pb_M_(DAVect(0.0)),
264 pbTransStiff_(0.0), pbAngStiff_(0.0){}
265 // parallel bond
266 int pb_state_; // Bond mode - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (B)
267 double pb_rmul_; // Radius multiplier
268 double pb_kn_; // normal stiffness
269 double pb_ks_; // shear stiffness
270 double pb_mcf_; // Moment contribution factor
271 double pb_ten_; // normal strength
272 double pb_coh_; // cohesion
273 double pb_fa_; // friction angle
274 DVect pb_F_; // Force in parallel bond
275 DAVect pb_M_; // moment in parallel bond
276 DVect2 pbTransStiff_; // (Normal,Shear) Translational stiffness of the parallel bond
277 DAVect pbAngStiff_; // (Normal,Shear) Rotational stiffness of the parallel bond
278 };
279
280 bool updateKn(const IContactMechanical *con);
281 bool updateKs(const IContactMechanical *con);
282 bool updateFric(const IContactMechanical *con);
283 double pbStrainEnergy() const; // Compute bond strain energy
284
285 void updateEffectiveStiffness(ContactModelMechanicalState *state);
286
287 DVect3 pbData(const IContactMechanical *con) const; // Bond area and inertia
288 DVect2 pbSMax(const IContactMechanical *con) const; // Maximum stress (tensile,shear) at bond periphery
289 double pbShearStrength(const double &pbArea) const; // Bond shear strength
290 void setDampCoefficients(const double &mass,double *vcn,double *vcs);
291
292 // inheritance fields
293 uint32 inheritanceField_;
294
295 // linear model
296 double kn_ = 0; // normal stiffness
297 double ks_ = 0; // shear stiffness
298 double fric_ = 0; // Coulomb friction coefficient
299 DVect lin_F_ = DVect(0.0); // Force carried in the linear model
300 bool lin_S_ = false; // the current sliding state
301 uint32 lin_mode_ = 0; // specifies absolute (0) or incremental (1) behavior for the the linear part
302 double rgap_ = 0; // reference gap for the linear part
303
304 dpProps * dpProps_ = nullptr; // The viscous properties
305 pbProps * pbProps_ = nullptr; // The parallel bond properties
306
307 double userArea_ = 0; // User specified area
308
309 Energies * energies_ = nullptr; // energies
310
311 DVect2 effectiveTranslationalStiffness_ = DVect2(0.0);
312
313 };
314} // namespace itascaxd
315// EoF
contactmodellinearpbond.cpp
1// contactmodellinear.cpp
2#include "contactmodellinearpbond.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 LINEARPBOND_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 "contactmodelmechanical3dlinearpbond";
29#else
30 return "contactmodelmechanical2dlinearpbond";
31#endif
32 }
33
34 extern "C" EXPORT_TAG unsigned getMajorVersion()
35 {
36 return MAJOR_VERSION;
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::ContactModelLinearPBond *m = NEWC(cmodelsxd::ContactModelLinearPBond());
47 return (void *)m;
48 }
49#endif // LINEARPBOND_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 ContactModelLinearPBond::index_ = -1;
59 uint32 ContactModelLinearPBond::getMinorVersion() const { return MINOR_VERSION;}
60
61 ContactModelLinearPBond::ContactModelLinearPBond() : inheritanceField_(linKnMask|linKsMask|linFricMask)
62 { }
63
64 ContactModelLinearPBond::ContactModelLinearPBond(const ContactModelLinearPBond &m) noexcept {
65 inheritanceField(linKnMask|linKsMask|linFricMask);
66 this->copy(&m);
67 }
68
69 const ContactModelLinearPBond & ContactModelLinearPBond::operator=(const ContactModelLinearPBond &m) {
70 inheritanceField(linKnMask|linKsMask|linFricMask);
71 this->copy(&m);
72 return *this;
73 }
74
75 void ContactModelLinearPBond::addToStorage(poly::vector<ContactModelMechanical> *s,int ww) {
76 s->addToStorage<ContactModelLinearPBond>(*this,ww);
77 }
78
79 ContactModelLinearPBond::~ContactModelLinearPBond() {
80 if (dpProps_)
81 delete dpProps_;
82 if (pbProps_)
83 delete pbProps_;
84 if (energies_)
85 delete energies_;
86 }
87
88 void ContactModelLinearPBond::archive(ArchiveStream &stream) {
89 stream & kn_;
90 stream & ks_;
91 stream & fric_;
92 stream & lin_F_;
93 stream & lin_S_;
94 stream & lin_mode_;
95 if (stream.getArchiveState()==ArchiveStream::Save) {
96 bool b = false;
97 if (dpProps_) {
98 b = true;
99 stream & b;
100 stream & dpProps_->dp_nratio_;
101 stream & dpProps_->dp_sratio_;
102 stream & dpProps_->dp_mode_;
103 stream & dpProps_->dp_F_;
104 }
105 else
106 stream & b;
107
108 b = false;
109 if (energies_) {
110 b = true;
111 stream & b;
112 stream & energies_->estrain_;
113 stream & energies_->eslip_;
114 stream & energies_->edashpot_;
115 stream & energies_->epbstrain_;
116 }
117 else
118 stream & b;
119
120 b = false;
121 if (pbProps_) {
122 b = true;
123 stream & b;
124 stream & pbProps_->pb_state_;
125 stream & pbProps_->pb_rmul_;
126 stream & pbProps_->pb_kn_;
127 stream & pbProps_->pb_ks_;
128 stream & pbProps_->pb_mcf_;
129 stream & pbProps_->pb_ten_;
130 stream & pbProps_->pb_coh_;
131 stream & pbProps_->pb_fa_;
132 stream & pbProps_->pb_F_;
133 stream & pbProps_->pb_M_;
134 }
135 else
136 stream & b;
137 } else {
138 bool b(false);
139 stream & b;
140 if (b) {
141 if (!dpProps_)
142 dpProps_ = NEWC(dpProps());
143 stream & dpProps_->dp_nratio_;
144 stream & dpProps_->dp_sratio_;
145 stream & dpProps_->dp_mode_;
146 stream & dpProps_->dp_F_;
147 }
148 stream & b;
149 if (b) {
150 if (!energies_)
151 energies_ = NEWC(Energies());
152 stream & energies_->estrain_;
153 stream & energies_->eslip_;
154 stream & energies_->edashpot_;
155 stream & energies_->epbstrain_;
156 }
157 stream & b;
158 if (b) {
159 if (!pbProps_)
160 pbProps_ = NEWC(pbProps());
161 stream & pbProps_->pb_state_;
162 stream & pbProps_->pb_rmul_;
163 stream & pbProps_->pb_kn_;
164 stream & pbProps_->pb_ks_;
165 stream & pbProps_->pb_mcf_;
166 stream & pbProps_->pb_ten_;
167 stream & pbProps_->pb_coh_;
168 stream & pbProps_->pb_fa_;
169 stream & pbProps_->pb_F_;
170 stream & pbProps_->pb_M_;
171 }
172 }
173
174 stream & inheritanceField_;
175 stream & effectiveTranslationalStiffness_;
176
177 if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() == getMinorVersion())
178 stream & rgap_;
179
180 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 1)
181 stream & userArea_;
182 }
183
184 void ContactModelLinearPBond::copy(const ContactModel *cm) {
185 ContactModelMechanical::copy(cm);
186 const ContactModelLinearPBond *in = dynamic_cast<const ContactModelLinearPBond*>(cm);
187 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
188 kn(in->kn());
189 ks(in->ks());
190 fric(in->fric());
191 lin_F(in->lin_F());
192 lin_S(in->lin_S());
193 lin_mode(in->lin_mode());
194 rgap(in->rgap());
195 if (in->hasDamping()) {
196 if (!dpProps_)
197 dpProps_ = NEWC(dpProps());
198 dp_nratio(in->dp_nratio());
199 dp_sratio(in->dp_sratio());
200 dp_mode(in->dp_mode());
201 dp_F(in->dp_F());
202 }
203 if (in->hasEnergies()) {
204 if (!energies_)
205 energies_ = NEWC(Energies());
206 estrain(in->estrain());
207 eslip(in->eslip());
208 edashpot(in->edashpot());
209 epbstrain(in->epbstrain());
210 }
211 if (in->hasPBond()) {
212 if (!pbProps_)
213 pbProps_ = NEWC(pbProps());
214 pbState(in->pbState());
215 pbRmul(in->pbRmul());
216 pbKn(in->pbKn());
217 pbKs(in->pbKs());
218 pbMCF(in->pbMCF());
219 pbTen(in->pbTen());
220 pbCoh(in->pbCoh());
221 pbFA(in->pbFA());
222 pbF(in->pbF());
223 pbM(in->pbM());
224 pbTransStiff(in->pbTransStiff());
225 pbAngStiff(in->pbAngStiff());
226 }
227 userArea_ = in->userArea_;
228 inheritanceField(in->inheritanceField());
229 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
230 }
231
232 QVariant ContactModelLinearPBond::getProperty(uint32 i,const IContact *con) const {
233 // The IContact pointer may be a nullptr!
234 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
235 QVariant var;
236 switch (i) {
237 case kwLinKn: return kn_;
238 case kwLinKs: return ks_;
239 case kwLinFric: return fric_;
240 case kwLinMode: return lin_mode_;
241 case kwLinF: var.setValue(lin_F_); return var;
242 case kwLinS: return lin_S_;
243 case kwRGap: return rgap_;
244 case kwEmod: {
245 if (c ==nullptr) return 0.0;
246 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
247 double rsum(0.0);
248 if (c->getEnd1Curvature().y())
249 rsum += 1.0/c->getEnd1Curvature().y();
250 if (c->getEnd2Curvature().y())
251 rsum += 1.0/c->getEnd2Curvature().y();
252 if (userArea_) {
253#ifdef THREED
254 rsq = std::sqrt(userArea_ / dPi);
255#else
256 rsq = userArea_ / 2.0;
257#endif
258 rsum = rsq + rsq;
259 rsq = 1. / rsq;
260 }
261#ifdef TWOD
262 return (kn_ * rsum * rsq / 2.0);
263#else
264 return (kn_ * rsum * rsq * rsq) / dPi;
265#endif
266 }
267 case kwKRatio: return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
268 case kwDpNRatio: return dpProps_ ? dpProps_->dp_nratio_ : 0;
269 case kwDpSRatio: return dpProps_ ? dpProps_->dp_sratio_ : 0;
270 case kwDpMode: return dpProps_ ? dpProps_->dp_mode_ : 0;
271 case kwUserArea: return userArea_;
272 case kwDpF: {
273 dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
274 return var;
275 }
276 case kwPbState: return pbProps_ ? pbProps_->pb_state_ : 0;
277 case kwPbRMul: return pbProps_ ? pbProps_->pb_rmul_ : 1.0;
278 case kwPbKn: return pbProps_ ? pbProps_->pb_kn_ : 0;
279 case kwPbKs: return pbProps_ ? pbProps_->pb_ks_ : 0;
280 case kwPbMcf: return pbProps_ ? pbProps_->pb_mcf_ : 1.0;
281 case kwPbTStrength: return pbProps_ ? pbProps_->pb_ten_ : 0.0;
282 case kwPbSStrength: {
283 if (!pbProps_ or not c) return 0.0;
284 double pbArea = pbData(c).x();
285 return pbShearStrength(pbArea);
286 }
287 case kwPbCoh: return pbProps_ ? pbProps_->pb_coh_ : 0;
288 case kwPbFa: return pbProps_ ? pbProps_->pb_fa_ : 0;
289 case kwPbSig: {
290 if (!pbProps_ or pbProps_->pb_state_ < 3 or not c) return 0.0;
291 return pbSMax(c).x();
292 }
293 case kwPbTau: {
294 if (!pbProps_ or pbProps_->pb_state_ < 3 or not c) return 0.0;
295 return pbSMax(c).y();
296 }
297 case kwPbF: {
298 pbProps_ ? var.setValue(pbProps_->pb_F_) : var.setValue(DVect(0.0));
299 return var;
300 }
301 case kwPbM: {
302 pbProps_ ? var.setValue(pbProps_->pb_M_) : var.setValue(DAVect(0.0));
303 return var;
304 }
305 case kwPbRadius: {
306 if (!pbProps_ or not c) return 0.0;
307 double Cmax1 = c->getEnd1Curvature().y();
308 double Cmax2 = c->getEnd2Curvature().y();
309 double br = pbProps_->pb_rmul_ * 1.0 / std::max(Cmax1,Cmax2);
310 if (userArea_)
311#ifdef THREED
312 br = std::sqrt(userArea_ / dPi);
313#else
314 br = userArea_ / 2.0;
315#endif
316 return br;
317 }
318 case kwPbEmod: {
319 if (!pbProps_ or not c) return 0.0;
320 double rsum(0.0);
321 if (c->getEnd1Curvature().y())
322 rsum += 1.0/c->getEnd1Curvature().y();
323 if (c->getEnd2Curvature().y())
324 rsum += 1.0/c->getEnd2Curvature().y();
325 if (userArea_) {
326#ifdef THREED
327 double rad = std::sqrt(userArea_ / dPi);
328#else
329 double rad = userArea_ / 2.0;
330#endif
331 rsum = 2 * rad;
332 }
333 double emod = pbProps_->pb_kn_ * rsum;
334 return emod;
335 }
336 case kwPbKRatio: {
337 if (!pbProps_) return 0.0;
338 return (pbProps_->pb_ks_ == 0.0) ? 0.0 : (pbProps_->pb_kn_/pbProps_->pb_ks_);
339 }
340 }
341 assert(0);
342 return QVariant();
343 }
344
345 bool ContactModelLinearPBond::getPropertyGlobal(uint32 i) const {
346 switch (i) {
347 case kwLinF:
348 case kwDpF:
349 case kwPbF:
350 return false;
351 }
352 return true;
353 }
354
355 bool ContactModelLinearPBond::setProperty(uint32 i,const QVariant &v,IContact *) {
356 pbProps pb;
357 dpProps dp;
358
359 switch (i) {
360 case kwLinKn: {
361 if (!v.canConvert<double>())
362 throw Exception("kn must be a double.");
363 double val(v.toDouble());
364 if (val<0.0)
365 throw Exception("Negative kn not allowed.");
366 kn_ = val;
367 return true;
368 }
369 case kwLinKs: {
370 if (!v.canConvert<double>())
371 throw Exception("ks must be a double.");
372 double val(v.toDouble());
373 if (val<0.0)
374 throw Exception("Negative ks not allowed.");
375 ks_ = val;
376 return true;
377 }
378 case kwLinFric: {
379 if (!v.canConvert<double>())
380 throw Exception("fric must be a double.");
381 double val(v.toDouble());
382 if (val<0.0)
383 throw Exception("Negative fric not allowed.");
384 fric_ = val;
385 return false;
386 }
387 case kwLinF: {
388 if (!v.canConvert<DVect>())
389 throw Exception("lin_force must be a vector.");
390 DVect val(v.value<DVect>());
391 lin_F_ = val;
392 return false;
393 }
394 case kwLinMode: {
395 if (!v.canConvert<uint32>())
396 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
397 uint32 val(v.toUInt());
398 if (val>1)
399 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
400 lin_mode_ = val;
401 return false;
402 }
403 case kwRGap: {
404 if (!v.canConvert<double>())
405 throw Exception("Reference gap must be a double.");
406 double val(v.toDouble());
407 rgap_ = val;
408 return false;
409 }
410 case kwPbRMul: {
411 if (!v.canConvert<double>())
412 throw Exception("pb_rmul must be a double.");
413 double val(v.toDouble());
414 if (val<=0.0)
415 throw Exception("pb_rmul must be positive.");
416 if (val == 1.0 && !pbProps_)
417 return false;
418 if (!pbProps_)
419 pbProps_ = NEWC(pbProps());
420 pbProps_->pb_rmul_ = val;
421 return false;
422 }
423 case kwPbKn: {
424 if (!v.canConvert<double>())
425 throw Exception("pb_kn must be a double.");
426 double val(v.toDouble());
427 if (val<0.0)
428 throw Exception("Negative pb_kn not allowed.");
429 if (val == 0.0 && !pbProps_)
430 return false;
431 if (!pbProps_)
432 pbProps_ = NEWC(pbProps());
433 pbProps_->pb_kn_ = val;
434 return true;
435 }
436 case kwPbKs: {
437 if (!v.canConvert<double>())
438 throw Exception("pb_ks must be a double.");
439 double val(v.toDouble());
440 if (val<0.0)
441 throw Exception("Negative pb_ks not allowed.");
442 if (val == 0.0 && !pbProps_)
443 return false;
444 if (!pbProps_)
445 pbProps_ = NEWC(pbProps());
446 pbProps_->pb_ks_ = val;
447 return true;
448 }
449 case kwPbMcf: {
450 if (!v.canConvert<double>())
451 throw Exception("pb_mcf must be a double.");
452 double val(v.toDouble());
453 if (val<0.0)
454 throw Exception("Negative pb_mcf not allowed.");
455 if (val > 1.0)
456 throw Exception("pb_mcf must be lower or equal to 1.0.");
457 if (val == 1.0 && !pbProps_)
458 return false;
459 if (!pbProps_)
460 pbProps_ = NEWC(pbProps());
461 pbProps_->pb_mcf_ = val;
462 return false;
463 }
464 case kwPbTStrength: {
465 if (!v.canConvert<double>())
466 throw Exception("pb_ten must be a double.");
467 double val(v.toDouble());
468 if (val < 0.0)
469 throw Exception("Negative pb_ten not allowed.");
470 if (val == 0.0 && !pbProps_)
471 return false;
472 if (!pbProps_)
473 pbProps_ = NEWC(pbProps());
474 pbProps_->pb_ten_ = val;
475 return false;
476 }
477 case kwPbCoh: {
478 if (!v.canConvert<double>())
479 throw Exception("pb_coh must be a double.");
480 double val(v.toDouble());
481 if (val<0.0)
482 throw Exception("Negative pb_coh not allowed.");
483 if (val == 0.0 && !pbProps_)
484 return false;
485 if (!pbProps_)
486 pbProps_ = NEWC(pbProps());
487 pbProps_->pb_coh_ = val;
488 return false;
489 }
490 case kwPbFa: {
491 if (!v.canConvert<double>())
492 throw Exception("pb_fa must be a double.");
493 double val(v.toDouble());
494 if (val<0.0)
495 throw Exception("Negative pb_fa not allowed.");
496 if (val>=90.0)
497 throw Exception("pb_fa must be lower than 90.0 degrees.");
498 if (val == 0.0 && !pbProps_)
499 return false;
500 if (!pbProps_)
501 pbProps_ = NEWC(pbProps());
502 pbProps_->pb_fa_ = val;
503 return false;
504 }
505 case kwPbF: {
506 if (!v.canConvert<DVect>())
507 throw Exception("pb_force must be a vector.");
508 DVect val(v.value<DVect>());
509 if (val.fsum() == 0.0 && !pbProps_)
510 return false;
511 if (!pbProps_)
512 pbProps_ = NEWC(pbProps());
513 pbProps_->pb_F_ = val;
514 return false;
515 }
516 case kwPbM: {
517 DAVect val(0.0);
518#ifdef TWOD
519 if (!v.canConvert<DAVect>() && !v.canConvert<double>())
520 throw Exception("pb_moment must be an angular vector.");
521 if (v.canConvert<DAVect>())
522 val = DAVect(v.value<DAVect>());
523 else
524 val = DAVect(v.toDouble());
525#else
526 if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
527 throw Exception("pb_moment must be an angular vector.");
528 if (v.canConvert<DAVect>())
529 val = DAVect(v.value<DAVect>());
530 else
531 val = DAVect(v.value<DVect>());
532#endif
533 if (val.fsum() == 0.0 && !pbProps_)
534 return false;
535 if (!pbProps_)
536 pbProps_ = NEWC(pbProps());
537 pbProps_->pb_M_ = val;
538 return false;
539 }
540 case kwDpNRatio: {
541 if (!v.canConvert<double>())
542 throw Exception("dp_nratio must be a double.");
543 double val(v.toDouble());
544 if (val<0.0)
545 throw Exception("Negative dp_nratio not allowed.");
546 if (val == 0.0 && !dpProps_)
547 return false;
548 if (!dpProps_)
549 dpProps_ = NEWC(dpProps());
550 dpProps_->dp_nratio_ = val;
551 return true;
552 }
553 case kwDpSRatio: {
554 if (!v.canConvert<double>())
555 throw Exception("dp_sratio must be a double.");
556 double val(v.toDouble());
557 if (val<0.0)
558 throw Exception("Negative dp_sratio not allowed.");
559 if (val == 0.0 && !dpProps_)
560 return false;
561 if (!dpProps_)
562 dpProps_ = NEWC(dpProps());
563 dpProps_->dp_sratio_ = val;
564 return true;
565 }
566 case kwDpMode: {
567 if (!v.canConvert<int>())
568 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
569 int val(v.toInt());
570 if (val == 0 && !dpProps_)
571 return false;
572 if (val < 0 || val > 3)
573 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
574 if (!dpProps_)
575 dpProps_ = NEWC(dpProps());
576 dpProps_->dp_mode_ = val;
577 return false;
578 }
579 case kwDpF: {
580 if (!v.canConvert<DVect>())
581 throw Exception("dp_force must be a vector.");
582 DVect val(v.value<DVect>());
583 if (val.fsum() == 0.0 && !dpProps_)
584 return false;
585 if (!dpProps_)
586 dpProps_ = NEWC(dpProps());
587 dpProps_->dp_F_ = val;
588 return false;
589 }
590 case kwUserArea: {
591 if (!v.canConvert<double>())
592 throw Exception("user_area must be a double.");
593 double val(v.toDouble());
594 if (val < 0.0)
595 throw Exception("Negative user_area not allowed.");
596 userArea_ = val;
597 return true;
598 }
599 }
600// assert(0);
601 return false;
602 }
603
604 bool ContactModelLinearPBond::getPropertyReadOnly(uint32 i) const {
605 switch (i) {
606 case kwDpF:
607 case kwLinS:
608 case kwEmod:
609 case kwKRatio:
610 case kwPbState:
611 case kwPbRadius:
612 case kwPbSStrength:
613 case kwPbSig:
614 case kwPbTau:
615 case kwPbEmod:
616 case kwPbKRatio:
617 return true;
618 default:
619 break;
620 }
621 return false;
622 }
623
624 bool ContactModelLinearPBond::supportsInheritance(uint32 i) const {
625 switch (i) {
626 case kwLinKn:
627 case kwLinKs:
628 case kwLinFric:
629 return true;
630 default:
631 break;
632 }
633 return false;
634 }
635
636 QString ContactModelLinearPBond::getMethodArguments(uint32 i) const {
637 QString def = QString();
638 switch (i) {
639 case kwDeformability:
640 return "emod,kratio";
641 case kwPbDeformability:
642 return "emod,kratio";
643 case kwPbBond:
644 return "gap";
645 case kwPbUnbond:
646 return "gap";
647 }
648 return def;
649 }
650
651 bool ContactModelLinearPBond::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
652 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
653 switch (i) {
654 case kwDeformability: {
655 double emod;
656 double krat;
657 if (vl.at(0).isNull())
658 throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
659 emod = vl.at(0).toDouble();
660 if (emod<0.0)
661 throw Exception("Negative emod not allowed in contact model %1.",getName());
662 if (vl.at(1).isNull())
663 throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
664 krat = vl.at(1).toDouble();
665 if (krat<0.0)
666 throw Exception("Negative linear stiffness ratio not allowed in contact model %1.",getName());
667 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
668 double rsum(0.0);
669 if (c->getEnd1Curvature().y())
670 rsum += 1.0/c->getEnd1Curvature().y();
671 if (c->getEnd2Curvature().y())
672 rsum += 1.0/c->getEnd2Curvature().y();
673 if (userArea_) {
674#ifdef THREED
675 rsq = std::sqrt(userArea_ / dPi);
676#else
677 rsq = userArea_ / 2.0;
678#endif
679 rsum = rsq + rsq;
680 rsq = 1. / rsq;
681 }
682#ifdef TWOD
683 kn_ = 2.0 * emod / (rsq * rsum);
684#else
685 kn_ = dPi * emod / (rsq * rsq * rsum);
686#endif
687 ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
688 setInheritance(1,false);
689 setInheritance(2,false);
690 return true;
691 }
692 case kwPbDeformability: {
693 //if (!pbProps_ || pbProps_->pb_state_ != 3) return false;
694 double emod;
695 double krat;
696 if (vl.at(0).isNull())
697 throw Exception("Argument emod must be specified with method pb_deformability in contact model %1.",getName());
698 emod = vl.at(0).toDouble();
699 if (emod<0.0)
700 throw Exception("Negative emod not allowed in contact model %1.",getName());
701 if (vl.at(1).isNull())
702 throw Exception("Argument kratio must be specified with method pb_deformability in contact model %1.",getName());
703 krat = vl.at(1).toDouble();
704 if (krat<0.0)
705 throw Exception("Negative parallel bond stiffness ratio not allowed in contact model %1.",getName());
706 double rsum(0.0);
707 if (c->getEnd1Curvature().y())
708 rsum += 1.0/c->getEnd1Curvature().y();
709 if (c->getEnd2Curvature().y())
710 rsum += 1.0/c->getEnd2Curvature().y();
711 if (!pbProps_)
712 pbProps_ = NEWC(pbProps());
713 if (userArea_)
714#ifdef THREED
715 rsum = 2 * std::sqrt(userArea_ / dPi);
716#else
717 rsum = 2 * userArea_ / 2.0;
718#endif
719 pbProps_->pb_kn_ = emod / rsum;
720 pbProps_->pb_ks_ = (krat == 0.0) ? 0.0 : pbProps_->pb_kn_ / krat;
721 return true;
722 }
723 case kwPbBond: {
724 if (pbProps_ && pbProps_->pb_state_ == 3) return false;
725 double mingap = -1.0 * limits<double>::max();
726 double maxgap = 0;
727 if (vl.at(0).canConvert<double>())
728 maxgap = vl.at(0).toDouble();
729 else if (vl.at(0).canConvert<DVect2>()) {
730 DVect2 value = vl.at(0).value<DVect2>();
731 mingap = value.minComp();
732 maxgap = value.maxComp();
733 } else if (!vl.at(0).isNull())
734 throw Exception("gap value %1 not recognized in method bond in contact model %2.",vl.at(0),getName());
735 double gap = c->getGap();
736 if ( gap >= mingap && gap <= maxgap) {
737 if (pbProps_)
738 pbProps_->pb_state_ = 3;
739 else {
740 pbProps_ = NEWC(pbProps());
741 pbProps_->pb_state_ = 3;
742 }
743 return true;
744 }
745 return false;
746 }
747 case kwPbUnbond: {
748 if (!pbProps_ || pbProps_->pb_state_ == 0) return false;
749 double mingap = -1.0 * limits<double>::max();
750 double maxgap = 0;
751 if (vl.at(0).canConvert<double>())
752 maxgap = vl.at(0).toDouble();
753 else if (vl.at(0).canConvert<DVect2>()) {
754 DVect2 value = vl.at(0).value<DVect2>();
755 mingap = value.minComp();
756 maxgap = value.maxComp();
757 }
758 else if (!vl.at(0).isNull())
759 throw Exception("gap value %1 not recognized in method unbond in contact model %2.",vl.at(0),getName());
760 double gap = c->getGap();
761 if ( gap >= mingap && gap <= maxgap) {
762 pbProps_->pb_state_ = 0;
763 return true;
764 }
765 return false;
766 }
767 case kwArea: {
768 if (!userArea_) {
769 double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
770#ifdef THREED
771 userArea_ = rsq * rsq * dPi;
772#else
773 userArea_ = rsq * 2.0;
774#endif
775 }
776 return true;
777 }
778 }
779 return false;
780 }
781
782 double ContactModelLinearPBond::getEnergy(uint32 i) const {
783 double ret(0.0);
784 if (!energies_)
785 return ret;
786 switch (i) {
787 case kwEStrain: return energies_->estrain_;
788 case kwESlip: return energies_->eslip_;
789 case kwEDashpot: return energies_->edashpot_;
790 case kwEPbStrain:return energies_->epbstrain_;
791 }
792 assert(0);
793 return ret;
794 }
795
796 bool ContactModelLinearPBond::getEnergyAccumulate(uint32 i) const {
797 switch (i) {
798 case kwEStrain: return false;
799 case kwESlip: return true;
800 case kwEDashpot: return true;
801 case kwEPbStrain: return false;
802 }
803 assert(0);
804 return false;
805 }
806
807 void ContactModelLinearPBond::setEnergy(uint32 i,const double &d) {
808 if (!energies_) return;
809 switch (i) {
810 case kwEStrain: energies_->estrain_ = d; return;
811 case kwESlip: energies_->eslip_ = d; return;
812 case kwEDashpot: energies_->edashpot_ = d; return;
813 case kwEPbStrain: energies_->epbstrain_= d; return;
814 }
815 assert(0);
816 return;
817 }
818
819 bool ContactModelLinearPBond::validate(ContactModelMechanicalState *state,const double &) {
820 assert(state);
821 const IContactMechanical *c = state->getMechanicalContact();
822 assert(c);
823
824 if (state->trackEnergy_)
825 activateEnergy();
826
827 if (inheritanceField_ & linKnMask)
828 updateKn(c);
829 if (inheritanceField_ & linKsMask)
830 updateKs(c);
831 if (inheritanceField_ & linFricMask)
832 updateFric(c);
833
834 updateEffectiveStiffness(state);
835 return checkActivity(state->gap_);
836 }
837
838 static const QString knstr("kn");
839 bool ContactModelLinearPBond::updateKn(const IContactMechanical *con) {
840 assert(con);
841 QVariant v1 = con->getEnd1()->getProperty(knstr);
842 QVariant v2 = con->getEnd2()->getProperty(knstr);
843 if (!v1.isValid() || !v2.isValid())
844 return false;
845 double kn1 = v1.toDouble();
846 double kn2 = v2.toDouble();
847 double val = kn_;
848 if (kn1 && kn2)
849 kn_ = kn1*kn2/(kn1+kn2);
850 else if (kn1)
851 kn_ = kn1;
852 else if (kn2)
853 kn_ = kn2;
854 return ( (kn_ != val) );
855 }
856
857 static const QString ksstr("ks");
858 bool ContactModelLinearPBond::updateKs(const IContactMechanical *con) {
859 assert(con);
860 QVariant v1 = con->getEnd1()->getProperty(ksstr);
861 QVariant v2 = con->getEnd2()->getProperty(ksstr);
862 if (!v1.isValid() || !v2.isValid())
863 return false;
864 double ks1 = v1.toDouble();
865 double ks2 = v2.toDouble();
866 double val = ks_;
867 if (ks1 && ks2)
868 ks_ = ks1*ks2/(ks1+ks2);
869 else if (ks1)
870 ks_ = ks1;
871 else if (ks2)
872 ks_ = ks2;
873 return ( (ks_ != val) );
874 }
875
876 static const QString fricstr("fric");
877 bool ContactModelLinearPBond::updateFric(const IContactMechanical *con) {
878 assert(con);
879 QVariant v1 = con->getEnd1()->getProperty(fricstr);
880 QVariant v2 = con->getEnd2()->getProperty(fricstr);
881 if (!v1.isValid() || !v2.isValid())
882 return false;
883 double fric1 = std::max(0.0,v1.toDouble());
884 double fric2 = std::max(0.0,v2.toDouble());
885 double val = fric_;
886 fric_ = std::min(fric1,fric2);
887 return ( (fric_ != val) );
888 }
889
890 bool ContactModelLinearPBond::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
891 assert(c);
892 QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
893 QRegExp rx(name,Qt::CaseInsensitive);
894 int idx = availableProperties.indexOf(rx)+1;
895 bool ret=false;
896
897 if (idx<=0)
898 return ret;
899
900 switch(idx) {
901 case kwLinKn: { //kn
902 if (inheritanceField_ & linKnMask)
903 ret = updateKn(c);
904 break;
905 }
906 case kwLinKs: { //ks
907 if (inheritanceField_ & linKsMask)
908 ret =updateKs(c);
909 break;
910 }
911 case kwLinFric: { //fric
912 if (inheritanceField_ & linFricMask)
913 updateFric(c);
914 break;
915 }
916 }
917 return ret;
918 }
919
920 void ContactModelLinearPBond::updateEffectiveStiffness(ContactModelMechanicalState *state) {
921 DVect2 ret(kn_,ks_);
922 // account for viscous damping
923 if (dpProps_) {
924 DVect2 correct(1.0);
925 if (dpProps_->dp_nratio_)
926 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
927 if (dpProps_->dp_sratio_)
928 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
929 ret /= (correct*correct);
930 }
931 effectiveTranslationalStiffness_ = ret;
932 if (isBonded()) {
933 double Cmin1 = state->end1Curvature_.x();
934 double Cmax1 = state->end1Curvature_.y();
935 double Cmax2 = state->end2Curvature_.y();
936 double dthick = (Cmin1 == 0.0) ? 1.0 : 0.0;
937 double br = pbProps_->pb_rmul_ * 1.0 / std::max(Cmax1,Cmax2);
938 if (userArea_)
939#ifdef THREED
940 br = std::sqrt(userArea_ / dPi);
941#else
942 br = userArea_ / 2.0;
943#endif
944 double br2 = br*br;
945 double pbArea = dthick <= 0.0 ? dPi*br2 : 2.0*br*dthick;
946 double bi = dthick <= 0.0 ? 0.25*pbArea*br2 : 2.0*br*br2*dthick/3.0;
947 pbProps_->pbTransStiff_.rx() = pbProps_->pb_kn_*pbArea;
948 pbProps_->pbTransStiff_.ry() = pbProps_->pb_ks_*pbArea;
949#if DIM==3
950 pbProps_->pbAngStiff_.rx() = pbProps_->pb_ks_* 2.0 * bi;
951 pbProps_->pbAngStiff_.ry() = pbProps_->pb_kn_* bi;
952#endif
953 pbProps_->pbAngStiff_.rz() = pbProps_->pb_kn_* bi;
954 }
955 }
956
957 double ContactModelLinearPBond::pbStrainEnergy() const {
958 double ret(0.0);
959 if (pbProps_->pb_kn_)
960 ret = 0.5 * pbProps_->pb_F_.x() * pbProps_->pb_F_.x() / pbProps_->pbTransStiff_.x();
961 if (pbProps_->pb_ks_) {
962 DVect tmp = pbProps_->pb_F_;
963 tmp.rx() = 0.0;
964 double smag2 = tmp.mag2();
965 ret += 0.5 * smag2 / pbProps_->pbTransStiff_.y();
966 }
967
968#ifdef THREED
969 if (pbProps_->pbAngStiff_.x())
970 ret += 0.5 * pbProps_->pb_M_.x() * pbProps_->pb_M_.x() / pbProps_->pbAngStiff_.x();
971#endif
972 if (pbProps_->pbAngStiff_.z()) {
973 DAVect tmp = pbProps_->pb_M_;
974#ifdef THREED
975 tmp.rx() = 0.0;
976 double smag2 = tmp.mag2();
977#else
978 double smag2 = tmp.z() * tmp.z();
979#endif
980 ret += 0.5 * smag2 / pbProps_->pbAngStiff_.z();
981 }
982 return ret;
983 }
984
985 bool ContactModelLinearPBond::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
986 assert(state);
987
988 double overlap = rgap_ - state->gap_;
989 DVect trans = state->relativeTranslationalIncrement_;
990 double correction = 1.0;
991
992 if (state->activated()) {
993 if (cmEvents_[fActivated] >= 0) {
994 auto c = state->getContact();
995 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
996 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
997 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
998 }
999 if (lin_mode_ == 0 && trans.x()) {
1000 correction = -1.0*overlap / trans.x();
1001 if (correction < 0)
1002 correction = 1.0;
1003 }
1004 }
1005
1006#ifdef THREED
1007 DVect norm(trans.x(),0.0,0.0);
1008#else
1009 DVect norm(trans.x(),0.0);
1010#endif
1011 DAVect ang = state->relativeAngularIncrement_;
1012 DVect ss_F_old = lin_F_;
1013
1014 if (lin_mode_==0)
1015 lin_F_.rx() = overlap * kn_;
1016 else
1017 lin_F_.rx() -= correction * norm.x() * kn_;
1018 lin_F_.rx() = std::max(0.0,lin_F_.x());
1019
1020 DVect u_s = trans;
1021 u_s.rx() = 0.0;
1022 DVect sforce = lin_F_ - u_s * ks_ * correction;
1023 sforce.rx() = 0.0;
1024 // Make sure that the shear force opposses the direction of translation - otherwise you can
1025 // have strange behavior
1026 //for (int j=1; j<dim; ++j)
1027 //{
1028 // qDebug()<<sforce.dof(j)<<trans.dof(j);
1029 // if (sign<double>(1,sforce.dof(j)) == sign<double>(1,trans.dof(j)))
1030 // sforce.rdof(j) = 0.0;
1031 //}
1032
1033 // Resolve sliding
1034 if (state->canFail_) {
1035 double crit = lin_F_.x() * fric_;
1036 double sfmag = sforce.mag();
1037 if (sfmag > crit) {
1038 double rat = crit / sfmag;
1039 sforce *= rat;
1040 if (!lin_S_ && cmEvents_[fSlipChange] >= 0) {
1041 auto c = state->getContact();
1042 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1043 fish::Parameter() };
1044 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1045 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1046 }
1047 lin_S_ = true;
1048 } else {
1049 if (lin_S_) {
1050 if (cmEvents_[fSlipChange] >= 0) {
1051 auto c = state->getContact();
1052 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1053 fish::Parameter((qint64)1) };
1054 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1055 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1056 }
1057 lin_S_ = false;
1058 }
1059 }
1060 }
1061
1062 sforce.rx() = lin_F_.x();
1063 lin_F_ = sforce; // total force in linear contact model
1064
1065 // Account for dashpot forces
1066 if (dpProps_) {
1067 dpProps_->dp_F_.fill(0.0);
1068 double vcn(0.0), vcs(0.0);
1069 setDampCoefficients(state->inertialMass_,&vcn,&vcs);
1070 // Need to change behavior based on the dp_mode
1071 if (dpProps_->dp_mode_ == 0) { // Damp all components
1072 dpProps_->dp_F_ = u_s * (-1.0* vcs) / timestep; // shear component
1073 dpProps_->dp_F_ -= norm * vcn / timestep; // normal component
1074 } else if (dpProps_->dp_mode_ == 1) { // limit the tensile
1075 dpProps_->dp_F_ -= norm * vcn / timestep; // normal component
1076 if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
1077 dpProps_->dp_F_.rx() = - lin_F_.rx();
1078 } else if (dpProps_->dp_mode_ == 2) { // limit the shear
1079 if (!lin_S_)
1080 dpProps_->dp_F_ = u_s * (-1.0* vcs) / timestep; // shear component
1081 } else {
1082 if (!lin_S_)
1083 dpProps_->dp_F_ = u_s * (-1.0* vcs) / timestep; // shear component
1084 dpProps_->dp_F_ -= norm * vcn / timestep; // normal component
1085 if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
1086 dpProps_->dp_F_.rx() = - lin_F_.rx();
1087 }
1088 }
1089
1090 // Account for parallel bonds
1091 bool doPb = false;
1092 if (pbProps_ && pbProps_->pb_state_ > 2) {
1093 doPb = true;
1094 // Parallel Bond Logic:
1095 // bond area and inertia
1096 // minimal curvature of end1
1097 double Cmin1 = state->end1Curvature_.x();
1098 double Cmax1 = state->end1Curvature_.y();
1099 double Cmax2 = state->end2Curvature_.y();
1100 double dthick = (Cmin1 == 0.0) ? 1.0 : 0.0;
1101 double br = pbProps_->pb_rmul_ * 1.0 / std::max(Cmax1,Cmax2);
1102 if (userArea_)
1103#ifdef THREED
1104 br = std::sqrt(userArea_ / dPi);
1105#else
1106 br = userArea_ / 2.0;
1107#endif
1108 double br2 = br*br;
1109 double pbArea = dthick <= 0.0 ? dPi*br2 : 2.0*br*dthick;
1110 double bi = dthick <= 0.0 ? 0.25*pbArea*br2 : 2.0*br*br2*dthick/3.0;
1111 pbProps_->pbTransStiff_.rx() = pbProps_->pb_kn_*pbArea;
1112 pbProps_->pbTransStiff_.ry() = pbProps_->pb_ks_*pbArea;
1113
1114 /* elastic force increments */
1115 pbProps_->pb_F_ -= norm *(pbProps_->pb_kn_*pbArea) + u_s * (pbProps_->pb_ks_*pbArea);
1116
1117 /* elastic moment increments */
1118 //DAVect pbMomentInc(0.0);
1119#if DIM==3
1120 pbProps_->pbAngStiff_.rx() = pbProps_->pb_ks_* 2.0 * bi;
1121 pbProps_->pbAngStiff_.ry() = pbProps_->pb_kn_* bi;
1122#endif
1123 pbProps_->pbAngStiff_.rz() = pbProps_->pb_kn_* bi;
1124
1125 DAVect pbMomentInc = ang * pbProps_->pbAngStiff_ *(-1.0);
1126 pbProps_->pb_M_ += pbMomentInc;
1127
1128 /* check bond failure */
1129 if (state->canFail_) {
1130 /* maximum stresses */
1131 double dbend = sqrt(pbProps_->pb_M_.y()*pbProps_->pb_M_.y() + pbProps_->pb_M_.z()*pbProps_->pb_M_.z());
1132 double dtwist = pbProps_->pb_M_.x();
1133 DVect bfs(pbProps_->pb_F_);
1134 bfs.rx() = 0.0;
1135 double dbfs = bfs.mag();
1136 double nsmax = -(pbProps_->pb_F_.x() / pbArea) + pbProps_->pb_mcf_ * dbend * br/bi;
1137 double ssmax = dbfs / pbArea + pbProps_->pb_mcf_ * std::abs(dtwist) * 0.5* br/bi;
1138 double ss ;
1139
1140 if (nsmax >= pbProps_->pb_ten_) {
1141 // Failed in tension
1142 double se = pbStrainEnergy(); // bond strain energy at the onset of failure
1143 pbProps_->pb_state_ = 1;
1144 pbProps_->pb_F_.fill(0.0);
1145 pbProps_->pb_M_.fill(0.0);
1146 if (cmEvents_[fBondBreak] >= 0) {
1147 auto c = state->getContact();
1148 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1149 fish::Parameter((qint64)pbProps_->pb_state_),
1150 fish::Parameter(pbProps_->pb_ten_),
1151 fish::Parameter(se) };
1152 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1153 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1154 }
1155 } else if ((ss = pbShearStrength(pbArea)) <= ssmax) {
1156 // Failed in shear
1157 double se = pbStrainEnergy(); // bond strain energy at the onset of failure
1158 pbProps_->pb_state_ = 2;
1159 pbProps_->pb_F_.fill(0.0);
1160 pbProps_->pb_M_.fill(0.0);
1161 if (cmEvents_[fBondBreak] >= 0) {
1162 auto c = state->getContact();
1163 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1164 fish::Parameter((qint64)pbProps_->pb_state_),
1165 fish::Parameter(ss),
1166 fish::Parameter(se) };
1167 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1168 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1169 }
1170 }
1171 }
1172 }
1173
1174 // Compute energies
1175 if (state->trackEnergy_) {
1176 assert(energies_);
1177 energies_->estrain_ = 0.0;
1178 energies_->epbstrain_ = 0.0;
1179 if (kn_)
1180 energies_->estrain_ = 0.5 * lin_F_.x()* lin_F_.x() /kn_;
1181 if (ks_) {
1182 DVect s = lin_F_;
1183 s.rx() = 0.0;
1184 double smag2 = s.mag2();
1185 energies_->estrain_ += 0.5* smag2 / ks_ ;
1186
1187 if (lin_S_) {
1188 ss_F_old.rx() = 0.0;
1189 DVect avg_F_s = (s + ss_F_old)*0.5;
1190 DVect u_s_el = (s - ss_F_old) / ks_;
1191 energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
1192 }
1193 }
1194 if (dpProps_)
1195 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1196 if (doPb)
1197 energies_->epbstrain_ = pbStrainEnergy();
1198 }
1199
1200 assert(lin_F_ == lin_F_);
1201 return checkActivity(state->gap_);
1202 }
1203
1204 bool ContactModelLinearPBond::thermalCoupling(ContactModelMechanicalState *ms,ContactModelThermalState *ts, IContactThermal *ct,const double &) {
1205 // Account for thermal expansion in linear group in incremental mode
1206 bool ret = false;
1207 if (lin_mode_ > 0 && ts->gapInc_ > 0.0) {
1208 DVect finc(0.0);
1209 finc.rx() = kn_ * ts->gapInc_;
1210 lin_F_ -= finc;
1211 ret = true;
1212 }
1213
1214 if (!pbProps_) return ret;
1215 if (pbProps_->pb_state_ < 3) return ret;
1216 int idx = ct->getModel()->getContactModel()->isProperty("thexp");
1217 if (idx<=0 ) return ret;
1218
1219 double thexp = (ct->getModel()->getContactModel()->getProperty(idx)).toDouble();
1220 double length = ts->length_;
1221 double delTemp =ts->tempInc_;
1222 double delUn = length * thexp * delTemp;
1223 if (delUn == 0.0) return ret;
1224
1225 double dthick = 0.0;
1226 double Cmin1 = ms->end1Curvature_.x();
1227 double Cmax1 = ms->end1Curvature_.y();
1228 double Cmin2 = ms->end2Curvature_.x();
1229 double Cmax2 = ms->end2Curvature_.y();
1230
1231 Cmin2;
1232 if (Cmin1 == 0.0)
1233 dthick = 1.0;
1234
1235 double br = pbProps_->pb_rmul_ * 1.0 / std::max(Cmax1,Cmax2);
1236 if (userArea_)
1237#ifdef THREED
1238 br = std::sqrt(userArea_ / dPi);
1239#else
1240 br = userArea_ / 2.0;
1241#endif
1242 double br2 = br*br;
1243 double pbArea = dthick <= 0.0 ? dPi*br2 : 2.0*br*dthick;
1244 //
1245 DVect finc(0.0);
1246 finc.rx() = pbProps_->pb_kn_ * pbArea * delUn;
1247 pbProps_->pb_F_ += finc;
1248
1249 ret = true;
1250 return ret;
1251 }
1252
1253 void ContactModelLinearPBond::setForce(const DVect &v,IContact *c) {
1254 lin_F(v);
1255 if (v.x() > 0)
1256 rgap_ = c->getGap() + v.x() / kn_;
1257 }
1258
1259 void ContactModelLinearPBond::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
1260 // Only do something if the contact model is of the same type
1261 if (old->getContactModel()->getName().compare("linearpbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
1262 ContactModelLinearPBond *oldCm = (ContactModelLinearPBond *)old;
1263#ifdef THREED
1264 // Need to rotate just the shear component from oldSystem to newSystem
1265
1266 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
1267 DVect axis = oldSystem.e1() & newSystem.e1();
1268 double c, ang, s;
1269 DVect re2;
1270 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1271 axis = axis.unit();
1272 c = oldSystem.e1()|newSystem.e1();
1273 if (c > 0)
1274 c = std::min(c,1.0);
1275 else
1276 c = std::max(c,-1.0);
1277 ang = acos(c);
1278 s = sin(ang);
1279 double t = 1. - c;
1280 DMatrix<3,3> rm;
1281 rm.get(0,0) = t*axis.x()*axis.x() + c;
1282 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
1283 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
1284 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
1285 rm.get(1,1) = t*axis.y()*axis.y() + c;
1286 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
1287 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
1288 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
1289 rm.get(2,2) = t*axis.z()*axis.z() + c;
1290 re2 = rm*oldSystem.e2();
1291 }
1292 else
1293 re2 = oldSystem.e2();
1294 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
1295 axis = re2 & newSystem.e2();
1296 DVect2 tpf;
1297 DMatrix<2,2> m;
1298 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1299 axis = axis.unit();
1300 c = re2|newSystem.e2();
1301 if (c > 0)
1302 c = std::min(c,1.0);
1303 else
1304 c = std::max(c,-1.0);
1305 ang = acos(c);
1306 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
1307 ang *= -1;
1308 s = sin(ang);
1309 m.get(0,0) = c;
1310 m.get(1,0) = s;
1311 m.get(0,1) = -m.get(1,0);
1312 m.get(1,1) = m.get(0,0);
1313 tpf = m*DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
1314 } else {
1315 m.get(0,0) = 1.;
1316 m.get(0,1) = 0.;
1317 m.get(1,0) = 0.;
1318 m.get(1,1) = 1.;
1319 tpf = DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
1320 }
1321 DVect pforce = DVect(0,tpf.x(),tpf.y());
1322#else
1323 oldSystem;
1324 newSystem;
1325 DVect pforce = DVect(0,oldCm->lin_F_.y());
1326#endif
1327 for (int i=1; i<dim; ++i)
1328 lin_F_.rdof(i) += pforce.dof(i);
1329 if (lin_mode_ && oldCm->lin_mode_)
1330 lin_F_.rx() = oldCm->lin_F_.x();
1331 oldCm->lin_F_ = DVect(0.0);
1332 if (dpProps_ && oldCm->dpProps_) {
1333#ifdef THREED
1334 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
1335 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
1336#else
1337 pforce = oldCm->dpProps_->dp_F_;
1338#endif
1339 dpProps_->dp_F_ += pforce;
1340 oldCm->dpProps_->dp_F_ = DVect(0.0);
1341 }
1342 if(oldCm->getEnergyActivated()) {
1343 activateEnergy();
1344 energies_->estrain_ = oldCm->energies_->estrain_;
1345 energies_->eslip_ = oldCm->energies_->eslip_;
1346 energies_->edashpot_ = oldCm->energies_->edashpot_;
1347 energies_->epbstrain_ = oldCm->energies_->epbstrain_;
1348 oldCm->energies_->estrain_ = 0.0;
1349 oldCm->energies_->edashpot_ = 0.0;
1350 oldCm->energies_->eslip_ = 0.0;
1351 oldCm->energies_->epbstrain_ = 0.0;
1352 }
1353 rgap_ = oldCm->rgap_;
1354 }
1355 assert(lin_F_ == lin_F_);
1356 }
1357
1358 void ContactModelLinearPBond::setNonForcePropsFrom(IContactModel *old) {
1359 // Only do something if the contact model is of the same type
1360 if (old->getName().compare("linearpbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
1361 ContactModelLinearPBond *oldCm = (ContactModelLinearPBond *)old;
1362 kn_ = oldCm->kn_;
1363 ks_ = oldCm->ks_;
1364 fric_ = oldCm->fric_;
1365 lin_mode_ = oldCm->lin_mode_;
1366 rgap_ = oldCm->rgap_;
1367 userArea_ = oldCm->userArea_;
1368
1369 if (oldCm->dpProps_) {
1370 if (!dpProps_)
1371 dpProps_ = NEWC(dpProps());
1372 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1373 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1374 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1375 }
1376
1377 if (oldCm->pbProps_) {
1378 if (!pbProps_)
1379 pbProps_ = NEWC(pbProps());
1380 pbProps_->pb_rmul_ = oldCm->pbProps_->pb_rmul_;
1381 pbProps_->pb_kn_ = oldCm->pbProps_->pb_kn_;
1382 pbProps_->pb_ks_ = oldCm->pbProps_->pb_ks_;
1383 pbProps_->pb_mcf_ = oldCm->pbProps_->pb_mcf_;
1384 pbProps_->pb_fa_ = oldCm->pbProps_->pb_fa_;
1385 pbProps_->pb_state_ = oldCm->pbProps_->pb_state_;
1386 pbProps_->pb_coh_ = oldCm->pbProps_->pb_coh_;
1387 pbProps_->pb_ten_ = oldCm->pbProps_->pb_ten_;
1388 pbProps_->pbTransStiff_ = oldCm->pbProps_->pbTransStiff_;
1389 pbProps_->pbAngStiff_ = oldCm->pbProps_->pbAngStiff_;
1390 }
1391 }
1392 }
1393
1394 DVect ContactModelLinearPBond::getForce() const {
1395 DVect ret(lin_F_);
1396 if (dpProps_)
1397 ret += dpProps_->dp_F_;
1398 if (pbProps_)
1399 ret += pbProps_->pb_F_;
1400 return ret;
1401 }
1402
1403 DAVect ContactModelLinearPBond::getMomentOn1(const IContactMechanical *c) const {
1404 DAVect ret(0.0);
1405 if (pbProps_)
1406 ret = pbProps_->pb_M_;
1407 if (c) {
1408 DVect force = getForce();
1409 c->updateResultingTorqueOn1Local(force,&ret);
1410 }
1411 return ret;
1412 }
1413
1414 DAVect ContactModelLinearPBond::getMomentOn2(const IContactMechanical *c) const {
1415 DAVect ret(0.0);
1416 if (pbProps_)
1417 ret = pbProps_->pb_M_;
1418 if (c) {
1419 DVect force = getForce();
1420 c->updateResultingTorqueOn2Local(force,&ret);
1421 }
1422 return ret;
1423 }
1424
1425 DAVect ContactModelLinearPBond::getModelMomentOn1() const {
1426 DAVect ret(0.0);
1427 if (pbProps_)
1428 ret = pbProps_->pb_M_;
1429 return ret;
1430 }
1431
1432 DAVect ContactModelLinearPBond::getModelMomentOn2() const {
1433 DAVect ret(0.0);
1434 if (pbProps_)
1435 ret = pbProps_->pb_M_;
1436 return ret;
1437 }
1438
1439 void ContactModelLinearPBond::objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *ret) const {
1440 ret->clear();
1441 ret->push_back({"kn",ScalarInfo});
1442 ret->push_back({"ks",ScalarInfo});
1443 ret->push_back({"fric",ScalarInfo});
1444 ret->push_back({"lin_force",VectorInfo});
1445 ret->push_back({"lin_slip",ScalarInfo});
1446 ret->push_back({"lin_mode",ScalarInfo});
1447 ret->push_back({"rgap",ScalarInfo});
1448 ret->push_back({"emod",ScalarInfo});
1449 ret->push_back({"kratio",ScalarInfo});
1450 ret->push_back({"dp_nratio",ScalarInfo});
1451 ret->push_back({"dp_sratio",ScalarInfo});
1452 ret->push_back({"dp_mode",ScalarInfo});
1453 ret->push_back({"dp_force",VectorInfo});
1454 ret->push_back({"pb_state",ScalarInfo});
1455 ret->push_back({"pb_rmul",ScalarInfo});
1456 ret->push_back({"pb_kn",ScalarInfo});
1457 ret->push_back({"pb_ks",ScalarInfo});
1458 ret->push_back({"pb_mcf",ScalarInfo});
1459 ret->push_back({"pb_ten",ScalarInfo});
1460 ret->push_back({"pb_shear",ScalarInfo});
1461 ret->push_back({"pb_coh",ScalarInfo});
1462 ret->push_back({"pb_fa",ScalarInfo});
1463 ret->push_back({"pb_sigma",ScalarInfo});
1464 ret->push_back({"pb_tau",ScalarInfo});
1465 ret->push_back({"pb_force",VectorInfo});
1466 ret->push_back({"pb_moment",VectorInfo});
1467 ret->push_back({"pb_radius",ScalarInfo});
1468 ret->push_back({"pb_emod",ScalarInfo});
1469 ret->push_back({"pb_kratio",ScalarInfo});
1470 ret->push_back({"user_area",ScalarInfo});
1471 }
1472
1473 void ContactModelLinearPBond::objectPropValues(std::vector<double> *ret,const IContact *con) const {
1474 ret->clear();
1475 ret->push_back(kn());
1476 ret->push_back(ks());
1477 ret->push_back(fric());
1478 ret->push_back(lin_F().mag());
1479 ret->push_back(lin_S());
1480 ret->push_back(lin_mode());
1481 ret->push_back(rgap());
1482 double d;
1483 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
1484 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
1485 double rsum(0.0);
1486 if (c->getEnd1Curvature().y())
1487 rsum += 1.0/c->getEnd1Curvature().y();
1488 if (c->getEnd2Curvature().y())
1489 rsum += 1.0/c->getEnd2Curvature().y();
1490 if (userArea_) {
1491#ifdef THREED
1492 rsq = std::sqrt(userArea_ / dPi);
1493#else
1494 rsq = userArea_ / 2.0;
1495#endif
1496 rsum = rsq + rsq;
1497 rsq = 1. / rsq;
1498 }
1499#ifdef TWOD
1500 d = (kn_ * rsum * rsq / 2.0);
1501#else
1502 d = (kn_ * rsum * rsq * rsq) / dPi;
1503#endif
1504 ret->push_back(d);
1505 ret->push_back((ks_ == 0.0) ? 0.0 : (kn_/ks_));
1506 ret->push_back(dp_nratio());
1507 ret->push_back(dp_sratio());
1508 ret->push_back(dp_mode());
1509 ret->push_back(dp_F().mag());
1510 ret->push_back(pbProps_ ? pbProps_->pb_state_ : 0);
1511 ret->push_back(pbProps_ ? pbProps_->pb_rmul_ : 1);
1512 ret->push_back(pbProps_ ? pbProps_->pb_kn_ : 0);
1513 ret->push_back(pbProps_ ? pbProps_->pb_ks_ : 0);
1514 ret->push_back(pbProps_ ? pbProps_->pb_mcf_ : 1);
1515 ret->push_back(pbProps_ ? pbProps_->pb_ten_ : 0);
1516 if (!pbProps_) d = 0; else d = pbShearStrength(pbData(c).x());
1517 ret->push_back(d);
1518 ret->push_back(pbProps_ ? pbProps_->pb_coh_ : 0);
1519 ret->push_back(pbProps_ ? pbProps_->pb_fa_ : 0);
1520 if (!pbProps_ || pbProps_->pb_state_ < 3) d =0.0; else pbSMax(c).x();
1521 ret->push_back(d);
1522 if (!pbProps_ || pbProps_->pb_state_ < 3) d =0.0; else pbSMax(c).y();
1523 ret->push_back(d);
1524 ret->push_back(pbProps_ ? (pbProps_->pb_F_).mag() : 0);
1525 ret->push_back(pbProps_ ? (pbProps_->pb_M_).mag() : 0);
1526 d = 0;
1527 if (pbProps_) {
1528 double Cmax1 = c->getEnd1Curvature().y();
1529 double Cmax2 = c->getEnd2Curvature().y();
1530 double br = pbProps_->pb_rmul_ * 1.0 / std::max(Cmax1,Cmax2);
1531 if (userArea_)
1532#ifdef THREED
1533 d = std::sqrt(userArea_ / dPi);
1534#else
1535 d = userArea_ / 2.0;
1536#endif
1537 }
1538 ret->push_back(d);
1539 d = 0;
1540 if (c->getEnd1Curvature().y())
1541 rsum += 1.0/c->getEnd1Curvature().y();
1542 if (c->getEnd2Curvature().y())
1543 rsum += 1.0/c->getEnd2Curvature().y();
1544 if (userArea_) {
1545#ifdef THREED
1546 double rad = std::sqrt(userArea_ / dPi);
1547#else
1548 double rad = userArea_ / 2.0;
1549#endif
1550 rsum = 2 * rad;
1551 }
1552 d = pbProps_ ? pbProps_->pb_kn_ * rsum : 0;
1553 ret->push_back(d);
1554 d = 0;
1555 if (pbProps_) d = (pbProps_->pb_ks_ == 0.0) ? 0.0 : (pbProps_->pb_kn_/pbProps_->pb_ks_);
1556 ret->push_back(d);
1557 ret->push_back(userArea_);
1558 }
1559
1560 DVect3 ContactModelLinearPBond::pbData(const IContactMechanical *c) const {
1561 double Cmax1 = c->getEnd1Curvature().y();
1562 double Cmax2 = c->getEnd2Curvature().y();
1563 double br = pbProps_->pb_rmul_ * 1.0 / std::max(Cmax1,Cmax2);
1564 if (userArea_)
1565#ifdef THREED
1566 br = std::sqrt(userArea_ / dPi);
1567#else
1568 br = userArea_ / 2.0;
1569#endif
1570 double br2 = br*br;
1571#ifdef TWOD
1572 double pbArea = 2.0*br;
1573 double bi = 2.0*br*br2/3.0;
1574#else
1575 double pbArea = dPi*br2;
1576 double bi = 0.25*pbArea*br2;
1577#endif
1578 return DVect3(pbArea,bi,br);
1579 }
1580
1581 DVect2 ContactModelLinearPBond::pbSMax(const IContactMechanical *c) const {
1582 DVect3 data = pbData(c);
1583 double pbArea = data.x();
1584 double bi = data.y();
1585 double br = data.z();
1586 /* maximum stresses */
1587 double dbend = sqrt(pbProps_->pb_M_.y()*pbProps_->pb_M_.y() + pbProps_->pb_M_.z()*pbProps_->pb_M_.z());
1588 double dtwist = pbProps_->pb_M_.x();
1589 DVect bfs(pbProps_->pb_F_);
1590 bfs.rx() = 0.0;
1591 double dbfs = bfs.mag();
1592 double nsmax = -(pbProps_->pb_F_.x() / pbArea) + pbProps_->pb_mcf_ * dbend * br/bi;
1593 double ssmax = dbfs / pbArea + pbProps_->pb_mcf_ * std::abs(dtwist) * 0.5* br/bi;
1594 return DVect2(nsmax,ssmax);
1595 }
1596
1597 double ContactModelLinearPBond::pbShearStrength(const double &pbArea) const {
1598 if (!pbProps_) return 0.0;
1599 double sig = -1.0*pbProps_->pb_F_.x() / pbArea;
1600 double nstr = pbProps_->pb_state_ > 2 ? pbProps_->pb_ten_ : 0.0;
1601 return sig <= nstr ? pbProps_->pb_coh_ - std::tan(dDegrad*pbProps_->pb_fa_)*sig
1602 : pbProps_->pb_coh_ - std::tan(dDegrad*pbProps_->pb_fa_)*nstr;
1603 }
1604
1605 void ContactModelLinearPBond::setDampCoefficients(const double &mass,double *vcn,double *vcs) {
1606 *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(kn_));
1607 *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(ks_));
1608 }
1609
1610} // namespace itascaxd
1611// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Dec 14, 2024 |
