#include "burmodel.h"
#include <math.h>

static CM_Burger bur(true);


CM_Burger::CM_Burger(bool bRegister) :
  ContactModel(ulModelBUR, bRegister),
    kn_k(0.0), cn_k(0.0), kn_m(0.0), cn_m(0.0),
    ks_k(0.0), cs_k(0.0), ks_m(0.0), cs_m(0.0),
    fric(0.0), bNtension(false),
    fn0(0.0),  u_n0(0.0), u_nk0(0.0),
    u_sk02(0.0),
    u_sk03(0.0,0.0,0.0),
		conAn(0.0), conB_An(0.0),	conCn(0.0), conDn(0.0),
		conAs(0.0), conB_As(0.0),	conCs(0.0), conDs(0.0),
		tdel(0.0),
    cstr_n(0.0), cstr_s(0.0) 	
{
}


const char *CM_Burger::Name(void) const {
  return("burger");
}


const char **CM_Burger::PropNames(void) const {
  static const char *strKey[] = {
  "bur_knk", "bur_cnk", "bur_knm", "bur_cnm",
  "bur_ksk", "bur_csk", "bur_ksm", "bur_csm",
  "bur_fric", "bur_notension",
	"bur_n_bond", "bur_s_bond",   
   0
  };
  return(strKey);
}


double CM_Burger::ReturnProp(int n) const {
  switch (n) {
    case  0: return kn_k;
    case  1: return cn_k;
    case  2: return kn_m;
    case  3: return cn_m;
    case  4: return ks_k;
    case  5: return cs_k;
    case  6: return ks_m;
    case  7: return cs_m;
    case  8: return fric;
    case  9: return bNtension ? 1.0 : 0.0;
    case 10: return cstr_n; 
    case 11: return cstr_s; 
    default: return 0.0;
  }
}


const char *CM_Burger::PreCycle(PropBlock &pb) {  
	cstr_n = pb.n_strength; 
  cstr_s = pb.s_strength; 
	if (cn_k==0.0) return("need to set bur_cnk non-zero value");
	if (cs_k==0.0) return("need to set bur_csk non-zero value");
	if (kn_m==0.0) return("need to set bur_knm non-zero value");
	if (ks_m==0.0) return("need to set bur_ksm non-zero value");
	if (cn_m==0.0) return("need to set bur_cnm non-zero value");
	if (cs_m==0.0) return("need to set bur_csm non-zero value");
	return(0);
}


void CM_Burger::AcceptProp(int n, double v) {
  switch (n) {
    case 0: kn_k = v; break;
    case 1: cn_k = v; break;
    case 2: kn_m = v; break;
    case 3: cn_m = v; break;
    case 4: ks_k = v; break;
    case 5: cs_k = v; break;
    case 6: ks_m = v; break;
    case 7: cs_m = v; break;
    case 8: fric = v; break;
    case 9: bNtension = v!=0.0 ? true : false; break;
    default:          break;
  }
  tdel = 0.0;
}


//------------ user-defined contact model "Burger's model" -----------
void CM_Burger::FDlaw(FdBlock &fb, char cDim) { // Force/displacement law
  if (fb.tdel!=tdel) { // re-calculated constants in eq.().
    tdel = fb.tdel;
		// need some protection devided by zero (k_k c_k k_m c_m = zero)
    conAn = A(kn_k, cn_k);
    double conBn = B(kn_k, cn_k);
	  conB_An = conBn / conAn;
	  conCn = C(kn_k, cn_k, kn_m, cn_m);
	  conDn = D(kn_k, cn_k, kn_m, cn_m);
    conAs = A(ks_k, cs_k);
    double conBs = B(ks_k, cs_k);
	  conB_As = conBs / conAs;
	  conCs = C(ks_k, cs_k, ks_m, cs_m);
	  conDs = D(ks_k, cs_k, ks_m, cs_m);
  }
	if (fb.u_n > 0.0 || fb.bflag) {   
    // normal force
    fb.n_force = 1.0/conCn*(fb.u_n-u_n0+(1.0-conB_An)*u_nk0-conDn*fn0);
		if(bNtension && fb.n_force<0.0) fb.n_force = 0.0;
    u_nk0 = conB_An*u_nk0+fb.tdel/(2.0*cn_k*conAn)*(fb.n_force+fn0);
    u_n0  = fb.u_n;
    fn0   = fb.n_force;
    // shear force
    if (cDim == 2) { // ... 2D section ...
      double f_old = fb.s_force;
      fb.s_force = -1.0/conCs*(fb.u_dot_s*fb.tdel+(1.0-conB_As)*u_sk02+conDs*f_old);
		  u_sk02 = conB_As*u_sk02-fb.tdel/(2.0*cs_k*conAs)*(fb.s_force+f_old);
			if (fb.bflag) { // bonding check 
        if (-fb.n_force >= cstr_n)           fb.bfishc_brokenn = true;
        else if (fabs(fb.s_force) >= cstr_s) fb.bfishc_brokens = true;
      } else { // sliding check
        double max_s_force = fb.n_force * fric;
        if (fabs(fb.s_force) > max_s_force) { // slip check
          fb.s_force = fb.s_force < 0.0 ? -fabs(max_s_force) : fabs(max_s_force);
          fb.bSliding = true;
        }
			}
    }
		else { // ... 3D section ...
      UMdvect f_old = fb_s_force;
      UMdvect v1 = fb_u_dot_s * fb.tdel;
      UMdvect v2 = u_sk03 * (1.0-conB_As);
      UMdvect v3 = f_old * conDs;
      fb_s_force = (v1+v2+v3) / conCs*(-1.0);
      double d1 = fb.tdel / (2.0*cs_k*conAs);
      v1 = fb_s_force + f_old;
      u_sk03 = u_sk03*conB_As-v1*d1;
      double dx = fb.pts_force->x;
      double dy = fb.pts_force->y;
      double dz = fb.pts_force->z;
      double sfmag = sqrt(dx*dx + dy*dy + dz*dz);
			if (fb.bflag) { // bonding check 
        if (-fb.n_force >= cstr_n) fb.bfishc_brokenn = true;
        else if (sfmag >= cstr_s)  fb.bfishc_brokens = true;
      } else { // sliding check
        double max_s_force = fb.n_force * fric;
        if (sfmag > max_s_force) { // slip check
          double rat = max_s_force / sfmag;
          fb_s_force= fb_s_force* rat;
				  fb.bSliding = true;
				}
      }
    }
    fb.knest = kn_k;
    fb.ksest = ks_k;
    fb.skip  = false;
  }
  else {
    fb.knest = 0.0;
    fb.ksest = 0.0;
    fb.skip  = true;
  }
}


const char *CM_Burger::SaveRestore(ModelSaveObject *mso) {
  const char *str = ContactModel::SaveRestore(mso);
  if (str) return(str);
  mso->Initialize(27,0,2);
  mso->Save( 0, delete_flag);
	mso->Save( 1, bNtension);
  mso->Save( 0, kn_k);
  mso->Save( 1, cn_k);
  mso->Save( 2, kn_m);
  mso->Save( 3, cn_m);
  mso->Save( 4, ks_k);
  mso->Save( 5, cs_k);
  mso->Save( 6, ks_m);
  mso->Save( 7, cs_m);
  mso->Save( 8, fric);
  mso->Save( 9, fn0);
  mso->Save(10, u_n0);
  mso->Save(11, u_nk0);
  mso->Save(12, conAn);
  mso->Save(13, conB_An);
  mso->Save(14, conCn);
  mso->Save(15, conDn);
  mso->Save(16, conAs);
  mso->Save(17, conB_As);
  mso->Save(18, conCs);
  mso->Save(19, conDs);
  mso->Save(20, tdel);
  mso->Save(21, u_sk02);
  mso->Save(22, u_sk03.x);
  mso->Save(23, u_sk03.y);
  mso->Save(24, u_sk03.z);
  mso->Save(25, cstr_n);
  mso->Save(26, cstr_s);
  return(0);
}



/* EoF */