Commit cd8d18dc authored by sjplimp's avatar sjplimp
Browse files 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15719 f3b2605a-c512-4ea7-a41b-209d697bcdaa
parent 5bc562b0

src/GRANULAR/fix_wall_gran.cpp

+450 −151
Original line number Diff line number Diff line
@@ -12,7 +12,8 @@ 
------------------------------------------------------------------------- */



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

   Contributing authors: Leo Silbert (SNL), Gary Grest (SNL)

   Contributing authors: Leo Silbert (SNL), Gary Grest (SNL),

                         Dan Bolintineanu (SNL)

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



#include <math.h>
@@ -34,8 +35,11 @@ using namespace LAMMPS_NS; 
using namespace FixConst;

using namespace MathConst;



enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER};    // XYZ PLANE need to be 0,1,2

enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY};

// XYZ PLANE need to be 0,1,2



enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION}; 

enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,BONDED_HISTORY};

enum{NONE,CONSTANT,EQUAL};



#define BIG 1.0e20
@@ -44,7 +48,7 @@ enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY}; 
FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :

  Fix(lmp, narg, arg)

{

  if (narg < 11) error->all(FLERR,"Illegal fix wall/gran command");

  if (narg < 4) error->all(FLERR,"Illegal fix wall/gran command");



  if (!atom->sphere_flag)

    error->all(FLERR,"Fix wall/gran requires atom style sphere");
@@ -53,16 +57,23 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : 
  create_attribute = 1;



  // set interaction style

  // disable bonded/history option for now



  if (strcmp(arg[3],"hooke") == 0) pairstyle = HOOKE;

  else if (strcmp(arg[3],"hooke/history") == 0) pairstyle = HOOKE_HISTORY;

  else if (strcmp(arg[3],"hertz/history") == 0) pairstyle = HERTZ_HISTORY;

  //else if (strcmp(arg[3],"bonded/history") == 0) pairstyle = BONDED_HISTORY;

  else error->all(FLERR,"Invalid fix wall/gran interaction style");



  history = 1;

  if (pairstyle == HOOKE) history = 0;



  // particle/wall coefficients

  // wall/particle coefficients



  int iarg;



  if (pairstyle != BONDED_HISTORY) {

    if (narg < 11) error->all(FLERR,"Illegal fix wall/gran command");



    kn = force->numeric(FLERR,arg[4]);

    if (strcmp(arg[5],"NULL") == 0) kt = kn * 2.0/7.0;
@@ -87,9 +98,30 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : 
      kt /= force->nktv2p;

    }



    iarg = 10;

  }



  else {

    if (narg < 10) error->all(FLERR,"Illegal fix wall/gran command");



    E = force->numeric(FLERR,arg[4]);

    G = force->numeric(FLERR,arg[5]);

    SurfEnergy = force->numeric(FLERR,arg[6]);

    // Note: this doesn't get used, check w/ Jeremy?

    gamman = force->numeric(FLERR,arg[7]);



    xmu = force->numeric(FLERR,arg[8]);

    // pois = E/(2.0*G) - 1.0;

    // kn = 2.0*E/(3.0*(1.0+pois)*(1.0-pois));

    // gammat=0.5*gamman;



    iarg = 9;

  }



  // wallstyle args

  

  int iarg = 10;

  idregion = NULL;



  if (strcmp(arg[iarg],"xplane") == 0) {

    if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/gran command");

    wallstyle = XPLANE;
@@ -118,11 +150,18 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : 
    if (narg < iarg+2) error->all(FLERR,"Illegal fix wall/gran command");

    wallstyle = ZCYLINDER;

    lo = hi = 0.0;

    cylradius = force->numeric(FLERR,arg[iarg+1]);

    cylradius = force->numeric(FLERR,arg[iarg+3]);

    iarg += 2;

  } else if (strcmp(arg[iarg],"region") == 0) {

    if (narg < iarg+2) error->all(FLERR,"Illegal fix wall/gran command");

    wallstyle = REGION;

    int n = strlen(arg[iarg+1]) + 1;

    idregion = new char[n];

    strcpy(idregion,arg[iarg+1]);

    iarg += 2;

  }

  

  // check for trailing keyword/values

  // optional args



  wiggle = 0;

  wshear = 0;
@@ -169,6 +208,8 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : 
    error->all(FLERR,"Invalid shear direction for fix wall/gran");

  if (wshear && wallstyle == ZPLANE && axis == 2)

    error->all(FLERR,"Invalid shear direction for fix wall/gran");

  if ((wiggle || wshear) && wallstyle == REGION)

    error->all(FLERR,"Cannot wiggle or shear with fix wall/gran/region");



  // setup oscillations
@@ -177,7 +218,10 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : 
  // perform initial allocation of atom-based arrays

  // register with Atom class



  shear = NULL;

  if (pairstyle == BONDED_HISTORY) sheardim = 7;

  else sheardim = 3;

  

  shearone = NULL;

  grow_arrays(atom->nmax);

  atom->add_callback(0);

  atom->add_callback(1);
@@ -185,12 +229,14 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : 
  nmax = 0;

  mass_rigid = NULL;



  // initialize as if particle is not touching wall

  // initialize shear history as if particle is not touching region

  // shearone will be NULL for wallstyle = REGION



  if (history) {

  if (history && shearone) {

    int nlocal = atom->nlocal;

    for (int i = 0; i < nlocal; i++)

      shear[i][0] = shear[i][1] = shear[i][2] = 0.0;

      for (int j = 0; j < sheardim; j++)

        shearone[i][j] = 0.0;

  }



  time_origin = update->ntimestep;
@@ -205,9 +251,10 @@ FixWallGran::~FixWallGran() 
  atom->delete_callback(id,0);

  atom->delete_callback(id,1);



  // delete locally stored arrays

  // delete local storage



  memory->destroy(shear);

  delete [] idregion;

  memory->destroy(shearone);

  memory->destroy(mass_rigid);

}
@@ -257,8 +304,8 @@ void FixWallGran::setup(int vflag) 


void FixWallGran::post_force(int vflag)

{

  int i;

  double dx,dy,dz,del1,del2,delxy,delr,rsq,meff;

  int i,j;

  double dx,dy,dz,del1,del2,delxy,delr,rsq,rwall,meff;

  double vwall[3];



  // do not update shear history during setup
@@ -322,6 +369,8 @@ void FixWallGran::post_force(int vflag) 
  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  

  rwall = 0.0;



  for (int i = 0; i < nlocal; i++) {

    if (mask[i] & groupbit) {
@@ -345,13 +394,17 @@ void FixWallGran::post_force(int vflag) 
      } else if (wallstyle == ZCYLINDER) {

        delxy = sqrt(x[i][0]*x[i][0] + x[i][1]*x[i][1]);

        delr = cylradius - delxy;

        if (delr > radius[i]) dz = cylradius;

        else {

        if (delr > radius[i]) {

          dz = cylradius;

          rwall = 0.0;

        } else {

          dx = -delr/delxy * x[i][0];

          dy = -delr/delxy * x[i][1];

          // rwall = -2r_c if inside cylinder, 2r_c outside

          rwall = 2*(1-2*(delxy < cylradius))*cylradius; 

          if (wshear && axis != 2) {

            vwall[0] = vshear * x[i][1]/delxy;

            vwall[1] = -vshear * x[i][0]/delxy;

            vwall[0] += vshear * x[i][1]/delxy;

            vwall[1] += -vshear * x[i][0]/delxy;

            vwall[2] = 0.0;

          }

        }
@@ -360,11 +413,10 @@ void FixWallGran::post_force(int vflag) 
      rsq = dx*dx + dy*dy + dz*dz;

      

      if (rsq > radius[i]*radius[i]) {

        if (pairstyle != HOOKE) {

          shear[i][0] = 0.0;

          shear[i][1] = 0.0;

          shear[i][2] = 0.0;

        }

        if (history)

          for (j = 0; j < sheardim; j++) 

            shearone[i][j] = 0.0;



      } else {



        // meff = effective mass of sphere
@@ -373,17 +425,20 @@ void FixWallGran::post_force(int vflag) 
        meff = rmass[i];

        if (fix_rigid && mass_rigid[i] > 0.0) meff = mass_rigid[i];



        // inovke sphere/wall interaction

        // invoke sphere/wall interaction



        if (pairstyle == HOOKE)

          hooke(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],

                radius[i],meff);

          hooke(rsq,dx,dy,dz,vwall,v[i],f[i],

                omega[i],torque[i],radius[i],meff);

        else if (pairstyle == HOOKE_HISTORY)

          hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],

                        radius[i],meff,shear[i]);

          hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],

                        omega[i],torque[i],radius[i],meff,shearone[i]);

        else if (pairstyle == HERTZ_HISTORY)

          hertz_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],

                        radius[i],meff,shear[i]);

          hertz_history(rsq,dx,dy,dz,vwall,rwall,v[i],f[i],

                        omega[i],torque[i],radius[i],meff,shearone[i]);

        else if (pairstyle == BONDED_HISTORY)

          bonded_history(rsq,dx,dy,dz,vwall,rwall,v[i],f[i],

                         omega[i],torque[i],radius[i],meff,shearone[i]);

      }

    }

  }
@@ -599,7 +654,7 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, 
/* ---------------------------------------------------------------------- */



void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,

                                double *vwall, double *v,

                                double *vwall, double rwall, double *v,

                                double *f, double *omega, double *torque,

                                double radius, double meff, double *shear)

{
@@ -638,10 +693,15 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, 
  wr3 = radius*omega[2] * rinv;

  

  // normal forces = Hertzian contact + normal velocity damping

  // rwall = 0 is flat wall case

  // rwall positive or negative is curved wall

  //   will break (as it should) if rwall is negative and 

  //   its absolute value < radius of particle

  

  damp = meff*gamman*vnnr*rsqinv;

  ccel = kn*(radius-r)*rinv - damp;

  polyhertz = sqrt((radius-r)*radius);

  if (rwall == 0.0) polyhertz = sqrt((radius-r)*radius);

  else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius));

  ccel *= polyhertz;

  

  // relative velocities
@@ -714,6 +774,247 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, 
  torque[2] -= radius*tor3;

}





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



void FixWallGran::bonded_history(double rsq, double dx, double dy, double dz,

                                 double *vwall, double rwall, double *v,

                                 double *f, double *omega, double *torque,

                                 double radius, double meff, double *shear)

{

  double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;

  double wr1,wr2,wr3,damp,ccel,vtr1,vtr2,vtr3,vrel;

  double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3;

  double shrmag,rsht,polyhertz,rinv,rsqinv;



  double pois,E_eff,G_eff,rad_eff;

  double a0,Fcrit,delcrit,delcritinv;

  double overlap,olapsq,olapcubed,sqrtterm,tmp,keyterm,keyterm2,keyterm3;

  double aovera0,foverFc;

  double gammatsuji;



  double ktwist,kroll,twistcrit,rollcrit;

  double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv;

  double magtwist,magtortwist;

  double magrollsq,magroll,magrollinv,magtorroll;



  r = sqrt(rsq);

  rinv = 1.0/r;

  rsqinv = 1.0/rsq;



  // relative translational velocity



  vr1 = v[0] - vwall[0];

  vr2 = v[1] - vwall[1];

  vr3 = v[2] - vwall[2];



  // normal component



  vnnr = vr1*dx + vr2*dy + vr3*dz;

  vn1 = dx*vnnr / rsq;

  vn2 = dy*vnnr / rsq;

  vn3 = dz*vnnr / rsq;



  // tangential component



  vt1 = vr1 - vn1;

  vt2 = vr2 - vn2;

  vt3 = vr3 - vn3;



  // relative rotational velocity



  wr1 = radius*omega[0] * rinv;

  wr2 = radius*omega[1] * rinv;

  wr3 = radius*omega[2] * rinv;



  // normal forces = Hertzian contact + normal velocity damping

  // material properties: currently assumes identical materials



  pois = E/(2.0*G) - 1.0;

  E_eff=0.5*E/(1.0-pois*pois);

  G_eff=G/(4.0-2.0*pois);



  // rwall = 0 is infinite wall radius of curvature (flat wall)



  if (rwall == 0) rad_eff = radius;

  else rad_eff = radius*rwall/(radius+rwall);



  Fcrit = rad_eff * (3.0 * M_PI * SurfEnergy);

  a0=pow(9.0*M_PI*SurfEnergy*rad_eff*rad_eff/E_eff,1.0/3.0);

  delcrit = 1.0/rad_eff*(0.5 * a0*a0/pow(6.0,1.0/3.0));

  delcritinv = 1.0/delcrit;



  overlap = (radius-r) * delcritinv;

  olapsq = overlap*overlap;

  olapcubed = olapsq*overlap;

  sqrtterm = sqrt(1.0 + olapcubed);

  tmp = 2.0 + olapcubed + 2.0*sqrtterm;

  keyterm = pow(tmp,THIRD);

  keyterm2 = olapsq/keyterm;

  keyterm3 = sqrt(overlap + keyterm2 + keyterm);

  aovera0 = pow(6.0,-TWOTHIRDS) * (keyterm3 +

            sqrt(2.0*overlap - keyterm2 - keyterm + 4.0/keyterm3));

  foverFc = 4.0*((aovera0*aovera0*aovera0) - pow(aovera0,1.5));

  ccel = Fcrit*foverFc*rinv;



  // damp = meff*gamman*vnnr*rsqinv;

  // ccel = kn*(radius-r)*rinv - damp;

  // polyhertz = sqrt((radius-r)*radius);

  // ccel *= polyhertz;



  // use Tsuji et al form



  polyhertz = 1.2728- 4.2783*0.9 + 11.087*0.9*0.9 - 22.348*0.9*0.9*0.9 + 

    27.467*0.9*0.9*0.9*0.9 - 18.022*0.9*0.9*0.9*0.9*0.9 + 

    4.8218*0.9*0.9*0.9*0.9*0.9*0.9;



  gammatsuji = 0.2*sqrt(meff*kn);

  damp = gammatsuji*vnnr/rsq;

  ccel = ccel - polyhertz * damp;



  // relative velocities



  vtr1 = vt1 - (dz*wr2-dy*wr3);

  vtr2 = vt2 - (dx*wr3-dz*wr1);

  vtr3 = vt3 - (dy*wr1-dx*wr2);

  vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;

  vrel = sqrt(vrel);



  // shear history effects



  if (shearupdate) {

    shear[0] += vtr1*dt;

    shear[1] += vtr2*dt;

    shear[2] += vtr3*dt;

  }

  shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]);



  // rotate shear displacements



  rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz;

  rsht = rsht*rsqinv;

  if (shearupdate) {

    shear[0] -= rsht*dx;

    shear[1] -= rsht*dy;

    shear[2] -= rsht*dz;

  }



  // tangential forces = shear + tangential velocity damping



  fs1 = -polyhertz * (kt*shear[0] + meff*gammat*vtr1);

  fs2 = -polyhertz * (kt*shear[1] + meff*gammat*vtr2);

  fs3 = -polyhertz * (kt*shear[2] + meff*gammat*vtr3);



  kt=8.0*G_eff*a0*aovera0;



  // shear damping uses Tsuji et al form also



  fs1 = -kt*shear[0] - polyhertz*gammatsuji*vtr1;

  fs2 = -kt*shear[1] - polyhertz*gammatsuji*vtr2;

  fs3 = -kt*shear[2] - polyhertz*gammatsuji*vtr3;



  // rescale frictional displacements and forces if needed



  fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);

  fn = xmu * fabs(ccel*r + 2.0*Fcrit);



  if (fs > fn) {

    if (shrmag != 0.0) {

      shear[0] = (fn/fs) * (shear[0] + polyhertz*gammatsuji*vtr1/kt) -

      polyhertz*gammatsuji*vtr1/kt;

      shear[1] = (fn/fs) * (shear[1] + polyhertz*gammatsuji*vtr2/kt) -

      polyhertz*gammatsuji*vtr2/kt;

      shear[2] = (fn/fs) * (shear[2] + polyhertz*gammatsuji*vtr3/kt) -

      polyhertz*gammatsuji*vtr3/kt;

      fs1 *= fn/fs ;

      fs2 *= fn/fs;

      fs3 *= fn/fs;

    } else fs1 = fs2 = fs3 = 0.0;

  }



  // calculate twisting and rolling components of torque

  // NOTE: this assumes spheres!



  relrot1 = omega[0];

  relrot2 = omega[1];

  relrot3 = omega[2];



  // rolling velocity

  // NOTE: this assumes mondisperse spheres!



  vrl1 = -rad_eff*rinv * (relrot2*dz - relrot3*dy);

  vrl2 = -rad_eff*rinv * (relrot3*dx - relrot1*dz);

  vrl3 = -rad_eff*rinv * (relrot1*dy - relrot2*dx);

  vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3);

  if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag;

  else vrlmaginv = 0.0;



  // bond history effects



  shear[3] += vrl1*dt;

  shear[4] += vrl2*dt;

  shear[5] += vrl3*dt;



  // rotate bonded displacements correctly



  double rlt = shear[3]*dx + shear[4]*dy + shear[5]*dz;

  rlt /= rsq;

  shear[3] -= rlt*dx;

  shear[4] -= rlt*dy;

  shear[5] -= rlt*dz;



  // twisting torque



  magtwist = rinv*(relrot1*dx + relrot2*dy + relrot3*dz);

  shear[6] += magtwist*dt;



  ktwist = 0.5*kt*(a0*aovera0)*(a0*aovera0);

  magtortwist = -ktwist*shear[6] - 

    0.5*polyhertz*gammatsuji*(a0*aovera0)*(a0*aovera0)*magtwist;



  twistcrit=TWOTHIRDS*a0*aovera0*Fcrit;

  if (fabs(magtortwist) > twistcrit)

    magtortwist = -twistcrit * magtwist/fabs(magtwist);



  // rolling torque



  magrollsq = shear[3]*shear[3] + shear[4]*shear[4] + shear[5]*shear[5];

  magroll = sqrt(magrollsq);

  if (magroll != 0.0) magrollinv = 1.0/magroll;

  else magrollinv = 0.0;



  kroll = 1.0*4.0*Fcrit*pow(aovera0,1.5);

  magtorroll = -kroll*magroll - 0.1*gammat*vrlmag;



  rollcrit = 0.01;

  if (magroll > rollcrit) magtorroll = -kroll*rollcrit;



  // forces & torques



  fx = dx*ccel + fs1;

  fy = dy*ccel + fs2;

  fz = dz*ccel + fs3;



  f[0] += fx;

  f[1] += fy;

  f[2] += fz;



  tor1 = rinv * (dy*fs3 - dz*fs2);

  tor2 = rinv * (dz*fs1 - dx*fs3);

  tor3 = rinv * (dx*fs2 - dy*fs1);

  torque[0] -= radius*tor1;

  torque[1] -= radius*tor2;

  torque[2] -= radius*tor3;



  torque[0] += magtortwist * dx*rinv;

  torque[1] += magtortwist * dy*rinv;

  torque[2] += magtortwist * dz*rinv;



  torque[0] += magtorroll * (shear[4]*dz - shear[5]*dy)*rinv*magrollinv;

  torque[1] += magtorroll * (shear[5]*dx - shear[3]*dz)*rinv*magrollinv;

  torque[2] += magtorroll * (shear[3]*dy - shear[4]*dx)*rinv*magrollinv;

}



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

   memory usage of local atom-based arrays

------------------------------------------------------------------------- */
@@ -722,7 +1023,7 @@ double FixWallGran::memory_usage() 
{

  int nmax = atom->nmax;

  double bytes = 0.0;

  if (history) bytes += 3*nmax * sizeof(double);   // shear history

  if (history) bytes += nmax*sheardim * sizeof(double);   // shear history

  if (fix_rigid) bytes += nmax * sizeof(int);             // mass_rigid

  return bytes;

}
@@ -733,7 +1034,7 @@ double FixWallGran::memory_usage() 


void FixWallGran::grow_arrays(int nmax)

{

  if (history) memory->grow(shear,nmax,3,"fix_wall_gran:shear");

  if (history) memory->grow(shearone,nmax,sheardim,"fix_wall_gran:shearone");

}



/* ----------------------------------------------------------------------
@@ -742,11 +1043,9 @@ void FixWallGran::grow_arrays(int nmax) 


void FixWallGran::copy_arrays(int i, int j, int delflag)

{

  if (history) {

    shear[j][0] = shear[i][0];

    shear[j][1] = shear[i][1];

    shear[j][2] = shear[i][2];

  }

  if (history)

    for (int m = 0; m < sheardim; m++)

      shearone[j][m] = shearone[i][m];

}



/* ----------------------------------------------------------------------
@@ -755,7 +1054,9 @@ void FixWallGran::copy_arrays(int i, int j, int delflag) 


void FixWallGran::set_arrays(int i)

{

  if (history) shear[i][0] = shear[i][1] = shear[i][2] = 0.0;

  if (history)

    for (int m = 0; m < sheardim; m++)

      shearone[i][m] = 0;

}



/* ----------------------------------------------------------------------
@@ -766,10 +1067,10 @@ int FixWallGran::pack_exchange(int i, double *buf) 
{

  if (!history) return 0;



  buf[0] = shear[i][0];

  buf[1] = shear[i][1];

  buf[2] = shear[i][2];

  return 3;

  int n = 0;

  for (int m = 0; m < sheardim; m++)

    buf[n++] = shearone[i][m];

  return n;

}



/* ----------------------------------------------------------------------
@@ -780,10 +1081,10 @@ int FixWallGran::unpack_exchange(int nlocal, double *buf) 
{

  if (!history) return 0;



  shear[nlocal][0] = buf[0];

  shear[nlocal][1] = buf[1];

  shear[nlocal][2] = buf[2];

  return 3;

  int n = 0;

  for (int m = 0; m < sheardim; m++)

    shearone[nlocal][m] = buf[n++];

  return n;

}



/* ----------------------------------------------------------------------
@@ -794,12 +1095,11 @@ int FixWallGran::pack_restart(int i, double *buf) 
{

  if (!history) return 0;



  int m = 0;

  buf[m++] = 4;

  buf[m++] = shear[i][0];

  buf[m++] = shear[i][1];

  buf[m++] = shear[i][2];

  return m;

  int n = 0;

  buf[n++] = sheardim + 1;

  for (int m = 0; m < sheardim; m++)

    buf[n++] = shearone[i][m];

  return n;

}



/* ----------------------------------------------------------------------
@@ -808,19 +1108,18 @@ int FixWallGran::pack_restart(int i, double *buf) 


void FixWallGran::unpack_restart(int nlocal, int nth)

{

  double **extra = atom->extra;



  if (!history) return;



  double **extra = atom->extra;



  // skip to Nth set of extra values

  

  int m = 0;

  for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);

  m++;

  

  shear[nlocal][0] = extra[nlocal][m++];

  shear[nlocal][1] = extra[nlocal][m++];

  shear[nlocal][2] = extra[nlocal][m++];

  for (int i = 0; i < sheardim; i++)

    shearone[nlocal][i] = extra[nlocal][m++];

}



/* ----------------------------------------------------------------------
@@ -830,7 +1129,7 @@ void FixWallGran::unpack_restart(int nlocal, int nth) 
int FixWallGran::maxsize_restart()

{

  if (!history) return 0;

  return 4;

  return 1 + sheardim;

}



/* ----------------------------------------------------------------------
@@ -840,7 +1139,7 @@ int FixWallGran::maxsize_restart() 
int FixWallGran::size_restart(int nlocal)

{

  if (!history) return 0;

  return 4;

  return 1 + sheardim;

}



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

src/GRANULAR/fix_wall_gran.h

+34 −26
Original line number Diff line number Diff line
@@ -29,51 +29,59 @@ class FixWallGran : public Fix { 
  FixWallGran(class LAMMPS *, int, char **);

  virtual ~FixWallGran();

  int setmask();

  void init();

  virtual void init();

  void setup(int);

  virtual void post_force(int);

  virtual void post_force_respa(int, int, int);



  double memory_usage();

  void grow_arrays(int);

  void copy_arrays(int, int, int);

  void set_arrays(int);

  int pack_exchange(int, double *);

  int unpack_exchange(int, double *);

  int pack_restart(int, double *);

  void unpack_restart(int, int);

  int size_restart(int);

  int maxsize_restart();

  virtual double memory_usage();

  virtual void grow_arrays(int);

  virtual void copy_arrays(int, int, int);

  virtual void set_arrays(int);

  virtual int pack_exchange(int, double *);

  virtual int unpack_exchange(int, double *);

  virtual int pack_restart(int, double *);

  virtual void unpack_restart(int, int);

  virtual int size_restart(int);

  virtual int maxsize_restart();

  void reset_dt();



  void hooke(double, double, double, double, double *,

             double *, double *, double *, double *, double, double);

  void hooke_history(double, double, double, double, double *,

                     double *, double *, double *, double *, double, double,

                     double *);

  void hertz_history(double, double, double, double, double *, double,

                     double *, double *, double *, double *, double, double,

                     double *);

  void bonded_history(double, double, double, double, double *, double,

                       double *, double *, double *, double *, double, double,

                       double *);



 protected:

  int wallstyle,pairstyle,history,wiggle,wshear,axis;

  int wallstyle,wiggle,wshear,axis;

  int pairstyle,nlevels_respa;

  bigint time_origin;

  double kn,kt,gamman,gammat,xmu;

  double E,G,SurfEnergy;

  double lo,hi,cylradius;

  double amplitude,period,omega,vshear;

  double dt;

  int nlevels_respa;

  int time_origin;

  char *idregion;



  // shear history values

  int history;       // if particle/wall interaction stores history

  int shearupdate;   // flag for whether shear history is updated

  int sheardim;      // # of shear history values per contact



  double **shear;

  int shearupdate;

  // shear history for single contact per particle



  double **shearone;



  // rigid body masses for use in granular interactions



  class Fix *fix_rigid;    // ptr to rigid body fix, NULL if none

  double *mass_rigid;      // rigid mass for owned+ghost atoms

  int nmax;                // allocated size of mass_rigid



  void hooke(double, double, double, double, double *,

             double *, double *, double *, double *, double, double);

  void hooke_history(double, double, double, double, double *,

                     double *, double *, double *, double *, double, double,

                     double *);

  void hertz_history(double, double, double, double, double *,

                     double *, double *, double *, double *, double, double,

                     double *);

};



}

src/GRANULAR/fix_wall_gran_region.cpp

0 → 100644
+548 −0

File added.

Preview size limit exceeded, changes collapsed.

src/GRANULAR/fix_wall_gran_region.h

0 → 100644
+106 −0
Original line number Diff line number Diff line 
/* -*- c++ -*- ----------------------------------------------------------

   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.

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



#ifdef FIX_CLASS



FixStyle(wall/gran/region,FixWallGranRegion)



#else



#ifndef LMP_FIX_WALL_GRAN_REGION_H

#define LMP_FIX_WALL_GRAN_REGION_H



#include "fix_wall_gran.h"



namespace LAMMPS_NS {



class FixWallGranRegion : public FixWallGran {

 public:

  FixWallGranRegion(class LAMMPS *, int, char **);

  ~FixWallGranRegion();

  void post_force(int);

  void write_restart(FILE *);

  void restart(char* );

  void init();

    

  double memory_usage();

  void grow_arrays(int);

  void copy_arrays(int, int, int);

  void set_arrays(int);

  int pack_exchange(int, double *);

  int unpack_exchange(int, double *);

  int pack_restart(int, double *);

  void unpack_restart(int, int);

  int size_restart(int);

  int maxsize_restart();



 private:

  class Region *region;

  char *region_style;

  int nregion;

  

  // shear history for multiple contacts per particle



  int tmax;              // max # of region walls one particle can touch 

  int *ncontact;         // # of shear contacts per particle

  int **walls;           // which wall each contact is with

  double ***shearmany;   // shear history per particle per contact

  int *c2r;              // contact to region mapping

                         // c2r[i] = index of Ith contact in

                         //   region-contact[] list of contacts

  int motion_resetflag;  // used by restart to indicate that region 

                         //    vel info is to be reset



  void update_contacts(int, int);

  int check_consistent_region(Region *, char*, char*, int);

};



}



#endif

#endif



/* ERROR/WARNING messages:



E: Illegal ... command



Self-explanatory.  Check the input script syntax and compare to the

documentation for the command.  You can use -echo screen as a

command-line option when running LAMMPS to see the offending line.



E: Fix wall/gran requires atom style sphere



Self-explanatory.



E: Cannot use wall in periodic dimension



Self-explanatory.



E: Cannot wiggle and shear fix wall/gran



Cannot specify both options at the same time.



E: Invalid wiggle direction for fix wall/gran



Self-explanatory.



E: Invalid shear direction for fix wall/gran



Self-explanatory.



E: Fix wall/gran is incompatible with Pair style



Must use a granular pair style to define the parameters needed for

this fix.



*/

CRA Git | Maintained and supported by SUSTech CRA and CCSE