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

src/GRANULAR/pair_gran_hertz_history.cpp

@@ -27,6 +27,7 @@
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "comm.h"
+#include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
@@ -60,12 +61,23 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
   int shearupdate = 1;
   if (update->setupflag) shearupdate = 0;
 
-  // update rigid body ptrs and values for ghost atoms if using FixRigid masses
+  // update rigid body info for owned & ghost atoms if using FixRigid masses
+  // body[i] = which body atom I is in, -1 if none
+  // mass_body = mass of each rigid body
 
   if (fix_rigid && neighbor->ago == 0) {
     int tmp;
-    body = (int *) fix_rigid->extract("body",tmp);
+    int *body = (int *) fix_rigid->extract("body",tmp);
-    mass_rigid = (double *) fix_rigid->extract("masstotal",tmp);
+    double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
+    if (atom->nmax > nmax) {
+      memory->destroy(mass_rigid);
+      nmax = atom->nmax;
+      memory->create(mass_rigid,nmax,"pair:mass_rigid");
+    }
+    int nlocal = atom->nlocal;
+    for (i = 0; i < nlocal; i++)
+      if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
+      else mass_rigid[i] = 0.0;
     comm->forward_comm_pair(this);
   }
 
@@ -164,8 +176,8 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
           mj = mass[type[j]];
         }
         if (fix_rigid) {
-          if (body[i] >= 0) mi = mass_rigid[body[i]];
+          if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
-          if (body[j] >= 0) mj = mass_rigid[body[j]];
+          if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
         }
 
         meff = mi*mj / (mi+mj);
@@ -374,13 +386,9 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype,
     mj = mass[type[j]];
   }
   if (fix_rigid) {
-    // NOTE: need to make sure ghost atoms have updated body?
+    // NOTE: insure mass_rigid is current for owned+ghost atoms?
-    // depends on where single() is called from
+    if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
-    int tmp;
+    if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
-    body = (int *) fix_rigid->extract("body",tmp);
-    mass_rigid = (double *) fix_rigid->extract("masstotal",tmp);
-    if (body[i] >= 0) mi = mass_rigid[body[i]];
-    if (body[j] >= 0) mj = mass_rigid[body[j]];
   }
 
   meff = mi*mj / (mi+mj);

src/GRANULAR/pair_gran_hooke.cpp

@@ -25,6 +25,7 @@
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "comm.h"
+#include "memory.h"
 
 using namespace LAMMPS_NS;
 
@@ -53,12 +54,23 @@ void PairGranHooke::compute(int eflag, int vflag)
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = 0;
 
-  // update rigid body ptrs and values for ghost atoms if using FixRigid masses
+  // update rigid body info for owned & ghost atoms if using FixRigid masses
+  // body[i] = which body atom I is in, -1 if none
+  // mass_body = mass of each rigid body
 
   if (fix_rigid && neighbor->ago == 0) {
     int tmp;
-    body = (int *) fix_rigid->extract("body",tmp);
+    int *body = (int *) fix_rigid->extract("body",tmp);
-    mass_rigid = (double *) fix_rigid->extract("masstotal",tmp);
+    double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
+    if (atom->nmax > nmax) {
+      memory->destroy(mass_rigid);
+      nmax = atom->nmax;
+      memory->create(mass_rigid,nmax,"pair:mass_rigid");
+    }
+    int nlocal = atom->nlocal;
+    for (i = 0; i < nlocal; i++)
+      if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
+      else mass_rigid[i] = 0.0;
     comm->forward_comm_pair(this);
   }
 
@@ -144,8 +156,8 @@ void PairGranHooke::compute(int eflag, int vflag)
           mj = mass[type[j]];
         }
         if (fix_rigid) {
-          if (body[i] >= 0) mi = mass_rigid[body[i]];
+          if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
-          if (body[j] >= 0) mj = mass_rigid[body[j]];
+          if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
         }
 
         meff = mi*mj / (mi+mj);
@@ -290,14 +302,11 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq,
     mi = mass[type[i]];
     mj = mass[type[j]];
   }
+
   if (fix_rigid) {
-    // NOTE: need to make sure ghost atoms have updated body?
+    // NOTE: insure mass_rigid is current for owned+ghost atoms?
-    // depends on where single() is called from
+    if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
-    int tmp;
+    if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
-    body = (int *) fix_rigid->extract("body",tmp);
-    mass_rigid = (double *) fix_rigid->extract("masstotal",tmp);
-    if (body[i] >= 0) mi = mass_rigid[body[i]];
-    if (body[j] >= 0) mj = mass_rigid[body[j]];
   }
 
   meff = mi*mj / (mi+mj);

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -53,6 +53,13 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
 
   computeflag = 0;
   neighprev = 0;
+
+  nmax = 0;
+  mass_rigid = NULL;
+
+  // set comm size needed by this Pair if used with fix rigid
+
+  comm_forward = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -72,6 +79,8 @@ PairGranHookeHistory::~PairGranHookeHistory()
     delete [] maxrad_dynamic;
     delete [] maxrad_frozen;
   }
+
+  memory->destroy(mass_rigid);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -98,12 +107,23 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
   int shearupdate = 1;
   if (update->setupflag) shearupdate = 0;
 
-  // update rigid body ptrs and values for ghost atoms if using FixRigid masses
+  // update rigid body info for owned & ghost atoms if using FixRigid masses
+  // body[i] = which body atom I is in, -1 if none
+  // mass_body = mass of each rigid body
 
   if (fix_rigid && neighbor->ago == 0) {
     int tmp;
-    body = (int *) fix_rigid->extract("body",tmp);
+    int *body = (int *) fix_rigid->extract("body",tmp);
-    mass_rigid = (double *) fix_rigid->extract("masstotal",tmp);
+    double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
+    if (atom->nmax > nmax) {
+      memory->destroy(mass_rigid);
+      nmax = atom->nmax;
+      memory->create(mass_rigid,nmax,"pair:mass_rigid");
+    }
+    int nlocal = atom->nlocal;
+    for (i = 0; i < nlocal; i++)
+      if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
+      else mass_rigid[i] = 0.0;
     comm->forward_comm_pair(this);
   }
 
@@ -202,8 +222,8 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
           mj = mass[type[j]];
         }
         if (fix_rigid) {
-          if (body[i] >= 0) mi = mass_rigid[body[i]];
+          if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
-          if (body[j] >= 0) mj = mass_rigid[body[j]];
+          if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
         }
 
         meff = mi*mj / (mi+mj);
@@ -653,13 +673,9 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype,
     mj = mass[type[j]];
   }
   if (fix_rigid) {
-    // NOTE: need to make sure ghost atoms have updated body?
+    // NOTE: insure mass_rigid is current for owned+ghost atoms?
-    // depends on where single() is called from
+    if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
-    int tmp;
+    if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
-    body = (int *) fix_rigid->extract("body",tmp);
-    mass_rigid = (double *) fix_rigid->extract("masstotal",tmp);
-    if (body[i] >= 0) mi = mass_rigid[body[i]];
-    if (body[j] >= 0) mj = mass_rigid[body[j]];
   }
 
   meff = mi*mj / (mi+mj);
@@ -745,7 +761,7 @@ int PairGranHookeHistory::pack_comm(int n, int *list,
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
-    buf[m++] = body[j];
+    buf[m++] = mass_rigid[j];
   }
   return 1;
 }
@@ -759,7 +775,7 @@ void PairGranHookeHistory::unpack_comm(int n, int first, double *buf)
   m = 0;
   last = first + n;
   for (i = first; i < last; i++)
-    body[i] = static_cast<int> (buf[m++]);
+    mass_rigid[i] = buf[m++];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -770,3 +786,13 @@ void *PairGranHookeHistory::extract(const char *str, int &dim)
   if (strcmp(str,"computeflag") == 0) return (void *) &computeflag;
   return NULL;
 }
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double PairGranHookeHistory::memory_usage()
+{
+  double bytes = nmax * sizeof(double);
+  return bytes;
+}

src/GRANULAR/pair_gran_hooke_history.h

@@ -45,6 +45,7 @@ class PairGranHookeHistory : public Pair {
   int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   void *extract(const char *, int &);
+  double memory_usage();
 
  protected:
   double kn,kt,gamman,gammat,xmu;
@@ -59,9 +60,12 @@ class PairGranHookeHistory : public Pair {
   double *maxrad_dynamic,*maxrad_frozen;
 
   class FixShearHistory *fix_history;
-  class Fix *fix_rigid;
+
-  int *body;
+  // storage of rigid body masses for use in granular interactions
-  double *mass_rigid;
+
+  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 allocate();
 };

src/RIGID/fix_rigid_small.cpp

@@ -151,6 +151,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
       bodyown[i] = nlocal_body++;
     } else bodyown[i] = -1;
 
+
   // bodysize = sizeof(Body) in doubles
 
   bodysize = sizeof(Body)/sizeof(double);
@@ -217,6 +218,11 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   random = NULL;
   if (langflag) random = new RanMars(lmp,seed + comm->me);
 
+  // mass vector for granular pair styles
+
+  mass_body = NULL;
+  nmax_mass = 0;
+
   // firstflag = 1 triggers one-time initialization of rigid body attributes
 
   firstflag = 1;
@@ -244,6 +250,7 @@ FixRigidSmall::~FixRigidSmall()
 
   delete random;
   memory->destroy(langextra);
+  memory->destroy(mass_body);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1289,10 +1296,7 @@ void FixRigidSmall::create_bodies()
   n = 0;
   for (i = 0; i < nlocal; i++) {
     if (!(mask[i] & groupbit)) continue;
-    if (hash->find(molecule[i]) == hash->end()) {
+    if (hash->find(molecule[i]) == hash->end()) (*hash)[molecule[i]] = n++;
-      hash->insert(std::pair<int,int> (molecule[i],n));
-      n++;
-    }
   }
 
   // bbox = bounding box of each rigid body my atoms are part of
@@ -1552,7 +1556,7 @@ void FixRigidSmall::setup_bodies()
       atom->tri_flag || atom->mu_flag) {
     int flag = 0;
     for (i = 0; i < nlocal; i++) {
-      if (atom2body[i] < 0) continue;
+      if (bodytag[i] == 0) continue;
       if (radius && radius[i] > 0.0) flag = 1;
       if (ellipsoid && ellipsoid[i] >= 0) flag = 1;
       if (line && line[i] >= 0) flag = 1;
@@ -1577,7 +1581,7 @@ void FixRigidSmall::setup_bodies()
 
     for (i = 0; i < nlocal; i++) {
       eflags[i] = 0;
-      if (atom2body[i] < 0) continue;
+      if (bodytag[i] == 0) continue;
 
       // set to POINT or SPHERE or ELLIPSOID or LINE
 
@@ -2518,6 +2522,37 @@ void FixRigidSmall::reset_dt()
   dtq = 0.5 * update->dt;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void *FixRigidSmall::extract(const char *str, int &dim)
+{
+  if (strcmp(str,"body") == 0) {
+    dim = 1;
+    return atom2body;
+  }
+
+  // return vector of rigid body masses, for owned+ghost bodies
+  // used by granular pair styles, indexed by atom2body
+
+  if (strcmp(str,"masstotal") == 0) {
+    dim = 1;
+
+    if (nmax_mass < nmax_body) {
+      memory->destroy(mass_body);
+      nmax_mass = nmax_body;
+      memory->create(mass_body,nmax_mass,"rigid:mass_body");
+    }
+
+    int n = nlocal_body + nghost_body;
+    for (int i = 0; i < n; i++)
+      mass_body[i] = body[i].mass;
+
+    return mass_body;
+  }
+
+  return NULL;
+}
+
 /* ----------------------------------------------------------------------
    return temperature of collection of rigid bodies
    non-active DOF are removed by fflag/tflag and in tfactor

src/RIGID/fix_rigid_small.h

@@ -59,6 +59,7 @@ class FixRigidSmall : public Fix {
   int dof(int);
   void deform(int);
   void reset_dt();
+  void *extract(const char*,int &);
   double compute_scalar();
   double memory_usage();
 
@@ -128,6 +129,11 @@ class FixRigidSmall : public Fix {
   int **counts;            // counts of atom types in bodies
   double **itensor;        // 6 space-frame components of inertia tensor
 
+  // mass per body, accessed by granular pair styles
+
+  double *mass_body;
+  int nmax_mass;
+
   // Langevin thermostatting
 
   int langflag;                     // 0/1 = no/yes Langevin thermostat