#include "cm3model.h"
#include <string.h>
#include <math.h>


CM3::CM3(bool bRegister) :
  ContactModel(ulModelCM3, bRegister),
    kn(0.0), ratks(0.0), ratfs(0.0), fric(0.0),
    FTmax(0.0), kduc(0.0), kun(0.0), broken(0.0), softened(0.0),
    ut(0.0), FT(0.0) {
  SetDelete(false);
}


const char *CM3::Name(void) const {
  return("ductile");
}


const char **CM3::PropNames(void) const {
  static const char *strKey[] = {
    "duc_kn",   "duc_ratks",  "duc_ratfs",
    "duc_ftmax", "duc_kduc",   "duc_fric",
    "duc_kun",  "duc_broken", "duc_softened",
     0
  };
  return(strKey);
}


double CM3::ReturnProp(int n) const {
  switch (n) {
    case 0:  return kn;
    case 1:  return ratks;
    case 2:  return ratfs;
    case 3:  return FTmax;
    case 4:  return kduc;
    case 5:  return fric;
    case 6:  return kun;
    case 7:  return broken;
    case 8:  return softened;
    default: return 0.0;
  }
}


void CM3::AcceptProp(int n, double v) {
  switch (n) {
    case 0:  kn       = v;        break;
    case 1:  ratks    = v;        break;
    case 2:  ratfs    = v;        break;
    case 3:  FTmax    = v;        break;
    case 4:  kduc     = v;        break;
    case 5:  fric     = v;        break;
    case 6:  kun      = v;        break;
    default: break;
  }
}


const char *CM3::PreCycle(PropBlock &) {
	bool bMes = false;
  if (broken == 0.0) {
    if (kun == 0.0)
      kun = kn;
    if (ut == 0.0) {
      if (kun==0.0)
				bMes = true;
      ut = -FTmax / kun;
    }
    if (FT == 0.0)
      FT = -FTmax;
  }
	if (bMes) return("Zero initial unloading modulus");
  else      return(0);
}


//------------ built-in, user-defined contact model "ductile" -----------
void CM3::FDlaw(FdBlock& fb, char cDim) {
	if (cDim==2) {  // ... 2D Section ...
    if (fb.u_n > 0.0) {          // COMPRESSIVE
      fb.n_force = kn * fb.u_n;
      double kst    = kn * ratks * fb.tdel;
      fb.s_force = fb.s_force - fb.u_dot_s * kst;
      if (broken) {
        double max_s_force = fb.n_force * fric; // Slip check
        double sfmag = fabs(fb.s_force);
        if (sfmag) {
          if (sfmag > max_s_force) {
            double rat    = max_s_force / sfmag;
            fb.s_force = fb.s_force * rat;
						fb.bSliding = true;
					}
				}
				else {
          fb.s_force = 0.0;
				}
			}
      fb.knest = kn;
      fb.ksest = kn * ratks;
      fb.skip  = false;
		}
		else {                     // EXTENSILE
      if (broken) {
        fb.n_force = 0.0;
        fb.s_force = 0.0;
			}
			else {
        double Fn =  kun * fb.u_n;
        if (Fn < FT) {
          FT = (ut - fb.u_n) * kduc - FTmax;
          if (FT > 0.0) {
            FT       = 0.0;
            broken   = 1.0;
            softened = 1.0;
            SetDelete(true);
            fb.bfishc_brokenn = true;
					}
          kun = FT / fb.u_n;
          Fn = FT;
          if (FTmax) {
            softened = 1.0 + FT / FTmax;
					}
					else {
            softened = 1.0;
					}
				}
        fb.n_force = Fn;
        if (broken == 0.0) {
          double kst         = kun * ratks * fb.tdel;
          fb.s_force      = fb.s_force - fb.u_dot_s * kst;
          double max_s_force = fabs(FT) * ratfs;
          double sfmag       = fabs(fb.s_force);
          if (sfmag) {
            if (sfmag > max_s_force) {
              double rat    = max_s_force / sfmag;
              fb.s_force = fb.s_force * rat;
							fb.bSliding = true;
						}
					}
					else {
            fb.s_force = 0.0;
					}
				}
				else {
          fb.s_force = 0.0;
				}
			}
      fb.knest = kun;
      fb.ksest = kun * ratks;
      fb.skip  = false;
		}
	}
	else {  // ... 3D section ...
    if (fb.u_n > 0.0) {          // COMPRESSIVE
      fb.n_force = kn * fb.u_n;
      double kst    = kn * ratks * fb.tdel;
      UMdvect vec     = fb_u_dot_s * kst;
      fb_s_force= fb_s_force- vec;
      if (broken) {
        double max_s_force = fb.n_force * fric; // Slip check
        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 (sfmag) {
          if (sfmag > max_s_force) {
            double rat    = max_s_force / sfmag;
            fb_s_force= fb_s_force* rat;
						fb.bSliding = true;
					}
				}
				else {
          fb_s_force.Fill(0.0);
				}
			}
      fb.knest = kn;
      fb.ksest = kn * ratks;
      fb.skip  = false;
		}
		else {                     // EXTENSILE
      if (broken) {
        fb.n_force = 0.0;
        fb_s_force.Fill(0.0);
			}
			else {
        double Fn =  kun * fb.u_n;
        if (Fn < FT) {
          FT = (ut - fb.u_n) * kduc - FTmax;
          if (FT > 0.0) {
            FT       = 0.0;
            broken   = 1.0;
            softened = 1.0;
            SetDelete(true);
            fb.bfishc_brokenn = true;
					}
          kun = FT / fb.u_n;
          Fn = FT;
          if (FTmax) {
            softened = 1.0 + FT / FTmax;
					}
					else {
            softened = 1.0;
					}
				}
        fb.n_force = Fn;
        if (broken == 0.0) {
          double kst         = kun * ratks * fb.tdel;
          UMdvect vec        = fb_u_dot_s * kst;
          fb_s_force     = fb_s_force- vec;
          double max_s_force = fabs(FT) * ratfs;
          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 (sfmag) {
            if (sfmag > max_s_force) {
              double rat    = max_s_force / sfmag;
              fb_s_force= fb_s_force* rat;
							fb.bSliding = true;
						}
					}
					else {
            fb_s_force.Fill(0.0);
					}
				}
				else {
          fb_s_force.Fill(0.0);
				}
			}
      fb.knest = kun;
      fb.ksest = kun * ratks;
      fb.skip  = false;
		}
  }
}


const char *CM3::SaveRestore(ModelSaveObject *mso) {
  const char *str = ContactModel::SaveRestore(mso);
  if (str) return(str);
  mso->Initialize(11,0,1);
  mso->Save( 0, delete_flag);
  mso->Save( 0, kn);
  mso->Save( 1, ratks);
  mso->Save( 2, ratfs);
  mso->Save( 3, fric);
  mso->Save( 4, FTmax);
  mso->Save( 5, kduc);
  mso->Save( 6, kun);
  mso->Save( 7, broken);
  mso->Save( 8, softened);
  mso->Save( 9, ut);
  mso->Save(10, FT);
  return(0);
}



/* EoF */
