/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

#include "stdlib.h"
#include "string.h"
#include "fix_wall_bounceback.h"
#include "atom.h"
#include "comm.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "lattice.h"
#include "input.h"
#include "variable.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "neigh_list.h"
#include "force.h"
#include "math.h"
#include "math_extra.h"
#include "math_vector.h"
#include "atom_vec_tri.h"
#include "memory.h"
#include "reader.h"
#include "error.h"



using namespace LAMMPS_NS;
using namespace FixConst;
enum{MOVE};

/* ---------------------------------------------------------------------- */

FixWallBounceback::FixWallBounceback(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg)
{
  if (narg < 3) error->all(FLERR,"Illegal fix wall/bounceback command");   // fix ID  group-ID  wall/bounceback  keyword(move)

  // parse args

  nwall = 0;
  int scaleflag = 0;    // 判断默认单位是box还是lattice 。0为box，1为lattice。我的文件中，细胞坐标以box为单位。

    int iarg= 2 ;  // move 关键词
  while (iarg < narg) 
  {
    if (strcmp(arg[iarg],"move") == 0) 
	{
	 if (iarg+2 > narg) error->all(FLERR,"Illegal fix wall/bounceback command");
    
      int newwall;
	  if (strcmp(arg[iarg],"move") == 0) newwall = MOVE;

     nwall++;
     iarg += 2;
    }
	else if (strcmp(arg[iarg],"units") == 0) 
	{
      if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect command");
      if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
      else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
      else error->all(FLERR,"Illegal fix wall/reflect command");
      iarg += 2;
    } 
	else error->all(FLERR,"Illegal fix wall/reflect command");	
  }

   // error check

  if (nwall == 0) error->all(FLERR,"Illegal fix wall command");

  // scale factors for walls
   for (int m = 0; m < nwall; m++)
   	
   	if (scaleflag) 
    {   // units 为lattice类型
      xscale = domain->lattice->xlattice;
      yscale = domain->lattice->ylattice;
      zscale = domain->lattice->zlattice;
    }
    else xscale = yscale = zscale = 1.0;    //units为box类型
 }
  /* ---------------------------------------------------------------------- */
  
 /* ---------------------------------------------------------------------- */
  
/* ---------------------------------------------------------------------- */

FixWallBounceback::~FixWallBounceback()
 {
 }

/* ---------------------------------------------------------------------- */

int FixWallBounceback::setmask()
{
  int mask = 0;
  mask |= POST_INTEGRATE;
  return mask;
}

/* ---------------------------------------------------------------------- */

void FixWallBounceback::init()
{
	neighbor->request(this);

}
/* ---------------------------------------------------------------------- */

void FixWallBounceback::init_list(int id, NeighList *ptr)
{
  list = ptr;
}

/* ---------------------------------------------------------------------- */


 void FixWallBounceback::compute()
{
  int i,j,ii,jj,inum,jnum,itype,jtype;
  int *ilist,*jlist,*numneigh,**firstneigh;
  double xtmp,ytmp,ztmp;
  double vxtmp,vytmp,vztmp;

  int atom1,atom2,atom3;
  int m,n;
  double x1tmp,y1tmp,z1tmp,x2tmp,y2tmp,z2tmp,x3tmp,y3tmp,z3tmp,pxtmp,pytmp,pztmp;
  double x1,y1,z1,x2,y2,z2,x3,y3,z3,px,py,pz;
  double vx1tmp,vy1tmp,vz1tmp,vx2tmp,vy2tmp,vz2tmp,vx3tmp,vy3tmp,vz3tmp,pvxtmp,pvytmp,pvztmp;
  double a1tmp[3],a2tmp[3],a1[3],a2[3],np[3],nptmp[3],ctmp[3],c[3],d1[3],d2[3],e1[3],q1[3],a1tmpcd2[3],d1ca2tmp[3],d1cd2[3];
  double a1new[3],a2new[3],vtri[3];
  double btmp,bnew,m0,m1,m2,m3,m4,m5,ft,dft,m6,m7,m8,m9,m10;;
  double t,t0,tol;
  double pnew[3],g[3],x1new,y1new,z1new,ksi,zita;
  int nt;
  double *ptr;
 
  double dt = update->dt;
  double **x = atom->x;
  double **v = atom->v;
  double **f = atom->f;
  int *type = atom->type;
  int **elementlist = neighbor->elementlist;
  int nelementlist = neighbor->nelementlist;
  int nlocal = atom->nlocal;

   int *num_element = atom->num_element;
   int **element_atom1 = atom->element_atom1;
   int **element_atom2 = atom->element_atom2;
   int **element_atom3 = atom->element_atom3;
   int **element_type = atom->element_type;
   int *tag = atom->tag;
   int *mask = atom->mask;
	 

  inum = list->inum;
  ilist = list->ilist;
  numneigh = list->numneigh;
  firstneigh = list->firstneigh;

  m = n = 0;
  for (atom2 = 0; atom2 < nlocal; atom2++) {
  	if (!(mask[atom2] & groupbit)) continue;
    for (i = 0; i < num_element[atom2]; i++) 
     { if (tag[atom2] != element_atom2[atom2][i]) continue;
      atom1 = atom->map(element_atom1[atom2][i]);
      if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
      atom3 = atom->map(element_atom3[atom2][i]);
      if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
      if (element_type[atom2][i] == 0) continue;

    //  if (flag) {
        
          x1tmp = x[atom1][0];
		  y1tmp = x[atom1][1];
		  z1tmp = x[atom1][2];
		  vx1tmp = v[atom1][0];
		  vy1tmp = v[atom1][1];
		  vz1tmp = v[atom1][2];
		  
		  x2tmp = x[atom2][0];
		  y2tmp = x[atom2][1];
		  z2tmp = x[atom2][2];
		  vx2tmp = v[atom2][0];
		  vy2tmp = v[atom2][1];
		  vz2tmp = v[atom2][2];
		  
		  x3tmp = x[atom3][0];
		  y3tmp = x[atom3][1];
		  z3tmp = x[atom3][2];
		  vx3tmp = v[atom3][0];
		  vy3tmp = v[atom3][1];
		  vz3tmp = v[atom3][2];
   
          x1 = x1tmp + vx1tmp*dt;
          y1 = y1tmp + vy1tmp*dt;
		  z1 = z1tmp + vz1tmp*dt;

		  x2 = x2tmp + vx2tmp*dt;
          y2 = y2tmp + vy2tmp*dt;
		  z2 = z2tmp + vz2tmp*dt;

		  x3 = x1tmp + vx3tmp*dt;
          y3 = y1tmp + vy3tmp*dt;
		  z3 = z1tmp + vz3tmp*dt;

			
          a1tmp[0] = x2tmp - x1tmp;
          a1tmp[1] = y2tmp - y1tmp;
          a1tmp[2] = z2tmp - z1tmp;
     
          a2tmp[0] = x3tmp - x1tmp;
          a2tmp[1] = y3tmp - y1tmp;
          a2tmp[2] = z3tmp - z1tmp;

		  d1[0] = vx2tmp - vx1tmp;
		  d1[1] = vy2tmp - vy1tmp;
		  d1[2] = vz2tmp - vz1tmp;

		  d2[0] = vx3tmp - vx1tmp;
		  d2[1] = vy3tmp - vy1tmp;
	      d2[2] = vz3tmp - vz1tmp;
		  	
	      a1[0] = x2 - x1;
		  a1[1] = y2 - y1;
		  a1[2] = z2 - z1;

		  a2[0] = x3 - x1;
		  a2[1] = y3 - y1;
		  a2[2] = z3 - z1;
   
        MathExtra::cross3(a1,a2,np);  // normal vector

        }
     }

  // loop over neighbors of my atoms

  for (ii = 0; ii < inum; ii++)
 {
    i = ilist[ii];
    xtmp = x[i][0];   // xtmp, ytmp, ztmp  实际上没用到
    ytmp = x[i][1];
    ztmp = x[i][2];
    vxtmp = v[i][0];
    vytmp = v[i][1];
    vztmp = v[i][2];
    itype = type[i];
    jlist = firstneigh[i];
    jnum = numneigh[i];


    x[i][0] += vxtmp*dt;
	x[i][1] += vytmp*dt;
	x[i][2] += vztmp*dt;
	
   for (jj = 0; jj < jnum; jj++)   // j 就是运动的点，判断j类原子和i类组成的三角形的位置关系
  {
      j = jlist[jj];
      j &= NEIGHMASK;
	  jtype = type[j];

     pxtmp = x[j][0];
	 pytmp = x[j][1];
	 pztmp = x[j][2];
	 pvxtmp = v[j][0];
	 pvytmp = v[j][1];
	 pvztmp = v[j][2];

	 px = pxtmp + pvxtmp*dt;
	 py = pytmp + pvytmp*dt;
	 pz = pztmp + pvztmp*dt;

     ctmp[0] = pxtmp - x1tmp;
	 ctmp[1] = pytmp - y1tmp;
	 ctmp[2] = pztmp - z1tmp;

	 e1[0] = pvxtmp - vx1tmp;
	 e1[1] = pvytmp - vy1tmp;
	 e1[2] = pvztmp - vz1tmp;

	 c[0] = px - x1;
	 c[1] = py - y1;
	 c[2] = pz - z1;

      MathExtra::cross3(a1tmp,a2tmp,nptmp);
	  btmp = MathExtra::dot3(nptmp,ctmp);
	  bnew = MathExtra::dot3(np,c);
      MathExtra::cross3(a1tmp,d2,a1tmpcd2);
	  MathExtra::cross3(d1,a2tmp,d1ca2tmp);
	  m0 = MathExtra::dot3(nptmp,ctmp);
	  m1 = MathExtra::dot3(nptmp,e1);
	  MathExtra::cross3(a1tmp,d2,a1tmpcd2);
	  MathExtra::cross3(d1,a2tmp,d1ca2tmp);
	  MathExtra::add3(a1tmpcd2,d1ca2tmp,q1);
      m2 = MathExtra::dot3(q1,ctmp);
	  m3 = MathExtra::dot3(q1,e1);
	  MathExtra::cross3(d1,d2,d1cd2);
	  m4 = MathExtra::dot3(d1cd2,ctmp);
	  m5 = MathExtra::dot3(d1cd2,e1);


	 double t=0,t0=0,tol=0.0001;
     int nt=20; 
	  if (!( btmp*bnew > 0 )) 
	 {
	    for(i=0;i<nt;i++)
	    {
	     ft = m5*t0*t0*t0+(m3+m4)*t0*t0+(m1+m2)*t0+m0;
	     dft = m5*3.0*t0*t0+(m3+m4)*2.0*t0+m1+m2;
		 t=t0-ft/dft;
		 if (fabs(t-t0)<=tol)
		  break;
		 t0=t;
        }
        pnew[0] = pxtmp + pvxtmp*t;
		pnew[1] = pytmp + pvytmp*t;
		pnew[2] = pztmp + pvztmp*t;

		x1new = x1tmp + vx1tmp*t;
		y1new = y1tmp + vy1tmp*t;
		z1new = z1tmp + vz1tmp*t;
		
        g[0] = pnew[0] - x1new;
		g[1] = pnew[1] - y1new;
		g[2] = pnew[2] - z1new;

		a1new[0] = a1tmp[0] + d1[0]*t;
		a1new[1] = a1tmp[1] + d1[1]*t;
		a1new[2] = a1tmp[2] + d1[2]*t;

        a2new[0] = a2tmp[0] + d2[0]*t;
		a2new[1] = a2tmp[1] + d2[1]*t;
		a2new[2] = a2tmp[2] + d2[2]*t;

		m6 = MathExtra::dot3(a1new,a2new);
		m7 = MathExtra::dot3(g,a1new);
		m8 = MathExtra::dot3(g,a2new);
		m9 = MathExtra::lensq3(a1new);
		m10 = MathExtra::lensq3(a2new);

		ksi = (m7*m10-m8*m6)/(m9*m10-m6*m6);
		zita = (m8*m9-m7*m6)/(m9*m10-m6*m6);
		if(ksi>=0 && zita>=0 && (ksi+zita)<=1 )
	      {
	       vtri[0] = (1-ksi-zita)*vx1tmp + ksi*vx2tmp + zita*vx3tmp;
           vtri[1] = (1-ksi-zita)*vy1tmp + ksi*vy2tmp + zita*vy3tmp;
		   vtri[2] = (1-ksi-zita)*vz1tmp + ksi*vz2tmp + zita*vz3tmp;

		   pvxtmp = 2*vtri[0] - pvxtmp;
		   pvytmp = 2*vtri[1] - pvytmp;
		   pvztmp = 2*vtri[2] - pvztmp;

           px = pnew[0] + (dt - t)*pvxtmp;
		   py = pnew[1] + (dt - t)*pvytmp;
		   pz = pnew[2] + (dt - t)*pvztmp;
		  }
	     }
     	}
      }
	 }
   
/* ---------------------------------------------------------------------- */

void FixWallBounceback::setup(int)
{
}



/* ---------------------------------------------------------------------- */
void FixWallBounceback::min_setup(int)
{
}


/* ---------------------------------------------------------------------- */
//double FixWallBounceback::compute_scalar()
//{
//}


/* ---------------------------------------------------------------------- */

//double FixWallBounceback::compute_vector()
//{
//}

/* ---------------------------------------------------------------------- */



/* ---------------------------------------------------------------------- */
/*
void FixWallBounceback::post_integrate()
{
  int i,dim,side;
  double coord;

  // coord = current position of wall
  // evaluate variable if necessary, wrap with clear/add

  double **x = atom->x;
  double **v = atom->v;
  int *mask = atom->mask;
  int nlocal = atom->nlocal;

  if (moveflag) modify->clearstep_compute();  //  ???????

  for (int m = 0; m < nwall; m++) 
   {
   /* if (wallstyle[m] == VARIABLE) 
	 {
      coord = input->variable->compute_equal(varindex[m]);
      if (wallwhich[m] < YLO) coord *= xscale;
      else if (wallwhich[m] < ZLO) coord *= yscale;
      else coord *= zscale;
     }
	else coord = coord0[m];
   */
   // dim = wallmove[m] / 2;
/*    side = wallnormalside[m] % 2;  //  positive normal=0   and negative normal side =1

    for (i = 0; i < nlocal; i++)
      if (mask[i] & groupbit) 
	   {
        if (side == 0) 
		{
          if (x[i][dim] < coord) 
		  {
            x[i][dim] = coord + (coord - x[i][dim]);
            v[i][dim] = -v[i][dim];
          }
        } 
		else 
		{
          if (x[i][dim] > coord) 
		  {
            x[i][dim] = coord - (x[i][dim] - coord);
            v[i][dim] = -v[i][dim];
          }
        }
      }
  }

  if (moveflag) modify->addstep_compute(update->ntimestep + 1);   
}

*/
/* ---------------------------------------------------------------------- */

void FixWallBounceback::post_force(int)
{
}


/* ---------------------------------------------------------------------- */