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