Rolling Resistance Linear Contact Model Implementation

See this page for the documentation of this contact model.

contactmodelrrlinear.h

  1#pragma once
  2// contactmodelrrlinear.h
  3
  4#include "contactmodel/src/contactmodelmechanical.h"
  5
  6#ifdef rrlinear_LIB
  7#  define rrlinear_EXPORT EXPORT_TAG
  8#elif defined(NO_MODEL_IMPORT)
  9#  define rrlinear_EXPORT
 10#else
 11#  define rrlinear_EXPORT IMPORT_TAG
 12#endif
 13
 14namespace cmodelsxd {
 15    using namespace itasca;
 16
 17    class ContactModelRRLinear : public ContactModelMechanical {
 18    public:
 19        // Constructor: Set default values for contact model properties.
 20        rrlinear_EXPORT ContactModelRRLinear();
 21        // Destructor, called when contact is deleted: free allocated memory, etc.
 22        rrlinear_EXPORT virtual ~ContactModelRRLinear();
 23        // Contact model name (used as keyword for commands and FISH).
 24        QString  getName() const override { return "rrlinear"; }
 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            , kwResFric
 56            , kwResMoment
 57            , kwResS
 58            , kwResKr
 59            , kwUserArea
 60        };
 61        // Contact model property names in a comma separated list. The order corresponds with
 62        // the order of the PropertyKeys enumerator above. One can visualize any of these 
 63        // properties in PFC automatically. 
 64        QString  getProperties() const override {
 65            return "kn"
 66                   ",ks"
 67                   ",fric"
 68                   ",lin_force"
 69                   ",lin_slip"
 70                   ",lin_mode"
 71                   ",rgap"
 72                   ",emod"
 73                   ",kratio"
 74                   ",dp_nratio"
 75                   ",dp_sratio"
 76                   ",dp_mode"
 77                   ",dp_force"
 78                   ",rr_fric"
 79                   ",rr_moment"
 80                   ",rr_slip"
 81                   ",rr_kr"
 82                   ",user_area" ;
 83        }
 84        // Enumerator for the energies.
 85        enum EnergyKeys { 
 86            kwEStrain=1
 87          , kwERRStrain
 88          , kwESlip
 89          , kwERRSlip
 90          , kwEDashpot
 91        };
 92        // Contact model energy names in a comma separated list. The order corresponds with
 93        // the order of the EnergyKeys enumerator above. 
 94        QString  getEnergies() const override {
 95            return "energy-strain"
 96                   ",energy-rrstrain"
 97                   ",energy-slip"
 98                   ",energy-rrslip"
 99                   ",energy-dashpot";
100        }
101        // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
102        double   getEnergy(uint32 i) const override;
103        // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1) 
104        // returns wther or not the estrain energy is accumulated which is false).
105        bool     getEnergyAccumulate(uint32 i) const override;
106        // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
107        void     setEnergy(uint32 i,const double &d) override; // Base 1
108        // Activate the energy. This is only called if the energy tracking is enabled. 
109        void     activateEnergy() override { if (energies_) return; energies_ = NEWC(Energies());}
110        // Returns whether or not the energy tracking has been enabled for this contact.
111        bool     getEnergyActivated() const override {return (energies_ != 0);}
112
113        // Enumerator for contact model related FISH callback events. 
114        enum FishCallEvents {
115            fActivated=0
116            ,fSlipChange
117        };
118        // Contact model FISH callback event names in a comma separated list. The order corresponds with
119        // the order of the FishCallEvents enumerator above. 
120        QString  getFishCallEvents() const override {
121            return 
122                "contact_activated"
123                ",slip_change"; 
124        }
125
126        // Return the specified contact model property.
127        QVariant getProperty(uint32 i,const IContact *) const override;
128        // The return value denotes whether or not the property corresponds to the global
129        // or local coordinate system (TRUE: global system, FALSE: local system). The
130        // local system is the contact-plane system (nst) defined as follows.
131        // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
132        // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
133        // and tc are unit vectors in directions of the nst axes.
134        // This is used when rendering contact model properties that are vectors.
135        bool     getPropertyGlobal(uint32 i) const override;
136        // Set the specified contact model property, ensuring that it is of the correct type
137        // and within the correct range --- if not, then throw an exception.
138        // The return value denotes whether or not the update has affected the timestep
139        // computation (by having modified the translational or rotational tangent stiffnesses).
140        // If true is returned, then the timestep will be recomputed.
141        bool     setProperty(uint32 i,const QVariant &v,IContact *) override;
142        // The return value denotes whether or not the property is read-only
143        // (TRUE: read-only, FALSE: read-write).
144        bool     getPropertyReadOnly(uint32 i) const override;
145
146        // The return value denotes whether or not the property is inheritable
147        // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
148        // the endPropertyUpdated method.
149        bool     supportsInheritance(uint32 i) const override;
150        // Return whether or not inheritance is enabled for the specified property.
151        bool     getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
152        // Set the inheritance flag for the specified property.
153        void     setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }
154
155        // Enumerator for contact model methods.
156        enum MethodKeys { kwDeformability=1,kwArea};
157        // Contact model methoid names in a comma separated list. The order corresponds with
158        // the order of the MethodKeys enumerator above.  
159        QString  getMethods() const override { return "deformability,area";}
160        // Return a comma seprated list of the contact model method arguments (base 1).
161        QString  getMethodArguments(uint32 i) const override;
162        // Set contact model method arguments (base 1). 
163        // The return value denotes whether or not the update has affected the timestep
164        // computation (by having modified the translational or rotational tangent stiffnesses).
165        // If true is returned, then the timestep will be recomputed.
166        bool     setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override;
167
168        // Prepare for entry into ForceDispLaw. The validate function is called when:
169        // (1) the contact is created, (2) a property of the contact that returns a true via
170        // the setProperty method has been modified and (3) when a set of cycles is executed
171        // via the {cycle N} command.
172        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
173        bool    validate(ContactModelMechanicalState *state,const double &timestep) override;
174        // The endPropertyUpdated method is called whenever a surface property (with a name
175        // that matches an inheritable contact model property name) of one of the contacting
176        // pieces is modified. This allows the contact model to update its associated
177        // properties. The return value denotes whether or not the update has affected
178        // the time step computation (by having modified the translational or rotational
179        // tangent stiffnesses). If true is returned, then the time step will be recomputed.  
180        bool    endPropertyUpdated(const QString &name,const IContactMechanical *c) override;
181        // The forceDisplacementLaw function is called during each cycle. Given the relative
182        // motion of the two contacting pieces (via
183        //   state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
184        //   state->relativeAngularIncrement_       (Dtt, Dtbs, Dtbt)
185        //     Ddn  : relative normal-displacement increment, Ddn > 0 is opening
186        //     Ddss : relative  shear-displacement increment (s-axis component)
187        //     Ddst : relative  shear-displacement increment (t-axis component)
188        //     Dtt  : relative twist-rotation increment
189        //     Dtbs : relative  bend-rotation increment (s-axis component)
190        //     Dtbt : relative  bend-rotation increment (t-axis component)
191        //       The relative displacement and rotation increments:
192        //         Dd = Ddn*nc + Ddss*sc + Ddst*tc
193        //         Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
194        //       where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
195        //       [see {Table 1: Contact State Variables} in PFC Model Components:
196        //       Contacts and Contact Models: Contact Resolution]
197        // ) and the contact properties, this function must update the contact force and
198        // moment.
199        //   The force_ is acting on piece 2, and is expressed in the local coordinate system
200        //   (defined in getPropertyGlobal) such that the first component positive denotes
201        //   compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
202        //   in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to 
203        //   get the total moment. 
204        // The return value indicates the contact activity status (TRUE: active, FALSE:
205        // inactive) during the next cycle.
206        // Additional information:
207        //   * If state->activated() is true, then the contact has just become active (it was
208        //     inactive during the previous time step).
209        //   * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
210        //     the forceDispLaw handle the case of {state->canFail_ == true}.
211        bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) override;
212        bool    thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
213        // The getEffectiveXStiffness functions return the translational and rotational
214        // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
215        // the translational tangent shear stiffness is zero (but this stiffness reduction
216        // is typically ignored when computing a stable time step). If the contact model
217        // includes a dashpot, then the translational stiffnesses must be increased (see
218        // Potyondy (2009)).
219        //   [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
220        //   Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
221        //   December 7, 2009.]
222        DVect2  getEffectiveTranslationalStiffness() const override { return effectiveTranslationalStiffness_; }
223        DAVect  getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness_;}
224
225        // Return a new instance of the contact model. This is used in the CMAT
226        // when a new contact is created. 
227        ContactModelRRLinear *clone() const override { return NEWC(ContactModelRRLinear()); }
228        // The getActivityDistance function is called by the contact-resolution logic when
229        // the CMAT is modified. Return value is the activity distance used by the
230        // checkActivity function.
231        double              getActivityDistance() const override {return rgap_;}
232        // The isOKToDelete function is called by the contact-resolution logic when...
233        // Return value indicates whether or not the contact may be deleted.
234        // If TRUE, then the contact may be deleted when it is inactive.
235        // If FALSE, then the contact may not be deleted (under any condition).
236        bool                isOKToDelete() const override { return !isBonded(); }
237        // Zero the forces and moments stored in the contact model. This function is called
238        // when the contact becomes inactive.
239        void                resetForcesAndMoments() override {
240            lin_F(DVect(0.0)); 
241            dp_F(DVect(0.0)); 
242            res_M(DAVect(0.0));
243            if (energies_) {
244                energies_->estrain_ = 0.0;
245                energies_->errstrain_ = 0.0;
246            }
247        }
248        void                setForce(const DVect &v,IContact *c) override;
249        void                setArea(const double &d) override { userArea_ = d; }
250        double              getArea() const override { return userArea_; }
251        // The checkActivity function is called by the contact-resolution logic when...
252        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
253        bool     checkActivity(const double &gap) override { return  gap <= rgap_; }
254        // Returns the sliding state (FALSE is returned if not implemented).
255        bool     isSliding() const override { return lin_S_; }
256        // Returns the bonding state (FALSE is returned if not implemented).
257        bool     isBonded() const override { return false; }
258        // Both of these methods are called only for contacts with facets where the wall 
259        // resolution scheme is set the full. In such cases one might wish to propagate 
260        // contact state information (e.g., shear force) from one active contact to another. 
261        // See the Faceted Wall section in the documentation. 
262        void     propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
263        void     setNonForcePropsFrom(IContactModel *oldCM) override;
264        /// Return the total force that the contact model holds.
265        DVect    getForce(const IContactMechanical *) const override;
266        /// Return the total moment on 1 that the contact model holds
267        DAVect   getMomentOn1(const IContactMechanical *) const override;
268        /// Return the total moment on 1 that the contact model holds
269        DAVect   getMomentOn2(const IContactMechanical *) const override;
270   
271        // Methods to get and set properties. 
272        const double & kn() const {return kn_;}
273        void           kn(const double &d) {kn_=d;}
274        const double & ks() const {return ks_;}
275        void           ks(const double &d) {ks_=d;}
276        const double & fric() const {return fric_;}
277        void           fric(const double &d) {fric_=d;}
278        const DVect &  lin_F() const {return lin_F_;}
279        void           lin_F(const DVect &f) { lin_F_=f;}
280        bool           lin_S() const {return lin_S_;}
281        void           lin_S(bool b) { lin_S_=b;}
282        uint32           lin_mode() const {return lin_mode_;}
283        void           lin_mode(uint32 i) { lin_mode_= i;}
284        const double & rgap() const {return rgap_;}
285        void           rgap(const double &d) {rgap_=d;}
286
287        bool     hasDamping() const {return dpProps_ ? true : false;}
288        double   dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
289        void     dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
290        double   dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
291        void     dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
292        int      dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
293        void     dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
294        DVect    dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
295        void     dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
296
297        bool    hasEnergies() const {return energies_ ? true:false;}
298        double  estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
299        void    estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
300        double  errstrain() const {return hasEnergies() ? energies_->errstrain_: 0.0;}
301        void    errstrain(const double &d) { if(!hasEnergies()) return; energies_->errstrain_=d;}
302        double  eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
303        void    eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
304        double  errslip() const {return hasEnergies() ? energies_->errslip_: 0.0;}
305        void    errslip(const double &d) { if(!hasEnergies()) return; energies_->errslip_=d;}
306        double  edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
307        void    edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
308
309        uint32 inheritanceField() const {return inheritanceField_;}
310        void inheritanceField(uint32 i) {inheritanceField_ = i;}
311
312        const DVect2 & effectiveTranslationalStiffness()  const          {return effectiveTranslationalStiffness_;}
313        void           effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
314        const DAVect & effectiveRotationalStiffness()  const             {return effectiveRotationalStiffness_;}
315        void           effectiveRotationalStiffness(const DAVect &v )    {effectiveRotationalStiffness_=v;}
316
317        // Rolling resistance methods
318        const double & res_fric() const {return res_fric_;}
319        void           res_fric(const double &d) {res_fric_=d;}
320        const DAVect & res_M() const               {return res_M_;}
321        void           res_M(const DAVect &f)       { res_M_=f;}
322        bool           res_S() const {return res_S_;}
323        void           res_S(bool b) { res_S_=b;}
324        const double & kr() const {return kr_;}
325        void           kr(const double &d) {kr_=d;}
326        const double & fr() const {return fr_;}
327        void           fr(const double &d) {fr_=d;}
328
329    private:
330        // Index - used internally by PFC. Should be set to -1 in the cpp file. 
331        static int index_;
332
333        // Structure to store the energies. 
334        struct Energies {
335            Energies() : estrain_(0.0),errstrain_(0.0),eslip_(0.0),errslip_(0.0),edashpot_(0.0) {}
336            double estrain_;   // elastic energy stored in linear group 
337            double errstrain_; // elastic energy stored in rolling resistance group
338            double eslip_;     // work dissipated by friction 
339            double errslip_;   // work dissipated by rolling resistance friction 
340            double edashpot_;  // work dissipated by dashpots
341        };
342
343        // Structure to store dashpot quantities. 
344        struct dpProps {
345            dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)) {}
346            double dp_nratio_;     // normal viscous critical damping ratio
347            double dp_sratio_;     // shear  viscous critical damping ratio
348            int    dp_mode_;      // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
349            DVect  dp_F_;  // Force in the dashpots
350        };
351
352        bool   updateKn(const IContactMechanical *con);
353        bool   updateKs(const IContactMechanical *con);
354        bool   updateFric(const IContactMechanical *con);
355        bool   updateResFric(const IContactMechanical *con);
356
357        void   updateStiffness(ContactModelMechanicalState *state);
358
359        void   setDampCoefficients(const double &mass,double *vcn,double *vcs);
360
361        // Contact model inheritance fields.
362        uint32 inheritanceField_;
363
364        // Effective translational stiffness.
365        DVect2  effectiveTranslationalStiffness_;
366        DAVect  effectiveRotationalStiffness_;      // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)
367
368        // linear model properties
369        double      kn_;        // Normal stiffness
370        double      ks_;        // Shear stiffness
371        double      fric_;      // Coulomb friction coefficient
372        DVect       lin_F_;     // Force carried in the linear model
373        bool        lin_S_;     // The current slip state
374        uint32        lin_mode_;  // Specifies absolute (0) or incremental (1) calculation mode 
375        double      rgap_;      // Reference gap 
376        dpProps *   dpProps_;   // The viscous properties
377
378        // rolling resistance properties
379        double res_fric_;       // rolling friction coefficient
380        DAVect res_M_;          // moment (bending only)         
381        bool   res_S_;          // The current rolling resistance slip state
382        double kr_;             // bending rotational stiffness (read-only, calculated internaly) 
383        double fr_;             // rolling friction coefficient (rbar*res_fric_) (calculated internaly, not a property) 
384
385        double      userArea_;  // Area as specified by the user 
386
387        Energies *   energies_; // The energies
388
389    };
390} // namespace cmodelsxd
391// EoF

Top

contactmodelrrlinear.cpp

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

Top