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

Browse Source
This commit is contained in:
sjplimp
2008-07-25 14:41:48 +00:00
parent 0afdb0b938
commit cc44a57275
10 changed files with 133 additions and 107 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/GRANULAR/atom_vec_granular.cpp
View File

@ -460,7 +460,7 @@ int AtomVecGranular::size_restart()
int i; int i;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int n = 19 * nlocal; int n = 16 * nlocal;
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++)

95
src/PERI/atom_vec_peri.cpp
View File

@ -37,13 +37,12 @@ AtomVec(lmp, narg, arg)
comm_x_only = 0; comm_x_only = 0;
size_comm = 4; size_comm = 4;
size_reverse = 3; size_reverse = 3;
size_border = 13; size_border = 11;
size_data_atom = 7; size_data_atom = 7;
size_data_vel = 4; size_data_vel = 4;
xcol_data = 3; xcol_data = 5;
atom->rmass_flag = atom->vfrac_flag = atom->s0_flag = 1; atom->vfrac_flag = atom->density_flag = atom->rmass_flag = 1;
atom->vinter_flag = 1;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -72,13 +71,13 @@ void AtomVecPeri::grow(int n)
vfrac = atom->vfrac = (double *) vfrac = atom->vfrac = (double *)
memory->srealloc(atom->vfrac,nmax*sizeof(double),"atom:vfrac"); memory->srealloc(atom->vfrac,nmax*sizeof(double),"atom:vfrac");
density = atom->density = (double *)
memory->srealloc(atom->density,nmax*sizeof(double),"atom:density");
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");
s0 = atom->s0 = (double *) s0 = atom->s0 = (double *)
memory->srealloc(atom->s0,nmax*sizeof(double),"atom:s0"); memory->srealloc(atom->s0,nmax*sizeof(double),"atom:s0");
x0 = atom->x0 = memory->grow_2d_double_array(atom->x0,nmax,3,"atom:x0"); x0 = atom->x0 = memory->grow_2d_double_array(atom->x0,nmax,3,"atom:x0");
vinter = atom->vinter = (double *)
memory->srealloc(atom->vinter,nmax*sizeof(double),"atom:vinter");
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++)
@ -101,12 +100,12 @@ void AtomVecPeri::copy(int i, int j)
v[j][2] = v[i][2]; v[j][2] = v[i][2];
vfrac[j] = vfrac[i]; vfrac[j] = vfrac[i];
density[j] = density[i];
rmass[j] = rmass[i]; rmass[j] = rmass[i];
s0[j] = s0[i]; s0[j] = s0[i];
x0[j][0] = x0[i][0]; x0[j][0] = x0[i][0];
x0[j][1] = x0[i][1]; x0[j][1] = x0[i][1];
x0[j][2] = x0[i][2]; x0[j][2] = x0[i][2];
vinter[j] = vinter[i];
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++)
@ -172,7 +171,7 @@ void AtomVecPeri::unpack_comm(int n, int first, double *buf)
x[i][0] = buf[m++]; x[i][0] = buf[m++];
x[i][1] = buf[m++]; x[i][1] = buf[m++];
x[i][2] = buf[m++]; x[i][2] = buf[m++];
s0[i] = buf[m++]; s0[i] = buf[m++];
} }
} }
@ -234,12 +233,10 @@ int AtomVecPeri::pack_border(int n, int *list, double *buf,
buf[m++] = type[j]; buf[m++] = type[j];
buf[m++] = mask[j]; buf[m++] = mask[j];
buf[m++] = vfrac[j]; buf[m++] = vfrac[j];
buf[m++] = rmass[j];
buf[m++] = s0[j]; buf[m++] = s0[j];
buf[m++] = x0[j][0]; buf[m++] = x0[j][0];
buf[m++] = x0[j][1]; buf[m++] = x0[j][1];
buf[m++] = x0[j][2]; buf[m++] = x0[j][2];
buf[m++] = vinter[j];
} }
} else { } else {
if (domain->triclinic == 0) { if (domain->triclinic == 0) {
@ -260,12 +257,10 @@ int AtomVecPeri::pack_border(int n, int *list, double *buf,
buf[m++] = type[j]; buf[m++] = type[j];
buf[m++] = mask[j]; buf[m++] = mask[j];
buf[m++] = vfrac[j]; buf[m++] = vfrac[j];
buf[m++] = rmass[j];
buf[m++] = s0[j]; buf[m++] = s0[j];
buf[m++] = x0[j][0]; buf[m++] = x0[j][0];
buf[m++] = x0[j][1]; buf[m++] = x0[j][1];
buf[m++] = x0[j][2]; buf[m++] = x0[j][2];
buf[m++] = vinter[j];
} }
} }
return m; return m;
@ -276,13 +271,11 @@ int AtomVecPeri::pack_border(int n, int *list, double *buf,
int AtomVecPeri::pack_border_one(int i, double *buf) int AtomVecPeri::pack_border_one(int i, double *buf)
{ {
buf[0] = vfrac[i]; buf[0] = vfrac[i];
buf[1] = rmass[i]; buf[1] = s0[i];
buf[2] = s0[i]; buf[2] = x0[i][0];
buf[3] = x0[i][0]; buf[3] = x0[i][1];
buf[4] = x0[i][1]; buf[4] = x0[i][2];
buf[5] = x0[i][2]; return 5;
buf[6] = vinter[i];
return 7;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -302,12 +295,10 @@ void AtomVecPeri::unpack_border(int n, int first, double *buf)
type[i] = static_cast<int> (buf[m++]); type[i] = static_cast<int> (buf[m++]);
mask[i] = static_cast<int> (buf[m++]); mask[i] = static_cast<int> (buf[m++]);
vfrac[i] = buf[m++]; vfrac[i] = buf[m++];
rmass[i] = buf[m++];
s0[i] = buf[m++]; s0[i] = buf[m++];
x0[i][0] = buf[m++]; x0[i][0] = buf[m++];
x0[i][1] = buf[m++]; x0[i][1] = buf[m++];
x0[i][2] = buf[m++]; x0[i][2] = buf[m++];
vinter[i] = buf[m++];
} }
} }
@ -316,13 +307,11 @@ void AtomVecPeri::unpack_border(int n, int first, double *buf)
int AtomVecPeri::unpack_border_one(int i, double *buf) int AtomVecPeri::unpack_border_one(int i, double *buf)
{ {
vfrac[i] = buf[0]; vfrac[i] = buf[0];
rmass[i] = buf[1]; s0[i] = buf[1];
s0[i] = buf[2]; x0[i][0] = buf[2];
x0[i][0] = buf[3]; x0[i][1] = buf[3];
x0[i][1] = buf[4]; x0[i][2] = buf[4];
x0[i][2] = buf[5]; return 5;
vinter[i]= buf[6];
return 7;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -345,12 +334,12 @@ int AtomVecPeri::pack_exchange(int i, double *buf)
buf[m++] = image[i]; buf[m++] = image[i];
buf[m++] = vfrac[i]; buf[m++] = vfrac[i];
buf[m++] = density[i];
buf[m++] = rmass[i]; buf[m++] = rmass[i];
buf[m++] = s0[i]; buf[m++] = s0[i];
buf[m++] = x0[i][0]; buf[m++] = x0[i][0];
buf[m++] = x0[i][1]; buf[m++] = x0[i][1];
buf[m++] = x0[i][2]; buf[m++] = x0[i][2];
buf[m++] = vinter[i];
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++)
@ -380,12 +369,12 @@ int AtomVecPeri::unpack_exchange(double *buf)
image[nlocal] = static_cast<int> (buf[m++]); image[nlocal] = static_cast<int> (buf[m++]);
vfrac[nlocal] = buf[m++]; vfrac[nlocal] = buf[m++];
density[nlocal] = buf[m++];
rmass[nlocal] = buf[m++]; rmass[nlocal] = buf[m++];
s0[nlocal] = buf[m++]; s0[nlocal] = buf[m++];
x0[nlocal][0] = buf[m++]; x0[nlocal][0] = buf[m++];
x0[nlocal][1] = buf[m++]; x0[nlocal][1] = buf[m++];
x0[nlocal][2] = buf[m++]; x0[nlocal][2] = buf[m++];
vinter[nlocal] = 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++)
@ -437,13 +426,13 @@ int AtomVecPeri::pack_restart(int i, double *buf)
buf[m++] = v[i][1]; buf[m++] = v[i][1];
buf[m++] = v[i][2]; buf[m++] = v[i][2];
buf[m++] = rmass[i];
buf[m++] = vfrac[i]; buf[m++] = vfrac[i];
buf[m++] = density[i];
buf[m++] = rmass[i];
buf[m++] = s0[i]; buf[m++] = s0[i];
buf[m++] = x0[i][0]; buf[m++] = x0[i][0];
buf[m++] = x0[i][1]; buf[m++] = x0[i][1];
buf[m++] = x0[i][2]; buf[m++] = x0[i][2];
buf[m++] = vinter[i];
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++)
@ -480,13 +469,13 @@ int AtomVecPeri::unpack_restart(double *buf)
v[nlocal][1] = buf[m++]; v[nlocal][1] = buf[m++];
v[nlocal][2] = buf[m++]; v[nlocal][2] = buf[m++];
rmass[nlocal] = buf[m++];
vfrac[nlocal] = buf[m++]; vfrac[nlocal] = buf[m++];
density[nlocal] = buf[m++];
rmass[nlocal] = buf[m++];
s0[nlocal] = buf[m++]; s0[nlocal] = buf[m++];
x0[nlocal][0] = buf[m++]; x0[nlocal][0] = buf[m++];
x0[nlocal][1] = buf[m++]; x0[nlocal][1] = buf[m++];
x0[nlocal][2] = buf[m++]; x0[nlocal][2] = buf[m++];
vinter[nlocal] = buf[m++];
double **extra = atom->extra; double **extra = atom->extra;
if (atom->nextra_store) { if (atom->nextra_store) {
@ -520,12 +509,12 @@ void AtomVecPeri::create_atom(int itype, double *coord)
v[nlocal][2] = 0.0; v[nlocal][2] = 0.0;
vfrac[nlocal] = 1.0; vfrac[nlocal] = 1.0;
rmass[nlocal] = 1.0; density[nlocal] = 1.0;
rmass[nlocal] = density[nlocal];
s0[nlocal] = MAXDOUBLE; s0[nlocal] = MAXDOUBLE;
x0[nlocal][0] = x[nlocal][0]; x0[nlocal][0] = coord[0];
x0[nlocal][1] = x[nlocal][1]; x0[nlocal][1] = coord[1];
x0[nlocal][2] = x[nlocal][2]; x0[nlocal][2] = coord[2];
vinter[nlocal] = 0.0;
atom->nlocal++; atom->nlocal++;
} }
@ -548,21 +537,25 @@ void AtomVecPeri::data_atom(double *coord, int imagetmp, char **values)
if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
error->one("Invalid atom type in Atoms section of data file"); error->one("Invalid atom type in Atoms section of data file");
vfrac[nlocal] = atof(values[5]); vfrac[nlocal] = atof(values[3]);
rmass[nlocal] = atof(values[6]); density[nlocal] = atof(values[4]);
s0[nlocal] = MAXDOUBLE; rmass[nlocal] = density[nlocal];
vinter[nlocal] = 0.0;
x[nlocal][0] = coord[0]; x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1]; x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2]; x[nlocal][2] = coord[2];
image[nlocal] = imagetmp; image[nlocal] = imagetmp;
mask[nlocal] = 1; mask[nlocal] = 1;
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;
s0[nlocal] = MAXDOUBLE;
x0[nlocal][0] = coord[0];
x0[nlocal][1] = coord[1];
x0[nlocal][2] = coord[2];
atom->nlocal++; atom->nlocal++;
} }
@ -576,13 +569,17 @@ void AtomVecPeri::data_atom(double *coord, int imagetmp, char **values)
int AtomVecPeri::data_atom_hybrid(int nlocal, char **values) int AtomVecPeri::data_atom_hybrid(int nlocal, char **values)
{ {
vfrac[nlocal] = atof(values[0]); vfrac[nlocal] = atof(values[0]);
rmass[nlocal] = atof(values[1]); density[nlocal] = atof(values[1]);
rmass[nlocal] = density[nlocal];
s0[nlocal] = MAXDOUBLE;
vinter[nlocal] = 0.0;
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;
s0[nlocal] = MAXDOUBLE;
x0[nlocal][0] = x[nlocal][0];
x0[nlocal][1] = x[nlocal][1];
x0[nlocal][2] = x[nlocal][2];
return 2; return 2;
} }
@ -604,10 +601,10 @@ double AtomVecPeri::memory_usage()
if (atom->memcheck("f")) bytes += nmax*3 * sizeof(double); if (atom->memcheck("f")) bytes += nmax*3 * sizeof(double);
if (atom->memcheck("vfrac")) bytes += nmax * sizeof(double); if (atom->memcheck("vfrac")) 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("s0")) bytes += nmax * sizeof(double); if (atom->memcheck("s0")) bytes += nmax * sizeof(double);
if (atom->memcheck("x0")) bytes += nmax*3 * sizeof(double); if (atom->memcheck("x0")) bytes += nmax*3 * sizeof(double);
if (atom->memcheck("vinter")) bytes += nmax * sizeof(double);
return bytes; return bytes;
} }

2
src/PERI/atom_vec_peri.h
View File

@ -46,7 +46,7 @@ class AtomVecPeri : public AtomVec {
private: private:
int *tag,*type,*mask,*image; int *tag,*type,*mask,*image;
double **x,**v,**f; double **x,**v,**f;
double *rmass,*vfrac,*s0,**x0,*vinter; double *vfrac,*density,*rmass,*s0,**x0;
}; };
} }

2
src/PERI/compute_damage_atom.cpp
View File

@ -87,8 +87,8 @@ void ComputeDamageAtom::compute_peratom()
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int *mask = atom->mask; int *mask = atom->mask;
double **x = atom->x; double **x = atom->x;
double *vinter = atom->vinter;
double *vfrac = atom->vfrac; double *vfrac = atom->vfrac;
double *vinter = ((FixPeriNeigh *) modify->fix[ifix_peri])->vinter;
int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
int i,j,jj,jnum; int i,j,jj,jnum;

29
src/PERI/fix_peri_neigh.cpp
View File

@ -50,6 +50,8 @@ FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) :
npartner = NULL; npartner = NULL;
partner = NULL; partner = NULL;
r0 = NULL; r0 = NULL;
vinter = NULL;
grow_arrays(atom->nmax); grow_arrays(atom->nmax);
atom->add_callback(0); atom->add_callback(0);
atom->add_callback(1); atom->add_callback(1);
@ -75,6 +77,7 @@ FixPeriNeigh::~FixPeriNeigh()
memory->sfree(npartner); memory->sfree(npartner);
memory->destroy_2d_int_array(partner); memory->destroy_2d_int_array(partner);
memory->destroy_2d_double_array(r0); memory->destroy_2d_double_array(r0);
memory->sfree(vinter);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -120,7 +123,6 @@ void FixPeriNeigh::setup(int vflag)
double **x = atom->x; double **x = atom->x;
double *vfrac = atom->vfrac; double *vfrac = atom->vfrac;
double *vinter = atom->vinter;
int *type = atom->type; int *type = atom->type;
int *tag = atom->tag; int *tag = atom->tag;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -144,8 +146,6 @@ void FixPeriNeigh::setup(int vflag)
Pair *anypair = force->pair_match("peri"); Pair *anypair = force->pair_match("peri");
double **cutsq = anypair->cutsq; double **cutsq = anypair->cutsq;
double *s0 = atom->s0;
for (ii = 0; ii < inum; ii++) { for (ii = 0; ii < inum; ii++) {
i = ilist[ii]; i = ilist[ii];
xtmp = x[i][0]; xtmp = x[i][0];
@ -173,9 +173,8 @@ void FixPeriNeigh::setup(int vflag)
int maxall; int maxall;
MPI_Allreduce(&maxpartner,&maxall,1,MPI_INT,MPI_MAX,world); MPI_Allreduce(&maxpartner,&maxall,1,MPI_INT,MPI_MAX,world);
// realloc arrays stored by fix with correct value for maxpartner // realloc arrays with correct value for maxpartner
memory->sfree(npartner);
memory->destroy_2d_int_array(partner); memory->destroy_2d_int_array(partner);
memory->destroy_2d_double_array(r0); memory->destroy_2d_double_array(r0);
npartner = NULL; npartner = NULL;
@ -184,8 +183,12 @@ void FixPeriNeigh::setup(int vflag)
grow_arrays(atom->nmax); grow_arrays(atom->nmax);
// create partner list and r0 values from neighbor list // create partner list and r0 values from neighbor list
// compute vinter for each atom
for (i = 0; i < nlocal; i++) npartner[i] = 0; for (i = 0; i < nlocal; i++) {
npartner[i] = 0;
vinter[i] = 0.0;
}
for (ii = 0; ii < inum; ii++) { for (ii = 0; ii < inum; ii++) {
i = ilist[ii]; i = ilist[ii];
@ -195,7 +198,6 @@ void FixPeriNeigh::setup(int vflag)
itype = type[i]; itype = type[i];
jlist = firstneigh[i]; jlist = firstneigh[i];
jnum = numneigh[i]; jnum = numneigh[i];
vinter[i] = 0.0;
for (jj = 0; jj < jnum; jj++) { for (jj = 0; jj < jnum; jj++) {
j = jlist[jj]; j = jlist[jj];
@ -246,6 +248,7 @@ double FixPeriNeigh::memory_usage()
int bytes = nmax * sizeof(int); int bytes = nmax * sizeof(int);
bytes += nmax*maxpartner * sizeof(int); bytes += nmax*maxpartner * sizeof(int);
bytes += nmax*maxpartner * sizeof(double); bytes += nmax*maxpartner * sizeof(double);
bytes += nmax * sizeof(double);
return bytes; return bytes;
} }
@ -260,6 +263,8 @@ void FixPeriNeigh::grow_arrays(int nmax)
partner = memory->grow_2d_int_array(partner,nmax,maxpartner, partner = memory->grow_2d_int_array(partner,nmax,maxpartner,
"peri_neigh:partner"); "peri_neigh:partner");
r0 = memory->grow_2d_double_array(r0,nmax,maxpartner,"peri_neigh:r0"); r0 = memory->grow_2d_double_array(r0,nmax,maxpartner,"peri_neigh:r0");
vinter = (double *) memory->srealloc(vinter,nmax*sizeof(double),
"peri_neigh:vinter");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -273,6 +278,7 @@ void FixPeriNeigh::copy_arrays(int i, int j)
partner[j][m] = partner[i][m]; partner[j][m] = partner[i][m];
r0[j][m] = r0[i][m]; r0[j][m] = r0[i][m];
} }
vinter[j] = vinter[i];
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -290,7 +296,9 @@ int FixPeriNeigh::pack_exchange(int i, double *buf)
buf[m++] = partner[i][n]; buf[m++] = partner[i][n];
buf[m++] = r0[i][n]; buf[m++] = r0[i][n];
} }
buf[0] = m/2; buf[0] = m/2;
buf[m++] = vinter[i];
return m; return m;
} }
@ -306,6 +314,7 @@ int FixPeriNeigh::unpack_exchange(int nlocal, double *buf)
partner[nlocal][n] = static_cast<int> (buf[m++]); partner[nlocal][n] = static_cast<int> (buf[m++]);
r0[nlocal][n] = buf[m++]; r0[nlocal][n] = buf[m++];
} }
vinter[nlocal] = buf[m++];
return m; return m;
} }
@ -322,6 +331,7 @@ int FixPeriNeigh::pack_restart(int i, double *buf)
buf[m++] = partner[i][n]; buf[m++] = partner[i][n];
buf[m++] = r0[i][n]; buf[m++] = r0[i][n];
} }
buf[m++] = vinter[i];
return m; return m;
} }
@ -344,6 +354,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth)
partner[nlocal][n] = static_cast<int> (extra[nlocal][m++]); partner[nlocal][n] = static_cast<int> (extra[nlocal][m++]);
r0[nlocal][n] = extra[nlocal][m++]; r0[nlocal][n] = extra[nlocal][m++];
} }
vinter[nlocal] = extra[nlocal][m++];
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -352,7 +363,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth)
int FixPeriNeigh::maxsize_restart() int FixPeriNeigh::maxsize_restart()
{ {
return 2*maxpartner + 2; return 2*maxpartner + 3;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -361,5 +372,5 @@ int FixPeriNeigh::maxsize_restart()
int FixPeriNeigh::size_restart(int nlocal) int FixPeriNeigh::size_restart(int nlocal)
{ {
return 2*npartner[nlocal] + 2; return 2*npartner[nlocal] + 3;
} }

3
src/PERI/fix_peri_neigh.h
View File

@ -41,11 +41,12 @@ class FixPeriNeigh : public Fix {
int maxsize_restart(); int maxsize_restart();
private: private:
int first; int first; // flag for first time initialization
int maxpartner; // max # of peridynamic neighs for any atom int maxpartner; // max # of peridynamic neighs for any atom
int *npartner; // # of neighbors for each atom int *npartner; // # of neighbors for each atom
int **partner; // neighs for each atom, stored as global IDs int **partner; // neighs for each atom, stored as global IDs
double **r0; // initial distance to partners double **r0; // initial distance to partners
double *vinter; // sum of vfrac for bonded neighbors
class NeighList *list; class NeighList *list;
}; };

57
src/PERI/pair_peri_pmb.cpp
View File

@ -96,7 +96,6 @@ void PairPeriPMB::compute(int eflag, int vflag)
double *vfrac = atom->vfrac; double *vfrac = atom->vfrac;
double *s0 = atom->s0; double *s0 = atom->s0;
double **x0 = atom->x0; double **x0 = atom->x0;
double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
@ -112,13 +111,15 @@ void PairPeriPMB::compute(int eflag, int vflag)
int nall = atom->nlocal + atom->nghost; int nall = atom->nlocal + atom->nghost;
int newton_pair = force->newton_pair; int newton_pair = force->newton_pair;
int nonperiodic = domain->nonperiodic;
inum = list->inum; inum = list->inum;
ilist = list->ilist; ilist = list->ilist;
numneigh = list->numneigh; numneigh = list->numneigh;
firstneigh = list->firstneigh; firstneigh = list->firstneigh;
// loop over neighbors of my atoms // loop over neighbors of my atoms
// need minimg() for x0 difference since not ghosted
for (ii = 0; ii < inum; ii++) { for (ii = 0; ii < inum; ii++) {
i = ilist[ii]; i = ilist[ii];
@ -139,24 +140,23 @@ void PairPeriPMB::compute(int eflag, int vflag)
delx = xtmp - x[j][0]; delx = xtmp - x[j][0];
dely = ytmp - x[j][1]; dely = ytmp - x[j][1];
delz = ztmp - x[j][2]; delz = ztmp - x[j][2];
domain->minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz; rsq = delx*delx + dely*dely + delz*delz;
delx0 = xtmp0 - x0[j][0]; delx0 = xtmp0 - x0[j][0];
dely0 = ytmp0 - x0[j][1]; dely0 = ytmp0 - x0[j][1];
delz0 = ztmp0 - x0[j][2]; delz0 = ztmp0 - x0[j][2];
domain->minimum_image(delx0,dely0,delz0); if (nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0);
rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
jtype = type[j]; jtype = type[j];
r = sqrt(rsq); r = sqrt(rsq);
// shortrange interaction distance based on initial particle position // short-range interaction distance based on initial particle position
// 0.9 and 1.35 are constants // 0.9 and 1.35 are constants
d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); d_ij = MIN(0.9*sqrt(rsq0),1.35*lc);
// apply short-range contact forces // short-range contact forces
// 15 is a constant taken from the EMU Theory Manual // 15 is constant taken from the EMU Theory Manual
// Silling, 12 May 2005, p 18 // Silling, 12 May 2005, p 18
if (r < d_ij) { if (r < d_ij) {
@ -191,13 +191,13 @@ void PairPeriPMB::compute(int eflag, int vflag)
s0_new = (double *) memory->smalloc(nmax*sizeof(double),"pair:s0_new"); s0_new = (double *) memory->smalloc(nmax*sizeof(double),"pair:s0_new");
} }
// first = flag indicating if this is first neighbor of particle i // loop over my particles and their partners
// partner list contains all bond partners, so I-J appears twice
// if bond already broken, skip this partner
// first = true if this is first neighbor of particle i
bool first; bool first;
// loop over my particles and their partners
// if bond already broken, skip this partner
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
xtmp = x[i][0]; xtmp = x[i][0];
ytmp = x[i][1]; ytmp = x[i][1];
@ -211,22 +211,22 @@ void PairPeriPMB::compute(int eflag, int vflag)
if (partner[i][jj] == 0) continue; if (partner[i][jj] == 0) continue;
j = atom->map(partner[i][jj]); j = atom->map(partner[i][jj]);
// check if we've lost a partner without first breaking bond // check if lost a partner without first breaking bond
if (j < 0) { if (j < 0) {
partner[i][jj] = 0; partner[i][jj] = 0;
continue; continue;
} }
// compute force density and add to PD equation of motion // compute force density, add to PD equation of motion
delx = xtmp - x[j][0]; delx = xtmp - x[j][0];
dely = ytmp - x[j][1]; dely = ytmp - x[j][1];
delz = ztmp - x[j][2]; delz = ztmp - x[j][2];
domain->minimum_image(delx,dely,delz); if (nonperiodic == 0) domain->minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz; rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j]; jtype = type[j];
delta = sqrt(cutsq[itype][jtype]); // the horizon delta = sqrt(cutsq[itype][jtype]);
r = sqrt(rsq); r = sqrt(rsq);
dr = r - r0[i][jj]; dr = r - r0[i][jj];
@ -234,7 +234,7 @@ void PairPeriPMB::compute(int eflag, int vflag)
if (fabs(dr) < 2.2204e-016) dr = 0.0; if (fabs(dr) < 2.2204e-016) dr = 0.0;
// apply scaling for vfrac[j] if particle j near the horizon // scale vfrac[j] if particle j near the horizon
if ((fabs(r0[i][jj] - delta)) <= half_lc) if ((fabs(r0[i][jj] - delta)) <= half_lc)
vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) +
@ -250,11 +250,13 @@ void PairPeriPMB::compute(int eflag, int vflag)
f[i][1] += dely*fbond; f[i][1] += dely*fbond;
f[i][2] += delz*fbond; f[i][2] += delz*fbond;
if (eflag) evdwl = rk*dr; // hardwire newton flag off and use 1/2 since I-J is double counted
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fbond,delx,dely,delz);
// find stretch in bond i-j and break if necessary if (eflag) evdwl = rk*dr;
//if (evflag) ev_tally(i,j,nlocal,0,
// 0.5*evdwl,0.0,0.5*fbond,delx,dely,delz);
// find stretch in bond I-J and break if necessary
// use s0 from previous timestep // use s0 from previous timestep
if (stretch > MIN(s0[i],s0[j])) partner[i][jj] = 0; if (stretch > MIN(s0[i],s0[j])) partner[i][jj] = 0;
@ -266,14 +268,13 @@ void PairPeriPMB::compute(int eflag, int vflag)
else else
s0_new[i] = MAX(s0_new[i], s0_new[i] = MAX(s0_new[i],
s00[itype][jtype] - (alpha[itype][jtype] * stretch)); s00[itype][jtype] - (alpha[itype][jtype] * stretch));
// first neighbor of particle i has now been examined
first = false; first = false;
} }
} }
// update with newly computed s0 if (vflag_fdotr) virial_compute();
// store new s0
for (i = 0; i < nlocal; i++) s0[i] = s0_new[i]; for (i = 0; i < nlocal; i++) s0[i] = s0_new[i];
} }
@ -353,6 +354,7 @@ double PairPeriPMB::init_one(int i, int j)
cutsq[j][i] = cutsq[i][j]; cutsq[j][i] = cutsq[i][j];
// set other j,i parameters // set other j,i parameters
kspring[j][i] = kspring[i][j]; kspring[j][i] = kspring[i][j];
alpha[j][i] = alpha[i][j]; alpha[j][i] = alpha[i][j];
s00[j][i] = s00[i][j]; s00[j][i] = s00[i][j];
@ -460,18 +462,17 @@ double PairPeriPMB::single(int i, int j, int itype, int jtype, double rsq,
double *vfrac = atom->vfrac; double *vfrac = atom->vfrac;
double **x0 = atom->x0; double **x0 = atom->x0;
double lc = domain->lattice->xlattice;
double half_lc = 0.5*lc;
double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
double lc = domain->lattice->xlattice;
double half_lc = 0.5*lc;
delx0 = x0[i][0] - x0[j][0]; delx0 = x0[i][0] - x0[j][0];
dely0 = x0[i][1] - x0[j][1]; dely0 = x0[i][1] - x0[j][1];
delz0 = x0[i][2] - x0[j][2]; delz0 = x0[i][2] - x0[j][2];
domain->minimum_image(delx0,dely0,delz0); if (domain->nonperiodic == 0) domain->minimum_image(delx0,dely0,delz0);
rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0;
d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); d_ij = MIN(0.9*sqrt(rsq0),1.35*lc);

8
src/atom.cpp
View File

@ -64,7 +64,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
q = NULL; q = NULL;
mu = NULL; mu = NULL;
quat = omega = angmom = torque = NULL; quat = omega = angmom = torque = NULL;
radius = density = rmass = vfrac = s0 = vinter = NULL; radius = density = rmass = NULL;
vfrac = s0 = NULL;
x0 = NULL; x0 = NULL;
maxspecial = 1; maxspecial = 1;
@ -91,7 +92,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
molecule_flag = 0; molecule_flag = 0;
q_flag = mu_flag = 0; q_flag = mu_flag = 0;
quat_flag = omega_flag = angmom_flag = torque_flag = 0; quat_flag = omega_flag = angmom_flag = torque_flag = 0;
radius_flag = density_flag = rmass_flag = vfrac_flag = s0_flag = vinter_flag = 0; radius_flag = density_flag = rmass_flag = vfrac_flag = 0;
// ntype-length arrays // ntype-length arrays
@ -164,7 +165,6 @@ Atom::~Atom()
memory->sfree(rmass); memory->sfree(rmass);
memory->sfree(vfrac); memory->sfree(vfrac);
memory->sfree(s0); memory->sfree(s0);
memory->sfree(vinter);
memory->destroy_2d_double_array(x0); memory->destroy_2d_double_array(x0);
memory->sfree(molecule); memory->sfree(molecule);
@ -293,7 +293,7 @@ void Atom::init()
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
return 1 if style matches atom style hybrid sub-style return 1 if style matches atom style or hybrid sub-style
else return 0 else return 0
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

5
src/atom.h
View File

@ -48,7 +48,8 @@ class Atom : protected Pointers {
int *molecule; int *molecule;
double *q,**mu; double *q,**mu;
double **quat,**omega,**angmom,**torque; double **quat,**omega,**angmom,**torque;
double *radius,*density,*rmass,*vfrac,*s0,**x0,*vinter; double *radius,*density,*rmass,*vfrac,*s0;
double **x0;
int maxspecial; int maxspecial;
int **nspecial,**special; int **nspecial,**special;
@ -75,7 +76,7 @@ class Atom : protected Pointers {
int molecule_flag; int molecule_flag;
int q_flag,mu_flag; int q_flag,mu_flag;
int quat_flag,omega_flag,angmom_flag,torque_flag; int quat_flag,omega_flag,angmom_flag,torque_flag;
int radius_flag,density_flag,rmass_flag,vfrac_flag,s0_flag,vinter_flag; int radius_flag,density_flag,rmass_flag,vfrac_flag;
// extra peratom info in restart file destined for fix & diag // extra peratom info in restart file destined for fix & diag

37
src/set.cpp
View File

@ -185,10 +185,8 @@ void Set::command(int narg, char **arg)
} else if (strcmp(arg[iarg],"density") == 0) { } else if (strcmp(arg[iarg],"density") == 0) {
if (iarg+2 > narg) error->all("Illegal set command"); if (iarg+2 > narg) error->all("Illegal set command");
dvalue = atof(arg[iarg+1]); dvalue = atof(arg[iarg+1]);
/*
if (!atom->density_flag) if (!atom->density_flag)
error->all("Cannot set this attribute for this atom style"); error->all("Cannot set this attribute for this atom style");
*/
set(DENSITY); set(DENSITY);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"volume") == 0) { } else if (strcmp(arg[iarg],"volume") == 0) {
@ -302,6 +300,13 @@ void Set::selection(int n)
void Set::set(int keyword) void Set::set(int keyword)
{ {
// set granular and peridynamic flags
int granflag = 0;
int periflag = 0;
if (atom->style_match("granular")) granflag = 1;
if (atom->style_match("peri")) periflag = 1;
if (keyword == DIPOLE) atom->check_dipole(); if (keyword == DIPOLE) atom->check_dipole();
selection(atom->nlocal); selection(atom->nlocal);
@ -318,6 +323,9 @@ void Set::set(int keyword)
else if (keyword == VY) atom->v[i][1] = dvalue; else if (keyword == VY) atom->v[i][1] = dvalue;
else if (keyword == VZ) atom->v[i][2] = dvalue; else if (keyword == VZ) atom->v[i][2] = dvalue;
else if (keyword == CHARGE) atom->q[i] = dvalue; else if (keyword == CHARGE) atom->q[i] = dvalue;
// diameter setting triggers setting of rmass
else if (keyword == DIAMETER) { else if (keyword == DIAMETER) {
atom->radius[i] = 0.5 * dvalue; atom->radius[i] = 0.5 * dvalue;
if (domain->dimension == 3) if (domain->dimension == 3)
@ -326,17 +334,24 @@ void Set::set(int keyword)
else else
atom->rmass[i] = PI * atom->rmass[i] = PI *
atom->radius[i]*atom->radius[i] * atom->density[i]; atom->radius[i]*atom->radius[i] * atom->density[i];
// density setting triggers setting of rmass
// for granular, rmass is function of diameter and density
// for peri, rmass stores density directly
} else if (keyword == DENSITY) { } else if (keyword == DENSITY) {
atom->rmass[i] = dvalue;
/*
atom->density[i] = dvalue; atom->density[i] = dvalue;
if (domain->dimension == 3) if (granflag) {
atom->rmass[i] = 4.0*PI/3.0 * if (domain->dimension == 3)
atom->radius[i]*atom->radius[i]*atom->radius[i] * atom->density[i]; atom->rmass[i] = 4.0*PI/3.0 *
else atom->radius[i]*atom->radius[i]*atom->radius[i] *
atom->rmass[i] = PI * atom->density[i];
atom->radius[i]*atom->radius[i] * atom->density[i]; else
*/ atom->rmass[i] = PI *
atom->radius[i]*atom->radius[i] * atom->density[i];
}
if (periflag) atom->rmass[i] = dvalue;
} else if (keyword == VOLUME) atom->vfrac[i] = dvalue; } else if (keyword == VOLUME) atom->vfrac[i] = dvalue;
else if (keyword == DIPOLE) { else if (keyword == DIPOLE) {