diff --git a/src/fix_shear_history.h b/src/fix_shear_history.h index d624c9efef..66ff9fc521 100644 --- a/src/fix_shear_history.h +++ b/src/fix_shear_history.h @@ -28,8 +28,8 @@ namespace LAMMPS_NS { class FixShearHistory : public Fix { friend class Neighbor; friend class PairGranHookeHistory; - friend class PairGranLine; - friend class PairGranTri; + friend class PairLineGranHooke; + friend class PairTriGranHooke; public: FixShearHistory(class LAMMPS *, int, char **); diff --git a/src/neighbor.cpp b/src/neighbor.cpp index f82a20acd9..8bd70ee9a8 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -24,6 +24,8 @@ #include "neigh_request.h" #include "atom.h" #include "atom_vec.h" +#include "atom_vec_line.h" +#include "atom_vec_tri.h" #include "comm.h" #include "force.h" #include "pair.h" @@ -35,6 +37,7 @@ #include "update.h" #include "respa.h" #include "output.h" +#include "math_extra.h" #include "citeme.h" #include "memory.h" #include "error.h" @@ -95,6 +98,8 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp) maxhold = 0; xhold = NULL; + line_hold = NULL; + tri_hold = NULL; lastcall = -1; // binning @@ -180,6 +185,8 @@ Neighbor::~Neighbor() delete [] fixchecklist; memory->destroy(xhold); + memory->destroy(line_hold); + memory->destroy(tri_hold); memory->destroy(binhead); memory->destroy(bins); @@ -241,6 +248,15 @@ void Neighbor::init() // ------------------------------------------------------------------ // settings + // linetri_flag = 1/2 if atom style allows for lines/tris + + avec_line = (AtomVecLine *) atom->style_match("line"); + avec_tri = (AtomVecTri *) atom->style_match("tri"); + + linetri_flag = 0; + if (avec_line) linetri_flag = 1; + if (avec_tri) linetri_flag = 2; + // bbox lo/hi = bounding box of entire domain, stored by Domain if (triclinic == 0) { @@ -379,6 +395,14 @@ void Neighbor::init() memory->destroy(xhold); maxhold = 0; xhold = NULL; + + if (linetri_flag == 1) { + memory->destroy(line_hold); + line_hold = NULL; + } else if (linetri_flag == 2) { + memory->destroy(tri_hold); + tri_hold = NULL; + } } if (style == NSQ) { @@ -389,6 +413,7 @@ void Neighbor::init() bins = NULL; // for USER-DPD Shardlow Splitting Algorithm (SSA) + memory->destroy(bins_ssa); memory->destroy(binhead_ssa); memory->destroy(gbinhead_ssa); @@ -399,11 +424,18 @@ void Neighbor::init() } // 1st time allocation of xhold and bins + // also line/tri hold if linetri_flag is set if (dist_check) { if (maxhold == 0) { maxhold = atom->nmax; memory->create(xhold,maxhold,3,"neigh:xhold"); + if (linetri_flag) { + if (linetri_flag == 1) + memory->create(line_hold,maxhold,4,"neigh:line_hold"); + else + memory->create(tri_hold,maxhold,9,"neigh:tri_hold"); + } } } @@ -1505,12 +1537,25 @@ int Neighbor::check_distance() if (includegroup) nlocal = atom->nfirst; int flag = 0; + for (int i = 0; i < nlocal; i++) { delx = x[i][0] - xhold[i][0]; dely = x[i][1] - xhold[i][1]; delz = x[i][2] - xhold[i][2]; rsq = delx*delx + dely*dely + delz*delz; - if (rsq > deltasq) flag = 1; + if (rsq > deltasq) { + flag = 1; + break; + } + } + + // if line or tri particles: + // also check distance moved by corner pts + // since rotation could mean corners move when x coord does not + + if (!flag && linetri_flag) { + if (linetri_flag == 1) flag = check_distance_line(deltasq); + else flag = check_distance_tri(deltasq); } int flagall; @@ -1519,10 +1564,154 @@ int Neighbor::check_distance() return flagall; } +/* ---------------------------------------------------------------------- + if any line end pt moved deltasq, return 1 +------------------------------------------------------------------------- */ + +int Neighbor::check_distance_line(double deltasq) +{ + double length,theta,dx,dy,rsq; + double endpts[4]; + + AtomVecLine::Bonus *bonus = avec_line->bonus; + double **x = atom->x; + int *line = atom->line; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (line[i] < 0) continue; + length = bonus[line[i]].length; + theta = bonus[line[i]].theta; + dx = 0.5*length*cos(theta); + dy = 0.5*length*sin(theta); + endpts[0] = x[i][0] - dx; + endpts[1] = x[i][1] - dy; + endpts[2] = x[i][0] + dx; + endpts[3] = x[i][1] + dy; + + dx = endpts[0] - line_hold[i][0]; + dy = endpts[1] - line_hold[i][1]; + rsq = dx*dx + dy*dy; + if (rsq > deltasq) return 1; + + dx = endpts[2] - line_hold[i][2]; + dy = endpts[3] - line_hold[i][3]; + rsq = dx*dx + dy*dy; + if (rsq > deltasq) return 1; + } + + return 0; +} + +/* ---------------------------------------------------------------------- + compute and store current line end pts in line_hold +------------------------------------------------------------------------- */ + +void Neighbor::calculate_endpts() +{ + double length,theta,dx,dy; + double *endpt; + + AtomVecLine::Bonus *bonus = avec_line->bonus; + double **x = atom->x; + int *line = atom->line; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (line[i] < 0) continue; + endpt = line_hold[i]; + length = bonus[line[i]].length; + theta = bonus[line[i]].theta; + dx = 0.5*length*cos(theta); + dy = 0.5*length*sin(theta); + endpt[0] = x[i][0] - dx; + endpt[1] = x[i][1] - dy; + endpt[2] = x[i][0] + dx; + endpt[3] = x[i][1] + dy; + } +} + +/* ---------------------------------------------------------------------- + if any tri corner pt moved deltasq, return 1 +------------------------------------------------------------------------- */ + +int Neighbor::check_distance_tri(double deltasq) +{ + int ibonus; + double dx,dy,dz,rsq; + double p[3][3],corner[9]; + + AtomVecTri::Bonus *bonus = avec_tri->bonus; + double **x = atom->x; + int *tri = atom->tri; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (tri[i] < 0) continue; + ibonus = tri[i]; + MathExtra::quat_to_mat(bonus[ibonus].quat,p); + MathExtra::matvec(p,bonus[ibonus].c1,&corner[0]); + MathExtra::add3(x[i],&corner[0],&corner[0]); + MathExtra::matvec(p,bonus[ibonus].c2,&corner[3]); + MathExtra::add3(x[i],&corner[3],&corner[3]); + MathExtra::matvec(p,bonus[ibonus].c3,&corner[6]); + MathExtra::add3(x[i],&corner[6],&corner[6]); + + dx = corner[0] - tri_hold[i][0]; + dy = corner[1] - tri_hold[i][1]; + dz = corner[2] - tri_hold[i][2]; + rsq = dx*dx + dy*dy + dz*dz; + if (rsq > deltasq) return 1; + + dx = corner[3] - tri_hold[i][3]; + dy = corner[4] - tri_hold[i][4]; + dz = corner[5] - tri_hold[i][5]; + rsq = dx*dx + dy*dy + dz*dz; + if (rsq > deltasq) return 1; + + dx = corner[6] - tri_hold[i][6]; + dy = corner[7] - tri_hold[i][7]; + dz = corner[8] - tri_hold[i][8]; + rsq = dx*dx + dy*dy + dz*dz; + if (rsq > deltasq) return 1; + } + + return 0; +} + +/* ---------------------------------------------------------------------- + compute and store current tri corner pts in tri_hold +------------------------------------------------------------------------- */ + +void Neighbor::calculate_corners() +{ + int ibonus; + double p[3][3]; + double *corner; + + AtomVecTri::Bonus *bonus = avec_tri->bonus; + double **x = atom->x; + int *tri = atom->tri; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (tri[i] < 0) continue; + ibonus = tri[i]; + corner = tri_hold[i]; + MathExtra::quat_to_mat(bonus[ibonus].quat,p); + MathExtra::matvec(p,bonus[ibonus].c1,&corner[0]); + MathExtra::add3(x[i],&corner[0],&corner[0]); + MathExtra::matvec(p,bonus[ibonus].c2,&corner[3]); + MathExtra::add3(x[i],&corner[3],&corner[3]); + MathExtra::matvec(p,bonus[ibonus].c3,&corner[6]); + MathExtra::add3(x[i],&corner[6],&corner[6]); + } +} + /* ---------------------------------------------------------------------- build perpetual neighbor lists called at setup and every few timesteps during run or minimization - topology lists also built if topoflag = 1, USER-CUDA calls with topoflag = 0 + topology lists also built if topoflag = 1, USER-CUDA called with tflag = 0 ------------------------------------------------------------------------- */ void Neighbor::build(int topoflag) @@ -1543,12 +1732,25 @@ void Neighbor::build(int topoflag) maxhold = atom->nmax; memory->destroy(xhold); memory->create(xhold,maxhold,3,"neigh:xhold"); + if (linetri_flag) { + if (linetri_flag == 1) { + memory->destroy(line_hold); + memory->create(line_hold,maxhold,4,"neigh:line_hold"); + } else { + memory->destroy(tri_hold); + memory->create(tri_hold,maxhold,4,"neigh:tri_hold"); + } + } } for (i = 0; i < nlocal; i++) { xhold[i][0] = x[i][0]; xhold[i][1] = x[i][1]; xhold[i][2] = x[i][2]; } + if (linetri_flag) { + if (linetri_flag == 1) calculate_endpts(); + else calculate_corners(); + } if (boxcheck) { if (triclinic == 0) { boxlo_hold[0] = bboxlo[0]; @@ -2204,6 +2406,8 @@ bigint Neighbor::memory_usage() { bigint bytes = 0; bytes += memory->usage(xhold,maxhold,3); + if (linetri_flag == 1) bytes += memory->usage(line_hold,maxhold,4); + if (linetri_flag == 2) bytes += memory->usage(tri_hold,maxhold,9); if (style != NSQ) { bytes += memory->usage(bins,maxbin); @@ -2232,4 +2436,3 @@ int Neighbor::exclude_setting() { return exclude; } - diff --git a/src/neighbor.h b/src/neighbor.h index b44e0fde00..167af590ea 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -116,6 +116,12 @@ class Neighbor : protected Pointers { double boxlo_hold[3],boxhi_hold[3]; // box size at last neighbor build double corners_hold[8][3]; // box corners at last neighbor build + int linetri_flag; // 1 if lines exist, 2 if tris exist + double **line_hold; // line corner pts at last neighbor build + double **tri_hold; // tri corner pts at last neighbor build + class AtomVecLine *avec_line; // used to extract line info + class AtomVecTri *avec_tri; // used to extract tri info + int binatomflag; // bin atoms or not when build neigh list // turned off by build_one() @@ -184,6 +190,11 @@ class Neighbor : protected Pointers { // methods + int check_distance_line(double); // check line move dist since last neigh + int check_distance_tri(double); // check tri move dist since last neigh + void calculate_endpts(); + void calculate_corners(); + void bin_atoms(); // bin all atoms double bin_distance(int, int, int); // distance between binx int coord2bin(double *); // mapping atom coord to a bin