Adhesive Rolling Resistance Linear Model

See this page for the documentation of this contact model.

contactmodelarrlinear.h

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

Top

contactmodelarrlinear.cpp

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

Top