Commit 201466f4 authored by sjplimp's avatar sjplimp
Browse files 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@930 f3b2605a-c512-4ea7-a41b-209d697bcdaa
parent 54a88a9b

src/GRANULAR/fix_gran_diag.cpp

+55 −29
Original line number Diff line number Diff line
@@ -24,6 +24,7 @@ 
#include "pair.h"

#include "atom.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "modify.h"

#include "comm.h"

#include "domain.h"
@@ -116,7 +117,7 @@ void FixGranDiag::init() 
  // insure use of granular pair_style

  // set local values from Pair values



  Pair *pair = force->pair_match("gran");

  pair = force->pair_match("gran");

  if (pair == NULL)

    error->all("Fix gran/diag is incompatible with Pair style");

  int dampflag;
@@ -151,6 +152,11 @@ void FixGranDiag::init() 


void FixGranDiag::setup()

{

  // get neighbor list from Pair class

  // can't do this until setup since lists are setup by neighbor::init()



  list = pair->list;



  if (first) end_of_step();

  first = 0;

}
@@ -462,7 +468,7 @@ void FixGranDiag::deallocate() 


void FixGranDiag::stress_no_history()

{

  int i,j,k,m,numneigh;

  int i,j,ii,jj,m,inum,jnum;

  double xtmp,ytmp,ztmp,delx,dely,delz;

  double radi,radj,radsum,rsq,r;

  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -470,7 +476,7 @@ void FixGranDiag::stress_no_history() 
  double vtr1,vtr2,vtr3,vrel;

  double xmeff,damp,ccel,ccelx,ccely,ccelz;

  double fn,fs,ft,fs1,fs2,fs3;

  int *neighs;

  int *ilist,*jlist,*numneigh,**firstneigh;



  double **x = atom->x;

  double **v = atom->v;
@@ -481,19 +487,25 @@ void FixGranDiag::stress_no_history() 
  int nlocal = atom->nlocal;

  int newton_pair = force->newton_pair;



  inum = list->inum;

  ilist = list->ilist;

  numneigh = list->numneigh;

  firstneigh = list->firstneigh;



  // loop over all neighbors of my atoms

  // store stress for both atoms i and j



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

  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];

    radi = radius[i];

    neighs = neighbor->firstneigh[i];

    numneigh = neighbor->numneigh[i];

    jlist = firstneigh[i];

    jnum = numneigh[i];



    for (k = 0; k < numneigh; k++) {

      j = neighs[k];

    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];



      // skip if neither atom is in fix group
@@ -613,7 +625,7 @@ void FixGranDiag::stress_no_history() 


void FixGranDiag::stress_history()

{

  int i,j,k,m,numneigh;

  int i,j,ii,jj,m,inum,jnum;

  double xtmp,ytmp,ztmp,delx,dely,delz;

  double radi,radj,radsum,rsq,r;

  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -622,8 +634,8 @@ void FixGranDiag::stress_history() 
  double xmeff,damp,ccel,ccelx,ccely,ccelz;

  double fn,fs,fs1,fs2,fs3;

  double shrmag,rsht;

  int *neighs;

  double *firstshear,*shear;

  int *ilist,*jlist,*numneigh,**firstneigh;

  double *shear,*allshear,**firstshear;



  double **x = atom->x;

  double **v = atom->v;
@@ -634,20 +646,27 @@ void FixGranDiag::stress_history() 
  int nlocal = atom->nlocal;

  int newton_pair = force->newton_pair;



  inum = list->inum;

  ilist = list->ilist;

  numneigh = list->numneigh;

  firstneigh = list->listgranhistory->firstneigh;

  firstshear = list->listgranhistory->firstdouble;



  // loop over all neighbors of my atoms

  // store stress on both atoms i and j



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

  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];

    radi = radius[i];

    neighs = neighbor->firstneigh[i];

    firstshear = neighbor->firstshear[i];

    numneigh = neighbor->numneigh[i];

    allshear = firstshear[i];

    jlist = firstneigh[i];

    jnum = numneigh[i];



    for (k = 0; k < numneigh; k++) {

      j = neighs[k];

    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];



      // skip if neither atom is in fix group
@@ -718,7 +737,7 @@ void FixGranDiag::stress_history() 
	// shrmag = magnitude of shear

	// do not update shear history since not timestepping



	shear = &firstshear[3*k];

	shear = &allshear[3*jj];

	shrx = shear[0] + vtr1;

	shry = shear[1] + vtr2;

	shrz = shear[2] + vtr3;
@@ -796,7 +815,7 @@ void FixGranDiag::stress_history() 


void FixGranDiag::stress_hertzian()

{

  int i,j,k,m,numneigh;

  int i,j,ii,jj,m,inum,jnum;

  double xtmp,ytmp,ztmp,delx,dely,delz;

  double radi,radj,radsum,rsq,r;

  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -805,8 +824,8 @@ void FixGranDiag::stress_hertzian() 
  double xmeff,damp,ccel,ccelx,ccely,ccelz;

  double fn,fs,fs1,fs2,fs3;

  double shrmag,rsht,rhertz;

  int *neighs;

  double *firstshear,*shear;

  int *ilist,*jlist,*numneigh,**firstneigh;

  double *shear,*allshear,**firstshear;



  double **x = atom->x;

  double **v = atom->v;
@@ -817,20 +836,27 @@ void FixGranDiag::stress_hertzian() 
  int nlocal = atom->nlocal;

  int newton_pair = force->newton_pair;



  inum = list->inum;

  ilist = list->ilist;

  numneigh = list->numneigh;

  firstneigh = list->firstneigh;

  firstshear = listgran->firstdouble;



  // loop over all neighbors of my atoms

  // store stress on both atoms i and j



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

  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];

    radi = radius[i];

    neighs = neighbor->firstneigh[i];

    firstshear = neighbor->firstshear[i];

    numneigh = neighbor->numneigh[i];

    allshear = firstshear[i];

    jlist = firstneigh[i];

    jnum = numneigh[i];



    for (k = 0; k < numneigh; k++) {

      j = neighs[k];

    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];



      // skip if neither atom is in fix group
@@ -903,7 +929,7 @@ void FixGranDiag::stress_hertzian() 
	// shrmag = magnitude of shear

	// do not update shear history since not timestepping



	shear = &firstshear[3*k];

	shear = &allshear[3*jj];

	shrx = shear[0] + vtr1;

	shry = shear[1] + vtr2;

	shrz = shear[2] + vtr3;

src/GRANULAR/fix_gran_diag.h

+3 −0
Original line number Diff line number Diff line
@@ -45,6 +45,9 @@ class FixGranDiag : public Fix { 
  double *velfxx,*velfyy,*velfzz,*velfxy,*velfxz,*velfyz;

  double *sigx2,*sigy2,*sigz2,*sigxy2,*sigxz2,*sigyz2;



  class Pair *pair;

  class NeighList *list,*listgran;



  void allocate();

  void deallocate();

  void stress_no_history();

src/GRANULAR/pair_gran_hertzian.cpp

+24 −15
Original line number Diff line number Diff line
@@ -21,7 +21,7 @@ 
#include "pair_gran_hertzian.h"

#include "atom.h"

#include "force.h"

#include "neighbor.h"

#include "neigh_list.h"



using namespace LAMMPS_NS;
@@ -39,7 +39,7 @@ PairGranHertzian::PairGranHertzian(LAMMPS *lmp) : PairGranHistory(lmp) 


void PairGranHertzian::compute(int eflag, int vflag)

{

  int i,j,k,numneigh;

  int i,j,ii,jj,inum,jnum;

  double xtmp,ytmp,ztmp,delx,dely,delz;

  double radi,radj,radsum,rsq,r,rinv;

  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -48,8 +48,9 @@ void PairGranHertzian::compute(int eflag, int vflag) 
  double xmeff,damp,ccel,ccelx,ccely,ccelz,tor1,tor2,tor3;

  double fn,fs,fs1,fs2,fs3;

  double shrmag,rsht,rhertz;

  int *neighs,*touch;

  double *firstshear,*shear;

  int *ilist,*jlist,*numneigh,**firstneigh;

  int *touch,**firsttouch;

  double *shear,*allshear,**firstshear;



  double **f = atom->f;

  double **x = atom->x;
@@ -62,20 +63,28 @@ void PairGranHertzian::compute(int eflag, int vflag) 
  int nlocal = atom->nlocal;

  int newton_pair = force->newton_pair;



  inum = list->inum;

  ilist = list->ilist;

  numneigh = list->numneigh;

  firstneigh = list->firstneigh;

  firsttouch = list->listgranhistory->firstneigh;

  firstshear = list->listgranhistory->firstdouble;



  // loop over neighbors of my atoms



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

  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];

    radi = radius[i];

    neighs = neighbor->firstneigh[i];

    touch = neighbor->firsttouch[i];

    firstshear = neighbor->firstshear[i];

    numneigh = neighbor->numneigh[i];

    touch = firsttouch[i];

    allshear = firstshear[i];

    jlist = firstneigh[i];

    jnum = numneigh[i];



    for (k = 0; k < numneigh; k++) {

      j = neighs[k];

    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];



      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];
@@ -88,8 +97,8 @@ void PairGranHertzian::compute(int eflag, int vflag) 


	// unset touching neighbors



        touch[k] = 0;

	shear = &firstshear[3*k];

        touch[jj] = 0;

	shear = &allshear[3*jj];

        shear[0] = 0.0;

        shear[1] = 0.0;

        shear[2] = 0.0;
@@ -152,8 +161,8 @@ void PairGranHertzian::compute(int eflag, int vflag) 
	// shear history effects

	// shrmag = magnitude of shear



	touch[k] = 1;

	shear = &firstshear[3*k];

	touch[jj] = 1;

	shear = &allshear[3*jj];

        shear[0] += vtr1;

        shear[1] += vtr2;

        shear[2] += vtr3;

src/GRANULAR/pair_gran_history.cpp

+70 −35
Original line number Diff line number Diff line
@@ -27,9 +27,12 @@ 
#include "modify.h"

#include "fix.h"

#include "fix_pour.h"

#include "fix_shear_history.h"

#include "comm.h"

#include "memory.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "neigh_request.h"

#include "memory.h"

#include "error.h"



using namespace LAMMPS_NS;
@@ -64,7 +67,7 @@ PairGranHistory::~PairGranHistory() 


void PairGranHistory::compute(int eflag, int vflag)

{

  int i,j,k,numneigh;

  int i,j,ii,jj,inum,jnum;

  double xtmp,ytmp,ztmp,delx,dely,delz;

  double radi,radj,radsum,rsq,r,rinv;

  double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@@ -73,8 +76,9 @@ void PairGranHistory::compute(int eflag, int vflag) 
  double xmeff,damp,ccel,ccelx,ccely,ccelz,tor1,tor2,tor3;

  double fn,fs,fs1,fs2,fs3;

  double shrmag,rsht;

  int *neighs,*touch;

  double *firstshear,*shear;

  int *ilist,*jlist,*numneigh,**firstneigh;

  int *touch,**firsttouch;

  double *shear,*allshear,**firstshear;



  double **f = atom->f;

  double **x = atom->x;
@@ -87,20 +91,28 @@ void PairGranHistory::compute(int eflag, int vflag) 
  int nlocal = atom->nlocal;

  int newton_pair = force->newton_pair;



  inum = list->inum;

  ilist = list->ilist;

  numneigh = list->numneigh;

  firstneigh = list->firstneigh;

  firsttouch = listgranhistory->firstneigh;

  firstshear = listgranhistory->firstdouble;

 

  // loop over neighbors of my atoms



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

  for (ii = 0; ii < inum; ii++) {

    i = ilist[ii];

    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];

    radi = radius[i];

    neighs = neighbor->firstneigh[i];

    touch = neighbor->firsttouch[i];

    firstshear = neighbor->firstshear[i];

    numneigh = neighbor->numneigh[i];

    touch = firsttouch[i];

    allshear = firstshear[i];

    jlist = firstneigh[i];

    jnum = numneigh[i];



    for (k = 0; k < numneigh; k++) {

      j = neighs[k];

    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];



      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];
@@ -113,8 +125,8 @@ void PairGranHistory::compute(int eflag, int vflag) 


	// unset touching neighbors



        touch[k] = 0;

	shear = &firstshear[3*k];

        touch[jj] = 0;

	shear = &allshear[3*jj];

        shear[0] = 0.0;

        shear[1] = 0.0;

        shear[2] = 0.0;
@@ -175,8 +187,8 @@ void PairGranHistory::compute(int eflag, int vflag) 
	// shear history effects

	// shrmag = magnitude of shear



	touch[k] = 1;

	shear = &firstshear[3*k];

	touch[jj] = 1;

	shear = &allshear[3*jj];

        shear[0] += vtr1;

        shear[1] += vtr2;

        shear[2] += vtr3;
@@ -274,8 +286,6 @@ void PairGranHistory::allocate() 


void PairGranHistory::settings(int narg, char **arg)

{

  if (domain->box_exist == 0)

    error->all("Pair_style granular command before simulation box is defined");

  if (narg != 4) error->all("Illegal pair_style command");



  xkk = atof(arg[0]);
@@ -302,21 +312,6 @@ void PairGranHistory::coeff(int narg, char **arg) 
  error->all("Granular pair styles do not use pair_coeff settings");

}



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

   init for one type pair i,j and corresponding j,i

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



double PairGranHistory::init_one(int i, int j)

{

  if (!allocated) allocate();



  // return dummy value used in neighbor setup,

  // but not in actual neighbor calculation

  // since particles have variable radius



  return 1.0;

}



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

   init specific to this pair style

------------------------------------------------------------------------- */
@@ -328,6 +323,19 @@ void PairGranHistory::init_style() 
  if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag)

    error->all("Pair granular requires atom attributes radius, omega, torque");



  // need a half neigh list and optionally a granular history neigh list



  int irequest = neighbor->request(this);

  neighbor->requests[irequest]->half = 0;

  neighbor->requests[irequest]->gran = 1;

  if (history) {

    irequest = neighbor->request(this);

    neighbor->requests[irequest]->id = 1;

    neighbor->requests[irequest]->half = 0;

    neighbor->requests[irequest]->granhistory = 1;

    neighbor->requests[irequest]->dnum = 3;

  }



  xkkt = xkk * 2.0/7.0;

  dt = update->dt;

  double gammas = 0.5*gamman;
@@ -344,12 +352,13 @@ void PairGranHistory::init_style() 


  if (history && fix_history == NULL) {

    char **fixarg = new char*[3];

    fixarg[0] = "SHEAR_HISTORY";

    fixarg[1] = "all";

    fixarg[2] = "SHEAR_HISTORY";

    fixarg[0] = (char *) "SHEAR_HISTORY";

    fixarg[1] = (char *) "all";

    fixarg[2] = (char *) "SHEAR_HISTORY";

    modify->add_fix(3,fixarg);

    delete [] fixarg;

    fix_history = (FixShearHistory *) modify->fix[modify->nfix-1];

    fix_history->pair = this;

  }



  // check for freeze Fix and set freeze_group_bit
@@ -393,6 +402,32 @@ void PairGranHistory::init_style() 
  cutforce = maxrad_dynamic + MAX(maxrad_dynamic,maxrad_frozen);

}



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

   neighbor callback to inform pair style of neighbor list to use

   optional granular history list

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



void PairGranHistory::init_list(int id, NeighList *ptr)

{

  if (id == 0) list = ptr;

  else if (id == 1) listgranhistory = ptr;

}



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

   init for one type pair i,j and corresponding j,i

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



double PairGranHistory::init_one(int i, int j)

{

  if (!allocated) allocate();



  // return dummy value used in neighbor setup,

  // but not in actual neighbor calculation

  // since particles have variable radius



  return 1.0;

}



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

  proc 0 writes to restart file

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

src/GRANULAR/pair_gran_history.h

+2 −1
Original line number Diff line number Diff line
@@ -25,8 +25,9 @@ class PairGranHistory : public Pair { 
  virtual void compute(int, int);

  void settings(int, char **);

  void coeff(int, char **);

  double init_one(int, int);

  void init_style();

  void init_list(int, class NeighList *);

  double init_one(int, int);

  void write_restart(FILE *);

  void read_restart(FILE *);

  void write_restart_settings(FILE *);

CRA Git | Maintained and supported by SUSTech CRA and CCSE