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