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