EEPA Model Implementation

See this page for the documentation of this contact model.

• contactmodeleepa.h
• contactmodeleepa.cpp

contactmodeleepa.h

#pragma once
// contactmodeleepa.h

#include "contactmodel/src/contactmodelmechanical.h"

#ifdef EEPA_LIB
#  define EEPA_EXPORT EXPORT_TAG
#elif defined(NO_MODEL_IMPORT)
#  define EEPA_EXPORT
#else
#  define EEPA_EXPORT IMPORT_TAG
#endif

namespace cmodelsxd {
	using namespace itasca;

	class ContactModelEEPA : public ContactModelMechanical {
	public:
		// Constructor: Set default values for contact model properties.
		EEPA_EXPORT ContactModelEEPA();
		// Destructor, called when contact is deleted: free allocated memory, etc.
		EEPA_EXPORT virtual ~ContactModelEEPA();
		// Contact model name (used as keyword for commands and FISH).
		virtual QString  getName() const { return "eepa"; }
		// The index provides a quick way to determine the type of contact model.
		// Each type of contact model in PFC must have a unique index; this is assigned
		// by PFC when the contact model is loaded. This index should be set to -1
		virtual void     setIndex(int i) { index_ = i; }
		virtual int      getIndex() const { return index_; }
		// Contact model version number (e.g., MyModel05_1). The version number can be
		// accessed during the save-restore operation (within the archive method,
		// testing {stream.getRestoreVersion() == getMinorVersion()} to allow for 
		// future modifications to the contact model data structure.
		virtual uint     getMinorVersion() const;
		// Copy the state information to a newly created contact model.
		// Provide access to state information, for use by copy method.
		virtual void     copy(const ContactModel *c);
		// Provide save-restore capability for the state information.
		virtual void     archive(ArchiveStream &);
		// Enumerator for the properties.
		enum PropertyKeys {
			kwShear = 1
			, kwPoiss
			, kwFric
			, kwEepaF
			, kwEepaS
			, kwRGap
			, kwDpNRatio
			, kwDpSRatio
			, kwDpMode
			, kwDpF
			, kwResFric
			, kwResMoment
			, kwResS
			, kwOverlapMax
			, kwPlasRat
			, kwLuExp
			, kwPullOff
			, kwAdhExp
			, kwSurfAdh
			, kwKsFac		
			, kwFmin
		};
		// Contact model property names in a comma separated list. The order corresponds with
		// the order of the PropertyKeys enumerator above. One can visualize any of these 
		// properties in PFC automatically. 
		virtual QString  getProperties() const {
			return "eepa_shear"
				",eepa_poiss"
				",fric"
				",eepa_force"
				",eepa_slip"
				",rgap"
				",dp_nratio"
				",dp_sratio"
				",dp_mode"
				",dp_force"
				",rr_fric"
				",rr_moment"
				",rr_slip"
				",overlap_max"
				",plas_ratio"
				",lu_exp"
				",pull_off"
				",adh_exp"
				",surf_adh"
				",ks_fac"
				",f_min";			
		}
		// Enumerator for the energies.
		enum EnergyKeys {
			kwEStrain = 1
			, kwERRStrain
			, kwESlip
			, kwERRSlip
			, kwEDashpot
		};
		// Contact model energy names in a comma separated list. The order corresponds with
		// the order of the EnergyKeys enumerator above. 
		virtual QString  getEnergies() const {
			return "energy-strain-eepa"
				",energy-rrstrain"
				",energy-slip"
				",energy-rrslip"
				",energy-dashpot";
		}
		// Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
		virtual double   getEnergy(uint i) const;
		// Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1) 
		// returns whether or not the estrain energy is accumulated which is false).
		virtual bool     getEnergyAccumulate(uint i) const;
		// Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
		virtual void     setEnergy(uint i, const double &d); // Base 1
		// Activate the energy. This is only called if the energy tracking is enabled. 
		virtual void     activateEnergy() { if (energies_) return; energies_ = NEWC(Energies()); }
		// Returns whether or not the energy tracking has been enabled for this contact.
		virtual bool     getEnergyActivated() const { return (energies_ != 0); }

		// Enumerator for contact model related FISH callback events. 
		enum FishCallEvents {
			fActivated = 0
			, fSlipChange
		};
		// Contact model FISH callback event names in a comma separated list. The order corresponds with
		// the order of the FishCallEvents enumerator above. 
		virtual QString  getFishCallEvents() const {
			return
				"contact_activated"
				",slip_change";
		}

		// Return the specified contact model property.
		virtual QVariant getProperty(uint i, const IContact *) const;
		// The return value denotes whether or not the property corresponds to the global
		// or local coordinate system (TRUE: global system, FALSE: local system). The
		// local system is the contact-plane system (nst) defined as follows.
		// If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
		// expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
		// and tc are unit vectors in directions of the nst axes.
		// This is used when rendering contact model properties that are vectors.
		virtual bool     getPropertyGlobal(uint i) const;
		// Set the specified contact model property, ensuring that it is of the correct type
		// and within the correct range --- if not, then throw an exception.
		// The return value denotes whether or not the update has affected the timestep
		// computation (by having modified the translational or rotational tangent stiffnesses).
		// If true is returned, then the timestep will be recomputed.
		virtual bool     setProperty(uint i, const QVariant &v, IContact *);
		// The return value denotes whether or not the property is read-only
		// (TRUE: read-only, FALSE: read-write).
		virtual bool     getPropertyReadOnly(uint i) const;

		// The return value denotes whether or not the property is inheritable
		// (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
		// the endPropertyUpdated method.
		virtual bool     supportsInheritance(uint i) const;
		// Return whether or not inheritance is enabled for the specified property.
		virtual bool     getInheritance(uint i) const { assert(i < 32); quint32 mask = to<quint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
		// Set the inheritance flag for the specified property.
		virtual void     setInheritance(uint i, bool b) { assert(i < 32); quint32 mask = to<quint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }

		// Prepare for entry into ForceDispLaw. The validate function is called when:
		// (1) the contact is created, (2) a property of the contact that returns a true via
		// the setProperty method has been modified and (3) when a set of cycles is executed
		// via the {cycle N} command.
		// Return value indicates contact activity (TRUE: active, FALSE: inactive).
		virtual bool    validate(ContactModelMechanicalState *state, const double &timestep);
		// The endPropertyUpdated method is called whenever a surface property (with a name
		// that matches an inheritable contact model property name) of one of the contacting
		// pieces is modified. This allows the contact model to update its associated
		// properties. The return value denotes whether or not the update has affected
		// the time step computation (by having modified the translational or rotational
		// tangent stiffnesses). If true is returned, then the time step will be recomputed.  
		virtual bool    endPropertyUpdated(const QString &name, const IContactMechanical *c);
		// The forceDisplacementLaw function is called during each cycle. Given the relative
		// motion of the two contacting pieces (via
		//   state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
		//   state->relativeAngularIncrement_       (Dtt, Dtbs, Dtbt)
		//     Ddn  : relative normal-displacement increment, Ddn > 0 is opening
		//     Ddss : relative  shear-displacement increment (s-axis component)
		//     Ddst : relative  shear-displacement increment (t-axis component)
		//     Dtt  : relative twist-rotation increment
		//     Dtbs : relative  bend-rotation increment (s-axis component)
		//     Dtbt : relative  bend-rotation increment (t-axis component)
		//       The relative displacement and rotation increments:
		//         Dd = Ddn*nc + Ddss*sc + Ddst*tc
		//         Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
		//       where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
		//       [see {Table 1: Contact State Variables} in PFC Model Components:
		//       Contacts and Contact Models: Contact Resolution]
		// ) and the contact properties, this function must update the contact force and
		// moment.
		//   The force_ is acting on piece 2, and is expressed in the local coordinate system
		//   (defined in getPropertyGlobal) such that the first component positive denotes
		//   compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
		//   in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to 
		// get the total moment. 
		// The return value indicates the contact activity status (TRUE: active, FALSE:
		// inactive) during the next cycle.
		// Additional information:
		//   * If state->activated() is true, then the contact has just become active (it was
		//     inactive during the previous time step).
		//   * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
		//     the forceDispLaw handle the case of {state->canFail_ == true}.
		virtual bool    forceDisplacementLaw(ContactModelMechanicalState *state, const double &timestep);
		// The getEffectiveXStiffness functions return the translational and rotational
		// tangent stiffnesses used to compute a stable time step. When a contact is sliding,
		// the translational tangent shear stiffness is zero (but this stiffness reduction
		// is typically ignored when computing a stable time step). If the contact model
		// includes a dashpot, then the translational stiffnesses must be increased (see
		// Potyondy (2009)).
		//   [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
		//   Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
		//   December 7, 2009.]
		virtual DVect2  getEffectiveTranslationalStiffness() const { return effectiveTranslationalStiffness_; }
		virtual DAVect  getEffectiveRotationalStiffness() const { return effectiveRotationalStiffness_; }

		// Return a new instance of the contact model. This is used in the CMAT
		// when a new contact is created. 
		virtual ContactModelEEPA *clone() const { return NEWC(ContactModelEEPA()); }
		// The getActivityDistance function is called by the contact-resolution logic when
		// the CMAT is modified. Return value is the activity distance used by the
		// checkActivity function.
		virtual double              getActivityDistance() const { return rgap_; }							
		// The isOKToDelete function is called by the contact-resolution logic when...
		// Return value indicates whether or not the contact may be deleted.
		// If TRUE, then the contact may be deleted when it is inactive.
		// If FALSE, then the contact may not be deleted (under any condition).
		virtual bool                isOKToDelete() const { return !isBonded(); }
		// Zero the forces and moments stored in the contact model. This function is called
		// when the contact becomes inactive.
		virtual void                resetForcesAndMoments() {
			eepa_F(DVect(0.0)); dp_F(DVect(0.0));
			res_M(DAVect(0.0));																			
			if (energies_) {
				//energies_->estrain_ = 0.0;
				energies_->errstrain_ = 0.0;
			}
			Moverlap_ = 0.0;																// reset the history dependend variables
			Poverlap_exp_ = 0.0;
			branch_ = 1;
		}
		virtual void     setForce(const DVect &v, IContact *c);
		virtual void     setArea(const double&) { throw Exception("The setArea method cannot be used with the EEPA contact model."); }
		virtual double   getArea() const { return 0.0; }
		// The checkActivity function is called by the contact-resolution logic when...
		// Return value indicates contact activity (TRUE: active, FALSE: inactive).
		// A contact with the arrlinear model is active if the surface gap is less than
		// or equal to the attraction range (a_d0_).
		virtual bool     checkActivity(const double &gap) { return  gap <= rgap_; }							
		
		// Returns the sliding state (FALSE is returned if not implemented).
		virtual bool     isSliding() const { return eepa_S_; }
		// Returns the bonding state (FALSE is returned if not implemented).
		virtual bool     isBonded() const { return false; }

		// Both of these methods are called only for contacts with facets where the wall 
		// resolution scheme is set the full. In such cases one might wish to propagate 
		// contact state information (e.g., shear force) from one active contact to another. 
		// See the Faceted Wall section in the documentation. 

		// Return the total force that the contact model holds.
		virtual DVect    getForce(const IContactMechanical *) const;

		// Return the total moment on 1 that the contact model holds
		virtual DAVect   getMomentOn1(const IContactMechanical *) const;

		// Return the total moment on 1 that the contact model holds
		virtual DAVect   getMomentOn2(const IContactMechanical *) const;

		// Methods to get and set properties. 
		const double & shear() const { return shear_; }
		void           shear(const double &d) { shear_ = d; }
		const double & poiss() const { return poiss_; }
		void           poiss(const double &d) { poiss_ = d; }
		const double & fric() const { return fric_; }
		void           fric(const double &d) { fric_ = d; }
		const DVect &  eepa_F() const { return eepa_F_; }
		void           eepa_F(const DVect &f) { eepa_F_ = f; }
		bool           eepa_S() const { return eepa_S_; }
		void           eepa_S(bool b) { eepa_S_ = b; }
		const double & rgap() const { return rgap_; }
		void           rgap(const double &d) { rgap_ = d; }
		const double & Moverlap() const { return Moverlap_; }
		void           Moverlap(const double &d) { Moverlap_ = d; }
		const double & Poverlap_exp() const { return Poverlap_exp_; }
		void           Poverlap_exp(const double &d) { Poverlap_exp_ = d; }
		const double & plas_ratio() const { return plas_ratio_; }
		void           plas_ratio(const double &d) { plas_ratio_ = d; }
		const double & lu_exp() const { return lu_exp_; }
		void		   lu_exp(const double &d) { lu_exp_ = d; }
		const double & pull_off() const { return pull_off_; }
		void           pull_off(const double &d) { pull_off_ = d; }
		const double & adh_exp() const { return adh_exp_; }
		void	       adh_exp(const double &d) { adh_exp_ = d; }
		const double & surf_adh() const { return surf_adh_; }
		void           surf_adh(const double &d) { surf_adh_ = d; }
		const double & lu_exp_inv() const { return lu_exp_inv_; }
		void           lu_exp_inv(const double &d) { lu_exp_inv_ = d; }
		const double & k1() const { return k1_; }
		void           k1(const double &d) { k1_ = d; }
		const double & k2() const { return k2_; }
		void           k2(const double &d) { k2_ = d; }
		const double & kadh() const { return kadh_; }
		void           kadh(const double &d) { kadh_ = d; }
		const double & ks() const { return ks_; }
		void           ks(const double &d) { ks_ = d; }
		const double & r_hertz() const { return r_hertz_; }
		void           r_hertz(const double &d) { r_hertz_ = d; }
		const double & ks_fac() const { return ks_fac_; }
		void           ks_fac(const double &d) { ks_fac_ = d; }
		const double & f_min() const { return f_min_; }
		void           f_min(const double &d) { f_min_ = d; }
		const int    & branch() const { return branch_; }
		void           branch(const int &d) { branch_ = d; }
		const double & rbar_square() const { return rbar_square_; }
		void           rbar_square(const double &d) { rbar_square_ = d; }

		bool     hasDamping() const { return dpProps_ ? true : false; }
		double   dp_nratio() const { return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0); }
		void     dp_nratio(const double &d) { if (!hasDamping()) return; dpProps_->dp_nratio_ = d; }
		double   dp_sratio() const { return hasDamping() ? dpProps_->dp_sratio_ : 0.0; }
		void     dp_sratio(const double &d) { if (!hasDamping()) return; dpProps_->dp_sratio_ = d; }
		int      dp_mode() const { return hasDamping() ? dpProps_->dp_mode_ : -1; }
		void     dp_mode(int i) { if (!hasDamping()) return; dpProps_->dp_mode_ = i; }
		DVect    dp_F() const { return hasDamping() ? dpProps_->dp_F_ : DVect(0.0); }
		void     dp_F(const DVect &f) { if (!hasDamping()) return; dpProps_->dp_F_ = f; }

		bool    hasEnergies() const { return energies_ ? true : false; }
		double  estrain() const { return hasEnergies() ? energies_->estrain_ : 0.0; }
		void    estrain(const double &d) { if (!hasEnergies()) return; energies_->estrain_ = d; }
		double  errstrain() const { return hasEnergies() ? energies_->errstrain_ : 0.0; }
		void    errstrain(const double &d) { if (!hasEnergies()) return; energies_->errstrain_ = d; }
		double  eslip() const { return hasEnergies() ? energies_->eslip_ : 0.0; }
		void    eslip(const double &d) { if (!hasEnergies()) return; energies_->eslip_ = d; }
		double  errslip() const { return hasEnergies() ? energies_->errslip_ : 0.0; }
		void    errslip(const double &d) { if (!hasEnergies()) return; energies_->errslip_ = d; }
		double  edashpot() const { return hasEnergies() ? energies_->edashpot_ : 0.0; }
		void    edashpot(const double &d) { if (!hasEnergies()) return; energies_->edashpot_ = d; }
		
		uint inheritanceField() const { return inheritanceField_; }
		void inheritanceField(uint i) { inheritanceField_ = i; }

		const DVect2 & effectiveTranslationalStiffness()  const { return effectiveTranslationalStiffness_; }
		void           effectiveTranslationalStiffness(const DVect2 &v) { effectiveTranslationalStiffness_ = v; }
		const DAVect & effectiveRotationalStiffness()  const { return effectiveRotationalStiffness_; }
		void           effectiveRotationalStiffness(const DAVect &v) { effectiveRotationalStiffness_ = v; }

		// Rolling resistance methods
		const double & res_fric() const { return res_fric_; }
		void           res_fric(const double &d) { res_fric_ = d; }
		const DAVect & res_M() const { return res_M_; }
		void           res_M(const DAVect &f) { res_M_ = f; }
		bool           res_S() const { return res_S_; }
		void           res_S(bool b) { res_S_ = b; }
		const double & kr() const { return kr_; }
		void           kr(const double &d) { kr_ = d; }
		const double & fr() const { return fr_; }
		void           fr(const double &d) { fr_ = d; }

	private:
		// Index - used internally by PFC. Should be set to -1 in the cpp file. 
		static int index_;

		// Structure to store the energies. 
		struct Energies {
			Energies() : estrain_(0.0), errstrain_(0.0), eslip_(0.0), errslip_(0.0), edashpot_(0.0) {}
			double estrain_;   // elastic energy stored in EEPA group 
			double errstrain_; // elastic energy stored in rolling resistance group
			double eslip_;     // work dissipated by friction 
			double errslip_;   // work dissipated by rolling resistance friction 
			double edashpot_;  // work dissipated by dashpots
		};

		// Structure to store dashpot quantities. 
		struct dpProps {
			dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)) {}
			double dp_nratio_;		// normal viscous critical damping ratio
			double dp_sratio_;		// shear  viscous critical damping ratio
			int    dp_mode_;		// for viscous mode (0-1) 0 = no cut-offs;    1 = shear damp cut-off when sliding
			DVect  dp_F_;			// Force in the dashpots
		};

		bool   updateEndStiffCoef(const IContactMechanical *con);
		bool   updateEndFric(const IContactMechanical *con);
		bool   updateEndResFric(const IContactMechanical *con);
		void   updateStiffCoef(const IContactMechanical *con);
		void   updateEffectiveStiffness(ContactModelMechanicalState *state);
		void   setDampCoefficients(const double &mass, double *kn_tangent, double *ks_tangent, double *vcn, double *vcs);

		// Contact model inheritance fields.
		quint32 inheritanceField_;

		// Effective translational stiffness.
		DVect2  effectiveTranslationalStiffness_;
		DAVect  effectiveRotationalStiffness_;      // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)

		// EEPA model properties 
		double      shear_;				// Shear modulus
		double      poiss_;				// Poisson's ratio
		double      fric_;				// Coulomb friction coefficient
		DVect       eepa_F_;			// Force carried in the eepa model
		bool        eepa_S_;			// The current slip state
		double      rgap_;				// Reference gap 
		dpProps *   dpProps_;			// The viscous properties
		double		Moverlap_;			// Maximum overlap - updated as it evolves
		double		Poverlap_exp_;	    // Plastic overlap raised to the power lu_exp - updated as it evolves
		double		plas_ratio_;		// plasticity ratio
		double		lu_exp_;			// load-unload exponent anlong the k1 and k2 brances
		double		pull_off_;			// constant pull-off force
		double		adh_exp_;			// adhesive branch exponent
		double		surf_adh_;			// surface adhesion energy
		double		lu_exp_inv_;		// inverse of load-unloading exponent
		double      k1_;				// k1 branch stiffness
		double		k2_;				// k2 branch stiffness
		double      kadh_;				// adhesive branch stiffness - updated as the maximum overlap (and plastic overlap) increases
		double      ks_;				// shear stiffness
		double      r_hertz_;    		// effective contact radius = (R1*R2)/(R1 + R2)
		double      ks_fac_;			// shear stiffness scaling factor
		double		f_min_;				// minimum adhesion force limit
		int         branch_;			// 1,2 or 3 for the current branch
		
		// rolling resistance properties
		double res_fric_;       // rolling friction coefficient
		DAVect res_M_;          // moment (bending only)         
		bool   res_S_;          // The current rolling resistance slip state
		double kr_;             // bending rotational stiffness (read-only, calculated internaly) 
		double fr_;             // rolling friction coefficient (rbar*res_fric_) (calculated internaly, not a property)
		double rbar_square_;	// curvature expression used in calculating the rolling stiffness - stays constant for  given contact where kr_ = rbar_square*ks_tangent

		Energies *   energies_; // The energies

	};
} // namespace cmodelsxd
// EoF

Top

contactmodeleepa.cpp

// contactmodeleepa.cpp
#include "contactmodeleepa.h"

#include "module/interface/icontactmechanical.h"
#include "module/interface/icontact.h"
#include "module/interface/ipiecemechanical.h"
#include "module/interface/ipiece.h"
#include "module/interface/ifishcalllist.h"

#include "../version.txt"

#include "utility/src/tptr.h"
#include "shared/src/mathutil.h"

#include "kernel/interface/iprogram.h"
#include "module/interface/icontactthermal.h"
#include "contactmodel/src/contactmodelthermal.h"
#include "fish/src/parameter.h"

#ifdef EEPA_LIB
int __stdcall DllMain(void *, unsigned, void *) {
	return 1;
}

extern "C" EXPORT_TAG const char *getName() {
#if DIM==3
	return "contactmodelmechanical3deepa";
#else
	return "contactmodelmechanical2deepa";
#endif
}

extern "C" EXPORT_TAG unsigned getMajorVersion() {
	return MAJOR_VERSION;
}

extern "C" EXPORT_TAG unsigned getMinorVersion() {
	return MINOR_VERSION;
}

extern "C" EXPORT_TAG void *createInstance() {
	cmodelsxd::ContactModelEEPA *m = new cmodelsxd::ContactModelEEPA();
	return (void *)m;
}
#endif 

namespace cmodelsxd {
	static const quint32 shearMask = 0x00000002; // Base 1!
	static const quint32 poissMask = 0x00000004;
	static const quint32 fricMask = 0x00000008;
	static const quint32 resFricMask = 0x00004000;

	using namespace itasca;

	int ContactModelEEPA::index_ = -1;
	UInt ContactModelEEPA::getMinorVersion() const { return MINOR_VERSION; }

	ContactModelEEPA::ContactModelEEPA() : inheritanceField_(shearMask | poissMask | fricMask | resFricMask)
		, effectiveTranslationalStiffness_(DVect2(0.0))
		, effectiveRotationalStiffness_(DAVect(0.0))
		, shear_(0.0)
		, poiss_(0.0)
		, fric_(0.0)
		, eepa_F_(DVect(0.0))
		, eepa_S_(false)
		, rgap_(0.0)
		, dpProps_(0)
		, res_fric_(0.0)
		, res_M_(DAVect(0.0))
		, res_S_(false)
		, kr_(0.0)
		, fr_(0.0)
		, Moverlap_(0.0)			
		, Poverlap_exp_(0.0)
		, plas_ratio_(0.5)
		, lu_exp_(1.5)
		, pull_off_(0.0)
		, adh_exp_(1.5)
		, surf_adh_(0.0)
		, lu_exp_inv_(0.0)
		, k1_(0.0)			
		, k2_(0.0)			
		, kadh_(0.0)
		, ks_(0.0)					
		, r_hertz_(0.0)
		, ks_fac_(1.0)
		, f_min_(0.0)
		, branch_(1)
		, rbar_square_(0.0)		
		, energies_(0) {
	}

	ContactModelEEPA::~ContactModelEEPA() {
		// Make sure to clean up after yourself!
		if (dpProps_)
			delete dpProps_;
		if (energies_)
			delete energies_;
	}

	void ContactModelEEPA::archive(ArchiveStream &stream) {
		// The stream allows one to archive the values of the contact model
		// so that it can be saved and restored. The minor version can be
		// used here to allow for incremental changes to the contact model too. 
		stream & shear_;
		stream & poiss_;
		stream & fric_;
		stream & eepa_F_;
		stream & eepa_S_;
		stream & rgap_;
		stream & res_fric_;
		stream & res_M_;
		stream & res_S_;
		stream & kr_;
		stream & fr_;
		stream & Moverlap_;									
		stream & Poverlap_exp_;
		stream & plas_ratio_;
		stream & lu_exp_;
		stream & pull_off_;
		stream & adh_exp_;
		stream & surf_adh_;
		stream & lu_exp_inv_;
		stream & k1_;			
		stream & k2_;			
		stream & kadh_;
		stream & ks_;										
		stream & r_hertz_;
		stream & ks_fac_;
		stream & f_min_;
		stream & branch_;
		stream & rbar_square_;							

		if (stream.getArchiveState() == ArchiveStream::Save) {
			bool b = false;
			if (dpProps_) {
				b = true;
				stream & b;
				stream & dpProps_->dp_nratio_;
				stream & dpProps_->dp_sratio_;
				stream & dpProps_->dp_mode_;
				stream & dpProps_->dp_F_;
			}
			else
				stream & b;

			b = false;
			if (energies_) {
				b = true;
				stream & b;
				stream & energies_->estrain_;
				stream & energies_->errstrain_;
				stream & energies_->eslip_;
				stream & energies_->errslip_;
				stream & energies_->edashpot_;
			}
			else
				stream & b;
		}
		else {
			bool b(false);
			stream & b;
			if (b) {
				if (!dpProps_)
					dpProps_ = NEWC(dpProps());
				stream & dpProps_->dp_nratio_;
				stream & dpProps_->dp_sratio_;
				stream & dpProps_->dp_mode_;
				stream & dpProps_->dp_F_;
			}
			stream & b;
			if (b) {
				if (!energies_)
					energies_ = NEWC(Energies());
				stream & energies_->estrain_;
				stream & energies_->errstrain_;
				stream & energies_->eslip_;
				stream & energies_->errslip_;
				stream & energies_->edashpot_;
			}
		}

		stream & inheritanceField_;
		stream & effectiveTranslationalStiffness_;
		stream & effectiveRotationalStiffness_;

	}

	void ContactModelEEPA::copy(const ContactModel *cm) {
		// Copy all of the contact model properties. Used in the CMAT 
		// when a new contact is created. 
		ContactModelMechanical::copy(cm);
		const ContactModelEEPA *in = dynamic_cast<const ContactModelEEPA*>(cm);
		if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
		shear(in->shear());
		poiss(in->poiss());
		fric(in->fric());
		eepa_F(in->eepa_F());
		eepa_S(in->eepa_S());
		rgap(in->rgap());
		res_fric(in->res_fric());
		res_M(in->res_M());
		res_S(in->res_S());
		kr(in->kr());
		fr(in->fr());
		Moverlap(in->Moverlap());					
		Poverlap_exp(in->Poverlap_exp());
		plas_ratio(in->plas_ratio());
		lu_exp(in->lu_exp());
		pull_off(in->pull_off());
		adh_exp(in->adh_exp());
		surf_adh(in->surf_adh());
		lu_exp_inv(in->lu_exp_inv());
		k1(in->k1());			
		k2(in->k2());			
		kadh(in->kadh());
		ks(in->ks());	
		r_hertz(in->r_hertz());
		ks_fac(in->ks_fac());
		f_min(in->f_min());
		branch(in->branch());
		rbar_square(in->rbar_square());				

		if (in->hasDamping()) {
			if (!dpProps_)
				dpProps_ = NEWC(dpProps());
			dp_nratio(in->dp_nratio());
			dp_sratio(in->dp_sratio());
			dp_mode(in->dp_mode());
			dp_F(in->dp_F());
		}
		if (in->hasEnergies()) {
			if (!energies_)
				energies_ = NEWC(Energies());
			estrain(in->estrain());
			errstrain(in->errstrain());
			eslip(in->eslip());
			errslip(in->errslip());
			edashpot(in->edashpot());
		}
		inheritanceField(in->inheritanceField());
		effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
		effectiveRotationalStiffness(in->effectiveRotationalStiffness());
	}

	QVariant ContactModelEEPA::getProperty(uint i, const IContact *con) const {
		// Return the property. The IContact pointer is provided so that 
		// more complicated properties, depending on contact characteristics,
		// can be calcualted. 
		QVariant var;
		switch (i) {
		case kwShear:     return shear_;
		case kwPoiss:     return poiss_;
		case kwFric:      return fric_;
		case kwEepaF:     var.setValue(eepa_F_); return var;
		case kwEepaS:     return eepa_S_;
		case kwRGap:      return rgap_;
		case kwDpNRatio:  return dpProps_ ? dpProps_->dp_nratio_ : 0;
		case kwDpSRatio:  return dpProps_ ? dpProps_->dp_sratio_ : 0;
		case kwDpMode:    return dpProps_ ? dpProps_->dp_mode_ : 0;
		case kwDpF: {
			dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
			return var;
		}
		case kwResFric:    return res_fric_;
		case kwResMoment:  var.setValue(res_M_); return var;
		case kwResS:       return res_S_;
		case kwOverlapMax: return Moverlap_;
		case kwPlasRat:	   return plas_ratio_;
		case kwLuExp:	   return lu_exp_;
		case kwPullOff:	   return pull_off_;
		case kwAdhExp:	   return adh_exp_;
		case kwSurfAdh:	   return surf_adh_;
		case kwKsFac:	   return ks_fac_;															
		case kwFmin:	   return f_min_;
		}
		assert(0);
		return QVariant();
	}

	bool ContactModelEEPA::getPropertyGlobal(uint i) const {
		// Returns whether or not a property is held in the global axis system (TRUE)
		// or the local system (FALSE). Used by the plotting logic.
		switch (i) {
		case kwEepaF:
		case kwDpF:
		case kwResMoment:
			return false;
		}
		return true;
	}

	bool ContactModelEEPA::setProperty(uint i, const QVariant &v, IContact *) {
		// Set a contact model property. Return value indicates that the timestep
		// should be recalculated. 
		dpProps dp;
		switch (i) {
		case kwShear: {
			if (!v.canConvert<double>())
				throw Exception("eepa_shear must be a double.");
			double val(v.toDouble());
			if (val <= 0.0)
				throw Exception("zero or negative eepa_shear not allowed.");
			shear_ = val;
			return true;
		}
		case kwPoiss: {
			if (!v.canConvert<double>())
				throw Exception("eepa_poiss must be a double.");
			double val(v.toDouble());
			if (val < 0.0 || val > 0.5)
				throw Exception("eepa_poiss must be in the range [0, 0.5].");
			poiss_ = val;
			return true;
		}
		case kwFric: {
			if (!v.canConvert<double>())
				throw Exception("fric must be a double.");
			double val(v.toDouble());
			if (val < 0.0)
				throw Exception("negative fric not allowed.");
			fric_ = val;
			return false;
		}
		case kwEepaF: {
			if (!v.canConvert<DVect>())
				throw Exception("eepa_force must be a vector.");
			DVect val(v.value<DVect>());
			eepa_F_ = val;
			return false;
		}
		case kwRGap: {
			if (!v.canConvert<double>())
				throw Exception("reference gap must be a double.");
			double val(v.toDouble());
			rgap_ = val;
			return false;
		}
		case kwDpNRatio: {
			if (!v.canConvert<double>())
				throw Exception("dp_nratio must be a double.");
			double val(v.toDouble());
			if (val < 0.0)
				throw Exception("negative dp_nratio not allowed.");
			if (val == 0.0 && !dpProps_)
				return false;
			if (!dpProps_)
				dpProps_ = NEWC(dpProps());
			dpProps_->dp_nratio_ = val;
			return true;
		}
		case kwDpSRatio: {
			if (!v.canConvert<double>())
				throw Exception("dp_sratio must be a double.");
			double val(v.toDouble());
			if (val < 0.0)
				throw Exception("negative dp_sratio not allowed.");
			if (val == 0.0 && !dpProps_)
				return false;
			if (!dpProps_)
				dpProps_ = NEWC(dpProps());
			dpProps_->dp_sratio_ = val;
			return true;
		}
		case kwDpMode: {
			if (!v.canConvert<int>())
				throw Exception("the viscous mode dp_mode must be 0 or 1");
			int val(v.toInt());
			if (val == 0 && !dpProps_)
				return false;
			if (val < 0 || val > 1)
				throw Exception("the viscous mode dp_mode must be 0 or 1.");
			if (!dpProps_)
				dpProps_ = NEWC(dpProps());
			dpProps_->dp_mode_ = val;
			return false;
		}
		case kwDpF: {
			if (!v.canConvert<DVect>())
				throw Exception("dp_force must be a vector.");
			DVect val(v.value<DVect>());
			if (val.fsum() == 0.0 && !dpProps_)
				return false;
			if (!dpProps_)
				dpProps_ = NEWC(dpProps());
			dpProps_->dp_F_ = val;
			return false;
		}
		case kwResFric: {
			if (!v.canConvert<double>())
				throw Exception("res_fric must be a double.");
			double val(v.toDouble());
			if (val < 0.0)
				throw Exception("negative res_fric not allowed.");
			res_fric_ = val;
			return true;
		}
		case kwResMoment: {
			DAVect val(0.0);
#ifdef TWOD               
			if (!v.canConvert<DAVect>() && !v.canConvert<double>())
				throw Exception("res_moment must be an angular vector.");
			if (v.canConvert<DAVect>())
				val = DAVect(v.value<DAVect>());
			else
				val = DAVect(v.toDouble());
#else
			if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
				throw Exception("res_moment must be an angular vector.");
			if (v.canConvert<DAVect>())
				val = DAVect(v.value<DAVect>());
			else
				val = DAVect(v.value<DVect>());
#endif
			res_M_ = val;
			return false;
		}				
		case kwPlasRat: {																				
			if (!v.canConvert<double>())
				throw Exception("plas_ratio must be a double.");
			 double val(v.toDouble());
			 if (val <= 0.0 || val >= 1.0)
			     throw Exception("plas_ratio must be larger than zero and smaller than 1 (0,1).");
			 plas_ratio_ = val;
			 return true;
		}
		case kwLuExp: {
			if (!v.canConvert<double>())
				throw Exception("loading-unloading branch exponent must be a double.");
			double val(v.toDouble());
			if (val < 1.0)
				throw Exception("lu_exp must be 1 or larger. For lu_exp = 1, use the linear contact model 'eepa_lin'.");
			lu_exp_ = val;
			return true;
		}
		case kwPullOff: {
			if (!v.canConvert<double>())
				throw Exception("constant pull-off force must be a double.");
			double val(v.toDouble());
			if (val > 0.0)
				throw Exception("pull_off must be a negative value or zero.");
			pull_off_ = val;
			return false;
		}
		case kwAdhExp: {
			if (!v.canConvert<double>())
				throw Exception("adhesion branch exponent must be a double.");
			double val(v.toDouble());
			if (val < 1.0)
				throw Exception("adh_exp must be equal to or greater than 1.");
			adh_exp_ = val;
			return true;
		}
		case kwSurfAdh: {
			if (!v.canConvert<double>())
				throw Exception("surface adhesion energy must be a double.");
			double val(v.toDouble());
			if (val < 0.0)
				throw Exception("negative surf_adh not allowed.");
			surf_adh_ = val;
			return false;
		}																								
		case kwKsFac: {
			if (!v.canConvert<double>())
				throw Exception("shear stiffness factor must be a double.");
			double val(v.toDouble());
			if (val <= 0.0)
				throw Exception("zero or negative ks_fac not allowed.");
			ks_fac_ = val;
			return true;
		}																								

		}//switch
		return false;
	}

	bool ContactModelEEPA::getPropertyReadOnly(uint i) const {
		// Returns TRUE if a property is read only or FALSE otherwise. 
		switch (i) {
		case kwDpF:
		case kwEepaS:
		case kwResS:
		case kwOverlapMax:
		case kwFmin:
			return true;
		default:
			break;
		}
		return false;
	}

	bool ContactModelEEPA::supportsInheritance(uint i) const {
		// Returns TRUE if a property supports inheritance or FALSE otherwise. 
		switch (i) {
		case kwShear:
		case kwPoiss:
		case kwFric:
		case kwResFric:
			return true;
		default:
			break;
		}
		return false;
	}

	double ContactModelEEPA::getEnergy(uint i) const {
		// Return an energy value. 
		double ret(0.0);
		if (!energies_)
			return ret;
		switch (i) {
		case kwEStrain:    return energies_->estrain_;
		case kwERRStrain:  return energies_->errstrain_;
		case kwESlip:      return energies_->eslip_;
		case kwERRSlip:    return energies_->errslip_;
		case kwEDashpot:   return energies_->edashpot_;
		}
		assert(0);
		return ret;
	}

	bool ContactModelEEPA::getEnergyAccumulate(uint i) const {
		// Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
		switch (i) {
		case kwEStrain:   return true;
		case kwERRStrain: return false;
		case kwESlip:     return true;
		case kwERRSlip:   return true;
		case kwEDashpot:  return true;
		}
		assert(0);
		return false;
	}

	void ContactModelEEPA::setEnergy(uint i, const double &d) {
		// Set an energy value. 
		if (!energies_) return;
		switch (i) {
		case kwEStrain:    energies_->estrain_ = d;   return;
		case kwERRStrain:  energies_->errstrain_ = d; return;
		case kwESlip:      energies_->eslip_ = d;     return;
		case kwERRSlip:    energies_->errslip_ = d;   return;
		case kwEDashpot:   energies_->edashpot_ = d;  return;
		}
		assert(0);
		return;
	}

	bool ContactModelEEPA::validate(ContactModelMechanicalState *state, const double &) {
		// Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
		// (1) the contact is created, (2) a property of the contact that returns a true via
		// the setProperty method has been modified and (3) when a set of cycles is executed
		// via the {cycle N} command.
		// Return value indicates contact activity (TRUE: active, FALSE: inactive).
		assert(state);
		const IContactMechanical *c = state->getMechanicalContact();
		assert(c);
		//
		if (state->trackEnergy_)
			activateEnergy();
		//		
		if ((inheritanceField_ & shearMask) || (inheritanceField_ & poissMask))
			updateEndStiffCoef(c);
		if (inheritanceField_ & fricMask)
			updateEndFric(c);
		if (inheritanceField_ & resFricMask)
			updateEndResFric(c);
		//
		if (shear_ <= 0.0) {
			throw Exception("'eepa_shear' must be specified using a value larger than zero in the 'eepa' model.");
		}
		//
		updateStiffCoef(c);												// calculate the stiffness values based on the material properties - specified directly or via inheritance
		updateEffectiveStiffness(state);								// effective stiffness for translation and rotation used in the time step estimation
		//
		return checkActivity(state->gap_);
	}

	void ContactModelEEPA::updateStiffCoef(const IContactMechanical *con) {
		//
		double c12 = con->getEnd1Curvature().y();
		double c22 = con->getEnd2Curvature().y();
		r_hertz_ = c12 + c22;
		if (r_hertz_ == 0.0)
			throw Exception("eepa contact model undefined for 2 non-curved surfaces");
		r_hertz_ = 1.0 / r_hertz_;										// R1*R2/(R1 + R2)
		k1_ = 4.0 / 3.0 * (shear_ / (1 - poiss_)) * sqrt(r_hertz_);	    // k1 loading branch 
		double shear_effective = shear_ / (4.0 - 2.0*poiss_);			// effective shear modulus for contact assuming the two contact pieces are identical in properties
		ks_ = ks_fac_ * 8.0*shear_effective*sqrt(r_hertz_);				// shear stiffness
		k2_ = k1_ / (1.0 - plas_ratio_);								// k2 loading branch
		lu_exp_inv_ = 1.0 / lu_exp_;									// inverse of load-unloading exponent - stored with contact
		// rolling resistance 
		kr_ = 0.0;
		fr_ = 0.0;
		if (res_fric_ > 0.0) {
			rbar_square_ = r_hertz_ * r_hertz_;							// store this value with contact - used again when 'ks_tangent' is updated during loading-unloading to calculate the rolling stiffness
			fr_ = res_fric_ * r_hertz_;									// store this value with contact - used in force-displacement calculation
			kr_ = ks_ * rbar_square_;									// based on 'ks_' for now, this is updated and based on the shear tangent stiffness during loading-unloading
		}
	}

	bool ContactModelEEPA::endPropertyUpdated(const QString &name, const IContactMechanical *c) {
		// The endPropertyUpdated method is called whenever a surface property (with a name
		// that matches an inheritable contact model property name) of one of the contacting
		// pieces is modified. This allows the contact model to update its associated
		// properties. The return value denotes whether or not the update has affected
		// the time step computation (by having modified the translational or rotational
		// tangent stiffnesses). If true is returned, then the time step will be recomputed.  
		assert(c);
		QStringList availableProperties = getProperties().simplified().replace(" ", "").split(",", QString::SkipEmptyParts);
		QRegExp rx(name, Qt::CaseInsensitive);
		int idx = availableProperties.indexOf(rx) + 1;
		bool ret = false;

		if (idx <= 0)
			return ret;

		switch (idx) {
		case kwShear: { //Shear
			if (inheritanceField_ & shearMask)
				ret = true;							// return 'true' to ensure 'validate()' is called to affect the change in Shear Modulus through 'updateEndStiffCoef()' and effective (tangent) stiffnesses
			break;
		}
		case kwPoiss: { //Poisson's
			if (inheritanceField_ & poissMask)
				ret =true;							// return 'true' to ensure 'validate()' is called to affect the change in Poisson's ratio through 'updateEndStiffCoef()' and effective (tangent) stiffnesses
			break;
		}
		case kwFric: { //fric
			if (inheritanceField_ & fricMask)
				updateEndFric(c);					// friction does not influence any of the other parameters or time step size, so directly update the contact model friction
			break;
		}
		case kwResFric: { //rr_fric
			if (inheritanceField_ & resFricMask)
				ret = true;							// return 'true' to ensure 'validate()' is called to affect the change in rolling friction through 'updateEndResFric()' and 'fr_'
			break;
		}
		//
		}
		return ret;
	}

	static const QString gstr("eepa_shear");
	static const QString nustr("eepa_poiss");
	bool ContactModelEEPA::updateEndStiffCoef(const IContactMechanical *con) {
		assert(con);
		double g1 = shear_;
		double g2 = shear_;
		double nu1 = poiss_;
		double nu2 = poiss_;
		QVariant vg1 = con->getEnd1()->getProperty(gstr);
		QVariant vg2 = con->getEnd2()->getProperty(gstr);
		QVariant vnu1 = con->getEnd1()->getProperty(nustr);
		QVariant vnu2 = con->getEnd2()->getProperty(nustr);
		if (vg1.isValid() && vg2.isValid()) {
			g1 = vg1.toDouble();
			g2 = vg2.toDouble();
			if (g1 <= 0.0 || g2 <= 0.0)
				throw Exception("Negative or zero shear modulus not allowed in eepa contact model");
		}
		if (vnu1.isValid() && vnu2.isValid()) {
			nu1 = vnu1.toDouble();
			nu2 = vnu2.toDouble();
			if (nu1 < 0.0 || nu1 > 0.5 || nu2 < 0.0 || nu2 > 0.5)
				throw Exception("Poisson ratio should be in range [0,0.5] in eepa contact model");
		}
		if (g1*g2 == 0.0) return false;
		double es = 1.0 / ((1.0 - nu1) / (2.0*g1) + (1.0 - nu2) / (2.0*g2));
		double gs = 1.0 / ((2.0 - nu1) / g1 + (2.0 - nu2) / g2);
		poiss_ = (4.0*gs - es) / (2.0*gs - es);
		shear_ = 2.0*gs*(2 - poiss_);
		if (shear_ <= 0.0)
			throw Exception("Negative or zero shear modulus not allowed in eepa contact model");
		if (poiss_ < 0.0 || poiss_ > 0.5)
			throw Exception("Poisson ratio should be in range [0,0.5] in eepa contact model");
		return true;
	}

	static const QString fricstr("fric");
	bool ContactModelEEPA::updateEndFric(const IContactMechanical *con) {
		assert(con);
		QVariant v1 = con->getEnd1()->getProperty(fricstr);
		QVariant v2 = con->getEnd2()->getProperty(fricstr);
		if (!v1.isValid() || !v2.isValid())
			return false;
		double fric1 = std::max(0.0, v1.toDouble());
		double fric2 = std::max(0.0, v2.toDouble());
		double val = fric_;
		fric_ = std::min(fric1, fric2);
		if (fric_ <= 0.0)
			throw Exception("Negative friction not allowed in eepa contact model");
		return ((fric_ != val));
	}

	static const QString rfricstr("rr_fric");
	bool ContactModelEEPA::updateEndResFric(const IContactMechanical *con) {
		assert(con);
		QVariant v1 = con->getEnd1()->getProperty(rfricstr);
		QVariant v2 = con->getEnd2()->getProperty(rfricstr);
		if (!v1.isValid() || !v2.isValid())
			return false;
		double rfric1 = std::max(0.0, v1.toDouble());
		double rfric2 = std::max(0.0, v2.toDouble());
		double val = res_fric_;
		res_fric_ = std::min(rfric1, rfric2);
		if (res_fric_ <= 0.0)
			throw Exception("Negative rolling friction not allowed in eepa contact model");
		return ((res_fric_ != val));
	}
	
	void ContactModelEEPA::updateEffectiveStiffness(ContactModelMechanicalState *state) {
		//
		double overlap = rgap_ - state->gap_;
		//
		double kn_tangent = 0.0;
		double ks_tangent = 0.0;
		//
		if (overlap <= 0.0) {													// inactive contact
			// The assumption below is based on the Hill contact model in PFC
			// 0.01% of Diameter = 0.0001*D			
			// For monodisperse system, r_hertz_ = 0.5*R where R is the particle radius
			// r_hertz_ = 0.5*(0.5*D) = 0.25*D where D is the particle diameter
			// D = 4*r_hertz_
			// 0.01%*D = 0.0001*D = 0.0001*4*r_hertz_ = 0.0004*r_hertz_
			overlap = 0.0004*r_hertz_;
			kn_tangent = lu_exp_ * k1_*pow(overlap, lu_exp_ - 1.0);				// k1 loading branch assumed for inactive contacts
			ks_tangent = ks_ * sqrt(overlap);									// tangent stiffness in the shear direction
		}
		else {																	// active contact
			switch (branch_) {
			case 1: {
				kn_tangent = lu_exp_ * k1_*pow(overlap, lu_exp_ - 1.0);
				break;
			}
			case 2: {
				kn_tangent = lu_exp_ * k2_*pow(overlap, lu_exp_ - 1.0);
				break;
			}
			case 3: {
				kn_tangent = adh_exp_ * kadh_*pow(overlap, adh_exp_ - 1.0);
				break;
			}
			}//switch
			ks_tangent = ks_ * sqrt(overlap);
		}
		//
		DVect2 retT(kn_tangent, ks_tangent);
		if (res_fric_ > 0.0) {													// rolling
			kr_ = ks_tangent * rbar_square_;									// based on the tangent shear stiffness (Wensrich and Katterfeld, 2012)
		}
		//
		// correction if viscous damping active
		if (dpProps_) {
			DVect2 correct(1.0);
			if (dpProps_->dp_nratio_)
				correct.rx() = sqrt(1.0 + dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
			if (dpProps_->dp_sratio_)
				correct.ry() = sqrt(1.0 + dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
			retT /= (correct*correct);
		}
		//
		effectiveTranslationalStiffness_ = retT;									// set
		effectiveRotationalStiffness_ = DAVect(kr_);								// Effective rotational stiffness (bending only)
#if DIM==3 
		effectiveRotationalStiffness_.rx() = 0.0;
#endif
	}

	bool ContactModelEEPA::forceDisplacementLaw(ContactModelMechanicalState *state, const double &timestep) {
		assert(state);
		//
		// Current overlap
		double overlap = rgap_ - state->gap_;
		// Relative translational increment
		DVect trans = state->relativeTranslationalIncrement_;
		//
		// The contact was just activated from an inactive state
		if (state->activated()) {
			// The contact just got activated, set the initial normal force equal to the pull-off force
			eepa_F_.rdof(0) = pull_off_;												
			//
			// Trigger the FISH callback if one is hooked up to the 
			// contact_activated event.
			if (cmEvents_[fActivated] >= 0) {
				// The contact was just activated from an inactive state
				// Trigger the FISH callback if one is hooked up to the 
				// contact_activated event.
				auto c = state->getContact();
				std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
				IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
				fi->setCMFishCallArguments(c, arg, cmEvents_[fActivated]);
			}
		}
		//
		// Angular dispacement increment.
		DAVect ang = state->relativeAngularIncrement_;
		//
		// EEPA NORMAL FORCE ==============================================================================================================	
		//
		double overlap_exp = pow(overlap, lu_exp_);										// current overlap raised to the power "lu_exp"
		double eepa_Fn = eepa_F_.dof(0);												// current normal force value
		DVect eepa_F_old = eepa_F_;														// store for energy computations
		//
		// The updates below only occur when the maximum overlap is exceeded
		if (overlap > Moverlap_) {														// the current overlap is exceeding the maximum it has ever been
			Moverlap_ = overlap;														// update the maximum overlap
			//
			double radius = sqrt(pow(Poverlap_exp_,lu_exp_inv_)*2.0*r_hertz_);			// contact patch radius at plastic overlap   
			f_min_ = pull_off_ - 1.5*M_PI*surf_adh_ *radius;							// update the minimum force based on the surface adhesion energy and the current contact surface area
			//
			double f_min_limit = eepa_Fn - k2_ * overlap_exp;							// the limit of Fmin that can be reached at zero overlap for the current force value
			if (f_min_ < f_min_limit) {
				f_min_ = 0.5*(pull_off_ + f_min_limit);									// if Fmin is less than the limit, set it to half the limit (Morissey, 2013)
			}
			//
			double overlap_min_exp = overlap_exp - (eepa_Fn - f_min_) / k2_;			// the minimum overlap raised to the power lu_exp - the overlap where the k2 unloading branch goes over into the adhesion branch
			kadh_ = (pull_off_ - f_min_) / (pow(overlap_min_exp, adh_exp_*lu_exp_inv_));// stiffness of the adhesion branch	
		}
		//
		double F1 = pull_off_ + k1_ * overlap_exp;								// k1 loading branch force for this overlap
		double F2 = pull_off_ + k2_ * (overlap_exp - Poverlap_exp_);			// k2 loading-unloading branch force for this overlap
		double Fadh = pull_off_ - kadh_ * pow(overlap, adh_exp_);				// adhesion branch force for this overlap
		//
		double kn_tangent = 0.0;												// tangent stiffness in the normal direction - used for time step calculation
		//		
		if (F2 >= F1) {															// the k1 loading branch should be followed
			eepa_Fn = F1;														// set the EEPA force to that of the k1 loading branch
			Poverlap_exp_ = overlap_exp - (F1 - pull_off_) / k2_;				// for the current force (F1), calculate the plastic overlap [not the usual equation based on maximum overlap, but the same result]  
			kn_tangent = lu_exp_ * k1_*pow(overlap, lu_exp_ - 1.0);				// set the tangent stiffness in the normal direction
			branch_ = 1;
		}
		else {																	// will only enter if F2 < F1
			if (F2 > Fadh) {													// on the k2 loading-unloading branch
				eepa_Fn = F2;													// set the EEPA force to that of the k2 loading-unloading branch
				kn_tangent = lu_exp_ * k2_*pow(overlap, lu_exp_ - 1.0); 		// set the tangent stiffness in the normal direction
				branch_ = 2;
			}
			else {																// only enter if F2 < F1 AND F2 < Fadh
				if (trans.x() >= 0.0) {											// contact is unloading along the adhesive branch (can not load along this branch)
					eepa_Fn = Fadh;												// set the EEPA force to that of the adhesive unloading branch
					kn_tangent = adh_exp_ * kadh_*pow(overlap, adh_exp_ - 1.0);	// set the tangent stiffness in the normal direction
					branch_ = 3;
				}
				else {															// loading on the adhesive branch which is not defined
					eepa_Fn += k2_ * pow(-trans.x(), lu_exp_);					// set the EEPA force by incrementing the force from the previous time-step along the k2 loading branch. 
																				// The next time step, the updated Poverlap_exp is used and the code will automatically follow the k2 loading branch since F2 is calculated from this updated Poverlap_exp 
																				// The negative of the increment should be used to ensure a positive value is rasied to the power "lu_exp". The sign is then corrected by adding the force increment to the previous value
					Poverlap_exp_ = overlap_exp + (pull_off_ - eepa_Fn) / k2_;	// update the plastic overlap (raised to the power "lu_exp") which means that the k2 loading-unloading branch has shifted in the unloading direction to start at this point
					kn_tangent = lu_exp_ * k2_*pow(overlap, lu_exp_ - 1.0);		// set the tangent stiffness in the normal direction
					branch_ = 2;
				}
			}
		}
		eepa_F_.rdof(0) = eepa_Fn;												//set the updated normal force component
		//
		// EEPA SHEAR FORCE ==============================================================================================================
		//
		double ks_tangent = ks_ * sqrt(overlap);						// tangent stiffness in the shear direction
		//
		DVect sforce(0.0);
		// dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
		// Loop over the shear components (note: the 0 component is the normal component) and calculate the shear force.
		for (int i = 1; i < dim; ++i)
			sforce.rdof(i) = eepa_F_.dof(i) - trans.dof(i) * ks_tangent;// shear force update
		//
		// The canFail flag corresponds to whether or not the contact can undergo non-linear
		// force-displacement response. If the SOLVE ELASTIC command is given then the 
		// canFail state is set to FALSE. Otherwise it is always TRUE. 
		if (state->canFail_) {
			// Resolve sliding. This is the normal force multiplied by the coefficient of friction.
			double crit = fric_ * std::abs(eepa_Fn - f_min_);														//take absolute value - contact with very small overlap might result in negative value here due to rounding
			// The is the magnitude of the shear force.
			double sfmag = sforce.mag();
			// Sliding occurs when the magnitude of the shear force is greater than the critical value.
			if (sfmag > crit) {
				// Lower the shear force to the critical value for sliding.
				double rat = crit / sfmag;
				sforce *= rat;
				// Handle the slip_change event if one has been hooked up. Sliding has commenced.  
				if (!eepa_S_ && cmEvents_[fSlipChange] >= 0) {
					auto c = state->getContact();
					std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
														 fish::Parameter((qint64)0) };
					IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
					fi->setCMFishCallArguments(c, arg, cmEvents_[fSlipChange]);
				}
				eepa_S_ = true;
			}
			else {
				// Handle the slip_change event if one has been hooked up and
				// the contact was previously sliding. Sliding has ceased.  
				if (eepa_S_) {
					if (cmEvents_[fSlipChange] >= 0) {
						auto c = state->getContact();
						std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
															 fish::Parameter((qint64)1) };
						IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
						fi->setCMFishCallArguments(c, arg, cmEvents_[fSlipChange]);
					}
					eepa_S_ = false;
				}
			}
		}
		//
		// Set the shear components of the total force.
		for (int i = 1; i < dim; ++i)
			eepa_F_.rdof(i) = sforce.dof(i);
		//
		// EEPA ROLLING RESISTANCE ==============================================================================================================
		//
		DAVect res_M_old = res_M_;													// used in energy calculation only
		if (fr_ == 0.0) {											
			res_M_.fill(0.0);
		}
		else {
			DAVect angStiff(0.0);
			DAVect MomentInc(0.0);
			kr_ = ks_tangent * rbar_square_;										// update rolling stiffness based on the tangent shear stiffness (Wensrich and Katterfeld, 2012)
#if DIM==3 
			angStiff.rx() = 0.0;
			angStiff.ry() = kr_;
#endif
			angStiff.rz() = kr_;
			MomentInc = ang * angStiff * (-1.0);
			res_M_ += MomentInc;
			if (state->canFail_) {
				// Account for bending strength
				double rmax = std::abs(fr_*(eepa_Fn - f_min_));						// Using the same normal force as in the shear limit
				double rmag = res_M_.mag();
				if (rmag > rmax) {
					double fac = rmax / rmag;
					res_M_ *= fac;
					res_S_ = true;
				}
				else {
					res_S_ = false;
				}
			}
		}
		//
		// EEPA TANGENT STIFFNESS FOR TIME STEP UPDATES ==============================================================================================================
		//
		// set effective or tangent stiffness in the normal and shear directions for timestep calculation
		effectiveTranslationalStiffness_ = DVect2(kn_tangent, ks_tangent);
		// set the effective rotational stiffness (bending only)
		effectiveRotationalStiffness_ = DAVect(kr_);
#if DIM==3 
		effectiveRotationalStiffness_.rx() = 0.0;
#endif
		// EEPA DAMPING FORCE ==============================================================================================================
		//
		// Account for dashpot forces if the dashpot structure has been defined. 
		if (dpProps_) {
			dpProps_->dp_F_.fill(0.0);
			double vcn(0.0), vcs(0.0);
			// Calculate the damping coefficients. 
			setDampCoefficients(state->inertialMass_, &kn_tangent, &ks_tangent, &vcn, &vcs);				// this is based on the current tangent loading stiffness
			// First damp the shear components
			for (int i = 1; i < dim; ++i)
				dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
			// Damp the normal component
			dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
			//
			if (eepa_S_ && dpProps_->dp_mode_ == 1) {
				// Limit in shear if sliding.
				double dfn = dpProps_->dp_F_.rx();
				dpProps_->dp_F_.fill(0.0);
				dpProps_->dp_F_.rx() = dfn;
			}
			// Correct effective translational stiffness
			DVect2 correct(1.0);
			if (dpProps_->dp_nratio_)
				correct.rx() = sqrt(1.0 + dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
			if (dpProps_->dp_sratio_)
				correct.ry() = sqrt(1.0 + dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
			effectiveTranslationalStiffness_ /= (correct*correct);
		}
		
		//Compute energies if energy tracking has been enabled. 
		if (state->trackEnergy_) {
			assert(energies_);
			// calculate the strain energy increment in the normal direction				 
			energies_->estrain_ -= 0.5*(eepa_F_old.dof(0) + eepa_Fn)*trans.dof(0);
			if (ks_tangent) {
				DVect u_s_elastic = trans;										// set the elastic displacement increment equal to the total displacement increment
				u_s_elastic.rx() = 0.0;											// set the normal component to zero: u_s_elatic = [0, trans_shear_1, trans_shear_2]
				DVect shearF = eepa_F_;											// set the shear force equal to the total EEPA force (including the normal component)
				shearF.rx() = 0.0;												// set the normal component to zero: shearF = [0, shear_force_1, shear_force_2]
				eepa_F_old.rx() = 0.0;											// set normal component of the previous force equal to zero
				DVect avg_F_s = (shearF + eepa_F_old)*0.5;						// average shear force vector
				if (eepa_S_) {													// if sliding, calculate the slip energy and accumulate it
					DVect u_s_total = u_s_elastic;								// total shear displacement increment
					u_s_elastic = (shearF - eepa_F_old) / ks_tangent;			// elastic shear displacement increment
					energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s_total + u_s_elastic))); // where (u_s_total + u_s_elatic) is the slip displacment increment (due to the sign convention, the terms are added up)
				}
				energies_->estrain_ -= avg_F_s | u_s_elastic;
			}
			// Add the rolling resistance energy contributions.							// done incrementally since the stiffness kr_ changes with the tangent shear stiffness (non-linearly)
			if (kr_) {
				DAVect t_s_elastic = ang;												// set the elastic rotation increment equal to the total angle increment
				DAVect avg_M = (res_M_ + res_M_old)*0.5;								// average moment from this time step and the previous
				if (res_S_) {															// if sliding, calculate the slip energy and accumulate it
					t_s_elastic = (res_M_ - res_M_old) / kr_;							// elastic angle increment
					energies_->errslip_ -= std::min(0.0, (avg_M | (ang + t_s_elastic)));// where (ang + t_s_elatic) is the slip rotation increment (due to the sign convention, the terms are added up)
				}
				energies_->errstrain_ -= avg_M | t_s_elastic;							// add the elastic component (if any - if previous force equals current force, the elastic incrment is zero)
			}
			// Calculate damping energy (accumulated) if the dashpots are active. 
			if (dpProps_) {
				energies_->edashpot_ -= dpProps_->dp_F_ | trans;
			}
		}

		// This is just a sanity check to ensure, in debug mode, that the force isn't wonky. 
		assert(eepa_F_ == eepa_F_);
		return true;
	}

	void ContactModelEEPA::setForce(const DVect &v, IContact *c) {
		//
	}

	DVect ContactModelEEPA::getForce(const IContactMechanical *) const {
		DVect ret(eepa_F_);																				
		if (dpProps_)
			ret += dpProps_->dp_F_;																		
		return ret;
	}

	DAVect ContactModelEEPA::getMomentOn1(const IContactMechanical *c) const {
		DVect force = getForce(c);
		DAVect ret(res_M_);
		c->updateResultingTorqueOn1Local(force, &ret);
		return ret;
	}

	DAVect ContactModelEEPA::getMomentOn2(const IContactMechanical *c) const {
		DVect force = getForce(c);
		DAVect ret(res_M_);
		c->updateResultingTorqueOn2Local(force, &ret);
		return ret;
	}

	void ContactModelEEPA::setDampCoefficients(const double &mass, double *kn_tangent, double *ks_tangent, double *vcn, double *vcs) {
		*vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(*kn_tangent));
		*vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(*ks_tangent));
	}

} // namespace cmodelsxd
// EoF

Top