From 2c8a7a318a4a1a0786c0ebacabbde82cb6d9c683 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 13 Oct 2016 16:55:53 -0600 Subject: [PATCH] bug fix for fix GCMC w/ fix shake, enhance of fix wall/gran/region with restarting --- doc/README | 3 - doc/src/compute_pe_atom.txt | 7 ++- doc/src/compute_stress_atom.txt | 5 ++ doc/src/fix_wall_gran_region.txt | 24 +++++--- doc/src/read_restart.txt | 29 +++++---- src/GRANULAR/fix_wall_gran.cpp | 5 +- src/GRANULAR/fix_wall_gran_region.cpp | 85 +++++++------------------- src/GRANULAR/fix_wall_gran_region.h | 1 - src/MC/fix_bond_create.cpp | 3 +- src/MC/fix_gcmc.cpp | 15 +++++ src/fix_recenter.cpp | 2 - src/lammps.cpp | 2 +- src/region.cpp | 87 +++++++++++++++++++++------ src/region.h | 9 ++- src/region_intersect.cpp | 60 +++++++++++++++--- src/region_intersect.h | 3 +- src/region_union.cpp | 55 ++++++++++++++--- src/region_union.h | 4 +- 18 files changed, 254 insertions(+), 145 deletions(-) diff --git a/doc/README b/doc/README index fa40e35910..ea0edc0356 100644 --- a/doc/README +++ b/doc/README @@ -91,6 +91,3 @@ This will install virtualenv from the Python Package Index. ---------------- Installing prerequisites for PDF build - - - diff --git a/doc/src/compute_pe_atom.txt b/doc/src/compute_pe_atom.txt index 319983a751..8e2fd17f24 100644 --- a/doc/src/compute_pe_atom.txt +++ b/doc/src/compute_pe_atom.txt @@ -72,9 +72,10 @@ compute peratom all pe/atom compute pe all reduce sum c_peratom thermo_style custom step temp etotal press pe c_pe :pre -NOTE: The per-atom energy does not any Lennard-Jones tail corrections -invoked by the "pair_modify tail yes"_pair_modify.html command, since -those are global contributions to the system energy. +NOTE: The per-atom energy does not include any Lennard-Jones tail +corrections to the energy added by the "pair_modify tail +yes"_pair_modify.html command, since those are contributions to the +global system energy. [Output info:] diff --git a/doc/src/compute_stress_atom.txt b/doc/src/compute_stress_atom.txt index 08e2dd0abe..dcfdee9f85 100644 --- a/doc/src/compute_stress_atom.txt +++ b/doc/src/compute_stress_atom.txt @@ -133,6 +133,11 @@ compute p all reduce sum c_peratom\[1\] c_peratom\[2\] c_peratom\[3\] variable press equal -(c_p\[1\]+c_p\[2\]+c_p\[3\])/(3*vol) thermo_style custom step temp etotal press v_press :pre +NOTE: The per-atom stress does not include any Lennard-Jones tail +corrections to the pressure added by the "pair_modify tail +yes"_pair_modify.html command, since those are contributions to the +global system pressure. + [Output info:] This compute calculates a per-atom array with 6 columns, which can be diff --git a/doc/src/fix_wall_gran_region.txt b/doc/src/fix_wall_gran_region.txt index 478f8b9b23..cfee0d9977 100644 --- a/doc/src/fix_wall_gran_region.txt +++ b/doc/src/fix_wall_gran_region.txt @@ -163,16 +163,20 @@ info on how to re-specify a fix in an input script that reads a restart file, so that the operation of the fix continues in an uninterrupted fashion. -Note that info about region definitions is NOT included in restart -files. So you must re-define your region and if it is a moving -region, define its motion attributes in a way that is consistent with -the simulation that wrote the restart file. In particular, if you -want to change its motion attributes (e.g. its velocity), then you -should insure the postition/orientation of the region at the initial -restart timestep is the same as it was on the timestep the restart -file was written. If this is not possible, then you may need to -ignore info in the restart file by defining a new fix wall/gran/region -command in your restart script (e.g. with a different fix ID). +NOTE: Information about region definitions is NOT included in restart +files, as discussed on the "read_restart"_read_restart.html doc page. +So you must re-define your region and if it is a moving region, define +its motion attributes in a way that is consistent with the simulation +that wrote the restart file. In particular, if you want to change the +region motion attributes (e.g. its velocity), then you should ensure +the postition/orientation of the region at the initial restart +timestep is the same as it was on the timestep the restart file was +written. If this is not possible, you may need to ignore info in the +restart file by defining a new fix wall/gran/region command in your +restart script, e.g. with a different fix ID. Or if you want to keep +the shear history info but discard the region motion information, you +can use the same fix ID for fix wall/gran/region, but assign it a +region with a different region ID. None of the "fix_modify"_fix_modify.html options are relevant to this fix. No global or per-atom quantities are stored by this fix for diff --git a/doc/src/read_restart.txt b/doc/src/read_restart.txt index dbb4e2d499..c79ca4460d 100644 --- a/doc/src/read_restart.txt +++ b/doc/src/read_restart.txt @@ -197,12 +197,13 @@ was used in the input script that wrote the restart file. If a match is found, LAMMPS prints a message indicating that the fix is being re-enabled. If no match is found before the first run or minimization is performed by the new script, the "state" information -for the saved fix is discarded. LAMMPS will also print a list of -fixes for which the information is being discarded. See the doc pages -for individual fixes for info on which ones can be restarted in this -manner. Note that fixes which are created internally by other LAMMPS -commands (computes, fixes, etc) will have style names which are -all-capitalized, and IDs which are generated internally. +for the saved fix is discarded. At the time the discard occurs, +LAMMPS will also print a list of fixes for which the information is +being discarded. See the doc pages for individual fixes for info on +which ones can be restarted in this manner. Note that fixes which are +created internally by other LAMMPS commands (computes, fixes, etc) +will have style names which are all-capitalized, and IDs which are +generated internally. Likewise, the "computes"_fix.html used for a simulation are not stored in the restart file. This means the new input script should specify @@ -219,12 +220,16 @@ described in the previous paragraph, so that the compute can continue its calculations in a consistent manner. NOTE: There are a handful of commands which can be used before or -between runs which require a system initialization. Examples include -the "balance", "displace_atoms", and "delete_atoms" commands. This is -because they may migrate atoms to new processors. Thus they will also -discard unused "state" information from fixes. This means that, if -desired, you must re-specify the relevant fixes and computes (which -create fixes) before those commands are used. +between runs which may require a system initialization. Examples +include the "balance", "displace_atoms", "delete_atoms", "set" (some +options), and "velocity" (some options) commands. This is because +they can migrate atoms to new processors. Thus they will also discard +unused "state" information from fixes. You will know the discard has +occurred because a list of discarded fixes will be printed to the +screen and log file, as explained above. This means that if you wish +to retain that info in a restarted run, you must re-specify the +relevant fixes and computes (which create fixes) before those commands +are used. Some pair styles, like the "granular pair styles"_pair_gran.html, also use a fix to store "state" information that persists from timestep to diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 8ba6a2402d..d53dadb807 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -53,7 +53,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : if (!atom->sphere_flag) error->all(FLERR,"Fix wall/gran requires atom style sphere"); - restart_peratom = 1; create_attribute = 1; // set interaction style @@ -65,8 +64,8 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : //else if (strcmp(arg[3],"bonded/history") == 0) pairstyle = BONDED_HISTORY; else error->all(FLERR,"Invalid fix wall/gran interaction style"); - history = 1; - if (pairstyle == HOOKE) history = 0; + history = restart_peratom = 1; + if (pairstyle == HOOKE) history = restart_peratom = 0; // wall/particle coefficients diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index 96b1f5071f..b3b21660d8 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -101,20 +101,27 @@ 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, " - "resetting its motion",idregion); - error->warning(FLERR,str); - } + 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(); + } +} /* ---------------------------------------------------------------------- */ @@ -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(®ion->nregion,sizeof(int),1,fp); - + int len = 0; + region->length_restart_string(len); + fwrite(&len, sizeof(int),1,fp); region->write_restart(fp); } @@ -505,44 +502,6 @@ 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; + int n = 0; + if (!region->restart(buf,n)) motion_resetflag = 1; } diff --git a/src/GRANULAR/fix_wall_gran_region.h b/src/GRANULAR/fix_wall_gran_region.h index d7bc8be53e..133542edac 100644 --- a/src/GRANULAR/fix_wall_gran_region.h +++ b/src/GRANULAR/fix_wall_gran_region.h @@ -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); }; } diff --git a/src/MC/fix_bond_create.cpp b/src/MC/fix_bond_create.cpp index 59e0e4cf75..81ec67fef1 100644 --- a/src/MC/fix_bond_create.cpp +++ b/src/MC/fix_bond_create.cpp @@ -40,7 +40,8 @@ using namespace FixConst; FixBondCreate::FixBondCreate(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - bondcount(NULL), partner(NULL), finalpartner(NULL), distsq(NULL), created(NULL), copy(NULL), random(NULL) + bondcount(NULL), partner(NULL), finalpartner(NULL), distsq(NULL), + created(NULL), copy(NULL), random(NULL) { if (narg < 8) error->all(FLERR,"Illegal fix bond/create command"); diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index ca93d7f993..8e19928e0a 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -2108,6 +2108,12 @@ double FixGCMC::energy_full() int eflag = 1; int vflag = 0; + // clear forces so they don't accumulate over multiple + // calls within fix gcmc timestep, e.g. for fix shake + + size_t nbytes = sizeof(double) * (atom->nlocal + atom->nghost); + if (nbytes) 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); @@ -2121,9 +2127,18 @@ double FixGCMC::energy_full() if (force->kspace) force->kspace->compute(eflag,vflag); + // unlike Verlet, not performing a reverse_comm() or forces here + // b/c GCMC does not care about forces + // don't think it will mess up energy due to any post_force() fixes + if (modify->n_post_force) modify->post_force(vflag); if (modify->n_end_of_step) modify->end_of_step(); + // NOTE: all fixes with THERMO_ENERGY mask set and which + // operate at pre_force() or post_force() or end_of_step() + // and which user has enable via fix_modify thermo yes, + // will contribute to total MC energy via pe->compute_scalar() + update->eflag_global = update->ntimestep; double total_energy = c_pe->compute_scalar(); diff --git a/src/fix_recenter.cpp b/src/fix_recenter.cpp index e10edcb2d9..499e1a4158 100644 --- a/src/fix_recenter.cpp +++ b/src/fix_recenter.cpp @@ -54,8 +54,6 @@ FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) : global_freq = 1; dynamic_group_allow = 1; -/* ---------------------------------------------------------------------- */ - if (strcmp(arg[3],"NULL") == 0) xflag = 0; else if (strcmp(arg[3],"INIT") == 0) xinitflag = 1; else xcom = force->numeric(FLERR,arg[3]); diff --git a/src/lammps.cpp b/src/lammps.cpp index 4c4e314879..d45a39311b 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -707,7 +707,7 @@ void LAMMPS::init() atom->init(); // atom must come after force and domain // atom deletes extra array // used by fix shear_history::unpack_restart() - // when force->pair->gran_history creates fix ?? + // when force->pair->gran_history creates fix // atom_vec init uses deform_vremap modify->init(); // modify must come after update, force, atom, domain neighbor->init(); // neighbor must come after force, modify diff --git a/src/region.cpp b/src/region.cpp index 08505d1ba5..c23cb031c2 100644 --- a/src/region.cpp +++ b/src/region.cpp @@ -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) && @@ -524,42 +519,95 @@ 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); + 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] += 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); +} + + +/* ---------------------------------------------------------------------- + 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 position/angle + 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 + 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++) { + 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; } /* ---------------------------------------------------------------------- @@ -568,6 +616,5 @@ int Region::restart(char *buf, int n) void Region::reset_vel() { - for (int i = 0; i < size_restart; i++) - prev[i] = 0; + for (int i = 0; i < size_restart; i++) prev[i] = 0; } diff --git a/src/region.h b/src/region.h index 2b920de496..c813772a82 100644 --- a/src/region.h +++ b/src/region.h @@ -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); diff --git a/src/region_intersect.cpp b/src/region_intersect.cpp index 39e8d90d59..02e3769ca6 100644 --- a/src/region_intersect.cpp +++ b/src/region_intersect.cpp @@ -91,11 +91,8 @@ RegIntersect::RegIntersect(LAMMPS *lmp, int narg, char **arg) : // for near contacts and touching contacts cmax = 0; - size_restart = 0; - for (int ilist = 0; ilist < nregion; ilist++){ + for (int ilist = 0; ilist < nregion; ilist++) cmax += regions[list[ilist]]->cmax; - size_restart += regions[list[ilist]]->size_restart; - } contact = new Contact[cmax]; tmax = 0; @@ -274,6 +271,7 @@ void RegIntersect::pretransform() regions[list[ilist]]->pretransform(); } + /* ---------------------------------------------------------------------- get translational/angular velocities of all subregions ------------------------------------------------------------------------- */ @@ -285,6 +283,19 @@ void RegIntersect::set_velocity() regions[list[ilist]]->set_velocity(); } +/* ---------------------------------------------------------------------- + increment length of restart buffer based on region info + used by restart of fix/wall/gran/region +------------------------------------------------------------------------- */ + +void RegIntersect::length_restart_string(int& n) +{ + n += sizeof(int) + strlen(id)+1 + + sizeof(int) + strlen(style)+1 + sizeof(int); + for (int ilist = 0; ilist < nregion; ilist++) + domain->regions[list[ilist]]->length_restart_string(n); + +} /* ---------------------------------------------------------------------- region writes its current position/angle needed by fix/wall/gran/region to compute velocity by differencing scheme @@ -292,8 +303,17 @@ void RegIntersect::set_velocity() void RegIntersect::write_restart(FILE *fp) { - for (int ilist = 0; ilist < nregion; ilist++) + 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); + + for (int ilist = 0; ilist < nregion; ilist++){ domain->regions[list[ilist]]->write_restart(fp); + } } /* ---------------------------------------------------------------------- @@ -301,11 +321,33 @@ void RegIntersect::write_restart(FILE *fp) needed by fix/wall/gran/region to compute velocity by differencing scheme ------------------------------------------------------------------------- */ -int RegIntersect::restart(char *buf, int n) +int RegIntersect::restart(char *buf, int& n) { - for (int ilist = 0; ilist < nregion; ilist++) - n = domain->regions[list[ilist]]->restart(buf, n); - return 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; + + for (int ilist = 0; ilist < nregion; ilist++){ + if (!domain->regions[list[ilist]]->restart(buf, n)){ + return 0; + } + } + return 1; } /* ---------------------------------------------------------------------- diff --git a/src/region_intersect.h b/src/region_intersect.h index 546baac265..2c4a2f2a0b 100644 --- a/src/region_intersect.h +++ b/src/region_intersect.h @@ -35,8 +35,9 @@ class RegIntersect : public Region { void shape_update(); void pretransform(); void set_velocity(); + void length_restart_string(int&); void write_restart(FILE *); - int restart(char *, int); + int restart(char *, int&); void reset_vel(); private: diff --git a/src/region_union.cpp b/src/region_union.cpp index b516af98d0..c506c811c3 100644 --- a/src/region_union.cpp +++ b/src/region_union.cpp @@ -85,11 +85,8 @@ RegUnion::RegUnion(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg) // for near contacts and touching contacts cmax = 0; - size_restart = 0; - for (int ilist = 0; ilist < nregion; ilist++){ + for (int ilist = 0; ilist < nregion; ilist++) cmax += regions[list[ilist]]->cmax; - size_restart += regions[list[ilist]]->size_restart; - } contact = new Contact[cmax]; tmax = 0; @@ -279,6 +276,20 @@ void RegUnion::set_velocity() for (int ilist = 0; ilist < nregion; ilist++) regions[list[ilist]]->set_velocity(); } + +/* ---------------------------------------------------------------------- + increment length of restart buffer based on region info + used by restart of fix/wall/gran/region +------------------------------------------------------------------------- */ + +void RegUnion::length_restart_string(int& n) +{ + n += sizeof(int) + strlen(id)+1 + + sizeof(int) + strlen(style)+1 + sizeof(int); + for (int ilist = 0; ilist < nregion; ilist++) + domain->regions[list[ilist]]->length_restart_string(n); + +} /* ---------------------------------------------------------------------- region writes its current position/angle needed by fix/wall/gran/region to compute velocity by differencing scheme @@ -286,8 +297,15 @@ void RegUnion::set_velocity() void RegUnion::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); for (int ilist = 0; ilist < nregion; ilist++) - domain->regions[list[ilist]]->write_restart(fp); + domain->regions[list[ilist]]->write_restart(fp); } /* ---------------------------------------------------------------------- @@ -295,11 +313,30 @@ void RegUnion::write_restart(FILE *fp) needed by fix/wall/gran/region to compute velocity by differencing scheme ------------------------------------------------------------------------- */ -int RegUnion::restart(char *buf, int n) +int RegUnion::restart(char *buf, int &n) { - for (int ilist = 0; ilist < nregion; ilist++) - n = domain->regions[list[ilist]]->restart(buf, n); - return 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; + + for (int ilist = 0; ilist < nregion; ilist++){ + if (!domain->regions[list[ilist]]->restart(buf, n)) return 0; + } + return 1; } /* ---------------------------------------------------------------------- diff --git a/src/region_union.h b/src/region_union.h index 091688d66e..cf481c324c 100644 --- a/src/region_union.h +++ b/src/region_union.h @@ -35,10 +35,10 @@ class RegUnion : public Region { void shape_update(); void pretransform(); void set_velocity(); + void length_restart_string(int&); void write_restart(FILE *); - int restart(char *, int); + int restart(char *, int&); void reset_vel(); - private: char **idsub; };