Linear Contact Model Implementation

See this page for the documentation of this contact model.

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

#include "contactmodel/src/contactmodelmechanical.h"

#ifdef LINEAR_LIB
#  define LINEAR_EXPORT EXPORT_TAG
#elif defined(NO_MODEL_IMPORT)
#  define LINEAR_EXPORT
#else
#  define LINEAR_EXPORT IMPORT_TAG
#endif

namespace cmodelsxd {
    using namespace itasca;

    class ContactModelLinear : public ContactModelMechanical {
    public:
        // Constructor: Set default values for contact model properties.
        LINEAR_EXPORT ContactModelLinear();
        // Destructor, called when contact is deleted: free allocated memory, etc.
        LINEAR_EXPORT virtual ~ContactModelLinear();
        // Contact model name (used as keyword for commands and FISH).
        virtual QString  getName() const { return "linear"; }
        // The index provides a quick way to determine the type of contact model.
        // Each type of contact model in PFC must have a unique index; this is assigned
        // by PFC when the contact model is loaded. This index should be set to -1
        virtual void     setIndex(int i) { index_=i;}
        virtual int      getIndex() const {return index_;}
        // Contact model version number (e.g., MyModel05_1). The version number can be
        // accessed during the save-restore operation (within the archive method,
        // testing {stream.getRestoreVersion() == getMinorVersion()} to allow for 
        // future modifications to the contact model data structure.
        virtual uint     getMinorVersion() const;
        // Copy the state information to a newly created contact model.
        // Provide access to state information, for use by copy method.
        virtual void     copy(const ContactModel *c);
        // Provide save-restore capability for the state information.
        virtual void     archive(ArchiveStream &); 
        // Enumerator for the properties.
        enum PropertyKeys { 
              kwKn=1
            , kwKs                            
            , kwFric   
            , kwLinF
            , kwLinS
            , kwLinMode
            , kwRGap
            , kwEmod
            , kwKRatio
            , kwDpNRatio 
            , kwDpSRatio
            , kwDpMode 
            , kwDpF
            , kwUserArea
        };
        // Contact model property names in a comma separated list. The order corresponds with
        // the order of the PropertyKeys enumerator above. One can visualize any of these 
        // properties in PFC automatically. 
        virtual QString  getProperties() const { 
            return "kn"
                   ",ks"
                   ",fric"
                   ",lin_force"
                   ",lin_slip"
                   ",lin_mode"
                   ",rgap"
                   ",emod"
                   ",kratio"
                   ",dp_nratio"
                   ",dp_sratio"
                   ",dp_mode"
                   ",dp_force"
                   ",user_area";
        }
        // Enumerator for the energies.
        enum EnergyKeys { 
            kwEStrain=1
          , kwESlip
          , kwEDashpot
        };
        // Contact model energy names in a comma separated list. The order corresponds with
        // the order of the EnergyKeys enumerator above. 
        virtual QString  getEnergies() const { 
            return "energy-strain"
                    ",energy-slip"
                    ",energy-dashpot";
        }
        // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
        virtual double   getEnergy(uint i) const; 
        // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1) 
        // returns wther or not the estrain energy is accumulated which is false).
        virtual bool     getEnergyAccumulate(uint i) const;
        // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
        virtual void     setEnergy(uint i,const double &d); // Base 1
        // Activate the energy. This is only called if the energy tracking is enabled. 
        virtual void     activateEnergy() { if (energies_) return; energies_ = NEWC(Energies());}
        // Returns whether or not the energy tracking has been enabled for this contact.
        virtual bool     getEnergyActivated() const {return (energies_ != 0);}

        // Enumerator for contact model related FISH callback events. 
        enum FishCallEvents {
            fActivated=0
            ,fSlipChange
        };
        // Contact model FISH callback event names in a comma separated list. The order corresponds with
        // the order of the FishCallEvents enumerator above. 
        virtual QString  getFishCallEvents() const { 
            return 
                "contact_activated"
                ",slip_change"; 
        }

        // Return the specified contact model property.
        virtual QVariant getProperty(uint i,const IContact *) const;
        // The return value denotes whether or not the property corresponds to the global
        // or local coordinate system (TRUE: global system, FALSE: local system). The
        // local system is the contact-plane system (nst) defined as follows.
        // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
        // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
        // and tc are unit vectors in directions of the nst axes.
        // This is used when rendering contact model properties that are vectors.
        virtual bool     getPropertyGlobal(uint i) const;
        // Set the specified contact model property, ensuring that it is of the correct type
        // and within the correct range --- if not, then throw an exception.
        // The return value denotes whether or not the update has affected the timestep
        // computation (by having modified the translational or rotational tangent stiffnesses).
        // If true is returned, then the timestep will be recomputed.
        virtual bool     setProperty(uint i,const QVariant &v,IContact *);
        // The return value denotes whether or not the property is read-only
        // (TRUE: read-only, FALSE: read-write).
        virtual bool     getPropertyReadOnly(uint i) const;

        // The return value denotes whether or not the property is inheritable
        // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
        // the endPropertyUpdated method.
        virtual bool     supportsInheritance(uint i) const; 
        // Return whether or not inheritance is enabled for the specified property.
        virtual bool     getInheritance(uint i) const { assert(i<32); quint32 mask = to<quint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
        // Set the inheritance flag for the specified property.
        virtual void     setInheritance(uint i,bool b) { assert(i<32); quint32 mask = to<quint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }

        // Enumerator for contact model methods.
        enum MethodKeys { kwDeformability=1, kwArea};
        // Contact model methoid names in a comma separated list. The order corresponds with
        // the order of the MethodKeys enumerator above.  
        virtual QString  getMethods() const { return "deformability,area";}
        // Return a comma seprated list of the contact model method arguments (base 1).
        virtual QString  getMethodArguments(uint i) const; 
        // Set contact model method arguments (base 1). 
        // The return value denotes whether or not the update has affected the timestep
        // computation (by having modified the translational or rotational tangent stiffnesses).
        // If true is returned, then the timestep will be recomputed.
        virtual bool     setMethod(uint i,const QVector<QVariant> &vl,IContact *con=0); 

        // Prepare for entry into ForceDispLaw. The validate function is called when:
        // (1) the contact is created, (2) a property of the contact that returns a true via
        // the setProperty method has been modified and (3) when a set of cycles is executed
        // via the {cycle N} command.
        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
        virtual bool    validate(ContactModelMechanicalState *state,const double &timestep);
        // The endPropertyUpdated method is called whenever a surface property (with a name
        // that matches an inheritable contact model property name) of one of the contacting
        // pieces is modified. This allows the contact model to update its associated
        // properties. The return value denotes whether or not the update has affected
        // the time step computation (by having modified the translational or rotational
        // tangent stiffnesses). If true is returned, then the time step will be recomputed.  
        virtual bool    endPropertyUpdated(const QString &name,const IContactMechanical *c);
        // The forceDisplacementLaw function is called during each cycle. Given the relative
        // motion of the two contacting pieces (via
        //   state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
        //   state->relativeAngularIncrement_       (Dtt, Dtbs, Dtbt)
        //     Ddn  : relative normal-displacement increment, Ddn > 0 is opening
        //     Ddss : relative  shear-displacement increment (s-axis component)
        //     Ddst : relative  shear-displacement increment (t-axis component)
        //     Dtt  : relative twist-rotation increment
        //     Dtbs : relative  bend-rotation increment (s-axis component)
        //     Dtbt : relative  bend-rotation increment (t-axis component)
        //       The relative displacement and rotation increments:
        //         Dd = Ddn*nc + Ddss*sc + Ddst*tc
        //         Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
        //       where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
        //       [see {Table 1: Contact State Variables} in PFC Model Components:
        //       Contacts and Contact Models: Contact Resolution]
        // ) and the contact properties, this function must update the contact force and
        // moment.
        //   The force_ is acting on piece 2, and is expressed in the local coordinate system
        //   (defined in getPropertyGlobal) such that the first component positive denotes
        //   compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
        //   in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to 
        //   get the total moment. 
        // The return value indicates the contact activity status (TRUE: active, FALSE:
        // inactive) during the next cycle.
        // Additional information:
        //   * If state->activated() is true, then the contact has just become active (it was
        //     inactive during the previous time step).
        //   * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
        //     the forceDispLaw handle the case of {state->canFail_ == true}.
        virtual bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep);
        // The getEffectiveXStiffness functions return the translational and rotational
        // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
        // the translational tangent shear stiffness is zero (but this stiffness reduction
        // is typically ignored when computing a stable time step). If the contact model
        // includes a dashpot, then the translational stiffnesses must be increased (see
        // Potyondy (2009)).
        //   [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
        //   Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
        //   December 7, 2009.]
        virtual DVect2  getEffectiveTranslationalStiffness() const { return effectiveTranslationalStiffness_; }
        virtual DAVect  getEffectiveRotationalStiffness() const { return DAVect(0.0);}

        // Return a new instance of the contact model. This is used in the CMAT
        // when a new contact is created. 
        virtual ContactModelLinear *clone() const { return NEWC(ContactModelLinear()); }
        // The getActivityDistance function is called by the contact-resolution logic when
        // the CMAT is modified. Return value is the activity distance used by the
        // checkActivity function. 
        virtual double              getActivityDistance() const {return rgap_;}
        // The isOKToDelete function is called by the contact-resolution logic when...
        // Return value indicates whether or not the contact may be deleted.
        // If TRUE, then the contact may be deleted when it is inactive.
        // If FALSE, then the contact may not be deleted (under any condition).
        virtual bool                isOKToDelete() const { return !isBonded(); }
        // Zero the forces and moments stored in the contact model. This function is called
        // when the contact becomes inactive.
        virtual void                resetForcesAndMoments() { lin_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 &d) { userArea_ = d; }
        // The checkActivity function is called by the contact-resolution logic when...
        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
        // A contact with the linear model is active if the contact gap is
        // less than or equal to zero.
        virtual bool     checkActivity(const double &gap) { return  gap <= rgap_; }

        // Returns the sliding state (FALSE is returned if not implemented).
        virtual bool     isSliding() const { return lin_S_; }
        // Returns the bonding state (FALSE is returned if not implemented).
        virtual bool     isBonded() const { return false; }

        // Both of these methods are called only for contacts with facets where the wall 
        // resolution scheme is set the full. In such cases one might wish to propagate 
        // contact state information (e.g., shear force) from one active contact to another. 
        // See the Faceted Wall section in the documentation. 
        virtual void     propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes());
        virtual void     setNonForcePropsFrom(IContactModel *oldCM);

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

        // Methods to get and set properties. 
        const double & kn() const {return kn_;}
        void           kn(const double &d) {kn_=d;}
        const double & ks() const {return ks_;}
        void           ks(const double &d) {ks_=d;}
        const double & fric() const {return fric_;}
        void           fric(const double &d) {fric_=d;}
        const DVect &  lin_F() const {return lin_F_;}
        void           lin_F(const DVect &f) { lin_F_=f;}
        bool           lin_S() const {return lin_S_;}
        void           lin_S(bool b) { lin_S_=b;}
        uint           lin_mode() const {return lin_mode_;}
        void           lin_mode(uint i) { lin_mode_= i;}
        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;}

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

    private:
        // Index - used internally by PFC. Should be set to -1 in the cpp file. 
        static int index_;

        // Structure to store the 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
        };

        // Structure to store dashpot quantities. 
        struct dpProps {
            dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(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
        };

        bool   updateKn(const IContactMechanical *con);
        bool   updateKs(const IContactMechanical *con);
        bool   updateFric(const IContactMechanical *con);

        void   updateEffectiveStiffness(ContactModelMechanicalState *state);

        void   setDampCoefficients(const double &mass,double *vcn,double *vcs);

        // Contact model inheritance fields.
        quint32 inheritanceField_;

        // Effective translational stiffness.
        DVect2  effectiveTranslationalStiffness_;

        // linear model properties
        double      kn_;        // Normal stiffness
        double      ks_;        // Shear stiffness
        double      fric_;      // Coulomb friction coefficient
        DVect       lin_F_;     // Force carried in the linear model
        bool        lin_S_;     // The current slip state
        uint        lin_mode_;  // Specifies absolute (0) or incremental (1) calculation mode 
        double      rgap_;      // Reference gap 
        dpProps *   dpProps_;   // The viscous properties
        double      userArea_;   // Area as specified by the user 

        Energies *   energies_; // The energies

         
    };
} // namespace cmodelsxd
// EoF

Top

contactmodellinear.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
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
// contactmodellinear.cpp
#include "contactmodellinear.h"

#include "module/interface/icontactmechanical.h"
#include "module/interface/icontact.h"
#include "module/interface/ipiecemechanical.h"
#include "module/interface/ipiece.h"
#include "module/interface/ifishcalllist.h"

#include "../version.txt"

#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 LINEAR_LIB
    int __stdcall DllMain(void *,unsigned, void *) {
        return 1;
    }

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

namespace cmodelsxd {
    static const quint32 linKnMask      = 0x00002; // Base 1!
    static const quint32 linKsMask      = 0x00004;
    static const quint32 linFricMask    = 0x00008;

    using namespace itasca;

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

    ContactModelLinear::ContactModelLinear() : inheritanceField_(linKnMask|linKsMask|linFricMask) 
                                             , effectiveTranslationalStiffness_(DVect2(0.0)) 
                                             , kn_(0.0)
                                             , ks_(0.0)
                                             , fric_(0.0)
                                             , lin_F_(DVect(0.0))
                                             , lin_S_(false)
                                             , lin_mode_(0)
                                             , rgap_(0.0)
                                             , dpProps_(0)
                                             , userArea_(0)
                                             , energies_(0) {
    }

    ContactModelLinear::~ContactModelLinear() {
        // Make sure to clean up after yourself!
        if (dpProps_)
            delete dpProps_;
        if (energies_)
            delete energies_;
    }

    void ContactModelLinear::archive(ArchiveStream &stream) {
        // The stream allows one to archive the values of the contact model
        // so that it can be saved and restored. The minor version can be
        // used here to allow for incremental changes to the contact model too. 
        stream & kn_;
        stream & ks_;
        stream & fric_;
        stream & lin_F_;
        stream & lin_S_;
        stream & lin_mode_;
 
        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_; 
            }
            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_; 
            }
            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() == getMinorVersion())
            stream & rgap_;

        if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 2) 
            stream & userArea_;
    }

    void ContactModelLinear::copy(const ContactModel *cm) {
        // Copy all of the contact model properties. Used in the CMAT 
        // when a new contact is created. 
        ContactModelMechanical::copy(cm);
        const ContactModelLinear *in = dynamic_cast<const ContactModelLinear*>(cm);
        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
        kn(in->kn());
        ks(in->ks());
        fric(in->fric());
        lin_F(in->lin_F());
        lin_S(in->lin_S());
        lin_mode(in->lin_mode());
        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()); 
        }
        if (in->hasEnergies()) {
            if (!energies_)
                energies_ = NEWC(Energies());
            estrain(in->estrain());
            eslip(in->eslip());
            edashpot(in->edashpot());
        }
        userArea_ = in->userArea_;
        inheritanceField(in->inheritanceField());
        effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
    }


    QVariant ContactModelLinear::getProperty(uint i,const IContact *con) const {
        // Return the property. The IContact pointer is provided so that 
        // more complicated properties, depending on contact characteristics,
        // can be calcualted. 
        QVariant var;
        switch (i) {
        case kwKn:        return kn_;
        case kwKs:        return ks_;
        case kwFric:      return fric_;
        case kwLinF:      var.setValue(lin_F_); return var;
        case kwLinS:      return lin_S_;
        case kwLinMode:   return lin_mode_;
        case kwRGap:      return rgap_;
        case kwEmod: {
                        const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
                        if (c ==nullptr) return 0.0;
                        double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
                        double rsum(0.0);
                        if (c->getEnd1Curvature().y())
                            rsum += 1.0/c->getEnd1Curvature().y();
                        if (c->getEnd2Curvature().y())
                            rsum += 1.0/c->getEnd2Curvature().y();
                        if (userArea_) {
#ifdef THREED
                            rsq = std::sqrt(userArea_ / dPi);
#else
                            rsq = userArea_ / 2.0;
#endif        
                            rsum = rsq + rsq;
                            rsq = 1. / rsq;
                        }
#ifdef TWOD               
                        return (kn_ * rsum * rsq / 2.0);
#else                   
                        return (kn_ * rsum * rsq * rsq) / dPi;
#endif                    
                    }
        case kwKRatio:    return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
        case kwDpNRatio:  return dpProps_ ? dpProps_->dp_nratio_ : 0;
        case kwDpSRatio:  return dpProps_ ? dpProps_->dp_sratio_ : 0;
        case kwDpMode:    return dpProps_ ? dpProps_->dp_mode_ : 0;
        case kwDpF: {
                dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
                return var;
            }
        case kwUserArea:    return userArea_;
        }
        assert(0);
        return QVariant();
    }

    bool ContactModelLinear::getPropertyGlobal(uint i) const {
        // Returns whether or not a property is held in the global axis system (TRUE)
        // or the local system (FALSE). Used by the plotting logic.
        switch (i) {
        case kwLinF:   
        case kwDpF:  
            return false;
        }
        return true;
    }

    bool ContactModelLinear::setProperty(uint i,const QVariant &v,IContact *) {
        // Set a contact model property. Return value indicates that the timestep
        // should be recalculated. 
        dpProps dp;
        switch (i) {
        case kwKn: {
                if (!v.canConvert<double>())
                    throw Exception("kn must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative kn not allowed.");
                kn_ = val;
                return true;
            }
        case kwKs: {
                if (!v.canConvert<double>())
                    throw Exception("ks must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative ks not allowed.");
                ks_ = 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 kwLinF: {
                if (!v.canConvert<DVect>())
                    throw Exception("lin_force must be a vector.");
                DVect val(v.value<DVect>());
                lin_F_ = val;
                return false;
            }
        case kwLinMode: {
                if (!v.canConvert<uint>())
                    throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
                uint val(v.toUInt());
                if (val >1)
                    throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
                lin_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 kwDpF: {
                if (!v.canConvert<DVect>())
                    throw Exception("dp_force must be a vector.");
                DVect val(v.value<DVect>());
                if (val.fsum() == 0.0 && !dpProps_)
                    return false;
                if (!dpProps_)
                    dpProps_ = NEWC(dpProps());
                dpProps_->dp_F_ = val;
                return false;
            }
        case kwUserArea: {
                if (!v.canConvert<double>())
                    throw Exception("user_area must be a double.");
                double val(v.toDouble());
                if (val < 0.0)
                    throw Exception("Negative user_area not allowed.");
                userArea_ = val;
                return true;
            }
        }
        return false;
    }

    bool ContactModelLinear::getPropertyReadOnly(uint i) const {
        // Returns TRUE if a property is read only or FALSE otherwise. 
        switch (i) {
        case kwDpF:
        case kwLinS:
        case kwEmod:
        case kwKRatio:
            return true;
        default:
            break;
        }
        return false;
    }

    bool ContactModelLinear::supportsInheritance(uint i) const {
        // Returns TRUE if a property supports inheritance or FALSE otherwise. 
        switch (i) {
        case kwKn:
        case kwKs:
        case kwFric:
            return true;
        default:
            break;
        }
        return false;
    }

    QString  ContactModelLinear::getMethodArguments(uint i) const {
        // Return a list of contact model method argument names. 
        switch (i) {
        case kwDeformability:
            return "emod,kratio";
        case kwArea:
            return QString();
        }
        assert(0);
        return QString();
    }

    bool ContactModelLinear::setMethod(uint i,const QVector<QVariant> &vl,IContact *con) {
        // Apply the specified method. 
        IContactMechanical *c(convert_getcast<IContactMechanical>(con));
        switch (i) {
        case kwDeformability: {
                double emod;
                double krat;
                if (vl.at(0).isNull()) 
                    throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
                emod = vl.at(0).toDouble();
                if (emod<0.0)
                    throw Exception("Negative emod not allowed in contact model %1.",getName());
                if (vl.at(1).isNull()) 
                    throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
                krat = vl.at(1).toDouble();
                if (krat<0.0)
                    throw Exception("Negative stiffness ratio not allowed in contact model %1.",getName());
                double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
                double rsum(0.0);
                if (c->getEnd1Curvature().y())
                    rsum += 1.0/c->getEnd1Curvature().y();
                if (c->getEnd2Curvature().y())
                    rsum += 1.0/c->getEnd2Curvature().y();
                if (userArea_) {
#ifdef THREED
                    rsq = std::sqrt(userArea_ / dPi);
#else
                    rsq = userArea_ / 2.0;
#endif        
                    rsum = rsq + rsq;
                    rsq = 1. / rsq;
                }
#ifdef TWOD
                kn_ = 2.0 * emod / (rsq * rsum);
#else
                kn_ = dPi * emod / (rsq * rsq * rsum);
#endif
                ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
                setInheritance(1,false);
                setInheritance(2,false);
                return true;
            }
        case kwArea: {
                if (!userArea_) {
                    double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
#ifdef THREED
                    userArea_ = rsq * rsq * dPi;
#else
                    userArea_ = rsq * 2.0;
#endif                            
                }
                return true;
            }
        }
        return false;
    }

    double ContactModelLinear::getEnergy(uint i) const {
        // Return an energy value. 
        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 ContactModelLinear::getEnergyAccumulate(uint i) const {
        // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
        switch (i) {
        case kwEStrain:  return false;
        case kwESlip:    return true;
        case kwEDashpot:    return true;
        }
        assert(0);
        return false;
    }

    void ContactModelLinear::setEnergy(uint i,const double &d) {
        // Set an energy value. 
        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 ContactModelLinear::validate(ContactModelMechanicalState *state,const double &) {
        // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
        // (1) the contact is created, (2) a property of the contact that returns a true via
        // the setProperty method has been modified and (3) when a set of cycles is executed
        // via the {cycle N} command.
        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
        assert(state);
        const IContactMechanical *c = state->getMechanicalContact(); 
        assert(c);

        if (state->trackEnergy_)
            activateEnergy();

        if (inheritanceField_ & linKnMask)
            updateKn(c);
        if (inheritanceField_ & linKsMask)
            updateKs(c);
        if (inheritanceField_ & linFricMask)
            updateFric(c);

        updateEffectiveStiffness(state);
        return checkActivity(state->gap_);
    }

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

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

    static const QString fricstr("fric");
    bool ContactModelLinear::updateFric(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 ContactModelLinear::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
        // The endPropertyUpdated method is called whenever a surface property (with a name
        // that matches an inheritable contact model property name) of one of the contacting
        // pieces is modified. This allows the contact model to update its associated
        // properties. The return value denotes whether or not the update has affected
        // the time step computation (by having modified the translational or rotational
        // tangent stiffnesses). If true is returned, then the time step will be recomputed.  
        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:  { //kn
                if (inheritanceField_ & linKnMask)
                    ret = updateKn(c);
                break;
            }
        case kwKs:  { //ks
                if (inheritanceField_ & linKsMask)
                    ret =updateKs(c);
                break;
            }
        case kwFric:  { //fric
                if (inheritanceField_ & linFricMask)
                    updateFric(c);
                break;
            }
        }
        return ret;
    }

    void ContactModelLinear::updateEffectiveStiffness(ContactModelMechanicalState *) {
        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 ContactModelLinear::forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) {
        assert(state);

        // Current overlap
        double overlap = rgap_ - state->gap_;
        // Relative translational increment
        DVect trans = state->relativeTranslationalIncrement_;
        // Correction factor to account for when the contact becomes newly active.
        // We estimate the time of activity during the timestep when the contact has first 
        // become active and scale the forces accordingly.
        double correction = 1.0;

        // The contact was just activated from an inactive state
        if (state->activated()) {
            // Trigger the FISH callback if one is hooked up to the 
            // contact_activated event.
            if (cmEvents_[fActivated] >= 0) {
                // An FArray of QVariant is returned and these will be passed
                // to the FISH function as an array of FISH symbols as the second
                // argument to the FISH callback function. 
                FArray<QVariant,2> arg;
                QVariant v;
                // Just put a pointer to the contact in the array and return it.
                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]);
            }
            // Calculate the correction factor.
            if (lin_mode_ == 0 && trans.x()) {
                correction = -1.0*overlap / trans.x();
                if (correction < 0)
                    correction = 1.0;
            }
        }

        // Angular dispacement increment.
        DAVect ang  = state->relativeAngularIncrement_;
        DVect lin_F_old = lin_F_;

        if (lin_mode_ == 0) 
            lin_F_.rx() = overlap * kn_; //  absolute mode for normal force calculation
        else
            lin_F_.rx() -= correction * trans.x() * kn_; // incremental  mode for normal force calculation

        // Normal force can only be positive.
        lin_F_.rx() = std::max(0.0,lin_F_.x());

        // Calculate the shear force.
        DVect sforce(0.0);
        // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
        // Loop over the shear components (note: the 0 component is the normal component)
        // and calculate the shear force.
        for (int i=1; i<dim; ++i)
            sforce.rdof(i) = lin_F_.dof(i) - trans.dof(i) * ks_ * correction;

        // The canFail flag corresponds to whether or not the contact can undergo non-linear
        // force-displacement response. If the SOLVE ELASTIC command is given then the 
        // canFail state is set to FALSE. Otherwise it is always TRUE. 
        if (state->canFail_) {
            // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
            double crit = lin_F_.x() * fric_;
            // The is the magnitude of the shear force.
            double sfmag = sforce.mag();
            // Sliding occurs when the magnitude of the shear force is greater than the 
            // critical value.
            if (sfmag > crit) {
                // Lower the shear force to the critical value for sliding.
                double rat = crit / sfmag;
                sforce *= rat;
                // Handle the slip_change event if one has been hooked up. Sliding has commenced.  
                if (!lin_S_ && cmEvents_[fSlipChange] >= 0) {
                    FArray<QVariant,3> arg;
                    QVariant p1;
                    // Put a pointer to the contact in the array plus 0 to indicate slip has initiated.
                    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]);
                }
                lin_S_ = true;
            } else {
                // Handle the slip_change event if one has been hooked up and
                // the contact was previously sliding. Sliding has ceased.  
                if (lin_S_) {
                    if (cmEvents_[fSlipChange] >= 0) {
                        FArray<QVariant,3> arg;
                        QVariant p1;
                        // Put a pointer to the contact in the array plus 1 to indicate slip has ceased.
                        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]);
                    }
                    lin_S_ = false;
                }
            }
        }
        
        // Set the shear components of the total force.
        for (int i=1; i<dim; ++i)
            lin_F_.rdof(i) = sforce.dof(i);
         
        // Account for dashpot forces if the dashpot structure has been defined. 
        if (dpProps_) {
            dpProps_->dp_F_.fill(0.0);
            double vcn(0.0), vcs(0.0);
            // Calculate the damping coefficients. 
            setDampCoefficients(state->inertialMass_,&vcn,&vcs);
            // First damp the shear components
            for (int i=1; i<dim; ++i)
                dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep; 
            // Damp the normal component
            dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;       
            // Need to change behavior based on the dp_mode.
            if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3))  { 
                // Limit in tension if not bonded.
                if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
                    dpProps_->dp_F_.rx() = - lin_F_.rx();
            }
            if (lin_S_ && dpProps_->dp_mode_ > 1)  { 
                // Limit in shear if not sliding.
                double dfn = dpProps_->dp_F_.rx();
                dpProps_->dp_F_.fill(0.0); 
                dpProps_->dp_F_.rx() = dfn; 
            }
        }

        //Compute energies if energy tracking has been enabled. 
        if (state->trackEnergy_) {
            assert(energies_);
            energies_->estrain_ =  0.0;
            if (kn_)
                // Calcualte the strain energy. 
                energies_->estrain_ = 0.5*lin_F_.x()*lin_F_.x()/kn_;
            if (ks_) {
                DVect s = lin_F_;
                s.rx() = 0.0;
                double smag2 = s.mag2();
                // Add the shear component of the strain energy.
                energies_->estrain_ += 0.5*smag2 / ks_;

                if (lin_S_) {
                    // If sliding calculate the slip energy and accumulate it.
                    lin_F_old.rx() = 0.0;
                    DVect avg_F_s = (s + lin_F_old)*0.5;
                    DVect u_s_el =  (s - lin_F_old) / ks_;
                    DVect u_s(0.0);
                    for (int i=1; i<dim; ++i)
                        u_s.rdof(i) = trans.dof(i);
                    energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
                }
            }
            if (dpProps_) {
                // Calculate damping energy (accumulated) if the dashpots are active. 
                energies_->edashpot_ -= dpProps_->dp_F_ | trans;
            }
        }

        // This is just a sanity check to ensure, in debug mode, that the force isn't wonky. 
        assert(lin_F_ == lin_F_);
        return true;
    }

    void ContactModelLinear::setForce(const DVect &v,IContact *c) { 
        lin_F(v); 
        if (v.x() > 0) 
            rgap_ = c->getGap() + v.x() / kn_; 
    } 

    void ContactModelLinear::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
        // Only called for contacts with wall facets when the wall resolution scheme
        // is set to full!
        // Only do something if the contact model is of the same type
        if (old->getContactModel()->getName().compare("linear",Qt::CaseInsensitive) == 0 && !isBonded()) {
            ContactModelLinear *oldCm = (ContactModelLinear *)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->lin_F_.y(),oldCm->lin_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->lin_F_.y(),oldCm->lin_F_.z());
            }
            DVect pforce = DVect(0,tpf.x(),tpf.y());
#else
            oldSystem;
            newSystem;
            DVect pforce = DVect(0,oldCm->lin_F_.y());
#endif
            for (int i=1; i<dim; ++i)
                lin_F_.rdof(i) += pforce.dof(i);
            oldCm->lin_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_->edashpot_ = oldCm->energies_->edashpot_;
                energies_->eslip_ = oldCm->energies_->eslip_;
                oldCm->energies_->estrain_ = 0.0;
                oldCm->energies_->edashpot_ = 0.0;
                oldCm->energies_->eslip_ = 0.0;
            }
            rgap_ = oldCm->rgap_;
        }
        assert(lin_F_ == lin_F_);
    }

    void ContactModelLinear::setNonForcePropsFrom(IContactModel *old) {
        // Only called for contacts with wall facets when the wall resolution scheme
        // is set to full!
        // Only do something if the contact model is of the same type
        if (old->getName().compare("linear",Qt::CaseInsensitive) == 0 && !isBonded()) {
            ContactModelLinear *oldCm = (ContactModelLinear *)old;
            kn_ = oldCm->kn_;
            ks_ = oldCm->ks_;
            fric_ = oldCm->fric_;
            lin_mode_ = oldCm->lin_mode_;
            rgap_ = oldCm->rgap_;
            userArea_ = oldCm->userArea_;

            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 ContactModelLinear::getForce(const IContactMechanical *) const {
        DVect ret(lin_F_);
        if (dpProps_)
            ret += dpProps_->dp_F_;
        return ret;
    }

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

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

    void ContactModelLinear::setDampCoefficients(const double &mass,double *vcn,double *vcs) {
        *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(kn_));
        *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(ks_));
    }

} // namespace cmodelsxd
// EoF

Top