diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index 0dbbf8e1bf..597c3faae6 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -2684,11 +2684,9 @@ void *FixShake::extract(const char *str, int &dim) return NULL; } - - /* ---------------------------------------------------------------------- - Add coordinate constraining forces; this method is called - at the end of a timestep. + add coordinate constraining forces + this method is called at the end of a timestep ------------------------------------------------------------------------- */ void FixShake::shake_end_of_step(int vflag) { @@ -2715,20 +2713,15 @@ void FixShake::shake_end_of_step(int vflag) { } } - - /* ---------------------------------------------------------------------- wrapper method for end_of_step fixes which modify velocities ------------------------------------------------------------------------- */ -void FixShake::correct_velocities() { - -} - +void FixShake::correct_velocities() {} /* ---------------------------------------------------------------------- - Calculate constraining forces based on the current configuration - and change coordinates. + calculate constraining forces based on the current configuration + change coordinates ------------------------------------------------------------------------- */ void FixShake::correct_coordinates(int vflag) { diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index bad1e56e43..afd3868cb4 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -50,13 +50,10 @@ class FixShake : public Fix { virtual int pack_forward_comm(int, int *, double *, int, int *); virtual void unpack_forward_comm(int, int, double *); - - virtual void shake_end_of_step(int vflag); virtual void correct_coordinates(int vflag); virtual void correct_velocities(); - int dof(int); virtual void reset_dt(); void *extract(const char *, int &); @@ -86,7 +83,7 @@ class FixShake : public Fix { double *step_respa; double **x,**v,**f; // local ptrs to atom class quantities - double **ftmp, **vtmp; // pointers to temporary arrays for forces and velocities + double **ftmp,**vtmp; // pointers to temporary arrays for f,v double *mass,*rmass; int *type; diff --git a/src/fix_shear_history.cpp b/src/fix_shear_history.cpp index a11def70ba..060466f8d8 100644 --- a/src/fix_shear_history.cpp +++ b/src/fix_shear_history.cpp @@ -45,6 +45,9 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : if (strcmp(id,"LINE_SHEAR_HISTORY") == 0) onesided = 1; if (strcmp(id,"TRI_SHEAR_HISTORY") == 0) onesided = 1; + surfglobal = 0; + if (strcmp(id,"SURFACE_GLOBAL_SHEAR_HISTORY") == 0) surfglobal = 1; + if (newton_pair) comm_reverse = 1; // just for single npartner value // variable-size history communicated via // reverse_comm_fix_variable() @@ -106,10 +109,15 @@ void FixShearHistory::init() { if (atom->tag_enable == 0) error->all(FLERR, - "Pair style granular with history requires atoms have IDs"); + "Granular interactions with history require atoms have IDs"); - int dim; - computeflag = (int *) pair->extract("computeflag",dim); + // extract pointers to neigh lists in FixSurfaceGlobal + + if (surfglobal) { + int dim; + fsg_list = (NeighList *) fsg->extract("list",dim); + fsg_listgranhistory = (NeighList *) fsg->extract("listgranhistory",dim); + } allocate_pages(); } @@ -161,6 +169,7 @@ void FixShearHistory::allocate_pages() void FixShearHistory::pre_exchange() { if (onesided) pre_exchange_onesided(); + else if (surfglobal) pre_exchange_surf(); else if (newton_pair) pre_exchange_newton(); else pre_exchange_no_newton(); } @@ -265,6 +274,103 @@ void FixShearHistory::pre_exchange_onesided() for (i = nlocal_neigh; i < nlocal; i++) npartner[i] = 0; } +/* ---------------------------------------------------------------------- + version for sphere contact with global line/tri particles + neighbor list has I = sphere, J = line/tri + only store history info with spheres +------------------------------------------------------------------------- */ + +void FixShearHistory::pre_exchange_surf() +{ + int i,j,ii,jj,m,n,inum,jnum; + int *ilist,*jlist,*numneigh,**firstneigh; + int *touch,**firsttouch; + double *shear,*allshear,**firstshear; + + // NOTE: all operations until very end are with nlocal_neigh <= current nlocal + // b/c previous neigh list was built with nlocal_neigh + // nlocal can be larger if other fixes added atoms at this pre_exchange() + + // zero npartner for owned atoms + // clear 2 page data structures + + for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0; + + ipage->reset(); + dpage->reset(); + + // 1st loop over neighbor list, I = sphere, J = global tri + // only calculate npartner for owned spheres + + inum = fsg_list->inum; + ilist = fsg_list->ilist; + numneigh = fsg_list->numneigh; + firstneigh = fsg_list->firstneigh; + firsttouch = fsg_listgranhistory->firstneigh; + firstshear = fsg_listgranhistory->firstdouble; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + touch = firsttouch[i]; + + for (jj = 0; jj < jnum; jj++) { + if (touch[jj]) npartner[i]++; + } + } + + // get page chunks to store atom IDs and shear history for owned atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + n = npartner[i]; + partner[i] = ipage->get(n); + shearpartner[i] = dpage->get(n); + if (partner[i] == NULL || shearpartner[i] == NULL) + error->one(FLERR,"Shear history overflow, boost neigh_modify one"); + } + + // 2nd loop over neighbor list, I = sphere, J = global tri + // store atom IDs and shear history for owned spheres + // re-zero npartner to use as counter + + for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + allshear = firstshear[i]; + jnum = numneigh[i]; + touch = firsttouch[i]; + + for (jj = 0; jj < jnum; jj++) { + if (touch[jj]) { + shear = &allshear[3*jj]; + j = jlist[jj]; + m = npartner[i]; + partner[i][m] = j; + shearpartner[i][m][0] = shear[0]; + shearpartner[i][m][1] = shear[1]; + shearpartner[i][m][2] = shear[2]; + npartner[i]++; + } + } + } + + // set maxtouch = max # of partners of any owned atom + // bump up comm->maxexchange_fix if necessary + + maxtouch = 0; + for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]); + comm->maxexchange_fix = MAX(comm->maxexchange_fix,4*maxtouch+1); + + // zero npartner values from previous nlocal_neigh to current nlocal + + int nlocal = atom->nlocal; + for (i = nlocal_neigh; i < nlocal; i++) npartner[i] = 0; +} + /* ---------------------------------------------------------------------- newton on version, for sphere/sphere contacts performs reverse comm to acquire shear partner info from ghost atoms diff --git a/src/fix_shear_history.h b/src/fix_shear_history.h index 0ef3154963..ddfc72a805 100644 --- a/src/fix_shear_history.h +++ b/src/fix_shear_history.h @@ -30,6 +30,7 @@ class FixShearHistory : public Fix { friend class PairGranHookeHistory; friend class PairLineGranHookeHistory; friend class PairTriGranHookeHistory; + friend class FixSurfaceGlobal; public: FixShearHistory(class LAMMPS *, int, char **); @@ -58,6 +59,7 @@ class FixShearHistory : public Fix { protected: int newton_pair; // same as force setting int onesided; // 1 for line/tri history, else 0 + int surfglobal; // 1 for surf/global history, else 0 int nlocal_neigh; // nlocal at last time neigh list was built int nall_neigh; // ditto for nlocal+nghost @@ -68,14 +70,17 @@ class FixShearHistory : public Fix { int commflag; // mode of reverse comm to get ghost info - class Pair *pair; - int *computeflag; // computeflag in PairGranHookeHistory + class Pair *pair; // pointers to callers and their data + class Fix *fsg; + class NeighList *fsg_list; + class NeighList *fsg_listgranhistory; int pgsize,oneatom; // copy of settings in Neighbor MyPage *ipage; // pages of partner atom IDs MyPage *dpage; // pages of shear history with partners void pre_exchange_onesided(); + void pre_exchange_surf(); void pre_exchange_newton(); void pre_exchange_no_newton(); void allocate_pages(); diff --git a/src/molecule.cpp b/src/molecule.cpp index 5e9aa292d0..a792c3674d 100644 --- a/src/molecule.cpp +++ b/src/molecule.cpp @@ -103,6 +103,9 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : sizescale = force->numeric(FLERR,arg[iarg+1]); if (sizescale <= 0.0) error->all(FLERR,"Illegal molecule command"); iarg += 2; + + // NOTE: add other geometric ops for points/lines/tris + } else break; } @@ -117,11 +120,14 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : initialize(); - // scan file for sizes of all fields and allocate them + // scan file for sizes of all fields if (me == 0) open(arg[ifile]); read(0); if (me == 0) fclose(fp); + + // allocate memory for all fields defined in file + allocate(); // read file again to populate all fields @@ -133,22 +139,34 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : // stats if (me == 0) { - if (screen) - fprintf(screen,"Read molecule %s:\n" - " %d atoms with %d types\n %d bonds with %d types\n" - " %d angles with %d types\n %d dihedrals with %d types\n" - " %d impropers with %d types\n", - id,natoms,ntypes, - nbonds,nbondtypes,nangles,nangletypes, - ndihedrals,ndihedraltypes,nimpropers,nimpropertypes); - if (logfile) - fprintf(logfile,"Read molecule %s:\n" - " %d atoms with %d types\n %d bonds with %d types\n" - " %d angles with %d types\n %d dihedrals with %d types\n" - " %d impropers with %d types\n", - id,natoms,ntypes, - nbonds,nbondtypes,nangles,nangletypes, - ndihedrals,ndihedraltypes,nimpropers,nimpropertypes); + if (screen) { + fprintf(screen,"Read molecule %s:\n",id); + if (natoms) + fprintf(screen," %d atoms with %d types\n %d bonds with %d types\n" + " %d angles with %d types\n %d dihedrals with %d types\n" + " %d impropers with %d types\n", + natoms,ntypes, + nbonds,nbondtypes,nangles,nangletypes, + ndihedrals,ndihedraltypes,nimpropers,nimpropertypes); + if (npoints && nlines) + fprintf(screen," %d lines with %d points\n",nlines,npoints); + if (npoints && ntris) + fprintf(screen," %d triangles with %d points\n",ntris,npoints); + } + if (logfile) { + fprintf(logfile,"Read molecule %s:\n",id); + if (natoms) + fprintf(logfile," %d atoms with %d types\n %d bonds with %d types\n" + " %d angles with %d types\n %d dihedrals with %d types\n" + " %d impropers with %d types\n", + natoms,ntypes, + nbonds,nbondtypes,nangles,nangletypes, + ndihedrals,ndihedraltypes,nimpropers,nimpropertypes); + if (npoints && nlines) + fprintf(logfile," %d lines with %d points\n",nlines,npoints); + if (npoints && ntris) + fprintf(logfile," %d triangles with %d points\n",ntris,npoints); + } } } @@ -417,14 +435,20 @@ void Molecule::read(int flag) if (strspn(line," \t\n\r") == strlen(line)) continue; // search line for header keywords and set corresponding variable + // check for triangles before angles so "triangles" not matched as "angles" if (strstr(line,"atoms")) sscanf(line,"%d",&natoms); + + else if (strstr(line,"points")) sscanf(line,"%d",&npoints); + else if (strstr(line,"lines")) sscanf(line,"%d",&nlines); + else if (strstr(line,"triangles")) sscanf(line,"%d",&ntris); + else if (strstr(line,"bonds")) sscanf(line,"%d",&nbonds); else if (strstr(line,"angles")) sscanf(line,"%d",&nangles); else if (strstr(line,"dihedrals")) sscanf(line,"%d",&ndihedrals); else if (strstr(line,"impropers")) sscanf(line,"%d",&nimpropers); - else if (strstr(line,"mass")) { + else if (strstr(line," mass")) { massflag = 1; sscanf(line,"%lg",&masstotal); masstotal *= sizescale*sizescale*sizescale; @@ -463,8 +487,11 @@ void Molecule::read(int flag) // error checks - if (natoms < 1) - error->all(FLERR,"No count or invalid atom count in molecule file"); + if (natoms < 0) error->all(FLERR,"Invalid atom count in molecule file"); + if (npoints < 0) error->all(FLERR,"Invalid point count in molecule file"); + if (nlines < 0) error->all(FLERR,"Invalid line count in molecule file"); + if (ntris < 0) error->all(FLERR,"Invalid triangle count in molecule file"); + if (nbonds < 0) error->all(FLERR,"Invalid bond count in molecule file"); if (nangles < 0) error->all(FLERR,"Invalid angle count in molecule file"); if (ndihedrals < 0) @@ -472,6 +499,21 @@ void Molecule::read(int flag) if (nimpropers < 0) error->all(FLERR,"Invalid improper count in molecule file"); + if (natoms == 0 && npoints == 0) + error->all(FLERR,"Molecule file must define either atoms or points"); + if (natoms && npoints) + error->all(FLERR,"Molecule file cannot define both atoms and points"); + + if (npoints && domain->dimension == 2 && nlines == 0) + error->all(FLERR,"Molecule file must define lines with points"); + if (npoints && domain->dimension == 2 && ntris) + error->all(FLERR,"Molecule file cannot define triangles for 2d model"); + + if (npoints && domain->dimension == 3 && ntris == 0) + error->all(FLERR,"Molecule file must define triangles with points"); + if (npoints && domain->dimension == 3 && nlines) + error->all(FLERR,"Molecule file cannot define lines for 3d model"); + // count = vector for tallying bonds,angles,etc per atom if (flag == 0) memory->create(count,natoms,"molecule:count"); @@ -506,6 +548,19 @@ void Molecule::read(int flag) if (flag) masses(line); else skip_lines(natoms,line); + } else if (strcmp(keyword,"Points") == 0) { + pointflag = 1; + if (flag) pts(line); + else skip_lines(npoints,line); + } else if (strcmp(keyword,"Lines") == 0) { + lineflag = 1; + if (flag) line_segments(line); + else skip_lines(nlines,line); + } else if (strcmp(keyword,"Triangles") == 0) { + triflag = 1; + if (flag) triangles(line); + else skip_lines(ntris,line); + } else if (strcmp(keyword,"Bonds") == 0) { if (nbonds == 0) error->all(FLERR,"Molecule file has bonds but no nbonds setting"); @@ -579,6 +634,15 @@ void Molecule::read(int flag) // error check if (flag == 0) { + if (natoms && !xflag) + error->all(FLERR,"Molecule file has no Coords section"); + if (npoints && !pointflag) + error->all(FLERR,"Molecule file has no Points section"); + if (nlines && !lineflag) + error->all(FLERR,"Molecule file has no Lines section"); + if (ntris && !triflag) + error->all(FLERR,"Molecule file has no Triangles section"); + if ((nspecialflag && !specialflag) || (!nspecialflag && specialflag)) error->all(FLERR,"Molecule file needs both Special Bond sections"); if (specialflag && !bondflag) @@ -739,6 +803,100 @@ void Molecule::masses(char *line) if (rmass[i] <= 0.0) error->all(FLERR,"Invalid atom mass in molecule file"); } +/* ---------------------------------------------------------------------- + read line/tri points from file +------------------------------------------------------------------------- */ + +void Molecule::pts(char *line) +{ + int tmp; + for (int i = 0; i < npoints; i++) { + readline(line); + if (i == 0) { + int nwords = atom->count_words(line); + if (nwords != 4) + error->all(FLERR,"Invalid Points section in molecule file"); + } + sscanf(line,"%d %lg %lg %lg",&tmp, + &points[i][0],&points[i][1],&points[i][2]); + + // apply geometric operations to each point + + points[i][0] *= sizescale; + points[i][1] *= sizescale; + points[i][2] *= sizescale; + } + + if (domain->dimension == 2) { + for (int i = 0; i < natoms; i++) + if (points[i][2] != 0.0) + error->all(FLERR,"Molecule file z point must be 0.0 for 2d"); + } +} + +/* ---------------------------------------------------------------------- + read line segments from file + do not subtract one from point indices +------------------------------------------------------------------------- */ + +void Molecule::line_segments(char *line) +{ + int tmp; + for (int i = 0; i < nlines; i++) { + readline(line); + if (i == 0) { + int nwords = atom->count_words(line); + if (nwords != 5) + error->all(FLERR,"Invalid Lines section in molecule file"); + } + sscanf(line,"%d %d %d %d %d",&tmp,&molline[i],&typeline[i], + &lines[i][0],&lines[i][1]); + typeline[i] += toffset; + } + + // check all types and point indices + + for (int i = 0; i < nlines; i++) { + if (typeline[i] <= 0) + error->all(FLERR,"Invalid line type in molecule file"); + if (lines[i][0] <= 0 || lines[i][0] > npoints || + lines[i][1] <= 0 || lines[i][1] > npoints) + error->all(FLERR,"Invalid point index for line in molecule file"); + } +} + +/* ---------------------------------------------------------------------- + read triangles from file + do not subtract one from point indices +------------------------------------------------------------------------- */ + +void Molecule::triangles(char *line) +{ + int tmp; + for (int i = 0; i < ntris; i++) { + readline(line); + if (i == 0) { + int nwords = atom->count_words(line); + if (nwords != 6) + error->all(FLERR,"Invalid Triangles section in molecule file"); + } + sscanf(line,"%d %d %d %d %d %d",&tmp,&moltri[i],&typetri[i], + &tris[i][0],&tris[i][1],&tris[i][2]); + typetri[i] += toffset; + } + + // check all types and point indices + + for (int i = 0; i < ntris; i++) { + if (typetri[i] <= 0) + error->all(FLERR,"Invalid triangle type in molecule file"); + if (tris[i][0] <= 0 || tris[i][0] > npoints || + tris[i][1] <= 0 || tris[i][1] > npoints || + tris[i][2] <= 0 || tris[i][2] > npoints) + error->all(FLERR,"Invalid point index for triangle in molecule file"); + } +} + /* ---------------------------------------------------------------------- read bonds from file set nbondtypes = max type of any bond @@ -1405,10 +1563,15 @@ void Molecule::initialize() nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0; nibody = ndbody = 0; + npoints = 0; + nlines = 0; + ntris = 0; + bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0; maxspecial = 0; xflag = typeflag = qflag = radiusflag = rmassflag = 0; + pointflag = lineflag = triflag = 0; bondflag = angleflag = dihedralflag = improperflag = 0; nspecialflag = specialflag = 0; shakeflag = shakeflagflag = shakeatomflag = shaketypeflag = 0; @@ -1423,6 +1586,14 @@ void Molecule::initialize() radius = NULL; rmass = NULL; + points = NULL; + lines = NULL; + molline = NULL; + typeline = NULL; + tris = NULL; + moltri = NULL; + typetri = NULL; + num_bond = NULL; bond_type = NULL; bond_atom = NULL; @@ -1467,6 +1638,20 @@ void Molecule::allocate() if (radiusflag) memory->create(radius,natoms,"molecule:radius"); if (rmassflag) memory->create(rmass,natoms,"molecule:rmass"); + // lines or triangles with corner points + + if (pointflag) memory->create(points,npoints,3,"molecule:points"); + if (lineflag) { + memory->create(lines,nlines,2,"molecule:lines"); + memory->create(molline,nlines,"molecule:molline"); + memory->create(typeline,nlines,"molecule:typeline"); + } + if (triflag) { + memory->create(tris,ntris,3,"molecule:tris"); + memory->create(moltri,ntris,"molecule:moltri"); + memory->create(typetri,ntris,"molecule:typetri"); + } + // always allocate num_bond,num_angle,etc and special+nspecial // even if not in molecule file, initialize to 0 // this is so methods that use these arrays don't have to check they exist @@ -1554,6 +1739,14 @@ void Molecule::deallocate() memory->destroy(radius); memory->destroy(rmass); + memory->destroy(points); + memory->destroy(lines); + memory->destroy(molline); + memory->destroy(typeline); + memory->destroy(tris); + memory->destroy(moltri); + memory->destroy(typetri); + memory->destroy(num_bond); memory->destroy(bond_type); memory->destroy(bond_atom); diff --git a/src/molecule.h b/src/molecule.h index 0bbe684636..7ec883ca62 100644 --- a/src/molecule.h +++ b/src/molecule.h @@ -34,6 +34,10 @@ class Molecule : protected Pointers { int nbondtypes,nangletypes,ndihedraltypes,nimpropertypes; int nibody,ndbody; + // for surface lines or tris with corner points + + int npoints,nlines,ntris; + // max bond,angle,etc per atom int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom; @@ -46,6 +50,7 @@ class Molecule : protected Pointers { int nspecialflag,specialflag; int shakeflag,shakeflagflag,shakeatomflag,shaketypeflag; int bodyflag,ibodyflag,dbodyflag; + int pointflag,lineflag,triflag; // 1 if attribute defined or computed, 0 if not @@ -63,6 +68,14 @@ class Molecule : protected Pointers { double *radius; // radius of each atom double *rmass; // mass of each atom + double **points; // end/corner pts of lines or tris + int **lines; // list of line indices into points + int **tris; // list of tri indices into points + int *molline; + int *typeline; + int *moltri; + int *typetri; + int *num_bond; // bonds, angles, dihedrals, impropers for each atom int **bond_type; tagint **bond_atom; @@ -134,6 +147,10 @@ class Molecule : protected Pointers { void charges(char *); void diameters(char *); void masses(char *); + void pts(char *); + void line_segments(char *); + void triangles(char *); + void bonds(int, char *); void angles(int, char *); void dihedrals(int, char *);