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

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);
     }
 

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;
 }
 
 /* ----------------------------------------------------------------------

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 dynamic;        // 1 if region position/orientation changes over time
-  int moveflag,rotateflag;
+  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 &);

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;
   }
 

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");
+}

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();
 };
 
 }

src/variable.cpp

@@ -2159,7 +2159,7 @@ double Variable::eval_tree(Tree *tree, int i)
   }
 
   if (tree->type == RMASK) {
-    if (domain->regions[tree->ivalue1]->inside(atom->x[i][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;
@@ -2167,7 +2167,7 @@ double Variable::eval_tree(Tree *tree, int i)
 
   if (tree->type == GRMASK) {
     if ((atom->mask[i] & tree->ivalue1) &&
-        (domain->regions[tree->ivalue2]->inside(atom->x[i][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;