Commit 65585e69 authored by sjplimp's avatar sjplimp
Browse files 

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

src/region.cpp

+175 −5
Original line number Diff line number Diff line
@@ -20,6 +20,7 @@ 
#include "lattice.h"

#include "input.h"

#include "variable.h"

#include "math_extra.h"

#include "error.h"

#include "force.h"
@@ -41,7 +42,15 @@ Region::Region(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) 
  xstr = ystr = zstr = tstr = NULL;

  dx = dy = dz = 0.0;



  // used by set_velocity() even if no rotation specified



  point[0] = point[1] = point[2] = 0.0;



  size_restart = 5;

  reset_vel();

  copymode = 0;

  list = NULL;

  nregion = 1;

}



/* ---------------------------------------------------------------------- */
@@ -87,6 +96,7 @@ void Region::init() 
    if (!input->variable->equalstyle(tvar))

      error->all(FLERR,"Variable for region is not equal style");

  }

  vel_timestep = -1;

}



/* ----------------------------------------------------------------------
@@ -129,6 +139,7 @@ void Region::prematch() 
int Region::match(double x, double y, double z)

{

  if (dynamic) inverse_transform(x,y,z);

  if (openflag) return 1;

  return !(inside(x,y,z) ^ interior);

}
@@ -170,8 +181,15 @@ int Region::surface(double x, double y, double z, double cutoff) 
  xnear[1] = y;

  xnear[2] = z;



  if (!openflag) {

    if (interior) ncontact = surface_interior(xnear,cutoff);

    else ncontact = surface_exterior(xnear,cutoff);

  }

  else{

    // one of surface_int/ext() will return 0

    // so no need to worry about offset of contact indices

    ncontact = surface_exterior(xnear,cutoff) + surface_interior(xnear,cutoff);

  }



  if (rotateflag && ncontact) {

    for (int i = 0; i < ncontact; i++) {
@@ -200,6 +218,7 @@ void Region::add_contact(int n, double *x, double xp, double yp, double zp) 
  double dely = x[1] - yp;

  double delz = x[2] - zp;

  contact[n].r = sqrt(delx*delx + dely*dely + delz*delz);

  contact[n].radius = 0;

  contact[n].delx = delx;

  contact[n].dely = dely;

  contact[n].delz = delz;
@@ -307,6 +326,9 @@ void Region::options(int narg, char **arg) 
  scaleflag = 1;

  moveflag = rotateflag = 0;



  openflag = 0;

  for (int i  = 0; i < 6; i++) open_faces[i] = 0;



  int iarg = 0;

  while (iarg < narg) {

    if (strcmp(arg[iarg],"units") == 0) {
@@ -363,9 +385,20 @@ void Region::options(int narg, char **arg) 
      axis[2] = force->numeric(FLERR,arg[iarg+7]);

      rotateflag = 1;

      iarg += 8;

    } else error->all(FLERR,"Illegal region command");



    } else if (strcmp(arg[iarg],"open") == 0) {

      if (iarg+2 > narg) error->all(FLERR,"Illegal region command");

      int iface = force->inumeric(FLERR,arg[iarg+1]);

      if (iface < 1 || iface > 6) error->all(FLERR,"Illegal region command");

      // additional checks on valid face index are done by region classes

      open_faces[iface-1] = 1;

      openflag = 1;

      iarg += 2;

    }

    else error->all(FLERR,"Illegal region command");

  }





  // error check



  if ((moveflag || rotateflag) &&
@@ -401,3 +434,140 @@ void Region::options(int narg, char **arg) 
  if (moveflag || rotateflag) dynamic = 1;

  else dynamic = 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 Region::point_on_line_segment(double *a, double *b,

                                   double *c, double *d)

{

  double ba[3],ca[3];



  MathExtra::sub3(b,a,ba);

  MathExtra::sub3(c,a,ca);

  double t = MathExtra::dot3(ca,ba) / MathExtra::dot3(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];

  }

}



/* ----------------------------------------------------------------------

   infer translational and angular velocity of region

   necessary b/c motion variables are for displacement & theta

     there is no analytic formula for v & omega

   prev[4] contains values of dx,dy,dz,theta at previous step

     used for difference, then updated to current step values

   dt is time elapsed since previous step

   rpoint = point updated by current displacement

   called by fix wall/gran/region every timestep

------------------------------------------------------------------------- */



void Region::set_velocity()

{

  if (vel_timestep == update->ntimestep) return;

  vel_timestep = update->ntimestep;

  if (moveflag) {

    if (update->ntimestep > 0) {

      v[0] = (dx - prev[0])/update->dt;

      v[1] = (dy - prev[1])/update->dt;

      v[2] = (dz - prev[2])/update->dt;

    }

    else v[0] = v[1] = v[2] = 0.0;

    prev[0] = dx;

    prev[1] = dy;

    prev[2] = dz;    

  }



  if (rotateflag) {

    rpoint[0] = point[0] + dx;

    rpoint[1] = point[1] + dy;

    rpoint[2] = point[2] + dz;

    if (update->ntimestep > 0) {

      double angvel = (theta-prev[3]) / update->dt;

      omega[0] = angvel*axis[0];

      omega[1] = angvel*axis[1];

      omega[2] = angvel*axis[2];

    }

    else omega[0] = omega[1] = omega[2] = 0.0;

    prev[3] = theta;

  }



  if (varshape){

    set_velocity_shape();

  }

}



/* ----------------------------------------------------------------------

   compute velocity of wall for given contact

   since contacts only store delx/y/z, need to pass particle coords

     to compute contact point

   called by fix/wall/gran/region every contact every timestep

------------------------------------------------------------------------- */



void Region::velocity_contact(double *vwall, double *x, int ic)

{

  Contact c = contact[ic];

  double xc[3];

  xc[0] = x[0] - contact[ic].delx;

  xc[1] = x[1] - contact[ic].dely;

  xc[2] = x[2] - contact[ic].delz;

  

  vwall[0] = v[0] + omega[1]*(xc[2] - rpoint[2]) - omega[2]*(xc[1] - rpoint[1]);

  vwall[1] = v[1] + omega[2]*(xc[0] - rpoint[0]) - omega[0]*(xc[2] - rpoint[2]);

  vwall[2] = v[2] + omega[0]*(xc[1] - rpoint[1]) - omega[1]*(xc[0] - rpoint[0]);



  if (varshape && contact[ic].varflag) velocity_contact_shape(vwall,xc);

}



/* ----------------------------------------------------------------------

   region writes its current position/angle

   needed by fix/wall/gran/region to compute velocity by differencing scheme

------------------------------------------------------------------------- */



void Region::write_restart(FILE *fp)

{

  fwrite(prev, sizeof(double), size_restart, fp);  

}



/* ----------------------------------------------------------------------

   region reads its previous position/angle

   needed by fix/wall/gran/region to compute velocity by differencing scheme

------------------------------------------------------------------------- */



int Region::restart(char *buf, int n)

{

  char *rlist = new char[size_restart*sizeof(double)];  

  for (int i = 0; i < size_restart*sizeof(double); i++)

    rlist[i] = buf[n++]; 

  for (int i = 0; i < size_restart; i++) {

    prev[i] = ((double *)rlist)[i];

  }

  delete [] rlist;

  return n;

}



/* ----------------------------------------------------------------------

   set prev vector to zero

------------------------------------------------------------------------- */



void Region::reset_vel()

{

  for (int i = 0; i < size_restart; i++)

    prev[i] = 0;

}

src/region.h

+40 −10
Original line number Diff line number Diff line
@@ -31,17 +31,44 @@ class Region : protected Pointers { 
  int bboxflag;                     // 1 if bounding box is computable

  int varshape;                     // 1 if region shape changes over time

  int dynamic;                      // 1 if position/orient changes over time

  int moveflag,rotateflag;          // 1 if position/orientation changes

  int openflag;			    // 1 if any face is open

  int open_faces[6];		    // flags for which faces are open



  int copymode;                     // 1 if copy of original class



  // contact = particle near region surface

  // contact = particle near region surface (for soft interactions)

  // touch = particle touching region surface (for granular interactions)



  struct Contact {

    double r;                 // distance between particle & surf, r > 0.0

    double delx,dely,delz;    // vector from surface pt to particle

    double radius;            // curvature of region at contact point

    int iwall;		      // unique id of wall for storing shear history

    int varflag;              // 1 if wall can be variable-controlled

  };

  Contact *contact;           // list of contacts

  int cmax;                   // max # of contacts possible with region

  int tmax;		      // max # of touching contacts possible



  // motion attributes of region

  // public so can be accessed by other classes



  double dx,dy,dz,theta;      // current displacement and orientation

  double v[3];	 	      // translational velocity

  double rpoint[3];	      // current origin of rotation axis

  double omega[3];	      // angular velocity

  double rprev;               // speed of time-dependent radius, if applicable

  double xcenter[3];          // translated/rotated center of cylinder/sphere

                              //   only used with varshape

  double prev[5];             // stores displacement (X3), angle and if 

                              //  necessary, region variable size (e.g. radius)

                              //  at previous time step

  int vel_timestep;           // store timestep at which set_velocity was called

                              //   prevents multiple fix/wall/gran/region calls

  int nregion;                // For union and intersect

  int size_restart;

  int *list;



  Region(class LAMMPS *, int, char **);

  virtual ~Region();
@@ -54,6 +81,12 @@ class Region : protected Pointers { 
  int match(double, double, double);

  int surface(double, double, double, double);



  virtual void set_velocity();

  void velocity_contact(double *, double *, int);

  virtual void write_restart(FILE *);

  virtual int restart(char *, int);

  virtual void reset_vel();



  // implemented by each region, not called by other classes



  virtual int inside(double, double, double) = 0;
@@ -61,6 +94,8 @@ class Region : protected Pointers { 
  virtual int surface_exterior(double *, double) = 0;

  virtual void shape_update() {}

  virtual void pretransform();

  virtual void set_velocity_shape() {}

  virtual void velocity_contact_shape(double*, double*) {}



  // Kokkos function, implemented by each Kokkos region
@@ -69,15 +104,14 @@ class Region : protected Pointers { 
 protected:

  void add_contact(int, double *, double, double, double);

  void options(int, char **);

  void point_on_line_segment(double *, double *, double *, double *);

  void forward_transform(double &, double &, double &);



  int moveflag,rotateflag;         // 1 if position/orientation changes



  double point[3],axis[3],runit[3];

 private:

  char *xstr,*ystr,*zstr,*tstr;

  int xvar,yvar,zvar,tvar;

  double dx,dy,dz,theta;

  double axis[3],point[3],runit[3];



  void forward_transform(double &, double &, double &);

  void inverse_transform(double &, double &, double &);

  void rotate(double &, double &, double &, double);

};
@@ -100,10 +134,6 @@ E: Variable for region is not equal style 


Self-explanatory.



E: Can only use Kokkos supported regions with Kokkos package



Self-explanatory.



E: Illegal ... command



Self-explanatory.  Check the input script syntax and compare to the

src/region_block.cpp

+224 −16
Original line number Diff line number Diff line
@@ -14,9 +14,10 @@ 
#include <stdlib.h>

#include <string.h>

#include "region_block.h"

#include "force.h"

#include "domain.h"

#include "math_extra.h"

#include "error.h"

#include "force.h"



using namespace LAMMPS_NS;
@@ -94,9 +95,91 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) 
  } else bboxflag = 0;



  // particle could be close to all 6 planes

  // particle can only touch 3 planes



  cmax = 6;

  contact = new Contact[cmax];

  if (interior) tmax = 3;

  else tmax = 1;



  // open face data structs



  face[0][0] = -1.0;

  face[0][1] = 0.0;

  face[0][2] = 0.0;

  face[1][0] = 1.0;

  face[1][1] = 0.0;

  face[1][2] = 0.0;

  face[2][0] = 0.0;

  face[2][1] = -1.0;

  face[2][2] = 0.0;

  face[3][0] = 0.0;

  face[3][1] = 1.0;

  face[3][2] = 0.0;

  face[4][0] = 0.0;

  face[4][1] = 0.0;

  face[4][2] = -1.0;

  face[5][0] = 0.0;

  face[5][1] = 0.0;

  face[5][2] = 1.0;



  // face[0]



  corners[0][0][0] = xlo;

  corners[0][0][1] = ylo;

  corners[0][0][2] = zlo;

  corners[0][1][0] = xlo;

  corners[0][1][1] = ylo;

  corners[0][1][2] = zhi;

  corners[0][2][0] = xlo;

  corners[0][2][1] = yhi;

  corners[0][2][2] = zhi;

  corners[0][3][0] = xlo;

  corners[0][3][1] = yhi;

  corners[0][3][2] = zlo;



  // face[1]



  corners[1][0][0] = xhi;

  corners[1][0][1] = ylo;

  corners[1][0][2] = zlo;

  corners[1][1][0] = xhi;

  corners[1][1][1] = ylo;

  corners[1][1][2] = zhi;

  corners[1][2][0] = xhi;

  corners[1][2][1] = yhi;

  corners[1][2][2] = zhi;

  corners[1][3][0] = xhi;

  corners[1][3][1] = yhi;

  corners[1][3][2] = zlo;



  // face[2]



  MathExtra::copy3(corners[2][0],corners[0][0]);

  MathExtra::copy3(corners[2][1],corners[1][0]);

  MathExtra::copy3(corners[2][2],corners[1][1]);

  MathExtra::copy3(corners[2][3],corners[0][1]);



  // face[3]



  MathExtra::copy3(corners[3][0],corners[0][3]);

  MathExtra::copy3(corners[3][1],corners[0][2]);

  MathExtra::copy3(corners[3][2],corners[1][2]);

  MathExtra::copy3(corners[3][3],corners[1][3]);



  // face[4]



  MathExtra::copy3(corners[4][0],corners[0][0]);

  MathExtra::copy3(corners[4][1],corners[0][3]);

  MathExtra::copy3(corners[4][2],corners[1][3]);

  MathExtra::copy3(corners[4][3],corners[1][0]);



  // face[5]



  MathExtra::copy3(corners[5][0],corners[0][1]);

  MathExtra::copy3(corners[5][1],corners[1][1]);

  MathExtra::copy3(corners[5][2],corners[1][2]);

  MathExtra::copy3(corners[5][3],corners[0][2]);

}



/* ---------------------------------------------------------------------- */
@@ -141,47 +224,59 @@ int RegBlock::surface_interior(double *x, double cutoff) 
  int n = 0;



  delta = x[0] - xlo;

  if (delta < cutoff) {

  if (delta < cutoff && !open_faces[0]) {

    contact[n].r = delta;

    contact[n].delx = delta;

    contact[n].dely = contact[n].delz = 0.0;

    contact[n].radius = 0;

    contact[n].iwall = 0;

    n++;

  }

  delta = xhi - x[0];

  if (delta < cutoff) {

  if (delta < cutoff && !open_faces[1]) {

    contact[n].r = delta;

    contact[n].delx = -delta;

    contact[n].dely = contact[n].delz = 0.0;

    contact[n].radius = 0;

    contact[n].iwall = 1;

    n++;

  }



  delta = x[1] - ylo;

  if (delta < cutoff) {

  if (delta < cutoff && !open_faces[2]) {

    contact[n].r = delta;

    contact[n].dely = delta;

    contact[n].delx = contact[n].delz = 0.0;

    contact[n].radius = 0;

    contact[n].iwall = 2;

    n++;

  }

  delta = yhi - x[1];

  if (delta < cutoff) {

  if (delta < cutoff && !open_faces[3]) {

    contact[n].r = delta;

    contact[n].dely = -delta;

    contact[n].delx = contact[n].delz = 0.0;

    contact[n].radius = 0;

    contact[n].iwall = 3;

    n++;

  }



  delta = x[2] - zlo;

  if (delta < cutoff) {

  if (delta < cutoff && !open_faces[4]) {

    contact[n].r = delta;

    contact[n].delz = delta;

    contact[n].delx = contact[n].dely = 0.0;

    contact[n].radius = 0;

    contact[n].iwall = 4;

    n++;

  }

  delta = zhi - x[2];

  if (delta < cutoff) {

  if (delta < cutoff && !open_faces[5]) {

    contact[n].r = delta;

    contact[n].delz = -delta;

    contact[n].delx = contact[n].dely = 0.0;

    contact[n].radius = 0;

    contact[n].iwall = 5;

    n++;

  }
@@ -197,6 +292,7 @@ int RegBlock::surface_interior(double *x, double cutoff) 
int RegBlock::surface_exterior(double *x, double cutoff)

{

  double xp,yp,zp;

  double xc,yc,zc,dist,mindist;



  // x is far enough from block that there is no contact

  // x is interior to block
@@ -212,6 +308,7 @@ int RegBlock::surface_exterior(double *x, double cutoff) 
  //            could be edge or corner pt of block

  // do not add contact point if r >= cutoff



  if (!openflag){

    if (x[0] < xlo) xp = xlo;

    else if (x[0] > xhi) xp = xhi;

    else xp = x[0];
@@ -221,8 +318,119 @@ int RegBlock::surface_exterior(double *x, double cutoff) 
    if (x[2] < zlo) zp = zlo;

    else if (x[2] > zhi) zp = zhi;

    else zp = x[2];

  }

  else{

    mindist = BIG;

    for (int i = 0; i < 6; i++){

      if (open_faces[i]) continue;

      dist = find_closest_point(i,x,xc,yc,zc);

      if (dist < mindist){

        xp = xc;

        yp = yc;

        zp = zc;

        mindist = dist;

      }

    }

  }



  add_contact(0,x,xp,yp,zp);

  contact[0].iwall = 0;

  if (contact[0].r < cutoff) return 1;

  return 0;

}



/*------------------------------------------------------------------------

  return distance to closest point on surface I of block region

  store closest point in xc,yc,zc

--------------------------------------------------------------------------*/



double RegBlock::find_closest_point(int i, double *x, 

                                    double &xc, double &yc, double &zc)

{

  double dot,d2,d2min;

  double xr[3],xproj[3],p[3];



  xr[0] = x[0] - corners[i][0][0];

  xr[1] = x[1] - corners[i][0][1];

  xr[2] = x[2] - corners[i][0][2];

  dot = face[i][0]*xr[0] + face[i][1]*xr[1] + face[i][2]*xr[2];

  xproj[0] = xr[0] - dot*face[i][0];

  xproj[1] = xr[1] - dot*face[i][1];

  xproj[2] = xr[2] - dot*face[i][2];



  d2min = BIG;



  // check if point projects inside of face



  if (inside_face(xproj, i)){

    d2 = d2min = dot*dot;

    xc = xproj[0] + corners[i][0][0];

    yc = xproj[1] + corners[i][0][1];

    zc = xproj[2] + corners[i][0][2];



 // check each edge



  } else {

    point_on_line_segment(corners[i][0],corners[i][1],x,p);

    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + 

      (p[2]-x[2])*(p[2]-x[2]);

    if (d2 < d2min) {

      d2min = d2;

      xc = p[0];

      yc = p[1];

      zc = p[2];

    }



    point_on_line_segment(corners[i][1],corners[i][2],x,p);

    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + 

      (p[2]-x[2])*(p[2]-x[2]);

    if (d2 < d2min) {

      d2min = d2;

      xc = p[0];

      yc = p[1];

      zc = p[2];

    }



    point_on_line_segment(corners[i][2],corners[i][3],x,p);

    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + 

      (p[2]-x[2])*(p[2]-x[2]);

    if (d2 < d2min) {

      d2min = d2;

      xc = p[0];

      yc = p[1];

      zc = p[2];

    }



    point_on_line_segment(corners[i][3],corners[i][4],x,p);

    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + 

      (p[2]-x[2])*(p[2]-x[2]);

    if (d2 < d2min) {

      d2min = d2;

      xc = p[0];

      yc = p[1];

      zc = p[2];

    }

  }



  return d2min;

}



/*------------------------------------------------------------------------

  determine if projected point is inside given face of the block

--------------------------------------------------------------------------*/



int RegBlock::inside_face(double *xproj, int iface)

{

  if (iface < 2) {

    if (xproj[1] > 0 && (xproj[1] < yhi-ylo) &&

        xproj[2] > 0 && (xproj[2] < zhi-zlo)) return 1;

  } else if (iface < 4) {

    if (xproj[0] > 0 && (xproj[0] < (xhi-xlo)) &&

        xproj[2] > 0 && (xproj[2] < (zhi-zlo))) return 1;

  } else {

    if (xproj[0] > 0 && xproj[0] < (xhi-xlo) &&

        xproj[1] > 0 && xproj[1] < (yhi-ylo)) return 1;

  }



  return 0;

}

src/region_block.h

+5 −0
Original line number Diff line number Diff line
@@ -36,6 +36,11 @@ class RegBlock : public Region { 


 protected:

  double xlo,xhi,ylo,yhi,zlo,zhi;

  double corners[6][4][3];

  double face[6][3];



  double find_closest_point(int, double *, double &, double &, double &);

  int inside_face(double *, int);

};



}

src/region_cone.cpp

+118 −89

File changed.

Preview size limit exceeded, changes collapsed.

CRA Git | Maintained and supported by SUSTech CRA and CCSE