diff --git a/src/GRANULAR/atom_vec_granular.cpp b/src/GRANULAR/atom_vec_granular.cpp index ad352f3a03..9bb5f05519 100644 --- a/src/GRANULAR/atom_vec_granular.cpp +++ b/src/GRANULAR/atom_vec_granular.cpp @@ -73,12 +73,12 @@ void AtomVecGranular::grow(int n) rmass = atom->rmass = (double *) memory->srealloc(atom->rmass,nmax*sizeof(double),"atom:rmass"); - phix = atom->phix = - memory->grow_2d_double_array(atom->phix,nmax,3,"atom:phix"); - phiv = atom->phiv = - memory->grow_2d_double_array(atom->phiv,nmax,3,"atom:phiv"); - phia = atom->phia = - memory->grow_2d_double_array(atom->phia,nmax,3,"atom:phia"); + xphi = atom->xphi = + memory->grow_2d_double_array(atom->xphi,nmax,3,"atom:xphi"); + omega = atom->omega = + memory->grow_2d_double_array(atom->omega,nmax,3,"atom:omega"); + torque = atom->torque = + memory->grow_2d_double_array(atom->torque,nmax,3,"atom:torque"); if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -100,9 +100,9 @@ void AtomVecGranular::reset_ptrs() radius = atom->radius; density = atom->density; rmass = atom->rmass; - phix = atom->phix; - phiv = atom->phiv; - phia = atom->phia; + xphi = atom->xphi; + omega = atom->omega; + torque = atom->torque; } /* ---------------------------------------------------------------------- @@ -115,12 +115,12 @@ void AtomVecGranular::zero_owned(int i) radius[i] = 0.0; density[i] = 0.0; rmass[i] = 0.0; - phix[i][0] = 0.0; - phix[i][1] = 0.0; - phix[i][2] = 0.0; - phiv[i][0] = 0.0; - phiv[i][1] = 0.0; - phiv[i][2] = 0.0; + xphi[i][0] = 0.0; + xphi[i][1] = 0.0; + xphi[i][2] = 0.0; + omega[i][0] = 0.0; + omega[i][1] = 0.0; + omega[i][2] = 0.0; } /* ---------------------------------------------------------------------- @@ -139,9 +139,9 @@ void AtomVecGranular::zero_ghost(int n, int first) v[i][2] = 0.0; radius[i] = 0.0; rmass[i] = 0.0; - phiv[i][0] = 0.0; - phiv[i][1] = 0.0; - phiv[i][2] = 0.0; + omega[i][0] = 0.0; + omega[i][1] = 0.0; + omega[i][2] = 0.0; } } @@ -163,12 +163,12 @@ void AtomVecGranular::copy(int i, int j) radius[j] = radius[i]; density[j] = density[i]; rmass[j] = rmass[i]; - phix[j][0] = phix[i][0]; - phix[j][1] = phix[i][1]; - phix[j][2] = phix[i][2]; - phiv[j][0] = phiv[i][0]; - phiv[j][1] = phiv[i][1]; - phiv[j][2] = phiv[i][2]; + xphi[j][0] = xphi[i][0]; + xphi[j][1] = xphi[i][1]; + xphi[j][2] = xphi[i][2]; + omega[j][0] = omega[i][0]; + omega[j][1] = omega[i][1]; + omega[j][2] = omega[i][2]; if (atom->nextra_grow) 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][1]; buf[m++] = v[j][2]; - buf[m++] = phiv[j][0]; - buf[m++] = phiv[j][1]; - buf[m++] = phiv[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; } } else { 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][1]; buf[m++] = v[j][2]; - buf[m++] = phiv[j][0]; - buf[m++] = phiv[j][1]; - buf[m++] = phiv[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; } } return m; @@ -230,9 +230,9 @@ int AtomVecGranular::pack_comm_one(int i, double *buf) buf[0] = v[i][0]; buf[1] = v[i][1]; buf[2] = v[i][2]; - buf[3] = phiv[i][0]; - buf[4] = phiv[i][1]; - buf[5] = phiv[i][2]; + buf[3] = omega[i][0]; + buf[4] = omega[i][1]; + buf[5] = omega[i][2]; return 6; } @@ -251,9 +251,9 @@ void AtomVecGranular::unpack_comm(int n, int first, double *buf) v[i][0] = buf[m++]; v[i][1] = buf[m++]; v[i][2] = buf[m++]; - phiv[i][0] = buf[m++]; - phiv[i][1] = buf[m++]; - phiv[i][2] = buf[m++]; + omega[i][0] = buf[m++]; + omega[i][1] = 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][1] = buf[1]; v[i][2] = buf[2]; - phiv[i][0] = buf[3]; - phiv[i][1] = buf[4]; - phiv[i][2] = buf[5]; + omega[i][0] = buf[3]; + omega[i][1] = buf[4]; + omega[i][2] = buf[5]; 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][1]; buf[m++] = f[i][2]; - buf[m++] = phia[i][0]; - buf[m++] = phia[i][1]; - buf[m++] = phia[i][2]; + buf[m++] = torque[i][0]; + buf[m++] = torque[i][1]; + buf[m++] = torque[i][2]; } 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) { - buf[0] = phia[i][0]; - buf[1] = phia[i][1]; - buf[2] = phia[i][2]; + buf[0] = torque[i][0]; + buf[1] = torque[i][1]; + buf[2] = torque[i][2]; return 3; } @@ -311,9 +311,9 @@ void AtomVecGranular::unpack_reverse(int n, int *list, double *buf) f[j][0] += buf[m++]; f[j][1] += buf[m++]; f[j][2] += buf[m++]; - phia[j][0] += buf[m++]; - phia[j][1] += buf[m++]; - phia[j][2] += buf[m++]; + torque[j][0] += buf[m++]; + torque[j][1] += 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) { - phia[i][0] = buf[0]; - phia[i][1] = buf[1]; - phia[i][2] = buf[2]; + torque[i][0] = buf[0]; + torque[i][1] = buf[1]; + torque[i][2] = buf[2]; return 3; } @@ -350,9 +350,9 @@ int AtomVecGranular::pack_border(int n, int *list, double *buf, buf[m++] = v[j][2]; buf[m++] = radius[j]; buf[m++] = rmass[j]; - buf[m++] = phiv[j][0]; - buf[m++] = phiv[j][1]; - buf[m++] = phiv[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; } } else { 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++] = radius[j]; buf[m++] = rmass[j]; - buf[m++] = phiv[j][0]; - buf[m++] = phiv[j][1]; - buf[m++] = phiv[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; } } return m; @@ -394,9 +394,9 @@ int AtomVecGranular::pack_border_one(int i, double *buf) buf[2] = v[i][2]; buf[3] = radius[i]; buf[4] = rmass[i]; - buf[5] = phiv[i][0]; - buf[6] = phiv[i][1]; - buf[7] = phiv[i][2]; + buf[5] = omega[i][0]; + buf[6] = omega[i][1]; + buf[7] = omega[i][2]; return 8; } @@ -421,9 +421,9 @@ void AtomVecGranular::unpack_border(int n, int first, double *buf) v[i][2] = buf[m++]; radius[i] = buf[m++]; rmass[i] = buf[m++]; - phiv[i][0] = buf[m++]; - phiv[i][1] = buf[m++]; - phiv[i][2] = buf[m++]; + omega[i][0] = buf[m++]; + omega[i][1] = 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]; radius[i] = buf[3]; rmass[i] = buf[4]; - phiv[i][0] = buf[5]; - phiv[i][1] = buf[6]; - phiv[i][2] = buf[7]; + omega[i][0] = buf[5]; + omega[i][1] = buf[6]; + omega[i][2] = buf[7]; return 8; } @@ -464,12 +464,12 @@ int AtomVecGranular::pack_exchange(int i, double *buf) buf[m++] = radius[i]; buf[m++] = density[i]; buf[m++] = rmass[i]; - buf[m++] = phix[i][0]; - buf[m++] = phix[i][1]; - buf[m++] = phix[i][2]; - buf[m++] = phiv[i][0]; - buf[m++] = phiv[i][1]; - buf[m++] = phiv[i][2]; + buf[m++] = xphi[i][0]; + buf[m++] = xphi[i][1]; + buf[m++] = xphi[i][2]; + buf[m++] = omega[i][0]; + buf[m++] = omega[i][1]; + buf[m++] = omega[i][2]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -501,12 +501,12 @@ int AtomVecGranular::unpack_exchange(double *buf) radius[nlocal] = buf[m++]; density[nlocal] = buf[m++]; rmass[nlocal] = buf[m++]; - phix[nlocal][0] = buf[m++]; - phix[nlocal][1] = buf[m++]; - phix[nlocal][2] = buf[m++]; - phiv[nlocal][0] = buf[m++]; - phiv[nlocal][1] = buf[m++]; - phiv[nlocal][2] = buf[m++]; + xphi[nlocal][0] = buf[m++]; + xphi[nlocal][1] = buf[m++]; + xphi[nlocal][2] = buf[m++]; + omega[nlocal][0] = buf[m++]; + omega[nlocal][1] = buf[m++]; + omega[nlocal][2] = buf[m++]; if (atom->nextra_grow) 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++] = density[i]; - buf[m++] = phix[i][0]; - buf[m++] = phix[i][1]; - buf[m++] = phix[i][2]; - buf[m++] = phiv[i][0]; - buf[m++] = phiv[i][1]; - buf[m++] = phiv[i][2]; + buf[m++] = xphi[i][0]; + buf[m++] = xphi[i][1]; + buf[m++] = xphi[i][2]; + buf[m++] = omega[i][0]; + buf[m++] = omega[i][1]; + buf[m++] = omega[i][2]; if (atom->nextra_restart) for (int iextra = 0; iextra < atom->nextra_restart; iextra++) @@ -619,12 +619,12 @@ int AtomVecGranular::unpack_restart(double *buf) else rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal]; - phix[nlocal][0] = buf[m++]; - phix[nlocal][1] = buf[m++]; - phix[nlocal][2] = buf[m++]; - phiv[nlocal][0] = buf[m++]; - phiv[nlocal][1] = buf[m++]; - phiv[nlocal][2] = buf[m++]; + xphi[nlocal][0] = buf[m++]; + xphi[nlocal][1] = buf[m++]; + xphi[nlocal][2] = buf[m++]; + omega[nlocal][0] = buf[m++]; + omega[nlocal][1] = buf[m++]; + omega[nlocal][2] = buf[m++]; double **extra = atom->extra; 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]; else rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal]; - phix[nlocal][0] = 0.0; - phix[nlocal][1] = 0.0; - phix[nlocal][2] = 0.0; - phiv[nlocal][0] = 0.0; - phiv[nlocal][1] = 0.0; - phiv[nlocal][2] = 0.0; + xphi[nlocal][0] = 0.0; + xphi[nlocal][1] = 0.0; + xphi[nlocal][2] = 0.0; + omega[nlocal][0] = 0.0; + omega[nlocal][1] = 0.0; + omega[nlocal][2] = 0.0; atom->nlocal++; } @@ -711,12 +711,12 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values, v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; v[nlocal][2] = 0.0; - phix[nlocal][0] = 0.0; - phix[nlocal][1] = 0.0; - phix[nlocal][2] = 0.0; - phiv[nlocal][0] = 0.0; - phiv[nlocal][1] = 0.0; - phiv[nlocal][2] = 0.0; + xphi[nlocal][0] = 0.0; + xphi[nlocal][1] = 0.0; + xphi[nlocal][2] = 0.0; + omega[nlocal][0] = 0.0; + omega[nlocal][1] = 0.0; + omega[nlocal][2] = 0.0; atom->nlocal++; } @@ -729,7 +729,8 @@ void AtomVecGranular::data_vel(int m, char *line, int ihybrid) { int tmp; 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("density")) bytes += nmax * sizeof(double); if (atom->memcheck("rmass")) bytes += nmax * sizeof(double); - if (atom->memcheck("phix")) bytes += nmax*3 * sizeof(double); - if (atom->memcheck("phiv")) bytes += nmax*3 * sizeof(double); - if (atom->memcheck("phia")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("xphi")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("omega")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("torque")) bytes += nmax*3 * sizeof(double); return bytes; } diff --git a/src/GRANULAR/atom_vec_granular.h b/src/GRANULAR/atom_vec_granular.h index 25c26a6755..61e6d7d8f4 100644 --- a/src/GRANULAR/atom_vec_granular.h +++ b/src/GRANULAR/atom_vec_granular.h @@ -55,7 +55,7 @@ class AtomVecGranular : public AtomVec { int *tag,*type,*mask,*image; double **x,**v,**f; double *radius,*density,*rmass; - double **phix,**phiv,**phia; + double **xphi,**omega,**torque; }; } diff --git a/src/GRANULAR/fix_freeze.cpp b/src/GRANULAR/fix_freeze.cpp index b264fe66bc..34e8a4014c 100644 --- a/src/GRANULAR/fix_freeze.cpp +++ b/src/GRANULAR/fix_freeze.cpp @@ -45,6 +45,7 @@ int FixFreeze::setmask() void FixFreeze::init() { // error if more than one freeze fix + // because accessed by pair style granular and fix gran/diag int count = 0; for (int i = 0; i < modify->nfix; i++) @@ -64,7 +65,7 @@ void FixFreeze::setup() void FixFreeze::post_force(int vflag) { double **f = atom->f; - double **phia = atom->phia; + double **torque = atom->torque; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -73,8 +74,8 @@ void FixFreeze::post_force(int vflag) f[i][0] = 0.0; f[i][1] = 0.0; f[i][2] = 0.0; - phia[i][0] = 0.0; - phia[i][1] = 0.0; - phia[i][2] = 0.0; + torque[i][0] = 0.0; + torque[i][1] = 0.0; + torque[i][2] = 0.0; } } diff --git a/src/GRANULAR/fix_gran_diag.cpp b/src/GRANULAR/fix_gran_diag.cpp index 8d7ce98cca..f606d5ee20 100644 --- a/src/GRANULAR/fix_gran_diag.cpp +++ b/src/GRANULAR/fix_gran_diag.cpp @@ -467,7 +467,7 @@ void FixGranDiag::stress_no_history() double **x = atom->x; double **v = atom->v; - double **phiv = atom->phiv; + double **omega = atom->omega; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; @@ -528,9 +528,9 @@ void FixGranDiag::stress_no_history() // relative rotational velocity - wr1 = radi*phiv[i][0] + radj*phiv[j][0]; - wr2 = radi*phiv[i][1] + radj*phiv[j][1]; - wr3 = radi*phiv[i][2] + radj*phiv[j][2]; + wr1 = radi*omega[i][0] + radj*omega[j][0]; + wr2 = radi*omega[i][1] + radj*omega[j][1]; + wr3 = radi*omega[i][2] + radj*omega[j][2]; wr1 *= dt/r; wr2 *= dt/r; @@ -620,7 +620,7 @@ void FixGranDiag::stress_history() double **x = atom->x; double **v = atom->v; - double **phiv = atom->phiv; + double **omega = atom->omega; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; @@ -682,9 +682,9 @@ void FixGranDiag::stress_history() // relative rotational velocity - wr1 = radi*phiv[i][0] + radj*phiv[j][0]; - wr2 = radi*phiv[i][1] + radj*phiv[j][1]; - wr3 = radi*phiv[i][2] + radj*phiv[j][2]; + wr1 = radi*omega[i][0] + radj*omega[j][0]; + wr2 = radi*omega[i][1] + radj*omega[j][1]; + wr3 = radi*omega[i][2] + radj*omega[j][2]; wr1 *= dt/r; wr2 *= dt/r; @@ -803,7 +803,7 @@ void FixGranDiag::stress_hertzian() double **x = atom->x; double **v = atom->v; - double **phiv = atom->phiv; + double **omega = atom->omega; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; @@ -865,9 +865,9 @@ void FixGranDiag::stress_hertzian() // relative rotational velocity - wr1 = radi*phiv[i][0] + radj*phiv[j][0]; - wr2 = radi*phiv[i][1] + radj*phiv[j][1]; - wr3 = radi*phiv[i][2] + radj*phiv[j][2]; + wr1 = radi*omega[i][0] + radj*omega[j][0]; + wr2 = radi*omega[i][1] + radj*omega[j][1]; + wr3 = radi*omega[i][2] + radj*omega[j][2]; wr1 *= dt/r; wr2 *= dt/r; diff --git a/src/GRANULAR/fix_nve_gran.cpp b/src/GRANULAR/fix_nve_gran.cpp index f98a9626ad..4a52af9ee9 100644 --- a/src/GRANULAR/fix_nve_gran.cpp +++ b/src/GRANULAR/fix_nve_gran.cpp @@ -53,10 +53,8 @@ void FixNVEGran::init() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; - if (force->dimension == 3) - dtfphi = 0.5 * update->dt * force->ftm2v / INERTIA3D; - else - dtfphi = 0.5 * update->dt * force->ftm2v / INERTIA2D; + if (force->dimension == 3) dtfphi = dtf / INERTIA3D; + else dtfphi = dtf / INERTIA2D; } /* ---------------------------------------------------------------------- */ @@ -68,9 +66,9 @@ void FixNVEGran::initial_integrate() double **x = atom->x; double **v = atom->v; double **f = atom->f; - double **phix = atom->phix; - double **phiv = atom->phiv; - double **phia = atom->phia; + double **xphi = atom->xphi; + double **omega = atom->omega; + double **torque = atom->torque; double *rmass = atom->rmass; double *radius = atom->radius; int *mask = atom->mask; @@ -86,12 +84,12 @@ void FixNVEGran::initial_integrate() x[i][1] += dtv * v[i][1]; x[i][2] += dtv * v[i][2]; dtfm = dtfphi / (radius[i]*radius[i]*rmass[i]); - phiv[i][0] += dtfm * phia[i][0]; - phiv[i][1] += dtfm * phia[i][1]; - phiv[i][2] += dtfm * phia[i][2]; - phix[i][0] += dtv * phiv[i][0]; - phix[i][1] += dtv * phiv[i][1]; - phix[i][2] += dtv * phiv[i][2]; + omega[i][0] += dtfm * torque[i][0]; + omega[i][1] += dtfm * torque[i][1]; + omega[i][2] += dtfm * torque[i][2]; + xphi[i][0] += dtv * omega[i][0]; + xphi[i][1] += dtv * omega[i][1]; + xphi[i][2] += dtv * omega[i][2]; } } } @@ -104,8 +102,8 @@ void FixNVEGran::final_integrate() double **v = atom->v; double **f = atom->f; - double **phiv = atom->phiv; - double **phia = atom->phia; + double **omega = atom->omega; + double **torque = atom->torque; double *rmass = atom->rmass; double *radius = atom->radius; int *mask = atom->mask; @@ -118,9 +116,9 @@ void FixNVEGran::final_integrate() v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; dtfm = dtfphi / (radius[i]*radius[i]*rmass[i]); - phiv[i][0] += dtfm * phia[i][0]; - phiv[i][1] += dtfm * phia[i][1]; - phiv[i][2] += dtfm * phia[i][2]; + omega[i][0] += dtfm * torque[i][0]; + omega[i][1] += dtfm * torque[i][1]; + omega[i][2] += dtfm * torque[i][2]; } } } diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 7d5116d1e1..a9d329f8d4 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -235,8 +235,8 @@ void FixWallGran::post_force(int vflag) double **x = atom->x; double **v = atom->v; double **f = atom->f; - double **phiv = atom->phiv; - double **phia = atom->phia; + double **omega = atom->omega; + double **torque = atom->torque; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; @@ -282,13 +282,13 @@ void FixWallGran::post_force(int vflag) } } else { 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]); 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]); 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]); } } @@ -299,7 +299,7 @@ void FixWallGran::post_force(int vflag) void FixWallGran::no_history(double rsq, double dx, double dy, double dz, double *vwall, double *v, - double *f, double *phiv, double *phia, + double *f, double *omega, double *torque, double radius, double mass) { 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 - wr1 = radius*phiv[0]; - wr2 = radius*phiv[1]; - wr3 = radius*phiv[2]; + wr1 = radius*omega[0]; + wr2 = radius*omega[1]; + wr3 = radius*omega[2]; wr1 *= 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; tor2 = dz*fs1 - dx*fs3; tor3 = dx*fs2 - dy*fs1; - phia[0] -= radius*tor1; - phia[1] -= radius*tor2; - phia[2] -= radius*tor3; + torque[0] -= radius*tor1; + torque[1] -= radius*tor2; + torque[2] -= radius*tor3; } /* ---------------------------------------------------------------------- */ void FixWallGran::history(double rsq, double dx, double dy, double dz, double *vwall, double *v, - double *f, double *phiv, double *phia, + double *f, double *omega, double *torque, double radius, double mass, double *shear) { 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 - wr1 = radius*phiv[0]; - wr2 = radius*phiv[1]; - wr3 = radius*phiv[2]; + wr1 = radius*omega[0]; + wr2 = radius*omega[1]; + wr3 = radius*omega[2]; wr1 *= 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); tor2 = rinv * (dz*fs1 - dx*fs3); tor3 = rinv * (dx*fs2 - dy*fs1); - phia[0] -= radius*tor1; - phia[1] -= radius*tor2; - phia[2] -= radius*tor3; + torque[0] -= radius*tor1; + torque[1] -= radius*tor2; + torque[2] -= radius*tor3; } /* ---------------------------------------------------------------------- */ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz, double *vwall, double *v, - double *f, double *phiv, double *phia, + double *f, double *omega, double *torque, double radius, double mass, double *shear) { 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 - wr1 = radius*phiv[0]; - wr2 = radius*phiv[1]; - wr3 = radius*phiv[2]; + wr1 = radius*omega[0]; + wr2 = radius*omega[1]; + wr3 = radius*omega[2]; wr1 *= 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; tor2 = dz*fs1 - dx*fs3; tor3 = dx*fs2 - dy*fs1; - phia[0] -= radius*tor1; - phia[1] -= radius*tor2; - phia[2] -= radius*tor3; + torque[0] -= radius*tor1; + torque[1] -= radius*tor2; + torque[2] -= radius*tor3; } /* ---------------------------------------------------------------------- diff --git a/src/GRANULAR/pair_gran_hertzian.cpp b/src/GRANULAR/pair_gran_hertzian.cpp index 87ba3d248b..2dad7b077e 100644 --- a/src/GRANULAR/pair_gran_hertzian.cpp +++ b/src/GRANULAR/pair_gran_hertzian.cpp @@ -54,8 +54,8 @@ void PairGranHertzian::compute(int eflag, int vflag) double **f = atom->f; double **x = atom->x; double **v = atom->v; - double **phiv = atom->phiv; - double **phia = atom->phia; + double **omega = atom->omega; + double **torque = atom->torque; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; @@ -122,9 +122,9 @@ void PairGranHertzian::compute(int eflag, int vflag) // relative rotational velocity - wr1 = radi*phiv[i][0] + radj*phiv[j][0]; - wr2 = radi*phiv[i][1] + radj*phiv[j][1]; - wr3 = radi*phiv[i][2] + radj*phiv[j][2]; + wr1 = radi*omega[i][0] + radj*omega[j][0]; + wr2 = radi*omega[i][1] + radj*omega[j][1]; + wr3 = radi*omega[i][2] + radj*omega[j][2]; wr1 *= dt/r; wr2 *= dt/r; @@ -211,17 +211,17 @@ void PairGranHertzian::compute(int eflag, int vflag) tor1 = rinv * (dely*fs3 - delz*fs2); tor2 = rinv * (delz*fs1 - delx*fs3); tor3 = rinv * (delx*fs2 - dely*fs1); - phia[i][0] -= radi*tor1; - phia[i][1] -= radi*tor2; - phia[i][2] -= radi*tor3; + torque[i][0] -= radi*tor1; + torque[i][1] -= radi*tor2; + torque[i][2] -= radi*tor3; if (newton_pair || j < nlocal) { f[j][0] -= ccelx; f[j][1] -= ccely; f[j][2] -= ccelz; - phia[j][0] -= radj*tor1; - phia[j][1] -= radj*tor2; - phia[j][2] -= radj*tor3; + torque[j][0] -= radj*tor1; + torque[j][1] -= radj*tor2; + torque[j][2] -= radj*tor3; } } } diff --git a/src/GRANULAR/pair_gran_history.cpp b/src/GRANULAR/pair_gran_history.cpp index 541b0e9950..8b6666212f 100644 --- a/src/GRANULAR/pair_gran_history.cpp +++ b/src/GRANULAR/pair_gran_history.cpp @@ -79,8 +79,8 @@ void PairGranHistory::compute(int eflag, int vflag) double **f = atom->f; double **x = atom->x; double **v = atom->v; - double **phiv = atom->phiv; - double **phia = atom->phia; + double **omega = atom->omega; + double **torque = atom->torque; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; @@ -147,9 +147,9 @@ void PairGranHistory::compute(int eflag, int vflag) // relative rotational velocity - wr1 = radi*phiv[i][0] + radj*phiv[j][0]; - wr2 = radi*phiv[i][1] + radj*phiv[j][1]; - wr3 = radi*phiv[i][2] + radj*phiv[j][2]; + wr1 = radi*omega[i][0] + radj*omega[j][0]; + wr2 = radi*omega[i][1] + radj*omega[j][1]; + wr3 = radi*omega[i][2] + radj*omega[j][2]; wr1 *= dt/r; wr2 *= dt/r; @@ -234,17 +234,17 @@ void PairGranHistory::compute(int eflag, int vflag) tor1 = rinv * (dely*fs3 - delz*fs2); tor2 = rinv * (delz*fs1 - delx*fs3); tor3 = rinv * (delx*fs2 - dely*fs1); - phia[i][0] -= radi*tor1; - phia[i][1] -= radi*tor2; - phia[i][2] -= radi*tor3; + torque[i][0] -= radi*tor1; + torque[i][1] -= radi*tor2; + torque[i][2] -= radi*tor3; if (newton_pair || j < nlocal) { f[j][0] -= ccelx; f[j][1] -= ccely; f[j][2] -= ccelz; - phia[j][0] -= radj*tor1; - phia[j][1] -= radj*tor2; - phia[j][2] -= radj*tor3; + torque[j][0] -= radj*tor1; + torque[j][1] -= radj*tor2; + torque[j][2] -= radj*tor3; } } } diff --git a/src/GRANULAR/pair_gran_no_history.cpp b/src/GRANULAR/pair_gran_no_history.cpp index 4542e26048..01e8a1f3a3 100644 --- a/src/GRANULAR/pair_gran_no_history.cpp +++ b/src/GRANULAR/pair_gran_no_history.cpp @@ -52,8 +52,8 @@ void PairGranNoHistory::compute(int eflag, int vflag) double **f = atom->f; double **x = atom->x; double **v = atom->v; - double **phiv = atom->phiv; - double **phia = atom->phia; + double **omega = atom->omega; + double **torque = atom->torque; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; @@ -108,9 +108,9 @@ void PairGranNoHistory::compute(int eflag, int vflag) // relative rotational velocity - wr1 = radi*phiv[i][0] + radj*phiv[j][0]; - wr2 = radi*phiv[i][1] + radj*phiv[j][1]; - wr3 = radi*phiv[i][2] + radj*phiv[j][2]; + wr1 = radi*omega[i][0] + radj*omega[j][0]; + wr2 = radi*omega[i][1] + radj*omega[j][1]; + wr3 = radi*omega[i][2] + radj*omega[j][2]; wr1 *= dt/r; wr2 *= dt/r; @@ -159,17 +159,17 @@ void PairGranNoHistory::compute(int eflag, int vflag) tor1 = rinv * (dely*fs3 - delz*fs2); tor2 = rinv * (delz*fs1 - delx*fs3); tor3 = rinv * (delx*fs2 - dely*fs1); - phia[i][0] -= radi*tor1; - phia[i][1] -= radi*tor2; - phia[i][2] -= radi*tor3; + torque[i][0] -= radi*tor1; + torque[i][1] -= radi*tor2; + torque[i][2] -= radi*tor3; if (newton_pair || j < nlocal) { f[j][0] -= ccelx; f[j][1] -= ccely; f[j][2] -= ccelz; - phia[j][0] -= radj*tor1; - phia[j][1] -= radj*tor2; - phia[j][2] -= radj*tor3; + torque[j][0] -= radj*tor1; + torque[j][1] -= radj*tor2; + torque[j][2] -= radj*tor3; } } }