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@8385 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2012-06-22 16:21:45 +00:00
parent b33ba405a3
commit 8f68628a20
5 changed files with 194 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

69
src/GRANULAR/pair_gran_hertz_history.cpp
View File

@ -23,7 +23,10 @@
#include "atom.h"
#include "update.h"
#include "force.h"
#include "fix_rigid.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "comm.h"
#include "error.h"
using namespace LAMMPS_NS;
@ -43,7 +46,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3;
double vtr1,vtr2,vtr3,vrel;
double meff,damp,ccel,tor1,tor2,tor3;
double mi,mj,meff,damp,ccel,tor1,tor2,tor3;
double fn,fs,fs1,fs2,fs3;
double shrmag,rsht,polyhertz;
int *ilist,*jlist,*numneigh,**firstneigh;
@ -57,6 +60,13 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
int shearupdate = 1;
if (update->setupflag) shearupdate = 0;
// update body ptr and values for ghost atoms if using FixRigid masses
if (fix_rigid && neighbor->ago == 0) {
body = fix_rigid->body;
comm->forward_comm_pair(this);
}
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
@ -140,19 +150,27 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
// normal force = Hertzian contact + normal velocity damping
// meff = effective mass of pair of particles
// if I or J part of rigid body, use body mass
// if I or J is frozen, meff is other particle
if (rmass) {
meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
if (mask[i] & freeze_group_bit) meff = rmass[j];
if (mask[j] & freeze_group_bit) meff = rmass[i];
mi = rmass[i];
mj = rmass[j];
} else {
itype = type[i];
jtype = type[j];
meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
if (mask[i] & freeze_group_bit) meff = mass[jtype];
if (mask[j] & freeze_group_bit) meff = mass[itype];
mi = mass[type[i]];
mj = mass[type[j]];
}
if (fix_rigid) {
if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]];
if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]];
}
meff = mi*mj / (mi+mj);
if (mask[i] & freeze_group_bit) meff = mj;
if (mask[j] & freeze_group_bit) meff = mi;
// normal force = Hertzian contact + normal velocity damping
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radsum-r)*rinv - damp;
@ -286,7 +304,7 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype,
double radi,radj,radsum;
double r,rinv,rsqinv,delx,dely,delz;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;
double meff,damp,ccel,polyhertz;
double mi,mj,meff,damp,ccel,polyhertz;
double vtr1,vtr2,vtr3,vrel,shrmag,rsht;
double fs1,fs2,fs3,fs,fn;
@ -337,21 +355,36 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype,
wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
// normal force = Hertzian contact + normal velocity damping
// meff = effective mass of pair of particles
// if I or J part of rigid body, use body mass
// if I or J is frozen, meff is other particle
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
if (rmass) {
meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
if (mask[i] & freeze_group_bit) meff = rmass[j];
if (mask[j] & freeze_group_bit) meff = rmass[i];
mi = rmass[i];
mj = rmass[j];
} else {
meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
if (mask[i] & freeze_group_bit) meff = mass[jtype];
if (mask[j] & freeze_group_bit) meff = mass[itype];
mi = mass[type[i]];
mj = mass[type[j]];
}
if (fix_rigid) {
// NOTE: need to make sure ghost atoms have updated body?
// depends on where single() is called from
body = fix_rigid->body;
if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]];
if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]];
}
meff = mi*mj / (mi+mj);
if (mask[i] & freeze_group_bit) meff = mj;
if (mask[j] & freeze_group_bit) meff = mi;
// normal force = Hertzian contact + normal velocity damping
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radsum-r)*rinv - damp;

68
src/GRANULAR/pair_gran_hooke.cpp
View File

@ -21,7 +21,10 @@
#include "pair_gran_hooke.h"
#include "atom.h"
#include "force.h"
#include "fix_rigid.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "comm.h"
using namespace LAMMPS_NS;
@ -43,13 +46,20 @@ void PairGranHooke::compute(int eflag, int vflag)
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3;
double vtr1,vtr2,vtr3,vrel;
double meff,damp,ccel,tor1,tor2,tor3;
double mi,mj,meff,damp,ccel,tor1,tor2,tor3;
double fn,fs,ft,fs1,fs2,fs3;
int *ilist,*jlist,*numneigh,**firstneigh;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
// update body ptr and values for ghost atoms if using FixRigid masses
if (fix_rigid && neighbor->ago == 0) {
body = fix_rigid->body;
comm->forward_comm_pair(this);
}
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
@ -120,19 +130,27 @@ void PairGranHooke::compute(int eflag, int vflag)
wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
// normal forces = Hookian contact + normal velocity damping
// meff = effective mass of pair of particles
// if I or J part of rigid body, use body mass
// if I or J is frozen, meff is other particle
if (rmass) {
meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
if (mask[i] & freeze_group_bit) meff = rmass[j];
if (mask[j] & freeze_group_bit) meff = rmass[i];
mi = rmass[i];
mj = rmass[j];
} else {
itype = type[i];
jtype = type[j];
meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
if (mask[i] & freeze_group_bit) meff = mass[jtype];
if (mask[j] & freeze_group_bit) meff = mass[itype];
mi = mass[type[i]];
mj = mass[type[j]];
}
if (fix_rigid) {
if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]];
if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]];
}
meff = mi*mj / (mi+mj);
if (mask[i] & freeze_group_bit) meff = mj;
if (mask[j] & freeze_group_bit) meff = mi;
// normal forces = Hookian contact + normal velocity damping
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radsum-r)*rinv - damp;
@ -202,7 +220,7 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq,
double delx,dely,delz;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;
double vtr1,vtr2,vtr3,vrel;
double meff,damp,ccel;
double mi,mj,meff,damp,ccel;
double fn,fs,ft,fs1,fs2,fs3;
double *radius = atom->radius;
@ -254,21 +272,35 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq,
wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
// normal forces = Hookian contact + normal velocity damping
// meff = effective mass of pair of particles
// if I or J part of rigid body, use body mass
// if I or J is frozen, meff is other particle
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
if (rmass) {
meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
if (mask[i] & freeze_group_bit) meff = rmass[j];
if (mask[j] & freeze_group_bit) meff = rmass[i];
mi = rmass[i];
mj = rmass[j];
} else {
meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
if (mask[i] & freeze_group_bit) meff = mass[jtype];
if (mask[j] & freeze_group_bit) meff = mass[itype];
mi = mass[type[i]];
mj = mass[type[j]];
}
if (fix_rigid) {
// NOTE: need to make sure ghost atoms have updated body?
// depends on where single() is called from
body = fix_rigid->body;
if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]];
if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]];
}
meff = mi*mj / (mi+mj);
if (mask[i] & freeze_group_bit) meff = mj;
if (mask[j] & freeze_group_bit) meff = mi;
// normal forces = Hookian contact + normal velocity damping
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radsum-r)*rinv - damp;

106
src/GRANULAR/pair_gran_hooke_history.cpp
View File

@ -29,6 +29,7 @@
#include "fix.h"
#include "fix_pour.h"
#include "fix_shear_history.h"
#include "fix_rigid.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
@ -84,7 +85,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3;
double vtr1,vtr2,vtr3,vrel;
double meff,damp,ccel,tor1,tor2,tor3;
double mi,mj,meff,damp,ccel,tor1,tor2,tor3;
double fn,fs,fs1,fs2,fs3;
double shrmag,rsht;
int *ilist,*jlist,*numneigh,**firstneigh;
@ -98,6 +99,13 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
int shearupdate = 1;
if (update->setupflag) shearupdate = 0;
// update body ptr and values for ghost atoms if using FixRigid masses
if (fix_rigid && neighbor->ago == 0) {
body = fix_rigid->body;
comm->forward_comm_pair(this);
}
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
@ -181,19 +189,27 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
// normal forces = Hookian contact + normal velocity damping
// meff = effective mass of pair of particles
// if I or J part of rigid body, use body mass
// if I or J is frozen, meff is other particle
if (rmass) {
meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
if (mask[i] & freeze_group_bit) meff = rmass[j];
if (mask[j] & freeze_group_bit) meff = rmass[i];
mi = rmass[i];
mj = rmass[j];
} else {
itype = type[i];
jtype = type[j];
meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
if (mask[i] & freeze_group_bit) meff = mass[jtype];
if (mask[j] & freeze_group_bit) meff = mass[itype];
mi = mass[type[i]];
mj = mass[type[j]];
}
if (fix_rigid) {
if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]];
if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]];
}
meff = mi*mj / (mi+mj);
if (mask[i] & freeze_group_bit) meff = mj;
if (mask[j] & freeze_group_bit) meff = mi;
// normal forces = Hookian contact + normal velocity damping
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radsum-r)*rinv - damp;
@ -406,14 +422,14 @@ void PairGranHookeHistory::init_style()
fix_history->pair = this;
}
// check for Fix freeze and set freeze_group_bit
// check for FixFreeze and set freeze_group_bit
for (i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"freeze") == 0) break;
if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit;
else freeze_group_bit = 0;
// check for Fix pour and set pour_type and pour_maxdiam
// check for FixPour and set pour_type and pour_maxdiam
int pour_type = 0;
double pour_maxrad = 0.0;
@ -424,6 +440,13 @@ void PairGranHookeHistory::init_style()
pour_maxrad = ((FixPour *) modify->fix[i])->radius_hi;
}
// check for FixRigid
fix_rigid = NULL;
for (i = 0; i < modify->nfix; i++)
if (strstr(modify->fix[i]->style,"rigid")) break;
if (i < modify->nfix) fix_rigid = (FixRigid *) modify->fix[i];
// set maxrad_dynamic and maxrad_frozen for each type
// include future Fix pour particles as dynamic
@ -561,10 +584,10 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype,
double radi,radj,radsum;
double r,rinv,rsqinv,delx,dely,delz;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;
double meff,damp,ccel,polyhertz;
double mi,mj,meff,damp,ccel,polyhertz;
double vtr1,vtr2,vtr3,vrel,shrmag,rsht;
double fs1,fs2,fs3,fs,fn;
double *radius = atom->radius;
radi = radius[i];
radj = radius[j];
@ -612,21 +635,35 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype,
wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
// normal force = Hertzian contact + normal velocity damping
// meff = effective mass of pair of particles
// if I or J part of rigid body, use body mass
// if I or J is frozen, meff is other particle
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
if (rmass) {
meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
if (mask[i] & freeze_group_bit) meff = rmass[j];
if (mask[j] & freeze_group_bit) meff = rmass[i];
mi = rmass[i];
mj = rmass[j];
} else {
meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]);
if (mask[i] & freeze_group_bit) meff = mass[jtype];
if (mask[j] & freeze_group_bit) meff = mass[itype];
mi = mass[type[i]];
mj = mass[type[j]];
}
if (fix_rigid) {
// NOTE: need to make sure ghost atoms have updated body?
// depends on where single() is called from
body = fix_rigid->body;
if (body[i] >= 0) mi = fix_rigid->masstotal[body[i]];
if (body[j] >= 0) mj = fix_rigid->masstotal[body[j]];
}
meff = mi*mj / (mi+mj);
if (mask[i] & freeze_group_bit) meff = mj;
if (mask[j] & freeze_group_bit) meff = mi;
// normal forces = Hookian contact + normal velocity damping
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radsum-r)*rinv - damp;
@ -697,6 +734,33 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype,
/* ---------------------------------------------------------------------- */
int PairGranHookeHistory::pack_comm(int n, int *list,
double *buf, int pbc_flag, int *pbc)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = body[j];
}
return 1;
}
/* ---------------------------------------------------------------------- */
void PairGranHookeHistory::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++)
body[i] = static_cast<int> (buf[m++]);
}
/* ---------------------------------------------------------------------- */
void *PairGranHookeHistory::extract(const char *str, int &dim)
{
dim = 0;

4
src/GRANULAR/pair_gran_hooke_history.h
View File

@ -42,6 +42,8 @@ class PairGranHookeHistory : public Pair {
void read_restart_settings(FILE *);
void reset_dt();
double single(int, int, int, int, double, double, double, double &);
int pack_comm(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
void *extract(const char *, int &);
protected:
@ -57,6 +59,8 @@ class PairGranHookeHistory : public Pair {
double *maxrad_dynamic,*maxrad_frozen;
class FixShearHistory *fix_history;
class FixRigid *fix_rigid;
int *body;
void allocate();
};

6
src/fix_rigid.h
View File

@ -26,6 +26,10 @@ namespace LAMMPS_NS {
class FixRigid : public Fix {
public:
// public so that granular pair styles can access them
int *body; // which body each atom is part of (-1 if none)
double *masstotal; // total mass of each rigid body
FixRigid(class LAMMPS *, int, char **);
virtual ~FixRigid();
virtual int setmask();
@ -68,7 +72,6 @@ class FixRigid : public Fix {
int *mol2body; // convert mol-ID to rigid body index
int maxmol; // size of mol2body = max mol-ID
double *masstotal; // total mass of each rigid body
double **xcm; // coords of center-of-mass of each rigid body
double **vcm; // velocity of center-of-mass of each
double **fcm; // force on center-of-mass of each
@ -84,7 +87,6 @@ class FixRigid : public Fix {
double **tflag; // flag for on/off of center-of-mass torque
double **langextra; // Langevin thermostat forces and torques
int *body; // which body each atom is part of (-1 if none)
double **displace; // displacement of each atom in body coords
double **sum,**all; // work vectors for each rigid body