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