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

src/COLLOID/fix_wall_colloid.cpp

@@ -83,7 +83,10 @@ void FixWallColloid::precompute(int m)
   offset[m] = coeff3[m]*r4inv*r4inv*rinv - coeff4[m]*r2inv*rinv;
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   interaction of all particles in group with all 6 walls (if defined)
+   error if any finite-size particle is touching or penetrating wall
+------------------------------------------------------------------------- */
 
 void FixWallColloid::wall_particle(int m, double coord)
 {
@@ -103,13 +106,19 @@ void FixWallColloid::wall_particle(int m, double coord)
   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 <= 0.0) continue;
+      if (delta >= cutoff[m]) continue;
-      if (delta > cutoff[m]) continue;
       rad = shape[type[i]][0];
+      if (rad >= delta) {
+	onflag = 1;
+	continue;
+      }
+
       new_coeff2 = coeff2[m]*rad*rad*rad;
       diam = 2.0*rad;
       rad2 = rad*rad;
@@ -141,4 +150,6 @@ void FixWallColloid::wall_particle(int m, double coord)
 		   (-rinv2)*rinv3) - offset[m];
       ewall[m+1] += fwall;
     }
+
+  if (onflag) error->one("Particle on or inside fix wall surface");
 }

src/Makefile.package

@@ -2,8 +2,8 @@
 # this file is auto-edited when those packages are included/excluded
 
 PKG_INC =  -I../../lib/reax -I../../lib/poems -I../../lib/meam 
-PKG_PATH = -L../../lib/gpu -L../../lib/reax -L../../lib/poems -L../../lib/meam 
+PKG_PATH = -L../../lib/reax -L../../lib/poems -L../../lib/meam 
-PKG_LIB =  -lgpu -lreax -lpoems -lmeam 
+PKG_LIB =  -lreax -lpoems -lmeam 
 
-PKG_SYSPATH = $(gpu_SYSPATH) $(reax_SYSPATH) $(meam_SYSPATH) 
+PKG_SYSPATH = $(reax_SYSPATH) $(meam_SYSPATH) 
-PKG_SYSLIB =  $(gpu_SYSLIB) $(reax_SYSLIB) $(meam_SYSLIB) 
+PKG_SYSLIB =  $(reax_SYSLIB) $(meam_SYSLIB) 

src/create_box.cpp

@@ -46,8 +46,8 @@ void CreateBox::command(int narg, char **arg)
 
   int iregion = domain->find_region(arg[1]);
   if (iregion == -1) error->all("Create_box region ID does not exist");
-  if (domain->regions[iregion]->interior == 0)
+  if (domain->regions[iregion]->bboxflag == 0)
-    error->all("Create_box region must be of type inside");
+    error->all("Create_box region does not support a bounding box");
 
   // if region not prism:
   //   setup orthogonal domain

src/fix_deposit.cpp

@@ -68,8 +68,8 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
   // error check on region and its extent being inside simulation box
 
   if (iregion == -1) error->all("Must specify a region in fix deposit");
-  if (domain->regions[iregion]->interior == 0)
+  if (domain->regions[iregion]->bboxflag == 0)
-    error->all("Must use region with side = in with fix deposit");
+    error->all("Fix deposit region does not support a bounding box");
 
   xlo = domain->regions[iregion]->extent_xlo;
   xhi = domain->regions[iregion]->extent_xhi;

src/fix_wall.cpp

@@ -101,6 +101,10 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
   for (int m = 0; m < 6; m++) if (wallflag[m]) flag = 1;
   if (!flag) error->all("Illegal fix wall command");
 
+  for (int m = 0; m < 6; m++)
+    if (wallflag[m] && cutoff[m] <= 0.0)
+      error->all("Fix wall cutoff <= 0.0");
+
   if (velflag && wigflag) 
     error->all("Cannot set both vel and wiggle in fix wall command");
 

src/fix_wall_lj126.cpp

@@ -14,6 +14,7 @@
 #include "math.h"
 #include "fix_wall_lj126.h"
 #include "atom.h"
+#include "error.h"
 
 using namespace LAMMPS_NS;
 
@@ -36,7 +37,10 @@ void FixWallLJ126::precompute(int m)
   offset[m] = r6inv*(coeff3[m]*r6inv - coeff4[m]);
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   interaction of all particles in group with all 6 walls (if defined)
+   error if any particle is on or behind wall
+------------------------------------------------------------------------- */
 
 void FixWallLJ126::wall_particle(int m, double coord)
 {
@@ -51,12 +55,17 @@ void FixWallLJ126::wall_particle(int m, double coord)
   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 <= 0.0) continue;
+      if (delta >= cutoff[m]) continue;
-      if (delta > cutoff[m]) continue;
+      if (delta <= 0.0) {
+	onflag = 1;
+	continue;
+      }
       rinv = 1.0/delta;
       r2inv = rinv*rinv;
       r6inv = r2inv*r2inv*r2inv;
@@ -65,4 +74,6 @@ void FixWallLJ126::wall_particle(int m, double coord)
       ewall[0] += r6inv*(coeff3[m]*r6inv - coeff4[m]) - offset[m];
       ewall[m+1] += fwall;
     }
+
+  if (onflag) error->one("Particle on or inside fix wall surface");
 }

src/fix_wall_lj93.cpp

@@ -14,6 +14,7 @@
 #include "math.h"
 #include "fix_wall_lj93.h"
 #include "atom.h"
+#include "error.h"
 
 using namespace LAMMPS_NS;
 
@@ -37,7 +38,10 @@ void FixWallLJ93::precompute(int m)
   offset[m] = coeff3[m]*r4inv*r4inv*rinv - coeff4[m]*r2inv*rinv;
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   interaction of all particles in group with all 6 walls (if defined)
+   error if any particle is on or behind wall
+------------------------------------------------------------------------- */
 
 void FixWallLJ93::wall_particle(int m, double coord)
 {
@@ -52,12 +56,17 @@ void FixWallLJ93::wall_particle(int m, double coord)
   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 <= 0.0) continue;
+      if (delta >= cutoff[m]) continue;
-      if (delta > cutoff[m]) continue;
+      if (delta <= 0.0) {
+	onflag = 1;
+	continue;
+      }
       rinv = 1.0/delta;
       r2inv = rinv*rinv;
       r4inv = r2inv*r2inv;
@@ -68,4 +77,6 @@ void FixWallLJ93::wall_particle(int m, double coord)
 	coeff4[m]*r2inv*rinv - offset[m];
       ewall[m+1] += fwall;
     }
+
+  if (onflag) error->one("Particle on or inside fix wall surface");
 }

src/fix_wall_region.cpp

@@ -0,0 +1,357 @@
+/* ----------------------------------------------------------------------
+   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 "fix_wall_region.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
+#include "region.h"
+#include "lattice.h"
+#include "update.h"
+#include "output.h"
+#include "respa.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{LJ93,LJ126,COLLOID,HARMONIC};
+
+/* ---------------------------------------------------------------------- */
+
+FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg != 8) error->all("Illegal fix wall/region command");
+
+  scalar_flag = 1;
+  vector_flag = 1;
+  size_vector = 3;
+  global_freq = 1;
+  extscalar = 1;
+  extvector = 1;
+
+  // parse args
+
+  iregion = domain->find_region(arg[3]);
+  if (iregion == -1) error->all("Fix wall/region region ID does not exist");
+
+  if (strcmp(arg[4],"lj93") == 0) style = LJ93;
+  else if (strcmp(arg[4],"lj126") == 0) style = LJ126;
+  else if (strcmp(arg[4],"colloid") == 0) style = COLLOID;
+  else if (strcmp(arg[4],"harmonic") == 0) style = HARMONIC;
+  else error->all("Illegal fix wall/region command");
+
+  epsilon = atof(arg[5]);
+  sigma = atof(arg[6]);
+  cutoff = atof(arg[7]);
+
+  if (cutoff <= 0.0) error->all("Fix wall/region cutoff <= 0.0");
+
+  eflag = 0;
+  ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
+
+  // set this when regions have time dependence
+  // time_depend = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixWallRegion::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+  mask |= THERMO_ENERGY;
+  mask |= POST_FORCE_RESPA;
+  mask |= MIN_POST_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallRegion::init()
+{
+  // error checks for style COLLOID
+  // insure all particle shapes are spherical
+  // can be polydisperse
+  // insure all particles in group are extended particles
+
+  if (style == COLLOID) {
+    if (!atom->avec->shape_type)
+      error->all("Fix wall/region colloid requires atom attribute shape");
+    if (atom->radius_flag)
+      error->all("Fix wall/region colloid cannot be used with "
+		 "atom attribute diameter");
+
+    for (int i = 1; i <= atom->ntypes; i++)
+      if ((atom->shape[i][0] != atom->shape[i][1]) || 
+	  (atom->shape[i][0] != atom->shape[i][2]) ||
+	  (atom->shape[i][1] != atom->shape[i][2]))
+	error->all("Fix wall/region colloid requires spherical particles");
+
+    double **shape = atom->shape;
+    int *type = atom->type;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+    
+    int flag = 0;
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit)
+	if (atom->shape[type[i]][0] == 0.0) flag = 1;
+    
+    int flagall;
+    MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+    if (flagall) 
+      error->all("Fix wall/region colloid requires extended particles");
+  }
+
+  // setup coefficients for each style
+
+  if (style == LJ93) {
+    coeff1 = 6.0/5.0 * epsilon * pow(sigma,9.0);
+    coeff2 = 3.0 * epsilon * pow(sigma,3.0);
+    coeff3 = 2.0/15.0 * epsilon * pow(sigma,9.0);
+    coeff4 = epsilon * pow(sigma,3.0);
+    double rinv = 1.0/cutoff;
+    double r2inv = rinv*rinv;
+    double r4inv = r2inv*r2inv;
+    offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv;
+  } else if (style == LJ126) {
+    coeff1 = 48.0 * epsilon * pow(sigma,12.0);
+    coeff2 = 24.0 * epsilon * pow(sigma,6.0);
+    coeff3 = 4.0 * epsilon * pow(sigma,12.0);
+    coeff4 = 4.0 * epsilon * pow(sigma,6.0);
+    double r2inv = 1.0/(cutoff*cutoff);
+    double r6inv = r2inv*r2inv*r2inv;
+    offset = r6inv*(coeff3*r6inv - coeff4);
+  } else if (style == COLLOID) {
+    coeff1 = -4.0/315.0 * epsilon * pow(sigma,6.0);
+    coeff2 = -2.0/3.0 * epsilon;
+    coeff3 = epsilon * pow(sigma,6.0)/7560.0;
+    coeff4 = epsilon/6.0;
+    double rinv = 1.0/cutoff;
+    double r2inv = rinv*rinv;
+    double r4inv = r2inv*r2inv;
+    offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv;
+  }
+
+  if (strcmp(update->integrate_style,"respa") == 0)
+    nlevels_respa = ((Respa *) update->integrate)->nlevels;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallRegion::setup(int vflag)
+{
+  if (strcmp(update->integrate_style,"verlet") == 0)
+    post_force(vflag);
+  else {
+    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+    post_force_respa(vflag,nlevels_respa-1,0);
+    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallRegion::min_setup(int vflag)
+{
+  post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallRegion::post_force(int vflag)
+{
+  int i,m,n;
+  double fx,fy,fz,tooclose;
+
+  eflag = 0;
+  ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  double **shape = atom->shape;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  Region *region = domain->regions[iregion];
+  int onflag = 0;
+
+  // region->match() insures particle is in region or on surface, else error
+  // if returned contact dist r = 0, is on surface, also an error
+  // in COLLOID case, r <= shape radius is an error
+  
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (!region->match(x[i][0],x[i][1],x[i][2])) {
+	onflag = 1;
+	continue;
+      }
+      if (style == COLLOID) tooclose = shape[type[i]][0];
+      else tooclose = 0.0;
+
+      n = region->surface(x[i],cutoff);
+
+      for (m = 0; m < n; m++) {
+	if (region->contact[m].r <= tooclose) {
+	  onflag = 1;
+	  continue;
+	}
+
+	if (style == LJ93) lj93(region->contact[m].r);
+	else if (style == LJ126) lj126(region->contact[m].r);
+	else if (style == COLLOID) 
+	  colloid(region->contact[m].r,shape[type[i]][0]);
+	else harmonic(region->contact[m].r);
+
+	ewall[0] += eng;
+	fx = fwall * region->contact[m].delx;
+	fy = fwall * region->contact[m].dely;
+	fz = fwall * region->contact[m].delz;
+	f[i][0] += fx;
+	f[i][1] += fy;
+	f[i][2] += fz;
+	ewall[1] -= fx;
+	ewall[2] -= fy;
+	ewall[3] -= fz;
+      }
+    }
+
+  if (onflag) error->one("Particle on or inside fix wall/region surface");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallRegion::post_force_respa(int vflag, int ilevel, int iloop)
+{
+  if (ilevel == nlevels_respa-1) post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallRegion::min_post_force(int vflag)
+{
+  post_force(vflag);
+}
+
+/* ----------------------------------------------------------------------
+   energy of wall interaction
+------------------------------------------------------------------------- */
+
+double FixWallRegion::compute_scalar()
+{
+  // only sum across procs one time
+
+  if (eflag == 0) {
+    MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world);
+    eflag = 1;
+  }
+  return ewall_all[0];
+}
+
+/* ----------------------------------------------------------------------
+   components of force on wall
+------------------------------------------------------------------------- */
+
+double FixWallRegion::compute_vector(int n)
+{
+  // only sum across procs one time
+
+  if (eflag == 0) {
+    MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world);
+    eflag = 1;
+  }
+  return ewall_all[n+1];
+}
+
+/* ----------------------------------------------------------------------
+   LJ 9/3 interaction for particle with wall
+   compute eng and fwall = magnitude of wall force
+------------------------------------------------------------------------- */
+
+void FixWallRegion::lj93(double r)
+{
+  double rinv = 1.0/r;
+  double r2inv = rinv*rinv;
+  double r4inv = r2inv*r2inv;
+  double r10inv = r4inv*r4inv*r2inv;
+  fwall = coeff1*r10inv - coeff2*r4inv;
+  eng = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv - offset;
+}
+
+/* ----------------------------------------------------------------------
+   LJ 12/6 interaction for particle with wall
+   compute eng and fwall = magnitude of wall force
+------------------------------------------------------------------------- */
+
+void FixWallRegion::lj126(double r)
+{
+  double rinv = 1.0/r;
+  double r2inv = rinv*rinv;
+  double r6inv = r2inv*r2inv*r2inv;
+  fwall = r6inv*(coeff1*r6inv - coeff2) * rinv;
+  eng = r6inv*(coeff3*r6inv - coeff4) - offset;
+}
+
+/* ----------------------------------------------------------------------
+   colloid interaction for finite-size particle of rad with wall
+   compute eng and fwall = magnitude of wall force
+------------------------------------------------------------------------- */
+
+void FixWallRegion::colloid(double r, double rad)
+{
+  double new_coeff2 = coeff2*rad*rad*rad;
+  double diam = 2.0*rad;
+
+  double rad2 = rad*rad;
+  double rad4 = rad2*rad2;
+  double rad8 = rad4*rad4;
+  double delta2 = rad2 - r*r;
+  double rinv = 1.0/delta2;
+  double r2inv = rinv*rinv;
+  double r4inv = r2inv*r2inv;
+  double r8inv = r4inv*r4inv;
+  fwall = coeff1*(rad8*rad + 27.0*rad4*rad2*rad*pow(r,2.0)
+		  + 63.0*rad4*rad*pow(r,4.0)
+		  + 21.0*rad2*rad*pow(r,6.0))*r8inv - new_coeff2*r2inv;
+
+  double r2 = 0.5*diam - r;
+  double rinv2 = 1.0/r2;
+  double r2inv2 = rinv2*rinv2;
+  double r4inv2 = r2inv2*r2inv2;
+  double r6inv2 = r4inv2*r2inv2;
+  double r3 = r + 0.5*diam;
+  double rinv3 = 1.0/r3;
+  double r2inv3 = rinv3*rinv3;
+  double r4inv3 = r2inv3*r2inv3;
+  double r6inv3 = r4inv3*r2inv3;
+  eng = coeff3*((-3.5*diam+r)*r4inv2*r2inv2*rinv2
+		+ (3.5*diam+r)*r4inv3*r2inv3*rinv3) -
+    coeff4*((-diam*r+r2*r3*(log(-r2)-log(r3)))*
+	    (-rinv2)*rinv3) - offset;
+}
+
+/* ----------------------------------------------------------------------
+   harmonic interaction for particle with wall
+   compute eng and fwall = magnitude of wall force
+------------------------------------------------------------------------- */
+
+void FixWallRegion::harmonic(double r)
+{
+  fwall = 2.0*epsilon*r;
+  eng = epsilon*r*r - offset;
+}

src/fix_wall_region.h

@@ -0,0 +1,54 @@
+/* ----------------------------------------------------------------------
+   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_REGION_H
+#define FIX_WALL_REGION_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixWallRegion : public Fix {
+ public:
+  FixWallRegion(class LAMMPS *, int, char **);
+  ~FixWallRegion() {}
+  int setmask();
+  void init();
+  void setup(int);
+  void min_setup(int);
+  void post_force(int);
+  void post_force_respa(int, int, int);
+  void min_post_force(int);
+  double compute_scalar();
+  double compute_vector(int);
+
+ private:
+  int style,iregion;
+  double epsilon,sigma,cutoff;
+  int eflag;
+  double ewall[4],ewall_all[4];
+  int nlevels_respa;
+  double dt;
+
+  double coeff1,coeff2,coeff3,coeff4,offset;
+  double eng,fwall;
+
+  void lj93(double);
+  void lj126(double);
+  void colloid(double, double);
+  void harmonic(double);
+};
+
+}
+
+#endif

src/region.cpp

@@ -11,6 +11,7 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+#include "math.h"
 #include "stdlib.h"
 #include "string.h"
 #include "region.h"
@@ -83,3 +84,30 @@ void Region::options(int narg, char **arg)
   }
   else xscale = yscale = zscale = 1.0;
 }
+
+/* ----------------------------------------------------------------------
+   generate list of contact points for interior or exterior regions
+------------------------------------------------------------------------- */
+
+int Region::surface(double *x, double cutoff)
+{
+  if (interior) return surface_interior(x,cutoff);
+  return surface_exterior(x,cutoff);
+}
+
+/* ----------------------------------------------------------------------
+   add a single contact at Nth location in contact array
+   x = particle position
+   xp,yp,zp = region surface point
+------------------------------------------------------------------------- */
+
+void Region::add_contact(int n, double *x, double xp, double yp, double zp)
+{
+  double delx = x[0] - xp;
+  double dely = x[1] - yp;
+  double delz = x[2] - zp;
+  contact[n].r = sqrt(delx*delx + dely*dely + delz*delz);
+  contact[n].delx = delx;
+  contact[n].dely = dely;
+  contact[n].delz = delz;
+}

src/region.h

@@ -27,13 +27,27 @@ class Region : protected Pointers {
   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
+
+  // 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 int match(double, double, double) = 0;
+  int surface(double *, double);
+  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);
 };
 
 }

src/region_block.cpp

@@ -21,6 +21,9 @@ using namespace LAMMPS_NS;
 
 #define BIG 1.0e20
 
+#define MIN(A,B) ((A) < (B)) ? (A) : (B)
+#define MAX(A,B) ((A) > (B)) ? (A) : (B)
+
 /* ---------------------------------------------------------------------- */
 
 RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
@@ -82,22 +85,147 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
 
   // extent of block
 
+  if (interior) {
+    bboxflag = 1;
     extent_xlo = xlo;
     extent_xhi = xhi;
     extent_ylo = ylo;
     extent_yhi = yhi;
     extent_zlo = zlo;
     extent_zhi = zhi;
+  } else bboxflag = 0;
+
+  // particle could be close to all 6 planes
+
+  cmax = 6;
+  contact = new Contact[cmax];
 }
 
 /* ---------------------------------------------------------------------- */
 
+RegBlock::~RegBlock()
+{
+  delete [] contact;
+}
+
+/* ----------------------------------------------------------------------
+   inside = 1 if x,y,z is inside or on surface
+   inside = 0 if x,y,z is outside and not on surface
+------------------------------------------------------------------------- */
+
 int RegBlock::match(double x, double y, double z)
 {
-  int inside;
+  int inside = 0;
   if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi)
     inside = 1;
-  else inside = 0;
 
   return !(inside ^ interior);         // 1 if same, 0 if different
 }
+
+/* ----------------------------------------------------------------------
+   contact if 0 <= x < cutoff from one or more inner surfaces of block
+   can be one contact for each of 6 faces
+   no contact if outside (possible if called from union/intersect)
+   delxyz = vector from nearest point on block to x
+------------------------------------------------------------------------- */
+
+int RegBlock::surface_interior(double *x, double cutoff)
+{
+  double delta;
+
+  // x is exterior to block
+
+  if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > zhi ||
+      x[2] < ylo || x[2] > zhi) return 0;
+
+  // x is interior to block or on its surface
+
+  int n = 0;
+
+  delta = x[0] - xlo;
+  if (delta < cutoff) {
+    contact[n].r = delta;
+    contact[n].delx = delta;
+    contact[n].dely = contact[n].delz = 0.0;
+    n++;
+  }
+  delta = xhi - x[0];
+  if (delta < cutoff) {
+    contact[n].r = delta;
+    contact[n].delx = -delta;
+    contact[n].dely = contact[n].delz = 0.0;
+    n++;
+  }
+
+  delta = x[1] - ylo;
+  if (delta < cutoff) {
+    contact[n].r = delta;
+    contact[n].dely = delta;
+    contact[n].delx = contact[n].delz = 0.0;
+    n++;
+  }
+  delta = yhi - x[1];
+  if (delta < cutoff) {
+    contact[n].r = delta;
+    contact[n].dely = -delta;
+    contact[n].delx = contact[n].delz = 0.0;
+    n++;
+  }
+
+  delta = x[2] - zlo;
+  if (delta < cutoff) {
+    contact[n].r = delta;
+    contact[n].delz = delta;
+    contact[n].delx = contact[n].dely = 0.0;
+    n++;
+  }
+  delta = zhi - x[2];
+  if (delta < cutoff) {
+    contact[n].r = delta;
+    contact[n].delz = -delta;
+    contact[n].delx = contact[n].dely = 0.0;
+    n++;
+  }
+
+  return n;
+}
+
+/* ----------------------------------------------------------------------
+   one contact if 0 <= x < cutoff from outer surface of block
+   no contact if inside (possible if called from union/intersect)
+   delxyz = vector from nearest point on block to x
+------------------------------------------------------------------------- */
+
+int RegBlock::surface_exterior(double *x, double cutoff)
+{
+  double xp,yp,zp;
+  double delta;
+
+  // x is far enough from block that there is no contact
+  // x is interior to block
+
+  if (x[0] <= xlo-cutoff || x[0] >= xhi+cutoff ||
+      x[1] <= ylo-cutoff || x[1] >= yhi+cutoff ||
+      x[2] <= zlo-cutoff || x[2] >= zhi+cutoff) return 0;
+  if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi &&
+      x[2] > zlo && x[2] < zhi) return 0;
+
+  // x is exterior to block or on its surface
+  // xp,yp,zp = point on surface of block that x is closest to
+  //            could be edge or corner pt of block
+  // do not add contact point if r >= cutoff
+
+  if (x[0] < xlo) xp = xlo;
+  else if (x[0] > xhi) xp = xhi;
+  else xp = x[0];
+  if (x[1] < ylo) yp = ylo;
+  else if (x[1] > yhi) yp = yhi;
+  else yp = x[1];
+  if (x[2] < zlo) zp = zlo;
+  else if (x[2] > zhi) zp = zhi;
+  else zp = x[2];
+
+  add_contact(0,x,xp,yp,zp);
+  if (contact[0].r < cutoff) return 1;
+  return 0;
+}

src/region_block.h

@@ -23,7 +23,10 @@ class RegBlock : public Region {
 
  public:
   RegBlock(class LAMMPS *, int, char **);
+  ~RegBlock();
   int match(double, double, double);
+  int surface_interior(double *, double);
+  int surface_exterior(double *, double);
 
  private:
   double xlo,xhi,ylo,yhi,zlo,zhi;

src/region_cone.cpp

@@ -40,15 +40,19 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) :
   if (axis == 'x') {
     c1 = yscale*atof(arg[3]);
     c2 = zscale*atof(arg[4]);
+    radiuslo = yscale*atof(arg[5]);
+    radiushi = yscale*atof(arg[6]);
   } else if (axis == 'y') {
     c1 = xscale*atof(arg[3]);
     c2 = zscale*atof(arg[4]);
+    radiuslo = xscale*atof(arg[5]);
+    radiushi = xscale*atof(arg[6]);
   } else if (axis == 'z') {
     c1 = xscale*atof(arg[3]);
     c2 = yscale*atof(arg[4]);
-  }
     radiuslo = xscale*atof(arg[5]);
     radiushi = xscale*atof(arg[6]);
+  }
 
   if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) {
     if (domain->box_exist == 0) 
@@ -102,47 +106,67 @@ RegCone::RegCone(LAMMPS *lmp, int narg, char **arg) :
 
   if (radiuslo < 0.0) error->all("Illegal radius in region cone command");
   if (radiushi < 0.0) error->all("Illegal radius in region cone command");
+  if (radiuslo == 0.0 && radiushi == 0.0)
+    error->all("Illegal radius in region cone command");
   if (hi == lo) error->all("Illegal cone length in region cone command");
 
   // extent of cone
 
-  double largestradius;
+  if (radiuslo > radiushi) maxradius = radiuslo;
-  if (radiuslo > radiushi) largestradius = radiuslo;
+  else maxradius = radiushi;
-  else largestradius = radiushi;
+
+  if (interior) {
+    bboxflag = 1;
 
     if (axis == 'x') {
       extent_xlo = lo;
       extent_xhi = hi;
-    extent_ylo = c1 - largestradius;
+      extent_ylo = c1 - maxradius;
-    extent_yhi = c1 + largestradius;
+      extent_yhi = c1 + maxradius;
-    extent_zlo = c2 - largestradius;
+      extent_zlo = c2 - maxradius;
-    extent_zhi = c2 + largestradius;
+      extent_zhi = c2 + maxradius;
     }
     if (axis == 'y') {
-    extent_xlo = c1 - largestradius;
+      extent_xlo = c1 - maxradius;
-    extent_xhi = c1 + largestradius;
+      extent_xhi = c1 + maxradius;
       extent_ylo = lo;
       extent_yhi = hi;
-    extent_zlo = c2 - largestradius;
+      extent_zlo = c2 - maxradius;
-    extent_zhi = c2 + largestradius;
+      extent_zhi = c2 + maxradius;
     }
     if (axis == 'z') {
-    extent_xlo = c1 - largestradius;
+      extent_xlo = c1 - maxradius;
-    extent_xhi = c1 + largestradius;
+      extent_xhi = c1 + maxradius;
-    extent_ylo = c2 - largestradius;
+      extent_ylo = c2 - maxradius;
-    extent_yhi = c2 + largestradius;
+      extent_yhi = c2 + maxradius;
       extent_zlo = lo;
       extent_zhi = hi;
     }
+  } else bboxflag = 0;
+
+  // particle could be contact cone surface and 2 ends
+
+  cmax = 3;
+  contact = new Contact[cmax];
 }
 
 /* ---------------------------------------------------------------------- */
 
+RegCone::~RegCone()
+{
+  delete [] contact;
+}
+
+/* ----------------------------------------------------------------------
+   inside = 1 if x,y,z is inside or on surface
+   inside = 0 if x,y,z is outside and not on surface
+------------------------------------------------------------------------- */
+
 int RegCone::match(double x, double y, double z)
 {
   double del1,del2,dist;
-  int inside;
   double currentradius;
+  int inside;
 
   if (axis == 'x') {
     del1 = y - c1;
@@ -171,3 +195,389 @@ int RegCone::match(double x, double y, double z)
 
   return !(inside ^ interior);         // 1 if same, 0 if different
 }
+
+/* ----------------------------------------------------------------------
+   contact if 0 <= x < cutoff from one or more inner surfaces of cone
+   can be one contact for each of 3 cone surfaces
+   no contact if outside (possible if called from union/intersect)
+   delxyz = vector from nearest point on cone to x
+   special case: no contact with curved surf if x is on center axis
+------------------------------------------------------------------------- */
+
+int RegCone::surface_interior(double *x, double cutoff)
+{
+  double del1,del2,r,currentradius,delx,dely,delz,dist,delta;
+  double surflo[3],surfhi[3],xs[3];
+
+  int n = 0;
+
+  if (axis == 'x') {
+    del1 = x[1] - c1;
+    del2 = x[2] - c2;
+    r = sqrt(del1*del1 + del2*del2);
+    currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo);
+
+    // x is exterior to cone
+    
+    if (r > currentradius || x[0] < lo || x[0] > hi) return 0;
+
+    // x is interior to cone or on its surface
+    // surflo = pt on outer circle of bottom end plane, same dir as x vs axis
+    // surfhi = pt on outer circle of top end plane, same dir as x vs axis
+
+    if (r > 0.0) {
+      surflo[0] = lo;
+      surflo[1] = c1 + del1*radiuslo/r;
+      surflo[2] = c2 + del2*radiuslo/r;
+      surfhi[0] = hi;
+      surfhi[1] = c1 + del1*radiushi/r;
+      surfhi[2] = c2 + del2*radiushi/r;
+      point_on_line_segment(surflo,surfhi,x,xs);
+      delx = x[0] - xs[0];
+      dely = x[1] - xs[1];
+      delz = x[2] - xs[2];
+      dist = sqrt(delx*delx + dely*dely + delz*delz);
+      if (dist < cutoff) {
+	contact[n].r = dist;
+	contact[n].delx = delx;
+	contact[n].dely = dely;
+	contact[n].delz = delz;
+	n++;
+      }
+    }
+
+    delta = x[0] - lo;
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].delx = delta;
+      contact[n].dely = contact[n].delz = 0.0;
+      n++;
+    }
+    delta = hi - x[0];
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].delx = -delta;
+      contact[n].dely = contact[n].delz = 0.0;
+      n++;
+    }
+
+  } else if (axis == 'y') {
+    delx = x[0] - c1;
+    delz = x[2] - c2;
+    r = sqrt(delx*delx + delz*delz);
+    currentradius = radiuslo + (x[1]-lo)*(radiushi-radiuslo)/(hi-lo);
+
+    // x is exterior to cone
+    
+    if (r > currentradius || x[1] < lo || x[1] > hi) return 0;
+
+    // x is interior to cone or on its surface
+    // surflo = pt on outer circle of bottom end plane, same dir as x vs axis
+    // surfhi = pt on outer circle of top end plane, same dir as x vs axis
+
+    if (r > 0.0) {
+      surflo[0] = c1 + del1*radiuslo/r;
+      surflo[1] = lo;
+      surflo[2] = c2 + del2*radiuslo/r;
+      surfhi[0] = c1 + del1*radiushi/r;
+      surfhi[1] = hi;
+      surfhi[2] = c2 + del2*radiushi/r;
+      point_on_line_segment(surflo,surfhi,x,xs);
+      delx = x[0] - xs[0];
+      dely = x[1] - xs[1];
+      delz = x[2] - xs[2];
+      dist = sqrt(delx*delx + dely*dely + delz*delz);
+      if (dist < cutoff) {
+	contact[n].r = dist;
+	contact[n].delx = delx;
+	contact[n].dely = dely;
+	contact[n].delz = delz;
+	n++;
+      }
+    }
+
+    delta = x[1] - lo;
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].delz = delta;
+      contact[n].delx = contact[n].dely = 0.0;
+      n++;
+    }
+    delta = hi - x[1];
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].delz = -delta;
+      contact[n].delx = contact[n].dely = 0.0;
+      n++;
+    }
+
+  } else {
+    delx = x[0] - c1;
+    dely = x[1] - c2;
+    r = sqrt(delx*delx + dely*dely);
+    currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo);
+
+    // x is exterior to cone
+    
+    if (r > currentradius || x[2] < lo || x[2] > hi) return 0;
+
+    // x is interior to cone or on its surface
+    // surflo = pt on outer circle of bottom end plane, same dir as x vs axis
+    // surfhi = pt on outer circle of top end plane, same dir as x vs axis
+
+    if (r > 0.0) {
+      surflo[0] = c1 + del1*radiuslo/r;
+      surflo[1] = c2 + del2*radiuslo/r;
+      surflo[2] = lo;
+      surfhi[0] = c1 + del1*radiushi/r;
+      surfhi[1] = c2 + del2*radiushi/r;
+      surfhi[2] = hi;
+      point_on_line_segment(surflo,surfhi,x,xs);
+      delx = x[0] - xs[0];
+      dely = x[1] - xs[1];
+      delz = x[2] - xs[2];
+      dist = sqrt(delx*delx + dely*dely + delz*delz);
+      if (dist < cutoff) {
+	contact[n].r = dist;
+	contact[n].delx = delx;
+	contact[n].dely = dely;
+	contact[n].delz = delz;
+	n++;
+      }
+    }
+
+    delta = x[2] - lo;
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].delz = delta;
+      contact[n].delx = contact[n].dely = 0.0;
+      n++;
+    }
+    delta = hi - x[2];
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].delz = -delta;
+      contact[n].delx = contact[n].dely = 0.0;
+      n++;
+    }
+  }
+
+  return n;
+}
+
+/* ----------------------------------------------------------------------
+   one contact if 0 <= x < cutoff from outer surface of cone
+   no contact if inside (possible if called from union/intersect)
+   delxyz = vector from nearest point on cone to x
+------------------------------------------------------------------------- */
+
+int RegCone::surface_exterior(double *x, double cutoff)
+{
+  double del1,del2,r,currentradius,distsq;
+  double corner1[3],corner2[3],corner3[3],corner4[3],xp[3],nearest[3];
+
+  if (axis == 'x') {
+    del1 = x[1] - c1;
+    del2 = x[2] - c2;
+    r = sqrt(del1*del1 + del2*del2);
+    currentradius = radiuslo + (x[0]-lo)*(radiushi-radiuslo)/(hi-lo);
+
+    // x is far enough from cone that there is no contact
+    // x is interior to cone
+
+    if (r >= maxradius+cutoff || 
+	x[0] <= lo-cutoff || x[0] >= hi+cutoff) return 0;
+    if (r < currentradius && x[0] > lo && x[0] < hi) return 0;
+
+    // x is exterior to cone or on its surface
+    // corner1234 = 4 corner pts of half trapezoid = cone surf in plane of x
+    // project x to 3 line segments in half trapezoid (4th is axis of cone)
+    // nearest = point on surface of cone that x is closest to
+    //           could be edge of cone
+    // do not add contact point if r >= cutoff
+
+    corner1[0] = lo;
+    corner1[1] = c1 + del1*radiuslo/r;
+    corner1[2] = c2 + del2*radiuslo/r;
+    corner2[0] = hi;
+    corner2[1] = c1 + del1*radiushi/r;
+    corner2[2] = c2 + del2*radiushi/r;
+    corner3[0] = lo;
+    corner3[1] = c1;
+    corner3[2] = c2;
+    corner4[0] = hi;
+    corner4[1] = c1;
+    corner4[2] = c2;
+
+    distsq = BIG;
+    point_on_line_segment(corner1,corner2,x,xp);
+    distsq = closest(x,xp,nearest,distsq);
+    point_on_line_segment(corner1,corner3,x,xp);
+    distsq = closest(x,xp,nearest,distsq);
+    point_on_line_segment(corner2,corner4,x,xp);
+    distsq = closest(x,xp,nearest,distsq);
+
+    add_contact(0,x,nearest[0],nearest[1],nearest[2]);
+    if (contact[0].r < cutoff) return 1;
+    return 0;
+
+  } else if (axis == 'y') {
+    del1 = x[0] - c1;
+    del2 = x[2] - c2;
+    r = sqrt(del1*del1 + del2*del2);
+    currentradius = radiuslo + (x[1]-lo)*(radiushi-radiuslo)/(hi-lo);
+
+    // x is far enough from cone that there is no contact
+    // x is interior to cone
+
+    if (r >= maxradius+cutoff || 
+	x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0;
+    if (r < currentradius && x[1] > lo && x[1] < hi) return 0;
+
+    // x is exterior to cone or on its surface
+    // corner1234 = 4 corner pts of half trapezoid = cone surf in plane of x
+    // project x to 3 line segments in half trapezoid (4th is axis of cone)
+    // nearest = point on surface of cone that x is closest to
+    //           could be edge of cone
+    // do not add contact point if r >= cutoff
+
+    corner1[0] = c1 + del1*radiuslo/r;
+    corner1[1] = lo;
+    corner1[2] = c2 + del2*radiuslo/r;
+    corner2[0] = c1 + del1*radiushi/r;
+    corner2[1] = hi;
+    corner2[2] = c2 + del2*radiushi/r;
+    corner3[0] = c1;
+    corner3[1] = lo;
+    corner3[2] = c2;
+    corner4[0] = c1;
+    corner4[1] = hi;
+    corner4[2] = c2;
+
+    distsq = BIG;
+    point_on_line_segment(corner1,corner2,x,xp);
+    distsq = closest(x,xp,nearest,distsq);
+    point_on_line_segment(corner1,corner3,x,xp);
+    distsq = closest(x,xp,nearest,distsq);
+    point_on_line_segment(corner2,corner4,x,xp);
+    distsq = closest(x,xp,nearest,distsq);
+
+    add_contact(0,x,nearest[0],nearest[1],nearest[2]);
+    if (contact[0].r < cutoff) return 1;
+    return 0;
+
+  } else {
+    del1 = x[0] - c1;
+    del2 = x[1] - c2;
+    r = sqrt(del1*del1 + del2*del2);
+    currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo);
+
+    // x is far enough from cone that there is no contact
+    // x is interior to cone
+
+    if (r >= maxradius+cutoff || 
+	x[2] <= lo-cutoff || x[2] >= hi+cutoff) return 0;
+    if (r < currentradius && x[2] > lo && x[2] < hi) return 0;
+
+    // x is exterior to cone or on its surface
+    // corner1234 = 4 corner pts of half trapezoid = cone surf in plane of x
+    // project x to 3 line segments in half trapezoid (4th is axis of cone)
+    // nearest = point on surface of cone that x is closest to
+    //           could be edge of cone
+    // do not add contact point if r >= cutoff
+
+    corner1[0] = c1 + del1*radiuslo/r;
+    corner1[1] = c2 + del2*radiuslo/r;
+    corner1[2] = lo;
+    corner2[0] = c1 + del1*radiushi/r;
+    corner2[1] = c2 + del2*radiushi/r;
+    corner2[2] = hi;
+    corner3[0] = c1;
+    corner3[1] = c2;
+    corner3[2] = lo;
+    corner4[0] = c1;
+    corner4[1] = c2;
+    corner4[2] = hi;
+
+    distsq = BIG;
+    point_on_line_segment(corner1,corner2,x,xp);
+    distsq = closest(x,xp,nearest,distsq);
+    point_on_line_segment(corner1,corner3,x,xp);
+    distsq = closest(x,xp,nearest,distsq);
+    point_on_line_segment(corner2,corner4,x,xp);
+    distsq = closest(x,xp,nearest,distsq);
+
+    add_contact(0,x,nearest[0],nearest[1],nearest[2]);
+    if (contact[0].r < cutoff) return 1;
+    return 0;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   find nearest point to C on line segment A,B and return it as D
+   project (C-A) onto (B-A)
+   t = length of that projection, normalized by length of (B-A)
+   t <= 0, C is closest to A
+   t >= 1, C is closest to B
+   else closest point is between A and B
+------------------------------------------------------------------------- */
+
+void RegCone::point_on_line_segment(double *a, double *b, 
+				    double *c, double *d)
+{
+  double ba[3],ca[3];
+
+  subtract(a,b,ba);
+  subtract(a,c,ca);
+  double t = dotproduct(ca,ba) / dotproduct(ba,ba);
+
+  if (t <= 0.0) {
+    d[0] = a[0];
+    d[1] = a[1];
+    d[2] = a[2];
+  } else if (t >= 1.0) {
+    d[0] = b[0];
+    d[1] = b[1];
+    d[2] = b[2];
+  } else {
+    d[0] = a[0] + t*ba[0];
+    d[1] = a[1] + t*ba[1];
+    d[2] = a[2] + t*ba[2];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double RegCone::closest(double *x, double *near, double *nearest, double dsq)
+{
+  double delx = x[0] - near[0];
+  double dely = x[1] - near[1];
+  double delz = x[2] - near[2];
+  double rsq = delx*delx + dely*dely + delz*delz;
+  if (rsq >= dsq) return dsq;
+
+  nearest[0] = near[0];
+  nearest[1] = near[1];
+  nearest[2] = near[2];
+  return rsq;
+}
+
+/* ----------------------------------------------------------------------
+   v3 = v2 - v1
+------------------------------------------------------------------------- */
+
+void RegCone::subtract(double *v1, double *v2, double *v3)
+{
+  v3[0] = v2[0] - v1[0];
+  v3[1] = v2[1] - v1[1];
+  v3[2] = v2[2] - v1[2];
+}
+
+/* ----------------------------------------------------------------------
+   return dotproduct = v1 dot v2
+------------------------------------------------------------------------- */
+
+double RegCone::dotproduct(double *v1, double *v2)
+{
+  return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]);
+}

src/region_cone.h

@@ -21,13 +21,23 @@ namespace LAMMPS_NS {
 class RegCone : public Region {
  public:
   RegCone(class LAMMPS *, int, char **);
+  ~RegCone();
   int match(double, double, double);
+  int surface_interior(double *, double);
+  int surface_exterior(double *, double);
 
  private:
   char axis;
   double c1,c2;
   double radiuslo,radiushi;
   double lo,hi;
+  double maxradius;
+
+  void point_on_line_segment(double *, double *, double *, double *);
+  double closest(double *, double *, double *, double);
+
+  void subtract(double *, double *, double *);
+  double dotproduct(double *, double *);
 };
 
 }

src/region_cylinder.cpp

@@ -22,6 +22,9 @@ using namespace LAMMPS_NS;
 
 #define BIG 1.0e20
 
+#define MIN(A,B) ((A) < (B)) ? (A) : (B)
+#define MAX(A,B) ((A) > (B)) ? (A) : (B)
+
 /* ---------------------------------------------------------------------- */
 
 RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
@@ -36,14 +39,16 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
   if (axis == 'x') {
     c1 = yscale*atof(arg[3]);
     c2 = zscale*atof(arg[4]);
+    radius = yscale*atof(arg[5]);
   } else if (axis == 'y') {
     c1 = xscale*atof(arg[3]);
     c2 = zscale*atof(arg[4]);
+    radius = xscale*atof(arg[5]);
   } else if (axis == 'z') {
     c1 = xscale*atof(arg[3]);
     c2 = yscale*atof(arg[4]);
-  }
     radius = xscale*atof(arg[5]);
+  }
 
   if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) {
     if (domain->box_exist == 0) 
@@ -95,10 +100,12 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
 
   // error check
 
-  if (radius < 0.0) error->all("Illegal region cylinder command");
+  if (radius <= 0.0) error->all("Illegal region cylinder command");
 
   // extent of cylinder
 
+  if (interior) {
+    bboxflag = 1;
     if (axis == 'x') {
       extent_xlo = lo;
       extent_xhi = hi;
@@ -123,10 +130,26 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
       extent_zlo = lo;
       extent_zhi = hi;
     }
+  } else bboxflag = 0;
+
+  // particle could be contact cylinder surface and 2 ends
+
+  cmax = 3;
+  contact = new Contact[cmax];
 }
 
 /* ---------------------------------------------------------------------- */
 
+RegCylinder::~RegCylinder()
+{
+  delete [] contact;
+}
+
+/* ----------------------------------------------------------------------
+   inside = 1 if x,y,z is inside or on surface
+   inside = 0 if x,y,z is outside and not on surface
+------------------------------------------------------------------------- */
+
 int RegCylinder::match(double x, double y, double z)
 {
   double del1,del2,dist;
@@ -138,15 +161,13 @@ int RegCylinder::match(double x, double y, double z)
     dist = sqrt(del1*del1 + del2*del2);
     if (dist <= radius && x >= lo && x <= hi) inside = 1;
     else inside = 0;
-  }
+  } else if (axis == 'y') {
-  if (axis == 'y') {
     del1 = x - c1;
     del2 = z - c2;
     dist = sqrt(del1*del1 + del2*del2);
     if (dist <= radius && y >= lo && y <= hi) inside = 1;
     else inside = 0;
-  }
+  } else {
-  if (axis == 'z') {
     del1 = x - c1;
     del2 = y - c2;
     dist = sqrt(del1*del1 + del2*del2);
@@ -156,3 +177,229 @@ int RegCylinder::match(double x, double y, double z)
 
   return !(inside ^ interior);         // 1 if same, 0 if different
 }
+
+/* ----------------------------------------------------------------------
+   contact if 0 <= x < cutoff from one or more inner surfaces of cylinder
+   can be one contact for each of 3 cylinder surfaces
+   no contact if outside (possible if called from union/intersect)
+   delxyz = vector from nearest point on cylinder to x
+   special case: no contact with curved surf if x is on center axis
+------------------------------------------------------------------------- */
+
+int RegCylinder::surface_interior(double *x, double cutoff)
+{
+  double del1,del2,r,delta;
+
+  int n = 0;
+
+  if (axis == 'x') {
+    del1 = x[1] - c1;
+    del2 = x[2] - c2;
+    r = sqrt(del1*del1 + del2*del2);
+
+    // x is exterior to cylinder
+    
+    if (r > radius || x[0] < lo || x[0] > hi) return 0;
+
+    // x is interior to cylinder or on its surface
+
+    delta = radius - r;
+    if (delta < cutoff && r > 0.0) {
+      contact[n].r = delta;
+      contact[n].delx = 0.0;
+      contact[n].dely = del1*(1.0-radius/r);
+      contact[n].delz = del2*(1.0-radius/r);
+      n++;
+    }
+    delta = x[0] - lo;
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].delx = delta;
+      contact[n].dely = contact[n].delz = 0.0;
+      n++;
+    }
+    delta = hi - x[0];
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].delx = -delta;
+      contact[n].dely = contact[n].delz = 0.0;
+      n++;
+    }
+
+  } else if (axis == 'y') {
+    del1 = x[0] - c1;
+    del2 = x[2] - c2;
+    r = sqrt(del1*del1 + del2*del2);
+
+    // x is exterior to cylinder
+    
+    if (r > radius || x[1] < lo || x[1] > hi) return 0;
+
+    // x is interior to cylinder or on its surface
+
+    delta = radius - r;
+    if (delta < cutoff && r > 0.0) {
+      contact[n].r = delta;
+      contact[n].delx = del1*(1.0-radius/r);
+      contact[n].dely = 0.0;
+      contact[n].delz = del2*(1.0-radius/r);
+      n++;
+    }
+    delta = x[1] - lo;
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].dely = delta;
+      contact[n].delx = contact[n].delz = 0.0;
+      n++;
+    }
+    delta = hi - x[1];
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].dely = -delta;
+      contact[n].delx = contact[n].delz = 0.0;
+      n++;
+    }
+
+  } else {
+    del1 = x[0] - c1;
+    del2 = x[1] - c2;
+    r = sqrt(del1*del1 + del2*del2);
+
+    // x is exterior to cylinder
+    
+    if (r > radius || x[2] < lo || x[2] > hi) return 0;
+
+    // x is interior to cylinder or on its surface
+
+    delta = radius - r;
+    if (delta < cutoff && r > 0.0) {
+      contact[n].r = delta;
+      contact[n].delx = del1*(1.0-radius/r);
+      contact[n].dely = del2*(1.0-radius/r);
+      contact[n].delz = 0.0;
+      n++;
+    }
+    delta = x[2] - lo;
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].delz = delta;
+      contact[n].delx = contact[n].dely = 0.0;
+      n++;
+    }
+    delta = hi - x[2];
+    if (delta < cutoff) {
+      contact[n].r = delta;
+      contact[n].delz = -delta;
+      contact[n].delx = contact[n].dely = 0.0;
+      n++;
+    }
+  }
+
+  return n;
+}
+
+/* ----------------------------------------------------------------------
+   one contact if 0 <= x < cutoff from outer surface of cylinder
+   no contact if inside (possible if called from union/intersect)
+   delxyz = vector from nearest point on cylinder to x
+------------------------------------------------------------------------- */
+
+int RegCylinder::surface_exterior(double *x, double cutoff)
+{
+  double del1,del2,r,delta;
+  double xp,yp,zp;
+
+  if (axis == 'x') {
+    del1 = x[1] - c1;
+    del2 = x[2] - c2;
+    r = sqrt(del1*del1 + del2*del2);
+
+    // x is far enough from cylinder that there is no contact
+    // x is interior to cylinder
+
+    if (r >= radius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff) return 0;
+    if (r < radius && x[0] > lo && x[0] < hi) return 0;
+
+    // x is exterior to cylinder or on its surface
+    // xp,yp,zp = point on surface of cylinder that x is closest to
+    //            could be edge of cylinder
+    // do not add contact point if r >= cutoff
+
+    if (r > radius) {
+      yp = c1 + del1*radius/r;
+      zp = c2 + del2*radius/r;
+    } else {
+      yp = x[1];
+      zp = x[2];
+    }
+    if (x[0] < lo) xp = lo;
+    else if (x[0] > hi) xp = hi;
+    else xp = x[0];
+
+    add_contact(0,x,xp,yp,zp);
+    if (contact[0].r < cutoff) return 1;
+    return 0;
+
+  } else if (axis == 'y') {
+    del1 = x[0] - c1;
+    del2 = x[2] - c2;
+    r = sqrt(del1*del1 + del2*del2);
+
+    // x is far enough from cylinder that there is no contact
+    // x is interior to cylinder
+
+    if (r >= radius+cutoff || x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0;
+    if (r < radius && x[1] > lo && x[1] < hi) return 0;
+
+    // x is exterior to cylinder or on its surface
+    // xp,yp,zp = point on surface of cylinder that x is closest to
+    //            could be edge of cylinder
+    // do not add contact point if r >= cutoff
+
+    if (r > radius) {
+      xp = c1 + del1*radius/r;
+      zp = c2 + del2*radius/r;
+    } else {
+      xp = x[0];
+      zp = x[2];
+    }
+    if (x[1] < lo) yp = lo;
+    else if (x[1] > hi) yp = hi;
+    else yp = x[1];
+
+    add_contact(0,x,xp,yp,zp);
+    if (contact[0].r < cutoff) return 1;
+    return 0;
+
+  } else {
+    del1 = x[0] - c1;
+    del2 = x[1] - c2;
+    r = sqrt(del1*del1 + del2*del2);
+
+    // x is far enough from cylinder that there is no contact
+    // x is interior to cylinder
+
+    if (r >= radius+cutoff || x[2] <= lo-cutoff || x[2] >= hi+cutoff) return 0;
+    if (r < radius && x[2] > lo && x[2] < hi) return 0;
+
+    // x is exterior to cylinder or on its surface
+    // xp,yp,zp = point on surface of cylinder that x is closest to
+    //            could be edge of cylinder
+    // do not add contact point if r >= cutoff
+
+    if (r > radius) {
+      xp = c1 + del1*radius/r;
+      yp = c2 + del2*radius/r;
+    } else {
+      xp = x[0];
+      yp = x[1];
+    }
+    if (x[2] < lo) zp = lo;
+    else if (x[2] > hi) zp = hi;
+    else zp = x[2];
+
+    add_contact(0,x,xp,yp,zp);
+    if (contact[0].r < cutoff) return 1;
+    return 0;
+  }
+}

src/region_cylinder.h

@@ -23,7 +23,10 @@ class RegCylinder : public Region {
 
  public:
   RegCylinder(class LAMMPS *, int, char **);
+  ~RegCylinder();
   int match(double, double, double);
+  int surface_interior(double *, double);
+  int surface_exterior(double *, double);
 
  private:
   char axis;

src/region_intersect.cpp

@@ -45,17 +45,29 @@ RegIntersect::RegIntersect(LAMMPS *lmp, int narg, char **arg) :
   }
 
   // extent of intersection of regions
+  // has bounding box if interior and any sub-region has bounding box
 
   Region **regions = domain->regions;
 
-  extent_xlo = regions[list[0]]->extent_xlo;
+  bboxflag = 0;
-  extent_ylo = regions[list[0]]->extent_ylo;
+  for (int ilist = 0; ilist < nregion; ilist++)
-  extent_zlo = regions[list[0]]->extent_zlo;
+    if (regions[list[ilist]]->bboxflag == 1) bboxflag = 1;
-  extent_xhi = regions[list[0]]->extent_xhi;
+  if (!interior) bboxflag = 0;
-  extent_yhi = regions[list[0]]->extent_yhi;
+
-  extent_zhi = regions[list[0]]->extent_zhi;
+  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;
+      }
 
-  for (int ilist = 1; ilist < nregion; ilist++) {
       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);
@@ -63,6 +75,14 @@ RegIntersect::RegIntersect(LAMMPS *lmp, int narg, char **arg) :
       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];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -70,9 +90,13 @@ RegIntersect::RegIntersect(LAMMPS *lmp, int narg, char **arg) :
 RegIntersect::~RegIntersect()
 {
   delete [] list;
+  delete [] contact;
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   inside = 1 if x,y,z is match() with all sub-regions
+   else inside = 0
+------------------------------------------------------------------------- */
 
 int RegIntersect::match(double x, double y, double z)
 {
@@ -81,9 +105,95 @@ int RegIntersect::match(double x, double y, double z)
   for (ilist = 0; ilist < nregion; ilist++)
     if (!regions[list[ilist]]->match(x,y,z)) break;
 
-  int inside;                          // inside if matched all regions
+  int inside = 0;                      // inside if matched all regions
   if (ilist == nregion) inside = 1;
-  else inside = 0;
 
   return !(inside ^ interior);         // 1 if same, 0 if different
 }
+
+/* ----------------------------------------------------------------------
+   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,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,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;
+}

src/region_intersect.h

@@ -23,6 +23,8 @@ class RegIntersect : public Region {
   RegIntersect(class LAMMPS *, int, char **);
   ~RegIntersect();
   int match(double, double, double);
+  int surface_interior(double *, double);
+  int surface_exterior(double *, double);
 
  private:
   int nregion;

src/region_prism.cpp

@@ -15,6 +15,7 @@
    Contributing author: Pieter in 't Veld (SNL)
 ------------------------------------------------------------------------- */
 
+#include "math.h"
 #include "stdlib.h"
 #include "string.h"
 #include "region_prism.h"
@@ -83,12 +84,30 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
 
   // error check
   // prism cannot be 0 thickness in any dim, else inverse blows up
+  // non-zero tilt values cannot be used if either dim is INF on both ends
 
   if (xlo >= xhi || ylo >= yhi || zlo >= zhi)
     error->all("Illegal region prism command");
 
+  if (xy != 0.0 && xlo == -BIG && xhi == BIG)
+    error->all("Illegal region prism command");
+  if (xy != 0.0 && ylo == -BIG && yhi == BIG)
+    error->all("Illegal region prism command");
+
+  if (xz != 0.0 && xlo == -BIG && xhi == BIG)
+    error->all("Illegal region prism command");
+  if (xz != 0.0 && zlo == -BIG && zhi == BIG)
+    error->all("Illegal region prism command");
+
+  if (yz != 0.0 && ylo == -BIG && yhi == BIG)
+    error->all("Illegal region prism command");
+  if (yz != 0.0 && zlo == -BIG && zhi == BIG)
+    error->all("Illegal region prism command");
+
   // extent of prism
   
+  if (interior) {
+    bboxflag = 1;
     extent_xlo = MIN(xlo,xlo+xy);
     extent_xlo = MIN(extent_xlo,extent_xlo+xz);
     extent_ylo = MIN(ylo,ylo+yz);
@@ -98,6 +117,12 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
     extent_xhi = MAX(extent_xhi,extent_xhi+xz);
     extent_yhi = MAX(yhi,yhi+yz);
     extent_zhi = zhi;
+  } else bboxflag = 0;
+
+  // particle could contact all 6 planes
+
+  cmax = 6;
+  contact = new Contact[cmax];
 
   // h = transformation matrix from tilt coords (0-1) to box coords (xyz)
   // columns of h are edge vectors of tilted box
@@ -119,10 +144,95 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
   hinv[1][1] = 1.0/h[1][1];
   hinv[1][2] = -h[1][2] / (h[1][1]*h[2][2]);
   hinv[2][2] = 1.0/h[2][2];
+
+  // corners = 8 corner points of prism
+  // order = x varies fastest, then y, finally z
+  // clo/chi = lo and hi corner pts of prism
+
+  a[0] = xhi-xlo;
+  a[1] = 0.0;
+  a[2] = 0.0;
+  b[0] = xy;
+  b[1] = yhi-ylo;
+  b[2] = 0.0;
+  c[0] = xz;
+  c[1] = yz;
+  c[2] = zhi-zlo;
+
+  clo[0] = corners[0][0] = xlo;
+  clo[1] = corners[0][1] = ylo;
+  clo[2] = corners[0][2] = zlo;
+
+  corners[1][0] = xlo + a[0];
+  corners[1][1] = ylo + a[1];
+  corners[1][2] = zlo + a[2];
+
+  corners[2][0] = xlo + b[0];
+  corners[2][1] = ylo + b[1];
+  corners[2][2] = zlo + b[2];
+
+  corners[3][0] = xlo + a[0] + b[0];
+  corners[3][1] = ylo + a[1] + b[1];
+  corners[3][2] = zlo + a[2] + b[2];
+
+  corners[4][0] = xlo + c[0];
+  corners[4][1] = ylo + c[1];
+  corners[4][2] = zlo + c[2];
+
+  corners[5][0] = xlo + a[0] + c[0];
+  corners[5][1] = ylo + a[1] + c[1];
+  corners[5][2] = zlo + a[2] + c[2];
+
+  corners[6][0] = xlo + b[0] + c[0];
+  corners[6][1] = ylo + b[1] + c[1];
+  corners[6][2] = zlo + b[2] + c[2];
+
+  chi[0] = corners[7][0] = xlo + a[0] + b[0] + c[0];
+  chi[1] = corners[7][1] = ylo + a[1] + b[1] + c[1];
+  chi[2] = corners[7][2] = zlo + a[2] + b[2] + c[2];
+
+  // face = 6 inward-facing unit normals to prism faces
+  // order = xy plane, xz plane, yz plane
+
+  cross(a,b,face[0]);
+  cross(b,a,face[1]);
+  cross(c,a,face[2]);
+  cross(a,c,face[3]);
+  cross(b,c,face[4]);
+  cross(c,b,face[5]);
+
+  for (int i = 0; i < 6; i++) normalize(face[i]);
+
+  // tri = 3 vertices (0-7) in each of 12 triangles on 6 faces
+  // verts in each tri are ordered so that right-hand rule gives inward norm
+  // order = xy plane, xz plane, yz plane
+
+  tri[0][0] =  0; tri[0][1] =  1; tri[0][2] =  3;
+  tri[1][0] =  0; tri[1][1] =  3; tri[1][2] =  2;
+  tri[2][0] =  4; tri[2][1] =  7; tri[2][2] =  5;
+  tri[3][0] =  4; tri[3][1] =  6; tri[3][2] =  7;
+
+  tri[4][0] =  0; tri[4][1] =  4; tri[4][2] =  5;
+  tri[5][0] =  0; tri[5][1] =  5; tri[5][2] =  1;
+  tri[6][0] =  2; tri[6][1] =  7; tri[6][2] =  6;
+  tri[7][0] =  2; tri[7][1] =  3; tri[7][2] =  7;
+
+  tri[8][0] =  2; tri[8][1] =  6; tri[8][2] =  4;
+  tri[9][0] =  2; tri[9][1] =  4; tri[9][2] =  0;
+  tri[10][0] = 1; tri[10][1] = 5; tri[10][2] = 7;
+  tri[11][0] = 1; tri[11][1] = 7; tri[11][2] = 3;
+}
+
+/* ---------------------------------------------------------------------- */
+
+RegPrism::~RegPrism()
+{
+  delete [] contact;
 }
 
 /* ----------------------------------------------------------------------
-   check xyz against prism
+   inside = 1 if x,y,z is inside or on surface
+   inside = 0 if x,y,z is outside and not on surface
    abc = Hinv * (xyz - xyz/lo)
    abc = tilt coords (0-1)
    Hinv = transformation matrix from box coords to tilt coords
@@ -136,10 +246,270 @@ int RegPrism::match(double x, double y, double z)
   double b = hinv[1][1]*(y-ylo) + hinv[1][2]*(z-zlo);
   double c = hinv[2][2]*(z-zlo);
 
-  int inside;
+  int inside = 0;
   if (a >= 0.0 && a <= 1.0 && b >= 0.0 && b <= 1.0 && c >= 0.0 && c <= 1.0)
     inside = 1;
-  else inside = 0;
 
   return !(inside ^ interior);         // 1 if same, 0 if different
 }
+
+/* ----------------------------------------------------------------------
+   contact if 0 <= x < cutoff from one or more inner surfaces of prism
+   can be one contact for each of 6 faces
+   no contact if outside (possible if called from union/intersect)
+   delxyz = vector from nearest point on prism to x
+------------------------------------------------------------------------- */
+
+int RegPrism::surface_interior(double *x, double cutoff)
+{
+  int i;
+  double dot;
+  double *corner;
+
+  // x is exterior to prism
+
+  for (i = 0; i < 6; i++) {
+    if (i % 2) corner = chi;
+    else corner = clo;
+    dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + 
+      (x[2]-corner[2])*face[i][2];
+    if (dot < 0.0) return 0;
+  }
+
+  // x is interior to prism or on its surface
+
+  int n = 0;
+
+  for (int i = 0; i < 6; i++) {
+    if (i % 2) corner = chi;
+    else corner = clo;
+    dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + 
+      (x[2]-corner[2])*face[i][2];
+    if (dot < cutoff) {
+      contact[n].r = dot;
+      contact[n].delx = dot*face[i][0];
+      contact[n].dely = dot*face[i][1];
+      contact[n].delz = dot*face[i][2];
+      n++;
+    }
+  }
+
+  return n;
+}
+
+/* ----------------------------------------------------------------------
+   one contact if 0 <= x < cutoff from outer surface of prism
+   no contact if inside (possible if called from union/intersect)
+   delxyz = vector from nearest point on prism to x
+------------------------------------------------------------------------- */
+
+int RegPrism::surface_exterior(double *x, double cutoff)
+{
+  int i;
+  double dot,delta;
+  double *corner;
+  double xp,yp,zp;
+
+  // x is far enough from prism that there is no contact
+
+  for (i = 0; i < 6; i++) {
+    if (i % 2) corner = chi;
+    else corner = clo;
+    dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + 
+      (x[2]-corner[2])*face[i][2];
+    if (dot <= -cutoff) return 0;
+  }
+
+  // x is interior to prism
+
+  for (i = 0; i < 6; i++) {
+    if (i % 2) corner = chi;
+    else corner = clo;
+    dot = (x[0]-corner[0])*face[i][0] + (x[1]-corner[1])*face[i][1] + 
+      (x[2]-corner[2])*face[i][2];
+    if (dot <= 0.0) break;
+  }
+
+  if (i == 6) return 0;
+
+  // x is exterior to prism or on its surface
+  // xp,yp,zp = point on surface of prism that x is closest to
+  //            could be edge or corner pt of prism
+  // do not add contact point if r >= cutoff
+
+  find_nearest(x,xp,yp,zp);
+  add_contact(0,x,xp,yp,zp);
+  if (contact[0].r < cutoff) return 1;
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   x is exterior to prism or on its surface
+   return (xp,yp,zp) = nearest pt to x that is on surface of prism
+------------------------------------------------------------------------- */
+
+void RegPrism::find_nearest(double *x, double &xp, double &yp, double &zp)
+{
+  int i,j,k,iface;
+  double xproj[3],xline[3],nearest[3];
+  double dot;
+
+  // generate successive xnear points, one nearest to x is (xp,yp,zp)
+  // loop over 6 faces and 2 triangles in each face
+  // xproj = x projected to plane of triangle
+  // if xproj is inside or on triangle boundary, that is xnear
+  // else: loop over 3 edges of triangle
+  //       compute distance to edge line
+  //       xnear = nearest point on line to xproj, bounded by segment end pts
+
+  double distsq = BIG;
+
+  for (int itri = 0; itri < 12; itri++) {
+    iface = itri/2;
+    i = tri[itri][0];
+    j = tri[itri][1];
+    k = tri[itri][2];
+    dot = (x[0]-corners[i][0])*face[iface][0] + 
+      (x[1]-corners[i][1])*face[iface][1] + 
+      (x[2]-corners[i][2])*face[iface][2];
+    xproj[0] = x[0] - dot*face[iface][0];
+    xproj[1] = x[1] - dot*face[iface][1];
+    xproj[2] = x[2] - dot*face[iface][2];
+    if (inside_tri(xproj,corners[i],corners[j],corners[k],face[iface]))
+      distsq = closest(x,xproj,nearest,distsq);
+    else {
+      point_on_line_segment(corners[i],corners[j],xproj,xline);
+      distsq = closest(x,xline,nearest,distsq);
+      point_on_line_segment(corners[j],corners[k],xproj,xline);
+      distsq = closest(x,xline,nearest,distsq);
+      point_on_line_segment(corners[i],corners[k],xproj,xline);
+      distsq = closest(x,xline,nearest,distsq);
+    }
+  }
+  
+  xp = nearest[0];
+  yp = nearest[1];
+  zp = nearest[2];
+}
+
+/* ----------------------------------------------------------------------
+   test if x is inside triangle with vertices v1,v2,v3
+   norm = normal to triangle, defined by right-hand rule for v1,v2,v3 ordering
+   edge = edge vector of triangle
+   pvec = vector from vertex to x
+   xproduct = cross product of edge with pvec
+   if xproduct dot norm < 0.0 for any of 3 edges, then x is outside triangle
+------------------------------------------------------------------------- */
+
+int RegPrism::inside_tri(double *x, double *v1, double *v2, double *v3,
+			 double *norm)
+{
+  double edge[3],pvec[3],xproduct[3];
+
+  subtract(v1,v2,edge);
+  subtract(v1,x,pvec);
+  cross(edge,pvec,xproduct);
+  if (dotproduct(xproduct,norm) < 0.0) return 0;
+
+  subtract(v2,v3,edge);
+  subtract(v2,x,pvec);
+  cross(edge,pvec,xproduct);
+  if (dotproduct(xproduct,norm) < 0.0) return 0;
+
+  subtract(v3,v1,edge);
+  subtract(v3,x,pvec);
+  cross(edge,pvec,xproduct);
+  if (dotproduct(xproduct,norm) < 0.0) return 0;
+
+  return 1;
+}
+
+/* ----------------------------------------------------------------------
+   find nearest point to C on line segment A,B and return it as D
+   project (C-A) onto (B-A)
+   t = length of that projection, normalized by length of (B-A)
+   t <= 0, C is closest to A
+   t >= 1, C is closest to B
+   else closest point is between A and B
+------------------------------------------------------------------------- */
+
+void RegPrism::point_on_line_segment(double *a, double *b, 
+				     double *c, double *d)
+{
+  double ba[3],ca[3];
+
+  subtract(a,b,ba);
+  subtract(a,c,ca);
+  double t = dotproduct(ca,ba) / dotproduct(ba,ba);
+
+  if (t <= 0.0) {
+    d[0] = a[0];
+    d[1] = a[1];
+    d[2] = a[2];
+  } else if (t >= 1.0) {
+    d[0] = b[0];
+    d[1] = b[1];
+    d[2] = b[2];
+  } else {
+    d[0] = a[0] + t*ba[0];
+    d[1] = a[1] + t*ba[1];
+    d[2] = a[2] + t*ba[2];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double RegPrism::closest(double *x, double *near, double *nearest, double dsq)
+{
+  double delx = x[0] - near[0];
+  double dely = x[1] - near[1];
+  double delz = x[2] - near[2];
+  double rsq = delx*delx + dely*dely + delz*delz;
+  if (rsq >= dsq) return dsq;
+
+  nearest[0] = near[0];
+  nearest[1] = near[1];
+  nearest[2] = near[2];
+  return rsq;
+}
+
+/* ----------------------------------------------------------------------
+   v3 = v2 - v1
+------------------------------------------------------------------------- */
+
+void RegPrism::subtract(double *v1, double *v2, double *v3)
+{
+  v3[0] = v2[0] - v1[0];
+  v3[1] = v2[1] - v1[1];
+  v3[2] = v2[2] - v1[2];
+}
+
+/* ----------------------------------------------------------------------
+   v3 = v1 x v2
+------------------------------------------------------------------------- */
+
+void RegPrism::cross(double *v1, double *v2, double *v3)
+{
+  v3[0] = v1[1]*v2[2] - v1[2]*v2[1];
+  v3[1] = v1[2]*v2[0] - v1[0]*v2[2];
+  v3[2] = v1[0]*v2[1] - v1[1]*v2[0];
+}
+
+/* ----------------------------------------------------------------------
+   return dotproduct = v1 dot v2
+------------------------------------------------------------------------- */
+
+double RegPrism::dotproduct(double *v1, double *v2)
+{
+  return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void RegPrism::normalize(double *x)
+{
+  double invlen = 1.0/sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
+  x[0] *= invlen;
+  x[1] *= invlen;
+  x[2] *= invlen;
+}

src/region_prism.h

@@ -23,12 +23,31 @@ class RegPrism : public Region {
 
  public:
   RegPrism(class LAMMPS *, int, char **);
+  ~RegPrism();
   int match(double, double, double);
+  int surface_interior(double *, double);
+  int surface_exterior(double *, double);
 
  private:
   double xlo,xhi,ylo,yhi,zlo,zhi;
   double xy,xz,yz;
   double h[3][3],hinv[3][3];
+  int dimension;
+  double a[3],b[3],c[3];       // edge vectors of region
+  double clo[3],chi[3];        // opposite corners of prism
+  double face[6][3];           // unit normals of 6 prism faces
+  double corners[8][3];        // 8 corner pts of prism
+  int tri[12][3];              // 3 corner pts of 12 triangles (2 per face)
+
+  void find_nearest(double *, double &, double &, double &);
+  int inside_tri(double *, double *, double *, double *, double *);
+  void point_on_line_segment(double *, double *, double *, double *);
+  double closest(double *, double *, double *, double);
+
+  void subtract(double *, double *, double *);
+  void cross(double *, double *, double *);
+  double dotproduct(double *, double *);
+  void normalize(double *);
 };
 
 }

src/region_sphere.cpp

@@ -37,26 +37,92 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
 
   // extent of sphere
 
+  if (interior) {
+    bboxflag = 1;
     extent_xlo = xc - radius;
     extent_xhi = xc + radius;
     extent_ylo = yc - radius;
     extent_yhi = yc + radius;
     extent_zlo = zc - radius;
     extent_zhi = zc + radius;
+  } else bboxflag = 0;
+
+  cmax = 1;
+  contact = new Contact[cmax];
 }
 
 /* ---------------------------------------------------------------------- */
 
+RegSphere::~RegSphere()
+{
+  delete [] contact;
+}
+
+/* ----------------------------------------------------------------------
+   inside = 1 if x,y,z is inside or on surface
+   inside = 0 if x,y,z is outside and not on surface
+------------------------------------------------------------------------- */
+
 int RegSphere::match(double x, double y, double z)
 {
   double delx = x - xc;
   double dely = y - yc;
   double delz = z - zc;
-  double dist = sqrt(delx*delx + dely*dely + delz*delz);
+  double r = sqrt(delx*delx + dely*dely + delz*delz);
 
-  int inside;
+  int inside = 0;
-  if (dist <= radius) inside = 1;
+  if (r <= radius) inside = 1;
-  else inside = 0;
 
   return !(inside ^ interior);         // 1 if same, 0 if different
 }
+
+/* ----------------------------------------------------------------------
+   one contact if 0 <= x < cutoff away 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
+------------------------------------------------------------------------- */
+
+int RegSphere::surface_interior(double *x, double cutoff)
+{
+  double delx = x[0] - xc;
+  double dely = x[1] - yc;
+  double delz = x[2] - zc;
+  double r = sqrt(delx*delx + dely*dely + delz*delz);
+  if (r > radius || r == 0.0) return 0;
+
+  double delta = radius - r;
+  if (delta < cutoff) {
+    contact[0].r = delta;
+    contact[0].delx = delx*(1.0-radius/r);
+    contact[0].dely = dely*(1.0-radius/r);
+    contact[0].delz = delz*(1.0-radius/r);
+    return 1;
+  }
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   one contact if 0 <= x < cutoff away from outer surface of sphere
+   no contact if inside (possible if called from union/intersect)
+   delxyz = vector from nearest point on sphere to x
+------------------------------------------------------------------------- */
+
+int RegSphere::surface_exterior(double *x, double cutoff)
+{
+  double delx = x[0] - xc;
+  double dely = x[1] - yc;
+  double delz = x[2] - zc;
+  double r = sqrt(delx*delx + dely*dely + delz*delz);
+  if (r < radius) return 0;
+
+  double delta = r - radius;
+  if (delta < cutoff) {
+    contact[0].r = delta;
+    contact[0].delx = delx*(1.0-radius/r);
+    contact[0].dely = dely*(1.0-radius/r);
+    contact[0].delz = delz*(1.0-radius/r);
+    return 1;
+  }
+  return 0;
+}

src/region_sphere.h

@@ -21,7 +21,10 @@ namespace LAMMPS_NS {
 class RegSphere : public Region {
  public:
   RegSphere(class LAMMPS *, int, char **);
+  ~RegSphere();
   int match(double, double, double);
+  int surface_interior(double *, double);
+  int surface_exterior(double *, double);
 
  private:
   double xc,yc,zc;

src/region_union.cpp

@@ -46,9 +46,16 @@ RegUnion::RegUnion(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
   }
 
   // extent of union of regions
+  // has bounding box if interior and all sub-regions have bounding box
 
   Region **regions = domain->regions;
 
+  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;
 
@@ -60,6 +67,14 @@ RegUnion::RegUnion(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
       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];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -67,9 +82,13 @@ RegUnion::RegUnion(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
 RegUnion::~RegUnion()
 {
   delete [] list;
+  delete [] contact;
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   inside = 1 if x,y,z is match() with any sub-region
+   else inside = 0
+------------------------------------------------------------------------- */
 
 int RegUnion::match(double x, double y, double z)
 {
@@ -78,9 +97,95 @@ int RegUnion::match(double x, double y, double z)
   for (ilist = 0; ilist < nregion; ilist++)
     if (regions[list[ilist]]->match(x,y,z)) break;
 
-  int inside;                          // inside if matched any region
+  int inside = 1;                      // inside if matched any region
   if (ilist == nregion) inside = 0;
-  else inside = 1;
 
   return !(inside ^ interior);         // 1 if same, 0 if different
 }
+
+/* ----------------------------------------------------------------------
+   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,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,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;
+}

src/region_union.h

@@ -23,6 +23,8 @@ class RegUnion : public Region {
   RegUnion(class LAMMPS *, int, char **);
   ~RegUnion();
   int match(double, double, double);
+  int surface_interior(double *, double);
+  int surface_exterior(double *, double);
 
  private:
   int nregion;

src/style.h

@@ -236,9 +236,11 @@ DumpStyle(xyz,DumpXYZ)
 #include "fix_ttm.h"
 #include "fix_viscosity.h"
 #include "fix_viscous.h"
+#include "fix_wall_harmonic.h"
 #include "fix_wall_lj126.h"
 #include "fix_wall_lj93.h"
 #include "fix_wall_reflect.h"
+#include "fix_wall_region.h"
 #endif
 
 #ifdef FixClass
@@ -296,9 +298,11 @@ FixStyle(tmd,FixTMD)
 FixStyle(ttm,FixTTM)
 FixStyle(viscosity,FixViscosity)
 FixStyle(viscous,FixViscous)
+FixStyle(wall/harmonic,FixWallHarmonic)
 FixStyle(wall/lj126,FixWallLJ126)
 FixStyle(wall/lj93,FixWallLJ93)
 FixStyle(wall/reflect,FixWallReflect)
+FixStyle(wall/region,FixWallRegion)
 #endif
 
 #ifdef ImproperInclude
@@ -384,6 +388,7 @@ PairStyle(yukawa,PairYukawa)
 #include "region_cone.h"
 #include "region_cylinder.h"
 #include "region_intersect.h"
+#include "region_plane.h"
 #include "region_prism.h"
 #include "region_sphere.h"
 #include "region_union.h"
@@ -394,6 +399,7 @@ RegionStyle(block,RegBlock)
 RegionStyle(cone,RegCone)
 RegionStyle(cylinder,RegCylinder)
 RegionStyle(intersect,RegIntersect)
+RegionStyle(plane,RegPlane)
 RegionStyle(prism,RegPrism)
 RegionStyle(sphere,RegSphere)
 RegionStyle(union,RegUnion)