Smooth-Joint Model Implementation

See this page for the documentation of this contact model.

contactmodelsmoothjoint.h

  1#pragma once
  2// contactmodelsmoothjoint.h
  3
  4#include "contactmodel/src/contactmodelmechanical.h"
  5
  6#ifdef SMOOTHJOINT_LIB
  7#  define SMOOTHJOINT_EXPORT EXPORT_TAG
  8#elif defined(NO_MODEL_IMPORT)
  9#  define SMOOTHJOINT_EXPORT
 10#else
 11#  define SMOOTHJOINT_EXPORT IMPORT_TAG
 12#endif
 13
 14namespace cmodelsxd {
 15    using namespace itasca;
 16
 17    class ContactModelSmoothJoint : public ContactModelMechanical {
 18    public:
 19        enum PropertyKeys { 
 20               kwKn=1    
 21             , kwKs
 22             , kwFric   
 23             , kwDA
 24             , kwState
 25             , kwTen    
 26             , kwBCoh    
 27             , kwBFA
 28             , kwShear
 29             , kwRMul   
 30             , kwLarge
 31             , kwFn
 32             , kwFs     
 33             , kwGap
 34             , kwNorm
 35             , kwArea 
 36             , kwRad
 37             , kwUn     
 38             , kwUs
 39             , kwSlip
 40             , kwDip    
 41             , kwDD    
 42        };
 43
 44        enum FishCallEvents { fBondBreak=0, fSlipChange };
 45
 46        SMOOTHJOINT_EXPORT ContactModelSmoothJoint();
 47        SMOOTHJOINT_EXPORT virtual ~ContactModelSmoothJoint();
 48        virtual void     copy(const ContactModel *c) override;
 49        virtual void     archive(ArchiveStream &); 
 50        virtual QString  getName() const { return "smoothjoint"; }
 51        virtual void     setIndex(int i) { index_=i;}
 52        virtual int      getIndex() const {return index_;}
 53
 54        virtual QString  getProperties() const { 
 55            return "sj_kn"
 56                   ",sj_ks"
 57                   ",sj_fric"
 58                   ",sj_da"
 59                   ",sj_state"
 60                   ",sj_ten"
 61                   ",sj_coh"
 62                   ",sj_fa"
 63                   ",sj_shear"
 64                   ",sj_rmul"
 65                   ",sj_large"
 66                   ",sj_fn"
 67                   ",sj_fs"
 68                   ",sj_gap"
 69                   ",sj_unorm"
 70                   ",sj_area"
 71                   ",sj_radius"
 72                   ",sj_un"
 73                   ",sj_us"
 74                   ",sj_slip"
 75                   ",dip"
 76#ifdef THREED
 77                   ",ddir"
 78#endif
 79                   ; 
 80        }
 81
 82        virtual QString  getFishCallEvents() const { return "bond_break,slip_change"; }
 83        virtual QVariant getProperty(uint32 i,const IContact *con=0) const;
 84        virtual bool     getPropertyGlobal(uint32 i) const;
 85        virtual bool     setProperty(uint32 i,const QVariant &v,IContact *con=0);
 86        virtual bool     getPropertyReadOnly(uint32 i) const;
 87        virtual bool     supportsInheritance(uint32 i) const; 
 88        virtual bool     getInheritance(uint32 i) const { assert(i<32);  uint32 mask = to<uint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
 89        virtual void     setInheritance(uint32 i,bool b) { assert(i<32);  uint32 mask = to<uint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }
 90
 91        enum MethodKeys { 
 92              kwSetForce=1,kwBond,kwUnbond 
 93        };
 94
 95        virtual QString  getMethods() const { 
 96            return  "sj_setforce,bond,unbond";
 97        }
 98
 99        virtual QString  getMethodArguments(uint32 i) const; 
100        virtual bool     setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0); // Base 1 - returns true if timestep contributions need to be updated
101
102        virtual uint32     getMinorVersion() const;
103
104        enum EnergyKeys { kwEStrain=1,kwESlip};
105        virtual QString  getEnergies() const { return "energy-strain,energy-slip";}
106        virtual double   getEnergy(uint32 i) const;  // Base 1
107        virtual bool     getEnergyAccumulate(uint32 i) const; // Base 1
108        virtual void     setEnergy(uint32 i,const double &d); // Base 1
109        virtual void     activateEnergy() { if (energies_) return; energies_ = NEWC(Energies());}
110        virtual bool     getEnergyActivated() const {return (energies_ !=0);}
111
112        virtual bool    validate(ContactModelMechanicalState *state,const double &timestep);
113        virtual bool    endPropertyUpdated(const QString &name,const IContactMechanical *c);
114        virtual bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep);
115        virtual DVect2  getEffectiveTranslationalStiffness() const {return effectiveTranslationalStiffness_;}
116        virtual DAVect  getEffectiveRotationalStiffness() const {return DAVect(0.0);}
117
118        virtual ContactModelSmoothJoint *clone() const override { return NEWC(ContactModelSmoothJoint()); }
119        virtual double  getActivityDistance() const {return 0.0;}
120        virtual void    resetForcesAndMoments() { sj_Fn(0.0); sj_Fs(DVect(0.0)); }
121        virtual void    setForce(const DVect &v,IContact *);
122        virtual void    setArea(const double &) { throw Exception("The setArea method cannot be used with this contact model."); }
123        virtual double  getArea() const { return 0.0; }
124
125        virtual bool    isOKToDelete() const { return (!isBonded() && sj_large_); }
126
127        virtual bool    checkActivity(const double &gap) { return ( !sj_large_ || gap <= 0.0 || isBonded()); }
128        virtual bool    isSliding() const { return isjSliding_; }
129        virtual bool    isBonded() const { return sj_state_ == 3; }
130        virtual void    unbond() { sj_state_ = 0; }
131        virtual bool    hasNormal() const { return true; }
132        virtual DVect3  getNormal() const { return toVect3(sj_unorm_); }
133
134        void sj_kn(const double &d)    {sj_kn_ = d;}
135        void sj_ks(const double &d)    {sj_ks_ = d;}
136        void sj_fric(const double &d)  {sj_fric_ = d;}
137        void sj_da(const double &d)    {sj_da_ = d;}
138        void sj_state(const double &d) {sj_state_ = d;}
139        void sj_bns(const double &d)   {sj_bns_ = d;}
140        void sj_bcoh(const double &d)  {sj_bcoh_ = d;}
141        void sj_bfa(const double &d)   {sj_bfa_ = d;}
142        void dip(const double &d)      {dip_ = d;}
143        void sj_rmul(const double &d)  {sj_rmul_ = d;}
144        void sj_large(bool b)          {sj_large_ = b;}
145
146        const double & sj_kn()    const  {return sj_kn_;}    
147        const double & sj_ks()    const  {return sj_ks_;}    
148        const double & sj_fric()  const  {return sj_fric_;}  
149        const double & sj_da()    const  {return sj_da_;}    
150        int            sj_state() const  {return sj_state_;} 
151        const double & sj_bns()   const  {return sj_bns_;}   
152        const double & sj_bcoh()  const  {return sj_bcoh_;}  
153        const double & sj_bfa()   const  {return sj_bfa_;}   
154        const double & dip()      const  {return dip_;}      
155        const double & sj_rmul()  const  {return sj_rmul_;}  
156        bool           sj_large() const  {return sj_large_;} 
157
158#ifdef THREED
159        const double & dd() const          {return dd_;}
160        void           dd(const double &d) {dd_ =d;}
161#endif
162
163        void sj_unorm(const DVect &v)  {sj_unorm_ = v;}
164        void sj_A(const double &d)     {sj_A_ = d;}
165        void sj_rad(const double &d)   {sj_rad_ = d;}
166        void sj_gap(const double &d)   {sj_gap_ = d;}
167        void sj_Un(const double &d)    {sj_Un_ = d;}
168        void sj_Us(const DVect &v)     {sj_Us_ = v;}
169        void sj_Fn(const double &d)    {sj_Fn_ = d;}
170        void sj_Fs(const DVect &v)     {sj_Fs_ = v;}
171        void isjFlip(bool b)           {isjFlip_ = b;}
172        void isjSliding(bool b)        {isjSliding_ = b;}
173
174        const DVect  & sj_unorm()   const {return sj_unorm_;}  
175        const double & sj_A()       const {return sj_A_;}      
176        const double & sj_rad()     const {return sj_rad_;}    
177        const double & sj_gap()     const {return sj_gap_;}   
178        const double & sj_Un()      const {return sj_Un_;}     
179        const DVect  & sj_Us()      const {return sj_Us_;}     
180        const double & sj_Fn()      const {return sj_Fn_;}     
181        const DVect  & sj_Fs()      const {return sj_Fs_;}     
182        bool           isjFlip()    const {return isjFlip_;}   
183        bool           isjSliding() const {return isjSliding_;}
184
185        bool    hasEnergies() const {return energies_ ? true:false;}
186        double  estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
187        void    estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
188        double  eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
189        void    eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
190
191        uint32 inheritanceField() const {return inheritanceField_;}
192        void inheritanceField(uint32 i) {inheritanceField_ = i;}
193
194        const DVect2 & effectiveTranslationalStiffness()  const          {return effectiveTranslationalStiffness_;}
195        void           effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
196
197        bool    geomRecomp()       const {return geomRecomp_ ;}
198        void    geomRecomp(bool b)       {geomRecomp_ = b;}
199
200        /// Return the total force that the contact model holds.
201        virtual DVect    getForce(const IContactMechanical *) const;
202
203        /// Return the total moment on 1 that the contact model holds
204        virtual DAVect   getMomentOn1(const IContactMechanical *) const;
205
206        /// Return the total moment on 1 that the contact model holds
207        virtual DAVect   getMomentOn2(const IContactMechanical *) const;
208
209    private:
210        static int index_;
211
212        struct Energies {
213            Energies() : estrain_(0.0), eslip_(0.0) {}
214            double estrain_;  // elastic energy stored in contact 
215            double eslip_;    // work dissipated by friction 
216        };
217
218        bool   updateKn(const IContactMechanical *con);
219        bool   updateKs(const IContactMechanical *con);
220        bool   updateFric(const IContactMechanical *con);
221        void   updateAreaAndNormal(ContactModelMechanicalState *state);
222        void   updateAreaAndNormalContact(const IContact *c);
223        double calcBSS() const;
224
225        void   updateEffectiveTranslationalStiffness();
226
227        // property set fields
228        uint32 inheritanceField_;
229
230        // smooth joint model - contact model properties
231        double sj_kn_;     // normal stiffness
232        double sj_ks_;     // shear stiffness
233        double sj_fric_;   // Coulomb friction coefficient
234        double sj_da_;     // Dilation angle (stored in radians, returned/set in degrees)
235        int    sj_state_;  // Bond state - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (bonded)
236        double sj_bns_;    // Bond normal (tensile) strength
237        double sj_bcoh_;   // Bonded system cohesion
238        double sj_bfa_;    // Bonded system friction angle (stored in radians, returned/set in degrees)
239        double dip_;       // Dip angle (stored in radians, returned/set in degrees)
240        double sj_rmul_;   // Radius multiplier
241        bool   sj_large_;  // Large strain indicator
242
243        // Internal state variables
244        DVect  sj_unorm_;  // Normal to the plane
245        double sj_A_;      // Cross-sectional area
246        double sj_rad_;    // Radius
247        double sj_gap_;    // Gap - this can be user modified
248        double sj_Un_;     // Normal displacement
249        DVect  sj_Us_;     // Shear displacement
250        double sj_Fn_;     // Normal force
251        DVect  sj_Fs_;     // Shear force
252        bool   isjFlip_;   // In order to flip the order
253        bool   isjSliding_;// True if this is sliding
254
255        DVect2  effectiveTranslationalStiffness_;
256        bool    geomRecomp_; // If true then there must be a validate before FD!
257#ifdef THREED
258        double dd_;        // Dip direction (stored in radians, returned/set in degrees)
259#endif
260        //DAVect  effectiveRotationalStiffness_;
261
262        Energies *energies_;
263    };
264} // namespace itascaxd
265// EoF

Top

contactmodelsmoothjoint.cpp

  1// contactmodelsmoothjoint.cpp
  2#include "contactmodelsmoothjoint.h"
  3
  4#include "../version.txt"
  5#include "fish/src/parameter.h"
  6#include "utility/src/tptr.h"
  7#include "base/src/basetoqt.h"
  8
  9#include "kernel/interface/iprogram.h"
 10#include "module/interface/icontact.h"
 11#include "module/interface/icontactmechanical.h"
 12#include "module/interface/ifishcalllist.h"
 13#include "module/interface/ipiece.h"
 14#include "module/interface/ipiecemechanical.h"
 15
 16#ifdef SMOOTHJOINT_LIB
 17#ifdef _WIN32
 18  int __stdcall DllMain(void *,unsigned, void *)
 19  {
 20    return 1;
 21  }
 22#endif
 23
 24  extern "C" EXPORT_TAG const char *getName() 
 25  {
 26#if DIM==3
 27    return "contactmodelmechanical3dsmoothjoint";
 28#else
 29    return "contactmodelmechanical2dsmoothjoint";
 30#endif
 31  }
 32
 33  extern "C" EXPORT_TAG unsigned getMajorVersion()
 34  {
 35    return MAJOR_VERSION;
 36  }
 37
 38  extern "C" EXPORT_TAG unsigned getMinorVersion()
 39  {
 40    return MINOR_VERSION;
 41  }
 42
 43  extern "C" EXPORT_TAG void *createInstance() 
 44  {
 45    cmodelsxd::ContactModelSmoothJoint *m = NEWC(cmodelsxd::ContactModelSmoothJoint());
 46    return (void *)m;
 47  }
 48#endif // SMOOTHJOINT_EXPORTS
 49
 50namespace cmodelsxd {
 51    static const uint32 sjKnMask      = 0x00002; // Base 1!
 52    static const uint32 sjKsMask      = 0x00004;
 53    static const uint32 sjFricMask    = 0x00008;
 54
 55    using namespace itasca;
 56
 57    int ContactModelSmoothJoint::index_ = -1;
 58    uint32 ContactModelSmoothJoint::getMinorVersion() const { return MINOR_VERSION;}
 59
 60    ContactModelSmoothJoint::ContactModelSmoothJoint() : inheritanceField_(sjKnMask|sjKsMask|sjFricMask)
 61                                                       , sj_kn_(0.0), sj_ks_(0.0), sj_fric_(0.0)
 62                                                       , sj_da_(0.0), sj_state_(0), sj_bns_(0.0)
 63                                                       , sj_bcoh_(0.0), sj_bfa_(0.0), dip_(0.0)
 64                                                       , sj_rmul_(1.0), sj_large_(false)  
 65                                                       , sj_unorm_(0.0), sj_A_(0.0), sj_rad_(0.0)
 66                                                       , sj_gap_(0.0), sj_Un_(0.0), sj_Us_(0.0)
 67                                                       , sj_Fn_(0.0), sj_Fs_(0.0), isjFlip_(false)
 68                                                       , isjSliding_(false)
 69                                                       , effectiveTranslationalStiffness_(0.0)
 70                                                       , geomRecomp_(true)
 71#ifdef THREED
 72                                                       , dd_(0.0)
 73#endif
 74                                                       , energies_(0) {
 75//    setFromParent(ContactModelMechanicalList::instance()->find(getName()));
 76    }
 77
 78    ContactModelSmoothJoint::~ContactModelSmoothJoint() {
 79        if (energies_)
 80            delete energies_;
 81    }
 82
 83    void ContactModelSmoothJoint::archive(ArchiveStream &stream) {
 84        stream & sj_kn_;
 85        stream & sj_ks_;
 86        stream & sj_fric_;
 87        stream & sj_da_;
 88        stream & sj_state_;
 89        stream & sj_bns_;  
 90        stream & sj_bcoh_; 
 91        stream & sj_bfa_;  
 92        stream & dip_;     
 93        stream & sj_rmul_; 
 94        stream & sj_large_;
 95        stream & sj_unorm_;  
 96        stream & sj_A_;      
 97        stream & sj_rad_;    
 98        stream & sj_gap_;    
 99        stream & sj_Un_;     
100        stream & sj_Us_;     
101        stream & sj_Fn_;     
102        stream & sj_Fs_;     
103        stream & isjFlip_;   
104        stream & isjSliding_;
105#ifdef THREED
106        stream & dd_;
107#endif
108        if (stream.getArchiveState()==ArchiveStream::Save) {
109            bool b = false;
110            if (energies_) {
111                b = true;
112                stream & b;
113                stream & energies_->estrain_;
114                stream & energies_->eslip_;
115            }
116            else
117                stream & b;
118        } else {
119            bool b(false);
120            stream & b;
121            if (b) {
122                if (!energies_)
123                    energies_ = NEWC(Energies());
124                stream & energies_->estrain_;
125                stream & energies_->eslip_;
126            }
127        }
128
129        stream & inheritanceField_;
130        stream & geomRecomp_;
131        stream & effectiveTranslationalStiffness_;
132
133    }
134
135    void ContactModelSmoothJoint::copy(const ContactModel *cm) {
136        ContactModelMechanical::copy(cm);
137        const ContactModelSmoothJoint *in = dynamic_cast<const ContactModelSmoothJoint*>(cm);
138        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
139        sj_kn(in->sj_kn());
140        sj_ks(in->sj_ks());
141        sj_fric(in->sj_fric());
142        sj_da(in->sj_da());
143        sj_state(in->sj_state());
144        sj_bns(in->sj_bns());
145        sj_bcoh(in->sj_bcoh());
146        sj_bfa(in->sj_bfa());
147        dip(in->dip());
148        sj_rmul(in->sj_rmul()); 
149        sj_large(in->sj_large());
150#ifdef THREED
151        dd(in->dd());
152#endif
153        sj_unorm(in->sj_unorm());  
154        sj_A(in->sj_A());      
155        sj_rad(in->sj_rad());    
156        sj_gap(in->sj_gap());    
157        sj_Un(in->sj_Un());     
158        sj_Us(in->sj_Us());     
159        sj_Fn(in->sj_Fn());     
160        sj_Fs(in->sj_Fs());     
161        isjFlip(in->isjFlip());   
162        isjSliding(in->isjSliding());
163
164        if (in->hasEnergies()) {
165            if (!energies_)
166                energies_ = NEWC(Energies());
167            estrain(in->estrain());
168            eslip(in->eslip());
169        }
170        inheritanceField(in->inheritanceField());
171        effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
172        geomRecomp(in->geomRecomp());
173    }
174
175    QVariant ContactModelSmoothJoint::getProperty(uint32 i,const IContact *) const {
176        QVariant var;
177        switch (i) {
178        case kwKn:      return sj_kn_;
179        case kwKs:      return sj_ks_;
180        case kwFric:    return sj_fric_;
181        case kwDA:      return sj_da_/dDegrad;  // Returned in degrees
182        case kwState:   return sj_state_;
183        case kwTen:     return sj_bns_;
184        case kwBCoh:    return sj_bcoh_;
185        case kwBFA:     return sj_bfa_/dDegrad; // Returned in degrees
186        case kwShear:   return calcBSS();
187        case kwDip:     return dip_/dDegrad; // Returned in degrees
188#ifdef THREED
189        case kwDD:      return dd_/dDegrad;  // Returned in degrees
190#endif
191        case kwRMul:    return sj_rmul_;
192        case kwLarge:   return sj_large_;
193        case kwFn:      return sj_Fn_;
194        case kwFs: {
195                var.setValue(sj_Fs_);
196                return var;
197            }
198        case kwGap:     return sj_gap_;
199        case kwNorm: {
200                var.setValue(sj_unorm_);
201                return var;
202            }
203        case kwArea:    return sj_A_;
204        case kwRad:     return sj_rad_;
205        case kwUn:      return sj_Un_;
206        case kwUs: {
207                var.setValue(sj_Us_);
208                return var;
209            }
210        case kwSlip: return isjSliding_;
211        }
212        assert(0);
213        return QVariant();
214    }
215
216    bool ContactModelSmoothJoint::getPropertyGlobal(uint32 ) const {
217        return true;
218    }
219
220    bool ContactModelSmoothJoint::setProperty(uint32 i,const QVariant &v,IContact *c) {
221        switch (i) {
222        case kwKn: {
223                if (!v.canConvert<double>())
224                    throw Exception("sj_kn must be a double.");
225                double val = v.toDouble();
226                if (val<0.0)
227                    throw Exception("Negative sj_kn not allowed.");
228                sj_kn_ = val;
229                updateEffectiveTranslationalStiffness();
230                return true;
231            }
232        case kwKs: {
233                if (!v.canConvert<double>())
234                    throw Exception("sj_ks must be a double.");
235                double val = v.toDouble();
236                if (val<0.0)
237                    throw Exception("Negative sj_ks not allowed.");
238                sj_ks_ = val;  
239                updateEffectiveTranslationalStiffness();
240                return true;
241            }
242        case kwFric: {
243                if (!v.canConvert<double>())
244                    throw Exception("sj_fric must be a double.");
245                double val = v.toDouble();
246                if (val<0.0)
247                    throw Exception("Negative friction not allowed.");
248                sj_fric_ = val;  
249                return false;
250            }
251        case kwDA: {
252                if (!v.canConvert<double>())
253                    throw Exception("sj_da must be a double.");
254                sj_da_ = v.toDouble()*dDegrad; // Given in degrees
255                return false;
256            }
257        case kwState: {
258                if (!v.canConvert<uint32>())
259                    throw Exception("sj_state must be must be an integer between 0 and 3.");
260                int val = v.toInt();
261                if (val<0 || val>3)
262                    throw Exception("The bond state must be an integer between 0 and 3.");
263                sj_state_ = val;  
264                return false;
265            }
266        case kwTen: {
267                if (!v.canConvert<double>())
268                    throw Exception("sj_ten must be a double.");
269                double val = v.toDouble();
270                if (val<0.0)
271                    throw Exception("Negative bond normal (tensile) strength not allowed.");
272                sj_bns_ = val;  
273                return false;
274            }
275        case kwBCoh: {
276                if (!v.canConvert<double>())
277                    throw Exception("sj_coh must be a double.");
278                double val = v.toDouble();
279                if (val<0.0)
280                    throw Exception("Negative bond system cohesion not allowed.");
281                sj_bcoh_ = val;  
282                return false;
283            }
284        case kwBFA: {
285                if (!v.canConvert<double>())
286                    throw Exception("sj_bfa must be a double.");
287                sj_bfa_ = v.toDouble()*dDegrad; // Given in degrees
288                return false;
289            }
290        case kwDip: {
291                if (!v.canConvert<double>())
292                    throw Exception("dip must be a double.");
293                dip_ = v.toDouble()*dDegrad; // Given in degrees
294                if (c)
295                    updateAreaAndNormalContact(c);
296                else 
297                    geomRecomp_ = true;
298                return true;
299            }
300        case kwDD: {
301#ifdef THREED
302                if (!v.canConvert<double>())
303                    throw Exception("ddir must be a double.");
304                dd_ = v.toDouble()*dDegrad; // Given in degrees
305                if (c)
306                    updateAreaAndNormalContact(c);
307                else 
308                    geomRecomp_ = true;
309#endif
310                return true;
311            }
312        case kwRMul: {
313                if (!v.canConvert<double>())
314                    throw Exception("sj_rmul must be a double.");
315                double val = v.toDouble();
316                if (val<0.0)
317                    throw Exception("Negative radius multiplier not allowed.");
318                sj_rmul_ = val;  
319                if (!geomRecomp_) geomRecomp_ = true;
320                return false;
321            }
322        case kwLarge: {
323                if (!v.canConvert<bool>())
324                    throw Exception("Large-strain flag must be a boolean.");
325                sj_large_ = v.toBool();  
326                return false;
327            }
328        case kwFn: {
329                 if (!v.canConvert<double>())
330                    throw Exception("sj_fn must be a double.");
331               sj_Fn_ = v.toDouble();
332                return false;
333            }
334        case kwFs: {
335                if (!v.canConvert<DVect>())
336                    throw Exception("sj_fs must be a vector.");
337                sj_Fs_ = v.value<DVect>();
338                return false;
339            }
340        case kwGap: {
341                if (!v.canConvert<double>())
342                    throw Exception("sj_gap must be a double.");
343                sj_gap_ = v.toDouble();
344                return false;
345            }
346
347         
348        // Following are read-only !
349        case kwNorm:
350        case kwRad:
351        case kwShear:
352        case kwUn:
353        case kwUs:
354        case kwSlip:
355            return false;
356        }
357        assert(0);
358        return false;
359    }
360
361    bool ContactModelSmoothJoint::getPropertyReadOnly(uint32 i) const {
362        switch (i) {
363        case kwNorm:
364        case kwRad:
365        case kwUn:
366        case kwUs:
367        case kwSlip:
368        case kwShear:
369        case kwArea:
370            return true;
371        default:
372            break;
373        }
374        return false;
375    }
376
377    bool ContactModelSmoothJoint::supportsInheritance(uint32 i) const {
378        switch (i) {
379        case kwKn:
380        case kwKs:
381        case kwFric:
382            return true;
383        default:
384            break;
385        }
386        return false;
387    }
388
389    QString  ContactModelSmoothJoint::getMethodArguments(uint32 i) const {
390        QString def = QString();
391        switch (i) {
392        case kwBond: 
393            return "gap";
394        case kwUnbond: 
395            return "gap";
396        }
397        return def;
398    }
399
400    bool ContactModelSmoothJoint::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
401        IContactMechanical *c(convert_getcast<IContactMechanical>(con));
402        switch (i) {
403        case kwSetForce: {
404                DVect gf = c->getGlobalForce();
405                if (isjFlip_)
406                    gf *= -1.0;
407                sj_Fn_ = gf | sj_unorm_;
408                return false;
409            }
410        case kwBond: {
411                if (sj_state_ == 3) return false;
412                double mingap = -1.0 * limits<double>::max();
413                double maxgap = 0;
414                if (vl.at(0).canConvert<double>()) 
415                    maxgap = vl.at(0).toDouble();
416                else if (vl.at(0).canConvert<DVect2>()) {
417                    DVect2 value = vl.at(0).value<DVect2>();
418                    mingap = value.minComp();
419                    maxgap = value.maxComp();
420                } else if (!vl.at(0).isNull())
421                    throw Exception("gap value %1 not recognized in method bond in contact model %2.",vl.at(0),getName());
422                if ( sj_gap_ >= mingap && sj_gap_ <= maxgap) {
423                    sj_state_ = 3;
424                    return true;
425                }
426                return false;
427        }
428        case kwUnbond: {
429                if (sj_state_ == 0) return false;
430                double mingap = -1.0 * limits<double>::max();
431                double maxgap = 0;
432                if (vl.at(0).canConvert<double>()) 
433                    maxgap = vl.at(0).toDouble();
434                else if (vl.at(0).canConvert<DVect2>()) {
435                    DVect2 value = vl.at(0).value<DVect2>();
436                    mingap = value.minComp();
437                    maxgap = value.maxComp();
438                }
439                else if (!vl.at(0).isNull())
440                    throw Exception("gap value %1 not recognized in method unbond in contact model %2.",vl.at(0),getName());
441                if ( sj_gap_ >= mingap && sj_gap_ <= maxgap) {
442                    sj_state_ = 0;
443                    return true;
444                }
445                return false;
446            }
447
448        }
449        return false;
450    }
451
452    double ContactModelSmoothJoint::getEnergy(uint32 i) const {
453        double ret(0.0);
454        if (!energies_)
455            return ret;
456        switch (i) {
457        case kwEStrain:  return energies_->estrain_;
458        case kwESlip:    return energies_->eslip_;
459        }
460        assert(0);
461        return ret;
462    }
463
464    bool ContactModelSmoothJoint::getEnergyAccumulate(uint32 i) const {
465        bool ret(false);
466        if (!energies_) return ret;
467        switch (i) {
468        case kwEStrain:  return false;
469        case kwESlip:    return true;
470        }
471        assert(0);
472        return ret;
473    }
474
475    void ContactModelSmoothJoint::setEnergy(uint32 i,const double &d) {
476        if (!energies_) return;
477        switch (i) {
478        case kwEStrain:  energies_->estrain_ = d; return;  
479        case kwESlip:    energies_->eslip_   = d; return;
480        }
481        assert(0);
482        return;
483    }
484
485    bool ContactModelSmoothJoint::validate(ContactModelMechanicalState *state,const double &) {
486        assert(state);
487        const IContactMechanical *c = state->getMechanicalContact(); 
488        assert(c);
489
490        if (state->trackEnergy_)
491            activateEnergy();
492
493        // The radius multiplier has been set previously so calculate the sj_A_
494        // Need to do this regardless of whether or not the radius multiplier has been updated as this is the 
495        // first place with the contact state to get the ball radii!
496        updateAreaAndNormal(state);
497
498        if (inheritanceField_ & sjKnMask)    
499            updateKn(c);
500        if (inheritanceField_ & sjKsMask)      
501            updateKs(c);
502        if (inheritanceField_ & sjFricMask)   
503            updateFric(c);
504
505        updateEffectiveTranslationalStiffness();
506
507       return checkActivity(state->gap_);
508    }
509
510    void ContactModelSmoothJoint::updateAreaAndNormal(ContactModelMechanicalState *state) {
511        updateAreaAndNormalContact(state->getContact());
512    }
513
514    void ContactModelSmoothJoint::updateAreaAndNormalContact(const IContact *c) {
515        assert(c);
516        if (geomRecomp_) geomRecomp_ = false;
517        // Use the maximum value of the curvature here - not the min!
518        auto *mc = convert_getcast<IContactMechanical>(c);
519        sj_rad_ = 1.0 / std::max(mc->getEnd1Curvature().y(),mc->getEnd2Curvature().y());
520        sj_rad_ = sj_rmul_ * sj_rad_;
521#ifdef THREED
522        sj_A_ = dPi * sj_rad_ * sj_rad_;
523        sj_unorm_.rx() = sin(dip_) * sin(dd_);
524        sj_unorm_.ry() = sin(dip_) * cos(dd_);
525        sj_unorm_.rz() = cos(dip_);
526#else
527        sj_A_ = 2.0 * sj_rad_; // Assumes thickness of 1 in 2D
528        sj_unorm_.rx() = sin(dip_);
529        sj_unorm_.ry() = cos(dip_);
530#endif
531        DVect nc = toVect(c->getNormal());
532        // Set the flip boolean so that the ordering is correct from side1 to side2
533        isjFlip_ = ( ((nc | sj_unorm_) >= 0.0 ) ? false : true );
534    }
535
536    static const QString kn("sj_kn");
537    bool ContactModelSmoothJoint::updateKn(const IContactMechanical *con) {
538        assert(con);
539        QVariant v1 = con->getEnd1()->getProperty(kn);
540        QVariant v2 = con->getEnd2()->getProperty(kn);
541        if (!v1.isValid() || !v2.isValid())
542            return false;
543        double kn1 = v1.toDouble();
544        double kn2 = v2.toDouble();
545        double val = sj_kn_;
546        if (kn1 && kn2)
547            sj_kn_ = kn1*kn2/(kn1+kn2);
548        else if (kn1)
549            sj_kn_ = kn1;
550        else if (kn2)
551            sj_kn_ = kn2;
552        return ( (sj_kn_ != val) );
553    }
554
555    static const QString ks("sj_ks");
556    bool ContactModelSmoothJoint::updateKs(const IContactMechanical *con) {
557        assert(con);
558        QVariant v1 = con->getEnd1()->getProperty(ks);
559        QVariant v2 = con->getEnd2()->getProperty(ks);
560        if (!v1.isValid() || !v2.isValid())
561            return false;
562        double kn1 = v1.toDouble();
563        double kn2 = v2.toDouble();
564        double val = sj_ks_;
565        if (kn1 && kn2)
566            sj_ks_ = kn1*kn2/(kn1+kn2);
567        else if (kn1)
568            sj_ks_ = kn1;
569        else if (kn2)
570            sj_ks_ = kn2;
571        return ( (sj_ks_ != val) );
572    }
573
574    static const QString fric("sj_fric");
575    bool ContactModelSmoothJoint::updateFric(const IContactMechanical *con) {
576        assert(con);
577        QVariant v1 = con->getEnd1()->getProperty(fric);
578        QVariant v2 = con->getEnd2()->getProperty(fric);
579        if (!v1.isValid() || !v2.isValid())
580            return false;
581        double fric1 = std::max(0.0,v1.toDouble());
582        double fric2 = std::max(0.0,v2.toDouble());
583        double val = sj_fric_;
584        sj_fric_ = std::min(fric1,fric2);
585        return ( (sj_fric_ != val) );
586    }
587
588    bool ContactModelSmoothJoint::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
589        assert(c);
590        QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
591        QRegularExpression rx(name, QRegularExpression::CaseInsensitiveOption);
592        int idx = availableProperties.indexOf(rx) + 1;
593        bool ret=false;
594
595        if (idx<=0)
596            return ret;
597         
598        switch(idx) {
599        case kwKn:  { //sj_kn_
600                if (inheritanceField_ & sjKnMask)
601                    ret = updateKn(c);
602                break;
603            }
604        case kwKs:  { //sj_ks_
605                if (inheritanceField_ & sjKsMask)
606                    ret =updateKs(c);
607                break;
608            }
609        case kwFric:  { //sj_fric_
610                if (inheritanceField_ & sjFricMask)
611                    updateFric(c);
612                break;
613            }
614        }
615
616        if (ret)
617            updateEffectiveTranslationalStiffness();
618        return ret;
619    }
620
621    void ContactModelSmoothJoint::updateEffectiveTranslationalStiffness() {
622        effectiveTranslationalStiffness_ = DVect2(sj_kn_,sj_ks_)*sj_A_;
623    }
624     
625    bool ContactModelSmoothJoint::forceDisplacementLaw(ContactModelMechanicalState *state,const double &) {
626        assert(state);
627         
628        if (geomRecomp_)
629            updateAreaAndNormal(state);
630
631        // Skip out if this should not be active
632        if (sj_large_ && state->gap_ > 0.0 && sj_state_ != 3) {
633            sj_Fn_ = 0.0;
634            sj_Fs_ = DVect(0.0);
635            sj_gap_ = 0.0;
636            sj_Un_ = 0.0;
637            sj_Us_ = DVect(0.0);
638            return false;
639        }
640
641        // Get the translational increment in the global system
642        CAxes localSys = state->getMechanicalContact()->getContact()->getLocalSystem();
643        DVect tInc = localSys.toGlobal(state->relativeTranslationalIncrement_);
644         
645        // Now have the translational increment and the normal in the global coordinate system
646        // Compute the normal/shear displacement increments along the sj
647        if (isjFlip_)
648            tInc *= -1.0;
649        double nInc = tInc | sj_unorm_;
650        DVect sInc = tInc - sj_unorm_ * nInc;
651        sj_Un_ -= nInc;
652        sj_Us_ += sInc;
653
654        double nElInc(0.0);
655        DVect sElInc(0.0);
656
657        double g0 = sj_gap_, g1 = sj_gap_ + nInc;
658        sj_gap_ = g1;
659         
660        if (!state->canFail_ || sj_state_ == 3 )  { // Bonded
661            nElInc  = nInc;
662            sElInc = sInc;
663        } else {
664            if ( g0 <= 0.0 ) {
665                if ( g1 <= 0.0 ) {
666                    nElInc  = nInc;
667                    sElInc = sInc;
668                } else { // g1 > 0.0
669                    double xi = -g0 / (g1 - g0);
670                    nElInc  = nInc  * xi;
671                    sElInc = sInc * xi;
672                }
673            } else { // g0 > 0.0
674                if ( g1 >= 0.0 ) {
675                    nElInc = 0.0;
676                    sElInc.fill(0.0);
677                } else { // g1 < 0.0
678                    double xi = -g0 / (g1 - g0);
679                    nElInc  = nInc  * (1.0 - xi);
680                    sElInc = sInc * (1.0 - xi);
681                }
682            }
683        }
684
685        double del_Fn = -sj_kn_ * sj_A_ * nElInc ;
686        sj_Fn_ += del_Fn ;
687        int slideChange(-1);
688        DVect sj_Fs_old = sj_Fs_;
689
690        if ( state->canFail_ && sj_state_ < 3 ) {  // Coulomb sliding with dilation
691            if ( sj_Fn_ <= 0.0 ) {
692                sj_Fn_ = 0.0; 
693                sj_Fs_.fill(0.0);
694            } else {
695                DVect del_Fs = sElInc * (-sj_ks_ * sj_A_);
696                sj_Fs_ += del_Fs;
697                double max_Fs = sj_Fn_ * sj_fric_;
698                double magFs = sj_Fs_.mag();
699                if ( magFs > max_Fs ) { // sliding
700                    if (!isjSliding_) {
701                        isjSliding_ = true;
702                        slideChange = 0;
703                    }
704                    if ( sj_ks_ > 0.0 )  // dilation
705                        sj_Fn_ += ( (magFs - max_Fs) / sj_ks_ ) * sj_kn_ * tan(sj_da_);
706                    double rat = max_Fs / magFs ;
707                    sj_Fs_ *= rat ;
708                } else {
709                    if (isjSliding_) {
710                        isjSliding_ = false;
711                        slideChange = 1;
712                    }
713                }
714            }
715        } else { // bonded behavior
716            if ( state->canFail_ && sj_Fn_ <= -sj_bns_ * sj_A_) {
717                sj_state_ = 1;  
718                if (cmEvents_[fBondBreak] >= 0) {
719                    auto c = state->getContact();
720                    std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
721                                                         fish::Parameter((qint64)sj_state_),
722                                                         fish::Parameter(sj_bns_)         };
723                    IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
724                    fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
725                }
726                sj_Fn_ = 0.0; 
727                sj_Fs_.fill(0.0);
728            } else {
729                DVect del_Fs = sElInc * (-sj_ks_ * sj_A_);
730                sj_Fs_ += del_Fs;
731                double magFs = sj_Fs_.mag();
732                double ss = calcBSS();
733                if ( state->canFail_ && magFs >= ss * sj_A_) { // break in shear
734                    sj_state_ = 2;  
735                    if (cmEvents_[fBondBreak] >= 0) {
736                        auto c = state->getContact();
737                        std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
738                                                             fish::Parameter((qint64)sj_state_),
739                                                             fish::Parameter(ss)                 };
740                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
741                        fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
742                    }
743
744                    if ( sj_Fn_ < 0.0 ) {
745                        sj_Fn_ = 0.0; 
746                        sj_Fs_.fill(0.0); 
747                    }  else { // was in compression // was in tension
748                        double max_Fs = sj_Fn_ * sj_fric_ ;
749                        if ( magFs > max_Fs) { // sliding, but no dilation
750                            if (!isjSliding_) {
751                                isjSliding_ = true;
752                                slideChange = 0;
753                            }
754                            double rat = max_Fs / magFs ;
755                            sj_Fs_ *= rat ;
756                        } else {
757                            if (isjSliding_ == true) {
758                                isjSliding_ = false;
759                                slideChange = 1;
760                            }
761                        }
762                    }
763                }
764            }
765        }
766        if (slideChange >= 0 && cmEvents_[fSlipChange] >= 0) {
767            auto c = state->getContact();
768            std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
769                                                 fish::Parameter((qint64)slideChange) };
770            IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
771            fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
772        }
773
774        // Have updated the normal and shear forces so need to put them into the contact local
775        // coordinate system and update the forces 
776        DVect Fj = sj_unorm_ * sj_Fn_ + sj_Fs_;   
777        if (isjFlip_)
778            Fj *= -1.0;
779
780        // Return the correct activity status
781        bool isactive = true;
782        if (sj_large_ && state->gap_ > 0.0 && sj_state_ != 3) {
783            sj_Fn_ = 0.0;
784            sj_Fs_ = DVect(0.0);
785            sj_gap_ = 0.0;
786            sj_Un_ = 0.0;
787            sj_Us_ = DVect(0.0);
788            isactive = false;
789        }
790
791        // update energies
792       if (state->trackEnergy_) {
793            assert(energies_);
794            energies_->estrain_ =  0.0;
795            if (sj_kn_)
796                energies_->estrain_ = 0.5*sj_Fn_*sj_Fn_/sj_kn_;
797            if (sj_ks_) {
798                double smag2 = sj_Fs_.mag2();
799                energies_->estrain_ += 0.5*smag2 / sj_ks_;
800                if (isjSliding_) {
801                    DVect avg_F_s = (sj_Fs_old + sj_Fs_)*0.5;
802                    DVect u_s_el =  (sj_Fs_ - sj_Fs_old) / sj_ks_;
803                    energies_->eslip_ -= std::min(0.0,(avg_F_s | (sInc + u_s_el)));
804                }
805            }
806        }
807
808        return isactive ;
809    }
810
811    void ContactModelSmoothJoint::setForce(const DVect &v,IContact *c) {
812        // this is in the local coordinate system
813        CAxes localSys = c->getLocalSystem();
814        DVect globForce = localSys.toGlobal(v);
815        if (isjFlip_)
816            globForce *= -1.0;
817        sj_Fn_ = (sj_unorm_ | globForce);
818        sj_Fs_ = globForce - sj_unorm_ * sj_Fn_;
819    }
820
821    DVect ContactModelSmoothJoint::getForce(const IContactMechanical *c) const {
822        CAxes localSys = c->getContact()->getLocalSystem();
823        DVect Fj = sj_unorm_ * sj_Fn_ + sj_Fs_;   
824        if (isjFlip_)
825            Fj *= -1.0;
826        DVect ret(localSys.toLocal(Fj));
827        return ret;
828    }
829
830    DAVect ContactModelSmoothJoint::getMomentOn1(const IContactMechanical *c) const {
831        DVect force = getForce(c);
832        DAVect ret(0.0);
833        c->updateResultingTorqueOn1Local(force,&ret);
834        return ret;
835    }
836
837    DAVect ContactModelSmoothJoint::getMomentOn2(const IContactMechanical *c) const {
838        DVect force = getForce(c);
839        DAVect ret(0.0);
840        c->updateResultingTorqueOn2Local(force,&ret);
841        return ret;
842    }
843
844    double ContactModelSmoothJoint::calcBSS() const {
845        if (sj_A_ > 0) {
846            double dSigma = sj_Fn_ / sj_A_;
847            return dSigma >= -sj_bns_ ? sj_bcoh_ + (tan(sj_bfa_) * dSigma)
848                                                    : sj_bcoh_ + (tan(sj_bfa_) * (-sj_bns_));
849        }
850        else 
851            return 0.0;
852    }
853
854
855} // namespace itascaxd
856// EoF

Top