From 201466f4f1ddebdccff6e315a7b1d85efe5b8cd7 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 3 Oct 2007 16:15:11 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@930 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GRANULAR/fix_gran_diag.cpp | 84 ++++++++++++++------- src/GRANULAR/fix_gran_diag.h | 3 + src/GRANULAR/pair_gran_hertzian.cpp | 39 ++++++---- src/GRANULAR/pair_gran_history.cpp | 105 +++++++++++++++++--------- src/GRANULAR/pair_gran_history.h | 3 +- src/GRANULAR/pair_gran_no_history.cpp | 22 ++++-- 6 files changed, 168 insertions(+), 88 deletions(-) diff --git a/src/GRANULAR/fix_gran_diag.cpp b/src/GRANULAR/fix_gran_diag.cpp index c541afe35a..d99ef04642 100644 --- a/src/GRANULAR/fix_gran_diag.cpp +++ b/src/GRANULAR/fix_gran_diag.cpp @@ -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,13 +117,13 @@ 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; double gamman; pair->extract_gran(&xkk,&gamman,&xmu,&dampflag); - + // same initialization as in pair_gran_history::init_style() xkkt = xkk * 2.0/7.0; @@ -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; diff --git a/src/GRANULAR/fix_gran_diag.h b/src/GRANULAR/fix_gran_diag.h index 8186eb4d3e..86fc8ec6b9 100644 --- a/src/GRANULAR/fix_gran_diag.h +++ b/src/GRANULAR/fix_gran_diag.h @@ -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(); diff --git a/src/GRANULAR/pair_gran_hertzian.cpp b/src/GRANULAR/pair_gran_hertzian.cpp index 2dad7b077e..84f87414b4 100644 --- a/src/GRANULAR/pair_gran_hertzian.cpp +++ b/src/GRANULAR/pair_gran_hertzian.cpp @@ -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; diff --git a/src/GRANULAR/pair_gran_history.cpp b/src/GRANULAR/pair_gran_history.cpp index b1f94cdf49..a9859650a5 100644 --- a/src/GRANULAR/pair_gran_history.cpp +++ b/src/GRANULAR/pair_gran_history.cpp @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/GRANULAR/pair_gran_history.h b/src/GRANULAR/pair_gran_history.h index 0ad8bb3b57..ddab1e7c9e 100644 --- a/src/GRANULAR/pair_gran_history.h +++ b/src/GRANULAR/pair_gran_history.h @@ -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 *); diff --git a/src/GRANULAR/pair_gran_no_history.cpp b/src/GRANULAR/pair_gran_no_history.cpp index 01e8a1f3a3..e1dd1831b5 100644 --- a/src/GRANULAR/pair_gran_no_history.cpp +++ b/src/GRANULAR/pair_gran_no_history.cpp @@ -21,7 +21,7 @@ #include "pair_gran_no_history.h" #include "atom.h" #include "force.h" -#include "neighbor.h" +#include "neigh_list.h" using namespace LAMMPS_NS; @@ -39,7 +39,7 @@ PairGranNoHistory::PairGranNoHistory(LAMMPS *lmp) : PairGranHistory(lmp) void PairGranNoHistory::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; @@ -47,7 +47,7 @@ void PairGranNoHistory::compute(int eflag, int vflag) double vtr1,vtr2,vtr3,vrel; double xmeff,damp,ccel,ccelx,ccely,ccelz,tor1,tor2,tor3; double fn,fs,ft,fs1,fs2,fs3; - int *neighs; + int *ilist,*jlist,*numneigh,**firstneigh; double **f = atom->f; double **x = atom->x; @@ -60,18 +60,24 @@ void PairGranNoHistory::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; + // 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]; - 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]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1];