Flat-Joint Model Implementation
See this page for the documentation of this contact model.
contactmodelflatjoint.h
1#pragma once
2// contactmodelflatjoint.h
3
4#include "contactmodel/src/contactmodelmechanical.h"
5
6#ifdef FLATJOINT_LIB
7# define FLATJOINT_EXPORT EXPORT_TAG
8#elif defined(NO_MODEL_IMPORT)
9# define FLATJOINT_EXPORT
10#else
11# define FLATJOINT_EXPORT IMPORT_TAG
12#endif
13
14namespace cmodelsxd {
15 using namespace itasca;
16
17 class ContactModelFlatJoint : public ContactModelMechanical {
18 public:
19 enum PropertyKeys {
20 kwFjNr=1
21 , kwFjElem
22 , kwFjKn
23 , kwFjKs
24 , kwFjFric
25 , kwFjEmod
26 , kwFjKRatio
27 , kwFjRmul
28 , kwFjRadius
29 , kwFjGap0
30 , kwFjTen
31 , kwFjCoh
32 , kwFjFa
33 , kwFjF
34 , kwFjM
35 , kwFjState
36 , kwFjSlip
37 , kwFjMType
38 , kwFjA
39 , kwFjEgap
40 , kwFjGap
41 , kwFjNstr
42 , kwFjSstr
43 , kwFjSs
44#ifdef THREED
45 , kwFjNa
46#endif
47 , kwFjRelBr
48 , kwFjCen
49 , kwFjTrack
50 , kwUserArea
51 , kwFjCohRes
52 , kwFjResMode
53 };
54
55 FLATJOINT_EXPORT ContactModelFlatJoint();
56 FLATJOINT_EXPORT virtual ~ContactModelFlatJoint();
57 void copy(const ContactModel *c) override;
58 void archive(ArchiveStream &) override;
59 QString getName() const override { return "flatjoint"; }
60 void setIndex(int i) override { index_=i;}
61 int getIndex() const override {return index_;}
62 QString getProperties() const override { return "fj_nr"
63 ",fj_elem"
64 ",fj_kn"
65 ",fj_ks"
66 ",fj_fric"
67 ",fj_emod"
68 ",fj_kratio"
69 ",fj_rmul"
70 ",fj_radius"
71 ",fj_gap0"
72 ",fj_ten"
73 ",fj_coh"
74 ",fj_fa"
75 ",fj_force"
76 ",fj_moment"
77 ",fj_state"
78 ",fj_slip"
79 ",fj_mtype"
80 ",fj_area"
81 ",fj_egap"
82 ",fj_gap"
83 ",fj_sigma"
84 ",fj_tau"
85 ",fj_shear"
86#ifdef THREED
87 ",fj_nal"
88#endif
89 ",fj_relbr"
90 ",fj_cen"
91 ",fj_track"
92 ",user_area"
93 ",fj_cohres"
94 ",fj_resmode"
95 ;}
96
97 enum EnergyKeys { kwEStrain=1,kwESlip};
98 QString getEnergies() const { return "energy-strain,energy-slip";}
99 double getEnergy(uint32 i) const override; // Base 1
100 bool getEnergyAccumulate(uint32 i) const override; // Base 1
101 void setEnergy(uint32 i,const double &d) override; // Base 1
102 void activateEnergy() override { if (energies_) return; energies_ = NEWC(Energies());}
103 bool getEnergyActivated() const override {return (energies_ !=0);}
104
105 enum FishCallEvents {fActivated=0,fBondBreak,fBroken,fSlipChange};
106 QString getFishCallEvents() const override { return "contact_activated,bond_break,broken,all_slip_change"; }
107 QVariant getProperty(uint32 i,const IContact *) const override;
108 bool getPropertyGlobal(uint32 i) const override;
109 bool setProperty(uint32 i,const QVariant &v,IContact *) override;
110 bool getPropertyReadOnly(uint32 i) const override;
111
112 bool supportsInheritance(uint32 ) const override { return false; }
113
114 enum MethodKeys { kwBond=1, kwUnbond, KwDeformability, KwUpdateGeom, kwArea, kwInitialize};
115
116 QString getMethods() const override { return "bond"
117 ",unbond"
118 ",deformability"
119 ",update_geometry"
120 ",area"
121 ",initialize"
122 ;}
123
124 QString getMethodArguments(uint32 i) const override;
125
126 bool setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override; // Base 1 - returns true if timestep contributions need to be updated
127
128 uint32 getMinorVersion() const override;
129
130 bool validate(ContactModelMechanicalState *state,const double ×tep) override;
131 bool endPropertyUpdated(const QString &,const IContactMechanical *) override { return false; }
132 bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) override;
133 bool thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
134 bool coupling(ContactModelMechanicalState*, ContactModelThermalState*,ContactModelFluidState*, const double&) override;
135 DVect2 getEffectiveTranslationalStiffness() const override {
136 return effectiveTranslationalStiffness();
137 }
138 DAVect getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness(); }
139
140 ContactModelFlatJoint *clone() const override { return NEWC(ContactModelFlatJoint()); }
141 double getActivityDistance() const override {return 0.0;}
142 bool isOKToDelete() const override { return !isBonded(); }
143 void resetForcesAndMoments() override { fj_f(DVect(0.0)); fj_m(DAVect(0.0)); for (int i=0; i<f_.size(); ++i) f_[i] = DVect(0.0); }
144 void setForce(const DVect &v,IContact *) override;
145 void setArea(const double &d) override { userArea_ = d; }
146 double getArea() const override { return userArea_; }
147 bool checkActivity(const double &inGap) override;
148
149 /// Return the total force that the contact model holds.
150 DVect getForce(const IContactMechanical*) const override;
151 /// Return the total moment on 1 that the contact model holds
152 DAVect getMomentOn1(const IContactMechanical*) const override;
153 /// Return the total moment on 1 that the contact model holds
154 DAVect getMomentOn2(const IContactMechanical*) const override;
155 //virtual bool isSliding() const { return fj_s_; }
156 bool isBonded() const override { FOR(it,bmode_) if ((*it) == 3) return true; return false; }
157 void unbond() override { FOR(it,bmode_) *it = 0; }
158
159 // For contact specific plotting
160 void getSphereList(const IContact* con, std::vector<DVect>* pos, std::vector<double>* rad, std::vector<double>* val) override;
161#ifdef THREED
162 void getDiskList(const IContact* con, std::vector<DVect>* pos, std::vector<DVect>* normal, std::vector<double>* radius, std::vector<double>* val) override;
163#endif
164 void getCylinderList(const IContact* con, std::vector<DVect>* bot, std::vector<DVect>* top, std::vector<double>* radlow, std::vector<double>* radhi, std::vector<double>* val) override;
165
166 int fj_nr() const {return fj_nr_;}
167 void fj_nr(int d) { fj_nr_= d;}
168#ifdef THREED
169 int fj_n() const { return fj_na_ * fj_nr_; }
170 int fj_na() const {return fj_na_;}
171 void fj_na(int d) { fj_na_= d;}
172#else
173 int fj_n() const { return fj_nr_; }
174#endif
175 int fj_elem() const {return fj_elem_;}
176 void fj_elem(int d) { fj_elem_= d;}
177 const double & fj_kn() const {return fj_kn_;}
178 void fj_kn(const double &d) { fj_kn_ = d;}
179 const double & fj_ks() const {return fj_ks_;}
180 void fj_ks(const double &d) { fj_ks_ = d;}
181 const double & fj_fric() const {return fj_fric_;}
182 void fj_fric(const double &d) { fj_fric_ = d;}
183 const double & fj_rmul() const {return fj_rmul_;}
184 void fj_rmul(const double &d) { fj_rmul_ = d;}
185 const double & fj_gap0() const {return fj_gap0_;}
186 void fj_gap0(const double &d) { fj_gap0_ = d;}
187 const double & fj_ten() const {return fj_ten_;}
188 void fj_ten(const double &d) { fj_ten_ = d;}
189 const double & fj_coh() const {return fj_coh_;}
190 void fj_coh(const double &d) { fj_coh_ = d;}
191 const double & fj_cohres() const {return fj_cohres_;}
192 void fj_cohres(const double &d) { fj_cohres_ = d;}
193 const double & fj_fa() const {return fj_fa_;}
194 void fj_fa(const double &d) { fj_fa_ = d;}
195 const DVect & fj_f() const {return fj_f_;}
196 void fj_f(const DVect &f) { fj_f_=f;}
197 const DAVect & fj_m() const {return fj_m_;}
198 void fj_m(const DAVect &f) { fj_m_=f;}
199 const DAVect & fj_m_set() const {return fj_m_set_;}
200 void fj_m_set(const DAVect &f) { fj_m_set_=f;}
201 const double & rmin() const {return rmin_;}
202 void rmin(const double &d) { rmin_ = d;}
203 const double & rbar() const {return rbar_;}
204 void rbar(const double &d) { rbar_ = d;}
205 const int & fj_resmode() const {return fj_resmode_;}
206 void fj_resmode(const int &i) { fj_resmode_ = i;}
207 const double & atot() const {return atot_;}
208 void atot(const double &d) { atot_ = d;}
209 bool propsFixed() const {return propsFixed_; }
210 void propsFixed(bool d) { propsFixed_ = d;}
211 int mType() const {return mType_; }
212 void mType(int d) { mType_ = d;}
213 const DVect & gap() const {return gap_; }
214 void gap(const DVect &d) { gap_ = d;}
215 const double & theta() const {return theta_; }
216 void theta(const double & d) { theta_ = d;}
217#ifdef THREED
218 const double & thetaM() const {return thetaM_; }
219 void thetaM(const double & d) { thetaM_ = d;}
220#else
221 double thetaM() const {return 0.0;}
222#endif
223
224 bool hasEnergies() const {return energies_ ? true:false;}
225 double estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
226 void estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
227 double eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
228 void eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
229
230 uint32 inheritanceField() const {return inheritanceField_;}
231 void inheritanceField(uint32 i) {inheritanceField_ = i;}
232
233 const DVect2 & effectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
234 void effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
235 const DAVect & effectiveRotationalStiffness() const {return effectiveRotationalStiffness_;}
236 void effectiveRotationalStiffness(const DAVect &v ) {effectiveRotationalStiffness_=v;}
237
238 private:
239 static int index_;
240
241 struct Energies {
242 Energies() : estrain_(0.0), eslip_(0.0) {}
243 double estrain_; // elastic energy stored in contact
244 double eslip_; // work dissipated by friction
245 };
246
247 void updateEffectiveStiffness(ContactModelMechanicalState *state);
248
249 // inheritance fields
250 uint32 inheritanceField_;
251
252 int fj_nr_; // radial number of elements >= 1 (total in 2D)
253#ifdef THREED
254 int fj_na_; // circumferential number of elements >= 3
255#endif
256 int fj_elem_; // Element to be queried
257 double fj_kn_; // normal stiffness
258 double fj_ks_; // shear stiffness
259 double fj_fric_; // Coulomb friction coefficient
260 double fj_rmul_; // Radius multiplier
261 double fj_gap0_; // Initial gap
262 double fj_ten_; // Tensile strength
263 double fj_coh_; // Cohesive strength
264 double fj_cohres_; // Residual cohesive strength
265 double fj_fa_; // Friction angle
266 DVect fj_f_; // Force carried in the model
267 DAVect fj_m_; // Moment carried in the model
268 DAVect fj_m_set_; // When initializing forces then need an extra moment term
269 // Area related quantities
270 double rmin_; // min(Ra,Rb) where Ra & Rb are particle radii
271 double rbar_; // flat-joint radius [m]
272 double atot_; // flat-joint area [m^2]
273 std::vector<double> a_; // cross-sectional area of elem[fj_elem-1] [m^2]
274#ifdef THREED
275 std::vector<DVect2> rBarl_; // centroid relative position of elem[fj_elem-1] [m] (3D)
276#else
277 std::vector<double> rBarl_; // centroid relative position of elem[fj_elem-1] [m] (2D)
278#endif
279 int fj_resmode_; // Residual mode
280 void setAreaQuantities(); // Set Rbar, Atot and A[]
281 DVect getRelElemPos(const IContact*,int ) const; // Return the relative location of element i
282 void setRelElemPos(const IContact*,int ,const DVect &); // Set the relative location of element i
283
284 bool propsFixed_; // {Rmul, N, G, bstate, mType} fixed, cannot reset
285 int mType_; // initial microstructural type
286 int getmType() const; // {1,2,3,4}={bonded, gapped, slit, other}
287
288 std::vector<int> bmode_; // bond mode - 0 unbonded, 1 failed in tension, 2 failed in shear, 3 bonded
289 std::vector<bool> smode_; // slip mode
290 bool Bonded(int e) const { return bmode_[e-1] == 3 ? true : false; }
291
292 // Set bstate and bmode (can only bond if fj_gap0_==0.0)
293 void bondElem(int iSeg,bool bBond);
294 // Set bstate & bmode
295 void breakBond(int iSeg,int fmode,ContactModelMechanicalState *state);
296 void slipChange(int iSeg,bool smode,ContactModelMechanicalState *state);
297
298 // For use in 2D only!
299 double tauC(const double &dSig,bool bBonded) const; // shear strength (positive) [N/m^2]
300
301 // INTERFACE RESPONSE QUANTITIES:
302 DVect gap_; // total relative displacement [m]
303 double theta_; // total relative rotation [rad]
304#ifdef THREED
305 double thetaM_; // total relative rotation [rad]
306 double thbMag() const { return sqrt(theta_*theta_ + thetaM_*thetaM_); }
307 // unit-vector xi of middle surface system xi-eta
308 // (If both thb_l and thb_m are zero, then xi is undefined
309 // and returns zero for both components.)
310 double xi(int comp /* component (l,m) = (1,2) */) const;
311#endif
312 std::vector<double> egap_; // gap at centroid of elem[fj_elem-1] [N]
313 std::vector<DVect> f_; // force on elem[fj_elem-1] [N]
314
315 void initVectors(); // Resize and zero all vector types based on current value of N
316#ifdef TWOD
317 double gap(const double &x) const; // Gap (g>0 is open) along the interface, x in [0, 2*Rbar]
318#else
319 double gap(const double &rl,const double &rm) const; // Gap (g>0 is open) gap at relative position (l,m) [m]
320 double sigBar( int e /* element, e = 1,2,...,Nel */ ) const; // normal stress at centroid of elem[eN-1] [N/m^2]
321 double tauBar( int e /* element, e = 1,2,...,Nel */ ) const; // shear stress at centroid of elem[eN-1] [N/m^2]
322#endif
323 double computeStrainEnergy(int e /* element, e = 1,2,...,Nel */) const; // strain energy in elem[eN-1]
324 // For use in 2D only! Segment normal stress
325 double computeSig(const double &g0, /* gap at left end */
326 const double &g1, /* gap at right end */
327 const double &rbar, /* length is 2*rbar */
328 const double &dA, /* area */
329 bool bBonded /* bond state */
330 ) const;
331 // For use in 2D only! Segment moment
332 double computeM(const double &g0, /* gap at left end */
333 const double &g1, /* gap at right end */
334 const double &rbar, /* length is 2*rbar */
335 bool bBonded /* bond state */
336 ) const;
337 // For use in 2D only! getCase used by ComputeSig and ComputeM
338 int getCase(const double &g0, /* gap at left end */
339 const double &g1 /* gap at right end */
340 ) const;
341 // Segment elastic shear-displacement increment, which is portion of
342 // increment that occurs while gap is negative.
343 double delUse(const double &gapStart, /* gap at start of FDlaw */
344 const double &gapEnd, /* gap at end of FDlaw */
345 bool bBonded, /* bond state */
346 const double &delUs /* shear displ. increment */
347 ) const;
348 double userArea_; // Area as specified by the user
349 Energies * energies_; // energies
350
351 DVect2 effectiveTranslationalStiffness_;
352 DAVect effectiveRotationalStiffness_;
353
354 struct orientProps {
355 orientProps() : origNormal_(DVect(0.0)) {}
356 Quat orient1_;
357 Quat orient2_;
358 DVect origNormal_;
359 };
360
361 orientProps *orientProps_;
362 };
363} // namespace itascaxd
364
365
366// EoF
contactmodelflatjoint.cpp
1// contactmodelflatjoint.cpp
2#include "contactmodelflatjoint.h"
3
4#include "../version.txt"
5#include "fish/src/parameter.h"
6#include "utility/src/tptr.h"
7#include "shared/src/mathutil.h"
8#include "base/src/basetoqt.h"
9#include "contactmodel/src/contactmodelthermal.h"
10#include "contactmodel/src/contactmodelfluid.h"
11
12#include "kernel/interface/iprogram.h"
13#include "module/interface/icontact.h"
14#include "module/interface/icontactmechanical.h"
15#include "module/interface/icontactthermal.h"
16#include "module/interface/icontactfluid.h"
17#include "module/interface/ifishcalllist.h"
18#include "module/interface/ipiece.h"
19#include "module/interface/ipiecemechanical.h"
20
21#ifdef FLATJOINT_LIB
22#ifdef _WIN32
23 int __stdcall DllMain(void *,unsigned, void *)
24 {
25 return 1;
26 }
27#endif
28
29 extern "C" EXPORT_TAG const char *getName()
30 {
31#if DIM==3
32 return "contactmodelmechanical3dflatjoint";
33#else
34 return "contactmodelmechanical2dflatjoint";
35#endif
36 }
37
38 extern "C" EXPORT_TAG unsigned getMajorVersion()
39 {
40 return MAJOR_VERSION;
41 }
42
43 extern "C" EXPORT_TAG unsigned getMinorVersion()
44 {
45 return MINOR_VERSION;
46 }
47
48 extern "C" EXPORT_TAG void *createInstance()
49 {
50 cmodelsxd::ContactModelFlatJoint *m = NEWC(cmodelsxd::ContactModelFlatJoint());
51 return (void *)m;
52 }
53#endif // FLATJOINT_LIB
54
55namespace cmodelsxd {
56 static const uint32 fjKnMask = 0x00002; // Base 1!
57 static const uint32 fjKsMask = 0x00004;
58 static const uint32 fjFricMask = 0x00008;
59
60 using namespace itasca;
61
62 int ContactModelFlatJoint::index_ = -1;
63 uint32 ContactModelFlatJoint::getMinorVersion() const { return MINOR_VERSION;}
64
65 ContactModelFlatJoint::ContactModelFlatJoint() : inheritanceField_(fjKnMask|fjKsMask|fjFricMask)
66 , fj_nr_(2)
67#ifdef THREED
68 , fj_na_(4)
69#endif
70 , fj_elem_(1)
71 , fj_kn_(0.0)
72 , fj_ks_(0.0)
73 , fj_fric_(0.0)
74 , fj_rmul_(1.0)
75 , fj_gap0_(0.0)
76 , fj_ten_(0.0)
77 , fj_coh_(0.0)
78 , fj_cohres_(0.0)
79 , fj_fa_(0.0)
80 , fj_f_(0.0)
81 , fj_m_(0.0)
82 , fj_m_set_(0.0)
83 , rmin_(1.0)
84 , rbar_(0.0)
85 , atot_(0.0)
86 , a_(2)
87 , rBarl_(2)
88 , fj_resmode_(0)
89 , propsFixed_(false)
90 , mType_(3)
91 , bmode_(2)
92 , smode_(2)
93 , gap_(0.0)
94 , theta_(0.0)
95#ifdef THREED
96 , thetaM_(0.0)
97#endif
98 , egap_(2)
99 , f_(2)
100 , userArea_(0)
101 , energies_(0)
102 , effectiveTranslationalStiffness_(DVect2(0.0))
103 , effectiveRotationalStiffness_(DAVect(0.0))
104 , orientProps_(0)
105 {
106 initVectors();
107 setAreaQuantities();
108 //setFromParent(ContactModelMechanicalList::instance()->find(getName()));
109 }
110
111 ContactModelFlatJoint::~ContactModelFlatJoint() {
112 if (orientProps_)
113 delete orientProps_;
114 if (energies_)
115 delete energies_;
116 }
117
118 void ContactModelFlatJoint::archive(ArchiveStream &stream) {
119 stream & fj_nr_;
120#ifdef THREED
121 stream & fj_na_;
122#endif
123 stream & fj_elem_;
124 stream & fj_kn_;
125 stream & fj_ks_;
126 stream & fj_fric_;
127 stream & fj_rmul_;
128 stream & fj_gap0_;
129 stream & fj_ten_;
130 stream & fj_coh_;
131 stream & fj_fa_;
132 stream & fj_f_;
133 stream & fj_m_;
134 stream & rmin_;
135 stream & rbar_;
136 stream & atot_;
137 stream & a_;
138 stream & rBarl_;
139 stream & propsFixed_;
140 stream & mType_;
141 stream & bmode_;
142 stream & smode_;
143 stream & gap_;
144 stream & theta_;
145#ifdef THREED
146 stream & thetaM_;
147#endif
148 stream & egap_;
149 stream & f_;
150
151 if (stream.getArchiveState()==ArchiveStream::Save) {
152 bool b = false;
153 if (orientProps_) {
154 b = true;
155 stream & b;
156 stream & orientProps_->orient1_;
157 stream & orientProps_->orient2_;
158 stream & orientProps_->origNormal_;
159 } else
160 stream & b;
161 b = false;
162 if (energies_) {
163 b = true;
164 stream & b;
165 stream & energies_->estrain_;
166 stream & energies_->eslip_;
167 } else
168 stream & b;
169 } else {
170 bool b(false);
171 stream & b;
172 if (b) {
173 if (!orientProps_)
174 orientProps_ = NEWC(orientProps());
175 stream & orientProps_->orient1_;
176 stream & orientProps_->orient2_;
177 stream & orientProps_->origNormal_;
178 }
179 stream & b;
180 if (b) {
181 if (!energies_)
182 energies_ = NEWC(Energies());
183 stream & energies_->estrain_;
184 stream & energies_->eslip_;
185 }
186 }
187
188 stream & inheritanceField_;
189 stream & effectiveTranslationalStiffness_;
190 stream & effectiveRotationalStiffness_;
191
192 if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() > 1)
193 stream & userArea_;
194
195 if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() > 2)
196 stream & fj_m_set_;
197
198 if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() > 3) {
199 stream & fj_cohres_;
200 stream & fj_resmode_;
201 }
202 }
203
204 void ContactModelFlatJoint::copy(const ContactModel *cm) {
205 ContactModelMechanical::copy(cm);
206 const ContactModelFlatJoint *in = dynamic_cast<const ContactModelFlatJoint*>(cm);
207 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
208 fj_nr(in->fj_nr());
209#ifdef THREED
210 fj_na(in->fj_na());
211#endif
212 fj_elem(in->fj_elem());
213 fj_kn(in->fj_kn());
214 fj_ks(in->fj_ks());
215 fj_fric(in->fj_fric());
216 fj_rmul(in->fj_rmul());
217 fj_gap0(in->fj_gap0());
218 fj_ten(in->fj_ten());
219 fj_coh(in->fj_coh());
220 fj_cohres(in->fj_cohres());
221 fj_fa(in->fj_fa());
222 fj_f(in->fj_f());
223 fj_m(in->fj_m());
224 fj_m_set(in->fj_m_set());
225 rmin(in->rmin());
226 rbar(in->rbar());
227 fj_resmode(in->fj_resmode());
228 atot(in->atot());
229 a_ = in->a_;
230 rBarl_ = in->rBarl_;
231 propsFixed(in->propsFixed());
232 mType(in->mType());
233 bmode_ = in->bmode_;
234 smode_ = in->smode_;
235 gap(in->gap());
236 theta(in->theta());
237#ifdef THREED
238 thetaM(in->thetaM());
239#endif
240 egap_ = in->egap_;
241 f_ = in->f_;
242 if (in->orientProps_) {
243 if (!orientProps_)
244 orientProps_ = NEWC(orientProps());
245 orientProps_->orient1_ = in->orientProps_->orient1_;
246 orientProps_->orient2_ = in->orientProps_->orient2_;
247 orientProps_->origNormal_ = in->orientProps_->origNormal_;
248 }
249 if (in->hasEnergies()) {
250 if (!energies_)
251 energies_ = NEWC(Energies());
252 estrain(in->estrain());
253 eslip(in->eslip());
254 }
255 userArea_ = in->userArea_;
256 inheritanceField(in->inheritanceField());
257 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
258 effectiveRotationalStiffness(in->effectiveRotationalStiffness());
259 }
260
261
262 QVariant ContactModelFlatJoint::getProperty(uint32 i,const IContact *con) const {
263 QVariant var;
264 switch (i) {
265 case kwFjNr : return fj_nr();
266 case kwFjElem : return fj_elem();
267 case kwFjKn : return fj_kn();
268 case kwFjKs : return fj_ks();
269 case kwFjFric : return fj_fric();
270 case kwFjEmod : {
271 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
272 if (c ==nullptr) return 0.0;
273 double rsum(0.0);
274 if (c->getEnd1Curvature().y())
275 rsum += 1.0/c->getEnd1Curvature().y();
276 if (c->getEnd2Curvature().y())
277 rsum += 1.0/c->getEnd2Curvature().y();
278 if (userArea_) {
279#ifdef THREED
280 rsum = std::sqrt(userArea_ / dPi);
281#else
282 rsum = userArea_ / 2.0;
283#endif
284 rsum += rsum;
285 }
286 return (fj_kn_ * rsum);
287 }
288 case kwFjKRatio : return (fj_ks_ == 0.0 ) ? 0.0 : (fj_kn_/fj_ks_);
289 case kwFjRmul : return fj_rmul();
290 case kwFjRadius : return rbar();
291 case kwFjGap0 : return fj_gap0();
292 case kwFjTen : return fj_ten();
293 case kwFjCoh : return fj_coh();
294 case kwFjFa : return fj_fa();
295 case kwFjF : var.setValue(fj_f()); return var;
296 case kwFjM : var.setValue(fj_m()); return var;
297 case kwFjState : return bmode_[fj_elem()-1];
298 case kwFjSlip : return smode_[fj_elem()-1];
299 case kwFjMType : return getmType();
300 case kwFjA : return a_[fj_elem()-1];
301 case kwFjEgap : return egap_[fj_elem()-1];
302 case kwFjGap : return gap().x();
303 case kwFjNstr : return -f_[fj_elem()-1].x() / a_[fj_elem()-1];
304 case kwFjSstr : return f_[fj_elem()-1].y() / a_[fj_elem()-1];
305 case kwFjSs : return tauC((-f_[fj_elem()-1].x() / a_[fj_elem()-1]),(bmode_[fj_elem()-1]==3));
306 case kwFjRelBr : var.setValue(DVect2(theta(),thetaM())); return var;
307 case kwFjCen : var.setValue(getRelElemPos(con,fj_elem()-1)); return var;
308#ifdef THREED
309 case kwFjNa : return fj_na();
310#endif
311 case kwFjTrack : var.setValue(orientProps_ ? true : false); return var;
312 case kwUserArea : return userArea_;
313 case kwFjCohRes : return fj_cohres();
314 case kwFjResMode: return fj_resmode();
315 }
316 assert(0);
317 return QVariant();
318 }
319
320 bool ContactModelFlatJoint::getPropertyGlobal(uint32 i) const {
321 switch (i) {
322 case kwFjF:
323 return false;
324 }
325 return true;
326 }
327
328 bool ContactModelFlatJoint::setProperty(uint32 i,const QVariant &v,IContact *c) {
329 bool ok(true);
330 switch (i) {
331 case kwFjNr: {
332 if (!propsFixed()) {
333 int val(v.toInt(&ok));
334 if (!ok || val < 1)
335 throw Exception("fj_nr must be an integer greater than 0.");
336 fj_nr(val);
337 if (fj_elem() > fj_n())
338 fj_elem(fj_n());
339 initVectors();
340 setAreaQuantities();
341 } else
342 throw Exception("fj_nr cannot be modified.");
343 return true;
344 }
345
346 case kwFjElem: {
347 int val(v.toInt(&ok));
348 if (!ok || val < 1 || val > fj_n())
349 throw Exception("fj_elem must be an integer between 1 and %1.",fj_n());
350 fj_elem(val);
351 return false;
352 }
353 case kwFjKn: {
354 double val(v.toDouble(&ok));
355 if (!ok || val<0.0)
356 throw Exception("fj_kn must be a positive double.");
357 fj_kn(val);
358 return true;
359 }
360 case kwFjKs: {
361 double val(v.toDouble(&ok));
362 if (!ok || val<0.0)
363 throw Exception("fj_ks must be a positive double.");
364 fj_ks(val);
365 return true;
366 }
367 case kwFjFric: {
368 double val(v.toDouble(&ok));
369 if (!ok || val<0.0)
370 throw Exception("fj_fric must be a positive double.");
371 fj_fric(val);
372 return false;
373 }
374 case kwFjRmul: {
375 if (!propsFixed()) {
376 double val(v.toDouble(&ok));
377 if (!ok || val<0.01)
378 throw Exception("fj_rmul must be a double greater than or equal to 0.01.");
379 fj_rmul(val);
380 setAreaQuantities();
381 return true;
382 } else
383 throw Exception("fj_rmul cannot be modified.");
384
385 return false;
386 }
387 case kwFjGap0: {
388 if (!propsFixed()) {
389 double val(v.toDouble(&ok));
390 if (!ok || val<0.0)
391 throw Exception("fj_gap0 must be a positive double.");
392 fj_gap0(val);
393 if (fj_gap0() > 0.0) {
394 for(int i=1; i<=fj_n(); ++i)
395 bondElem(i,false);
396 // surfaces are parallel w/ gap G
397 DVect temp(0.0);
398 temp.rx() = fj_gap0();
399 gap(temp);
400 theta(0.0);
401 }
402 } else
403 throw Exception("fj_gap0 cannot be modified.");
404 return true;
405 }
406 case kwFjTen: {
407 double val(v.toDouble(&ok));
408 if (!ok || val<0.0)
409 throw Exception("fj_ten must be a positive double.");
410 fj_ten(val);
411 return false;
412 }
413 case kwFjFa: {
414 double val(v.toDouble(&ok));
415 if (!ok || val<0.0)
416 throw Exception("fj_fa must be a positive double.");
417 fj_fa(val);
418 return false;
419 }
420 case kwFjCoh: {
421 double val(v.toDouble(&ok));
422 if (!ok || val<0.0)
423 throw Exception("fj_coh must be a positive double.");
424 fj_coh(val);
425 return false;
426 }
427 case kwFjA: {
428 double val(v.toDouble(&ok));
429 if (!ok || val<0.0)
430 throw Exception("fj_area must be a positive double.");
431 a_[fj_elem()-1] = val;
432 return false;
433 }
434 case kwFjNstr: {
435 double val(v.toDouble(&ok));
436 if (!ok || val<0.0)
437 throw Exception("fj_sigma must be a positive double.");
438 f_[fj_elem()-1].rx() = -val * a_[fj_elem()-1];
439 return false;
440 }
441 case kwFjSstr: {
442 double val(v.toDouble(&ok));
443 if (!ok || val<0.0)
444 throw Exception("fj_tau must be a positive double.");
445 f_[fj_elem()-1].ry() = val * a_[fj_elem()-1];
446 return false;
447 }
448#ifdef THREED
449 case kwFjNa: {
450 if (!propsFixed()) {
451 int val(v.toInt(&ok));
452 if (!ok || val < 1)
453 throw Exception("fj_na must be an integer greater than 0.");
454 fj_na(val);
455 if (fj_elem() > fj_n())
456 fj_elem(fj_n());
457 initVectors();
458 setAreaQuantities();
459 } else
460 throw Exception("fj_na cannot be modified.");
461 return true;
462 }
463#endif
464 case kwFjCen: {
465 if (!v.canConvert<DVect>())
466 throw Exception("fj_cen cannot be modified.");
467 DVect val(v.value<DVect>());
468 int el = fj_elem()-1;
469 setRelElemPos(c,el,val);
470 return false;
471 }
472 case kwFjTrack: {
473 if (!v.canConvert<bool>())
474 throw Exception("fj_track must be a boolean.");
475 bool b = v.toBool();
476 if (b) {
477 if (!orientProps_)
478 orientProps_ = NEWC(orientProps());
479 } else {
480 if (orientProps_) {
481 delete orientProps_;
482 orientProps_ = 0;
483 }
484 }
485 return true;
486 }
487 case kwUserArea: {
488 if (!v.canConvert<double>())
489 throw Exception("user_area must be a double.");
490 double val(v.toDouble());
491 if (val < 0.0)
492 throw Exception("Negative user_area not allowed.");
493 userArea_ = val;
494 propsFixed_ = false;
495 return true;
496 }
497 case kwFjCohRes: {
498 double val(v.toDouble(&ok));
499 if (!ok || val<0.0)
500 throw Exception("fj_cohres must be a positive double.");
501 fj_cohres(val);
502 return false;
503 }
504 case kwFjResMode: {
505 int val(v.toInt(&ok));
506 if (!ok || (val != 0 && val != 1))
507 throw Exception("fj_resmode must be 0 or 1.");
508 fj_resmode(val);
509 return false;
510 }
511 }
512 return false;
513 }
514
515 bool ContactModelFlatJoint::getPropertyReadOnly(uint32 i) const {
516 switch (i) {
517 case kwFjF:
518 case kwFjM:
519 case kwFjGap:
520 case kwFjRelBr:
521 case kwFjState:
522 case kwFjSlip:
523 case kwFjEgap:
524 case kwFjNstr:
525 case kwFjSstr:
526 case kwFjSs:
527 case kwFjRadius:
528 return true;
529 default:
530 break;
531 }
532 return false;
533 }
534
535 QString ContactModelFlatJoint::getMethodArguments(uint32 i) const {
536 switch (i) {
537 case kwBond:
538 case kwUnbond:
539 return "gap,element";
540 case KwDeformability:
541 return "emod,kratio";
542 case kwInitialize:
543 return "force,moment";
544 }
545 return QString();
546 }
547
548 bool ContactModelFlatJoint::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
549 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
550 bool bond(false);
551 switch (i) {
552 case kwBond:
553 bond = true;
554 [[fallthrough]];
555 case kwUnbond: {
556 int seg(0);
557 double mingap = -1.0 * limits<double>::max();
558 double maxgap = 0;
559 if (vl.size()==2) {
560 // The first is the gap
561 QVariant arg = vl.at(0);
562 if (!arg.isNull()) {
563 if (arg.canConvert<double>())
564 maxgap = vl.at(0).toDouble();
565 else if (arg.canConvert<DVect2>()) {
566 DVect2 value = vl.at(0).value<DVect2>();
567 mingap = value.minComp();
568 maxgap = value.maxComp();
569 } else
570 throw Exception("Argument %1 not recognized in method %2 in contact model %3.",vl.at(0),bond ? "bond":"unbond",getName());
571 }
572 arg = vl.at(1);
573 if (!arg.isNull()) {
574 seg = vl.at(1).toUInt();
575 if (seg < 1)
576 throw Exception("Element indices start at 1 in method %1 in contact model %2.",bond ? "bond":"unbond",getName());
577 if (seg > fj_n())
578 throw Exception("Element index %1 exceeds segments number (%2) in method %3 in contact model %4.",seg,fj_n(),bond ? "bond":"unbond",getName());
579 }
580 }
581 double gap = c->getGap();
582 if (gap >= mingap && gap <= maxgap) {
583 if (!seg) {
584 for(int iSeg=1; iSeg<=fj_n(); ++iSeg)
585 bondElem(iSeg,bond);
586 } else {
587 bondElem(seg,bond);
588 }
589 // If have installed bonds and tracking is enabled then set the contact normal appropriately
590 if (orientProps_) {
591 orientProps_->orient1_ = Quat::identity();
592 orientProps_->orient2_ = Quat::identity();
593 orientProps_->origNormal_ = toVect(con->getNormal());
594 }
595 }
596 return true;
597 }
598 case KwDeformability:
599 {
600 double emod;
601 double krat;
602 if (vl.at(0).isNull())
603 throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
604 emod = vl.at(0).toDouble();
605 if (emod<0.0)
606 throw Exception("Negative emod not allowed in contact model %1.",getName());
607 if (vl.at(1).isNull())
608 throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
609 krat = vl.at(1).toDouble();
610 if (krat<0.0)
611 throw Exception("Negative stiffness ratio not allowed in contact model %1.",getName());
612 double rsum(0.0);
613 if (c->getEnd1Curvature().y())
614 rsum += 1.0/c->getEnd1Curvature().y();
615 if (c->getEnd2Curvature().y())
616 rsum += 1.0/c->getEnd2Curvature().y();
617 if (userArea_) {
618#ifdef THREED
619 rsum = std::sqrt(userArea_ / dPi);
620#else
621 rsum = userArea_ / 2.0;
622#endif
623 rsum += rsum;
624 }
625 fj_kn_ = emod / rsum;
626 fj_ks_ = (krat == 0.0) ? 0.0 : fj_kn_ / krat;
627 return true;
628 }
629 case KwUpdateGeom: {
630 // go through and update the total area (atot) and the
631 // radius rbar
632 double at = 0.0;
633 for (int i=1; i<=fj_n(); ++i)
634 at += a_[i-1];
635 atot(at);
636 //get the equivalent radius
637#ifdef THREED
638 rbar(sqrt(at/dPi));
639#else
640 rbar(at/2.0);
641#endif
642 return true;
643 }
644 case kwArea: {
645 if (!userArea_) {
646 double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
647#ifdef THREED
648 userArea_ = rsq * rsq * dPi;
649#else
650 userArea_ = rsq * 2.0;
651#endif
652 }
653 return true;
654 }
655 case kwInitialize: {
656 DVect force;
657 DAVect moment;
658 if (vl.at(0).isNull())
659 throw Exception("Argument force must be specified with method initialize in contact model %1.",getName());
660 force = vl.at(0).value<DVect>();
661 if (vl.at(1).isNull())
662 throw Exception("Argument moment must be specified with method initialize in contact model %1.",getName());
663#ifdef THREED
664 moment = vl.at(1).value<DVect>();
665#else
666 moment.rz() = vl.at(1).toDouble();
667#endif
668 // Set the gap accordingly to get the correct force
669 setForce(force,con);
670 fj_m_set(moment);
671 return true;
672 }
673 }
674 return false;
675 }
676
677 double ContactModelFlatJoint::getEnergy(uint32 i) const {
678 double ret(0.0);
679 if (!energies_)
680 return ret;
681 switch (i) {
682 case kwEStrain: return energies_->estrain_;
683 case kwESlip: return energies_->eslip_;
684 }
685 assert(0);
686 return ret;
687 }
688
689 bool ContactModelFlatJoint::getEnergyAccumulate(uint32 i) const {
690 switch (i) {
691 case kwEStrain: return false;
692 case kwESlip: return true;
693 }
694 assert(0);
695 return false;
696 }
697
698 void ContactModelFlatJoint::setEnergy(uint32 i,const double &d) {
699 if (!energies_) return;
700 switch (i) {
701 case kwEStrain: energies_->estrain_ = d; return;
702 case kwESlip: energies_->eslip_ = d; return;
703 }
704 assert(0);
705 return;
706 }
707
708 bool ContactModelFlatJoint::validate(ContactModelMechanicalState *state,const double &) {
709 assert(state);
710 const IContactMechanical *c = state->getMechanicalContact();
711 assert(c);
712 // This presumes that one of the ends has a non-zero curvature
713 rmin(1.0/std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
714 if (userArea_) {
715#ifdef THREED
716 rmin(std::sqrt(userArea_ / dPi));
717#else
718 rmin(userArea_ / 2.0);
719#endif
720 }
721 if (!propsFixed()) {
722 setAreaQuantities();
723 mType(getmType());
724 }
725
726 // Initialize the tracking if not initialized
727 if (orientProps_ && orientProps_->origNormal_ == DVect(0.0)) {
728 orientProps_->origNormal_ = toVect(c->getContact()->getNormal());
729 orientProps_->orient1_ = Quat::identity();
730 orientProps_->orient2_ = Quat::identity();
731 }
732
733 if (state->trackEnergy_)
734 activateEnergy();
735
736 updateEffectiveStiffness(state);
737 return checkActivity(state->gap_);
738 }
739
740 void ContactModelFlatJoint::updateEffectiveStiffness(ContactModelMechanicalState *) {
741 DVect2 ret(fj_kn_,fj_ks_);
742 ret *= atot();
743 effectiveTranslationalStiffness(ret);
744#ifdef TWOD
745 effectiveRotationalStiffness(DAVect(fj_kn() * (2.0/3.0)*rbar()*rbar()*rbar()));
746#else
747 double piR4 = dPi * rbar() * rbar() * rbar() * rbar();
748 double t = fj_kn() * 0.25 * piR4;
749 effectiveRotationalStiffness(DAVect(fj_ks() * 0.50 * piR4,t,t));
750#endif
751 }
752
753 bool ContactModelFlatJoint::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
754 if (!propsFixed())
755 propsFixed(true);
756 timestep;
757 assert(state);
758
759 if (state->activated()) {
760 if (cmEvents_[fActivated] >= 0) {
761 auto c = state->getContact();
762 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
763 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
764 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
765 }
766 }
767
768 // Update the orientations
769 if (orientProps_) {
770 orientProps_->orient1_.increment(state->getMechanicalContact()->getEnd1Mechanical()->getAngVelocity()*timestep);
771 orientProps_->orient2_.increment(state->getMechanicalContact()->getEnd2Mechanical()->getAngVelocity()*timestep);
772 }
773
774#ifdef TWOD
775 // Translational increment in local coordinates
776 DVect del_U = state->relativeTranslationalIncrement_;
777 double del_theta = state->relativeAngularIncrement_.z();
778 gap(gap() + del_U); // in normal and shear direction in 2D
779 theta(theta() + del_theta);
780 double dSig=0.0, dTau=0.0;
781 double delX = 2*rbar() / fj_n();
782 double rbar2 = rbar() / fj_n();
783 DVect dFSum(0.0);
784 double dMSum = 0.0;
785 if (state->trackEnergy_) {
786 assert(energies_);
787 energies_->estrain_ = 0.0;
788 }
789 bool oneBonded = false;
790 for(int i=0; i<fj_n(); ++i) {
791 double g0 = gap((i )*delX);
792 double g1 = gap((i+1)*delX);
793 double gMid = 0.5*(g0 + g1);
794 if (bmode_[i] != 3 && gMid > 0) {
795 egap_[i] = gMid;
796 f_[i] = DVect(0.0);
797 continue;
798 }
799 dSig = computeSig(g0,g1,rbar2,a_[i],(bmode_[i]==3));
800 bool tensileBreak = false;
801 if (bmode_[i]==3) {
802 if (state->canFail_ && dSig >= fj_ten()) {
803 breakBond(i+1,1,state);
804 dSig = dTau = 0.0;
805 tensileBreak = true;
806 }
807 }
808 if (!tensileBreak) {
809 dTau = f_[i].y() / a_[i];
810 double dUse = delUse(egap_[i],gMid,(bmode_[i]==3),del_U.y());
811 double dtauP = dTau - fj_ks()*dUse;
812 double dtauPabs = abs(dtauP);
813 if (bmode_[i]==3) { // bonded
814 if (dtauPabs < tauC(dSig,true))
815 dTau = dtauP;
816 else {
817 if (state->canFail_) {
818 breakBond(i+1,2,state);
819 if (fj_resmode() == 0)
820 dSig = dTau = 0.0;
821 else
822 dTau = fj_cohres() - dSig * fj_fric();
823 }
824 }
825 } else { // unbonded
826 double dtauC = tauC(dSig,false);
827 if (dtauPabs <= dtauC) {
828 dTau = dtauP;
829 slipChange(i+1,false,state);
830 } else {
831 dTau = dtauP * ( dtauC / dtauPabs );
832 slipChange(i+1,true,state);
833 if (state->trackEnergy_) { energies_->eslip_ += dtauC*a_[i]*abs(dUse);}
834 }
835 }
836 }
837 oneBonded = true;
838 egap_[i] = gMid;
839 f_[i] = DVect(-dSig*a_[i],dTau*a_[i]);
840 dFSum += f_[i];
841 double m = computeM(g0,g1,rbar2,(bmode_[i]==3)) + fj_m_set().z()/fj_n();
842 dMSum += m - rBarl_[i]*f_[i].x();
843 if (state->trackEnergy_) {
844 if (fj_kn_) {
845 double ie = 2.0*rBarl_[i]*rBarl_[i]*rBarl_[i] / 3.0;
846 energies_->estrain_ += 0.5*(dSig*dSig*a_[i] + m*m/ie) / fj_kn_;
847 }
848 if (fj_ks_) {
849 energies_->estrain_ += 0.5 * dTau*dTau*a_[i] / fj_ks_;
850 }
851 }
852 }
853 //
854 fj_f(dFSum);
855 fj_m(DAVect(dMSum));
856 if (!oneBonded)
857 fj_m_set(DAVect(0.0));
858#else
859 CAxes localSys = state->getMechanicalContact()->getContact()->getLocalSystem();
860 DVect trans = state->relativeTranslationalIncrement_; // translation increment in local coordinates
861 DAVect ang = state->relativeAngularIncrement_; // rotational increment in local coordinates
862 DVect shear(0.0,trans.y(),trans.z());
863 DVect del_Us = localSys.toGlobal(shear); // In global coordinates
864 // What is the twist in global coordinates?
865 DVect del_Theta_t = localSys.e1()*ang.x();
866 theta_ += ang.y();
867 thetaM_ += ang.z();
868
869 gap(gap() + trans);
870 if (state->trackEnergy_) {
871 assert(energies_);
872 energies_->estrain_ = 0.0;
873 }
874 DVect force(0.0);
875 DAVect mom(0.0);
876 bool oneBonded = false;
877 for (int e=1,i=0; e<=fj_n(); ++e, ++i) {
878 double gBar1 = gap( rBarl_[i].x(),rBarl_[i].y());
879 if (!Bonded(e) && gBar1 > 0) {
880 egap_[i] = gBar1;
881 f_[i] = DVect(0.0);
882 continue;
883 }
884 DVect r = localSys.e2()*rBarl_[i].x() + localSys.e3()*rBarl_[i].y(); // location of element point
885 double sigBar_e = sigBar(e);
886 f_[i].rx() = -sigBar_e * a_[i]; // Step 1...
887 if (Bonded(e) && (sigBar_e >= fj_ten())) { // break bond in tension
888 if (state->canFail_) {
889 breakBond(e,1,state);
890 f_[i] = DVect(0.0);
891 }
892 } else {
893 DVect del_us = del_Us + (del_Theta_t & r); // In global - has the twist in there
894 double del_usl = delUse(egap_[i],gBar1,Bonded(e),(del_us | localSys.e2()));
895 double del_usm = delUse(egap_[i],gBar1,Bonded(e),(del_us | localSys.e3()));
896 double Fs_lP = f_[i].y() - fj_ks() * a_[i] * del_usl;
897 double Fs_mP = f_[i].z() - fj_ks() * a_[i] * del_usm;
898 double FsPMag = sqrt( Fs_lP*Fs_lP + Fs_mP*Fs_mP );
899 double tauBarP = FsPMag / a_[i];
900 if ( !Bonded(e) ) {
901 double tau_c = sigBar_e < 0.0 ? fj_cohres()-fj_fric()*sigBar_e : 0.0;
902 if ( tauBarP <= tau_c ) {
903 f_[i].ry() = Fs_lP;
904 f_[i].rz() = Fs_mP;
905 slipChange(e,false,state);
906 } else { // enforce sliding
907 double sFac = tau_c * a_[i] / FsPMag;
908 f_[i].ry() = Fs_lP * sFac;
909 f_[i].rz() = Fs_mP * sFac;
910 slipChange(e,true,state);
911 if (state->trackEnergy_) { energies_->eslip_ += tau_c*a_[i]*sqrt(del_usl*del_usl+del_usm*del_usm);}
912 }
913 } else { // Bonded(e)
914 double tau_c = fj_coh() - sigBar_e * tan(dDegrad*fj_fa());
915 if ( tauBarP <= tau_c ) {
916 f_[i].ry() = Fs_lP;
917 f_[i].rz() = Fs_mP;
918 } else { // break bond in shear
919 if (state->canFail_) {
920 breakBond(e,2,state);
921 if (fj_resmode() == 0)
922 f_[i] = DVect(0.0);
923 else {
924 double newForce = fj_cohres() - sigBar_e * fj_fric();
925 if (!userArea_)
926 newForce *= a_[i];
927 else
928 newForce *= userArea_ / fj_n();
929 newForce /= std::sqrt(f_[i].y()*f_[i].y() + f_[i].z()*f_[i].z());
930 f_[i].ry() *= newForce;
931 f_[i].rz() *= newForce;
932 }
933 }
934 }
935 }
936 }
937 oneBonded = true;
938 force += f_[i];
939 mom += (localSys.toLocal(r) & f_[i]) + fj_m_set()/fj_n();
940 egap_[i] = gBar1;
941 if (state->trackEnergy_) {
942 energies_->estrain_ += computeStrainEnergy(e);
943 }
944 }
945 fj_f(force);
946 fj_m(mom);
947 if (!oneBonded)
948 fj_m_set(DAVect(0.0));
949#endif
950 assert(fj_f_ == fj_f_);
951 return checkActivity(0.0);
952 }
953
954 bool ContactModelFlatJoint::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
955 // Account for thermal expansion in incremental mode
956 if (ts->gapInc_ == 0.0) return false;
957 DVect dg(0.0);
958 dg.rx() = ts->gapInc_;
959 gap(gap() + dg);
960 return true;
961 }
962
963 bool ContactModelFlatJoint::coupling(ContactModelMechanicalState* , ContactModelThermalState*,ContactModelFluidState* fs, const double&) {
964 // Account for strain-induced pressure increase in fluid contact
965 if (!fs)
966 return false;
967 double dp = fs->dp_geom_;
968 DVect force = DVect(0.0);
969 for (int e = 1, i = 0; e <= fj_n(); ++e, ++i) {
970 DVect finc(0.0);
971 finc.rx() = dp * a_[i];
972 f_[i] -= finc;
973 force += f_[i];
974 }
975 fj_f(force);
976 return true;
977 }
978
979
980 void ContactModelFlatJoint::setAreaQuantities() {
981 rbar(fj_rmul() * rmin());
982#ifdef TWOD
983 atot(2.0 * rbar());
984 double v = atot()/fj_n();
985 for (int i=1; i<=fj_n(); ++i) {
986 a_[i-1] = v;
987 rBarl_[i-1] = rbar() * (double(-2*i + 1 + fj_n()) / fj_n());
988 }
989#else
990 atot(dPi * rbar() * rbar());
991 double del_r = rbar() / fj_nr();
992 double del_al = 2.0*dPi / fj_na();
993 double fac = 2.0/3.0;
994 for (int i=0; i < fj_n(); ++i) {
995 double dVal = i / fj_na();
996 int I = (int)floor( dVal );
997 int J = i - I*fj_na();
998 double r1 = I * del_r;
999 double r2 = (I + 1) * del_r;
1000 double al1 = J * del_al;
1001 double al2 = (J + 1) * del_al;
1002 a_[i] = 0.5 * (al2 - al1) * (r2*r2 - r1*r1);
1003 rBarl_[i] = DVect2(((sin(al2) - sin(al1)) / (al2 - al1))*((r2*r2*r2 - r1*r1*r1)/(r2*r2 - r1*r1)),
1004 ((cos(al1) - cos(al2)) / (al2 - al1))*((r2*r2*r2 - r1*r1*r1)/(r2*r2 - r1*r1)))*fac;
1005 }
1006#endif
1007 updateEffectiveStiffness(0);
1008 }
1009
1010 DVect ContactModelFlatJoint::getRelElemPos(const IContact* c,int i) const {
1011 DVect ret(0.0);
1012 if (c) {
1013 ret = c->getPosition();
1014 CAxes localSys = c->getLocalSystem();
1015#ifdef TWOD
1016 ret += localSys.e2()*rBarl_[i];
1017#else
1018 ret += localSys.e2()*rBarl_[i].x() + localSys.e3()*rBarl_[i].y();
1019#endif
1020 }
1021 return ret;
1022 }
1023
1024 void ContactModelFlatJoint::setRelElemPos(const IContact* c,int i,const DVect &pos) {
1025 // pos is a position in space in global coordinates
1026 propsFixed(true);
1027 if (c) {
1028 // project onto the plane
1029 DVect cp = c->getPosition();
1030 DVect norm = toVect(c->getNormal());
1031 double sd = norm|(cp - pos);
1032 // np is the point on the plane
1033 DVect np = pos+norm*sd;
1034 np = np-cp;
1035 CAxes localSys = c->getLocalSystem();
1036 np = localSys.toLocal(np);
1037#ifdef TWOD
1038 rBarl_[i] = np.y();
1039#else
1040 rBarl_[i] = DVect2(np.y(),np.z());
1041#endif
1042 }
1043 }
1044
1045 int ContactModelFlatJoint::getmType() const {
1046 if (propsFixed()) return mType();
1047 //
1048 if (fj_gap0() > 0.0) return 2;
1049 //
1050 // If we get to here, then G == 0.0.
1051 bool AllBonded = true;
1052 bool AllSlit = true;
1053 for(int i=0; i<fj_n(); ++i) {
1054 if (bmode_[i] != 3) AllBonded = false;
1055 else AllSlit = false;
1056 }
1057 if (AllBonded) return 1;
1058 if (AllSlit) return 3;
1059 //
1060 return 4;
1061 }
1062
1063 void ContactModelFlatJoint::bondElem(int iSeg,bool bBond ) {
1064 if (bBond) {
1065 if (fj_gap0() == 0.0) {
1066 bmode_[iSeg-1] = 3;
1067 } else
1068 bmode_[iSeg-1] = 0;
1069 } else
1070 bmode_[iSeg-1] = 0;
1071 }
1072
1073 void ContactModelFlatJoint::breakBond(int iSeg,int fmode,ContactModelMechanicalState *state) {
1074 bmode_[iSeg-1] = fmode;
1075 if (cmEvents_[fBondBreak] >= 0) {
1076 auto c = state->getContact();
1077 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1078 fish::Parameter((qint64)iSeg),
1079 fish::Parameter((qint64)fmode),
1080 fish::Parameter(computeStrainEnergy(iSeg))
1081 };
1082 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1083 fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1084 }
1085 if (!isBonded() && cmEvents_[fBroken] >= 0) {
1086 auto c = state->getContact();
1087 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
1088 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1089 fi->setCMFishCallArguments(c,arg,cmEvents_[fBroken]);
1090 }
1091 }
1092
1093 void ContactModelFlatJoint::slipChange(int iSeg,bool smode,ContactModelMechanicalState *state) {
1094 bool emitEvent = false;
1095 if (smode) {
1096 if (!smode_[iSeg-1]) {
1097 emitEvent = true;
1098 smode_[iSeg-1] = smode;
1099 }
1100 } else {
1101 if (smode_[iSeg-1]) {
1102 emitEvent = true;
1103 smode_[iSeg-1] = smode;
1104 }
1105 }
1106 if (emitEvent && cmEvents_[fSlipChange] >= 0) {
1107 auto c = state->getContact();
1108 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1109 fish::Parameter((qint64)iSeg),
1110 fish::Parameter(smode) };
1111 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1112 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1113 }
1114 }
1115
1116 double ContactModelFlatJoint::tauC(const double &dSig,bool bBonded) const {
1117 if (bBonded) return (fj_coh() + (tan(dDegrad*fj_fa()) * (-dSig)) );
1118 else
1119 return (dSig < 0.0 ? fj_cohres() - fj_fric() * dSig : 0.0 );
1120 }
1121
1122#ifdef THREED
1123 double ContactModelFlatJoint::xi(int comp) const {
1124 if (comp == 1) return abs(theta_) <= 1e-12 ? 0.0 : theta_/thbMag();
1125 else return abs(thetaM_) <= 1e-12 ? 0.0 : thetaM_/thbMag();
1126 }
1127#endif
1128
1129 void ContactModelFlatJoint::initVectors() {
1130 a_.resize(fj_n());
1131 rBarl_.resize(fj_n());
1132 bmode_.resize(fj_n());
1133 smode_.resize(fj_n());
1134 egap_.resize(fj_n());
1135 f_.resize(fj_n());
1136 for (int i=0; i<fj_n(); ++i) {
1137 a_[i] = egap_[i] = 0.0;
1138#ifdef THREED
1139 rBarl_[i] = DVect2(0.0);
1140#else
1141 rBarl_[i] = 0.0;
1142#endif
1143 f_[i] = DVect(0.0);
1144 bmode_[i] = 0;
1145 smode_[i] = false;
1146 }
1147 }
1148
1149#ifdef TWOD
1150 double ContactModelFlatJoint::gap(const double &x) const {
1151 return gap().x() + theta()*(x - rbar());
1152 }
1153#else
1154 double ContactModelFlatJoint::gap(const double &r_l,const double &r_m ) const {
1155 return gap().x() + ( r_m*xi(1) - r_l*xi(2) ) * thbMag();
1156 }
1157
1158 double ContactModelFlatJoint::sigBar(int e) const {
1159 if (!Bonded(e)&& gap(rBarl_[e-1].x(),rBarl_[e-1].y()) >= 0.0)
1160 return 0.0;
1161 else
1162 return fj_kn() * gap(rBarl_[e-1].x(),rBarl_[e-1].y());
1163 }
1164
1165 double ContactModelFlatJoint::tauBar(int e) const {
1166 return a_[e-1] <= 1e-12 ?
1167 0.0 : sqrt(f_[e-1].y()*f_[e-1].y() + f_[e-1].z()*f_[e-1].z())/a_[e-1] ;
1168 }
1169
1170#endif
1171
1172 double ContactModelFlatJoint::computeStrainEnergy(int e) const {
1173 double ret(0.0);
1174 int i = e - 1;
1175#ifdef TWOD
1176 double delX = 2 * rbar() / fj_n();
1177 double g0 = gap((i)*delX);
1178 double g1 = gap((i + 1)*delX);
1179 double rbar2 = rbar() / fj_n();
1180 double dSig = computeSig(g0, g1, rbar2, a_[i], (bmode_[i] == 3));
1181 double m = computeM(g0, g1, rbar2, (bmode_[i] == 3));
1182 double dTau = f_[i].y() / a_[i]; // only valid before failure
1183 if (fj_kn_) {
1184 double ie = 2.0*rBarl_[i] * rBarl_[i] * rBarl_[i] / 3.0;
1185 ret += 0.5*(dSig*dSig*a_[i] + m * m / ie) / fj_kn_;
1186 }
1187 if (fj_ks_) {
1188 ret += 0.5 * dTau*dTau*a_[i] / fj_ks_;
1189 }
1190#else
1191 if (fj_kn_) {
1192 ret += 0.5*(sigBar(e)*sigBar(e)*a_[i]) / fj_kn_;
1193 }
1194 if (fj_ks_) {
1195 ret += 0.5 * (f_[i].y()*f_[i].y() + f_[i].z()*f_[i].z()) / (fj_ks_*a_[i]);
1196 }
1197#endif
1198 return ret;
1199 }
1200
1201 double ContactModelFlatJoint::computeSig(const double &g0,const double &g1,const double &rbar,
1202 const double &dA,bool bBonded ) const {
1203 double gTerm;
1204 switch (getCase(g0, g1)) {
1205 case 1:
1206 if (bBonded) gTerm = (g0 + g1);
1207 else if (g0 < 0.0) gTerm = -( g0*g0 / (g1 - g0) );
1208 else gTerm = ( g1*g1 / (g1 - g0) );
1209 break;
1210 case 2:
1211 if (bBonded) gTerm = (g0 + g1);
1212 else gTerm = 0.0;
1213 break;
1214 case 3:
1215 gTerm = (g0 + g1);
1216 break;
1217 default: gTerm = 0.0; break;
1218 }
1219 return (fj_kn() * gTerm * rbar) / dA;
1220 }
1221
1222 double ContactModelFlatJoint::computeM(const double &g0,const double &g1,const double &rbar,
1223 bool bBonded) const {
1224 double gTerm;
1225 switch (getCase(g0,g1)) {
1226 case 1:
1227 if (bBonded) gTerm = -((g1 - g0) / 3.0);
1228 else if (g0 < 0.0) gTerm = g0*g0*(g0 - 3.0*g1) / (3.0*(g1-g0)*(g1-g0));
1229 else gTerm = -(((g1-g0)*(g1-g0)*(g1-g0) + g0*g0*(g0 - 3.0*g1))
1230 / (3.0*(g1-g0)*(g1-g0)));
1231 break;
1232 case 2:
1233 if (bBonded) gTerm = -((g1 - g0) / 3.0);
1234 else gTerm = 0.0;
1235 break;
1236 case 3:
1237 gTerm = -((g1 - g0) / 3.0);
1238 break;
1239 default: gTerm = 0.0; break;
1240 }
1241 return fj_kn() * gTerm * rbar*rbar;
1242 }
1243
1244 int ContactModelFlatJoint::getCase(const double &g0,const double &g1) const {
1245 if (g0 * g1 < 0.0) // Case 1: gap changes sign
1246 return 1;
1247 else if (g0 >= 0.0 && g1 >= 0.0) // Case 2: gap remains positive or zero
1248 return 2;
1249 else // Case 3: gap remains negative
1250 return 3;
1251 }
1252
1253 double ContactModelFlatJoint::delUse(const double &gapStart,const double &gapEnd,bool bBonded,
1254 const double &delUs) const {
1255 if ( bBonded ) return delUs;
1256 if ( gapStart <= 0.0 ) {
1257 if ( gapEnd <= 0.0 )
1258 return delUs;
1259 else { // gapEnd > 0.0
1260 double xi = -gapStart / (gapEnd - gapStart);
1261 return delUs * xi;
1262 }
1263 } else { // gapStart > 0.0
1264 if ( gapEnd >= 0.0 )
1265 return 0.0;
1266 else { // gapEnd < 0.0
1267 double xi = -gapStart / (gapEnd - gapStart);
1268 return delUs * (1.0 - xi);
1269 }
1270 }
1271 }
1272
1273 bool ContactModelFlatJoint::checkActivity(const double &inGap) {
1274 // If any subcontact is bonded return true
1275 FOR(it,bmode_) if ((*it) == 3)
1276 return true;
1277 // If the normal gap is less than 2*rbar return true
1278 if (gap().x() < 2.0*rbar())
1279 return true;
1280 // check to see if there is overlap (based on the initial gap) to reset activity if the contact has been previously deactivated
1281 if (inGap < 0) {
1282 // reset the relative rotation
1283 theta(0.0);
1284#ifdef THREED
1285 thetaM(0.0);
1286#endif
1287 // set the gap to be the current gap, removing the shear displacement
1288 DVect inp(inGap,0.0);
1289 gap(inp);
1290 return true;
1291 }
1292 return false;
1293 }
1294
1295 void ContactModelFlatJoint::setForce(const DVect &v,IContact *) {
1296 fj_f_ = v;
1297 DVect df = v / f_.size();
1298 for (int i=0; i<f_.size(); ++i)
1299 f_[i] = df;
1300 // Set gap consistent with normal force
1301 double at = userArea_;
1302 if (!userArea_) {
1303 for (int i = 1; i <= fj_n(); ++i)
1304 at += a_[i - 1];
1305 }
1306 gap_.rx() = -1.0 * v.x() / (fj_kn_ * at);
1307 }
1308
1309 void ContactModelFlatJoint::getSphereList(const IContact *con,std::vector<DVect> *pos,std::vector<double> *rad,std::vector<double> *val) {
1310 assert(pos);
1311 assert(rad);
1312 assert(val);
1313 if (!orientProps_)
1314 return;
1315 // find minimal radii for end1
1316 double br = convert_getcast<IContactMechanical>(con)->getEnd1Curvature().y();
1317 if (br) {
1318 const IPiece *p = con->getEnd1();
1319 FArray<const IContact*> arr;
1320 p->getContactList(&arr);
1321 double maxgap = 0.0;
1322 FOR(ic,arr) {
1323 const IContactMechanical *mc = convert_getcast<IContactMechanical>(*ic);
1324 const IContactModelMechanical *mcm = mc->getModelMechanical();
1325 if (mcm->getContactModel()->getIndex() == ContactModelFlatJoint::getIndex()) {
1326 const ContactModelFlatJoint *in = dynamic_cast<const ContactModelFlatJoint*>(mcm);
1327 maxgap = std::max<double>(maxgap,in->gap().x()- mc->getGap());
1328 }
1329 }
1330 br = 1.0 / br - 0.5*maxgap;
1331 const IContactMechanical *mc = convert_getcast<IContactMechanical>(con);
1332 pos->push_back(convert_getcast<IPieceMechanical>(mc->getEnd1())->getPosition());
1333 rad->push_back(br);
1334 val->push_back(mc->getEnd1()->getIThing()->getID());
1335 }
1336
1337 // Give the end2 sphere - bummer
1338 br = convert_getcast<IContactMechanical>(con)->getEnd2Curvature().y();
1339 if (br) {
1340 const IPiece *p = con->getEnd2();
1341 FArray<const IContact*> arr;
1342 p->getContactList(&arr);
1343 double maxgap = 0.0;
1344 FOR(ic,arr) {
1345 const IContactMechanical *mc = convert_getcast<IContactMechanical>(*ic);
1346 const IContactModelMechanical *mcm = mc->getModelMechanical();
1347 if (mcm->getContactModel()->getIndex() == ContactModelFlatJoint::getIndex()) {
1348 const ContactModelFlatJoint *in = dynamic_cast<const ContactModelFlatJoint*>(mcm);
1349 maxgap = std::max<double>(maxgap,in->gap().x()- mc->getGap());
1350 }
1351 }
1352 br = 1.0 / br - 0.5*maxgap;
1353 const IContactMechanical *mc = convert_getcast<IContactMechanical>(con);
1354 pos->push_back(convert_getcast<IPieceMechanical>(mc->getEnd2())->getPosition());
1355 rad->push_back(br);
1356 val->push_back(mc->getEnd2()->getIThing()->getID());
1357 }
1358 }
1359
1360#ifdef THREED
1361
1362 void ContactModelFlatJoint::getDiskList(const IContact *con,std::vector<DVect> *pos,std::vector<DVect> *normal,std::vector<double> *radius,std::vector<double> *val) {
1363 assert(pos);
1364 assert(normal);
1365 assert(radius);
1366 assert(val);
1367 if (!orientProps_)
1368 return;
1369 // plot the contact plane right in the middle of the 2 normals
1370 double rad = fj_rmul()*rmin();
1371 DVect axis1 = orientProps_->orient1_.rotate(orientProps_->origNormal_);
1372 DVect axis2 = orientProps_->orient2_.rotate(orientProps_->origNormal_);
1373 DVect norm = ((axis1.unit()+axis2.unit())*0.5).unit();
1374 pos->push_back(con->getPosition());
1375 normal->push_back(norm);
1376 radius->push_back(rad);
1377 const IContactMechanical *mc = convert_getcast<IContactMechanical>(con);
1378 val->push_back(mc->getLocalForce().mag());
1379 }
1380
1381#endif
1382
1383 void ContactModelFlatJoint::getCylinderList(const IContact *con,std::vector<DVect> *bot,std::vector<DVect> *top,std::vector<double> *radlow,std::vector<double> *radhi,std::vector<double> *val) {
1384 assert(bot);
1385 assert(top);
1386 assert(radlow);
1387 assert(radhi);
1388 assert(val);
1389 if (!orientProps_)
1390 return;
1391 const IContactMechanical *mc = convert_getcast<IContactMechanical>(con);
1392 double br = mc->getEnd1Curvature().y(), br2 = mc->getEnd2Curvature().y();
1393 if (userArea_) {
1394#ifdef THREED
1395 br = std::sqrt(userArea_ / dPi);
1396#else
1397 br = userArea_ / 2.0;
1398#endif
1399 br = 1. / br;
1400 br2 = br;
1401 }
1402
1403 double cgap = mc->getGap();
1404 if (br > 0 && br2 > 0) {
1405 br = 1.0 / br;
1406 br2 = 1.0 / br2;
1407 double rad = fj_rmul()*rmin();
1408 DVect bp = convert_getcast<IPieceMechanical>(mc->getEnd1())->getPosition();
1409 DVect axis = orientProps_->orient1_.rotate(orientProps_->origNormal_);
1410 bot->push_back(axis.unit()*(br-0.5*(gap().x()- cgap))+bp);
1411 top->push_back(bp);
1412 radlow->push_back(rad);
1413 radhi->push_back(0.0);
1414 val->push_back(mc->getEnd1()->getIThing()->getID());
1415 bp = convert_getcast<IPieceMechanical>(mc->getEnd2())->getPosition();
1416 axis = orientProps_->orient2_.rotate(orientProps_->origNormal_);
1417 bot->push_back(axis.unit()*(br2-0.5*(gap().x()-cgap))*(-1.0)+bp);
1418 top->push_back(bp);
1419 radlow->push_back(rad);
1420 radhi->push_back(0.0);
1421 val->push_back(mc->getEnd2()->getIThing()->getID());
1422 }
1423 }
1424
1425 DVect ContactModelFlatJoint::getForce(const IContactMechanical *) const {
1426 DVect ret(fj_f_);
1427 return ret;
1428 }
1429
1430 DAVect ContactModelFlatJoint::getMomentOn1(const IContactMechanical *c) const {
1431 DVect force = getForce(c);
1432 DAVect ret(fj_m_);
1433 c->updateResultingTorqueOn1Local(force,&ret);
1434 return ret;
1435 }
1436
1437 DAVect ContactModelFlatJoint::getMomentOn2(const IContactMechanical *c) const {
1438 DVect force = getForce(c);
1439 DAVect ret(fj_m_);
1440 c->updateResultingTorqueOn2Local(force,&ret);
1441 return ret;
1442 }
1443
1444} // namespace itascaxd
1445
1446// EoF
Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Aug 13, 2024 |