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

This commit is contained in:
sjplimp 2009-08-08 18:45:22 +00:00
parent 39c25bd61a
commit efea2a1835
3 changed files with 80 additions and 22 deletions

View File

@ -676,6 +676,52 @@ void Domain::remap(double *x, int &image)
if (triclinic) lamda2x(coord,x);
}
/* ----------------------------------------------------------------------
remap the point into the periodic box no matter how far away
resulting coord must satisfy lo <= coord < hi
MAX is important since coord - prd < lo can happen when coord = hi
for triclinic, point is converted to lamda coords (0-1) before doing remap
------------------------------------------------------------------------- */
void Domain::remap(double *x)
{
double *lo,*hi,*period,*coord;
double lamda[3];
if (triclinic == 0) {
lo = boxlo;
hi = boxhi;
period = prd;
coord = x;
} else {
lo = boxlo_lamda;
hi = boxhi_lamda;
period = prd_lamda;
x2lamda(x,lamda);
coord = lamda;
}
if (xperiodic) {
while (coord[0] < lo[0]) coord[0] += period[0];
while (coord[0] >= hi[0]) coord[0] -= period[0];
coord[0] = MAX(coord[0],lo[0]);
}
if (yperiodic) {
while (coord[1] < lo[1]) coord[1] += period[1];
while (coord[1] >= hi[1]) coord[1] -= period[1];
coord[1] = MAX(coord[1],lo[1]);
}
if (zperiodic) {
while (coord[2] < lo[2]) coord[2] += period[2];
while (coord[2] >= hi[2]) coord[2] -= period[2];
coord[2] = MAX(coord[2],lo[2]);
}
if (triclinic) lamda2x(coord,x);
}
/* ----------------------------------------------------------------------
unmap the point via image flags
x overwritten with result, don't reset image flag

View File

@ -90,6 +90,7 @@ class Domain : protected Pointers {
void reset_box();
void pbc();
void remap(double *, int &);
void remap(double *);
void unmap(double *, int);
void unmap(double *, int, double *);
void minimum_image(double &, double &, double &);

View File

@ -171,12 +171,15 @@ void FixIndent::post_force(int vflag)
if (istyle == SPHERE) {
// x1,y1,z1 = current position of indenter from original x0,y0,z0
// ctr = current indenter center from original x0,y0,z0
// remap into periodic box
double ctr[3];
double delta = (update->ntimestep - update->beginstep) * update->dt;
double x1 = x0 + delta*vx;
double y1 = y0 + delta*vy;
double z1 = z0 + delta*vz;
ctr[0] = x0 + delta*vx;
ctr[1] = y0 + delta*vy;
ctr[2] = z0 + delta*vz;
domain->remap(ctr);
double **x = atom->x;
double **f = atom->f;
@ -187,9 +190,10 @@ void FixIndent::post_force(int vflag)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
delx = x[i][0] - x1;
dely = x[i][1] - y1;
delz = x[i][2] - z1;
delx = x[i][0] - ctr[0];
dely = x[i][1] - ctr[1];
delz = x[i][2] - ctr[2];
domain->minimum_image(delx,dely,delz);
r = sqrt(delx*delx + dely*dely + delz*delz);
if (side == OUTSIDE) {
dr = r - r0;
@ -215,21 +219,27 @@ void FixIndent::post_force(int vflag)
} else if (istyle == CYLINDER) {
// c1new,c2new = current coords of indenter axis from original c1,c2
// ctr = current indenter axis from original c1,c2
// remap into periodic box
// 3rd coord is just near box for remap(), since isn't used
double ctr[3];
double delta = (update->ntimestep - update->beginstep) * update->dt;
double c1new,c2new;
if (cdim == 0) {
c1new = c1 + delta*vy;
c2new = c2 + delta*vz;
ctr[0] = domain->boxlo[0];
ctr[1] = c1 + delta*vy;
ctr[2] = c2 + delta*vz;
} else if (cdim == 1) {
c1new = c1 + delta*vx;
c2new = c2 + delta*vz;
ctr[0] = c1 + delta*vx;
ctr[1] = domain->boxlo[1];
ctr[2] = c2 + delta*vz;
} else {
c1new = c1 + delta*vx;
c2new = c2 + delta*vy;
ctr[0] = c1 + delta*vx;
ctr[1] = c2 + delta*vy;
ctr[2] = domain->boxlo[2];
}
domain->remap(ctr);
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
@ -241,17 +251,18 @@ void FixIndent::post_force(int vflag)
if (mask[i] & groupbit) {
if (cdim == 0) {
delx = 0;
dely = x[i][1] - c1new;
delz = x[i][2] - c2new;
dely = x[i][1] - ctr[1];
delz = x[i][2] - ctr[2];
} else if (cdim == 1) {
delx = x[i][0] - c1new;
delx = x[i][0] - ctr[0];
dely = 0;
delz = x[i][2] - c2new;
delz = x[i][2] - ctr[2];
} else {
delx = x[i][0] - c1new;
dely = x[i][1] - c2new;
delx = x[i][0] - ctr[0];
dely = x[i][1] - ctr[1];
delz = 0;
}
domain->minimum_image(delx,dely,delz);
r = sqrt(delx*delx + dely*dely + delz*delz);
if (side == OUTSIDE) {
dr = r - r0;