diff --git a/src/region.h b/src/region.h index 8f7959025..0a1d888a9 100644 --- a/src/region.h +++ b/src/region.h @@ -1,116 +1,115 @@ /* ---------------------------------------------------------------------- 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 LMP_REGION_H #define LMP_REGION_H #include "pointers.h" namespace LAMMPS_NS { class Region : protected Pointers { public: char *id,*style; int interior; // 1 for interior, 0 for exterior int scaleflag; // 1 for lattice, 0 for box double xscale,yscale,zscale; // scale factors for box/lattice units double extent_xlo,extent_xhi; // bounding box on region double extent_ylo,extent_yhi; double extent_zlo,extent_zhi; int bboxflag; // 1 if bounding box is computable + int varshape; // 1 if region shape changes over time // contact = particle near region surface struct Contact { double r; // distance between particle & surf, r > 0.0 double delx,dely,delz; // vector from surface pt to particle }; Contact *contact; // list of contacts int cmax; // max # of contacts possible with region Region(class LAMMPS *, int, char **); virtual ~Region(); virtual void init(); virtual int dynamic_check(); // called by other classes to check point versus region int match(double, double, double); int surface(double, double, double, double); // implemented by each region, not called by other classes virtual int inside(double, double, double) = 0; virtual int surface_interior(double *, double) = 0; virtual int surface_exterior(double *, double) = 0; virtual void shape_update() {} protected: - int varshape; // 1 if region shape changes over time - void add_contact(int, double *, double, double, double); void options(int, char **); private: int dynamic; // 1 if region position/orientation changes over time int moveflag,rotateflag; // 1 if position/orientation changes double point[3],axis[3],runit[3]; char *xstr,*ystr,*zstr,*tstr; int xvar,yvar,zvar,tvar; double dx,dy,dz,theta; bigint lastshape,lastdynamic; void forward_transform(double &, double &, double &); void inverse_transform(double &, double &, double &); void rotate(double &, double &, double &, double); }; } #endif /* ERROR/WARNING messages: E: Variable name for region does not exist Self-explanatory. E: Variable for region is invalid style Only equal-style variables can be used. E: Variable for region is not equal style Self-explanatory. E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Region union or intersect cannot be dynamic The sub-regions can be dynamic, but not the combined region. E: Use of region with undefined lattice If units = lattice (the default) for the region command, then a lattice must first be defined via the lattice command. E: Region cannot have 0 length rotation vector Self-explanatory. */ diff --git a/src/region_intersect.cpp b/src/region_intersect.cpp index 736ec3c6e..bae73952d 100644 --- a/src/region_intersect.cpp +++ b/src/region_intersect.cpp @@ -1,208 +1,233 @@ /* ---------------------------------------------------------------------- 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 "stdlib.h" #include "string.h" #include "region_intersect.h" #include "domain.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ RegIntersect::RegIntersect(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) { if (narg < 5) error->all(FLERR,"Illegal region command"); int n = atoi(arg[2]); if (n < 2) error->all(FLERR,"Illegal region command"); options(narg-(n+3),&arg[n+3]); // build list of regions to intersect list = new int[n]; nregion = 0; int iregion; for (int iarg = 0; iarg < n; iarg++) { iregion = domain->find_region(arg[iarg+3]); if (iregion == -1) error->all(FLERR,"Region intersect region ID does not exist"); list[nregion++] = iregion; } - // extent of intersection of regions - // has bounding box if interior and any sub-region has bounding box + // this region is variable shape if any of sub-regions are Region **regions = domain->regions; + for (int ilist = 0; ilist < nregion; ilist++) + if (regions[list[ilist]]->varshape) varshape = 1; + + // extent of intersection of regions + // has bounding box if interior and any sub-region has bounding box bboxflag = 0; for (int ilist = 0; ilist < nregion; ilist++) if (regions[list[ilist]]->bboxflag == 1) bboxflag = 1; if (!interior) bboxflag = 0; if (bboxflag) { int first = 1; for (int ilist = 0; ilist < nregion; ilist++) { if (regions[list[ilist]]->bboxflag == 0) continue; if (first) { extent_xlo = regions[list[ilist]]->extent_xlo; extent_ylo = regions[list[ilist]]->extent_ylo; extent_zlo = regions[list[ilist]]->extent_zlo; extent_xhi = regions[list[ilist]]->extent_xhi; extent_yhi = regions[list[ilist]]->extent_yhi; extent_zhi = regions[list[ilist]]->extent_zhi; first = 0; } extent_xlo = MAX(extent_xlo,regions[list[ilist]]->extent_xlo); extent_ylo = MAX(extent_ylo,regions[list[ilist]]->extent_ylo); extent_zlo = MAX(extent_zlo,regions[list[ilist]]->extent_zlo); extent_xhi = MIN(extent_xhi,regions[list[ilist]]->extent_xhi); extent_yhi = MIN(extent_yhi,regions[list[ilist]]->extent_yhi); extent_zhi = MIN(extent_zhi,regions[list[ilist]]->extent_zhi); } } // possible contacts = sum of possible contacts in all sub-regions cmax = 0; for (int ilist = 0; ilist < nregion; ilist++) cmax += regions[list[ilist]]->cmax; contact = new Contact[cmax]; } /* ---------------------------------------------------------------------- */ RegIntersect::~RegIntersect() { delete [] list; delete [] contact; } +/* ---------------------------------------------------------------------- */ + +void RegIntersect::init() +{ + Region::init(); + Region **regions = domain->regions; + for (int ilist = 0; ilist < nregion; ilist++) + regions[list[ilist]]->init(); +} + /* ---------------------------------------------------------------------- return 1 if region is dynamic, 0 if static dynamic if any sub-region is dynamic, else static ------------------------------------------------------------------------- */ int RegIntersect::dynamic_check() { Region **regions = domain->regions; for (int ilist = 0; ilist < nregion; ilist++) if (regions[list[ilist]]->dynamic_check()) return 1; return 0; } /* ---------------------------------------------------------------------- inside = 1 if x,y,z is match() with all sub-regions else inside = 0 ------------------------------------------------------------------------- */ 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; if (ilist == nregion) return 1; return 0; } /* ---------------------------------------------------------------------- compute contacts with interior of intersection of sub-regions (1) compute contacts in each sub-region (2) only keep a contact if surface point is match() to all other regions ------------------------------------------------------------------------- */ int RegIntersect::surface_interior(double *x, double cutoff) { int m,ilist,jlist,iregion,jregion,ncontacts; double xs,ys,zs; Region **regions = domain->regions; int n = 0; for (ilist = 0; ilist < nregion; ilist++) { iregion = list[ilist]; 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; zs = x[2] - regions[iregion]->contact[m].delz; for (jlist = 0; jlist < nregion; jlist++) { if (jlist == ilist) continue; jregion = list[jlist]; if (!regions[jregion]->match(xs,ys,zs)) break; } if (jlist == nregion) { contact[n].r = regions[iregion]->contact[m].r; contact[n].delx = regions[iregion]->contact[m].delx; contact[n].dely = regions[iregion]->contact[m].dely; contact[n].delz = regions[iregion]->contact[m].delz; n++; } } } return n; } /* ---------------------------------------------------------------------- compute contacts with interior of intersection of sub-regions (1) flip interior/exterior flag of each sub-region (2) compute contacts in each sub-region (3) only keep a contact if surface point is not match() to all other regions (4) flip interior/exterior flags back to original settings this is effectively same algorithm as surface_interior() for RegUnion ------------------------------------------------------------------------- */ int RegIntersect::surface_exterior(double *x, double cutoff) { int m,ilist,jlist,iregion,jregion,ncontacts; double xs,ys,zs; Region **regions = domain->regions; int n = 0; for (ilist = 0; ilist < nregion; ilist++) regions[list[ilist]]->interior ^= 1; for (ilist = 0; ilist < nregion; ilist++) { iregion = list[ilist]; 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; zs = x[2] - regions[iregion]->contact[m].delz; for (jlist = 0; jlist < nregion; jlist++) { if (jlist == ilist) continue; jregion = list[jlist]; if (regions[jregion]->match(xs,ys,zs)) break; } if (jlist == nregion) { contact[n].r = regions[iregion]->contact[m].r; contact[n].delx = regions[iregion]->contact[m].delx; contact[n].dely = regions[iregion]->contact[m].dely; contact[n].delz = regions[iregion]->contact[m].delz; n++; } } } for (ilist = 0; ilist < nregion; ilist++) regions[list[ilist]]->interior ^= 1; return n; } + +/* ---------------------------------------------------------------------- + change region shape of all sub-regions +------------------------------------------------------------------------- */ + +void RegIntersect::shape_update() +{ + Region **regions = domain->regions; + for (int ilist = 0; ilist < nregion; ilist++) + regions[list[ilist]]->shape_update(); +} diff --git a/src/region_intersect.h b/src/region_intersect.h index a8bd460f3..62946a65c 100644 --- a/src/region_intersect.h +++ b/src/region_intersect.h @@ -1,58 +1,60 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef REGION_CLASS RegionStyle(intersect,RegIntersect) #else #ifndef LMP_REGION_INTERSECT_H #define LMP_REGION_INTERSECT_H #include "region.h" namespace LAMMPS_NS { class RegIntersect : public Region { public: RegIntersect(class LAMMPS *, int, char **); ~RegIntersect(); + void init(); int dynamic_check(); int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); + void shape_update(); private: int nregion; int *list; }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Region intersect region ID does not exist Self-explanatory. */ diff --git a/src/region_union.cpp b/src/region_union.cpp index ca1e693f3..3741e8fdd 100644 --- a/src/region_union.cpp +++ b/src/region_union.cpp @@ -1,199 +1,225 @@ /* ---------------------------------------------------------------------- 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 "stdlib.h" #include "string.h" #include "region_union.h" #include "domain.h" #include "error.h" using namespace LAMMPS_NS; #define BIG 1.0e20 /* ---------------------------------------------------------------------- */ RegUnion::RegUnion(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) { if (narg < 5) error->all(FLERR,"Illegal region command"); int n = atoi(arg[2]); if (n < 2) error->all(FLERR,"Illegal region command"); options(narg-(n+3),&arg[n+3]); // build list of regions to union list = new int[n]; nregion = 0; int iregion; for (int iarg = 0; iarg < n; iarg++) { iregion = domain->find_region(arg[iarg+3]); - if (iregion == -1) error->all(FLERR,"Region union region ID does not exist"); + if (iregion == -1) + error->all(FLERR,"Region union region ID does not exist"); list[nregion++] = iregion; } - // extent of union of regions - // has bounding box if interior and all sub-regions have bounding box + // this region is variable shape if any of sub-regions are Region **regions = domain->regions; + for (int ilist = 0; ilist < nregion; ilist++) + if (regions[list[ilist]]->varshape) varshape = 1; + + // extent of union of regions + // has bounding box if interior and all sub-regions have bounding box bboxflag = 1; for (int ilist = 0; ilist < nregion; ilist++) if (regions[list[ilist]]->bboxflag == 0) bboxflag = 0; if (!interior) bboxflag = 0; if (bboxflag) { extent_xlo = extent_ylo = extent_zlo = BIG; extent_xhi = extent_yhi = extent_zhi = -BIG; for (int ilist = 0; ilist < nregion; ilist++) { extent_xlo = MIN(extent_xlo,regions[list[ilist]]->extent_xlo); extent_ylo = MIN(extent_ylo,regions[list[ilist]]->extent_ylo); extent_zlo = MIN(extent_zlo,regions[list[ilist]]->extent_zlo); extent_xhi = MAX(extent_xhi,regions[list[ilist]]->extent_xhi); extent_yhi = MAX(extent_yhi,regions[list[ilist]]->extent_yhi); extent_zhi = MAX(extent_zhi,regions[list[ilist]]->extent_zhi); } } // possible contacts = sum of possible contacts in all sub-regions cmax = 0; for (int ilist = 0; ilist < nregion; ilist++) cmax += regions[list[ilist]]->cmax; contact = new Contact[cmax]; } /* ---------------------------------------------------------------------- */ RegUnion::~RegUnion() { delete [] list; delete [] contact; } +/* ---------------------------------------------------------------------- */ + +void RegUnion::init() +{ + Region::init(); + Region **regions = domain->regions; + for (int ilist = 0; ilist < nregion; ilist++) + regions[list[ilist]]->init(); +} + /* ---------------------------------------------------------------------- return 1 if region is dynamic, 0 if static dynamic if any sub-region is dynamic, else static ------------------------------------------------------------------------- */ int RegUnion::dynamic_check() { Region **regions = domain->regions; for (int ilist = 0; ilist < nregion; ilist++) if (regions[list[ilist]]->dynamic_check()) return 1; return 0; } /* ---------------------------------------------------------------------- inside = 1 if x,y,z is match() with any sub-region else inside = 0 ------------------------------------------------------------------------- */ 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; if (ilist == nregion) return 0; return 1; } /* ---------------------------------------------------------------------- compute contacts with interior of union of sub-regions (1) compute contacts in each sub-region (2) only keep a contact if surface point is not match() to all other regions ------------------------------------------------------------------------- */ int RegUnion::surface_interior(double *x, double cutoff) { int m,ilist,jlist,iregion,jregion,ncontacts; double xs,ys,zs; Region **regions = domain->regions; int n = 0; for (ilist = 0; ilist < nregion; ilist++) { iregion = list[ilist]; 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; zs = x[2] - regions[iregion]->contact[m].delz; for (jlist = 0; jlist < nregion; jlist++) { if (jlist == ilist) continue; jregion = list[jlist]; if (regions[jregion]->match(xs,ys,zs)) break; } if (jlist == nregion) { contact[n].r = regions[iregion]->contact[m].r; contact[n].delx = regions[iregion]->contact[m].delx; contact[n].dely = regions[iregion]->contact[m].dely; contact[n].delz = regions[iregion]->contact[m].delz; n++; } } } return n; } /* ---------------------------------------------------------------------- compute contacts with exterior of union of sub-regions (1) flip interior/exterior flag of each sub-region (2) compute contacts in each sub-region (3) only keep a contact if surface point is match() to all other regions (4) flip interior/exterior flags back to original settings this is effectively same algorithm as surface_interior() for RegIntersect ------------------------------------------------------------------------- */ int RegUnion::surface_exterior(double *x, double cutoff) { int m,ilist,jlist,iregion,jregion,ncontacts; double xs,ys,zs; Region **regions = domain->regions; int n = 0; for (ilist = 0; ilist < nregion; ilist++) regions[list[ilist]]->interior ^= 1; for (ilist = 0; ilist < nregion; ilist++) { iregion = list[ilist]; 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; zs = x[2] - regions[iregion]->contact[m].delz; for (jlist = 0; jlist < nregion; jlist++) { if (jlist == ilist) continue; jregion = list[jlist]; if (!regions[jregion]->match(xs,ys,zs)) break; } if (jlist == nregion) { contact[n].r = regions[iregion]->contact[m].r; contact[n].delx = regions[iregion]->contact[m].delx; contact[n].dely = regions[iregion]->contact[m].dely; contact[n].delz = regions[iregion]->contact[m].delz; n++; } } } for (ilist = 0; ilist < nregion; ilist++) regions[list[ilist]]->interior ^= 1; return n; } + +/* ---------------------------------------------------------------------- + change region shape of all sub-regions +------------------------------------------------------------------------- */ + +void RegUnion::shape_update() +{ + Region **regions = domain->regions; + for (int ilist = 0; ilist < nregion; ilist++) + regions[list[ilist]]->shape_update(); +} diff --git a/src/region_union.h b/src/region_union.h index 44d99375c..6d0ec063c 100644 --- a/src/region_union.h +++ b/src/region_union.h @@ -1,59 +1,61 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef REGION_CLASS RegionStyle(union,RegUnion) #else #ifndef LMP_REGION_UNION_H #define LMP_REGION_UNION_H #include "region.h" namespace LAMMPS_NS { class RegUnion : public Region { public: RegUnion(class LAMMPS *, int, char **); ~RegUnion(); + void init(); int dynamic_check(); int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); + void shape_update(); private: int nregion; int *list; }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Region union region ID does not exist One or more of the region IDs specified by the region union command does not exist. */