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

src/BODY/body_nparticle.cpp

@@ -98,25 +98,16 @@ int BodyNparticle::unpack_border_body(AtomVecBody::Bonus *bonus, double *buf)
 ------------------------------------------------------------------------- */
 
 void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
-                              char **ifile, char **dfile)
+                              int *ifile, double *dfile)
 {
   AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
 
-  // error in data file if any values are NULL
-
-  for (int i = 0; i < ninteger; i++)
-    if (ifile[i] == NULL)
-      error->one(FLERR,"Invalid format in Bodies section of data file");
-  for (int i = 0; i < ndouble; i++)
-    if (dfile[i] == NULL)
-      error->one(FLERR,"Invalid format in Bodies section of data file");
-
   // set ninteger, ndouble in bonus and allocate 2 vectors of ints, doubles
 
   if (ninteger != 1)
     error->one(FLERR,"Incorrect # of integer values in "
                "Bodies section of data file");
-  int nsub = atoi(ifile[0]);
+  int nsub = ifile[0];
   if (nsub < 1)
     error->one(FLERR,"Incorrect integer value in "
                "Bodies section of data file");
@@ -133,12 +124,12 @@ void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
   // diagonalize inertia tensor
 
   double tensor[3][3];
-  tensor[0][0] = atof(dfile[0]);
-  tensor[1][1] = atof(dfile[1]);
-  tensor[2][2] = atof(dfile[2]);
-  tensor[0][1] = tensor[1][0] = atof(dfile[3]);
-  tensor[0][2] = tensor[2][0] = atof(dfile[4]);
-  tensor[1][2] = tensor[2][1] = atof(dfile[5]);
+  tensor[0][0] = dfile[0];
+  tensor[1][1] = dfile[1];
+  tensor[2][2] = dfile[2];
+  tensor[0][1] = tensor[1][0] = dfile[3];
+  tensor[0][2] = tensor[2][0] = dfile[4];
+  tensor[1][2] = tensor[2][1] = dfile[5];
 
   double *inertia = bonus->inertia;
   double evectors[3][3];
@@ -188,9 +179,9 @@ void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
   int j = 6;
   int k = 0;
   for (int i = 0; i < nsub; i++) {
-    delta[0] = atof(dfile[j]);
-    delta[1] = atof(dfile[j+1]);
-    delta[2] = atof(dfile[j+2]);
+    delta[0] = dfile[j];
+    delta[1] = dfile[j+1];
+    delta[2] = dfile[j+2];
     MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
                                 delta,&bonus->dvalue[k]);
     j += 3;
@@ -198,6 +189,43 @@ void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
   }
 }
 
+/* ----------------------------------------------------------------------
+   return radius of body particle defined by ifile/dfile params
+   params are ordered as in data file
+   called by Molecule class which needs single body size
+------------------------------------------------------------------------- */
+
+double BodyNparticle::radius_body(int ninteger, int ndouble,
+				  int *ifile, double *dfile)
+{
+  int nsub = ifile[0];
+  if (nsub < 1)
+    error->one(FLERR,"Incorrect integer value in "
+               "Bodies section of data file");
+  if (ndouble != 6 + 3*nsub)
+    error->one(FLERR,"Incorrect # of floating-point values in "
+               "Bodies section of data file");
+
+  // sub-particle coords are relative to body center at (0,0,0)
+  // offset = 6 for sub-particle coords
+
+  double onerad;
+  double maxrad = 0.0;
+  double delta[3];
+  
+  int offset = 6;          
+  for (int i = 0; i < nsub; i++) {
+    delta[0] = dfile[offset];
+    delta[1] = dfile[offset+1];
+    delta[2] = dfile[offset+2];
+    offset += 3;
+    onerad = MathExtra::len3(delta);
+    maxrad = MAX(maxrad,onerad);
+  }
+
+  return maxrad;
+}
+
 /* ---------------------------------------------------------------------- */
 
 int BodyNparticle::noutcol()

src/BODY/body_nparticle.h

@@ -34,7 +34,8 @@ class BodyNparticle : public Body {
 
   int pack_border_body(struct AtomVecBody::Bonus *, double *);
   int unpack_border_body(struct AtomVecBody::Bonus *, double *);
-  void data_body(int, int, int, char **, char **);
+  void data_body(int, int, int, int *, double *);
+  double radius_body(int, int, int *, double *);
 
   int noutrow(int);
   int noutcol();

src/GRANULAR/fix_pour.cpp

@@ -378,6 +378,13 @@ void FixPour::pre_exchange()
 
   if (next_reneighbor != update->ntimestep) return;
 
+  // clear ghost count and any ghost bonus data internal to AtomVec
+  // same logic as beginning of Comm::exchange()
+  // do it now b/c inserting atoms will overwrite ghost atoms
+
+  atom->nghost = 0;
+  atom->avec->clear_bonus();
+
   // find current max atom and molecule IDs if necessary
 
   if (!idnext) find_maxid();
@@ -638,7 +645,11 @@ void FixPour::pre_exchange()
           radtmp = newcoord[3];
           atom->radius[n] = radtmp;
           atom->rmass[n] = 4.0*MY_PI/3.0 * radtmp*radtmp*radtmp * denstmp;
-        } else atom->add_molecule_atom(onemols[imol],m,n,maxtag_all);
+        } else {
+	  onemols[imol]->quat_external = quat;
+	  atom->add_molecule_atom(onemols[imol],m,n,maxtag_all);
+	}
+	
         modify->create_attribute(n);
       }
     }
@@ -667,7 +678,8 @@ void FixPour::pre_exchange()
   // reset global natoms,nbonds,etc
   // increment maxtag_all and maxmol_all if necessary
   // if global map exists, reset it now instead of waiting for comm
-  // since adding atoms messes up ghosts
+  //   since other pre-exchange fixes may use it
+  //   invoke map_init() b/c atom count has grown
 
   if (ninserted_atoms) {
     atom->natoms += ninserted_atoms;
@@ -682,7 +694,6 @@ void FixPour::pre_exchange()
     if (maxtag_all >= MAXTAGINT)
       error->all(FLERR,"New atom IDs exceed maximum allowed ID");
     if (atom->map_style) {
-      atom->nghost = 0;
       atom->map_init();
       atom->map_set();
     }
@@ -1008,23 +1019,24 @@ void *FixPour::extract(const char *str, int &itype)
 
     } else {
 
-      // find a molecule in template with matching type
+      // loop over onemols molecules
+      // skip a molecule with no atoms as large as itype
 
+      oneradius = 0.0;
       for (int i = 0; i < nmol; i++) {
-        if (itype-ntype > onemols[i]->ntypes) continue;
+        if (itype > ntype+onemols[i]->ntypes) continue;
         double *radius = onemols[i]->radius;
         int *type = onemols[i]->type;
         int natoms = onemols[i]->natoms;
 
-        // check radii of matching types in Molecule
+        // check radii of atoms in Molecule with matching types
         // default to 0.5, if radii not defined in Molecule
         //   same as atom->avec->create_atom(), invoked in pre_exchange()
 
-        oneradius = 0.0;
         for (int i = 0; i < natoms; i++)
-          if (type[i] == itype-ntype) {
+          if (type[i]+ntype == itype) {
             if (radius) oneradius = MAX(oneradius,radius[i]);
-            else oneradius = 0.5;
+            else oneradius = MAX(oneradius,0.5);
           }
       }
     }