Rolling Resistance Linear Contact Model Implementation
See this page for the documentation of this contact model.
contactmodelrrlinear.h
1#pragma once
2// contactmodelrrlinear.h
3
4#include "contactmodel/src/contactmodelmechanical.h"
5
6#ifdef rrlinear_LIB
7# define rrlinear_EXPORT EXPORT_TAG
8#elif defined(NO_MODEL_IMPORT)
9# define rrlinear_EXPORT
10#else
11# define rrlinear_EXPORT IMPORT_TAG
12#endif
13
14namespace cmodelsxd {
15 using namespace itasca;
16
17 class ContactModelRRLinear : public ContactModelMechanical {
18 public:
19 // Constructor: Set default values for contact model properties.
20 rrlinear_EXPORT ContactModelRRLinear();
21 rrlinear_EXPORT ContactModelRRLinear(const ContactModelRRLinear &) noexcept;
22 rrlinear_EXPORT const ContactModelRRLinear & operator=(const ContactModelRRLinear &);
23 rrlinear_EXPORT void addToStorage(poly::vector<ContactModelMechanical> *s,int ww = -1) override;
24 // Destructor, called when contact is deleted: free allocated memory, etc.
25 rrlinear_EXPORT virtual ~ContactModelRRLinear();
26 // Contact model name (used as keyword for commands and FISH).
27 string getName() const override { return "rrlinear"; }
28 // The index provides a quick way to determine the type of contact model.
29 // Each type of contact model in PFC must have a unique index; this is assigned
30 // by PFC when the contact model is loaded. This index should be set to -1
31 void setIndex(int i) override { index_=i;}
32 int getIndex() const override {return index_;}
33 // Contact model version number (e.g., MyModel05_1). The version number can be
34 // accessed during the save-restore operation (within the archive method,
35 // testing {stream.getRestoreVersion() == getMinorVersion()} to allow for
36 // future modifications to the contact model data structure.
37 uint32 getMinorVersion() const override;
38 // Copy the state information to a newly created contact model.
39 // Provide access to state information, for use by copy method.
40 void copy(const ContactModel *c) override;
41 // Provide save-restore capability for the state information.
42 void archive(ArchiveStream &) override;
43 // Enumerator for the properties.
44 enum PropertyKeys {
45 kwKn=1
46 , kwKs
47 , kwFric
48 , kwLinF
49 , kwLinS
50 , kwLinMode
51 , kwRGap
52 , kwEmod
53 , kwKRatio
54 , kwDpNRatio
55 , kwDpSRatio
56 , kwDpMode
57 , kwDpF
58 , kwResFric
59 , kwResMoment
60 , kwResS
61 , kwResKr
62 , kwUserArea
63 };
64 // Contact model property names in a comma separated list. The order corresponds with
65 // the order of the PropertyKeys enumerator above. One can visualize any of these
66 // properties in PFC automatically.
67 string getProperties() const override {
68 return "kn"
69 ",ks"
70 ",fric"
71 ",lin_force"
72 ",lin_slip"
73 ",lin_mode"
74 ",rgap"
75 ",emod"
76 ",kratio"
77 ",dp_nratio"
78 ",dp_sratio"
79 ",dp_mode"
80 ",dp_force"
81 ",rr_fric"
82 ",rr_moment"
83 ",rr_slip"
84 ",rr_kr"
85 ",user_area" ;
86 }
87 // Enumerator for the energies.
88 enum EnergyKeys {
89 kwEStrain=1
90 , kwERRStrain
91 , kwESlip
92 , kwERRSlip
93 , kwEDashpot
94 };
95 // Contact model energy names in a comma separated list. The order corresponds with
96 // the order of the EnergyKeys enumerator above.
97 string getEnergies() const override {
98 return "energy-strain"
99 ",energy-rrstrain"
100 ",energy-slip"
101 ",energy-rrslip"
102 ",energy-dashpot";
103 }
104 // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
105 double getEnergy(uint32 i) const override;
106 // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1)
107 // returns wther or not the estrain energy is accumulated which is false).
108 bool getEnergyAccumulate(uint32 i) const override;
109 // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
110 void setEnergy(uint32 i,const double &d) override; // Base 1
111 // Activate the energy. This is only called if the energy tracking is enabled.
112 void activateEnergy() override { if (energies_) return; energies_ = NEW Energies();}
113 // Returns whether or not the energy tracking has been enabled for this contact.
114 bool getEnergyActivated() const override {return (energies_ != 0);}
115
116 // Enumerator for contact model related FISH callback events.
117 enum FishCallEvents {
118 fActivated=0
119 ,fSlipChange
120 };
121 // Contact model FISH callback event names in a comma separated list. The order corresponds with
122 // the order of the FishCallEvents enumerator above.
123 string getFishCallEvents() const override {
124 return
125 "contact_activated"
126 ",slip_change";
127 }
128
129 // Return the specified contact model property.
130 base::Property getProperty(uint32 i,const IContact *) const override;
131 // The return value denotes whether or not the property corresponds to the global
132 // or local coordinate system (TRUE: global system, FALSE: local system). The
133 // local system is the contact-plane system (nst) defined as follows.
134 // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
135 // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
136 // and tc are unit vectors in directions of the nst axes.
137 // This is used when rendering contact model properties that are vectors.
138 bool getPropertyGlobal(uint32 i) const override;
139 // Set the specified contact model property, ensuring that it is of the correct type
140 // and within the correct range --- if not, then throw an exception.
141 // The return value denotes whether or not the update has affected the timestep
142 // computation (by having modified the translational or rotational tangent stiffnesses).
143 // If true is returned, then the timestep will be recomputed.
144 bool setProperty(uint32 i,const base::Property &v,IContact *) override;
145 // The return value denotes whether or not the property is read-only
146 // (TRUE: read-only, FALSE: read-write).
147 bool getPropertyReadOnly(uint32 i) const override;
148
149 // The return value denotes whether or not the property is inheritable
150 // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
151 // the endPropertyUpdated method.
152 bool supportsInheritance(uint32 i) const override;
153 // Return whether or not inheritance is enabled for the specified property.
154 bool getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i); return (inheritanceField_ & mask) ? true : false; }
155 // Set the inheritance flag for the specified property.
156 void setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i); if (b) inheritanceField_ |= mask; else inheritanceField_ &= ~mask; }
157
158 // Enumerator for contact model methods.
159 enum MethodKeys { kwDeformability=1,kwArea};
160 // Contact model methoid names in a comma separated list. The order corresponds with
161 // the order of the MethodKeys enumerator above.
162 string getMethods() const override { return "deformability,area";}
163 // Return a comma seprated list of the contact model method arguments (base 1).
164 string getMethodArguments(uint32 i) const override;
165 // Set contact model method arguments (base 1).
166 // The return value denotes whether or not the update has affected the timestep
167 // computation (by having modified the translational or rotational tangent stiffnesses).
168 // If true is returned, then the timestep will be recomputed.
169 bool setMethod(uint32 i,const std::vector<base::Property> &vl,IContact *con=0) override;
170
171 // Prepare for entry into ForceDispLaw. The validate function is called when:
172 // (1) the contact is created, (2) a property of the contact that returns a true via
173 // the setProperty method has been modified and (3) when a set of cycles is executed
174 // via the {cycle N} command.
175 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
176 bool validate(ContactModelMechanicalState *state,const double ×tep) override;
177 // The endPropertyUpdated method is called whenever a surface property (with a name
178 // that matches an inheritable contact model property name) of one of the contacting
179 // pieces is modified. This allows the contact model to update its associated
180 // properties. The return value denotes whether or not the update has affected
181 // the time step computation (by having modified the translational or rotational
182 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
183 bool endPropertyUpdated(const string &name,const IContactMechanical *c) override;
184 // The forceDisplacementLaw function is called during each cycle. Given the relative
185 // motion of the two contacting pieces (via
186 // state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
187 // state->relativeAngularIncrement_ (Dtt, Dtbs, Dtbt)
188 // Ddn : relative normal-displacement increment, Ddn > 0 is opening
189 // Ddss : relative shear-displacement increment (s-axis component)
190 // Ddst : relative shear-displacement increment (t-axis component)
191 // Dtt : relative twist-rotation increment
192 // Dtbs : relative bend-rotation increment (s-axis component)
193 // Dtbt : relative bend-rotation increment (t-axis component)
194 // The relative displacement and rotation increments:
195 // Dd = Ddn*nc + Ddss*sc + Ddst*tc
196 // Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
197 // where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
198 // [see {Table 1: Contact State Variables} in PFC Model Components:
199 // Contacts and Contact Models: Contact Resolution]
200 // ) and the contact properties, this function must update the contact force and
201 // moment.
202 // The force_ is acting on piece 2, and is expressed in the local coordinate system
203 // (defined in getPropertyGlobal) such that the first component positive denotes
204 // compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
205 // in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to
206 // get the total moment.
207 // The return value indicates the contact activity status (TRUE: active, FALSE:
208 // inactive) during the next cycle.
209 // Additional information:
210 // * If state->activated() is true, then the contact has just become active (it was
211 // inactive during the previous time step).
212 // * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
213 // the forceDispLaw handle the case of {state->canFail_ == true}.
214 bool forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) override;
215 bool thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
216 // The getEffectiveXStiffness functions return the translational and rotational
217 // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
218 // the translational tangent shear stiffness is zero (but this stiffness reduction
219 // is typically ignored when computing a stable time step). If the contact model
220 // includes a dashpot, then the translational stiffnesses must be increased (see
221 // Potyondy (2009)).
222 // [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
223 // Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
224 // December 7, 2009.]
225 DVect2 getEffectiveTranslationalStiffness() const override { return effectiveTranslationalStiffness_; }
226 DAVect getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness_;}
227
228 // Return a new instance of the contact model. This is used in the CMAT
229 // when a new contact is created.
230 ContactModelRRLinear *clone() const override { return NEW ContactModelRRLinear(); }
231 // The getActivityDistance function is called by the contact-resolution logic when
232 // the CMAT is modified. Return value is the activity distance used by the
233 // checkActivity function.
234 double getActivityDistance() const override {return rgap_;}
235 // The isOKToDelete function is called by the contact-resolution logic when...
236 // Return value indicates whether or not the contact may be deleted.
237 // If TRUE, then the contact may be deleted when it is inactive.
238 // If FALSE, then the contact may not be deleted (under any condition).
239 bool isOKToDelete() const override { return !isBonded(); }
240 // Zero the forces and moments stored in the contact model. This function is called
241 // when the contact becomes inactive.
242 void resetForcesAndMoments() override {
243 lin_F(DVect(0.0));
244 dp_F(DVect(0.0));
245 res_M(DAVect(0.0));
246 if (energies_) {
247 energies_->estrain_ = 0.0;
248 energies_->errstrain_ = 0.0;
249 }
250 }
251 void setForce(const DVect &v,IContact *c) override;
252 void setArea(const double &d) override { userArea_ = d; }
253 double getArea() const override { return userArea_; }
254 // The checkActivity function is called by the contact-resolution logic when...
255 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
256 bool checkActivity(const double &gap) override { return gap <= rgap_; }
257 // Returns the sliding state (FALSE is returned if not implemented).
258 bool isSliding() const override { return lin_S_; }
259 // Returns the bonding state (FALSE is returned if not implemented).
260 bool isBonded() const override { return false; }
261 // Both of these methods are called only for contacts with facets where the wall
262 // resolution scheme is set the full. In such cases one might wish to propagate
263 // contact state information (e.g., shear force) from one active contact to another.
264 // See the Faceted Wall section in the documentation.
265 void propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
266 void setNonForcePropsFrom(IContactModel *oldCM) override;
267 /// Return the total force that the contact model holds.
268 DVect getForce() const override;
269 /// Return the total moment on 1 that the contact model holds
270 DAVect getMomentOn1(const IContactMechanical *) const override;
271 /// Return the total moment on 1 that the contact model holds
272 DAVect getMomentOn2(const IContactMechanical *) const override;
273
274 DAVect getModelMomentOn1() const override;
275 DAVect getModelMomentOn2() const override;
276
277 // Used to efficiently get properties from the contact model for the object pane.
278 // List of properties for the object pane, comma separated.
279 // All properties will be cast to doubles for comparison. No other comparisons
280 // are supported. This may not be the same as the entire property list.
281 // Return property name and type for plotting.
282 void objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *) const override;
283 // All properties cast to doubles - this is what can be compared.
284 void objectPropValues(std::vector<double> *,const IContact *) const override;
285
286 // Methods to get and set properties.
287 const double & kn() const {return kn_;}
288 void kn(const double &d) {kn_=d;}
289 const double & ks() const {return ks_;}
290 void ks(const double &d) {ks_=d;}
291 const double & fric() const {return fric_;}
292 void fric(const double &d) {fric_=d;}
293 const DVect & lin_F() const {return lin_F_;}
294 void lin_F(const DVect &f) { lin_F_=f;}
295 bool lin_S() const {return lin_S_;}
296 void lin_S(bool b) { lin_S_=b;}
297 uint32 lin_mode() const {return lin_mode_;}
298 void lin_mode(uint32 i) { lin_mode_= i;}
299 const double & rgap() const {return rgap_;}
300 void rgap(const double &d) {rgap_=d;}
301
302 bool hasDamping() const {return dpProps_ ? true : false;}
303 double dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
304 void dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
305 double dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
306 void dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
307 int dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
308 void dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
309 DVect dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
310 void dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
311
312 bool hasEnergies() const {return energies_ ? true:false;}
313 double estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
314 void estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
315 double errstrain() const {return hasEnergies() ? energies_->errstrain_: 0.0;}
316 void errstrain(const double &d) { if(!hasEnergies()) return; energies_->errstrain_=d;}
317 double eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
318 void eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
319 double errslip() const {return hasEnergies() ? energies_->errslip_: 0.0;}
320 void errslip(const double &d) { if(!hasEnergies()) return; energies_->errslip_=d;}
321 double edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
322 void edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
323
324 uint32 inheritanceField() const {return inheritanceField_;}
325 void inheritanceField(uint32 i) {inheritanceField_ = i;}
326
327 const DVect2 & effectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
328 void effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
329 const DAVect & effectiveRotationalStiffness() const {return effectiveRotationalStiffness_;}
330 void effectiveRotationalStiffness(const DAVect &v ) {effectiveRotationalStiffness_=v;}
331
332 // Rolling resistance methods
333 const double & res_fric() const {return res_fric_;}
334 void res_fric(const double &d) {res_fric_=d;}
335 const DAVect & res_M() const {return res_M_;}
336 void res_M(const DAVect &f) { res_M_=f;}
337 bool res_S() const {return res_S_;}
338 void res_S(bool b) { res_S_=b;}
339 const double & kr() const {return kr_;}
340 void kr(const double &d) {kr_=d;}
341 const double & fr() const {return fr_;}
342 void fr(const double &d) {fr_=d;}
343
344 private:
345 // Index - used internally by PFC. Should be set to -1 in the cpp file.
346 static int index_;
347
348 // Structure to store the energies.
349 struct Energies {
350 Energies() : estrain_(0.0),errstrain_(0.0),eslip_(0.0),errslip_(0.0),edashpot_(0.0) {}
351 double estrain_; // elastic energy stored in linear group
352 double errstrain_; // elastic energy stored in rolling resistance group
353 double eslip_; // work dissipated by friction
354 double errslip_; // work dissipated by rolling resistance friction
355 double edashpot_; // work dissipated by dashpots
356 };
357
358 // Structure to store dashpot quantities.
359 struct dpProps {
360 dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)) {}
361 double dp_nratio_; // normal viscous critical damping ratio
362 double dp_sratio_; // shear viscous critical damping ratio
363 int dp_mode_; // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
364 DVect dp_F_; // Force in the dashpots
365 };
366
367 bool updateKn(const IContactMechanical *con);
368 bool updateKs(const IContactMechanical *con);
369 bool updateFric(const IContactMechanical *con);
370 bool updateResFric(const IContactMechanical *con);
371
372 void updateStiffness(ContactModelMechanicalState *state);
373
374 void setDampCoefficients(const double &mass,double *vcn,double *vcs);
375
376 // Contact model inheritance fields.
377 uint32 inheritanceField_;
378
379 // Effective translational stiffness.
380 DVect2 effectiveTranslationalStiffness_ = DVect2(0.0);
381 DAVect effectiveRotationalStiffness_ = DAVect(0.0); // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)
382
383 // linear model properties
384 double kn_ = 0.0; // Normal stiffness
385 double ks_ = 0.0; // Shear stiffness
386 double fric_ = 0.0; // Coulomb friction coefficient
387 DVect lin_F_ = DVect(0.0); // Force carried in the linear model
388 bool lin_S_ = false; // The current slip state
389 uint32 lin_mode_ = 0; // Specifies absolute (0) or incremental (1) calculation mode
390 double rgap_ = 0.0; // Reference gap
391 dpProps * dpProps_ = nullptr; // The viscous properties
392
393 // rolling resistance properties
394 double res_fric_ = 0.0; // rolling friction coefficient
395 DAVect res_M_ = DAVect(0.0); // moment (bending only)
396 bool res_S_ = false; // The current rolling resistance slip state
397 double kr_ = 0.0; // bending rotational stiffness (read-only, calculated internaly)
398 double fr_ = 0.0; // rolling friction coefficient (rbar*res_fric_) (calculated internaly, not a property)
399
400 double userArea_ = 0.0; // Area as specified by the user
401
402 Energies * energies_ = nullptr; // The energies
403
404 };
405} // namespace cmodelsxd
406// EoF
contactmodelrrlinear.cpp
1// contactmodelrrlinear.cpp
2#include "contactmodelrrlinear.h"
3
4#include "module/interface/icontactmechanical.h"
5#include "module/interface/icontact.h"
6#include "module/interface/ipiece.h"
7#include "module/interface/ifishcalllist.h"
8#include "fish/src/parameter.h"
9
10#include "../version.txt"
11
12#include "utility/src/tptr.h"
13#include "shared/src/mathutil.h"
14
15#include "kernel/interface/iprogram.h"
16#include "module/interface/icontactthermal.h"
17#include "contactmodel/src/contactmodelthermal.h"
18
19#ifdef rrlinear_LIB
20#ifdef _WIN32
21 int __stdcall DllMain(void *,unsigned, void *) {
22 return 1;
23 }
24#endif
25 extern "C" EXPORT_TAG const char *getName() {
26#if DIM==3
27 return "contactmodelmechanical3drrlinear";
28#else
29 return "contactmodelmechanical2drrlinear";
30#endif
31 }
32
33 extern "C" EXPORT_TAG unsigned getMajorVersion() {
34 return MAJOR_VERSION;
35 }
36
37 extern "C" EXPORT_TAG unsigned getMinorVersion() {
38 return MINOR_VERSION;
39 }
40
41 extern "C" EXPORT_TAG void *createInstance() {
42 cmodelsxd::ContactModelRRLinear *m = NEW cmodelsxd::ContactModelRRLinear();
43 return (void *)m;
44 }
45#endif
46
47namespace cmodelsxd {
48 static const uint32 linKnMask = 0x00000002; // Base 1!
49 static const uint32 linKsMask = 0x00000004;
50 static const uint32 linFricMask = 0x00000008;
51 static const uint32 resFricMask = 0x00004000;
52
53 using namespace itasca;
54
55 int ContactModelRRLinear::index_ = -1;
56 uint32 ContactModelRRLinear::getMinorVersion() const { return MINOR_VERSION; }
57
58 ContactModelRRLinear::ContactModelRRLinear() : inheritanceField_(linKnMask|linKsMask|linFricMask|resFricMask) {
59 }
60
61 ContactModelRRLinear::ContactModelRRLinear(const ContactModelRRLinear &m) noexcept {
62 inheritanceField(linKnMask|linKsMask|linFricMask|resFricMask);
63 this->copy(&m);
64 }
65
66 const ContactModelRRLinear & ContactModelRRLinear::operator=(const ContactModelRRLinear &m) {
67 inheritanceField(linKnMask|linKsMask|linFricMask|resFricMask);
68 this->copy(&m);
69 return *this;
70 }
71
72 void ContactModelRRLinear::addToStorage(poly::vector<ContactModelMechanical> *s,int ww) {
73 s->addToStorage<ContactModelRRLinear>(*this,ww);
74 }
75
76 ContactModelRRLinear::~ContactModelRRLinear() {
77 // Make sure to clean up after yourself!
78 if (dpProps_)
79 delete dpProps_;
80 if (energies_)
81 delete energies_;
82 }
83
84 void ContactModelRRLinear::archive(ArchiveStream &stream) {
85 // The stream allows one to archive the values of the contact model
86 // so that it can be saved and restored. The minor version can be
87 // used here to allow for incremental changes to the contact model too.
88 stream & kn_;
89 stream & ks_;
90 stream & fric_;
91 stream & lin_F_;
92 stream & lin_S_;
93 stream & lin_mode_;
94 stream & rgap_;
95 stream & res_fric_;
96 stream & res_M_;
97 stream & res_S_;
98 stream & kr_;
99 stream & fr_;
100
101 if (stream.getArchiveState()==ArchiveStream::Save) {
102 bool b = false;
103 if (dpProps_) {
104 b = true;
105 stream & b;
106 stream & dpProps_->dp_nratio_;
107 stream & dpProps_->dp_sratio_;
108 stream & dpProps_->dp_mode_;
109 stream & dpProps_->dp_F_;
110 }
111 else
112 stream & b;
113
114 b = false;
115 if (energies_) {
116 b = true;
117 stream & b;
118 stream & energies_->estrain_;
119 stream & energies_->errstrain_;
120 stream & energies_->eslip_;
121 stream & energies_->errslip_;
122 stream & energies_->edashpot_;
123 }
124 else
125 stream & b;
126 } else {
127 bool b(false);
128 stream & b;
129 if (b) {
130 if (!dpProps_)
131 dpProps_ = NEW dpProps();
132 stream & dpProps_->dp_nratio_;
133 stream & dpProps_->dp_sratio_;
134 stream & dpProps_->dp_mode_;
135 stream & dpProps_->dp_F_;
136 }
137 stream & b;
138 if (b) {
139 if (!energies_)
140 energies_ = NEW Energies();
141 stream & energies_->estrain_;
142 stream & energies_->errstrain_;
143 stream & energies_->eslip_;
144 stream & energies_->errslip_;
145 stream & energies_->edashpot_;
146 }
147 }
148
149 stream & inheritanceField_;
150 stream & effectiveTranslationalStiffness_;
151 stream & effectiveRotationalStiffness_;
152
153 if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() >= 20)
154 stream & userArea_;
155 }
156
157 void ContactModelRRLinear::copy(const ContactModel *cm) {
158 // Copy all of the contact model properties. Used in the CMAT
159 // when a new contact is created.
160 ContactModelMechanical::copy(cm);
161 const ContactModelRRLinear *in = dynamic_cast<const ContactModelRRLinear*>(cm);
162 if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
163 kn(in->kn());
164 ks(in->ks());
165 fric(in->fric());
166 lin_F(in->lin_F());
167 lin_S(in->lin_S());
168 lin_mode(in->lin_mode());
169 rgap(in->rgap());
170 res_fric(in->res_fric());
171 res_M(in->res_M());
172 res_S(in->res_S());
173 kr(in->kr());
174 fr(in->fr());
175
176 if (in->hasDamping()) {
177 if (!dpProps_)
178 dpProps_ = NEW dpProps();
179 dp_nratio(in->dp_nratio());
180 dp_sratio(in->dp_sratio());
181 dp_mode(in->dp_mode());
182 dp_F(in->dp_F());
183 }
184 if (in->hasEnergies()) {
185 if (!energies_)
186 energies_ = NEW Energies();
187 estrain(in->estrain());
188 errstrain(in->errstrain());
189 eslip(in->eslip());
190 errslip(in->errslip());
191 edashpot(in->edashpot());
192 }
193 userArea_ = in->userArea_;
194 inheritanceField(in->inheritanceField());
195 effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
196 effectiveRotationalStiffness(in->effectiveRotationalStiffness());
197 }
198
199
200 base::Property ContactModelRRLinear::getProperty(uint32 i,const IContact *con) const {
201 // Return the property. The IContact pointer is provided so that
202 // more complicated properties, depending on contact characteristics,
203 // can be calcualted.
204 // The IContact pointer may be a nullptr!
205 base::Property var;
206 switch (i) {
207 case kwKn: return kn_;
208 case kwKs: return ks_;
209 case kwFric: return fric_;
210 case kwLinF: var.setValue(lin_F_); return var;
211 case kwLinS: return lin_S_;
212 case kwLinMode: return lin_mode_;
213 case kwRGap: return rgap_;
214 case kwEmod: {
215 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
216 if (c ==nullptr) return 0.0;
217 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
218 double rsum = calcRSum(c);
219 if (userArea_) {
220#ifdef THREED
221 rsq = std::sqrt(userArea_ / dPi);
222#else
223 rsq = userArea_ / 2.0;
224#endif
225 rsum = rsq + rsq;
226 rsq = 1. / rsq;
227 }
228#ifdef TWOD
229 return (kn_ * rsum * rsq / 2.0);
230#else
231 return (kn_ * rsum * rsq * rsq) / dPi;
232#endif
233 }
234 case kwKRatio: return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
235 case kwDpNRatio: return dpProps_ ? dpProps_->dp_nratio_ : 0;
236 case kwDpSRatio: return dpProps_ ? dpProps_->dp_sratio_ : 0;
237 case kwDpMode: return dpProps_ ? dpProps_->dp_mode_ : 0;
238 case kwDpF: {
239 dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
240 return var;
241 }
242 case kwResFric: return res_fric_;
243 case kwResMoment: var.setValue(res_M_); return var;
244 case kwResS: return res_S_;
245 case kwResKr: return kr_;
246 case kwUserArea: return userArea_;
247 }
248 assert(0);
249 return base::Property();
250 }
251
252 bool ContactModelRRLinear::getPropertyGlobal(uint32 i) const {
253 // Returns whether or not a property is held in the global axis system (TRUE)
254 // or the local system (FALSE). Used by the plotting logic.
255 switch (i) {
256 case kwLinF:
257 case kwDpF:
258 case kwResMoment:
259 return false;
260 }
261 return true;
262 }
263
264 bool ContactModelRRLinear::setProperty(uint32 i,const base::Property &v,IContact *) {
265 // Set a contact model property. Return value indicates that the timestep
266 // should be recalculated.
267 dpProps dp;
268 switch (i) {
269 case kwKn: {
270 if (!v.canConvert<double>())
271 throw Exception("kn must be a double.");
272 double val(v.toDouble());
273 if (val<0.0)
274 throw Exception("Negative kn not allowed.");
275 kn_ = val;
276 return true;
277 }
278 case kwKs: {
279 if (!v.canConvert<double>())
280 throw Exception("ks must be a double.");
281 double val(v.toDouble());
282 if (val<0.0)
283 throw Exception("Negative ks not allowed.");
284 ks_ = val;
285 return true;
286 }
287 case kwFric: {
288 if (!v.canConvert<double>())
289 throw Exception("fric must be a double.");
290 double val(v.toDouble());
291 if (val<0.0)
292 throw Exception("Negative fric not allowed.");
293 fric_ = val;
294 return false;
295 }
296 case kwLinF: {
297 if (!v.canConvert<DVect>())
298 throw Exception("lin_force must be a vector.");
299 DVect val(v.value<DVect>());
300 lin_F_ = val;
301 return false;
302 }
303 case kwLinMode: {
304 if (!v.canConvert<int64>())
305 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
306 uint32 val(v.toUInt());
307 if (val >1)
308 throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
309 lin_mode_ = val;
310 return false;
311 }
312 case kwRGap: {
313 if (!v.canConvert<double>())
314 throw Exception("Reference gap must be a double.");
315 double val(v.toDouble());
316 rgap_ = val;
317 return false;
318 }
319 case kwDpNRatio: {
320 if (!v.canConvert<double>())
321 throw Exception("dp_nratio must be a double.");
322 double val(v.toDouble());
323 if (val<0.0)
324 throw Exception("Negative dp_nratio not allowed.");
325 if (val == 0.0 && !dpProps_)
326 return false;
327 if (!dpProps_)
328 dpProps_ = NEW dpProps();
329 dpProps_->dp_nratio_ = val;
330 return true;
331 }
332 case kwDpSRatio: {
333 if (!v.canConvert<double>())
334 throw Exception("dp_sratio must be a double.");
335 double val(v.toDouble());
336 if (val<0.0)
337 throw Exception("Negative dp_sratio not allowed.");
338 if (val == 0.0 && !dpProps_)
339 return false;
340 if (!dpProps_)
341 dpProps_ = NEW dpProps();
342 dpProps_->dp_sratio_ = val;
343 return true;
344 }
345 case kwDpMode: {
346 if (!v.canConvert<int64>())
347 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
348 int val(v.toInt());
349 if (val == 0 && !dpProps_)
350 return false;
351 if (val < 0 || val > 3)
352 throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
353 if (!dpProps_)
354 dpProps_ = NEW dpProps();
355 dpProps_->dp_mode_ = val;
356 return false;
357 }
358 case kwDpF: {
359 if (!v.canConvert<DVect>())
360 throw Exception("dp_force must be a vector.");
361 DVect val(v.value<DVect>());
362 if (val.fsum() == 0.0 && !dpProps_)
363 return false;
364 if (!dpProps_)
365 dpProps_ = NEW dpProps();
366 dpProps_->dp_F_ = val;
367 return false;
368 }
369 case kwResFric: {
370 if (!v.canConvert<double>())
371 throw Exception("res_fric must be a double.");
372 double val(v.toDouble());
373 if (val<0.0)
374 throw Exception("Negative res_fric not allowed.");
375 res_fric_ = val;
376 return false;
377 }
378 case kwResMoment: {
379 DAVect val(0.0);
380#ifdef TWOD
381 if (!v.canConvert<DAVect>() && !v.canConvert<double>())
382 throw Exception("res_moment must be an angular vector.");
383 if (v.canConvert<DAVect>())
384 val = DAVect(v.value<DAVect>());
385 else
386 val = DAVect(v.toDouble());
387#else
388 if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
389 throw Exception("res_moment must be an angular vector.");
390 if (v.canConvert<DAVect>())
391 val = DAVect(v.value<DAVect>());
392 else
393 val = DAVect(v.value<DVect>());
394#endif
395 res_M_ = val;
396 return false;
397 }
398 case kwUserArea: {
399 if (!v.canConvert<double>())
400 throw Exception("user_area must be a double.");
401 double val(v.toDouble());
402 if (val < 0.0)
403 throw Exception("Negative user_area not allowed.");
404 userArea_ = val;
405 return true;
406 }
407 }
408 return false;
409 }
410
411 bool ContactModelRRLinear::getPropertyReadOnly(uint32 i) const {
412 // Returns TRUE if a property is read only or FALSE otherwise.
413 switch (i) {
414 case kwDpF:
415 case kwLinS:
416 case kwEmod:
417 case kwKRatio:
418 case kwResS:
419 case kwResKr:
420 return true;
421 default:
422 break;
423 }
424 return false;
425 }
426
427 bool ContactModelRRLinear::supportsInheritance(uint32 i) const {
428 // Returns TRUE if a property supports inheritance or FALSE otherwise.
429 switch (i) {
430 case kwKn:
431 case kwKs:
432 case kwFric:
433 case kwResFric:
434 return true;
435 default:
436 break;
437 }
438 return false;
439 }
440
441 string ContactModelRRLinear::getMethodArguments(uint32 i) const {
442 // Return a list of contact model method argument names.
443 switch (i) {
444 case kwDeformability:
445 return "emod,kratio";
446 case kwArea:
447 return string();
448 }
449 assert(0);
450 return string();
451 }
452
453 bool ContactModelRRLinear::setMethod(uint32 i,const std::vector<base::Property> &vl,IContact *con) {
454 FP_S;
455 // Apply the specified method.
456 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
457 switch (i) {
458 case kwDeformability: {
459 double emod;
460 double krat;
461 if (vl.at(0).isNull())
462 throw Exception("Argument emod must be specified with method deformability in contact model {0}.",getName());
463 emod = vl.at(0).toDouble();
464 if (emod<0.0)
465 throw Exception("Negative emod not allowed in contact model {0}.",getName());
466 if (vl.at(1).isNull())
467 throw Exception("Argument kratio must be specified with method deformability in contact model {0}.",getName());
468 krat = vl.at(1).toDouble();
469 if (krat<0.0)
470 throw Exception("Negative stiffness ratio not allowed in contact model {0}.",getName());
471 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
472 double rsum = calcRSum(c);
473 if (userArea_) {
474#ifdef THREED
475 rsq = std::sqrt(userArea_ / dPi);
476#else
477 rsq = userArea_ / 2.0;
478#endif
479 rsum = rsq + rsq;
480 rsq = safeDiv(1. , rsq);
481 }
482#ifdef TWOD
483 kn_ = safeDiv(2.0 * emod ,rsq * rsum);
484#else
485 kn_ = safeDiv(dPi * emod ,rsq * rsq * rsum);
486#endif
487 ks_ = safeDiv(kn_ , krat);
488 setInheritance(1,false);
489 setInheritance(2,false);
490 FP_S;
491 return true;
492 }
493 case kwArea: {
494 if (!userArea_) {
495 double rsq = calcRSQ(c);
496#ifdef THREED
497 userArea_ = rsq * rsq * dPi;
498#else
499 userArea_ = rsq * 2.0;
500#endif
501 }
502 FP_S;
503 return true;
504 }
505
506 }
507 FP_S;
508 return false;
509 }
510
511 double ContactModelRRLinear::getEnergy(uint32 i) const {
512 // Return an energy value.
513 double ret(0.0);
514 if (!energies_)
515 return ret;
516 switch (i) {
517 case kwEStrain: return energies_->estrain_;
518 case kwERRStrain: return energies_->errstrain_;
519 case kwESlip: return energies_->eslip_;
520 case kwERRSlip: return energies_->errslip_;
521 case kwEDashpot: return energies_->edashpot_;
522 }
523 assert(0);
524 return ret;
525 }
526
527 bool ContactModelRRLinear::getEnergyAccumulate(uint32 i) const {
528 // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
529 switch (i) {
530 case kwEStrain: return false;
531 case kwERRStrain: return false;
532 case kwESlip: return true;
533 case kwERRSlip: return true;
534 case kwEDashpot: return true;
535 }
536 assert(0);
537 return false;
538 }
539
540 void ContactModelRRLinear::setEnergy(uint32 i,const double &d) {
541 // Set an energy value.
542 if (!energies_) return;
543 switch (i) {
544 case kwEStrain: energies_->estrain_ = d; return;
545 case kwERRStrain: energies_->errstrain_ = d; return;
546 case kwESlip: energies_->eslip_ = d; return;
547 case kwERRSlip: energies_->errslip_ = d; return;
548 case kwEDashpot: energies_->edashpot_= d; return;
549 }
550 assert(0);
551 return;
552 }
553
554 bool ContactModelRRLinear::validate(ContactModelMechanicalState *state,const double &) {
555 // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
556 // (1) the contact is created, (2) a property of the contact that returns a true via
557 // the setProperty method has been modified and (3) when a set of cycles is executed
558 // via the {cycle N} command.
559 // Return value indicates contact activity (TRUE: active, FALSE: inactive).
560 FP_S;
561 assert(state);
562 const IContactMechanical *c = state->getMechanicalContact();
563 assert(c);
564
565 FP_S;
566 if (state->trackEnergy_)
567 activateEnergy();
568
569 FP_S;
570 if (inheritanceField_ & linKnMask)
571 updateKn(c);
572 FP_S;
573 if (inheritanceField_ & linKsMask)
574 updateKs(c);
575 FP_S;
576 if (inheritanceField_ & linFricMask)
577 updateFric(c);
578 FP_S;
579 if (inheritanceField_ & resFricMask)
580 updateResFric(c);
581 FP_S;
582
583 updateStiffness(state);
584 FP_S;
585 auto ret = checkActivity(state->gap_);
586 FP_S;
587 return ret;
588 }
589
590 static const string knstr("kn");
591 bool ContactModelRRLinear::updateKn(const IContactMechanical *con) {
592 assert(con);
593 base::Property v1 = con->getEnd1()->getProperty(knstr);
594 base::Property v2 = con->getEnd2()->getProperty(knstr);
595 if (!v1.isValid() || !v2.isValid())
596 return false;
597 double kn1 = v1.toDouble();
598 double kn2 = v2.toDouble();
599 double val = kn_;
600 if (kn1 && kn2)
601 kn_ = kn1*kn2/(kn1+kn2);
602 else if (kn1)
603 kn_ = kn1;
604 else if (kn2)
605 kn_ = kn2;
606 return ( (kn_ != val) );
607 }
608
609 static const string ksstr("ks");
610 bool ContactModelRRLinear::updateKs(const IContactMechanical *con) {
611 assert(con);
612 base::Property v1 = con->getEnd1()->getProperty(ksstr);
613 base::Property v2 = con->getEnd2()->getProperty(ksstr);
614 if (!v1.isValid() || !v2.isValid())
615 return false;
616 double ks1 = v1.toDouble();
617 double ks2 = v2.toDouble();
618 double val = ks_;
619 if (ks1 && ks2)
620 ks_ = ks1*ks2/(ks1+ks2);
621 else if (ks1)
622 ks_ = ks1;
623 else if (ks2)
624 ks_ = ks2;
625 return ( (ks_ != val) );
626 }
627
628 static const string fricstr("fric");
629 bool ContactModelRRLinear::updateFric(const IContactMechanical *con) {
630 assert(con);
631 base::Property v1 = con->getEnd1()->getProperty(fricstr);
632 base::Property v2 = con->getEnd2()->getProperty(fricstr);
633 if (!v1.isValid() || !v2.isValid())
634 return false;
635 double fric1 = std::max(0.0,v1.toDouble());
636 double fric2 = std::max(0.0,v2.toDouble());
637 double val = fric_;
638 fric_ = std::min(fric1,fric2);
639 return ( (fric_ != val) );
640 }
641
642 static const string rfricstr("rr_fric");
643 bool ContactModelRRLinear::updateResFric(const IContactMechanical *con) {
644 assert(con);
645 base::Property v1 = con->getEnd1()->getProperty(rfricstr);
646 base::Property v2 = con->getEnd2()->getProperty(rfricstr);
647 if (!v1.isValid() || !v2.isValid())
648 return false;
649 double rfric1 = std::max(0.0,v1.toDouble());
650 double rfric2 = std::max(0.0,v2.toDouble());
651 double val = res_fric_;
652 res_fric_ = std::min(rfric1,rfric2);
653 return ( (res_fric_ != val) );
654 }
655
656 bool ContactModelRRLinear::endPropertyUpdated(const string &name,const IContactMechanical *c) {
657 // The endPropertyUpdated method is called whenever a surface property (with a name
658 // that matches an inheritable contact model property name) of one of the contacting
659 // pieces is modified. This allows the contact model to update its associated
660 // properties. The return value denotes whether or not the update has affected
661 // the time step computation (by having modified the translational or rotational
662 // tangent stiffnesses). If true is returned, then the time step will be recomputed.
663 assert(c);
664 StringList availableProperties = split(replace(simplified(getProperties())," ",""),",");
665 auto idx = findRegex(availableProperties,name);
666 if (idx==string::npos) return false;
667 idx += 1;
668 bool ret=false;
669 switch(idx) {
670 case kwKn: { //kn
671 if (inheritanceField_ & linKnMask)
672 ret = updateKn(c);
673 break;
674 }
675 case kwKs: { //ks
676 if (inheritanceField_ & linKsMask)
677 ret =updateKs(c);
678 break;
679 }
680 case kwFric: { //fric
681 if (inheritanceField_ & linFricMask)
682 updateFric(c);
683 break;
684 }
685 case kwResFric: { //rr_fric
686 if (inheritanceField_ & resFricMask)
687 ret = updateResFric(c);
688 break;
689 }
690 }
691 return ret;
692 }
693
694 void ContactModelRRLinear::updateStiffness(ContactModelMechanicalState *state) {
695 FP_S;
696 // first compute rolling resistance stiffness
697 kr_ = 0.0;
698 if (res_fric_ > 0.0) {
699 FP_S;
700 double rbar = 0.0;
701 double r1 = safeDiv(1.0,state->end1Curvature_.y());
702 rbar = r1;
703 double r2 = 0.0;
704 if (state->end2Curvature_.y()) {
705 r2 = safeDiv(1.0 , state->end2Curvature_.y());
706 rbar = (r1*r2) / (r1+r2);
707 }
708 if (userArea_) {
709#ifdef THREED
710 r1 = std::sqrt(userArea_ / dPi);
711#else
712 r1 = userArea_ / 2.0;
713#endif
714 rbar = r1/2.0;
715 }
716 kr_ = ks_*rbar*rbar;
717 fr_ = res_fric_*rbar;
718 }
719 // Now calculate effective stiffness
720 DVect2 retT(kn_,ks_);
721 // correction if viscous damping active
722 if (dpProps_) {
723 DVect2 correct(1.0);
724 if (dpProps_->dp_nratio_)
725 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
726 if (dpProps_->dp_sratio_)
727 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
728 retT /= (correct*correct);
729 }
730 effectiveTranslationalStiffness_ = retT;
731 // Effective rotational stiffness (bending only)
732 effectiveRotationalStiffness_ = DAVect(kr_);
733#if DIM==3
734 effectiveRotationalStiffness_.rx() = 0.0;
735#endif
736 }
737
738 bool ContactModelRRLinear::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
739 assert(state);
740
741 // Current overlap
742 double overlap = rgap_ - state->gap_;
743 // Relative translational increment
744 DVect trans = state->relativeTranslationalIncrement_;
745 // Correction factor to account for when the contact becomes newly active.
746 // We estimate the time of activity during the timestep when the contact has first
747 // become active and scale the forces accordingly.
748 double correction = 1.0;
749
750 // The contact was just activated from an inactive state
751 if (state->activated()) {
752 // Trigger the FISH callback if one is hooked up to the
753 // contact_activated event.
754 if (cmEvents_[fActivated] >= 0) {
755 // An FArray of base::Property is returned and these will be passed
756 // to the FISH function as an array of FISH symbols as the second
757 // argument to the FISH callback function.
758 auto c = state->getContact();
759 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
760 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
761 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
762 }
763 // Calculate the correction factor.
764 if (lin_mode_ == 0 && trans.x()) {
765 correction = -1.0*overlap / trans.x();
766 if (correction < 0)
767 correction = 1.0;
768 }
769 }
770
771 // Angular dispacement increment.
772 DAVect ang = state->relativeAngularIncrement_;
773 DVect lin_F_old = lin_F_;
774
775 if (lin_mode_ == 0)
776 lin_F_.rx() = overlap * kn_; // incremental mode for normal force calculation
777 else
778 lin_F_.rx() -= correction * trans.x() * kn_; // absolute mode for normal force calculation
779
780 // Normal force can only be positive.
781 lin_F_.rx() = std::max(0.0,lin_F_.x());
782
783 // Calculate the shear force.
784 DVect sforce(0.0);
785 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
786 // Loop over the shear components (note: the 0 component is the normal component)
787 // and calculate the shear force.
788 for (int i=1; i<dim; ++i)
789 sforce.rdof(i) = lin_F_.dof(i) - trans.dof(i) * ks_ * correction;
790
791 // The canFail flag corresponds to whether or not the contact can undergo non-linear
792 // force-displacement response. If the SOLVE ELASTIC command is given then the
793 // canFail state is set to FALSE. Otherwise it is always TRUE.
794 if (state->canFail_) {
795 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
796 double crit = lin_F_.x() * fric_;
797 // The is the magnitude of the shear force.
798 double sfmag = sforce.mag();
799 // Sliding occurs when the magnitude of the shear force is greater than the
800 // critical value.
801 if (sfmag > crit) {
802 // Lower the shear force to the critical value for sliding.
803 double rat = crit / sfmag;
804 sforce *= rat;
805 // Handle the slip_change event if one has been hooked up. Sliding has commenced.
806 if (!lin_S_ && cmEvents_[fSlipChange] >= 0) {
807 auto c = state->getContact();
808 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
809 fish::Parameter() };
810 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
811 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
812 }
813 lin_S_ = true;
814 } else {
815 // Handle the slip_change event if one has been hooked up and
816 // the contact was previously sliding. Sliding has ceased.
817 if (lin_S_) {
818 if (cmEvents_[fSlipChange] >= 0) {
819 auto c = state->getContact();
820 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
821 fish::Parameter((int64)1) };
822 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
823 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
824 }
825 lin_S_ = false;
826 }
827 }
828 }
829
830 // Set the shear components of the total force.
831 for (int i=1; i<dim; ++i)
832 lin_F_.rdof(i) = sforce.dof(i);
833
834 // Rolling resistance
835 DAVect res_M_old = res_M_;
836 if ((fr_ == 0.0) || (kr_==0.0)) {
837 res_M_.fill(0.0);
838 } else {
839 DAVect angStiff(0.0);
840 DAVect MomentInc(0.0);
841#if DIM==3
842 angStiff.rx() = 0.0;
843 angStiff.ry() = kr_;
844#endif
845 angStiff.rz() = kr_;
846 MomentInc = ang * angStiff * (-1.0);
847 res_M_ += MomentInc;
848 if (state->canFail_) {
849 // Account for bending strength
850 double rmax = std::abs(fr_*lin_F_.x());
851 double rmag = res_M_.mag();
852 if (rmag > rmax) {
853 double fac = rmax/rmag;
854 res_M_ *= fac;
855 res_S_ = true;
856 } else {
857 res_S_ = false;
858 }
859 }
860 }
861
862 // Account for dashpot forces if the dashpot structure has been defined.
863 if (dpProps_) {
864 dpProps_->dp_F_.fill(0.0);
865 double vcn(0.0), vcs(0.0);
866 // Calculate the damping coefficients.
867 setDampCoefficients(state->inertialMass_,&vcn,&vcs);
868 // First damp the shear components
869 for (int i=1; i<dim; ++i)
870 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
871 // Damp the normal component
872 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
873 // Need to change behavior based on the dp_mode.
874 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
875 // Limit in tension if not bonded.
876 if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
877 dpProps_->dp_F_.rx() = - lin_F_.rx();
878 }
879 if (lin_S_ && dpProps_->dp_mode_ > 1) {
880 // Limit in shear if not sliding.
881 double dfn = dpProps_->dp_F_.rx();
882 dpProps_->dp_F_.fill(0.0);
883 dpProps_->dp_F_.rx() = dfn;
884 }
885 }
886
887 //Compute energies if energy tracking has been enabled.
888 if (state->trackEnergy_) {
889 assert(energies_);
890 energies_->estrain_ = 0.0;
891 if (kn_)
892 // Calculate the strain energy.
893 energies_->estrain_ = 0.5*lin_F_.x()*lin_F_.x()/kn_;
894 if (ks_) {
895 DVect s = lin_F_;
896 s.rx() = 0.0;
897 double smag2 = s.mag2();
898 // Add the shear component of the strain energy.
899 energies_->estrain_ += 0.5*smag2 / ks_;
900
901 if (lin_S_) {
902 // If sliding calculate the slip energy and accumulate it.
903 lin_F_old.rx() = 0.0;
904 DVect avg_F_s = (s + lin_F_old)*0.5;
905 DVect u_s_el = (s - lin_F_old) / ks_;
906 DVect u_s(0.0);
907 for (int i=1; i<dim; ++i)
908 u_s.rdof(i) = trans.dof(i);
909 energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
910 }
911 }
912 // Add the rolling resistance energy contributions.
913 energies_->errstrain_ = 0.0;
914 if (kr_) {
915 energies_->errstrain_ = 0.5*res_M_.mag2() / kr_;
916 if (res_S_) {
917 // If sliding calculate the slip energy and accumulate it.
918 DAVect avg_M = (res_M_ + res_M_old)*0.5;
919 DAVect t_s_el = (res_M_ - res_M_old) / kr_;
920 energies_->errslip_ -= std::min(0.0,(avg_M | (ang + t_s_el)));
921 }
922
923 }
924 if (dpProps_) {
925 // Calculate damping energy (accumulated) if the dashpots are active.
926 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
927 }
928 }
929
930 // This is just a sanity check to ensure, in debug mode, that the force isn't wonky.
931 assert(lin_F_ == lin_F_);
932 return true;
933 }
934
935 bool ContactModelRRLinear::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
936 // Account for thermal expansion in incremental mode
937 if (lin_mode_ == 0 || ts->gapInc_ == 0.0) return false;
938 DVect finc(0.0);
939 finc.rx() = kn_ * ts->gapInc_;
940 lin_F_ -= finc;
941 return true;
942 }
943
944 void ContactModelRRLinear::setForce(const DVect &v,IContact *c) {
945 lin_F(v);
946 if (v.x() > 0)
947 rgap_ = c->getGap() + v.x() / kn_;
948 }
949
950 void ContactModelRRLinear::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
951 // Only called for contacts with wall facets when the wall resolution scheme
952 // is set to full!
953 // Only do something if the contact model is of the same type
954 if (equal(old->getContactModel()->getName(),"rrlinear") && !isBonded()) {
955 ContactModelRRLinear *oldCm = (ContactModelRRLinear *)old;
956#ifdef THREED
957 // Need to rotate just the shear component from oldSystem to newSystem
958
959 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
960 DVect axis = oldSystem.e1() & newSystem.e1();
961 double c, ang, s;
962 DVect re2;
963 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
964 axis = axis.unit();
965 c = oldSystem.e1()|newSystem.e1();
966 if (c > 0)
967 c = std::min(c,1.0);
968 else
969 c = std::max(c,-1.0);
970 ang = acos(c);
971 s = sin(ang);
972 double t = 1. - c;
973 DMatrix<3,3> rm;
974 rm.get(0,0) = t*axis.x()*axis.x() + c;
975 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
976 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
977 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
978 rm.get(1,1) = t*axis.y()*axis.y() + c;
979 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
980 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
981 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
982 rm.get(2,2) = t*axis.z()*axis.z() + c;
983 re2 = rm*oldSystem.e2();
984 }
985 else
986 re2 = oldSystem.e2();
987 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
988 axis = re2 & newSystem.e2();
989 DVect2 tpf;
990 DVect2 tpm;
991 DMatrix<2,2> m;
992 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
993 axis = axis.unit();
994 c = re2|newSystem.e2();
995 if (c > 0)
996 c = std::min(c,1.0);
997 else
998 c = std::max(c,-1.0);
999 ang = acos(c);
1000 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
1001 ang *= -1;
1002 s = sin(ang);
1003 m.get(0,0) = c;
1004 m.get(1,0) = s;
1005 m.get(0,1) = -m.get(1,0);
1006 m.get(1,1) = m.get(0,0);
1007 tpf = m*DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
1008 tpm = m*DVect2(oldCm->res_M_.y(),oldCm->res_M_.z());
1009 } else {
1010 m.get(0,0) = 1.;
1011 m.get(0,1) = 0.;
1012 m.get(1,0) = 0.;
1013 m.get(1,1) = 1.;
1014 tpf = DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
1015 tpm = DVect2(oldCm->res_M_.y(),oldCm->res_M_.z());
1016 }
1017 DVect pforce = DVect(0,tpf.x(),tpf.y());
1018 //DVect pm = DVect(0,tpm.x(),tpm.y());
1019#else
1020 oldSystem;
1021 newSystem;
1022 DVect pforce = DVect(0,oldCm->lin_F_.y());
1023 //DVect pm = DVect(0,oldCm->res_M_.y());
1024#endif
1025 for (int i=1; i<dim; ++i)
1026 lin_F_.rdof(i) += pforce.dof(i);
1027 if (lin_mode_ && oldCm->lin_mode_)
1028 lin_F_.rx() = oldCm->lin_F_.x();
1029 oldCm->lin_F_ = DVect(0.0);
1030 oldCm->res_M_ = DAVect(0.0);
1031 if (dpProps_ && oldCm->dpProps_) {
1032#ifdef THREED
1033 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
1034 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
1035#else
1036 pforce = oldCm->dpProps_->dp_F_;
1037#endif
1038 dpProps_->dp_F_ += pforce;
1039 oldCm->dpProps_->dp_F_ = DVect(0.0);
1040 }
1041 if(oldCm->getEnergyActivated()) {
1042 activateEnergy();
1043 energies_->estrain_ = oldCm->energies_->estrain_;
1044 energies_->edashpot_ = oldCm->energies_->edashpot_;
1045 energies_->eslip_ = oldCm->energies_->eslip_;
1046 oldCm->energies_->estrain_ = 0.0;
1047 oldCm->energies_->edashpot_ = 0.0;
1048 oldCm->energies_->eslip_ = 0.0;
1049 }
1050 }
1051 assert(lin_F_ == lin_F_);
1052 }
1053
1054 void ContactModelRRLinear::setNonForcePropsFrom(IContactModel *old) {
1055 // Only called for contacts with wall facets when the wall resolution scheme
1056 // is set to full!
1057 // Only do something if the contact model is of the same type
1058 if (equal(old->getName(),"rrlinear") && !isBonded()) {
1059 ContactModelRRLinear *oldCm = (ContactModelRRLinear *)old;
1060 kn_ = oldCm->kn_;
1061 ks_ = oldCm->ks_;
1062 fric_ = oldCm->fric_;
1063 lin_mode_ = oldCm->lin_mode_;
1064 rgap_ = oldCm->rgap_;
1065 res_fric_ = oldCm->res_fric_;
1066 res_S_ = oldCm->res_S_;
1067 kr_ = oldCm->kr_;
1068 fr_ = oldCm->fr_;
1069 userArea_ = oldCm->userArea_;
1070
1071 if (oldCm->dpProps_) {
1072 if (!dpProps_)
1073 dpProps_ = NEW dpProps();
1074 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1075 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1076 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1077 }
1078 }
1079 }
1080
1081 DVect ContactModelRRLinear::getForce() const {
1082 DVect ret(lin_F_);
1083 if (dpProps_)
1084 ret += dpProps_->dp_F_;
1085 return ret;
1086 }
1087
1088 DAVect ContactModelRRLinear::getMomentOn1(const IContactMechanical *c) const {
1089 DAVect ret(res_M_);
1090 if (c) {
1091 DVect force = getForce();
1092 c->updateResultingTorqueOn1Local(force,&ret);
1093 }
1094 return ret;
1095 }
1096
1097 DAVect ContactModelRRLinear::getMomentOn2(const IContactMechanical *c) const {
1098 DAVect ret(res_M_);
1099 if (c) {
1100 DVect force = getForce();
1101 c->updateResultingTorqueOn2Local(force,&ret);
1102 }
1103 return ret;
1104 }
1105
1106 DAVect ContactModelRRLinear::getModelMomentOn1() const {
1107 DAVect ret(res_M_);
1108 return ret;
1109 }
1110
1111 DAVect ContactModelRRLinear::getModelMomentOn2() const {
1112 DAVect ret(res_M_);
1113 return ret;
1114 }
1115
1116 void ContactModelRRLinear::objectPropsTypes(std::vector<std::pair<string,InfoTypes>> *ret) const {
1117 ret->clear();
1118 ret->push_back({"kn",ScalarInfo});
1119 ret->push_back({"ks",ScalarInfo});
1120 ret->push_back({"fric",ScalarInfo});
1121 ret->push_back({"lin_force",VectorInfo});
1122 ret->push_back({"lin_slip",ScalarInfo});
1123 ret->push_back({"lin_mode",ScalarInfo});
1124 ret->push_back({"rgap",ScalarInfo});
1125 ret->push_back({"emod",ScalarInfo});
1126 ret->push_back({"kratio",ScalarInfo});
1127 ret->push_back({"dp_nratio",ScalarInfo});
1128 ret->push_back({"dp_sratio",ScalarInfo});
1129 ret->push_back({"dp_mode",ScalarInfo});
1130 ret->push_back({"dp_force",ScalarInfo});
1131 ret->push_back({"rr_fric",ScalarInfo});
1132 ret->push_back({"rr_moment",VectorInfo});
1133 ret->push_back({"rr_slip",ScalarInfo});
1134 ret->push_back({"rr_kr",ScalarInfo});
1135 ret->push_back({"user_area",ScalarInfo});
1136 }
1137
1138 void ContactModelRRLinear::objectPropValues(std::vector<double> *ret,const IContact *con) const {
1139 FP_S;
1140 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
1141 ret->clear();
1142 ret->push_back(kn_);
1143 ret->push_back(ks_);
1144 ret->push_back(fric_);
1145 ret->push_back(lin_F_.mag());
1146 ret->push_back(lin_S_);
1147 ret->push_back(lin_mode_);
1148 ret->push_back(rgap_);
1149 double d;
1150 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
1151 double rsum = calcRSum(c);
1152 if (userArea_) {
1153#ifdef THREED
1154 rsq = std::sqrt(userArea_ / dPi);
1155#else
1156 rsq = userArea_ / 2.0;
1157#endif
1158 rsum = rsq + rsq;
1159 rsq = safeDiv(1. , rsq);
1160 }
1161#ifdef TWOD
1162 d = (kn_ * rsum * rsq / 2.0);
1163#else
1164 d = (kn_ * rsum * rsq * rsq) / dPi;
1165#endif
1166 ret->push_back(d);
1167 ret->push_back(safeDiv(kn_,ks_));
1168 ret->push_back(dpProps_ ? dpProps_->dp_nratio_ : 0);
1169 ret->push_back(dpProps_ ? dpProps_->dp_sratio_ : 0);
1170 ret->push_back(dpProps_ ? dpProps_->dp_mode_ : 0);
1171 ret->push_back(dpProps_ ? safeMag(dpProps_->dp_F_) : 0);
1172 ret->push_back(res_fric_);
1173 ret->push_back(safeMag(res_M_));
1174 ret->push_back(res_S_);
1175 ret->push_back(kr_);
1176 ret->push_back(userArea_);
1177 FP_S;
1178 }
1179
1180 void ContactModelRRLinear::setDampCoefficients(const double &mass,double *vcn,double *vcs) {
1181 *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(kn_));
1182 *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(ks_));
1183 }
1184
1185} // namespace cmodelsxd
1186// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Jun 07, 2025 |
