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 QString 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 QString 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 QString 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_ = NEWC(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 QString getFishCallEvents() const override {
112 return
113 "contact_activated"
114 ",slip_change";
115 }
116
117 // Return the specified contact model property.
118 QVariant 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 QVariant &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 QString getMethods() const override { return "deformability,area";}
151 // Return a comma seprated list of the contact model method arguments (base 1).
152 QString 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 QVector<QVariant> &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 QString &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 NEWC(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<QString,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 = NEWC(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_ = NEWC(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_ = NEWC(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_ = NEWC(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_ = NEWC(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 QVariant 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 QVariant 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 QVariant();
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 QVariant &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<uint32>())
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_ = NEWC(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_ = NEWC(dpProps());
296 dpProps_->dp_sratio_ = val;
297 return true;
298 }
299 case kwDpMode: {
300 if (!v.canConvert<int>())
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_ = NEWC(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_ = NEWC(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 QString 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 QString();
370 }
371 assert(0);
372 return QString();
373 }
374
375 bool ContactModelLinear::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
376 // Apply the specified method.
377 IContactMechanical *c(convert_getcast<IContactMechanical>(con));
378 switch (i) {
379 case kwDeformability: {
380 double emod;
381 double krat;
382 if (vl.at(0).isNull())
383 throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
384 emod = vl.at(0).toDouble();
385 if (emod<0.0)
386 throw Exception("Negative emod not allowed in contact model %1.",getName());
387 if (vl.at(1).isNull())
388 throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
389 krat = vl.at(1).toDouble();
390 if (krat<0.0)
391 throw Exception("Negative stiffness ratio not allowed in contact model %1.",getName());
392 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
393 double rsum(0.0);
394 if (c->getEnd1Curvature().y())
395 rsum += 1.0/c->getEnd1Curvature().y();
396 if (c->getEnd2Curvature().y())
397 rsum += 1.0/c->getEnd2Curvature().y();
398 if (userArea_) {
399#ifdef THREED
400 rsq = std::sqrt(userArea_ / dPi);
401#else
402 rsq = userArea_ / 2.0;
403#endif
404 rsum = rsq + rsq;
405 rsq = 1. / rsq;
406 }
407#ifdef TWOD
408 kn_ = 2.0 * emod / (rsq * rsum);
409#else
410 kn_ = dPi * emod / (rsq * rsq * rsum);
411#endif
412 ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
413 setInheritance(1,false);
414 setInheritance(2,false);
415 return true;
416 }
417 case kwArea: {
418 if (!userArea_) {
419 double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
420#ifdef THREED
421 userArea_ = rsq * rsq * dPi;
422#else
423 userArea_ = rsq * 2.0;
424#endif
425 }
426 return true;
427 }
428 }
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 QString knstr("kn");
494 bool ContactModelLinear::updateKn(const IContactMechanical *con) {
495 assert(con);
496 QVariant v1 = con->getEnd1()->getProperty(knstr);
497 QVariant 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 QString ksstr("ks");
513 bool ContactModelLinear::updateKs(const IContactMechanical *con) {
514 assert(con);
515 QVariant v1 = con->getEnd1()->getProperty(ksstr);
516 QVariant 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 QString fricstr("fric");
532 bool ContactModelLinear::updateFric(const IContactMechanical *con) {
533 assert(con);
534 QVariant v1 = con->getEnd1()->getProperty(fricstr);
535 QVariant 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 QString &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 QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
554 QRegExp rx(name,Qt::CaseInsensitive);
555 int idx = availableProperties.indexOf(rx)+1;
556 bool ret=false;
557
558 if (idx<=0)
559 return ret;
560
561 switch(idx) {
562 case kwKn: { //kn
563 if (inheritanceField_ & linKnMask)
564 ret = updateKn(c);
565 break;
566 }
567 case kwKs: { //ks
568 if (inheritanceField_ & linKsMask)
569 ret =updateKs(c);
570 break;
571 }
572 case kwFric: { //fric
573 if (inheritanceField_ & linFricMask)
574 updateFric(c);
575 break;
576 }
577 }
578 return ret;
579 }
580
581 void ContactModelLinear::updateEffectiveStiffness(ContactModelMechanicalState *) {
582 DVect2 ret(kn_,ks_);
583 // correction if viscous damping active
584 if (dpProps_) {
585 DVect2 correct(1.0);
586 if (dpProps_->dp_nratio_)
587 correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
588 if (dpProps_->dp_sratio_)
589 correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
590 ret /= (correct*correct);
591 }
592 effectiveTranslationalStiffness_ = ret;
593 }
594
595 bool ContactModelLinear::forceDisplacementLaw(ContactModelMechanicalState *state,const double ×tep) {
596 assert(state);
597
598 // Current overlap
599 double overlap = rgap_ - state->gap_;
600 // Relative translational increment
601 DVect trans = state->relativeTranslationalIncrement_;
602 // Correction factor to account for when the contact becomes newly active.
603 // We estimate the time of activity during the timestep when the contact has first
604 // become active and scale the forces accordingly.
605 double correction = 1.0;
606
607 // The contact was just activated from an inactive state
608 if (state->activated()) {
609 // Trigger the FISH callback if one is hooked up to the
610 // contact_activated event.
611 if (cmEvents_[fActivated] >= 0) {
612 // An FArray of QVariant is returned and these will be passed
613 // to the FISH function as an array of FISH symbols as the second
614 // argument to the FISH callback function.
615 auto c = state->getContact();
616 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
617 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
618 fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
619 }
620 // Calculate the correction factor.
621 if (lin_mode_ == 0 && trans.x()) {
622 correction = -1.0*overlap / trans.x();
623 if (correction < 0)
624 correction = 1.0;
625 }
626 }
627
628 DVect lin_F_old = lin_F_;
629
630 if (lin_mode_ == 0)
631 lin_F_.rx() = overlap * kn_; // absolute mode for normal force calculation
632 else
633 lin_F_.rx() -= correction * trans.x() * kn_; // incremental mode for normal force calculation
634
635 // Normal force can only be positive.
636 lin_F_.rx() = std::max(0.0,lin_F_.x());
637
638 // Calculate the shear force.
639 DVect sforce(0.0);
640 // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
641 // Loop over the shear components (note: the 0 component is the normal component)
642 // and calculate the shear force.
643 for (int i=1; i<dim; ++i)
644 sforce.rdof(i) = lin_F_.dof(i) - trans.dof(i) * ks_ * correction;
645
646 // The canFail flag corresponds to whether or not the contact can undergo non-linear
647 // force-displacement response. If the SOLVE ELASTIC command is given then the
648 // canFail state is set to FALSE. Otherwise it is always TRUE.
649 if (state->canFail_) {
650 // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
651 double crit = lin_F_.x() * fric_;
652 // The is the magnitude of the shear force.
653 double sfmag = sforce.mag();
654 // Sliding occurs when the magnitude of the shear force is greater than the
655 // critical value.
656 if (sfmag > crit) {
657 // Lower the shear force to the critical value for sliding.
658 double rat = crit / sfmag;
659 sforce *= rat;
660 // Handle the slip_change event if one has been hooked up. Sliding has commenced.
661 if (!lin_S_ && cmEvents_[fSlipChange] >= 0) {
662 auto c = state->getContact();
663 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
664 fish::Parameter() };
665 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
666 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
667 }
668 lin_S_ = true;
669 } else {
670 // Handle the slip_change event if one has been hooked up and
671 // the contact was previously sliding. Sliding has ceased.
672 if (lin_S_) {
673 if (cmEvents_[fSlipChange] >= 0) {
674 auto c = state->getContact();
675 std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
676 fish::Parameter((qint64)1) };
677 IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
678 fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
679 }
680 lin_S_ = false;
681 }
682 }
683 }
684
685 // Set the shear components of the total force.
686 for (int i=1; i<dim; ++i)
687 lin_F_.rdof(i) = sforce.dof(i);
688
689 // Account for dashpot forces if the dashpot structure has been defined.
690 if (dpProps_) {
691 dpProps_->dp_F_.fill(0.0);
692 double vcn(0.0), vcs(0.0);
693 // Calculate the damping coefficients.
694 setDampCoefficients(state->inertialMass_,&vcn,&vcs);
695 // First damp the shear components
696 for (int i=1; i<dim; ++i)
697 dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
698 // Damp the normal component
699 dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
700 // Need to change behavior based on the dp_mode.
701 if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
702 // Limit in tension if not bonded.
703 if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
704 dpProps_->dp_F_.rx() = - lin_F_.rx();
705 }
706 if (lin_S_ && dpProps_->dp_mode_ > 1) {
707 // Limit in shear if not sliding.
708 double dfn = dpProps_->dp_F_.rx();
709 dpProps_->dp_F_.fill(0.0);
710 dpProps_->dp_F_.rx() = dfn;
711 }
712 }
713
714 //Compute energies if energy tracking has been enabled.
715 if (state->trackEnergy_) {
716 assert(energies_);
717 energies_->estrain_ = 0.0;
718 if (kn_)
719 // Calculate the strain energy.
720 energies_->estrain_ = 0.5*lin_F_.x()*lin_F_.x()/kn_;
721 if (ks_) {
722 DVect s = lin_F_;
723 s.rx() = 0.0;
724 double smag2 = s.mag2();
725 // Add the shear component of the strain energy.
726 energies_->estrain_ += 0.5*smag2 / ks_;
727
728 if (lin_S_) {
729 // If sliding calculate the slip energy and accumulate it.
730 lin_F_old.rx() = 0.0;
731 DVect avg_F_s = (s + lin_F_old)*0.5;
732 DVect u_s_el = (s - lin_F_old) / ks_;
733 DVect u_s(0.0);
734 for (int i=1; i<dim; ++i)
735 u_s.rdof(i) = trans.dof(i);
736 energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
737 }
738 }
739 if (dpProps_) {
740 // Calculate damping energy (accumulated) if the dashpots are active.
741 energies_->edashpot_ -= dpProps_->dp_F_ | trans;
742 }
743 }
744
745 // This is just a sanity check to ensure, in debug mode, that the force isn't wonky.
746 assert(lin_F_ == lin_F_);
747 return true;
748 }
749
750 bool ContactModelLinear::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
751 // Account for thermal expansion in incremental mode
752 if (lin_mode_ == 0 || ts->gapInc_ == 0.0) return false;
753 DVect finc(0.0);
754 finc.rx() = kn_ * ts->gapInc_;
755 lin_F_ -= finc;
756 return true;
757 }
758
759 void ContactModelLinear::setForce(const DVect &v,IContact *c) {
760 lin_F(v);
761 if (v.x() > 0)
762 rgap_ = c->getGap() + v.x() / kn_;
763 }
764
765 void ContactModelLinear::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
766 // Only called for contacts with wall facets when the wall resolution scheme
767 // is set to full!
768 // Only do something if the contact model is of the same type
769 if (old->getContactModel()->getName().compare("linear",Qt::CaseInsensitive) == 0 && !isBonded()) {
770 ContactModelLinear *oldCm = (ContactModelLinear *)old;
771#ifdef THREED
772 // Need to rotate just the shear component from oldSystem to newSystem
773
774 // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
775 DVect axis = oldSystem.e1() & newSystem.e1();
776 double c, ang, s;
777 DVect re2;
778 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
779 axis = axis.unit();
780 c = oldSystem.e1()|newSystem.e1();
781 if (c > 0)
782 c = std::min(c,1.0);
783 else
784 c = std::max(c,-1.0);
785 ang = acos(c);
786 s = sin(ang);
787 double t = 1. - c;
788 DMatrix<3,3> rm;
789 rm.get(0,0) = t*axis.x()*axis.x() + c;
790 rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
791 rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
792 rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
793 rm.get(1,1) = t*axis.y()*axis.y() + c;
794 rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
795 rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
796 rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
797 rm.get(2,2) = t*axis.z()*axis.z() + c;
798 re2 = rm*oldSystem.e2();
799 }
800 else
801 re2 = oldSystem.e2();
802 // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
803 axis = re2 & newSystem.e2();
804 DVect2 tpf;
805 DMatrix<2,2> m;
806 if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
807 axis = axis.unit();
808 c = re2|newSystem.e2();
809 if (c > 0)
810 c = std::min(c,1.0);
811 else
812 c = std::max(c,-1.0);
813 ang = acos(c);
814 if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
815 ang *= -1;
816 s = sin(ang);
817 m.get(0,0) = c;
818 m.get(1,0) = s;
819 m.get(0,1) = -m.get(1,0);
820 m.get(1,1) = m.get(0,0);
821 tpf = m*DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
822 } else {
823 m.get(0,0) = 1.;
824 m.get(0,1) = 0.;
825 m.get(1,0) = 0.;
826 m.get(1,1) = 1.;
827 tpf = DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
828 }
829 DVect pforce = DVect(0,tpf.x(),tpf.y());
830#else
831 oldSystem;
832 newSystem;
833 DVect pforce = DVect(0,oldCm->lin_F_.y());
834#endif
835 for (int i=1; i<dim; ++i)
836 lin_F_.rdof(i) += pforce.dof(i);
837 if (lin_mode_ && oldCm->lin_mode_)
838 lin_F_.rx() = oldCm->lin_F_.x();
839 oldCm->lin_F_ = DVect(0.0);
840 if (dpProps_ && oldCm->dpProps_) {
841#ifdef THREED
842 tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
843 pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
844#else
845 pforce = oldCm->dpProps_->dp_F_;
846#endif
847 dpProps_->dp_F_ += pforce;
848 oldCm->dpProps_->dp_F_ = DVect(0.0);
849 }
850 if(oldCm->getEnergyActivated()) {
851 activateEnergy();
852 energies_->estrain_ = oldCm->energies_->estrain_;
853 energies_->edashpot_ = oldCm->energies_->edashpot_;
854 energies_->eslip_ = oldCm->energies_->eslip_;
855 oldCm->energies_->estrain_ = 0.0;
856 oldCm->energies_->edashpot_ = 0.0;
857 oldCm->energies_->eslip_ = 0.0;
858 }
859 rgap_ = oldCm->rgap_;
860 }
861 assert(lin_F_ == lin_F_);
862 }
863
864 void ContactModelLinear::setNonForcePropsFrom(IContactModel *old) {
865 // Only called for contacts with wall facets when the wall resolution scheme
866 // is set to full!
867 // Only do something if the contact model is of the same type
868 if (old->getName().compare("linear",Qt::CaseInsensitive) == 0 && !isBonded()) {
869 ContactModelLinear *oldCm = (ContactModelLinear *)old;
870 kn_ = oldCm->kn_;
871 ks_ = oldCm->ks_;
872 fric_ = oldCm->fric_;
873 lin_mode_ = oldCm->lin_mode_;
874 rgap_ = oldCm->rgap_;
875 userArea_ = oldCm->userArea_;
876
877 if (oldCm->dpProps_) {
878 if (!dpProps_)
879 dpProps_ = NEWC(dpProps());
880 dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
881 dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
882 dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
883 }
884 }
885 }
886
887 DVect ContactModelLinear::getForce() const {
888 DVect ret(lin_F_);
889 if (dpProps_)
890 ret += dpProps_->dp_F_;
891 return ret;
892 }
893
894 DAVect ContactModelLinear::getMomentOn1(const IContactMechanical *c) const {
895 if (!c)
896 return DAVect(0.0);
897 DVect force = getForce();
898 DAVect ret(0.0);
899 c->updateResultingTorqueOn1Local(force,&ret);
900 return ret;
901 }
902
903 DAVect ContactModelLinear::getMomentOn2(const IContactMechanical *c) const {
904 if (!c)
905 return DAVect(0.0);
906 DVect force = getForce();
907 DAVect ret(0.0);
908 c->updateResultingTorqueOn2Local(force,&ret);
909 return ret;
910 }
911
912 DAVect ContactModelLinear::getModelMomentOn1() const {
913 DAVect ret(0.0);
914 return ret;
915 }
916
917 DAVect ContactModelLinear::getModelMomentOn2() const {
918 DAVect ret(0.0);
919 return ret;
920 }
921
922 void ContactModelLinear::objectPropsTypes(std::vector<std::pair<QString,InfoTypes>> *ret) const {
923 ret->clear();
924 ret->push_back({"kn",ScalarInfo});
925 ret->push_back({"ks",ScalarInfo});
926 ret->push_back({"fric",ScalarInfo});
927 ret->push_back({"lin_force",VectorInfo});
928 ret->push_back({"lin_slip",ScalarInfo});
929 ret->push_back({"lin_mode",ScalarInfo});
930 ret->push_back({"rgap",ScalarInfo});
931 ret->push_back({"emod",ScalarInfo});
932 ret->push_back({"kratio",ScalarInfo});
933 ret->push_back({"dp_nratio",ScalarInfo});
934 ret->push_back({"dp_sratio",ScalarInfo});
935 ret->push_back({"dp_mode",ScalarInfo});
936 ret->push_back({"dp_force",VectorInfo});
937 ret->push_back({"user_area",ScalarInfo});
938 }
939
940 void ContactModelLinear::objectPropValues(std::vector<double> *ret,const IContact *c) const {
941 ret->clear();
942 ret->push_back(kn());
943 ret->push_back(ks());
944 ret->push_back(fric());
945 ret->push_back(lin_F().mag());
946 ret->push_back(lin_S());
947 ret->push_back(lin_mode());
948 ret->push_back(rgap());
949 ret->push_back(emod(c));
950 ret->push_back(kratio());
951 ret->push_back(dp_nratio());
952 ret->push_back(dp_sratio());
953 ret->push_back(dp_mode());
954 ret->push_back(dp_F().mag());
955 ret->push_back(getArea());
956 }
957
958 double ContactModelLinear::emod(const IContact *con) const {
959 const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
960 if (not c) return 0;
961 double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
962 double rsum(0.0);
963 if (c->getEnd1Curvature().y())
964 rsum += 1.0/c->getEnd1Curvature().y();
965 if (c->getEnd2Curvature().y())
966 rsum += 1.0/c->getEnd2Curvature().y();
967 if (userArea_) {
968#ifdef THREED
969 rsq = std::sqrt(userArea_ / dPi);
970#else
971 rsq = userArea_ / 2.0;
972#endif
973 rsum = rsq + rsq;
974 rsq = 1. / rsq;
975 }
976#ifdef TWOD
977 return (kn_ * rsum * rsq / 2.0);
978#else
979 return (kn_ * rsum * rsq * rsq) / dPi;
980#endif
981 }
982
983 void ContactModelLinear::setDampCoefficients(const double &mass,double *vcn,double *vcs) {
984 *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(kn_));
985 *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(ks_));
986 }
987
988} // namespace cmodelsxd
989// EoF
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Oct 31, 2024 |
