Burger’s Contact Model Implementation

See this page for the documentation of this contact model.

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

#include "contactmodel/src/contactmodelmechanical.h"

#ifdef BURGER_LIB
#  define BURGER_EXPORT EXPORT_TAG
#elif defined(NO_MODEL_IMPORT)
#  define BURGER_EXPORT
#else
#  define BURGER_EXPORT IMPORT_TAG
#endif

namespace cmodelsxd {
    using namespace itasca;

    class ContactModelBurger : public ContactModelMechanical {
    public:
        // Constructor: Set default values for contact model properties.
        BURGER_EXPORT ContactModelBurger();
        // Destructor, called when contact is deleted: free allocated memory, etc.
        BURGER_EXPORT virtual ~ContactModelBurger();
        // Contact model name (used as keyword for commands and FISH).
        virtual QString  getName() const { return "burger"; }
        // 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 { 
              kwKnK=1
            , kwCnK                            
            , kwKnM
            , kwCnM                            
            , kwKsK                            
            , kwCsK                            
            , kwKsM
            , kwCsM                            
            , kwMode
            , kwFric   
            , kwF
            , kwS
        };
        // 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 "bur_knk"
                   ",bur_cnk"
                   ",bur_knm"
                   ",bur_cnm"
                   ",bur_ksk"
                   ",bur_csk"
                   ",bur_ksm"
                   ",bur_csm"
                   ",bur_mode"
                   ",bur_fric"
                   ",bur_force"
                   ",bur_slip";
        }
        // 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;

        // 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 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 DVect2(knk_,ksk_); }

        // Return a new instance of the contact model. This is used in the CMAT
        // when a new contact is created. 
        virtual ContactModelBurger *clone() const { return NEWC(ContactModelBurger()); }
        // 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 0.0;}
        // 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() { force_.fill(0.0);}
        virtual void                setForce(const DVect &v,IContact *) { force_ = v; }
        virtual void                setArea(const double &) { throw Exception("The setArea method cannot be used with this contact model."); }
        // The checkActivity function is called by the contact-resolution logic when...
        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
        virtual bool     checkActivity(const double &gap) { return  gap <= 0.0; }

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

        virtual bool    endPropertyUpdated(const QString &,const IContactMechanical *) {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 & knk() const {return knk_;}
        void           knk(const double &d) {knk_=d;}
        const double & cnk() const {return cnk_;}
        void           cnk(const double &d) {cnk_=d;}
        const double & knm() const {return knm_;}
        void           knm(const double &d) {knm_=d;}
        const double & cnm() const {return cnm_;}
        void           cnm(const double &d) {cnm_=d;}
        const double & ksk() const {return ksk_;}
        void           ksk(const double &d) {ksk_=d;}
        const double & csk() const {return csk_;}
        void           csk(const double &d) {csk_=d;}
        const double & ksm() const {return ksm_;}
        void           ksm(const double &d) {ksm_=d;}
        const double & csm() const {return csm_;}
        void           csm(const double &d) {csm_=d;}

        const double & fric() const { return fric_;}
        void           fric(const double &d) {fric_=d;}
        uint           bmode() const {return bmode_;}
        void           bmode(uint b) {bmode_= b;}

        const double & fn0() const { return fn0_;}
        void           fn0(const double &d) {fn0_=d;}
        const double & u_n0() const { return u_n0_;}
        void           u_n0(const double &d) {u_n0_=d;}
        const double & u_nk0() const { return u_nk0_;}
        void           u_nk0(const double &d) {u_nk0_=d;}
        const DVect &  u_sk() const { return u_sk_;}
        void           u_sk(const DVect &v) {u_sk_=v;}

        const double & conAn() const { return conAn_;}
        void           conAn(const double &d) {conAn_=d;}
        const double & conB_An() const { return conB_An_;}
        void           conB_An(const double &d) {conB_An_=d;}
        const double & conCn() const { return conCn_;}
        void           conCn(const double &d) {conCn_=d;}
        const double & conDn() const { return conDn_;}
        void           conDn(const double &d) {conDn_=d;}

        const double & conAs() const { return conAs_;}
        void           conAs(const double &d) {conAs_=d;}
        const double & conB_As() const { return conB_As_;}
        void           conB_As(const double &d) {conB_As_=d;}
        const double & conCs() const { return conCs_;}
        void           conCs(const double &d) {conCs_=d;}
        const double & conDs() const { return conDs_;}
        void           conDs(const double &d) {conDs_=d;}

        const double & tdel() const { return tdel_;}
        void           tdel(const double &d) {tdel_=d;}
        const DVect &  force() const { return force_;}
        void           force(const DVect &v) {force_=v;}
        bool           s() const {return s_;}
        void           s(bool b) {s_= b;}

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

        // Burger model properties
        double  knk_;     // normal stiffness for Kelvin section  (bur_knk)
        double  cnk_;     // normal viscosity for Kelvin section  (bur_cnk)
        double  knm_;     // normal stiffness for Maxwell section (bur_knm)
        double  cnm_;     // normal viscosity for Maxwell section (bur_cnm)
        double  ksk_;     // shear stiffness for Kelvin section   (bur_ksk)
        double  csk_;     // shear viscosity for Kelvin section   (bur_csk)
        double  ksm_;     // shear stiffness for Maxwell section  (bur_ksm)
        double  csm_;     // shear viscosity for Maxwell section  (bur_csm)
        double  fric_;    // friction coefficient                 (bur_fric)
        uint    bmode_;   // FDLaw option, with or without tensile force, default false (with tensile)
        double  fn0_;     // normal contact force 1 step before
        double  u_n0_;    // normal total overlap 1 step before
        double  u_nk0_;   // normal overlap of Kelvin part 1step before
        DVect   u_sk_;    // shear relative displacement of Kelvin part 1step before
        double  conAn_;   // constant A in eq.(), normal
        double  conB_An_; // constant B/A in eq.(), normal
        double  conCn_;   // constant C in eq.(), normal
        double  conDn_;   // constant D in eq.(), normal
        double  conAs_;   // constant A in eq.(), shear
        double  conB_As_; // constant B/A in eq.(), shear
        double  conCs_;   // constant C in eq.(), shear
        double  conDs_;   // constant D in eq.(), shear
        double  tdel_;    // current timestep
        DVect   force_;   // current total force
        bool    s_;       // current sliding state

        // Constants
        inline double A(const double k_k, const double c_k) {    return(1.0 + k_k*tdel_/(2.0*c_k)); }
        inline double B(const double k_k, const double c_k) { return(1.0 - k_k*tdel_/(2.0*c_k)); }
        inline double B_A(const double k_k, const double c_k) { return(B(k_k,c_k)/A(k_k,c_k)); }
        inline double C(const double k_k, const double c_k, const double k_m, const double c_m) {
          return(tdel_/(2.0*c_k*A(k_k,c_k)) + 1.0/k_m + tdel_/(2.0*c_m)); }
        inline double D(const double k_k, const double c_k, const double k_m, const double c_m) {
          return(tdel_/(2.0*c_k*A(k_k,c_k)) - 1.0/k_m + tdel_/(2.0*c_m)); }

    };
} // namespace cmodelsxd
// EoF

Top

contactmodelburger.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
// contactmodelburger.cpp
#include "contactmodelburger.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 BURGER_LIB
    int __stdcall DllMain(void *,unsigned, void *) {
        return 1;
    }

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

namespace cmodelsxd {
    using namespace itasca;

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

    ContactModelBurger::ContactModelBurger() :knk_(0.0)      
                                             ,cnk_(0.0)      
                                             ,knm_(0.0)     
                                             ,cnm_(0.0)     
                                             ,ksk_(0.0)      
                                             ,csk_(0.0)   
                                             ,ksm_(0.0)  
                                             ,csm_(0.0)  
                                             ,fric_(0.0)  
                                             ,bmode_(0)  
                                             ,fn0_(0.0)  
                                             ,u_n0_(0.0)  
                                             ,u_nk0_(0.0)  
                                             ,u_sk_(0.0) 
                                             ,conAn_(0.0) 
                                             ,conB_An_(0.0)
                                             ,conCn_(0.0)   
                                             ,conDn_(0.0)   
                                             ,conAs_(0.0) 
                                             ,conB_As_(0.0)
                                             ,conCs_(0.0)  
                                             ,conDs_(0.0) 
                                             ,tdel_(0.0) 
                                             ,force_(0.0) 
                                             ,s_(false) 
    {  

    }

    ContactModelBurger::~ContactModelBurger() {
    }

    void ContactModelBurger::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 & knk_;     
        stream & cnk_;     
        stream & knm_;     
        stream & cnm_;     
        stream & ksk_;     
        stream & csk_;     
        stream & ksm_;     
        stream & csm_;       
        stream & fric_;    
        stream & bmode_;   
        stream & fn0_;     
        stream & u_n0_;    
        stream & u_nk0_;   
        stream & u_sk_;    
        stream & conAn_;   
        stream & conB_An_; 
        stream & conCn_;   
        stream & conDn_;   
        stream & conAs_;   
        stream & conB_As_; 
        stream & conCs_;   
        stream & conDs_;   
        stream & tdel_;    
        stream & force_;    
        stream & s_;    
    }

    void ContactModelBurger::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 ContactModelBurger *in = dynamic_cast<const ContactModelBurger*>(cm);
        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
        
        knk(in->knk());     
        cnk(in->cnk());     
        knm(in->knm());     
        cnm(in->cnm());     
        ksk(in->ksk());     
        csk(in->csk());     
        ksm(in->ksm());     
        csm(in->csm());     
        fric(in->fric());    
        bmode(in->bmode());   
        fn0(in->fn0());     
        u_n0(in->u_n0());    
        u_nk0(in->u_nk0());   
        u_sk(in->u_sk());    
        conAn(in->conAn());   
        conB_An(in->conB_An()); 
        conCn(in->conCn());   
        conDn(in->conDn());   
        conAs(in->conAs());   
        conB_As(in->conB_As()); 
        conCs(in->conCs());   
        conDs(in->conDs());   
        tdel(in->tdel()); 
        force(in->force());
        s(in->s());
    }


    QVariant ContactModelBurger::getProperty(uint i,const IContact *) const {
        // Return the property. The IContact pointer is provided so that 
        // more complicated properties, depending on contact characteristics,
        // can be calculated. Not used with the Burger model.
        QVariant var;
        switch (i) {
        case kwKnK   : return knk_;
        case kwCnK   : return cnk_;
        case kwKnM   : return knm_;
        case kwCnM   : return cnm_;
        case kwKsK   : return ksk_;
        case kwCsK   : return csk_;
        case kwKsM   : return ksm_;
        case kwCsM   : return csm_;
        case kwMode  : return bmode_;
        case kwFric  : return fric_;
        case kwF     : var.setValue(force_); return var;
        case kwS     : return s_;
        }
        assert(0);
        return QVariant();
    }

    bool ContactModelBurger::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 kwF:   
            return false;
        }
        return true;
    }

    bool ContactModelBurger::setProperty(uint i,const QVariant &v,IContact *) {
        // Set a contact model property. Return value indicates that the timestep
        // should be recalculated. 
        switch (i) {
        case kwKnK: {
                if (!v.canConvert<double>())
                    throw Exception("bur_knk must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative bur_knk not allowed.");
                knk_ = val;
                return true;
            }
        case kwCnK: {
                if (!v.canConvert<double>())
                    throw Exception("bur_cnk must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative bur_cnk not allowed.");
                cnk_ = val;
                return true;
            }
        case kwKnM: {
                if (!v.canConvert<double>())
                    throw Exception("bur_knm must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative bur_knm not allowed.");
                knm_ = val;
                return true;
            }
        case kwCnM: {
                if (!v.canConvert<double>())
                    throw Exception("bur_cnm must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative bur_cnm not allowed.");
                cnm_ = val;
                return true;
            }
        case kwKsK: {
                if (!v.canConvert<double>())
                    throw Exception("bur_ksk must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative bur_ksk not allowed.");
                ksk_ = val;
                return true;
            }
        case kwCsK: {
                if (!v.canConvert<double>())
                    throw Exception("bur_csk must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative bur_csk not allowed.");
                csk_ = val;
                return true;
            }
        case kwKsM: {
                if (!v.canConvert<double>())
                    throw Exception("bur_ksm must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative bur_ksm not allowed.");
                ksm_ = val;
                return true;
            }
        case kwCsM: {
                if (!v.canConvert<double>())
                    throw Exception("bur_csm must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative bur_csm not allowed.");
                csm_ = val;
                return true;
            }
        case kwMode: {
                if (!v.canConvert<uint>())
                    throw Exception("bur_mode must be 0 (tensile) or 1 (no-tension).");
                uint val(v.toUInt());
                if (val >1)
                    throw Exception("bur_mode must be 0 (tensile) or 1 (no-tension).");
                bmode_ = val;
                return false;
            }
        case kwFric: {
                if (!v.canConvert<double>())
                    throw Exception("bur_fric must be a double.");
                double val(v.toDouble());
                if (val<0.0)
                    throw Exception("Negative bur_fric not allowed.");
                fric_ = val;  
                return false;
            }
        }
        return false;
    }

    bool ContactModelBurger::getPropertyReadOnly(uint i) const {
        // Returns TRUE if a property is read only or FALSE otherwise. 
        switch (i) {
        case kwF:
        case kwS:
            return true;
        default:
            break;
        }
        return false;
    }

    bool ContactModelBurger::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).
        return checkActivity(state->gap_);
    }
     
    bool ContactModelBurger::forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) {
        assert(state);

        // Current overlap
        double overlap = - 1.0*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 (trans.x()) {
                correction = -1.0*overlap / trans.x();
                if (correction < 0)
                    correction = 1.0;
            }
        }

        trans*=correction;

        if (timestep!=tdel_) { // re-calculated constants.
            tdel_ = timestep;
            // need some protection for divided by zero (k_k c_k k_m c_m = zero)
            conAn_ = A(knk_, cnk_);
            double conBn = B(knk_, cnk_);
            conB_An_ = conBn / conAn_;
            conCn_ = C(knk_, cnk_, knm_, cnm_);
            conDn_ = D(knk_, cnk_, knm_, cnm_);
            conAs_ = A(ksk_, csk_);
            double conBs = B(ksk_, csk_);
            conB_As_ = conBs / conAs_;
            conCs_ = C(ksk_, csk_, ksm_, csm_);
            conDs_ = D(ksk_, csk_, ksm_, csm_);
        }


        // normal force 
        force_.rx() = 1.0/conCn_*(overlap-u_n0_+(1.0-conB_An_)*u_nk0_-conDn_*fn0_);
        if (bmode_ && force_.x()<0.0) force_.rx() = 0.0;
        u_nk0_ = conB_An_*u_nk0_+timestep/(2.0*cnk_*conAn_)*(force_.x()+fn0_);
        u_n0_  = overlap;
        fn0_   = force_.x();

        // Calculate the shear force.
        DVect sforce(0.0);
        DVect sforce_old = force_;
        sforce_old.rx()=0.0;
        DVect v1 = trans;
        DVect v2 = u_sk_ * (1.0-conB_As_);
        DVect v3 = sforce_old * conDs_;
        sforce = (v1+v2+v3) / conCs_ * (-1.0);

        double d1 = timestep / (2.0*csk_*conAs_);
        sforce.rx() = 0.0;
        v1 = sforce + sforce_old;
        u_sk_ = u_sk_*conB_As_-v1*d1;

        // 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 = force_.x() * fric_;
            crit = max(0.0,crit);
            // This 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 (!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]);
                }
                s_ = true;
            } else {
                // Handle the slip_change event if one has been hooked up and
                // the contact was previously sliding. Sliding has ceased.  
                if (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]);
                    }
                    s_ = false;
                }
            }
        }
        
        // Set the shear components of the total force.
        for (int i=1; i<dim; ++i)
            force_.rdof(i) = sforce.dof(i);
        // This is just a sanity check to ensure, in debug mode, that the force isn't wonky. 
        assert(force_ == force_);
        return true;
    }

    void ContactModelBurger::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("burger",Qt::CaseInsensitive) == 0 && !isBonded()) {
            ContactModelBurger *oldCm = (ContactModelBurger *)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->force_.y(),oldCm->force_.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->force_.y(),oldCm->force_.z());
            }
            DVect pforce = DVect(0,tpf.x(),tpf.y());
#else
            oldSystem;
            newSystem;
            DVect pforce = DVect(0,oldCm->force_.y());
#endif
            for (int i=1; i<dim; ++i)
                force_.rdof(i) += pforce.dof(i);
            oldCm->force_ = DVect(0.0);
        }
        assert(force_ == force_);
    }

    void ContactModelBurger::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("burger",Qt::CaseInsensitive)) {
            ContactModelBurger *oldCm = (ContactModelBurger *)old;
            knk_ = oldCm->knk_;
            cnk_ = oldCm->cnk_;
            knm_ = oldCm->knm_;
            cnm_ = oldCm->cnm_;
            ksk_ = oldCm->ksk_;
            csk_ = oldCm->csk_;
            ksm_ = oldCm->ksm_;
            csm_ = oldCm->csm_;
            fric_ = oldCm->fric_;
            bmode_ = bmode_;
        }
    }

    DVect ContactModelBurger::getForce(const IContactMechanical *) const {
        DVect ret(force_);
        return ret;
    }

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

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

} // namespace cmodelsxd
// EoF

Top