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