Linear Contact Model Implementation

See this page for the documentation of this contact model.

contactmodellinear.h

  1#pragma once
  2// contactmodellinear.h
  3
  4#include "contactmodel/src/contactmodelmechanical.h"
  5
  6#ifdef LINEAR_LIB
  7#  define LINEAR_EXPORT EXPORT_TAG
  8#elif defined(NO_MODEL_IMPORT)
  9#  define LINEAR_EXPORT
 10#else
 11#  define LINEAR_EXPORT IMPORT_TAG
 12#endif
 13
 14namespace cmodelsxd {
 15    using namespace itasca;
 16
 17    class ContactModelLinear : public ContactModelMechanical {
 18    public:
 19        // Constructor: Set default values for contact model properties.
 20        LINEAR_EXPORT ContactModelLinear();
 21        // Destructor, called when contact is deleted: free allocated memory, etc.
 22        LINEAR_EXPORT virtual ~ContactModelLinear();
 23        // Contact model name (used as keyword for commands and FISH).
 24        QString  getName() const override { return "linear"; }
 25        // The index provides a quick way to determine the type of contact model.
 26        // Each type of contact model in PFC must have a unique index; this is assigned
 27        // by PFC when the contact model is loaded. This index should be set to -1
 28        void     setIndex(int i) override { index_=i;}
 29        int      getIndex() const override {return index_;}
 30        // Contact model version number (e.g., MyModel05_1). The version number can be
 31        // accessed during the save-restore operation (within the archive method,
 32        // testing {stream.getRestoreVersion() == getMinorVersion()} to allow for 
 33        // future modifications to the contact model data structure.
 34        uint32     getMinorVersion() const override;
 35        // Copy the state information to a newly created contact model.
 36        // Provide access to state information, for use by copy method.
 37        void     copy(const ContactModel *c) override;
 38        // Provide save-restore capability for the state information.
 39        void     archive(ArchiveStream &) override;
 40        // Enumerator for the properties.
 41        enum PropertyKeys { 
 42              kwKn=1
 43            , kwKs                            
 44            , kwFric   
 45            , kwLinF
 46            , kwLinS
 47            , kwLinMode
 48            , kwRGap
 49            , kwEmod
 50            , kwKRatio
 51            , kwDpNRatio 
 52            , kwDpSRatio
 53            , kwDpMode 
 54            , kwDpF
 55            , kwUserArea
 56        };
 57        // Contact model property names in a comma separated list. The order corresponds with
 58        // the order of the PropertyKeys enumerator above. One can visualize any of these 
 59        // properties in PFC automatically. 
 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                   ",user_area";
 75        }
 76        // Enumerator for the energies.
 77        enum EnergyKeys { 
 78            kwEStrain=1
 79          , kwESlip
 80          , kwEDashpot
 81        };
 82        // Contact model energy names in a comma separated list. The order corresponds with
 83        // the order of the EnergyKeys enumerator above. 
 84        QString  getEnergies() const override {
 85            return "energy-strain"
 86                    ",energy-slip"
 87                    ",energy-dashpot";
 88        }
 89        // Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
 90        double   getEnergy(uint32 i) const override;
 91        // Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1) 
 92        // returns wther or not the estrain energy is accumulated which is false).
 93        bool     getEnergyAccumulate(uint32 i) const override;
 94        // Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
 95        void     setEnergy(uint32 i,const double &d) override; // Base 1
 96        // Activate the energy. This is only called if the energy tracking is enabled. 
 97        void     activateEnergy() override { if (energies_) return; energies_ = NEWC(Energies());}
 98        // Returns whether or not the energy tracking has been enabled for this contact.
 99        bool     getEnergyActivated() const override {return (energies_ != 0);}
100
101        // Enumerator for contact model related FISH callback events. 
102        enum FishCallEvents {
103            fActivated=0
104            ,fSlipChange
105        };
106        // Contact model FISH callback event names in a comma separated list. The order corresponds with
107        // the order of the FishCallEvents enumerator above. 
108        QString  getFishCallEvents() const override {
109            return 
110                "contact_activated"
111                ",slip_change"; 
112        }
113
114        // Return the specified contact model property.
115        QVariant getProperty(uint32 i,const IContact *) const override;
116        // The return value denotes whether or not the property corresponds to the global
117        // or local coordinate system (TRUE: global system, FALSE: local system). The
118        // local system is the contact-plane system (nst) defined as follows.
119        // If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
120        // expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
121        // and tc are unit vectors in directions of the nst axes.
122        // This is used when rendering contact model properties that are vectors.
123        bool     getPropertyGlobal(uint32 i) const override;
124        // Set the specified contact model property, ensuring that it is of the correct type
125        // and within the correct range --- if not, then throw an exception.
126        // The return value denotes whether or not the update has affected the timestep
127        // computation (by having modified the translational or rotational tangent stiffnesses).
128        // If true is returned, then the timestep will be recomputed.
129        bool     setProperty(uint32 i,const QVariant &v,IContact *) override;
130        // The return value denotes whether or not the property is read-only
131        // (TRUE: read-only, FALSE: read-write).
132        bool     getPropertyReadOnly(uint32 i) const override;
133
134        // The return value denotes whether or not the property is inheritable
135        // (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
136        // the endPropertyUpdated method.
137        bool     supportsInheritance(uint32 i) const override;
138        // Return whether or not inheritance is enabled for the specified property.
139        bool     getInheritance(uint32 i) const override { assert(i<32); uint32 mask = to<uint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
140        // Set the inheritance flag for the specified property.
141        void     setInheritance(uint32 i,bool b) override { assert(i<32); uint32 mask = to<uint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }
142
143        // Enumerator for contact model methods.
144        enum MethodKeys { kwDeformability=1, kwArea};
145        // Contact model methoid names in a comma separated list. The order corresponds with
146        // the order of the MethodKeys enumerator above.  
147        QString  getMethods() const override { return "deformability,area";}
148        // Return a comma seprated list of the contact model method arguments (base 1).
149        QString  getMethodArguments(uint32 i) const override;
150        // Set contact model method arguments (base 1). 
151        // The return value denotes whether or not the update has affected the timestep
152        // computation (by having modified the translational or rotational tangent stiffnesses).
153        // If true is returned, then the timestep will be recomputed.
154        bool     setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con=0) override;
155
156        // Prepare for entry into ForceDispLaw. The validate function is called when:
157        // (1) the contact is created, (2) a property of the contact that returns a true via
158        // the setProperty method has been modified and (3) when a set of cycles is executed
159        // via the {cycle N} command.
160        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
161        bool    validate(ContactModelMechanicalState *state,const double &timestep) override;
162        // The endPropertyUpdated method is called whenever a surface property (with a name
163        // that matches an inheritable contact model property name) of one of the contacting
164        // pieces is modified. This allows the contact model to update its associated
165        // properties. The return value denotes whether or not the update has affected
166        // the time step computation (by having modified the translational or rotational
167        // tangent stiffnesses). If true is returned, then the time step will be recomputed.  
168        bool    endPropertyUpdated(const QString &name,const IContactMechanical *c) override;
169        // The forceDisplacementLaw function is called during each cycle. Given the relative
170        // motion of the two contacting pieces (via
171        //   state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
172        //   state->relativeAngularIncrement_       (Dtt, Dtbs, Dtbt)
173        //     Ddn  : relative normal-displacement increment, Ddn > 0 is opening
174        //     Ddss : relative  shear-displacement increment (s-axis component)
175        //     Ddst : relative  shear-displacement increment (t-axis component)
176        //     Dtt  : relative twist-rotation increment
177        //     Dtbs : relative  bend-rotation increment (s-axis component)
178        //     Dtbt : relative  bend-rotation increment (t-axis component)
179        //       The relative displacement and rotation increments:
180        //         Dd = Ddn*nc + Ddss*sc + Ddst*tc
181        //         Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
182        //       where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
183        //       [see {Table 1: Contact State Variables} in PFC Model Components:
184        //       Contacts and Contact Models: Contact Resolution]
185        // ) and the contact properties, this function must update the contact force and
186        // moment.
187        //   The force_ is acting on piece 2, and is expressed in the local coordinate system
188        //   (defined in getPropertyGlobal) such that the first component positive denotes
189        //   compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
190        //   in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to 
191        //   get the total moment. 
192        // The return value indicates the contact activity status (TRUE: active, FALSE:
193        // inactive) during the next cycle.
194        // Additional information:
195        //   * If state->activated() is true, then the contact has just become active (it was
196        //     inactive during the previous time step).
197        //   * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
198        //     the forceDispLaw handle the case of {state->canFail_ == true}.
199        bool    forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) override;
200        // Perform thermal coupling
201        bool    thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState*, IContactThermal*, const double&) override;
202
203        // The getEffectiveXStiffness functions return the translational and rotational
204        // tangent stiffnesses used to compute a stable time step. When a contact is sliding,
205        // the translational tangent shear stiffness is zero (but this stiffness reduction
206        // is typically ignored when computing a stable time step). If the contact model
207        // includes a dashpot, then the translational stiffnesses must be increased (see
208        // Potyondy (2009)).
209        //   [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
210        //   Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
211        //   December 7, 2009.]
212        DVect2  getEffectiveTranslationalStiffness() const override { return effectiveTranslationalStiffness_; }
213        DAVect  getEffectiveRotationalStiffness() const override { return DAVect(0.0);}
214
215        // Return a new instance of the contact model. This is used in the CMAT
216        // when a new contact is created. 
217        ContactModelLinear *clone() const override { return NEWC(ContactModelLinear()); }
218        // The getActivityDistance function is called by the contact-resolution logic when
219        // the CMAT is modified. Return value is the activity distance used by the
220        // checkActivity function. 
221        double              getActivityDistance() const override {return rgap_;}
222        // The isOKToDelete function is called by the contact-resolution logic when...
223        // Return value indicates whether or not the contact may be deleted.
224        // If TRUE, then the contact may be deleted when it is inactive.
225        // If FALSE, then the contact may not be deleted (under any condition).
226        bool                isOKToDelete() const override { return !isBonded(); }
227        // Zero the forces and moments stored in the contact model. This function is called
228        // when the contact becomes inactive.
229        void                resetForcesAndMoments() override { lin_F(DVect(0.0)); dp_F(DVect(0.0)); if (energies_) energies_->estrain_ = 0.0;}
230        void                setForce(const DVect &v,IContact *c) override;
231        void                setArea(const double &d) override { userArea_ = d; }
232        double              getArea() const override { return userArea_; }
233
234        // The checkActivity function is called by the contact-resolution logic when...
235        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
236        // A contact with the linear model is active if the contact gap is
237        // less than or equal to zero.
238        bool     checkActivity(const double &gap) override { return  gap <= rgap_; }
239
240        // Returns the sliding state (FALSE is returned if not implemented).
241        bool     isSliding() const override { return lin_S_; }
242        // Returns the bonding state (FALSE is returned if not implemented).
243        bool     isBonded() const override { return false; }
244
245        // Both of these methods are called only for contacts with facets where the wall 
246        // resolution scheme is set the full. In such cases one might wish to propagate 
247        // contact state information (e.g., shear force) from one active contact to another. 
248        // See the Faceted Wall section in the documentation. 
249        void     propagateStateInformation(IContactModelMechanical* oldCm,const CAxes &oldSystem=CAxes(),const CAxes &newSystem=CAxes()) override;
250        void     setNonForcePropsFrom(IContactModel *oldCM) override;
251
252        /// Return the total force that the contact model holds.
253        DVect    getForce(const IContactMechanical *) const override;
254
255        /// Return the total moment on 1 that the contact model holds
256        DAVect   getMomentOn1(const IContactMechanical *) const override;
257
258        /// Return the total moment on 1 that the contact model holds
259        DAVect   getMomentOn2(const IContactMechanical *) const override;
260
261        // Methods to get and set properties. 
262        const double & kn() const {return kn_;}
263        void           kn(const double &d) {kn_=d;}
264        const double & ks() const {return ks_;}
265        void           ks(const double &d) {ks_=d;}
266        const double & fric() const {return fric_;}
267        void           fric(const double &d) {fric_=d;}
268        const DVect &  lin_F() const {return lin_F_;}
269        void           lin_F(const DVect &f) { lin_F_=f;}
270        bool           lin_S() const {return lin_S_;}
271        void           lin_S(bool b) { lin_S_=b;}
272        uint32           lin_mode() const {return lin_mode_;}
273        void           lin_mode(uint32 i) { lin_mode_= i;}
274        const double & rgap() const {return rgap_;}
275        void           rgap(const double &d) {rgap_=d;}
276
277        bool     hasDamping() const {return dpProps_ ? true : false;}
278        double   dp_nratio() const {return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0);}
279        void     dp_nratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_nratio_=d;}
280        double   dp_sratio() const {return hasDamping() ? dpProps_->dp_sratio_: 0.0;}
281        void     dp_sratio(const double &d) { if(!hasDamping()) return; dpProps_->dp_sratio_=d;}
282        int      dp_mode() const {return hasDamping() ? dpProps_->dp_mode_: -1;}
283        void     dp_mode(int i) { if(!hasDamping()) return; dpProps_->dp_mode_=i;}
284        DVect    dp_F() const {return hasDamping() ? dpProps_->dp_F_: DVect(0.0);}
285        void     dp_F(const DVect &f) { if(!hasDamping()) return; dpProps_->dp_F_=f;}
286
287        bool    hasEnergies() const {return energies_ ? true:false;}
288        double  estrain() const {return hasEnergies() ? energies_->estrain_: 0.0;}
289        void    estrain(const double &d) { if(!hasEnergies()) return; energies_->estrain_=d;}
290        double  eslip() const {return hasEnergies() ? energies_->eslip_: 0.0;}
291        void    eslip(const double &d) { if(!hasEnergies()) return; energies_->eslip_=d;}
292        double  edashpot() const {return hasEnergies() ? energies_->edashpot_: 0.0;}
293        void    edashpot(const double &d) { if(!hasEnergies()) return; energies_->edashpot_=d;}
294
295        uint32 inheritanceField() const {return inheritanceField_;}
296        void inheritanceField(uint32 i) {inheritanceField_ = i;}
297
298        const DVect2 & effectiveTranslationalStiffness()  const          {return effectiveTranslationalStiffness_;}
299        void           effectiveTranslationalStiffness(const DVect2 &v ) {effectiveTranslationalStiffness_=v;}
300
301    private:
302        // Index - used internally by PFC. Should be set to -1 in the cpp file. 
303        static int index_;
304
305        // Structure to store the energies. 
306        struct Energies {
307            Energies() : estrain_(0.0), eslip_(0.0),edashpot_(0.0) {}
308            double estrain_;  // elastic energy stored in contact 
309            double eslip_;    // work dissipated by friction 
310            double edashpot_; // work dissipated by dashpots
311        };
312
313        // Structure to store dashpot quantities. 
314        struct dpProps {
315            dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)) {}
316            double dp_nratio_;     // normal viscous critical damping ratio
317            double dp_sratio_;     // shear  viscous critical damping ratio
318            int    dp_mode_;      // for viscous mode (0-4) 0 = dashpots, 1 = tensile limit, 2 = shear limit, 3 = limit both
319            DVect  dp_F_;  // Force in the dashpots
320        };
321
322        bool   updateKn(const IContactMechanical *con);
323        bool   updateKs(const IContactMechanical *con);
324        bool   updateFric(const IContactMechanical *con);
325
326        void   updateEffectiveStiffness(ContactModelMechanicalState *state);
327
328        void   setDampCoefficients(const double &mass,double *vcn,double *vcs);
329
330        // Contact model inheritance fields.
331        uint32 inheritanceField_;
332
333        // Effective translational stiffness.
334        DVect2  effectiveTranslationalStiffness_;
335
336        // linear model properties
337        double      kn_;        // Normal stiffness
338        double      ks_;        // Shear stiffness
339        double      fric_;      // Coulomb friction coefficient
340        DVect       lin_F_;     // Force carried in the linear model
341        bool        lin_S_;     // The current slip state
342        uint32        lin_mode_;  // Specifies absolute (0) or incremental (1) calculation mode 
343        double      rgap_;      // Reference gap 
344        dpProps *   dpProps_;   // The viscous properties
345        double      userArea_;   // Area as specified by the user 
346
347        Energies *   energies_; // The energies
348
349         
350    };
351} // namespace cmodelsxd
352// EoF

Top

contactmodellinear.cpp

  1// contactmodellinear.cpp
  2#include "contactmodellinear.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
 10#include "kernel/interface/iprogram.h"
 11#include "module/interface/icontact.h"
 12#include "module/interface/icontactmechanical.h"
 13#include "module/interface/icontactthermal.h"
 14#include "module/interface/ifishcalllist.h"
 15#include "module/interface/ipiece.h"
 16#include "module/interface/ipiecemechanical.h"
 17
 18#ifdef LINEAR_LIB
 19#ifdef _WIN32
 20    int __stdcall DllMain(void *,unsigned, void *) {
 21        return 1;
 22    }
 23#endif
 24
 25    extern "C" EXPORT_TAG const char *getName() {
 26#if DIM==3
 27        return "contactmodelmechanical3dlinear";
 28#else
 29        return "contactmodelmechanical2dlinear";
 30#endif
 31    }
 32
 33    extern "C" EXPORT_TAG unsigned getMajorVersion() {
 34        return MAJOR_VERSION
 35    }
 36
 37    extern "C" EXPORT_TAG unsigned getMinorVersion() {
 38        return MINOR_VERSION;
 39    }
 40
 41    extern "C" EXPORT_TAG void *createInstance() {
 42        cmodelsxd::ContactModelLinear *m = NEWC(cmodelsxd::ContactModelLinear());
 43        return (void *)m;
 44    }
 45#endif 
 46
 47namespace cmodelsxd {
 48    static const uint32 linKnMask      = 0x00002; // Base 1!
 49    static const uint32 linKsMask      = 0x00004;
 50    static const uint32 linFricMask    = 0x00008;
 51
 52    using namespace itasca;
 53
 54    int ContactModelLinear::index_ = -1;
 55    uint32 ContactModelLinear::getMinorVersion() const { return MINOR_VERSION;}
 56
 57    ContactModelLinear::ContactModelLinear() : inheritanceField_(linKnMask|linKsMask|linFricMask) 
 58                                             , effectiveTranslationalStiffness_(DVect2(0.0)) 
 59                                             , kn_(0.0)
 60                                             , ks_(0.0)
 61                                             , fric_(0.0)
 62                                             , lin_F_(DVect(0.0))
 63                                             , lin_S_(false)
 64                                             , lin_mode_(0)
 65                                             , rgap_(0.0)
 66                                             , dpProps_(0)
 67                                             , userArea_(0)
 68                                             , energies_(0) {
 69    }
 70
 71    ContactModelLinear::~ContactModelLinear() {
 72        // Make sure to clean up after yourself!
 73        if (dpProps_)
 74            delete dpProps_;
 75        if (energies_)
 76            delete energies_;
 77    }
 78
 79    void ContactModelLinear::archive(ArchiveStream &stream) {
 80        // The stream allows one to archive the values of the contact model
 81        // so that it can be saved and restored. The minor version can be
 82        // used here to allow for incremental changes to the contact model too. 
 83        stream & kn_;
 84        stream & ks_;
 85        stream & fric_;
 86        stream & lin_F_;
 87        stream & lin_S_;
 88        stream & lin_mode_;
 89 
 90        if (stream.getArchiveState()==ArchiveStream::Save) {
 91            bool b = false;
 92            if (dpProps_) {
 93                b = true;
 94                stream & b;
 95                stream & dpProps_->dp_nratio_; 
 96                stream & dpProps_->dp_sratio_; 
 97                stream & dpProps_->dp_mode_; 
 98                stream & dpProps_->dp_F_; 
 99            }
100            else
101                stream & b;
102
103            b = false;
104            if (energies_) {
105                b = true;
106                stream & b;
107                stream & energies_->estrain_;
108                stream & energies_->eslip_;
109                stream & energies_->edashpot_;
110            }
111            else
112                stream & b;
113        } else {
114            bool b(false);
115            stream & b;
116            if (b) {
117                if (!dpProps_)
118                    dpProps_ = NEWC(dpProps());
119                stream & dpProps_->dp_nratio_; 
120                stream & dpProps_->dp_sratio_; 
121                stream & dpProps_->dp_mode_; 
122                stream & dpProps_->dp_F_; 
123            }
124            stream & b;
125            if (b) {
126                if (!energies_)
127                    energies_ = NEWC(Energies());
128                stream & energies_->estrain_;
129                stream & energies_->eslip_;
130                stream & energies_->edashpot_;
131            }
132        }
133
134        stream & inheritanceField_;
135        stream & effectiveTranslationalStiffness_;
136        
137        if (stream.getArchiveState()==ArchiveStream::Save || stream.getRestoreVersion() == getMinorVersion())
138            stream & rgap_;
139
140        if (stream.getArchiveState() == ArchiveStream::Save || stream.getRestoreVersion() > 2) 
141            stream & userArea_;
142    }
143
144    void ContactModelLinear::copy(const ContactModel *cm) {
145        // Copy all of the contact model properties. Used in the CMAT 
146        // when a new contact is created. 
147        ContactModelMechanical::copy(cm);
148        const ContactModelLinear *in = dynamic_cast<const ContactModelLinear*>(cm);
149        if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
150        kn(in->kn());
151        ks(in->ks());
152        fric(in->fric());
153        lin_F(in->lin_F());
154        lin_S(in->lin_S());
155        lin_mode(in->lin_mode());
156        rgap(in->rgap());
157        if (in->hasDamping()) {
158            if (!dpProps_)
159                dpProps_ = NEWC(dpProps());
160            dp_nratio(in->dp_nratio()); 
161            dp_sratio(in->dp_sratio()); 
162            dp_mode(in->dp_mode()); 
163            dp_F(in->dp_F()); 
164        }
165        if (in->hasEnergies()) {
166            if (!energies_)
167                energies_ = NEWC(Energies());
168            estrain(in->estrain());
169            eslip(in->eslip());
170            edashpot(in->edashpot());
171        }
172        userArea_ = in->userArea_;
173        inheritanceField(in->inheritanceField());
174        effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
175    }
176
177
178    QVariant ContactModelLinear::getProperty(uint32 i,const IContact *con) const {
179        // Return the property. The IContact pointer is provided so that 
180        // more complicated properties, depending on contact characteristics,
181        // can be calcualted. 
182        QVariant var;
183        switch (i) {
184        case kwKn:        return kn_;
185        case kwKs:        return ks_;
186        case kwFric:      return fric_;
187        case kwLinF:      var.setValue(lin_F_); return var;
188        case kwLinS:      return lin_S_;
189        case kwLinMode:   return lin_mode_;
190        case kwRGap:      return rgap_;
191        case kwEmod: {
192                        const IContactMechanical *c(convert_getcast<IContactMechanical>(con));
193                        if (c ==nullptr) return 0.0;
194                        double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
195                        double rsum(0.0);
196                        if (c->getEnd1Curvature().y())
197                            rsum += 1.0/c->getEnd1Curvature().y();
198                        if (c->getEnd2Curvature().y())
199                            rsum += 1.0/c->getEnd2Curvature().y();
200                        if (userArea_) {
201#ifdef THREED
202                            rsq = std::sqrt(userArea_ / dPi);
203#else
204                            rsq = userArea_ / 2.0;
205#endif        
206                            rsum = rsq + rsq;
207                            rsq = 1. / rsq;
208                        }
209#ifdef TWOD               
210                        return (kn_ * rsum * rsq / 2.0);
211#else                   
212                        return (kn_ * rsum * rsq * rsq) / dPi;
213#endif                    
214                    }
215        case kwKRatio:    return (ks_ == 0.0) ? 0.0 : (kn_/ks_);
216        case kwDpNRatio:  return dpProps_ ? dpProps_->dp_nratio_ : 0;
217        case kwDpSRatio:  return dpProps_ ? dpProps_->dp_sratio_ : 0;
218        case kwDpMode:    return dpProps_ ? dpProps_->dp_mode_ : 0;
219        case kwDpF: {
220                dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
221                return var;
222            }
223        case kwUserArea:    return userArea_;
224        }
225        assert(0);
226        return QVariant();
227    }
228
229    bool ContactModelLinear::getPropertyGlobal(uint32 i) const {
230        // Returns whether or not a property is held in the global axis system (TRUE)
231        // or the local system (FALSE). Used by the plotting logic.
232        switch (i) {
233        case kwLinF:   
234        case kwDpF:  
235            return false;
236        }
237        return true;
238    }
239
240    bool ContactModelLinear::setProperty(uint32 i,const QVariant &v,IContact *) {
241        // Set a contact model property. Return value indicates that the timestep
242        // should be recalculated. 
243        dpProps dp;
244        switch (i) {
245        case kwKn: {
246                if (!v.canConvert<double>())
247                    throw Exception("kn must be a double.");
248                double val(v.toDouble());
249                if (val<0.0)
250                    throw Exception("Negative kn not allowed.");
251                kn_ = val;
252                return true;
253            }
254        case kwKs: {
255                if (!v.canConvert<double>())
256                    throw Exception("ks must be a double.");
257                double val(v.toDouble());
258                if (val<0.0)
259                    throw Exception("Negative ks not allowed.");
260                ks_ = val;  
261                return true;
262            }
263        case kwFric: {
264                if (!v.canConvert<double>())
265                    throw Exception("fric must be a double.");
266                double val(v.toDouble());
267                if (val<0.0)
268                    throw Exception("Negative fric not allowed.");
269                fric_ = val;  
270                return false;
271            }
272        case kwLinF: {
273                if (!v.canConvert<DVect>())
274                    throw Exception("lin_force must be a vector.");
275                DVect val(v.value<DVect>());
276                lin_F_ = val;
277                return false;
278            }
279        case kwLinMode: {
280                if (!v.canConvert<uint32>())
281                    throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
282                uint32 val(v.toUInt());
283                if (val >1)
284                    throw Exception("lin_mode must be 0 (absolute) or 1 (incremental).");
285                lin_mode_ = val;
286                return false;
287            }
288        case kwRGap: {
289                if (!v.canConvert<double>())
290                    throw Exception("Reference gap must be a double.");
291                double val(v.toDouble());
292                rgap_ = val;  
293                return false;
294            }
295        case kwDpNRatio: {
296                if (!v.canConvert<double>())
297                    throw Exception("dp_nratio must be a double.");
298                double val(v.toDouble());
299                if (val<0.0)
300                    throw Exception("Negative dp_nratio not allowed.");
301                if (val == 0.0 && !dpProps_)
302                    return false;
303                if (!dpProps_)
304                    dpProps_ = NEWC(dpProps());
305                dpProps_->dp_nratio_ = val; 
306                return true;
307            }
308        case kwDpSRatio: {
309                if (!v.canConvert<double>())
310                    throw Exception("dp_sratio must be a double.");
311                double val(v.toDouble());
312                if (val<0.0)
313                    throw Exception("Negative dp_sratio not allowed.");
314                if (val == 0.0 && !dpProps_)
315                    return false;
316                if (!dpProps_)
317                    dpProps_ = NEWC(dpProps());
318                dpProps_->dp_sratio_ = val;
319                return true;
320            }
321        case kwDpMode: {
322                if (!v.canConvert<int>())
323                    throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
324                int val(v.toInt());
325                if (val == 0 && !dpProps_)
326                    return false;
327                if (val < 0 || val > 3)
328                    throw Exception("The viscous mode dp_mode must be 0, 1, 2, or 3.");
329                if (!dpProps_)
330                    dpProps_ = NEWC(dpProps());
331                dpProps_->dp_mode_ = val;
332                return false;
333            }
334        case kwDpF: {
335                if (!v.canConvert<DVect>())
336                    throw Exception("dp_force must be a vector.");
337                DVect val(v.value<DVect>());
338                if (val.fsum() == 0.0 && !dpProps_)
339                    return false;
340                if (!dpProps_)
341                    dpProps_ = NEWC(dpProps());
342                dpProps_->dp_F_ = val;
343                return false;
344            }
345        case kwUserArea: {
346                if (!v.canConvert<double>())
347                    throw Exception("user_area must be a double.");
348                double val(v.toDouble());
349                if (val < 0.0)
350                    throw Exception("Negative user_area not allowed.");
351                userArea_ = val;
352                return true;
353            }
354        }
355        return false;
356    }
357
358    bool ContactModelLinear::getPropertyReadOnly(uint32 i) const {
359        // Returns TRUE if a property is read only or FALSE otherwise. 
360        switch (i) {
361        case kwDpF:
362        case kwLinS:
363        case kwEmod:
364        case kwKRatio:
365            return true;
366        default:
367            break;
368        }
369        return false;
370    }
371
372    bool ContactModelLinear::supportsInheritance(uint32 i) const {
373        // Returns TRUE if a property supports inheritance or FALSE otherwise. 
374        switch (i) {
375        case kwKn:
376        case kwKs:
377        case kwFric:
378            return true;
379        default:
380            break;
381        }
382        return false;
383    }
384
385    QString  ContactModelLinear::getMethodArguments(uint32 i) const {
386        // Return a list of contact model method argument names. 
387        switch (i) {
388        case kwDeformability:
389            return "emod,kratio";
390        case kwArea:
391            return QString();
392        }
393        assert(0);
394        return QString();
395    }
396
397    bool ContactModelLinear::setMethod(uint32 i,const QVector<QVariant> &vl,IContact *con) {
398        // Apply the specified method. 
399        IContactMechanical *c(convert_getcast<IContactMechanical>(con));
400        switch (i) {
401        case kwDeformability: {
402                double emod;
403                double krat;
404                if (vl.at(0).isNull()) 
405                    throw Exception("Argument emod must be specified with method deformability in contact model %1.",getName());
406                emod = vl.at(0).toDouble();
407                if (emod<0.0)
408                    throw Exception("Negative emod not allowed in contact model %1.",getName());
409                if (vl.at(1).isNull()) 
410                    throw Exception("Argument kratio must be specified with method deformability in contact model %1.",getName());
411                krat = vl.at(1).toDouble();
412                if (krat<0.0)
413                    throw Exception("Negative stiffness ratio not allowed in contact model %1.",getName());
414                double rsq(std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
415                double rsum(0.0);
416                if (c->getEnd1Curvature().y())
417                    rsum += 1.0/c->getEnd1Curvature().y();
418                if (c->getEnd2Curvature().y())
419                    rsum += 1.0/c->getEnd2Curvature().y();
420                if (userArea_) {
421#ifdef THREED
422                    rsq = std::sqrt(userArea_ / dPi);
423#else
424                    rsq = userArea_ / 2.0;
425#endif        
426                    rsum = rsq + rsq;
427                    rsq = 1. / rsq;
428                }
429#ifdef TWOD
430                kn_ = 2.0 * emod / (rsq * rsum);
431#else
432                kn_ = dPi * emod / (rsq * rsq * rsum);
433#endif
434                ks_ = (krat == 0.0) ? 0.0 : kn_ / krat;
435                setInheritance(1,false);
436                setInheritance(2,false);
437                return true;
438            }
439        case kwArea: {
440                if (!userArea_) {
441                    double rsq(1./std::max(c->getEnd1Curvature().y(),c->getEnd2Curvature().y()));
442#ifdef THREED
443                    userArea_ = rsq * rsq * dPi;
444#else
445                    userArea_ = rsq * 2.0;
446#endif                            
447                }
448                return true;
449            }
450        }
451        return false;
452    }
453
454    double ContactModelLinear::getEnergy(uint32 i) const {
455        // Return an energy value. 
456        double ret(0.0);
457        if (!energies_)
458            return ret;
459        switch (i) {
460        case kwEStrain:  return energies_->estrain_;
461        case kwESlip:    return energies_->eslip_;
462        case kwEDashpot: return energies_->edashpot_;
463        }
464        assert(0);
465        return ret;
466    }
467
468    bool ContactModelLinear::getEnergyAccumulate(uint32 i) const {
469        // Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
470        switch (i) {
471        case kwEStrain:  return false;
472        case kwESlip:    return true;
473        case kwEDashpot:    return true;
474        }
475        assert(0);
476        return false;
477    }
478
479    void ContactModelLinear::setEnergy(uint32 i,const double &d) {
480        // Set an energy value. 
481        if (!energies_) return;
482        switch (i) {
483        case kwEStrain:  energies_->estrain_ = d; return;  
484        case kwESlip:    energies_->eslip_   = d; return;
485        case kwEDashpot: energies_->edashpot_= d; return;
486        }
487        assert(0);
488        return;
489    }
490
491    bool ContactModelLinear::validate(ContactModelMechanicalState *state,const double &) {
492        // Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
493        // (1) the contact is created, (2) a property of the contact that returns a true via
494        // the setProperty method has been modified and (3) when a set of cycles is executed
495        // via the {cycle N} command.
496        // Return value indicates contact activity (TRUE: active, FALSE: inactive).
497        assert(state);
498        const IContactMechanical *c = state->getMechanicalContact(); 
499        assert(c);
500
501        if (state->trackEnergy_)
502            activateEnergy();
503
504        if (inheritanceField_ & linKnMask)
505            updateKn(c);
506        if (inheritanceField_ & linKsMask)
507            updateKs(c);
508        if (inheritanceField_ & linFricMask)
509            updateFric(c);
510
511        updateEffectiveStiffness(state);
512        return checkActivity(state->gap_);
513    }
514
515    static const QString knstr("kn");
516    bool ContactModelLinear::updateKn(const IContactMechanical *con) {
517        assert(con);
518        QVariant v1 = con->getEnd1()->getProperty(knstr);
519        QVariant v2 = con->getEnd2()->getProperty(knstr);
520        if (!v1.isValid() || !v2.isValid())
521            return false;
522        double kn1 = v1.toDouble();
523        double kn2 = v2.toDouble();
524        double val = kn_;
525        if (kn1 && kn2)
526            kn_ = kn1*kn2/(kn1+kn2);
527        else if (kn1)
528            kn_ = kn1;
529        else if (kn2)
530            kn_ = kn2;
531        return ( (kn_ != val) );
532    }
533
534    static const QString ksstr("ks");
535    bool ContactModelLinear::updateKs(const IContactMechanical *con) {
536        assert(con);
537        QVariant v1 = con->getEnd1()->getProperty(ksstr);
538        QVariant v2 = con->getEnd2()->getProperty(ksstr);
539        if (!v1.isValid() || !v2.isValid())
540            return false;
541        double ks1 = v1.toDouble();
542        double ks2 = v2.toDouble();
543        double val = ks_;
544        if (ks1 && ks2)
545            ks_ = ks1*ks2/(ks1+ks2);
546        else if (ks1)
547            ks_ = ks1;
548        else if (ks2)
549            ks_ = ks2;
550        return ( (ks_ != val) );
551    }
552
553    static const QString fricstr("fric");
554    bool ContactModelLinear::updateFric(const IContactMechanical *con) {
555        assert(con);
556        QVariant v1 = con->getEnd1()->getProperty(fricstr);
557        QVariant v2 = con->getEnd2()->getProperty(fricstr);
558        if (!v1.isValid() || !v2.isValid())
559            return false;
560        double fric1 = std::max(0.0,v1.toDouble());
561        double fric2 = std::max(0.0,v2.toDouble());
562        double val = fric_;
563        fric_ = std::min(fric1,fric2);
564        return ( (fric_ != val) );
565    }
566
567    bool ContactModelLinear::endPropertyUpdated(const QString &name,const IContactMechanical *c) {
568        // The endPropertyUpdated method is called whenever a surface property (with a name
569        // that matches an inheritable contact model property name) of one of the contacting
570        // pieces is modified. This allows the contact model to update its associated
571        // properties. The return value denotes whether or not the update has affected
572        // the time step computation (by having modified the translational or rotational
573        // tangent stiffnesses). If true is returned, then the time step will be recomputed.  
574        assert(c);
575        QStringList availableProperties = getProperties().simplified().replace(" ","").split(",",Qt::SkipEmptyParts);
576        QRegularExpression rx(name, QRegularExpression::CaseInsensitiveOption);
577        int idx = availableProperties.indexOf(rx)+1;
578        bool ret=false;
579
580        if (idx<=0)
581            return ret;
582         
583        switch(idx) {
584        case kwKn:  { //kn
585                if (inheritanceField_ & linKnMask)
586                    ret = updateKn(c);
587                break;
588            }
589        case kwKs:  { //ks
590                if (inheritanceField_ & linKsMask)
591                    ret =updateKs(c);
592                break;
593            }
594        case kwFric:  { //fric
595                if (inheritanceField_ & linFricMask)
596                    updateFric(c);
597                break;
598            }
599        }
600        return ret;
601    }
602
603    void ContactModelLinear::updateEffectiveStiffness(ContactModelMechanicalState *) {
604        DVect2 ret(kn_,ks_);
605        // correction if viscous damping active
606        if (dpProps_) {
607            DVect2 correct(1.0);
608            if (dpProps_->dp_nratio_)
609                correct.rx() = sqrt(1.0+dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
610            if (dpProps_->dp_sratio_)
611                correct.ry() = sqrt(1.0+dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
612            ret /= (correct*correct);
613        }
614        effectiveTranslationalStiffness_ = ret;
615    }
616     
617    bool ContactModelLinear::forceDisplacementLaw(ContactModelMechanicalState *state,const double &timestep) {
618        assert(state);
619
620        // Current overlap
621        double overlap = rgap_ - state->gap_;
622        // Relative translational increment
623        DVect trans = state->relativeTranslationalIncrement_;
624        // Correction factor to account for when the contact becomes newly active.
625        // We estimate the time of activity during the timestep when the contact has first 
626        // become active and scale the forces accordingly.
627        double correction = 1.0;
628
629        // The contact was just activated from an inactive state
630        if (state->activated()) {
631            // Trigger the FISH callback if one is hooked up to the 
632            // contact_activated event.
633            if (cmEvents_[fActivated] >= 0) {
634                // An FArray of QVariant is returned and these will be passed
635                // to the FISH function as an array of FISH symbols as the second
636                // argument to the FISH callback function. 
637                auto c = state->getContact();
638                std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
639                IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
640                fi->setCMFishCallArguments(c,arg,cmEvents_[fActivated]);
641            }
642            // Calculate the correction factor.
643            if (lin_mode_ == 0 && trans.x()) {
644                correction = -1.0*overlap / trans.x();
645                if (correction < 0)
646                    correction = 1.0;
647            }
648        }
649
650        DVect lin_F_old = lin_F_;
651
652        if (lin_mode_ == 0) 
653            lin_F_.rx() = overlap * kn_; //  absolute mode for normal force calculation
654        else
655            lin_F_.rx() -= correction * trans.x() * kn_; // incremental  mode for normal force calculation
656
657        // Normal force can only be positive.
658        lin_F_.rx() = std::max(0.0,lin_F_.x());
659
660        // Calculate the shear force.
661        DVect sforce(0.0);
662        // dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
663        // Loop over the shear components (note: the 0 component is the normal component)
664        // and calculate the shear force.
665        for (int i=1; i<dim; ++i)
666            sforce.rdof(i) = lin_F_.dof(i) - trans.dof(i) * ks_ * correction;
667
668        // The canFail flag corresponds to whether or not the contact can undergo non-linear
669        // force-displacement response. If the SOLVE ELASTIC command is given then the 
670        // canFail state is set to FALSE. Otherwise it is always TRUE. 
671        if (state->canFail_) {
672            // Resolve sliding. This is the normal force multiplied by the coefficient of friction.
673            double crit = lin_F_.x() * fric_;
674            // The is the magnitude of the shear force.
675            double sfmag = sforce.mag();
676            // Sliding occurs when the magnitude of the shear force is greater than the 
677            // critical value.
678            if (sfmag > crit) {
679                // Lower the shear force to the critical value for sliding.
680                double rat = crit / sfmag;
681                sforce *= rat;
682                // Handle the slip_change event if one has been hooked up. Sliding has commenced.  
683                if (!lin_S_ && cmEvents_[fSlipChange] >= 0) {
684                    auto c = state->getContact();
685                    std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
686                                                         fish::Parameter()                  };
687                    IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
688                    fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
689                }
690                lin_S_ = true;
691            } else {
692                // Handle the slip_change event if one has been hooked up and
693                // the contact was previously sliding. Sliding has ceased.  
694                if (lin_S_) {
695                    if (cmEvents_[fSlipChange] >= 0) {
696                        auto c = state->getContact();
697                        std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
698                                                             fish::Parameter((qint64)1)      };
699                        IFishCallList *fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
700                        fi->setCMFishCallArguments(c,arg,cmEvents_[fSlipChange]);
701                    }
702                    lin_S_ = false;
703                }
704            }
705        }
706        
707        // Set the shear components of the total force.
708        for (int i=1; i<dim; ++i)
709            lin_F_.rdof(i) = sforce.dof(i);
710         
711        // Account for dashpot forces if the dashpot structure has been defined. 
712        if (dpProps_) {
713            dpProps_->dp_F_.fill(0.0);
714            double vcn(0.0), vcs(0.0);
715            // Calculate the damping coefficients. 
716            setDampCoefficients(state->inertialMass_,&vcn,&vcs);
717            // First damp the shear components
718            for (int i=1; i<dim; ++i)
719                dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep; 
720            // Damp the normal component
721            dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;       
722            // Need to change behavior based on the dp_mode.
723            if ((dpProps_->dp_mode_ == 1 || dpProps_->dp_mode_ == 3))  { 
724                // Limit in tension if not bonded.
725                if (dpProps_->dp_F_.x() + lin_F_.x() < 0)
726                    dpProps_->dp_F_.rx() = - lin_F_.rx();
727            }
728            if (lin_S_ && dpProps_->dp_mode_ > 1)  { 
729                // Limit in shear if not sliding.
730                double dfn = dpProps_->dp_F_.rx();
731                dpProps_->dp_F_.fill(0.0); 
732                dpProps_->dp_F_.rx() = dfn; 
733            }
734        }
735
736        //Compute energies if energy tracking has been enabled. 
737        if (state->trackEnergy_) {
738            assert(energies_);
739            energies_->estrain_ =  0.0;
740            if (kn_)
741                // Calculate the strain energy. 
742                energies_->estrain_ = 0.5*lin_F_.x()*lin_F_.x()/kn_;
743            if (ks_) {
744                DVect s = lin_F_;
745                s.rx() = 0.0;
746                double smag2 = s.mag2();
747                // Add the shear component of the strain energy.
748                energies_->estrain_ += 0.5*smag2 / ks_;
749
750                if (lin_S_) {
751                    // If sliding calculate the slip energy and accumulate it.
752                    lin_F_old.rx() = 0.0;
753                    DVect avg_F_s = (s + lin_F_old)*0.5;
754                    DVect u_s_el =  (s - lin_F_old) / ks_;
755                    DVect u_s(0.0);
756                    for (int i=1; i<dim; ++i)
757                        u_s.rdof(i) = trans.dof(i);
758                    energies_->eslip_ -= std::min(0.0,(avg_F_s | (u_s + u_s_el)));
759                }
760            }
761            if (dpProps_) {
762                // Calculate damping energy (accumulated) if the dashpots are active. 
763                energies_->edashpot_ -= dpProps_->dp_F_ | trans;
764            }
765        }
766
767        // This is just a sanity check to ensure, in debug mode, that the force isn't wonky. 
768        assert(lin_F_ == lin_F_);
769        return true;
770    }
771
772    bool ContactModelLinear::thermalCoupling(ContactModelMechanicalState*, ContactModelThermalState* ts, IContactThermal*, const double&) {
773        // Account for thermal expansion in incremental mode
774        if (lin_mode_ == 0 || ts->gapInc_ == 0.0) return false;
775        DVect finc(0.0);
776        finc.rx() = kn_ * ts->gapInc_;
777        lin_F_ -= finc;
778        return true;
779    }
780
781    void ContactModelLinear::setForce(const DVect &v,IContact *c) { 
782        lin_F(v); 
783        if (v.x() > 0) 
784            rgap_ = c->getGap() + v.x() / kn_; 
785    } 
786
787    void ContactModelLinear::propagateStateInformation(IContactModelMechanical* old,const CAxes &oldSystem,const CAxes &newSystem) {
788        // Only called for contacts with wall facets when the wall resolution scheme
789        // is set to full!
790        // Only do something if the contact model is of the same type
791        if (old->getContactModel()->getName().compare("linear",Qt::CaseInsensitive) == 0 && !isBonded()) {
792            ContactModelLinear *oldCm = (ContactModelLinear *)old;
793#ifdef THREED
794            // Need to rotate just the shear component from oldSystem to newSystem
795
796            // Step 1 - rotate oldSystem so that the normal is the same as the normal of newSystem
797            DVect axis = oldSystem.e1() & newSystem.e1();
798            double c, ang, s;
799            DVect re2;
800            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
801                axis = axis.unit();
802                c = oldSystem.e1()|newSystem.e1();
803                if (c > 0)
804                    c = std::min(c,1.0);
805                else
806                    c = std::max(c,-1.0);
807                ang = acos(c);
808                s = sin(ang);
809                double t = 1. - c;
810                DMatrix<3,3> rm;
811                rm.get(0,0) = t*axis.x()*axis.x() + c;
812                rm.get(0,1) = t*axis.x()*axis.y() - axis.z()*s;
813                rm.get(0,2) = t*axis.x()*axis.z() + axis.y()*s;
814                rm.get(1,0) = t*axis.x()*axis.y() + axis.z()*s;
815                rm.get(1,1) = t*axis.y()*axis.y() + c;
816                rm.get(1,2) = t*axis.y()*axis.z() - axis.x()*s;
817                rm.get(2,0) = t*axis.x()*axis.z() - axis.y()*s;
818                rm.get(2,1) = t*axis.y()*axis.z() + axis.x()*s;
819                rm.get(2,2) = t*axis.z()*axis.z() + c;
820                re2 = rm*oldSystem.e2();
821            }
822            else
823                re2 = oldSystem.e2();
824            // Step 2 - get the angle between the oldSystem rotated shear and newSystem shear
825            axis = re2 & newSystem.e2();
826            DVect2 tpf;
827            DMatrix<2,2> m;
828            if (!checktol(axis.abs().maxComp(),0.0,1.0,1000)) {
829                axis = axis.unit();
830                c = re2|newSystem.e2();
831                if (c > 0)
832                    c = std::min(c,1.0);
833                else
834                    c = std::max(c,-1.0);
835                ang = acos(c);
836                if (!checktol(axis.x(),newSystem.e1().x(),1.0,100))
837                    ang *= -1;
838                s = sin(ang);
839                m.get(0,0) = c;
840                m.get(1,0) = s;
841                m.get(0,1) = -m.get(1,0);
842                m.get(1,1) = m.get(0,0);
843                tpf = m*DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
844            } else {
845                m.get(0,0) = 1.;
846                m.get(0,1) = 0.;
847                m.get(1,0) = 0.;
848                m.get(1,1) = 1.;
849                tpf = DVect2(oldCm->lin_F_.y(),oldCm->lin_F_.z());
850            }
851            DVect pforce = DVect(0,tpf.x(),tpf.y());
852#else
853            oldSystem;
854            newSystem;
855            DVect pforce = DVect(0,oldCm->lin_F_.y());
856#endif
857            for (int i=1; i<dim; ++i)
858                lin_F_.rdof(i) += pforce.dof(i);
859            if (lin_mode_ && oldCm->lin_mode_)
860                lin_F_.rx() = oldCm->lin_F_.x();
861            oldCm->lin_F_ = DVect(0.0);
862            if (dpProps_ && oldCm->dpProps_) {
863#ifdef THREED
864                tpf = m*DVect2(oldCm->dpProps_->dp_F_.y(),oldCm->dpProps_->dp_F_.z());
865                pforce = DVect(oldCm->dpProps_->dp_F_.x(),tpf.x(),tpf.y());
866#else
867                pforce = oldCm->dpProps_->dp_F_;
868#endif
869                dpProps_->dp_F_ += pforce;
870                oldCm->dpProps_->dp_F_ = DVect(0.0);
871            }
872            if(oldCm->getEnergyActivated()) {
873                activateEnergy();
874                energies_->estrain_ = oldCm->energies_->estrain_;
875                energies_->edashpot_ = oldCm->energies_->edashpot_;
876                energies_->eslip_ = oldCm->energies_->eslip_;
877                oldCm->energies_->estrain_ = 0.0;
878                oldCm->energies_->edashpot_ = 0.0;
879                oldCm->energies_->eslip_ = 0.0;
880            }
881            rgap_ = oldCm->rgap_;
882        }
883        assert(lin_F_ == lin_F_);
884    }
885
886    void ContactModelLinear::setNonForcePropsFrom(IContactModel *old) {
887        // Only called for contacts with wall facets when the wall resolution scheme
888        // is set to full!
889        // Only do something if the contact model is of the same type
890        if (old->getName().compare("linear",Qt::CaseInsensitive) == 0 && !isBonded()) {
891            ContactModelLinear *oldCm = (ContactModelLinear *)old;
892            kn_ = oldCm->kn_;
893            ks_ = oldCm->ks_;
894            fric_ = oldCm->fric_;
895            lin_mode_ = oldCm->lin_mode_;
896            rgap_ = oldCm->rgap_;
897            userArea_ = oldCm->userArea_;
898
899            if (oldCm->dpProps_) {
900                if (!dpProps_)
901                    dpProps_ = NEWC(dpProps());
902                dpProps_->dp_nratio_ = oldCm->dpProps_->dp_nratio_;
903                dpProps_->dp_sratio_ = oldCm->dpProps_->dp_sratio_;
904                dpProps_->dp_mode_ = oldCm->dpProps_->dp_mode_;
905            }
906        }
907    }
908
909    DVect ContactModelLinear::getForce(const IContactMechanical *) const {
910        DVect ret(lin_F_);
911        if (dpProps_)
912            ret += dpProps_->dp_F_;
913        return ret;
914    }
915
916    DAVect ContactModelLinear::getMomentOn1(const IContactMechanical *c) const {
917        DVect force = getForce(c);
918        DAVect ret(0.0);
919        c->updateResultingTorqueOn1Local(force,&ret);
920        return ret;
921    }
922
923    DAVect ContactModelLinear::getMomentOn2(const IContactMechanical *c) const {
924        DVect force = getForce(c);
925        DAVect ret(0.0);
926        c->updateResultingTorqueOn2Local(force,&ret);
927        return ret;
928    }
929
930    void ContactModelLinear::setDampCoefficients(const double &mass,double *vcn,double *vcs) {
931        *vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(kn_));
932        *vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(ks_));
933    }
934
935} // namespace cmodelsxd
936// EoF

Top