Hysteretic Model Implementation

See this page for the documentation of this contact model.

contactmodelhysteretic.h

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
#pragma once
// contactmodelhysteretic.h

#include "contactmodel/src/contactmodelmechanical.h"

#ifdef HYSTERETIC_LIB
#  define HYSTERETIC_EXPORT EXPORT_TAG
#elif defined(NO_MODEL_IMPORT)
#  define HYSTERETIC_EXPORT
#else
#  define HYSTERETIC_EXPORT IMPORT_TAG
#endif

namespace cmodelsxd {
    using namespace itasca;

    class ContactModelHysteretic : public ContactModelMechanical {
    public:
        enum PropertyKeys { kwHzShear=1
                          , kwHzPoiss                            
                          , kwFric   
                          , kwHzF
                          , kwHzS
                          , kwHzSd
                          , kwHzAlpha
                          , kwDpMode
                          , kwDpEn
                          , kwDpEnMin
                          , kwDpF 
                         };
       
        HYSTERETIC_EXPORT ContactModelHysteretic();
        HYSTERETIC_EXPORT virtual ~ContactModelHysteretic();
        virtual void                copy(const ContactModel *c);
        virtual void                archive(ArchiveStream &); 
  
        virtual QString  getName() const { return "hysteretic"; }
        virtual void     setIndex(int i) { index_=i;}
        virtual int      getIndex() const {return index_;}
      
        virtual QString  getProperties() const { 
            return "hz_shear"
                   ",hz_poiss"
                   ",fric"
                   ",hz_force"
                   ",hz_slip"
                   ",hz_mode"
                   ",hz_alpha"
                   ",dp_mode"
                   ",dp_en"
                   ",dp_enmin"
                   ",dp_force"
            ;
        }
  
        enum EnergyKeys { kwEStrain=1,kwESlip,kwEDashpot};
        virtual QString  getEnergies() const { return "energy-strain,energy-slip,energy-dashpot";}
        virtual double   getEnergy(uint i) const;  // Base 1
        virtual bool     getEnergyAccumulate(uint i) const; // Base 1
        virtual void     setEnergy(uint i,const double &d); // Base 1
        virtual void     activateEnergy() { if (energies_) return; energies_ = NEWC(Energies());}
        virtual bool     getEnergyActivated() const {return (energies_ !=0);}
        
        enum FishCallEvents { fActivated=0, fSlipChange};
        virtual QString  getFishCallEvents() const { return "contact_activated,slip_change"; }
        virtual QVariant getProperty(uint i,const IContact *) const;
        virtual bool     getPropertyGlobal(uint i) const;
        virtual bool     setProperty(uint i,const QVariant &v,IContact *);
        virtual bool     getPropertyReadOnly(uint i) const;
        
        virtual bool     supportsInheritance(uint i) const; 
        virtual bool     getInheritance(uint i) const { assert(i<32); quint32 mask = to<quint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
        virtual void     setInheritance(uint i,bool b) { assert(i<32); quint32 mask = to<quint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }
                
        virtual uint     getMinorVersion() const;
        
        virtual bool    validate(ContactModelMechanicalState *state,const double &timestep);
        virtual bool    endPropertyUpdated(const QString &name,const IContactMechanical *c);
        virtual bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep);
        virtual DVect2  getEffectiveTranslationalStiffness() const { return effectiveTranslationalStiffness_;}
        virtual DAVect  getEffectiveRotationalStiffness() const { return DAVect(0.0);}
        
        virtual ContactModelHysteretic *clone() const { return NEWC(ContactModelHysteretic()); }
        virtual double              getActivityDistance() const {return 0.0;}
        virtual bool                isOKToDelete() const { return !isBonded(); }
        virtual void                resetForcesAndMoments() { hz_F(DVect(0.0)); dp_F(DVect(0.0));  if (energies_) energies_->estrain_ = 0.0;}
        virtual void                setForce(const DVect &v,IContact *) { hz_F(v); }
        virtual void                setArea(const double &) { throw Exception("The setArea method cannot be used with this contact model."); }

        virtual bool     checkActivity(const double &gap) { return gap <= 0.0; }
        
        virtual bool     isSliding() const { return hz_slip_; }
        virtual bool     isBonded() const { return false; }
        virtual void     propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes());
        virtual void     setNonForcePropsFrom(IContactModel *oldCM);
        
        const double & hz_shear() const {return hz_shear_;}
        void           hz_shear(const double &d) {hz_shear_=d;}
        const double & hz_poiss() const {return hz_poiss_;}
        void           hz_poiss(const double &d) {hz_poiss_=d;}
        const double & fric() const {return fric_;}
        void           fric(const double &d) {fric_=d;}
        uint           hz_mode() const {return hz_mode_;}
        void           hz_mode(uint i) {hz_mode_=i;}
        const DVect &  hz_F() const {return hz_F_;}
        void           hz_F(const DVect &f) { hz_F_=f;}
        bool           hz_S() const {return hz_slip_;}
        void           hz_S(bool b) { hz_slip_=b;}
        const double & hz_alpha() const {return hz_alpha_;}
        void           hz_alpha(const double &d) {hz_alpha_=d;}
        int            dp_mode() const {return dp_mode_;}
        void           dp_mode(int i) {dp_mode_=i;}
        const double & dp_en() const {return dp_en_;}
        void           dp_en(const double &d) {dp_en_=d;}
        const double & dp_enmin() const {return dp_enmin_;}
        void           dp_enmin(const double &d) {dp_enmin_=d;}
        const DVect &  dp_F() const {return dp_F_;}
        void           dp_F(const DVect &f) { dp_F_=f;}
        const double & hn() const {return hn_;}
        void           hn(const double &d) {hn_=d;}
        const double & hs() const {return hs_;}
        void           hs(const double &d) {hs_=d;}
        const double & vni() const {return vni_;}
        void           vni(const double &d) {vni_=d;}
        double         pfac() const {return pfac_;}
        void           pfac(int i) {pfac_=i;}
                
        bool    hasEnergies() const {return energies_ ? true:false;}
        double  estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
        void    estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
        double  eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
        void    eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
        double  edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
        void    edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
        
        uint inheritanceField() const {return inheritanceField_;}
        void inheritanceField(uint i) {inheritanceField_ = i;}
        
        const DVect2 & effectiveTranslationalStiffness()  const          {return effectiveTranslationalStiffness_;}
        void           effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
  
        /// Return the total force that the contact model holds.
        virtual DVect    getForce(const IContactMechanical *) const;

        /// Return the total moment on 1 that the contact model holds
        virtual DAVect   getMomentOn1(const IContactMechanical *) const;

        /// Return the total moment on 1 that the contact model holds
        virtual DAVect   getMomentOn2(const IContactMechanical *) const;

    private:
        static int index_;
        
        bool   updateStiffCoef(const IContactMechanical *con);
        bool   updateEndStiffCoef(const IContactMechanical *con);
        bool   updateEndFric(const IContactMechanical *con);
        void   updateEffectiveStiffness(ContactModelMechanicalState *state);
        // inheritance fields
        quint32 inheritanceField_;
        
        // hertz model
        double      hz_shear_;  // Shear modulus
        double      hz_poiss_;  // Poisson ratio
        double      fric_;       // Coulomb friction coefficient
        DVect       hz_F_;      // Force carried in the hertz model
        bool        hz_slip_;   // the current sliding state
        uint        hz_mode_;   // specifies down-scaling of the shear force when normal unloading occurs 
        double      hz_alpha_;  // alpha exponent
        int         dp_mode_;   // Damping scheme mode
        double      dp_en_;     // normal restitution coefficient
        double      dp_enmin_;  // minimal normal restitution coefficient in default mode
        DVect       dp_F_;      // Damping Force
                
        // energies
        struct Energies {
            Energies() : estrain_(0.0), eslip_(0.0),edashpot_(0.0) {}
            double estrain_;  // elastic energy stored in contact 
            double eslip_;    // work dissipated by friction 
            double edashpot_;    // work dissipated by dashpots
        };
        Energies *   energies_;    
                
        double hn_;                               // normal stiffness coefficient
        double hs_;                               // shear stiffness coefficient
        DVect2 effectiveTranslationalStiffness_;  // effective stiffness
        double vni_;                              // normal impact velocity
        double pfac_;                             // previous damping factor

    };

} // namespace cmodelsxd
// EoF

Top

contactmodelhysteretic.cpp

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
// contactmodelhysteretic.cpp
#include "contactmodelhysteretic.h"

#include "module/interface/icontactmechanical.h"
#include "module/interface/icontact.h"

#include "module/interface/ipiecemechanical.h"
#include "module/interface/ipiece.h"
#include "../version.txt"

#include "module/interface/ifishcalllist.h"
#include "utility/src/tptr.h"
#include "shared/src/mathutil.h"

#include "kernel/interface/iprogram.h"
#include "module/interface/icontactthermal.h"
#include "contactmodel/src/contactmodelthermal.h"


#ifdef HYSTERETIC_LIB
    int __stdcall DllMain(void *,unsigned, void *) {
        return 1;
    }

    extern "C" EXPORT_TAG const char *getName() {
#if DIM==3
        return "contactmodelmechanical3dhysteretic";
#else
        return "contactmodelmechanical2dhysteretic";
#endif
    }

    extern "C" EXPORT_TAG unsigned getMajorVersion() {
        return MAJOR_VERSION;
    }

    extern "C" EXPORT_TAG unsigned getMinorVersion() {
        return MINOR_VERSION;
    }

    extern "C" EXPORT_TAG void *createInstance() {
        cmodelsxd::ContactModelHysteretic *m = NEWC(cmodelsxd::ContactModelHysteretic());
        return (void *)m;
    }
#endif // HYSTERETIC_LIB

namespace cmodelsxd {

    static const quint32 shearMask   = 0x00002;
    static const quint32 poissMask   = 0x00004;
    static const quint32 fricMask    = 0x00008;
  
    using namespace itasca;
  
    int ContactModelHysteretic::index_ = -1;
    UInt ContactModelHysteretic::getMinorVersion() const { return MINOR_VERSION;}
  
    ContactModelHysteretic::ContactModelHysteretic() : inheritanceField_(shearMask|poissMask|fricMask) 
                                            , hz_shear_(0.0)
                                            , hz_poiss_(0.0)
                                            , fric_(0.0)
                                            , hz_F_(DVect(0.0))
                                            , hz_slip_(false)
                                            , hz_mode_(0)
                                            , hz_alpha_(1.5)
                                            , dp_mode_(0)
                                            , dp_en_(1.0)
                                            , dp_enmin_(0.0)
                                            , dp_F_(DVect(0.0))
                                            , energies_(0)
                                            , hn_(0.0)
                                            , hs_(0.0)
                                            , effectiveTranslationalStiffness_(DVect2(0.0)) 
                                            , vni_(0.0) 
                                            , pfac_(0.0) 
    {
    }
  
    ContactModelHysteretic::~ContactModelHysteretic() {
        if (energies_)
          delete energies_;
    }
  
    void ContactModelHysteretic::archive(ArchiveStream &stream) {
        stream & hz_shear_;
        stream & hz_poiss_;
        stream & fric_;
        stream & hz_F_;
        stream & hz_slip_;
        stream & hz_mode_;
        stream & hz_alpha_;
        stream & dp_mode_;
        stream & dp_en_;
        stream & dp_enmin_;
        stream & dp_F_;
        stream & hn_;
        stream & hs_;
        stream & vni_;
        stream & pfac_;
    
        if (stream.getArchiveState()==ArchiveStream::Save) {
            bool b = false;
            if (energies_) {
                b = true;
                stream & b;
                stream & energies_->estrain_;
                stream & energies_->eslip_;
                stream & energies_->edashpot_;
            } else
                stream & b;
        } else {
            bool b(false);
            stream & b;
            if (b) {
                if (!energies_)
                    energies_ = NEWC(Energies());
                stream & energies_->estrain_;
                stream & energies_->eslip_;
                stream & energies_->edashpot_;
            }
        }
    
        stream & inheritanceField_;
        stream & effectiveTranslationalStiffness_;
    }
  
    void ContactModelHysteretic::copy(const ContactModel *cm) {
        ContactModelMechanical::copy(cm);
        const ContactModelHysteretic *in = dynamic_cast<const ContactModelHysteretic*>(cm);
        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
        
        hz_shear(in->hz_shear());
        hz_poiss(in->hz_poiss());
        fric(in->fric());
        hz_F(in->hz_F());
        hz_S(in->hz_S());
        hz_mode(in->hz_mode());
        hz_alpha(in->hz_alpha());
        dp_mode(in->dp_mode());
        dp_en(in->dp_en());
        dp_enmin(in->dp_enmin());
        hn(in->hn());
        hs(in->hs()); 
        vni(in->vni()); 
        pfac(in->pfac()); 
        if (in->hasEnergies()) {
            if (!energies_)
                energies_ = NEWC(Energies());
            estrain(in->estrain());
            eslip(in->eslip());
            edashpot(in->edashpot());
        }      
        inheritanceField(in->inheritanceField());
        effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
    }
  
    QVariant ContactModelHysteretic::getProperty(uint i,const IContact *) const {
      QVariant var;
      switch (i) {
          case kwHzShear:  return hz_shear_;
          case kwHzPoiss:  return hz_poiss_;
          case kwFric:      return fric_;
          case kwHzF:      var.setValue(hz_F_); return var;
          case kwHzS:      return hz_slip_;
          case kwHzSd:     return hz_mode_;
          case kwHzAlpha:  return hz_alpha_;
          case kwDpMode:    return dp_mode_;
          case kwDpEn:      return dp_en_;
          case kwDpEnMin:      return dp_enmin_;
          case kwDpF:       var.setValue(dp_F_); return var;
      }
      assert(0);
      return QVariant();
    }
  
    bool ContactModelHysteretic::getPropertyGlobal(uint i) const {
        switch (i) {
            case kwHzF: return false;   
            case kwDpF: return false;   
        }
        return true;
    }
  
    bool ContactModelHysteretic::setProperty(uint i,const QVariant &v,IContact *) {
        switch (i) {
            case kwHzShear: {
                if (!v.canConvert<double>())
                    throw Exception("hz_shear must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative shear modulus (hz_shear) not allowed.");
                hz_shear_ = val;  
                return true;
            }
            case kwHzPoiss: {
                if (!v.canConvert<double>())
                    throw Exception("hz_poiss must be a double.");
                double val(v.toDouble());
                if (val<=-1.0 || val>0.5)
                    throw Exception("Poisson ratio (hz_poiss) must be in range (-1.0,0.5] ");
                hz_poiss_ = val;  
                return true;
            }
            case kwFric: {
                if (!v.canConvert<double>())
                    throw Exception("fric must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative fric not allowed.");
                fric_ = val;  
                return false;
            }
            case kwHzSd: {
               if (!v.canConvert<uint>())
                    throw Exception("hz_mode must be 0 or 1.");
                uint val(v.toUInt());
                if (val >1)
                    throw Exception("hz_mode must be 0 or 1.");
                hz_mode_ = val;
                return false;
            }
            case kwHzAlpha: {
                if (!v.canConvert<double>())
                    throw Exception("hz_alpha must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Alpha exponent (hz_alpha) must be positive.");
                hz_alpha_ = val;
                return false;
            }
            case kwDpMode: {
               if (!v.canConvert<int>())
                    throw Exception("dp_mode must be an integer.");
                int val(v.toInt());
                dp_mode_ = val;
                return false;
            }
            case kwDpEn: {
                if (!v.canConvert<double>())
                    throw Exception("dp_en must be a double.");
                double val(v.toDouble());
                if (val<0.0 || val>1.0)
                    throw Exception("Restitution coefficient (dp_en) must be in range [0.0,1.0] ");
                dp_en_ = val;
                return false;
            }
            case kwDpEnMin: {
                if (!v.canConvert<double>())
                    throw Exception("dp_enmin must be a double.");
                double val(v.toDouble());
                if (val<0.0 || val>1.0)
                    throw Exception("Minimal restitution coefficient (dp_enmin) must be in range [0.0,1.0] ");
                dp_enmin_ = val;
                return false;
            }
        }
        return false;
    }
  
    bool ContactModelHysteretic::getPropertyReadOnly(uint i) const {
        switch (i) {
            case kwHzF:
            case kwHzS:
            case kwDpF:
                return true;
            default:
                break;
        }
        return false;
    }
  
    bool ContactModelHysteretic::supportsInheritance(uint i) const {
        switch (i) {
            case kwHzShear:
            case kwHzPoiss:
            case kwFric:
                return true;
            default:
                break;
        }
        return false;
    }
  
    double ContactModelHysteretic::getEnergy(uint i) const {
        double ret(0.0);
        if (!energies_)
            return ret;
        switch (i) {
            case kwEStrain:  return energies_->estrain_;
            case kwESlip:    return energies_->eslip_;
            case kwEDashpot: return energies_->edashpot_;
        }
        assert(0);
        return ret;
    }
  
    bool ContactModelHysteretic::getEnergyAccumulate(uint i) const {
        switch (i) {
            case kwEStrain:  return false;
            case kwESlip:    return true;
            case kwEDashpot: return true;
        }
        assert(0);
        return false;
    }
  
    void ContactModelHysteretic::setEnergy(uint i,const double &d) {
        if (!energies_) return;
        switch (i) {
            case kwEStrain:  energies_->estrain_ = d; return;  
            case kwESlip:    energies_->eslip_   = d; return;
            case kwEDashpot: energies_->edashpot_   = d; return;
        }
        assert(0);
        return;
    }
  
    bool ContactModelHysteretic::validate(ContactModelMechanicalState *state,const double &) {
        assert(state);
        const IContactMechanical *c = state->getMechanicalContact(); 
        assert(c);
    
        if (state->trackEnergy_)
            activateEnergy();
    
        updateStiffCoef(c);
        if ((inheritanceField_ & shearMask) || (inheritanceField_ & poissMask))
            updateEndStiffCoef(c);
    
        if (inheritanceField_ & fricMask)
            updateEndFric(c);
    
        updateEffectiveStiffness(state);
        return checkActivity(state->gap_);
    }
  
    bool ContactModelHysteretic::updateStiffCoef(const IContactMechanical *con) {
        double hnold = hn_;
        double hsold = hs_;
        double c12 = con->getEnd1Curvature().y();
        double c22 = con->getEnd2Curvature().y();
        double reff = c12+c22;
        if (reff == 0.0) 
            throw Exception("Hysteretic contact model undefined for 2 non-curved surfaces");
        reff = 2.0 /reff;
        hn_ = 2.0/3.0 * (hz_shear_/(1 -hz_poiss_)) * sqrt(2.0*reff);
        hs_ = (2.0*pow(hz_shear_*hz_shear_*3.0*(1-hz_poiss_)*(reff),1.0/3.0)) / (2.0- hz_poiss_);
        return ( (hn_ != hnold) || (hs_ != hsold) );
    }

    static const QString gstr("hz_shear");
    static const QString nustr("hz_poiss");
    bool ContactModelHysteretic::updateEndStiffCoef(const IContactMechanical *con) {
        assert(con);
        double g1 = hz_shear_;
        double g2 = hz_shear_;
        double nu1 = hz_poiss_;
        double nu2 = hz_poiss_;
        QVariant vg1 = con->getEnd1()->getProperty(gstr);
        QVariant vg2 = con->getEnd2()->getProperty(gstr);
        QVariant vnu1 = con->getEnd1()->getProperty(nustr);
        QVariant vnu2 = con->getEnd2()->getProperty(nustr);
        if (vg1.isValid() && vg2.isValid()) {
            g1 = vg1.toDouble();
            g2 = vg2.toDouble();
            if (g1 < 0.0 || g2 < 0.0) 
                throw Exception("Negative shear modulus not allowed in Hysteretic contact model");
        }    
        if (vnu1.isValid() && vnu2.isValid()) {
            nu1 = vnu1.toDouble();
            nu2 = vnu2.toDouble();
            if (nu1 <= -1.0 || nu1 > 0.5 || nu2 <= -1.0 || nu2 > 0.5) 
                throw Exception("Poisson ratio should be in range (-1.0,0.5] in Hysteretic contact model");
        }
        if (g1*g2 == 0.0) return false;
        double es = 1.0 / ((1.0-nu1) / (2.0*g1) + (1.0-nu2) / (2.0*g2));
        double gs = 1.0 / ((2.0-nu1) / g1 + (2.0-nu2) /g2);
        hz_poiss_ = (4.0*gs-es)/(2.0*gs-es);
        hz_shear_ = 2.0*gs*(2-hz_poiss_);
        if (hz_shear_ < 0.0) 
            throw Exception("Negative shear modulus not allowed in Hysteretic contact model");
        if (hz_poiss_ <= -1.0 || hz_poiss_ > 0.5) 
            throw Exception("Poisson ratio should be in range (-1.0,0.5] in Hysteretic contact model");
        return updateStiffCoef(con);
    }
  
    static const QString fricstr("fric");
    bool ContactModelHysteretic::updateEndFric(const IContactMechanical *con) {
        assert(con);
        QVariant v1 = con->getEnd1()->getProperty(fricstr);
        QVariant v2 = con->getEnd2()->getProperty(fricstr);
        if (!v1.isValid() || !v2.isValid())
            return false;
        double fric1 = std::max(0.0,v1.toDouble());
        double fric2 = std::max(0.0,v2.toDouble());
        double val = fric_;
        fric_ = std::min(fric1,fric2);
        return ( (fric_ != val) );
    }
  
    bool ContactModelHysteretic::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
        assert(c);
        QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",QString::SkipEmptyParts);
        QRegExp rx(name,Qt::CaseInsensitive);
        int idx = availableProperties.indexOf(rx)+1;
        bool ret=false;
    
        if (idx<=0)
            return ret;
        
        switch(idx) {
            case kwHzShear: {
                if (inheritanceField_ & shearMask)
                    ret = updateEndStiffCoef(c);
                break;
            }
            case kwHzPoiss: {
                if (inheritanceField_ & poissMask)
                    ret = updateEndStiffCoef(c);
                break;
            }
            case kwFric: {
                if (inheritanceField_ & fricMask)
                    ret = updateEndFric(c);
                break;
            }
        }
        return ret;
    }
  
    void ContactModelHysteretic::updateEffectiveStiffness(ContactModelMechanicalState *state) {
        effectiveTranslationalStiffness_ = DVect2(hn_,hs_);
        if (state->gap_ >= 0.0) return;
        double overlap = - state->gap_;
        double kn = 1.5*hn_*sqrt(overlap);
        double ks = hs_ * pow(hz_F_.x(),(1.0/3.0));
        DVect2 ret(kn,ks);
        effectiveTranslationalStiffness_ = ret;
    }
    
    bool ContactModelHysteretic::forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) {
        assert(state);
        bool firstActive = false;
        if (state->activated()) {
            if (cmEvents_[fActivated] >= 0) {
                FArray<QVariant,2> arg;
                QVariant v;
                IContact * c = const_cast<IContact*>(state->getContact());
                TPtr<IThing> t(c->getIThing());
                v.setValue(t);
                arg.push_back(v);
                IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
                fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
            }
            firstActive = true; 
        }

        double overlap = - state->gap_;
        DVect trans = state->relativeTranslationalIncrement_;
#ifdef THREED
        DVect norm(trans.x(),0.0,0.0);
#else
        DVect norm(trans.x(),0.0);
#endif

        DAVect ang  = state->relativeAngularIncrement_;
        DVect hz_F_old = hz_F_;
        hz_F_ = DVect(0.0);
        dp_F_ = DVect(0.0);
        double ks_old = hs_ * pow(hz_F_old.x(),(1.0/3.0));
        DVect fs_old = hz_F_old;
        fs_old.rx() = 0.0;
        if (overlap > 0) {
    
            double kn = 1.5 * hn_ * sqrt(overlap);
    
            double vn = norm.x() / timestep;
            if (firstActive) {
                if (vn <= 0.0)
                    vni_= vn;
                else
                    vni_ = 0.0;
                fs_old = DVect(0.0);
                ks_old = 0.0;
                hz_F_old = DVect(0.0);
                pfac_ = 0.0;
            }
            hz_F_.rx() = hn_*pow(overlap,hz_alpha_);
    
            double ks = hs_ * pow(hz_F_.x(),(1.0/3.0));
            effectiveTranslationalStiffness_ = DVect2(kn,ks);

            //damping force
            if (std::abs(vni_) <= limits<double>::epsilon()*1000.0)
                dp_F_.rx() = 0.0;
            else {
                double ratio = vn/vni_;
                double fac(0.0);
                switch (dp_mode_) {
                default:
                    {
                        double used = max(dp_en_,dp_enmin_);
                        fac = max(0.0,(1.0-used*used)/used)*ratio; //Gonthier et al.
                    }
                    break;
                case 1: 
                    fac = 0.75*(1.0-dp_en_*dp_en_)*ratio; // Lankarani et al (1989). 
                    break;
                case 2:
                    fac = 0.75*(1.0-dp_en_*dp_en_)*exp(2.0*(1-dp_en_))*ratio; // Zhiying et al.
                    break;
                }

                if (fac <= -1.0) { // sucking
                    if (pfac_ >= 0.0) { //switch in one timestep from pushing to sucking - instability
                        fac = 0.0;
                        pfac_ = 0.0;
                    } else
                        pfac_ = fac;
                } else if (pfac_ != 0.0 && abs(fac) > 10.0*abs(pfac_)) { // growing too fast - instability
                    fac = 0.0;
                    pfac_ = 0.0;
                } else
                    pfac_ = fac;

                dp_F_.rx() = hz_F_.x()*fac;
            }
    
            DVect u_s = trans;
            u_s.rx() = 0.0;
            DVect vec = u_s * ks;
            if (hz_mode_ && (hz_F_.x() < hz_F_old.x())) {
                double rat = ks / ks_old;
                fs_old *= rat;  
            }
            DVect fs = fs_old - vec;
        
            if (state->canFail_) {
                // resolve sliding
                double crit = hz_F_.x() * fric_;
                double sfmag = fs.mag();
                if (sfmag > crit) {
                    double rat = crit / sfmag;
                    fs *= rat;
                    if (!hz_slip_ && cmEvents_[fSlipChange] >= 0) {
                        FArray<QVariant,3> arg;
                        QVariant p1;
                        IContact * c = const_cast<IContact*>(state->getContact());
                        TPtr<IThing> t(c->getIThing());
                        p1.setValue(t);
                        arg.push_back(p1);
                        p1.setValue(0);
                        arg.push_back(p1);
                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
                        fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
                    }
                    hz_slip_ = true;
                } else {
                    if (hz_slip_) {
                        if (cmEvents_[fSlipChange] >= 0) {
                            FArray<QVariant,3> arg;
                            QVariant p1;
                            IContact * c = const_cast<IContact*>(state->getContact());
                            TPtr<IThing> t(c->getIThing());
                            p1.setValue(t);
                            arg.push_back(p1);
                            p1.setValue(1);
                            arg.push_back(p1);
                            IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
                            fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
                        }
                        hz_slip_ = false;
                    }
                }
            }
        
            fs.rx() = hz_F_.x();
            hz_F_ = fs;          // total force in hertz part
            
            // 5) Compute energies
            if (state->trackEnergy_) {
                assert(energies_);
                energies_->estrain_ =  0.0;
                if (kn)
                    energies_->estrain_ = hz_alpha_*hz_F_.x()*hz_F_.x()/((hz_alpha_+1)*kn);
                if (ks) {
                    DVect s = hz_F_;
                    s.rx() = 0.0;
                    double smag2 = s.mag2();
                    energies_->estrain_ += 0.5*smag2 / ks;
                    if (hz_slip_) {
                        hz_F_old.rx() = 0.0;
                        DVect avg_F_s = (s + hz_F_old)*0.5;
                        DVect u_s_el =  (s - hz_F_old) / ks;
                        energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
                    }
                }
                energies_->edashpot_ -= dp_F_ | trans;
            }
        
        } else {
            hz_F_ = DVect(0.0);
            dp_F_ = DVect(0.0);
            pfac_ = 0.0;
        }


        return true;
    }
  
    void ContactModelHysteretic::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
        // Only do something if the contact model is of the same type
        if (old->getContactModel()->getName().compare("hysteretic",Qt::CaseInsensitive) == 0) {
            ContactModelHysteretic *oldCm = (ContactModelHysteretic *)old;
#ifdef THREED
            // Need to rotate just the shear component from oldSystem to newSystem
            
            // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
            DVect axis = oldSystem.e1() & newSystem.e1();
            double c, ang, s;
            DVect re2;
            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
                axis = axis.unit();
                c = oldSystem.e1()|newSystem.e1();
                if (c > 0)
                  c = std::min(c,1.0);
                else
                  c = std::max(c,-1.0);
                ang = acos(c);
                s = sin(ang);
                double t = 1. - c;
                DMatrix<3,3> rm;
                rm.get(0,0) = t*axis.x()*axis.x() + c;
                rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
                rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
                rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
                rm.get(1,1) = t*axis.y()*axis.y() + c;
                rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
                rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
                rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
                rm.get(2,2) = t*axis.z()*axis.z() + c;
                re2 = rm*oldSystem.e2();
            } else
                re2 = oldSystem.e2();
            
            // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
            axis = re2 & newSystem.e2();
            DVect2 tpf;
            DMatrix<2,2> m;
            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
                axis = axis.unit();
                c = re2|newSystem.e2();
                if (c > 0)
                    c = std::min(c,1.0);
                else
                    c = std::max(c,-1.0);
                ang = acos(c);
                if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
                    ang *= -1;
                s = sin(ang);
                m.get(0,0) = c;
                m.get(1,0) = s;
                m.get(0,1) = -m.get(1,0);
                m.get(1,1) = m.get(0,0);
                tpf = m*DVect2(oldCm->hz_F_.y(),oldCm->hz_F_.z());
            } else {
                m.get(0,0) = 1.;
                m.get(0,1) = 0.;
                m.get(1,0) = 0.;
                m.get(1,1) = 1.;
                tpf = DVect2(oldCm->hz_F_.y(),oldCm->hz_F_.z());
            }
            DVect pforce = DVect(0,tpf.x(),tpf.y());
#else
            oldSystem;
            newSystem;
            DVect pforce = DVect(0,oldCm->hz_F_.y());
#endif
            for (int i=1; i<dim; ++i)
                hz_F_.rdof(i) += pforce.dof(i);
            oldCm->hz_F_ = DVect(0.0);

            if(oldCm->getEnergyActivated()) {
                activateEnergy();
                energies_->estrain_  = oldCm->energies_->estrain_;
                energies_->eslip_    = oldCm->energies_->eslip_;
                energies_->edashpot_ = oldCm->energies_->edashpot_;
                oldCm->energies_->estrain_  = 0.0;
                oldCm->energies_->eslip_    = 0.0;
                oldCm->energies_->edashpot_ = 0.0;
            }
        }
    }
  
    void ContactModelHysteretic::setNonForcePropsFrom(IContactModel *old) {
        // Only do something if the contact model is of the same type
        if (old->getName().compare("hysteretic",Qt::CaseInsensitive) == 0 && !isBonded()) {
            ContactModelHysteretic *oldCm = (ContactModelHysteretic *)old;
            hn_ = oldCm->hn_;
            hs_ = oldCm->hs_;
            fric_ = oldCm->fric_;
            vni_ = oldCm->vni_;
            pfac_ = oldCm->pfac_;
          }
    }

    DVect ContactModelHysteretic::getForce(const IContactMechanical *) const {
        DVect ret(hz_F_);
        ret += dp_F_;
        return ret;
    }

    DAVect ContactModelHysteretic::getMomentOn1(const IContactMechanical *c) const {
        DVect force = getForce(c);
        DAVect ret(0.0);
        c->updateResultingTorqueOn1Local(force,&ret);
        return ret;
    }

    DAVect ContactModelHysteretic::getMomentOn2(const IContactMechanical *c) const {
        DVect force = getForce(c);
        DAVect ret(0.0);
        c->updateResultingTorqueOn2Local(force,&ret);
        return ret;
    }
  
} // namespace cmodelsxd
// EoF

Top