Soft-Bond Model Implementation

See this page for the documentation of this contact model.

contactmodelsoftbond.h

  1#pragma once
  2// contactmodelsoftbond.h
  3
  4#include "contactmodel/src/contactmodelmechanical.h"
  5
  6#ifdef SOFTBOND_LIB
  7#  define SOFTBOND_EXPORT EXPORT_TAG
  8#elif defined(NO_MODEL_IMPORT)
  9#  define SOFTBOND_EXPORT
 10#else
 11#  define SOFTBOND_EXPORT IMPORT_TAG
 12#endif
 13
 14namespace cmodelsxd {
 15    using namespace itasca;
 16
 17    class ContactModelSoftBond : public ContactModelMechanical {
 18    public:
 19        // Constructor: Set default values for contact model properties.
 20        SOFTBOND_EXPORT ContactModelSoftBond();
 21        // Destructor, called when contact is deleted: free allocated memory, etc.
 22        SOFTBOND_EXPORT virtual ~ContactModelSoftBond();
 23        // Contact model name (used as keyword for commands and FISH).
 24        QString  getName() const override { return "softbond"; }
 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            , kwFric 
 45            , kwBMul
 46            , kwTMul
 47            , kwSBMode
 48            , kwSBF
 49            , kwSBM
 50            , kwSBS
 51            , kwSBBS
 52            , kwSBTS
 53            , kwSBRMul
 54            , kwSBRadius
 55            , kwEmod
 56            , kwKRatio
 57            , kwDpNRatio 
 58            , kwDpSRatio
 59            , kwDpMode 
 60            , kwDpF
 61            , kwSBState
 62            , kwSBTStr
 63            , kwSBSStr
 64            , kwSBCoh
 65            , kwSBFa
 66            , kwSBMCF
 67            , kwSBSig
 68            , kwSBTau
 69            , kwSBSoft
 70            , kwSBCut
 71            , kwSBArea
 72            , kwUserArea
 73            , kwRGap
 74        };
 75        // Contact model property names in a comma separated list. The order corresponds with
 76        // the order of the PropertyKeys enumerator above. One can visualize any of these 
 77        // properties in PFC automatically. 
 78        QString  getProperties() const override {
 79            return "kn"
 80                   ",ks"
 81                   ",fric"
 82                   ",sb_bmul"
 83                   ",sb_tmul"
 84                   ",sb_mode"
 85                   ",sb_force"
 86                   ",sb_moment"
 87                   ",sb_slip"
 88                   ",sb_slipb"
 89                   ",sb_slipt"
 90                   ",sb_rmul"
 91                   ",sb_radius"
 92                   ",emod"
 93                   ",kratio"
 94                   ",dp_nratio"
 95                   ",dp_sratio"
 96                   ",dp_mode"
 97                   ",dp_force"
 98                   ",sb_state"
 99                   ",sb_ten"
100                   ",sb_shear"
101                   ",sb_coh"
102                   ",sb_fa"
103                   ",sb_mcf"
104                   ",sb_sigma"
105                   ",sb_tau"
106                   ",sb_soft"
107                   ",sb_cut"
108                   ",sb_area"
109                   ",user_area"
110                   ",rgap"
111                ;
112        }
113        // Enumerator for the energies.
114        enum EnergyKeys { 
115            kwEStrain=1
116          , kwESlip
117          , kwEDashpot
118        };
119        // Contact model energy names in a comma separated list. The order corresponds with
120        // the order of the EnergyKeys enumerator above. 
121        QString  getEnergies() const override {
122            return "energy-strain"
123                   ",energy-slip"
124                   ",energy-dashpot";
125        }
126        // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
127        double   getEnergy(uint32 i) const override;
128        // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1) 
129        // returns wther or not the estrain energy is accumulated which is false).
130        bool     getEnergyAccumulate(uint32 i) const override;
131        // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
132        void     setEnergy(uint32 i,const double &d) override; // Base 1
133        // Activate the energy. This is only called if the energy tracking is enabled. 
134        void     activateEnergy() override { if (energies_) return; energies_ = NEW Energies();}
135        // Returns whether or not the energy tracking has been enabled for this contact.
136        bool     getEnergyActivated() const override {return (energies_ != 0);}
137
138        // Enumerator for contact model related FISH callback events. 
139        enum FishCallEvents {
140             fActivated=0
141            ,fSlipChange
142            ,fBondBreak
143        };
144        // Contact model FISH callback event names in a comma separated list. The order corresponds with
145        // the order of the FishCallEvents enumerator above. 
146        QString  getFishCallEvents() const override {
147            return 
148                "contact_activated"
149                ",slip_change"
150                ",bond_break"; 
151        }
152
153        // Return the specified contact model property.
154        QVariant getProperty(uint32 i,const IContact *) const override;
155        // The return value denotes whether or not the property corresponds to the global
156        // or local coordinate system (TRUE: global system, FALSE: local system). The
157        // local system is the contact-plane system (nst) defined as follows.
158        // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
159        // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
160        // and tc are unit vectors in directions of the nst axes.
161        // This is used when rendering contact model properties that are vectors.
162        bool     getPropertyGlobal(uint32 i) const override;
163        // Set the specified contact model property, ensuring that it is of the correct type
164        // and within the correct range --- if not, then throw an exception.
165        // The return value denotes whether or not the update has affected the timestep
166        // computation (by having modified the translational or rotational tangent stiffnesses).
167        // If true is returned, then the timestep will be recomputed.
168        bool     setProperty(uint32 i,const QVariant &v,IContact *) override;
169        // The return value denotes whether or not the property is read-only
170        // (TRUE: read-only, FALSE: read-write).
171        bool     getPropertyReadOnly(uint32 i) const override;
172
173        // The return value denotes whether or not the property is inheritable
174        // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
175        // the endPropertyUpdated method.
176        bool     supportsInheritance(uint32 i) const override;
177        // Return whether or not inheritance is enabled for the specified property.
178        bool     getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
179        // Set the inheritance flag for the specified property.
180        void     setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }
181
182        // Enumerator for contact model methods.
183        enum MethodKeys { kwDeformability=1, kwBond, kwUnbond, kwArea};
184        // Contact model methoid names in a comma separated list. The order corresponds with
185        // the order of the MethodKeys enumerator above.  
186        QString  getMethods() const override { return "deformability,bond,unbond,area";}
187        // Return a comma seprated list of the contact model method arguments (base 1).
188        QString  getMethodArguments(uint32 i) const override;
189        // Set contact model method arguments (base 1). 
190        // The return value denotes whether or not the update has affected the timestep
191        // computation (by having modified the translational or rotational tangent stiffnesses).
192        // If true is returned, then the timestep will be recomputed.
193        bool     setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override;
194
195        // Prepare for entry into ForceDispLaw. The validate function is called when:
196        // (1) the contact is created, (2) a property of the contact that returns a true via
197        // the setProperty method has been modified and (3) when a set of cycles is executed
198        // via the {cycle N} command.
199        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
200        bool    validate(ContactModelMechanicalState *state,const double &timestep) override;
201        // The endPropertyUpdated method is called whenever a surface property (with a name
202        // that matches an inheritable contact model property name) of one of the contacting
203        // pieces is modified. This allows the contact model to update its associated
204        // properties. The return value denotes whether or not the update has affected
205        // the time step computation (by having modified the translational or rotational
206        // tangent stiffnesses). If true is returned, then the time step will be recomputed.  
207        bool    endPropertyUpdated(const QString &name,const IContactMechanical *c) override;
208        // The forceDisplacementLaw function is called during each cycle. Given the relative
209        // motion of the two contacting pieces (via
210        //   state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
211        //   state->relativeAngularIncrement_       (Dtt, Dtbs, Dtbt)
212        //     Ddn  : relative normal-displacement increment, Ddn > 0 is opening
213        //     Ddss : relative  shear-displacement increment (s-axis component)
214        //     Ddst : relative  shear-displacement increment (t-axis component)
215        //     Dtt  : relative twist-rotation increment
216        //     Dtbs : relative  bend-rotation increment (s-axis component)
217        //     Dtbt : relative  bend-rotation increment (t-axis component)
218        //       The relative displacement and rotation increments:
219        //         Dd = Ddn*nc + Ddss*sc + Ddst*tc
220        //         Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
221        //       where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
222        //       [see {Table 1: Contact State Variables} in PFC Model Components:
223        //       Contacts and Contact Models: Contact Resolution]
224        // ) and the contact properties, this function must update the contact force and
225        // moment.
226        //   The force_ is acting on piece 2, and is expressed in the local coordinate system
227        //   (defined in getPropertyGlobal) such that the first component positive denotes
228        //   compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
229        //   in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to 
230        //   get the total moment. 
231        // The return value indicates the contact activity status (TRUE: active, FALSE:
232        // inactive) during the next cycle.
233        // Additional information:
234        //   * If state->activated() is true, then the contact has just become active (it was
235        //     inactive during the previous time step).
236        //   * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
237        //     the forceDispLaw handle the case of {state->canFail_ == true}.
238        bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) override;
239        bool    thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
240        // The getEffectiveXStiffness functions return the translational and rotational
241        // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
242        // the translational tangent shear stiffness is zero (but this stiffness reduction
243        // is typically ignored when computing a stable time step). If the contact model
244        // includes a dashpot, then the translational stiffnesses must be increased (see
245        // Potyondy (2009)).
246        //   [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
247        //   Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
248        //   December 7, 2009.]
249        DVect2  getEffectiveTranslationalStiffness() const override { return effectiveTranslationalStiffness_; }
250        DAVect  getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness_;}
251
252        // Return a new instance of the contact model. This is used in the CMAT
253        // when a new contact is created. 
254        ContactModelSoftBond *clone() const override { return NEW ContactModelSoftBond(); }
255        // The getActivityDistance function is called by the contact-resolution logic when
256        // the CMAT is modified. Return value is the activity distance used by the
257        // checkActivity function.
258        double              getActivityDistance() const override {return rgap_;}
259        // The isOKToDelete function is called by the contact-resolution logic when...
260        // Return value indicates whether or not the contact may be deleted.
261        // If TRUE, then the contact may be deleted when it is inactive.
262        // If FALSE, then the contact may not be deleted (under any condition).
263        bool                isOKToDelete() const override { return !isBonded(); }
264        // Zero the forces and moments stored in the contact model. This function is called
265        // when the contact becomes inactive.
266        void                resetForcesAndMoments() override {
267            sb_F(DVect(0.0)); 
268            dp_F(DVect(0.0)); 
269            sb_M(DAVect(0.0));
270            if (energies_) {
271                energies_->estrain_ = 0.0;
272            }
273        }
274        void     setForce(const DVect &v,IContact *c) override;
275        void     setArea(const double &d) override { userArea_ = d; }
276        double   getArea() const override { return userArea_; }
277
278        // The checkActivity function is called by the contact-resolution logic when...
279        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
280        bool     checkActivity(const double &gap) override { return  gap <= rgap_ || isBonded();}
281
282        // Returns the sliding state (FALSE is returned if not implemented).
283        bool     isSliding() const override { return sb_S_; }
284        // Returns the bonding state (FALSE is returned if not implemented).
285        bool     isBonded() const override { return bProps_ ? (bProps_->sb_state_ >= 3) : false; }
286        void     unbond() override { if (bProps_) bProps_->sb_state_ = 0; }
287
288        // Both of these methods are called only for contacts with facets where the wall 
289        // resolution scheme is set the full. In such cases one might wish to propagate 
290        // contact state information (e.g., shear force) from one active contact to another. 
291        // See the Faceted Wall section in the documentation. 
292        void     propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
293        void     setNonForcePropsFrom(IContactModel *oldCM) override;   
294        /// Return the total force that the contact model holds.
295        DVect    getForce(const IContactMechanical *) const override;
296        /// Return the total moment on 1 that the contact model holds
297        DAVect   getMomentOn1(const IContactMechanical *) const override;
298        /// Return the total moment on 1 that the contact model holds
299        DAVect   getMomentOn2(const IContactMechanical *) const override;
300
301        // Methods to get and set properties. 
302        const double & kn() const {return kn_;}
303        void           kn(const double &d) {kn_=d;}
304        const double & ks() const {return ks_;}
305        void           ks(const double &d) {ks_=d;}
306        const double & fric() const {return fric_;}
307        void           fric(const double &d) {fric_=d;}
308        const double & sb_bmul() const { return   sb_bmul_; }
309        void           sb_bmul(const double &d) { sb_bmul_ = d; }
310        const double & sb_tmul() const { return   sb_tmul_; }
311        void           sb_tmul(const double &d) { sb_tmul_ = d; }
312        const DVect &  sb_F() const {return sb_F_;}
313        void           sb_F(const DVect &f) { sb_F_=f;}
314        const DAVect & sb_M() const { return sb_M_; }
315        void           sb_M(const DAVect &f) { sb_M_ = f; }
316        bool           sb_S() const {return sb_S_;}
317        void           sb_S(bool b) { sb_S_=b;}
318        bool           sb_BS() const { return sb_BS_; }
319        void           sb_BS(bool b) { sb_BS_ = b; }
320        bool           sb_TS() const { return sb_TS_; }
321        void           sb_TS(bool b) { sb_TS_ = b; }
322        const double & sb_rmul() const { return   sb_rmul_;        }
323        void           sb_rmul(const double &d) { sb_rmul_ = d; }
324        uint32           sb_mode() const {return sb_mode_;}
325        void           sb_mode(uint32 i) { sb_mode_=i;}
326
327        bool           hasBond() const { return bProps_ ? true : false; }
328        int            sb_state()   const { return (hasBond() ? bProps_->sb_state_ : 0); }
329        void           sb_state(int i) { if (!hasBond()) return; bProps_->sb_state_ = i; }
330        double         sb_Ten() const { return (hasBond() ? (bProps_->sb_ten_) : 0.0); }
331        void           sb_Ten(const double &d) { if (!hasBond()) return; bProps_->sb_ten_ = d; }
332        double         sb_Coh() const { return (hasBond() ? (bProps_->sb_coh_) : 0.0); }
333        void           sb_Coh(const double &d) { if (!hasBond()) return; bProps_->sb_coh_ = d; }
334        double         sb_FA() const { return (hasBond() ? (bProps_->sb_fa_) : 0.0); }
335        void           sb_FA(const double &d) { if (!hasBond()) return; bProps_->sb_fa_ = d; }
336        double         sb_MCF() const {return (hasBond() ? (bProps_->sb_mcf_) : 0.0);}
337        void           sb_MCF(const double &d) { if(!hasBond()) return; bProps_->sb_mcf_=d;}
338        double         sb_soft() const { return (hasBond() ? (bProps_->sb_soft_) : 0.0); }
339        void           sb_soft(const double &d) { if (!hasBond()) return; bProps_->sb_soft_ = d; }
340        double         sb_cut() const { return (hasBond() ? (bProps_->sb_cut_) : 0.0); }
341        void           sb_cut(const double &d) { if (!hasBond()) return; bProps_->sb_cut_ = d; }
342        double         sb_maxTen() const { return (hasBond() ? (bProps_->sb_maxTen_) : 0.0); }
343        void           sb_maxTen(const double &d) { if (!hasBond()) return; bProps_->sb_maxTen_ = d; }
344        double         sb_delu() const { return (hasBond() ? (bProps_->sb_delu_) : 0.0); }
345        void           sb_delu(const double &d) { if (!hasBond()) return; bProps_->sb_delu_ = d; }
346        Quat           sb_delo() const { return (hasBond() ? (bProps_->sb_delo_) : Quat::identity()); }
347        void           sb_delo(const Quat &d) { if (!hasBond()) return; bProps_->sb_delo_ = d; }
348        double         sb_maxu() const { return (hasBond() ? (bProps_->sb_maxu_) : 0.0); }
349        void           sb_maxu(const double &d) { if (!hasBond()) return; bProps_->sb_maxu_ = d; }
350        double         sb_critu() const { return (hasBond() ? (bProps_->sb_critu_) : 0.0); }
351        void           sb_critu(const double &d) { if (!hasBond()) return; bProps_->sb_critu_ = d; }
352
353        bool     hasDamping() const {return dpProps_ ? true : false;}
354        double   dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
355        void     dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
356        double   dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
357        void     dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
358        int      dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
359        void     dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
360        DVect    dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
361        void     dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
362
363        bool    hasEnergies() const {return energies_ ? true:false;}
364        double  estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
365        void    estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
366        double  eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
367        void    eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
368        double  edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
369        void    edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
370
371        uint32 inheritanceField() const {return inheritanceField_;}
372        void inheritanceField(uint32 i) {inheritanceField_ = i;}
373
374        const DVect2 & effectiveTranslationalStiffness()  const          {return effectiveTranslationalStiffness_;}
375        void           effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
376        const DAVect & effectiveRotationalStiffness()  const             {return effectiveRotationalStiffness_;}
377        void           effectiveRotationalStiffness(const DAVect &v )    {effectiveRotationalStiffness_=v;}
378
379    private:
380        // Index - used internally by PFC. Should be set to -1 in the cpp file. 
381        static int index_;
382
383        bool  FDLawBonded(ContactModelMechanicalState *state, const double &timestep);
384        bool  FDLawUnBonded(ContactModelMechanicalState *state, const double &timestep);
385
386        // Structure to compute stiffness
387        struct StiffData {
388            DVect2 trans_ = DVect2(0.0);
389            DAVect ang_   = DAVect(0.0);
390            double reff_ = 0.0;
391        };
392
393        // Structure to store the energies. 
394        struct Energies {
395            double estrain_  = 0.0;   // elastic energy  
396            double eslip_    = 0.0;   // work dissipated by friction 
397            double edashpot_ = 0.0;   // work dissipated by dashpots
398        };
399
400        // Structure to store dashpot quantities. 
401        struct dpProps {
402            double dp_nratio_ = 0.0;         // normal viscous critical damping ratio
403            double dp_sratio_ = 0.0;         // shear  viscous critical damping ratio
404            int    dp_mode_   = 0;           // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
405            DVect  dp_F_      = DVect(0.0);  // Force in the dashpots
406        };
407
408        // Structure to store bond-related quantities. 
409        struct bProps {
410            int     sb_state_    = 0;                // bond mode - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (B), 4 (B-Softening), 5 (B-Compression from Softening)
411            double  sb_ten_      = 0.0;              // normal strength 
412            double  sb_coh_      = 0.0;              // cohesion
413            double  sb_fa_       = 0.0;              // friction angle
414            double  sb_mcf_      = 1.0;              // moment contribution factor
415            double  sb_soft_     = 0.0;              // softening factor 
416            double  sb_cut_      = 1.0;              // critical bond length
417            double  sb_maxTen_   = 0.0;              // tensile strength one needs to reach for softening
418            double  sb_delu_     = 0.0;              // incremental elongation in softening
419            Quat    sb_delo_     = Quat::identity(); // incremental orientation in softening
420            double  sb_maxu_     = 0.0;              // max elongation for softening
421            double  sb_critu_    = 0.0;              // critical elongation for softening
422        };
423
424
425        bool   updateKn(const IContactMechanical *con);
426        bool   updateKs(const IContactMechanical *con);
427        bool   updateFric(const IContactMechanical *con);
428
429        StiffData computeStiffData(ContactModelMechanicalState *state) const;
430        DVect3    computeGeomData(const IContactMechanical *c) const;
431        DVect2    SMax(const IContactMechanical *con) const; // Maximum stress (tensile,shear) at bond periphery
432        double    shearStrength(const double &pbArea) const;      // Bond shear strength
433        double    strainEnergy(double kn, double ks, double kb, double kt) const;
434
435        void      updateStiffness(ContactModelMechanicalState *state);
436
437        // Contact model inheritance fields.
438        uint32 inheritanceField_;
439
440        // Effective translational stiffness.
441        DVect2  effectiveTranslationalStiffness_;
442        DAVect  effectiveRotationalStiffness_;      // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)
443
444        // linear model properties
445        double      kn_;          // Normal stiffness
446        double      ks_;          // Shear stiffness
447        double      fric_;        // Coulomb friction coefficient
448        double      sb_bmul_;     // Bending friction multiplier
449        double      sb_tmul_;     // Twisting friction  multiplier
450        uint32        sb_mode_;     // specifies absolute (0) or incremental (1) behavior for the the normal force 
451        DVect       sb_F_;        // Force carried in the model
452        DAVect      sb_M_;        // moment (bending + twisting in 3D)         
453        bool        sb_S_;        // The current slip state
454        bool        sb_BS_;       // The bending  slip state
455        bool        sb_TS_;       // The twisting slip state
456        double      sb_rmul_;     // Radius multiplier
457        double      userArea_;    // Area as specified by the user 
458        double      rgap_;        // Reference gap
459
460        dpProps *   dpProps_;   // The viscous properties
461        bProps *    bProps_;     // The bond properties
462
463        Energies *   energies_; // The energies
464
465    };
466} // namespace cmodelsxd
467// EoF

Top

contactmodelsoftbond.cpp

   1// contactmodelsoftbond.cpp
   2#include "contactmodelsoftbond.h"
   3
   4#include "module/interface/icontactmechanical.h"
   5#include "module/interface/icontact.h"
   6#include "module/interface/ipiecemechanical.h"
   7#include "module/interface/ipiece.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 "module/interface/icontactfluid.h"
  16#include "contactmodel/src/contactmodelthermal.h"
  17#include "contactmodel/src/contactmodelfluid.h"
  18#include "../version.txt"
  19#include "fish/src/parameter.h"
  20
  21#ifdef SOFTBOND_LIB
  22#ifdef _WIN32
  23    int __stdcall DllMain(void *,unsigned, void *) {
  24        return 1;
  25    }
  26#endif
  27    extern "C" EXPORT_TAG const char *getName() {
  28#if DIM==3
  29        return "contactmodelmechanical3dsoftbond";
  30#else
  31        return "contactmodelmechanical2dsoftbond";
  32#endif
  33    }
  34
  35    extern "C" EXPORT_TAG unsigned getMajorVersion() {
  36        return MAJOR_VERSION;
  37    }
  38
  39    extern "C" EXPORT_TAG unsigned getMinorVersion() {
  40        return MINOR_VERSION;
  41    }
  42
  43    extern "C" EXPORT_TAG void *createInstance() {
  44        cmodelsxd::ContactModelSoftBond *m = NEW cmodelsxd::ContactModelSoftBond();
  45        return (void *)m;
  46    }
  47#endif
  48
  49namespace cmodelsxd {
  50    static const uint32 KnMask      = 0x00000002; // Base 1!
  51    static const uint32 KsMask      = 0x00000004;
  52    static const uint32 FricMask    = 0x00000008;
  53
  54    using namespace itasca;
  55
  56    int ContactModelSoftBond::index_ = -1;
  57    uint32 ContactModelSoftBond::getMinorVersion() const { return MINOR_VERSION;  }
  58
  59    ContactModelSoftBond::ContactModelSoftBond() : inheritanceField_(KnMask|KsMask|FricMask)
  60                                             , effectiveTranslationalStiffness_(DVect2(0.0))
  61                                             , effectiveRotationalStiffness_(DAVect(0.0))
  62                                             , kn_(0.0)
  63                                             , ks_(0.0)
  64                                             , fric_(0.0)
  65                                             , sb_bmul_(1.0)
  66                                             , sb_tmul_(1.0)
  67                                             , sb_mode_(0)
  68                                             , sb_F_(DVect(0.0))
  69                                             , sb_M_(DAVect(0.0))
  70                                             , sb_S_(false)
  71                                             , sb_BS_(false)
  72                                             , sb_TS_(false)
  73                                             , sb_rmul_(1.0)
  74                                             , userArea_(0.0)
  75                                             , rgap_(0.0)
  76                                             , dpProps_(nullptr)
  77                                             , bProps_(nullptr)
  78                                             , energies_(nullptr) {
  79    }
  80
  81    ContactModelSoftBond::~ContactModelSoftBond() {
  82        // Make sure to clean up after yourself!
  83        if (dpProps_)
  84            delete dpProps_;
  85        if (bProps_)
  86            delete bProps_;
  87        if (energies_)
  88            delete energies_;
  89    }
  90
  91    void ContactModelSoftBond::archive(ArchiveStream &stream) {
  92        // The stream allows one to archive the values of the contact model
  93        // so that it can be saved and restored. The minor version can be
  94        // used here to allow for incremental changes to the contact model too.
  95        stream & kn_;
  96        stream & ks_;
  97        stream & fric_;
  98        stream & sb_mode_;
  99        stream & sb_F_;
 100        stream & sb_M_;
 101        stream & sb_S_;
 102        stream & sb_BS_;
 103        stream & sb_TS_;
 104        stream & sb_rmul_;
 105
 106        if (stream.getArchiveState()==ArchiveStream::Save) {
 107            bool b = false;
 108            if (dpProps_) {
 109                b = true;
 110                stream & b;
 111                stream & dpProps_->dp_nratio_;
 112                stream & dpProps_->dp_sratio_;
 113                stream & dpProps_->dp_mode_;
 114                stream & dpProps_->dp_F_;
 115            }
 116            else
 117                stream & b;
 118
 119            b = false;
 120            if (energies_) {
 121                b = true;
 122                stream & b;
 123                stream & energies_->estrain_;
 124                stream & energies_->eslip_;
 125                stream & energies_->edashpot_;
 126            }
 127            else
 128                stream & b;
 129
 130            b = false;
 131            if (bProps_) {
 132                b = true;
 133                stream & b;
 134                stream & bProps_->sb_state_;
 135                stream & bProps_->sb_ten_;
 136                stream & bProps_->sb_coh_;
 137                stream & bProps_->sb_fa_;
 138                stream & bProps_->sb_mcf_;
 139                stream & bProps_->sb_soft_;
 140                stream & bProps_->sb_cut_;
 141                stream & bProps_->sb_maxTen_;
 142                stream & bProps_->sb_delu_;
 143                stream & bProps_->sb_delo_;
 144                stream & bProps_->sb_maxu_;
 145                stream & bProps_->sb_critu_;
 146            }
 147            else
 148                stream & b;
 149
 150        } else {
 151            bool b(false);
 152            stream & b;
 153            if (b) {
 154                if (!dpProps_)
 155                    dpProps_ = NEW dpProps();
 156                stream & dpProps_->dp_nratio_;
 157                stream & dpProps_->dp_sratio_;
 158                stream & dpProps_->dp_mode_;
 159                stream & dpProps_->dp_F_;
 160            }
 161            stream & b;
 162            if (b) {
 163                if (!energies_)
 164                    energies_ = NEW Energies();
 165                stream & energies_->estrain_;
 166                stream & energies_->eslip_;
 167                stream & energies_->edashpot_;
 168            }
 169            stream & b;
 170            if (b) {
 171                if (!bProps_)
 172                    bProps_ = NEW bProps();
 173                stream & bProps_->sb_state_;
 174                stream & bProps_->sb_ten_;
 175                stream & bProps_->sb_coh_;
 176                stream & bProps_->sb_fa_;
 177                stream & bProps_->sb_mcf_;
 178                stream & bProps_->sb_soft_;
 179                stream & bProps_->sb_cut_;
 180                stream & bProps_->sb_maxTen_;
 181                stream & bProps_->sb_delu_;
 182                stream & bProps_->sb_delo_;
 183                stream & bProps_->sb_maxu_;
 184                stream & bProps_->sb_critu_;
 185            }
 186
 187        }
 188
 189        stream & inheritanceField_;
 190        stream & effectiveTranslationalStiffness_;
 191        stream & effectiveRotationalStiffness_;
 192
 193        if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 1) {
 194            stream & sb_bmul_;
 195            stream & sb_tmul_;
 196        }
 197
 198        if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 2)
 199            stream & userArea_;
 200
 201        if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 3)
 202            stream & rgap_;
 203
 204    }
 205
 206    void ContactModelSoftBond::copy(const ContactModel *cm) {
 207        // Copy all of the contact model properties. Used in the CMAT
 208        // when a new contact is created.
 209        ContactModelMechanical::copy(cm);
 210        const ContactModelSoftBond *in = dynamic_cast<const ContactModelSoftBond*>(cm);
 211        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
 212        kn(in->kn());
 213        ks(in->ks());
 214        fric(in->fric());
 215        sb_bmul(in->sb_bmul());
 216        sb_tmul(in->sb_tmul());
 217        sb_mode(in->sb_mode());
 218        sb_F(in->sb_F());
 219        sb_S(in->sb_S());
 220        sb_BS(in->sb_BS());
 221        sb_TS(in->sb_TS());
 222        sb_rmul(in->sb_rmul());
 223        sb_M(in->sb_M());
 224
 225        if (in->hasDamping()) {
 226            if (!dpProps_)
 227                dpProps_ = NEW dpProps();
 228            dp_nratio(in->dp_nratio());
 229            dp_sratio(in->dp_sratio());
 230            dp_mode(in->dp_mode());
 231            dp_F(in->dp_F());
 232        }
 233        if (in->hasEnergies()) {
 234            if (!energies_)
 235                energies_ = NEW Energies();
 236            estrain(in->estrain());
 237            eslip(in->eslip());
 238            edashpot(in->edashpot());
 239        }
 240        if (in->hasBond()) {
 241            if (!bProps_)
 242                bProps_ = NEW bProps();
 243            sb_state(in->sb_state());
 244            sb_Ten(in->sb_Ten());
 245            sb_Coh(in->sb_Coh());
 246            sb_FA(in->sb_FA());
 247            sb_MCF(in->sb_MCF());
 248            sb_soft(in->sb_soft());
 249            sb_cut(in->sb_cut());
 250            sb_maxTen(in->sb_maxTen());
 251            sb_delu(in->sb_delu());
 252            sb_delo(in->sb_delo());
 253            sb_maxu(in->sb_maxu());
 254            sb_critu(in->sb_critu());
 255        }
 256        userArea_ = in->userArea_;
 257        rgap_ = in->rgap_;
 258        inheritanceField(in->inheritanceField());
 259        effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
 260        effectiveRotationalStiffness(in->effectiveRotationalStiffness());
 261    }
 262
 263
 264    QVariant ContactModelSoftBond::getProperty(uint32 i,const IContact *con) const {
 265        // Return the property. The IContact pointer is provided so that
 266        // more complicated properties, depending on contact characteristics,
 267        // can be calcualted.
 268        QVariant var;
 269        switch (i) {
 270        case kwKn:       return kn_;
 271        case kwKs:       return ks_;
 272        case kwFric:     return fric_;
 273        case kwBMul:     return sb_bmul_;
 274        case kwTMul:     return sb_tmul_;
 275        case kwSBMode:   return sb_mode_;
 276        case kwSBF:      var.setValue(sb_F_); return var;
 277        case kwSBM:      var.setValue(sb_M_); return var;
 278        case kwSBS:      return sb_S_;
 279        case kwSBBS:     return sb_BS_;
 280        case kwSBTS:     return sb_TS_;
 281        case kwSBRMul:   return sb_rmul_;
 282        case kwSBRadius: {
 283            const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 284            if (!c) return 0.0;
 285            double Cmax1 = c->getEnd1Curvature().y();
 286            double Cmax2 = c->getEnd2Curvature().y();
 287            if (!userArea_)
 288                return sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
 289            else {
 290#ifdef THREED
 291                double rad = std::sqrt(userArea_ / dPi);
 292#else
 293                double rad = userArea_ / 2.0;
 294#endif
 295                return rad;
 296            }
 297
 298        }
 299        case kwEmod: {
 300                const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 301                if (!c) return 0.0;
 302                double rsum(0.0);
 303                if (c->getEnd1Curvature().y())
 304                    rsum += 1.0/c->getEnd1Curvature().y();
 305                if (c->getEnd2Curvature().y())
 306                    rsum += 1.0/c->getEnd2Curvature().y();
 307                if (userArea_)
 308#ifdef THREED
 309                    rsum = 2.0 * std::sqrt(userArea_ / dPi);
 310#else
 311                    rsum = userArea_;
 312#endif
 313                return kn_ * rsum;
 314        }
 315        case kwKRatio:    return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
 316        case kwDpNRatio:  return dpProps_ ? dpProps_->dp_nratio_ : 0;
 317        case kwDpSRatio:  return dpProps_ ? dpProps_->dp_sratio_ : 0;
 318        case kwDpMode:    return dpProps_ ? dpProps_->dp_mode_ : 0;
 319        case kwDpF: {
 320                dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
 321                return var;
 322            }
 323        case kwSBState:     return bProps_ ? bProps_->sb_state_ : 0;
 324        case kwSBTStr:      return bProps_ ? bProps_->sb_ten_ : 0.0;
 325        case kwSBSStr: {
 326            if (!bProps_) return 0.0;
 327            const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 328            double area = computeGeomData(c).x();
 329            return shearStrength(area);
 330        }
 331        case kwSBCoh:       return bProps_ ? bProps_->sb_coh_ : 0;
 332        case kwSBFa:        return bProps_ ? bProps_->sb_fa_ : 0;
 333        case kwSBMCF:       return bProps_ ? bProps_->sb_mcf_ : 0;
 334        case kwSBSig: {
 335            if (!bProps_ || bProps_->sb_state_ < 3) return 0.0;
 336            const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 337            return SMax(c).x();
 338        }
 339        case kwSBTau: {
 340            if (!bProps_ || bProps_->sb_state_ < 3) return 0.0;
 341            const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 342            return SMax(c).y();
 343        }
 344        case kwSBSoft:
 345            if (!bProps_) return 0.0;
 346            return bProps_->sb_soft_;
 347        case kwSBCut:
 348            if (!bProps_) return 0.0;
 349            return bProps_->sb_cut_;
 350        case kwSBArea: {
 351                if (userArea_) return userArea_;
 352                //if (!bProps_) return 0.0;
 353                const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 354                if (!c)
 355                    return 0.0;
 356                return computeGeomData(c).x();
 357            }
 358        case kwUserArea:
 359            return userArea_;
 360        case kwRGap:
 361            return rgap_;
 362        }
 363        assert(0);
 364        return QVariant();
 365    }
 366
 367    bool ContactModelSoftBond::getPropertyGlobal(uint32 i) const {
 368        // Returns whether or not a property is held in the global axis system (TRUE)
 369        // or the local system (FALSE). Used by the plotting logic.
 370        switch (i) {
 371        case kwSBF:
 372        case kwSBM:
 373        case kwDpF:
 374            return false;
 375        }
 376        return true;
 377    }
 378
 379    bool ContactModelSoftBond::setProperty(uint32 i,const QVariant &v,IContact *) {
 380        // Set a contact model property. Return value indicates that the timestep
 381        // should be recalculated.
 382        dpProps dp;
 383        switch (i) {
 384        case kwKn: {
 385                if (!v.canConvert<double>())
 386                    throw Exception("kn must be a double.");
 387                double val(v.toDouble());
 388                if (val<0.0)
 389                    throw Exception("Negative kn not allowed.");
 390                kn_ = val;
 391                return true;
 392            }
 393        case kwKs: {
 394                if (!v.canConvert<double>())
 395                    throw Exception("ks must be a double.");
 396                double val(v.toDouble());
 397                if (val<0.0)
 398                    throw Exception("Negative ks not allowed.");
 399                ks_ = val;
 400                return true;
 401            }
 402        case kwFric: {
 403                if (!v.canConvert<double>())
 404                    throw Exception("fric must be a double.");
 405                double val(v.toDouble());
 406                if (val<0.0)
 407                    throw Exception("Negative fric not allowed.");
 408                fric_ = val;
 409                return false;
 410            }
 411        case kwBMul: {
 412                if (!v.canConvert<double>())
 413                    throw Exception("sb_bmul must be a double.");
 414                double val(v.toDouble());
 415                if (val<0.0)
 416                    throw Exception("Negative sb_bmul not allowed.");
 417                sb_bmul_ = val;
 418                return false;
 419            }
 420        case kwTMul: {
 421                if (!v.canConvert<double>())
 422                    throw Exception("sb_tmul must be a double.");
 423                double val(v.toDouble());
 424                if (val<0.0)
 425                    throw Exception("Negative st_bmul not allowed.");
 426                sb_tmul_ = val;
 427                return false;
 428            }
 429        case kwSBMode: {
 430                if (!v.canConvert<uint32>())
 431                    throw Exception("sb_mode must be 0 (absolute) or 1 (incremental).");
 432                double val(v.toUInt());
 433                if (val>1)
 434                    throw Exception("sb_mode must be 0 (absolute) or 1 (incremental).");
 435                sb_mode_ = val;
 436                return false;
 437            }
 438        case kwSBRMul: {
 439                if (!v.canConvert<double>())
 440                    throw Exception("rmul must be a double.");
 441                double val(v.toDouble());
 442                if (val<0.0)
 443                    throw Exception("Negative rmul not allowed.");
 444                sb_rmul_ = val;
 445                return false;
 446            }
 447        case kwSBF: {
 448                if (!v.canConvert<DVect>())
 449                    throw Exception("sb_force must be a vector.");
 450                DVect val(v.value<DVect>());
 451                sb_F_ = val;
 452                return false;
 453            }
 454        case kwSBM: {
 455                DAVect val(0.0);
 456#ifdef TWOD
 457                if (!v.canConvert<DAVect>() && !v.canConvert<double>())
 458                    throw Exception("res_moment must be an angular vector.");
 459                if (v.canConvert<DAVect>())
 460                    val = DAVect(v.value<DAVect>());
 461                else
 462                    val = DAVect(v.toDouble());
 463#else
 464                if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
 465                    throw Exception("res_moment must be an angular vector.");
 466                if (v.canConvert<DAVect>())
 467                    val = DAVect(v.value<DAVect>());
 468                else
 469                    val = DAVect(v.value<DVect>());
 470#endif
 471                sb_M_ = val;
 472                return false;
 473            }
 474        case kwDpNRatio: {
 475                if (!v.canConvert<double>())
 476                    throw Exception("dp_nratio must be a double.");
 477                double val(v.toDouble());
 478                if (val<0.0)
 479                    throw Exception("Negative dp_nratio not allowed.");
 480                if (val == 0.0 && !dpProps_)
 481                    return false;
 482                if (!dpProps_)
 483                    dpProps_ = NEW dpProps();
 484                dpProps_->dp_nratio_ = val;
 485                return true;
 486            }
 487        case kwDpSRatio: {
 488                if (!v.canConvert<double>())
 489                    throw Exception("dp_sratio must be a double.");
 490                double val(v.toDouble());
 491                if (val<0.0)
 492                    throw Exception("Negative dp_sratio not allowed.");
 493                if (val == 0.0 && !dpProps_)
 494                    return false;
 495                if (!dpProps_)
 496                    dpProps_ = NEW dpProps();
 497                dpProps_->dp_sratio_ = val;
 498                return true;
 499            }
 500        case kwDpMode: {
 501                if (!v.canConvert<int>())
 502                    throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
 503                int val(v.toInt());
 504                if (val == 0 && !dpProps_)
 505                    return false;
 506                if (val < 0 || val > 3)
 507                    throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
 508                if (!dpProps_)
 509                    dpProps_ = NEW dpProps();
 510                dpProps_->dp_mode_ = val;
 511                return false;
 512            }
 513        case kwDpF: {
 514                if (!v.canConvert<DVect>())
 515                    throw Exception("dp_force must be a vector.");
 516                DVect val(v.value<DVect>());
 517                if (val.fsum() == 0.0 && !dpProps_)
 518                    return false;
 519                if (!dpProps_)
 520                    dpProps_ = NEW dpProps();
 521                dpProps_->dp_F_ = val;
 522                return false;
 523            }
 524        case kwSBTStr: {
 525                if (!v.canConvert<double>())
 526                    throw Exception("sb_ten must be a double.");
 527                double val(v.toDouble());
 528                if (val < 0.0)
 529                    throw Exception("Negative sb_ten not allowed.");
 530                if (val == 0.0 && !bProps_)
 531                    return false;
 532                if (!bProps_)
 533                    bProps_ = NEW bProps();
 534                bProps_->sb_ten_ = val;
 535                return false;
 536            }
 537        case kwSBCoh: {
 538                if (!v.canConvert<double>())
 539                    throw Exception("sb_coh must be a double.");
 540                double val(v.toDouble());
 541                if (val<0.0)
 542                    throw Exception("Negative pb_coh not allowed.");
 543                if (val == 0.0 && !bProps_)
 544                    return false;
 545                if (!bProps_)
 546                    bProps_ = NEW bProps();
 547                bProps_->sb_coh_ = val;
 548                return false;
 549            }
 550        case kwSBFa: {
 551                if (!v.canConvert<double>())
 552                    throw Exception("sb_fa must be a double.");
 553                double val(v.toDouble());
 554                if (val<0.0)
 555                    throw Exception("Negative sb_fa not allowed.");
 556                if (val >= 90.0)
 557                    throw Exception("sb_fa must be lower than 90.0 degrees.");
 558                if (val == 0.0 && !bProps_)
 559                    return false;
 560                if (!bProps_)
 561                    bProps_ = NEW bProps();
 562                bProps_->sb_fa_ = val;
 563                return false;
 564            }
 565        case kwSBMCF: {
 566                if (!v.canConvert<double>())
 567                    throw Exception("sb_mcf must be a double.");
 568                double val(v.toDouble());
 569                if (val<0.0)
 570                    throw Exception("Negative sb_mcf not allowed.");
 571                if (val > 1.0)
 572                    throw Exception("sb_mcf must be lower or equal to 1.0.");
 573                if (val == 1.0 && !bProps_)
 574                    return false;
 575                if (!bProps_)
 576                    bProps_ = NEW bProps();
 577                bProps_->sb_mcf_ = val;
 578                return false;
 579            }
 580        case kwSBSoft: {
 581                if (!v.canConvert<double>())
 582                    throw Exception("sb_soft must be a double.");
 583                double val(v.toDouble());
 584                if (val < 0.0)
 585                    throw Exception("Negative pb_soft not allowed.");
 586                if (!bProps_)
 587                    bProps_ = NEW bProps();
 588                bProps_->sb_soft_ = val;
 589                return false;
 590            }
 591        case kwSBCut: {
 592                if (!v.canConvert<double>())
 593                    throw Exception("sb_cut must be a double.");
 594                double val(v.toDouble());
 595                if (val < 0.0)
 596                    throw Exception("Negative sb_cut not allowed.");
 597                if (!bProps_)
 598                    bProps_ = NEW bProps();
 599                bProps_->sb_cut_ = val;
 600                return false;
 601            }
 602        case kwSBArea:
 603        case kwUserArea: {
 604                if (!v.canConvert<double>())
 605                    throw Exception("area must be a double.");
 606                double val(v.toDouble());
 607                if (val < 0.0)
 608                    throw Exception("Negative area not allowed.");
 609                userArea_ = val;
 610                return true;
 611            }
 612        case kwRGap: {
 613                if (!v.canConvert<double>())
 614                    throw Exception("Reference gap must be a double.");
 615                double val(v.toDouble());
 616                rgap_ = val;
 617                return false;
 618            }
 619        }
 620        return false;
 621    }
 622
 623    bool ContactModelSoftBond::getPropertyReadOnly(uint32 i) const {
 624        // Returns TRUE if a property is read only or FALSE otherwise.
 625        switch (i) {
 626        case kwDpF:
 627        case kwSBS:
 628        case kwSBBS:
 629        case kwSBTS:
 630        case kwEmod:
 631        case kwKRatio:
 632        case kwSBState:
 633        case kwSBRadius:
 634        case kwSBSStr:
 635        case kwSBSig:
 636        case kwSBTau:
 637            return true;
 638        default:
 639            break;
 640        }
 641        return false;
 642    }
 643
 644    bool ContactModelSoftBond::supportsInheritance(uint32 i) const {
 645        // Returns TRUE if a property supports inheritance or FALSE otherwise.
 646        switch (i) {
 647        case kwKn:
 648        case kwKs:
 649        case kwFric:
 650            return true;
 651        default:
 652            break;
 653        }
 654        return false;
 655    }
 656
 657    QString  ContactModelSoftBond::getMethodArguments(uint32 i) const {
 658        // Return a list of contact model method argument names.
 659        switch (i) {
 660        case kwDeformability:
 661            return "emod,kratio";
 662        case kwBond:
 663            return "gap,soft,cut";
 664        case kwUnbond:
 665            return "gap";
 666        case kwArea:
 667            return QString();
 668        }
 669        assert(0);
 670        return QString();
 671    }
 672
 673    bool ContactModelSoftBond::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
 674        // Apply the specified method.
 675        IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 676        switch (i) {
 677        case kwDeformability: {
 678                double emod;
 679                double krat;
 680                if (vl.at(0).isNull())
 681                    throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
 682                emod = vl.at(0).toDouble();
 683                if (emod<0.0)
 684                    throw Exception("Negative emod not allowed in contact model %1.",getName());
 685                if (vl.at(1).isNull())
 686                    throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
 687                krat = vl.at(1).toDouble();
 688                if (krat<0.0)
 689                    throw Exception("Negative stiffness ratio not allowed in contact model %1.",getName());
 690                double rsum(0.0);
 691                if (c->getEnd1Curvature().y())
 692                    rsum += 1.0 / c->getEnd1Curvature().y();
 693                if (c->getEnd2Curvature().y())
 694                    rsum += 1.0 / c->getEnd2Curvature().y();
 695                if (userArea_)
 696#ifdef THREED
 697                    rsum = 2.0 * std::sqrt(userArea_ / dPi);
 698#else
 699                    rsum = userArea_;
 700#endif
 701                kn_ = emod / rsum;
 702                ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
 703                setInheritance(1,false);
 704                setInheritance(2,false);
 705                return true;
 706            }
 707        case kwBond: {
 708                if (bProps_ && bProps_->sb_state_ >= 3) return false;
 709                double mingap = -1.0 * limits<double>::max();
 710                double maxgap = 0;
 711                if (vl.at(0).canConvert<double>())
 712                    maxgap = vl.at(0).toDouble();
 713                else if (vl.at(0).canConvert<DVect2>()) {
 714                    DVect2 value = vl.at(0).value<DVect2>();
 715                    mingap = value.minComp();
 716                    maxgap = value.maxComp();
 717                }
 718                else if (!vl.at(0).isNull())
 719                    throw Exception("gap value %1 not recognized in method bond in contact model %2.", vl.at(1), getName());
 720                double soft = bProps_ ? bProps_->sb_soft_ : 0.0;
 721                if (!vl.at(1).isNull()) {
 722                    soft = vl.at(1).toDouble();
 723                if (soft < 0.0)
 724                    throw Exception("Negative soft not allowed in contact model %1.", getName());
 725                }
 726                double cut = bProps_ ? bProps_->sb_cut_ : 1.0;
 727                if (!vl.at(2).isNull()) {
 728                    if (vl.at(2).canConvert<double>())
 729                        cut = vl.at(2).toDouble();
 730                    if (cut < 0.0)
 731                        throw Exception("cut value %1 is negative, or not recognized in method bond in contact model %2.", vl.at(2), getName());
 732                    if (cut > 1.0)
 733                        throw Exception("cut value %1 must be in range [0,1] in method bond in contact model %2.", vl.at(2), getName());
 734                }
 735                double gap = c->getGap();
 736                if (gap >= mingap && gap <= maxgap) {
 737                    if (!bProps_)
 738                        bProps_ = NEW bProps();
 739                    bProps_->sb_state_ = 3;
 740                    bProps_->sb_soft_ = soft;
 741                    // Update the critical distance
 742                    if (cut != -1)
 743                        bProps_->sb_cut_ = cut;
 744                    // seet to incremental normal force calculation
 745                    sb_mode_ = 1;
 746                    return true;
 747                }
 748                return false;
 749            }
 750        case kwUnbond: {
 751                if (!bProps_ || bProps_->sb_state_ == 0) return false;
 752                double mingap = -1.0 * limits<double>::max();
 753                double maxgap = 0;
 754                if (vl.at(0).canConvert<double>())
 755                    maxgap = vl.at(0).toDouble();
 756                else if (vl.at(0).canConvert<DVect2>()) {
 757                    DVect2 value = vl.at(0).value<DVect2>();
 758                    mingap = value.minComp();
 759                    maxgap = value.maxComp();
 760                }
 761                else if (!vl.at(0).isNull())
 762                    throw Exception("gap value %1 not recognized in method unbond in contact model %2.", vl.at(0), getName());
 763                double gap = c->getGap();
 764                if (gap >= mingap && gap <= maxgap) {
 765                    bProps_->sb_state_ = 0;
 766                    return true;
 767                }
 768                return false;
 769            }
 770        case kwArea: {
 771                if (!userArea_) {
 772                    double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
 773#ifdef THREED
 774                    userArea_ = rsq * rsq * dPi;
 775#else
 776                    userArea_ = rsq * 2.0;
 777#endif
 778                }
 779                return true;
 780            }
 781
 782        }
 783        return false;
 784    }
 785
 786    double ContactModelSoftBond::getEnergy(uint32 i) const {
 787        // Return an energy value.
 788        double ret(0.0);
 789        if (!energies_)
 790            return ret;
 791        switch (i) {
 792        case kwEStrain:    return energies_->estrain_;
 793        case kwESlip:      return energies_->eslip_;
 794        case kwEDashpot:   return energies_->edashpot_;
 795        }
 796        assert(0);
 797        return ret;
 798    }
 799
 800    bool ContactModelSoftBond::getEnergyAccumulate(uint32 i) const {
 801        // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
 802        switch (i) {
 803        case kwEStrain:   return false;
 804        case kwESlip:     return true;
 805        case kwEDashpot:  return true;
 806        }
 807        assert(0);
 808        return false;
 809    }
 810
 811    void ContactModelSoftBond::setEnergy(uint32 i,const double &d) {
 812        // Set an energy value.
 813        if (!energies_) return;
 814        switch (i) {
 815        case kwEStrain:    energies_->estrain_ = d;   return;
 816        case kwESlip:      energies_->eslip_   = d;   return;
 817        case kwEDashpot:   energies_->edashpot_= d;   return;
 818        }
 819        assert(0);
 820        return;
 821    }
 822
 823    bool ContactModelSoftBond::validate(ContactModelMechanicalState *state,const double &) {
 824        // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
 825        // (1) the contact is created, (2) a property of the contact that returns a true via
 826        // the setProperty method has been modified and (3) when a set of cycles is executed
 827        // via the {cycle N} command.
 828        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
 829        assert(state);
 830        const IContactMechanical *c = state->getMechanicalContact();
 831        assert(c);
 832
 833        if (state->trackEnergy_)
 834            activateEnergy();
 835
 836        if (inheritanceField_ & KnMask)
 837            updateKn(c);
 838        if (inheritanceField_ & KsMask)
 839            updateKs(c);
 840        if (inheritanceField_ & FricMask)
 841            updateFric(c);
 842
 843        updateStiffness(state);
 844        return checkActivity(state->gap_);
 845    }
 846
 847    static const QString knstr("kn");
 848    bool ContactModelSoftBond::updateKn(const IContactMechanical *con) {
 849        assert(con);
 850        QVariant v1 = con->getEnd1()->getProperty(knstr);
 851        QVariant v2 = con->getEnd2()->getProperty(knstr);
 852        if (!v1.isValid() || !v2.isValid())
 853            return false;
 854        double kn1 = v1.toDouble();
 855        double kn2 = v2.toDouble();
 856        double val = kn_;
 857        if (kn1 && kn2)
 858            kn_ = kn1*kn2/(kn1+kn2);
 859        else if (kn1)
 860            kn_ = kn1;
 861        else if (kn2)
 862            kn_ = kn2;
 863        return ( (kn_ != val) );
 864    }
 865
 866    static const QString ksstr("ks");
 867    bool ContactModelSoftBond::updateKs(const IContactMechanical *con) {
 868        assert(con);
 869        QVariant v1 = con->getEnd1()->getProperty(ksstr);
 870        QVariant v2 = con->getEnd2()->getProperty(ksstr);
 871        if (!v1.isValid() || !v2.isValid())
 872            return false;
 873        double ks1 = v1.toDouble();
 874        double ks2 = v2.toDouble();
 875        double val = ks_;
 876        if (ks1 && ks2)
 877            ks_ = ks1*ks2/(ks1+ks2);
 878        else if (ks1)
 879            ks_ = ks1;
 880        else if (ks2)
 881            ks_ = ks2;
 882        return ( (ks_ != val) );
 883    }
 884
 885    static const QString fricstr("fric");
 886    bool ContactModelSoftBond::updateFric(const IContactMechanical *con) {
 887        assert(con);
 888        QVariant v1 = con->getEnd1()->getProperty(fricstr);
 889        QVariant v2 = con->getEnd2()->getProperty(fricstr);
 890        if (!v1.isValid() || !v2.isValid())
 891            return false;
 892        double fric1 = std::max(0.0,v1.toDouble());
 893        double fric2 = std::max(0.0,v2.toDouble());
 894        double val = fric_;
 895        fric_ = std::min(fric1,fric2);
 896        return ( (fric_ != val) );
 897    }
 898
 899    bool ContactModelSoftBond::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
 900        // The endPropertyUpdated method is called whenever a surface property (with a name
 901        // that matches an inheritable contact model property name) of one of the contacting
 902        // pieces is modified. This allows the contact model to update its associated
 903        // properties. The return value denotes whether or not the update has affected
 904        // the time step computation (by having modified the translational or rotational
 905        // tangent stiffnesses). If true is returned, then the time step will be recomputed.
 906        assert(c);
 907        QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
 908        QRegularExpression rx(name, QRegularExpression::CaseInsensitiveOption);
 909        int idx = availableProperties.indexOf(rx)+1;
 910        bool ret=false;
 911
 912        if (idx<=0)
 913            return ret;
 914
 915        switch(idx) {
 916        case kwKn:  { //kn
 917                if (inheritanceField_ & KnMask)
 918                    ret = updateKn(c);
 919                break;
 920            }
 921        case kwKs:  { //ks
 922                if (inheritanceField_ & KsMask)
 923                    ret =updateKs(c);
 924                break;
 925            }
 926        case kwFric:  { //fric
 927                if (inheritanceField_ & FricMask)
 928                    updateFric(c);
 929                break;
 930            }
 931        }
 932        return ret;
 933    }
 934
 935    ContactModelSoftBond::StiffData ContactModelSoftBond::computeStiffData(ContactModelMechanicalState *state) const {
 936        // Update contact data
 937        double Cmin1 = state->end1Curvature_.x();
 938        double Cmax1 = state->end1Curvature_.y();
 939        double Cmax2 = state->end2Curvature_.y();
 940        double dthick = (Cmin1 == 0.0) ? 1.0 : 0.0;
 941        double br = sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
 942        if (userArea_)
 943#ifdef THREED
 944            br = std::sqrt(userArea_ / dPi);
 945#else
 946            br = userArea_ / 2.0;
 947#endif
 948        double br2 = br * br;
 949        double area = dthick <= 0.0 ? dPi * br2 : 2.0*br*dthick;
 950        double bi = dthick <= 0.0 ? 0.25*area*br2 : 2.0*br*br2*dthick / 3.0;
 951        StiffData ret;
 952        ret.reff_ = br;
 953        ret.trans_ = DVect2(kn_ * area , ks_ * area);
 954        ret.ang_ = DAVect(kn_ * bi);
 955#if DIM==3
 956        ret.ang_.rx() = ks_ * 2.0*bi;
 957#endif
 958        return ret;
 959    }
 960
 961    void ContactModelSoftBond::updateStiffness(ContactModelMechanicalState *state) {
 962        // first compute stiffness data
 963        StiffData stiff = computeStiffData(state);
 964        // Now calculate effective stiffness
 965        DVect2 retT = stiff.trans_;
 966        // correction if viscous damping active
 967        if (dpProps_) {
 968            DVect2 correct(1.0);
 969            if (dpProps_->dp_nratio_)
 970                correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
 971            if (dpProps_->dp_sratio_)
 972                correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
 973            retT /= (correct*correct);
 974        }
 975        effectiveTranslationalStiffness_ = retT;
 976        // Effective rotational stiffness (bending and twisting)
 977        effectiveRotationalStiffness_ = stiff.ang_;
 978    }
 979
 980    bool ContactModelSoftBond::forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) {
 981        assert(state);
 982
 983        if (state->activated()) {
 984            // The contact was just activated from an inactive state
 985            // Trigger the FISH callback if one is hooked up to the
 986            // contact_activated event.
 987            if (cmEvents_[fActivated] >= 0) {
 988                auto c = state->getContact();
 989                std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
 990                IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
 991                fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
 992            }
 993        }
 994        updateStiffness(state);
 995        if (isBonded()) return FDLawBonded(state, timestep);
 996        else return FDLawUnBonded(state, timestep);
 997
 998    }
 999
1000    bool ContactModelSoftBond::FDLawBonded(ContactModelMechanicalState *state, const double &timestep) {
1001        // Relative translational/rotational displacement increments
1002        DVect  trans = state->relativeTranslationalIncrement_;
1003        DAVect ang = state->relativeAngularIncrement_;
1004
1005        // Store previous force and moment
1006        DVect  sb_F_old = sb_F_;
1007        DAVect sb_M_old = sb_M_;
1008
1009        // Update stiffness data
1010        StiffData stiff = computeStiffData(state);
1011        DVect3 geom = computeGeomData(state->getMechanicalContact());
1012        double area = geom.x();
1013        double bi = geom.y();
1014        double br = geom.z();
1015        double kn = stiff.trans_.x();
1016        double ks = stiff.trans_.y();
1017        double kb = stiff.ang_.z();
1018#if DIM==3
1019        double kt = stiff.ang_.x();
1020#else
1021        double kt = 0.0;
1022#endif
1023
1024        double nsmax0 = -(sb_F_.x() / area) + bProps_->sb_mcf_* sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z()) * br / bi;
1025
1026        // incremental normal force calculation
1027        sb_F_.rx() -= trans.x() * kn;
1028
1029        // shear force calculation
1030        // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1031        // Loop over the shear components (note: the 0 component is the normal component)
1032        // and calculate the shear force.
1033        for (int i = 1; i<dim; ++i)
1034            sb_F_.rdof(i) -= trans.dof(i) * ks;
1035
1036        // moment calculation
1037        sb_M_ -= ang * stiff.ang_;
1038        double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1039
1040        // maximum tensile stress at bond periphery
1041        double nsmax = -(sb_F_.x() / area) + bProps_->sb_mcf_* dbend * br / bi;
1042
1043        bool softened = false;
1044        // Mode check
1045        if (state->canFail_) {
1046            if (bProps_->sb_state_ == 3 || bProps_->sb_state_ == 5) {
1047                double compVal = bProps_->sb_state_ == 3 ? bProps_->sb_ten_ : bProps_->sb_maxTen_;
1048                if (nsmax >= compVal ) {
1049                    // enter softening regime
1050                    // current bond elongation when softening starts
1051                    // This is the elongation at the bond periphery
1052                    double ls = - sb_F_.x() / kn + bProps_->sb_mcf_*dbend* br / kb;
1053                    bProps_->sb_maxTen_ = compVal;
1054                    if (bProps_->sb_state_ == 3)
1055                        bProps_->sb_critu_ = ls /**(1.0+bProps_->sb_soft_)*/;
1056                    bProps_->sb_delu_ = 0.0;
1057                    bProps_->sb_delo_ = Quat::identity();
1058                    if (bProps_->sb_state_ == 5 && nsmax < bProps_->sb_maxTen_)
1059                        softened = true;
1060                    bProps_->sb_state_ = 4;
1061                }
1062            }
1063        }
1064
1065        if (bProps_->sb_state_ == 4 && !softened && !checktol(bProps_->sb_soft_,0.0,1.0,100.0)) {
1066            double ls = bProps_->sb_critu_;
1067            double lc = ls * (1.0+bProps_->sb_soft_);
1068            DVect normal(0.0);
1069            normal.rx() = 1.0;
1070            DVect backNormal = (bProps_->sb_delo_.getConj().rotate(normal)).unit();
1071            double bend = acos(qBound(-1.0,normal|backNormal,1.0));
1072            double l0 = ls + bProps_->sb_maxu_ + bProps_->sb_delu_ + br*abs(bend);
1073            bProps_->sb_delu_ += trans.x();
1074            bProps_->sb_delo_.increment(ang);
1075            // Take the current contact normal and rotate it in the opposite direction of
1076            // the orientation - get the angle of bend from there
1077            backNormal = (bProps_->sb_delo_.getConj().rotate(normal)).unit();
1078            bend = acos(qBound(-1.0,normal|backNormal,1.0));
1079            double l = ls + bProps_->sb_maxu_ + bProps_->sb_delu_ + br*abs(bend);
1080            // target tensile stress
1081            double ns = bProps_->sb_ten_*(lc-l) / (bProps_->sb_soft_*ls);
1082            if (ns > 0) {
1083                if (nsmax >= ns) {
1084                    double fac = ns / nsmax;
1085                    sb_F_.rx() = fac*sb_F_.x();
1086#if DIM==3
1087                    sb_M_.ry() = fac*sb_M_.y();
1088#endif
1089                    sb_M_.rz() = fac*sb_M_.z();
1090                } else {
1091                    bProps_->sb_state_ = 5;
1092                    bProps_->sb_maxTen_ = nsmax0;
1093                    bProps_->sb_maxu_ = (l0-ls);
1094                }
1095            } else {
1096                sb_F_.rx() = 0.0;
1097#if DIM==3
1098                sb_M_.ry() = 0.0;
1099#endif
1100                sb_M_.rz() = 0.0;
1101            }
1102        }
1103
1104        if (state->canFail_) {
1105            /* check for normal failure */
1106            bool failed = false;
1107            if (bProps_->sb_state_ == 4) {
1108                double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1109                double nsmax = -(sb_F_.x() / area) + bProps_->sb_mcf_*dbend * br / bi;
1110                if (nsmax <= bProps_->sb_ten_*bProps_->sb_cut_ || checktol(bProps_->sb_soft_,0.0,1.0,100.0)) {
1111                    // Failed in tension
1112                    double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1113                    bProps_->sb_state_ = 1;
1114                    sb_F_.fill(0.0);
1115                    sb_M_.fill(0.0);
1116                    failed = true;
1117                    if (cmEvents_[fBondBreak] >= 0) {
1118                        auto c = state->getContact();
1119                        std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1120                                                             fish::Parameter((qint64)bProps_->sb_state_),
1121                                                             fish::Parameter(nsmax),
1122                                                             fish::Parameter(se)
1123                                                           };
1124                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1125                        fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1126                    }
1127                }
1128            }
1129
1130            if (!failed) {
1131                /* check for shear failure */
1132                double dtwist = sb_M_.x();
1133                DVect bfs(sb_F_);
1134                bfs.rx() = 0.0;
1135                double dbfs = bfs.mag();
1136                double ssmax = dbfs / area + bProps_->sb_mcf_*std::abs(dtwist) * 0.5* br / bi;
1137                double ss = shearStrength(area);
1138                if (ss < 0)
1139                    ss = 0;
1140                if (ss <= ssmax) {
1141                    // Failed in shear
1142                    double se = strainEnergy(kn, ks, kb, kt); // bond strain energy at the onset of failure
1143                    bProps_->sb_state_ = 2;
1144                    if (cmEvents_[fBondBreak] >= 0) {
1145                        auto c = state->getContact();
1146                        std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1147                                                             fish::Parameter((qint64)bProps_->sb_state_),
1148                                                             fish::Parameter(ss),
1149                                                             fish::Parameter(se)
1150                                                           };
1151                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1152                        fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1153                    }
1154                    // Resolve sliding.
1155                    double crit = sb_F_.x() * fric_;
1156                    if (crit < 0)
1157                        crit = 0;
1158                    DVect sforce = sb_F_; sforce.rx() = 0.0;
1159                    // The is the magnitude of the shear force.
1160                    double sfmag = sforce.mag();
1161                    // Sliding occurs when the magnitude of the shear force is greater than the
1162                    // critical value.
1163                    if (sfmag > crit) {
1164                        // Lower the shear force to the critical value for sliding.
1165                        double rat = crit / sfmag;
1166                        sforce *= rat;
1167                        sforce.rx() = sb_F_.x();
1168                        sb_F_ = sforce;
1169                        sb_S_ = true;
1170                    }
1171
1172                    // Resolve bending
1173                    crit = sb_bmul_*2.1*0.25*stiff.reff_*std::abs(sb_F_.x()); // Jiang 2015
1174                    DAVect test = sb_M_;
1175#if DIM==3
1176                    test.rx() = 0.0;
1177#endif
1178                    double tmag = test.mag();
1179                    if (tmag > crit) {
1180                        // Lower the bending moment to the critical value for sliding.
1181                        double rat = crit / tmag;
1182                        test *= rat;
1183                        sb_BS_ = true;
1184                    }
1185                    sb_M_.rz() = test.z();
1186#if DIM==3
1187                    sb_M_.ry() = test.y();
1188                    // Resolve twisting
1189                    crit = sb_tmul_ * 0.65*fric_* stiff.reff_*std::abs(sb_F_.x()) ; // Jiang 2015
1190                    tmag = std::abs(sb_M_.x());
1191                    if (tmag > crit) {
1192                        // Lower the shear force to the critical value for sliding.
1193                        double rat = crit / tmag;
1194                        tmag = sb_M_.x() * rat;
1195                        sb_TS_ = true;
1196                    } else
1197                        tmag = sb_M_.x();
1198                    sb_M_.rx() = tmag;
1199#endif
1200                }
1201            }
1202        }
1203
1204        // Account for dashpot forces if the dashpot structure has been defined.
1205        if (dpProps_) {
1206            dpProps_->dp_F_.fill(0.0);
1207            double vcn(0.0), vcs(0.0);
1208            // Calculate the damping coefficients.
1209            vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1210            vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1211            // First damp the shear components
1212            for (int i = 1; i<dim; ++i)
1213                dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1214            // Damp the normal component
1215            dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1216            // Need to change behavior based on the dp_mode.
1217            if (bProps_->sb_state_ < 3 && (dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1218                // Limit in tension if not bonded.
1219                if (dpProps_->dp_F_.x() + sb_F_.x() < 0)
1220                    dpProps_->dp_F_.rx() = -sb_F_.rx();
1221            }
1222            if (sb_S_ && dpProps_->dp_mode_ > 1) {
1223                // Limit in shear if sliding.
1224                double dfn = dpProps_->dp_F_.rx();
1225                dpProps_->dp_F_.fill(0.0);
1226                dpProps_->dp_F_.rx() = dfn;
1227            }
1228        }
1229
1230        //Compute energies if energy tracking has been enabled.
1231        if (state->trackEnergy_) {
1232            assert(energies_);
1233            energies_->estrain_ = 0.0;
1234            if (kn)
1235                // Calculate the strain energy.
1236                energies_->estrain_ = 0.5*sb_F_.x()*sb_F_.x() / kn;
1237            if (ks) {
1238                DVect s = sb_F_;
1239                s.rx() = 0.0;
1240                double smag2 = s.mag2();
1241                // Add the shear component of the strain energy.
1242                energies_->estrain_ += 0.5*smag2 / ks;
1243
1244                if (sb_S_) {
1245                    // If sliding calculate the slip energy and accumulate it.
1246                    sb_F_old.rx() = 0.0;
1247                    DVect avg_F_s = (s + sb_F_old)*0.5;
1248                    DVect u_s_el = (s - sb_F_old) / ks;
1249                    DVect u_s(0.0);
1250                    for (int i = 1; i < dim; ++i)
1251                        u_s.rdof(i) = trans.dof(i);
1252                    energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1253                }
1254            }
1255            // Add the bending/twisting resistance energy contributions.
1256            if (kb) {
1257                DAVect tmp = sb_M_;
1258#ifdef THREED
1259                tmp.rx() = 0.0;
1260#endif
1261                energies_->estrain_ += 0.5*tmp.mag2() / kb;
1262                if (sb_BS_) {
1263                    //  accumulate bending slip energy.
1264                    DAVect tmp_old = sb_M_old;
1265#ifdef THREED
1266                    tmp_old.rx() = 0.0;
1267#endif
1268                    DAVect avg_M = (tmp + tmp_old)*0.5;
1269                    DAVect t_s_el = (tmp - tmp_old) / kb;
1270                    energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1271                }
1272            }
1273#ifdef THREED
1274            if (kt) {
1275                double mt = std::abs(sb_M_.x());
1276                energies_->estrain_ += 0.5*mt*mt / kt;
1277                if (sb_TS_) {
1278                    //  accumulate twisting slip energy.
1279                    DAVect tmp(0.0);
1280                    DAVect tmp_old(0.0);
1281                    tmp.rx() = sb_M_.x();
1282                    tmp_old.rx() = sb_M_old.x();
1283                    DAVect avg_M = (tmp + tmp_old)*0.5;
1284                    DAVect t_s_el = (tmp - tmp_old) / kt;
1285                    energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1286                }
1287            }
1288#endif
1289
1290            if (dpProps_) {
1291                // Calculate damping energy (accumulated) if the dashpots are active.
1292                energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1293            }
1294        }
1295
1296        // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1297        assert(sb_F_ == sb_F_);
1298        assert(sb_M_ == sb_M_);
1299        return true;
1300    }
1301
1302    bool ContactModelSoftBond::FDLawUnBonded(ContactModelMechanicalState *state, const double &timestep) {
1303
1304        // Relative translational/rotational displacement increments
1305        DVect  trans = state->relativeTranslationalIncrement_;
1306        DAVect ang   = state->relativeAngularIncrement_;
1307        double overlap = rgap_ - state->gap_;
1308        double correction = 1.0;
1309        if (state->activated() && sb_mode_ == 0 && trans.x()) {
1310                correction = -1.0*overlap / trans.x();
1311                if (correction < 0)
1312                    correction = 1.0;
1313        }
1314
1315        // Store previous force and moment
1316        DVect  sb_F_old = sb_F_;
1317        DAVect sb_M_old = sb_M_;
1318
1319        // Update stiffness data
1320        StiffData stiff = computeStiffData(state);
1321        double kn = stiff.trans_.x();
1322        double ks = stiff.trans_.y();
1323        double kb = stiff.ang_.z();
1324#if DIM==3
1325        double kt = stiff.ang_.x();
1326#endif
1327        // absolute/incremental normal force calculation
1328        if (sb_mode_==0)
1329            sb_F_.rx() = overlap * kn;
1330        else
1331            sb_F_.rx() -= trans.x() * kn;
1332        // Normal force can only be positive if unbonded
1333        sb_F_.rx() = std::max(0.0, sb_F_.x());
1334
1335        // Calculate the trial shear force.
1336        DVect sforce(0.0);
1337        // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
1338        // Loop over the shear components (note: the 0 component is the normal component)
1339        // and calculate the shear force.
1340        for (int i = 1; i<dim; ++i)
1341            sforce.rdof(i) = sb_F_.dof(i) - trans.dof(i) * ks;
1342
1343        // Calculate the trial moment.
1344        DAVect mom = sb_M_ - ang*stiff.ang_;
1345
1346        // If the SOLVE ELASTIC command is given then the
1347        // canFail state is set to FALSE. Otherwise it is always TRUE.
1348        if (state->canFail_) {
1349            bool changed = false;
1350            // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
1351            bool slip_changed = false;
1352            double crit = sb_F_.x() * fric_;
1353            // The is the magnitude of the shear force.
1354            double sfmag = sforce.mag();
1355            // Sliding occurs when the magnitude of the shear force is greater than the
1356            // critical value.
1357            if (sfmag > crit) {
1358                // Lower the shear force to the critical value for sliding.
1359                double rat = crit / sfmag;
1360                sforce *= rat;
1361                if (!sb_S_) {
1362                    slip_changed = true;
1363                    changed = true;
1364                }
1365                sb_S_ = true;
1366            }
1367            else {
1368                if (sb_S_) {
1369                    slip_changed = true;
1370                    changed = true;
1371                }
1372                sb_S_ = false;
1373            }
1374
1375            // Resolve bending
1376            bool bslip_changed = false;
1377            crit = sb_bmul_ * 2.1*0.25*sb_F_.x() * stiff.reff_; // Jiang 2015
1378            DAVect test = mom;
1379#if DIM==3
1380            test.rx() = 0.0;
1381#endif
1382            double tmag = test.mag();
1383            if (tmag > crit) {
1384                // Lower the bending moment to the critical value for sliding.
1385                double rat = crit / tmag;
1386                test *= rat;
1387                if (!sb_BS_) {
1388                    bslip_changed = true;
1389                    changed = true;
1390                }
1391                sb_BS_ = true;
1392            }
1393            else {
1394                if (sb_BS_) {
1395                    bslip_changed = true;
1396                    changed = true;
1397                }
1398                sb_BS_ = false;
1399            }
1400            mom.rz() = test.z();
1401#if DIM==3
1402            mom.ry() = test.y();
1403            // Resolve twisting
1404            bool tslip_changed = false;
1405            crit = sb_tmul_ * 0.65*fric_*sb_F_.x() * stiff.reff_; // Jiang 2015
1406            tmag = std::abs(mom.x());
1407            if (tmag > crit) {
1408                // Lower the twisting moment to the critical value for sliding.
1409                double rat = crit / tmag;
1410                mom.rx() *= rat;
1411                if (!sb_TS_) {
1412                    tslip_changed = true;
1413                    changed = true;
1414                }
1415                sb_TS_ = true;
1416            } else {
1417                if (sb_TS_) {
1418                    tslip_changed = true;
1419                    changed = true;
1420                }
1421                sb_TS_ = false;
1422            }
1423#endif
1424            if (changed && cmEvents_[fSlipChange] >= 0) {
1425                qint64 code = 0;
1426                if (slip_changed) {
1427                    code = 1;
1428                    if (bslip_changed) {
1429                        code = 4;
1430#if DIM==3
1431                        if (tslip_changed)
1432                            code = 7;
1433#endif
1434                    }
1435                }
1436                else if (bslip_changed) {
1437                    code = 2;
1438#if DIM==3
1439                    if (tslip_changed)
1440                        code = 6;
1441#endif
1442                }
1443#if DIM==3
1444                else if (tslip_changed) {
1445                    code = 3;
1446                    if (slip_changed)
1447                        code = 5;
1448                }
1449#endif
1450                auto c = state->getContact();
1451                std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1452                                                     fish::Parameter(code),
1453                                                     fish::Parameter(sb_S_),
1454                                                     fish::Parameter(sb_BS_)
1455#ifdef THREED
1456                                                     ,fish::Parameter(sb_TS_)
1457#endif
1458                                                   };
1459                IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1460                fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1461            }
1462        }
1463
1464        // Set the shear components of the total force.
1465        for (int i = 1; i<dim; ++i)
1466            sb_F_.rdof(i) = sforce.dof(i);
1467
1468        // Set the moment.
1469        sb_M_ = mom;
1470
1471        // Account for dashpot forces if the dashpot structure has been defined.
1472        if (dpProps_) {
1473            dpProps_->dp_F_.fill(0.0);
1474            double vcn(0.0), vcs(0.0);
1475            // Calculate the damping coefficients.
1476            vcn = dpProps_->dp_nratio_ * 2.0 * sqrt((state->inertialMass_*(kn)));
1477            vcs = dpProps_->dp_sratio_ * 2.0 * sqrt((state->inertialMass_*(ks)));
1478            // First damp the shear components
1479            for (int i = 1; i<dim; ++i)
1480                dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
1481            // Damp the normal component
1482            dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
1483            // Need to change behavior based on the dp_mode.
1484            if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3)) {
1485                // Limit in tension if not bonded.
1486                if (dpProps_->dp_F_.x() + sb_F_.x() < 0)
1487                    dpProps_->dp_F_.rx() = -sb_F_.rx();
1488            }
1489            if (sb_S_ && dpProps_->dp_mode_ > 1) {
1490                // Limit in shear if not sliding.
1491                double dfn = dpProps_->dp_F_.rx();
1492                dpProps_->dp_F_.fill(0.0);
1493                dpProps_->dp_F_.rx() = dfn;
1494            }
1495        }
1496
1497        //Compute energies if energy tracking has been enabled.
1498        if (state->trackEnergy_) {
1499            assert(energies_);
1500            energies_->estrain_ = 0.0;
1501            if (kn_)
1502                // Calculate the strain energy.
1503                energies_->estrain_ = 0.5*sb_F_.x()*sb_F_.x() / kn;
1504            if (ks_) {
1505                DVect s = sb_F_;
1506                s.rx() = 0.0;
1507                double smag2 = s.mag2();
1508                // Add the shear component of the strain energy.
1509                energies_->estrain_ += 0.5*smag2 / ks;
1510
1511                if (sb_S_) {
1512                    // If sliding calculate the slip energy and accumulate it.
1513                    sb_F_old.rx() = 0.0;
1514                    DVect avg_F_s = (s + sb_F_old)*0.5;
1515                    DVect u_s_el = (s - sb_F_old) / ks;
1516                    DVect u_s(0.0);
1517                    for (int i = 1; i < dim; ++i)
1518                        u_s.rdof(i) = trans.dof(i);
1519                    energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s + u_s_el)));
1520                }
1521            }
1522            // Add the bending/twisting resistance energy contributions.
1523            if (kb) {
1524                DAVect tmp = sb_M_;
1525#ifdef THREED
1526                tmp.rx() = 0.0;
1527#endif
1528                energies_->estrain_ += 0.5*tmp.mag2() / kb;
1529                if (sb_BS_) {
1530                    //  accumulate bending slip energy.
1531                    DAVect tmp_old = sb_M_old;
1532#ifdef THREED
1533                    tmp_old.rx() = 0.0;
1534#endif
1535                    DAVect avg_M = (tmp + tmp_old)*0.5;
1536                    DAVect t_s_el = (tmp - tmp_old) / kb;
1537                    energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1538                }
1539            }
1540#ifdef THREED
1541            if (kt) {
1542                double mt = std::abs(sb_M_.x());
1543                energies_->estrain_ += 0.5*mt*mt / kt;
1544                if (sb_TS_) {
1545                    //  accumulate twisting slip energy.
1546                    DAVect tmp(0.0);
1547                    DAVect tmp_old(0.0);
1548                    tmp.rx() = sb_M_.x();
1549                    tmp_old.rx() = sb_M_old.x();
1550                    DAVect avg_M = (tmp + tmp_old)*0.5;
1551                    DAVect t_s_el = (tmp - tmp_old) / kt;
1552                    energies_->eslip_ -= std::min(0.0, (avg_M | (ang + t_s_el)));
1553                }
1554            }
1555#endif
1556
1557            if (dpProps_) {
1558                // Calculate damping energy (accumulated) if the dashpots are active.
1559                energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1560            }
1561        }
1562
1563        // This is just a sanity check to ensure, in debug mode, that the force/moment aren't wonky.
1564        assert(sb_F_ == sb_F_);
1565        assert(sb_M_ == sb_M_);
1566        return true;
1567    }
1568
1569    bool ContactModelSoftBond::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
1570        // Account for thermal expansion in incremental mode
1571        if (sb_mode_ == 0 || ts->gapInc_ == 0.0) return false;
1572        DVect finc(0.0);
1573        finc.rx() = kn_ * ts->gapInc_;
1574        sb_F_ -= finc;
1575        return true;
1576    }
1577
1578    void ContactModelSoftBond::setForce(const DVect &v,IContact *c) {
1579        sb_F(v);
1580        if (v.x() > 0)
1581            rgap_ = c->getGap() + v.x() / (kn_ * computeGeomData(convert_getcast<IContactMechanical>(c)).x());
1582    }
1583
1584    void ContactModelSoftBond::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
1585        // Only called for contacts with wall facets when the wall resolution scheme
1586        // is set to full!
1587        // Only do something if the contact model is of the same type
1588        if (old->getContactModel()->getName().compare("softbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
1589            ContactModelSoftBond *oldCm = (ContactModelSoftBond *)old;
1590#ifdef THREED
1591            // Need to rotate just the shear component from oldSystem to newSystem
1592
1593            // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
1594            DVect axis = oldSystem.e1() & newSystem.e1();
1595            double c, ang, s;
1596            DVect re2;
1597            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1598                axis = axis.unit();
1599                c = oldSystem.e1()|newSystem.e1();
1600                if (c > 0)
1601                    c = std::min(c,1.0);
1602                else
1603                    c = std::max(c,-1.0);
1604                ang = acos(c);
1605                s = sin(ang);
1606                double t = 1. - c;
1607                DMatrix<3,3> rm;
1608                rm.get(0,0) = t*axis.x()*axis.x() + c;
1609                rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
1610                rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
1611                rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
1612                rm.get(1,1) = t*axis.y()*axis.y() + c;
1613                rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
1614                rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
1615                rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
1616                rm.get(2,2) = t*axis.z()*axis.z() + c;
1617                re2 = rm*oldSystem.e2();
1618            }
1619            else
1620                re2 = oldSystem.e2();
1621            // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
1622            axis = re2 & newSystem.e2();
1623            DVect2 tpf;
1624            DVect2 tpm;
1625            DMatrix<2,2> m;
1626            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1627                axis = axis.unit();
1628                c = re2|newSystem.e2();
1629                if (c > 0)
1630                    c = std::min(c,1.0);
1631                else
1632                    c = std::max(c,-1.0);
1633                ang = acos(c);
1634                if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
1635                    ang *= -1;
1636                s = sin(ang);
1637                m.get(0,0) = c;
1638                m.get(1,0) = s;
1639                m.get(0,1) = -m.get(1,0);
1640                m.get(1,1) = m.get(0,0);
1641                tpf = m*DVect2(oldCm->sb_F_.y(),oldCm->sb_F_.z());
1642                tpm = m*DVect2(oldCm->sb_M_.y(),oldCm->sb_M_.z());
1643            } else {
1644                m.get(0,0) = 1.;
1645                m.get(0,1) = 0.;
1646                m.get(1,0) = 0.;
1647                m.get(1,1) = 1.;
1648                tpf = DVect2(oldCm->sb_F_.y(),oldCm->sb_F_.z());
1649                tpm = DVect2(oldCm->sb_M_.y(),oldCm->sb_M_.z());
1650            }
1651            DVect pforce = DVect(0,tpf.x(),tpf.y());
1652            //DVect pm     = DVect(0,tpm.x(),tpm.y());
1653#else
1654            oldSystem;
1655            newSystem;
1656            DVect pforce = DVect(0,oldCm->sb_F_.y());
1657            //DVect pm     = DVect(0,oldCm->sb_M_.y());
1658#endif
1659            for (int i=1; i<dim; ++i)
1660                sb_F_.rdof(i) += pforce.dof(i);
1661            if (sb_mode_ && oldCm->sb_mode_)
1662                sb_F_.rx() = oldCm->sb_F_.x();
1663            oldCm->sb_F_ = DVect(0.0);
1664            oldCm->sb_M_ = DAVect(0.0);
1665            if (dpProps_ && oldCm->dpProps_) {
1666#ifdef THREED
1667                tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
1668                pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
1669#else
1670                pforce = oldCm->dpProps_->dp_F_;
1671#endif
1672                dpProps_->dp_F_ += pforce;
1673                oldCm->dpProps_->dp_F_ = DVect(0.0);
1674            }
1675            if(oldCm->getEnergyActivated()) {
1676                activateEnergy();
1677                energies_->estrain_ = oldCm->energies_->estrain_;
1678                energies_->edashpot_ = oldCm->energies_->edashpot_;
1679                energies_->eslip_ = oldCm->energies_->eslip_;
1680                oldCm->energies_->estrain_ = 0.0;
1681                oldCm->energies_->edashpot_ = 0.0;
1682                oldCm->energies_->eslip_ = 0.0;
1683            }
1684            rgap_ = oldCm->rgap_;
1685        }
1686        assert(sb_F_ == sb_F_);
1687    }
1688
1689    void ContactModelSoftBond::setNonForcePropsFrom(IContactModel *old) {
1690        // Only called for contacts with wall facets when the wall resolution scheme
1691        // is set to full!
1692        // Only do something if the contact model is of the same type
1693        if (old->getName().compare("softbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
1694            ContactModelSoftBond *oldCm = (ContactModelSoftBond *)old;
1695            kn_       = oldCm->kn_;
1696            ks_       = oldCm->ks_;
1697            fric_     = oldCm->fric_;
1698            sb_bmul_  = oldCm->sb_bmul_;
1699            sb_tmul_  = oldCm->sb_tmul_;
1700            sb_mode_  = oldCm->sb_mode_;
1701            sb_rmul_  = oldCm->sb_rmul_;
1702            sb_S_     = oldCm->sb_S_;
1703            sb_BS_    = oldCm->sb_BS_;
1704            sb_TS_    = oldCm->sb_TS_;
1705            rgap_     = oldCm->rgap_;
1706            userArea_ = oldCm->userArea_;
1707
1708            if (oldCm->dpProps_) {
1709                if (!dpProps_)
1710                    dpProps_ = NEW dpProps();
1711                dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1712                dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1713                dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1714            }
1715            if (oldCm->bProps_) {
1716                if (!bProps_)
1717                    bProps_ = NEW bProps();
1718                bProps_->sb_mcf_    = oldCm->bProps_->sb_mcf_;
1719                bProps_->sb_fa_     = oldCm->bProps_->sb_fa_;
1720                bProps_->sb_state_  = oldCm->bProps_->sb_state_;
1721                bProps_->sb_coh_    = oldCm->bProps_->sb_coh_;
1722                bProps_->sb_ten_    = oldCm->bProps_->sb_ten_;
1723                bProps_->sb_maxTen_ = oldCm->bProps_->sb_maxTen_;
1724                bProps_->sb_cut_    = oldCm->bProps_->sb_cut_;
1725                bProps_->sb_delu_   = oldCm->bProps_->sb_delu_;
1726                bProps_->sb_delo_   = oldCm->bProps_->sb_delo_;
1727                bProps_->sb_maxu_   = oldCm->bProps_->sb_maxu_;
1728                bProps_->sb_critu_  = oldCm->bProps_->sb_critu_;
1729            }
1730
1731        }
1732    }
1733
1734    DVect ContactModelSoftBond::getForce(const IContactMechanical *) const {
1735        DVect ret(sb_F_);
1736        if (dpProps_)
1737            ret += dpProps_->dp_F_;
1738        return ret;
1739    }
1740
1741    DAVect ContactModelSoftBond::getMomentOn1(const IContactMechanical *c) const {
1742        DVect force = getForce(c);
1743        DAVect ret(sb_M_);
1744        c->updateResultingTorqueOn1Local(force,&ret);
1745        return ret;
1746    }
1747
1748    DAVect ContactModelSoftBond::getMomentOn2(const IContactMechanical *c) const {
1749        DVect force = getForce(c);
1750        DAVect ret(sb_M_);
1751        c->updateResultingTorqueOn2Local(force,&ret);
1752        return ret;
1753    }
1754
1755    DVect3 ContactModelSoftBond::computeGeomData(const IContactMechanical *c) const {
1756        double Cmax1 = c->getEnd1Curvature().y();
1757        double Cmax2 = c->getEnd2Curvature().y();
1758        double br = sb_rmul_ * 1.0 / std::max(Cmax1, Cmax2);
1759        if (userArea_)
1760#ifdef THREED
1761            br = std::sqrt(userArea_ / dPi);
1762#else
1763            br = userArea_ / 2.0;
1764#endif
1765        double br2 = br * br;
1766#ifdef TWOD
1767        double area = 2.0*br;
1768        double bi = 2.0*br*br2 / 3.0;
1769#else
1770        double area = dPi * br2;
1771        double bi = 0.25*area*br2;
1772#endif
1773        return DVect3(area, bi, br);
1774    }
1775
1776    DVect2 ContactModelSoftBond::SMax(const IContactMechanical *c) const {
1777        DVect3 data = computeGeomData(c);
1778        double area = data.x();
1779        double bi = data.y();
1780        double br = data.z();
1781        /* maximum stresses */
1782        double dbend = sqrt(sb_M_.y()*sb_M_.y() + sb_M_.z()*sb_M_.z());
1783        double dtwist = sb_M_.x();
1784        DVect bfs(sb_F_);
1785        bfs.rx() = 0.0;
1786        double dbfs = bfs.mag();
1787        double nsmax = -(sb_F_.x() / area) + dbend * br / bi;
1788        double ssmax = dbfs / area + std::abs(dtwist) * 0.5* br / bi;
1789        return DVect2(nsmax, ssmax);
1790    }
1791
1792    double ContactModelSoftBond::shearStrength(const double &area) const {
1793        if (!bProps_) return 0.0;
1794        double sig = -1.0*sb_F_.x() / area;
1795        double nstr = bProps_->sb_state_ > 2 ? bProps_->sb_ten_ : 0.0;
1796        return sig <= nstr ? bProps_->sb_coh_ - std::tan(dDegrad*bProps_->sb_fa_)*sig
1797            : bProps_->sb_coh_ - std::tan(dDegrad*bProps_->sb_fa_)*nstr;
1798    }
1799
1800
1801    double ContactModelSoftBond::strainEnergy(double kn,double ks,double kb,double kt) const {
1802        double ret(0.0);
1803        if (kn)
1804            ret = 0.5 * sb_F_.x() * sb_F_.x() / kn;
1805        if (ks) {
1806            DVect tmp = sb_F_;
1807            tmp.rx() = 0.0;
1808            double smag2 = tmp.mag2();
1809            ret += 0.5 * smag2 / ks;
1810        }
1811
1812        if (kt)
1813            ret += 0.5 * sb_M_.x() * sb_M_.x() / kt;
1814        if (kb) {
1815            DAVect tmp = sb_M_;
1816#ifdef THREED
1817            tmp.rx() = 0.0;
1818            double smag2 = tmp.mag2();
1819#else
1820            double smag2 = tmp.z() * tmp.z();
1821#endif
1822            ret += 0.5 * smag2 / kb;
1823        }
1824        return ret;
1825    }
1826
1827} // namespace cmodelsxd
1828// EoF

Top