lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2007-03-08 00:49:32 +00:00
parent ab24117893
commit e1b5b447f1
5 changed files with 73 additions and 65 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
src/GRANULAR/atom_vec_granular.cpp
View File

@ -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;

8
src/GRANULAR/atom_vec_granular.h
View File

@ -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();

68
src/GRANULAR/fix_pour.cpp
View File

@ -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;
}
}

2
src/GRANULAR/fix_pour.h
View File

@ -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 *);
};
}

10
src/GRANULAR/fix_wall_gran.cpp
View File

@ -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");