From 2d6d2de507ff921827d4a23f168a819f0abe2a60 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 18 Jan 2013 20:39:12 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9299 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/fix_wall.cpp | 10 ++++++-- src/region.cpp | 34 ++++++++++++++++++++------- src/region.h | 17 ++++++++++---- src/region_intersect.cpp | 3 ++- src/region_sphere.cpp | 50 +++++++++++++++++++++++++++++++++++++--- src/region_sphere.h | 6 +++++ src/variable.cpp | 12 +++++----- 7 files changed, 108 insertions(+), 24 deletions(-) diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp index 70fcbbea2c..297f22741c 100644 --- a/src/fix_wall.cpp +++ b/src/fix_wall.cpp @@ -314,10 +314,16 @@ void FixWall::post_force(int vflag) else coord *= zscale; } else coord = coord0[m]; if (wstyle[m] == VARIABLE) { - if (estyle[m] == VARIABLE) + if (estyle[m] == VARIABLE) { epsilon[m] = input->variable->compute_equal(eindex[m]); - if (sstyle[m] == VARIABLE) + if (epsilon[m] < 0.0) + error->all(FLERR,"Variable evaluation in fix wall gave bad value"); + } + if (sstyle[m] == VARIABLE) { sigma[m] = input->variable->compute_equal(sindex[m]); + if (sigma[m] < 0.0) + error->all(FLERR,"Variable evaluation in fix wall gave bad value"); + } precompute(m); } diff --git a/src/region.cpp b/src/region.cpp index a9d27e1a38..14f7155e4a 100644 --- a/src/region.cpp +++ b/src/region.cpp @@ -36,9 +36,10 @@ Region::Region(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) style = new char[n]; strcpy(style,arg[1]); + varshape = 0; xstr = ystr = zstr = tstr = NULL; dx = dy = dz = 0.0; - laststep = -1; + lastshape = lastdynamic = -1; } /* ---------------------------------------------------------------------- */ @@ -100,23 +101,35 @@ int Region::dynamic_check() 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 has variable shape, invoke shape_update() once per timestep if region is dynamic, apply inverse transform to x,y,z unmove first, then unrotate, so don't have to change rotation point + caller is responsible for wrapping this call with + modify->clearstep_compute() and modify->addstep_compute() if needed ------------------------------------------------------------------------- */ int Region::match(double x, double y, double z) { - if (dynamic) inverse_transform(x,y,z); + if (varshape && update->ntimestep != lastshape) { + shape_update(); + lastshape = update->ntimestep; + } + + if (dynamic && update->ntimestep) inverse_transform(x,y,z); + return !(inside(x,y,z) ^ interior); } /* ---------------------------------------------------------------------- generate list of contact points for interior or exterior regions + if region has variable shape, invoke shape_update() once per timestep if region is dynamic: before: inverse transform x,y,z (unmove, then unrotate) after: forward transform contact point xs,yx,zs (rotate, then move), then reset contact delx,dely,delz based on new contact point no need to do this if no rotation since delxyz doesn't change + caller is responsible for wrapping this call with + modify->clearstep_compute() and modify->addstep_compute() if needed ------------------------------------------------------------------------- */ int Region::surface(double x, double y, double z, double cutoff) @@ -125,6 +138,11 @@ 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; @@ -179,13 +197,13 @@ 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 != laststep) + if (update->ntimestep != lastdynamic) theta = input->variable->compute_equal(tvar); rotate(x,y,z,theta); } if (moveflag) { - if (update->ntimestep != laststep) { + 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); @@ -195,7 +213,7 @@ void Region::forward_transform(double &x, double &y, double &z) z += dz; } - laststep = update->ntimestep; + lastdynamic = update->ntimestep; } /* ---------------------------------------------------------------------- @@ -206,7 +224,7 @@ 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 != laststep) { + 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); @@ -217,12 +235,12 @@ void Region::inverse_transform(double &x, double &y, double &z) } if (rotateflag) { - if (update->ntimestep != laststep) + if (update->ntimestep != lastdynamic) theta = input->variable->compute_equal(tvar); rotate(x,y,z,-theta); } - laststep = update->ntimestep; + lastdynamic = update->ntimestep; } /* ---------------------------------------------------------------------- diff --git a/src/region.h b/src/region.h index 2406dc186f..8f79590258 100644 --- a/src/region.h +++ b/src/region.h @@ -40,27 +40,36 @@ class Region : protected Pointers { Region(class LAMMPS *, int, char **); virtual ~Region(); - void init(); + 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 changes over time - int moveflag,rotateflag; + 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 laststep; + bigint lastshape,lastdynamic; void forward_transform(double &, double &, double &); void inverse_transform(double &, double &, double &); diff --git a/src/region_intersect.cpp b/src/region_intersect.cpp index 1d3d78d186..736ec3c6e1 100644 --- a/src/region_intersect.cpp +++ b/src/region_intersect.cpp @@ -37,7 +37,8 @@ RegIntersect::RegIntersect(LAMMPS *lmp, int narg, char **arg) : 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"); + if (iregion == -1) + error->all(FLERR,"Region intersect region ID does not exist"); list[nregion++] = iregion; } diff --git a/src/region_sphere.cpp b/src/region_sphere.cpp index e90704dd74..773fc4cb16 100644 --- a/src/region_sphere.cpp +++ b/src/region_sphere.cpp @@ -15,10 +15,15 @@ #include "stdlib.h" #include "string.h" #include "region_sphere.h" +#include "input.h" +#include "variable.h" +#include "update.h" #include "error.h" using namespace LAMMPS_NS; +enum{CONSTANT,VARIABLE}; + /* ---------------------------------------------------------------------- */ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : @@ -29,13 +34,23 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : xc = xscale*atof(arg[2]); yc = yscale*atof(arg[3]); zc = zscale*atof(arg[4]); - radius = xscale*atof(arg[5]); - // error check + rstr = NULL; + if (strstr(arg[5],"v_") == arg[5]) { + int n = strlen(&arg[5][2]) + 1; + rstr = new char[n]; + strcpy(rstr,&arg[5][2]); + radius = 0.0; + rstyle = VARIABLE; + varshape = 1; + } else { + radius = xscale*atof(arg[5]); + rstyle = CONSTANT; + } if (radius < 0.0) error->all(FLERR,"Illegal region sphere command"); - // extent of sphere + // extent of sphere, will be 0.0 for variable radius if (interior) { bboxflag = 1; @@ -55,9 +70,27 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : RegSphere::~RegSphere() { + delete [] rstr; delete [] contact; } +/* ---------------------------------------------------------------------- */ + +void RegSphere::init() +{ + Region::init(); + + // check variable + + if (rstr) { + rvar = input->variable->find(rstr); + if (rvar < 0) + error->all(FLERR,"Variable name for region sphere does not exist"); + if (!input->variable->equalstyle(rvar)) + error->all(FLERR,"Variable for region sphere is invalid style"); + } +} + /* ---------------------------------------------------------------------- inside = 1 if x,y,z is inside or on surface inside = 0 if x,y,z is outside and not on surface @@ -124,3 +157,14 @@ int RegSphere::surface_exterior(double *x, double cutoff) } return 0; } + +/* ---------------------------------------------------------------------- + change region shape via variable evaluation +------------------------------------------------------------------------- */ + +void RegSphere::shape_update() +{ + radius = xscale * input->variable->compute_equal(rvar); + if (radius < 0.0) + error->one(FLERR,"Variable evaluation in region gave bad value"); +} diff --git a/src/region_sphere.h b/src/region_sphere.h index 73556d07c3..ed34369221 100644 --- a/src/region_sphere.h +++ b/src/region_sphere.h @@ -28,13 +28,19 @@ class RegSphere : public Region { public: RegSphere(class LAMMPS *, int, char **); ~RegSphere(); + void init(); int inside(double, double, double); int surface_interior(double *, double); int surface_exterior(double *, double); + void shape_update(); private: double xc,yc,zc; double radius; + int rstyle,rvar; + char *rstr; + + void variable_check(); }; } diff --git a/src/variable.cpp b/src/variable.cpp index 5b9402a731..d60acccc58 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -2159,17 +2159,17 @@ double Variable::eval_tree(Tree *tree, int i) } if (tree->type == RMASK) { - if (domain->regions[tree->ivalue1]->inside(atom->x[i][0], - atom->x[i][1], - atom->x[i][2])) return 1.0; + if (domain->regions[tree->ivalue1]->match(atom->x[i][0], + atom->x[i][1], + atom->x[i][2])) return 1.0; else return 0.0; } if (tree->type == GRMASK) { if ((atom->mask[i] & tree->ivalue1) && - (domain->regions[tree->ivalue2]->inside(atom->x[i][0], - atom->x[i][1], - atom->x[i][2]))) return 1.0; + (domain->regions[tree->ivalue2]->match(atom->x[i][0], + atom->x[i][1], + atom->x[i][2]))) return 1.0; else return 0.0; }