#ifndef __CMODEL_BURMODEL_H
#define __CMODEL_BURMODEL_H

#ifndef __CMODEL_CONTMODEL_H
#include "base\contmodel.h"
#endif


//------------ user-defined contact model "Burger's model" -----------
const unsigned long ulModelBUR = 1001;
class CM_Burger : public ContactModel {
  private:
    double  kn_k;    // normal stiffness for Kelvin section  (bur_knk)
    double  cn_k;    // normal viscosity for Kelvin section  (bur_cnk)
    double  kn_m;    // normal stiffness for Maxwell section (bur_knm)
    double  cn_m;		 // normal viscosity for Maxwell section (bur_cnm)
    double  ks_k;		 // shear stiffness for Kelvin section   (bur_ksk)
    double  cs_k;		 // shear viscosity for Kelvin section   (bur_csk)
    double  ks_m;		 // shear stiffness for Maxwell section  (bur_ksm)
    double  cs_m;		 // shear viscosity for Maxwell section  (bur_csm)
    double  fric;		 // friction coefficient                (bur_fric)
    bool bNtension;	 // FDlaw option, with or without tensile force, default bFalse (with tensile)
    double  fn0;     // normal contact force 1 step before
    double  u_n0;    // normal total overlap 1 step before
    double  u_nk0;   // normal overlap of Kelvin part 1step before
    double  u_sk02;  // shear relative displacement of Kelvin part 1step before (2D)
		UMdvect u_sk03;  // shear relative displacement of Kelvin part 1step before (3D)
		double  conAn;   // constant A in eq.(), normal
		double  conB_An; // constant B/A in eq.(), normal
		double  conCn;   // constant C in eq.(), normal
		double  conDn;   // constant D in eq.(), normal
		double  conAs;   // constant A in eq.(), shear
		double  conB_As; // constant B/A in eq.(), shear
		double  conCs;   // constant C in eq.(), shear
		double  conDs;   // constant D in eq.(), shear
    double  tdel;    // current time
		double  cstr_n;  // contact bond normal strength  
		double  cstr_s;  // contact bond shear strength  
  public:
    // Creators
    EXPORT CM_Burger(bool bRegister=false);
    EXPORT ~CM_Burger(void) { }

    // Accessors
    EXPORT const char *Name(void) const;
    EXPORT const char **PropNames(void) const;
    EXPORT ContactModel *Clone(PropBlock * =0) const { return(new CM_Burger); }
    EXPORT double  ReturnProp(int n) const;
    EXPORT double  KnEstimate(void) const { return kn_m; }
    EXPORT double  KsEstimate(void) const { return ks_m; }
    EXPORT unsigned long   Version(void) const { return 4; }

    // Manipulators
    EXPORT void AcceptProp(int n, double v);
    EXPORT const char *PreCycle(PropBlock &pb);
    EXPORT void FDlaw(FdBlock &fb, char cDim);
    EXPORT const char *SaveRestore(ModelSaveObject *mso);
    EXPORT void remap(PropBlock &) {}

    inline double A(const double k_k, const double c_k) {	return(1.0 + k_k*tdel/(2.0*c_k)); }
    inline double B(const double k_k, const double c_k) { return(1.0 - k_k*tdel/(2.0*c_k)); }
    inline double B_A(const double k_k, const double c_k) { return(B(k_k,c_k)/A(k_k,c_k)); }
    inline double C(const double k_k, const double c_k, const double k_m, const double c_m) {
      return(tdel/(2.0*c_k*A(k_k,c_k)) + 1.0/k_m + tdel/(2.0*c_m)); }
    inline double D(const double k_k, const double c_k, const double k_m, const double c_m) {
      return(tdel/(2.0*c_k*A(k_k,c_k)) - 1.0/k_m + tdel/(2.0*c_m)); }
};



#endif
