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

src/GRANULAR/pair_gran_hertz_history.cpp

@@ -89,6 +89,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
   double *rmass = atom->rmass;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
 
   inum = list->inum;
   ilist = list->ilist;
@@ -254,7 +255,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
         torque[i][1] -= radi*tor2;
         torque[i][2] -= radi*tor3;
 
-        if (j < nlocal) {
+        if (newton_pair || j < nlocal) {
           f[j][0] -= fx;
           f[j][1] -= fy;
           f[j][2] -= fz;
@@ -268,6 +269,8 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
       }
     }
   }
+
+  if (vflag_fdotr) virial_fdotr_compute();
 }
 
 /* ----------------------------------------------------------------------
@@ -321,7 +324,7 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype,
 
   if (rsq >= radsum*radsum) {
     fforce = 0.0;
-    svector[0] = svector[1] = svector[2] = svector[3] = 0.0;
+    for (int m = 0; m < single_extra; m++) svector[m] = 0.0;
     return 0.0;
   }
 
@@ -441,12 +444,22 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype,
     } else fs1 = fs2 = fs3 = fs = 0.0;
   }
 
-  // set all forces and return no energy
+  // set force and return no energy
 
   fforce = ccel;
+
+  // set single_extra quantities
+
   svector[0] = fs1;
   svector[1] = fs2;
   svector[2] = fs3;
   svector[3] = fs;
+  svector[4] = vn1;
+  svector[5] = vn2;
+  svector[6] = vn3;
+  svector[7] = vt1;
+  svector[8] = vt2;
+  svector[9] = vt3;
+
   return 0.0;
 }

src/GRANULAR/pair_gran_hooke.cpp

@@ -239,7 +239,7 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq,
 
   if (rsq >= radsum*radsum) {
     fforce = 0.0;
-    svector[0] = svector[1] = svector[2] = svector[3] = 0.0;
+    for (int m = 0; m < single_extra; m++) svector[m] = 0.0;
     return 0.0;
   }
 
@@ -318,14 +318,25 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq,
   if (vrel != 0.0) ft = MIN(fn,fs) / vrel;
   else ft = 0.0;
 
-  // set all forces and return no energy
+  // set force and return no energy
 
   fforce = ccel;
+
+
+  // set single_extra quantities
+
   svector[0] = -ft*vtr1;
   svector[1] = -ft*vtr2;
   svector[2] = -ft*vtr3;
   svector[3] = sqrt(svector[0]*svector[0] +
                     svector[1]*svector[1] +
                     svector[2]*svector[2]);
+  svector[4] = vn1;
+  svector[5] = vn2;
+  svector[6] = vn3;
+  svector[7] = vt1;
+  svector[8] = vt2;
+  svector[9] = vt3;
+
   return 0.0;
 }

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -46,8 +46,8 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
   history = 1;
   fix_history = NULL;
 
-  single_extra = 4;
+  single_extra = 10;
-  svector = new double[4];
+  svector = new double[10];
 
   neighprev = 0;
 
@@ -131,6 +131,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
   double *rmass = atom->rmass;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
 
   inum = list->inum;
   ilist = list->ilist;
@@ -295,7 +296,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
         torque[i][1] -= radi*tor2;
         torque[i][2] -= radi*tor3;
 
-        if (j < nlocal) {
+        if (newton_pair || j < nlocal) {
           f[j][0] -= fx;
           f[j][1] -= fy;
           f[j][2] -= fz;
@@ -309,6 +310,8 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
       }
     }
   }
+
+  if (vflag_fdotr) virial_fdotr_compute();
 }
 
 /* ----------------------------------------------------------------------
@@ -413,13 +416,8 @@ void PairGranHookeHistory::init_style()
   dt = update->dt;
 
   // if shear history is stored:
-  // check if newton flag is valid
   // if first init, create Fix needed for storing shear history
 
-  if (history && force->newton_pair == 1)
-    error->all(FLERR,
-               "Pair granular with shear history requires newton pair off");
-
   if (history && fix_history == NULL) {
     char **fixarg = new char*[3];
     fixarg[0] = (char *) "SHEAR_HISTORY";
@@ -624,7 +622,7 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype,
 
   if (rsq >= radsum*radsum) {
     fforce = 0.0;
-    svector[0] = svector[1] = svector[2] = svector[3] = 0.0;
+    for (int m = 0; m < single_extra; m++) svector[m] = 0.0;
     return 0.0;
   }
 
@@ -741,13 +739,23 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype,
     } else fs1 = fs2 = fs3 = fs = 0.0;
   }
 
-  // set all forces and return no energy
+  // set force and return no energy
 
   fforce = ccel;
+
+  // set single_extra quantities
+
   svector[0] = fs1;
   svector[1] = fs2;
   svector[2] = fs3;
   svector[3] = fs;
+  svector[4] = vn1;
+  svector[5] = vn2;
+  svector[6] = vn3;
+  svector[7] = vt1;
+  svector[8] = vt2;
+  svector[9] = vt3;
+
   return 0.0;
 }
 

src/comm.h

@@ -80,6 +80,7 @@ class Comm : protected Pointers {
   virtual void reverse_comm_pair(class Pair *) = 0;
   virtual void forward_comm_fix(class Fix *, int size=0) = 0;
   virtual void reverse_comm_fix(class Fix *, int size=0) = 0;
+  virtual void reverse_comm_fix_variable(class Fix *) = 0;
   virtual void forward_comm_compute(class Compute *) = 0;
   virtual void reverse_comm_compute(class Compute *) = 0;
   virtual void forward_comm_dump(class Dump *) = 0;

src/comm_brick.cpp

@@ -1037,6 +1037,50 @@ void CommBrick::reverse_comm_fix(Fix *fix, int size)
   }
 }
 
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Fix with variable size data
+   query fix for pack size to insure buf_send is big enough
+   handshake sizes before each Irecv/Send to insure buf_recv is big enough
+------------------------------------------------------------------------- */
+
+void CommBrick::reverse_comm_fix_variable(Fix *fix)
+{
+  int iswap,nsend,nrecv;
+  double *buf;
+  MPI_Request request;
+
+  for (iswap = nswap-1; iswap >= 0; iswap--) {
+
+    // pack buffer
+
+    nsend = fix->pack_reverse_comm_size(recvnum[iswap],firstrecv[iswap]);
+    if (nsend > maxsend) grow_send(nsend,0);
+    nsend = fix->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send);
+
+    // exchange with another proc
+    // if self, set recv buffer to send buffer
+
+    if (sendproc[iswap] != me) {
+      MPI_Sendrecv(&nsend,1,MPI_INT,recvproc[iswap],0,
+                   &nrecv,1,MPI_INT,sendproc[iswap],0,world,
+                   MPI_STATUS_IGNORE);
+      if (sendnum[iswap]) {
+        if (nrecv > maxrecv) grow_recv(nrecv);
+        MPI_Irecv(buf_recv,maxrecv,MPI_DOUBLE,sendproc[iswap],0,
+                  world,&request);
+      }
+      if (recvnum[iswap])
+        MPI_Send(buf_send,nsend,MPI_DOUBLE,recvproc[iswap],0,world);
+      if (sendnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
+      buf = buf_recv;
+    } else buf = buf_send;
+
+    // unpack buffer
+
+    fix->unpack_reverse_comm(sendnum[iswap],sendlist[iswap],buf);
+  }
+}
+
 /* ----------------------------------------------------------------------
    forward communication invoked by a Compute
    nsize used only to set recv buffer limit
@@ -1055,7 +1099,7 @@ void CommBrick::forward_comm_compute(Compute *compute)
     // pack buffer
 
     n = compute->pack_forward_comm(sendnum[iswap],sendlist[iswap],
-                           buf_send,pbc_flag[iswap],pbc[iswap]);
+                                   buf_send,pbc_flag[iswap],pbc[iswap]);
 
     // exchange with another proc
     // if self, set recv buffer to send buffer
@@ -1276,7 +1320,8 @@ int CommBrick::exchange_variable(int n, double *inbuf, double *&outbuf)
     nrecv += nrecv1;
     if (procgrid[dim] > 2) {
       MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][1],0,
-                   &nrecv2,1,MPI_INT,procneigh[dim][0],0,world,MPI_STATUS_IGNORE);
+                   &nrecv2,1,MPI_INT,procneigh[dim][0],0,world,
+                   MPI_STATUS_IGNORE);
       nrecv += nrecv2;
     } else nrecv2 = 0;
 

src/comm_brick.h

@@ -37,6 +37,8 @@ class CommBrick : public Comm {
                                                    // forward comm from a Fix
   virtual void reverse_comm_fix(class Fix *, int size=0);
                                                    // reverse comm from a Fix
+  virtual void reverse_comm_fix_variable(class Fix *);
+                                     // variable size reverse comm from a Fix
   virtual void forward_comm_compute(class Compute *);  // forward from a Compute
   virtual void reverse_comm_compute(class Compute *);  // reverse from a Compute
   virtual void forward_comm_dump(class Dump *);    // forward comm from a Dump

src/comm_tiled.cpp

@@ -1180,6 +1180,18 @@ void CommTiled::reverse_comm_fix(Fix *fix, int size)
   }
 }
 
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Fix with variable size data
+   query fix for all pack sizes to insure buf_send is big enough
+   handshake sizes before irregular comm to insure buf_recv is big enough
+   NOTE: how to setup one big buf recv with correct offsets ??
+------------------------------------------------------------------------- */
+
+void CommTiled::reverse_comm_fix_variable(Fix *fix)
+{
+  error->all(FLERR,"Reverse comm fix variable not yet supported by CommTiled");
+}
+
 /* ----------------------------------------------------------------------
    forward communication invoked by a Compute
    nsize used only to set recv buffer limit

src/comm_tiled.h

@@ -37,6 +37,8 @@ class CommTiled : public Comm {
                                                    // forward comm from a Fix
   virtual void reverse_comm_fix(class Fix *, int size=0);
                                                    // reverse comm from a Fix
+  virtual void reverse_comm_fix_variable(class Fix *);
+                                     // variable size reverse comm from a Fix
   void forward_comm_compute(class Compute *);  // forward from a Compute
   void reverse_comm_compute(class Compute *);  // reverse from a Compute
   void forward_comm_dump(class Dump *);    // forward comm from a Dump

src/fix.h

@@ -170,6 +170,7 @@ class Fix : protected Pointers {
 
   virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;}
   virtual void unpack_forward_comm(int, int, double *) {}
+  virtual int pack_reverse_comm_size(int, int) {return 0;}
   virtual int pack_reverse_comm(int, int, double *) {return 0;}
   virtual void unpack_reverse_comm(int, int *, double *) {}
 

src/fix_shear_history.cpp

@@ -29,6 +29,8 @@
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
+enum{NPARTNER,PERPARTNER};
+
 /* ---------------------------------------------------------------------- */
 
 FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) :
@@ -37,6 +39,12 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) :
   restart_peratom = 1;
   create_attribute = 1;
 
+  newton_pair = force->newton_pair;
+
+  if (newton_pair) comm_reverse = 1;   // just for single npartner value
+                                       // variable-size history communicated via
+                                       // reverse_comm_fix_variable()
+
   // perform initial allocation of atom-based arrays
   // register with atom class
 
@@ -56,6 +64,8 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) :
   int nlocal = atom->nlocal;
   for (int i = 0; i < nlocal; i++) npartner[i] = 0;
   maxtouch = 0;
+
+  nlocal_neigh = nall_neigh = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -131,40 +141,51 @@ void FixShearHistory::allocate_pages()
 /* ----------------------------------------------------------------------
    copy shear partner info from neighbor lists to atom arrays
    should be called whenever neighbor list stores current history info
-     and need to have atoms store the info
+     and need to store the info with owned atoms
    e.g. so atoms can migrate to new procs or between runs
      when atoms may be added or deleted (neighbor list becomes out-of-date)
    the next granular neigh list build will put this info back into neigh list
-   called during run before atom exchanges
+   called during run before atom exchanges, including for restart files
    called at end of run via post_run()
    do not call during setup of run (setup_pre_exchange)
      b/c there is no guarantee of a current neigh list (even on continued run)
    if run command does a 2nd run with pre = no, then no neigh list
      will be built, but old neigh list will still have the info
+   newton ON and newton OFF versions
+     newton ON does reverse comm to acquire shear partner info from ghost atoms
+     newton OFF works with smaller vectors that don't include ghost info
 ------------------------------------------------------------------------- */
 
 void FixShearHistory::pre_exchange()
+{
+  if (newton_pair) pre_exchange_newton();
+  else pre_exchange_no_newton();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixShearHistory::pre_exchange_newton()
 {
   int i,j,ii,jj,m,n,inum,jnum;
   int *ilist,*jlist,*numneigh,**firstneigh;
   int *touch,**firsttouch;
   double *shear,*allshear,**firstshear;
 
-  // nlocal may include atoms added since last neigh build
+  // NOTE: all operations until very end are with 
+  //   nlocal_neigh  <= current nlocal and nall_neigh
+  // b/c previous neigh list was built with nlocal_neigh,nghost_neigh
+  // nlocal can be larger if other fixes added atoms at this pre_exchange()
 
-  int nlocal = atom->nlocal;
+  // zero npartner for owned+ghost atoms
-
-  // zero npartner for all current atoms
   // clear 2 page data structures
 
-  for (i = 0; i < nlocal; i++) npartner[i] = 0;
+  for (i = 0; i < nall_neigh; i++) npartner[i] = 0;
 
   ipage->reset();
   dpage->reset();
 
   // 1st loop over neighbor list
-  // calculate npartner for each owned atom
+  // calculate npartner for owned+ghost atoms
-  // nlocal_neigh = nlocal when neigh list was built, may be smaller than nlocal
 
   tagint *tag = atom->tag;
   NeighList *list = pair->list;
@@ -175,8 +196,133 @@ void FixShearHistory::pre_exchange()
   firsttouch = list->listgranhistory->firstneigh;
   firstshear = list->listgranhistory->firstdouble;
 
-  int nlocal_neigh = 0;
+  for (ii = 0; ii < inum; ii++) {
-  if (inum) nlocal_neigh = ilist[inum-1] + 1;
+    i = ilist[ii];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    touch = firsttouch[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      if (touch[jj]) {
+        npartner[i]++;
+        j = jlist[jj];
+        j &= NEIGHMASK;
+        npartner[j]++;
+      }
+    }
+  }
+
+  // perform reverse comm to augment owned npartner counts with ghost counts
+
+  commflag = NPARTNER;
+  comm->reverse_comm_fix(this,0);
+
+  // get page chunks to store atom IDs and shear history for owned+ghost atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    n = npartner[i];
+    partner[i] = ipage->get(n);
+    shearpartner[i] = dpage->get(n);
+    if (partner[i] == NULL || shearpartner[i] == NULL) {
+      error->one(FLERR,"Shear history overflow, boost neigh_modify one");
+    }
+  }
+
+  for (i = nlocal_neigh; i < nall_neigh; i++) {
+    n = npartner[i];
+    partner[i] = ipage->get(n);
+    shearpartner[i] = dpage->get(n);
+    if (partner[i] == NULL || shearpartner[i] == NULL) {
+      error->one(FLERR,"Shear history overflow, boost neigh_modify one");
+    }
+  }
+
+  // 2nd loop over neighbor list
+  // store atom IDs and shear history for owned+ghost atoms
+  // re-zero npartner to use as counter
+
+  for (i = 0; i < nall_neigh; i++) npartner[i] = 0;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    jlist = firstneigh[i];
+    allshear = firstshear[i];
+    jnum = numneigh[i];
+    touch = firsttouch[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      if (touch[jj]) {
+        shear = &allshear[3*jj];
+        j = jlist[jj];
+        j &= NEIGHMASK;
+        m = npartner[i]++;
+        partner[i][m] = tag[j];
+        shearpartner[i][m][0] = shear[0];
+        shearpartner[i][m][1] = shear[1];
+        shearpartner[i][m][2] = shear[2];
+        m = npartner[j]++;
+        partner[j][m] = tag[i];
+        shearpartner[j][m][0] = -shear[0];
+        shearpartner[j][m][1] = -shear[1];
+        shearpartner[j][m][2] = -shear[2];
+      }
+    }
+  }
+
+  // perform reverse comm to augment
+  // owned atom partner/shearpartner with ghost info
+  // use variable variant b/c size of packed data can be arbitrarily large
+  //   if many touching neighbors for large particle
+
+  commflag = PERPARTNER;
+  comm->reverse_comm_fix_variable(this);
+
+  // set maxtouch = max # of partners of any owned atom
+  // bump up comm->maxexchange_fix if necessary
+
+  maxtouch = 0;
+  for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]);
+  comm->maxexchange_fix = MAX(comm->maxexchange_fix,4*maxtouch+1);
+
+  // zero npartner values from previous nlocal_neigh to current nlocal
+
+  int nlocal = atom->nlocal;
+  for (i = nlocal_neigh; i < nlocal; i++) npartner[i] = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixShearHistory::pre_exchange_no_newton()
+{
+  int i,j,ii,jj,m,n,inum,jnum;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int *touch,**firsttouch;
+  double *shear,*allshear,**firstshear;
+
+  // NOTE: all operations until very end are with nlocal_neigh <= current nlocal
+  // b/c previous neigh list was built with nlocal_neigh
+  // nlocal can be larger if other fixes added atoms at this pre_exchange()
+
+  // zero npartner for owned atoms
+  // clear 2 page data structures
+
+  for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0;
+
+  ipage->reset();
+  dpage->reset();
+
+  // 1st loop over neighbor list
+  // calculate npartner for owned atoms
+
+  tagint *tag = atom->tag;
+  NeighList *list = pair->list;
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+  firsttouch = list->listgranhistory->firstneigh;
+  firstshear = list->listgranhistory->firstdouble;
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
@@ -194,7 +340,7 @@ void FixShearHistory::pre_exchange()
     }
   }
 
-  // get page chunks to store atom IDs and shear history for my atoms
+  // get page chunks to store atom IDs and shear history for owned atoms
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
@@ -206,10 +352,10 @@ void FixShearHistory::pre_exchange()
   }
 
   // 2nd loop over neighbor list
-  // store atom IDs and shear history for my atoms
+  // store atom IDs and shear history for owned atoms
-  // re-zero npartner to use as counter for all my atoms
+  // re-zero npartner to use as counter
 
-  for (i = 0; i < nlocal; i++) npartner[i] = 0;
+  for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0;
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
@@ -243,10 +389,15 @@ void FixShearHistory::pre_exchange()
 
   // set maxtouch = max # of partners of any owned atom
   // bump up comm->maxexchange_fix if necessary
-
+  
   maxtouch = 0;
-  for (i = 0; i < nlocal; i++) maxtouch = MAX(maxtouch,npartner[i]);
+  for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]);
   comm->maxexchange_fix = MAX(comm->maxexchange_fix,4*maxtouch+1);
+
+  // zero npartner values from previous nlocal_neigh to current nlocal
+
+  int nlocal = atom->nlocal;
+  for (i = nlocal_neigh; i < nlocal; i++) npartner[i] = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -324,6 +475,83 @@ void FixShearHistory::set_arrays(int i)
   npartner[i] = 0;
 }
 
+/* ----------------------------------------------------------------------
+   only called by Comm::reverse_comm_fix_variable for PERPARTNER mode
+------------------------------------------------------------------------- */
+
+int FixShearHistory::pack_reverse_comm_size(int n, int first)
+{
+  int i,j,k,last;
+
+  int m = 0;
+  last = first + n;
+
+  for (i = first; i < last; i++)
+    m += 1 + 4*npartner[i];
+
+  return m;
+}
+
+/* ----------------------------------------------------------------------
+   two modes: NPARTNER and PERPARTNER
+------------------------------------------------------------------------- */
+
+int FixShearHistory::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i,j,k,last;
+
+  int m = 0;
+  last = first + n;
+
+  if (commflag == NPARTNER) {
+    for (i = first; i < last; i++) {
+      buf[m++] = npartner[i];
+    }
+  } else if (commflag == PERPARTNER) {
+    for (i = first; i < last; i++) {
+      buf[m++] = npartner[i];
+      for (int k = 0; k < npartner[i]; k++) {
+        buf[m++] = partner[i][k];
+        buf[m++] = shearpartner[i][k][0];
+        buf[m++] = shearpartner[i][k][1];
+        buf[m++] = shearpartner[i][k][2];
+      }
+    }
+  }
+
+  return m;
+}
+
+/* ----------------------------------------------------------------------
+   two modes: NPARTNER and PERPARTNER
+------------------------------------------------------------------------- */
+
+void FixShearHistory::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i,j,k,kk,ncount;
+
+  int m = 0;
+
+  if (commflag == NPARTNER) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      npartner[j] += static_cast<int> (buf[m++]);
+    }
+  } else if (commflag == PERPARTNER) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      ncount = static_cast<int> (buf[m++]);
+      for (int k = 0; k < ncount; k++) {
+        kk = npartner[j]++;
+        partner[j][kk] = static_cast<tagint> (buf[m++]);
+        shearpartner[j][kk][0] = buf[m++];
+        shearpartner[j][kk][1] = buf[m++];
+        shearpartner[j][kk][2] = buf[m++];
+      }
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    pack values in local atom-based arrays for exchange with another proc
 ------------------------------------------------------------------------- */

src/fix_shear_history.h

@@ -44,6 +44,10 @@ class FixShearHistory : public Fix {
   void grow_arrays(int);
   void copy_arrays(int, int, int);
   void set_arrays(int);
+
+  int pack_reverse_comm_size(int, int);
+  int pack_reverse_comm(int, int, double *);
+  void unpack_reverse_comm(int, int *, double *);
   int pack_exchange(int, double *);
   int unpack_exchange(int, double *);
   int pack_restart(int, double *);
@@ -52,11 +56,17 @@ class FixShearHistory : public Fix {
   int maxsize_restart();
 
  protected:
+  int newton_pair;
+  int nlocal_neigh;             // nlocal at last time neigh list was built
+  int nall_neigh;               // ditto for nlocal+nghost
+
   int *npartner;                // # of touching partners of each atom
   tagint **partner;             // global atom IDs for the partners
   double (**shearpartner)[3];   // 3 shear values with the partner
   int maxtouch;                 // max # of touching partners for my atoms
 
+  int commflag;                 // mode of reverse comm to get ghost info
+
   class Pair *pair;
   int *computeflag;             // computeflag in PairGranHookeHistory
 
@@ -64,6 +74,8 @@ class FixShearHistory : public Fix {
   MyPage<tagint> *ipage;        // pages of partner atom IDs
   MyPage<double[3]> *dpage;     // pages of shear history with partners
 
+  void pre_exchange_newton();
+  void pre_exchange_no_newton();
   void allocate_pages();
 };
 

src/neigh_gran.cpp

@@ -65,6 +65,8 @@ void Neighbor::granular_nsq_no_newton(NeighList *list)
 
   FixShearHistory *fix_history = list->fix_history;
   if (fix_history) {
+    fix_history->nlocal_neigh = nlocal;
+    fix_history->nall_neigh = nall;
     npartner = fix_history->npartner;
     partner = fix_history->partner;
     shearpartner = fix_history->shearpartner;
@@ -160,7 +162,6 @@ void Neighbor::granular_nsq_no_newton(NeighList *list)
 /* ----------------------------------------------------------------------
    granular particles
    N^2 / 2 search for neighbor pairs with full Newton's 3rd law
-   no shear history is allowed for this option
    pair added to list if atoms i and j are both owned and i < j
    if j is ghost only me or other proc adds pair
    decision based on itag,jtag tests
@@ -168,10 +169,20 @@ void Neighbor::granular_nsq_no_newton(NeighList *list)
 
 void Neighbor::granular_nsq_newton(NeighList *list)
 {
-  int i,j,n,itag,jtag,bitmask;
+  int i,j,m,n,nn,itag,jtag,bitmask;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
-  int *neighptr;
+  int *neighptr,*touchptr;
+  double *shearptr;
+
+  NeighList *listgranhistory;
+  int *npartner;
+  tagint **partner;
+  double (**shearpartner)[3];
+  int **firsttouch;
+  double **firstshear;
+  MyPage<int> *ipage_touch;
+  MyPage<double> *dpage_shear;
 
   double **x = atom->x;
   double *radius = atom->radius;
@@ -191,12 +202,35 @@ void Neighbor::granular_nsq_newton(NeighList *list)
   int **firstneigh = list->firstneigh;
   MyPage<int> *ipage = list->ipage;
 
+  FixShearHistory *fix_history = list->fix_history;
+  if (fix_history) {
+    fix_history->nlocal_neigh = nlocal;
+    fix_history->nall_neigh = nall;
+    npartner = fix_history->npartner;
+    partner = fix_history->partner;
+    shearpartner = fix_history->shearpartner;
+    listgranhistory = list->listgranhistory;
+    firsttouch = listgranhistory->firstneigh;
+    firstshear = listgranhistory->firstdouble;
+    ipage_touch = listgranhistory->ipage;
+    dpage_shear = listgranhistory->dpage;
+  }
+
   int inum = 0;
   ipage->reset();
+  if (fix_history) {
+    ipage_touch->reset();
+    dpage_shear->reset();
+  }
 
   for (i = 0; i < nlocal; i++) {
     n = 0;
     neighptr = ipage->vget();
+    if (fix_history) {
+      nn = 0;
+      touchptr = ipage_touch->vget();
+      shearptr = dpage_shear->vget();
+    }
 
     itag = tag[i];
     xtmp = x[i][0];
@@ -233,7 +267,34 @@ void Neighbor::granular_nsq_newton(NeighList *list)
       radsum = radi + radius[j];
       cutsq = (radsum+skin) * (radsum+skin);
 
-      if (rsq <= cutsq) neighptr[n++] = j;
+      if (rsq <= cutsq) {
+        neighptr[n++] = j;
+
+        if (fix_history) {
+          if (rsq < radsum*radsum) {
+            for (m = 0; m < npartner[i]; m++)
+              if (partner[i][m] == tag[j]) break;
+            if (m < npartner[i]) {
+              touchptr[n] = 1;
+              shearptr[nn++] = shearpartner[i][m][0];
+              shearptr[nn++] = shearpartner[i][m][1];
+              shearptr[nn++] = shearpartner[i][m][2];
+            } else {
+              touchptr[n] = 0;
+              shearptr[nn++] = 0.0;
+              shearptr[nn++] = 0.0;
+              shearptr[nn++] = 0.0;
+            }
+          } else {
+            touchptr[n] = 0;
+            shearptr[nn++] = 0.0;
+            shearptr[nn++] = 0.0;
+            shearptr[nn++] = 0.0;
+          }
+        }
+
+        n++;
+      }
     }
 
     ilist[inum++] = i;
@@ -242,6 +303,13 @@ void Neighbor::granular_nsq_newton(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+
+    if (fix_history) {
+      firsttouch[i] = touchptr;
+      firstshear[i] = shearptr;
+      ipage_touch->vgot(n);
+      dpage_shear->vgot(nn);
+    }
   }
 
   list->inum = inum;
@@ -297,6 +365,8 @@ void Neighbor::granular_bin_no_newton(NeighList *list)
 
   FixShearHistory *fix_history = list->fix_history;
   if (fix_history) {
+    fix_history->nlocal_neigh = nlocal;
+    fix_history->nall_neigh = nlocal + atom->nghost;
     npartner = fix_history->npartner;
     partner = fix_history->partner;
     shearpartner = fix_history->shearpartner;
@@ -398,17 +468,26 @@ void Neighbor::granular_bin_no_newton(NeighList *list)
 /* ----------------------------------------------------------------------
    granular particles
    binned neighbor list construction with full Newton's 3rd law
-   no shear history is allowed for this option
    each owned atom i checks its own bin and other bins in Newton stencil
    every pair stored exactly once by some processor
 ------------------------------------------------------------------------- */
 
 void Neighbor::granular_bin_newton(NeighList *list)
 {
-  int i,j,k,n,ibin;
+  int i,j,k,m,n,nn,ibin;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
-  int *neighptr;
+  int *neighptr,*touchptr;
+  double *shearptr;
+
+  NeighList *listgranhistory;
+  int *npartner;
+  tagint **partner;
+  double (**shearpartner)[3];
+  int **firsttouch;
+  double **firstshear;
+  MyPage<int> *ipage_touch;
+  MyPage<double> *dpage_shear;
 
   // bin local & ghost atoms
 
@@ -418,6 +497,7 @@ void Neighbor::granular_bin_newton(NeighList *list)
 
   double **x = atom->x;
   double *radius = atom->radius;
+  tagint *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
   tagint *molecule = atom->molecule;
@@ -431,12 +511,35 @@ void Neighbor::granular_bin_newton(NeighList *list)
   int *stencil = list->stencil;
   MyPage<int> *ipage = list->ipage;
 
+  FixShearHistory *fix_history = list->fix_history;
+  if (fix_history) {
+    fix_history->nlocal_neigh = nlocal;
+    fix_history->nall_neigh = nlocal + atom->nghost;
+    npartner = fix_history->npartner;
+    partner = fix_history->partner;
+    shearpartner = fix_history->shearpartner;
+    listgranhistory = list->listgranhistory;
+    firsttouch = listgranhistory->firstneigh;
+    firstshear = listgranhistory->firstdouble;
+    ipage_touch = listgranhistory->ipage;
+    dpage_shear = listgranhistory->dpage;
+  }
+
   int inum = 0;
   ipage->reset();
+  if (fix_history) {
+    ipage_touch->reset();
+    dpage_shear->reset();
+  }
 
   for (i = 0; i < nlocal; i++) {
     n = 0;
     neighptr = ipage->vget();
+    if (fix_history) {
+      nn = 0;
+      touchptr = ipage_touch->vget();
+      shearptr = dpage_shear->vget();
+    }
 
     xtmp = x[i][0];
     ytmp = x[i][1];
@@ -465,7 +568,34 @@ void Neighbor::granular_bin_newton(NeighList *list)
       radsum = radi + radius[j];
       cutsq = (radsum+skin) * (radsum+skin);
 
-      if (rsq <= cutsq) neighptr[n++] = j;
+      if (rsq <= cutsq) {
+        neighptr[n] = j;
+
+        if (fix_history) {
+          if (rsq < radsum*radsum) {
+            for (m = 0; m < npartner[i]; m++)
+              if (partner[i][m] == tag[j]) break;
+            if (m < npartner[i]) {
+              touchptr[n] = 1;
+              shearptr[nn++] = shearpartner[i][m][0];
+              shearptr[nn++] = shearpartner[i][m][1];
+              shearptr[nn++] = shearpartner[i][m][2];
+            } else {
+              touchptr[n] = 0;
+              shearptr[nn++] = 0.0;
+              shearptr[nn++] = 0.0;
+              shearptr[nn++] = 0.0;
+            }
+          } else {
+            touchptr[n] = 0;
+            shearptr[nn++] = 0.0;
+            shearptr[nn++] = 0.0;
+            shearptr[nn++] = 0.0;
+          }
+        }
+
+        n++;
+      }
     }
 
     // loop over all atoms in other bins in stencil, store every pair
@@ -482,7 +612,34 @@ void Neighbor::granular_bin_newton(NeighList *list)
         radsum = radi + radius[j];
         cutsq = (radsum+skin) * (radsum+skin);
 
-        if (rsq <= cutsq) neighptr[n++] = j;
+        if (rsq <= cutsq) {
+          neighptr[n] = j;
+
+          if (fix_history) {
+            if (rsq < radsum*radsum) {
+              for (m = 0; m < npartner[i]; m++)
+                if (partner[i][m] == tag[j]) break;
+              if (m < npartner[i]) {
+                touchptr[n] = 1;
+                shearptr[nn++] = shearpartner[i][m][0];
+                shearptr[nn++] = shearpartner[i][m][1];
+                shearptr[nn++] = shearpartner[i][m][2];
+              } else {
+                touchptr[n] = 0;
+                shearptr[nn++] = 0.0;
+                shearptr[nn++] = 0.0;
+                shearptr[nn++] = 0.0;
+              }
+            } else {
+              touchptr[n] = 0;
+              shearptr[nn++] = 0.0;
+              shearptr[nn++] = 0.0;
+              shearptr[nn++] = 0.0;
+            }
+          }
+          
+          n++;
+        }
       }
     }
 
@@ -492,6 +649,13 @@ void Neighbor::granular_bin_newton(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+
+    if (fix_history) {
+      firsttouch[i] = touchptr;
+      firstshear[i] = shearptr;
+      ipage_touch->vgot(n);
+      dpage_shear->vgot(nn);
+    }
   }
 
   list->inum = inum;
@@ -500,17 +664,26 @@ void Neighbor::granular_bin_newton(NeighList *list)
 /* ----------------------------------------------------------------------
    granular particles
    binned neighbor list construction with Newton's 3rd law for triclinic
-   no shear history is allowed for this option
    each owned atom i checks its own bin and other bins in triclinic stencil
    every pair stored exactly once by some processor
 ------------------------------------------------------------------------- */
 
 void Neighbor::granular_bin_newton_tri(NeighList *list)
 {
-  int i,j,k,n,ibin;
+  int i,j,k,m,n,nn,ibin;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   double radi,radsum,cutsq;
-  int *neighptr;
+  int *neighptr,*touchptr;
+  double *shearptr;
+
+  NeighList *listgranhistory;
+  int *npartner;
+  tagint **partner;
+  double (**shearpartner)[3];
+  int **firsttouch;
+  double **firstshear;
+  MyPage<int> *ipage_touch;
+  MyPage<double> *dpage_shear;
 
   // bin local & ghost atoms
 
@@ -520,6 +693,7 @@ void Neighbor::granular_bin_newton_tri(NeighList *list)
 
   double **x = atom->x;
   double *radius = atom->radius;
+  tagint *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
   tagint *molecule = atom->molecule;
@@ -533,12 +707,35 @@ void Neighbor::granular_bin_newton_tri(NeighList *list)
   int *stencil = list->stencil;
   MyPage<int> *ipage = list->ipage;
 
+  FixShearHistory *fix_history = list->fix_history;
+  if (fix_history) {
+    fix_history->nlocal_neigh = nlocal;
+    fix_history->nall_neigh = nlocal + atom->nghost;
+    npartner = fix_history->npartner;
+    partner = fix_history->partner;
+    shearpartner = fix_history->shearpartner;
+    listgranhistory = list->listgranhistory;
+    firsttouch = listgranhistory->firstneigh;
+    firstshear = listgranhistory->firstdouble;
+    ipage_touch = listgranhistory->ipage;
+    dpage_shear = listgranhistory->dpage;
+  }
+
   int inum = 0;
   ipage->reset();
+  if (fix_history) {
+    ipage_touch->reset();
+    dpage_shear->reset();
+  }
 
   for (i = 0; i < nlocal; i++) {
     n = 0;
     neighptr = ipage->vget();
+    if (fix_history) {
+      nn = 0;
+      touchptr = ipage_touch->vget();
+      shearptr = dpage_shear->vget();
+    }
 
     xtmp = x[i][0];
     ytmp = x[i][1];
@@ -572,7 +769,34 @@ void Neighbor::granular_bin_newton_tri(NeighList *list)
         radsum = radi + radius[j];
         cutsq = (radsum+skin) * (radsum+skin);
 
-        if (rsq <= cutsq) neighptr[n++] = j;
+        if (rsq <= cutsq) {
+          neighptr[n++] = j;
+
+          if (fix_history) {
+            if (rsq < radsum*radsum) {
+              for (m = 0; m < npartner[i]; m++)
+                if (partner[i][m] == tag[j]) break;
+              if (m < npartner[i]) {
+                touchptr[n] = 1;
+                shearptr[nn++] = shearpartner[i][m][0];
+                shearptr[nn++] = shearpartner[i][m][1];
+                shearptr[nn++] = shearpartner[i][m][2];
+              } else {
+                touchptr[n] = 0;
+                shearptr[nn++] = 0.0;
+                shearptr[nn++] = 0.0;
+                shearptr[nn++] = 0.0;
+              }
+            } else {
+              touchptr[n] = 0;
+              shearptr[nn++] = 0.0;
+              shearptr[nn++] = 0.0;
+              shearptr[nn++] = 0.0;
+            }
+          }
+
+          n++;
+        }
       }
     }
 
@@ -582,6 +806,13 @@ void Neighbor::granular_bin_newton_tri(NeighList *list)
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+
+    if (fix_history) {
+      firsttouch[i] = touchptr;
+      firstshear[i] = shearptr;
+      ipage_touch->vgot(n);
+      dpage_shear->vgot(nn);
+    }
   }
 
   list->inum = inum;