Flat-Joint Model Implementation

See this page for the documentation of this contact model.

contactmodelflatjoint.h

  1#pragma once
  2// contactmodelflatjoint.h
  3
  4#include "contactmodel/src/contactmodelmechanical.h"
  5
  6#ifdef FLATJOINT_LIB
  7#  define FLATJOINT_EXPORT EXPORT_TAG
  8#elif defined(NO_MODEL_IMPORT)
  9#  define FLATJOINT_EXPORT
 10#else
 11#  define FLATJOINT_EXPORT IMPORT_TAG
 12#endif
 13
 14namespace cmodelsxd {
 15    using namespace itasca;
 16
 17    class ContactModelFlatJoint : public ContactModelMechanical {
 18    public:
 19        enum PropertyKeys { 
 20              kwFjNr=1
 21            , kwFjElem
 22            , kwFjKn
 23            , kwFjKs                            
 24            , kwFjFric   
 25            , kwFjEmod
 26            , kwFjKRatio                            
 27            , kwFjRmul
 28            , kwFjRadius
 29            , kwFjGap0
 30            , kwFjTen 
 31            , kwFjCoh
 32            , kwFjFa 
 33            , kwFjF
 34            , kwFjM
 35            , kwFjState
 36            , kwFjSlip
 37            , kwFjMType
 38            , kwFjA
 39            , kwFjEgap
 40            , kwFjGap
 41            , kwFjNstr
 42            , kwFjSstr
 43            , kwFjSs
 44#ifdef THREED
 45            , kwFjNa
 46#endif
 47            , kwFjRelBr
 48            , kwFjCen
 49            , kwFjTrack
 50            , kwUserArea
 51            , kwFjCohRes
 52            , kwFjResMode
 53        };
 54         
 55        FLATJOINT_EXPORT ContactModelFlatJoint();
 56        FLATJOINT_EXPORT virtual ~ContactModelFlatJoint();
 57        void     copy(const ContactModel *c) override;
 58        void     archive(ArchiveStream &) override;
 59        QString  getName() const override { return "flatjoint"; }
 60        void     setIndex(int i) override { index_=i;}
 61        int      getIndex() const override {return index_;}
 62        QString  getProperties() const override { return "fj_nr"
 63                                                        ",fj_elem"
 64                                                        ",fj_kn"
 65                                                        ",fj_ks"
 66                                                        ",fj_fric"
 67                                                        ",fj_emod"
 68                                                        ",fj_kratio"
 69                                                        ",fj_rmul"
 70                                                        ",fj_radius"
 71                                                        ",fj_gap0"
 72                                                        ",fj_ten"
 73                                                        ",fj_coh"
 74                                                        ",fj_fa"
 75                                                        ",fj_force"
 76                                                        ",fj_moment"
 77                                                        ",fj_state"
 78                                                        ",fj_slip"
 79                                                        ",fj_mtype"
 80                                                        ",fj_area"
 81                                                        ",fj_egap"
 82                                                        ",fj_gap"
 83                                                        ",fj_sigma"
 84                                                        ",fj_tau"
 85                                                        ",fj_shear"
 86#ifdef THREED
 87                                                        ",fj_nal"
 88#endif
 89                                                        ",fj_relbr"
 90                                                        ",fj_cen"
 91                                                        ",fj_track"
 92                                                        ",user_area"
 93                                                        ",fj_cohres"
 94                                                        ",fj_resmode"
 95                                                        ;}
 96
 97        enum EnergyKeys { kwEStrain=1,kwESlip};
 98        QString  getEnergies() const { return "energy-strain,energy-slip";}
 99        double   getEnergy(uint32 i) const override;  // Base 1
100        bool     getEnergyAccumulate(uint32 i) const override; // Base 1
101        void     setEnergy(uint32 i,const double &d) override; // Base 1
102        void     activateEnergy() override { if (energies_) return; energies_ = NEWC(Energies());}
103        bool     getEnergyActivated() const override {return (energies_ !=0);}
104
105        enum FishCallEvents {fActivated=0,fBondBreak,fBroken,fSlipChange};
106        QString  getFishCallEvents() const override { return "contact_activated,bond_break,broken,all_slip_change"; }
107        QVariant getProperty(uint32 i,const IContact *) const override;
108        bool     getPropertyGlobal(uint32 i) const override;
109        bool     setProperty(uint32 i,const QVariant &v,IContact *) override;
110        bool     getPropertyReadOnly(uint32 i) const override;
111
112        bool     supportsInheritance(uint32 ) const override { return false; }
113
114        enum MethodKeys { kwBond=1, kwUnbond, KwDeformability, KwUpdateGeom, kwArea, kwInitialize};
115
116        QString  getMethods() const override { return "bond"
117                                                     ",unbond"
118                                                     ",deformability"
119                                                     ",update_geometry"
120                                                     ",area"
121                                                     ",initialize"
122                                            ;}
123        
124        QString  getMethodArguments(uint32 i) const override;
125        
126        bool     setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override; // Base 1 - returns true if timestep contributions need to be updated
127
128        uint32     getMinorVersion() const override;
129
130        bool    validate(ContactModelMechanicalState *state,const double &timestep) override;
131        bool    endPropertyUpdated(const QString &,const IContactMechanical *) override { return false; }
132        bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) override;
133        bool    thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
134        bool    coupling(ContactModelMechanicalState*, ContactModelThermalState*,ContactModelFluidState*, const double&) override;
135        DVect2 getEffectiveTranslationalStiffness() const override {
136              return effectiveTranslationalStiffness();
137        }
138        DAVect  getEffectiveRotationalStiffness() const override { return effectiveRotationalStiffness(); }
139
140        ContactModelFlatJoint *clone() const override { return NEWC(ContactModelFlatJoint()); }
141        double  getActivityDistance() const override {return 0.0;}
142        bool    isOKToDelete() const override { return !isBonded(); }
143        void    resetForcesAndMoments() override { fj_f(DVect(0.0)); fj_m(DAVect(0.0)); for (int i=0; i<f_.size(); ++i) f_[i] = DVect(0.0); }
144        void    setForce(const DVect &v,IContact *) override;
145        void    setArea(const double &d) override { userArea_ = d; }
146        double  getArea() const override { return userArea_; }
147        bool    checkActivity(const double &inGap) override;
148
149        /// Return the total force that the contact model holds.
150        DVect    getForce(const IContactMechanical*) const override;
151        /// Return the total moment on 1 that the contact model holds
152        DAVect   getMomentOn1(const IContactMechanical*) const override;
153        /// Return the total moment on 1 that the contact model holds
154        DAVect   getMomentOn2(const IContactMechanical*) const override;
155        //virtual bool     isSliding() const { return fj_s_; }
156        bool    isBonded() const override { FOR(it,bmode_) if ((*it) == 3) return true; return false; }
157        void    unbond() override { FOR(it,bmode_) *it = 0; }
158
159        // For contact specific plotting
160        void getSphereList(const IContact* con, std::vector<DVect>* pos, std::vector<double>* rad, std::vector<double>* val) override;
161#ifdef THREED
162        void getDiskList(const IContact* con, std::vector<DVect>* pos, std::vector<DVect>* normal, std::vector<double>* radius, std::vector<double>* val) override;
163#endif
164        void getCylinderList(const IContact* con, std::vector<DVect>* bot, std::vector<DVect>* top, std::vector<double>* radlow, std::vector<double>* radhi, std::vector<double>* val) override;
165
166        int             fj_nr() const               {return fj_nr_;}
167        void            fj_nr(int d)                {       fj_nr_= d;}
168#ifdef THREED
169        int             fj_n() const                { return fj_na_ * fj_nr_; }
170        int             fj_na() const               {return fj_na_;}
171        void            fj_na(int d)                {       fj_na_= d;}
172#else
173        int             fj_n() const                { return fj_nr_; }
174#endif
175        int             fj_elem() const             {return fj_elem_;}
176        void            fj_elem(int d)              {       fj_elem_= d;}
177        const double &  fj_kn() const               {return fj_kn_;}
178        void            fj_kn(const double &d)      {       fj_kn_ = d;}
179        const double &  fj_ks() const               {return fj_ks_;}
180        void            fj_ks(const double &d)      {       fj_ks_ = d;}
181        const double &  fj_fric() const             {return fj_fric_;}
182        void            fj_fric(const double &d)    {       fj_fric_ = d;}
183        const double &  fj_rmul() const             {return fj_rmul_;}
184        void            fj_rmul(const double &d)    {       fj_rmul_ = d;}
185        const double &  fj_gap0() const             {return fj_gap0_;}
186        void            fj_gap0(const double &d)    {       fj_gap0_ = d;}
187        const double &  fj_ten() const              {return fj_ten_;}
188        void            fj_ten(const double &d)     {       fj_ten_ = d;}
189        const double &  fj_coh() const              {return fj_coh_;}
190        void            fj_coh(const double &d)     {       fj_coh_ = d;}
191        const double &  fj_cohres() const           {return fj_cohres_;}
192        void            fj_cohres(const double &d)  {       fj_cohres_ = d;}
193        const double &  fj_fa() const               {return fj_fa_;}
194        void            fj_fa(const double &d)      {       fj_fa_ = d;}
195        const DVect &   fj_f() const                {return fj_f_;}
196        void            fj_f(const DVect &f)        {       fj_f_=f;}
197        const DAVect &  fj_m() const                {return fj_m_;}
198        void            fj_m(const DAVect &f)       {       fj_m_=f;}
199        const DAVect &  fj_m_set() const            {return fj_m_set_;}
200        void            fj_m_set(const DAVect &f)   {       fj_m_set_=f;}
201        const double &  rmin() const                {return rmin_;}
202        void            rmin(const double &d)       {       rmin_ = d;}
203        const double &  rbar() const                {return rbar_;}
204        void            rbar(const double &d)       {       rbar_ = d;}
205        const int &     fj_resmode() const          {return fj_resmode_;}
206        void            fj_resmode(const int &i)    {       fj_resmode_ = i;}
207        const double &  atot() const                {return atot_;}
208        void            atot(const double &d)       {       atot_ = d;}
209        bool            propsFixed() const          {return propsFixed_; }
210        void            propsFixed(bool d)          {       propsFixed_ = d;}
211        int             mType() const               {return mType_; }
212        void            mType(int d)                {       mType_ = d;}
213        const DVect &   gap() const                 {return gap_; }
214        void            gap(const DVect &d)         {       gap_ = d;}
215        const double &  theta() const               {return theta_; }
216        void            theta(const double & d)     {       theta_ = d;}
217#ifdef THREED
218        const double &  thetaM() const              {return thetaM_; }
219        void            thetaM(const double & d)    {       thetaM_ = d;}
220#else
221        double thetaM() const                       {return 0.0;}
222#endif
223
224        bool    hasEnergies() const {return energies_ ? true:false;}
225        double  estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
226        void    estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
227        double  eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
228        void    eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
229
230        uint32 inheritanceField() const {return inheritanceField_;}
231        void inheritanceField(uint32 i) {inheritanceField_ = i;}
232
233        const DVect2 & effectiveTranslationalStiffness()  const             {return effectiveTranslationalStiffness_;}
234        void           effectiveTranslationalStiffness(const DVect2 &v )    {effectiveTranslationalStiffness_=v;}
235        const DAVect & effectiveRotationalStiffness()  const                {return effectiveRotationalStiffness_;}
236        void           effectiveRotationalStiffness(const DAVect &v )       {effectiveRotationalStiffness_=v;}
237
238    private:
239        static int index_;
240
241        struct Energies {
242            Energies() : estrain_(0.0), eslip_(0.0) {}
243            double estrain_;  // elastic energy stored in contact 
244            double eslip_;    // work dissipated by friction 
245        };
246
247        void   updateEffectiveStiffness(ContactModelMechanicalState *state);
248
249        // inheritance fields
250        uint32 inheritanceField_;
251
252        int                     fj_nr_;             // radial number of elements >= 1 (total in 2D)
253#ifdef THREED
254        int                     fj_na_;             // circumferential number of elements >= 3
255#endif
256        int                     fj_elem_;           // Element to be queried
257        double                  fj_kn_;             // normal stiffness
258        double                  fj_ks_;             // shear stiffness
259        double                  fj_fric_;           // Coulomb friction coefficient
260        double                  fj_rmul_;           // Radius multiplier
261        double                  fj_gap0_;           // Initial gap
262        double                  fj_ten_;            // Tensile strength 
263        double                  fj_coh_;            // Cohesive strength
264        double                  fj_cohres_;         // Residual cohesive strength
265        double                  fj_fa_;             // Friction angle 
266        DVect                   fj_f_;              // Force carried in the model
267        DAVect                  fj_m_;              // Moment carried in the model
268        DAVect                  fj_m_set_;          // When initializing forces then need an extra moment term
269        // Area related quantities
270        double                  rmin_;              // min(Ra,Rb) where Ra & Rb are particle radii
271        double                  rbar_;              // flat-joint radius [m]
272        double                  atot_;              // flat-joint area [m^2]
273        std::vector<double>     a_;                 // cross-sectional area of elem[fj_elem-1] [m^2]
274#ifdef THREED
275        std::vector<DVect2>     rBarl_;             // centroid relative position of elem[fj_elem-1] [m] (3D)
276#else
277        std::vector<double>     rBarl_;             // centroid relative position of elem[fj_elem-1] [m] (2D)
278#endif
279        int                     fj_resmode_;         // Residual mode
280        void setAreaQuantities();                   // Set Rbar, Atot and A[]
281        DVect getRelElemPos(const IContact*,int ) const;   // Return the relative location of element i
282        void setRelElemPos(const IContact*,int ,const DVect &);   // Set the relative location of element i
283
284        bool                    propsFixed_;        // {Rmul, N, G, bstate, mType} fixed, cannot reset
285        int                     mType_;             // initial microstructural type
286        int getmType() const;                       // {1,2,3,4}={bonded, gapped, slit, other}
287        
288        std::vector<int>        bmode_;             // bond mode - 0 unbonded, 1 failed in tension, 2 failed in shear, 3 bonded
289        std::vector<bool>       smode_;             // slip mode
290        bool Bonded(int e) const { return bmode_[e-1] == 3 ? true : false; }
291
292        // Set bstate and bmode (can only bond if fj_gap0_==0.0)
293        void bondElem(int iSeg,bool bBond);
294        // Set bstate & bmode 
295        void breakBond(int iSeg,int fmode,ContactModelMechanicalState *state);
296        void slipChange(int iSeg,bool smode,ContactModelMechanicalState *state);
297
298        // For use in 2D only!
299        double tauC(const double &dSig,bool bBonded) const; // shear strength (positive) [N/m^2]
300
301        // INTERFACE RESPONSE QUANTITIES:
302        DVect                   gap_;               // total relative displacement [m]
303        double                  theta_;             // total relative rotation [rad]
304#ifdef THREED
305        double                  thetaM_;            // total relative rotation [rad]
306        double thbMag() const   { return sqrt(theta_*theta_ + thetaM_*thetaM_); }
307        // unit-vector xi of middle surface system xi-eta
308        // (If both thb_l and thb_m are zero, then xi is undefined
309        // and returns zero for both components.)
310        double xi(int comp /* component (l,m) = (1,2) */) const;
311#endif
312        std::vector<double>     egap_;          // gap at centroid of elem[fj_elem-1] [N]
313        std::vector<DVect>      f_;             // force on elem[fj_elem-1] [N]
314
315        void   initVectors();                   // Resize and zero all vector types based on current value of N
316#ifdef TWOD
317        double gap(const double &x) const;      // Gap (g>0 is open) along the interface, x in [0, 2*Rbar]
318#else
319        double gap(const double &rl,const double &rm) const; // Gap (g>0 is open) gap at relative position (l,m) [m]
320        double sigBar( int e /* element, e = 1,2,...,Nel */ ) const; // normal stress at centroid of elem[eN-1] [N/m^2]
321        double tauBar( int e /* element, e = 1,2,...,Nel */ ) const; // shear  stress at centroid of elem[eN-1] [N/m^2]
322#endif
323        double computeStrainEnergy(int e /* element, e = 1,2,...,Nel */) const; // strain energy in elem[eN-1]
324        // For use in 2D only! Segment normal stress
325        double computeSig(const double &g0,   /* gap at left end  */
326                          const double &g1,   /* gap at right end */
327                          const double &rbar, /* length is 2*rbar */
328                          const double &dA,   /* area             */
329                          bool bBonded        /* bond state       */
330                          ) const;
331        // For use in 2D only! Segment moment
332        double computeM(const double &g0,   /* gap at left end  */ 
333                        const double &g1,   /* gap at right end */ 
334                        const double &rbar, /* length is 2*rbar */
335                        bool bBonded        /* bond state       */
336                        ) const;
337        // For use in 2D only! getCase used by ComputeSig and ComputeM
338        int getCase(const double &g0, /* gap at left end  */ 
339                 const double &g1  /* gap at right end */ 
340                 ) const;
341        // Segment elastic shear-displacement increment, which is portion of
342        // increment that occurs while gap is negative.
343        double delUse(const double &gapStart, /* gap at start of FDlaw  */
344                      const double &gapEnd,   /* gap at end of FDlaw    */
345                      bool bBonded,           /* bond state             */
346                      const double &delUs     /* shear displ. increment */
347                     ) const;
348        double      userArea_;   // Area as specified by the user 
349        Energies *   energies_;    // energies
350
351        DVect2  effectiveTranslationalStiffness_;
352        DAVect  effectiveRotationalStiffness_;
353
354        struct orientProps {
355            orientProps() : origNormal_(DVect(0.0)) {}
356            Quat    orient1_;
357            Quat    orient2_;
358            DVect   origNormal_;
359        };
360
361        orientProps *orientProps_;
362    };
363} // namespace itascaxd
364
365
366// EoF

Top

contactmodelflatjoint.cpp

   1// contactmodelflatjoint.cpp
   2#include "contactmodelflatjoint.h"
   3
   4#include "../version.txt"
   5#include "fish/src/parameter.h"
   6#include "utility/src/tptr.h"
   7#include "shared/src/mathutil.h"
   8#include "base/src/basetoqt.h"
   9#include "contactmodel/src/contactmodelthermal.h"
  10#include "contactmodel/src/contactmodelfluid.h"
  11
  12#include "kernel/interface/iprogram.h"
  13#include "module/interface/icontact.h"
  14#include "module/interface/icontactmechanical.h"
  15#include "module/interface/icontactthermal.h"
  16#include "module/interface/icontactfluid.h"
  17#include "module/interface/ifishcalllist.h"
  18#include "module/interface/ipiece.h"
  19#include "module/interface/ipiecemechanical.h"
  20
  21#ifdef FLATJOINT_LIB
  22#ifdef _WIN32
  23  int __stdcall DllMain(void *,unsigned, void *)
  24  {
  25    return 1;
  26  }
  27#endif
  28
  29  extern "C" EXPORT_TAG const char *getName()
  30  {
  31#if DIM==3
  32    return "contactmodelmechanical3dflatjoint";
  33#else
  34    return "contactmodelmechanical2dflatjoint";
  35#endif
  36  }
  37
  38  extern "C" EXPORT_TAG unsigned getMajorVersion()
  39  {
  40    return MAJOR_VERSION;
  41  }
  42
  43  extern "C" EXPORT_TAG unsigned getMinorVersion()
  44  {
  45    return MINOR_VERSION;
  46  }
  47
  48  extern "C" EXPORT_TAG void *createInstance()
  49  {
  50    cmodelsxd::ContactModelFlatJoint *m = NEWC(cmodelsxd::ContactModelFlatJoint());
  51    return (void *)m;
  52  }
  53#endif // FLATJOINT_LIB
  54
  55namespace cmodelsxd {
  56    static const uint32 fjKnMask      = 0x00002; // Base 1!
  57    static const uint32 fjKsMask      = 0x00004;
  58    static const uint32 fjFricMask    = 0x00008;
  59
  60    using namespace itasca;
  61
  62    int ContactModelFlatJoint::index_ = -1;
  63    uint32 ContactModelFlatJoint::getMinorVersion() const { return MINOR_VERSION;}
  64
  65    ContactModelFlatJoint::ContactModelFlatJoint() : inheritanceField_(fjKnMask|fjKsMask|fjFricMask)
  66                                            , fj_nr_(2)
  67#ifdef THREED
  68                                            , fj_na_(4)
  69#endif
  70                                            , fj_elem_(1)
  71                                            , fj_kn_(0.0)
  72                                            , fj_ks_(0.0)
  73                                            , fj_fric_(0.0)
  74                                            , fj_rmul_(1.0)
  75                                            , fj_gap0_(0.0)
  76                                            , fj_ten_(0.0)
  77                                            , fj_coh_(0.0)
  78                                            , fj_cohres_(0.0)
  79                                            , fj_fa_(0.0)
  80                                            , fj_f_(0.0)
  81                                            , fj_m_(0.0)
  82                                            , fj_m_set_(0.0)
  83                                            , rmin_(1.0)
  84                                            , rbar_(0.0)
  85                                            , atot_(0.0)
  86                                            , a_(2)
  87                                            , rBarl_(2)
  88                                            , fj_resmode_(0)
  89                                            , propsFixed_(false)
  90                                            , mType_(3)
  91                                            , bmode_(2)
  92                                            , smode_(2)
  93                                            , gap_(0.0)
  94                                            , theta_(0.0)
  95#ifdef THREED
  96                                            , thetaM_(0.0)
  97#endif
  98                                            , egap_(2)
  99                                            , f_(2)
 100                                            , userArea_(0)
 101                                            , energies_(0)
 102                                            , effectiveTranslationalStiffness_(DVect2(0.0))
 103                                            , effectiveRotationalStiffness_(DAVect(0.0))
 104                                            , orientProps_(0)
 105    {
 106        initVectors();
 107        setAreaQuantities();
 108        //setFromParent(ContactModelMechanicalList::instance()->find(getName()));
 109    }
 110
 111    ContactModelFlatJoint::~ContactModelFlatJoint() {
 112        if (orientProps_)
 113            delete orientProps_;
 114        if (energies_)
 115            delete energies_;
 116    }
 117
 118    void ContactModelFlatJoint::archive(ArchiveStream &stream) {
 119        stream & fj_nr_;
 120#ifdef THREED
 121        stream & fj_na_;
 122#endif
 123        stream & fj_elem_;
 124        stream & fj_kn_;
 125        stream & fj_ks_;
 126        stream & fj_fric_;
 127        stream & fj_rmul_;
 128        stream & fj_gap0_;
 129        stream & fj_ten_;
 130        stream & fj_coh_;
 131        stream & fj_fa_;
 132        stream & fj_f_;
 133        stream & fj_m_;
 134        stream & rmin_;
 135        stream & rbar_;
 136        stream & atot_;
 137        stream & a_;
 138        stream & rBarl_;
 139        stream & propsFixed_;
 140        stream & mType_;
 141        stream & bmode_;
 142        stream & smode_;
 143        stream & gap_;
 144        stream & theta_;
 145#ifdef THREED
 146        stream & thetaM_;
 147#endif
 148        stream & egap_;
 149        stream & f_;
 150
 151        if (stream.getArchiveState()==ArchiveStream::Save) {
 152            bool b = false;
 153            if (orientProps_) {
 154                b = true;
 155                stream & b;
 156                stream & orientProps_->orient1_;
 157                stream & orientProps_->orient2_;
 158                stream & orientProps_->origNormal_;
 159            } else
 160                stream & b;
 161            b = false;
 162            if (energies_) {
 163                b = true;
 164                stream & b;
 165                stream & energies_->estrain_;
 166                stream & energies_->eslip_;
 167            } else
 168                stream & b;
 169        } else {
 170            bool b(false);
 171            stream & b;
 172            if (b) {
 173                if (!orientProps_)
 174                    orientProps_ = NEWC(orientProps());
 175                stream & orientProps_->orient1_;
 176                stream & orientProps_->orient2_;
 177                stream & orientProps_->origNormal_;
 178            }
 179            stream & b;
 180            if (b) {
 181                if (!energies_)
 182                    energies_ = NEWC(Energies());
 183                stream & energies_->estrain_;
 184                stream & energies_->eslip_;
 185            }
 186        }
 187
 188        stream & inheritanceField_;
 189        stream & effectiveTranslationalStiffness_;
 190        stream & effectiveRotationalStiffness_;
 191
 192        if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() > 1)
 193            stream & userArea_;
 194
 195        if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() > 2)
 196            stream & fj_m_set_;
 197
 198        if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() > 3) {
 199            stream & fj_cohres_;
 200            stream & fj_resmode_;
 201        }
 202    }
 203
 204    void ContactModelFlatJoint::copy(const ContactModel *cm) {
 205        ContactModelMechanical::copy(cm);
 206        const ContactModelFlatJoint *in = dynamic_cast<const ContactModelFlatJoint*>(cm);
 207        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
 208        fj_nr(in->fj_nr());
 209#ifdef THREED
 210        fj_na(in->fj_na());
 211#endif
 212        fj_elem(in->fj_elem());
 213        fj_kn(in->fj_kn());
 214        fj_ks(in->fj_ks());
 215        fj_fric(in->fj_fric());
 216        fj_rmul(in->fj_rmul());
 217        fj_gap0(in->fj_gap0());
 218        fj_ten(in->fj_ten());
 219        fj_coh(in->fj_coh());
 220        fj_cohres(in->fj_cohres());
 221        fj_fa(in->fj_fa());
 222        fj_f(in->fj_f());
 223        fj_m(in->fj_m());
 224        fj_m_set(in->fj_m_set());
 225        rmin(in->rmin());
 226        rbar(in->rbar());
 227        fj_resmode(in->fj_resmode());
 228        atot(in->atot());
 229        a_ = in->a_;
 230        rBarl_ = in->rBarl_;
 231        propsFixed(in->propsFixed());
 232        mType(in->mType());
 233        bmode_ = in->bmode_;
 234        smode_ = in->smode_;
 235        gap(in->gap());
 236        theta(in->theta());
 237#ifdef THREED
 238        thetaM(in->thetaM());
 239#endif
 240        egap_ = in->egap_;
 241        f_ = in->f_;
 242        if (in->orientProps_) {
 243            if (!orientProps_)
 244                orientProps_ = NEWC(orientProps());
 245            orientProps_->orient1_ = in->orientProps_->orient1_;
 246            orientProps_->orient2_ = in->orientProps_->orient2_;
 247            orientProps_->origNormal_ =  in->orientProps_->origNormal_;
 248        }
 249        if (in->hasEnergies()) {
 250            if (!energies_)
 251                energies_ = NEWC(Energies());
 252            estrain(in->estrain());
 253            eslip(in->eslip());
 254        }
 255        userArea_ = in->userArea_;
 256        inheritanceField(in->inheritanceField());
 257        effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
 258        effectiveRotationalStiffness(in->effectiveRotationalStiffness());
 259    }
 260
 261
 262    QVariant ContactModelFlatJoint::getProperty(uint32 i,const IContact *con) const {
 263        QVariant var;
 264        switch (i) {
 265        case kwFjNr     :   return fj_nr();
 266        case kwFjElem   :   return fj_elem();
 267        case kwFjKn     :   return fj_kn();
 268        case kwFjKs     :   return fj_ks();
 269        case kwFjFric   :   return fj_fric();
 270        case kwFjEmod   :  {
 271                                const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 272                                if (c ==nullptr) return 0.0;
 273                                double rsum(0.0);
 274                                if (c->getEnd1Curvature().y())
 275                                    rsum += 1.0/c->getEnd1Curvature().y();
 276                                if (c->getEnd2Curvature().y())
 277                                    rsum += 1.0/c->getEnd2Curvature().y();
 278                                if (userArea_) {
 279#ifdef THREED
 280                                    rsum = std::sqrt(userArea_ / dPi);
 281#else
 282                                    rsum = userArea_ / 2.0;
 283#endif
 284                                    rsum += rsum;
 285                                }
 286                                return (fj_kn_ * rsum);
 287                           }
 288        case kwFjKRatio :  return (fj_ks_ == 0.0 ) ? 0.0 : (fj_kn_/fj_ks_);
 289        case kwFjRmul   :   return fj_rmul();
 290        case kwFjRadius :   return rbar();
 291        case kwFjGap0   :   return fj_gap0();
 292        case kwFjTen    :   return fj_ten();
 293        case kwFjCoh    :   return fj_coh();
 294        case kwFjFa     :   return fj_fa();
 295        case kwFjF      :   var.setValue(fj_f()); return var;
 296        case kwFjM      :   var.setValue(fj_m()); return var;
 297        case kwFjState  :   return bmode_[fj_elem()-1];
 298        case kwFjSlip   :   return smode_[fj_elem()-1];
 299        case kwFjMType  :   return getmType();
 300        case kwFjA      :   return a_[fj_elem()-1];
 301        case kwFjEgap   :   return egap_[fj_elem()-1];
 302        case kwFjGap    :   return gap().x();
 303        case kwFjNstr   :   return -f_[fj_elem()-1].x() / a_[fj_elem()-1];
 304        case kwFjSstr   :   return f_[fj_elem()-1].y() / a_[fj_elem()-1];
 305        case kwFjSs     :   return tauC((-f_[fj_elem()-1].x() / a_[fj_elem()-1]),(bmode_[fj_elem()-1]==3));
 306        case kwFjRelBr  :   var.setValue(DVect2(theta(),thetaM())); return var;
 307        case kwFjCen    :   var.setValue(getRelElemPos(con,fj_elem()-1)); return var;
 308#ifdef THREED
 309        case kwFjNa     :   return fj_na();
 310#endif
 311        case kwFjTrack  :   var.setValue(orientProps_ ? true : false); return var;
 312        case kwUserArea :   return userArea_;
 313        case kwFjCohRes :   return fj_cohres();
 314        case kwFjResMode:   return fj_resmode();
 315        }
 316        assert(0);
 317        return QVariant();
 318    }
 319
 320    bool ContactModelFlatJoint::getPropertyGlobal(uint32 i) const {
 321        switch (i) {
 322        case kwFjF:
 323            return false;
 324        }
 325        return true;
 326    }
 327
 328    bool ContactModelFlatJoint::setProperty(uint32 i,const QVariant &v,IContact *c) {
 329        bool ok(true);
 330        switch (i) {
 331        case kwFjNr: {
 332                if (!propsFixed()) {
 333                    int val(v.toInt(&ok));
 334                    if (!ok || val < 1)
 335                        throw Exception("fj_nr must be an integer greater than 0.");
 336                    fj_nr(val);
 337                    if (fj_elem() > fj_n())
 338                        fj_elem(fj_n());
 339                    initVectors();
 340                    setAreaQuantities();
 341                } else
 342                    throw Exception("fj_nr cannot be modified.");
 343                return true;
 344            }
 345
 346        case kwFjElem: {
 347               int val(v.toInt(&ok));
 348               if (!ok || val < 1 || val > fj_n())
 349                   throw Exception("fj_elem must be an integer between 1 and %1.",fj_n());
 350               fj_elem(val);
 351               return false;
 352           }
 353        case kwFjKn: {
 354                double val(v.toDouble(&ok));
 355                if (!ok || val<0.0)
 356                    throw Exception("fj_kn must be a positive double.");
 357                fj_kn(val);
 358                return true;
 359            }
 360        case kwFjKs: {
 361                double val(v.toDouble(&ok));
 362                if (!ok || val<0.0)
 363                    throw Exception("fj_ks must be a positive double.");
 364                fj_ks(val);
 365                return true;
 366            }
 367        case kwFjFric: {
 368                double val(v.toDouble(&ok));
 369                if (!ok || val<0.0)
 370                    throw Exception("fj_fric must be a positive double.");
 371                fj_fric(val);
 372                return false;
 373            }
 374        case kwFjRmul: {
 375                if (!propsFixed()) {
 376                    double val(v.toDouble(&ok));
 377                    if (!ok || val<0.01)
 378                        throw Exception("fj_rmul must be a double greater than or equal to 0.01.");
 379                    fj_rmul(val);
 380                    setAreaQuantities();
 381                    return true;
 382                } else
 383                    throw Exception("fj_rmul cannot be modified.");
 384
 385                return false;
 386            }
 387        case kwFjGap0: {
 388                if (!propsFixed()) {
 389                    double val(v.toDouble(&ok));
 390                    if (!ok || val<0.0)
 391                        throw Exception("fj_gap0 must be a positive double.");
 392                    fj_gap0(val);
 393                    if (fj_gap0() > 0.0) {
 394                        for(int i=1; i<=fj_n(); ++i)
 395                            bondElem(i,false);
 396                        // surfaces are parallel w/ gap G
 397                        DVect temp(0.0);
 398                        temp.rx() = fj_gap0();
 399                        gap(temp);
 400                        theta(0.0);
 401                    }
 402                } else
 403                    throw Exception("fj_gap0 cannot be modified.");
 404                return true;
 405            }
 406        case kwFjTen: {
 407                double val(v.toDouble(&ok));
 408                if (!ok || val<0.0)
 409                    throw Exception("fj_ten must be a positive double.");
 410                fj_ten(val);
 411                return false;
 412            }
 413        case kwFjFa: {
 414                double val(v.toDouble(&ok));
 415                if (!ok || val<0.0)
 416                    throw Exception("fj_fa must be a positive double.");
 417                fj_fa(val);
 418                return false;
 419            }
 420        case kwFjCoh: {
 421                double val(v.toDouble(&ok));
 422                if (!ok || val<0.0)
 423                    throw Exception("fj_coh must be a positive double.");
 424                fj_coh(val);
 425                return false;
 426            }
 427        case kwFjA: {
 428                double val(v.toDouble(&ok));
 429                if (!ok || val<0.0)
 430                    throw Exception("fj_area must be a positive double.");
 431                a_[fj_elem()-1] = val;
 432                return false;
 433            }
 434        case kwFjNstr: {
 435                double val(v.toDouble(&ok));
 436                if (!ok || val<0.0)
 437                    throw Exception("fj_sigma must be a positive double.");
 438                f_[fj_elem()-1].rx() = -val * a_[fj_elem()-1];
 439                return false;
 440            }
 441        case kwFjSstr: {
 442                double val(v.toDouble(&ok));
 443                if (!ok || val<0.0)
 444                    throw Exception("fj_tau must be a positive double.");
 445                f_[fj_elem()-1].ry() = val * a_[fj_elem()-1];
 446                return false;
 447            }
 448#ifdef THREED
 449        case kwFjNa: {
 450                if (!propsFixed()) {
 451                    int val(v.toInt(&ok));
 452                    if (!ok || val < 1)
 453                        throw Exception("fj_na must be an integer greater than 0.");
 454                    fj_na(val);
 455                    if (fj_elem() > fj_n())
 456                        fj_elem(fj_n());
 457                    initVectors();
 458                    setAreaQuantities();
 459                } else
 460                    throw Exception("fj_na cannot be modified.");
 461                return true;
 462            }
 463#endif
 464        case kwFjCen: {
 465                if (!v.canConvert<DVect>())
 466                    throw Exception("fj_cen cannot be modified.");
 467                DVect val(v.value<DVect>());
 468                int el = fj_elem()-1;
 469                setRelElemPos(c,el,val);
 470                return false;
 471            }
 472        case kwFjTrack: {
 473                if (!v.canConvert<bool>())
 474                    throw Exception("fj_track must be a boolean.");
 475                bool b = v.toBool();
 476                if (b) {
 477                    if (!orientProps_)
 478                        orientProps_ = NEWC(orientProps());
 479                } else {
 480                    if (orientProps_) {
 481                        delete orientProps_;
 482                        orientProps_ = 0;
 483                    }
 484                }
 485                return true;
 486            }
 487        case kwUserArea: {
 488                if (!v.canConvert<double>())
 489                    throw Exception("user_area must be a double.");
 490                double val(v.toDouble());
 491                if (val < 0.0)
 492                    throw Exception("Negative user_area not allowed.");
 493                userArea_ = val;
 494                propsFixed_ = false;
 495                return true;
 496            }
 497        case kwFjCohRes: {
 498                double val(v.toDouble(&ok));
 499                if (!ok || val<0.0)
 500                    throw Exception("fj_cohres must be a positive double.");
 501                fj_cohres(val);
 502                return false;
 503            }
 504        case kwFjResMode: {
 505                int val(v.toInt(&ok));
 506                if (!ok || (val != 0 && val != 1))
 507                    throw Exception("fj_resmode must be 0 or 1.");
 508                fj_resmode(val);
 509                return false;
 510            }
 511        }
 512        return false;
 513    }
 514
 515    bool ContactModelFlatJoint::getPropertyReadOnly(uint32 i) const {
 516        switch (i) {
 517        case kwFjF:
 518        case kwFjM:
 519        case kwFjGap:
 520        case kwFjRelBr:
 521        case kwFjState:
 522        case kwFjSlip:
 523        case kwFjEgap:
 524        case kwFjNstr:
 525        case kwFjSstr:
 526        case kwFjSs:
 527        case kwFjRadius:
 528            return true;
 529        default:
 530            break;
 531        }
 532        return false;
 533    }
 534
 535    QString  ContactModelFlatJoint::getMethodArguments(uint32 i) const {
 536        switch (i) {
 537        case kwBond:
 538        case kwUnbond:
 539            return "gap,element";
 540        case KwDeformability:
 541            return "emod,kratio";
 542        case kwInitialize:
 543            return "force,moment";
 544        }
 545        return QString();
 546    }
 547
 548    bool ContactModelFlatJoint::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
 549        IContactMechanical *c(convert_getcast<IContactMechanical>(con));
 550        bool bond(false);
 551        switch (i) {
 552        case kwBond:
 553            bond = true;
 554            [[fallthrough]];
 555        case kwUnbond: {
 556                int seg(0);
 557                double mingap = -1.0 * limits<double>::max();
 558                double maxgap = 0;
 559                if (vl.size()==2) {
 560                    // The first is the gap
 561                    QVariant arg = vl.at(0);
 562                    if (!arg.isNull()) {
 563                        if (arg.canConvert<double>())
 564                            maxgap = vl.at(0).toDouble();
 565                        else if (arg.canConvert<DVect2>()) {
 566                            DVect2 value = vl.at(0).value<DVect2>();
 567                            mingap = value.minComp();
 568                            maxgap = value.maxComp();
 569                        } else
 570                            throw Exception("Argument %1 not recognized in method %2 in contact model %3.",vl.at(0),bond ? "bond":"unbond",getName());
 571                    }
 572                    arg = vl.at(1);
 573                    if (!arg.isNull()) {
 574                        seg = vl.at(1).toUInt();
 575                        if (seg < 1)
 576                            throw Exception("Element indices start at 1 in method %1 in contact model %2.",bond ? "bond":"unbond",getName());
 577                        if (seg > fj_n())
 578                            throw Exception("Element index %1 exceeds segments number (%2) in method %3 in contact model %4.",seg,fj_n(),bond ? "bond":"unbond",getName());
 579                    }
 580                }
 581                double gap = c->getGap();
 582                if (gap >= mingap && gap <= maxgap) {
 583                    if (!seg) {
 584                        for(int iSeg=1; iSeg<=fj_n(); ++iSeg)
 585                            bondElem(iSeg,bond);
 586                    } else {
 587                        bondElem(seg,bond);
 588                    }
 589                    // If have installed bonds and tracking is enabled then set the contact normal appropriately
 590                    if (orientProps_) {
 591                        orientProps_->orient1_ = Quat::identity();
 592                        orientProps_->orient2_ = Quat::identity();
 593                        orientProps_->origNormal_ = toVect(con->getNormal());
 594                    }
 595                }
 596                return true;
 597             }
 598        case KwDeformability:
 599            {
 600                double emod;
 601                double krat;
 602                if (vl.at(0).isNull())
 603                    throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
 604                emod = vl.at(0).toDouble();
 605                if (emod<0.0)
 606                    throw Exception("Negative emod not allowed in contact model %1.",getName());
 607                if (vl.at(1).isNull())
 608                    throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
 609                krat = vl.at(1).toDouble();
 610                if (krat<0.0)
 611                    throw Exception("Negative stiffness ratio not allowed in contact model %1.",getName());
 612                double rsum(0.0);
 613                if (c->getEnd1Curvature().y())
 614                    rsum += 1.0/c->getEnd1Curvature().y();
 615                if (c->getEnd2Curvature().y())
 616                    rsum += 1.0/c->getEnd2Curvature().y();
 617                if (userArea_) {
 618#ifdef THREED
 619                    rsum = std::sqrt(userArea_ / dPi);
 620#else
 621                    rsum = userArea_ / 2.0;
 622#endif
 623                    rsum += rsum;
 624                }
 625                fj_kn_ = emod / rsum;
 626                fj_ks_ = (krat == 0.0) ? 0.0 : fj_kn_ / krat;
 627                return true;
 628            }
 629        case KwUpdateGeom: {
 630                // go through and update the total area (atot) and the
 631                // radius rbar
 632                double at = 0.0;
 633                for (int i=1; i<=fj_n(); ++i)
 634                    at += a_[i-1];
 635                atot(at);
 636                //get the equivalent radius
 637#ifdef THREED
 638                rbar(sqrt(at/dPi));
 639#else
 640                rbar(at/2.0);
 641#endif
 642                return true;
 643            }
 644        case kwArea: {
 645                if (!userArea_) {
 646                    double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
 647#ifdef THREED
 648                    userArea_ = rsq * rsq * dPi;
 649#else
 650                    userArea_ = rsq * 2.0;
 651#endif
 652                }
 653                return true;
 654            }
 655        case kwInitialize: {
 656                DVect force;
 657                DAVect moment;
 658                if (vl.at(0).isNull())
 659                    throw Exception("Argument force must be specified with method initialize in contact model %1.",getName());
 660                force = vl.at(0).value<DVect>();
 661                if (vl.at(1).isNull())
 662                    throw Exception("Argument moment must be specified with method initialize in contact model %1.",getName());
 663#ifdef THREED
 664                moment = vl.at(1).value<DVect>();
 665#else
 666                moment.rz() = vl.at(1).toDouble();
 667#endif
 668                // Set the gap accordingly to get the correct force
 669                setForce(force,con);
 670                fj_m_set(moment);
 671                return true;
 672            }
 673        }
 674        return false;
 675    }
 676
 677    double ContactModelFlatJoint::getEnergy(uint32 i) const {
 678        double ret(0.0);
 679        if (!energies_)
 680            return ret;
 681        switch (i) {
 682        case kwEStrain:  return energies_->estrain_;
 683        case kwESlip:    return energies_->eslip_;
 684        }
 685        assert(0);
 686        return ret;
 687    }
 688
 689    bool ContactModelFlatJoint::getEnergyAccumulate(uint32 i) const {
 690        switch (i) {
 691        case kwEStrain:  return false;
 692        case kwESlip:    return true;
 693        }
 694        assert(0);
 695        return false;
 696    }
 697
 698    void ContactModelFlatJoint::setEnergy(uint32 i,const double &d) {
 699        if (!energies_) return;
 700        switch (i) {
 701        case kwEStrain:  energies_->estrain_ = d; return;
 702        case kwESlip:    energies_->eslip_   = d; return;
 703        }
 704        assert(0);
 705        return;
 706    }
 707
 708    bool ContactModelFlatJoint::validate(ContactModelMechanicalState *state,const double &) {
 709        assert(state);
 710        const IContactMechanical *c = state->getMechanicalContact();
 711        assert(c);
 712        // This presumes that one of the ends has a non-zero curvature
 713        rmin(1.0/std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
 714        if (userArea_) {
 715#ifdef THREED
 716            rmin(std::sqrt(userArea_ / dPi));
 717#else
 718            rmin(userArea_ / 2.0);
 719#endif
 720        }
 721        if (!propsFixed()) {
 722            setAreaQuantities();
 723            mType(getmType());
 724        }
 725
 726        // Initialize the tracking if not initialized
 727        if (orientProps_ && orientProps_->origNormal_ == DVect(0.0)) {
 728            orientProps_->origNormal_ = toVect(c->getContact()->getNormal());
 729            orientProps_->orient1_ = Quat::identity();
 730            orientProps_->orient2_ = Quat::identity();
 731        }
 732
 733        if (state->trackEnergy_)
 734            activateEnergy();
 735
 736        updateEffectiveStiffness(state);
 737        return checkActivity(state->gap_);
 738    }
 739
 740    void ContactModelFlatJoint::updateEffectiveStiffness(ContactModelMechanicalState *) {
 741        DVect2 ret(fj_kn_,fj_ks_);
 742        ret *= atot();
 743        effectiveTranslationalStiffness(ret);
 744#ifdef TWOD
 745        effectiveRotationalStiffness(DAVect(fj_kn() * (2.0/3.0)*rbar()*rbar()*rbar()));
 746#else
 747        double piR4 = dPi * rbar() * rbar() * rbar() * rbar();
 748        double t = fj_kn() * 0.25 * piR4;
 749        effectiveRotationalStiffness(DAVect(fj_ks() * 0.50 * piR4,t,t));
 750#endif
 751    }
 752
 753    bool ContactModelFlatJoint::forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) {
 754        if (!propsFixed())
 755            propsFixed(true);
 756        timestep;
 757        assert(state);
 758
 759        if (state->activated()) {
 760            if (cmEvents_[fActivated] >= 0) {
 761                auto c = state->getContact();
 762                std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
 763                IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
 764                fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
 765            }
 766        }
 767
 768        // Update the orientations
 769        if (orientProps_) {
 770            orientProps_->orient1_.increment(state->getMechanicalContact()->getEnd1Mechanical()->getAngVelocity()*timestep);
 771            orientProps_->orient2_.increment(state->getMechanicalContact()->getEnd2Mechanical()->getAngVelocity()*timestep);
 772        }
 773
 774#ifdef TWOD
 775        // Translational increment in local coordinates
 776        DVect del_U = state->relativeTranslationalIncrement_;
 777        double del_theta  = state->relativeAngularIncrement_.z();
 778        gap(gap() + del_U); // in normal and shear direction in 2D
 779        theta(theta() + del_theta);
 780        double dSig=0.0, dTau=0.0;
 781        double delX = 2*rbar() / fj_n();
 782        double rbar2 = rbar() / fj_n();
 783        DVect dFSum(0.0);
 784        double dMSum = 0.0;
 785        if (state->trackEnergy_) {
 786            assert(energies_);
 787            energies_->estrain_ =  0.0;
 788        }
 789        bool oneBonded = false;
 790        for(int i=0; i<fj_n(); ++i) {
 791            double g0 = gap((i  )*delX);
 792            double g1 = gap((i+1)*delX);
 793            double gMid = 0.5*(g0 + g1);
 794            if (bmode_[i] != 3 && gMid > 0) {
 795                egap_[i] = gMid;
 796                f_[i] = DVect(0.0);
 797                continue;
 798            }
 799            dSig = computeSig(g0,g1,rbar2,a_[i],(bmode_[i]==3));
 800            bool tensileBreak = false;
 801            if (bmode_[i]==3) {
 802                if (state->canFail_ && dSig >= fj_ten()) {
 803                    breakBond(i+1,1,state);
 804                    dSig = dTau = 0.0;
 805                    tensileBreak = true;
 806                }
 807            }
 808            if (!tensileBreak) {
 809                dTau = f_[i].y() / a_[i];
 810                double dUse = delUse(egap_[i],gMid,(bmode_[i]==3),del_U.y());
 811                double dtauP = dTau - fj_ks()*dUse;
 812                double dtauPabs = abs(dtauP);
 813                if (bmode_[i]==3) { // bonded
 814                    if (dtauPabs < tauC(dSig,true))
 815                        dTau = dtauP;
 816                    else {
 817                        if (state->canFail_) {
 818                            breakBond(i+1,2,state);
 819                            if (fj_resmode() == 0)
 820                                dSig = dTau = 0.0;
 821                            else
 822                                dTau = fj_cohres() - dSig * fj_fric();
 823                        }
 824                    }
 825                } else {             // unbonded
 826                    double dtauC = tauC(dSig,false);
 827                    if (dtauPabs <= dtauC) {
 828                        dTau = dtauP;
 829                        slipChange(i+1,false,state);
 830                    } else {
 831                        dTau = dtauP * ( dtauC / dtauPabs );
 832                        slipChange(i+1,true,state);
 833                        if (state->trackEnergy_) { energies_->eslip_ += dtauC*a_[i]*abs(dUse);}
 834                    }
 835                }
 836            }
 837            oneBonded = true;
 838            egap_[i] = gMid;
 839            f_[i] = DVect(-dSig*a_[i],dTau*a_[i]);
 840            dFSum += f_[i];
 841            double m = computeM(g0,g1,rbar2,(bmode_[i]==3)) + fj_m_set().z()/fj_n();
 842            dMSum  += m - rBarl_[i]*f_[i].x();
 843            if (state->trackEnergy_) {
 844                if (fj_kn_) {
 845                    double ie = 2.0*rBarl_[i]*rBarl_[i]*rBarl_[i] / 3.0;
 846                    energies_->estrain_ += 0.5*(dSig*dSig*a_[i] + m*m/ie) / fj_kn_;
 847                }
 848                if (fj_ks_) {
 849                    energies_->estrain_ += 0.5 * dTau*dTau*a_[i] / fj_ks_;
 850                }
 851            }
 852        }
 853        //
 854        fj_f(dFSum);
 855        fj_m(DAVect(dMSum));
 856        if (!oneBonded)
 857            fj_m_set(DAVect(0.0));
 858#else
 859        CAxes localSys = state->getMechanicalContact()->getContact()->getLocalSystem();
 860        DVect trans = state->relativeTranslationalIncrement_; // translation increment in local coordinates
 861        DAVect ang = state->relativeAngularIncrement_; // rotational increment in local coordinates
 862        DVect shear(0.0,trans.y(),trans.z());
 863        DVect del_Us = localSys.toGlobal(shear); // In global coordinates
 864        // What is the twist in global coordinates?
 865        DVect del_Theta_t = localSys.e1()*ang.x();
 866        theta_ += ang.y();
 867        thetaM_ += ang.z();
 868
 869        gap(gap() + trans);
 870        if (state->trackEnergy_) {
 871            assert(energies_);
 872            energies_->estrain_ =  0.0;
 873        }
 874        DVect force(0.0);
 875        DAVect mom(0.0);
 876        bool oneBonded = false;
 877        for (int e=1,i=0; e<=fj_n(); ++e, ++i) {
 878            double gBar1 = gap( rBarl_[i].x(),rBarl_[i].y());
 879            if (!Bonded(e) && gBar1 > 0) {
 880                egap_[i] = gBar1;
 881                f_[i] = DVect(0.0);
 882                continue;
 883            }
 884            DVect r = localSys.e2()*rBarl_[i].x() + localSys.e3()*rBarl_[i].y(); // location of element point
 885            double sigBar_e = sigBar(e);
 886            f_[i].rx() = -sigBar_e * a_[i]; // Step 1...
 887            if (Bonded(e) && (sigBar_e >= fj_ten())) { // break bond in tension
 888                if (state->canFail_) {
 889                    breakBond(e,1,state);
 890                    f_[i] = DVect(0.0);
 891                }
 892            } else {
 893                DVect del_us  = del_Us + (del_Theta_t & r); // In global - has the twist in there
 894                double  del_usl = delUse(egap_[i],gBar1,Bonded(e),(del_us | localSys.e2()));
 895                double  del_usm = delUse(egap_[i],gBar1,Bonded(e),(del_us | localSys.e3()));
 896                double Fs_lP = f_[i].y() - fj_ks() * a_[i] * del_usl;
 897                double Fs_mP = f_[i].z() - fj_ks() * a_[i] * del_usm;
 898                double FsPMag = sqrt( Fs_lP*Fs_lP + Fs_mP*Fs_mP );
 899                double tauBarP = FsPMag / a_[i];
 900                if ( !Bonded(e) ) {
 901                    double tau_c = sigBar_e < 0.0 ? fj_cohres()-fj_fric()*sigBar_e : 0.0;
 902                    if ( tauBarP <= tau_c ) {
 903                        f_[i].ry() = Fs_lP;
 904                        f_[i].rz() = Fs_mP;
 905                        slipChange(e,false,state);
 906                    } else { // enforce sliding
 907                        double sFac = tau_c * a_[i] / FsPMag;
 908                        f_[i].ry() = Fs_lP * sFac;
 909                        f_[i].rz() = Fs_mP * sFac;
 910                        slipChange(e,true,state);
 911                        if (state->trackEnergy_) { energies_->eslip_ += tau_c*a_[i]*sqrt(del_usl*del_usl+del_usm*del_usm);}
 912                    }
 913                } else { // Bonded(e)
 914                    double tau_c = fj_coh() - sigBar_e * tan(dDegrad*fj_fa());
 915                    if ( tauBarP <= tau_c ) {
 916                        f_[i].ry() = Fs_lP;
 917                        f_[i].rz() = Fs_mP;
 918                    } else { // break bond in shear
 919                        if (state->canFail_) {
 920                            breakBond(e,2,state);
 921                            if (fj_resmode() == 0)
 922                                f_[i] = DVect(0.0);
 923                            else {
 924                                double newForce = fj_cohres() - sigBar_e * fj_fric();
 925                                if (!userArea_)
 926                                    newForce *= a_[i];
 927                                else
 928                                    newForce *= userArea_ / fj_n();
 929                                newForce /= std::sqrt(f_[i].y()*f_[i].y() + f_[i].z()*f_[i].z());
 930                                f_[i].ry() *= newForce;
 931                                f_[i].rz() *= newForce;
 932                            }
 933                        }
 934                    }
 935                }
 936            }
 937            oneBonded = true;
 938            force += f_[i];
 939            mom += (localSys.toLocal(r) & f_[i]) + fj_m_set()/fj_n();
 940            egap_[i] = gBar1;
 941            if (state->trackEnergy_) {
 942                energies_->estrain_ += computeStrainEnergy(e);
 943            }
 944        }
 945        fj_f(force);
 946        fj_m(mom);
 947        if (!oneBonded)
 948            fj_m_set(DAVect(0.0));
 949#endif
 950        assert(fj_f_ == fj_f_);
 951        return checkActivity(0.0);
 952    }
 953
 954    bool ContactModelFlatJoint::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
 955        // Account for thermal expansion in incremental mode
 956        if (ts->gapInc_ == 0.0) return false;
 957        DVect dg(0.0);
 958        dg.rx() = ts->gapInc_;
 959        gap(gap() + dg);
 960        return true;
 961    }
 962
 963    bool ContactModelFlatJoint::coupling(ContactModelMechanicalState* , ContactModelThermalState*,ContactModelFluidState* fs, const double&) {
 964        // Account for strain-induced pressure increase in fluid contact
 965        if (!fs)
 966            return false;
 967        double dp = fs->dp_geom_;
 968        DVect force = DVect(0.0);
 969        for (int e = 1, i = 0; e <= fj_n(); ++e, ++i) {
 970            DVect finc(0.0);
 971            finc.rx() = dp * a_[i];
 972            f_[i] -= finc;
 973            force += f_[i];
 974        }
 975        fj_f(force);
 976        return true;
 977    }
 978
 979
 980    void ContactModelFlatJoint::setAreaQuantities() {
 981        rbar(fj_rmul() * rmin());
 982#ifdef TWOD
 983        atot(2.0 * rbar());
 984        double v = atot()/fj_n();
 985        for (int i=1; i<=fj_n(); ++i) {
 986            a_[i-1] = v;
 987            rBarl_[i-1] = rbar() * (double(-2*i + 1 + fj_n()) / fj_n());
 988        }
 989#else
 990        atot(dPi * rbar() * rbar());
 991        double del_r  = rbar() / fj_nr();
 992        double del_al = 2.0*dPi / fj_na();
 993        double fac = 2.0/3.0;
 994        for (int i=0; i < fj_n(); ++i) {
 995            double dVal = i / fj_na();
 996            int I = (int)floor( dVal );
 997            int J = i - I*fj_na();
 998            double r1  =  I      * del_r;
 999            double r2  = (I + 1) * del_r;
1000            double al1 =  J      * del_al;
1001            double al2 = (J + 1) * del_al;
1002            a_[i] = 0.5 * (al2 - al1) * (r2*r2 - r1*r1);
1003            rBarl_[i] = DVect2(((sin(al2) - sin(al1)) / (al2 - al1))*((r2*r2*r2 - r1*r1*r1)/(r2*r2 - r1*r1)),
1004                               ((cos(al1) - cos(al2)) / (al2 - al1))*((r2*r2*r2 - r1*r1*r1)/(r2*r2 - r1*r1)))*fac;
1005        }
1006#endif
1007        updateEffectiveStiffness(0);
1008    }
1009
1010    DVect ContactModelFlatJoint::getRelElemPos(const IContact* c,int i) const {
1011        DVect ret(0.0);
1012        if (c) {
1013            ret = c->getPosition();
1014            CAxes localSys = c->getLocalSystem();
1015#ifdef TWOD
1016            ret += localSys.e2()*rBarl_[i];
1017#else
1018            ret += localSys.e2()*rBarl_[i].x() + localSys.e3()*rBarl_[i].y();
1019#endif
1020        }
1021        return ret;
1022    }
1023
1024    void ContactModelFlatJoint::setRelElemPos(const IContact* c,int i,const DVect &pos) {
1025        // pos is a position in space in global coordinates
1026        propsFixed(true);
1027        if (c) {
1028            // project onto the plane
1029            DVect cp = c->getPosition();
1030            DVect norm = toVect(c->getNormal());
1031            double sd = norm|(cp - pos);
1032            // np is the point on the plane
1033            DVect np = pos+norm*sd;
1034            np = np-cp;
1035            CAxes localSys = c->getLocalSystem();
1036            np = localSys.toLocal(np);
1037#ifdef TWOD
1038            rBarl_[i] = np.y();
1039#else
1040            rBarl_[i] = DVect2(np.y(),np.z());
1041#endif
1042        }
1043    }
1044
1045    int ContactModelFlatJoint::getmType() const {
1046        if (propsFixed()) return mType();
1047        //
1048        if (fj_gap0() > 0.0)   return 2;
1049        //
1050        // If we get to here, then G == 0.0.
1051        bool AllBonded = true;
1052        bool AllSlit = true;
1053        for(int i=0; i<fj_n(); ++i) {
1054            if (bmode_[i] != 3) AllBonded = false;
1055            else AllSlit = false;
1056        }
1057        if (AllBonded) return 1;
1058        if (AllSlit)   return 3;
1059        //
1060        return 4;
1061    }
1062
1063    void ContactModelFlatJoint::bondElem(int iSeg,bool bBond ) {
1064        if (bBond) {
1065            if (fj_gap0() == 0.0) {
1066                bmode_[iSeg-1]  = 3;
1067            } else
1068                bmode_[iSeg-1] = 0;
1069        } else
1070            bmode_[iSeg-1] = 0;
1071    }
1072
1073    void ContactModelFlatJoint::breakBond(int iSeg,int fmode,ContactModelMechanicalState *state) {
1074        bmode_[iSeg-1]  = fmode;
1075        if (cmEvents_[fBondBreak] >= 0) {
1076            auto c = state->getContact();
1077            std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1078                                                 fish::Parameter((qint64)iSeg),
1079                                                 fish::Parameter((qint64)fmode),
1080                                                 fish::Parameter(computeStrainEnergy(iSeg))
1081                                               };
1082            IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1083            fi->setCMFishCallArguments(c,arg,cmEvents_[fBondBreak]);
1084        }
1085        if (!isBonded() && cmEvents_[fBroken] >= 0) {
1086            auto c = state->getContact();
1087            std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
1088            IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1089            fi->setCMFishCallArguments(c,arg,cmEvents_[fBroken]);
1090        }
1091    }
1092
1093    void ContactModelFlatJoint::slipChange(int iSeg,bool smode,ContactModelMechanicalState *state) {
1094        bool emitEvent = false;
1095        if (smode) {
1096            if (!smode_[iSeg-1]) {
1097                emitEvent = true;
1098                smode_[iSeg-1] = smode;
1099            }
1100        } else {
1101            if (smode_[iSeg-1]) {
1102                emitEvent = true;
1103                smode_[iSeg-1] = smode;
1104            }
1105        }
1106        if (emitEvent && cmEvents_[fSlipChange] >= 0) {
1107            auto c = state->getContact();
1108            std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
1109                                                 fish::Parameter((qint64)iSeg),
1110                                                 fish::Parameter(smode) };
1111            IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
1112            fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
1113        }
1114    }
1115
1116    double ContactModelFlatJoint::tauC(const double &dSig,bool bBonded) const {
1117        if (bBonded) return (fj_coh() + (tan(dDegrad*fj_fa()) * (-dSig)) );
1118        else
1119            return (dSig < 0.0 ? fj_cohres() - fj_fric() * dSig : 0.0 );
1120    }
1121
1122#ifdef THREED
1123    double ContactModelFlatJoint::xi(int comp) const {
1124        if (comp == 1) return abs(theta_) <= 1e-12 ? 0.0 : theta_/thbMag();
1125        else           return abs(thetaM_) <= 1e-12 ? 0.0 : thetaM_/thbMag();
1126    }
1127#endif
1128
1129    void ContactModelFlatJoint::initVectors() {
1130        a_.resize(fj_n());
1131        rBarl_.resize(fj_n());
1132        bmode_.resize(fj_n());
1133        smode_.resize(fj_n());
1134        egap_.resize(fj_n());
1135        f_.resize(fj_n());
1136        for (int i=0; i<fj_n(); ++i) {
1137            a_[i] = egap_[i] = 0.0;
1138#ifdef THREED
1139            rBarl_[i] = DVect2(0.0);
1140#else
1141            rBarl_[i] = 0.0;
1142#endif
1143            f_[i] = DVect(0.0);
1144            bmode_[i] = 0;
1145            smode_[i] = false;
1146        }
1147    }
1148
1149#ifdef TWOD
1150    double ContactModelFlatJoint::gap(const double &x) const {
1151        return gap().x() + theta()*(x - rbar());
1152    }
1153#else
1154    double ContactModelFlatJoint::gap(const double &r_l,const double &r_m ) const {
1155       return gap().x() + ( r_m*xi(1) - r_l*xi(2) ) * thbMag();
1156    }
1157
1158    double ContactModelFlatJoint::sigBar(int e) const {
1159        if (!Bonded(e)&& gap(rBarl_[e-1].x(),rBarl_[e-1].y()) >= 0.0)
1160            return 0.0;
1161        else
1162            return fj_kn() * gap(rBarl_[e-1].x(),rBarl_[e-1].y());
1163    }
1164
1165    double ContactModelFlatJoint::tauBar(int e) const {
1166        return a_[e-1] <= 1e-12 ?
1167        0.0 : sqrt(f_[e-1].y()*f_[e-1].y() + f_[e-1].z()*f_[e-1].z())/a_[e-1] ;
1168    }
1169
1170#endif
1171
1172    double ContactModelFlatJoint::computeStrainEnergy(int e) const {
1173        double ret(0.0);
1174        int i = e - 1;
1175#ifdef TWOD
1176        double delX = 2 * rbar() / fj_n();
1177        double g0 = gap((i)*delX);
1178        double g1 = gap((i + 1)*delX);
1179        double rbar2 = rbar() / fj_n();
1180        double dSig = computeSig(g0, g1, rbar2, a_[i], (bmode_[i] == 3));
1181        double m = computeM(g0, g1, rbar2, (bmode_[i] == 3));
1182        double dTau = f_[i].y() / a_[i]; // only valid before failure
1183        if (fj_kn_) {
1184            double ie = 2.0*rBarl_[i] * rBarl_[i] * rBarl_[i] / 3.0;
1185            ret += 0.5*(dSig*dSig*a_[i] + m * m / ie) / fj_kn_;
1186        }
1187        if (fj_ks_) {
1188            ret += 0.5 * dTau*dTau*a_[i] / fj_ks_;
1189        }
1190#else
1191        if (fj_kn_) {
1192            ret += 0.5*(sigBar(e)*sigBar(e)*a_[i]) / fj_kn_;
1193        }
1194        if (fj_ks_) {
1195            ret += 0.5 * (f_[i].y()*f_[i].y() + f_[i].z()*f_[i].z()) / (fj_ks_*a_[i]);
1196        }
1197#endif
1198        return ret;
1199    }
1200
1201    double ContactModelFlatJoint::computeSig(const double &g0,const double &g1,const double &rbar,
1202                                             const double &dA,bool bBonded ) const {
1203        double gTerm;
1204        switch (getCase(g0, g1)) {
1205        case 1:
1206            if (bBonded)       gTerm =  (g0 + g1);
1207            else if (g0 < 0.0) gTerm = -( g0*g0 / (g1 - g0) );
1208            else               gTerm =  ( g1*g1 / (g1 - g0) );
1209            break;
1210        case 2:
1211            if (bBonded) gTerm = (g0 + g1);
1212            else         gTerm = 0.0;
1213            break;
1214        case 3:
1215            gTerm = (g0 + g1);
1216            break;
1217        default:  gTerm = 0.0;  break;
1218        }
1219        return (fj_kn() * gTerm * rbar) / dA;
1220    }
1221
1222    double ContactModelFlatJoint::computeM(const double &g0,const double &g1,const double &rbar,
1223                                           bool bBonded) const {
1224        double gTerm;
1225        switch (getCase(g0,g1)) {
1226            case 1:
1227                if (bBonded)       gTerm = -((g1 - g0) / 3.0);
1228                else if (g0 < 0.0) gTerm = g0*g0*(g0 - 3.0*g1) / (3.0*(g1-g0)*(g1-g0));
1229                else               gTerm = -(((g1-g0)*(g1-g0)*(g1-g0) + g0*g0*(g0 - 3.0*g1))
1230                                            / (3.0*(g1-g0)*(g1-g0)));
1231            break;
1232          case 2:
1233                if (bBonded) gTerm = -((g1 - g0) / 3.0);
1234                else         gTerm = 0.0;
1235                break;
1236          case 3:
1237                gTerm = -((g1 - g0) / 3.0);
1238                break;
1239          default:  gTerm = 0.0;  break;
1240        }
1241        return fj_kn() * gTerm * rbar*rbar;
1242    }
1243
1244    int ContactModelFlatJoint::getCase(const double &g0,const double &g1) const {
1245        if (g0 * g1 < 0.0) // Case 1: gap changes sign
1246            return 1;
1247        else if (g0 >= 0.0 && g1 >= 0.0) // Case 2: gap remains positive or zero
1248            return 2;
1249        else // Case 3: gap remains negative
1250            return 3;
1251    }
1252
1253    double ContactModelFlatJoint::delUse(const double &gapStart,const double &gapEnd,bool bBonded,
1254                                         const double &delUs) const {
1255        if ( bBonded ) return delUs;
1256        if ( gapStart <= 0.0 ) {
1257            if ( gapEnd <= 0.0 )
1258                return delUs;
1259            else { // gapEnd > 0.0
1260                double xi = -gapStart / (gapEnd - gapStart);
1261                return delUs * xi;
1262            }
1263        } else { // gapStart > 0.0
1264            if ( gapEnd >= 0.0 )
1265                return 0.0;
1266            else { // gapEnd < 0.0
1267                double xi = -gapStart / (gapEnd - gapStart);
1268                return delUs * (1.0 - xi);
1269            }
1270        }
1271    }
1272
1273    bool ContactModelFlatJoint::checkActivity(const double &inGap) {
1274        // If any subcontact is bonded return true
1275        FOR(it,bmode_) if ((*it) == 3)
1276            return true;
1277        // If the normal gap is less than 2*rbar return true
1278        if (gap().x() < 2.0*rbar())
1279            return true;
1280        // check to see if there is overlap (based on the initial gap) to reset activity if the contact has been previously deactivated
1281        if (inGap < 0) {
1282            // reset the relative rotation
1283            theta(0.0);
1284#ifdef THREED
1285            thetaM(0.0);
1286#endif
1287            // set the gap to be the current gap, removing the shear displacement
1288            DVect inp(inGap,0.0);
1289            gap(inp);
1290            return true;
1291        }
1292        return false;
1293    }
1294
1295    void ContactModelFlatJoint::setForce(const DVect &v,IContact *) {
1296        fj_f_ = v;
1297        DVect df = v / f_.size();
1298        for (int i=0; i<f_.size(); ++i)
1299            f_[i] = df;
1300        // Set gap consistent with normal force
1301        double at = userArea_;
1302        if (!userArea_) {
1303            for (int i = 1; i <= fj_n(); ++i)
1304                at += a_[i - 1];
1305        }
1306        gap_.rx() = -1.0 * v.x() / (fj_kn_ * at);
1307    }
1308
1309    void ContactModelFlatJoint::getSphereList(const IContact *con,std::vector<DVect> *pos,std::vector<double> *rad,std::vector<double> *val) {
1310        assert(pos);
1311        assert(rad);
1312        assert(val);
1313        if (!orientProps_)
1314            return;
1315        // find minimal radii for end1
1316        double br = convert_getcast<IContactMechanical>(con)->getEnd1Curvature().y();
1317        if (br) {
1318            const IPiece *p = con->getEnd1();
1319            FArray<const IContact*> arr;
1320            p->getContactList(&arr);
1321            double maxgap = 0.0;
1322            FOR(ic,arr) {
1323                const IContactMechanical *mc = convert_getcast<IContactMechanical>(*ic);
1324                const IContactModelMechanical *mcm = mc->getModelMechanical();
1325                if (mcm->getContactModel()->getIndex() == ContactModelFlatJoint::getIndex()) {
1326                    const ContactModelFlatJoint *in = dynamic_cast<const ContactModelFlatJoint*>(mcm);
1327                    maxgap = std::max<double>(maxgap,in->gap().x()- mc->getGap());
1328                }
1329            }
1330            br = 1.0 / br - 0.5*maxgap;
1331            const IContactMechanical *mc = convert_getcast<IContactMechanical>(con);
1332            pos->push_back(convert_getcast<IPieceMechanical>(mc->getEnd1())->getPosition());
1333            rad->push_back(br);
1334            val->push_back(mc->getEnd1()->getIThing()->getID());
1335        }
1336
1337        // Give the end2 sphere - bummer
1338        br = convert_getcast<IContactMechanical>(con)->getEnd2Curvature().y();
1339        if (br) {
1340            const IPiece *p = con->getEnd2();
1341            FArray<const IContact*> arr;
1342            p->getContactList(&arr);
1343            double maxgap = 0.0;
1344            FOR(ic,arr) {
1345                const IContactMechanical *mc = convert_getcast<IContactMechanical>(*ic);
1346                const IContactModelMechanical *mcm = mc->getModelMechanical();
1347                if (mcm->getContactModel()->getIndex() == ContactModelFlatJoint::getIndex()) {
1348                    const ContactModelFlatJoint *in = dynamic_cast<const ContactModelFlatJoint*>(mcm);
1349                    maxgap = std::max<double>(maxgap,in->gap().x()- mc->getGap());
1350                }
1351            }
1352            br = 1.0 / br - 0.5*maxgap;
1353            const IContactMechanical *mc = convert_getcast<IContactMechanical>(con);
1354            pos->push_back(convert_getcast<IPieceMechanical>(mc->getEnd2())->getPosition());
1355            rad->push_back(br);
1356            val->push_back(mc->getEnd2()->getIThing()->getID());
1357        }
1358    }
1359
1360#ifdef THREED
1361
1362    void ContactModelFlatJoint::getDiskList(const IContact *con,std::vector<DVect> *pos,std::vector<DVect> *normal,std::vector<double> *radius,std::vector<double> *val) {
1363        assert(pos);
1364        assert(normal);
1365        assert(radius);
1366        assert(val);
1367        if (!orientProps_)
1368            return;
1369        // plot the contact plane right in the middle of the 2 normals
1370        double rad = fj_rmul()*rmin();
1371        DVect axis1 = orientProps_->orient1_.rotate(orientProps_->origNormal_);
1372        DVect axis2 = orientProps_->orient2_.rotate(orientProps_->origNormal_);
1373        DVect norm = ((axis1.unit()+axis2.unit())*0.5).unit();
1374        pos->push_back(con->getPosition());
1375        normal->push_back(norm);
1376        radius->push_back(rad);
1377        const IContactMechanical *mc = convert_getcast<IContactMechanical>(con);
1378        val->push_back(mc->getLocalForce().mag());
1379    }
1380
1381#endif
1382
1383    void ContactModelFlatJoint::getCylinderList(const IContact *con,std::vector<DVect> *bot,std::vector<DVect> *top,std::vector<double> *radlow,std::vector<double> *radhi,std::vector<double> *val) {
1384        assert(bot);
1385        assert(top);
1386        assert(radlow);
1387        assert(radhi);
1388        assert(val);
1389        if (!orientProps_)
1390            return;
1391        const IContactMechanical *mc = convert_getcast<IContactMechanical>(con);
1392        double br = mc->getEnd1Curvature().y(), br2 = mc->getEnd2Curvature().y();
1393        if (userArea_) {
1394#ifdef THREED
1395            br = std::sqrt(userArea_ / dPi);
1396#else
1397            br = userArea_ / 2.0;
1398#endif
1399            br = 1. / br;
1400            br2 = br;
1401        }
1402
1403        double cgap = mc->getGap();
1404        if (br > 0 && br2 > 0) {
1405            br = 1.0 / br;
1406            br2 = 1.0 / br2;
1407            double rad = fj_rmul()*rmin();
1408            DVect bp = convert_getcast<IPieceMechanical>(mc->getEnd1())->getPosition();
1409            DVect axis = orientProps_->orient1_.rotate(orientProps_->origNormal_);
1410            bot->push_back(axis.unit()*(br-0.5*(gap().x()- cgap))+bp);
1411            top->push_back(bp);
1412            radlow->push_back(rad);
1413            radhi->push_back(0.0);
1414            val->push_back(mc->getEnd1()->getIThing()->getID());
1415            bp = convert_getcast<IPieceMechanical>(mc->getEnd2())->getPosition();
1416            axis = orientProps_->orient2_.rotate(orientProps_->origNormal_);
1417            bot->push_back(axis.unit()*(br2-0.5*(gap().x()-cgap))*(-1.0)+bp);
1418            top->push_back(bp);
1419            radlow->push_back(rad);
1420            radhi->push_back(0.0);
1421            val->push_back(mc->getEnd2()->getIThing()->getID());
1422        }
1423    }
1424
1425    DVect ContactModelFlatJoint::getForce(const IContactMechanical *) const {
1426        DVect ret(fj_f_);
1427        return ret;
1428    }
1429
1430    DAVect ContactModelFlatJoint::getMomentOn1(const IContactMechanical *c) const {
1431        DVect force = getForce(c);
1432        DAVect ret(fj_m_);
1433        c->updateResultingTorqueOn1Local(force,&ret);
1434        return ret;
1435    }
1436
1437    DAVect ContactModelFlatJoint::getMomentOn2(const IContactMechanical *c) const {
1438        DVect force = getForce(c);
1439        DAVect ret(fj_m_);
1440        c->updateResultingTorqueOn2Local(force,&ret);
1441        return ret;
1442    }
1443
1444} // namespace itascaxd
1445
1446// EoF

Top