Smooth-Joint Model Implementation

See this page for the documentation of this contact model.

contactmodelsmoothjoint.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
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
#pragma once
// contactmodelsmoothjoint.h

#include "contactmodel/src/contactmodelmechanical.h"

#ifdef SMOOTHJOINT_LIB
#  define SMOOTHJOINT_EXPORT EXPORT_TAG
#elif defined(NO_MODEL_IMPORT)
#  define SMOOTHJOINT_EXPORT
#else
#  define SMOOTHJOINT_EXPORT IMPORT_TAG
#endif

namespace cmodelsxd {
    using namespace itasca;

    class ContactModelSmoothJoint : public ContactModelMechanical {
    public:
        enum PropertyKeys { 
               kwKn=1    
             , kwKs
             , kwFric   
             , kwDA
             , kwState
             , kwTen    
             , kwBCoh    
             , kwBFA
             , kwShear
             , kwRMul   
             , kwLarge
             , kwFn
             , kwFs     
             , kwGap
             , kwNorm
             , kwArea 
             , kwRad
             , kwUn     
             , kwUs
             , kwSlip
             , kwDip    
             , kwDD    
        };

        enum FishCallEvents { fBondBreak=0, fSlipChange };

        SMOOTHJOINT_EXPORT ContactModelSmoothJoint();
        SMOOTHJOINT_EXPORT virtual ~ContactModelSmoothJoint();
        virtual void     copy(const ContactModel *c);
        virtual void     archive(ArchiveStream &); 
        virtual QString  getName() const { return "smoothjoint"; }
        virtual void     setIndex(int i) { index_=i;}
        virtual int      getIndex() const {return index_;}

        virtual QString  getProperties() const { 
            return "sj_kn"
                   ",sj_ks"
                   ",sj_fric"
                   ",sj_da"
                   ",sj_state"
                   ",sj_ten"
                   ",sj_coh"
                   ",sj_fa"
                   ",sj_shear"
                   ",sj_rmul"
                   ",sj_large"
                   ",sj_fn"
                   ",sj_fs"
                   ",sj_gap"
                   ",sj_unorm"
                   ",sj_area"
                   ",sj_radius"
                   ",sj_un"
                   ",sj_us"
                   ",sj_slip"
                   ",dip"
#ifdef THREED
                   ",ddir"
#endif
                   ; 
        }

        virtual QString  getFishCallEvents() const { return "bond_break,slip_change"; }
        virtual QVariant getProperty(uint i,const IContact *con=0) const;
        virtual bool     getPropertyGlobal(uint i) const;
        virtual bool     setProperty(uint i,const QVariant &v,IContact *con=0);
        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; }

        enum MethodKeys { 
              kwSetForce=1,kwBond,kwUnbond 
        };

        virtual QString  getMethods() const { 
            return  "sj_setforce,bond,unbond";
        }

        virtual QString  getMethodArguments(uint i) const; 
        virtual bool     setMethod(uint i,const QVector<QVariant> &vl,IContact *con=0); // Base 1 - returns true if timestep contributions need to be updated

        virtual uint     getMinorVersion() const;

        enum EnergyKeys { kwEStrain=1,kwESlip};
        virtual QString  getEnergies() const { return "energy-strain,energy-slip";}
        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);}

        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 ContactModelSmoothJoint *clone() const { return NEWC(ContactModelSmoothJoint()); }
        virtual double  getActivityDistance() const {return 0.0;}
        virtual void    resetForcesAndMoments() { sj_Fn(0.0); sj_Fs(DVect(0.0)); }
        virtual void    setForce(const DVect &v,IContact *);
        virtual void    setArea(const double &) { throw Exception("The setArea method cannot be used with this contact model."); }

        virtual bool    isOKToDelete() const { return (!isBonded() && sj_large_); }

        virtual bool    checkActivity(const double &gap) { return ( !sj_large_ || gap <= 0.0 || isBonded()); }
        virtual bool    isSliding() const { return isjSliding_; }
        virtual bool    isBonded() const { return sj_state_ == 3; }
        virtual bool    hasNormal() const { return true; }
        virtual DVect   getNormal() const { return sj_unorm_; }

        void sj_kn(const double &d)    {sj_kn_ = d;}
        void sj_ks(const double &d)    {sj_ks_ = d;}
        void sj_fric(const double &d)  {sj_fric_ = d;}
        void sj_da(const double &d)    {sj_da_ = d;}
        void sj_state(const double &d) {sj_state_ = d;}
        void sj_bns(const double &d)   {sj_bns_ = d;}
        void sj_bcoh(const double &d)  {sj_bcoh_ = d;}
        void sj_bfa(const double &d)   {sj_bfa_ = d;}
        void dip(const double &d)      {dip_ = d;}
        void sj_rmul(const double &d)  {sj_rmul_ = d;}
        void sj_large(bool b)          {sj_large_ = b;}

        const double & sj_kn()    const  {return sj_kn_;}    
        const double & sj_ks()    const  {return sj_ks_;}    
        const double & sj_fric()  const  {return sj_fric_;}  
        const double & sj_da()    const  {return sj_da_;}    
        int            sj_state() const  {return sj_state_;} 
        const double & sj_bns()   const  {return sj_bns_;}   
        const double & sj_bcoh()  const  {return sj_bcoh_;}  
        const double & sj_bfa()   const  {return sj_bfa_;}   
        const double & dip()      const  {return dip_;}      
        const double & sj_rmul()  const  {return sj_rmul_;}  
        bool           sj_large() const  {return sj_large_;} 

#ifdef THREED
        const double & dd() const          {return dd_;}
        void           dd(const double &d) {dd_ =d;}
#endif

        void sj_unorm(const DVect &v)  {sj_unorm_ = v;}
        void sj_A(const double &d)     {sj_A_ = d;}
        void sj_rad(const double &d)   {sj_rad_ = d;}
        void sj_gap(const double &d)   {sj_gap_ = d;}
        void sj_Un(const double &d)    {sj_Un_ = d;}
        void sj_Us(const DVect &v)     {sj_Us_ = v;}
        void sj_Fn(const double &d)    {sj_Fn_ = d;}
        void sj_Fs(const DVect &v)     {sj_Fs_ = v;}
        void isjFlip(bool b)           {isjFlip_ = b;}
        void isjSliding(bool b)        {isjSliding_ = b;}

        const DVect  & sj_unorm()   const {return sj_unorm_;}  
        const double & sj_A()       const {return sj_A_;}      
        const double & sj_rad()     const {return sj_rad_;}    
        const double & sj_gap()     const {return sj_gap_;}   
        const double & sj_Un()      const {return sj_Un_;}     
        const DVect  & sj_Us()      const {return sj_Us_;}     
        const double & sj_Fn()      const {return sj_Fn_;}     
        const DVect  & sj_Fs()      const {return sj_Fs_;}     
        bool           isjFlip()    const {return isjFlip_;}   
        bool           isjSliding() const {return isjSliding_;}

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

        uint inheritanceField() const {return inheritanceField_;}
        void inheritanceField(uint i) {inheritanceField_ = i;}

        const DVect2 & effectiveTranslationalStiffness()  const          {return effectiveTranslationalStiffness_;}
        void           effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}

        bool    geomRecomp()       const {return geomRecomp_ ;}
        void    geomRecomp(bool b)       {geomRecomp_ = b;}

        /// 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_;

        struct Energies {
            Energies() : estrain_(0.0), eslip_(0.0) {}
            double estrain_;  // elastic energy stored in contact 
            double eslip_;    // work dissipated by friction 
        };

        bool   updateKn(const IContactMechanical *con);
        bool   updateKs(const IContactMechanical *con);
        bool   updateFric(const IContactMechanical *con);
        void   updateAreaAndNormal(ContactModelMechanicalState *state);
        double calcBSS() const;

        void   updateEffectiveTranslationalStiffness();

        // property set fields
        quint32 inheritanceField_;

        // smooth joint model - contact model properties
        double sj_kn_;     // normal stiffness
        double sj_ks_;     // shear stiffness
        double sj_fric_;   // Coulomb friction coefficient
        double sj_da_;     // Dilation angle (stored in radians, returned/set in degrees)
        int    sj_state_;  // Bond state - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (bonded)
        double sj_bns_;    // Bond normal (tensile) strength
        double sj_bcoh_;   // Bonded system cohesion
        double sj_bfa_;    // Bonded system friction angle (stored in radians, returned/set in degrees)
        double dip_;       // Dip angle (stored in radians, returned/set in degrees)
        double sj_rmul_;   // Radius multiplier
        bool   sj_large_;  // Large strain indicator

        // Internal state variables
        DVect  sj_unorm_;  // Normal to the plane
        double sj_A_;      // Cross-sectional area
        double sj_rad_;    // Radius
        double sj_gap_;    // Gap - this can be user modified
        double sj_Un_;     // Normal displacement
        DVect  sj_Us_;     // Shear displacement
        double sj_Fn_;     // Normal force
        DVect  sj_Fs_;     // Shear force
        bool   isjFlip_;   // In order to flip the order
        bool   isjSliding_;// True if this is sliding

        DVect2  effectiveTranslationalStiffness_;
        bool    geomRecomp_; // If true then there must be a validate before FD!
#ifdef THREED
        double dd_;        // Dip direction (stored in radians, returned/set in degrees)
#endif
        //DAVect  effectiveRotationalStiffness_;

        Energies *energies_;
    };
} // namespace itascaxd
// EoF

Top

contactmodelsmoothjoint.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
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
// contactmodelsmoothjoint.cpp
#include "contactmodelsmoothjoint.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 "base/src/basetoqt.h"

#include "module/interface/ifishcalllist.h"
#include "kernel/interface/iprogram.h"
#include "utility/src/tptr.h"

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

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

namespace cmodelsxd {
    static const quint32 sjKnMask      = 0x00002; // Base 1!
    static const quint32 sjKsMask      = 0x00004;
    static const quint32 sjFricMask    = 0x00008;

    using namespace itasca;

    int ContactModelSmoothJoint::index_ = -1;
    UInt ContactModelSmoothJoint::getMinorVersion() const { return MINOR_VERSION;}

    ContactModelSmoothJoint::ContactModelSmoothJoint() : inheritanceField_(sjKnMask|sjKsMask|sjFricMask)
                                                       , sj_kn_(0.0), sj_ks_(0.0), sj_fric_(0.0)
                                                       , sj_da_(0.0), sj_state_(0), sj_bns_(0.0)
                                                       , sj_bcoh_(0.0), sj_bfa_(0.0), dip_(0.0)
                                                       , sj_rmul_(1.0), sj_large_(false)  
                                                       , sj_unorm_(0.0), sj_A_(0.0), sj_rad_(0.0)
                                                       , sj_gap_(0.0), sj_Un_(0.0), sj_Us_(0.0)
                                                       , sj_Fn_(0.0), sj_Fs_(0.0), isjFlip_(false)
                                                       , isjSliding_(false)
                                                       , effectiveTranslationalStiffness_(0.0)
                                                       , geomRecomp_(true)
#ifdef THREED
                                                       , dd_(0.0)
#endif
                                                       , energies_(0) {
//    setFromParent(ContactModelMechanicalList::instance()->find(getName()));
    }

    ContactModelSmoothJoint::~ContactModelSmoothJoint() {
        if (energies_)
            delete energies_;
    }

    void ContactModelSmoothJoint::archive(ArchiveStream &stream) {
        stream & sj_kn_;
        stream & sj_ks_;
        stream & sj_fric_;
        stream & sj_da_;
        stream & sj_state_;
        stream & sj_bns_;  
        stream & sj_bcoh_; 
        stream & sj_bfa_;  
        stream & dip_;     
        stream & sj_rmul_; 
        stream & sj_large_;
        stream & sj_unorm_;  
        stream & sj_A_;      
        stream & sj_rad_;    
        stream & sj_gap_;    
        stream & sj_Un_;     
        stream & sj_Us_;     
        stream & sj_Fn_;     
        stream & sj_Fs_;     
        stream & isjFlip_;   
        stream & isjSliding_;
#ifdef THREED
        stream & dd_;
#endif
        if (stream.getArchiveState()==ArchiveStream::Save) {
            bool b = false;
            if (energies_) {
                b = true;
                stream & b;
                stream & energies_->estrain_;
                stream & energies_->eslip_;
            }
            else
                stream & b;
        } else {
            bool b(false);
            stream & b;
            if (b) {
                if (!energies_)
                    energies_ = NEWC(Energies());
                stream & energies_->estrain_;
                stream & energies_->eslip_;
            }
        }

        stream & inheritanceField_;
        stream & geomRecomp_;
        stream & effectiveTranslationalStiffness_;

    }

    void ContactModelSmoothJoint::copy(const ContactModel *cm) {
        ContactModelMechanical::copy(cm);
        const ContactModelSmoothJoint *in = dynamic_cast<const ContactModelSmoothJoint*>(cm);
        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
        sj_kn(in->sj_kn());
        sj_ks(in->sj_ks());
        sj_fric(in->sj_fric());
        sj_da(in->sj_da());
        sj_state(in->sj_state());
        sj_bns(in->sj_bns());
        sj_bcoh(in->sj_bcoh());
        sj_bfa(in->sj_bfa());
        dip(in->dip());
        sj_rmul(in->sj_rmul()); 
        sj_large(in->sj_large());
#ifdef THREED
        dd(in->dd());
#endif
        sj_unorm(in->sj_unorm());  
        sj_A(in->sj_A());      
        sj_rad(in->sj_rad());    
        sj_gap(in->sj_gap());    
        sj_Un(in->sj_Un());     
        sj_Us(in->sj_Us());     
        sj_Fn(in->sj_Fn());     
        sj_Fs(in->sj_Fs());     
        isjFlip(in->isjFlip());   
        isjSliding(in->isjSliding());

        if (in->hasEnergies()) {
            if (!energies_)
                energies_ = NEWC(Energies());
            estrain(in->estrain());
            eslip(in->eslip());
        }
        inheritanceField(in->inheritanceField());
        effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
        geomRecomp(in->geomRecomp());
    }

    QVariant ContactModelSmoothJoint::getProperty(uint i,const IContact *) const {
        QVariant var;
        switch (i) {
        case kwKn:      return sj_kn_;
        case kwKs:      return sj_ks_;
        case kwFric:    return sj_fric_;
        case kwDA:      return sj_da_/dDegrad;  // Returned in degrees
        case kwState:   return sj_state_;
        case kwTen:     return sj_bns_;
        case kwBCoh:    return sj_bcoh_;
        case kwBFA:     return sj_bfa_/dDegrad; // Returned in degrees
        case kwShear:   return calcBSS();
        case kwDip:     return dip_/dDegrad; // Returned in degrees
#ifdef THREED
        case kwDD:      return dd_/dDegrad;  // Returned in degrees
#endif
        case kwRMul:    return sj_rmul_;
        case kwLarge:   return sj_large_;
        case kwFn:      return sj_Fn_;
        case kwFs: {
                var.setValue(sj_Fs_);
                return var;
            }
        case kwGap:     return sj_gap_;
        case kwNorm: {
                var.setValue(sj_unorm_);
                return var;
            }
        case kwArea:    return sj_A_;
        case kwRad:     return sj_rad_;
        case kwUn:      return sj_Un_;
        case kwUs: {
                var.setValue(sj_Us_);
                return var;
            }
        case kwSlip: return isjSliding_;
        }
        assert(0);
        return QVariant();
    }

    bool ContactModelSmoothJoint::getPropertyGlobal(uint ) const {
        return true;
    }

    bool ContactModelSmoothJoint::setProperty(uint i,const QVariant &v,IContact *) {
        switch (i) {
        case kwKn: {
                if (!v.canConvert<double>())
                    throw Exception("sj_kn must be a double.");
                double val = v.toDouble();
                if (val<0.0)
                    throw Exception("Negative sj_kn not allowed.");
                sj_kn_ = val;
                updateEffectiveTranslationalStiffness();
                return true;
            }
        case kwKs: {
                if (!v.canConvert<double>())
                    throw Exception("sj_ks must be a double.");
                double val = v.toDouble();
                if (val<0.0)
                    throw Exception("Negative sj_ks not allowed.");
                sj_ks_ = val;  
                updateEffectiveTranslationalStiffness();
                return true;
            }
        case kwFric: {
                if (!v.canConvert<double>())
                    throw Exception("sj_fric must be a double.");
                double val = v.toDouble();
                if (val<0.0)
                    throw Exception("Negative friction not allowed.");
                sj_fric_ = val;  
                return false;
            }
        case kwDA: {
                if (!v.canConvert<double>())
                    throw Exception("sj_da must be a double.");
                sj_da_ = v.toDouble()*dDegrad; // Given in degrees
                return false;
            }
        case kwState: {
                if (!v.canConvert<uint>())
                    throw Exception("sj_state must be must be an integer between 0 and 3.");
                int val = v.toInt();
                if (val<0 || val>3)
                    throw Exception("The bond state must be an integer between 0 and 3.");
                sj_state_ = val;  
                return false;
            }
        case kwTen: {
                if (!v.canConvert<double>())
                    throw Exception("sj_ten must be a double.");
                double val = v.toDouble();
                if (val<0.0)
                    throw Exception("Negative bond normal (tensile) strength not allowed.");
                sj_bns_ = val;  
                return false;
            }
        case kwBCoh: {
                if (!v.canConvert<double>())
                    throw Exception("sj_coh must be a double.");
                double val = v.toDouble();
                if (val<0.0)
                    throw Exception("Negative bond system cohesion not allowed.");
                sj_bcoh_ = val;  
                return false;
            }
        case kwBFA: {
                if (!v.canConvert<double>())
                    throw Exception("sj_bfa must be a double.");
                sj_bfa_ = v.toDouble()*dDegrad; // Given in degrees
                return false;
            }
        case kwDip: {
                if (!v.canConvert<double>())
                    throw Exception("dip must be a double.");
                dip_ = v.toDouble()*dDegrad; // Given in degrees
                if (!geomRecomp_) geomRecomp_ = true;
                return true;
            }
        case kwDD: {
#ifdef THREED
                if (!v.canConvert<double>())
                    throw Exception("ddir must be a double.");
                dd_ = v.toDouble()*dDegrad; // Given in degrees
                if (!geomRecomp_) geomRecomp_ = true;
#endif
                return true;
            }
        case kwRMul: {
                if (!v.canConvert<double>())
                    throw Exception("sj_rmul must be a double.");
                double val = v.toDouble();
                if (val<0.0)
                    throw Exception("Negative radius multiplier not allowed.");
                sj_rmul_ = val;  
                if (!geomRecomp_) geomRecomp_ = true;
                return false;
            }
        case kwLarge: {
                if (!v.canConvert<bool>())
                    throw Exception("Large-strain flag must be a boolean.");
                sj_large_ = v.toBool();  
                return false;
            }
        case kwFn: {
                 if (!v.canConvert<double>())
                    throw Exception("sj_fn must be a double.");
               sj_Fn_ = v.toDouble();
                return false;
            }
        case kwFs: {
                if (!v.canConvert<DVect>())
                    throw Exception("sj_fs must be a vector.");
                sj_Fs_ = v.value<DVect>();
                return false;
            }
        case kwGap: {
                if (!v.canConvert<double>())
                    throw Exception("sj_gap must be a double.");
                sj_gap_ = v.toDouble();
                return false;
            }

         
        // Following are read-only !
        case kwNorm:
        case kwRad:
        case kwShear:
        case kwUn:
        case kwUs:
        case kwSlip:
            return false;
        }
        assert(0);
        return false;
    }

    bool ContactModelSmoothJoint::getPropertyReadOnly(uint i) const {
        switch (i) {
        case kwNorm:
        case kwRad:
        case kwUn:
        case kwUs:
        case kwSlip:
        case kwShear:
        case kwArea:
            return true;
        default:
            break;
        }
        return false;
    }

    bool ContactModelSmoothJoint::supportsInheritance(uint i) const {
        switch (i) {
        case kwKn:
        case kwKs:
        case kwFric:
            return true;
        default:
            break;
        }
        return false;
    }

    QString  ContactModelSmoothJoint::getMethodArguments(uint i) const {
        QString def = QString();
        switch (i) {
        case kwBond: 
            return "gap";
        case kwUnbond: 
            return "gap";
        }
        return def;
    }

    bool ContactModelSmoothJoint::setMethod(uint i,const QVector<QVariant> &vl,IContact *con) {
        IContactMechanical *c(convert_getcast<IContactMechanical>(con));
        switch (i) {
        case kwSetForce: {
                DVect gf = c->getGlobalForce();
                if (isjFlip_)
                    gf *= -1.0;
                sj_Fn_ = gf | sj_unorm_;
                return false;
            }
        case kwBond: {
                if (sj_state_ == 3) return false;
                double mingap = -1.0 * limits<double>::max();
                double maxgap = 0;
                if (vl.at(0).canConvert<Double>()) 
                    maxgap = vl.at(0).toDouble();
                else if (vl.at(0).canConvert<DVect2>()) {
                    DVect2 value = vl.at(0).value<DVect2>();
                    mingap = value.minComp();
                    maxgap = value.maxComp();
                } else if (!vl.at(0).isNull())
                    throw Exception("gap value %1 not recognized in method bond in contact model %2.",vl.at(0),getName());
                if ( sj_gap_ >= mingap && sj_gap_ <= maxgap) {
                    sj_state_ = 3;
                    return true;
                }
                return false;
        }
        case kwUnbond: {
                if (sj_state_ == 0) return false;
                double mingap = -1.0 * limits<double>::max();
                double maxgap = 0;
                if (vl.at(0).canConvert<double>()) 
                    maxgap = vl.at(0).toDouble();
                else if (vl.at(0).canConvert<DVect2>()) {
                    DVect2 value = vl.at(0).value<DVect2>();
                    mingap = value.minComp();
                    maxgap = value.maxComp();
                }
                else if (!vl.at(0).isNull())
                    throw Exception("gap value %1 not recognized in method unbond in contact model %2.",vl.at(0),getName());
                if ( sj_gap_ >= mingap && sj_gap_ <= maxgap) {
                    sj_state_ = 0;
                    return true;
                }
                return false;
            }

        }
        return false;
    }

    double ContactModelSmoothJoint::getEnergy(uint i) const {
        double ret(0.0);
        if (!energies_)
            return ret;
        switch (i) {
        case kwEStrain:  return energies_->estrain_;
        case kwESlip:    return energies_->eslip_;
        }
        assert(0);
        return ret;
    }

    bool ContactModelSmoothJoint::getEnergyAccumulate(uint i) const {
        bool ret(false);
        if (!energies_) return ret;
        switch (i) {
        case kwEStrain:  return false;
        case kwESlip:    return true;
        }
        assert(0);
        return ret;
    }

    void ContactModelSmoothJoint::setEnergy(uint i,const double &d) {
        if (!energies_) return;
        switch (i) {
        case kwEStrain:  energies_->estrain_ = d; return;  
        case kwESlip:    energies_->eslip_   = d; return;
        }
        assert(0);
        return;
    }

    bool ContactModelSmoothJoint::validate(ContactModelMechanicalState *state,const double &) {
        assert(state);
        const IContactMechanical *c = state->getMechanicalContact(); 
        assert(c);

        if (state->trackEnergy_)
            activateEnergy();

        // The radius multiplier has been set previously so calculate the sj_A_
        // Need to do this regardless of whether or not the radius multiplier has been updated as this is the 
        // first place with the contact state to get the ball radii!
        updateAreaAndNormal(state);

        if (inheritanceField_ & sjKnMask)    
            updateKn(c);
        if (inheritanceField_ & sjKsMask)      
            updateKs(c);
        if (inheritanceField_ & sjFricMask)   
            updateFric(c);

        updateEffectiveTranslationalStiffness();

       return checkActivity(state->gap_);
    }

    void ContactModelSmoothJoint::updateAreaAndNormal(ContactModelMechanicalState *state) {
        assert(state);
        if (geomRecomp_) geomRecomp_ = false;
        // Use the maximum value of the curvature here - not the min!
        sj_rad_ = 1.0 / std::max(state->end1Curvature_.y(),state->end2Curvature_.y());
        sj_rad_ = sj_rmul_ * sj_rad_;
#ifdef THREED
        sj_A_ = dPi * sj_rad_ * sj_rad_;
        sj_unorm_.rx() = sin(dip_) * sin(dd_);
        sj_unorm_.ry() = sin(dip_) * cos(dd_);
        sj_unorm_.rz() = cos(dip_);
#else
        sj_A_ = 2.0 * sj_rad_; // Assumes thickness of 1 in 2D
        sj_unorm_.rx() = sin(dip_);
        sj_unorm_.ry() = cos(dip_);
#endif
        DVect nc = state->getMechanicalContact()->getContact()->getNormal();
        // Set the flip boolean so that the ordering is correct from side1 to side2
        isjFlip_ = ( ((nc | sj_unorm_) >= 0.0 ) ? false : true );
    }

    static const QString kn("sj_kn");
    bool ContactModelSmoothJoint::updateKn(const IContactMechanical *con) {
        assert(con);
        QVariant v1 = con->getEnd1()->getProperty(kn);
        QVariant v2 = con->getEnd2()->getProperty(kn);
        if (!v1.isValid() || !v2.isValid())
            return false;
        double kn1 = v1.toDouble();
        double kn2 = v2.toDouble();
        double val = sj_kn_;
        if (kn1 && kn2)
            sj_kn_ = kn1*kn2/(kn1+kn2);
        else if (kn1)
            sj_kn_ = kn1;
        else if (kn2)
            sj_kn_ = kn2;
        return ( (sj_kn_ != val) );
    }

    static const QString ks("sj_ks");
    bool ContactModelSmoothJoint::updateKs(const IContactMechanical *con) {
        assert(con);
        QVariant v1 = con->getEnd1()->getProperty(ks);
        QVariant v2 = con->getEnd2()->getProperty(ks);
        if (!v1.isValid() || !v2.isValid())
            return false;
        double kn1 = v1.toDouble();
        double kn2 = v2.toDouble();
        double val = sj_ks_;
        if (kn1 && kn2)
            sj_ks_ = kn1*kn2/(kn1+kn2);
        else if (kn1)
            sj_ks_ = kn1;
        else if (kn2)
            sj_ks_ = kn2;
        return ( (sj_ks_ != val) );
    }

    static const QString fric("sj_fric");
    bool ContactModelSmoothJoint::updateFric(const IContactMechanical *con) {
        assert(con);
        QVariant v1 = con->getEnd1()->getProperty(fric);
        QVariant v2 = con->getEnd2()->getProperty(fric);
        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 = sj_fric_;
        sj_fric_ = std::min(fric1,fric2);
        return ( (sj_fric_ != val) );
    }

    bool ContactModelSmoothJoint::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 kwKn:  { //sj_kn_
                if (inheritanceField_ & sjKnMask)
                    ret = updateKn(c);
                break;
            }
        case kwKs:  { //sj_ks_
                if (inheritanceField_ & sjKsMask)
                    ret =updateKs(c);
                break;
            }
        case kwFric:  { //sj_fric_
                if (inheritanceField_ & sjFricMask)
                    updateFric(c);
                break;
            }
        }

        if (ret)
            updateEffectiveTranslationalStiffness();
        return ret;
    }

    void ContactModelSmoothJoint::updateEffectiveTranslationalStiffness() {
        effectiveTranslationalStiffness_ = DVect2(sj_kn_,sj_ks_)*sj_A_;
    }
     
    bool ContactModelSmoothJoint::forceDisplacementLaw(ContactModelMechanicalState *state,const double &) {
        assert(state);
         
        if (geomRecomp_)
            updateAreaAndNormal(state);

        // Skip out if this should not be active
        if (sj_large_ && state->gap_ > 0.0 && sj_state_ != 3) {
            sj_Fn_ = 0.0;
            sj_Fs_ = DVect(0.0);
            sj_gap_ = 0.0;
            sj_Un_ = 0.0;
            sj_Us_ = DVect(0.0);
            return false;
        }

        // Get the translational increment in the global system
        CAxes localSys = state->getMechanicalContact()->getContact()->getLocalSystem();
        DVect tInc = localSys.toGlobal(state->relativeTranslationalIncrement_);
         
        // Now have the translational increment and the normal in the global coordinate system
        // Compute the normal/shear displacement increments along the sj
        if (isjFlip_)
            tInc *= -1.0;
        double nInc = tInc | sj_unorm_;
        DVect sInc = tInc - sj_unorm_ * nInc;
        sj_Un_ -= nInc;
        sj_Us_ += sInc;

        double nElInc(0.0);
        DVect sElInc(0.0);

        double g0 = sj_gap_, g1 = sj_gap_ + nInc;
        sj_gap_ = g1;
         
        if (!state->canFail_ || sj_state_ == 3 )  { // Bonded
            nElInc  = nInc;
            sElInc = sInc;
        } else {
            if ( g0 <= 0.0 ) {
                if ( g1 <= 0.0 ) {
                    nElInc  = nInc;
                    sElInc = sInc;
                } else { // g1 > 0.0
                    double xi = -g0 / (g1 - g0);
                    nElInc  = nInc  * xi;
                    sElInc = sInc * xi;
                }
            } else { // g0 > 0.0
                if ( g1 >= 0.0 ) {
                    nElInc = 0.0;
                    sElInc.fill(0.0);
                } else { // g1 < 0.0
                    double xi = -g0 / (g1 - g0);
                    nElInc  = nInc  * (1.0 - xi);
                    sElInc = sInc * (1.0 - xi);
                }
            }
        }

        double del_Fn = -sj_kn_ * sj_A_ * nElInc ;
        sj_Fn_ += del_Fn ;
        int slideChange(-1);
        DVect sj_Fs_old = sj_Fs_;

        if ( state->canFail_ && sj_state_ < 3 ) {  // Coulomb sliding with dilation
            if ( sj_Fn_ <= 0.0 ) {
                sj_Fn_ = 0.0; 
                sj_Fs_.fill(0.0);
            } else {
                DVect del_Fs = sElInc * (-sj_ks_ * sj_A_);
                sj_Fs_ += del_Fs;
                double max_Fs = sj_Fn_ * sj_fric_;
                double magFs = sj_Fs_.mag();
                if ( magFs > max_Fs ) { // sliding
                    if (!isjSliding_) {
                        isjSliding_ = true;
                        slideChange = 0;
                    }
                    if ( sj_ks_ > 0.0 )  // dilation
                        sj_Fn_ += ( (magFs - max_Fs) / sj_ks_ ) * sj_kn_ * tan(sj_da_);
                    double rat = max_Fs / magFs ;
                    sj_Fs_ *= rat ;
                } else {
                    if (isjSliding_) {
                        isjSliding_ = false;
                        slideChange = 1;
                    }
                }
            }
        } else { // bonded behavior
            if ( state->canFail_ && sj_Fn_ <= -sj_bns_ * sj_A_) {
                sj_state_ = 1;  
                if (cmEvents_[fBondBreak] >= 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(sj_state_);
                    arg.push_back(p1);
                    p1.setValue(sj_bns_);
                    arg.push_back(p1);
                    IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
                    fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
                }
                sj_Fn_ = 0.0; 
                sj_Fs_.fill(0.0);
            } else {
                DVect del_Fs = sElInc * (-sj_ks_ * sj_A_);
                sj_Fs_ += del_Fs;
                double magFs = sj_Fs_.mag();
                double ss = calcBSS();
                if ( state->canFail_ && magFs >= ss * sj_A_) { // break in shear
                    sj_state_ = 2;  
                    if (cmEvents_[fBondBreak] >= 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(sj_state_);
                        arg.push_back(p1);
                        p1.setValue(ss);
                        arg.push_back(p1);
                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
                        fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
                    }

                    if ( sj_Fn_ < 0.0 ) {
                        sj_Fn_ = 0.0; 
                        sj_Fs_.fill(0.0); 
                    }  else { // was in compression // was in tension
                        double max_Fs = sj_Fn_ * sj_fric_ ;
                        if ( magFs > max_Fs) { // sliding, but no dilation
                            if (!isjSliding_) {
                                isjSliding_ = true;
                                slideChange = 0;
                            }
                            double rat = max_Fs / magFs ;
                            sj_Fs_ *= rat ;
                        } else {
                            if (isjSliding_ == true) {
                                isjSliding_ = false;
                                slideChange = 1;
                            }
                        }
                    }
                }
            }
        }
        if (slideChange >= 0 && 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(slideChange);
            arg.push_back(p1);
            IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
            fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
        }

        // Have updated the normal and shear forces so need to put them into the contact local
        // coordinate system and update the forces 
        DVect Fj = sj_unorm_ * sj_Fn_ + sj_Fs_;   
        if (isjFlip_)
            Fj *= -1.0;

        // Return the correct activity status
        bool isactive = true;
        if (sj_large_ && state->gap_ > 0.0 && sj_state_ != 3) {
            sj_Fn_ = 0.0;
            sj_Fs_ = DVect(0.0);
            sj_gap_ = 0.0;
            sj_Un_ = 0.0;
            sj_Us_ = DVect(0.0);
            isactive = false;
        }

        // update energies
       if (state->trackEnergy_) {
            assert(energies_);
            energies_->estrain_ =  0.0;
            if (sj_kn_)
                energies_->estrain_ = 0.5*sj_Fn_*sj_Fn_/sj_kn_;
            if (sj_ks_) {
                double smag2 = sj_Fs_.mag2();
                energies_->estrain_ += 0.5*smag2 / sj_ks_;
                if (isjSliding_) {
                    DVect avg_F_s = (sj_Fs_old + sj_Fs_)*0.5;
                    DVect u_s_el =  (sj_Fs_ - sj_Fs_old) / sj_ks_;
                    energies_->eslip_ -= std::min(0.0,(avg_F_s | (sInc + u_s_el)));
                }
            }
        }

        return isactive ;
    }

    void ContactModelSmoothJoint::setForce(const DVect &v,IContact *c) {
        // this is in the local coordinate system
        CAxes localSys = c->getLocalSystem();
        DVect globForce = localSys.toGlobal(v);
        if (isjFlip_)
            globForce *= -1.0;
        sj_Fn_ = (sj_unorm_ | globForce);
        sj_Fs_ = globForce - sj_unorm_ * sj_Fn_;
    }

    DVect ContactModelSmoothJoint::getForce(const IContactMechanical *c) const {
        CAxes localSys = c->getContact()->getLocalSystem();
        DVect Fj = sj_unorm_ * sj_Fn_ + sj_Fs_;   
        if (isjFlip_)
            Fj *= -1.0;
        DVect ret(localSys.toLocal(Fj));
        return ret;
    }

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

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

    double ContactModelSmoothJoint::calcBSS() const {
        if (sj_A_ > 0) {
            double dSigma = sj_Fn_ / sj_A_;
            return dSigma >= -sj_bns_ ? sj_bcoh_ + (tan(sj_bfa_) * dSigma)
                                                    : sj_bcoh_ + (tan(sj_bfa_) * (-sj_bns_));
        }
        else 
            return 0.0;
    }


} // namespace itascaxd
// EoF

Top