jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2007-04-13 15:23:13 +00:00
parent c501cef961
commit d15f2ee3b6
9 changed files with 200 additions and 200 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

213
src/GRANULAR/atom_vec_granular.cpp
View File

@ -73,12 +73,12 @@ void AtomVecGranular::grow(int n)
rmass = atom->rmass = (double *) rmass = atom->rmass = (double *)
memory->srealloc(atom->rmass,nmax*sizeof(double),"atom:rmass"); memory->srealloc(atom->rmass,nmax*sizeof(double),"atom:rmass");
phix = atom->phix = xphi = atom->xphi =
memory->grow_2d_double_array(atom->phix,nmax,3,"atom:phix"); memory->grow_2d_double_array(atom->xphi,nmax,3,"atom:xphi");
phiv = atom->phiv = omega = atom->omega =
memory->grow_2d_double_array(atom->phiv,nmax,3,"atom:phiv"); memory->grow_2d_double_array(atom->omega,nmax,3,"atom:omega");
phia = atom->phia = torque = atom->torque =
memory->grow_2d_double_array(atom->phia,nmax,3,"atom:phia"); memory->grow_2d_double_array(atom->torque,nmax,3,"atom:torque");
if (atom->nextra_grow) if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++) for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@ -100,9 +100,9 @@ void AtomVecGranular::reset_ptrs()
radius = atom->radius; radius = atom->radius;
density = atom->density; density = atom->density;
rmass = atom->rmass; rmass = atom->rmass;
phix = atom->phix; xphi = atom->xphi;
phiv = atom->phiv; omega = atom->omega;
phia = atom->phia; torque = atom->torque;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -115,12 +115,12 @@ void AtomVecGranular::zero_owned(int i)
radius[i] = 0.0; radius[i] = 0.0;
density[i] = 0.0; density[i] = 0.0;
rmass[i] = 0.0; rmass[i] = 0.0;
phix[i][0] = 0.0; xphi[i][0] = 0.0;
phix[i][1] = 0.0; xphi[i][1] = 0.0;
phix[i][2] = 0.0; xphi[i][2] = 0.0;
phiv[i][0] = 0.0; omega[i][0] = 0.0;
phiv[i][1] = 0.0; omega[i][1] = 0.0;
phiv[i][2] = 0.0; omega[i][2] = 0.0;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -139,9 +139,9 @@ void AtomVecGranular::zero_ghost(int n, int first)
v[i][2] = 0.0; v[i][2] = 0.0;
radius[i] = 0.0; radius[i] = 0.0;
rmass[i] = 0.0; rmass[i] = 0.0;
phiv[i][0] = 0.0; omega[i][0] = 0.0;
phiv[i][1] = 0.0; omega[i][1] = 0.0;
phiv[i][2] = 0.0; omega[i][2] = 0.0;
} }
} }
@ -163,12 +163,12 @@ void AtomVecGranular::copy(int i, int j)
radius[j] = radius[i]; radius[j] = radius[i];
density[j] = density[i]; density[j] = density[i];
rmass[j] = rmass[i]; rmass[j] = rmass[i];
phix[j][0] = phix[i][0]; xphi[j][0] = xphi[i][0];
phix[j][1] = phix[i][1]; xphi[j][1] = xphi[i][1];
phix[j][2] = phix[i][2]; xphi[j][2] = xphi[i][2];
phiv[j][0] = phiv[i][0]; omega[j][0] = omega[i][0];
phiv[j][1] = phiv[i][1]; omega[j][1] = omega[i][1];
phiv[j][2] = phiv[i][2]; omega[j][2] = omega[i][2];
if (atom->nextra_grow) if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++) for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@ -193,9 +193,9 @@ int AtomVecGranular::pack_comm(int n, int *list, double *buf,
buf[m++] = v[j][0]; buf[m++] = v[j][0];
buf[m++] = v[j][1]; buf[m++] = v[j][1];
buf[m++] = v[j][2]; buf[m++] = v[j][2];
buf[m++] = phiv[j][0]; buf[m++] = omega[j][0];
buf[m++] = phiv[j][1]; buf[m++] = omega[j][1];
buf[m++] = phiv[j][2]; buf[m++] = omega[j][2];
} }
} else { } else {
if (domain->triclinic == 0) { if (domain->triclinic == 0) {
@ -215,9 +215,9 @@ int AtomVecGranular::pack_comm(int n, int *list, double *buf,
buf[m++] = v[j][0]; buf[m++] = v[j][0];
buf[m++] = v[j][1]; buf[m++] = v[j][1];
buf[m++] = v[j][2]; buf[m++] = v[j][2];
buf[m++] = phiv[j][0]; buf[m++] = omega[j][0];
buf[m++] = phiv[j][1]; buf[m++] = omega[j][1];
buf[m++] = phiv[j][2]; buf[m++] = omega[j][2];
} }
} }
return m; return m;
@ -230,9 +230,9 @@ int AtomVecGranular::pack_comm_one(int i, double *buf)
buf[0] = v[i][0]; buf[0] = v[i][0];
buf[1] = v[i][1]; buf[1] = v[i][1];
buf[2] = v[i][2]; buf[2] = v[i][2];
buf[3] = phiv[i][0]; buf[3] = omega[i][0];
buf[4] = phiv[i][1]; buf[4] = omega[i][1];
buf[5] = phiv[i][2]; buf[5] = omega[i][2];
return 6; return 6;
} }
@ -251,9 +251,9 @@ void AtomVecGranular::unpack_comm(int n, int first, double *buf)
v[i][0] = buf[m++]; v[i][0] = buf[m++];
v[i][1] = buf[m++]; v[i][1] = buf[m++];
v[i][2] = buf[m++]; v[i][2] = buf[m++];
phiv[i][0] = buf[m++]; omega[i][0] = buf[m++];
phiv[i][1] = buf[m++]; omega[i][1] = buf[m++];
phiv[i][2] = buf[m++]; omega[i][2] = buf[m++];
} }
} }
@ -264,9 +264,9 @@ int AtomVecGranular::unpack_comm_one(int i, double *buf)
v[i][0] = buf[0]; v[i][0] = buf[0];
v[i][1] = buf[1]; v[i][1] = buf[1];
v[i][2] = buf[2]; v[i][2] = buf[2];
phiv[i][0] = buf[3]; omega[i][0] = buf[3];
phiv[i][1] = buf[4]; omega[i][1] = buf[4];
phiv[i][2] = buf[5]; omega[i][2] = buf[5];
return 6; return 6;
} }
@ -282,9 +282,9 @@ int AtomVecGranular::pack_reverse(int n, int first, double *buf)
buf[m++] = f[i][0]; buf[m++] = f[i][0];
buf[m++] = f[i][1]; buf[m++] = f[i][1];
buf[m++] = f[i][2]; buf[m++] = f[i][2];
buf[m++] = phia[i][0]; buf[m++] = torque[i][0];
buf[m++] = phia[i][1]; buf[m++] = torque[i][1];
buf[m++] = phia[i][2]; buf[m++] = torque[i][2];
} }
return m; return m;
} }
@ -293,9 +293,9 @@ int AtomVecGranular::pack_reverse(int n, int first, double *buf)
int AtomVecGranular::pack_reverse_one(int i, double *buf) int AtomVecGranular::pack_reverse_one(int i, double *buf)
{ {
buf[0] = phia[i][0]; buf[0] = torque[i][0];
buf[1] = phia[i][1]; buf[1] = torque[i][1];
buf[2] = phia[i][2]; buf[2] = torque[i][2];
return 3; return 3;
} }
@ -311,9 +311,9 @@ void AtomVecGranular::unpack_reverse(int n, int *list, double *buf)
f[j][0] += buf[m++]; f[j][0] += buf[m++];
f[j][1] += buf[m++]; f[j][1] += buf[m++];
f[j][2] += buf[m++]; f[j][2] += buf[m++];
phia[j][0] += buf[m++]; torque[j][0] += buf[m++];
phia[j][1] += buf[m++]; torque[j][1] += buf[m++];
phia[j][2] += buf[m++]; torque[j][2] += buf[m++];
} }
} }
@ -321,9 +321,9 @@ void AtomVecGranular::unpack_reverse(int n, int *list, double *buf)
int AtomVecGranular::unpack_reverse_one(int i, double *buf) int AtomVecGranular::unpack_reverse_one(int i, double *buf)
{ {
phia[i][0] = buf[0]; torque[i][0] = buf[0];
phia[i][1] = buf[1]; torque[i][1] = buf[1];
phia[i][2] = buf[2]; torque[i][2] = buf[2];
return 3; return 3;
} }
@ -350,9 +350,9 @@ int AtomVecGranular::pack_border(int n, int *list, double *buf,
buf[m++] = v[j][2]; buf[m++] = v[j][2];
buf[m++] = radius[j]; buf[m++] = radius[j];
buf[m++] = rmass[j]; buf[m++] = rmass[j];
buf[m++] = phiv[j][0]; buf[m++] = omega[j][0];
buf[m++] = phiv[j][1]; buf[m++] = omega[j][1];
buf[m++] = phiv[j][2]; buf[m++] = omega[j][2];
} }
} else { } else {
if (domain->triclinic == 0) { if (domain->triclinic == 0) {
@ -377,9 +377,9 @@ int AtomVecGranular::pack_border(int n, int *list, double *buf,
buf[m++] = v[j][2]; buf[m++] = v[j][2];
buf[m++] = radius[j]; buf[m++] = radius[j];
buf[m++] = rmass[j]; buf[m++] = rmass[j];
buf[m++] = phiv[j][0]; buf[m++] = omega[j][0];
buf[m++] = phiv[j][1]; buf[m++] = omega[j][1];
buf[m++] = phiv[j][2]; buf[m++] = omega[j][2];
} }
} }
return m; return m;
@ -394,9 +394,9 @@ int AtomVecGranular::pack_border_one(int i, double *buf)
buf[2] = v[i][2]; buf[2] = v[i][2];
buf[3] = radius[i]; buf[3] = radius[i];
buf[4] = rmass[i]; buf[4] = rmass[i];
buf[5] = phiv[i][0]; buf[5] = omega[i][0];
buf[6] = phiv[i][1]; buf[6] = omega[i][1];
buf[7] = phiv[i][2]; buf[7] = omega[i][2];
return 8; return 8;
} }
@ -421,9 +421,9 @@ void AtomVecGranular::unpack_border(int n, int first, double *buf)
v[i][2] = buf[m++]; v[i][2] = buf[m++];
radius[i] = buf[m++]; radius[i] = buf[m++];
rmass[i] = buf[m++]; rmass[i] = buf[m++];
phiv[i][0] = buf[m++]; omega[i][0] = buf[m++];
phiv[i][1] = buf[m++]; omega[i][1] = buf[m++];
phiv[i][2] = buf[m++]; omega[i][2] = buf[m++];
} }
} }
@ -436,9 +436,9 @@ int AtomVecGranular::unpack_border_one(int i, double *buf)
v[i][2] = buf[2]; v[i][2] = buf[2];
radius[i] = buf[3]; radius[i] = buf[3];
rmass[i] = buf[4]; rmass[i] = buf[4];
phiv[i][0] = buf[5]; omega[i][0] = buf[5];
phiv[i][1] = buf[6]; omega[i][1] = buf[6];
phiv[i][2] = buf[7]; omega[i][2] = buf[7];
return 8; return 8;
} }
@ -464,12 +464,12 @@ int AtomVecGranular::pack_exchange(int i, double *buf)
buf[m++] = radius[i]; buf[m++] = radius[i];
buf[m++] = density[i]; buf[m++] = density[i];
buf[m++] = rmass[i]; buf[m++] = rmass[i];
buf[m++] = phix[i][0]; buf[m++] = xphi[i][0];
buf[m++] = phix[i][1]; buf[m++] = xphi[i][1];
buf[m++] = phix[i][2]; buf[m++] = xphi[i][2];
buf[m++] = phiv[i][0]; buf[m++] = omega[i][0];
buf[m++] = phiv[i][1]; buf[m++] = omega[i][1];
buf[m++] = phiv[i][2]; buf[m++] = omega[i][2];
if (atom->nextra_grow) if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++) for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@ -501,12 +501,12 @@ int AtomVecGranular::unpack_exchange(double *buf)
radius[nlocal] = buf[m++]; radius[nlocal] = buf[m++];
density[nlocal] = buf[m++]; density[nlocal] = buf[m++];
rmass[nlocal] = buf[m++]; rmass[nlocal] = buf[m++];
phix[nlocal][0] = buf[m++]; xphi[nlocal][0] = buf[m++];
phix[nlocal][1] = buf[m++]; xphi[nlocal][1] = buf[m++];
phix[nlocal][2] = buf[m++]; xphi[nlocal][2] = buf[m++];
phiv[nlocal][0] = buf[m++]; omega[nlocal][0] = buf[m++];
phiv[nlocal][1] = buf[m++]; omega[nlocal][1] = buf[m++];
phiv[nlocal][2] = buf[m++]; omega[nlocal][2] = buf[m++];
if (atom->nextra_grow) if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++) for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@ -569,12 +569,12 @@ int AtomVecGranular::pack_restart(int i, double *buf)
buf[m++] = radius[i]; buf[m++] = radius[i];
buf[m++] = density[i]; buf[m++] = density[i];
buf[m++] = phix[i][0]; buf[m++] = xphi[i][0];
buf[m++] = phix[i][1]; buf[m++] = xphi[i][1];
buf[m++] = phix[i][2]; buf[m++] = xphi[i][2];
buf[m++] = phiv[i][0]; buf[m++] = omega[i][0];
buf[m++] = phiv[i][1]; buf[m++] = omega[i][1];
buf[m++] = phiv[i][2]; buf[m++] = omega[i][2];
if (atom->nextra_restart) if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++) for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
@ -619,12 +619,12 @@ int AtomVecGranular::unpack_restart(double *buf)
else else
rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal]; rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal];
phix[nlocal][0] = buf[m++]; xphi[nlocal][0] = buf[m++];
phix[nlocal][1] = buf[m++]; xphi[nlocal][1] = buf[m++];
phix[nlocal][2] = buf[m++]; xphi[nlocal][2] = buf[m++];
phiv[nlocal][0] = buf[m++]; omega[nlocal][0] = buf[m++];
phiv[nlocal][1] = buf[m++]; omega[nlocal][1] = buf[m++];
phiv[nlocal][2] = buf[m++]; omega[nlocal][2] = buf[m++];
double **extra = atom->extra; double **extra = atom->extra;
if (atom->nextra_store) { if (atom->nextra_store) {
@ -664,12 +664,12 @@ void AtomVecGranular::create_atom(int itype, double *coord, int ihybrid)
radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];
else else
rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal]; rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal];
phix[nlocal][0] = 0.0; xphi[nlocal][0] = 0.0;
phix[nlocal][1] = 0.0; xphi[nlocal][1] = 0.0;
phix[nlocal][2] = 0.0; xphi[nlocal][2] = 0.0;
phiv[nlocal][0] = 0.0; omega[nlocal][0] = 0.0;
phiv[nlocal][1] = 0.0; omega[nlocal][1] = 0.0;
phiv[nlocal][2] = 0.0; omega[nlocal][2] = 0.0;
atom->nlocal++; atom->nlocal++;
} }
@ -711,12 +711,12 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values,
v[nlocal][0] = 0.0; v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0; v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0; v[nlocal][2] = 0.0;
phix[nlocal][0] = 0.0; xphi[nlocal][0] = 0.0;
phix[nlocal][1] = 0.0; xphi[nlocal][1] = 0.0;
phix[nlocal][2] = 0.0; xphi[nlocal][2] = 0.0;
phiv[nlocal][0] = 0.0; omega[nlocal][0] = 0.0;
phiv[nlocal][1] = 0.0; omega[nlocal][1] = 0.0;
phiv[nlocal][2] = 0.0; omega[nlocal][2] = 0.0;
atom->nlocal++; atom->nlocal++;
} }
@ -729,7 +729,8 @@ void AtomVecGranular::data_vel(int m, char *line, int ihybrid)
{ {
int tmp; int tmp;
sscanf(line,"%d %lg %lg %lg %lg %lg %lg", sscanf(line,"%d %lg %lg %lg %lg %lg %lg",
&tmp,&v[m][0],&v[m][1],&v[m][2],&phiv[m][0],&phiv[m][1],&phiv[m][2]); &tmp,&v[m][0],&v[m][1],&v[m][2],
&omega[m][0],&omega[m][1],&omega[m][2]);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -751,9 +752,9 @@ int AtomVecGranular::memory_usage()
if (atom->memcheck("radius")) bytes += nmax * sizeof(double); if (atom->memcheck("radius")) bytes += nmax * sizeof(double);
if (atom->memcheck("density")) bytes += nmax * sizeof(double); if (atom->memcheck("density")) bytes += nmax * sizeof(double);
if (atom->memcheck("rmass")) bytes += nmax * sizeof(double); if (atom->memcheck("rmass")) bytes += nmax * sizeof(double);
if (atom->memcheck("phix")) bytes += nmax*3 * sizeof(double); if (atom->memcheck("xphi")) bytes += nmax*3 * sizeof(double);
if (atom->memcheck("phiv")) bytes += nmax*3 * sizeof(double); if (atom->memcheck("omega")) bytes += nmax*3 * sizeof(double);
if (atom->memcheck("phia")) bytes += nmax*3 * sizeof(double); if (atom->memcheck("torque")) bytes += nmax*3 * sizeof(double);
return bytes; return bytes;
} }

2
src/GRANULAR/atom_vec_granular.h
View File

@ -55,7 +55,7 @@ class AtomVecGranular : public AtomVec {
int *tag,*type,*mask,*image; int *tag,*type,*mask,*image;
double **x,**v,**f; double **x,**v,**f;
double *radius,*density,*rmass; double *radius,*density,*rmass;
double **phix,**phiv,**phia; double **xphi,**omega,**torque;
}; };
} }

9
src/GRANULAR/fix_freeze.cpp
View File

@ -45,6 +45,7 @@ int FixFreeze::setmask()
void FixFreeze::init() void FixFreeze::init()
{ {
// error if more than one freeze fix // error if more than one freeze fix
// because accessed by pair style granular and fix gran/diag
int count = 0; int count = 0;
for (int i = 0; i < modify->nfix; i++) for (int i = 0; i < modify->nfix; i++)
@ -64,7 +65,7 @@ void FixFreeze::setup()
void FixFreeze::post_force(int vflag) void FixFreeze::post_force(int vflag)
{ {
double **f = atom->f; double **f = atom->f;
double **phia = atom->phia; double **torque = atom->torque;
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -73,8 +74,8 @@ void FixFreeze::post_force(int vflag)
f[i][0] = 0.0; f[i][0] = 0.0;
f[i][1] = 0.0; f[i][1] = 0.0;
f[i][2] = 0.0; f[i][2] = 0.0;
phia[i][0] = 0.0; torque[i][0] = 0.0;
phia[i][1] = 0.0; torque[i][1] = 0.0;
phia[i][2] = 0.0; torque[i][2] = 0.0;
} }
} }

24
src/GRANULAR/fix_gran_diag.cpp
View File

@ -467,7 +467,7 @@ void FixGranDiag::stress_no_history()
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **phiv = atom->phiv; double **omega = atom->omega;
double *radius = atom->radius; double *radius = atom->radius;
double *rmass = atom->rmass; double *rmass = atom->rmass;
int *mask = atom->mask; int *mask = atom->mask;
@ -528,9 +528,9 @@ void FixGranDiag::stress_no_history()
// relative rotational velocity // relative rotational velocity
wr1 = radi*phiv[i][0] + radj*phiv[j][0]; wr1 = radi*omega[i][0] + radj*omega[j][0];
wr2 = radi*phiv[i][1] + radj*phiv[j][1]; wr2 = radi*omega[i][1] + radj*omega[j][1];
wr3 = radi*phiv[i][2] + radj*phiv[j][2]; wr3 = radi*omega[i][2] + radj*omega[j][2];
wr1 *= dt/r; wr1 *= dt/r;
wr2 *= dt/r; wr2 *= dt/r;
@ -620,7 +620,7 @@ void FixGranDiag::stress_history()
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **phiv = atom->phiv; double **omega = atom->omega;
double *radius = atom->radius; double *radius = atom->radius;
double *rmass = atom->rmass; double *rmass = atom->rmass;
int *mask = atom->mask; int *mask = atom->mask;
@ -682,9 +682,9 @@ void FixGranDiag::stress_history()
// relative rotational velocity // relative rotational velocity
wr1 = radi*phiv[i][0] + radj*phiv[j][0]; wr1 = radi*omega[i][0] + radj*omega[j][0];
wr2 = radi*phiv[i][1] + radj*phiv[j][1]; wr2 = radi*omega[i][1] + radj*omega[j][1];
wr3 = radi*phiv[i][2] + radj*phiv[j][2]; wr3 = radi*omega[i][2] + radj*omega[j][2];
wr1 *= dt/r; wr1 *= dt/r;
wr2 *= dt/r; wr2 *= dt/r;
@ -803,7 +803,7 @@ void FixGranDiag::stress_hertzian()
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **phiv = atom->phiv; double **omega = atom->omega;
double *radius = atom->radius; double *radius = atom->radius;
double *rmass = atom->rmass; double *rmass = atom->rmass;
int *mask = atom->mask; int *mask = atom->mask;
@ -865,9 +865,9 @@ void FixGranDiag::stress_hertzian()
// relative rotational velocity // relative rotational velocity
wr1 = radi*phiv[i][0] + radj*phiv[j][0]; wr1 = radi*omega[i][0] + radj*omega[j][0];
wr2 = radi*phiv[i][1] + radj*phiv[j][1]; wr2 = radi*omega[i][1] + radj*omega[j][1];
wr3 = radi*phiv[i][2] + radj*phiv[j][2]; wr3 = radi*omega[i][2] + radj*omega[j][2];
wr1 *= dt/r; wr1 *= dt/r;
wr2 *= dt/r; wr2 *= dt/r;

34
src/GRANULAR/fix_nve_gran.cpp
View File

@ -53,10 +53,8 @@ void FixNVEGran::init()
{ {
dtv = update->dt; dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v; dtf = 0.5 * update->dt * force->ftm2v;
if (force->dimension == 3) if (force->dimension == 3) dtfphi = dtf / INERTIA3D;
dtfphi = 0.5 * update->dt * force->ftm2v / INERTIA3D; else dtfphi = dtf / INERTIA2D;
else
dtfphi = 0.5 * update->dt * force->ftm2v / INERTIA2D;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -68,9 +66,9 @@ void FixNVEGran::initial_integrate()
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
double **phix = atom->phix; double **xphi = atom->xphi;
double **phiv = atom->phiv; double **omega = atom->omega;
double **phia = atom->phia; double **torque = atom->torque;
double *rmass = atom->rmass; double *rmass = atom->rmass;
double *radius = atom->radius; double *radius = atom->radius;
int *mask = atom->mask; int *mask = atom->mask;
@ -86,12 +84,12 @@ void FixNVEGran::initial_integrate()
x[i][1] += dtv * v[i][1]; x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2]; x[i][2] += dtv * v[i][2];
dtfm = dtfphi / (radius[i]*radius[i]*rmass[i]); dtfm = dtfphi / (radius[i]*radius[i]*rmass[i]);
phiv[i][0] += dtfm * phia[i][0]; omega[i][0] += dtfm * torque[i][0];
phiv[i][1] += dtfm * phia[i][1]; omega[i][1] += dtfm * torque[i][1];
phiv[i][2] += dtfm * phia[i][2]; omega[i][2] += dtfm * torque[i][2];
phix[i][0] += dtv * phiv[i][0]; xphi[i][0] += dtv * omega[i][0];
phix[i][1] += dtv * phiv[i][1]; xphi[i][1] += dtv * omega[i][1];
phix[i][2] += dtv * phiv[i][2]; xphi[i][2] += dtv * omega[i][2];
} }
} }
} }
@ -104,8 +102,8 @@ void FixNVEGran::final_integrate()
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
double **phiv = atom->phiv; double **omega = atom->omega;
double **phia = atom->phia; double **torque = atom->torque;
double *rmass = atom->rmass; double *rmass = atom->rmass;
double *radius = atom->radius; double *radius = atom->radius;
int *mask = atom->mask; int *mask = atom->mask;
@ -118,9 +116,9 @@ void FixNVEGran::final_integrate()
v[i][1] += dtfm * f[i][1]; v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2]; v[i][2] += dtfm * f[i][2];
dtfm = dtfphi / (radius[i]*radius[i]*rmass[i]); dtfm = dtfphi / (radius[i]*radius[i]*rmass[i]);
phiv[i][0] += dtfm * phia[i][0]; omega[i][0] += dtfm * torque[i][0];
phiv[i][1] += dtfm * phia[i][1]; omega[i][1] += dtfm * torque[i][1];
phiv[i][2] += dtfm * phia[i][2]; omega[i][2] += dtfm * torque[i][2];
} }
} }
} }

52
src/GRANULAR/fix_wall_gran.cpp
View File

@ -235,8 +235,8 @@ void FixWallGran::post_force(int vflag)
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
double **phiv = atom->phiv; double **omega = atom->omega;
double **phia = atom->phia; double **torque = atom->torque;
double *radius = atom->radius; double *radius = atom->radius;
double *rmass = atom->rmass; double *rmass = atom->rmass;
int *mask = atom->mask; int *mask = atom->mask;
@ -282,13 +282,13 @@ void FixWallGran::post_force(int vflag)
} }
} else { } else {
if (pairstyle == NO_HISTORY) if (pairstyle == NO_HISTORY)
no_history(rsq,dx,dy,dz,vwall,v[i],f[i],phiv[i],phia[i], no_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
radius[i],rmass[i]); radius[i],rmass[i]);
else if (pairstyle == HISTORY) else if (pairstyle == HISTORY)
history(rsq,dx,dy,dz,vwall,v[i],f[i],phiv[i],phia[i], history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
radius[i],rmass[i],shear[i]); radius[i],rmass[i],shear[i]);
else if (pairstyle == HERTZIAN) else if (pairstyle == HERTZIAN)
hertzian(rsq,dx,dy,dz,vwall,v[i],f[i],phiv[i],phia[i], hertzian(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
radius[i],rmass[i],shear[i]); radius[i],rmass[i],shear[i]);
} }
} }
@ -299,7 +299,7 @@ void FixWallGran::post_force(int vflag)
void FixWallGran::no_history(double rsq, double dx, double dy, double dz, void FixWallGran::no_history(double rsq, double dx, double dy, double dz,
double *vwall, double *v, double *vwall, double *v,
double *f, double *phiv, double *phia, double *f, double *omega, double *torque,
double radius, double mass) double radius, double mass)
{ {
double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@ -333,9 +333,9 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz,
// relative rotational velocity // relative rotational velocity
wr1 = radius*phiv[0]; wr1 = radius*omega[0];
wr2 = radius*phiv[1]; wr2 = radius*omega[1];
wr3 = radius*phiv[2]; wr3 = radius*omega[2];
wr1 *= dt/r; wr1 *= dt/r;
wr2 *= dt/r; wr2 *= dt/r;
@ -386,16 +386,16 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz,
tor1 = dy*fs3 - dz*fs2; tor1 = dy*fs3 - dz*fs2;
tor2 = dz*fs1 - dx*fs3; tor2 = dz*fs1 - dx*fs3;
tor3 = dx*fs2 - dy*fs1; tor3 = dx*fs2 - dy*fs1;
phia[0] -= radius*tor1; torque[0] -= radius*tor1;
phia[1] -= radius*tor2; torque[1] -= radius*tor2;
phia[2] -= radius*tor3; torque[2] -= radius*tor3;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixWallGran::history(double rsq, double dx, double dy, double dz, void FixWallGran::history(double rsq, double dx, double dy, double dz,
double *vwall, double *v, double *vwall, double *v,
double *f, double *phiv, double *phia, double *f, double *omega, double *torque,
double radius, double mass, double *shear) double radius, double mass, double *shear)
{ {
double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@ -430,9 +430,9 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz,
// relative rotational velocity // relative rotational velocity
wr1 = radius*phiv[0]; wr1 = radius*omega[0];
wr2 = radius*phiv[1]; wr2 = radius*omega[1];
wr3 = radius*phiv[2]; wr3 = radius*omega[2];
wr1 *= dt/r; wr1 *= dt/r;
wr2 *= dt/r; wr2 *= dt/r;
@ -512,16 +512,16 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz,
tor1 = rinv * (dy*fs3 - dz*fs2); tor1 = rinv * (dy*fs3 - dz*fs2);
tor2 = rinv * (dz*fs1 - dx*fs3); tor2 = rinv * (dz*fs1 - dx*fs3);
tor3 = rinv * (dx*fs2 - dy*fs1); tor3 = rinv * (dx*fs2 - dy*fs1);
phia[0] -= radius*tor1; torque[0] -= radius*tor1;
phia[1] -= radius*tor2; torque[1] -= radius*tor2;
phia[2] -= radius*tor3; torque[2] -= radius*tor3;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixWallGran::hertzian(double rsq, double dx, double dy, double dz, void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,
double *vwall, double *v, double *vwall, double *v,
double *f, double *phiv, double *phia, double *f, double *omega, double *torque,
double radius, double mass, double *shear) double radius, double mass, double *shear)
{ {
double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
@ -556,9 +556,9 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,
// relative rotational velocity // relative rotational velocity
wr1 = radius*phiv[0]; wr1 = radius*omega[0];
wr2 = radius*phiv[1]; wr2 = radius*omega[1];
wr3 = radius*phiv[2]; wr3 = radius*omega[2];
wr1 *= dt/r; wr1 *= dt/r;
wr2 *= dt/r; wr2 *= dt/r;
@ -639,9 +639,9 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,
tor1 = dy*fs3 - dz*fs2; tor1 = dy*fs3 - dz*fs2;
tor2 = dz*fs1 - dx*fs3; tor2 = dz*fs1 - dx*fs3;
tor3 = dx*fs2 - dy*fs1; tor3 = dx*fs2 - dy*fs1;
phia[0] -= radius*tor1; torque[0] -= radius*tor1;
phia[1] -= radius*tor2; torque[1] -= radius*tor2;
phia[2] -= radius*tor3; torque[2] -= radius*tor3;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

22
src/GRANULAR/pair_gran_hertzian.cpp
View File

@ -54,8 +54,8 @@ void PairGranHertzian::compute(int eflag, int vflag)
double **f = atom->f; double **f = atom->f;
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **phiv = atom->phiv; double **omega = atom->omega;
double **phia = atom->phia; double **torque = atom->torque;
double *radius = atom->radius; double *radius = atom->radius;
double *rmass = atom->rmass; double *rmass = atom->rmass;
int *mask = atom->mask; int *mask = atom->mask;
@ -122,9 +122,9 @@ void PairGranHertzian::compute(int eflag, int vflag)
// relative rotational velocity // relative rotational velocity
wr1 = radi*phiv[i][0] + radj*phiv[j][0]; wr1 = radi*omega[i][0] + radj*omega[j][0];
wr2 = radi*phiv[i][1] + radj*phiv[j][1]; wr2 = radi*omega[i][1] + radj*omega[j][1];
wr3 = radi*phiv[i][2] + radj*phiv[j][2]; wr3 = radi*omega[i][2] + radj*omega[j][2];
wr1 *= dt/r; wr1 *= dt/r;
wr2 *= dt/r; wr2 *= dt/r;
@ -211,17 +211,17 @@ void PairGranHertzian::compute(int eflag, int vflag)
tor1 = rinv * (dely*fs3 - delz*fs2); tor1 = rinv * (dely*fs3 - delz*fs2);
tor2 = rinv * (delz*fs1 - delx*fs3); tor2 = rinv * (delz*fs1 - delx*fs3);
tor3 = rinv * (delx*fs2 - dely*fs1); tor3 = rinv * (delx*fs2 - dely*fs1);
phia[i][0] -= radi*tor1; torque[i][0] -= radi*tor1;
phia[i][1] -= radi*tor2; torque[i][1] -= radi*tor2;
phia[i][2] -= radi*tor3; torque[i][2] -= radi*tor3;
if (newton_pair || j < nlocal) { if (newton_pair || j < nlocal) {
f[j][0] -= ccelx; f[j][0] -= ccelx;
f[j][1] -= ccely; f[j][1] -= ccely;
f[j][2] -= ccelz; f[j][2] -= ccelz;
phia[j][0] -= radj*tor1; torque[j][0] -= radj*tor1;
phia[j][1] -= radj*tor2; torque[j][1] -= radj*tor2;
phia[j][2] -= radj*tor3; torque[j][2] -= radj*tor3;
} }
} }
} }

22
src/GRANULAR/pair_gran_history.cpp
View File

@ -79,8 +79,8 @@ void PairGranHistory::compute(int eflag, int vflag)
double **f = atom->f; double **f = atom->f;
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **phiv = atom->phiv; double **omega = atom->omega;
double **phia = atom->phia; double **torque = atom->torque;
double *radius = atom->radius; double *radius = atom->radius;
double *rmass = atom->rmass; double *rmass = atom->rmass;
int *mask = atom->mask; int *mask = atom->mask;
@ -147,9 +147,9 @@ void PairGranHistory::compute(int eflag, int vflag)
// relative rotational velocity // relative rotational velocity
wr1 = radi*phiv[i][0] + radj*phiv[j][0]; wr1 = radi*omega[i][0] + radj*omega[j][0];
wr2 = radi*phiv[i][1] + radj*phiv[j][1]; wr2 = radi*omega[i][1] + radj*omega[j][1];
wr3 = radi*phiv[i][2] + radj*phiv[j][2]; wr3 = radi*omega[i][2] + radj*omega[j][2];
wr1 *= dt/r; wr1 *= dt/r;
wr2 *= dt/r; wr2 *= dt/r;
@ -234,17 +234,17 @@ void PairGranHistory::compute(int eflag, int vflag)
tor1 = rinv * (dely*fs3 - delz*fs2); tor1 = rinv * (dely*fs3 - delz*fs2);
tor2 = rinv * (delz*fs1 - delx*fs3); tor2 = rinv * (delz*fs1 - delx*fs3);
tor3 = rinv * (delx*fs2 - dely*fs1); tor3 = rinv * (delx*fs2 - dely*fs1);
phia[i][0] -= radi*tor1; torque[i][0] -= radi*tor1;
phia[i][1] -= radi*tor2; torque[i][1] -= radi*tor2;
phia[i][2] -= radi*tor3; torque[i][2] -= radi*tor3;
if (newton_pair || j < nlocal) { if (newton_pair || j < nlocal) {
f[j][0] -= ccelx; f[j][0] -= ccelx;
f[j][1] -= ccely; f[j][1] -= ccely;
f[j][2] -= ccelz; f[j][2] -= ccelz;
phia[j][0] -= radj*tor1; torque[j][0] -= radj*tor1;
phia[j][1] -= radj*tor2; torque[j][1] -= radj*tor2;
phia[j][2] -= radj*tor3; torque[j][2] -= radj*tor3;
} }
} }
} }

22
src/GRANULAR/pair_gran_no_history.cpp
View File

@ -52,8 +52,8 @@ void PairGranNoHistory::compute(int eflag, int vflag)
double **f = atom->f; double **f = atom->f;
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **phiv = atom->phiv; double **omega = atom->omega;
double **phia = atom->phia; double **torque = atom->torque;
double *radius = atom->radius; double *radius = atom->radius;
double *rmass = atom->rmass; double *rmass = atom->rmass;
int *mask = atom->mask; int *mask = atom->mask;
@ -108,9 +108,9 @@ void PairGranNoHistory::compute(int eflag, int vflag)
// relative rotational velocity // relative rotational velocity
wr1 = radi*phiv[i][0] + radj*phiv[j][0]; wr1 = radi*omega[i][0] + radj*omega[j][0];
wr2 = radi*phiv[i][1] + radj*phiv[j][1]; wr2 = radi*omega[i][1] + radj*omega[j][1];
wr3 = radi*phiv[i][2] + radj*phiv[j][2]; wr3 = radi*omega[i][2] + radj*omega[j][2];
wr1 *= dt/r; wr1 *= dt/r;
wr2 *= dt/r; wr2 *= dt/r;
@ -159,17 +159,17 @@ void PairGranNoHistory::compute(int eflag, int vflag)
tor1 = rinv * (dely*fs3 - delz*fs2); tor1 = rinv * (dely*fs3 - delz*fs2);
tor2 = rinv * (delz*fs1 - delx*fs3); tor2 = rinv * (delz*fs1 - delx*fs3);
tor3 = rinv * (delx*fs2 - dely*fs1); tor3 = rinv * (delx*fs2 - dely*fs1);
phia[i][0] -= radi*tor1; torque[i][0] -= radi*tor1;
phia[i][1] -= radi*tor2; torque[i][1] -= radi*tor2;
phia[i][2] -= radi*tor3; torque[i][2] -= radi*tor3;
if (newton_pair || j < nlocal) { if (newton_pair || j < nlocal) {
f[j][0] -= ccelx; f[j][0] -= ccelx;
f[j][1] -= ccely; f[j][1] -= ccely;
f[j][2] -= ccelz; f[j][2] -= ccelz;
phia[j][0] -= radj*tor1; torque[j][0] -= radj*tor1;
phia[j][1] -= radj*tor2; torque[j][1] -= radj*tor2;
phia[j][2] -= radj*tor3; torque[j][2] -= radj*tor3;
} }
} }
} }