support variables for more region properties

add code to allow the center of a spherical or cylindrical
region to be determined as variables and thus change over time.

doc/src/region.txt

@@ -26,7 +26,7 @@ style = {delete} or {block} or {cone} or {cylinder} or {plane} or {prism} or {sp
     dim = {x} or {y} or {z} = axis of cylinder
     c1,c2 = coords of cylinder axis in other 2 dimensions (distance units)
     radius = cylinder radius (distance units)
-      radius can be a variable (see below)
+      c1,c2, and radius can be a variable (see below)
     lo,hi = bounds of cylinder in dim (distance units)
   {plane} args = px py pz nx ny nz
     px,py,pz = point on the plane (distance units)
@@ -39,7 +39,7 @@ style = {delete} or {block} or {cone} or {cylinder} or {plane} or {prism} or {sp
   {sphere} args = x y z radius
     x,y,z = center of sphere (distance units)
     radius = radius of sphere (distance units)
-      radius can be a variable (see below)
+      x,y,z, and radius can be a variable (see below)
   {union} args = N reg-ID1 reg-ID2 ...
     N = # of regions to follow, must be 2 or greater
     reg-ID1,reg-ID2, ... = IDs of regions to join together
@@ -179,12 +179,17 @@ The {radius} value for style {sphere} and {cylinder} can be specified
 as an equal-style "variable"_variable.html.  If the value is a
 variable, it should be specified as v_name, where name is the variable
 name.  In this case, the variable will be evaluated each timestep, and
-its value used to determine the radius of the region.
+its value used to determine the radius of the region. For style {sphere}
+also the x-, y-, and z- coordinate of the center of the sphere and for
+style {cylinder} the two center positions c1 and c2 for the location of
+the cylinder axes can be a variable with the same kind of effect and
+requirements than for the radius.
 
 Equal-style variables can specify formulas with various mathematical
 functions, and include "thermo_style"_thermo_style.html command
 keywords for the simulation box parameters and timestep and elapsed
-time.  Thus it is easy to specify a time-dependent radius.
+time.  Thus it is easy to specify a time-dependent radius or have
+a time dependent position of the sphere or cylinder region.
 
 See the "Howto tricilinc"_Howto_triclinic.html doc page for a
 geometric description of triclinic boxes, as defined by LAMMPS, and

src/region_cylinder.cpp

@@ -30,7 +30,7 @@ enum{CONSTANT,VARIABLE};
 /* ---------------------------------------------------------------------- */
 
 RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
-  Region(lmp, narg, arg), rstr(NULL)
+  Region(lmp, narg, arg), c1str(NULL), c2str(NULL), rstr(NULL)
 {
   options(narg-8,&arg[8]);
 
@@ -44,17 +44,76 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
   axis = arg[2][0];
 
   if (axis == 'x') {
-    c1 = yscale*force->numeric(FLERR,arg[3]);
+    if (strstr(arg[3],"v_") == arg[3]) {
-    c2 = zscale*force->numeric(FLERR,arg[4]);
+      int n = strlen(arg[3]+2) + 1;
+      c1str = new char[n];
+      strcpy(c1str,arg[3]+2);
+      c1 = 0.0;
+      c1style = VARIABLE;
+      varshape = 1;
+    } else {
+      c1 = yscale*force->numeric(FLERR,arg[3]);
+      c1style = CONSTANT;
+    }
+    if (strstr(arg[4],"v_") == arg[4]) {
+      int n = strlen(arg[4]+2) + 1;
+      c2str = new char[n];
+      strcpy(c2str,arg[4]+2);
+      c2 = 0.0;
+      c2style = VARIABLE;
+      varshape = 1;
+    } else {
+      c2 = zscale*force->numeric(FLERR,arg[4]);
+      c2style = CONSTANT;
+    }
   } else if (axis == 'y') {
-    c1 = xscale*force->numeric(FLERR,arg[3]);
+    if (strstr(arg[3],"v_") == arg[3]) {
-    c2 = zscale*force->numeric(FLERR,arg[4]);
+      int n = strlen(arg[3]+2) + 1;
+      c1str = new char[n];
+      strcpy(c1str,arg[3]+2);
+      c1 = 0.0;
+      c1style = VARIABLE;
+      varshape = 1;
+    } else {
+      c1 = xscale*force->numeric(FLERR,arg[3]);
+      c1style = CONSTANT;
+    }
+    if (strstr(arg[4],"v_") == arg[4]) {
+      int n = strlen(arg[4]+2) + 1;
+      c2str = new char[n];
+      strcpy(c2str,arg[4]+2);
+      c2 = 0.0;
+      c2style = VARIABLE;
+      varshape = 1;
+    } else {
+      c2 = zscale*force->numeric(FLERR,arg[4]);
+      c2style = CONSTANT;
+    }
   } else if (axis == 'z') {
-    c1 = xscale*force->numeric(FLERR,arg[3]);
+    if (strstr(arg[3],"v_") == arg[3]) {
-    c2 = yscale*force->numeric(FLERR,arg[4]);
+      int n = strlen(arg[3]+2) + 1;
+      c1str = new char[n];
+      strcpy(c1str,arg[3]+2);
+      c1 = 0.0;
+      c1style = VARIABLE;
+      varshape = 1;
+    } else {
+      c1 = xscale*force->numeric(FLERR,arg[3]);
+      c1style = CONSTANT;
+    }
+    if (strstr(arg[4],"v_") == arg[4]) {
+      int n = strlen(arg[4]+2) + 1;
+      c2str = new char[n];
+      strcpy(c2str,arg[4]+2);
+      c2 = 0.0;
+      c2style = VARIABLE;
+      varshape = 1;
+    } else {
+      c2 = yscale*force->numeric(FLERR,arg[4]);
+      c2style = CONSTANT;
+    }
   }
 
-  rstr = NULL;
   if (strstr(arg[5],"v_") == arg[5]) {
     int n = strlen(&arg[5][2]) + 1;
     rstr = new char[n];
@@ -62,8 +121,6 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
     radius = 0.0;
     rstyle = VARIABLE;
     varshape = 1;
-    variable_check();
-    shape_update();
   } else {
     radius = force->numeric(FLERR,arg[5]);
     if (axis == 'x') radius *= yscale;
@@ -71,6 +128,11 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
     rstyle = CONSTANT;
   }
 
+  if (varshape) {
+    variable_check();
+    shape_update();
+  }
+
   if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) {
     if (domain->box_exist == 0)
       error->all(FLERR,"Cannot use region INF or EDGE when box does not exist");
@@ -167,6 +229,8 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
 
 RegCylinder::~RegCylinder()
 {
+  delete [] c1str;
+  delete [] c2str;
   delete [] rstr;
   delete [] contact;
 }
@@ -176,7 +240,7 @@ RegCylinder::~RegCylinder()
 void RegCylinder::init()
 {
   Region::init();
-  if (rstr) variable_check();
+  if (varshape) variable_check();
 }
 
 /* ----------------------------------------------------------------------
@@ -667,12 +731,27 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
 
 void RegCylinder::shape_update()
 {
-  radius = input->variable->compute_equal(rvar);
+  if (c1style == VARIABLE) c1 = input->variable->compute_equal(c1var);
-  if (radius < 0.0)
+  if (c2style == VARIABLE) c2 = input->variable->compute_equal(c2var);
-    error->one(FLERR,"Variable evaluation in region gave bad value");
+  if (rstyle == VARIABLE) {
-  if (axis == 'x') radius *= xscale;
+    radius = input->variable->compute_equal(rvar);
-  else if (axis == 'y') radius*= yscale;
+    if (radius < 0.0)
-  else radius *= zscale;
+      error->one(FLERR,"Variable evaluation in region gave bad value");
+  }
+
+  if (axis == 'x') {
+    if (c1style == VARIABLE) c1 *= yscale;
+    if (c2style == VARIABLE) c2 *= zscale;
+    if (rstyle == VARIABLE)  radius *= yscale;
+  } else if (axis == 'y') {
+    if (c1style == VARIABLE) c1 *= xscale;
+    if (c2style == VARIABLE) c2 *= zscale;
+    if (rstyle == VARIABLE)  radius *= xscale;
+  } else { // axis == 'z'
+    if (c1style == VARIABLE) c1 *= xscale;
+    if (c2style == VARIABLE) c2 *= yscale;
+    if (rstyle == VARIABLE)  radius *= xscale;
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -681,11 +760,29 @@ void RegCylinder::shape_update()
 
 void RegCylinder::variable_check()
 {
-  rvar = input->variable->find(rstr);
+  if (c1style == VARIABLE) {
-  if (rvar < 0)
+    c1var = input->variable->find(c1str);
-    error->all(FLERR,"Variable name for region cylinder does not exist");
+    if (c1var < 0)
-  if (!input->variable->equalstyle(rvar))
+      error->all(FLERR,"Variable name for region cylinder does not exist");
-    error->all(FLERR,"Variable for region cylinder is invalid style");
+    if (!input->variable->equalstyle(c1var))
+      error->all(FLERR,"Variable for region cylinder is invalid style");
+  }
+
+  if (c2style == VARIABLE) {
+    c2var = input->variable->find(c2str);
+    if (c2var < 0)
+      error->all(FLERR,"Variable name for region cylinder does not exist");
+    if (!input->variable->equalstyle(c2var))
+      error->all(FLERR,"Variable for region cylinder is invalid style");
+  }
+
+  if (rstyle == VARIABLE) {
+    rvar = input->variable->find(rstr);
+    if (rvar < 0)
+      error->all(FLERR,"Variable name for region cylinder does not exist");
+    if (!input->variable->equalstyle(rvar))
+      error->all(FLERR,"Variable for region cylinder is invalid style");
+  }
 }
 
 

src/region_cylinder.h

@@ -43,8 +43,10 @@ class RegCylinder : public Region {
   double c1,c2;
   double radius;
   double lo,hi;
+  int c1style,c1var;
+  int c2style,c2var;
   int rstyle,rvar;
-  char *rstr;
+  char *c1str,*c2str,*rstr;
 
   void variable_check();
 

src/region_sphere.cpp

@@ -28,15 +28,46 @@ enum{CONSTANT,VARIABLE};
 /* ---------------------------------------------------------------------- */
 
 RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
-  Region(lmp, narg, arg)
+  Region(lmp, narg, arg), xstr(NULL), ystr(NULL), zstr(NULL), rstr(NULL)
 {
   options(narg-6,&arg[6]);
 
-  xc = xscale*force->numeric(FLERR,arg[2]);
+  if (strstr(arg[2],"v_") == arg[2]) {
-  yc = yscale*force->numeric(FLERR,arg[3]);
+    int n = strlen(arg[2]+2) + 1;
-  zc = zscale*force->numeric(FLERR,arg[4]);
+    xstr = new char[n];
+    strcpy(xstr,arg[2]+2);
+    xc = 0.0;
+    xstyle = VARIABLE;
+    varshape = 1;
+  } else {
+    xc = xscale*force->numeric(FLERR,arg[2]);
+    xstyle = CONSTANT;
+  }
+
+  if (strstr(arg[3],"v_") == arg[3]) {
+    int n = strlen(arg[3]+2) + 1;
+    ystr = new char[n];
+    strcpy(ystr,arg[3]+2);
+    yc = 0.0;
+    ystyle = VARIABLE;
+    varshape = 1;
+  } else {
+    yc = yscale*force->numeric(FLERR,arg[3]);
+    ystyle = CONSTANT;
+  }
+
+  if (strstr(arg[4],"v_") == arg[4]) {
+    int n = strlen(arg[4]+2) + 1;
+    zstr = new char[n];
+    strcpy(zstr,arg[4]+2);
+    zc = 0.0;
+    zstyle = VARIABLE;
+    varshape = 1;
+  } else {
+    zc = zscale*force->numeric(FLERR,arg[4]);
+    zstyle = CONSTANT;
+  }
 
-  rstr = NULL;
   if (strstr(arg[5],"v_") == arg[5]) {
     int n = strlen(&arg[5][2]) + 1;
     rstr = new char[n];
@@ -44,19 +75,22 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
     radius = 0.0;
     rstyle = VARIABLE;
     varshape = 1;
-    variable_check();
-    shape_update();
   } else {
     radius = xscale*force->numeric(FLERR,arg[5]);
     rstyle = CONSTANT;
   }
 
+  if (varshape) {
+    variable_check();
+    shape_update();
+  }
+
   // error check
 
   if (radius < 0.0) error->all(FLERR,"Illegal region sphere command");
 
   // extent of sphere
-  // for variable radius, uses initial radius
+  // for variable radius, uses initial radius and origin for variable center
 
   if (interior) {
     bboxflag = 1;
@@ -77,6 +111,9 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
 
 RegSphere::~RegSphere()
 {
+  delete [] xstr;
+  delete [] ystr;
+  delete [] zstr;
   delete [] rstr;
   delete [] contact;
 }
@@ -86,7 +123,7 @@ RegSphere::~RegSphere()
 void RegSphere::init()
 {
   Region::init();
-  if (rstr) variable_check();
+  if (varshape) variable_check();
 }
 
 /* ----------------------------------------------------------------------
@@ -168,9 +205,20 @@ int RegSphere::surface_exterior(double *x, double cutoff)
 
 void RegSphere::shape_update()
 {
-  radius = xscale * input->variable->compute_equal(rvar);
+  if (xstyle == VARIABLE)
-  if (radius < 0.0)
+    xc = xscale * input->variable->compute_equal(xvar);
-    error->one(FLERR,"Variable evaluation in region gave bad value");
+
+  if (ystyle == VARIABLE)
+    yc = yscale * input->variable->compute_equal(yvar);
+
+  if (zstyle == VARIABLE)
+    zc = zscale * input->variable->compute_equal(zvar);
+
+  if (rstyle == VARIABLE) {
+    radius = xscale * input->variable->compute_equal(rvar);
+    if (radius < 0.0)
+      error->one(FLERR,"Variable evaluation in region gave bad value");
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -179,13 +227,38 @@ void RegSphere::shape_update()
 
 void RegSphere::variable_check()
 {
-  rvar = input->variable->find(rstr);
+  if (xstyle == VARIABLE) {
-  if (rvar < 0)
+    xvar = input->variable->find(xstr);
-    error->all(FLERR,"Variable name for region sphere does not exist");
+    if (xvar < 0)
-  if (!input->variable->equalstyle(rvar))
+      error->all(FLERR,"Variable name for region sphere does not exist");
-    error->all(FLERR,"Variable for region sphere is invalid style");
+    if (!input->variable->equalstyle(xvar))
-}
+      error->all(FLERR,"Variable for region sphere is invalid style");
+  }
 
+  if (ystyle == VARIABLE) {
+    yvar = input->variable->find(ystr);
+    if (yvar < 0)
+      error->all(FLERR,"Variable name for region sphere does not exist");
+    if (!input->variable->equalstyle(yvar))
+      error->all(FLERR,"Variable for region sphere is invalid style");
+  }
+
+  if (zstyle == VARIABLE) {
+    zvar = input->variable->find(zstr);
+    if (zvar < 0)
+      error->all(FLERR,"Variable name for region sphere does not exist");
+    if (!input->variable->equalstyle(zvar))
+      error->all(FLERR,"Variable for region sphere is invalid style");
+  }
+
+  if (rstyle == VARIABLE) {
+    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");
+  }
+}
 
 /* ----------------------------------------------------------------------
    Set values needed to calculate velocity due to shape changes.

src/region_sphere.h

@@ -40,8 +40,11 @@ class RegSphere : public Region {
  private:
   double xc,yc,zc;
   double radius;
+  int xstyle,xvar;
+  int ystyle,yvar;
+  int zstyle,zvar;
   int rstyle,rvar;
-  char *rstr;
+  char *xstr,*ystr,*zstr,*rstr;
 
   void variable_check();
 };