From e1b5b447f1dc31c3aabfdd836b8e508e06e1e71b Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 8 Mar 2007 00:49:32 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@365 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GRANULAR/atom_vec_granular.cpp | 50 ++++++++++------------ src/GRANULAR/atom_vec_granular.h | 8 ++-- src/GRANULAR/fix_pour.cpp | 68 ++++++++++++++++-------------- src/GRANULAR/fix_pour.h | 2 +- src/GRANULAR/fix_wall_gran.cpp | 10 +++++ 5 files changed, 73 insertions(+), 65 deletions(-) diff --git a/src/GRANULAR/atom_vec_granular.cpp b/src/GRANULAR/atom_vec_granular.cpp index 8744d79c8f..2e5a5a9c81 100644 --- a/src/GRANULAR/atom_vec_granular.cpp +++ b/src/GRANULAR/atom_vec_granular.cpp @@ -15,7 +15,6 @@ #include "stdlib.h" #include "atom_vec_granular.h" #include "atom.h" -#include "domain.h" #include "force.h" #include "modify.h" #include "fix.h" @@ -177,12 +176,13 @@ void AtomVecGranular::copy(int i, int j) /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_comm(int n, int *list, double *buf, int *pbc_flags) +int AtomVecGranular::pack_comm(int n, int *list, double *buf, + int pbc_flag, double *pbc_dist) { int i,j,m; m = 0; - if (pbc_flags[0] == 0) { + if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; @@ -196,14 +196,11 @@ int AtomVecGranular::pack_comm(int n, int *list, double *buf, int *pbc_flags) buf[m++] = phiv[j][2]; } } else { - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = x[j][0] + pbc_flags[1]*xprd; - buf[m++] = x[j][1] + pbc_flags[2]*yprd; - buf[m++] = x[j][2] + pbc_flags[3]*zprd; + buf[m++] = x[j][0] + pbc_dist[0]; + buf[m++] = x[j][1] + pbc_dist[1]; + buf[m++] = x[j][2] + pbc_dist[2]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; @@ -321,12 +318,13 @@ int AtomVecGranular::unpack_reverse_one(int i, double *buf) /* ---------------------------------------------------------------------- */ -int AtomVecGranular::pack_border(int n, int *list, double *buf, int *pbc_flags) +int AtomVecGranular::pack_border(int n, int *list, double *buf, + int pbc_flag, double *pbc_dist) { int i,j,m; m = 0; - if (pbc_flags[0] == 0) { + if (pbc_flag == 0) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = x[j][0]; @@ -345,14 +343,11 @@ int AtomVecGranular::pack_border(int n, int *list, double *buf, int *pbc_flags) buf[m++] = phiv[j][2]; } } else { - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = x[j][0] + pbc_flags[1]*xprd; - buf[m++] = x[j][1] + pbc_flags[2]*yprd; - buf[m++] = x[j][2] + pbc_flags[3]*zprd; + buf[m++] = x[j][0] + pbc_dist[0]; + buf[m++] = x[j][1] + pbc_dist[1]; + buf[m++] = x[j][2] + pbc_dist[2]; buf[m++] = tag[j]; buf[m++] = type[j]; buf[m++] = mask[j]; @@ -621,21 +616,20 @@ int AtomVecGranular::unpack_restart(double *buf) } /* ---------------------------------------------------------------------- - create one atom of type itype at x0,y0,z0 + create one atom of itype at coord set other values to defaults ------------------------------------------------------------------------- */ -void AtomVecGranular::create_atom(int itype, double x0, double y0, double z0, - int ihybrid) +void AtomVecGranular::create_atom(int itype, double *coord, int ihybrid) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); tag[nlocal] = 0; type[nlocal] = itype; - x[nlocal][0] = x0; - x[nlocal][1] = y0; - x[nlocal][2] = z0; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; mask[nlocal] = 1; image[nlocal] = (512 << 20) | (512 << 10) | 512; v[nlocal][0] = 0.0; @@ -664,8 +658,8 @@ void AtomVecGranular::create_atom(int itype, double x0, double y0, double z0, initialize other atom quantities ------------------------------------------------------------------------- */ -void AtomVecGranular::data_atom(double xtmp, double ytmp, double ztmp, - int imagetmp, char **values, int ihybrid) +void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values, + int ihybrid) { int nlocal = atom->nlocal; if (nlocal == nmax) grow(0); @@ -686,9 +680,9 @@ void AtomVecGranular::data_atom(double xtmp, double ytmp, double ztmp, else rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal]; - x[nlocal][0] = xtmp; - x[nlocal][1] = ytmp; - x[nlocal][2] = ztmp; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; image[nlocal] = imagetmp; diff --git a/src/GRANULAR/atom_vec_granular.h b/src/GRANULAR/atom_vec_granular.h index 30987ef710..15d2887a3a 100644 --- a/src/GRANULAR/atom_vec_granular.h +++ b/src/GRANULAR/atom_vec_granular.h @@ -27,7 +27,7 @@ class AtomVecGranular : public AtomVec { void zero_owned(int); void zero_ghost(int, int); void copy(int, int); - int pack_comm(int, int *, double *, int *); + int pack_comm(int, int *, double *, int, double *); int pack_comm_one(int, double *); void unpack_comm(int, int, double *); int unpack_comm_one(int, double *); @@ -35,7 +35,7 @@ class AtomVecGranular : public AtomVec { int pack_reverse_one(int, double *); void unpack_reverse(int, int *, double *); int unpack_reverse_one(int, double *); - int pack_border(int, int *, double *, int *); + int pack_border(int, int *, double *, int, double *); int pack_border_one(int, double *); void unpack_border(int, int, double *); int unpack_border_one(int, double *); @@ -45,8 +45,8 @@ class AtomVecGranular : public AtomVec { int size_restart_one(int); int pack_restart(int, double *); int unpack_restart(double *); - void create_atom(int, double, double, double, int); - void data_atom(double, double, double, int, char **, int); + void create_atom(int, double *, int); + void data_atom(double *, int, char **, int); void data_vel(int, char *, int); int memory_usage(); diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 843b85be94..1381e15a5a 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -45,6 +45,8 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : if (atom->check_style("granular") == 0) error->all("Must use fix pour with atom style granular"); + if (domain->triclinic) error->all("Cannot use fix pour with triclinic box"); + // required args ninsert = atoi(arg[3]); @@ -365,7 +367,7 @@ void FixPour::pre_exchange() // h = height, biased to give uniform distribution in time of insertion int success; - double xtmp,ytmp,ztmp,radtmp,delx,dely,delz,rsq,radsum,rn,h; + double coord[3],radtmp,delx,dely,delz,rsq,radsum,rn,h; int attempt = 0; int max = nnew * maxattempt; @@ -378,11 +380,11 @@ void FixPour::pre_exchange() success = 0; while (attempt < max) { attempt++; - xyz_random(h,xtmp,ytmp,ztmp); + xyz_random(h,coord); for (i = 0; i < nnear; i++) { - delx = xtmp - xnear[i][0]; - dely = ytmp - xnear[i][1]; - delz = ztmp - xnear[i][2]; + delx = coord[0] - xnear[i][0]; + dely = coord[1] - xnear[i][1]; + delz = coord[2] - xnear[i][2]; rsq = delx*delx + dely*dely + delz*delz; radsum = radtmp + xnear[i][3]; if (rsq <= radsum*radsum) break; @@ -393,9 +395,9 @@ void FixPour::pre_exchange() } } if (success) { - xnear[nnear][0] = xtmp; - xnear[nnear][1] = ytmp; - xnear[nnear][2] = ztmp; + xnear[nnear][0] = coord[0]; + xnear[nnear][1] = coord[1]; + xnear[nnear][2] = coord[2]; xnear[nnear][3] = radtmp; nnear++; } else break; @@ -421,37 +423,39 @@ void FixPour::pre_exchange() int m,flag; double denstmp,vxtmp,vytmp,vztmp; double g = 1.0; + double *sublo = domain->sublo; + double *subhi = domain->subhi; for (i = nprevious; i < nnear; i++) { - xtmp = xnear[i][0]; - ytmp = xnear[i][1]; - ztmp = xnear[i][2]; + coord[0] = xnear[i][0]; + coord[1] = xnear[i][1]; + coord[2] = xnear[i][2]; radtmp = xnear[i][3]; denstmp = density_lo + random->uniform() * (density_hi-density_lo); if (force->dimension == 3) { vxtmp = vxlo + random->uniform() * (vxhi-vxlo); vytmp = vylo + random->uniform() * (vyhi-vylo); - vztmp = vz - sqrt(2.0*g*(hi_current-ztmp)); + vztmp = vz - sqrt(2.0*g*(hi_current-coord[2])); } else { vxtmp = vxlo + random->uniform() * (vxhi-vxlo); - vytmp = vy - sqrt(2.0*g*(hi_current-ytmp)); + vytmp = vy - sqrt(2.0*g*(hi_current-coord[1])); vztmp = 0.0; } flag = 0; - if (xtmp >= domain->subxlo && xtmp < domain->subxhi && - ytmp >= domain->subylo && ytmp < domain->subyhi && - ztmp >= domain->subzlo && ztmp < domain->subzhi) flag = 1; - else if (force->dimension == 3 && ztmp >= domain->boxzhi && + if (coord[0] >= sublo[0] && coord[0] < subhi[0] && + coord[1] >= sublo[1] && coord[1] < subhi[1] && + coord[2] >= sublo[2] && coord[2] < subhi[2]) flag = 1; + else if (force->dimension == 3 && coord[2] >= domain->boxzhi && comm->myloc[2] == comm->procgrid[2]-1 && - xtmp >= domain->subxlo && xtmp < domain->subxhi && - ytmp >= domain->subylo && ytmp < domain->subyhi) flag = 1; - else if (force->dimension == 2 && ytmp >= domain->boxyhi && + coord[0] >= sublo[0] && coord[0] < subhi[0] && + coord[1] >= sublo[1] && coord[1] < subhi[1]) flag = 1; + else if (force->dimension == 2 && coord[1] >= domain->boxyhi && comm->myloc[1] == comm->procgrid[1]-1 && - xtmp >= domain->subxlo && xtmp < domain->subxhi) flag = 1; + coord[0] >= sublo[0] && coord[0] < subhi[0]) flag = 1; if (flag) { - avec->create_atom(ntype,xtmp,ytmp,ztmp,0); + avec->create_atom(ntype,coord,0); m = atom->nlocal - 1; atom->type[m] = ntype; atom->radius[m] = radtmp; @@ -526,13 +530,13 @@ int FixPour::overlap(int i) /* ---------------------------------------------------------------------- */ -void FixPour::xyz_random(double h, double &x, double &y, double &z) +void FixPour::xyz_random(double h, double *coord) { if (force->dimension == 3) { if (region_style == 1) { - x = xlo + random->uniform() * (xhi-xlo); - y = ylo + random->uniform() * (yhi-ylo); - z = h; + coord[0] = xlo + random->uniform() * (xhi-xlo); + coord[1] = ylo + random->uniform() * (yhi-ylo); + coord[2] = h; } else { double r1,r2; while (1) { @@ -540,13 +544,13 @@ void FixPour::xyz_random(double h, double &x, double &y, double &z) r2 = random->uniform() - 0.5; if (r1*r1 + r2*r2 < 0.25) break; } - x = xc + 2.0*r1*rc; - y = yc + 2.0*r2*rc; - z = h; + coord[0] = xc + 2.0*r1*rc; + coord[1] = yc + 2.0*r2*rc; + coord[2] = h; } } else { - x = xlo + random->uniform() * (xhi-xlo); - y = h; - z = 0.0; + coord[0] = xlo + random->uniform() * (xhi-xlo); + coord[1] = h; + coord[2] = 0.0; } } diff --git a/src/GRANULAR/fix_pour.h b/src/GRANULAR/fix_pour.h index c2a011218d..0155ad6e14 100644 --- a/src/GRANULAR/fix_pour.h +++ b/src/GRANULAR/fix_pour.h @@ -51,7 +51,7 @@ class FixPour : public Fix { class RanPark *random; int overlap(int); - void xyz_random(double, double &, double &, double &); + void xyz_random(double, double *); }; } diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 5850d4ade0..7d5116d1e1 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -20,6 +20,7 @@ #include "string.h" #include "fix_wall_gran.h" #include "atom.h" +#include "domain.h" #include "update.h" #include "force.h" #include "pair.h" @@ -108,6 +109,15 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : } else error->all("Illegal fix wall/gran command"); } + if (wallstyle == XPLANE && domain->xperiodic) + error->all("Cannot use wall in periodic dimension"); + if (wallstyle == YPLANE && domain->yperiodic) + error->all("Cannot use wall in periodic dimension"); + if (wallstyle == ZPLANE && domain->zperiodic) + error->all("Cannot use wall in periodic dimension"); + if (wallstyle == ZCYLINDER && (domain->xperiodic || domain->yperiodic)) + error->all("Cannot use wall in periodic dimension"); + if (wallstyle == ZCYLINDER && wiggle) if (axis != 2) error->all("Can only wiggle zcylinder wall in z dim");