Spring Network Model Implementation

See this page for the documentation of this contact model.

contactmodelrbsn.h

  1#pragma once
  2// contactmodelrbsn.h
  3
  4#include "contactmodel/src/contactmodelmechanical.h"
  5
  6#ifdef RBSN_LIB
  7#  define RBSN_EXPORT EXPORT_TAG
  8#elif defined(NO_MODEL_IMPORT)
  9#  define RBSN_EXPORT
 10#else
 11#  define RBSN_EXPORT IMPORT_TAG
 12#endif
 13
 14namespace cmodelsxd {
 15    using namespace itasca;
 16
 17    class ContactModelRBSN : public ContactModelMechanical {
 18    public:
 19        // Constructor: Set default values for contact model properties.
 20        RBSN_EXPORT ContactModelRBSN();
 21        // Destructor, called when contact is deleted: free allocated memory, etc.
 22        RBSN_EXPORT virtual ~ContactModelRBSN();
 23        // Contact model name (used as keyword for commands and FISH).
 24        QString  getName() const override { return "springnetwork"; }
 25        // The index provides a quick way to determine the type of contact model.
 26        // Each type of contact model in PFC must have a unique index; this is assigned
 27        // by PFC when the contact model is loaded. This index should be set to -1
 28        void     setIndex(int i) override { index_=i;}
 29        int      getIndex() const override {return index_;}
 30        // Contact model version number (e.g., MyModel05_1). The version number can be
 31        // accessed during the save-restore operation (within the archive method,
 32        // testing {stream.getRestoreVersion() == getMinorVersion()} to allow for 
 33        // future modifications to the contact model data structure.
 34        uint32     getMinorVersion() const override;
 35        // Copy the state information to a newly created contact model.
 36        // Provide access to state information, for use by copy method.
 37        void     copy(const ContactModel *c) override;
 38        // Provide save-restore capability for the state information.
 39        void     archive(ArchiveStream &) override; 
 40        // Enumerator for the properties.
 41        enum PropertyKeys { 
 42              kwKn=1
 43            , kwKs    
 44            , kwKrot
 45            , kwFric 
 46            , kwBMul
 47            , kwTMul
 48            , kwSNMode
 49            , kwSNF
 50            , kwSNM
 51            , kwSNS
 52            , kwSNBS
 53            , kwSNTS
 54            , kwSNRMul
 55            , kwSNRadius
 56            , kwEmod
 57            , kwKRatio
 58            , kwDpNRatio 
 59            , kwDpSRatio
 60            , kwDpMode 
 61            , kwDpF
 62            , kwSNState
 63            , kwSNTStr
 64            , kwSNSStr
 65            , kwSNCoh
 66            , kwSNFa
 67            , kwSNMCF
 68            , kwSNSig
 69            , kwSNTau
 70            , kwSNArea
 71            , kwUserArea
 72            , kwRGap
 73            , kwPForce
 74            , kwPois
 75            , kwPoisDiag
 76            , kwSnCohRes
 77            , kwSnDil
 78            , kwSnDilZ
 79            , kwSnNormDisp
 80            , kwSnShearDisp
 81            , kwSnCohDc
 82            , kwSnHeal
 83            , kwSnTenDc
 84            , kwTenTable
 85            , kwCohTable
 86            , kwTablePos
 87            , kwPorP
 88            , kwStressNorm
 89        };
 90        // Contact model property names in a comma separated list. The order corresponds with
 91        // the order of the PropertyKeys enumerator above. One can visualize any of these 
 92        // properties in PFC automatically. 
 93        QString  getProperties() const override { 
 94            return "kn"
 95                   ",ks"
 96                   ",krot"
 97                   ",fric"
 98                   ",sn_bmul"
 99                   ",sn_tmul"
100                   ",sn_mode"
101                   ",sn_force"
102                   ",sn_moment"
103                   ",sn_slip"
104                   ",sn_slipb"
105                   ",sn_slipt"
106                   ",sn_rmul"
107                   ",sn_radius"
108                   ",emod"
109                   ",kratio"
110                   ",dp_nratio"
111                   ",dp_sratio"
112                   ",dp_mode"
113                   ",dp_force"
114                   ",sn_state"
115                   ",sn_ten"
116                   ",sn_shear"
117                   ",sn_coh"
118                   ",sn_fa"
119                   ",sn_mcf"
120                   ",sn_sigma"
121                   ",sn_tau"
122                   ",sn_area"
123                   ",user_area"
124                   ",rgap"
125                   ",sn_pois_force"
126                   ",sn_pois"
127                   ",sn_non_diag"
128                   ",sn_cohres"
129                   ",sn_dil"
130                   ",sn_dilzero"
131                   ",sn_ndisp"
132                   ",sn_sdisp"
133                   ",sn_cohdc"
134                   ",sn_heal"
135                   ",sn_tendc"
136                   ",sn_tentab"
137                   ",sn_cohtab"
138                   ",sn_tabpos"
139                   ",sn_porp"
140                   ",sn_esigma"
141                ;
142        }
143        // Enumerator for the energies.
144        enum EnergyKeys { 
145            kwEStrain=1
146          , kwESlip
147          , kwEDashpot
148        };
149        // Contact model energy names in a comma separated list. The order corresponds with
150        // the order of the EnergyKeys enumerator above. 
151        QString  getEnergies() const override { 
152            return "energy-strain"
153                   ",energy-slip"
154                   ",energy-dashpot";
155        }
156        // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
157        double   getEnergy(uint32 i) const override; 
158        // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1) 
159        // returns wther or not the estrain energy is accumulated which is false).
160        bool     getEnergyAccumulate(uint32 i) const override;
161        // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
162        void     setEnergy(uint32 i,const double &d) override; // Base 1
163        // Activate the energy. This is only called if the energy tracking is enabled. 
164        void     activateEnergy()  override { if (energies_) return; energies_ = NEWC(Energies());}
165        // Returns whether or not the energy tracking has been enabled for this contact.
166        bool     getEnergyActivated() const override {return (energies_ != 0);}
167
168        // Enumerator for contact model related FISH callback events. 
169        enum FishCallEvents {
170             fActivated=0
171            ,fSlipChange
172            ,fBondBreak
173        };
174        // Contact model FISH callback event names in a comma separated list. The order corresponds with
175        // the order of the FishCallEvents enumerator above. 
176        QString  getFishCallEvents() const override { 
177            return 
178                "contact_activated"
179                ",slip_change"
180                ",bond_break"; 
181        }
182
183        // Return the specified contact model property.
184        QVariant getProperty(uint32 i,const IContact *) const override;
185        // The return value denotes whether or not the property corresponds to the global
186        // or local coordinate system (TRUE: global system, FALSE: local system). The
187        // local system is the contact-plane system (nst) defined as follows.
188        // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
189        // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
190        // and tc are unit vectors in directions of the nst axes.
191        // This is used when rendering contact model properties that are vectors.
192        bool     getPropertyGlobal(uint32 i) const override;
193        // Set the specified contact model property, ensuring that it is of the correct type
194        // and within the correct range --- if not, then throw an exception.
195        // The return value denotes whether or not the update has affected the timestep
196        // computation (by having modified the translational or rotational tangent stiffnesses).
197        // If true is returned, then the timestep will be recomputed.
198        bool     setProperty(uint32 i,const QVariant &v,IContact *) override;
199        // The return value denotes whether or not the property is read-only
200        // (TRUE: read-only, FALSE: read-write).
201        bool     getPropertyReadOnly(uint32 i) const override;
202
203        // The return value denotes whether or not the property is inheritable
204        // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
205        // the endPropertyUpdated method.
206        bool     supportsInheritance(uint32 i) const override; 
207        // Return whether or not inheritance is enabled for the specified property.
208        bool     getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
209        // Set the inheritance flag for the specified property.
210        void     setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }
211
212        // Enumerator for contact model methods.
213        enum MethodKeys { kwAssignStiffness=1, kwStiffness, kwBond, kwUnbond, kwArea, kwResetDisp};
214        // Contact model methoid names in a comma separated list. The order corresponds with
215        // the order of the MethodKeys enumerator above.  
216        QString  getMethods() const override { return "assign-stiffness,compute-stiffness,bond,unbond,area,reset-disp";}
217        // Return a comma seprated list of the contact model method arguments (base 1).
218        QString  getMethodArguments(uint32 i) const override; 
219        // Set contact model method arguments (base 1). 
220        // The return value denotes whether or not the update has affected the timestep
221        // computation (by having modified the translational or rotational tangent stiffnesses).
222        // If true is returned, then the timestep will be recomputed.
223        bool     setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override; 
224
225        // Prepare for entry into ForceDispLaw. The validate function is called when:
226        // (1) the contact is created, (2) a property of the contact that returns a true via
227        // the setProperty method has been modified and (3) when a set of cycles is executed
228        // via the {cycle N} command.
229        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
230        bool    validate(ContactModelMechanicalState *state,const double &timestep) override;
231        // The endPropertyUpdated method is called whenever a surface property (with a name
232        // that matches an inheritable contact model property name) of one of the contacting
233        // pieces is modified. This allows the contact model to update its associated
234        // properties. The return value denotes whether or not the update has affected
235        // the time step computation (by having modified the translational or rotational
236        // tangent stiffnesses). If true is returned, then the time step will be recomputed.  
237        bool    endPropertyUpdated(const QString &name,const IContactMechanical *c) override;
238        // The forceDisplacementLaw function is called during each cycle. Given the relative
239        // motion of the two contacting pieces (via
240        //   state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
241        //   state->relativeAngularIncrement_       (Dtt, Dtbs, Dtbt)
242        //     Ddn  : relative normal-displacement increment, Ddn > 0 is opening
243        //     Ddss : relative  shear-displacement increment (s-axis component)
244        //     Ddst : relative  shear-displacement increment (t-axis component)
245        //     Dtt  : relative twist-rotation increment
246        //     Dtbs : relative  bend-rotation increment (s-axis component)
247        //     Dtbt : relative  bend-rotation increment (t-axis component)
248        //       The relative displacement and rotation increments:
249        //         Dd = Ddn*nc + Ddss*sc + Ddst*tc
250        //         Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
251        //       where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
252        //       [see {Table 1: Contact State Variables} in PFC Model Components:
253        //       Contacts and Contact Models: Contact Resolution]
254        // ) and the contact properties, this function must update the contact force and
255        // moment.
256        //   The force_ is acting on piece 2, and is expressed in the local coordinate system
257        //   (defined in getPropertyGlobal) such that the first component positive denotes
258        //   compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
259        //   in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to 
260        //   get the total moment. 
261        // The return value indicates the contact activity status (TRUE: active, FALSE:
262        // inactive) during the next cycle.
263        // Additional information:
264        //   * If state->activated() is true, then the contact has just become active (it was
265        //     inactive during the previous time step).
266        //   * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
267        //     the forceDispLaw handle the case of {state->canFail_ == true}.
268        bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) override;
269        bool    thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
270        
271        // The getEffectiveXStiffness functions return the translational and rotational
272        // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
273        // the translational tangent shear stiffness is zero (but this stiffness reduction
274        // is typically ignored when computing a stable time step). If the contact model
275        // includes a dashpot, then the translational stiffnesses must be increased (see
276        // Potyondy (2009)).
277        //   [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
278        //   Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
279        //   December 7, 2009.]
280        DVect2  getEffectiveTranslationalStiffness() const override { return effectiveTranslationalStiffness_; }
281        DAVect  getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness_;}
282
283        // Return a new instance of the contact model. This is used in the CMAT
284        // when a new contact is created. 
285        ContactModelRBSN *clone() const override { return NEWC(ContactModelRBSN()); }
286        // The getActivityDistance function is called by the contact-resolution logic when
287        // the CMAT is modified. Return value is the activity distance used by the
288        // checkActivity function.
289        double              getActivityDistance() const override {return rgap_;}
290        // The isOKToDelete function is called by the contact-resolution logic when...
291        // Return value indicates whether or not the contact may be deleted.
292        // If TRUE, then the contact may be deleted when it is inactive.
293        // If FALSE, then the contact may not be deleted (under any condition).
294        bool                isOKToDelete() const override { return !isBonded(); }
295        // Zero the forces and moments stored in the contact model. This function is called
296        // when the contact becomes inactive.
297        void                resetForcesAndMoments() override { 
298            sn_F_ = DVect(0.0);
299            fictForce_ = DVect(0.0);
300            sn_M_ = DAVect(0.0);
301            dp_F(DVect(0.0)); 
302            if (energies_) {
303                energies_->estrain_ = 0.0;
304            }
305        }
306        void     setForce(const DVect &v,IContact *c) override;
307        void     setArea(const double &d) override { userArea_ = d; }
308        double   getArea() const override { return userArea_; }
309
310        // The checkActivity function is called by the contact-resolution logic when...
311        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
312        bool     checkActivity(const double &gap) override { return  gap <= rgap_ || isBonded();}
313
314        // Returns the sliding state (FALSE is returned if not implemented).
315        bool     isSliding() const override { return sn_S_; }
316        // Returns the bonding state (FALSE is returned if not implemented).
317        bool     isBonded() const override { return sn_state_ >= 3; }
318        void     unbond() override { sn_state_ = 0; }
319
320        // Both of these methods are called only for contacts with facets where the wall 
321        // resolution scheme is set the full. In such cases one might wish to propagate 
322        // contact state information (e.g., shear force) from one active contact to another. 
323        // See the Faceted Wall section in the documentation. 
324        void     propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
325        void     setNonForcePropsFrom(IContactModel *oldCM) override;
326           
327        /// Return the total force that the contact model holds.
328        DVect    getForce(const IContactMechanical *) const override;
329
330        /// Return the total moment on 1 that the contact model holds
331        DAVect   getMomentOn1(const IContactMechanical *) const override;
332
333        /// Return the total moment on 1 that the contact model holds
334        DAVect   getMomentOn2(const IContactMechanical *) const override;
335
336        // Methods to get and set properties. 
337        double         sn_Ten() const { return tenTable_[0].x(); }
338        void           sn_Ten(const double &d) { tenTable_[0].rx() = d; }
339        double         sn_Coh() const { return cohTable_[0].x(); }
340        void           sn_Coh(const double &d) { cohTable_[0].rx() = d; }
341        void           sn_MCF(const double &d) { sn_mcf_=d;}
342        double         sn_cohdc() const           {return cohTable_.back().y(); }
343        double         sn_tendc() const           {return tenTable_.back().y(); }
344
345        bool     hasDamping() const {return dpProps_ ? true : false;}
346        double   dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
347        void     dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
348        double   dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
349        void     dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
350        int      dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
351        void     dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
352        DVect    dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
353        void     dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
354
355        bool    hasEnergies() const {return energies_ ? true:false;}
356        double  estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
357        void    estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
358        double  eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
359        void    eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
360        double  edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
361        void    edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
362
363        uint32 inheritanceField() const {return inheritanceField_;}
364        void inheritanceField(uint32 i) {inheritanceField_ = i;}
365
366        const DVect2 & effectiveTranslationalStiffness()  const          {return effectiveTranslationalStiffness_;}
367        void           effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
368        const DAVect & effectiveRotationalStiffness()  const             {return effectiveRotationalStiffness_;}
369        void           effectiveRotationalStiffness(const DAVect &v )    {effectiveRotationalStiffness_=v;}
370
371    private:
372        // Index - used internally by PFC. Should be set to -1 in the cpp file. 
373        static int index_;
374
375        bool  FDLawBonded(ContactModelMechanicalState *state, const double &timestep);
376        bool  FDLawUnBonded(ContactModelMechanicalState *state, const double &timestep);
377
378        // Structure to compute stiffness
379        struct StiffData {
380            DVect2 trans_ = DVect2(0.0);
381            DAVect ang_   = DAVect(0.0);
382            double reff_ = 0.0;
383        };
384
385        // Structure to store the energies. 
386        struct Energies {
387            double estrain_  = 0.0;   // elastic energy  
388            double eslip_    = 0.0;   // work dissipated by friction 
389            double edashpot_ = 0.0;   // work dissipated by dashpots
390        };
391
392        // Structure to store dashpot quantities. 
393        struct dpProps {
394            double dp_nratio_ = 0.0;         // normal viscous critical damping ratio
395            double dp_sratio_ = 0.0;         // shear  viscous critical damping ratio
396            int    dp_mode_   = 0;           // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
397            DVect  dp_F_      = DVect(0.0);  // Force in the dashpots
398        };
399
400        bool   updateKn(const IContactMechanical *con);
401        bool   updateKs(const IContactMechanical *con);
402        bool   updateFric(const IContactMechanical *con);
403
404        StiffData computeStiffData(ContactModelMechanicalState *state) const;
405        DVect3    computeGeomData(const IContactMechanical *c) const;
406        DVect2    SMax(const IContactMechanical *con) const; // Maximum stress (tensile,shear) at bond periphery
407        double    shearStrength(const double &pbArea) const;      // Bond shear strength
408        double    strainEnergy(double kn, double ks, double kb, double kt) const;
409
410        void      updateStiffness(ContactModelMechanicalState *state);
411
412        // Contact model inheritance fields.
413        uint32 inheritanceField_;
414
415        // Effective translational stiffness.
416        DVect2  effectiveTranslationalStiffness_ = DVect2(0.0); 
417        DAVect  effectiveRotationalStiffness_ = DAVect(0.0);      // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)
418
419        // linear model properties
420        DVect       fictForce_ = DVect(0.0);// Ficticous force to be added 
421        DVect       sn_F_ = DVect(0.0);     // Force carried in the model
422        DVect2      sn_sdisp_ = DVect2(0.0); // Accumulated total shear displacement
423                                            // The x component holds the current slip
424        DAVect      sn_M_ = DAVect(0.0);    // moment (bending + twisting in 3D)         
425        DAVect      kRot_ = DAVect(0.0);    // Translational degrees of freedom
426        double      kTran_ = 0.0;           // Translational degrees of freedom
427        double      kRatio_ = 1.0;          // Ratio of normal to shear stiffness
428        double      E_ = 0.0;               // Young's modulus
429        double      poisson_ = 0.0;         // Poisson ratio
430        double      fric_ = 0.0;            // Coulomb friction coefficient
431        double      sn_bmul_ = 1.0;         // Bending friction multiplier
432        double      sn_tmul_ = 1.0;         // Twisting friction  multiplier
433        double      sn_rmul_ = 1.0;         // Radius multiplier
434        double      userArea_ = 0.0;        // Area as specified by the user 
435        double      rgap_ = 0.0;            // Reference gap
436        double      sn_fa_  = 0.0;          // friction angle (stored as tan(dDegrad*fa))
437        double      sn_mcf_ = 1.0;          // moment contribution factor
438        double      sn_dil_ = 0.0;          // Dilation (stored as tan(dDegrad*dil))
439        double      sn_dilzero_ = 0.0;      // Dilation zero
440        double      transTen_ = 0.0;        // Force for transition from tensile to compression
441        double      sn_elong_ = 0.0;        // Elongation (or normal displacement since softening)
442        double      sn_ndisp_ = 0.0;        // Accumulated normal displacement
443        double      sn_cohres_ = 0.0;       // Residual cohesion
444        double      sn_por_ = 0.0;          // Pore Pressure
445        uint        sn_mode_ = 0;           // specifies absolute (0) or incremental (1) behavior for the the normal force
446        uint        sn_state_  = 0;         // bond mode - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (B)
447        int         sn_tabPos_ = 0;         // Position in the table for query
448        bool        poisOffDiag_ = false;   // Add the off diagonal terms 
449        bool        sn_S_ = false;          // The current slip state
450        bool        sn_BS_ = false;         // The bending  slip state
451        bool        sn_TS_ = false;         // The twisting slip state
452        bool        forceSet_ = false;      // About setting the force
453        bool        sn_heal_ = false;       // Healing behavior
454
455        std::vector<DVect2> tenTable_ = { DVect2(0) };     
456        std::vector<DVect2> cohTable_ = { DVect2(0) };
457
458        dpProps *   dpProps_ = nullptr;     // The viscous properties
459
460        Energies *   energies_ = nullptr;   // The energies
461
462    };
463} // namespace cmodelsxd
464// EoF

Top

contactmodelrbsn.cpp

   1// contactmodelrbsn.cpp
   2#include "contactmodelrbsn.h"
   3
   4#include "module/interface/icontactmechanical.h"
   5#include "module/interface/icontact.h"
   6#include "module/interface/ipiece.h"
   7#include "module/interface/ibody.h"
   8#include "module/interface/ifishcalllist.h"
   9
  10#include "utility/src/tptr.h"
  11#include "shared/src/mathutil.h"
  12
  13#include "kernel/interface/iprogram.h"
  14#include "module/interface/icontactthermal.h"
  15#include "contactmodel/src/contactmodelthermal.h"
  16#include "../version.txt"
  17#include "fish/src/parameter.h"
  18
  19#ifdef RBSN_LIB
  20    int __stdcall DllMain(void *,unsigned, void *) {
  21        return 1;
  22    }
  23
  24    extern "C" EXPORT_TAG const char *getName() {
  25#if DIM==3
  26        return "contactmodelmechanical3drbsn";
  27#else
  28        return "contactmodelmechanical2drbsn";
  29#endif
  30    }
  31
  32    extern "C" EXPORT_TAG unsigned getMajorVersion() {
  33        return MAJOR_VERSION;
  34    }
  35
  36    extern "C" EXPORT_TAG unsigned getMinorVersion() {
  37        return MINOR_VERSION;
  38    }
  39
  40    extern "C" EXPORT_TAG void *createInstance() {
  41        cmodelsxd::ContactModelRBSN *m = NEWC(cmodelsxd::ContactModelRBSN());
  42        return (void *)m;
  43    }
  44#endif
  45
  46namespace cmodelsxd {
  47    static const uint32 KnMask      = 0x00000002; // Base 1!
  48    static const uint32 KsMask      = 0x00000004;
  49    static const uint32 FricMask    = 0x00000008;
  50
  51    using namespace itasca;
  52
  53    int ContactModelRBSN::index_ = -1;
  54    uint32 ContactModelRBSN::getMinorVersion() const { return MINOR_VERSION; }
  55
  56    ContactModelRBSN::ContactModelRBSN() : inheritanceField_(KnMask|KsMask|FricMask) {
  57    }
  58
  59    ContactModelRBSN::~ContactModelRBSN() {
  60        // Make sure to clean up after yourself!
  61        if (dpProps_)
  62            delete dpProps_;
  63        if (energies_)
  64            delete energies_;
  65    }
  66
  67    void ContactModelRBSN::archive(ArchiveStream &stream) {
  68        if (stream.getRestoreVersion() > 4) {
  69            // New version
  70            stream & fictForce_;
  71            stream & sn_F_;
  72            stream & sn_sdisp_;
  73            stream & sn_M_;
  74            stream & kRot_;
  75            stream & kTran_;
  76            stream & kRatio_;
  77            stream & E_;
  78            stream & poisson_;
  79            stream & fric_;
  80            stream & sn_bmul_;
  81            stream & sn_tmul_;
  82            stream & sn_rmul_;
  83            stream & userArea_;
  84            stream & rgap_;
  85            stream & sn_fa_;
  86            stream & sn_mcf_;
  87            stream & sn_dil_;
  88            stream & sn_dilzero_;
  89            stream & transTen_;
  90            stream & sn_elong_;
  91            stream & sn_ndisp_;
  92            stream & sn_mode_;
  93            stream & sn_state_;
  94            stream & poisOffDiag_;
  95            stream & sn_S_;
  96            stream & sn_BS_;
  97            stream & sn_TS_;
  98            stream & forceSet_;
  99            stream & sn_heal_;
 100            stream & tenTable_;
 101            stream & cohTable_;
 102            stream & inheritanceField_;
 103            stream & effectiveTranslationalStiffness_;
 104            stream & effectiveRotationalStiffness_;
 105            if (stream.getArchiveState()==ArchiveStream::Save) {
 106                bool b = false;
 107                if (dpProps_) {
 108                    b = true;
 109                    stream & b;
 110                    stream & dpProps_->dp_nratio_;
 111                    stream & dpProps_->dp_sratio_;
 112                    stream & dpProps_->dp_mode_;
 113                    stream & dpProps_->dp_F_;
 114                }
 115                else
 116                    stream & b;
 117                b = false;
 118                if (energies_) {
 119                    b = true;
 120                    stream & b;
 121                    stream & energies_->estrain_;
 122                    stream & energies_->eslip_;
 123                    stream & energies_->edashpot_;
 124                }
 125                else
 126                    stream & b;
 127            } else {
 128                bool b(false);
 129                stream & b;
 130                if (b) {
 131                    if (!dpProps_)
 132                        dpProps_ = NEWC(dpProps());
 133                    stream & dpProps_->dp_nratio_;
 134                    stream & dpProps_->dp_sratio_;
 135                    stream & dpProps_->dp_mode_;
 136                    stream & dpProps_->dp_F_;
 137                }
 138                stream & b;
 139                if (b) {
 140                    if (!energies_)
 141                        energies_ = NEWC(Energies());
 142                    stream & energies_->estrain_;
 143                    stream & energies_->eslip_;
 144                    stream & energies_->edashpot_;
 145                }
 146            }
 147
 148            if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 5) {
 149                stream & sn_tabPos_;
 150                stream & sn_cohres_;
 151            }
 152
 153            if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 6)
 154                stream & sn_por_;
 155                
 156        } else {
 157            // Backward compatibility
 158            stream & kTran_;
 159            stream & E_;
 160            stream & kRot_;
 161            stream & fictForce_;
 162            stream & poisson_;
 163            stream & fric_;
 164            stream & sn_mode_;
 165            stream & sn_F_;
 166            stream & sn_M_;
 167            stream & sn_S_;
 168            stream & sn_BS_;
 169            stream & sn_TS_;
 170            stream & sn_rmul_;
 171
 172            bool b(false);
 173            stream & b;
 174            if (b) {
 175                if (!dpProps_)
 176                    dpProps_ = NEWC(dpProps());
 177                stream & dpProps_->dp_nratio_;
 178                stream & dpProps_->dp_sratio_;
 179                stream & dpProps_->dp_mode_;
 180                stream & dpProps_->dp_F_;
 181            }
 182            stream & b;
 183            if (b) {
 184                if (!energies_)
 185                    energies_ = NEWC(Energies());
 186                stream & energies_->estrain_;
 187                stream & energies_->eslip_;
 188                stream & energies_->edashpot_;
 189            }
 190            stream & b;
 191            if (b) {
 192                int vi = 0;
 193                stream & vi;
 194                sn_state_ = abs(vi);
 195                double val = 0.0;
 196                stream & val;
 197                tenTable_[0].rx() = val;
 198                val = 0.0;
 199                stream & val;
 200                cohTable_[0].rx() = val;
 201                stream & sn_fa_;
 202                stream & sn_mcf_;
 203                stream & val;
 204                stream & val;
 205                stream & val;
 206                stream & val;
 207                Quat q;
 208                stream & q;
 209                stream & val;
 210                stream & val;
 211            }
 212
 213            stream & inheritanceField_;
 214            stream & effectiveTranslationalStiffness_;
 215            stream & effectiveRotationalStiffness_;
 216
 217            stream & sn_bmul_;
 218            stream & sn_tmul_;
 219            stream & userArea_;
 220            stream & rgap_;
 221
 222            if (stream.getRestoreVersion() > 1)
 223                stream & kRatio_;
 224
 225            if (stream.getRestoreVersion() > 2) {
 226                uint32 v;
 227                stream & v;
 228                poisOffDiag_ = v == 0 ? false : true;
 229            }
 230
 231            if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 3)
 232                stream & forceSet_;
 233
 234        }
 235    }
 236
 237    void ContactModelRBSN::copy(const ContactModel *cm) {
 238        // Copy all of the contact model properties. Used in the CMAT
 239        // when a new contact is created.
 240        ContactModelMechanical::copy(cm);
 241        const ContactModelRBSN *in = dynamic_cast<const ContactModelRBSN*>(cm);
 242        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
 243        fictForce_ = in->fictForce_;
 244        sn_F_ = in->sn_F_;
 245        sn_sdisp_ = in->sn_sdisp_;
 246        sn_M_ = in->sn_M_;
 247        kRot_ = in->kRot_;
 248        kTran_ = in->kTran_;
 249        kRatio_ = in->kRatio_;
 250        E_ = in->E_;
 251        poisson_ = in->poisson_;
 252        fric_ = in->fric_;
 253        sn_bmul_ = in->sn_bmul_;
 254        sn_tmul_ = in->sn_tmul_;
 255        sn_rmul_ = in->sn_rmul_;
 256        userArea_ = in->userArea_;
 257        rgap_ = in->rgap_;
 258        sn_fa_ = in->sn_fa_;
 259        sn_mcf_ = in->sn_mcf_;
 260        sn_dil_ = in->sn_dil_;
 261        sn_dilzero_ = in->sn_dilzero_;
 262        transTen_ = in->transTen_;
 263        sn_elong_ = in->sn_elong_;
 264        sn_ndisp_ = in->sn_ndisp_;
 265        sn_mode_ = in->sn_mode_;
 266        sn_state_ = in->sn_state_;
 267        poisOffDiag_ = in->poisOffDiag_;
 268        sn_S_ = in->sn_S_;
 269        sn_BS_ = in->sn_BS_;
 270        sn_TS_ = in->sn_TS_;
 271        forceSet_ = in->forceSet_;
 272        sn_heal_ = in->sn_heal_;
 273        tenTable_ = in->tenTable_;
 274        cohTable_ = in->cohTable_;
 275        sn_por_ = in->sn_por_;
 276        if (in->hasDamping()) {
 277            if (!dpProps_)
 278                dpProps_ = NEWC(dpProps());
 279            dp_nratio(in->dp_nratio());
 280            dp_sratio(in->dp_sratio());
 281            dp_mode(in->dp_mode());
 282            dp_F(in->dp_F());
 283        }
 284        if (in->hasEnergies()) {
 285            if (!energies_)
 286                energies_ = NEWC(Energies());
 287            estrain(in->estrain());
 288            eslip(in->eslip());
 289            edashpot(in->edashpot());
 290        }
 291        inheritanceField(in->inheritanceField());
 292        effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
 293        effectiveRotationalStiffness(in->effectiveRotationalStiffness());
 294    }
 295
 296
 297    QVariant ContactModelRBSN::getProperty(uint32 i,const IContact *con) const {
 298        // Return the property. The IContact pointer is provided so that
 299        // more complicated properties, depending on contact characteristics,
 300        // can be calcualted.
 301        QVariant var;
 302        switch (i) {
 303        case kwKn:       return kTran_;
 304        case kwKs:       return kTran_ / kRatio_;
 305        case kwKrot:     var.setValue(kRot_); return var;
 306        case kwFric:     return fric_;
 307        case kwBMul:     return sn_bmul_;
 308        case kwTMul:     return sn_tmul_;
 309        case kwSNMode:   return sn_mode_;
 310        case kwSNF:      var.setValue(sn_F_); return var;
 311        case kwSNM:      var.setValue(sn_M_); return var;
 312        case kwSNS:      return sn_S_;
 313        case kwSNBS:     return sn_BS_;
 314        case kwSNTS:     return sn_TS_;
 315        case kwPoisDiag: return poisOffDiag_ == false ? 0 : 1;
 316        case kwSNRMul:   return sn_rmul_;
 317        case kwSNRadius: {
 318            const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 319            if (!c) return 0.0;
 320            double Cmax1 = c->getEnd1Curvature().y();
 321            double Cmax2 = c->getEnd2Curvature().y();
 322            if (!userArea_)
 323                return sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
 324            else {
 325#ifdef THREED
 326                double rad = std::sqrt(userArea_ / dPi);
 327#else
 328                double rad = userArea_ / 2.0;
 329#endif
 330                return rad;
 331            }
 332
 333        }
 334        case kwEmod:      return E_;
 335        case kwKRatio:    return 1.0;
 336        case kwDpNRatio:  return dpProps_ ? dpProps_->dp_nratio_ : 0;
 337        case kwDpSRatio:  return dpProps_ ? dpProps_->dp_sratio_ : 0;
 338        case kwDpMode:    return dpProps_ ? dpProps_->dp_mode_ : 0;
 339        case kwDpF: {
 340                dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
 341                return var;
 342            }
 343        case kwSNState:     return sn_state_;
 344        case kwSNTStr:      return sn_Ten();
 345        case kwSNSStr: {
 346            const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 347            if (!c) return 0.0;
 348            double area = computeGeomData(c).x();
 349            return shearStrength(area);
 350        }
 351        case kwSNCoh:       return cohTable_[0].x();
 352        case kwSNFa:        return std::atan(sn_fa_)/dDegrad;
 353        case kwSNMCF:       return sn_mcf_;
 354        case kwSNSig: {
 355            const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 356            if (!c) return 0.0;
 357            return SMax(c).x();
 358        }
 359        case kwSNTau: {
 360            const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 361            if (!c) return 0.0;
 362            return SMax(c).y();
 363        }
 364        case kwSNArea: {
 365                if (userArea_) return userArea_;
 366                const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 367                if (!c)
 368                    return 0.0;
 369                return computeGeomData(c).x();
 370            }
 371        case kwUserArea:
 372            return userArea_;
 373        case kwRGap:
 374            return rgap_;
 375        case kwPForce:  var.setValue(fictForce_); return var;
 376        case kwPois: return poisson_;
 377        case kwSnCohRes :   return sn_cohres_;
 378        //case kwSnTenRes :   return sn_tenres();
 379        case kwSnDil :      return std::atan(sn_dil_)/dDegrad;
 380        case kwSnDilZ :     return sn_dilzero_;
 381        case kwSnNormDisp:  return sn_ndisp_;
 382        case kwSnShearDisp:  var.setValue(sn_sdisp_); return var;
 383        case kwSnCohDc :    return sn_cohdc();
 384        case kwSnTenDc :    return sn_tendc();
 385        case kwSnHeal :     return sn_heal_;
 386        case kwTenTable:
 387            if (sn_tabPos_ < tenTable_.size())
 388                if (sn_tabPos_ == 0)
 389                    var.setValue(DVect2(1,0));
 390                else
 391                    var.setValue(tenTable_[sn_tabPos_]);
 392            else
 393                var.setValue(DVect2(0,0));
 394            return var;
 395        case kwCohTable:
 396            if (sn_tabPos_ < cohTable_.size())
 397                if (sn_tabPos_ == 0)
 398                    var.setValue(DVect2(1,0));
 399                else
 400                    var.setValue(cohTable_[sn_tabPos_]);
 401            else
 402                var.setValue(DVect2(0,0));
 403            return var;
 404        case kwTablePos: return sn_tabPos_+1;
 405        case kwPorP: return sn_por_;
 406        case kwStressNorm: {
 407                const IContactMechanical* c(convert_getcast<IContactMechanical>(con));
 408                if (!c)
 409                    return 0.0;
 410                return SMax(c).x() + sn_por_;
 411            }
 412        }
 413        assert(0);
 414        return QVariant();
 415    }
 416
 417    bool ContactModelRBSN::getPropertyGlobal(uint32 i) const {
 418        // Returns whether or not a property is held in the global axis system (TRUE)
 419        // or the local system (FALSE). Used by the plotting logic.
 420        switch (i) {
 421        case kwSNF:
 422        case kwSNM:
 423        case kwDpF:
 424            return false;
 425        }
 426        return true;
 427    }
 428
 429    bool ContactModelRBSN::setProperty(uint32 i,const QVariant &v,IContact *) {
 430        // Set a contact model property. Return value indicates that the timestep
 431        // should be recalculated.
 432        dpProps dp;
 433        switch (i) {
 434        case kwKn: {
 435                if (!v.canConvert<double>())
 436                    throw Exception("kn must be a double.");
 437                double val(v.toDouble());
 438                if (val<0.0)
 439                    throw Exception("Negative kn not allowed.");
 440                kTran_ = val;
 441                return true;
 442            }
 443        case kwKs: {
 444                if (!v.canConvert<double>())
 445                    throw Exception("ks must be a double.");
 446                double val(v.toDouble());
 447                if (val<0.0)
 448                    throw Exception("Negative ks not allowed.");
 449                if (!kTran_)
 450                    kTran_ = val;
 451                else
 452                    kRatio_ = kTran_ / val;
 453                return true;
 454            }
 455        case kwKrot: {
 456                DAVect val(0.0);
 457#ifdef TWOD
 458                if (!v.canConvert<DAVect>() && !v.canConvert<double>())
 459                    throw Exception("krot must be an angular vector.");
 460                if (v.canConvert<DAVect>())
 461                    val = DAVect(v.value<DAVect>());
 462                else
 463                    val = DAVect(v.toDouble());
 464#else
 465                if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
 466                    throw Exception("krot must be an angular vector.");
 467                if (v.canConvert<DAVect>())
 468                    val = DAVect(v.value<DAVect>());
 469                else
 470                    val = DAVect(v.value<DVect>());
 471#endif
 472                kRot_ = val;
 473                return false;
 474            }
 475        case kwFric: {
 476                if (!v.canConvert<double>())
 477                    throw Exception("fric must be a double.");
 478                double val(v.toDouble());
 479                if (val<0.0)
 480                    throw Exception("Negative fric not allowed.");
 481                fric_ = val;
 482                //if (!sn_fa_)
 483                //    sn_fa_ = fric_;
 484                return false;
 485            }
 486        case kwBMul: {
 487                if (!v.canConvert<double>())
 488                    throw Exception("sn_bmul must be a double.");
 489                double val(v.toDouble());
 490                if (val<0.0)
 491                    throw Exception("Negative sn_bmul not allowed.");
 492                sn_bmul_ = val;
 493                return false;
 494            }
 495        case kwTMul: {
 496                if (!v.canConvert<double>())
 497                    throw Exception("sn_tmul must be a double.");
 498                double val(v.toDouble());
 499                if (val<0.0)
 500                    throw Exception("Negative st_bmul not allowed.");
 501                sn_tmul_ = val;
 502                return false;
 503            }
 504        case kwSNMode: {
 505                if (!v.canConvert<uint32>())
 506                    throw Exception("sn_mode must be 0 (absolute) or 1 (incremental).");
 507                double val(v.toUInt());
 508                if (val>1)
 509                    throw Exception("sn_mode must be 0 (absolute) or 1 (incremental).");
 510                sn_mode_ = val;
 511                return false;
 512            }
 513        case kwSNRMul: {
 514                if (!v.canConvert<double>())
 515                    throw Exception("rmul must be a double.");
 516                double val(v.toDouble());
 517                if (val<0.0)
 518                    throw Exception("Negative rmul not allowed.");
 519                sn_rmul_ = val;
 520                return false;
 521            }
 522        case kwSNF: {
 523                if (!v.canConvert<DVect>())
 524                    throw Exception("sn_force must be a vector.");
 525                DVect val(v.value<DVect>());
 526                sn_F_ = val;
 527                return false;
 528            }
 529        case kwSNM: {
 530                DAVect val(0.0);
 531#ifdef TWOD
 532                if (!v.canConvert<DAVect>() && !v.canConvert<double>())
 533                    throw Exception("res_moment must be an angular vector.");
 534                if (v.canConvert<DAVect>())
 535                    val = DAVect(v.value<DAVect>());
 536                else
 537                    val = DAVect(v.toDouble());
 538#else
 539                if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
 540                    throw Exception("res_moment must be an angular vector.");
 541                if (v.canConvert<DAVect>())
 542                    val = DAVect(v.value<DAVect>());
 543                else
 544                    val = DAVect(v.value<DVect>());
 545#endif
 546                sn_M_ = val;
 547                return false;
 548            }
 549        case kwDpNRatio: {
 550                if (!v.canConvert<double>())
 551                    throw Exception("dp_nratio must be a double.");
 552                double val(v.toDouble());
 553                if (val<0.0)
 554                    throw Exception("Negative dp_nratio not allowed.");
 555                if (val == 0.0 && !dpProps_)
 556                    return false;
 557                if (!dpProps_)
 558                    dpProps_ = NEWC(dpProps());
 559                dpProps_->dp_nratio_ = val;
 560                return true;
 561            }
 562        case kwDpSRatio: {
 563                if (!v.canConvert<double>())
 564                    throw Exception("dp_sratio must be a double.");
 565                double val(v.toDouble());
 566                if (val<0.0)
 567                    throw Exception("Negative dp_sratio not allowed.");
 568                if (val == 0.0 && !dpProps_)
 569                    return false;
 570                if (!dpProps_)
 571                    dpProps_ = NEWC(dpProps());
 572                dpProps_->dp_sratio_ = val;
 573                return true;
 574            }
 575        case kwDpMode: {
 576                if (!v.canConvert<int>())
 577                    throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
 578                int val(v.toInt());
 579                if (val == 0 && !dpProps_)
 580                    return false;
 581                if (val < 0 || val > 3)
 582                    throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
 583                if (!dpProps_)
 584                    dpProps_ = NEWC(dpProps());
 585                dpProps_->dp_mode_ = val;
 586                return false;
 587            }
 588        case kwDpF: {
 589                if (!v.canConvert<DVect>())
 590                    throw Exception("dp_force must be a vector.");
 591                DVect val(v.value<DVect>());
 592                if (val.fsum() == 0.0 && !dpProps_)
 593                    return false;
 594                if (!dpProps_)
 595                    dpProps_ = NEWC(dpProps());
 596                dpProps_->dp_F_ = val;
 597                return false;
 598            }
 599        case kwSNTStr: {
 600                if (!v.canConvert<double>())
 601                    throw Exception("sn_ten must be a double.");
 602                double val(v.toDouble());
 603                if (val < 0.0)
 604                    throw Exception("Negative sn_ten not allowed.");
 605                tenTable_[0].rx() = val;
 606                return false;
 607            }
 608        case kwSNCoh: {
 609                if (!v.canConvert<double>())
 610                    throw Exception("sn_coh must be a double.");
 611                double val(v.toDouble());
 612                if (val<0.0)
 613                    throw Exception("Negative pb_coh not allowed.");
 614                cohTable_[0].rx() = val;
 615                return false;
 616            }
 617        case kwSNFa: {
 618                if (!v.canConvert<double>())
 619                    throw Exception("sn_fa must be a double.");
 620                double val(v.toDouble());
 621                if (val<0.0)
 622                    throw Exception("Negative sn_fa not allowed.");
 623                if (val >= 90.0)
 624                    throw Exception("sn_fa must be lower than 90.0 degrees.");
 625                sn_fa_ = std::tan(val*dDegrad);
 626                if (!fric_)
 627                    fric_ = sn_fa_;
 628                return false;
 629            }
 630        case kwSNMCF: {
 631                if (!v.canConvert<double>())
 632                    throw Exception("sn_mcf must be a double.");
 633                double val(v.toDouble());
 634                if (val<0.0)
 635                    throw Exception("Negative sn_mcf not allowed.");
 636                if (val > 1.0)
 637                    throw Exception("sn_mcf must be lower or equal to 1.0.");
 638                sn_mcf_ = val;
 639                return false;
 640            }
 641        case kwSNArea:
 642        case kwUserArea: {
 643                if (!v.canConvert<double>())
 644                    throw Exception("area must be a double.");
 645                double val(v.toDouble());
 646                if (val < 0.0)
 647                    throw Exception("Negative area not allowed.");
 648                if (userArea_ && val) {
 649                    double rat = userArea_ / val;
 650                    kTran_ *=  rat;
 651                    kRot_ *= rat;
 652                }
 653                userArea_ = val;
 654                return true;
 655            }
 656        case kwRGap: {
 657                if (!v.canConvert<double>())
 658                    throw Exception("Reference gap must be a double.");
 659                double val(v.toDouble());
 660                rgap_ = val;
 661                return false;
 662            }
 663        case kwPois: {
 664                if (!v.canConvert<double>())
 665                    throw Exception("Reference poisson must be a double.");
 666                double val(v.toDouble());
 667                poisson_ = val;
 668                return false;
 669            }
 670        case kwPoisDiag: {
 671                if (!v.canConvert<uint32>())
 672                    throw Exception("Reference diagonal must be an integer.");
 673                uint32 b(v.toUInt());
 674                if (b > 1)
 675                    throw Exception("diagonal must be 0 (diagonal terms only) or 1 (all terms).");
 676                poisOffDiag_ = b == 0 ? false : true;
 677                return false;
 678            }
 679        case kwSnCohRes: {
 680                bool ok;
 681                double val(v.toDouble(&ok));
 682                if (!ok || val<0.0)
 683                    throw Exception("sn_cohres must be a positive double.");
 684                sn_cohres_ = val;
 685                return false;
 686            }
 687        case kwSnDil: {
 688                bool ok;
 689                double val(v.toDouble(&ok));
 690                if (!ok || val<0.0)
 691                    throw Exception("sn_dil must be a positive double.");
 692                sn_dil_ = std::tan(val*dDegrad);
 693                return false;
 694            }
 695        case kwSnDilZ: {
 696                bool ok;
 697                double val(v.toDouble(&ok));
 698                if (!ok || val<0.0)
 699                    throw Exception("sn_dil_zero must be a positive double.");
 700                sn_dilzero_ = val;
 701                return false;
 702            }
 703        case kwSnNormDisp: {
 704                bool ok;
 705                double val(v.toDouble(&ok));
 706                if (!ok)
 707                    throw Exception("sn_ndisp must be a positive double.");
 708                sn_ndisp_ = val;
 709                return false;
 710            }
 711        case kwSnShearDisp: {
 712                if (!v.canConvert<DVect2>())
 713                    throw Exception("sn_sdisp must be a vector.");
 714                DVect2 val(v.value<DVect2>());
 715                sn_sdisp_ = val;
 716                return false;
 717            }
 718        case kwSnCohDc: {
 719                bool ok;
 720                double val(v.toDouble(&ok));
 721                if (!ok || val<0.0)
 722                    throw Exception("sn_cohdc must be a positive double.");
 723                if (cohTable_.size() == 1)
 724                    cohTable_.push_back(DVect2(0,val));
 725                else {
 726                    cohTable_.back().ry() = val;
 727                    std::sort(cohTable_.begin(),cohTable_.end(), [](const DVect2& a, const DVect2& b) { return a.y() < b.y(); } );
 728                    while (cohTable_.back().y() > val)
 729                        cohTable_.pop_back();
 730                }
 731                if (sn_state_ == 0)
 732                    sn_state_ = 6;
 733                return false;
 734            }
 735        case kwSnTenDc: {
 736                bool ok;
 737                double val(v.toDouble(&ok));
 738                if (!ok || val<0.0)
 739                    throw Exception("sn_tendc must be a positive double.");
 740                if (tenTable_.size() == 1)
 741                    tenTable_.push_back(DVect2(0,val));
 742                else {
 743                    tenTable_.back().ry() = val;
 744                    std::sort(tenTable_.begin(),tenTable_.end(), [](const DVect2& a, const DVect2& b) { return a.y() < b.y(); } );
 745                    while (tenTable_.back().y() > val)
 746                        tenTable_.pop_back();
 747                }
 748                return false;
 749            }
 750        case kwSnHeal: {
 751                bool ok;
 752                int val(v.toInt(&ok));
 753                if (!ok || (val != 0 && val != 1))
 754                    throw Exception("sn_heal must be 0 or 1.");
 755                sn_heal_ = val == 0 ? false : true;
 756                return false;
 757            }
 758        case kwTenTable: {
 759                if (!v.canConvert<DVect2>())
 760                    throw Exception("sn_tentab entry must be a strength and distance.");
 761                DVect2 vl(v.value<DVect2>());
 762                if (vl.x() < 0 || vl.y() < 0)
 763                    throw Exception("The sn_tentab entries must be positive.");
 764                if (vl.y() == 0)
 765                    throw Exception("Use sn_ten to set the tensile strength at 0 elongation.");
 766                tenTable_.push_back(vl);
 767                std::sort(tenTable_.begin(),tenTable_.end(), [](const DVect2& a, const DVect2& b) { return a.y() < b.y(); } );
 768            }
 769            return false;
 770        case kwCohTable: {
 771                if (!v.canConvert<DVect2>())
 772                    throw Exception("sn_cohtab entry must be a strength and distance.");
 773                DVect2 vl(v.value<DVect2>());
 774                if (vl.x() < 0 || vl.y() < 0)
 775                    throw Exception("The sn_cohtab entries must be positive.");
 776                if (vl.y() == 0)
 777                    throw Exception("Use sn_coh to set the cohesive strength at 0 elongation.");
 778                cohTable_.push_back(vl);
 779                std::sort(cohTable_.begin(),cohTable_.end(), [](const DVect2& a, const DVect2& b) { return a.y() < b.y(); } );
 780            }
 781            return false;
 782        case kwTablePos: {
 783                bool ok;
 784                int val(v.toInt(&ok));
 785                if (!ok || val < 1)
 786                    throw Exception("sn_tabpos must be 1 or greater.");
 787                sn_tabPos_ = val - 1;
 788                return false;
 789            }
 790        case kwPorP: {
 791                if (!v.canConvert<double>())
 792                    throw Exception("sn_porp must be a double.");
 793                double val = v.toDouble();
 794                sn_por_ = val;
 795                return true;
 796            }
 797        }
 798        return false;
 799    }
 800
 801    bool ContactModelRBSN::getPropertyReadOnly(uint32 i) const {
 802        // Returns TRUE if a property is read only or FALSE otherwise.
 803        switch (i) {
 804        case kwDpF:
 805        case kwSNS:
 806        case kwSNBS:
 807        case kwSNTS:
 808        case kwEmod:
 809        case kwKRatio:
 810        case kwSNState:
 811        case kwSNRadius:
 812        case kwSNSStr:
 813        case kwSNSig:
 814        case kwSNTau:
 815        case kwPForce:
 816        case kwStressNorm:
 817            return true;
 818        default:
 819            break;
 820        }
 821        return false;
 822    }
 823
 824    bool ContactModelRBSN::supportsInheritance(uint32 i) const {
 825        // Returns TRUE if a property supports inheritance or FALSE otherwise.
 826        switch (i) {
 827        case kwKn:
 828        case kwKs:
 829        case kwFric:
 830            return true;
 831        default:
 832            break;
 833        }
 834        return false;
 835    }
 836
 837    QString  ContactModelRBSN::getMethodArguments(uint32 i) const {
 838        // Return a list of contact model method argument names.
 839        switch (i) {
 840        case kwAssignStiffness:
 841            return "kn,kratio";
 842        case kwStiffness:
 843            return "emod,poisson";
 844        case kwBond:
 845            return "gap";
 846        case kwUnbond:
 847            return "gap";
 848        case kwArea:
 849        case kwResetDisp:
 850            return QString();
 851        }
 852        assert(0);
 853        return QString();
 854    }
 855
 856    bool ContactModelRBSN::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
 857        // Apply the specified method.
 858        IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 859        switch (i) {
 860        case kwAssignStiffness: {
 861                poisson_ = 0.0;
 862                if (vl.at(0).isNull())
 863                    throw Exception("Argument kn must be specified with method assign-stiffness in contact model %1.",getName());
 864                double knpa = vl.at(0).toDouble();
 865                if (knpa<=0.0)
 866                    throw Exception("Negative or zero kn not allowed in contact model %1.",getName());
 867                if (vl.at(1).isNull())
 868                    throw Exception("Argument kratio must be specified with method assign-stiffness in contact model %1.",getName());
 869                kRatio_ = vl.at(1).toDouble();
 870                if (kRatio_<0.0) {
 871                    kRatio_ = 0.0;
 872                    throw Exception("Negative kratio not allowed in contact model %1.",getName());
 873                }
 874                const IContactMechanical *mc = convert_getcast<IContactMechanical>(con);
 875                assert(mc);
 876                std::vector<DVect> pts;
 877                mc->getJointGeometry(&pts);
 878                double area = 0.0;
 879#ifdef THREED
 880                // Step 1: get centroid and area
 881                for (int i=1; i<pts.size()-1; ++i) {
 882                    double a = (pts[0] - pts[i]).mag();
 883                    double b = (pts[0] - pts[i+1]).mag();
 884                    double c = (pts[i] - pts[i+1]).mag();
 885                    double la = 0.0;
 886                    if (b > a)
 887                        std::swap(a,b);
 888                    if (c > a)
 889                        std::swap(a,c);
 890                    if (c > b)
 891                        std::swap(b,c);
 892                    if (c - (a - b) >= 0)
 893                        la = 0.25 * sqrt((a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c)));
 894                    area += la;
 895                }
 896#else
 897                // Assume unit thickness in the out of plane direction
 898                area = (pts[1] - pts[0]).mag();
 899#endif
 900                userArea_ = area;
 901                kTran_ = knpa * area;
 902                E_ = kTran_ / area;
 903                kRot_ = DAVect(0.0);
 904                setInheritance(1,false);
 905                setInheritance(2,false);
 906                sn_mode_ = 1.0;
 907                return true;
 908            }
 909        case kwStiffness: {
 910                FP_S;
 911                poisson_ = 0.0;
 912                if (vl.at(0).isNull())
 913                    throw Exception("Argument emod must be specified with method compute-stiffness in contact model %1.",getName());
 914                E_ = vl.at(0).toDouble();
 915                if (E_<=0.0)
 916                    throw Exception("Negative or zero emod not allowed in contact model %1.",getName());
 917                if (vl.at(1).isNull())
 918                    throw Exception("Argument poisson must be specified with method compute-stiffness in contact model %1.",getName());
 919                poisson_ = vl.at(1).toDouble();
 920                if (poisson_ < 0.0) {
 921                    poisson_ = 0.0;
 922                    throw Exception("Negative poisson not allowed in contact model %1.",getName());
 923                }
 924                const IBody *b1 = con->getEnd1()->getIBody();
 925                const IBody *b2 = con->getEnd2()->getIBody();
 926                //double vol1 = b1->getVolume();
 927                //double vol2 = b2->getVolume();
 928                //if (std::max(vol1,vol2) > 10.*std::min(vol1,vol2))
 929                //    poisson_ = 0.0;
 930                DVect pos1 = toDVect(b1->getIThing()->getLocation());
 931                DVect pos2 = toDVect(b2->getIThing()->getLocation()) + con->getOffSet();
 932                double dist = (pos1-pos2).mag();
 933                if (con->withWall())
 934                    dist = (pos1 - con->getPosition()).mag();
 935                double tol = std::max(pos1.abs().maxComp(),pos2.abs().maxComp())*limits<double>::epsilon()*1000;
 936                if (dist < tol) {
 937                    poisson_ = 0;
 938                    userArea_ = 0;
 939                    kTran_ = 0;
 940                    kRot_.fill(0);
 941                    return true;
 942                }
 943                const IContactMechanical *mc = convert_getcast<IContactMechanical>(con);
 944                assert(mc);
 945                std::vector<DVect> pts;
 946                FP_S;
 947                mc->getJointGeometry(&pts);
 948                FP_S;
 949                double area = 0.0;
 950                DAVect inertia(0.0);
 951#ifdef THREED
 952                // Step 1: get centroid and area
 953                DVect cm(0.0);
 954                for (int i=1; i<pts.size()-1; ++i) {
 955                    DVect lcm = (pts[0] + pts[i] + pts[i+1])/3.0;
 956                    double a = (pts[0] - pts[i]).mag();
 957                    double b = (pts[0] - pts[i+1]).mag();
 958                    double c = (pts[i] - pts[i+1]).mag();
 959                    double la = 0.0;
 960                    if (b > a)
 961                        std::swap(a,b);
 962                    if (c > a)
 963                        std::swap(a,c);
 964                    if (c > b)
 965                        std::swap(b,c);
 966                    if (c - (a - b) >= 0)
 967                        la = 0.25 * sqrt((a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c)));
 968                    cm += lcm * la;
 969                    area += la;
 970                }
 971                FP_S;
 972                if (area == 0.0) {
 973                    poisson_ = 0;
 974                    userArea_ = 0;
 975                    kTran_ = 0;
 976                    kRot_.fill(0);
 977                    return true;
 978                }
 979                cm /= area;
 980                FP_S;
 981                // Step 2 - center it and put in the local system
 982                for (int i=0; i<pts.size(); ++i) {
 983                    pts[i] -= cm;
 984                    pts[i] = con->getLocalSystem().toLocal(pts[i]);
 985                }
 986                // Step 3: compute the polar inertia
 987                for (int i=0; i<pts.size(); ++i) {
 988                    int j = i < pts.size()-1 ? i+1 : 0;
 989                    double xi = pts[i].y();
 990                    double xip1 = pts[j].y();
 991                    double yi = pts[i].z();
 992                    double yip1 = pts[j].z();
 993                    double frnt = (xi*yip1-xip1*yi);
 994                    inertia.ry() += frnt*(xi*xi+xi*xip1+xip1*xip1);
 995                    inertia.rz() += frnt*(yi*yi+yi*yip1+yip1*yip1);
 996                }
 997                inertia.ry() = std::abs(inertia.y() / 12.);
 998                inertia.rz() = std::abs(inertia.z() / 12.);
 999                inertia.rx() = inertia.y() + inertia.z();
1000
1001#else
1002                // Assume unit thickness in the out of plane direction
1003                area = (pts[1] - pts[0]).mag();
1004                inertia.rz() = area*area*area/12.;
1005#endif
1006                userArea_ = area;
1007                kTran_ = E_ * area / dist;
1008                kRot_ = inertia *E_ / dist;
1009                setInheritance(1,false);
1010                setInheritance(2,false);
1011                sn_mode_ = 1.0;
1012                return true;
1013            }
1014        case kwBond: {
1015                if (sn_state_ == 3) return false;
1016                double mingap = -1.0 * limits<double>::max();
1017                double maxgap = 0;
1018                if (vl.at(0).canConvert<double>())
1019                    maxgap = vl.at(0).toDouble();
1020                else if (vl.at(0).canConvert<DVect2>()) {
1021                    DVect2 value = vl.at(0).value<DVect2>();
1022                    mingap = value.minComp();
1023                    maxgap = value.maxComp();
1024                }
1025                else if (!vl.at(0).isNull())
1026                    throw Exception("gap value %1 not recognized in method bond in contact model %2.", vl.at(1), getName());
1027                double gap = c->getGap();
1028                if (gap >= mingap && gap <= maxgap) {
1029                    sn_state_ = 3;
1030                    sn_mode_ = 1;
1031                    return true;
1032                }
1033                return false;
1034            }
1035        case kwUnbond: {
1036                if (sn_state_ == 0) return false;
1037                double mingap = -1.0 * limits<double>::max();
1038                double maxgap = 0;
1039                if (vl.at(0).canConvert<double>())
1040                    maxgap = vl.at(0).toDouble();
1041                else if (vl.at(0).canConvert<DVect2>()) {
1042                    DVect2 value = vl.at(0).value<DVect2>();
1043                    mingap = value.minComp();
1044                    maxgap = value.maxComp();
1045                }
1046                else if (!vl.at(0).isNull())
1047                    throw Exception("gap value %1 not recognized in method unbond in contact model %2.", vl.at(0), getName());
1048                double gap = c->getGap();
1049                if (gap >= mingap && gap <= maxgap) {
1050                    sn_state_ = 0;
1051                    return true;
1052                }
1053                return false;
1054            }
1055        case kwArea: {
1056                if (!userArea_) {
1057                    double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
1058#ifdef THREED
1059                    userArea_ = rsq * rsq * dPi;
1060#else
1061                    userArea_ = rsq * 2.0;
1062#endif
1063                }
1064                return true;
1065            }
1066        case kwResetDisp:
1067            sn_ndisp_ = 0.0;
1068            for (int i=1; i<dim; ++i)
1069                sn_sdisp_.rdof(i) = 0;
1070            break;
1071        }
1072        return false;
1073    }
1074
1075    double ContactModelRBSN::getEnergy(uint32 i) const {
1076        // Return an energy value.
1077        double ret(0.0);
1078        if (!energies_)
1079            return ret;
1080        switch (i) {
1081        case kwEStrain:    return energies_->estrain_;
1082        case kwESlip:      return energies_->eslip_;
1083        case kwEDashpot:   return energies_->edashpot_;
1084        }
1085        assert(0);
1086        return ret;
1087    }
1088
1089    bool ContactModelRBSN::getEnergyAccumulate(uint32 i) const {
1090        // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
1091        switch (i) {
1092        case kwEStrain:   return false;
1093        case kwESlip:     return true;
1094        case kwEDashpot:  return true;
1095        }
1096        assert(0);
1097        return false;
1098    }
1099
1100    void ContactModelRBSN::setEnergy(uint32 i,const double &d) {
1101        // Set an energy value.
1102        if (!energies_) return;
1103        switch (i) {
1104        case kwEStrain:    energies_->estrain_ = d;   return;
1105        case kwESlip:      energies_->eslip_   = d;   return;
1106        case kwEDashpot:   energies_->edashpot_= d;   return;
1107        }
1108        assert(0);
1109        return;
1110    }
1111
1112    bool ContactModelRBSN::validate(ContactModelMechanicalState *state,const double &) {
1113        // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
1114        // (1) the contact is created, (2) a property of the contact that returns a true via
1115        // the setProperty method has been modified and (3) when a set of cycles is executed
1116        // via the {cycle N} command.
1117        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
1118        assert(state);
1119        const IContactMechanical *c = state->getMechanicalContact();
1120        assert(c);
1121
1122        if (state->trackEnergy_)
1123            activateEnergy();
1124
1125        if (inheritanceField_ & KnMask)
1126            updateKn(c);
1127        if (inheritanceField_ & KsMask)
1128            updateKs(c);
1129        if (inheritanceField_ & FricMask)
1130            updateFric(c);
1131
1132        updateStiffness(state);
1133        return checkActivity(state->gap_);
1134    }
1135
1136    static const QString knstr("kn");
1137    bool ContactModelRBSN::updateKn(const IContactMechanical *con) {
1138        assert(con);
1139        QVariant v1 = con->getEnd1()->getProperty(knstr);
1140        QVariant v2 = con->getEnd2()->getProperty(knstr);
1141        if (!v1.isValid() || !v2.isValid())
1142            return false;
1143        double kn1 = v1.toDouble();
1144        double kn2 = v2.toDouble();
1145        double val = kTran_;
1146        if (kn1 && kn2)
1147            kTran_ = kn1*kn2/(kn1+kn2);
1148        else if (kn1)
1149            kTran_ = kn1;
1150        else if (kn2)
1151            kTran_ = kn2;
1152        return ( (kTran_ != val) );
1153    }
1154
1155    static const QString ksstr("ks");
1156    bool ContactModelRBSN::updateKs(const IContactMechanical *con) {
1157        assert(con);
1158        QVariant v1 = con->getEnd1()->getProperty(ksstr);
1159        QVariant v2 = con->getEnd2()->getProperty(ksstr);
1160        if (!v1.isValid() || !v2.isValid())
1161            return false;
1162        double ks1 = v1.toDouble();
1163        double ks2 = v2.toDouble();
1164        double val = kTran_;
1165        if (ks1 && ks2)
1166            kTran_ = ks1*ks2/(ks1+ks2);
1167        else if (ks1)
1168            kTran_ = ks1;
1169        else if (ks2)
1170            kTran_ = ks2;
1171        return ( (kTran_ != val) );
1172    }
1173
1174    static const QString fricstr("fric");
1175    bool ContactModelRBSN::updateFric(const IContactMechanical *con) {
1176        assert(con);
1177        QVariant v1 = con->getEnd1()->getProperty(fricstr);
1178        QVariant v2 = con->getEnd2()->getProperty(fricstr);
1179        if (!v1.isValid() || !v2.isValid())
1180            return false;
1181        double fric1 = std::max(0.0,v1.toDouble());
1182        double fric2 = std::max(0.0,v2.toDouble());
1183        double val = fric_;
1184        fric_ = std::min(fric1,fric2);
1185        return ( (fric_ != val) );
1186    }
1187
1188    bool ContactModelRBSN::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
1189        // The endPropertyUpdated method is called whenever a surface property (with a name
1190        // that matches an inheritable contact model property name) of one of the contacting
1191        // pieces is modified. This allows the contact model to update its associated
1192        // properties. The return value denotes whether or not the update has affected
1193        // the time step computation (by having modified the translational or rotational
1194        // tangent stiffnesses). If true is returned, then the time step will be recomputed.
1195        assert(c);
1196        QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
1197        QRegularExpression rx(name, QRegularExpression::CaseInsensitiveOption);
1198        int idx = availableProperties.indexOf(rx)+1;
1199        bool ret=false;
1200
1201        if (idx<=0)
1202            return ret;
1203
1204        switch(idx) {
1205        case kwKn:  { //kn
1206                if (inheritanceField_ & KnMask)
1207                    ret = updateKn(c);
1208                break;
1209            }
1210        case kwKs:  { //ks
1211                if (inheritanceField_ & KsMask)
1212                    ret =updateKs(c);
1213                break;
1214            }
1215        case kwFric:  { //fric
1216                if (inheritanceField_ & FricMask)
1217                    updateFric(c);
1218                break;
1219            }
1220        }
1221        return ret;
1222    }
1223
1224    ContactModelRBSN::StiffData ContactModelRBSN::computeStiffData(ContactModelMechanicalState *state) const {
1225        // Update contact data
1226        //double Cmin1 = state->end1Curvature_.x();
1227        double Cmax1 = state->end1Curvature_.y();
1228        double Cmax2 = state->end2Curvature_.y();
1229        double br = sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
1230        if (userArea_)
1231#ifdef THREED
1232            br = std::sqrt(userArea_ / dPi);
1233#else
1234            br = userArea_ / 2.0;
1235#endif
1236        StiffData ret;
1237        ret.reff_ = br;
1238        ret.trans_ = DVect2(kTran_,kTran_/kRatio_);
1239        ret.ang_ = kRot_;
1240        return ret;
1241    }
1242
1243    void ContactModelRBSN::updateStiffness(ContactModelMechanicalState *state) {
1244        // first compute stiffness data
1245        StiffData stiff = computeStiffData(state);
1246        // Now calculate effective stiffness
1247        DVect2 retT = stiff.trans_;
1248        // correction if viscous damping active
1249        if (dpProps_) {
1250            DVect2 correct(1.0);
1251            if (dpProps_->dp_nratio_)
1252                correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
1253            if (dpProps_->dp_sratio_)
1254                correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
1255            retT /= (correct*correct);
1256        }
1257        effectiveTranslationalStiffness_ = retT;
1258        // Effective rotational stiffness (bending and twisting)
1259        effectiveRotationalStiffness_ = stiff.ang_;
1260    }
1261
1262    bool ContactModelRBSN::forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) {
1263        assert(state);
1264
1265        if (state->activated()) {
1266            // The contact was just activated from an inactive state
1267            // Trigger the FISH callback if one is hooked up to the
1268            // contact_activated event.
1269            if (cmEvents_[fActivated] >= 0) {
1270                auto c = state->getContact();
1271                std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
1272                IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1273                fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
1274            }
1275        }
1276        updateStiffness(state);
1277        // accumulate shear displacement for dilation
1278        sn_ndisp_ += state->relativeTranslationalIncrement_.x();
1279        DVect shearInc =  state->relativeTranslationalIncrement_;
1280        shearInc.rx() = 0;
1281        sn_sdisp_.ry() += shearInc.mag();
1282
1283        if (isBonded()) return FDLawBonded(state, timestep);
1284        else return FDLawUnBonded(state, timestep);
1285
1286    }
1287    
1288    bool ContactModelRBSN::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
1289        // Account for thermal expansion in incremental mode
1290        if (sn_mode_ == 0 || ts->gapInc_ == 0.0) return false;
1291        DVect finc(0.0);
1292        finc.rx() = kTran_ * ts->gapInc_;
1293        sn_F_ -= finc;
1294        return true;
1295    }
1296
1297    bool ContactModelRBSN::FDLawBonded(ContactModelMechanicalState *state, const double &timestep) {
1298        // initialize ... get the geometry information
1299        DVect3 geom = computeGeomData(state->getMechanicalContact());
1300        double area = geom.x();
1301        double bi = geom.y();
1302        double br = geom.z();
1303        double kn = kTran_;
1304        double ks = kn / kRatio_;
1305        double kb = kRot_.z();
1306#if DIM==3
1307        kb = std::sqrt(kb*kb + kRot_.y()*kRot_.y());
1308        double kt = kRot_.x();
1309#else
1310        double kt = 0.0;
1311#endif
1312
1313        DVect totalForce = sn_F_ + fictForce_;
1314
1315        //double nsmax0 = -(totalForce.x() / area) + sn_mcf_* sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z()) * br / bi;
1316
1317        // Relative translational/rotational displacement increments
1318        DVect  trans = state->relativeTranslationalIncrement_;
1319        DAVect ang = state->relativeAngularIncrement_;
1320
1321        // Store previous force and moment
1322        DVect  sn_F_old = totalForce;
1323        sn_F_old.rx() -= sn_por_ * area;
1324        DAVect sn_M_old = sn_M_;
1325
1326        DVect theStiff(ks);
1327        theStiff.rx() = kn;
1328        sn_F_ -= trans * theStiff;
1329        sn_M_ -= ang * kRot_;
1330        if (poisson_ != 0.0) {
1331            const IContact *con = state->getContact();
1332            const IBody *b1 = con->getEnd1()->getIBody();
1333            const IBody *b2 = con->getEnd2()->getIBody();
1334#ifdef THREED
1335            std::array<double,6> stress11 = {0,0,0,0,0,0};
1336            std::array<double,6> stress22 = {0,0,0,0,0,0};
1337#else
1338            std::array<double,3> stress11 = {0,0,0};
1339            std::array<double,3> stress22 = {0,0,0};
1340#endif
1341
1342            double vol1 = b1->getVolume();
1343            double vol2 = b2->getVolume();
1344            b1->getOldStress(stress11);
1345            b2->getOldStress(stress22);
1346            double ms = 0.0;
1347            for (int i=0; i<stress11.size(); ++i)  {
1348                stress11[i] = (stress11[i]*vol1 + stress22[i]*vol2)/(vol1 + vol2);
1349                ms = std::max(ms,abs(stress11[i]));
1350            }
1351            DMatrix<dim,dim> stresst(0.0);
1352#ifdef THREED
1353            stresst.get(0,0) = -poisson_*stress11[1] - poisson_*stress11[2];
1354            stresst.get(1,1) = -poisson_*stress11[0] - poisson_*stress11[2];
1355#else
1356            stresst.get(0,0) = -poisson_*stress11[1];
1357            stresst.get(1,1) = -poisson_*stress11[0];
1358#endif
1359#ifdef THREED
1360            stresst.get(2,2) = -poisson_*stress11[0] - poisson_*stress11[1];
1361#endif
1362            if (poisOffDiag_) {
1363#ifdef THREED
1364                double sxy = stress11[3];
1365                double szx = stress11[4];
1366                double syz = stress11[5];
1367                stresst.get(0,1) = poisson_ * sxy;
1368                stresst.get(1,0) = stresst.get(0,1);
1369                stresst.get(0,2) = poisson_ * szx;
1370                stresst.get(2,0) = stresst.get(0,2);
1371                stresst.get(1,2) = poisson_ * syz;
1372                stresst.get(2,1) = stresst.get(1,2);
1373#else
1374                double sxy = stress11[2];
1375                stresst.get(0,1) = poisson_ * sxy;
1376                stresst.get(1,0) = stresst.get(0,1);
1377#endif
1378            }
1379
1380
1381            DVect norm = toVect(con->getNormal());
1382            DVect traction = stresst * norm * userArea_;
1383            DVect shear(0.0);
1384            shear.ry() = 1.0;
1385            shear = con->getLocalSystem().toGlobal(shear);
1386#ifdef THREED
1387            DVect ns = con->getLocalSystem().toGlobal(DVect(0.,0.,1.));
1388            fictForce_ = DVect(norm|traction,shear|traction,ns|traction);
1389#else
1390            fictForce_ = DVect(norm|traction,shear|traction);
1391#endif
1392            if (forceSet_ && ms) {
1393                forceSet_ = false;
1394                sn_F_ -= fictForce_;
1395            }
1396        }
1397        FP_S;
1398        double porForce = sn_por_ * area;
1399        sn_F_.rx() -= porForce;
1400        totalForce = sn_F_ + fictForce_;
1401
1402        if (state->canFail_) {
1403            double dbend = sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z());
1404            double nsmax = -(totalForce.x() / area) + sn_mcf_*dbend * br / bi;
1405            bool failed = false;
1406            if (sn_state_ == 3 || sn_state_ == 5) {
1407                double compVal = sn_state_ == 3 ? tenTable_[0].x() : transTen_;
1408                if (nsmax >= compVal ) {
1409                    if (tenTable_.back().y() < limits<double>::epsilon()*100)
1410                        failed = true;
1411                    else {
1412                        if (sn_state_ == 3)
1413                            sn_elong_ = 0;
1414                        transTen_ = compVal;
1415                        sn_state_ = 4;
1416                    }
1417                }
1418            }
1419            FP_S;
1420
1421            if (sn_state_ == 4) {
1422                sn_elong_ += trans.x();
1423                sn_elong_ = std::max(0.0,sn_elong_);
1424                int ww = -1;
1425                if (sn_elong_ <= tenTable_.back().y()) {
1426                    for (int i=0; i<tenTable_.size(); ++i) {
1427                        if (tenTable_[i].y() >= sn_elong_) {
1428                            ww = i;
1429                            break;
1430                        }
1431                    }
1432                } else
1433                    ww = tenTable_.size() - 1;
1434                if (ww > 0) {
1435                    //double factor = std::min(1.0, sn_elong_ / tenTable_[ww].y());
1436                    double prevVal = ww == 1 ? 1 : tenTable_[ww-1].x();
1437                    double curVal =  tenTable_[ww].x();
1438                    double prevElong = tenTable_[ww-1].y();
1439                    double curElong = tenTable_[ww].y();
1440                    double slope = (curVal - prevVal)/(curElong - prevElong);
1441                    FP_S;
1442                    //y-y0 = m(x-x0)
1443                    double nstren = slope * (sn_elong_ - prevElong) + prevVal;
1444                    if (nstren <= 0)
1445                        failed = true;
1446                    else {
1447                        nstren *= tenTable_[0].x();
1448                        if (nsmax >= nstren || slope > 0) {
1449                            double fac = nstren / nsmax;
1450                            sn_F_.rx() *= fac;
1451#if DIM==3
1452                            sn_M_.ry() *= fac;
1453#endif
1454                            sn_M_.rz() *= fac;
1455                            fictForce_.rx() *= fac;
1456                        } else {
1457                            sn_state_ = 5;
1458                            transTen_ = -(sn_F_old.x() / area) + sn_mcf_* sqrt(sn_M_old.y()*sn_M_old.y() + sn_M_old.z()*sn_M_old.z()) * br / bi;
1459                        }
1460                    }
1461                }
1462            }
1463            if (sn_state_ == 6 && nsmax >= 0)
1464                failed = true;
1465            FP_S;
1466            if (failed) {
1467                // Failed in tension
1468                double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1469                sn_state_ = 1;
1470                sn_F_.fill(0.0);
1471                sn_M_.fill(0.0);
1472                failed = true;
1473                fictForce_ = DVect(0.0);
1474                //sn_F_.rx() = -sn_tenres_ * area;
1475                if (cmEvents_[fBondBreak] >= 0) {
1476                    auto c = state->getContact();
1477                    std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1478                                                         fish::Parameter((qint64)sn_state_),
1479                                                         fish::Parameter(nsmax),
1480                                                         fish::Parameter(se)
1481                                                       };
1482                    IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1483                    fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1484                }
1485            }
1486            FP_S;
1487
1488            if (!failed) {
1489                /* check for shear failure */
1490                double dtwist = sn_M_.x();
1491                DVect bfs(totalForce);
1492                bfs.rx() = 0.0;
1493                double dbfs = bfs.mag();
1494                double ssmax = dbfs / area + sn_mcf_*std::abs(dtwist) * 0.5* br / bi;
1495                double ss = shearStrength(area);
1496                FP_S;
1497                if (ss < 0)
1498                    ss = 0;
1499                if (ss <= ssmax) {
1500                    // strength when it breaks for
1501                    // Failed in shear
1502                    double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1503                    sn_state_ = 2;
1504                    fictForce_ = DVect(0.0);
1505                    FP_S;
1506                    sn_F_ = totalForce;
1507                    if (cmEvents_[fBondBreak] >= 0) {
1508                        auto c = state->getContact();
1509                        std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1510                                                             fish::Parameter((qint64)sn_state_),
1511                                                             fish::Parameter(ss),
1512                                                             fish::Parameter(se)
1513                                                           };
1514                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1515                        fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1516                    }
1517                    double mm = 0.0;
1518                    for (int i=1; i<dim; ++i)
1519                        mm += trans[i]*trans[i];
1520                    sn_sdisp_.rx() = std::sqrt(mm);
1521                    double crit = sn_F_.x() * fric_ + sn_cohres_ * area;
1522                    if (sn_cohdc() > std::numeric_limits<double>::epsilon()) {
1523                        int ww = -1;
1524                        if (sn_sdisp_.x() <= sn_cohdc()) {
1525                            for (int i=0; i<cohTable_.size(); ++i) {
1526                                if (cohTable_[i].y() >= sn_sdisp_.x()) {
1527                                    ww = i;
1528                                    break;
1529                                }
1530                            }
1531                        } else
1532                            ww = cohTable_.size() - 1;
1533                        if (ww > 0) {
1534                            double prevVal = ww == 1 ? 1: cohTable_[ww-1].x() ;//critBreak_ : cohTable_[ww-1].x() + sn_F_.x() * fric_;
1535                            double curVal =  cohTable_[ww].x();//cohTable_[ww].x() + sn_F_.x() * fric_;
1536                            double prevShear = cohTable_[ww-1].y();
1537                            double curShear = cohTable_[ww].y();
1538                            double slope = (curVal - prevVal)/(curShear - prevShear);
1539                            // should go from 0 to 1
1540                            double mval = slope * (sn_sdisp_.x() - prevShear) + prevVal;
1541                            double fact = std::max(0.,std::min(1.0,1.0 - mval));
1542                            assert(fact >= 0 && fact <= 1.0);
1543                            double theFric = fric_;
1544                            if (sn_fa_)
1545                                theFric = sn_fa_ + (fric_ - sn_fa_)*fact;
1546                            double theCoh = sn_Coh();
1547                            if (theCoh)
1548                                theCoh = sn_Coh() + (sn_cohres_ - sn_Coh())*fact;
1549                            crit = sn_F_.x() * theFric + theCoh * geom.x();
1550                        }
1551                    }
1552
1553                    // Resolve sliding.
1554                    FP_S;
1555                    if (crit < 0)
1556                        crit = 0;
1557                    DVect sforce = sn_F_; sforce.rx() = 0.0;
1558                    // The is the magnitude of the shear force.
1559                    double sfmag = sforce.mag();
1560                    // Sliding occurs when the magnitude of the shear force is greater than the
1561                    // critical value.
1562                    if (sfmag > crit) {
1563                        // Lower the shear force to the critical value for sliding.
1564                        double rat = crit / sfmag;
1565                        sforce *= rat;
1566                        sforce.rx() = sn_F_.x();
1567                        sn_F_ = sforce;
1568                        sn_S_ = true;
1569                    }
1570                    if (sn_S_) {
1571                        if (sn_dil_ > 0) {
1572                            double zdd = sn_dilzero_ != 0 ? sn_dilzero_ : limits<double>::max();
1573                            double usm = sn_sdisp_.y();
1574                            if (usm < zdd) {
1575                                double sInc = 0.0;
1576                                for (int i=1; i<dim; ++i)
1577                                    sInc += trans.dof(i)*trans.dof(i);
1578                                sInc = std::sqrt(sInc);
1579                                sn_F_.rx() += kTran_ * sn_dil_ * sInc;
1580                            }
1581                        }
1582                    }
1583
1584                    // Resolve bending
1585                    crit = sn_bmul_*2.1*0.25*br*std::abs(sn_F_.x()); // Jiang 2015
1586                    FP_S;
1587                    DAVect test = sn_M_;
1588#if DIM==3
1589                    test.rx() = 0.0;
1590#endif
1591                    double tmag = test.mag();
1592                    if (tmag > crit) {
1593                        // Lower the bending moment to the critical value for sliding.
1594                        double rat = crit / tmag;
1595                        test *= rat;
1596                        sn_BS_ = true;
1597                    }
1598                    sn_M_.rz() = test.z();
1599#if DIM==3
1600                    sn_M_.ry() = test.y();
1601                    // Resolve twisting
1602                    crit = sn_tmul_ * 0.65*fric_* br*std::abs(sn_F_.x()) ; // Jiang 2015
1603                    tmag = std::abs(sn_M_.x());
1604                    if (tmag > crit) {
1605                        // Lower the shear force to the critical value for sliding.
1606                        double rat = crit / tmag;
1607                        tmag = sn_M_.x() * rat;
1608                        sn_TS_ = true;
1609                    }
1610                    sn_M_.rx() = tmag;
1611                    FP_S;
1612#endif
1613                }
1614            }
1615        }
1616        sn_F_old.rx() += porForce;
1617        sn_F_.rx() += porForce;
1618        totalForce = sn_F_ + fictForce_;
1619        FP_S;
1620
1621        // Account for dashpot forces if the dashpot structure has been defined.
1622        if (dpProps_) {
1623            dpProps_->dp_F_.fill(0.0);
1624            double vcn(0.0), vcs(0.0);
1625            // Calculate the damping coefficients.
1626            vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1627            vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1628            // First damp the shear components
1629            for (int i = 1; i<dim; ++i)
1630                dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1631            // Damp the normal component
1632            dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1633            // Need to change behavior based on the dp_mode.
1634            if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1635                // Limit in tension if not bonded.
1636                if (sn_state_ < 3 && (dpProps_->dp_F_.x() + totalForce.x() < 0))
1637                    dpProps_->dp_F_.rx() = -totalForce.rx();
1638            }
1639            if (sn_S_ && dpProps_->dp_mode_ > 1) {
1640                // Limit in shear if sliding.
1641                double dfn = dpProps_->dp_F_.rx();
1642                dpProps_->dp_F_.fill(0.0);
1643                dpProps_->dp_F_.rx() = dfn;
1644            }
1645        }
1646        FP_S;
1647
1648        //Compute energies if energy tracking has been enabled.
1649        if (state->trackEnergy_) {
1650            assert(energies_);
1651            energies_->estrain_ = 0.0;
1652            if (kn)
1653                // Calculate the strain energy.
1654                energies_->estrain_ = 0.5*totalForce.x()*totalForce.x() / kn;
1655            if (ks) {
1656                DVect s = totalForce;
1657                s.rx() = 0.0;
1658                double smag2 = s.mag2();
1659                // Add the shear component of the strain energy.
1660                energies_->estrain_ += 0.5*smag2 / ks;
1661
1662                if (sn_S_) {
1663                    // If sliding calculate the slip energy and accumulate it.
1664                    sn_F_old.rx() = 0.0;
1665                    DVect avg_F_s = (s + sn_F_old)*0.5;
1666                    DVect u_s_el = (s - sn_F_old) / ks;
1667                    DVect u_s(0.0);
1668                    for (int i = 1; i < dim; ++i)
1669                        u_s.rdof(i) = trans.dof(i);
1670                    energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1671                }
1672            }
1673            // Add the bending/twisting resistance energy contributions.
1674            if (kb) {
1675                DAVect tmp = sn_M_;
1676#ifdef THREED
1677                tmp.rx() = 0.0;
1678#endif
1679                energies_->estrain_ += 0.5*tmp.mag2() / kb;
1680                if (sn_BS_) {
1681                    //  accumulate bending slip energy.
1682                    DAVect tmp_old = sn_M_old;
1683#ifdef THREED
1684                    tmp_old.rx() = 0.0;
1685#endif
1686                    DAVect avg_M = (tmp + tmp_old)*0.5;
1687                    DAVect t_s_el = (tmp - tmp_old) / kb;
1688                    energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1689                }
1690            }
1691#ifdef THREED
1692            if (kt) {
1693                double mt = std::abs(sn_M_.x());
1694                energies_->estrain_ += 0.5*mt*mt / kt;
1695                if (sn_TS_) {
1696                    //  accumulate twisting slip energy.
1697                    DAVect tmp(0.0);
1698                    DAVect tmp_old(0.0);
1699                    tmp.rx() = sn_M_.x();
1700                    tmp_old.rx() = sn_M_old.x();
1701                    DAVect avg_M = (tmp + tmp_old)*0.5;
1702                    DAVect t_s_el = (tmp - tmp_old) / kt;
1703                    energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1704                }
1705            }
1706#endif
1707
1708            if (dpProps_) {
1709                // Calculate damping energy (accumulated) if the dashpots are active.
1710                energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1711            }
1712        }
1713
1714        // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1715        assert(sn_F_ == sn_F_);
1716        assert(sn_M_ == sn_M_);
1717        assert(fictForce_ == fictForce_);
1718        FP_S;
1719        return true;
1720    }
1721
1722    bool ContactModelRBSN::FDLawUnBonded(ContactModelMechanicalState *state, const double &timestep) {
1723        DVect3 geom = computeGeomData(state->getMechanicalContact());
1724        double br = geom.z();
1725        // Relative translational/rotational displacement increments
1726        DVect  trans = state->relativeTranslationalIncrement_;
1727        DAVect ang   = state->relativeAngularIncrement_;
1728        double overlap = rgap_ - state->gap_;
1729        double correction = 1.0;
1730        if (state->activated() && sn_mode_ == 0 && trans.x()) {
1731                correction = -1.0*overlap / trans.x();
1732                if (correction < 0)
1733                    correction = 1.0;
1734        }
1735
1736        // Store previous force and moment
1737        DVect  sn_F_old = sn_F_;
1738        double porForce = sn_por_ * geom.x();
1739        sn_F_old.rx() -= porForce;
1740        DAVect sn_M_old = sn_M_;
1741
1742        double kb = kRot_.z();
1743#if DIM==3
1744        double kt = kRot_.x();
1745        //kb = std::sqrt(kb * kb + kRot_.y() * kRot_.y());
1746#endif
1747        // absolute/incremental normal force calculation
1748        DVect theStiff(kTran_/kRatio_);
1749        theStiff.rx() = kTran_;
1750
1751        if (sn_mode_==0)
1752            sn_F_.rx() = overlap * theStiff.x();
1753        else
1754            sn_F_.rx() -= trans.x() * theStiff.x();
1755        // Normal force can only be positive if unbonded
1756        sn_F_.rx() = std::max(0.0, sn_F_.x()) - porForce;
1757
1758        // Calculate the trial shear force.
1759        DVect sforce(0.0);
1760        // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1761        // Loop over the shear components (note: the 0 component is the normal component)
1762        // and calculate the shear force.
1763        for (int i = 1; i<dim; ++i)
1764            sforce.rdof(i) = sn_F_.dof(i) - trans.dof(i) * theStiff.dof(i);
1765
1766        // Calculate the trial moment.
1767        DAVect mom = sn_M_ - ang*kRot_;
1768
1769        // If the SOLVE ELASTIC command is given then the
1770        // canFail state is set to FALSE. Otherwise it is always TRUE.
1771        if (state->canFail_) {
1772            bool changed = false;
1773            // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
1774            bool slip_changed = false;
1775
1776            double crit = sn_F_.x() * fric_ + sn_cohres_ * geom.x();
1777            if (sn_state_ != 0) {
1778                bool doComp = true; 
1779                if (!sn_S_) {
1780                    if (sn_heal_) {
1781                        sn_sdisp_.rx() = 0;
1782                        crit = sn_F_.x() * sn_fa_ + cohTable_[0].x() * geom.x();
1783                        doComp = false;
1784                    }
1785                }
1786                if (doComp) {
1787                    double mm = 0.0;
1788                    for (int i=1; i<dim; ++i)
1789                        mm += trans[i]*trans[i];
1790                    sn_sdisp_.rx() += std::sqrt(mm);
1791                    if (sn_cohdc() > std::numeric_limits<double>::epsilon()) {
1792                        int ww = -1;
1793                        if (sn_sdisp_.x() <= sn_cohdc()) {
1794                            for (int i=0; i<cohTable_.size(); ++i) {
1795                                if (cohTable_[i].y() >= sn_sdisp_.x()) {
1796                                    ww = i;
1797                                    break;
1798                                }
1799                            }
1800                        } else
1801                            ww = cohTable_.size() - 1;
1802                        if (ww > 0) {
1803                            double prevVal = ww == 1 ? 1: cohTable_[ww-1].x() ;//critBreak_ : cohTable_[ww-1].x() + sn_F_.x() * fric_;
1804                            double curVal =  cohTable_[ww].x();//cohTable_[ww].x() + sn_F_.x() * fric_;
1805                            double prevShear = cohTable_[ww-1].y();
1806                            double curShear = cohTable_[ww].y();
1807                            double slope = (curVal - prevVal)/(curShear - prevShear);
1808                            // should go from 0 to 1
1809                            double mval = slope * (sn_sdisp_.x() - prevShear) + prevVal;
1810                            double fact = std::max(0.,std::min(1.0,1.0 - mval));
1811                            assert(fact >= 0 && fact <= 1.0);
1812                            double theFric = fric_;
1813                            if (sn_fa_)
1814                                theFric = sn_fa_ + (fric_ - sn_fa_)*fact;
1815                            double theCoh = sn_Coh();
1816                            if (theCoh)
1817                                theCoh = sn_Coh() + (sn_cohres_ - sn_Coh())*fact;
1818                            crit = sn_F_.x() * theFric + theCoh * geom.x();
1819                        }
1820                    }
1821                }
1822            }
1823
1824            // Resolve sliding.
1825            if (crit < 0)
1826                crit = 0.0;
1827            // The is the magnitude of the shear force.
1828            double sfmag = sforce.mag();
1829            if (sfmag > crit) {
1830                // Lower the shear force to the critical value for sliding.
1831                double rat = crit / sfmag;
1832                sforce *= rat;
1833                if (!sn_S_) {
1834                    slip_changed = true;
1835                    changed = true;
1836                }
1837                sn_S_ = true;
1838            } else {
1839                if (sn_S_) {
1840                    slip_changed = true;
1841                    changed = true;
1842                }
1843                sn_S_ = false;
1844            }
1845            if (sn_S_) {
1846                if (sn_dil_ > 0) {
1847                    double zdd = sn_dilzero_ != 0 ? sn_dilzero_ : limits<double>::max();
1848                    double usm = sn_sdisp_.y();
1849                    if (usm < zdd) {
1850                        double sInc = 0.0;
1851                        for (int i=1; i<dim; ++i)
1852                            sInc += trans.dof(i)*trans.dof(i);
1853                        sInc = std::sqrt(sInc);
1854                        sn_F_.rx() += kTran_ * sn_dil_ * sInc;
1855                    }
1856                }
1857            } else {
1858                if (sn_heal_) {
1859                    if (shearStrength(geom.x()))
1860                        sn_state_ = 6;
1861                }
1862            }
1863            // Resolve bending
1864            bool bslip_changed = false;
1865            crit = sn_bmul_ * 2.1*0.25*sn_F_.x() * br; // Jiang 2015
1866            DAVect test = mom;
1867#if DIM==3
1868            test.rx() = 0.0;
1869#endif
1870            double tmag = test.mag();
1871            if (tmag > crit) {
1872                // Lower the bending moment to the critical value for sliding.
1873                double rat = crit / tmag;
1874                test *= rat;
1875                if (!sn_BS_) {
1876                    bslip_changed = true;
1877                    changed = true;
1878                }
1879                sn_BS_ = true;
1880            }
1881            else {
1882                if (sn_BS_) {
1883                    bslip_changed = true;
1884                    changed = true;
1885                }
1886                sn_BS_ = false;
1887            }
1888            mom.rz() = test.z();
1889#if DIM==3
1890            mom.ry() = test.y();
1891            // Resolve twisting
1892            bool tslip_changed = false;
1893            crit = sn_tmul_ * 0.65*fric_*sn_F_.x() * br; // Jiang 2015
1894            tmag = std::abs(mom.x());
1895            if (tmag > crit) {
1896                // Lower the twisting moment to the critical value for sliding.
1897                double rat = crit / tmag;
1898                mom.rx() *= rat;
1899                if (!sn_TS_) {
1900                    tslip_changed = true;
1901                    changed = true;
1902                }
1903                sn_TS_ = true;
1904            }
1905            else {
1906                if (sn_TS_) {
1907                    tslip_changed = true;
1908                    changed = true;
1909                }
1910                sn_TS_ = false;
1911            }
1912#endif
1913            if (changed && cmEvents_[fSlipChange] >= 0) {
1914                qint64 code = 0;
1915                if (slip_changed) {
1916                    code = 1;
1917                    if (bslip_changed) {
1918                        code = 4;
1919#if DIM==3
1920                        if (tslip_changed)
1921                            code = 7;
1922#endif
1923                    }
1924                }
1925                else if (bslip_changed) {
1926                    code = 2;
1927#if DIM==3
1928                    if (tslip_changed)
1929                        code = 6;
1930#endif
1931                }
1932#if DIM==3
1933                else if (tslip_changed) {
1934                    code = 3;
1935                    if (slip_changed)
1936                        code = 5;
1937                }
1938#endif
1939                auto c = state->getContact();
1940                std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1941                                                     fish::Parameter(code),
1942                                                     fish::Parameter(sn_S_),
1943                                                     fish::Parameter(sn_BS_)
1944#ifdef THREED
1945                                                     ,fish::Parameter(sn_TS_)
1946#endif
1947                                                   };
1948                IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1949                fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1950            }
1951        }
1952
1953        sn_F_.rx() += porForce;
1954        sn_F_old.rx() += porForce;
1955
1956        // Set the shear components of the total force.
1957        for (int i = 1; i<dim; ++i)
1958            sn_F_.rdof(i) = sforce.dof(i);
1959
1960        // Set the moment.
1961        sn_M_ = mom;
1962
1963        // Account for dashpot forces if the dashpot structure has been defined.
1964        if (dpProps_) {
1965            dpProps_->dp_F_.fill(0.0);
1966            double vcn(0.0), vcs(0.0);
1967            // Calculate the damping coefficients.
1968            vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kTran_)));
1969            vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(kTran_/kRatio_)));
1970            // First damp the shear components
1971            for (int i = 1; i<dim; ++i)
1972                dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1973            // Damp the normal component
1974            dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1975            // Need to change behavior based on the dp_mode.
1976            if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1977                // Limit in tension if not bonded.
1978                if (dpProps_->dp_F_.x() + sn_F_.x() < 0)
1979                    dpProps_->dp_F_.rx() = -sn_F_.rx();
1980            }
1981            if (sn_S_ && dpProps_->dp_mode_ > 1) {
1982                // Limit in shear if not sliding.
1983                double dfn = dpProps_->dp_F_.rx();
1984                dpProps_->dp_F_.fill(0.0);
1985                dpProps_->dp_F_.rx() = dfn;
1986            }
1987        }
1988
1989        //Compute energies if energy tracking has been enabled.
1990        if (state->trackEnergy_) {
1991            assert(energies_);
1992            energies_->estrain_ = 0.0;
1993            if (kTran_)
1994                // Calculate the strain energy.
1995                energies_->estrain_ = 0.5*sn_F_.x()*sn_F_.x() / kTran_;
1996            if (kTran_) {
1997                DVect s = sn_F_;
1998                s.rx() = 0.0;
1999                double smag2 = s.mag2();
2000                // Add the shear component of the strain energy.
2001                energies_->estrain_ += 0.5*smag2 / (kTran_/kRatio_);
2002
2003                if (sn_S_) {
2004                    // If sliding calculate the slip energy and accumulate it.
2005                    sn_F_old.rx() = 0.0;
2006                    DVect avg_F_s = (s + sn_F_old)*0.5;
2007                    DVect u_s_el = (s - sn_F_old) / (kTran_/kRatio_);
2008                    DVect u_s(0.0);
2009                    for (int i = 1; i < dim; ++i)
2010                        u_s.rdof(i) = trans.dof(i);
2011                    energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
2012                }
2013            }
2014            // Add the bending/twisting resistance energy contributions.
2015            if (kb) {
2016                DAVect tmp = sn_M_;
2017#ifdef THREED
2018                tmp.rx() = 0.0;
2019#endif
2020                energies_->estrain_ += 0.5*tmp.mag2() / kb;
2021                if (sn_BS_) {
2022                    //  accumulate bending slip energy.
2023                    DAVect tmp_old = sn_M_old;
2024#ifdef THREED
2025                    tmp_old.rx() = 0.0;
2026#endif
2027                    DAVect avg_M = (tmp + tmp_old)*0.5;
2028                    DAVect t_s_el = (tmp - tmp_old) / kb;
2029                    energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
2030                }
2031            }
2032#ifdef THREED
2033            if (kt) {
2034                double mt = std::abs(sn_M_.x());
2035                energies_->estrain_ += 0.5*mt*mt / kt;
2036                if (sn_TS_) {
2037                    //  accumulate twisting slip energy.
2038                    DAVect tmp(0.0);
2039                    DAVect tmp_old(0.0);
2040                    tmp.rx() = sn_M_.x();
2041                    tmp_old.rx() = sn_M_old.x();
2042                    DAVect avg_M = (tmp + tmp_old)*0.5;
2043                    DAVect t_s_el = (tmp - tmp_old) / kt;
2044                    energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
2045                }
2046            }
2047#endif
2048
2049            if (dpProps_) {
2050                // Calculate damping energy (accumulated) if the dashpots are active.
2051                energies_->edashpot_ -= dpProps_->dp_F_ | trans;
2052            }
2053        }
2054
2055        // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
2056        assert(sn_F_ == sn_F_);
2057        assert(fictForce_ == fictForce_);
2058        assert(sn_M_ == sn_M_);
2059        return true;
2060    }
2061
2062    void ContactModelRBSN::setForce(const DVect &v,IContact *c) {
2063        sn_F_ = v;
2064        fictForce_ = DVect(0.0);
2065        forceSet_ = true;
2066        if (v.x() > 0) 
2067            rgap_ = c->getGap() + v.x() / kTran_;
2068    }
2069
2070    void ContactModelRBSN::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
2071        // Only called for contacts with wall facets when the wall resolution scheme
2072        // is set to full!
2073        // Only do something if the contact model is of the same type
2074        if (old->getContactModel()->getName().compare("springnetwork",Qt::CaseInsensitive) == 0 && !isBonded()) {
2075            ContactModelRBSN *oldCm = (ContactModelRBSN *)old;
2076#ifdef THREED
2077            // Need to rotate just the shear component from oldSystem to newSystem
2078
2079            // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
2080            DVect axis = oldSystem.e1() & newSystem.e1();
2081            double c, ang, s;
2082            DVect re2;
2083            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
2084                axis = axis.unit();
2085                c = oldSystem.e1()|newSystem.e1();
2086                if (c > 0)
2087                    c = std::min(c,1.0);
2088                else
2089                    c = std::max(c,-1.0);
2090                ang = acos(c);
2091                s = sin(ang);
2092                double t = 1. - c;
2093                DMatrix<3,3> rm;
2094                rm.get(0,0) = t*axis.x()*axis.x() + c;
2095                rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
2096                rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
2097                rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
2098                rm.get(1,1) = t*axis.y()*axis.y() + c;
2099                rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
2100                rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
2101                rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
2102                rm.get(2,2) = t*axis.z()*axis.z() + c;
2103                re2 = rm*oldSystem.e2();
2104            }
2105            else
2106                re2 = oldSystem.e2();
2107            // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
2108            axis = re2 & newSystem.e2();
2109            DVect2 tpf;
2110            DVect2 tpm;
2111            DMatrix<2,2> m;
2112            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
2113                axis = axis.unit();
2114                c = re2|newSystem.e2();
2115                if (c > 0)
2116                    c = std::min(c,1.0);
2117                else
2118                    c = std::max(c,-1.0);
2119                ang = acos(c);
2120                if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
2121                    ang *= -1;
2122                s = sin(ang);
2123                m.get(0,0) = c;
2124                m.get(1,0) = s;
2125                m.get(0,1) = -m.get(1,0);
2126                m.get(1,1) = m.get(0,0);
2127                tpf = m*DVect2(oldCm->sn_F_.y(),oldCm->sn_F_.z());
2128                tpm = m*DVect2(oldCm->sn_M_.y(),oldCm->sn_M_.z());
2129            } else {
2130                m.get(0,0) = 1.;
2131                m.get(0,1) = 0.;
2132                m.get(1,0) = 0.;
2133                m.get(1,1) = 1.;
2134                tpf = DVect2(oldCm->sn_F_.y(),oldCm->sn_F_.z());
2135                tpm = DVect2(oldCm->sn_M_.y(),oldCm->sn_M_.z());
2136            }
2137            DVect pforce = DVect(0,tpf.x(),tpf.y());
2138            //DVect pm     = DVect(0,tpm.x(),tpm.y());
2139#else
2140            oldSystem;
2141            newSystem;
2142            DVect pforce = DVect(0,oldCm->sn_F_.y());
2143            //DVect pm     = DVect(0,oldCm->sn_M_.y());
2144#endif
2145            for (int i=1; i<dim; ++i)
2146                sn_F_.rdof(i) += pforce.dof(i);
2147            if (sn_mode_ && oldCm->sn_mode_)
2148                sn_F_.rx() = oldCm->sn_F_.x();
2149            oldCm->sn_F_ = DVect(0.0);
2150            oldCm->sn_M_ = DAVect(0.0);
2151            if (dpProps_ && oldCm->dpProps_) {
2152#ifdef THREED
2153                tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
2154                pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
2155#else
2156                pforce = oldCm->dpProps_->dp_F_;
2157#endif
2158                dpProps_->dp_F_ += pforce;
2159                oldCm->dpProps_->dp_F_ = DVect(0.0);
2160            }
2161            if(oldCm->getEnergyActivated()) {
2162                activateEnergy();
2163                energies_->estrain_ = oldCm->energies_->estrain_;
2164                energies_->edashpot_ = oldCm->energies_->edashpot_;
2165                energies_->eslip_ = oldCm->energies_->eslip_;
2166                oldCm->energies_->estrain_ = 0.0;
2167                oldCm->energies_->edashpot_ = 0.0;
2168                oldCm->energies_->eslip_ = 0.0;
2169            }
2170            rgap_ = oldCm->rgap_;
2171        }
2172        assert(sn_F_ == sn_F_);
2173    }
2174
2175    void ContactModelRBSN::setNonForcePropsFrom(IContactModel *old) {
2176        // Only called for contacts with wall facets when the wall resolution scheme
2177        // is set to full!
2178        // Only do something if the contact model is of the same type
2179        if (old->getName().compare("springnetwork",Qt::CaseInsensitive) == 0 && !isBonded()) {
2180            ContactModelRBSN *oldCm = (ContactModelRBSN *)old;
2181
2182            fictForce_ = oldCm->fictForce_;
2183            sn_F_ = oldCm->sn_F_;
2184            sn_sdisp_ = oldCm->sn_sdisp_;
2185            sn_M_ = oldCm->sn_M_;
2186            kRot_ = oldCm->kRot_;
2187            kTran_ = oldCm->kTran_;
2188            kRatio_ = oldCm->kRatio_;
2189            E_ = oldCm->E_;
2190            poisson_ = oldCm->poisson_;
2191            fric_ = oldCm->fric_;
2192            sn_bmul_ = oldCm->sn_bmul_;
2193            sn_tmul_ = oldCm->sn_tmul_;
2194            sn_rmul_ = oldCm->sn_rmul_;
2195            userArea_ = oldCm->userArea_;
2196            rgap_ = oldCm->rgap_;
2197            sn_fa_ = oldCm->sn_fa_;
2198            sn_mcf_ = oldCm->sn_mcf_;
2199            sn_dil_ = oldCm->sn_dil_;
2200            sn_dilzero_ = oldCm->sn_dilzero_;
2201            transTen_ = oldCm->transTen_;
2202            sn_elong_ = oldCm->sn_elong_;
2203            sn_ndisp_ = oldCm->sn_ndisp_;
2204            sn_mode_ = oldCm->sn_mode_;
2205            sn_state_ = oldCm->sn_state_;
2206            poisOffDiag_ = oldCm->poisOffDiag_;
2207            sn_S_ = oldCm->sn_S_;
2208            sn_BS_ = oldCm->sn_BS_;
2209            sn_TS_ = oldCm->sn_TS_;
2210            forceSet_ = oldCm->forceSet_;
2211            sn_heal_ = oldCm->sn_heal_;
2212            tenTable_ = oldCm->tenTable_;
2213            cohTable_ = oldCm->cohTable_;
2214            if (oldCm->hasDamping()) {
2215                if (!dpProps_)
2216                    dpProps_ = NEWC(dpProps());
2217                dp_nratio(oldCm->dp_nratio());
2218                dp_sratio(oldCm->dp_sratio());
2219                dp_mode(oldCm->dp_mode());
2220                dp_F(oldCm->dp_F());
2221            }
2222        }
2223    }
2224
2225    DVect ContactModelRBSN::getForce(const IContactMechanical *) const {
2226        DVect ret(sn_F_);
2227        ret += fictForce_;
2228        if (dpProps_)
2229            ret += dpProps_->dp_F_;
2230        return ret;
2231    }
2232
2233    DAVect ContactModelRBSN::getMomentOn1(const IContactMechanical *c) const {
2234        DVect force = getForce(c);
2235        DAVect ret(sn_M_);
2236        c->updateResultingTorqueOn1Local(force,&ret);
2237        return ret;
2238    }
2239
2240    DAVect ContactModelRBSN::getMomentOn2(const IContactMechanical *c) const {
2241        DVect force = getForce(c);
2242        DAVect ret(sn_M_);
2243        c->updateResultingTorqueOn2Local(force,&ret);
2244        return ret;
2245    }
2246
2247    DVect3 ContactModelRBSN::computeGeomData(const IContactMechanical *c) const {
2248        double Cmax1 = c->getEnd1Curvature().y();
2249        double Cmax2 = c->getEnd2Curvature().y();
2250        double br = sn_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
2251        if (userArea_)
2252#ifdef THREED
2253            br = std::sqrt(userArea_ / dPi);
2254#else
2255            br = userArea_ / 2.0;
2256#endif
2257        double br2 = br * br;
2258#ifdef THREED
2259        double area = dPi * br2;
2260        double bi = 0.25*area*br2;
2261#else
2262        double area = 2.0*br;
2263        double bi = 2.0*br*br2 / 3.0;
2264#endif
2265        return DVect3(area, bi, br);
2266    }
2267
2268    DVect2 ContactModelRBSN::SMax(const IContactMechanical *c) const {
2269        DVect3 data = computeGeomData(c);
2270        double area = data.x();
2271        double bi = data.y();
2272        double br = data.z();
2273        /* maximum stresses */
2274        //double nsmax0 = -(totalForce.x() / area) + sn_mcf_* sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z()) * br / bi;
2275        double dbend = sqrt(sn_M_.y()*sn_M_.y() + sn_M_.z()*sn_M_.z());
2276        dbend *= sn_mcf_;
2277        double dtwist = sn_M_.x();
2278        DVect bfs(sn_F_);
2279        bfs.rx() = 0.0;
2280        double dbfs = bfs.mag();
2281        double nsmax = -((sn_F_.x()+fictForce_.x()) / area) + dbend * br / bi;
2282        double ssmax = dbfs / area + std::abs(dtwist) * 0.5* br / bi;
2283        return DVect2(nsmax, ssmax);
2284    }
2285
2286    double ContactModelRBSN::shearStrength(const double &area) const {
2287        double sig = -1.0*(sn_F_.x() + fictForce_.x()) / area;
2288        double nstr = (sn_state_ > 2 && sn_state_ != 6) ? tenTable_[0].x() : 0.0;
2289        return sig <= nstr ? cohTable_[0].x() - sn_fa_*sig
2290            : cohTable_[0].x() - sn_fa_*nstr;
2291    }
2292
2293
2294    double ContactModelRBSN::strainEnergy(double kn,double ks,double kb,double kt) const {
2295        double ret(0.0);
2296        if (kn)
2297            ret = 0.5 * (sn_F_.x()+fictForce_.x()) * (sn_F_.x()+fictForce_.x()) / kn;
2298        if (ks) {
2299            DVect tmp = sn_F_ + fictForce_;
2300            tmp.rx() = 0.0;
2301            double smag2 = tmp.mag2();
2302            ret += 0.5 * smag2 / ks;
2303        }
2304
2305        if (kt)
2306            ret += 0.5 * sn_M_.x() * sn_M_.x() / kt;
2307        if (kb) {
2308            DAVect tmp = sn_M_;
2309#ifdef THREED
2310            tmp.rx() = 0.0;
2311            double smag2 = tmp.mag2();
2312#else
2313            double smag2 = tmp.z() * tmp.z();
2314#endif
2315            ret += 0.5 * smag2 / kb;
2316        }
2317        return ret;
2318    }
2319
2320} // namespace cmodelsxd
2321// EoF

Top