diff --git a/src/COLLOID/fix_wall_colloid.cpp b/src/COLLOID/fix_wall_colloid.cpp index 713235f465..8fcef726c6 100644 --- a/src/COLLOID/fix_wall_colloid.cpp +++ b/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"); } diff --git a/src/Makefile.package b/src/Makefile.package index e6ba0bbcff..0c8aa9d850 100644 --- a/src/Makefile.package +++ b/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_LIB = -lgpu -lreax -lpoems -lmeam +PKG_PATH = -L../../lib/reax -L../../lib/poems -L../../lib/meam +PKG_LIB = -lreax -lpoems -lmeam -PKG_SYSPATH = $(gpu_SYSPATH) $(reax_SYSPATH) $(meam_SYSPATH) -PKG_SYSLIB = $(gpu_SYSLIB) $(reax_SYSLIB) $(meam_SYSLIB) +PKG_SYSPATH = $(reax_SYSPATH) $(meam_SYSPATH) +PKG_SYSLIB = $(reax_SYSLIB) $(meam_SYSLIB) diff --git a/src/create_box.cpp b/src/create_box.cpp index 7d4b9f7cb8..cc9ad9513e 100644 --- a/src/create_box.cpp +++ b/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) - error->all("Create_box region must be of type inside"); + if (domain->regions[iregion]->bboxflag == 0) + error->all("Create_box region does not support a bounding box"); // if region not prism: // setup orthogonal domain diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp index 8e0dc9ec86..4b446f2da8 100644 --- a/src/fix_deposit.cpp +++ b/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) - error->all("Must use region with side = in with fix deposit"); + if (domain->regions[iregion]->bboxflag == 0) + error->all("Fix deposit region does not support a bounding box"); xlo = domain->regions[iregion]->extent_xlo; xhi = domain->regions[iregion]->extent_xhi; diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp index 4d9def3847..fee8d27fa2 100644 --- a/src/fix_wall.cpp +++ b/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"); diff --git a/src/fix_wall_lj126.cpp b/src/fix_wall_lj126.cpp index f70c23f9e9..1379676a49 100644 --- a/src/fix_wall_lj126.cpp +++ b/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"); } diff --git a/src/fix_wall_lj93.cpp b/src/fix_wall_lj93.cpp index ea180d6448..99f5f90903 100644 --- a/src/fix_wall_lj93.cpp +++ b/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"); } diff --git a/src/fix_wall_region.cpp b/src/fix_wall_region.cpp new file mode 100644 index 0000000000..b738890b94 --- /dev/null +++ b/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; +} diff --git a/src/fix_wall_region.h b/src/fix_wall_region.h new file mode 100644 index 0000000000..e79ea7ef79 --- /dev/null +++ b/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 diff --git a/src/region.cpp b/src/region.cpp index 5a20cc10a9..b034e1a2f0 100644 --- a/src/region.cpp +++ b/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; +} diff --git a/src/region.h b/src/region.h index 804c39797e..1bad9faf84 100644 --- a/src/region.h +++ b/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); }; } diff --git a/src/region_block.cpp b/src/region_block.cpp index 63daa72d4e..f63dc66883 100644 --- a/src/region_block.cpp +++ b/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) @@ -81,23 +84,148 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) error->all("Illegal region block command"); // extent of block - - extent_xlo = xlo; - extent_xhi = xhi; - extent_ylo = ylo; - extent_yhi = yhi; - extent_zlo = zlo; - extent_zhi = zhi; + + 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; +} diff --git a/src/region_block.h b/src/region_block.h index 5c7de84769..5f61764398 100644 --- a/src/region_block.h +++ b/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; diff --git a/src/region_cone.cpp b/src/region_cone.cpp index f33d30c436..dc53f6ddfb 100644 --- a/src/region_cone.cpp +++ b/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]); } - 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) largestradius = radiuslo; - else largestradius = radiushi; + if (radiuslo > radiushi) maxradius = radiuslo; + else maxradius = radiushi; - if (axis == 'x') { - extent_xlo = lo; - extent_xhi = hi; - extent_ylo = c1 - largestradius; - extent_yhi = c1 + largestradius; - extent_zlo = c2 - largestradius; - extent_zhi = c2 + largestradius; - } - if (axis == 'y') { - extent_xlo = c1 - largestradius; - extent_xhi = c1 + largestradius; - extent_ylo = lo; - extent_yhi = hi; - extent_zlo = c2 - largestradius; - extent_zhi = c2 + largestradius; - } - if (axis == 'z') { - extent_xlo = c1 - largestradius; - extent_xhi = c1 + largestradius; - extent_ylo = c2 - largestradius; - extent_yhi = c2 + largestradius; - extent_zlo = lo; - extent_zhi = hi; - } + if (interior) { + bboxflag = 1; + + if (axis == 'x') { + extent_xlo = lo; + extent_xhi = hi; + extent_ylo = c1 - maxradius; + extent_yhi = c1 + maxradius; + extent_zlo = c2 - maxradius; + extent_zhi = c2 + maxradius; + } + if (axis == 'y') { + extent_xlo = c1 - maxradius; + extent_xhi = c1 + maxradius; + extent_ylo = lo; + extent_yhi = hi; + extent_zlo = c2 - maxradius; + extent_zhi = c2 + maxradius; + } + if (axis == 'z') { + extent_xlo = c1 - maxradius; + extent_xhi = c1 + maxradius; + extent_ylo = c2 - maxradius; + 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]); +} diff --git a/src/region_cone.h b/src/region_cone.h index bbcca14a75..791a6d3f56 100644 --- a/src/region_cone.h +++ b/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 *); }; } diff --git a/src/region_cylinder.cpp b/src/region_cylinder.cpp index 9f8ac5b8e0..c81954c137 100644 --- a/src/region_cylinder.cpp +++ b/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]); } - radius = xscale*atof(arg[5]); if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) { if (domain->box_exist == 0) @@ -95,38 +100,56 @@ 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 (axis == 'x') { - extent_xlo = lo; - extent_xhi = hi; - extent_ylo = c1 - radius; - extent_yhi = c1 + radius; - extent_zlo = c2 - radius; - extent_zhi = c2 + radius; - } - if (axis == 'y') { - extent_xlo = c1 - radius; - extent_xhi = c1 + radius; - extent_ylo = lo; - extent_yhi = hi; - extent_zlo = c2 - radius; - extent_zhi = c2 + radius; - } - if (axis == 'z') { - extent_xlo = c1 - radius; - extent_xhi = c1 + radius; - extent_ylo = c2 - radius; - extent_yhi = c2 + radius; - extent_zlo = lo; - extent_zhi = hi; - } + if (interior) { + bboxflag = 1; + if (axis == 'x') { + extent_xlo = lo; + extent_xhi = hi; + extent_ylo = c1 - radius; + extent_yhi = c1 + radius; + extent_zlo = c2 - radius; + extent_zhi = c2 + radius; + } + if (axis == 'y') { + extent_xlo = c1 - radius; + extent_xhi = c1 + radius; + extent_ylo = lo; + extent_yhi = hi; + extent_zlo = c2 - radius; + extent_zhi = c2 + radius; + } + if (axis == 'z') { + extent_xlo = c1 - radius; + extent_xhi = c1 + radius; + extent_ylo = c2 - radius; + extent_yhi = c2 + radius; + 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; - } - if (axis == 'y') { + } else 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; - } - if (axis == 'z') { + } else { 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; + } +} diff --git a/src/region_cylinder.h b/src/region_cylinder.h index 782710599a..25273c114a 100644 --- a/src/region_cylinder.h +++ b/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; diff --git a/src/region_intersect.cpp b/src/region_intersect.cpp index 872cc3cc04..2ecfe70713 100644 --- a/src/region_intersect.cpp +++ b/src/region_intersect.cpp @@ -45,24 +45,44 @@ 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; - extent_ylo = regions[list[0]]->extent_ylo; - extent_zlo = regions[list[0]]->extent_zlo; - extent_xhi = regions[list[0]]->extent_xhi; - extent_yhi = regions[list[0]]->extent_yhi; - extent_zhi = regions[list[0]]->extent_zhi; + bboxflag = 0; + for (int ilist = 0; ilist < nregion; ilist++) + if (regions[list[ilist]]->bboxflag == 1) bboxflag = 1; + if (!interior) bboxflag = 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); - extent_xhi = MIN(extent_xhi,regions[list[ilist]]->extent_xhi); - extent_yhi = MIN(extent_yhi,regions[list[ilist]]->extent_yhi); - extent_zhi = MIN(extent_zhi,regions[list[ilist]]->extent_zhi); + if (bboxflag) { + int first = 1; + for (int ilist = 0; ilist < nregion; ilist++) { + if (regions[list[ilist]]->bboxflag == 0) continue; + if (first) { + extent_xlo = regions[list[ilist]]->extent_xlo; + extent_ylo = regions[list[ilist]]->extent_ylo; + extent_zlo = regions[list[ilist]]->extent_zlo; + extent_xhi = regions[list[ilist]]->extent_xhi; + extent_yhi = regions[list[ilist]]->extent_yhi; + extent_zhi = regions[list[ilist]]->extent_zhi; + first = 0; + } + + extent_xlo = MAX(extent_xlo,regions[list[ilist]]->extent_xlo); + extent_ylo = MAX(extent_ylo,regions[list[ilist]]->extent_ylo); + extent_zlo = MAX(extent_zlo,regions[list[ilist]]->extent_zlo); + extent_xhi = MIN(extent_xhi,regions[list[ilist]]->extent_xhi); + extent_yhi = MIN(extent_yhi,regions[list[ilist]]->extent_yhi); + extent_zhi = MIN(extent_zhi,regions[list[ilist]]->extent_zhi); + } } + + // possible contacts = sum of possible contacts in all sub-regions + + cmax = 0; + for (int ilist = 0; ilist < nregion; ilist++) + cmax += regions[list[ilist]]->cmax; + contact = new Contact[cmax]; } /* ---------------------------------------------------------------------- */ @@ -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; +} diff --git a/src/region_intersect.h b/src/region_intersect.h index 01ec72648d..8d7d76bf23 100644 --- a/src/region_intersect.h +++ b/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; diff --git a/src/region_prism.cpp b/src/region_prism.cpp index 6b157176c4..3ccb3ea181 100644 --- a/src/region_prism.cpp +++ b/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,21 +84,45 @@ 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 - extent_xlo = MIN(xlo,xlo+xy); - extent_xlo = MIN(extent_xlo,extent_xlo+xz); - extent_ylo = MIN(ylo,ylo+yz); - extent_zlo = zlo; + if (interior) { + bboxflag = 1; + extent_xlo = MIN(xlo,xlo+xy); + extent_xlo = MIN(extent_xlo,extent_xlo+xz); + extent_ylo = MIN(ylo,ylo+yz); + extent_zlo = zlo; + + extent_xhi = MAX(xhi,xhi+xy); + extent_xhi = MAX(extent_xhi,extent_xhi+xz); + extent_yhi = MAX(yhi,yhi+yz); + extent_zhi = zhi; + } else bboxflag = 0; - extent_xhi = MAX(xhi,xhi+xy); - extent_xhi = MAX(extent_xhi,extent_xhi+xz); - extent_yhi = MAX(yhi,yhi+yz); - extent_zhi = zhi; + // 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; +} diff --git a/src/region_prism.h b/src/region_prism.h index 2bdae8c9cd..39c87452ff 100644 --- a/src/region_prism.h +++ b/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 *); }; } diff --git a/src/region_sphere.cpp b/src/region_sphere.cpp index 8248fc1d48..95f6fdfe15 100644 --- a/src/region_sphere.cpp +++ b/src/region_sphere.cpp @@ -37,26 +37,92 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) : // extent of sphere - extent_xlo = xc - radius; - extent_xhi = xc + radius; - extent_ylo = yc - radius; - extent_yhi = yc + radius; - extent_zlo = zc - radius; - extent_zhi = zc + radius; + 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; - if (dist <= radius) inside = 1; - else inside = 0; + int inside = 0; + if (r <= radius) inside = 1; 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; +} diff --git a/src/region_sphere.h b/src/region_sphere.h index 0055378258..09d9760873 100644 --- a/src/region_sphere.h +++ b/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; diff --git a/src/region_union.cpp b/src/region_union.cpp index a0ff63e523..4d339cecca 100644 --- a/src/region_union.cpp +++ b/src/region_union.cpp @@ -46,20 +46,35 @@ 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; - extent_xlo = extent_ylo = extent_zlo = BIG; - extent_xhi = extent_yhi = extent_zhi = -BIG; + bboxflag = 1; + for (int ilist = 0; ilist < nregion; ilist++) + if (regions[list[ilist]]->bboxflag == 0) bboxflag = 0; + if (!interior) bboxflag = 0; - for (int ilist = 0; ilist < nregion; ilist++) { - extent_xlo = MIN(extent_xlo,regions[list[ilist]]->extent_xlo); - extent_ylo = MIN(extent_ylo,regions[list[ilist]]->extent_ylo); - extent_zlo = MIN(extent_zlo,regions[list[ilist]]->extent_zlo); - extent_xhi = MAX(extent_xhi,regions[list[ilist]]->extent_xhi); - extent_yhi = MAX(extent_yhi,regions[list[ilist]]->extent_yhi); - extent_zhi = MAX(extent_zhi,regions[list[ilist]]->extent_zhi); + if (bboxflag) { + extent_xlo = extent_ylo = extent_zlo = BIG; + extent_xhi = extent_yhi = extent_zhi = -BIG; + + for (int ilist = 0; ilist < nregion; ilist++) { + extent_xlo = MIN(extent_xlo,regions[list[ilist]]->extent_xlo); + extent_ylo = MIN(extent_ylo,regions[list[ilist]]->extent_ylo); + extent_zlo = MIN(extent_zlo,regions[list[ilist]]->extent_zlo); + extent_xhi = MAX(extent_xhi,regions[list[ilist]]->extent_xhi); + extent_yhi = MAX(extent_yhi,regions[list[ilist]]->extent_yhi); + extent_zhi = MAX(extent_zhi,regions[list[ilist]]->extent_zhi); + } } + + // possible contacts = sum of possible contacts in all sub-regions + + cmax = 0; + for (int ilist = 0; ilist < nregion; ilist++) + cmax += regions[list[ilist]]->cmax; + contact = new Contact[cmax]; } /* ---------------------------------------------------------------------- */ @@ -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; +} diff --git a/src/region_union.h b/src/region_union.h index 9be04ab193..01ebe79ce5 100644 --- a/src/region_union.h +++ b/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; diff --git a/src/style.h b/src/style.h index d9988be610..e8d150b787 100644 --- a/src/style.h +++ b/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)