ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2010-09-10 17:08:44 +00:00
parent edf05df900
commit 66178b66ab
4 changed files with 235 additions and 74 deletions
Download Patch File Download Diff File Expand all files Collapse all files

152
src/GRANULAR/atom_vec_granular.cpp
View File

@ -13,12 +13,14 @@
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "atom_vec_granular.h"
#include "atom.h"
#include "domain.h"
#include "modify.h"
#include "force.h"
#include "fix.h"
#include "fix_adapt.h"
#include "memory.h"
#include "error.h"
@ -46,9 +48,27 @@ AtomVecGranular::AtomVecGranular(LAMMPS *lmp, int narg, char **arg) :
atom->radius_flag = atom->density_flag = atom->rmass_flag = 1;
atom->omega_flag = atom->torque_flag = 1;
radvary = 0;
PI = 4.0*atan(1.0);
}
/* ---------------------------------------------------------------------- */
void AtomVecGranular::init()
{
// set radvary if particle diameters are time-varying due to fix adapt
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"adapt") == 0) {
FixAdapt *fix = (FixAdapt *) modify->fix[i];
if (fix->diamflag) {
radvary = 1;
size_forward = 5;
}
}
}
/* ----------------------------------------------------------------------
grow atom arrays
n = 0 grows arrays by DELTA
@ -138,6 +158,7 @@ int AtomVecGranular::pack_comm(int n, int *list, double *buf,
int i,j,m;
double dx,dy,dz;
if (radvary == 0) {
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
@ -163,6 +184,39 @@ int AtomVecGranular::pack_comm(int n, int *list, double *buf,
buf[m++] = x[j][2] + dz;
}
}
} else {
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
buf[m++] = radius[i];
buf[m++] = rmass[i];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
dz = pbc[2]*domain->zprd;
}
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = radius[i];
buf[m++] = rmass[i];
}
}
}
return m;
}
@ -174,6 +228,7 @@ int AtomVecGranular::pack_comm_vel(int n, int *list, double *buf,
int i,j,m;
double dx,dy,dz;
if (radvary == 0) {
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
@ -211,15 +266,72 @@ int AtomVecGranular::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = omega[j][2];
}
}
} else {
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
buf[m++] = radius[i];
buf[m++] = rmass[i];
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
buf[m++] = omega[j][0];
buf[m++] = omega[j][1];
buf[m++] = omega[j][2];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
dz = pbc[2]*domain->zprd;
}
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = radius[i];
buf[m++] = rmass[i];
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
buf[m++] = omega[j][0];
buf[m++] = omega[j][1];
buf[m++] = omega[j][2];
}
}
}
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecGranular::pack_comm_one(int i, double *buf)
{
if (radvary == 0) return 0;
buf[0] = radius[i];
buf[1] = rmass[i];
return 2;
}
/* ---------------------------------------------------------------------- */
void AtomVecGranular::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
if (radvary == 0) {
m = 0;
last = first + n;
for (i = first; i < last; i++) {
@ -227,6 +339,17 @@ void AtomVecGranular::unpack_comm(int n, int first, double *buf)
x[i][1] = buf[m++];
x[i][2] = buf[m++];
}
} else {
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
radius[i] = buf[m++];
rmass[i] = buf[m++];
}
}
}
/* ---------------------------------------------------------------------- */
@ -235,6 +358,7 @@ void AtomVecGranular::unpack_comm_vel(int n, int first, double *buf)
{
int i,m,last;
if (radvary == 0) {
m = 0;
last = first + n;
for (i = first; i < last; i++) {
@ -248,6 +372,34 @@ void AtomVecGranular::unpack_comm_vel(int n, int first, double *buf)
omega[i][1] = buf[m++];
omega[i][2] = buf[m++];
}
} else {
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
radius[i] = buf[m++];
rmass[i] = buf[m++];
v[i][0] = buf[m++];
v[i][1] = buf[m++];
v[i][2] = buf[m++];
omega[i][0] = buf[m++];
omega[i][1] = buf[m++];
omega[i][2] = buf[m++];
}
}
}
/* ---------------------------------------------------------------------- */
int AtomVecGranular::unpack_comm_one(int i, double *buf)
{
if (radvary == 0) return 0;
radius[i] = buf[0];
rmass[i] = buf[1];
return 2;
}
/* ---------------------------------------------------------------------- */

4
src/GRANULAR/atom_vec_granular.h
View File

@ -28,13 +28,16 @@ class AtomVecGranular : public AtomVec {
public:
AtomVecGranular(class LAMMPS *, int, char **);
~AtomVecGranular() {}
void init();
void grow(int);
void grow_reset();
void copy(int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
int pack_comm_one(int, double *);
void unpack_comm(int, int, double *);
void unpack_comm_vel(int, int, double *);
int unpack_comm_one(int, double *);
int pack_reverse(int, int, double *);
int pack_reverse_one(int, double *);
void unpack_reverse(int, int *, double *);
@ -63,6 +66,7 @@ class AtomVecGranular : public AtomVec {
double **x,**v,**f;
double *radius,*density,*rmass;
double **omega,**torque;
int radvary;
};
}

7
src/fix_adapt.cpp
View File

@ -71,6 +71,7 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
// parse keywords
diamflag = 0;
nadapt = 0;
iarg = 4;
@ -99,6 +100,7 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
int n = strlen(arg[iarg+1]) + 1;
param[nadapt] = new char[n];
strcpy(param[nadapt],arg[iarg+1]);
if (strcmp(param[nadapt],"diameter") == 0) diamflag = 1;
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
n = strlen(&arg[iarg+2][2]) + 1;
var[nadapt] = new char[n];
@ -201,7 +203,7 @@ void FixAdapt::pre_force(int vflag)
else if (which[m] == ATOM) {
// set radius from diameter
// set rmass if both rmass and density are defined
// also set rmass if both rmass and density are defined
if (awhich[m] == DIAMETER) {
int mflag = 0;
@ -217,7 +219,8 @@ void FixAdapt::pre_force(int vflag)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
radius[i] = 0.5*value;
rmass[i] = 4.0*PI/3.0 * radius[i]*radius[i]*radius[i] * density[i];
if (mflag) rmass[i] = 4.0*PI/3.0 *
radius[i]*radius[i]*radius[i] * density[i];
}
}
}

2
src/fix_adapt.h
View File

@ -26,6 +26,8 @@ namespace LAMMPS_NS {
class FixAdapt : public Fix {
public:
int diamflag;
FixAdapt(class LAMMPS *, int, char **);
~FixAdapt();
int setmask();