From af5d91115f4bf90d1c557277b95ce2f9590472ac Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 2 May 2014 14:28:56 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11875 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_reduce_region.cpp | 4 +- src/compute_temp_region.cpp | 4 ++ src/create_atoms.cpp | 2 + src/delete_atoms.cpp | 3 ++ src/dump_custom.cpp | 1 + src/fix_addforce.cpp | 18 +++++---- src/fix_ave_spatial.cpp | 4 ++ src/fix_aveforce.cpp | 34 ++++++++-------- src/fix_heat.cpp | 13 +++--- src/fix_setforce.cpp | 26 ++++++++---- src/fix_wall_region.cpp | 2 + src/group.cpp | 15 ++++++- src/region.cpp | 75 ++++++++++++++++------------------- src/region.h | 3 +- src/set.cpp | 1 + 15 files changed, 125 insertions(+), 80 deletions(-) diff --git a/src/compute_reduce_region.cpp b/src/compute_reduce_region.cpp index 3f54dd5bbe..fe33fc3c79 100644 --- a/src/compute_reduce_region.cpp +++ b/src/compute_reduce_region.cpp @@ -58,6 +58,7 @@ double ComputeReduceRegion::compute_one(int m, int flag) int i; Region *region = domain->regions[iregion]; + region->prematch(); // invoke the appropriate attribute,compute,fix,variable // compute scalar quantity by summing over atom scalars @@ -156,7 +157,8 @@ double ComputeReduceRegion::compute_one(int m, int flag) } else if (which[m] == FIX) { if (update->ntimestep % modify->fix[n]->peratom_freq) - error->all(FLERR,"Fix used in compute reduce not computed at compatible time"); + error->all(FLERR,"Fix used in compute reduce not computed at " + "compatible time"); Fix *fix = modify->fix[n]; if (flavor[m] == PERATOM) { diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp index c89c0a10ae..cc7f79111b 100644 --- a/src/compute_temp_region.cpp +++ b/src/compute_temp_region.cpp @@ -102,6 +102,8 @@ double ComputeTempRegion::compute_scalar() int nlocal = atom->nlocal; Region *region = domain->regions[iregion]; + region->prematch(); + int count = 0; double t = 0.0; @@ -147,6 +149,8 @@ void ComputeTempRegion::compute_vector() int nlocal = atom->nlocal; Region *region = domain->regions[iregion]; + region->prematch(); + double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 351d3656b1..a9393729b7 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -72,6 +72,7 @@ void CreateAtoms::command(int narg, char **arg) if (nregion == -1) error->all(FLERR, "Create_atoms region ID does not exist"); domain->regions[nregion]->init(); + domain->regions[nregion]->prematch(); iarg = 3;; } else if (strcmp(arg[1],"single") == 0) { style = SINGLE; @@ -91,6 +92,7 @@ void CreateAtoms::command(int narg, char **arg) if (nregion == -1) error->all(FLERR, "Create_atoms region ID does not exist"); domain->regions[nregion]->init(); + domain->regions[nregion]->prematch(); } iarg = 5; } else error->all(FLERR,"Illegal create_atoms command"); diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp index e75e2e5abe..916fcadada 100644 --- a/src/delete_atoms.cpp +++ b/src/delete_atoms.cpp @@ -153,6 +153,8 @@ void DeleteAtoms::delete_region(int narg, char **arg) int iregion = domain->find_region(arg[1]); if (iregion == -1) error->all(FLERR,"Could not find delete_atoms region ID"); + domain->regions[iregion]->prematch(); + options(narg-2,&arg[2]); // allocate and initialize deletion list @@ -332,6 +334,7 @@ void DeleteAtoms::delete_porosity(int narg, char **arg) int iregion = domain->find_region(arg[1]); if (iregion == -1) error->all(FLERR,"Could not find delete_atoms region ID"); + domain->regions[iregion]->prematch(); double porosity_fraction = force->numeric(FLERR,arg[2]); int seed = force->inumeric(FLERR,arg[3]); diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 1c7bf6fcc6..bafe82882a 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -419,6 +419,7 @@ int DumpCustom::count() if (iregion >= 0) { Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; for (i = 0; i < nlocal; i++) if (choose[i] && region->match(x[i][0],x[i][1],x[i][2]) == 0) diff --git a/src/fix_addforce.cpp b/src/fix_addforce.cpp index dcac27a012..fe15268376 100644 --- a/src/fix_addforce.cpp +++ b/src/fix_addforce.cpp @@ -224,6 +224,14 @@ void FixAddForce::post_force(int vflag) imageint *image = atom->image; int nlocal = atom->nlocal; + // update region if necessary + + Region *region = NULL; + if (iregion >= 0) { + region = domain->regions[iregion]; + region->prematch(); + } + // reallocate sforce array if necessary if ((varflag == ATOM || estyle == ATOM) && nlocal > maxatom) { @@ -245,10 +253,7 @@ void FixAddForce::post_force(int vflag) double unwrap[3]; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (iregion >= 0 && - !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) - continue; - + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; domain->unmap(x[i],image[i],unwrap); foriginal[0] -= xvalue*unwrap[0] + yvalue*unwrap[1] + zvalue*unwrap[2]; foriginal[1] += f[i][0]; @@ -291,10 +296,7 @@ void FixAddForce::post_force(int vflag) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (iregion >= 0 && - !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) - continue; - + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; if (estyle == ATOM) foriginal[0] += sforce[i][3]; foriginal[1] += f[i][0]; foriginal[2] += f[i][1]; diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp index c9defd4078..c55ffd0e7f 100644 --- a/src/fix_ave_spatial.cpp +++ b/src/fix_ave_spatial.cpp @@ -508,6 +508,10 @@ void FixAveSpatial::end_of_step() bigint ntimestep = update->ntimestep; if (ntimestep != nvalid) return; + // update region if necessary + + if (regionflag) region->prematch(); + // zero out arrays that accumulate over many samples // if box changes, first re-setup bins diff --git a/src/fix_aveforce.cpp b/src/fix_aveforce.cpp index 6b4bf40bc8..7cd44412eb 100644 --- a/src/fix_aveforce.cpp +++ b/src/fix_aveforce.cpp @@ -189,6 +189,14 @@ void FixAveForce::min_setup(int vflag) void FixAveForce::post_force(int vflag) { + // update region if necessary + + Region *region = NULL; + if (iregion >= 0) { + region = domain->regions[iregion]; + region->prematch(); + } + // sum forces on participating atoms double **x = atom->x; @@ -201,10 +209,7 @@ void FixAveForce::post_force(int vflag) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (iregion >= 0 && - !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) - continue; - + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; foriginal[0] += f[i][0]; foriginal[1] += f[i][1]; foriginal[2] += f[i][2]; @@ -238,10 +243,7 @@ void FixAveForce::post_force(int vflag) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (iregion >= 0 && - !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) - continue; - + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; if (xstyle) f[i][0] = fave[0]; if (ystyle) f[i][1] = fave[1]; if (zstyle) f[i][2] = fave[2]; @@ -257,6 +259,12 @@ void FixAveForce::post_force_respa(int vflag, int ilevel, int iloop) if (ilevel == nlevels_respa-1) post_force(vflag); else { + Region *region = NULL; + if (iregion >= 0) { + region = domain->regions[iregion]; + region->prematch(); + } + double **x = atom->x; double **f = atom->f; int *mask = atom->mask; @@ -267,10 +275,7 @@ void FixAveForce::post_force_respa(int vflag, int ilevel, int iloop) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (iregion >= 0 && - !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) - continue; - + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; foriginal[0] += f[i][0]; foriginal[1] += f[i][1]; foriginal[2] += f[i][2]; @@ -289,10 +294,7 @@ void FixAveForce::post_force_respa(int vflag, int ilevel, int iloop) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (iregion >= 0 && - !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) - continue; - + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; if (xstyle) f[i][0] = fave[0]; if (ystyle) f[i][1] = fave[1]; if (zstyle) f[i][2] = fave[2]; diff --git a/src/fix_heat.cpp b/src/fix_heat.cpp index 82d16f235d..de23927fe9 100644 --- a/src/fix_heat.cpp +++ b/src/fix_heat.cpp @@ -142,7 +142,6 @@ void FixHeat::end_of_step() int i; double heat,ke,massone; double vsub[3],vcm[3]; - Region *region; double **x = atom->x; double **v = atom->v; @@ -189,6 +188,12 @@ void FixHeat::end_of_step() // vsub = velocity subtracted from each atom to preserve momentum // overall KE cannot go negative + Region *region = NULL; + if (iregion >= 0) { + region = domain->regions[iregion]; + region->prematch(); + } + if (hstyle != ATOM) { heat = heat_input*nevery*update->dt*force->ftm2v; double escale = @@ -207,7 +212,6 @@ void FixHeat::end_of_step() v[i][2] = scale*v[i][2] - vsub[2]; } } else { - region = domain->regions[iregion]; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { v[i][0] = scale*v[i][0] - vsub[0]; @@ -254,7 +258,6 @@ void FixHeat::end_of_step() } } else { - region = domain->regions[iregion]; for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { heat = vheat[i]*nevery*update->dt*force->ftm2v; @@ -306,8 +309,8 @@ double FixHeat::compute_scalar() } } } else { - Region *region; - region = domain->regions[iregion]; + Region *region = domain->regions[iregion]; + region->prematch(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { scale_sum += sqrt(vscale[i]); diff --git a/src/fix_setforce.cpp b/src/fix_setforce.cpp index 848c4a003b..fd21ada90f 100644 --- a/src/fix_setforce.cpp +++ b/src/fix_setforce.cpp @@ -218,6 +218,14 @@ void FixSetForce::post_force(int vflag) int *mask = atom->mask; int nlocal = atom->nlocal; + // update region if necessary + + Region *region = NULL; + if (iregion >= 0) { + region = domain->regions[iregion]; + region->prematch(); + } + // reallocate sforce array if necessary if (varflag == ATOM && nlocal > maxatom) { @@ -232,10 +240,7 @@ void FixSetForce::post_force(int vflag) if (varflag == CONSTANT) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (iregion >= 0 && - !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) - continue; - + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; foriginal[0] += f[i][0]; foriginal[1] += f[i][1]; foriginal[2] += f[i][2]; @@ -270,10 +275,7 @@ void FixSetForce::post_force(int vflag) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (iregion >= 0 && - !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) - continue; - + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; foriginal[0] += f[i][0]; foriginal[1] += f[i][1]; foriginal[2] += f[i][2]; @@ -295,12 +297,20 @@ void FixSetForce::post_force_respa(int vflag, int ilevel, int iloop) if (ilevel == nlevels_respa-1) post_force(vflag); else { + Region *region = NULL; + if (iregion >= 0) { + region = domain->regions[iregion]; + region->prematch(); + } + + double **x = atom->x; double **f = atom->f; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { + if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; if (xstyle) f[i][0] = 0.0; if (ystyle) f[i][1] = 0.0; if (zstyle) f[i][2] = 0.0; diff --git a/src/fix_wall_region.cpp b/src/fix_wall_region.cpp index 315c2a0924..de8e894c4d 100644 --- a/src/fix_wall_region.cpp +++ b/src/fix_wall_region.cpp @@ -192,6 +192,8 @@ void FixWallRegion::post_force(int vflag) int nlocal = atom->nlocal; Region *region = domain->regions[iregion]; + region->prematch(); + int onflag = 0; // region->match() insures particle is in region or on surface, else error diff --git a/src/group.cpp b/src/group.cpp index 3b494a0b42..17a79fefb1 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -162,6 +162,7 @@ void Group::assign(int narg, char **arg) int iregion = domain->find_region(arg[2]); if (iregion == -1) error->all(FLERR,"Group region ID does not exist"); domain->regions[iregion]->init(); + domain->regions[iregion]->prematch(); for (i = 0; i < nlocal; i++) if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) @@ -661,6 +662,7 @@ bigint Group::count(int igroup, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; int *mask = atom->mask; @@ -715,6 +717,7 @@ double Group::mass(int igroup, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; double *mass = atom->mass; @@ -769,6 +772,7 @@ double Group::charge(int igroup, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; double *q = atom->q; @@ -837,6 +841,7 @@ void Group::bounds(int igroup, double *minmax, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double extent[6]; extent[0] = extent[2] = extent[4] = BIG; @@ -936,6 +941,7 @@ void Group::xcm(int igroup, double masstotal, double *cm, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; int *mask = atom->mask; @@ -1035,6 +1041,7 @@ void Group::vcm(int igroup, double masstotal, double *cm, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; double **v = atom->v; @@ -1106,6 +1113,7 @@ void Group::fcm(int igroup, double *cm, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; double **f = atom->f; @@ -1168,6 +1176,7 @@ double Group::ke(int igroup, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; double **v = atom->v; @@ -1246,6 +1255,7 @@ double Group::gyration(int igroup, double masstotal, double *cm, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; int *mask = atom->mask; @@ -1327,6 +1337,7 @@ void Group::angmom(int igroup, double *cm, double *lmom, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; double **v = atom->v; @@ -1405,6 +1416,7 @@ void Group::torque(int igroup, double *cm, double *tq, int iregion) { int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; double **f = atom->f; @@ -1492,6 +1504,7 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion) int groupbit = bitmask[igroup]; Region *region = domain->regions[iregion]; + region->prematch(); double **x = atom->x; int *mask = atom->mask; @@ -1533,7 +1546,7 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion) /* ---------------------------------------------------------------------- compute angular velocity omega from L = Iw, inverting I to solve for w - really not a group operation, but L and I were computed for a group + really not a group/region operation, but L,I were computed for a group/region ------------------------------------------------------------------------- */ void Group::omega(double *angmom, double inertia[3][3], double *w) diff --git a/src/region.cpp b/src/region.cpp index 2d26369980..3b836d0cfc 100644 --- a/src/region.cpp +++ b/src/region.cpp @@ -40,7 +40,6 @@ Region::Region(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) varshape = 0; xstr = ystr = zstr = tstr = NULL; dx = dy = dz = 0.0; - lastshape = lastdynamic = -1; } /* ---------------------------------------------------------------------- */ @@ -87,14 +86,30 @@ void Region::init() } /* ---------------------------------------------------------------------- - return 1 if region is dynamic, 0 if static + return 1 if region is dynamic (moves/rotates) or has variable shape + else return 0 if static only primitive regions define it here union/intersect regions have their own dynamic_check() ------------------------------------------------------------------------- */ int Region::dynamic_check() { - return dynamic; + if (dynamic || varshape) return 1; + return 0; +} + +/* ---------------------------------------------------------------------- + called before looping over atoms with match() or surface() + this insures any variables used by region are invoked once per timestep + also insures variables are invoked by all procs even those w/out atoms + necessary if equal-style variable invokes global operation + with MPI_Allreduce, e.g. xcm() or count() +------------------------------------------------------------------------- */ + +void Region::prematch() +{ + if (varshape) shape_update(); + if (dynamic) pretransform(); } /* ---------------------------------------------------------------------- @@ -111,13 +126,7 @@ int Region::dynamic_check() int Region::match(double x, double y, double z) { - if (varshape && update->ntimestep != lastshape) { - shape_update(); - lastshape = update->ntimestep; - } - if (dynamic) inverse_transform(x,y,z); - return !(inside(x,y,z) ^ interior); } @@ -139,11 +148,6 @@ int Region::surface(double x, double y, double z, double cutoff) double xs,ys,zs; double xnear[3],xorig[3]; - if (varshape && update->ntimestep != lastshape) { - shape_update(); - lastshape = update->ntimestep; - } - if (dynamic) { xorig[0] = x; xorig[1] = y; @@ -190,6 +194,21 @@ void Region::add_contact(int n, double *x, double xp, double yp, double zp) contact[n].delz = delz; } +/* ---------------------------------------------------------------------- + pre-compute dx,dy,dz and theta for a moving/rotating region + called once for the region before per-atom loop, via prematch() +------------------------------------------------------------------------- */ + +void Region::pretransform() +{ + if (moveflag) { + if (xstr) dx = input->variable->compute_equal(xvar); + if (ystr) dy = input->variable->compute_equal(yvar); + if (zstr) dz = input->variable->compute_equal(zvar); + } + if (rotateflag) theta = input->variable->compute_equal(tvar); +} + /* ---------------------------------------------------------------------- transform a point x,y,z in region space to moved space rotate first (around original P), then displace @@ -197,24 +216,12 @@ void Region::add_contact(int n, double *x, double xp, double yp, double zp) void Region::forward_transform(double &x, double &y, double &z) { - if (rotateflag) { - if (update->ntimestep != lastdynamic) - theta = input->variable->compute_equal(tvar); - rotate(x,y,z,theta); - } - + if (rotateflag) rotate(x,y,z,theta); if (moveflag) { - if (update->ntimestep != lastdynamic) { - if (xstr) dx = input->variable->compute_equal(xvar); - if (ystr) dy = input->variable->compute_equal(yvar); - if (zstr) dz = input->variable->compute_equal(zvar); - } x += dx; y += dy; z += dz; } - - lastdynamic = update->ntimestep; } /* ---------------------------------------------------------------------- @@ -225,23 +232,11 @@ void Region::forward_transform(double &x, double &y, double &z) void Region::inverse_transform(double &x, double &y, double &z) { if (moveflag) { - if (update->ntimestep != lastdynamic) { - if (xstr) dx = input->variable->compute_equal(xvar); - if (ystr) dy = input->variable->compute_equal(yvar); - if (zstr) dz = input->variable->compute_equal(zvar); - } x -= dx; y -= dy; z -= dz; } - - if (rotateflag) { - if (update->ntimestep != lastdynamic) - theta = input->variable->compute_equal(tvar); - rotate(x,y,z,-theta); - } - - lastdynamic = update->ntimestep; + if (rotateflag) rotate(x,y,z,-theta); } /* ---------------------------------------------------------------------- diff --git a/src/region.h b/src/region.h index dc39eff5da..f9f1841f35 100644 --- a/src/region.h +++ b/src/region.h @@ -46,6 +46,7 @@ class Region : protected Pointers { // called by other classes to check point versus region + void prematch(); int match(double, double, double); int surface(double, double, double, double); @@ -68,8 +69,8 @@ class Region : protected Pointers { char *xstr,*ystr,*zstr,*tstr; int xvar,yvar,zvar,tvar; double dx,dy,dz,theta; - bigint lastshape,lastdynamic; + void pretransform(); void forward_transform(double &, double &, double &); void inverse_transform(double &, double &, double &); void rotate(double &, double &, double &, double); diff --git a/src/set.cpp b/src/set.cpp index 23a654af67..5f3af82122 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -480,6 +480,7 @@ void Set::selection(int n) } else if (style == REGION_SELECT) { int iregion = domain->find_region(id); if (iregion == -1) error->all(FLERR,"Set region ID does not exist"); + domain->regions[iregion]->prematch(); double **x = atom->x; for (int i = 0; i < n; i++)