Commit 0f1c56d0 authored by sjplimp's avatar sjplimp
Browse files 

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

src/GRANULAR/fix_wall_gran_region.cpp

+22 −63
Original line number Diff line number Diff line
@@ -101,21 +101,28 @@ void FixWallGranRegion::init() 
    error->all(FLERR,"Region ID for fix wall/gran/region does not exist");

  region = domain->regions[iregion];



  // region displacement and orientation theta at previous step

  // check if region properties changed between runs

  // reset if restart info was inconsistent



  if (motion_resetflag) {

    if (comm->me == 0) {

      char str[128];

      sprintf(str,"Properties for region %s do not match restart file, "

  if (strcmp(idregion,region->id) != 0 ||

      strcmp(region_style,region->style) != 0 ||

      nregion != region->nregion) {

    char str[256];

    sprintf(str,"Region properties for region %s changed between runs, "

            "resetting its motion",idregion);

    error->warning(FLERR,str);

    region->reset_vel();

  }



  if (motion_resetflag){    

    char str[256];

    sprintf(str,"Region properties for region %s are inconsistent "

            "with restart file, resetting its motion",idregion);

    error->warning(FLERR,str);

    region->reset_vel();

  }

}





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



void FixWallGranRegion::post_force(int vflag)
@@ -483,19 +490,9 @@ int FixWallGranRegion::size_restart(int nlocal) 
void FixWallGranRegion::write_restart(FILE *fp)

{

  if (comm->me) return;

  int size_id_str = (strlen(region->id) + 1) * sizeof(char);

  int size_style_str = (strlen(region->style) + 1) * sizeof(char);

  int size_tot = sizeof(int) + size_id_str + 

    sizeof(int) + size_style_str + sizeof(int) + 

    region->size_restart*sizeof(double);



  fwrite(&size_tot,sizeof(int),1,fp);

  fwrite(&size_id_str,sizeof(int),1,fp);

  fwrite(region->id,sizeof(char),size_id_str,fp);

  fwrite(&size_style_str,sizeof(int),1,fp);

  fwrite(region->style,sizeof(char),size_style_str,fp);  

  fwrite(&region->nregion,sizeof(int),1,fp);



  int len = 0;

  region->length_restart_string(len);

  fwrite(&len, sizeof(int),1,fp);

  region->write_restart(fp);

}
@@ -506,43 +503,5 @@ void FixWallGranRegion::write_restart(FILE *fp) 
void FixWallGranRegion::restart(char *buf)

{

  int n = 0;

  int size_id_str = buf[n];

  n += sizeof(int);

  char *region_id_restart = new char[size_id_str];

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

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

  }

  

  int size_style_str = buf[n];

  n += sizeof(int);

  char *region_style_restart = new char[size_style_str];

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

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

  

  int nregion_restart = buf[n];

  n += sizeof(int);



  if (check_consistent_region(region,region_id_restart,

                              region_style_restart,nregion_restart))

    region->restart(buf,n);

  else motion_resetflag = 1;



  delete [] region_id_restart;

  delete [] region_style_restart;

}





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

   check that region id/style/number of sub-regions are consistent

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



int FixWallGranRegion::check_consistent_region(Region *region, 

                                               char* region_id, 

                                               char* region_style, int nregion)

{

  if (strcmp(region_id, region->id) != 0 ||

      strcmp(region_style, region->style) != 0 ||

      nregion != region->nregion)

    return 0;       

  return 1;

  if (!region->restart(buf,n)) motion_resetflag = 1;  

}

src/GRANULAR/fix_wall_gran_region.h

+0 −1
Original line number Diff line number Diff line
@@ -62,7 +62,6 @@ class FixWallGranRegion : public FixWallGran { 
                         //    vel info is to be reset



  void update_contacts(int, int);

  int check_consistent_region(Region *, char*, char*, int);

};



}

src/MC/fix_gcmc.cpp

+0 −5
Original line number Diff line number Diff line
@@ -2108,11 +2108,6 @@ double FixGCMC::energy_full() 
  int eflag = 1;

  int vflag = 0;



  // clear forces because they are used by fix shake 

  

  int nbytes = sizeof(double) * (atom->nlocal + atom->nghost);

  memset(&atom->f[0][0],0,3*nbytes);



  if (modify->n_pre_force) modify->pre_force(vflag);



  if (force->pair) force->pair->compute(eflag,vflag);

src/region.cpp

+70 −27
Original line number Diff line number Diff line
@@ -42,10 +42,6 @@ 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;
@@ -398,7 +394,6 @@ void Region::options(int narg, char **arg) 
    else error->all(FLERR,"Illegal region command");

  }





  // error check



  if ((moveflag || rotateflag) &&
@@ -514,58 +509,106 @@ void Region::set_velocity() 
}



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

   compute velocity of wall for given contact

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

   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

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



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



  if (moveflag){

    vwall[0] = v[0];

    vwall[1] = v[1];

    vwall[2] = v[2];

  }

  if (rotateflag){

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

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

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

    vwall[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

   Increment length of restart buffer based on region info

   Used by restart of fix/wall/gran/region

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

void Region::length_restart_string(int& n)

{

  n += sizeof(int) + strlen(id)+1 + 

    sizeof(int) + strlen(style)+1 + sizeof(int)+

    size_restart*sizeof(double);    

}



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

   Region writes its current style, id, number of sub-regions and position/angle

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

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

void Region::write_restart(FILE *fp)

{

  int sizeid = (strlen(id)+1);

  int sizestyle = (strlen(style)+1);

  fwrite(&sizeid, sizeof(int), 1, fp);

  fwrite(id, 1, sizeid, fp);

  fwrite(&sizestyle, sizeof(int), 1, fp);

  fwrite(style, 1, sizestyle, fp);  

  fwrite(&nregion,sizeof(int),1,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

   Region reads style, id, number of sub-regions from restart file. If they

    match current region, also read previous position/angle

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

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



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

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

{

  int sizeid = buf[n];

  n += sizeof(int);

  char *restart_id = new char[sizeid];

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

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

  if (strcmp(restart_id,id) != 0) return 0;



  int sizestyle = buf[n];

  n += sizeof(int);

  char *restart_style = new char[sizestyle];

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

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

  if (strcmp(restart_style,style) != 0) return 0;    



  int restart_nreg = buf[n];

  n += sizeof(int);

  if (restart_nreg != nregion) return 0;



  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;

  delete [] restart_id;

  delete [] restart_style;

  return 1;

}



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

   set prev vector to zero

   Set prev vector to zero

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



void Region::reset_vel()

{

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

src/region.h

+4 −5
Original line number Diff line number Diff line
@@ -59,8 +59,7 @@ class Region : protected Pointers { 
  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 xcenter[3];          // translated/rotated center of cylinder/sphere (only used if varshape)

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

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

                              //  at previous time step
@@ -84,7 +83,8 @@ class Region : protected Pointers { 
  virtual void set_velocity();

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

  virtual void write_restart(FILE *);

  virtual int restart(char *, int);

  virtual int restart(char *, int&);

  virtual void length_restart_string(int&);

  virtual void reset_vel();



  // implemented by each region, not called by other classes
@@ -106,12 +106,11 @@ class Region : protected Pointers { 
  void options(int, char **);

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

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

  double point[3],runit[3];



 private:

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

  int xvar,yvar,zvar,tvar;

  double axis[3];

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



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

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

CRA Git | Maintained and supported by SUSTech CRA and CCSE