Hertz Model Implementation

See this page for the documentation of this contact model.

contactmodelhertz.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
193
194
195
196
197
198
199
200
201
202
203
204
#pragma once
// contactmodelhertz.h

#include "contactmodel/src/contactmodelmechanical.h"

#ifdef HERTZ_LIB
#  define HERTZ_EXPORT EXPORT_TAG
#elif defined(NO_MODEL_IMPORT)
#  define HERTZ_EXPORT
#else
#  define HERTZ_EXPORT IMPORT_TAG
#endif

namespace cmodelsxd {
    using namespace itasca;

    class ContactModelHertz : public ContactModelMechanical {
    public:
        enum PropertyKeys { kwHzShear=1
                          , kwHzPoiss                            
                          , kwFric   
                          , kwHzAlpha
                          , kwHzS
                          , kwHzSd
                          , kwHzF
                          , kwDpNRatio 
                          , kwDpSRatio
                          , kwDpMode 
                          , kwDpF 
                          , kwDpAlpha 
                          , kwRGap
                          };
       
        HERTZ_EXPORT ContactModelHertz();
        HERTZ_EXPORT virtual ~ContactModelHertz();
        virtual void                copy(const ContactModel *c);
        virtual void                archive(ArchiveStream &); 
  
        virtual QString  getName() const { return "hertz"; }
        virtual void     setIndex(int i) { index_=i;}
        virtual int      getIndex() const {return index_;}
      
        virtual QString  getProperties() const { 
            return "hz_shear"
                   ",hz_poiss"
                   ",fric"
                   ",hz_alpha"
                   ",hz_slip"
                   ",hz_mode"
                   ",hz_force"
                   ",dp_nratio"
                   ",dp_sratio"
                   ",dp_mode"
                   ",dp_force"
                   ",dp_alpha"
                   ",rgap"
            ;
        }
  
        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 ContactModelHertz *clone() const { return NEWC(ContactModelHertz()); }
        virtual double              getActivityDistance() const {return rgap_;}
        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 *c);
        virtual void                setArea(const double &) { throw Exception("The setArea method cannot be used with the Hertz contact model."); }
        
        virtual bool     checkActivity(const double &gap) { return gap <= rgap_; }
        
        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 double & hz_alpha() const {return hz_alpha_;}
        void           hz_alpha(const double &d) {hz_alpha_=d;}
        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 & 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 & rgap() const {return rgap_;}
        void           rgap(const double &d) {rgap_=d;}
        
        bool     hasDamping() const {return dpProps_ ? true : false;}
        double   dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
        void     dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
        double   dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
        void     dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
        int      dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
        void     dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
        DVect    dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
        void     dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
        double   dp_alpha() const {return hasDamping() ? dpProps_->dp_alpha_: 0.0;}
        void     dp_alpha(const double &d) { if(!hasDamping()) return; dpProps_->dp_alpha_=d;}
        
        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);
        void   setDampCoefficients(const ContactModelMechanicalState &state,double *vcn,double *vcs);
        // inheritance fields
        quint32 inheritanceField_;
        
        // hertz model
        double      hz_shear_;  // Shear modulus
        double      hz_poiss_;  // Poisson ratio
        double      fric_;      // Coulomb friction coefficient
        double      hz_alpha_;  // Exponent
        bool        hz_slip_;      // the current sliding state
        uint        hz_mode_;     // specifies down-scaling of the shear force when normal unloading occurs 
        DVect       hz_F_;      // Force carried in the hertz model
        double      rgap_;      // Reference gap 
       
        //viscous model
        struct dpProps {
            dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)),dp_alpha_(0.0) {}
            double dp_nratio_;     // normal viscous critical damping ratio
            double dp_sratio_;     // shear  viscous critical damping ratio
            int    dp_mode_;       // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
            DVect  dp_F_;          // Force in the dashpots
            double dp_alpha_;      // exponent
        };
        dpProps *   dpProps_;  
        
        // 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
    };

} // namespace cmodelsxd
// EoF

Top

contactmodelhertz.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
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
// contactmodelhertz.cpp
#include "contactmodelhertz.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 HERTZ_LIB
    int __stdcall DllMain(void *,unsigned, void *) {
        return 1;
    }

    extern "C" EXPORT_TAG const char *getName() {
#if DIM==3
        return "contactmodelmechanical3dhertz";
#else
        return "contactmodelmechanical2dhertz";
#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::ContactModelHertz *m = NEWC(cmodelsxd::ContactModelHertz());
        return (void *)m;
    }
#endif // HERTZ_LIB

namespace cmodelsxd {

    static const quint32 shearMask   = 0x00002;
    static const quint32 poissMask   = 0x00004;
    static const quint32 fricMask    = 0x00008;
  
    using namespace itasca;
  
    int ContactModelHertz::index_ = -1;
    UInt ContactModelHertz::getMinorVersion() const { return MINOR_VERSION;}
  
    ContactModelHertz::ContactModelHertz() : inheritanceField_(shearMask|poissMask|fricMask) 
                                            , hz_shear_(0.0)
                                            , hz_poiss_(0.0)
                                            , fric_(0.0)
                                            , hz_alpha_(1.5)
                                            , hz_slip_(false)
                                            , hz_mode_(0)
                                            , hz_F_(DVect(0.0))
                                            , rgap_(0.0)
                                            , dpProps_(0)
                                            , energies_(0)
                                            , hn_(0.0)
                                            , hs_(0.0)
                                            , effectiveTranslationalStiffness_(DVect2(0.0)) 
    {
    }
  
    ContactModelHertz::~ContactModelHertz() {
        if (dpProps_)
          delete dpProps_;
        if (energies_)
          delete energies_;
    }
  
    void ContactModelHertz::archive(ArchiveStream &stream) {
        stream & hz_shear_;
        stream & hz_poiss_;
        stream & fric_;
        stream & hz_alpha_;
        stream & hz_slip_;
        stream & hz_mode_;
        stream & hz_F_;
        stream & hn_;
        stream & hs_;

        if (stream.getArchiveState()==ArchiveStream::Save) {
            bool b = false;
            if (dpProps_) {
                b = true;
                stream & b;
                stream & dpProps_->dp_nratio_; 
                stream & dpProps_->dp_sratio_; 
                stream & dpProps_->dp_mode_; 
                stream & dpProps_->dp_F_; 
                stream & dpProps_->dp_alpha_; 
            } else
                stream & b;
            
            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 (!dpProps_)
                    dpProps_ = NEWC(dpProps());
                stream & dpProps_->dp_nratio_; 
                stream & dpProps_->dp_sratio_; 
                stream & dpProps_->dp_mode_; 
                stream & dpProps_->dp_F_; 
                if (stream.getRestoreVersion() >= 2)
                    stream & dpProps_->dp_alpha_;
            }
            stream & b;
            if (b) {
                if (!energies_)
                    energies_ = NEWC(Energies());
                stream & energies_->estrain_;
                stream & energies_->eslip_;
                stream & energies_->edashpot_;
            }
        }
    
        stream & inheritanceField_;
        stream & effectiveTranslationalStiffness_;
        if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() >= 2)
            stream & rgap_;

    }
  
    void ContactModelHertz::copy(const ContactModel *cm) {
        ContactModelMechanical::copy(cm);
        const ContactModelHertz *in = dynamic_cast<const ContactModelHertz*>(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_alpha(in->hz_alpha());
        hz_S(in->hz_S());
        hz_mode(in->hz_mode());
        hz_F(in->hz_F());
        hn(in->hn());
        hs(in->hs());
        rgap(in->rgap());
        if (in->hasDamping()) {
            if (!dpProps_)
                dpProps_ = NEWC(dpProps());
            dp_nratio(in->dp_nratio()); 
            dp_sratio(in->dp_sratio()); 
            dp_mode(in->dp_mode()); 
            dp_F(in->dp_F()); 
            dp_alpha(in->dp_alpha()); 
        }
        if (in->hasEnergies()) {
            if (!energies_)
                energies_ = NEWC(Energies());
            estrain(in->estrain());
            eslip(in->eslip());
            edashpot(in->edashpot());
        }
        inheritanceField(in->inheritanceField());
        effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
    }
  
    QVariant ContactModelHertz::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 kwHzAlpha:   return hz_alpha_;
          case kwHzS:       return hz_slip_;
          case kwHzSd:      return hz_mode_;
          case kwHzF:       var.setValue(hz_F_); return var;
          case kwRGap:      return rgap_;
          case kwDpNRatio:  return dpProps_ ? dpProps_->dp_nratio_ : 0.0;
          case kwDpSRatio:  return dpProps_ ? dpProps_->dp_sratio_ : 0.0;
          case kwDpMode:    return dpProps_ ? dpProps_->dp_mode_   : 0;
          case kwDpAlpha:   return dpProps_ ? dpProps_->dp_alpha_  : 0.0;
          case kwDpF:{
              dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
              return var;
            }
      }
      assert(0);
      return QVariant();
    }
  
    bool ContactModelHertz::getPropertyGlobal(uint i) const {
        switch (i) {
            case kwHzF: // fall through   
            case kwDpF: return false;
        }
        return true;
    }
  
    bool ContactModelHertz::setProperty(uint i,const QVariant &v,IContact *) {
        dpProps dp;
        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 kwHzAlpha: {
                if (!v.canConvert<double>())
                    throw Exception("hz_alpha must be a double.");
                double val(v.toDouble());
                if (val<=0.0)
                    throw Exception("Negative exponent value not allowed.");
                hz_alpha_ = 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 kwRGap: {
                if (!v.canConvert<double>())
                    throw Exception("Reference gap must be a double.");
                double val(v.toDouble());
                rgap_ = val;  
                return false;
            }
            case kwDpNRatio: {
                if (!v.canConvert<double>())
                    throw Exception("dp_nratio must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative dp_nratio not allowed.");
                if (val == 0.0 && !dpProps_)
                    return false;
                if (!dpProps_)
                    dpProps_ = NEWC(dpProps());
                dpProps_->dp_nratio_ = val; 
                return true;
            }
            case kwDpSRatio: {
                if (!v.canConvert<double>())
                    throw Exception("dp_sratio must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative dp_sratio not allowed.");
                if (val == 0.0 && !dpProps_)
                    return false;
                if (!dpProps_)
                    dpProps_ = NEWC(dpProps());
                dpProps_->dp_sratio_ = val;
                return true;
            }
            case kwDpMode: {
                if (!v.canConvert<int>())
                    throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
                int val(v.toInt());
                if (val == 0 && !dpProps_)
                    return false;
                if (val < 0 || val > 3)
                    throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
                if (!dpProps_)
                    dpProps_ = NEWC(dpProps());
                dpProps_->dp_mode_ = val;
                return false;
            }
            case kwDpAlpha: {
                if (!v.canConvert<double>())
                    throw Exception("dp_alpha must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative dp_alpha not allowed.");
                if (val == 0.0 && !dpProps_)
                    return false;
                if (!dpProps_)
                    dpProps_ = NEWC(dpProps());
                dpProps_->dp_alpha_ = val; 
                return true;
            }
            case kwHzF: {
                if (!v.canConvert<DVect>())
                    throw Exception("hz_force must be a vector.");
                DVect val(v.value<DVect>());
                hz_F_ = val;
                return false;
            }
        }
        return false;
    }
  
    bool ContactModelHertz::getPropertyReadOnly(uint i) const {
        switch (i) {
//            case kwHzF:
            case kwDpF:
            case kwHzS:
                return true;
            default:
                break;
        }
        return false;
    }
  
    bool ContactModelHertz::supportsInheritance(uint i) const {
        switch (i) {
            case kwHzShear:
            case kwHzPoiss:
            case kwFric:
                return true;
            default:
                break;
        }
        return false;
    }
  
    double ContactModelHertz::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 ContactModelHertz::getEnergyAccumulate(uint i) const {
        switch (i) {
            case kwEStrain:  return false;
            case kwESlip:    return true;
            case kwEDashpot: return true;
        }
        assert(0);
        return false;
    }
  
    void ContactModelHertz::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 ContactModelHertz::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 ContactModelHertz::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("Hertz 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*(1-hz_poiss_)/(2.0- hz_poiss_))*hz_alpha_*pow(hn_,1.0/hz_alpha_);
        return ( (hn_ != hnold) || (hs_ != hsold) );
    }


    static const QString gstr("hz_shear");
    static const QString nustr("hz_poiss");
    bool ContactModelHertz::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 Hertz 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 Hertz 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 Hertz contact model");
        if (hz_poiss_ <= -1.0 || hz_poiss_ > 0.5) 
            throw Exception("Poisson ratio should be in range (-1.0,0.5] in Hertz contact model");
        return updateStiffCoef(con);
    }
  
    static const QString fricstr("fric");
    bool ContactModelHertz::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 ContactModelHertz::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 ContactModelHertz::updateEffectiveStiffness(ContactModelMechanicalState *state) {
        effectiveTranslationalStiffness_ = DVect2(hn_,hs_);
        double overlap = rgap_ - state->gap_;
        if (overlap <= 0.0) return;
        double kn = hz_alpha_*hn_*pow(overlap,hz_alpha_-1.0);
        double ks = hs_ * pow(hz_F_.x(),(hz_alpha_-1.0)/hz_alpha_);
        DVect2 ret(kn,ks);
        // correction if viscous damping active
        if (dpProps_) {
            DVect2 correct(1.0);
            if (dpProps_->dp_nratio_)
                correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
            if (dpProps_->dp_sratio_)
                correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
            ret /= (correct*correct);
        }
        effectiveTranslationalStiffness_ = ret;
    }
    
    bool ContactModelHertz::forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) {
        assert(state);

        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]);
            }
        }

        double overlap = rgap_ - 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_;
        // normal force in Hertz part
        double fn = hn_*pow(overlap,hz_alpha_);
        // tangent normal stiffness
        double kn = hz_alpha_ * hn_ * pow(overlap,hz_alpha_-1.0);
        // initial tangent shear stiffness 
        double ks = hs_ * pow(fn,(hz_alpha_- 1.0)/hz_alpha_);

        DVect fs_old = hz_F_;
        fs_old.rx() = 0.0;

        if (hz_mode_ && fn < hz_F_.x()) {
            double ks_old = hs_ * pow(hz_F_.x(),(hz_alpha_- 1.0)/hz_alpha_);
            double rat = ks/ks_old;
            fs_old *= rat; 
        }
        
        DVect u_s = trans;
        u_s.rx() = 0.0;
        DVect vec = u_s * ks;

        DVect fs = fs_old - vec;
    
        if (state->canFail_) {
            // resolve sliding
            double crit = fn * 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;
                }
            }
        }
    
        hz_F_ = fs ;          // total force in hertz part
        hz_F_.rx() += fn;
        effectiveTranslationalStiffness_ = DVect2(kn,ks);
        // 3) Account for dashpot forces
        if (dpProps_) {
            dpProps_->dp_F_.fill(0.0);
            double vcn(0.0), vcs(0.0);
            setDampCoefficients(*state,&vcn,&vcs);
            double fac = 1.0;
            if (dpProps_->dp_alpha_ > 0.0) fac = pow(overlap,dpProps_->dp_alpha_);
            // First damp all components
            dpProps_->dp_F_  = u_s * (-1.0* vcs*fac) / timestep; // shear component   
            dpProps_->dp_F_ -= norm * vcn*fac / timestep;       // normal component
            // Need to change behavior based on the dp_mode
            if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
                // limit the tensile if not bonded
                if (dpProps_->dp_F_.x() + hz_F_.x() < 0)
                    dpProps_->dp_F_.rx() = - hz_F_.rx();
            }
            if (hz_slip_ && dpProps_->dp_mode_ > 1) {
                // limit the shear if not sliding
                double dfn = dpProps_->dp_F_.rx();
                dpProps_->dp_F_.fill(0.0); 
                dpProps_->dp_F_.rx() = dfn; 
            }
            // Correct effective translational stiffness
            DVect2 correct(1.0);
            if (dpProps_->dp_nratio_)
                correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
            if (dpProps_->dp_sratio_)
                correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
            effectiveTranslationalStiffness_ /= (correct*correct);       
        }
    
        // 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.0)*kn);
            if (ks) {
                double smag2 = fs.mag2();
                energies_->estrain_ += 0.5*smag2 / ks;
                
                if (hz_slip_) {
                    DVect avg_F_s = (fs + fs_old)*0.5;
                    DVect u_s_el =  (fs - fs_old) / ks;
                    energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
                }
            }
            if (dpProps_) {
                energies_->edashpot_ -= dpProps_->dp_F_ | trans;
            }
        }
    
        return true;
    }

    void ContactModelHertz::setForce(const DVect &v,IContact *c) { 
        hz_F(v); 
        if (v.x() > 0) 
            rgap_ = c->getGap() + pow(v.x()/hn_,1./hz_alpha_); 
    }
  
    void ContactModelHertz::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("hertz",Qt::CaseInsensitive) == 0) {
            ContactModelHertz *oldCm = (ContactModelHertz *)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 (dpProps_ && oldCm->dpProps_) {
#ifdef THREED
                tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
                pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
#else
                pforce = oldCm->dpProps_->dp_F_;
#endif
                dpProps_->dp_F_ += pforce;
                oldCm->dpProps_->dp_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;
            }
            rgap_ = oldCm->rgap_;
        }
    }
  
    void ContactModelHertz::setNonForcePropsFrom(IContactModel *old) {
        // Only do something if the contact model is of the same type
        if (old->getName().compare("hertz",Qt::CaseInsensitive) == 0 && !isBonded()) {
            ContactModelHertz *oldCm = (ContactModelHertz *)old;
            hn_ = oldCm->hn_;
            hs_ = oldCm->hs_;
            fric_ = oldCm->fric_;
            rgap_ = oldCm->rgap_;
  
            if (oldCm->dpProps_) {
                if (!dpProps_)
                    dpProps_ = NEWC(dpProps());
                dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
                dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
                dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
            }
        }
    }

    DVect ContactModelHertz::getForce(const IContactMechanical *) const {
        DVect ret(hz_F_);
        if (dpProps_)
            ret += dpProps_->dp_F_;
        return ret;
    }

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

    DAVect ContactModelHertz::getMomentOn2(const IContactMechanical *c) const {
        DVect force = getForce(c);
        DAVect ret(0.0);
        c->updateResultingTorqueOn2Local(force,&ret);
        return ret;
    }
  
    void ContactModelHertz::setDampCoefficients(const ContactModelMechanicalState &state,double *vcn,double *vcs) {
        double overlap = rgap_ - state.gap_;
        double kn = hz_alpha_*hn_*pow(overlap,hz_alpha_-1.0);
        double ks = hs_ * pow(hz_F_.x(),(hz_alpha_-1.0)/hz_alpha_);
        *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(state.inertialMass_*(kn));
        *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(state.inertialMass_*(ks));
    }

} // namespace cmodelsxd
// EoF

Top