diff --git a/src/domain.cpp b/src/domain.cpp index 1180dfb6c6..660c53eace 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -119,6 +119,10 @@ void Domain::init() deform_groupbit = modify->fix[i]->groupbit; } else deform_remap = 0; } + + // region inits + + for (int i = 0; i < nregion; i++) regions[i]->init(); } /* ---------------------------------------------------------------------- diff --git a/src/fix_wall_harmonic.cpp b/src/fix_wall_harmonic.cpp new file mode 100644 index 0000000000..d8b18aebd5 --- /dev/null +++ b/src/fix_wall_harmonic.cpp @@ -0,0 +1,63 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "fix_wall_harmonic.h" +#include "atom.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixWallHarmonic::FixWallHarmonic(LAMMPS *lmp, int narg, char **arg) : + FixWall(lmp, narg, arg) {} + +/* ---------------------------------------------------------------------- + interaction of all particles in group with all 6 walls (if defined) + error if any particle is on or behind wall +------------------------------------------------------------------------- */ + +void FixWallHarmonic::wall_particle(int m, double coord) +{ + double delta,dr,fwall; + + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int dim = m/2; + int side = m % 2; + if (side == 0) side = -1; + + int onflag = 0; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (side < 0) delta = x[i][dim] - coord; + else delta = coord - x[i][dim]; + if (delta >= cutoff[m]) continue; + if (delta <= 0.0) { + onflag = 1; + continue; + } + dr = cutoff[m]-delta; + fwall = side * 2.0*epsilon[m]*dr; + f[i][dim] -= fwall; + ewall[0] += epsilon[m]*dr*dr; + ewall[m+1] += fwall; + } + + if (onflag) error->one("Particle on or inside fix wall surface"); +} diff --git a/src/fix_wall_harmonic.h b/src/fix_wall_harmonic.h new file mode 100644 index 0000000000..debfd64bd1 --- /dev/null +++ b/src/fix_wall_harmonic.h @@ -0,0 +1,33 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef FIX_WALL_HARMONIC_H +#define FIX_WALL_HARMONIC_H + +#include "fix_wall.h" + +namespace LAMMPS_NS { + +class FixWallHarmonic : public FixWall { + public: + FixWallHarmonic(class LAMMPS *, int, char **); + void precompute(int) {} + void wall_particle(int, double); + + private: + double offset[6]; +}; + +} + +#endif diff --git a/src/fix_wall_region.cpp b/src/fix_wall_region.cpp index b738890b94..e977d557b9 100644 --- a/src/fix_wall_region.cpp +++ b/src/fix_wall_region.cpp @@ -204,7 +204,7 @@ void FixWallRegion::post_force(int vflag) if (style == COLLOID) tooclose = shape[type[i]][0]; else tooclose = 0.0; - n = region->surface(x[i],cutoff); + n = region->surface(x[i][0],x[i][1],x[i][2],cutoff); for (m = 0; m < n; m++) { if (region->contact[m].r <= tooclose) { diff --git a/src/region.cpp b/src/region.cpp index b034e1a2f0..fa1f5a893d 100644 --- a/src/region.cpp +++ b/src/region.cpp @@ -15,12 +15,15 @@ #include "stdlib.h" #include "string.h" #include "region.h" +#include "update.h" #include "domain.h" #include "lattice.h" #include "error.h" using namespace LAMMPS_NS; +enum{NONE,LINEAR,WIGGLE,ROTATE,VARIABLE}; + /* ---------------------------------------------------------------------- */ Region::Region(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) @@ -32,6 +35,8 @@ Region::Region(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) n = strlen(arg[1]) + 1; style = new char[n]; strcpy(style,arg[1]); + + time_origin = update->ntimestep; } /* ---------------------------------------------------------------------- */ @@ -42,6 +47,13 @@ Region::~Region() delete [] style; } +/* ---------------------------------------------------------------------- */ + +void Region::init() +{ + dt = update->dt; +} + /* ---------------------------------------------------------------------- parse optional parameters at end of region input line ------------------------------------------------------------------------- */ @@ -54,6 +66,7 @@ void Region::options(int narg, char **arg) interior = 1; scaleflag = 1; + dynamic = NONE; int iarg = 0; while (iarg < narg) { @@ -69,9 +82,41 @@ void Region::options(int narg, char **arg) else if (strcmp(arg[iarg+1],"out") == 0) interior = 0; else error->all("Illegal region command"); iarg += 2; + } else if (strcmp(arg[iarg],"linear") == 0) { + if (iarg+4 > narg) error->all("Illegal region command"); + vx = atof(arg[iarg+1]); + vy = atof(arg[iarg+2]); + vz = atof(arg[iarg+3]); + dynamic = LINEAR; + iarg += 4; + } else if (strcmp(arg[iarg],"wiggle") == 0) { + if (iarg+5 > narg) error->all("Illegal region command"); + ax = atof(arg[iarg+1]); + ay = atof(arg[iarg+2]); + az = atof(arg[iarg+3]); + period = atof(arg[iarg+4]); + dynamic = WIGGLE; + iarg += 5; + } else if (strcmp(arg[iarg],"rotate") == 0) { + if (iarg+8 > narg) error->all("Illegal region command"); + point[0] = atof(arg[iarg+1]); + point[1] = atof(arg[iarg+2]); + point[2] = atof(arg[iarg+3]); + axis[0] = atof(arg[iarg+4]); + axis[1] = atof(arg[iarg+5]); + axis[2] = atof(arg[iarg+6]); + period = atof(arg[iarg+7]); + dynamic = ROTATE; + iarg += 8; } else error->all("Illegal region command"); } + // error check + + if (dynamic && + (strcmp(style,"union") == 0 || strcmp(style,"intersect") == 0)) + error->all("Region union or intersect cannot be dynamic"); + // setup scaling if (scaleflag && domain->lattice == NULL) @@ -83,16 +128,128 @@ void Region::options(int narg, char **arg) zscale = domain->lattice->zlattice; } else xscale = yscale = zscale = 1.0; + + if (dynamic == LINEAR) { + vx *= xscale; + vy *= yscale; + vz *= zscale; + } else if (dynamic == WIGGLE) { + ax *= xscale; + ay *= yscale; + az *= zscale; + } else if (dynamic == ROTATE) { + point[0] *= xscale; + point[1] *= yscale; + point[2] *= zscale; + } + + if (dynamic == WIGGLE || dynamic == ROTATE) { + double PI = 4.0 * atan(1.0); + omega_rotate = 2.0*PI / period; + } + + // runit = unit vector along rotation axis + + if (dynamic == ROTATE) { + double len = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]); + if (len == 0.0) + error->all("Region cannot have 0 length rotation vector"); + runit[0] = axis[0]/len; + runit[1] = axis[1]/len; + runit[2] = axis[2]/len; + } +} + +/* ---------------------------------------------------------------------- + determine if point x,y,z is a match to region volume + XOR computes 0 if 2 args are the same, 1 if different + note that inside() returns 1 for points on surface of region + thus point on surface of exterior region will not match + if region is dynamic, apply inverse of region change to x +------------------------------------------------------------------------- */ + +int Region::match(double x, double y, double z) +{ + double a[3],b[3],c[3],d[3],disp[0]; + + if (dynamic) { + double delta = (update->ntimestep - time_origin) * dt; + if (dynamic == LINEAR) { + x -= vx*delta; + y -= vy*delta; + z -= vz*delta; + } else if (dynamic == WIGGLE) { + double arg = omega_rotate * delta; + double sine = sin(arg); + x -= ax*sine; + y -= ay*sine; + z -= az*sine; + } else if (dynamic == ROTATE) { + double angle = -omega_rotate*delta; + rotate(x,y,z,angle); + } + } + + return !(inside(x,y,z) ^ interior); } /* ---------------------------------------------------------------------- generate list of contact points for interior or exterior regions + if region is dynamic: + change x by inverse of region change + change contact point by region change ------------------------------------------------------------------------- */ -int Region::surface(double *x, double cutoff) +int Region::surface(double x, double y, double z, double cutoff) { - if (interior) return surface_interior(x,cutoff); - return surface_exterior(x,cutoff); + int ncontact; + double xnear[3],xhold[3]; + + if (dynamic) { + double delta = (update->ntimestep - time_origin) * dt; + if (dynamic == LINEAR) { + x -= vx*delta; + y -= vy*delta; + z -= vz*delta; + } else if (dynamic == WIGGLE) { + double arg = omega_rotate * delta; + double sine = sin(arg); + x -= ax*sine; + y -= ay*sine; + z -= az*sine; + } else if (dynamic == ROTATE) { + xhold[0] = x; + xhold[1] = y; + xhold[2] = z; + double angle = -omega_rotate*delta; + rotate(x,y,z,angle); + } + } + + xnear[0] = x; + xnear[1] = y; + xnear[2] = z; + + if (interior) ncontact = surface_interior(xnear,cutoff); + else ncontact = surface_exterior(xnear,cutoff); + + if (dynamic && ncontact) { + double delta = (update->ntimestep - time_origin) * dt; + if (dynamic == ROTATE) { + for (int i = 0; i < ncontact; i++) { + x -= contact[i].delx; + y -= contact[i].dely; + z -= contact[i].delz; + double angle = omega_rotate*delta; + rotate(x,y,z,angle); + contact[i].delx = xhold[0] - x; + contact[i].dely = xhold[1] - y; + contact[i].delz = xhold[2] - z; + } + } + } + + return ncontact; } /* ---------------------------------------------------------------------- @@ -111,3 +268,47 @@ void Region::add_contact(int n, double *x, double xp, double yp, double zp) contact[n].dely = dely; contact[n].delz = delz; } + +/* ---------------------------------------------------------------------- + rotate x,y,z by angle via right-hand rule around point and runit normal + sign of angle determines whether rotating forward/backward in time + return updated x,y,z + P = point = vector = point of rotation + R = vector = axis of rotation + w = omega of rotation (from period) + X0 = x,y,z = initial coord of atom + R0 = runit = unit vector for R + C = (X0 dot R0) R0 = projection of atom coord onto R + D = X0 - P = vector from P to X0 + A = D - C = vector from R line to X0 + B = R0 cross A = vector perp to A in plane of rotation + A,B define plane of circular rotation around R line + x,y,z = P + C + A cos(w*dt) + B sin(w*dt) +------------------------------------------------------------------------- */ + +void Region::rotate(double &x, double &y, double &z, double angle) +{ + double a[3],b[3],c[3],d[3],disp[0]; + + double sine = sin(angle); + double cosine = cos(angle); + double x0dotr = x*runit[0] + y*runit[1] + z*runit[2]; + c[0] = x0dotr * runit[0]; + c[1] = x0dotr * runit[1]; + c[2] = x0dotr * runit[2]; + d[0] = x - point[0]; + d[1] = y - point[1]; + d[2] = z - point[2]; + a[0] = d[0] - c[0]; + a[1] = d[1] - c[1]; + a[2] = d[2] - c[2]; + b[0] = runit[1]*a[2] - runit[2]*a[1]; + b[1] = runit[2]*a[0] - runit[0]*a[2]; + b[2] = runit[0]*a[1] - runit[1]*a[0]; + disp[0] = a[0]*cosine + b[0]*sine; + disp[1] = a[1]*cosine + b[1]*sine; + disp[2] = a[2]*cosine + b[2]*sine; + x = point[0] + c[0] + disp[0]; + y = point[1] + c[1] + disp[1]; + z = point[2] + c[2] + disp[2]; +} diff --git a/src/region.h b/src/region.h index 1bad9faf84..cd796398dc 100644 --- a/src/region.h +++ b/src/region.h @@ -28,6 +28,7 @@ class Region : protected Pointers { double extent_ylo,extent_yhi; double extent_zlo,extent_zhi; int bboxflag; // 1 if bounding box is computable + int dynamic; // 1 if region changes over time // contact = particle near region surface @@ -40,14 +41,26 @@ class Region : protected Pointers { Region(class LAMMPS *, int, char **); virtual ~Region(); - virtual int match(double, double, double) = 0; - int surface(double *, double); + void init(); + int match(double, double, double); + int surface(double, double, double, double); + + virtual int inside(double, double, double) = 0; virtual int surface_interior(double *, double) = 0; virtual int surface_exterior(double *, double) = 0; protected: void options(int, char **); void add_contact(int, double *, double, double, double); + + private: + int time_origin; + double dt,period,omega_rotate; + double vx,vy,vz; + double ax,ay,az; + double point[3],axis[3],runit[3]; + + void rotate(double &, double &, double &, double); }; } diff --git a/src/region_block.cpp b/src/region_block.cpp index f63dc66883..21321d25c5 100644 --- a/src/region_block.cpp +++ b/src/region_block.cpp @@ -29,7 +29,7 @@ using namespace LAMMPS_NS; RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) { options(narg-8,&arg[8]); - + if (strcmp(arg[2],"INF") == 0 || strcmp(arg[2],"EDGE") == 0) { if (domain->box_exist == 0) error->all("Cannot use region INF or EDGE when box does not exist"); @@ -113,13 +113,11 @@ RegBlock::~RegBlock() inside = 0 if x,y,z is outside and not on surface ------------------------------------------------------------------------- */ -int RegBlock::match(double x, double y, double z) +int RegBlock::inside(double x, double y, double z) { - int inside = 0; if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi) - inside = 1; - - return !(inside ^ interior); // 1 if same, 0 if different + return 1; + return 0; } /* ---------------------------------------------------------------------- diff --git a/src/region_block.h b/src/region_block.h index 5f61764398..1cd05ae00b 100644 --- a/src/region_block.h +++ b/src/region_block.h @@ -24,7 +24,7 @@ class RegBlock : public Region { public: RegBlock(class LAMMPS *, int, char **); ~RegBlock(); - int match(double, double, double); + int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); diff --git a/src/region_cone.cpp b/src/region_cone.cpp index dc53f6ddfb..0c45de6adb 100644 --- a/src/region_cone.cpp +++ b/src/region_cone.cpp @@ -162,7 +162,7 @@ RegCone::~RegCone() inside = 0 if x,y,z is outside and not on surface ------------------------------------------------------------------------- */ -int RegCone::match(double x, double y, double z) +int RegCone::inside(double x, double y, double z) { double del1,del2,dist; double currentradius; @@ -193,7 +193,7 @@ int RegCone::match(double x, double y, double z) else inside = 0; } - return !(inside ^ interior); // 1 if same, 0 if different + return inside; } /* ---------------------------------------------------------------------- diff --git a/src/region_cone.h b/src/region_cone.h index 791a6d3f56..bd15eb3eb8 100644 --- a/src/region_cone.h +++ b/src/region_cone.h @@ -22,7 +22,7 @@ class RegCone : public Region { public: RegCone(class LAMMPS *, int, char **); ~RegCone(); - int match(double, double, double); + int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); diff --git a/src/region_cylinder.cpp b/src/region_cylinder.cpp index c81954c137..1c342052e1 100644 --- a/src/region_cylinder.cpp +++ b/src/region_cylinder.cpp @@ -150,7 +150,7 @@ RegCylinder::~RegCylinder() inside = 0 if x,y,z is outside and not on surface ------------------------------------------------------------------------- */ -int RegCylinder::match(double x, double y, double z) +int RegCylinder::inside(double x, double y, double z) { double del1,del2,dist; int inside; @@ -175,7 +175,7 @@ int RegCylinder::match(double x, double y, double z) else inside = 0; } - return !(inside ^ interior); // 1 if same, 0 if different + return inside; } /* ---------------------------------------------------------------------- diff --git a/src/region_cylinder.h b/src/region_cylinder.h index 25273c114a..5d033b077d 100644 --- a/src/region_cylinder.h +++ b/src/region_cylinder.h @@ -24,7 +24,7 @@ class RegCylinder : public Region { public: RegCylinder(class LAMMPS *, int, char **); ~RegCylinder(); - int match(double, double, double); + int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); diff --git a/src/region_intersect.cpp b/src/region_intersect.cpp index 2ecfe70713..ac773165a9 100644 --- a/src/region_intersect.cpp +++ b/src/region_intersect.cpp @@ -98,17 +98,15 @@ RegIntersect::~RegIntersect() else inside = 0 ------------------------------------------------------------------------- */ -int RegIntersect::match(double x, double y, double z) +int RegIntersect::inside(double x, double y, double z) { int ilist; Region **regions = domain->regions; for (ilist = 0; ilist < nregion; ilist++) if (!regions[list[ilist]]->match(x,y,z)) break; - int inside = 0; // inside if matched all regions - if (ilist == nregion) inside = 1; - - return !(inside ^ interior); // 1 if same, 0 if different + if (ilist == nregion) return 1; + return 0; } /* ---------------------------------------------------------------------- @@ -127,7 +125,7 @@ int RegIntersect::surface_interior(double *x, double cutoff) for (ilist = 0; ilist < nregion; ilist++) { iregion = list[ilist]; - ncontacts = regions[iregion]->surface(x,cutoff); + ncontacts = regions[iregion]->surface(x[0],x[1],x[2],cutoff); for (m = 0; m < ncontacts; m++) { xs = x[0] - regions[iregion]->contact[m].delx; ys = x[1] - regions[iregion]->contact[m].dely; @@ -172,7 +170,7 @@ int RegIntersect::surface_exterior(double *x, double cutoff) for (ilist = 0; ilist < nregion; ilist++) { iregion = list[ilist]; - ncontacts = regions[iregion]->surface(x,cutoff); + ncontacts = regions[iregion]->surface(x[0],x[1],x[2],cutoff); for (m = 0; m < ncontacts; m++) { xs = x[0] - regions[iregion]->contact[m].delx; ys = x[1] - regions[iregion]->contact[m].dely; diff --git a/src/region_intersect.h b/src/region_intersect.h index 8d7d76bf23..8c09277d57 100644 --- a/src/region_intersect.h +++ b/src/region_intersect.h @@ -22,7 +22,7 @@ class RegIntersect : public Region { public: RegIntersect(class LAMMPS *, int, char **); ~RegIntersect(); - int match(double, double, double); + int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); diff --git a/src/region_plane.cpp b/src/region_plane.cpp new file mode 100644 index 0000000000..45b6f103ab --- /dev/null +++ b/src/region_plane.cpp @@ -0,0 +1,110 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "region_plane.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) : + Region(lmp, narg, arg) +{ + options(narg-8,&arg[8]); + + xp = xscale*atof(arg[2]); + yp = yscale*atof(arg[3]); + zp = zscale*atof(arg[4]); + normal[0] = xscale*atof(arg[5]); + normal[1] = yscale*atof(arg[6]); + normal[2] = zscale*atof(arg[7]); + + // enforce unit normal + + double rsq = normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]; + if (rsq == 0.0) error->all("Illegal region plane command"); + normal[0] /= sqrt(rsq); + normal[1] /= sqrt(rsq); + normal[2] /= sqrt(rsq); + + // plane has no bounding box + + bboxflag = 0; + + cmax = 1; + contact = new Contact[cmax]; +} + +/* ---------------------------------------------------------------------- */ + +RegPlane::~RegPlane() +{ + delete [] contact; +} + +/* ---------------------------------------------------------------------- + inside = 1 if x,y,z is on normal side of plane or on plane + inside = 0 if x,y,z is on non-normal side of plane and not on plane + x,y,z is inside if (x-xp) dot normal >= 0 +------------------------------------------------------------------------- */ + +int RegPlane::inside(double x, double y, double z) +{ + double dot = (x-xp)*normal[0] + (y-yp)*normal[1] + (z-zp)*normal[2]; + + if (dot >= 0.0) return 1; + return 0; +} + +/* ---------------------------------------------------------------------- + one contact if 0 <= x < cutoff from normal side of plane + no contact if on other side (possible if called from union/intersect) + delxyz = vector from nearest projected point on plane to x +------------------------------------------------------------------------- */ + +int RegPlane::surface_interior(double *x, double cutoff) +{ + double dot = (x[0]-xp)*normal[0] + (x[1]-yp)*normal[1] + (x[2]-zp)*normal[2]; + if (dot < cutoff && dot >= 0.0) { + contact[0].r = dot; + contact[0].delx = dot*normal[0]; + contact[0].dely = dot*normal[1]; + contact[0].delz = dot*normal[2]; + return 1; + } + return 0; +} + +/* ---------------------------------------------------------------------- + one contact if 0 <= x < cutoff from non-normal side of plane + no contact if on other side (possible if called from union/intersect) + delxyz = vector from nearest projected point on plane to x +------------------------------------------------------------------------- */ + +int RegPlane::surface_exterior(double *x, double cutoff) +{ + double dot = (x[0]-xp)*normal[0] + (x[1]-yp)*normal[1] + (x[2]-zp)*normal[2]; + dot = -dot; + if (dot < cutoff && dot >= 0.0) { + contact[0].r = dot; + contact[0].delx = dot*normal[0]; + contact[0].dely = dot*normal[1]; + contact[0].delz = dot*normal[2]; + return 1; + } + return 0; +} diff --git a/src/region_plane.h b/src/region_plane.h new file mode 100644 index 0000000000..0cea928058 --- /dev/null +++ b/src/region_plane.h @@ -0,0 +1,36 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef REGION_PLANE_H +#define REGION_PLANE_H + +#include "region.h" + +namespace LAMMPS_NS { + +class RegPlane : public Region { + public: + RegPlane(class LAMMPS *, int, char **); + ~RegPlane(); + int inside(double, double, double); + int surface_interior(double *, double); + int surface_exterior(double *, double); + + private: + double xp,yp,zp; + double normal[3]; +}; + +} + +#endif diff --git a/src/region_prism.cpp b/src/region_prism.cpp index 3ccb3ea181..0be3a83b94 100644 --- a/src/region_prism.cpp +++ b/src/region_prism.cpp @@ -240,17 +240,15 @@ RegPrism::~RegPrism() xyz/lo = lower-left corner of prism ------------------------------------------------------------------------- */ -int RegPrism::match(double x, double y, double z) +int RegPrism::inside(double x, double y, double z) { double a = hinv[0][0]*(x-xlo) + hinv[0][1]*(y-ylo) + hinv[0][2]*(z-zlo); double b = hinv[1][1]*(y-ylo) + hinv[1][2]*(z-zlo); double c = hinv[2][2]*(z-zlo); - int inside = 0; if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0 && c >= 0.0 && c <= 1.0) - inside = 1; - - return !(inside ^ interior); // 1 if same, 0 if different + return 1; + return 0; } /* ---------------------------------------------------------------------- diff --git a/src/region_prism.h b/src/region_prism.h index 39c87452ff..9b9229c88a 100644 --- a/src/region_prism.h +++ b/src/region_prism.h @@ -24,7 +24,7 @@ class RegPrism : public Region { public: RegPrism(class LAMMPS *, int, char **); ~RegPrism(); - int match(double, double, double); + int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); diff --git a/src/region_sphere.cpp b/src/region_sphere.cpp index 95f6fdfe15..351421f94d 100644 --- a/src/region_sphere.cpp +++ b/src/region_sphere.cpp @@ -63,21 +63,19 @@ RegSphere::~RegSphere() inside = 0 if x,y,z is outside and not on surface ------------------------------------------------------------------------- */ -int RegSphere::match(double x, double y, double z) +int RegSphere:: inside(double x, double y, double z) { double delx = x - xc; double dely = y - yc; double delz = z - zc; double r = sqrt(delx*delx + dely*dely + delz*delz); - int inside = 0; - if (r <= radius) inside = 1; - - return !(inside ^ interior); // 1 if same, 0 if different + if (r <= radius) return 1; + return 0; } /* ---------------------------------------------------------------------- - one contact if 0 <= x < cutoff away from inner surface of sphere + one contact if 0 <= x < cutoff from inner surface of sphere no contact if outside (possible if called from union/intersect) delxyz = vector from nearest point on sphere to x special case: no contact if x is at center of sphere @@ -103,7 +101,7 @@ int RegSphere::surface_interior(double *x, double cutoff) } /* ---------------------------------------------------------------------- - one contact if 0 <= x < cutoff away from outer surface of sphere + one contact if 0 <= x < cutoff from outer surface of sphere no contact if inside (possible if called from union/intersect) delxyz = vector from nearest point on sphere to x ------------------------------------------------------------------------- */ diff --git a/src/region_sphere.h b/src/region_sphere.h index 09d9760873..2d06f12403 100644 --- a/src/region_sphere.h +++ b/src/region_sphere.h @@ -22,7 +22,7 @@ class RegSphere : public Region { public: RegSphere(class LAMMPS *, int, char **); ~RegSphere(); - int match(double, double, double); + int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); diff --git a/src/region_union.cpp b/src/region_union.cpp index 4d339cecca..89947acab4 100644 --- a/src/region_union.cpp +++ b/src/region_union.cpp @@ -90,17 +90,15 @@ RegUnion::~RegUnion() else inside = 0 ------------------------------------------------------------------------- */ -int RegUnion::match(double x, double y, double z) +int RegUnion::inside(double x, double y, double z) { int ilist; Region **regions = domain->regions; for (ilist = 0; ilist < nregion; ilist++) if (regions[list[ilist]]->match(x,y,z)) break; - int inside = 1; // inside if matched any region - if (ilist == nregion) inside = 0; - - return !(inside ^ interior); // 1 if same, 0 if different + if (ilist == nregion) return 0; + return 1; } /* ---------------------------------------------------------------------- @@ -119,7 +117,7 @@ int RegUnion::surface_interior(double *x, double cutoff) for (ilist = 0; ilist < nregion; ilist++) { iregion = list[ilist]; - ncontacts = regions[iregion]->surface(x,cutoff); + ncontacts = regions[iregion]->surface(x[0],x[1],x[2],cutoff); for (m = 0; m < ncontacts; m++) { xs = x[0] - regions[iregion]->contact[m].delx; ys = x[1] - regions[iregion]->contact[m].dely; @@ -164,7 +162,7 @@ int RegUnion::surface_exterior(double *x, double cutoff) for (ilist = 0; ilist < nregion; ilist++) { iregion = list[ilist]; - ncontacts = regions[iregion]->surface(x,cutoff); + ncontacts = regions[iregion]->surface(x[0],x[1],x[2],cutoff); for (m = 0; m < ncontacts; m++) { xs = x[0] - regions[iregion]->contact[m].delx; ys = x[1] - regions[iregion]->contact[m].dely; diff --git a/src/region_union.h b/src/region_union.h index 01ebe79ce5..d73ca85ef1 100644 --- a/src/region_union.h +++ b/src/region_union.h @@ -22,7 +22,7 @@ class RegUnion : public Region { public: RegUnion(class LAMMPS *, int, char **); ~RegUnion(); - int match(double, double, double); + int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); diff --git a/src/update.cpp b/src/update.cpp index 87cb472b94..0665031684 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -18,6 +18,8 @@ #include "force.h" #include "modify.h" #include "fix.h" +#include "domain.h" +#include "region.h" #include "compute.h" #include "output.h" #include "memory.h" @@ -222,6 +224,7 @@ void Update::create_minimize(int narg, char **arg) reset timestep from input script do not allow dump files or a restart to be defined do not allow any timestep-dependent fixes to be defined + do not allow any dynamic regions to be defined reset eflag/vflag global so nothing will think eng/virial are current reset invoked flags of computes, so nothing will think they are current clear timestep list of computes that store future invocation times @@ -241,6 +244,10 @@ void Update::reset_timestep(int narg, char **arg) if (modify->fix[i]->time_depend) error->all("Cannot reset timestep with a time-dependent fix defined"); + for (int i = 0; i < domain->nregion; i++) + if (domain->regions[i]->dynamic) + error->all("Cannot reset timestep with a dynamic region defined"); + eflag_global = vflag_global = -1; for (int i = 0; i < modify->ncompute; i++) {