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

static CM_VisDamp cm_visdamp(true);

CM_VisDamp::CM_VisDamp(bool bRegister) :
  ContactModel(ulModelCM_VisDamp, bRegister),
            kn(0.0), ks(0.0), fric(0.0),
						cn(0.0), cs(0.0),
						cstr_n(0.0), cstr_s(0.0),
						udotn(0.0), udots(0.0), tdot(0.0),
						edn(0.0), eds(0.0), edf(0.0),
						nt(false), mode(0),
						cs_force2(0.0)
{
	cs_force3.Fill(0.0);
  SetDelete(true);
}

const char *CM_VisDamp::Name(void) const
{
  return("visdamp");
}

const char **CM_VisDamp::PropNames(void) const
{
  static const char *strKey[] = {
    "visd_kn",  "visd_ks",  "visd_fric",
    "visd_cn",  "visd_cs",
		"visd_n_bond",  "visd_s_bond",
		"visd_udotn", "visd_udots", "visd_tdot",
		"visd_edn", "visd_eds", "visd_edf",
		"visd_nt",  "visd_mode",
     0
  };
  return(strKey);
}

double CM_VisDamp::ReturnProp(int n) const
{
  switch (n) {
    case  0:  return kn;
    case  1:  return ks;
    case  2:  return fric;
    case  3:  return cn;
    case  4:  return cs;
    case  5:  return cstr_n;
    case  6:  return cstr_s;
    case  7:  return udotn;		
    case  8:  return udots;
    case  9:  return tdot;
    case 10:  return edn;
    case 11:  return eds;
    case 12:  return edf;
    case 13:  return nt==false ? 0 : 1;
    case 14:  return mode;			
    default:  return 0.0;
  }
}

void CM_VisDamp::AcceptProp(int n, double v)
{
  switch (n) {
    case  0: kn   = v; break;
    case  1: ks   = v; break;
    case  2: fric = v; break;
    case  3: cn   = v; break;
    case  4: cs   = v; break;			
    case  5: break; // prohibitted to set
    case  6: break;	//        "		
    case  7: break; //        "  
    case  8: break; //        "
    case  9: break; //        " 
		case 10: break; //        "
		case 11: break; //        "
		case 12: break; //        "
    case 13: nt   = (int)v==0 ? false : true;
    case 14: mode = (int)v;
    default: break; 
  }
}

const char *CM_VisDamp::PreCycle(PropBlock &pb)
{  
	if (mode!=0 && mode!=1) return("visd_mode should be ZERO or ONE");
	cstr_n = pb.n_strength; 
  cstr_s = pb.s_strength; 
	return(0);
}

//------------ user-defined contact model "CM_VisDamp" -----------
void CM_VisDamp::FDlaw(FdBlock &fb, char cDim)
{
  if (cDim==2) {  // ... 2d section ...
		double s_force0 = fb.s_force;
		double cs_force0 = cs_force2;
		if (fb.u_n > 0.0 || fb.bflag) {   
		  // elastic part
      fb.n_force  = kn * fb.u_n;
			fb.s_force -= ks * fb.u_dot_s * fb.tdel;
			switch (mode) {
			  case 0: {
          // mode 0 : first check bond, slip, then add viscous term     
			    // check bond, slip
					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
							fb.s_force = fb.s_force < 0.0 ? -fabs(max_s_force) : fabs(max_s_force);
							fb.bSliding = true;
						}
					}
          // add viscous damping
			    fb.n_force -= cn * fb.u_dot_n;
          if (nt==true)
            fb.n_force = fb.n_force > 0.0 ? fb.n_force : 0.0;
          fb.cs_force = fb.s_force;
          if (!fb.bSliding)
				    fb.cs_force -= fb.u_dot_s * cs;   // in case without slip
					break;
				}
        case 1: {
          // mode 1 : first add viscous term, then check bond, slip
          // add viscous damping
			    fb.n_force -= cn * fb.u_dot_n;
          if (nt==true)
            fb.n_force = fb.n_force > 0.0 ? fb.n_force : 0.0;
			    fb.cs_force = fb.s_force - fb.u_dot_s * cs;   
			    // check bond, slip
					if (fb.bflag) { // bonding check 
				    if (-fb.n_force >= cstr_n)            fb.bfishc_brokenn = true;
						else if (fabs(fb.cs_force) >= cstr_s) fb.bfishc_brokens = true;
					} else { // sliding check
				    double max_s_force = fb.n_force * fric;
		        if (fabs(fb.cs_force) > max_s_force) {        // slip
			        fb.cs_force = fb.cs_force < 0.0 ? -fabs(max_s_force) : fabs(max_s_force);
							fb.s_force  = fb.cs_force + fb.u_dot_s * cs;
					    fb.bSliding = true;
						}
					}
					break;
				}
        default:
					break;
      }
			cs_force2 = fb.cs_force;
      fb.knest  = kn;
      fb.ksest  = ks;
      fb.skip   = false;
	    // dissipated energy
			edn = (nt==true && fb.n_force==0.0) ? 0.0 : cn * fb.u_dot_n * fb.u_dot_n * fb.tdel;
			eds = fb.bSliding ? 0.0 : cs * fb.u_dot_s * fb.u_dot_s * fb.tdel;
			double avF_s = mode!=1 ? 0.5*(fb.s_force+s_force0):0.5*(fb.cs_force+cs_force0);
			double delU_s    = fb.u_dot_s * fb.tdel ;
			double delU_elas = (fb.s_force - s_force0) / ks;
			double delU_slip = delU_s + delU_elas ;
			edf = -avF_s * delU_slip;
			edf = 0.0 > edf ? 0.0 : edf;
		} else {
			cs_force2   = 0.0;
      fb.n_force  = 0.0;
      fb.s_force  = 0.0;
			fb.cs_force = 0.0;
      fb.knest = 0.0;
      fb.ksest = 0.0;
      fb.skip  = true;
		  edn = 0.0;
			eds = 0.0;
			edf = 0.0;
		}
	} else {   // ... 3d section ... 
   	UMdvect s_force0 = fb_s_force;
		UMdvect cs_force0 = cs_force3;
    if (fb.u_n > 0.0 || fb.bflag) {
		  // elastic part
      fb.n_force  = kn * fb.u_n;
      double kst = ks * fb.tdel;
      UMdvect vec = fb_u_dot_s * kst;
      fb_s_force = fb_s_force - vec;
			switch (mode) {
			  case 0: {
          // mode 0 : first check bond, slip, then add viscous term     
			    // check bond, slip
          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) { // bond 
						if (-fb.n_force >= cstr_n) fb.bfishc_brokenn = true;
						else if (sfmag >= cstr_s)  fb.bfishc_brokens = true;
					} else { // slip
						double max_s_force = fb.n_force * fric;
						if (sfmag > max_s_force) {        
							double rat = max_s_force / sfmag;
							fb_s_force = fb_s_force* rat;
							fb.bSliding = true;
						}
					}
          // add viscous damping
			    fb.n_force -= cn * fb.u_dot_n;
          if (nt==true)
            fb.n_force = fb.n_force > 0.0 ? fb.n_force : 0.0;
          fb_cs_force = fb_s_force;
          if (!fb.bSliding)
				    fb_cs_force -= fb_u_dot_s * cs;   // in case without slip
					break;
				}
        case 1: { 
          // mode 1 : first add viscous term, then check bond, slip
          // add viscous damping
			    fb.n_force -= cn * fb.u_dot_n;
          if (nt==true)
            fb.n_force = fb.n_force > 0.0 ? fb.n_force : 0.0;
			    fb_cs_force = fb_s_force - fb_u_dot_s * cs;   
			    // check bond, slip
          double dX = fb.ptcs_force->x;
          double dY = fb.ptcs_force->y;
          double dZ = fb.ptcs_force->z;
          double sfmag = sqrt(dX*dX + dY*dY + dZ*dZ);
					if (fb.bflag) { // bond 
						if (-fb.n_force >= cstr_n) fb.bfishc_brokenn = true;
						else if (sfmag >= cstr_s)  fb.bfishc_brokens = true;
					} else { // slip
	          double max_s_force = fb.n_force * fric;
		        if (sfmag > max_s_force) {        
							double rat = max_s_force / sfmag;
							fb_cs_force = fb_cs_force * rat;
							fb_s_force  = fb_cs_force + fb_u_dot_s * cs;
							fb.bSliding = true;
						}
					}
					break;
				}
        default:
					break;
      }
			cs_force3 = fb_cs_force;
      fb.knest = kn;
      fb.ksest = ks;
      fb.skip  = false;
	    // dissipated energy
			edn = (nt==true && fb.n_force==0.0) ? 0.0 : cn * fb.u_dot_n * fb.u_dot_n * fb.tdel;
			if (!fb.bSliding) {
				double dX = fb.ptu_dot_s->x;
				double dY = fb.ptu_dot_s->y;
				double dZ = fb.ptu_dot_s->z;
				double sfmag2 = dX*dX + dY*dY + dZ*dZ;
				eds = cs * sfmag2 * fb.tdel;
			} else
				eds = 0.0;
			UMdvect avF_s = mode!=1 ? fb_s_force + s_force0:fb_cs_force + cs_force0;
			avF_s        *= 0.5;
			UMdvect delU_s    = fb_u_dot_s * fb.tdel;
			UMdvect delU_elas = fb_s_force - s_force0;
			delU_elas      /= ks;
			UMdvect delU_slip = delU_s + delU_elas;
			edf             = -(avF_s | delU_slip);
			edf             = 0.0 > edf ? 0.0 : edf;

		} else {
			cs_force3.Fill(0.0);
      fb.n_force = 0.0;
      fb_s_force.Fill(0.0);
			fb_cs_force.Fill(0.0);
      fb.knest = 0.0;
      fb.ksest = 0.0;
      fb.skip  = true;
		  edn = 0.0;
			eds = 0.0;
			edf = 0.0;
		}
	}
}

const char *CM_VisDamp::SaveRestore(ModelSaveObject *mso)
{
  const char *str = ContactModel::SaveRestore(mso);
  if (str) return(str);
  mso->Initialize(17,1,2);
  mso->Save( 0, kn);
  mso->Save( 1, ks);
  mso->Save( 2, fric);
  mso->Save( 3, cn);
  mso->Save( 4, cs);
  mso->Save( 5, cstr_n);
  mso->Save( 6, cstr_s);
  mso->Save( 7, udotn);
  mso->Save( 8, udots);
  mso->Save( 9, tdot);
  mso->Save(10, edn);
  mso->Save(11, eds);
  mso->Save(12, edf);
  mso->Save(13, cs_force2);
  mso->Save(14, cs_force3.x);
  mso->Save(15, cs_force3.y);
  mso->Save(16, cs_force3.z);
	mso->Save( 0, mode);
  mso->Save( 0, delete_flag);
  mso->Save( 1, nt);
  return(0);
}



/* EoF */
