Linear Parallel Bond Model Implementation

See this page for the documentation of this contact model.

contactmodellinearpbond.h

  1#pragma once
  2// contactmodellinearpbond.h
  3
  4#include "contactmodel/src/contactmodelmechanical.h"
  5
  6#ifdef LINEARPBOND_LIB
  7#  define LINEARPBOND_EXPORT EXPORT_TAG
  8#elif defined(NO_MODEL_IMPORT)
  9#  define LINEARPBOND_EXPORT
 10#else
 11#  define LINEARPBOND_EXPORT IMPORT_TAG
 12#endif
 13
 14namespace cmodelsxd {
 15    using namespace itasca;
 16 
 17    class ContactModelLinearPBond : public ContactModelMechanical {
 18    public:
 19        LINEARPBOND_EXPORT ContactModelLinearPBond();
 20        LINEARPBOND_EXPORT virtual ~ContactModelLinearPBond();
 21        void     copy(const ContactModel *c) override;
 22        void     archive(ArchiveStream &) override;
 23        QString  getName() const override { return "linearpbond"; }
 24        void     setIndex(int i) override { index_=i;}
 25        int      getIndex() const override {return index_;}
 26
 27        enum PropertyKeys { 
 28              kwLinKn=1
 29            , kwLinKs                            
 30            , kwLinFric   
 31            , kwLinF
 32            , kwLinS 
 33            , kwLinMode
 34            , kwRGap
 35            , kwEmod
 36            , kwKRatio
 37            , kwDpNRatio 
 38            , kwDpSRatio
 39            , kwDpMode 
 40            , kwDpF
 41            , kwPbState
 42            , kwPbRMul                        
 43            , kwPbKn  
 44            , kwPbKs
 45            , kwPbMcf
 46            , kwPbTStrength 
 47            , kwPbSStrength 
 48            , kwPbCoh
 49            , kwPbFa 
 50            , kwPbSig
 51            , kwPbTau 
 52            , kwPbF
 53            , kwPbM 
 54            , kwPbRadius 
 55            , kwPbEmod
 56            , kwPbKRatio
 57            , kwUserArea
 58        };
 59         
 60        QString  getProperties() const override {
 61            return "kn"
 62                   ",ks"
 63                   ",fric"
 64                   ",lin_force"
 65                   ",lin_slip"
 66                   ",lin_mode"
 67                   ",rgap"
 68                   ",emod"
 69                   ",kratio"
 70                   ",dp_nratio"
 71                   ",dp_sratio"
 72                   ",dp_mode"
 73                   ",dp_force"
 74                   ",pb_state"
 75                   ",pb_rmul"
 76                   ",pb_kn"
 77                   ",pb_ks"
 78                   ",pb_mcf"
 79                   ",pb_ten"
 80                   ",pb_shear"
 81                   ",pb_coh"
 82                   ",pb_fa"
 83                   ",pb_sigma"
 84                   ",pb_tau"
 85                   ",pb_force"
 86                   ",pb_moment"
 87                   ",pb_radius"
 88                   ",pb_emod"
 89                   ",pb_kratio"
 90                   ",user_area";
 91        }
 92
 93        enum EnergyKeys { kwEStrain=1,kwESlip,kwEDashpot,kwEPbStrain};
 94        QString  getEnergies() const override { return "energy-strain,energy-slip,energy-dashpot,energy-pbstrain";}
 95        double   getEnergy(uint32 i) const override;  // Base 1
 96        bool     getEnergyAccumulate(uint32 i) const override; // Base 1
 97        void     setEnergy(uint32 i,const double &d) override; // Base 1
 98        void     activateEnergy() override { if (energies_) return; energies_ = NEW Energies();}
 99        bool     getEnergyActivated() const override {return (energies_ !=0);}
100
101        enum FishCallEvents {fActivated=0,fBondBreak,fSlipChange};
102        QString  getFishCallEvents() const override { return "contact_activated,bond_break,slip_change"; }
103        QVariant getProperty(uint32 i,const IContact *con=0) const override;
104        bool     getPropertyGlobal(uint32 i) const override;
105        bool     setProperty(uint32 i,const QVariant &v,IContact *con=0);
106        bool     getPropertyReadOnly(uint32 i) const override;
107        bool     supportsInheritance(uint32 i) const override;
108        bool     getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
109        void     setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }
110
111        enum MethodKeys { kwDeformability=1
112                        , kwPbDeformability
113                        , kwPbBond 
114                        , kwPbUnbond 
115                        , kwArea
116        };
117
118        QString  getMethods() const override {
119            return "deformability"
120                   ",pb_deformability"
121                   ",bond"
122                   ",unbond"
123                   ",area";
124        }
125
126        QString  getMethodArguments(uint32 i) const override;
127
128        bool     setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override; // Base 1 - returns true if timestep contributions need to be updated
129
130        uint32     getMinorVersion() const override;
131
132        bool    validate(ContactModelMechanicalState *state,const double &timestep) override;
133        bool    endPropertyUpdated(const QString &name,const IContactMechanical *c) override;
134        bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) override;
135        DVect2  getEffectiveTranslationalStiffness() const override { DVect2 ret = effectiveTranslationalStiffness_; if(pbProps_) ret+= pbProps_->pbTransStiff_ ;return ret;}
136        DAVect  getEffectiveRotationalStiffness() const override {if (!pbProps_) return DAVect(0.0); return pbProps_->pbAngStiff_;}
137
138        bool thermalCoupling(ContactModelMechanicalState *, ContactModelThermalState * , IContactThermal *,const double &) override;
139
140        ContactModelLinearPBond *clone() const override { return NEW ContactModelLinearPBond(); }
141        double   getActivityDistance() const override {return rgap_;}
142        bool     isOKToDelete() const override { return !isBonded(); }
143        void     resetForcesAndMoments() override { lin_F(DVect(0.0)); dp_F(DVect(0.0)); pbF(DVect(0.0)); pbM(DAVect(0.0));  if (energies_) { energies_->estrain_ = 0.0;  if (energies_) energies_->epbstrain_ = 0.0;}}
144        void     setForce(const DVect &v,IContact *c) override;
145        void     setArea(const double &d) override { userArea_ = d; }
146        double   getArea() const override { return userArea_; }
147        bool     checkActivity(const double &gap) override { return (gap <= rgap_ || isBonded()); }
148        bool     isSliding() const override { return lin_S_; }
149        bool     isBonded() const override { return pbProps_ ? (pbProps_->pb_state_==3) : false; }
150        void     unbond() override { if (pbProps_) pbProps_->pb_state_= 0; }
151        void     propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
152        void     setNonForcePropsFrom(IContactModel *oldCM) override;
153        /// Return the total force that the contact model holds.
154        DVect    getForce(const IContactMechanical*) const override;
155        /// Return the total moment on 1 that the contact model holds
156        DAVect   getMomentOn1(const IContactMechanical*) const override;
157        /// Return the total moment on 1 that the contact model holds
158        DAVect   getMomentOn2(const IContactMechanical*) const override;
159
160        const double & kn() const {return kn_;}
161        void           kn(const double &d) {kn_=d;}
162        const double & ks() const {return ks_;}
163        void           ks(const double &d) {ks_=d;}
164        const double & fric() const {return fric_;}
165        void           fric(const double &d) {fric_=d;}
166        const DVect &  lin_F() const {return lin_F_;}
167        void           lin_F(const DVect &f) { lin_F_=f;}
168        bool           lin_S() const {return lin_S_;}
169        void           lin_S(bool b) { lin_S_=b;}
170        uint32           lin_mode() const {return lin_mode_;}
171        void           lin_mode(uint32 i) { lin_mode_=i;}
172        const double & rgap() const {return rgap_;}
173        void           rgap(const double &d) {rgap_=d;}
174
175        bool     hasDamping() const {return dpProps_ ? true : false;}
176        double   dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
177        void     dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
178        double   dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
179        void     dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
180        int      dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
181        void     dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
182        DVect    dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
183        void     dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
184
185        bool    hasEnergies() const {return energies_ ? true:false;}
186        double  estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
187        void    estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
188        double  eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
189        void    eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
190        double  edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
191        void    edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
192        double  epbstrain() const {return hasEnergies() ? energies_->epbstrain_: 0.0;}
193        void    epbstrain(const double &d) { if(!hasEnergies()) return; energies_->epbstrain_=d;}
194
195        bool     hasPBond() const {return pbProps_ ? true:false;}
196        int      pbState()   const {return hasPBond() ? pbProps_->pb_state_: 0;}   
197        void     pbState(int i) { if(!hasPBond()) return; pbProps_->pb_state_=i;}
198        double   pbRmul() const {return (hasPBond() ? (pbProps_->pb_rmul_) : 0.0);}
199        void     pbRmul(const double &d) { if(!hasPBond()) return; pbProps_->pb_rmul_=d;}
200        double   pbKn() const {return (hasPBond() ? (pbProps_->pb_kn_) : 0.0);}
201        void     pbKn(const double &d) { if(!hasPBond()) return; pbProps_->pb_kn_=d;}
202        double   pbKs() const {return (hasPBond() ? (pbProps_->pb_ks_) : 0.0);}
203        void     pbKs(const double &d) { if(!hasPBond()) return; pbProps_->pb_ks_=d;}
204        double   pbMCF() const {return (hasPBond() ? (pbProps_->pb_mcf_) : 0.0);}
205        void     pbMCF(const double &d) { if(!hasPBond()) return; pbProps_->pb_mcf_=d;}
206        double   pbTen() const {return (hasPBond() ? (pbProps_->pb_ten_) : 0.0);}
207        void     pbTen(const double &d) { if(!hasPBond()) return; pbProps_->pb_ten_=d;}
208        double   pbCoh() const {return (hasPBond() ? (pbProps_->pb_coh_) : 0.0);}
209        void     pbCoh(const double &d) { if(!hasPBond()) return; pbProps_->pb_coh_=d;}
210        double   pbFA() const {return (hasPBond() ? (pbProps_->pb_fa_) : 0.0);}
211        void     pbFA(const double &d) { if(!hasPBond()) return; pbProps_->pb_fa_=d;}
212        DVect    pbF() const {return hasPBond() ? pbProps_->pb_F_: DVect(0.0);}
213        void     pbF(const DVect &f) { if(!hasPBond()) return; pbProps_->pb_F_=f;}
214        DAVect   pbM() const {return hasPBond() ? pbProps_->pb_M_: DAVect(0.0);}
215        void     pbM(const DAVect &m) { if(!hasPBond()) return; pbProps_->pb_M_=m;}
216        DVect2   pbTransStiff() const {return hasPBond() ? pbProps_->pbTransStiff_: DVect2(0.0);}
217        void     pbTransStiff(const DVect2 &f) { if(!hasPBond()) return; pbProps_->pbTransStiff_=f;}
218        DAVect   pbAngStiff() const {return hasPBond() ? pbProps_->pbAngStiff_: DAVect(0.0);}
219        void     pbAngStiff(const DAVect &m) { if(!hasPBond()) return; pbProps_->pbAngStiff_=m;}
220
221        uint32 inheritanceField() const {return inheritanceField_;}
222        void inheritanceField(uint32 i) {inheritanceField_ = i;}
223
224        const DVect2 & effectiveTranslationalStiffness()  const          {return effectiveTranslationalStiffness_;}
225        void           effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
226
227    private:
228        static int index_;
229
230        struct Energies {
231            Energies() : estrain_(0.0), eslip_(0.0), edashpot_(0.0), epbstrain_(0.0) {}
232            double estrain_;  // elastic energy stored in contact 
233            double eslip_;    // work dissipated by friction 
234            double edashpot_;    // work dissipated by dashpots
235            double epbstrain_; // parallel bond strain energy
236        };
237
238        struct dpProps {
239            dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)) {}
240            double dp_nratio_;     // normal viscous critical damping ratio
241            double dp_sratio_;     // shear  viscous critical damping ratio
242            int    dp_mode_;      // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
243            DVect  dp_F_;  // Force in the dashpots
244        };
245
246        struct pbProps {
247            pbProps() : pb_state_(0), pb_rmul_(1.0), pb_kn_(0.0), pb_ks_(0.0), 
248                        pb_mcf_(1.0), pb_ten_(0.0), pb_coh_(0.0), pb_fa_(0.0), pb_F_(DVect(0.0)), pb_M_(DAVect(0.0)),
249                        pbTransStiff_(0.0), pbAngStiff_(0.0){}
250            // parallel bond
251            int     pb_state_;         // Bond mode - 0 (NBNF), 1 (NBFT), 2 (NBFS), 3 (B)
252            double  pb_rmul_;         // Radius multiplier
253            double  pb_kn_;           // normal stiffness
254            double  pb_ks_;           // shear stiffness
255            double  pb_mcf_;          // Moment contribution factor 
256            double  pb_ten_;          // normal strength 
257            double  pb_coh_;          // cohesion
258            double  pb_fa_;           // friction angle
259            DVect   pb_F_;            // Force in parallel bond
260            DAVect  pb_M_;            // moment in parallel bond
261            DVect2  pbTransStiff_;    // (Normal,Shear) Translational stiffness of the parallel bond
262            DAVect  pbAngStiff_;      // (Normal,Shear) Rotational stiffness of the parallel bond
263        };
264
265        bool   updateKn(const IContactMechanical *con);
266        bool   updateKs(const IContactMechanical *con);
267        bool   updateFric(const IContactMechanical *con);
268        double pbStrainEnergy() const;                     // Compute  bond strain energy 
269
270        void   updateEffectiveStiffness(ContactModelMechanicalState *state);
271
272        DVect3 pbData(const IContactMechanical *con) const; // Bond area and inertia
273        DVect2 pbSMax(const IContactMechanical *con) const; // Maximum stress (tensile,shear) at bond periphery
274        double pbShearStrength(const double &pbArea) const;      // Bond shear strength
275        void   setDampCoefficients(const double &mass,double *vcn,double *vcs);
276
277        // inheritance fields
278        uint32 inheritanceField_;
279
280        // linear model
281        double      kn_;        // normal stiffness
282        double      ks_;        // shear stiffness
283        double      fric_;      // Coulomb friction coefficient
284        DVect       lin_F_;     // Force carried in the linear model
285        bool        lin_S_;     // the current sliding state
286        uint32        lin_mode_;  // specifies absolute (0) or incremental (1) behavior for the the linear part 
287        double      rgap_;      // reference gap for the linear part
288
289        dpProps *   dpProps_;     // The viscous properties
290        pbProps *   pbProps_;     // The parallel bond properties
291
292        double      userArea_;    // User specified area 
293
294        Energies *   energies_;   // energies
295
296        DVect2  effectiveTranslationalStiffness_;
297         
298    };
299} // namespace itascaxd
300// EoF

Top

contactmodellinearpbond.cpp

   1// contactmodellinear.cpp
   2#include "contactmodellinearpbond.h"
   3
   4#include "../version.txt"
   5#include "contactmodel/src/contactmodelthermal.h"
   6#include "fish/src/parameter.h"
   7#include "utility/src/tptr.h"
   8#include "shared/src/mathutil.h"
   9#include "base/src/basetoqt.h"
  10
  11#include "kernel/interface/iprogram.h"
  12#include "module/interface/icontact.h"
  13#include "module/interface/icontactmechanical.h"
  14#include "module/interface/icontactthermal.h"
  15#include "module/interface/ifishcalllist.h"
  16#include "module/interface/ipiece.h"
  17#include "module/interface/ipiecemechanical.h"
  18
  19#ifdef LINEARPBOND_LIB
  20#ifdef _WIN32
  21  int __stdcall DllMain(void *,unsigned, void *)
  22  {
  23    return 1;
  24  }
  25#endif
  26
  27  extern "C" EXPORT_TAG const char *getName() 
  28  {
  29#if DIM==3
  30    return "contactmodelmechanical3dlinearpbond";
  31#else
  32    return "contactmodelmechanical2dlinearpbond";
  33#endif
  34  }
  35
  36  extern "C" EXPORT_TAG unsigned getMajorVersion()
  37  {
  38    return MAJOR_VERSION;
  39  }
  40
  41  extern "C" EXPORT_TAG unsigned getMinorVersion()
  42  {
  43    return MINOR_VERSION;
  44  }
  45
  46  extern "C" EXPORT_TAG void *createInstance() 
  47  {
  48    cmodelsxd::ContactModelLinearPBond *m = NEW cmodelsxd::ContactModelLinearPBond();
  49    return (void *)m;
  50  }
  51#endif // LINEARPBOND_EXPORTS
  52
  53namespace cmodelsxd {
  54    static const uint32 linKnMask      = 0x00002; // Base 1!
  55    static const uint32 linKsMask      = 0x00004;
  56    static const uint32 linFricMask    = 0x00008;
  57
  58    using namespace itasca;
  59
  60    int ContactModelLinearPBond::index_ = -1;
  61    uint32 ContactModelLinearPBond::getMinorVersion() const { return MINOR_VERSION; }
  62
  63    ContactModelLinearPBond::ContactModelLinearPBond() : inheritanceField_(linKnMask|linKsMask|linFricMask) 
  64                                                       , kn_(0.0)
  65                                                       , ks_(0.0)
  66                                                       , fric_(0.0)
  67                                                       , lin_F_(DVect(0.0))
  68                                                       , lin_S_(false)
  69                                                       , lin_mode_(0)
  70                                                       , rgap_(0.0)
  71                                                       , dpProps_(0)
  72                                                       , pbProps_(0)
  73                                                       , userArea_(0)
  74                                                       , energies_(0)
  75                                                       , effectiveTranslationalStiffness_(DVect2(0.0)) {
  76//    setFromParent(ContactModelMechanicalList::instance()->find(getName()));
  77    }
  78
  79    ContactModelLinearPBond::~ContactModelLinearPBond() {
  80        if (dpProps_)
  81            delete dpProps_;
  82        if (pbProps_)
  83            delete pbProps_;
  84        if (energies_)
  85            delete energies_;
  86    }
  87
  88    void ContactModelLinearPBond::archive(ArchiveStream &stream) {
  89        stream & kn_;
  90        stream & ks_;
  91        stream & fric_;
  92        stream & lin_F_;
  93        stream & lin_S_;
  94        stream & lin_mode_;
  95        if (stream.getArchiveState()==ArchiveStream::Save) {
  96            bool b = false;
  97            if (dpProps_) {
  98                b = true;
  99                stream & b;
 100                stream & dpProps_->dp_nratio_; 
 101                stream & dpProps_->dp_sratio_; 
 102                stream & dpProps_->dp_mode_; 
 103                stream & dpProps_->dp_F_; 
 104            }
 105            else
 106                stream & b;
 107
 108            b = false;
 109            if (energies_) {
 110                b = true;
 111                stream & b;
 112                stream & energies_->estrain_;
 113                stream & energies_->eslip_;
 114                stream & energies_->edashpot_;
 115                stream & energies_->epbstrain_;
 116            }
 117            else
 118                stream & b;
 119
 120            b = false;
 121            if (pbProps_) {
 122                b = true;
 123                stream & b;
 124                stream & pbProps_->pb_state_;
 125                stream & pbProps_->pb_rmul_;
 126                stream & pbProps_->pb_kn_;
 127                stream & pbProps_->pb_ks_;
 128                stream & pbProps_->pb_mcf_;
 129                stream & pbProps_->pb_ten_;
 130                stream & pbProps_->pb_coh_;
 131                stream & pbProps_->pb_fa_;
 132                stream & pbProps_->pb_F_;
 133                stream & pbProps_->pb_M_;
 134            }
 135            else
 136                stream & b;
 137        } else {
 138            bool b(false);
 139            stream & b;
 140            if (b) {
 141                if (!dpProps_)
 142                    dpProps_ = NEW dpProps();
 143                stream & dpProps_->dp_nratio_; 
 144                stream & dpProps_->dp_sratio_; 
 145                stream & dpProps_->dp_mode_; 
 146                stream & dpProps_->dp_F_; 
 147            }
 148            stream & b;
 149            if (b) {
 150                if (!energies_)
 151                    energies_ = NEW Energies();
 152                stream & energies_->estrain_;
 153                stream & energies_->eslip_;
 154                stream & energies_->edashpot_;
 155                stream & energies_->epbstrain_;
 156            }
 157            stream & b;
 158            if (b) {
 159                if (!pbProps_)
 160                    pbProps_ = NEW pbProps();
 161                stream & pbProps_->pb_state_;
 162                stream & pbProps_->pb_rmul_;
 163                stream & pbProps_->pb_kn_;
 164                stream & pbProps_->pb_ks_;
 165                stream & pbProps_->pb_mcf_;
 166                stream & pbProps_->pb_ten_;
 167                stream & pbProps_->pb_coh_;
 168                stream & pbProps_->pb_fa_;
 169                stream & pbProps_->pb_F_;
 170                stream & pbProps_->pb_M_;
 171            }
 172        }
 173
 174        stream & inheritanceField_;
 175        stream & effectiveTranslationalStiffness_;
 176
 177        if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() == getMinorVersion())
 178            stream & rgap_;
 179
 180        if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 1) 
 181            stream & userArea_;
 182    }
 183
 184    void ContactModelLinearPBond::copy(const ContactModel *cm) {
 185        ContactModelMechanical::copy(cm);
 186        const ContactModelLinearPBond *in = dynamic_cast<const ContactModelLinearPBond*>(cm);
 187        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
 188        kn(in->kn());
 189        ks(in->ks());
 190        fric(in->fric());
 191        lin_F(in->lin_F());
 192        lin_S(in->lin_S());
 193        lin_mode(in->lin_mode());
 194        rgap(in->rgap());
 195        if (in->hasDamping()) {
 196            if (!dpProps_)
 197                dpProps_ = NEW dpProps();
 198            dp_nratio(in->dp_nratio()); 
 199            dp_sratio(in->dp_sratio()); 
 200            dp_mode(in->dp_mode()); 
 201            dp_F(in->dp_F()); 
 202        }
 203        if (in->hasEnergies()) {
 204            if (!energies_)
 205                energies_ = NEW Energies();
 206            estrain(in->estrain());
 207            eslip(in->eslip());
 208            edashpot(in->edashpot());
 209            epbstrain(in->epbstrain());
 210        }
 211        if (in->hasPBond()) {
 212            if (!pbProps_)
 213                pbProps_ = NEW pbProps();
 214            pbState(in->pbState());
 215            pbRmul(in->pbRmul());
 216            pbKn(in->pbKn());
 217            pbKs(in->pbKs());
 218            pbMCF(in->pbMCF());
 219            pbTen(in->pbTen());
 220            pbCoh(in->pbCoh());
 221            pbFA(in->pbFA());
 222            pbF(in->pbF());
 223            pbM(in->pbM());
 224            pbTransStiff(in->pbTransStiff());
 225            pbAngStiff(in->pbAngStiff());
 226        }
 227        userArea_ = in->userArea_;
 228        inheritanceField(in->inheritanceField());
 229        effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
 230    }
 231
 232    QVariant ContactModelLinearPBond::getProperty(uint32 i,const IContact *con) const {
 233        QVariant var;
 234        switch (i) {
 235        case kwLinKn:        return kn_;
 236        case kwLinKs:        return ks_;
 237        case kwLinFric:      return fric_;
 238        case kwLinMode:      return lin_mode_;
 239        case kwLinF:         var.setValue(lin_F_); return var;
 240        case kwLinS:   return lin_S_;
 241        case kwRGap:       return rgap_;
 242        case kwEmod: {
 243                        const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 244                        if (c ==nullptr) return 0.0;
 245                        double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
 246                        double rsum(0.0);
 247                        if (c->getEnd1Curvature().y())
 248                            rsum += 1.0/c->getEnd1Curvature().y();
 249                        if (c->getEnd2Curvature().y())
 250                            rsum += 1.0/c->getEnd2Curvature().y();
 251                        if (userArea_) {
 252#ifdef THREED
 253                            rsq = std::sqrt(userArea_ / dPi);
 254#else
 255                            rsq = userArea_ / 2.0;
 256#endif        
 257                            rsum = rsq + rsq;
 258                            rsq = 1. / rsq;
 259                        }
 260#ifdef TWOD             
 261                        return (kn_ * rsum * rsq / 2.0);
 262#else                     
 263                        return (kn_ * rsum * rsq * rsq) / dPi;
 264#endif                    
 265                    }
 266        case kwKRatio:      return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
 267        case kwDpNRatio:    return dpProps_ ? dpProps_->dp_nratio_ : 0;
 268        case kwDpSRatio:    return dpProps_ ? dpProps_->dp_sratio_ : 0;
 269        case kwDpMode:      return dpProps_ ? dpProps_->dp_mode_ : 0;
 270        case kwUserArea:    return userArea_;
 271        case kwDpF: {
 272                dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
 273                return var;
 274            }
 275        case kwPbState:     return pbProps_ ? pbProps_->pb_state_ : 0;
 276        case kwPbRMul:      return pbProps_ ? pbProps_->pb_rmul_ : 1.0;
 277        case kwPbKn:        return pbProps_ ? pbProps_->pb_kn_ : 0;
 278        case kwPbKs:        return pbProps_ ? pbProps_->pb_ks_ : 0;
 279        case kwPbMcf:       return pbProps_ ? pbProps_->pb_mcf_ : 1.0;
 280        case kwPbTStrength: return pbProps_ ? pbProps_->pb_ten_ : 0.0;
 281        case kwPbSStrength: {
 282                if (!pbProps_) return 0.0;
 283                const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 284                double pbArea = pbData(c).x();
 285                return pbShearStrength(pbArea);
 286            }
 287        case kwPbCoh:       return pbProps_ ? pbProps_->pb_coh_ : 0;
 288        case kwPbFa:        return pbProps_ ? pbProps_->pb_fa_ : 0;
 289        case kwPbSig: {
 290                if (!pbProps_ || pbProps_->pb_state_ < 3) return 0.0;
 291                const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 292                return pbSMax(c).x();
 293            }
 294        case kwPbTau: {
 295                if (!pbProps_ || pbProps_->pb_state_ < 3) return 0.0;
 296                const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 297                return pbSMax(c).y();
 298            }
 299        case kwPbF: {
 300                pbProps_ ? var.setValue(pbProps_->pb_F_) : var.setValue(DVect(0.0));
 301                return var;
 302            }
 303        case kwPbM: {
 304                pbProps_ ? var.setValue(pbProps_->pb_M_) : var.setValue(DAVect(0.0));
 305                return var;
 306            }
 307        case kwPbRadius: {
 308                if (!pbProps_) return 0.0;
 309                const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 310                double Cmax1 = c->getEnd1Curvature().y();
 311                double Cmax2 = c->getEnd2Curvature().y();
 312                double br = pbProps_->pb_rmul_ * 1.0 / std::max(Cmax1,Cmax2);
 313                if (userArea_)
 314#ifdef THREED
 315                    br = std::sqrt(userArea_ / dPi);
 316#else
 317                    br = userArea_ / 2.0;
 318#endif
 319                return br;
 320            }
 321        case kwPbEmod: {
 322                if (!pbProps_) return 0.0;
 323                const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 324                double rsum(0.0);
 325                if (c->getEnd1Curvature().y())
 326                    rsum += 1.0/c->getEnd1Curvature().y();
 327                if (c->getEnd2Curvature().y())
 328                    rsum += 1.0/c->getEnd2Curvature().y();
 329                if (userArea_) {
 330#ifdef THREED
 331                    double rad = std::sqrt(userArea_ / dPi);
 332#else
 333                    double rad = userArea_ / 2.0;
 334#endif        
 335                    rsum = 2 * rad;
 336                }
 337                double emod = pbProps_->pb_kn_ * rsum;
 338                return emod;
 339            }
 340        case kwPbKRatio: {
 341                if (!pbProps_) return 0.0;
 342                return (pbProps_->pb_ks_ == 0.0) ? 0.0 : (pbProps_->pb_kn_/pbProps_->pb_ks_);
 343            }
 344        }
 345        assert(0);
 346        return QVariant();
 347    }
 348
 349    bool ContactModelLinearPBond::getPropertyGlobal(uint32 i) const {
 350        switch (i) {
 351        case kwLinF:   
 352        case kwDpF:  
 353        case kwPbF:
 354            return false;
 355        }
 356        return true;
 357    }
 358
 359    bool ContactModelLinearPBond::setProperty(uint32 i,const QVariant &v,IContact *) {
 360        pbProps pb;
 361        dpProps dp;
 362
 363        switch (i) {
 364        case kwLinKn: {
 365                if (!v.canConvert<double>())
 366                    throw Exception("kn must be a double.");
 367                double val(v.toDouble());
 368                if (val<0.0)
 369                    throw Exception("Negative kn not allowed.");
 370                kn_ = val;
 371                return true;
 372            }
 373        case kwLinKs: {
 374                if (!v.canConvert<double>())
 375                    throw Exception("ks must be a double.");
 376                double val(v.toDouble());
 377                if (val<0.0)
 378                    throw Exception("Negative ks not allowed.");
 379                ks_ = val;  
 380                return true;
 381            }
 382        case kwLinFric: {
 383                if (!v.canConvert<double>())
 384                    throw Exception("fric must be a double.");
 385                double val(v.toDouble());
 386                if (val<0.0)
 387                    throw Exception("Negative fric not allowed.");
 388                fric_ = val;  
 389                return false;
 390            }
 391        case kwLinF: {
 392                if (!v.canConvert<DVect>())
 393                    throw Exception("lin_force must be a vector.");
 394                DVect val(v.value<DVect>());
 395                lin_F_ = val;
 396                return false;
 397            }
 398        case kwLinMode: {
 399                if (!v.canConvert<uint32>())
 400                    throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
 401                uint32 val(v.toUInt());
 402                if (val>1)
 403                    throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
 404                lin_mode_ = val;
 405                return false;
 406            }
 407        case kwRGap: {
 408                if (!v.canConvert<double>())
 409                    throw Exception("Reference gap must be a double.");
 410                double val(v.toDouble());
 411                rgap_ = val;  
 412                return false;
 413            }
 414        case kwPbRMul: {
 415                if (!v.canConvert<double>())
 416                    throw Exception("pb_rmul must be a double.");
 417                double val(v.toDouble());
 418                if (val<=0.0)
 419                    throw Exception("pb_rmul must be positive.");
 420                if (val == 1.0 && !pbProps_)
 421                    return false;
 422                if (!pbProps_)
 423                    pbProps_ = NEW pbProps();
 424                pbProps_->pb_rmul_ = val;
 425                return false;
 426            }
 427        case kwPbKn: {
 428                if (!v.canConvert<double>())
 429                    throw Exception("pb_kn must be a double.");
 430                double val(v.toDouble());
 431                if (val<0.0)
 432                    throw Exception("Negative pb_kn not allowed.");
 433                if (val == 0.0 && !pbProps_)
 434                    return false;
 435                if (!pbProps_)
 436                    pbProps_ = NEW pbProps();
 437                pbProps_->pb_kn_ = val;
 438                return true;
 439            }
 440        case kwPbKs: {
 441                if (!v.canConvert<double>())
 442                    throw Exception("pb_ks must be a double.");
 443                double val(v.toDouble());
 444                if (val<0.0)
 445                    throw Exception("Negative pb_ks not allowed.");
 446                if (val == 0.0 && !pbProps_)
 447                    return false;
 448                if (!pbProps_)
 449                    pbProps_ = NEW pbProps();
 450                pbProps_->pb_ks_ = val;
 451                return true;
 452            }
 453        case kwPbMcf: {
 454                if (!v.canConvert<double>())
 455                    throw Exception("pb_mcf must be a double.");
 456                double val(v.toDouble());
 457                if (val<0.0)
 458                    throw Exception("Negative pb_mcf not allowed.");
 459                if (val > 1.0)
 460                    throw Exception("pb_mcf must be lower or equal to 1.0.");
 461                if (val == 1.0 && !pbProps_)
 462                    return false;
 463                if (!pbProps_)
 464                    pbProps_ = NEW pbProps();
 465                pbProps_->pb_mcf_ = val;  
 466                return false;
 467            }
 468        case kwPbTStrength: {
 469                if (!v.canConvert<double>())
 470                    throw Exception("pb_ten must be a double.");
 471                double val(v.toDouble());
 472                if (val < 0.0)
 473                    throw Exception("Negative pb_ten not allowed.");
 474                if (val == 0.0 && !pbProps_)
 475                    return false;
 476                if (!pbProps_)
 477                    pbProps_ = NEW pbProps();
 478                pbProps_->pb_ten_ = val;  
 479                return false;
 480            }
 481        case kwPbCoh: {
 482                if (!v.canConvert<double>())
 483                    throw Exception("pb_coh must be a double.");
 484                double val(v.toDouble());
 485                if (val<0.0)
 486                    throw Exception("Negative pb_coh not allowed.");
 487                if (val == 0.0 && !pbProps_)
 488                    return false;
 489                if (!pbProps_)
 490                    pbProps_ = NEW pbProps();
 491                pbProps_->pb_coh_ = val;  
 492                return false;
 493            }
 494        case kwPbFa: {
 495                if (!v.canConvert<double>())
 496                    throw Exception("pb_fa must be a double.");
 497                double val(v.toDouble());
 498                if (val<0.0)
 499                    throw Exception("Negative pb_fa not allowed.");
 500                if (val>=90.0)
 501                    throw Exception("pb_fa must be lower than 90.0 degrees.");
 502                if (val == 0.0 && !pbProps_)
 503                    return false;
 504                if (!pbProps_)
 505                    pbProps_ = NEW pbProps();
 506                pbProps_->pb_fa_ = val;  
 507                return false;
 508            }
 509        case kwPbF: {
 510                if (!v.canConvert<DVect>())
 511                    throw Exception("pb_force must be a vector.");
 512                DVect val(v.value<DVect>());
 513                if (val.fsum() == 0.0 && !pbProps_)
 514                    return false;
 515                if (!pbProps_)
 516                    pbProps_ = NEW pbProps();
 517                pbProps_->pb_F_ = val;
 518                return false;
 519            }
 520        case kwPbM: {
 521                DAVect val(0.0);
 522#ifdef TWOD               
 523                if (!v.canConvert<DAVect>() && !v.canConvert<double>())
 524                    throw Exception("pb_moment must be an angular vector.");
 525                if (v.canConvert<DAVect>())
 526                    val = DAVect(v.value<DAVect>());
 527                else
 528                    val = DAVect(v.toDouble());
 529#else
 530                if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
 531                    throw Exception("pb_moment must be an angular vector.");
 532                if (v.canConvert<DAVect>())
 533                    val = DAVect(v.value<DAVect>());
 534                else
 535                    val = DAVect(v.value<DVect>());
 536#endif
 537                if (val.fsum() == 0.0 && !pbProps_)
 538                    return false;
 539                if (!pbProps_)
 540                    pbProps_ = NEW pbProps();
 541                pbProps_->pb_M_ = val;
 542                return false;
 543            }      
 544        case kwDpNRatio: {
 545                if (!v.canConvert<double>())
 546                    throw Exception("dp_nratio must be a double.");
 547                double val(v.toDouble());
 548                if (val<0.0)
 549                    throw Exception("Negative dp_nratio not allowed.");
 550                if (val == 0.0 && !dpProps_)
 551                    return false;
 552                if (!dpProps_)
 553                    dpProps_ = NEW dpProps();
 554                dpProps_->dp_nratio_ = val; 
 555                return true;
 556            }
 557        case kwDpSRatio: {
 558                if (!v.canConvert<double>())
 559                    throw Exception("dp_sratio must be a double.");
 560                double val(v.toDouble());
 561                if (val<0.0)
 562                    throw Exception("Negative dp_sratio not allowed.");
 563                if (val == 0.0 && !dpProps_)
 564                    return false;
 565                if (!dpProps_)
 566                    dpProps_ = NEW dpProps();
 567                dpProps_->dp_sratio_ = val;
 568                return true;
 569            }
 570        case kwDpMode: {
 571                if (!v.canConvert<int>())
 572                    throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
 573                int val(v.toInt());
 574                if (val == 0 && !dpProps_)
 575                    return false;
 576                if (val < 0 || val > 3)
 577                    throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
 578                if (!dpProps_)
 579                    dpProps_ = NEW dpProps();
 580                dpProps_->dp_mode_ = val;
 581                return false;
 582            }
 583        case kwDpF: {
 584                if (!v.canConvert<DVect>())
 585                    throw Exception("dp_force must be a vector.");
 586                DVect val(v.value<DVect>());
 587                if (val.fsum() == 0.0 && !dpProps_)
 588                    return false;
 589                if (!dpProps_)
 590                    dpProps_ = NEW dpProps();
 591                dpProps_->dp_F_ = val;
 592                return false;
 593            }
 594        case kwUserArea: {
 595                if (!v.canConvert<double>())
 596                    throw Exception("user_area must be a double.");
 597                double val(v.toDouble());
 598                if (val < 0.0)
 599                    throw Exception("Negative user_area not allowed.");
 600                userArea_ = val;
 601                return true;
 602            }
 603        }
 604//    assert(0);
 605        return false;
 606    }
 607
 608    bool ContactModelLinearPBond::getPropertyReadOnly(uint32 i) const {
 609        switch (i) {
 610        case kwDpF:
 611        case kwLinS:
 612        case kwEmod:
 613        case kwKRatio:
 614        case kwPbState:
 615        case kwPbRadius:
 616        case kwPbSStrength:
 617        case kwPbSig:
 618        case kwPbTau:
 619        case kwPbEmod:
 620        case kwPbKRatio:
 621            return true;
 622        default:
 623            break;
 624        }
 625        return false;
 626    }
 627
 628    bool ContactModelLinearPBond::supportsInheritance(uint32 i) const {
 629        switch (i) {
 630        case kwLinKn:
 631        case kwLinKs:
 632        case kwLinFric:
 633            return true;
 634        default:
 635            break;
 636        }
 637        return false;
 638    }
 639
 640    QString  ContactModelLinearPBond::getMethodArguments(uint32 i) const {
 641        QString def = QString();
 642        switch (i) {
 643        case kwDeformability:
 644            return "emod,kratio";
 645        case kwPbDeformability:
 646            return "emod,kratio";
 647        case kwPbBond: 
 648            return "gap";
 649        case kwPbUnbond: 
 650            return "gap";
 651        }
 652        return def;
 653    }
 654
 655    bool ContactModelLinearPBond::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
 656        IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 657        switch (i) {
 658        case kwDeformability: {
 659                double emod;
 660                double krat;
 661                if (vl.at(0).isNull()) 
 662                    throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
 663                emod = vl.at(0).toDouble();
 664                if (emod<0.0)
 665                    throw Exception("Negative emod not allowed in contact model %1.",getName());
 666                if (vl.at(1).isNull()) 
 667                    throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
 668                krat = vl.at(1).toDouble();
 669                if (krat<0.0)
 670                    throw Exception("Negative linear stiffness ratio not allowed in contact model %1.",getName());
 671                double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
 672                double rsum(0.0);
 673                if (c->getEnd1Curvature().y())
 674                    rsum += 1.0/c->getEnd1Curvature().y();
 675                if (c->getEnd2Curvature().y())
 676                    rsum += 1.0/c->getEnd2Curvature().y();
 677                if (userArea_) {
 678#ifdef THREED
 679                    rsq = std::sqrt(userArea_ / dPi);
 680#else
 681                    rsq = userArea_ / 2.0;
 682#endif        
 683                    rsum = rsq + rsq;
 684                    rsq = 1. / rsq;
 685                }
 686#ifdef TWOD
 687                kn_ = 2.0 * emod / (rsq * rsum);
 688#else
 689                kn_ = dPi * emod / (rsq * rsq * rsum);
 690#endif
 691                ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
 692                setInheritance(1,false);
 693                setInheritance(2,false);
 694                return true;
 695            }
 696        case kwPbDeformability: {
 697                //if (!pbProps_ || pbProps_->pb_state_ != 3) return false;
 698                double emod;
 699                double krat;
 700                if (vl.at(0).isNull()) 
 701                    throw Exception("Argument emod must be specified with method pb_deformability in contact model %1.",getName());
 702                emod = vl.at(0).toDouble();
 703                if (emod<0.0)
 704                    throw Exception("Negative emod not allowed in contact model %1.",getName());
 705                if (vl.at(1).isNull()) 
 706                    throw Exception("Argument kratio must be specified with method pb_deformability in contact model %1.",getName());
 707                krat = vl.at(1).toDouble();
 708                if (krat<0.0)
 709                    throw Exception("Negative parallel bond stiffness ratio not allowed in contact model %1.",getName());
 710                double rsum(0.0);
 711                if (c->getEnd1Curvature().y())
 712                    rsum += 1.0/c->getEnd1Curvature().y();
 713                if (c->getEnd2Curvature().y())
 714                    rsum += 1.0/c->getEnd2Curvature().y();
 715                if (!pbProps_)
 716                    pbProps_ = NEW pbProps();
 717                if (userArea_) 
 718#ifdef THREED
 719                    rsum = 2 * std::sqrt(userArea_ / dPi);
 720#else
 721                    rsum = 2 * userArea_ / 2.0;
 722#endif
 723                pbProps_->pb_kn_ = emod / rsum;
 724                pbProps_->pb_ks_ = (krat == 0.0) ? 0.0 : pbProps_->pb_kn_ / krat;
 725                return true;
 726            }
 727        case kwPbBond: {
 728                if (pbProps_ && pbProps_->pb_state_ == 3) return false;
 729                double mingap = -1.0 * limits<double>::max();
 730                double maxgap = 0;
 731                if (vl.at(0).canConvert<double>()) 
 732                    maxgap = vl.at(0).toDouble();
 733                else if (vl.at(0).canConvert<DVect2>()) {
 734                    DVect2 value = vl.at(0).value<DVect2>();
 735                    mingap = value.minComp();
 736                    maxgap = value.maxComp();
 737                } else if (!vl.at(0).isNull())
 738                    throw Exception("gap value %1 not recognized in method bond in contact model %2.",vl.at(0),getName());
 739                double gap = c->getGap(); 
 740                if ( gap >= mingap && gap <= maxgap) {
 741                    if (pbProps_)
 742                        pbProps_->pb_state_ = 3;
 743                    else {
 744                        pbProps_ = NEW pbProps();
 745                        pbProps_->pb_state_ = 3;
 746                    }
 747                    return true;
 748                }
 749                return false;
 750            }
 751        case kwPbUnbond: {
 752                if (!pbProps_ || pbProps_->pb_state_ == 0) return false;
 753                double mingap = -1.0 * limits<double>::max();
 754                double maxgap = 0;
 755                if (vl.at(0).canConvert<double>()) 
 756                    maxgap = vl.at(0).toDouble();
 757                else if (vl.at(0).canConvert<DVect2>()) {
 758                    DVect2 value = vl.at(0).value<DVect2>();
 759                    mingap = value.minComp();
 760                    maxgap = value.maxComp();
 761                }
 762                else if (!vl.at(0).isNull())
 763                    throw Exception("gap value %1 not recognized in method unbond in contact model %2.",vl.at(0),getName());
 764                double gap = c->getGap(); 
 765                if ( gap >= mingap && gap <= maxgap) {
 766                    pbProps_->pb_state_ = 0;
 767                    return true;
 768                }
 769                return false;
 770            }
 771        case kwArea: {
 772                if (!userArea_) {
 773                    double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
 774#ifdef THREED
 775                    userArea_ = rsq * rsq * dPi;
 776#else
 777                    userArea_ = rsq * 2.0;
 778#endif                            
 779                }
 780                return true;
 781            }
 782        }
 783        return false;
 784    }
 785
 786    double ContactModelLinearPBond::getEnergy(uint32 i) const {
 787        double ret(0.0);
 788        if (!energies_)
 789            return ret;
 790        switch (i) {
 791        case kwEStrain:  return energies_->estrain_;
 792        case kwESlip:    return energies_->eslip_;
 793        case kwEDashpot: return energies_->edashpot_;
 794        case kwEPbStrain:return energies_->epbstrain_;
 795        }
 796        assert(0);
 797        return ret;
 798    }
 799
 800    bool ContactModelLinearPBond::getEnergyAccumulate(uint32 i) const {
 801        switch (i) {
 802        case kwEStrain:   return false;
 803        case kwESlip:     return true;
 804        case kwEDashpot:  return true;
 805        case kwEPbStrain: return false;
 806        }
 807        assert(0);
 808        return false;
 809    }
 810
 811    void ContactModelLinearPBond::setEnergy(uint32 i,const double &d) {
 812        if (!energies_) return;
 813        switch (i) {
 814        case kwEStrain:   energies_->estrain_  = d; return;  
 815        case kwESlip:     energies_->eslip_    = d; return;
 816        case kwEDashpot:  energies_->edashpot_ = d; return;
 817        case kwEPbStrain: energies_->epbstrain_= d; return;
 818        }
 819        assert(0);
 820        return;
 821    }
 822
 823    bool ContactModelLinearPBond::validate(ContactModelMechanicalState *state,const double &) {
 824        assert(state);
 825        const IContactMechanical *c = state->getMechanicalContact(); 
 826        assert(c);
 827
 828        if (state->trackEnergy_)
 829            activateEnergy();
 830
 831        if (inheritanceField_ & linKnMask)
 832            updateKn(c);
 833        if (inheritanceField_ & linKsMask)
 834            updateKs(c);
 835        if (inheritanceField_ & linFricMask)
 836            updateFric(c);
 837
 838        updateEffectiveStiffness(state);
 839        return checkActivity(state->gap_);
 840    }
 841
 842    static const QString knstr("kn");
 843    bool ContactModelLinearPBond::updateKn(const IContactMechanical *con) {
 844        assert(con);
 845        QVariant v1 = con->getEnd1()->getProperty(knstr);
 846        QVariant v2 = con->getEnd2()->getProperty(knstr);
 847        if (!v1.isValid() || !v2.isValid())
 848            return false;
 849        double kn1 = v1.toDouble();
 850        double kn2 = v2.toDouble();
 851        double val = kn_;
 852        if (kn1 && kn2)
 853            kn_ = kn1*kn2/(kn1+kn2);
 854        else if (kn1)
 855            kn_ = kn1;
 856        else if (kn2)
 857            kn_ = kn2;
 858        return ( (kn_ != val) );
 859    }
 860
 861    static const QString ksstr("ks");
 862    bool ContactModelLinearPBond::updateKs(const IContactMechanical *con) {
 863        assert(con);
 864        QVariant v1 = con->getEnd1()->getProperty(ksstr);
 865        QVariant v2 = con->getEnd2()->getProperty(ksstr);
 866        if (!v1.isValid() || !v2.isValid())
 867            return false;
 868        double ks1 = v1.toDouble();
 869        double ks2 = v2.toDouble();
 870        double val = ks_;
 871        if (ks1 && ks2)
 872            ks_ = ks1*ks2/(ks1+ks2);
 873        else if (ks1)
 874            ks_ = ks1;
 875        else if (ks2)
 876            ks_ = ks2;
 877        return ( (ks_ != val) );
 878    }
 879
 880    static const QString fricstr("fric");
 881    bool ContactModelLinearPBond::updateFric(const IContactMechanical *con) {
 882        assert(con);
 883        QVariant v1 = con->getEnd1()->getProperty(fricstr);
 884        QVariant v2 = con->getEnd2()->getProperty(fricstr);
 885        if (!v1.isValid() || !v2.isValid())
 886            return false;
 887        double fric1 = std::max(0.0,v1.toDouble());
 888        double fric2 = std::max(0.0,v2.toDouble());
 889        double val = fric_;
 890        fric_ = std::min(fric1,fric2);
 891        return ( (fric_ != val) );
 892    }
 893
 894    bool ContactModelLinearPBond::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
 895        assert(c);
 896        QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
 897        QRegularExpression rx(name, QRegularExpression::CaseInsensitiveOption);
 898        int idx = availableProperties.indexOf(rx)+1;
 899        bool ret=false;
 900
 901        if (idx<=0)
 902            return ret;
 903         
 904        switch(idx) {
 905        case kwLinKn:  { //kn
 906                if (inheritanceField_ & linKnMask)
 907                    ret = updateKn(c);
 908                break;
 909            }
 910        case kwLinKs:  { //ks
 911                if (inheritanceField_ & linKsMask)
 912                    ret =updateKs(c);
 913                break;
 914            }
 915        case kwLinFric:  { //fric
 916                if (inheritanceField_ & linFricMask)
 917                    updateFric(c);
 918                break;
 919            }
 920        }
 921        return ret;
 922    }
 923
 924    void ContactModelLinearPBond::updateEffectiveStiffness(ContactModelMechanicalState *state) {
 925        DVect2 ret(kn_,ks_);
 926        // account for viscous damping
 927        if (dpProps_) {
 928            DVect2 correct(1.0);
 929            if (dpProps_->dp_nratio_)
 930                correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
 931            if (dpProps_->dp_sratio_)
 932                correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
 933            ret /= (correct*correct);
 934        }
 935        effectiveTranslationalStiffness_ = ret;
 936        if (isBonded()) {
 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 = pbProps_->pb_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 pbArea = dthick <= 0.0 ? dPi*br2 : 2.0*br*dthick;
 950            double bi = dthick <= 0.0 ? 0.25*pbArea*br2 : 2.0*br*br2*dthick/3.0;
 951            pbProps_->pbTransStiff_.rx() = pbProps_->pb_kn_*pbArea;
 952            pbProps_->pbTransStiff_.ry() = pbProps_->pb_ks_*pbArea;
 953#if DIM==3 
 954            pbProps_->pbAngStiff_.rx() = pbProps_->pb_ks_* 2.0 * bi;
 955            pbProps_->pbAngStiff_.ry() = pbProps_->pb_kn_* bi;
 956#endif
 957            pbProps_->pbAngStiff_.rz() = pbProps_->pb_kn_* bi;
 958        }
 959    }
 960     
 961    double ContactModelLinearPBond::pbStrainEnergy() const {
 962        double ret(0.0);
 963        if (pbProps_->pb_kn_)
 964            ret = 0.5 * pbProps_->pb_F_.x() * pbProps_->pb_F_.x() / pbProps_->pbTransStiff_.x();
 965        if (pbProps_->pb_ks_) {
 966            DVect tmp = pbProps_->pb_F_;
 967            tmp.rx() = 0.0;
 968            double smag2 = tmp.mag2();
 969            ret += 0.5 * smag2 / pbProps_->pbTransStiff_.y();
 970        }
 971
 972#ifdef THREED
 973        if (pbProps_->pbAngStiff_.x())
 974            ret += 0.5 * pbProps_->pb_M_.x() * pbProps_->pb_M_.x() / pbProps_->pbAngStiff_.x();
 975#endif
 976        if (pbProps_->pbAngStiff_.z()) {
 977            DAVect tmp = pbProps_->pb_M_;
 978#ifdef THREED
 979            tmp.rx() = 0.0;
 980            double smag2 = tmp.mag2();
 981#else
 982            double smag2 = tmp.z() * tmp.z();
 983#endif
 984            ret += 0.5 * smag2 / pbProps_->pbAngStiff_.z();
 985        }
 986        return ret;
 987    }
 988
 989    bool ContactModelLinearPBond::forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) {
 990        assert(state);
 991
 992        double overlap = rgap_ - state->gap_;
 993        DVect trans = state->relativeTranslationalIncrement_;
 994        double correction = 1.0;
 995
 996        if (state->activated()) {
 997            if (cmEvents_[fActivated] >= 0) {
 998                auto c = state->getContact();
 999                std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
1000                IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1001                fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
1002            }
1003            if (lin_mode_ == 0 && trans.x()) {
1004                correction = -1.0*overlap / trans.x();
1005                if (correction < 0)
1006                    correction = 1.0;
1007            }
1008        }
1009
1010#ifdef THREED
1011        DVect norm(trans.x(),0.0,0.0);
1012#else
1013        DVect norm(trans.x(),0.0);
1014#endif
1015        DAVect ang  = state->relativeAngularIncrement_;
1016        DVect ss_F_old = lin_F_;
1017
1018        if (lin_mode_==0)
1019            lin_F_.rx() = overlap * kn_;
1020        else
1021            lin_F_.rx() -= correction * norm.x() * kn_;
1022        lin_F_.rx() = std::max(0.0,lin_F_.x());
1023
1024        DVect u_s = trans;
1025        u_s.rx() = 0.0;
1026        DVect sforce = lin_F_ - u_s * ks_ * correction;
1027        sforce.rx() = 0.0;
1028        // Make sure that the shear force opposses the direction of translation - otherwise you can
1029        // have strange behavior
1030        //for (int j=1; j<dim; ++j)
1031        //{
1032        //  qDebug()<<sforce.dof(j)<<trans.dof(j);
1033        //  if (sign<double>(1,sforce.dof(j)) == sign<double>(1,trans.dof(j)))
1034        //    sforce.rdof(j) = 0.0;
1035        //}
1036
1037        // Resolve sliding
1038        if (state->canFail_) {
1039            double crit = lin_F_.x() * fric_;
1040            double sfmag = sforce.mag();
1041            if (sfmag > crit) {
1042                double rat = crit / sfmag;
1043                sforce *= rat;
1044                if (!lin_S_ && cmEvents_[fSlipChange] >= 0) {
1045                    auto c = state->getContact();
1046                    std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1047                                                         fish::Parameter()                };
1048                    IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1049                    fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1050                }
1051                lin_S_ = true;
1052            } else {
1053                if (lin_S_) {
1054                    if (cmEvents_[fSlipChange] >= 0) {
1055                        auto c = state->getContact();
1056                        std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1057                                                             fish::Parameter((qint64)1)        };
1058                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1059                        fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1060                    }
1061                    lin_S_ = false;
1062                }
1063            }
1064        }
1065
1066        sforce.rx() = lin_F_.x();
1067        lin_F_ = sforce;          // total force in linear contact model
1068         
1069        // Account for dashpot forces
1070        if (dpProps_) {
1071            dpProps_->dp_F_.fill(0.0);
1072            double vcn(0.0), vcs(0.0);
1073            setDampCoefficients(state->inertialMass_,&vcn,&vcs);
1074            // Need to change behavior based on the dp_mode
1075            if (dpProps_->dp_mode_ == 0)  { // Damp all components
1076                dpProps_->dp_F_ = u_s * (-1.0* vcs) / timestep; // shear component   
1077                dpProps_->dp_F_ -= norm * vcn / timestep;       // normal component
1078            } else if (dpProps_->dp_mode_ == 1)  { // limit the tensile
1079                dpProps_->dp_F_ -= norm * vcn / timestep;       // normal component
1080                if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
1081                    dpProps_->dp_F_.rx() = - lin_F_.rx();
1082            } else if (dpProps_->dp_mode_ == 2)  { // limit the shear
1083                if (!lin_S_)
1084                    dpProps_->dp_F_ = u_s * (-1.0* vcs) / timestep; // shear component 
1085            } else {
1086                if (!lin_S_)
1087                    dpProps_->dp_F_ = u_s * (-1.0* vcs) / timestep; // shear component   
1088                dpProps_->dp_F_ -= norm * vcn / timestep;       // normal component
1089                if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
1090                    dpProps_->dp_F_.rx() = - lin_F_.rx();
1091            }
1092        }
1093
1094        // Account for parallel bonds
1095        bool doPb = false;
1096        if (pbProps_ && pbProps_->pb_state_ > 2) {
1097            doPb = true;
1098            // Parallel Bond Logic:
1099            // bond area and inertia
1100            // minimal curvature of end1
1101            double Cmin1 = state->end1Curvature_.x();
1102            double Cmax1 = state->end1Curvature_.y();
1103            double Cmax2 = state->end2Curvature_.y();
1104            double dthick = (Cmin1 == 0.0) ? 1.0 : 0.0;    
1105            double br = pbProps_->pb_rmul_ * 1.0 / std::max(Cmax1,Cmax2);
1106            if (userArea_)
1107#ifdef THREED
1108                br = std::sqrt(userArea_ / dPi);
1109#else
1110                br = userArea_ / 2.0;
1111#endif
1112            double br2 = br*br;
1113            double pbArea = dthick <= 0.0 ? dPi*br2 : 2.0*br*dthick;
1114            double bi = dthick <= 0.0 ? 0.25*pbArea*br2 : 2.0*br*br2*dthick/3.0;
1115            pbProps_->pbTransStiff_.rx() = pbProps_->pb_kn_*pbArea;
1116            pbProps_->pbTransStiff_.ry() = pbProps_->pb_ks_*pbArea;
1117     
1118            /* elastic force increments */
1119            pbProps_->pb_F_ -= norm *(pbProps_->pb_kn_*pbArea) + u_s * (pbProps_->pb_ks_*pbArea);
1120
1121            /* elastic moment increments */
1122            //DAVect pbMomentInc(0.0);
1123#if DIM==3 
1124            pbProps_->pbAngStiff_.rx() = pbProps_->pb_ks_* 2.0 * bi;
1125            pbProps_->pbAngStiff_.ry() = pbProps_->pb_kn_* bi;
1126#endif
1127            pbProps_->pbAngStiff_.rz() = pbProps_->pb_kn_* bi;
1128
1129            DAVect pbMomentInc = ang * pbProps_->pbAngStiff_ *(-1.0);
1130            pbProps_->pb_M_ += pbMomentInc;
1131
1132            /* check bond failure */
1133            if (state->canFail_) {
1134                /* maximum stresses */
1135                double dbend = sqrt(pbProps_->pb_M_.y()*pbProps_->pb_M_.y() + pbProps_->pb_M_.z()*pbProps_->pb_M_.z());
1136                double dtwist = pbProps_->pb_M_.x();
1137                DVect bfs(pbProps_->pb_F_);
1138                bfs.rx() = 0.0;
1139                double dbfs = bfs.mag();
1140                double nsmax = -(pbProps_->pb_F_.x() / pbArea)  +  pbProps_->pb_mcf_ * dbend * br/bi;
1141                double ssmax = dbfs / pbArea  +  pbProps_->pb_mcf_ * std::abs(dtwist) * 0.5* br/bi;
1142                double ss ;
1143
1144                if (nsmax >= pbProps_->pb_ten_) {
1145                    // Failed in tension
1146                    double se = pbStrainEnergy(); // bond strain energy at the onset of failure
1147                    pbProps_->pb_state_ = 1;
1148                    pbProps_->pb_F_.fill(0.0);
1149                    pbProps_->pb_M_.fill(0.0);
1150                    if (cmEvents_[fBondBreak] >= 0) {
1151                        auto c = state->getContact();
1152                        std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1153                                                             fish::Parameter((qint64)pbProps_->pb_state_),
1154                                                             fish::Parameter(pbProps_->pb_ten_),
1155                                                             fish::Parameter(se) };
1156                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1157                        fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1158                    }
1159                } else if ((ss = pbShearStrength(pbArea)) <= ssmax) {
1160                    // Failed in shear
1161                    double se = pbStrainEnergy(); // bond strain energy at the onset of failure
1162                    pbProps_->pb_state_ = 2;
1163                    pbProps_->pb_F_.fill(0.0);
1164                    pbProps_->pb_M_.fill(0.0);
1165                    if (cmEvents_[fBondBreak] >= 0) {
1166                        auto c = state->getContact();
1167                        std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1168                                                             fish::Parameter((qint64)pbProps_->pb_state_),
1169                                                             fish::Parameter(ss),
1170                                                             fish::Parameter(se)                      };
1171                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1172                        fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1173                    }
1174                }
1175            }
1176        }
1177
1178        // Compute energies
1179        if (state->trackEnergy_) {
1180            assert(energies_);
1181            energies_->estrain_ = 0.0;
1182            energies_->epbstrain_ = 0.0;
1183            if (kn_)
1184                energies_->estrain_ = 0.5 * lin_F_.x()* lin_F_.x() /kn_;
1185            if (ks_) {
1186                DVect s = lin_F_;
1187                s.rx() = 0.0;
1188                double smag2 = s.mag2();
1189                energies_->estrain_ += 0.5* smag2 / ks_ ;
1190
1191                if (lin_S_) {
1192                    ss_F_old.rx() = 0.0;
1193                    DVect avg_F_s = (s + ss_F_old)*0.5;
1194                    DVect u_s_el =  (s - ss_F_old) / ks_;
1195                    energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
1196                }
1197            }
1198            if (dpProps_)
1199                energies_->edashpot_ -= dpProps_->dp_F_ | trans;
1200            if (doPb)
1201                energies_->epbstrain_ = pbStrainEnergy();
1202        }
1203
1204        assert(lin_F_ == lin_F_);
1205        return checkActivity(state->gap_);
1206    }
1207
1208    bool ContactModelLinearPBond::thermalCoupling(ContactModelMechanicalState *ms,ContactModelThermalState *ts, IContactThermal *ct,const double &) {
1209        // Account for thermal expansion in linear group in incremental mode
1210        bool ret = false;
1211        if (lin_mode_ > 0 && ts->gapInc_ > 0.0) {
1212            DVect finc(0.0);
1213            finc.rx() = kn_ * ts->gapInc_;
1214            lin_F_ -= finc;
1215            ret = true; 
1216        }
1217
1218        if (!pbProps_) return ret;
1219        if (pbProps_->pb_state_ < 3) return ret;
1220        int idx = ct->getModel()->getContactModel()->isProperty("thexp");
1221        if (idx<=0 ) return ret; 
1222
1223        double thexp = (ct->getModel()->getContactModel()->getProperty(idx)).toDouble();
1224        double length = ts->length_;
1225        double delTemp =ts->tempInc_;
1226        double delUn = length * thexp * delTemp;
1227        if (delUn == 0.0) return ret;
1228         
1229        double dthick = 0.0;
1230        double Cmin1 = ms->end1Curvature_.x();
1231        double Cmax1 = ms->end1Curvature_.y();
1232        double Cmin2 = ms->end2Curvature_.x();
1233        double Cmax2 = ms->end2Curvature_.y();
1234
1235        Cmin2;
1236        if (Cmin1 == 0.0)
1237            dthick = 1.0; 
1238
1239        double br = pbProps_->pb_rmul_ * 1.0 / std::max(Cmax1,Cmax2);
1240        if (userArea_)
1241#ifdef THREED
1242            br = std::sqrt(userArea_ / dPi);
1243#else
1244            br = userArea_ / 2.0;
1245#endif
1246        double br2 = br*br;
1247        double pbArea = dthick <= 0.0 ? dPi*br2 : 2.0*br*dthick;
1248        //
1249        DVect finc(0.0);
1250        finc.rx() = pbProps_->pb_kn_ * pbArea * delUn;
1251        pbProps_->pb_F_ += finc;
1252         
1253        ret = true;
1254        return ret;
1255    }
1256
1257    void ContactModelLinearPBond::setForce(const DVect &v,IContact *c) { 
1258        lin_F(v); 
1259        if (v.x() > 0) 
1260            rgap_ = c->getGap() + v.x() / kn_; 
1261    } 
1262
1263    void ContactModelLinearPBond::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
1264        // Only do something if the contact model is of the same type
1265        if (old->getContactModel()->getName().compare("linearpbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
1266            ContactModelLinearPBond *oldCm = (ContactModelLinearPBond *)old;
1267#ifdef THREED
1268            // Need to rotate just the shear component from oldSystem to newSystem
1269
1270            // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
1271            DVect axis = oldSystem.e1() & newSystem.e1();
1272            double c, ang, s;
1273            DVect re2;
1274            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1275                axis = axis.unit();
1276                c = oldSystem.e1()|newSystem.e1();
1277                if (c > 0)
1278                    c = std::min(c,1.0);
1279                else
1280                    c = std::max(c,-1.0);
1281                ang = acos(c);
1282                s = sin(ang);
1283                double t = 1. - c;
1284                DMatrix<3,3> rm;
1285                rm.get(0,0) = t*axis.x()*axis.x() + c;
1286                rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
1287                rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
1288                rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
1289                rm.get(1,1) = t*axis.y()*axis.y() + c;
1290                rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
1291                rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
1292                rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
1293                rm.get(2,2) = t*axis.z()*axis.z() + c;
1294                re2 = rm*oldSystem.e2();
1295            }
1296            else
1297                re2 = oldSystem.e2();
1298            // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
1299            axis = re2 & newSystem.e2();
1300            DVect2 tpf;
1301            DMatrix<2,2> m;
1302            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
1303                axis = axis.unit();
1304                c = re2|newSystem.e2();
1305                if (c > 0)
1306                    c = std::min(c,1.0);
1307                else
1308                    c = std::max(c,-1.0);
1309                ang = acos(c);
1310                if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
1311                    ang *= -1;
1312                s = sin(ang);
1313                m.get(0,0) = c;
1314                m.get(1,0) = s;
1315                m.get(0,1) = -m.get(1,0);
1316                m.get(1,1) = m.get(0,0);
1317                tpf = m*DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
1318            } else {
1319                m.get(0,0) = 1.;
1320                m.get(0,1) = 0.;
1321                m.get(1,0) = 0.;
1322                m.get(1,1) = 1.;
1323                tpf = DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
1324            }
1325            DVect pforce = DVect(0,tpf.x(),tpf.y());
1326#else
1327            oldSystem;
1328            newSystem;
1329            DVect pforce = DVect(0,oldCm->lin_F_.y());
1330#endif
1331            for (int i=1; i<dim; ++i)
1332                lin_F_.rdof(i) += pforce.dof(i);
1333            if (lin_mode_ && oldCm->lin_mode_)
1334                lin_F_.rx() = oldCm->lin_F_.x();
1335            oldCm->lin_F_ = DVect(0.0);
1336            if (dpProps_ && oldCm->dpProps_) {
1337#ifdef THREED
1338                tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
1339                pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
1340#else
1341                pforce = oldCm->dpProps_->dp_F_;
1342#endif
1343                dpProps_->dp_F_ += pforce;
1344                oldCm->dpProps_->dp_F_ = DVect(0.0);
1345            }
1346            if(oldCm->getEnergyActivated()) {
1347                activateEnergy();
1348                energies_->estrain_ = oldCm->energies_->estrain_;
1349                energies_->eslip_ = oldCm->energies_->eslip_;
1350                energies_->edashpot_ = oldCm->energies_->edashpot_;
1351                energies_->epbstrain_ = oldCm->energies_->epbstrain_;
1352                oldCm->energies_->estrain_ = 0.0;
1353                oldCm->energies_->edashpot_ = 0.0;
1354                oldCm->energies_->eslip_ = 0.0;
1355                oldCm->energies_->epbstrain_ = 0.0;
1356            }
1357            rgap_ = oldCm->rgap_;
1358        }
1359        assert(lin_F_ == lin_F_);
1360    }
1361
1362    void ContactModelLinearPBond::setNonForcePropsFrom(IContactModel *old) {
1363        // Only do something if the contact model is of the same type
1364        if (old->getName().compare("linearpbond",Qt::CaseInsensitive) == 0 && !isBonded()) {
1365            ContactModelLinearPBond *oldCm = (ContactModelLinearPBond *)old;
1366            kn_ = oldCm->kn_;
1367            ks_ = oldCm->ks_;
1368            fric_ = oldCm->fric_;
1369            lin_mode_ = oldCm->lin_mode_;
1370            rgap_ = oldCm->rgap_;
1371            userArea_ = oldCm->userArea_;
1372
1373            if (oldCm->dpProps_) {
1374                if (!dpProps_)
1375                    dpProps_ = NEW dpProps();
1376                dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
1377                dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
1378                dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
1379            }
1380
1381            if (oldCm->pbProps_) {
1382                if (!pbProps_)
1383                    pbProps_ = NEW pbProps();
1384                pbProps_->pb_rmul_ = oldCm->pbProps_->pb_rmul_;
1385                pbProps_->pb_kn_ = oldCm->pbProps_->pb_kn_;
1386                pbProps_->pb_ks_ = oldCm->pbProps_->pb_ks_; 
1387                pbProps_->pb_mcf_ = oldCm->pbProps_->pb_mcf_; 
1388                pbProps_->pb_fa_ = oldCm->pbProps_->pb_fa_; 
1389                pbProps_->pb_state_ = oldCm->pbProps_->pb_state_; 
1390                pbProps_->pb_coh_ = oldCm->pbProps_->pb_coh_; 
1391                pbProps_->pb_ten_ = oldCm->pbProps_->pb_ten_; 
1392                pbProps_->pbTransStiff_ = oldCm->pbProps_->pbTransStiff_; 
1393                pbProps_->pbAngStiff_ = oldCm->pbProps_->pbAngStiff_; 
1394            }
1395        }
1396    }
1397
1398    DVect ContactModelLinearPBond::getForce(const IContactMechanical *) const {
1399        DVect ret(lin_F_);
1400        if (dpProps_)
1401            ret += dpProps_->dp_F_;
1402        if (pbProps_)
1403            ret += pbProps_->pb_F_;
1404        return ret;
1405    }
1406
1407    DAVect ContactModelLinearPBond::getMomentOn1(const IContactMechanical *c) const {
1408        DVect force = getForce(c);
1409        DAVect ret(0.0);
1410        if (pbProps_)
1411            ret = pbProps_->pb_M_;
1412        c->updateResultingTorqueOn1Local(force,&ret);
1413        return ret;
1414    }
1415
1416    DAVect ContactModelLinearPBond::getMomentOn2(const IContactMechanical *c) const {
1417        DVect force = getForce(c);
1418        DAVect ret(0.0);
1419        if (pbProps_)
1420            ret = pbProps_->pb_M_;
1421        c->updateResultingTorqueOn2Local(force,&ret);
1422        return ret;
1423    }
1424
1425    DVect3 ContactModelLinearPBond::pbData(const IContactMechanical *c) const {
1426        double Cmax1 = c->getEnd1Curvature().y();
1427        double Cmax2 = c->getEnd2Curvature().y();
1428        double br = pbProps_->pb_rmul_ * 1.0 / std::max(Cmax1,Cmax2);
1429        if (userArea_)
1430#ifdef THREED
1431            br = std::sqrt(userArea_ / dPi);
1432#else
1433            br = userArea_ / 2.0;
1434#endif
1435        double br2 = br*br;
1436#ifdef TWOD
1437        double pbArea = 2.0*br;
1438        double bi = 2.0*br*br2/3.0;
1439#else
1440        double pbArea = dPi*br2;
1441        double bi = 0.25*pbArea*br2;
1442#endif
1443        return DVect3(pbArea,bi,br);
1444    }
1445
1446    DVect2 ContactModelLinearPBond::pbSMax(const IContactMechanical *c) const {
1447        DVect3 data = pbData(c);
1448        double pbArea = data.x();
1449        double bi = data.y();
1450        double br = data.z();
1451        /* maximum stresses */
1452        double dbend = sqrt(pbProps_->pb_M_.y()*pbProps_->pb_M_.y() + pbProps_->pb_M_.z()*pbProps_->pb_M_.z());
1453        double dtwist = pbProps_->pb_M_.x();
1454        DVect bfs(pbProps_->pb_F_);
1455        bfs.rx() = 0.0;
1456        double dbfs = bfs.mag();
1457        double nsmax = -(pbProps_->pb_F_.x() / pbArea)  +  pbProps_->pb_mcf_ * dbend * br/bi;
1458        double ssmax = dbfs / pbArea  +  pbProps_->pb_mcf_ * std::abs(dtwist) * 0.5* br/bi;
1459        return DVect2(nsmax,ssmax);
1460    }
1461
1462    double ContactModelLinearPBond::pbShearStrength(const double &pbArea) const {
1463        if (!pbProps_) return 0.0;
1464        double sig = -1.0*pbProps_->pb_F_.x() / pbArea;
1465        double nstr = pbProps_->pb_state_ > 2 ? pbProps_->pb_ten_ : 0.0;
1466        return sig <= nstr ? pbProps_->pb_coh_ - std::tan(dDegrad*pbProps_->pb_fa_)*sig
1467                            : pbProps_->pb_coh_ - std::tan(dDegrad*pbProps_->pb_fa_)*nstr;
1468    }
1469
1470    void ContactModelLinearPBond::setDampCoefficients(const double &mass,double *vcn,double *vcs) {
1471        *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(kn_));
1472        *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(ks_));
1473    }
1474
1475} // namespace itascaxd
1476// EoF

Top