diff --git a/src/BODY/body_rounded_polygon.cpp b/src/BODY/body_rounded_polygon.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d848a8fa959adcf2fc1e723ba96a2b0b2f0c613c
--- /dev/null
+++ b/src/BODY/body_rounded_polygon.cpp
@@ -0,0 +1,452 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)
+------------------------------------------------------------------------- */
+
+#include <cstdlib>
+#include "body_rounded_polygon.h"
+#include "atom_vec_body.h"
+#include "atom.h"
+#include "force.h"
+#include "domain.h"
+#include "math_extra.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define EPSILON 1.0e-7
+enum{SPHERE,LINE}; // also in DumpImage
+
+/* ---------------------------------------------------------------------- */
+
+BodyRoundedPolygon::BodyRoundedPolygon(LAMMPS *lmp, int narg, char **arg) :
+ Body(lmp, narg, arg)
+{
+ if (narg != 3) error->all(FLERR,"Invalid body rounded/polygon command");
+
+ if (domain->dimension != 2)
+ error->all(FLERR,"Atom_style body rounded/polygon "
+ "can only be used in 2d simulations");
+
+ // nmin and nmax are minimum and maximum number of vertices
+
+ int nmin = force->inumeric(FLERR,arg[1]);
+ int nmax = force->inumeric(FLERR,arg[2]);
+ if (nmin <= 0 || nmin > nmax)
+ error->all(FLERR,"Invalid body rounded/polygon command");
+
+ size_forward = 0;
+
+ // 1 integer for number of vertices,
+ // 3*nmax doubles for vertex coordinates + 2*nmax doubles for edge ends
+ // 1 double for the enclosing radius
+ // 1 double for the rounded radius
+
+ size_border = 1 + 3*nmax + 2*nmax + 1 + 1;
+
+ // NOTE: need to set appropriate nnbin param for dcp
+
+ icp = new MyPoolChunk<int>(1,1);
+ dcp = new MyPoolChunk<double>(3*nmin+2*nmin+1+1,3*nmax+2*nmax+1+1);
+
+ memory->create(imflag,nmax,"body/rounded/polygon:imflag");
+ memory->create(imdata,nmax,7,"body/nparticle:imdata");
+}
+
+/* ---------------------------------------------------------------------- */
+
+BodyRoundedPolygon::~BodyRoundedPolygon()
+{
+ delete icp;
+ delete dcp;
+ memory->destroy(imflag);
+ memory->destroy(imdata);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolygon::nsub(AtomVecBody::Bonus *bonus)
+{
+ return bonus->ivalue[0];
+}
+
+/* ---------------------------------------------------------------------- */
+
+double *BodyRoundedPolygon::coords(AtomVecBody::Bonus *bonus)
+{
+ return bonus->dvalue;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolygon::nedges(AtomVecBody::Bonus *bonus)
+{
+ int nvertices = bonus->ivalue[0];
+ if (nvertices == 1) return 0;
+ else if (nvertices == 2) return 1;
+ return nvertices;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double *BodyRoundedPolygon::edges(AtomVecBody::Bonus *bonus)
+{
+ return bonus->dvalue+3*nsub(bonus);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double BodyRoundedPolygon::enclosing_radius(struct AtomVecBody::Bonus *bonus)
+{
+ int nvertices = bonus->ivalue[0];
+ if (nvertices == 1 || nvertices == 2)
+ return *(bonus->dvalue+3*nsub(bonus)+2);
+ return *(bonus->dvalue+3*nsub(bonus)+2*nsub(bonus));
+}
+
+/* ---------------------------------------------------------------------- */
+
+double BodyRoundedPolygon::rounded_radius(struct AtomVecBody::Bonus *bonus)
+{
+ int nvertices = bonus->ivalue[0];
+ if (nvertices == 1 || nvertices == 2)
+ return *(bonus->dvalue+3*nsub(bonus)+2+1);
+ return *(bonus->dvalue+3*nsub(bonus)+2*nsub(bonus)+1);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolygon::pack_border_body(AtomVecBody::Bonus *bonus, double *buf)
+{
+ int nsub = bonus->ivalue[0];
+ buf[0] = nsub;
+ memcpy(&buf[1],bonus->dvalue,(3*nsub+2*nsub+1+1)*sizeof(double));
+ return 1+(3*nsub+2*nsub+1+1);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolygon::unpack_border_body(AtomVecBody::Bonus *bonus,
+ double *buf)
+{
+ int nsub = static_cast<int> (buf[0]);
+ bonus->ivalue[0] = nsub;
+ memcpy(bonus->dvalue,&buf[1],(3*nsub+2*nsub+1+1)*sizeof(double));
+ return 1+(3*nsub+2*nsub+1+1);
+}
+
+/* ----------------------------------------------------------------------
+ populate bonus data structure with data file values
+------------------------------------------------------------------------- */
+
+void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
+ int *ifile, double *dfile)
+{
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+
+ // set ninteger, ndouble in bonus and allocate 2 vectors of ints, doubles
+
+ if (ninteger != 1)
+ error->one(FLERR,"Incorrect # of integer values in "
+ "Bodies section of data file");
+ int nsub = ifile[0];
+ if (nsub < 1)
+ error->one(FLERR,"Incorrect integer value in "
+ "Bodies section of data file");
+
+ // nentries = number of double entries to be read from Body section:
+ // 6 for inertia + 3*nsub for vertex coords + 1 for rounded radius
+
+ int nentries = 6 + 3*nsub + 1;
+ if (ndouble != nentries)
+ error->one(FLERR,"Incorrect # of floating-point values in "
+ "Bodies section of data file");
+
+ bonus->ninteger = 1;
+ bonus->ivalue = icp->get(bonus->iindex);
+ bonus->ivalue[0] = nsub;
+ bonus->ndouble = 3*nsub + 2*nsub + 1 + 1;
+ bonus->dvalue = dcp->get(bonus->ndouble,bonus->dindex);
+
+ // diagonalize inertia tensor
+
+ double tensor[3][3];
+ tensor[0][0] = dfile[0];
+ tensor[1][1] = dfile[1];
+ tensor[2][2] = dfile[2];
+ tensor[0][1] = tensor[1][0] = dfile[3];
+ tensor[0][2] = tensor[2][0] = dfile[4];
+ tensor[1][2] = tensor[2][1] = dfile[5];
+
+ double *inertia = bonus->inertia;
+ double evectors[3][3];
+ int ierror = MathExtra::jacobi(tensor,inertia,evectors);
+ if (ierror) error->one(FLERR,
+ "Insufficient Jacobi rotations for body nparticle");
+
+ // if any principal moment < scaled EPSILON, set to 0.0
+
+ double max;
+ max = MAX(inertia[0],inertia[1]);
+ max = MAX(max,inertia[2]);
+
+ if (inertia[0] < EPSILON*max) inertia[0] = 0.0;
+ if (inertia[1] < EPSILON*max) inertia[1] = 0.0;
+ if (inertia[2] < EPSILON*max) inertia[2] = 0.0;
+
+ // exyz_space = principal axes in space frame
+
+ double ex_space[3],ey_space[3],ez_space[3];
+
+ ex_space[0] = evectors[0][0];
+ ex_space[1] = evectors[1][0];
+ ex_space[2] = evectors[2][0];
+ ey_space[0] = evectors[0][1];
+ ey_space[1] = evectors[1][1];
+ ey_space[2] = evectors[2][1];
+ ez_space[0] = evectors[0][2];
+ ez_space[1] = evectors[1][2];
+ ez_space[2] = evectors[2][2];
+
+ // enforce 3 evectors as a right-handed coordinate system
+ // flip 3rd vector if needed
+
+ double cross[3];
+ MathExtra::cross3(ex_space,ey_space,cross);
+ if (MathExtra::dot3(cross,ez_space) < 0.0) MathExtra::negate3(ez_space);
+
+ // create initial quaternion
+
+ MathExtra::exyz_to_q(ex_space,ey_space,ez_space,bonus->quat);
+
+ // bonus->dvalue = the first 3*nsub elements are sub-particle displacements
+ // find the enclosing radius of the body from the maximum displacement
+
+ int i,m;
+ double delta[3], rsq, erad, rrad;
+ double erad2 = 0;
+ int j = 6;
+ int k = 0;
+ for (i = 0; i < nsub; i++) {
+ delta[0] = dfile[j];
+ delta[1] = dfile[j+1];
+ delta[2] = dfile[j+2];
+ MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
+ delta,&bonus->dvalue[k]);
+ rsq = delta[0] * delta[0] + delta[1] * delta[1] +
+ delta[2] * delta[2];
+ if (rsq > erad2) erad2 = rsq;
+ j += 3;
+ k += 3;
+ }
+
+ // .. the next 2*nsub elements are edge ends
+
+ int nedges;
+ if (nsub == 1) { // spheres
+ nedges = 0;
+ bonus->dvalue[k] = 0;
+ *(&bonus->dvalue[k]+1) = 0;
+ k += 2;
+
+ // the last element of bonus->dvalue is the rounded & enclosing radius
+
+ rrad = 0.5 * dfile[j];
+ bonus->dvalue[k] = rrad;
+ erad = rrad;
+
+ k++;
+ bonus->dvalue[k] = rrad;
+
+ atom->radius[bonus->ilocal] = erad;
+
+ } else if (nsub == 2) { // rods
+ nedges = 1;
+ for (i = 0; i < nedges; i++) {
+ bonus->dvalue[k] = 0;
+ *(&bonus->dvalue[k]+1) = 1;
+ k += 2;
+ }
+
+ erad = sqrt(erad2);
+ bonus->dvalue[k] = erad;
+
+ // the last element of bonus->dvalue is the rounded radius
+
+ rrad = 0.5 * dfile[j];
+ k++;
+ bonus->dvalue[k] = rrad;
+
+ atom->radius[bonus->ilocal] = erad + rrad;
+
+ } else { // polygons
+ nedges = nsub;
+ for (i = 0; i < nedges; i++) {
+ bonus->dvalue[k] = i;
+ m = i+1;
+ if (m == nedges) m = 0;
+ *(&bonus->dvalue[k]+1) = m;
+ k += 2;
+ }
+
+ // the next to last element is the enclosing radius
+
+ erad = sqrt(erad2);
+ bonus->dvalue[k] = erad;
+
+ // the last element of bonus->dvalue is the rounded radius
+
+ rrad = 0.5 * dfile[j];
+ k++;
+ bonus->dvalue[k] = rrad;
+
+ atom->radius[bonus->ilocal] = erad + rrad;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ return radius of body particle defined by ifile/dfile params
+ params are ordered as in data file
+ called by Molecule class which needs single body size
+------------------------------------------------------------------------- */
+
+double BodyRoundedPolygon::radius_body(int ninteger, int ndouble,
+ int *ifile, double *dfile)
+{
+ int nsub = ifile[0];
+ if (nsub < 1)
+ error->one(FLERR,"Incorrect integer value in "
+ "Bodies section of data file");
+ if (ndouble != 6 + 3*nsub + 1)
+ error->one(FLERR,"Incorrect # of floating-point values in "
+ "Bodies section of data file");
+
+ // sub-particle coords are relative to body center at (0,0,0)
+ // offset = 6 for sub-particle coords
+
+ double onerad;
+ double maxrad = 0.0;
+ double delta[3];
+
+ int offset = 6;
+ for (int i = 0; i < nsub; i++) {
+ delta[0] = dfile[offset];
+ delta[1] = dfile[offset+1];
+ delta[2] = dfile[offset+2];
+ offset += 3;
+ onerad = MathExtra::len3(delta);
+ maxrad = MAX(maxrad,onerad);
+ }
+
+ // add in radius of rounded corners
+
+ return maxrad + 0.5*dfile[offset];
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolygon::noutcol()
+{
+ // the number of columns for the vertex coordinates
+
+ return 3;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolygon::noutrow(int ibonus)
+{
+ // only return the first nsub rows for the vertex coordinates
+
+ return avec->bonus[ibonus].ivalue[0];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BodyRoundedPolygon::output(int ibonus, int m, double *values)
+{
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+
+ double p[3][3];
+ MathExtra::quat_to_mat(bonus->quat,p);
+ MathExtra::matvec(p,&bonus->dvalue[3*m],values);
+
+ double *x = atom->x[bonus->ilocal];
+ values[0] += x[0];
+ values[1] += x[1];
+ values[2] += x[2];
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolygon::image(int ibonus, double flag1, double flag2,
+ int *&ivec, double **&darray)
+{
+ int j;
+ double p[3][3];
+ double *x, rrad;
+
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+ int n = bonus->ivalue[0];
+
+ if (n == 1) {
+ for (int i = 0; i < n; i++) {
+ imflag[i] = SPHERE;
+ MathExtra::quat_to_mat(bonus->quat,p);
+ MathExtra::matvec(p,&bonus->dvalue[3*i],imdata[i]);
+
+ rrad = enclosing_radius(bonus);
+ x = atom->x[bonus->ilocal];
+ imdata[i][0] += x[0];
+ imdata[i][1] += x[1];
+ imdata[i][2] += x[2];
+ if (flag1 <= 0) imdata[i][3] = 2*rrad;
+ else imdata[i][3] = flag1;
+ }
+
+ } else {
+
+ // first end pt of each line
+
+ for (int i = 0; i < n; i++) {
+ imflag[i] = LINE;
+ MathExtra::quat_to_mat(bonus->quat,p);
+ MathExtra::matvec(p,&bonus->dvalue[3*i],imdata[i]);
+
+ rrad = rounded_radius(bonus);
+ x = atom->x[bonus->ilocal];
+ imdata[i][0] += x[0];
+ imdata[i][1] += x[1];
+ imdata[i][2] += x[2];
+ if (flag1 <= 0) imdata[i][6] = 2*rrad;
+ else imdata[i][6] = flag1;
+ }
+
+ // second end pt of each line
+
+ for (int i = 0; i < n; i++) {
+ j = i+1;
+ if (j == n) j = 0;
+ imdata[i][3] = imdata[j][0];
+ imdata[i][4] = imdata[j][1];
+ imdata[i][5] = imdata[j][2];
+ }
+ }
+
+ ivec = imflag;
+ darray = imdata;
+ return n;
+}
diff --git a/src/BODY/body_rounded_polygon.h b/src/BODY/body_rounded_polygon.h
new file mode 100644
index 0000000000000000000000000000000000000000..b6f45c5cf551cc35fc1ff2d817684f5250cf23f4
--- /dev/null
+++ b/src/BODY/body_rounded_polygon.h
@@ -0,0 +1,86 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef BODY_CLASS
+
+BodyStyle(rounded/polygon,BodyRoundedPolygon)
+
+#else
+
+#ifndef LMP_BODY_ROUNDED_POLYGON_H
+#define LMP_BODY_ROUNDED_POLYGON_H
+
+#include "body.h"
+#include "atom_vec_body.h"
+
+namespace LAMMPS_NS {
+
+class BodyRoundedPolygon : public Body {
+ public:
+ BodyRoundedPolygon(class LAMMPS *, int, char **);
+ ~BodyRoundedPolygon();
+ int nsub(struct AtomVecBody::Bonus *);
+ double *coords(struct AtomVecBody::Bonus *);
+ int nedges(struct AtomVecBody::Bonus *);
+ double *edges(struct AtomVecBody::Bonus *);
+ double enclosing_radius(struct AtomVecBody::Bonus *);
+ double rounded_radius(struct AtomVecBody::Bonus *);
+
+ int pack_border_body(struct AtomVecBody::Bonus *, double *);
+ int unpack_border_body(struct AtomVecBody::Bonus *, double *);
+ void data_body(int, int, int, int *, double *);
+ double radius_body(int, int, int *, double *);
+
+ int noutrow(int);
+ int noutcol();
+ void output(int, int, double *);
+ int image(int, double, double, int *&, double **&);
+
+ private:
+ int *imflag;
+ double **imdata;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Invalid body rounded/polygon command
+
+Arguments in atom-style command are not correct.
+
+E: Invalid format in Bodies section of data file
+
+The specified number of integer or floating point values does not
+appear.
+
+E: Incorrect # of integer values in Bodies section of data file
+
+See doc page for body style.
+
+E: Incorrect integer value in Bodies section of data file
+
+See doc page for body style.
+
+E: Incorrect # of floating-point values in Bodies section of data file
+
+See doc page for body style.
+
+E: Insufficient Jacobi rotations for body nparticle
+
+Eigensolve for rigid body was not sufficiently accurate.
+
+*/
diff --git a/src/BODY/body_rounded_polyhedron.cpp b/src/BODY/body_rounded_polyhedron.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a26b6d0cbd4ef0ddc0f1294a8c6d6a57349eafad
--- /dev/null
+++ b/src/BODY/body_rounded_polyhedron.cpp
@@ -0,0 +1,523 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)
+------------------------------------------------------------------------- */
+
+#include <cstdlib>
+#include "body_rounded_polyhedron.h"
+#include "atom_vec_body.h"
+#include "atom.h"
+#include "force.h"
+#include "domain.h"
+#include "math_extra.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define EPSILON 1.0e-7
+#define MAX_FACE_SIZE 4 // maximum number of vertices per face (for now)
+
+enum{SPHERE,LINE}; // also in DumpImage
+
+/* ---------------------------------------------------------------------- */
+
+BodyRoundedPolyhedron::BodyRoundedPolyhedron(LAMMPS *lmp, int narg, char **arg) :
+ Body(lmp, narg, arg)
+{
+ if (narg != 3) error->all(FLERR,"Invalid body rounded/polygon command");
+
+ // nmin and nmax are minimum and maximum number of vertices
+
+ int nmin = force->inumeric(FLERR,arg[1]);
+ int nmax = force->inumeric(FLERR,arg[2]);
+ if (nmin <= 0 || nmin > nmax)
+ error->all(FLERR,"Invalid body rounded/polyhedron command");
+
+ size_forward = 0;
+
+ // 1 integer for number of vertices,
+ // 3*nmax doubles for vertex coordinates + 2*nmax doubles for edge ends +
+ // (MAX_FACE_SIZE+1)*nmax for faces
+ // 1 double for the enclosing radius
+ // 1 double for the rounded radius
+
+ size_border = 1 + 3*nmax + 2*nmax + MAX_FACE_SIZE*nmax + 1 + 1;
+
+ // NOTE: need to set appropriate nnbin param for dcp
+
+ icp = new MyPoolChunk<int>(1,3);
+ dcp = new MyPoolChunk<double>(3*nmin+2+1+1,
+ 3*nmax+2*nmax+MAX_FACE_SIZE*nmax+1+1);
+
+ memory->create(imflag,2*nmax,"body/rounded/polyhedron:imflag");
+ memory->create(imdata,2*nmax,7,"body/polyhedron:imdata");
+}
+
+/* ---------------------------------------------------------------------- */
+
+BodyRoundedPolyhedron::~BodyRoundedPolyhedron()
+{
+ delete icp;
+ delete dcp;
+ memory->destroy(imflag);
+ memory->destroy(imdata);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolyhedron::nsub(AtomVecBody::Bonus *bonus)
+{
+ return bonus->ivalue[0];
+}
+
+/* ---------------------------------------------------------------------- */
+
+double *BodyRoundedPolyhedron::coords(AtomVecBody::Bonus *bonus)
+{
+ return bonus->dvalue;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolyhedron::nedges(AtomVecBody::Bonus *bonus)
+{
+ int nvertices = bonus->ivalue[0];
+ int nedges = bonus->ivalue[1];
+ int nfaces = bonus->ivalue[2];
+ if (nvertices == 1) return 0;
+ else if (nvertices == 2) return 1;
+ return nedges; //(nvertices+nfaces-2); // Euler's polyon formula: V-E+F=2
+}
+
+/* ---------------------------------------------------------------------- */
+
+double *BodyRoundedPolyhedron::edges(AtomVecBody::Bonus *bonus)
+{
+ return bonus->dvalue+3*nsub(bonus);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolyhedron::nfaces(AtomVecBody::Bonus *bonus)
+{
+ return bonus->ivalue[2];
+}
+
+/* ---------------------------------------------------------------------- */
+
+double *BodyRoundedPolyhedron::faces(AtomVecBody::Bonus *bonus)
+{
+ int nvertices = bonus->ivalue[0];
+ if (nvertices == 1 || nvertices == 2) return NULL;
+ return bonus->dvalue+3*nsub(bonus)+2*nedges(bonus);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double BodyRoundedPolyhedron::enclosing_radius(struct AtomVecBody::Bonus *bonus)
+{
+ int nvertices = bonus->ivalue[0];
+ if (nvertices == 1 || nvertices == 2)
+ return *(bonus->dvalue+3*nsub(bonus)+2);
+ return *(bonus->dvalue+3*nsub(bonus)+2*nedges(bonus)+MAX_FACE_SIZE*nfaces(bonus));
+}
+
+/* ---------------------------------------------------------------------- */
+
+double BodyRoundedPolyhedron::rounded_radius(struct AtomVecBody::Bonus *bonus)
+{
+ int nvertices = bonus->ivalue[0];
+ if (nvertices == 1 || nvertices == 2)
+ return *(bonus->dvalue+3*nsub(bonus)+2+1);
+ return *(bonus->dvalue+3*nsub(bonus)+2*nedges(bonus)+MAX_FACE_SIZE*nfaces(bonus)+1);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolyhedron::pack_border_body(AtomVecBody::Bonus *bonus, double *buf)
+{
+ int nsub = bonus->ivalue[0];
+ int ned = bonus->ivalue[1];
+ int nfac = bonus->ivalue[2];
+ buf[0] = nsub;
+ buf[1] = ned;
+ buf[2] = nfac;
+ int ndouble;
+ if (nsub == 1 || nsub == 2) ndouble = 3*nsub+2+MAX_FACE_SIZE*nfac+1+1;
+ else ndouble = 3*nsub+2*nedges(bonus)+MAX_FACE_SIZE*nfac+1+1;
+ memcpy(&buf[3],bonus->dvalue,ndouble*sizeof(double));
+ return 3+ndouble;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolyhedron::unpack_border_body(AtomVecBody::Bonus *bonus,
+ double *buf)
+{
+ int nsub = static_cast<int> (buf[0]);
+ int ned = static_cast<int> (buf[1]);
+ int nfac = static_cast<int> (buf[2]);
+ bonus->ivalue[0] = nsub;
+ bonus->ivalue[1] = ned;
+ bonus->ivalue[2] = nfac;
+ int ndouble;
+ if (nsub == 1 || nsub == 2) ndouble = 3*nsub+2+MAX_FACE_SIZE*nfac+1+1;
+ else ndouble = 3*nsub+2*nedges(bonus)+MAX_FACE_SIZE*nfac+1+1;
+ memcpy(bonus->dvalue,&buf[3],ndouble*sizeof(double));
+ return 3+ndouble;
+}
+
+/* ----------------------------------------------------------------------
+ populate bonus data structure with data file values
+------------------------------------------------------------------------- */
+
+void BodyRoundedPolyhedron::data_body(int ibonus, int ninteger, int ndouble,
+ int *ifile, double *dfile)
+{
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+
+ // set ninteger, ndouble in bonus and allocate 2 vectors of ints, doubles
+
+ if (ninteger != 3)
+ error->one(FLERR,"Incorrect # of integer values in "
+ "Bodies section of data file");
+ int nsub = ifile[0];
+ int ned = ifile[1];
+ int nfac = ifile[2];
+ if (nsub < 1)
+ error->one(FLERR,"Incorrect integer value in "
+ "Bodies section of data file");
+
+ // nentries = number of double entries to be read from Body section:
+ // nsub == 1 || nsub == 2 || nsub == 3:
+ // 6 for inertia + 3*nsub for vertex coords + 1 for rounded radius
+ // nsub > 3:
+ // 6 for inertia + 3*nsub for vertex coords + 2*nsub for edges + 3*nfaces + 1 for rounded radius
+
+ int nedges,nentries;
+ if (nsub == 1 || nsub == 2) {
+ nentries = 6 + 3*nsub + 1;
+ } else {
+ nedges = ned; //nsub + nfac - 2;
+ nentries = 6 + 3*nsub + 2*nedges + MAX_FACE_SIZE*nfac + 1;
+ }
+ if (ndouble != nentries)
+ error->one(FLERR,"Incorrect # of floating-point values in "
+ "Bodies section of data file");
+
+ bonus->ninteger = 3;
+ bonus->ivalue = icp->get(bonus->iindex);
+ bonus->ivalue[0] = nsub;
+ bonus->ivalue[1] = ned;
+ bonus->ivalue[2] = nfac;
+ if (nsub == 1 || nsub == 2) bonus->ndouble = 3*nsub + 2*nsub + 1 + 1;
+ else bonus->ndouble = 3*nsub + 2*nedges + MAX_FACE_SIZE*nfac + 1 + 1;
+ bonus->dvalue = dcp->get(bonus->ndouble,bonus->dindex);
+
+ // diagonalize inertia tensor
+
+ double tensor[3][3];
+ tensor[0][0] = dfile[0];
+ tensor[1][1] = dfile[1];
+ tensor[2][2] = dfile[2];
+ tensor[0][1] = tensor[1][0] = dfile[3];
+ tensor[0][2] = tensor[2][0] = dfile[4];
+ tensor[1][2] = tensor[2][1] = dfile[5];
+
+ double *inertia = bonus->inertia;
+ double evectors[3][3];
+ int ierror = MathExtra::jacobi(tensor,inertia,evectors);
+ if (ierror) error->one(FLERR,
+ "Insufficient Jacobi rotations for body nparticle");
+
+ // if any principal moment < scaled EPSILON, set to 0.0
+
+ double max;
+ max = MAX(inertia[0],inertia[1]);
+ max = MAX(max,inertia[2]);
+
+ if (inertia[0] < EPSILON*max) inertia[0] = 0.0;
+ if (inertia[1] < EPSILON*max) inertia[1] = 0.0;
+ if (inertia[2] < EPSILON*max) inertia[2] = 0.0;
+
+ // exyz_space = principal axes in space frame
+
+ double ex_space[3],ey_space[3],ez_space[3];
+
+ ex_space[0] = evectors[0][0];
+ ex_space[1] = evectors[1][0];
+ ex_space[2] = evectors[2][0];
+ ey_space[0] = evectors[0][1];
+ ey_space[1] = evectors[1][1];
+ ey_space[2] = evectors[2][1];
+ ez_space[0] = evectors[0][2];
+ ez_space[1] = evectors[1][2];
+ ez_space[2] = evectors[2][2];
+
+ // enforce 3 evectors as a right-handed coordinate system
+ // flip 3rd vector if needed
+
+ double cross[3];
+ MathExtra::cross3(ex_space,ey_space,cross);
+ if (MathExtra::dot3(cross,ez_space) < 0.0) MathExtra::negate3(ez_space);
+
+ // create initial quaternion
+
+ MathExtra::exyz_to_q(ex_space,ey_space,ez_space,bonus->quat);
+
+ // bonus->dvalue = the first 3*nsub elements are sub-particle displacements
+ // find the enclosing radius of the body from the maximum displacement
+
+ int i,m;
+ double delta[3], rsq, erad, rrad;
+ double erad2 = 0;
+ int j = 6;
+ int k = 0;
+ for (i = 0; i < nsub; i++) {
+ delta[0] = dfile[j];
+ delta[1] = dfile[j+1];
+ delta[2] = dfile[j+2];
+ MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
+ delta,&bonus->dvalue[k]);
+ rsq = delta[0] * delta[0] + delta[1] * delta[1] +
+ delta[2] * delta[2];
+ if (rsq > erad2) erad2 = rsq;
+ j += 3;
+ k += 3;
+ }
+
+ // .. the next 2*nsub elements are edge ends
+
+ if (nsub == 1) { // spheres
+ nedges = 0;
+ bonus->dvalue[k] = 0;
+ *(&bonus->dvalue[k]+1) = 0;
+ k += 2;
+
+ rrad = 0.5 * dfile[j];
+ bonus->dvalue[k] = rrad;
+ erad = rrad; // enclosing radius = rounded_radius
+
+ // the last element of bonus->dvalue is the rounded radius
+
+ k++;
+ bonus->dvalue[k] = rrad;
+
+ atom->radius[bonus->ilocal] = erad;
+
+ } else if (nsub == 2) { // rods
+ nedges = 1;
+ for (i = 0; i < nedges; i++) {
+ bonus->dvalue[k] = 0;
+ *(&bonus->dvalue[k]+1) = 1;
+ k += 2;
+ }
+
+ erad = sqrt(erad2);
+ bonus->dvalue[k] = erad;
+
+ // the last element of bonus->dvalue is the rounded radius
+
+ rrad = 0.5 * dfile[j];
+ k++;
+ bonus->dvalue[k] = rrad;
+
+ atom->radius[bonus->ilocal] = erad + rrad;
+
+ } else { // polyhedra
+
+ // edges
+
+ for (i = 0; i < nedges; i++) {
+ bonus->dvalue[k] = dfile[j];
+ *(&bonus->dvalue[k]+1) = dfile[j+1];
+ k += 2;
+ j += 2;
+ }
+
+ // faces
+
+ for (i = 0; i < nfac; i++) {
+ for (m = 0; m < MAX_FACE_SIZE; m++)
+ *(&bonus->dvalue[k]+m) = dfile[j+m];
+ k += MAX_FACE_SIZE;
+ j += MAX_FACE_SIZE;
+ }
+
+ // the next to last element is the enclosing radius
+
+ erad = sqrt(erad2);
+ bonus->dvalue[k] = erad;
+
+ // the last element bonus-> dvalue is the rounded radius
+
+ rrad = 0.5 * dfile[j];
+ k++;
+ bonus->dvalue[k] = rrad;
+
+ atom->radius[bonus->ilocal] = erad + rrad;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ return radius of body particle defined by ifile/dfile params
+ params are ordered as in data file
+ called by Molecule class which needs single body size
+------------------------------------------------------------------------- */
+
+double BodyRoundedPolyhedron::radius_body(int ninteger, int ndouble,
+ int *ifile, double *dfile)
+{
+ int nsub = ifile[0];
+ int ned = ifile[1];
+ int nfac = ifile[2];
+ int nedges = ned; //nsub + nfac - 2;
+
+ int nentries;
+ if (nsub == 1 || nsub == 2) nentries = 6 + 3*nsub + 1;
+ else nentries = 6 + 3*nsub + 2*nedges + MAX_FACE_SIZE*nfac + 1;
+
+ if (nsub < 1)
+ error->one(FLERR,"Incorrect integer value in "
+ "Bodies section of data file");
+ if (ndouble != nentries)
+ error->one(FLERR,"Incorrect # of floating-point values in "
+ "Bodies section of data file");
+
+ // sub-particle coords are relative to body center at (0,0,0)
+ // offset = 6 for sub-particle coords
+
+ double onerad;
+ double maxrad = 0.0;
+ double delta[3];
+
+ int offset = 6;
+ for (int i = 0; i < nsub; i++) {
+ delta[0] = dfile[offset];
+ delta[1] = dfile[offset+1];
+ delta[2] = dfile[offset+2];
+ offset += 3;
+ onerad = MathExtra::len3(delta);
+ maxrad = MAX(maxrad,onerad);
+ }
+
+ if (nsub > 2) offset += (2*nedges+MAX_FACE_SIZE*nfac);
+
+ // add in radius of rounded corners
+
+ return maxrad + 0.5*dfile[offset];
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolyhedron::noutcol()
+{
+ // the number of columns for the vertex coordinates
+
+ return 3;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolyhedron::noutrow(int ibonus)
+{
+ // only return the first nsub rows for the vertex coordinates
+
+ return avec->bonus[ibonus].ivalue[0];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BodyRoundedPolyhedron::output(int ibonus, int m, double *values)
+{
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+
+ double p[3][3];
+ MathExtra::quat_to_mat(bonus->quat,p);
+ MathExtra::matvec(p,&bonus->dvalue[3*m],values);
+
+ double *x = atom->x[bonus->ilocal];
+ values[0] += x[0];
+ values[1] += x[1];
+ values[2] += x[2];
+}
+
+/* ---------------------------------------------------------------------- */
+
+int BodyRoundedPolyhedron::image(int ibonus, double flag1, double flag2,
+ int *&ivec, double **&darray)
+{
+ int j, nelements;
+ double p[3][3];
+ double *x, rrad;
+
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+ int nvertices = bonus->ivalue[0];
+
+ if (nvertices == 1) { // spheres
+
+ for (int i = 0; i < nvertices; i++) {
+ imflag[i] = SPHERE;
+ MathExtra::quat_to_mat(bonus->quat,p);
+ MathExtra::matvec(p,&bonus->dvalue[3*i],imdata[i]);
+
+ rrad = enclosing_radius(bonus);
+ x = atom->x[bonus->ilocal];
+ imdata[i][0] += x[0];
+ imdata[i][1] += x[1];
+ imdata[i][2] += x[2];
+ if (flag1 <= 0) imdata[i][3] = 2*rrad;
+ else imdata[i][3] = flag1;
+ }
+
+ nelements = nvertices;
+ } else {
+ int nfaces = bonus->ivalue[2];
+ int nedges = bonus->ivalue[1]; //nvertices + nfaces - 2;
+ if (nvertices == 2) nedges = 1; // special case: rods
+ double* edge_ends = &bonus->dvalue[3*nvertices];
+ int pt1, pt2;
+
+ for (int i = 0; i < nedges; i++) {
+ imflag[i] = LINE;
+
+ pt1 = static_cast<int>(edge_ends[2*i]);
+ pt2 = static_cast<int>(edge_ends[2*i+1]);
+
+ MathExtra::quat_to_mat(bonus->quat,p);
+ MathExtra::matvec(p,&bonus->dvalue[3*pt1],imdata[i]);
+ MathExtra::matvec(p,&bonus->dvalue[3*pt2],&imdata[i][3]);
+
+ rrad = rounded_radius(bonus);
+ x = atom->x[bonus->ilocal];
+ imdata[i][0] += x[0];
+ imdata[i][1] += x[1];
+ imdata[i][2] += x[2];
+ imdata[i][3] += x[0];
+ imdata[i][4] += x[1];
+ imdata[i][5] += x[2];
+
+ if (flag1 <= 0) imdata[i][6] = 2*rrad;
+ else imdata[i][6] = flag1;
+ }
+
+ nelements = nedges;
+ }
+
+ ivec = imflag;
+ darray = imdata;
+ return nelements;
+}
diff --git a/src/BODY/body_rounded_polyhedron.h b/src/BODY/body_rounded_polyhedron.h
new file mode 100644
index 0000000000000000000000000000000000000000..e5b15fd8f93f6c0a65d896262a52dd85d7569c41
--- /dev/null
+++ b/src/BODY/body_rounded_polyhedron.h
@@ -0,0 +1,88 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef BODY_CLASS
+
+BodyStyle(rounded/polyhedron,BodyRoundedPolyhedron)
+
+#else
+
+#ifndef LMP_BODY_ROUNDED_POLYHEDRON_H
+#define LMP_BODY_ROUNDED_POLYHEDRON_H
+
+#include "body.h"
+#include "atom_vec_body.h"
+
+namespace LAMMPS_NS {
+
+class BodyRoundedPolyhedron : public Body {
+ public:
+ BodyRoundedPolyhedron(class LAMMPS *, int, char **);
+ ~BodyRoundedPolyhedron();
+ int nsub(struct AtomVecBody::Bonus *);
+ double *coords(struct AtomVecBody::Bonus *);
+ int nedges(struct AtomVecBody::Bonus *);
+ double *edges(struct AtomVecBody::Bonus *);
+ int nfaces(struct AtomVecBody::Bonus *);
+ double *faces(struct AtomVecBody::Bonus *);
+ double enclosing_radius(struct AtomVecBody::Bonus *);
+ double rounded_radius(struct AtomVecBody::Bonus *);
+
+ int pack_border_body(struct AtomVecBody::Bonus *, double *);
+ int unpack_border_body(struct AtomVecBody::Bonus *, double *);
+ void data_body(int, int, int, int *, double *);
+ double radius_body(int, int, int *, double *);
+
+ int noutrow(int);
+ int noutcol();
+ void output(int, int, double *);
+ int image(int, double, double, int *&, double **&);
+
+ private:
+ int *imflag;
+ double **imdata;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Invalid body rounded/polyhedron command
+
+Arguments in atom-style command are not correct.
+
+E: Invalid format in Bodies section of data file
+
+The specified number of integer or floating point values does not
+appear.
+
+E: Incorrect # of integer values in Bodies section of data file
+
+See doc page for body style.
+
+E: Incorrect integer value in Bodies section of data file
+
+See doc page for body style.
+
+E: Incorrect # of floating-point values in Bodies section of data file
+
+See doc page for body style.
+
+E: Insufficient Jacobi rotations for body rounded/polyhedron
+
+Eigensolve for rigid body was not sufficiently accurate.
+
+*/
diff --git a/src/BODY/fix_wall_body_polygon.cpp b/src/BODY/fix_wall_body_polygon.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ea81ae26dfce881840769636e7a28c2d218ab43e
--- /dev/null
+++ b/src/BODY/fix_wall_body_polygon.cpp
@@ -0,0 +1,831 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing authors: Trung Dac Nguyen (ndactrung@gmail.com)
+------------------------------------------------------------------------- */
+
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+#include "fix_wall_body_polygon.h"
+#include "atom.h"
+#include "atom_vec_body.h"
+#include "body_rounded_polygon.h"
+#include "domain.h"
+#include "update.h"
+#include "force.h"
+#include "pair.h"
+#include "modify.h"
+#include "respa.h"
+#include "math_const.h"
+#include "math_extra.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+using namespace MathConst;
+
+enum{XPLANE=0,YPLANE=1,ZCYLINDER}; // XYZ PLANE need to be 0,1,2
+enum{HOOKE,HOOKE_HISTORY};
+
+enum {INVALID=0,NONE=1,VERTEX=2};
+enum {FAR=0,XLO,XHI,YLO,YHI};
+
+//#define _POLYGON_DEBUG
+#define DELTA 10000
+#define EPSILON 1e-2
+#define BIG 1.0e20
+#define MAX_CONTACTS 4 // maximum number of contacts for 2D models
+#define EFF_CONTACTS 2 // effective contacts for 2D models
+
+/* ---------------------------------------------------------------------- */
+
+FixWallBodyPolygon::FixWallBodyPolygon(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg)
+{
+ if (narg < 7) error->all(FLERR,"Illegal fix wall/body/polygon command");
+
+ if (!atom->body_flag)
+ error->all(FLERR,"Fix wall/body/polygon requires atom style body/rounded/polygon");
+
+ restart_peratom = 1;
+ create_attribute = 1;
+
+ // wall/particle coefficients
+
+ kn = force->numeric(FLERR,arg[3]);
+
+ c_n = force->numeric(FLERR,arg[4]);
+ if (strcmp(arg[5],"NULL") == 0) c_t = 0.5 * c_n;
+ else c_t = force->numeric(FLERR,arg[5]);
+
+ if (kn < 0.0 || c_n < 0.0 || c_t < 0.0)
+ error->all(FLERR,"Illegal fix wall/body/polygon command");
+
+ // wallstyle args
+
+ int iarg = 6;
+ if (strcmp(arg[iarg],"xplane") == 0) {
+ if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/body/polygon command");
+ wallstyle = XPLANE;
+ if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
+ else lo = force->numeric(FLERR,arg[iarg+1]);
+ if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
+ else hi = force->numeric(FLERR,arg[iarg+2]);
+ iarg += 3;
+ } else if (strcmp(arg[iarg],"yplane") == 0) {
+ if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/body/polygon command");
+ wallstyle = YPLANE;
+ if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
+ else lo = force->numeric(FLERR,arg[iarg+1]);
+ if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
+ else hi = force->numeric(FLERR,arg[iarg+2]);
+ iarg += 3;
+ } else if (strcmp(arg[iarg],"zcylinder") == 0) {
+ if (narg < iarg+2) error->all(FLERR,"Illegal fix wall/body/polygon command");
+ wallstyle = ZCYLINDER;
+ lo = hi = 0.0;
+ cylradius = force->numeric(FLERR,arg[iarg+1]);
+ iarg += 2;
+ }
+
+ // check for trailing keyword/values
+
+ wiggle = 0;
+
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"wiggle") == 0) {
+ if (iarg+4 > narg) error->all(FLERR,"Illegal fix wall/body/polygon command");
+ if (strcmp(arg[iarg+1],"x") == 0) axis = 0;
+ else if (strcmp(arg[iarg+1],"y") == 0) axis = 1;
+ else if (strcmp(arg[iarg+1],"z") == 0) axis = 2;
+ else error->all(FLERR,"Illegal fix wall/body/polygon command");
+ amplitude = force->numeric(FLERR,arg[iarg+2]);
+ period = force->numeric(FLERR,arg[iarg+3]);
+ wiggle = 1;
+ iarg += 4;
+ } else error->all(FLERR,"Illegal fix wall/body/polygon command");
+ }
+
+ if (wallstyle == XPLANE && domain->xperiodic)
+ error->all(FLERR,"Cannot use wall in periodic dimension");
+ if (wallstyle == YPLANE && domain->yperiodic)
+ error->all(FLERR,"Cannot use wall in periodic dimension");
+ if (wallstyle == ZCYLINDER && (domain->xperiodic || domain->yperiodic))
+ error->all(FLERR,"Cannot use wall in periodic dimension");
+
+ if (wiggle && wallstyle == ZCYLINDER && axis != 2)
+ error->all(FLERR,"Invalid wiggle direction for fix wall/body/polygon");
+
+ // setup oscillations
+
+ if (wiggle) omega = 2.0*MY_PI / period;
+
+ time_origin = update->ntimestep;
+
+ dmax = nmax = 0;
+ discrete = NULL;
+ dnum = dfirst = NULL;
+
+ edmax = ednummax = 0;
+ edge = NULL;
+ ednum = edfirst = NULL;
+
+ enclosing_radius = NULL;
+ rounded_radius = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixWallBodyPolygon::~FixWallBodyPolygon()
+{
+ memory->destroy(discrete);
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+
+ memory->destroy(edge);
+ memory->destroy(ednum);
+ memory->destroy(edfirst);
+
+ memory->destroy(enclosing_radius);
+ memory->destroy(rounded_radius);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixWallBodyPolygon::setmask()
+{
+ int mask = 0;
+ mask |= POST_FORCE;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallBodyPolygon::init()
+{
+ dt = update->dt;
+
+ avec = (AtomVecBody *) atom->style_match("body");
+ if (!avec)
+ error->all(FLERR,"Pair body/rounded/polygon requires atom style body");
+ if (strcmp(avec->bptr->style,"rounded/polygon") != 0)
+ error->all(FLERR,"Pair body/rounded/polygon requires body style rounded/polygon");
+ bptr = (BodyRoundedPolygon *) avec->bptr;
+
+ // set pairstyle from body/polygonular pair style
+
+ if (force->pair_match("body/rounded/polygon",1))
+ pairstyle = HOOKE;
+ else error->all(FLERR,"Fix wall/body/polygon is incompatible with Pair style");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallBodyPolygon::setup(int vflag)
+{
+ if (strstr(update->integrate_style,"verlet"))
+ post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallBodyPolygon::post_force(int vflag)
+{
+ double vwall[3],dx,dy,dz,del1,del2,delxy,delr,rsq,eradi,rradi,wall_pos;
+ int i,ni,npi,ifirst,nei,iefirst,side;
+ double facc[3];
+
+ // set position of wall to initial settings and velocity to 0.0
+ // if wiggle, set wall position and velocity accordingly
+
+ double wlo = lo;
+ double whi = hi;
+ vwall[0] = vwall[1] = vwall[2] = 0.0;
+ if (wiggle) {
+ double arg = omega * (update->ntimestep - time_origin) * dt;
+ if (wallstyle == axis) {
+ wlo = lo + amplitude - amplitude*cos(arg);
+ whi = hi + amplitude - amplitude*cos(arg);
+ }
+ vwall[axis] = amplitude*omega*sin(arg);
+ }
+
+ // loop over all my atoms
+ // rsq = distance from wall
+ // dx,dy,dz = signed distance from wall
+ // for rotating cylinder, reset vwall based on particle position
+ // skip atom if not close enough to wall
+ // if wall was set to NULL, it's skipped since lo/hi are infinity
+ // compute force and torque on atom if close enough to wall
+ // via wall potential matched to pair potential
+
+ double **x = atom->x;
+ double **v = atom->v;
+ double **f = atom->f;
+ int *body = atom->body;
+ double *radius = atom->radius;
+ double **torque = atom->torque;
+ double **angmom = atom->angmom;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ // grow the per-atom lists if necessary and initialize
+
+ if (atom->nmax > nmax) {
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+ memory->destroy(ednum);
+ memory->destroy(edfirst);
+ memory->destroy(enclosing_radius);
+ memory->destroy(rounded_radius);
+ nmax = atom->nmax;
+ memory->create(dnum,nmax,"fix:dnum");
+ memory->create(dfirst,nmax,"fix:dfirst");
+ memory->create(ednum,nmax,"fix:ednum");
+ memory->create(edfirst,nmax,"fix:edfirst");
+ memory->create(enclosing_radius,nmax,"fix:enclosing_radius");
+ memory->create(rounded_radius,nmax,"fix:rounded_radius");
+ }
+
+ ndiscrete = nedge = 0;
+ for (i = 0; i < nlocal; i++)
+ dnum[i] = ednum[i] = 0;
+
+ for (i = 0; i < nlocal; i++) {
+ if (mask[i] & groupbit) {
+
+ if (body[i] < 0) continue;
+
+ dx = dy = dz = 0.0;
+ side = FAR;
+ if (wallstyle == XPLANE) {
+ del1 = x[i][0] - wlo;
+ del2 = whi - x[i][0];
+ if (del1 < del2) {
+ dx = del1;
+ wall_pos = wlo;
+ side = XLO;
+ } else {
+ dx = -del2;
+ wall_pos = whi;
+ side = XHI;
+ }
+ } else if (wallstyle == YPLANE) {
+ del1 = x[i][1] - wlo;
+ del2 = whi - x[i][1];
+ if (del1 < del2) {
+ dy = del1;
+ wall_pos = wlo;
+ side = YLO;
+ } else {
+ dy = -del2;
+ wall_pos = whi;
+ side = YHI;
+ }
+ } else if (wallstyle == ZCYLINDER) {
+ delxy = sqrt(x[i][0]*x[i][0] + x[i][1]*x[i][1]);
+ delr = cylradius - delxy;
+ if (delr > eradi) dz = cylradius;
+ else {
+ dx = -delr/delxy * x[i][0];
+ dy = -delr/delxy * x[i][1];
+ }
+ }
+
+ rsq = dx*dx + dy*dy + dz*dz;
+ if (rsq > radius[i]*radius[i]) continue;
+
+ double r = sqrt(rsq);
+ double rsqinv = 1.0 / rsq;
+
+ if (dnum[i] == 0) body2space(i);
+ npi = dnum[i];
+ ifirst = dfirst[i];
+ nei = ednum[i];
+ iefirst = edfirst[i];
+ eradi = enclosing_radius[i];
+ rradi = rounded_radius[i];
+
+ // reset vertex and edge forces
+
+ for (ni = 0; ni < npi; ni++) {
+ discrete[ifirst+ni][3] = 0;
+ discrete[ifirst+ni][4] = 0;
+ discrete[ifirst+ni][5] = 0;
+ }
+
+ for (ni = 0; ni < nei; ni++) {
+ edge[iefirst+ni][2] = 0;
+ edge[iefirst+ni][3] = 0;
+ edge[iefirst+ni][4] = 0;
+ }
+
+ int interact, num_contacts, done;
+ double delta_a, delta_ua, j_a;
+ Contact contact_list[MAX_CONTACTS];
+
+ num_contacts = 0;
+ facc[0] = facc[1] = facc[2] = 0;
+ interact = vertex_against_wall(i, wall_pos, x, f, torque, side,
+ contact_list, num_contacts, facc);
+
+ if (num_contacts >= 2) {
+
+ // find the first two distinct contacts
+
+ done = 0;
+ for (int m = 0; m < num_contacts-1; m++) {
+ for (int n = m+1; n < num_contacts; n++) {
+ delta_a = contact_separation(contact_list[m], contact_list[n]);
+ if (delta_a > 0) {
+ delta_ua = 1.0;
+ j_a = delta_a / (EFF_CONTACTS * delta_ua);
+ if (j_a < 1.0) j_a = 1.0;
+
+ // scale the force at both contacts
+
+ contact_forces(contact_list[m], j_a, x, v, angmom, f, torque,
+ vwall, facc);
+ contact_forces(contact_list[n], j_a, x, v, angmom, f, torque,
+ vwall, facc);
+ done = 1;
+ break;
+ }
+ }
+ if (done == 1) break;
+ }
+
+ } else if (num_contacts == 1) {
+
+ // if there's only one contact, it should be handled here
+ // since forces/torques have not been accumulated from vertex2wall()
+
+ contact_forces(contact_list[0], 1.0, x, v, angmom, f, torque,
+ vwall, facc);
+ }
+ } // group bit
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallBodyPolygon::reset_dt()
+{
+ dt = update->dt;
+}
+
+/* ----------------------------------------------------------------------
+ convert N sub-particles in body I to space frame using current quaternion
+ store sub-particle space-frame displacements from COM in discrete list
+------------------------------------------------------------------------- */
+
+void FixWallBodyPolygon::body2space(int i)
+{
+ int ibonus = atom->body[i];
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+ int nsub = bptr->nsub(bonus);
+ double *coords = bptr->coords(bonus);
+ int body_num_edges = bptr->nedges(bonus);
+ double* vertices = bptr->edges(bonus);
+ double eradius = bptr->enclosing_radius(bonus);
+ double rradius = bptr->rounded_radius(bonus);
+
+ // get the number of sub-particles (vertices)
+ // and the index of the first vertex of my body in the list
+
+ dnum[i] = nsub;
+ dfirst[i] = ndiscrete;
+
+ // grow the vertex list if necessary
+ // the first 3 columns are for coords, the last 3 for forces
+
+ if (ndiscrete + nsub > dmax) {
+ dmax += DELTA;
+ memory->grow(discrete,dmax,6,"fix:discrete");
+ }
+
+ double p[3][3];
+ MathExtra::quat_to_mat(bonus->quat,p);
+
+ for (int m = 0; m < nsub; m++) {
+ MathExtra::matvec(p,&coords[3*m],discrete[ndiscrete]);
+ discrete[ndiscrete][3] = 0;
+ discrete[ndiscrete][4] = 0;
+ discrete[ndiscrete][5] = 0;
+ ndiscrete++;
+ }
+
+ // get the number of edges (vertices)
+ // and the index of the first edge of my body in the list
+
+ ednum[i] = body_num_edges;
+ edfirst[i] = nedge;
+
+ // grow the edge list if necessary
+ // the first 2 columns are for vertex indices within body,
+ // the last 3 for forces
+
+ if (nedge + body_num_edges > edmax) {
+ edmax += DELTA;
+ memory->grow(edge,edmax,5,"fix:edge");
+ }
+
+ for (int m = 0; m < body_num_edges; m++) {
+ edge[nedge][0] = static_cast<int>(vertices[2*m+0]);
+ edge[nedge][1] = static_cast<int>(vertices[2*m+1]);
+ edge[nedge][2] = 0;
+ edge[nedge][3] = 0;
+ edge[nedge][4] = 0;
+ nedge++;
+ }
+
+ enclosing_radius[i] = eradius;
+ rounded_radius[i] = rradius;
+}
+
+/* ----------------------------------------------------------------------
+ Determine the interaction mode between i's vertices against the wall
+
+ i = atom i (body i)
+ x = atoms' coordinates
+ f = atoms' forces
+ torque = atoms' torques
+ Return:
+ contact_list = list of contacts between i and the wall
+ num_contacts = number of contacts between i's vertices and the wall
+ interact = 0 no interaction with the wall
+ 1 there's at least one vertex of i interacts
+ with the wall
+---------------------------------------------------------------------- */
+
+int FixWallBodyPolygon::vertex_against_wall(int i, double wall_pos,
+ double** x, double** f, double** torque, int side,
+ Contact* contact_list, int &num_contacts, double* facc)
+{
+ int ni, npi, ifirst, interact;
+ double xpi[3], xpj[3], dist, eradi, rradi;
+ double fx, fy, fz, rx, ry, rz;
+ int nlocal = atom->nlocal;
+
+ npi = dnum[i];
+ ifirst = dfirst[i];
+ eradi = enclosing_radius[i];
+ rradi = rounded_radius[i];
+
+ interact = 0;
+
+ // loop through body i's vertices
+
+ for (ni = 0; ni < npi; ni++) {
+
+ // convert body-fixed coordinates to space-fixed, xi
+
+ xpi[0] = x[i][0] + discrete[ifirst+ni][0];
+ xpi[1] = x[i][1] + discrete[ifirst+ni][1];
+ xpi[2] = x[i][2] + discrete[ifirst+ni][2];
+
+ int mode, contact, p2vertex;
+ double d, R, hi[3], t, delx, dely, delz, fpair, shift;
+ double xj[3], rij;
+
+ // compute the distance from the vertex xpi to the wall
+
+ mode = compute_distance_to_wall(xpi, rradi, wall_pos, side,
+ d, hi, contact);
+
+ if (mode == INVALID || mode == NONE) continue;
+
+ if (mode == VERTEX) {
+
+ interact = 1;
+
+ // vertex i interacts with the wall
+
+ delx = xpi[0] - hi[0];
+ dely = xpi[1] - hi[1];
+ delz = xpi[2] - hi[2];
+
+ // R = surface separation = d shifted by the rounded radius
+ // R = d - p1.rounded_radius;
+ // note: the force is defined for R, not for d
+ // R > 0: no interaction
+ // R <= 0: deformation between vertex i and the wall
+
+ rij = sqrt(delx*delx + dely*dely + delz*delz);
+ R = rij - rradi;
+
+ // the normal frictional term -c_n * vn will be added later
+
+ if (R <= 0) { // deformation occurs
+ fpair = -kn * R;
+ } else fpair = 0.0;
+
+ fx = delx*fpair/rij;
+ fy = dely*fpair/rij;
+ fz = delz*fpair/rij;
+
+ #ifdef _POLYGON_DEBUG
+ printf(" Interaction between vertex %d of %d and wall:", ni);
+ printf(" mode = %d; contact = %d; d = %f; rij = %f\n",
+ mode, contact, d, rij);
+ printf(" R = %f\n", R);
+ printf(" fpair = %f\n", fpair);
+ #endif
+
+ if (contact == 1) {
+
+ // vertex ni of body i contacts with edge nj of body j
+
+ contact_list[num_contacts].ibody = i;
+ contact_list[num_contacts].jbody = -1;
+ contact_list[num_contacts].vertex = ni;
+ contact_list[num_contacts].edge = -1;
+ contact_list[num_contacts].xv[0] = xpi[0];
+ contact_list[num_contacts].xv[1] = xpi[1];
+ contact_list[num_contacts].xv[2] = xpi[2];
+ contact_list[num_contacts].xe[0] = hi[0];
+ contact_list[num_contacts].xe[1] = hi[1];
+ contact_list[num_contacts].xe[2] = hi[2];
+ contact_list[num_contacts].separation = R;
+ num_contacts++;
+
+ // store forces to vertex ni to be rescaled later,
+ // if there are 2 contacts
+
+ discrete[ifirst+ni][3] = fx;
+ discrete[ifirst+ni][4] = fy;
+ discrete[ifirst+ni][5] = fz;
+
+ #ifdef _POLYGON_DEBUG
+ printf(" Stored forces at vertex and edge for accumulating later.\n");
+ #endif
+
+ } else { // no contact
+
+ // accumulate force and torque to the body directly
+
+ f[i][0] += fx;
+ f[i][1] += fy;
+ f[i][2] += fz;
+ sum_torque(x[i], xpi, fx, fy, fz, torque[i]);
+
+ } // end if contact
+
+ } // end if mode
+
+ } // end for looping through the vertices of body i
+
+ return interact;
+}
+
+/* -------------------------------------------------------------------------
+ Compute the distance between a vertex to the wall
+ another body
+ Input:
+ x0 = coordinate of the tested vertex
+ rradi = rounded radius of the vertex
+ wall_pos = position of the wall
+ Output:
+ d = Distance from a point x0 to an wall
+ hi = coordinates of the projection of x0 on the wall
+ contact = 0 no contact between the queried vertex and the wall
+ 1 contact detected
+ return NONE if there is no interaction
+ EDGE if the tested vertex interacts with the wall
+------------------------------------------------------------------------- */
+
+int FixWallBodyPolygon::compute_distance_to_wall(double* x0, double rradi,
+ double wall_pos, int side, double &d, double hi[3], int &contact)
+{
+ int mode;
+ double delxy;
+
+ // h0 = position of the projection of x0 on the wall
+ if (wallstyle == XPLANE) {
+ hi[0] = wall_pos;
+ hi[1] = x0[1];
+ hi[2] = x0[2];
+ } else if (wallstyle == YPLANE) {
+ hi[0] = x0[0];
+ hi[1] = wall_pos;
+ hi[2] = x0[2];
+ } else if (wallstyle == ZCYLINDER) {
+ delxy = sqrt(x0[0]*x0[0] + x0[1]*x0[1]);
+ hi[0] = x0[0]*cylradius/delxy;
+ hi[1] = x0[1]*cylradius/delxy;
+ hi[2] = x0[2];
+ }
+
+ // distance from x0 to the wall = distance from x0 to hi
+
+ distance(hi, x0, d);
+
+ // determine the interaction mode
+
+ if (d < rradi) {
+ mode = VERTEX;
+ contact = 1;
+ } else {
+ if (side == XLO) {
+ if (x0[0] < wall_pos) mode = VERTEX;
+ else mode = NONE;
+ } else if (side == XHI) {
+ if (x0[0] > wall_pos) mode = VERTEX;
+ else mode = NONE;
+ } else if (side == YLO) {
+ if (x0[1] < wall_pos) mode = VERTEX;
+ else mode = NONE;
+ } else if (side == YHI) {
+ if (x0[1] > wall_pos) mode = VERTEX;
+ else mode = NONE;
+ }
+ }
+
+ if (mode == NONE) contact = 0;
+ else contact = 1;
+
+ return mode;
+}
+
+/* ----------------------------------------------------------------------
+ Compute the contact forces between two bodies
+ modify the force stored at the vertex and edge in contact by j_a
+ sum forces and torque to the corresponding bodies
+ fn = normal friction component
+ ft = tangential friction component (-c_t * vrt)
+------------------------------------------------------------------------- */
+
+void FixWallBodyPolygon::contact_forces(Contact& contact, double j_a,
+ double** x, double** v, double** angmom, double** f,
+ double** torque, double* vwall, double* facc)
+{
+ int ibody,ibonus,ifirst, jefirst, ni;
+ double fx,fy,fz,delx,dely,delz,rsq,rsqinv;
+ double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
+ double fn[3],ft[3],vi[3];
+ double *quat, *inertia;
+ AtomVecBody::Bonus *bonus;
+
+ ibody = contact.ibody;
+
+ // compute the velocity of the vertex in the space-fixed frame
+
+ ibonus = atom->body[ibody];
+ bonus = &avec->bonus[ibonus];
+ quat = bonus->quat;
+ inertia = bonus->inertia;
+ total_velocity(contact.xv, x[ibody], v[ibody], angmom[ibody],
+ inertia, quat, vi);
+
+ // vector pointing from the vertex to the point on the wall
+
+ delx = contact.xv[0] - contact.xe[0];
+ dely = contact.xv[1] - contact.xe[1];
+ delz = contact.xv[2] - contact.xe[2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ rsqinv = 1.0/rsq;
+
+ // relative translational velocity
+
+ vr1 = vi[0] - vwall[0];
+ vr2 = vi[1] - vwall[1];
+ vr3 = vi[2] - vwall[2];
+
+ // normal component
+
+ vnnr = vr1*delx + vr2*dely + vr3*delz;
+ vn1 = delx*vnnr * rsqinv;
+ vn2 = dely*vnnr * rsqinv;
+ vn3 = delz*vnnr * rsqinv;
+
+ // tangential component
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // normal friction term at contact
+
+ fn[0] = -c_n * vn1;
+ fn[1] = -c_n * vn2;
+ fn[2] = -c_n * vn3;
+
+ // tangential friction term at contact
+ // excluding the tangential deformation term for now
+
+ ft[0] = -c_t * vt1;
+ ft[1] = -c_t * vt2;
+ ft[2] = -c_t * vt3;
+
+ // only the cohesive force is scaled by j_a
+
+ ifirst = dfirst[ibody];
+ ni = contact.vertex;
+
+ fx = discrete[ifirst+ni][3] * j_a + fn[0] + ft[0];
+ fy = discrete[ifirst+ni][4] * j_a + fn[1] + ft[1];
+ fz = discrete[ifirst+ni][5] * j_a + fn[2] + ft[2];
+ f[ibody][0] += fx;
+ f[ibody][1] += fy;
+ f[ibody][2] += fz;
+ sum_torque(x[ibody], contact.xv, fx, fy, fz, torque[ibody]);
+
+ // accumulate forces to the vertex only
+
+ facc[0] += fx; facc[1] += fy; facc[2] += fz;
+
+ #ifdef _POLYGON_DEBUG
+ printf("From contact forces: vertex fx %f fy %f fz %f\n"
+ " torque body %d: %f %f %f\n",
+ discrete[ifirst+ni][3], discrete[ifirst+ni][4], discrete[ifirst+ni][5],
+ atom->tag[ibody],torque[ibody][0],torque[ibody][1],torque[ibody][2]);
+ #endif
+}
+
+/* ----------------------------------------------------------------------
+ Determine the length of the contact segment, i.e. the separation between
+ 2 contacts
+------------------------------------------------------------------------- */
+
+double FixWallBodyPolygon::contact_separation(const Contact& c1,
+ const Contact& c2)
+{
+ double x1 = c1.xv[0];
+ double y1 = c1.xv[1];
+ double x2 = c1.xe[0];
+ double y2 = c1.xe[1];
+ double x3 = c2.xv[0];
+ double y3 = c2.xv[1];
+
+ double delta_a = 0.0;
+ if (fabs(x2 - x1) > EPSILON) {
+ double A = (y2 - y1) / (x2 - x1);
+ delta_a = fabs(y1 - A * x1 - y3 + A * x3) / sqrt(1 + A * A);
+ } else {
+ delta_a = fabs(x1 - x3);
+ }
+
+ return delta_a;
+}
+
+/* ----------------------------------------------------------------------
+ Accumulate torque to body from the force f=(fx,fy,fz) acting at point x
+------------------------------------------------------------------------- */
+
+void FixWallBodyPolygon::sum_torque(double* xm, double *x, double fx,
+ double fy, double fz, double* torque)
+{
+ double rx = x[0] - xm[0];
+ double ry = x[1] - xm[1];
+ double rz = x[2] - xm[2];
+ double tx = ry * fz - rz * fy;
+ double ty = rz * fx - rx * fz;
+ double tz = rx * fy - ry * fx;
+ torque[0] += tx;
+ torque[1] += ty;
+ torque[2] += tz;
+}
+
+/* ----------------------------------------------------------------------
+ Calculate the total velocity of a point (vertex, a point on an edge):
+ vi = vcm + omega ^ (p - xcm)
+------------------------------------------------------------------------- */
+
+void FixWallBodyPolygon::total_velocity(double* p, double *xcm,
+ double* vcm, double *angmom, double *inertia,
+ double *quat, double* vi)
+{
+ double r[3],omega[3],ex_space[3],ey_space[3],ez_space[3];
+ r[0] = p[0] - xcm[0];
+ r[1] = p[1] - xcm[1];
+ r[2] = p[2] - xcm[2];
+ MathExtra::q_to_exyz(quat,ex_space,ey_space,ez_space);
+ MathExtra::angmom_to_omega(angmom,ex_space,ey_space,ez_space,
+ inertia,omega);
+ vi[0] = omega[1]*r[2] - omega[2]*r[1] + vcm[0];
+ vi[1] = omega[2]*r[0] - omega[0]*r[2] + vcm[1];
+ vi[2] = omega[0]*r[1] - omega[1]*r[0] + vcm[2];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallBodyPolygon::distance(const double* x2, const double* x1,
+ double& r) {
+ r = sqrt((x2[0] - x1[0]) * (x2[0] - x1[0])
+ + (x2[1] - x1[1]) * (x2[1] - x1[1])
+ + (x2[2] - x1[2]) * (x2[2] - x1[2]));
+}
+
diff --git a/src/BODY/fix_wall_body_polygon.h b/src/BODY/fix_wall_body_polygon.h
new file mode 100644
index 0000000000000000000000000000000000000000..3521c6c797afe4169a0dea4a01cc3839cced02f2
--- /dev/null
+++ b/src/BODY/fix_wall_body_polygon.h
@@ -0,0 +1,129 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(wall/body/polygon,FixWallBodyPolygon)
+
+#else
+
+#ifndef LMP_FIX_WALL_BODY_POLYGON_H
+#define LMP_FIX_WALL_BODY_POLYGON_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixWallBodyPolygon : public Fix {
+ public:
+ FixWallBodyPolygon(class LAMMPS *, int, char **);
+ virtual ~FixWallBodyPolygon();
+ int setmask();
+ void init();
+ void setup(int);
+ virtual void post_force(int);
+ void reset_dt();
+
+ struct Contact {
+ int ibody, jbody; // body (i.e. atom) indices (not tags)
+ int vertex; // vertex of the first polygon
+ int edge; // edge of the second polygon
+ double xv[3]; // coordinates of the vertex
+ double xe[3]; // coordinates of the projection of the vertex on the edge
+ double separation;// separation at contact
+ };
+
+ protected:
+ int wallstyle,pairstyle,wiggle,axis;
+ double kn,c_n,c_t;
+ double lo,hi,cylradius;
+ double amplitude,period,omega;
+ double dt;
+ int time_origin;
+
+ class AtomVecBody *avec;
+ class BodyRoundedPolygon *bptr;
+
+ double **discrete; // list of all sub-particles for all bodies
+ int ndiscrete; // number of discretes in list
+ int dmax; // allocated size of discrete list
+ int *dnum; // number of discretes per line, 0 if uninit
+ int *dfirst; // index of first discrete per each line
+ int nmax; // allocated size of dnum,dfirst vectors
+
+ double **edge; // list of all edge for all bodies
+ int nedge; // number of edge in list
+ int edmax; // allocated size of edge list
+ int *ednum; // number of edges per line, 0 if uninit
+ int *edfirst; // index of first edge per each line
+ int ednummax; // allocated size of ednum,edfirst vectors
+
+ double *enclosing_radius; // enclosing radii for all bodies
+ double *rounded_radius; // rounded radii for all bodies
+
+ void body2space(int);
+
+ int vertex_against_wall(int ibody, double wall_pos, double** x,
+ double** f, double** torque, int side,
+ Contact* contact_list, int &num_contacts,
+ double* facc);
+
+ int compute_distance_to_wall(double* x0, double rradi, double wall_pos,
+ int side, double &d, double hi[3], int &contact);
+ double contact_separation(const Contact& c1, const Contact& c2);
+ void contact_forces(Contact& contact, double j_a, double** x,
+ double** v, double** angmom, double** f, double** torque,
+ double* vwall, double* facc);
+ void sum_torque(double* xm, double *x, double fx,
+ double fy, double fz, double* torque);
+ void total_velocity(double* p, double *xcm, double* vcm, double *angmom,
+ double *inertia, double *quat, double* vi);
+ void distance(const double* x2, const double* x1, double& r);
+
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Fix wall/body/polygon requires atom style body rounded/polygon
+
+Self-explanatory.
+
+E: Cannot use wall in periodic dimension
+
+Self-explanatory.
+
+E: Cannot wiggle and shear fix wall/body/polygon
+
+Cannot specify both options at the same time.
+
+E: Invalid wiggle direction for fix wall/body/polygon
+
+Self-explanatory.
+
+E: Fix wall/body/polygon is incompatible with Pair style
+
+Must use a body pair style to define the parameters needed for
+this fix.
+
+*/
diff --git a/src/BODY/fix_wall_body_polyhedron.cpp b/src/BODY/fix_wall_body_polyhedron.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4806da9673ecbeaddd06a4d80bf564c3ef740858
--- /dev/null
+++ b/src/BODY/fix_wall_body_polyhedron.cpp
@@ -0,0 +1,944 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing authors: Trung Dac Nguyen (ndactrung@gmail.com)
+------------------------------------------------------------------------- */
+
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+#include "fix_wall_body_polyhedron.h"
+#include "atom.h"
+#include "atom_vec_body.h"
+#include "body_rounded_polyhedron.h"
+#include "domain.h"
+#include "update.h"
+#include "force.h"
+#include "pair.h"
+#include "modify.h"
+#include "respa.h"
+#include "math_const.h"
+#include "math_extra.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+using namespace MathConst;
+
+enum{XPLANE=0,YPLANE=1,ZPLANE}; // XYZ PLANE need to be 0,1,2
+enum{HOOKE,HOOKE_HISTORY};
+
+enum {INVALID=0,NONE=1,VERTEX=2};
+enum {FAR=0,XLO,XHI,YLO,YHI,ZLO,ZHI};
+
+//#define _POLYHEDRON_DEBUG
+#define DELTA 10000
+#define EPSILON 1e-2
+#define BIG 1.0e20
+#define MAX_CONTACTS 4 // maximum number of contacts for 2D models
+#define EFF_CONTACTS 2 // effective contacts for 2D models
+
+/* ---------------------------------------------------------------------- */
+
+FixWallBodyPolyhedron::FixWallBodyPolyhedron(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg)
+{
+ if (narg < 7) error->all(FLERR,"Illegal fix wall/body/polyhedron command");
+
+ if (!atom->body_flag)
+ error->all(FLERR,"Fix wall/body/polyhedron requires atom style body/rounded/polyhedron");
+
+ restart_peratom = 1;
+ create_attribute = 1;
+
+ // wall/particle coefficients
+
+ kn = force->numeric(FLERR,arg[3]);
+
+ c_n = force->numeric(FLERR,arg[4]);
+ if (strcmp(arg[5],"NULL") == 0) c_t = 0.5 * c_n;
+ else c_t = force->numeric(FLERR,arg[5]);
+
+ if (kn < 0.0 || c_n < 0.0 || c_t < 0.0)
+ error->all(FLERR,"Illegal fix wall/body/polyhedron command");
+
+ // wallstyle args
+
+ int iarg = 6;
+ if (strcmp(arg[iarg],"xplane") == 0) {
+ if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/body/polyhedron command");
+ wallstyle = XPLANE;
+ if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
+ else lo = force->numeric(FLERR,arg[iarg+1]);
+ if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
+ else hi = force->numeric(FLERR,arg[iarg+2]);
+ iarg += 3;
+ } else if (strcmp(arg[iarg],"yplane") == 0) {
+ if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/body/polyhedron command");
+ wallstyle = YPLANE;
+ if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
+ else lo = force->numeric(FLERR,arg[iarg+1]);
+ if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
+ else hi = force->numeric(FLERR,arg[iarg+2]);
+ iarg += 3;
+ } else if (strcmp(arg[iarg],"zplane") == 0) {
+ if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/body/polyhedron command");
+ wallstyle = ZPLANE;
+ if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
+ else lo = force->numeric(FLERR,arg[iarg+1]);
+ if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
+ else hi = force->numeric(FLERR,arg[iarg+2]);
+ iarg += 3;
+ }
+
+ // check for trailing keyword/values
+
+ wiggle = 0;
+
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"wiggle") == 0) {
+ if (iarg+4 > narg) error->all(FLERR,"Illegal fix wall/body/polyhedron command");
+ if (strcmp(arg[iarg+1],"x") == 0) axis = 0;
+ else if (strcmp(arg[iarg+1],"y") == 0) axis = 1;
+ else if (strcmp(arg[iarg+1],"z") == 0) axis = 2;
+ else error->all(FLERR,"Illegal fix wall/body/polyhedron command");
+ amplitude = force->numeric(FLERR,arg[iarg+2]);
+ period = force->numeric(FLERR,arg[iarg+3]);
+ wiggle = 1;
+ iarg += 4;
+ } else error->all(FLERR,"Illegal fix wall/body/polyhedron command");
+ }
+
+ if (wallstyle == XPLANE && domain->xperiodic)
+ error->all(FLERR,"Cannot use wall in periodic dimension");
+ if (wallstyle == YPLANE && domain->yperiodic)
+ error->all(FLERR,"Cannot use wall in periodic dimension");
+ if (wallstyle == ZPLANE && domain->zperiodic)
+ error->all(FLERR,"Cannot use wall in periodic dimension");
+
+ // setup oscillations
+
+ if (wiggle) omega = 2.0*MY_PI / period;
+
+ time_origin = update->ntimestep;
+
+ dmax = nmax = 0;
+ discrete = NULL;
+ dnum = dfirst = NULL;
+
+ edmax = ednummax = 0;
+ edge = NULL;
+ ednum = edfirst = NULL;
+
+ facmax = facnummax = 0;
+ face = NULL;
+ facnum = facfirst = NULL;
+
+ enclosing_radius = NULL;
+ rounded_radius = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixWallBodyPolyhedron::~FixWallBodyPolyhedron()
+{
+ memory->destroy(discrete);
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+
+ memory->destroy(edge);
+ memory->destroy(ednum);
+ memory->destroy(edfirst);
+
+ memory->destroy(face);
+ memory->destroy(facnum);
+ memory->destroy(facfirst);
+
+ memory->destroy(enclosing_radius);
+ memory->destroy(rounded_radius);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixWallBodyPolyhedron::setmask()
+{
+ int mask = 0;
+ mask |= POST_FORCE;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallBodyPolyhedron::init()
+{
+ dt = update->dt;
+
+ avec = (AtomVecBody *) atom->style_match("body");
+ if (!avec)
+ error->all(FLERR,"Pair body/rounded/polyhedron requires atom style body");
+ if (strcmp(avec->bptr->style,"rounded/polyhedron") != 0)
+ error->all(FLERR,"Pair body/rounded/polyhedron requires body style rounded/polyhedron");
+ bptr = (BodyRoundedPolyhedron *) avec->bptr;
+
+ // set pairstyle from body/polyhedronular pair style
+
+ if (force->pair_match("body/rounded/polyhedron",1))
+ pairstyle = HOOKE;
+ else error->all(FLERR,"Fix wall/body/polyhedron is incompatible with Pair style");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallBodyPolyhedron::setup(int vflag)
+{
+ if (strstr(update->integrate_style,"verlet"))
+ post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallBodyPolyhedron::post_force(int vflag)
+{
+ double vwall[3],dx,dy,dz,del1,del2,delxy,delr,rsq,eradi,rradi,wall_pos;
+ int i,ni,npi,ifirst,nei,iefirst,nfi,iffirst,side;
+ double facc[3];
+
+ // set position of wall to initial settings and velocity to 0.0
+ // if wiggle, set wall position and velocity accordingly
+
+ double wlo = lo;
+ double whi = hi;
+ vwall[0] = vwall[1] = vwall[2] = 0.0;
+ if (wiggle) {
+ double arg = omega * (update->ntimestep - time_origin) * dt;
+ if (wallstyle == axis) {
+ wlo = lo + amplitude - amplitude*cos(arg);
+ whi = hi + amplitude - amplitude*cos(arg);
+ }
+ vwall[axis] = amplitude*omega*sin(arg);
+ }
+
+ // loop over all my atoms
+ // rsq = distance from wall
+ // dx,dy,dz = signed distance from wall
+ // for rotating cylinder, reset vwall based on particle position
+ // skip atom if not close enough to wall
+ // if wall was set to NULL, it's skipped since lo/hi are infinity
+ // compute force and torque on atom if close enough to wall
+ // via wall potential matched to pair potential
+
+ double **x = atom->x;
+ double **v = atom->v;
+ double **f = atom->f;
+ int *body = atom->body;
+ double *radius = atom->radius;
+ double **torque = atom->torque;
+ double **angmom = atom->angmom;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ // grow the per-atom lists if necessary and initialize
+
+ if (atom->nmax > nmax) {
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+ memory->destroy(ednum);
+ memory->destroy(edfirst);
+ memory->destroy(facnum);
+ memory->destroy(facfirst);
+ memory->destroy(enclosing_radius);
+ memory->destroy(rounded_radius);
+ nmax = atom->nmax;
+ memory->create(dnum,nmax,"fix:dnum");
+ memory->create(dfirst,nmax,"fix:dfirst");
+ memory->create(ednum,nmax,"fix:ednum");
+ memory->create(edfirst,nmax,"fix:edfirst");
+ memory->create(facnum,nmax,"fix:facnum");
+ memory->create(facfirst,nmax,"fix:facfirst");
+ memory->create(enclosing_radius,nmax,"fix:enclosing_radius");
+ memory->create(rounded_radius,nmax,"fix:rounded_radius");
+ }
+
+ ndiscrete = nedge = nface = 0;
+ for (i = 0; i < nlocal; i++)
+ dnum[i] = ednum[i] = facnum[i] = 0;
+
+ for (i = 0; i < nlocal; i++) {
+ if (mask[i] & groupbit) {
+
+ if (body[i] < 0) continue;
+
+ dx = dy = dz = 0.0;
+ side = FAR;
+ if (wallstyle == XPLANE) {
+ del1 = x[i][0] - wlo;
+ del2 = whi - x[i][0];
+ if (del1 < del2) {
+ dx = del1;
+ wall_pos = wlo;
+ side = XLO;
+ } else {
+ dx = -del2;
+ wall_pos = whi;
+ side = XHI;
+ }
+ } else if (wallstyle == YPLANE) {
+ del1 = x[i][1] - wlo;
+ del2 = whi - x[i][1];
+ if (del1 < del2) {
+ dy = del1;
+ wall_pos = wlo;
+ side = YLO;
+ } else {
+ dy = -del2;
+ wall_pos = whi;
+ side = YHI;
+ }
+ } else if (wallstyle == ZPLANE) {
+ del1 = x[i][2] - wlo;
+ del2 = whi - x[i][2];
+ if (del1 < del2) {
+ dy = del1;
+ wall_pos = wlo;
+ side = ZLO;
+ } else {
+ dy = -del2;
+ wall_pos = whi;
+ side = ZHI;
+ }
+ }
+
+ rsq = dx*dx + dy*dy + dz*dz;
+ if (rsq > radius[i]*radius[i]) continue;
+
+ double r = sqrt(rsq);
+ double rsqinv = 1.0 / rsq;
+
+ if (dnum[i] == 0) body2space(i);
+ npi = dnum[i];
+ ifirst = dfirst[i];
+ nei = ednum[i];
+ iefirst = edfirst[i];
+ nfi = facnum[i];
+ iffirst = facfirst[i];
+ eradi = enclosing_radius[i];
+ rradi = rounded_radius[i];
+
+ if (npi == 1) {
+ sphere_against_wall(i, wall_pos, side, vwall, x, v, f, angmom, torque);
+ continue;
+ }
+
+ // reset vertex and edge forces
+
+ for (ni = 0; ni < npi; ni++) {
+ discrete[ifirst+ni][3] = 0;
+ discrete[ifirst+ni][4] = 0;
+ discrete[ifirst+ni][5] = 0;
+ discrete[ifirst+ni][6] = 0;
+ }
+
+ for (ni = 0; ni < nei; ni++) {
+ edge[iefirst+ni][2] = 0;
+ edge[iefirst+ni][3] = 0;
+ edge[iefirst+ni][4] = 0;
+ edge[iefirst+ni][5] = 0;
+ }
+
+ int interact, num_contacts, done;
+ double delta_a, delta_ua, j_a;
+ Contact contact_list[MAX_CONTACTS];
+
+ num_contacts = 0;
+ facc[0] = facc[1] = facc[2] = 0;
+ interact = edge_against_wall(i, wall_pos, side, vwall, x, f, torque,
+ contact_list, num_contacts, facc);
+
+ } // group bit
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallBodyPolyhedron::reset_dt()
+{
+ dt = update->dt;
+}
+
+/* ----------------------------------------------------------------------
+ convert N sub-particles in body I to space frame using current quaternion
+ store sub-particle space-frame displacements from COM in discrete list
+------------------------------------------------------------------------- */
+
+void FixWallBodyPolyhedron::body2space(int i)
+{
+ int ibonus = atom->body[i];
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+ int nsub = bptr->nsub(bonus);
+ double *coords = bptr->coords(bonus);
+ int body_num_edges = bptr->nedges(bonus);
+ double* edge_ends = bptr->edges(bonus);
+ int body_num_faces = bptr->nfaces(bonus);
+ double* face_pts = bptr->faces(bonus);
+ double eradius = bptr->enclosing_radius(bonus);
+ double rradius = bptr->rounded_radius(bonus);
+
+ // get the number of sub-particles (vertices)
+ // and the index of the first vertex of my body in the list
+
+ dnum[i] = nsub;
+ dfirst[i] = ndiscrete;
+
+ // grow the vertex list if necessary
+ // the first 3 columns are for coords, the last 3 for forces
+
+ if (ndiscrete + nsub > dmax) {
+ dmax += DELTA;
+ memory->grow(discrete,dmax,7,"fix:discrete");
+ }
+
+ double p[3][3];
+ MathExtra::quat_to_mat(bonus->quat,p);
+
+ for (int m = 0; m < nsub; m++) {
+ MathExtra::matvec(p,&coords[3*m],discrete[ndiscrete]);
+ discrete[ndiscrete][3] = 0;
+ discrete[ndiscrete][4] = 0;
+ discrete[ndiscrete][5] = 0;
+ discrete[ndiscrete][6] = 0;
+ ndiscrete++;
+ }
+
+ // get the number of edges (vertices)
+ // and the index of the first edge of my body in the list
+
+ ednum[i] = body_num_edges;
+ edfirst[i] = nedge;
+
+ // grow the edge list if necessary
+ // the first 2 columns are for vertex indices within body,
+ // the last 3 for forces
+
+ if (nedge + body_num_edges > edmax) {
+ edmax += DELTA;
+ memory->grow(edge,edmax,6,"fix:edge");
+ }
+
+ for (int m = 0; m < body_num_edges; m++) {
+ edge[nedge][0] = static_cast<int>(edge_ends[2*m+0]);
+ edge[nedge][1] = static_cast<int>(edge_ends[2*m+1]);
+ edge[nedge][2] = 0;
+ edge[nedge][3] = 0;
+ edge[nedge][4] = 0;
+ edge[nedge][5] = 0;
+ nedge++;
+ }
+
+ // get the number of faces and the index of the first face
+
+ facnum[i] = body_num_faces;
+ facfirst[i] = nface;
+
+ // grow the face list if necessary
+ // the first 3 columns are for vertex indices within body, the last 3 for forces
+
+ if (nface + body_num_faces > facmax) {
+ facmax += DELTA;
+ memory->grow(face,facmax,6,"pair:face");
+ }
+
+ for (int m = 0; m < body_num_faces; m++) {
+ face[nface][0] = static_cast<int>(face_pts[3*m+0]);
+ face[nface][1] = static_cast<int>(face_pts[3*m+1]);
+ face[nface][2] = static_cast<int>(face_pts[3*m+2]);
+ face[nface][3] = 0;
+ face[nface][4] = 0;
+ face[nface][5] = 0;
+ nface++;
+ }
+
+ enclosing_radius[i] = eradius;
+ rounded_radius[i] = rradius;
+}
+
+/* ----------------------------------------------------------------------
+ Determine the interaction mode between a sphere against the wall
+
+ i = atom i (body i)
+ x = atoms' coordinates
+ f = atoms' forces
+ torque = atoms' torques
+---------------------------------------------------------------------- */
+
+int FixWallBodyPolyhedron::sphere_against_wall(int i, double wall_pos,
+ int side, double* vwall, double** x, double** v, double** f,
+ double** angmom, double** torque)
+{
+ int mode;
+ double rradi,hi[3],d,delx,dely,delz,R,fpair,fx,fy,fz;
+
+ rradi = rounded_radius[i];
+ mode = NONE;
+
+ if (wallstyle == XPLANE) {
+ hi[0] = wall_pos;
+ hi[1] = x[i][1];
+ hi[2] = x[i][2];
+ } else if (wallstyle == YPLANE) {
+ hi[0] = x[i][0];
+ hi[1] = wall_pos;
+ hi[2] = x[i][2];
+ } else if (wallstyle == ZPLANE) {
+ hi[0] = x[i][0];
+ hi[1] = x[i][1];
+ hi[2] = wall_pos;
+ }
+
+ distance(hi, x[i], d);
+
+ if (d <= rradi) {
+ delx = x[i][0] - hi[0];
+ dely = x[i][1] - hi[1];
+ delz = x[i][2] - hi[2];
+ R = d - rradi;
+
+ fpair = -kn * R;
+
+ fx = delx*fpair/d;
+ fy = dely*fpair/d;
+ fz = delz*fpair/d;
+
+ contact_forces(i, 1.0, x[i], hi, delx, dely, delz,
+ fx, fy, fz, x, v, angmom, f, torque, vwall);
+ mode = VERTEX;
+ }
+
+ return mode;
+}
+
+/* ----------------------------------------------------------------------
+ Determine the interaction mode between i's vertices against the wall
+
+ i = atom i (body i)
+ x = atoms' coordinates
+ f = atoms' forces
+ torque = atoms' torques
+ Output:
+ contact_list = list of contacts between i and the wall
+ num_contacts = number of contacts between i's vertices and the wall
+ Return:
+ number of contacts of the edge to the wall (0, 1 or 2)
+---------------------------------------------------------------------- */
+
+int FixWallBodyPolyhedron::edge_against_wall(int i, double wall_pos,
+ int side, double* vwall, double** x, double** f, double** torque,
+ Contact* contact_list, int &num_contacts, double* facc)
+{
+ int ni, nei, mode, contact;
+ double rradi;
+ int nlocal = atom->nlocal;
+
+ nei = ednum[i];
+ rradi = rounded_radius[i];
+
+ contact = 0;
+
+ // loop through body i's edges
+
+ for (ni = 0; ni < nei; ni++)
+ mode = compute_distance_to_wall(i, ni, x[i], rradi, wall_pos, side, vwall,
+ contact);
+
+ return contact;
+}
+
+/* -------------------------------------------------------------------------
+ Compute the distance between a vertex to the wall
+ another body
+ Input:
+ x0 = coordinate of the tested vertex
+ rradi = rounded radius of the vertex
+ wall_pos = position of the wall
+ Output:
+ d = Distance from a point x0 to an wall
+ hi = coordinates of the projection of x0 on the wall
+ contact = 0 no contact between the queried vertex and the wall
+ 1 contact detected
+ return NONE if there is no interaction
+ VERTEX if the tested vertex interacts with the wall
+------------------------------------------------------------------------- */
+
+int FixWallBodyPolyhedron::compute_distance_to_wall(int ibody, int edge_index,
+ double *xmi, double rounded_radius_i, double wall_pos,
+ int side, double* vwall, int &contact)
+{
+ int mode,ifirst,iefirst,npi1,npi2;
+ double d1,d2,xpi1[3],xpi2[3],hi[3];
+ double fx,fy,fz,fpair,delx,dely,delz,R;
+
+ double** x = atom->x;
+ double** v = atom->v;
+ double** f = atom->f;
+ double** torque = atom->torque;
+ double** angmom = atom->angmom;
+
+ // two ends of the edge from body i
+
+ ifirst = dfirst[ibody];
+ iefirst = edfirst[ibody];
+ npi1 = static_cast<int>(edge[iefirst+edge_index][0]);
+ npi2 = static_cast<int>(edge[iefirst+edge_index][1]);
+
+ xpi1[0] = xmi[0] + discrete[ifirst+npi1][0];
+ xpi1[1] = xmi[1] + discrete[ifirst+npi1][1];
+ xpi1[2] = xmi[2] + discrete[ifirst+npi1][2];
+
+ xpi2[0] = xmi[0] + discrete[ifirst+npi2][0];
+ xpi2[1] = xmi[1] + discrete[ifirst+npi2][1];
+ xpi2[2] = xmi[2] + discrete[ifirst+npi2][2];
+
+ // determine the intersection of the edge to the wall
+
+ mode = NONE;
+ double j_a = 1.0;
+
+ if (wallstyle == XPLANE) {
+ hi[0] = wall_pos;
+ hi[1] = xpi1[1];
+ hi[2] = xpi1[2];
+ } else if (wallstyle == YPLANE) {
+ hi[0] = xpi1[0];
+ hi[1] = wall_pos;
+ hi[2] = xpi1[2];
+ } else if (wallstyle == ZPLANE) {
+ hi[0] = xpi1[0];
+ hi[1] = xpi1[1];
+ hi[2] = wall_pos;
+ }
+
+ distance(hi, xpi1, d1);
+
+ if (d1 <= rounded_radius_i && static_cast<int>(discrete[ifirst+npi1][6]) == 0) {
+ delx = xpi1[0] - hi[0];
+ dely = xpi1[1] - hi[1];
+ delz = xpi1[2] - hi[2];
+ R = d1 - rounded_radius_i;
+
+ fpair = -kn * R;
+
+ fx = delx*fpair/d1;
+ fy = dely*fpair/d1;
+ fz = delz*fpair/d1;
+
+ contact_forces(ibody, j_a, xpi1, hi, delx, dely, delz,
+ fx, fy, fz, x, v, angmom, f, torque, vwall);
+ discrete[ifirst+npi1][6] = 1;
+ contact++;
+ mode = VERTEX;
+ }
+
+ if (wallstyle == XPLANE) {
+ hi[0] = wall_pos;
+ hi[1] = xpi2[1];
+ hi[2] = xpi2[2];
+ } else if (wallstyle == YPLANE) {
+ hi[0] = xpi2[0];
+ hi[1] = wall_pos;
+ hi[2] = xpi2[2];
+ } else if (wallstyle == ZPLANE) {
+ hi[0] = xpi2[0];
+ hi[1] = xpi2[1];
+ hi[2] = wall_pos;
+ }
+
+ distance(hi, xpi2, d2);
+
+ if (d2 <= rounded_radius_i && static_cast<int>(discrete[ifirst+npi2][6]) == 0) {
+ delx = xpi2[0] - hi[0];
+ dely = xpi2[1] - hi[1];
+ delz = xpi2[2] - hi[2];
+ R = d2 - rounded_radius_i;
+
+ fpair = -kn * R;
+
+ fx = delx*fpair/d2;
+ fy = dely*fpair/d2;
+ fz = delz*fpair/d2;
+
+ contact_forces(ibody, j_a, xpi2, hi, delx, dely, delz,
+ fx, fy, fz, x, v, angmom, f, torque, vwall);
+ discrete[ifirst+npi2][6] = 1;
+ contact++;
+ mode = VERTEX;
+ }
+
+ return mode;
+}
+
+/* ----------------------------------------------------------------------
+ Compute contact forces between two bodies
+ modify the force stored at the vertex and edge in contact by j_a
+ sum forces and torque to the corresponding bodies
+ fn = normal friction component
+ ft = tangential friction component (-c_t * v_t)
+------------------------------------------------------------------------- */
+
+void FixWallBodyPolyhedron::contact_forces(int ibody,
+ double j_a, double *xi, double *xj, double delx, double dely, double delz,
+ double fx, double fy, double fz, double** x, double** v, double** angmom,
+ double** f, double** torque, double* vwall)
+{
+ int ibonus,jbonus;
+ double fxt,fyt,fzt,rsq,rsqinv;
+ double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
+ double fn[3],ft[3],vi[3],vj[3];
+ double *quat, *inertia;
+ AtomVecBody::Bonus *bonus;
+
+ // compute the velocity of the vertex in the space-fixed frame
+
+ ibonus = atom->body[ibody];
+ bonus = &avec->bonus[ibonus];
+ quat = bonus->quat;
+ inertia = bonus->inertia;
+ total_velocity(xi, x[ibody], v[ibody], angmom[ibody],
+ inertia, quat, vi);
+
+ // vector pointing from the contact point on ibody to the wall
+
+ rsq = delx*delx + dely*dely + delz*delz;
+ rsqinv = 1.0/rsq;
+
+ // relative translational velocity
+
+ vr1 = vi[0] - vwall[0];
+ vr2 = vi[1] - vwall[1];
+ vr3 = vi[2] - vwall[2];
+
+ // normal component
+
+ vnnr = vr1*delx + vr2*dely + vr3*delz;
+ vn1 = delx*vnnr * rsqinv;
+ vn2 = dely*vnnr * rsqinv;
+ vn3 = delz*vnnr * rsqinv;
+
+ // tangential component
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // normal friction term at contact
+
+ fn[0] = -c_n * vn1;
+ fn[1] = -c_n * vn2;
+ fn[2] = -c_n * vn3;
+
+ // tangential friction term at contact
+ // excluding the tangential deformation term for now
+
+ ft[0] = -c_t * vt1;
+ ft[1] = -c_t * vt2;
+ ft[2] = -c_t * vt3;
+
+ fxt = fx; fyt = fy; fzt = fz;
+ fx = fxt * j_a + fn[0] + ft[0];
+ fy = fyt * j_a + fn[1] + ft[1];
+ fz = fzt * j_a + fn[2] + ft[2];
+
+ f[ibody][0] += fx;
+ f[ibody][1] += fy;
+ f[ibody][2] += fz;
+ sum_torque(x[ibody], xi, fx, fy, fz, torque[ibody]);
+
+ #ifdef _POLYHEDRON_DEBUG
+ printf("From contact forces: vertex fx %f fy %f fz %f\n"
+ " torque body %d: %f %f %f\n"
+ " torque body %d: %f %f %f\n",
+ fxt, fyt, fzt,
+ atom->tag[ibody],torque[ibody][0],torque[ibody][1],torque[ibody][2],
+ atom->tag[jbody],torque[jbody][0],torque[jbody][1],torque[jbody][2]);
+ #endif
+}
+
+/* ----------------------------------------------------------------------
+ Compute the contact forces between two bodies
+ modify the force stored at the vertex and edge in contact by j_a
+ sum forces and torque to the corresponding bodies
+ fn = normal friction component
+ ft = tangential friction component (-c_t * vrt)
+------------------------------------------------------------------------- */
+
+void FixWallBodyPolyhedron::contact_forces(Contact& contact, double j_a,
+ double** x, double** v, double** angmom, double** f,
+ double** torque, double* vwall, double* facc)
+{
+ int ibody,ibonus,ifirst, jefirst, ni;
+ double fx,fy,fz,delx,dely,delz,rsq,rsqinv;
+ double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
+ double fn[3],ft[3],vi[3];
+ double *quat, *inertia;
+ AtomVecBody::Bonus *bonus;
+
+ ibody = contact.ibody;
+
+ // compute the velocity of the vertex in the space-fixed frame
+
+ ibonus = atom->body[ibody];
+ bonus = &avec->bonus[ibonus];
+ quat = bonus->quat;
+ inertia = bonus->inertia;
+ total_velocity(contact.xv, x[ibody], v[ibody], angmom[ibody],
+ inertia, quat, vi);
+
+ // vector pointing from the vertex to the point on the wall
+
+ delx = contact.xv[0] - contact.xe[0];
+ dely = contact.xv[1] - contact.xe[1];
+ delz = contact.xv[2] - contact.xe[2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ rsqinv = 1.0/rsq;
+
+ // relative translational velocity
+
+ vr1 = vi[0] - vwall[0];
+ vr2 = vi[1] - vwall[1];
+ vr3 = vi[2] - vwall[2];
+
+ // normal component
+
+ vnnr = vr1*delx + vr2*dely + vr3*delz;
+ vn1 = delx*vnnr * rsqinv;
+ vn2 = dely*vnnr * rsqinv;
+ vn3 = delz*vnnr * rsqinv;
+
+ // tangential component
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // normal friction term at contact
+
+ fn[0] = -c_n * vn1;
+ fn[1] = -c_n * vn2;
+ fn[2] = -c_n * vn3;
+
+ // tangential friction term at contact
+ // excluding the tangential deformation term for now
+
+ ft[0] = -c_t * vt1;
+ ft[1] = -c_t * vt2;
+ ft[2] = -c_t * vt3;
+
+ // only the cohesive force is scaled by j_a
+
+ ifirst = dfirst[ibody];
+ ni = contact.vertex;
+
+ fx = discrete[ifirst+ni][3] * j_a + fn[0] + ft[0];
+ fy = discrete[ifirst+ni][4] * j_a + fn[1] + ft[1];
+ fz = discrete[ifirst+ni][5] * j_a + fn[2] + ft[2];
+ f[ibody][0] += fx;
+ f[ibody][1] += fy;
+ f[ibody][2] += fz;
+ sum_torque(x[ibody], contact.xv, fx, fy, fz, torque[ibody]);
+
+ // accumulate forces to the vertex only
+
+ facc[0] += fx; facc[1] += fy; facc[2] += fz;
+
+ #ifdef _POLYHEDRON_DEBUG
+ printf("From contact forces: vertex fx %f fy %f fz %f\n"
+ " torque body %d: %f %f %f\n",
+ discrete[ifirst+ni][3], discrete[ifirst+ni][4], discrete[ifirst+ni][5],
+ atom->tag[ibody],torque[ibody][0],torque[ibody][1],torque[ibody][2]);
+ #endif
+}
+
+/* ----------------------------------------------------------------------
+ Determine the length of the contact segment, i.e. the separation between
+ 2 contacts
+------------------------------------------------------------------------- */
+
+double FixWallBodyPolyhedron::contact_separation(const Contact& c1,
+ const Contact& c2)
+{
+ double x1 = c1.xv[0];
+ double y1 = c1.xv[1];
+ double x2 = c1.xe[0];
+ double y2 = c1.xe[1];
+ double x3 = c2.xv[0];
+ double y3 = c2.xv[1];
+
+ double delta_a = 0.0;
+ if (fabs(x2 - x1) > EPSILON) {
+ double A = (y2 - y1) / (x2 - x1);
+ delta_a = fabs(y1 - A * x1 - y3 + A * x3) / sqrt(1 + A * A);
+ } else {
+ delta_a = fabs(x1 - x3);
+ }
+
+ return delta_a;
+}
+
+/* ----------------------------------------------------------------------
+ Accumulate torque to body from the force f=(fx,fy,fz) acting at point x
+------------------------------------------------------------------------- */
+
+void FixWallBodyPolyhedron::sum_torque(double* xm, double *x, double fx,
+ double fy, double fz, double* torque)
+{
+ double rx = x[0] - xm[0];
+ double ry = x[1] - xm[1];
+ double rz = x[2] - xm[2];
+ double tx = ry * fz - rz * fy;
+ double ty = rz * fx - rx * fz;
+ double tz = rx * fy - ry * fx;
+ torque[0] += tx;
+ torque[1] += ty;
+ torque[2] += tz;
+}
+
+/* ----------------------------------------------------------------------
+ Calculate the total velocity of a point (vertex, a point on an edge):
+ vi = vcm + omega ^ (p - xcm)
+------------------------------------------------------------------------- */
+
+void FixWallBodyPolyhedron::total_velocity(double* p, double *xcm,
+ double* vcm, double *angmom, double *inertia,
+ double *quat, double* vi)
+{
+ double r[3],omega[3],ex_space[3],ey_space[3],ez_space[3];
+ r[0] = p[0] - xcm[0];
+ r[1] = p[1] - xcm[1];
+ r[2] = p[2] - xcm[2];
+ MathExtra::q_to_exyz(quat,ex_space,ey_space,ez_space);
+ MathExtra::angmom_to_omega(angmom,ex_space,ey_space,ez_space,
+ inertia,omega);
+ vi[0] = omega[1]*r[2] - omega[2]*r[1] + vcm[0];
+ vi[1] = omega[2]*r[0] - omega[0]*r[2] + vcm[1];
+ vi[2] = omega[0]*r[1] - omega[1]*r[0] + vcm[2];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallBodyPolyhedron::distance(const double* x2, const double* x1,
+ double& r) {
+ r = sqrt((x2[0] - x1[0]) * (x2[0] - x1[0])
+ + (x2[1] - x1[1]) * (x2[1] - x1[1])
+ + (x2[2] - x1[2]) * (x2[2] - x1[2]));
+}
+
diff --git a/src/BODY/fix_wall_body_polyhedron.h b/src/BODY/fix_wall_body_polyhedron.h
new file mode 100644
index 0000000000000000000000000000000000000000..ff7b7ca7cfd0806de2d5f4fcd16c79b48c980954
--- /dev/null
+++ b/src/BODY/fix_wall_body_polyhedron.h
@@ -0,0 +1,143 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(wall/body/polyhedron,FixWallBodyPolyhedron)
+
+#else
+
+#ifndef LMP_FIX_WALL_BODY_POLYHERON_H
+#define LMP_FIX_WALL_BODY_POLYHERON_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixWallBodyPolyhedron : public Fix {
+ public:
+ FixWallBodyPolyhedron(class LAMMPS *, int, char **);
+ virtual ~FixWallBodyPolyhedron();
+ int setmask();
+ void init();
+ void setup(int);
+ virtual void post_force(int);
+ void reset_dt();
+
+ struct Contact {
+ int ibody, jbody; // body (i.e. atom) indices (not tags)
+ int vertex; // vertex of the first polygon
+ int edge; // edge of the second polygon
+ double xv[3]; // coordinates of the vertex
+ double xe[3]; // coordinates of the projection of the vertex on the edge
+ double separation;// separation at contact
+ };
+
+ protected:
+ int wallstyle,pairstyle,wiggle,axis;
+ double kn,c_n,c_t;
+ double lo,hi,cylradius;
+ double amplitude,period,omega;
+ double dt;
+ int time_origin;
+
+ class AtomVecBody *avec;
+ class BodyRoundedPolyhedron *bptr;
+
+ double **discrete; // list of all sub-particles for all bodies
+ int ndiscrete; // number of discretes in list
+ int dmax; // allocated size of discrete list
+ int *dnum; // number of discretes per line, 0 if uninit
+ int *dfirst; // index of first discrete per each line
+ int nmax; // allocated size of dnum,dfirst vectors
+
+ double **edge; // list of all edge for all bodies
+ int nedge; // number of edge in list
+ int edmax; // allocated size of edge list
+ int *ednum; // number of edges per line, 0 if uninit
+ int *edfirst; // index of first edge per each line
+ int ednummax; // allocated size of ednum,edfirst vectors
+
+ double **face; // list of all edge for all bodies
+ int nface; // number of faces in list
+ int facmax; // allocated size of face list
+ int *facnum; // number of faces per line, 0 if uninit
+ int *facfirst; // index of first face per each line
+ int facnummax; // allocated size of facnum,facfirst vectors
+
+ double *enclosing_radius; // enclosing radii for all bodies
+ double *rounded_radius; // rounded radii for all bodies
+
+ void body2space(int);
+
+ int edge_against_wall(int ibody, double wall_pos, int side, double* vwall,
+ double** x, double** f, double** torque, Contact* contact_list,
+ int &num_contacts, double* facc);
+ int sphere_against_wall(int i, double wall_pos, int side, double* vwall,
+ double** x, double** v, double** f, double** angmom, double** torque);
+
+ int compute_distance_to_wall(int ibody, int edge_index, double *xmi,
+ double rounded_radius_i, double wall_pos, int side,
+ double* vwall, int &contact);
+ double contact_separation(const Contact& c1, const Contact& c2);
+ void contact_forces(int ibody, double j_a, double *xi, double *xj,
+ double delx, double dely, double delz,
+ double fx, double fy, double fz, double** x, double** v,
+ double** angmom, double** f, double** torque, double* vwall);
+
+ void contact_forces(Contact& contact, double j_a, double** x,
+ double** v, double** angmom, double** f, double** torque,
+ double* vwall, double* facc);
+ void sum_torque(double* xm, double *x, double fx,
+ double fy, double fz, double* torque);
+ void total_velocity(double* p, double *xcm, double* vcm, double *angmom,
+ double *inertia, double *quat, double* vi);
+ void distance(const double* x2, const double* x1, double& r);
+
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Fix wall/body/polyhedron requires atom style body rounded/polyhedron
+
+Self-explanatory.
+
+E: Cannot use wall in periodic dimension
+
+Self-explanatory.
+
+E: Cannot wiggle and shear fix wall/body/polygon
+
+Cannot specify both options at the same time.
+
+E: Invalid wiggle direction for fix wall/body/polygon
+
+Self-explanatory.
+
+E: Fix wall/body/polygon is incompatible with Pair style
+
+Must use a body pair style to define the parameters needed for
+this fix.
+
+*/
diff --git a/src/BODY/pair_body_rounded_polygon.cpp b/src/BODY/pair_body_rounded_polygon.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..72591d2bd04e3d2534ff5c87e4765c6a921d9580
--- /dev/null
+++ b/src/BODY/pair_body_rounded_polygon.cpp
@@ -0,0 +1,1359 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)
+ Ref: Fraige, Langston, Matchett and Dodds, Particuology 2008, 6:455-466
+ Note: The current implementation has not taken into account
+ the contact history for friction forces.
+------------------------------------------------------------------------- */
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include "pair_body_rounded_polygon.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "atom_vec_body.h"
+#include "body_rounded_polygon.h"
+#include "comm.h"
+#include "force.h"
+#include "fix.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+//#define _POLYGON_DEBUG
+#define DELTA 10000
+#define EPSILON 1e-3
+#define MAX_CONTACTS 4 // maximum number of contacts for 2D models
+#define EFF_CONTACTS 2 // effective contacts for 2D models
+
+enum {INVALID=0,NONE=1,VERTEXI=2,VERTEXJ=3,EDGE=4};
+
+/* ---------------------------------------------------------------------- */
+
+PairBodyRoundedPolygon::PairBodyRoundedPolygon(LAMMPS *lmp) : Pair(lmp)
+{
+ dmax = nmax = 0;
+ discrete = NULL;
+ dnum = dfirst = NULL;
+
+ edmax = ednummax = 0;
+ edge = NULL;
+ ednum = edfirst = NULL;
+
+ enclosing_radius = NULL;
+ rounded_radius = NULL;
+ maxerad = NULL;
+
+ single_enable = 0;
+ restartinfo = 0;
+
+ c_n = 0.1;
+ c_t = 0.2;
+ mu = 0.0;
+ delta_ua = 1.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairBodyRoundedPolygon::~PairBodyRoundedPolygon()
+{
+ memory->destroy(discrete);
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+
+ memory->destroy(edge);
+ memory->destroy(ednum);
+ memory->destroy(edfirst);
+
+ memory->destroy(enclosing_radius);
+ memory->destroy(rounded_radius);
+
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+
+ memory->destroy(k_n);
+ memory->destroy(k_na);
+ memory->destroy(maxerad);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ int ni,nj,npi,npj,ifirst,jfirst;
+ int nei,nej,iefirst,jefirst;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fx,fy,fz;
+ double rsq,rsqinv,r,radi,radj,eradi,eradj,rradi,rradj,k_nij,k_naij;
+ double xi[3],xj[3],fi[3],fj[3],ti[3],tj[3],facc[3];
+ double *dxi,*dxj;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **v = atom->v;
+ double **f = atom->f;
+ double **torque = atom->torque;
+ double **angmom = atom->angmom;
+ double *radius = atom->radius;
+ tagint* tag = atom->tag;
+ int *body = atom->body;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ int nall = nlocal + atom->nghost;
+ int newton_pair = force->newton_pair;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // grow the per-atom lists if necessary and initialize
+
+ if (atom->nmax > nmax) {
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+ memory->destroy(ednum);
+ memory->destroy(edfirst);
+ memory->destroy(enclosing_radius);
+ memory->destroy(rounded_radius);
+ nmax = atom->nmax;
+ memory->create(dnum,nmax,"pair:dnum");
+ memory->create(dfirst,nmax,"pair:dfirst");
+ memory->create(ednum,nmax,"pair:ednum");
+ memory->create(edfirst,nmax,"pair:edfirst");
+ memory->create(enclosing_radius,nmax,"pair:enclosing_radius");
+ memory->create(rounded_radius,nmax,"pair:rounded_radius");
+ }
+
+ ndiscrete = nedge = 0;
+ for (i = 0; i < nall; i++)
+ dnum[i] = ednum[i] = 0;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ radi = radius[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ if (body[i] >= 0) {
+ if (dnum[i] == 0) body2space(i);
+ npi = dnum[i];
+ ifirst = dfirst[i];
+ nei = ednum[i];
+ iefirst = edfirst[i];
+ eradi = enclosing_radius[i];
+ rradi = rounded_radius[i];
+ }
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+ radj = radius[j];
+
+ // body/body interactions
+
+ evdwl = 0.0;
+ facc[0] = facc[1] = facc[2] = 0;
+
+ if (body[i] < 0 || body[j] < 0) continue;
+
+ if (dnum[j] == 0) body2space(j);
+ npj = dnum[j];
+ jfirst = dfirst[j];
+ nej = ednum[j];
+ jefirst = edfirst[j];
+ eradj = enclosing_radius[j];
+ rradj = rounded_radius[j];
+
+ k_nij = k_n[itype][jtype];
+ k_naij = k_na[itype][jtype];
+
+ // no interaction
+
+ r = sqrt(rsq);
+ if (r > radi + radj + cut_inner) continue;
+ rsqinv = 1.0 / rsq;
+
+ if (npi == 1 && npj == 1) {
+ sphere_against_sphere(i, j, delx, dely, delz, rsq,
+ k_nij, k_naij, x, v, f, evflag);
+ continue;
+ }
+
+ // reset vertex and edge forces
+
+ for (ni = 0; ni < npi; ni++) {
+ discrete[ifirst+ni][3] = 0;
+ discrete[ifirst+ni][4] = 0;
+ discrete[ifirst+ni][5] = 0;
+ }
+
+ for (nj = 0; nj < npj; nj++) {
+ discrete[jfirst+nj][3] = 0;
+ discrete[jfirst+nj][4] = 0;
+ discrete[jfirst+nj][5] = 0;
+ }
+
+ for (ni = 0; ni < nei; ni++) {
+ edge[iefirst+ni][2] = 0;
+ edge[iefirst+ni][3] = 0;
+ edge[iefirst+ni][4] = 0;
+ }
+
+ for (nj = 0; nj < nej; nj++) {
+ edge[jefirst+nj][2] = 0;
+ edge[jefirst+nj][3] = 0;
+ edge[jefirst+nj][4] = 0;
+ }
+
+ int interact, num_contacts, done;
+ double delta_a, j_a;
+ Contact contact_list[MAX_CONTACTS];
+
+ num_contacts = 0;
+
+ // check interaction between i's vertices and j' edges
+
+ interact = vertex_against_edge(i, j, k_nij, k_naij,
+ x, f, torque, tag, contact_list,
+ num_contacts, evdwl, facc);
+
+ // check interaction between j's vertices and i' edges
+
+ interact = vertex_against_edge(j, i, k_nij, k_naij,
+ x, f, torque, tag, contact_list,
+ num_contacts, evdwl, facc);
+
+ if (num_contacts >= 2) {
+
+ // find the first two distinct contacts
+
+ done = 0;
+ for (int m = 0; m < num_contacts-1; m++) {
+ for (int n = m+1; n < num_contacts; n++) {
+ delta_a = contact_separation(contact_list[m], contact_list[n]);
+ if (delta_a > 0) {
+ j_a = delta_a / (EFF_CONTACTS * delta_ua);
+ if (j_a < 1.0) j_a = 1.0;
+
+ // scale the force at both contacts
+
+ contact_forces(contact_list[m], j_a, x, v, angmom, f, torque, evdwl, facc);
+ contact_forces(contact_list[n], j_a, x, v, angmom, f, torque, evdwl, facc);
+ done = 1;
+
+ #ifdef _POLYGON_DEBUG
+ printf(" Two separate contacts %d and %d: delta_a = %f; j_a = %f\n",
+ m, n, delta_a, j_a);
+ printf(" %d: vertex %d of body %d and edge %d of body %d; "
+ "xv = %f %f %f; xe = %f %f %f\n",
+ m, contact_list[m].vertex, contact_list[m].ibody,
+ contact_list[m].edge, contact_list[m].jbody,
+ contact_list[m].xv[0], contact_list[m].xv[1], contact_list[m].xv[2],
+ contact_list[m].xe[0], contact_list[m].xe[1], contact_list[m].xe[2]);
+ printf(" %d: vertex %d of body %d and edge %d of body %d; "
+ "xv = %f %f %f; xe = %f %f %f\n",
+ n, contact_list[n].vertex, contact_list[n].ibody,
+ contact_list[n].edge, contact_list[n].jbody,
+ contact_list[n].xv[0], contact_list[n].xv[1], contact_list[n].xv[2],
+ contact_list[n].xe[0], contact_list[n].xe[1], contact_list[n].xe[2]);
+ #endif
+
+ break;
+ }
+ }
+ if (done == 1) break;
+ }
+
+
+ } else if (num_contacts == 1) {
+
+ // if there's only one contact, it should be handled here
+ // since forces/torques have not been accumulated from vertex2edge()
+
+ contact_forces(contact_list[0], 1.0, x, v, angmom, f, torque, evdwl, facc);
+
+ #ifdef _POLYGON_DEBUG
+ printf("One contact between vertex %d of body %d and edge %d of body %d:\n",
+ contact_list[0].vertex, tag[contact_list[0].ibody],
+ contact_list[0].edge, tag[contact_list[0].jbody]);
+ printf("xv = %f %f %f; xe = %f %f %f\n",
+ contact_list[0].xv[0], contact_list[0].xv[1], contact_list[0].xv[2],
+ contact_list[0].xe[0], contact_list[0].xe[1], contact_list[0].xe[2]);
+ #endif
+ }
+
+ #ifdef _POLYGON_DEBUG
+ int num_overlapping_contacts = 0;
+ for (int m = 0; m < num_contacts-1; m++) {
+ for (int n = m+1; n < num_contacts; n++) {
+ double l = contact_separation(contact_list[m], contact_list[n]);
+ if (l < EPSILON) num_overlapping_contacts++;
+ }
+ }
+ printf("There are %d contacts detected, %d of which overlap.\n",
+ num_contacts, num_overlapping_contacts);
+ #endif
+
+ if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0.0,
+ facc[0],facc[1],facc[2],delx,dely,delz);
+
+ } // end for jj
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+ memory->create(k_n,n+1,n+1,"pair:k_n");
+ memory->create(k_na,n+1,n+1,"pair:k_na");
+ memory->create(maxerad,n+1,"pair:maxerad");
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::settings(int narg, char **arg)
+{
+ if (narg < 5) error->all(FLERR,"Illegal pair_style command");
+
+ c_n = force->numeric(FLERR,arg[0]);
+ c_t = force->numeric(FLERR,arg[1]);
+ mu = force->numeric(FLERR,arg[2]);
+ delta_ua = force->numeric(FLERR,arg[3]);
+ cut_inner = force->numeric(FLERR,arg[4]);
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::coeff(int narg, char **arg)
+{
+ if (narg < 4 || narg > 5)
+ error->all(FLERR,"Incorrect args for pair coefficients");
+ if (!allocated) allocate();
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
+
+ double k_n_one = force->numeric(FLERR,arg[2]);
+ double k_na_one = force->numeric(FLERR,arg[3]);
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ k_n[i][j] = k_n_one;
+ k_na[i][j] = k_na_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::init_style()
+{
+ avec = (AtomVecBody *) atom->style_match("body");
+ if (!avec) error->all(FLERR,"Pair body/rounded/polygon requires atom style body");
+ if (strcmp(avec->bptr->style,"rounded/polygon") != 0)
+ error->all(FLERR,"Pair body/rounded/polygon requires body style rounded/polygon");
+ bptr = (BodyRoundedPolygon *) avec->bptr;
+
+ if (comm->ghost_velocity == 0)
+ error->all(FLERR,"Pair body/rounded/polygon requires ghost atoms store velocity");
+
+ neighbor->request(this);
+
+ // find the maximum enclosing radius for each atom type
+
+ int i, itype;
+ double eradi;
+ int* body = atom->body;
+ int* type = atom->type;
+ int ntypes = atom->ntypes;
+ int nlocal = atom->nlocal;
+
+ if (atom->nmax > nmax) {
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+ memory->destroy(ednum);
+ memory->destroy(edfirst);
+ memory->destroy(enclosing_radius);
+ memory->destroy(rounded_radius);
+ nmax = atom->nmax;
+ memory->create(dnum,nmax,"pair:dnum");
+ memory->create(dfirst,nmax,"pair:dfirst");
+ memory->create(ednum,nmax,"pair:ednum");
+ memory->create(edfirst,nmax,"pair:edfirst");
+ memory->create(enclosing_radius,nmax,"pair:enclosing_radius");
+ memory->create(rounded_radius,nmax,"pair:rounded_radius");
+ }
+
+ ndiscrete = nedge = 0;
+ for (i = 0; i < nlocal; i++)
+ dnum[i] = ednum[i] = 0;
+
+ double *merad = NULL;
+ memory->create(merad,ntypes+1,"pair:merad");
+ for (i = 1; i <= ntypes; i++)
+ maxerad[i] = merad[i] = 0;
+
+ int ipour;
+ for (ipour = 0; ipour < modify->nfix; ipour++)
+ if (strcmp(modify->fix[ipour]->style,"pour") == 0) break;
+ if (ipour == modify->nfix) ipour = -1;
+
+ int idep;
+ for (idep = 0; idep < modify->nfix; idep++)
+ if (strcmp(modify->fix[idep]->style,"deposit") == 0) break;
+ if (idep == modify->nfix) idep = -1;
+
+ for (i = 1; i <= ntypes; i++) {
+ merad[i] = 0.0;
+ if (ipour >= 0) {
+ itype = i;
+ merad[i] =
+ *((double *) modify->fix[ipour]->extract("radius",itype));
+ }
+ if (idep >= 0) {
+ itype = i;
+ merad[i] =
+ *((double *) modify->fix[idep]->extract("radius",itype));
+ }
+ }
+
+ for (i = 0; i < nlocal; i++) {
+ itype = type[i];
+ if (body[i] >= 0) {
+ if (dnum[i] == 0) body2space(i);
+ eradi = enclosing_radius[i];
+ if (eradi > merad[itype]) merad[itype] = eradi;
+ } else
+ merad[itype] = 0;
+ }
+
+ MPI_Allreduce(&merad[1],&maxerad[1],ntypes,MPI_DOUBLE,MPI_MAX,world);
+
+ memory->destroy(merad);
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairBodyRoundedPolygon::init_one(int i, int j)
+{
+ k_n[j][i] = k_n[i][j];
+ k_na[j][i] = k_na[i][j];
+
+ return (maxerad[i]+maxerad[j]);
+}
+
+/* ----------------------------------------------------------------------
+ convert N sub-particles in body I to space frame using current quaternion
+ store sub-particle space-frame displacements from COM in discrete list
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::body2space(int i)
+{
+ int ibonus = atom->body[i];
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+ int nsub = bptr->nsub(bonus);
+ double *coords = bptr->coords(bonus);
+ int body_num_edges = bptr->nedges(bonus);
+ double* edge_ends = bptr->edges(bonus);
+ double eradius = bptr->enclosing_radius(bonus);
+ double rradius = bptr->rounded_radius(bonus);
+
+ // get the number of sub-particles (vertices)
+ // and the index of the first vertex of my body in the list
+
+ dnum[i] = nsub;
+ dfirst[i] = ndiscrete;
+
+ // grow the vertex list if necessary
+ // the first 3 columns are for coords, the last 3 for forces
+
+ if (ndiscrete + nsub > dmax) {
+ dmax += DELTA;
+ memory->grow(discrete,dmax,6,"pair:discrete");
+ }
+
+ double p[3][3];
+ MathExtra::quat_to_mat(bonus->quat,p);
+
+ for (int m = 0; m < nsub; m++) {
+ MathExtra::matvec(p,&coords[3*m],discrete[ndiscrete]);
+ discrete[ndiscrete][3] = 0;
+ discrete[ndiscrete][4] = 0;
+ discrete[ndiscrete][5] = 0;
+ ndiscrete++;
+ }
+
+ // get the number of edges (vertices)
+ // and the index of the first edge of my body in the list
+
+ ednum[i] = body_num_edges;
+ edfirst[i] = nedge;
+
+ // grow the edge list if necessary
+ // the first 2 columns are for vertex indices within body, the last 3 for forces
+
+ if (nedge + body_num_edges > edmax) {
+ edmax += DELTA;
+ memory->grow(edge,edmax,5,"pair:edge");
+ }
+
+ for (int m = 0; m < body_num_edges; m++) {
+ edge[nedge][0] = static_cast<int>(edge_ends[2*m+0]);
+ edge[nedge][1] = static_cast<int>(edge_ends[2*m+1]);
+ edge[nedge][2] = 0;
+ edge[nedge][3] = 0;
+ edge[nedge][4] = 0;
+ nedge++;
+ }
+
+ enclosing_radius[i] = eradius;
+ rounded_radius[i] = rradius;
+}
+
+/* ----------------------------------------------------------------------
+ Interaction between two spheres with different radii
+ according to the 2D model from Fraige et al.
+---------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::sphere_against_sphere(int i, int j,
+ double delx, double dely, double delz, double rsq,
+ double k_n, double k_na, double** x, double** v,
+ double** f, int evflag)
+{
+ double eradi,eradj,rradi,rradj;
+ double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
+ double rij,rsqinv,R,fx,fy,fz,fn[3],ft[3],fpair,shift,energy;
+ int nlocal = atom->nlocal;
+ int newton_pair = force->newton_pair;
+
+ eradi = enclosing_radius[i];
+ rradi = rounded_radius[i];
+
+ eradj = enclosing_radius[j];
+ rradj = rounded_radius[j];
+
+ rsqinv = 1.0/rsq;
+ rij = sqrt(rsq);
+ R = rij - (rradi + rradj);
+ shift = k_na * cut_inner;
+
+ energy = 0;
+
+ if (R <= 0) { // deformation occurs
+ fpair = -k_n * R - shift;
+ energy = (0.5 * k_n * R + shift) * R;
+ } else if (R <= cut_inner) { // not deforming but cohesive ranges overlap
+ fpair = k_na * R - shift;
+ energy = (-0.5 * k_na * R + shift) * R;
+ } else fpair = 0.0;
+
+ fx = delx*fpair/rij;
+ fy = dely*fpair/rij;
+ fz = delz*fpair/rij;
+
+ if (R <= EPSILON) { // in contact
+
+ // relative translational velocity
+
+ vr1 = v[i][0] - v[j][0];
+ vr2 = v[i][1] - v[j][1];
+ vr3 = v[i][2] - v[j][2];
+
+ // normal component
+
+ vnnr = vr1*delx + vr2*dely + vr3*delz;
+ vn1 = delx*vnnr * rsqinv;
+ vn2 = dely*vnnr * rsqinv;
+ vn3 = delz*vnnr * rsqinv;
+
+ // tangential component
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // normal friction term at contact
+
+ fn[0] = -c_n * vn1;
+ fn[1] = -c_n * vn2;
+ fn[2] = -c_n * vn3;
+
+ // tangential friction term at contact
+ // excluding the tangential deformation term
+
+ ft[0] = -c_t * vt1;
+ ft[1] = -c_t * vt2;
+ ft[2] = -c_t * vt3;
+ }
+
+ f[i][0] += fx;
+ f[i][1] += fy;
+ f[i][2] += fz;
+
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= fx;
+ f[j][1] -= fy;
+ f[j][2] -= fz;
+ }
+
+ if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
+ energy,0.0,fx,fy,fz,delx,dely,delz);
+}
+
+/* ----------------------------------------------------------------------
+ Determine the interaction mode between i's vertices against j's edges
+
+ i = atom i (body i)
+ j = atom j (body j)
+ x = atoms' coordinates
+ f = atoms' forces
+ torque = atoms' torques
+ tag = atoms' tags
+ contact_list = list of contacts
+ num_contacts = number of contacts between i's vertices and j's edges
+ Return:
+ interact = 0 no interaction at all
+ 1 there's at least one case where i's vertices interacts
+ with j's edges
+---------------------------------------------------------------------- */
+
+int PairBodyRoundedPolygon::vertex_against_edge(int i, int j,
+ double k_n, double k_na,
+ double** x, double** f,
+ double** torque, tagint* tag,
+ Contact* contact_list,
+ int &num_contacts,
+ double &evdwl, double* facc)
+{
+ int ni, npi, ifirst, nei, iefirst;
+ int nj, npj, jfirst, nej, jefirst;
+ double xpi[3], xpj[3], dist, eradi, eradj, rradi, rradj;
+ double fx, fy, fz, rx, ry, rz, energy;
+ int interact;
+ int nlocal = atom->nlocal;
+
+ npi = dnum[i];
+ ifirst = dfirst[i];
+ nei = ednum[i];
+ iefirst = edfirst[i];
+ eradi = enclosing_radius[i];
+ rradi = rounded_radius[i];
+
+ npj = dnum[j];
+ jfirst = dfirst[j];
+ nej = ednum[j];
+ jefirst = edfirst[j];
+ eradj = enclosing_radius[j];
+ rradj = rounded_radius[j];
+
+ energy = 0;
+ interact = 0;
+
+ // loop through body i's vertices
+
+ for (ni = 0; ni < npi; ni++) {
+
+ // convert body-fixed coordinates to space-fixed, xi
+
+ xpi[0] = x[i][0] + discrete[ifirst+ni][0];
+ xpi[1] = x[i][1] + discrete[ifirst+ni][1];
+ xpi[2] = x[i][2] + discrete[ifirst+ni][2];
+
+ // compute the distance from the vertex to the COM of body j
+
+ distance(xpi, x[j], dist);
+
+ #ifdef _POLYGON_DEBUG
+ printf("Distance between vertex %d of body %d (%0.1f %0.1f %0.1f) "
+ "to body %d's COM: %f (cut = %0.1f)\n",
+ ni, xpi[0], xpi[1], xpi[2], atom->tag[i], atom->tag[j], dist,
+ eradj + rradi + rradj + cut_inner);
+ #endif
+
+ // the vertex is within the enclosing circle (sphere) of body j,
+ // possibly interacting
+
+ if (dist > eradj + rradj + rradi + cut_inner) continue;
+
+ int mode, contact, p2vertex;
+ double d, R, hi[3], t, delx, dely, delz, fpair, shift;
+ double xj[3], rij;
+
+ // loop through body j's edges
+
+ for (nj = 0; nj < nej; nj++) {
+
+ // compute the distance between the edge nj to the vertex xpi
+
+ mode = compute_distance_to_vertex(j, nj, x[j], rradj,
+ xpi, rradi, cut_inner,
+ d, hi, t, contact);
+
+ if (mode == INVALID || mode == NONE) continue;
+
+ if (mode == VERTEXI || mode == VERTEXJ) {
+
+ interact = 1;
+
+ // vertex i interacts with a vertex of the edge, but does not contact
+
+ if (mode == VERTEXI) p2vertex = edge[jefirst+nj][0];
+ else if (mode == VERTEXJ) p2vertex = edge[jefirst+nj][1];
+
+ // p2.body2space(p2vertex, xj);
+ xpj[0] = x[j][0] + discrete[jfirst+p2vertex][0];
+ xpj[1] = x[j][1] + discrete[jfirst+p2vertex][1];
+ xpj[2] = x[j][2] + discrete[jfirst+p2vertex][2];
+
+ delx = xpi[0] - xpj[0];
+ dely = xpi[1] - xpj[1];
+ delz = xpi[2] - xpj[2];
+
+ // R = surface separation = rij shifted by the rounded radii
+ // R = rij - (p1.rounded_radius + p2.rounded_radius);
+ // note: the force is defined for R, not for rij
+ // R > rc: no interaction between vertex ni and p2vertex
+ // 0 < R < rc: cohesion between vertex ni and p2vertex
+ // R < 0: deformation between vertex ni and p2vertex
+
+ rij = sqrt(delx*delx + dely*dely + delz*delz);
+ R = rij - (rradi + rradj);
+ shift = k_na * cut_inner;
+
+ // the normal frictional term -c_n * vn will be added later
+
+ if (R <= 0) { // deformation occurs
+ fpair = -k_n * R - shift;
+ energy += (0.5 * k_n * R + shift) * R;
+ } else if (R <= cut_inner) { // not deforming but cohesive ranges overlap
+ fpair = k_na * R - shift;
+ energy += (-0.5 * k_na * R + shift) * R;
+ } else fpair = 0.0;
+
+ fx = delx*fpair/rij;
+ fy = dely*fpair/rij;
+ fz = delz*fpair/rij;
+
+ #ifdef _POLYGON_DEBUG
+ printf(" Interaction between vertex %d of %d and vertex %d of %d:",
+ ni, tag[i], p2vertex, tag[j]);
+ printf(" mode = %d; contact = %d; d = %f; rij = %f, t = %f\n",
+ mode, contact, d, rij, t);
+ printf(" R = %f; cut_inner = %f\n", R, cut_inner);
+ printf(" fpair = %f\n", fpair);
+ #endif
+
+ // add forces to body i and body j directly
+ // avoid double counts this pair of vertices
+ // i and j can be either local or ghost atoms (bodies)
+ // probably need more work here when the vertices' interaction
+ // are not symmetric, e.g. j interacts with the edge
+ // consisting of i but in mode = EDGE instead of VERTEX*.
+ // OR, for the time being assume that the edge length is
+ // sufficiently greater than the rounded radius to distinguish
+ // vertex-vertex from vertex-edge contact modes.
+ // Special case: when i is a sphere, also accumulate
+
+ if (tag[i] < tag[j] || npi == 1) {
+
+ f[i][0] += fx;
+ f[i][1] += fy;
+ f[i][2] += fz;
+ sum_torque(x[i], xpi, fx, fy, fz, torque[i]);
+
+ f[j][0] -= fx;
+ f[j][1] -= fy;
+ f[j][2] -= fz;
+ sum_torque(x[j], xpj, -fx, -fy, -fz, torque[j]);
+
+ facc[0] += fx; facc[1] += fy; facc[2] += fz;
+
+ #ifdef _POLYGON_DEBUG
+ printf(" from vertex-vertex: "
+ "force on vertex %d of body %d: fx %f fy %f fz %f\n"
+ " torque body %d: %f %f %f\n"
+ " torque body %d: %f %f %f\n", ni, tag[i], fx, fy, fz,
+ tag[i],torque[i][0],torque[i][1],torque[i][2],
+ tag[j],torque[j][0],torque[j][1],torque[j][2]);
+ #endif
+ }
+
+ // done with the edges from body j,
+ // given that vertex ni interacts with only one vertex from one edge of body j
+
+// break;
+
+ } else if (mode == EDGE) {
+
+ interact = 1;
+
+ // vertex i interacts with the edge
+
+ delx = xpi[0] - hi[0];
+ dely = xpi[1] - hi[1];
+ delz = xpi[2] - hi[2];
+
+ // R = surface separation = d shifted by the rounded radii
+ // R = d - (p1.rounded_radius + p2.rounded_radius);
+ // Note: the force is defined for R, not for d
+ // R > rc: no interaction between vertex i and edge j
+ // 0 < R < rc: cohesion between vertex i and edge j
+ // R < 0: deformation between vertex i and edge j
+ // rij = sqrt(delx*delx + dely*dely + delz*delz);
+
+ R = d - (rradi + rradj);
+ shift = k_na * cut_inner;
+
+ // the normal frictional term -c_n * vn will be added later
+
+ if (R <= 0) { // deformation occurs
+ fpair = -k_n * R - shift;
+ energy += (0.5 * k_n * R + shift) * R;
+ } else if (R <= cut_inner) { // not deforming but cohesive ranges overlap
+ fpair = k_na * R - shift;
+ energy += (-0.5 * k_na * R + shift) * R;
+ } else fpair = 0.0;
+
+ fx = delx*fpair/d;
+ fy = dely*fpair/d;
+ fz = delz*fpair/d;
+
+ #ifdef _POLYGON_DEBUG
+ printf(" Interaction between vertex %d of %d and edge %d of %d:",
+ ni, tag[i], nj, tag[j]);
+ printf(" mode = %d; contact = %d; d = %f; t = %f\n",
+ mode, contact, d, t);
+ printf(" R = %f; cut_inner = %f\n", R, cut_inner);
+ printf(" fpair = %f\n", fpair);
+ #endif
+
+ if (contact == 1) {
+
+ // vertex ni of body i contacts with edge nj of body j
+
+ contact_list[num_contacts].ibody = i;
+ contact_list[num_contacts].jbody = j;
+ contact_list[num_contacts].vertex = ni;
+ contact_list[num_contacts].edge = nj;
+ contact_list[num_contacts].xv[0] = xpi[0];
+ contact_list[num_contacts].xv[1] = xpi[1];
+ contact_list[num_contacts].xv[2] = xpi[2];
+ contact_list[num_contacts].xe[0] = hi[0];
+ contact_list[num_contacts].xe[1] = hi[1];
+ contact_list[num_contacts].xe[2] = hi[2];
+ contact_list[num_contacts].separation = R;
+ num_contacts++;
+
+ // store forces to vertex ni and the edge nj
+ // to be rescaled later
+
+ discrete[ifirst+ni][3] = fx;
+ discrete[ifirst+ni][4] = fy;
+ discrete[ifirst+ni][5] = fz;
+
+ edge[jefirst+nj][2] = -fx;
+ edge[jefirst+nj][3] = -fy;
+ edge[jefirst+nj][4] = -fz;
+
+ #ifdef _POLYGON_DEBUG
+ printf(" Stored forces at vertex and edge for accumulating later.\n");
+ #endif
+
+ } else { // no contact
+
+ // accumulate force and torque to both bodies directly
+
+ f[i][0] += fx;
+ f[i][1] += fy;
+ f[i][2] += fz;
+ sum_torque(x[i], xpi, fx, fy, fz, torque[i]);
+
+ f[j][0] -= fx;
+ f[j][1] -= fy;
+ f[j][2] -= fz;
+ sum_torque(x[j], hi, -fx, -fy, -fz, torque[j]);
+
+ facc[0] += fx; facc[1] += fy; facc[2] += fz;
+
+ #ifdef _POLYGON_DEBUG
+ printf(" from vertex-edge, no contact: "
+ "force on vertex %d of body %d: fx %f fy %f fz %f\n"
+ " torque body %d: %f %f %f\n"
+ " torque body %d: %f %f %f\n", ni, tag[i], fx, fy, fz,
+ tag[i],torque[i][0],torque[i][1],torque[i][2],
+ tag[j],torque[j][0],torque[j][1],torque[j][2]);
+ #endif
+ } // end if contact
+
+ // done with the edges from body j,
+ // given that vertex ni interacts with only one edge from body j
+
+// break;
+
+ } // end if mode
+
+ } // end for looping through the edges of body j
+
+ } // end for looping through the vertices of body i
+
+ evdwl += energy;
+
+ return interact;
+}
+
+/* -------------------------------------------------------------------------
+ Compute the distance between an edge of body i and a vertex from
+ another body
+ Input:
+ ibody = body i (i.e. atom i)
+ edge_index = edge index of body i
+ xmi = atom i's coordinates (body i's center of mass)
+ x0 = coordinate of the tested vertex from another body
+ x0_rounded_radius = rounded radius of the tested vertex
+ cut_inner = cutoff for vertex-vertex and vertex-edge interaction
+ Output:
+ d = Distance from a point x0 to an edge
+ hi = coordinates of the projection of x0 on the edge
+ t = ratio to determine the relative position of hi
+ wrt xi and xj on the segment
+ contact = 0 no contact between the queried vertex and the edge
+ 1 contact detected
+ return
+ INVALID if the edge index is invalid
+ NONE if there is no interaction
+ VERTEXI if the tested vertex interacts with the first vertex of the edge
+ VERTEXJ if the tested vertex interacts with the second vertex of the edge
+ EDGE if the tested vertex interacts with the edge
+------------------------------------------------------------------------- */
+
+int PairBodyRoundedPolygon::compute_distance_to_vertex(int ibody,
+ int edge_index,
+ double *xmi,
+ double rounded_radius,
+ double* x0,
+ double x0_rounded_radius,
+ double cut_inner,
+ double &d,
+ double hi[3],
+ double &t,
+ int &contact)
+{
+ if (edge_index >= ednum[ibody]) return INVALID;
+
+ int mode,ifirst,iefirst,npi1,npi2;
+ double xi1[3],xi2[3],u[3],v[3],uij[3];
+ double udotv, magv, magucostheta;
+ double delx,dely,delz;
+
+ ifirst = dfirst[ibody];
+ iefirst = edfirst[ibody];
+ npi1 = static_cast<int>(edge[iefirst+edge_index][0]);
+ npi2 = static_cast<int>(edge[iefirst+edge_index][1]);
+
+ // compute the space-fixed coordinates for the vertices of the edge
+
+ xi1[0] = xmi[0] + discrete[ifirst+npi1][0];
+ xi1[1] = xmi[1] + discrete[ifirst+npi1][1];
+ xi1[2] = xmi[2] + discrete[ifirst+npi1][2];
+
+ xi2[0] = xmi[0] + discrete[ifirst+npi2][0];
+ xi2[1] = xmi[1] + discrete[ifirst+npi2][1];
+ xi2[2] = xmi[2] + discrete[ifirst+npi2][2];
+
+ // u = x0 - xi1
+
+ u[0] = x0[0] - xi1[0];
+ u[1] = x0[1] - xi1[1];
+ u[2] = x0[2] - xi1[2];
+
+ // v = xi2 - xi1
+
+ v[0] = xi2[0] - xi1[0];
+ v[1] = xi2[1] - xi1[1];
+ v[2] = xi2[2] - xi1[2];
+
+ // dot product between u and v = magu * magv * costheta
+
+ udotv = u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
+ magv = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+ magucostheta = udotv / magv;
+
+ // uij is the unit vector pointing from xi to xj
+
+ uij[0] = v[0] / magv;
+ uij[1] = v[1] / magv;
+ uij[2] = v[2] / magv;
+
+ // position of the projection of x0 on the line (xi, xj)
+
+ hi[0] = xi1[0] + magucostheta * uij[0];
+ hi[1] = xi1[1] + magucostheta * uij[1];
+ hi[2] = xi1[2] + magucostheta * uij[2];
+
+ // distance from x0 to the line (xi, xj) = distance from x0 to hi
+
+ distance(hi, x0, d);
+
+ // determine the interaction mode
+ // for 2D: a vertex can interact with one edge at most
+ // for 3D: a vertex can interact with one face at most
+
+ mode = NONE;
+ contact = 0;
+
+ if (d > rounded_radius + x0_rounded_radius + cut_inner) {
+
+ // if the vertex is far away from the edge
+
+ mode = NONE;
+
+ } else {
+
+ // check if x0 (the queried vertex) and xmi (the body's center of mass)
+ // are on the different sides of the edge
+
+ int m = opposite_sides(xi1, xi2, x0, xmi);
+
+ if (m == 0) {
+
+ // x0 and xmi are on not the opposite sides of the edge
+ // leave xpi for another edge to detect
+
+ mode = NONE;
+
+ } else {
+
+ // x0 and xmi are on the different sides
+ // t is the ratio to detect if x0 is closer to the vertices xi or xj
+
+ if (fabs(xi2[0] - xi1[0]) > EPSILON)
+ t = (hi[0] - xi1[0]) / (xi2[0] - xi1[0]);
+ else if (fabs(xi2[1] - xi1[1]) > EPSILON)
+ t = (hi[1] - xi1[1]) / (xi2[1] - xi1[1]);
+ else if (fabs(xi2[2] - xi1[2]) > EPSILON)
+ t = (hi[2] - xi1[2]) / (xi2[2] - xi1[2]);
+
+ double contact_dist = rounded_radius + x0_rounded_radius;
+ if (t >= 0 && t <= 1) {
+ mode = EDGE;
+ if (d < contact_dist + EPSILON)
+ contact = 1;
+
+ } else { // t < 0 || t > 1: closer to either vertices of the edge
+
+ if (t < 0) {
+ // measure the distance from x0 to xi1
+ delx = x0[0] - xi1[0];
+ dely = x0[1] - xi1[1];
+ delz = x0[2] - xi1[2];
+ double dx0xi1 = sqrt(delx*delx + dely*dely + delz*delz);
+ if (dx0xi1 > contact_dist + cut_inner)
+ mode = NONE;
+ else
+ mode = VERTEXI;
+ } else {
+ // measure the distance from x0 to xi2
+ delx = x0[0] - xi2[0];
+ dely = x0[1] - xi2[1];
+ delz = x0[2] - xi2[2];
+ double dx0xi2 = sqrt(delx*delx + dely*dely + delz*delz);
+ if (dx0xi2 > contact_dist + cut_inner)
+ mode = NONE;
+ else
+ mode = VERTEXJ;
+ }
+ } // end if t >= 0 && t <= 1
+ } // end if x0 and xmi are on the same side of the edge
+ }
+
+ return mode;
+}
+
+/* ----------------------------------------------------------------------
+ Compute contact forces between two bodies
+ modify the force stored at the vertex and edge in contact by j_a
+ sum forces and torque to the corresponding bodies
+ fn = normal friction component
+ ft = tangential friction component (-c_t * v_t)
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::contact_forces(Contact& contact, double j_a,
+ double** x, double** v, double** angmom, double** f,
+ double** torque, double &evdwl, double* facc)
+{
+ int ibody,jbody,ibonus,jbonus,ifirst,jefirst,ni,nj;
+ double fx,fy,fz,delx,dely,delz,rsq,rsqinv;
+ double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
+ double fn[3],ft[3],vi[3],vj[3];
+ double *quat, *inertia;
+ AtomVecBody::Bonus *bonus;
+
+ ibody = contact.ibody;
+ jbody = contact.jbody;
+
+ // compute the velocity of the vertex in the space-fixed frame
+
+ ibonus = atom->body[ibody];
+ bonus = &avec->bonus[ibonus];
+ quat = bonus->quat;
+ inertia = bonus->inertia;
+ total_velocity(contact.xv, x[ibody], v[ibody], angmom[ibody],
+ inertia, quat, vi);
+
+ // compute the velocity of the point on the edge in the space-fixed frame
+
+ jbonus = atom->body[jbody];
+ bonus = &avec->bonus[jbonus];
+ quat = bonus->quat;
+ inertia = bonus->inertia;
+ total_velocity(contact.xe, x[jbody], v[jbody], angmom[jbody],
+ inertia, quat, vj);
+
+ // vector pointing from the vertex to the point on the edge
+
+ delx = contact.xv[0] - contact.xe[0];
+ dely = contact.xv[1] - contact.xe[1];
+ delz = contact.xv[2] - contact.xe[2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ rsqinv = 1.0/rsq;
+
+ // relative translational velocity
+
+ vr1 = vi[0] - vj[0];
+ vr2 = vi[1] - vj[1];
+ vr3 = vi[2] - vj[2];
+
+ // normal component
+
+ vnnr = vr1*delx + vr2*dely + vr3*delz;
+ vn1 = delx*vnnr * rsqinv;
+ vn2 = dely*vnnr * rsqinv;
+ vn3 = delz*vnnr * rsqinv;
+
+ // tangential component
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // normal friction term at contact
+
+ fn[0] = -c_n * vn1;
+ fn[1] = -c_n * vn2;
+ fn[2] = -c_n * vn3;
+
+ // tangential friction term at contact
+ // excluding the tangential deformation term for now
+
+ ft[0] = -c_t * vt1;
+ ft[1] = -c_t * vt2;
+ ft[2] = -c_t * vt3;
+
+ // only the cohesive force is scaled by j_a
+ // mu * fne = tangential friction deformation during gross sliding
+ // see Eq. 4, Fraige et al.
+
+ ifirst = dfirst[ibody];
+ ni = contact.vertex;
+
+ fx = discrete[ifirst+ni][3] * j_a + fn[0] + ft[0] +
+ mu * discrete[ifirst+ni][3];
+ fy = discrete[ifirst+ni][4] * j_a + fn[1] + ft[1] +
+ mu * discrete[ifirst+ni][4];
+ fz = discrete[ifirst+ni][5] * j_a + fn[2] + ft[2] +
+ mu * discrete[ifirst+ni][5];
+ f[ibody][0] += fx;
+ f[ibody][1] += fy;
+ f[ibody][2] += fz;
+ sum_torque(x[ibody], contact.xv, fx, fy, fz, torque[ibody]);
+
+ // accumulate forces to the vertex only
+
+ facc[0] += fx; facc[1] += fy; facc[2] += fz;
+
+ // only the cohesive force is scaled by j_a
+ // mu * fne = tangential friction deformation during gross sliding
+ // Eq. 4, Fraige et al.
+
+ jefirst = edfirst[jbody];
+ nj = contact.edge;
+
+ fx = edge[jefirst+nj][2] * j_a - fn[0] - ft[0] +
+ mu * edge[jefirst+nj][2];
+ fy = edge[jefirst+nj][3] * j_a - fn[1] - ft[1] +
+ mu * edge[jefirst+nj][3];
+ fz = edge[jefirst+nj][4] * j_a - fn[2] - ft[2] +
+ mu * edge[jefirst+nj][4];
+ f[jbody][0] += fx;
+ f[jbody][1] += fy;
+ f[jbody][2] += fz;
+ sum_torque(x[jbody], contact.xe, fx, fy, fz, torque[jbody]);
+
+ #ifdef _POLYGON_DEBUG
+ printf("From contact forces: vertex fx %f fy %f fz %f\n"
+ " torque body %d: %f %f %f\n"
+ " torque body %d: %f %f %f\n",
+ discrete[ifirst+ni][3], discrete[ifirst+ni][4], discrete[ifirst+ni][5],
+ atom->tag[ibody],torque[ibody][0],torque[ibody][1],torque[ibody][2],
+ atom->tag[jbody],torque[jbody][0],torque[jbody][1],torque[jbody][2]);
+ #endif
+}
+
+/* ----------------------------------------------------------------------
+ Determine the length of the contact segment, i.e. the separation between
+ 2 contacts, should be extended for 3D models.
+------------------------------------------------------------------------- */
+
+double PairBodyRoundedPolygon::contact_separation(const Contact& c1,
+ const Contact& c2)
+{
+ double x1 = c1.xv[0];
+ double y1 = c1.xv[1];
+ double x2 = c1.xe[0];
+ double y2 = c1.xe[1];
+ double x3 = c2.xv[0];
+ double y3 = c2.xv[1];
+
+ double delta_a = 0.0;
+ if (fabs(x2 - x1) > EPSILON) {
+ double A = (y2 - y1) / (x2 - x1);
+ delta_a = fabs(y1 - A * x1 - y3 + A * x3) / sqrt(1 + A * A);
+ } else {
+ delta_a = fabs(x1 - x3);
+ }
+
+ return delta_a;
+}
+
+/* ----------------------------------------------------------------------
+ Accumulate torque to body from the force f=(fx,fy,fz) acting at point x
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::sum_torque(double* xm, double *x, double fx,
+ double fy, double fz, double* torque)
+{
+ double rx = x[0] - xm[0];
+ double ry = x[1] - xm[1];
+ double rz = x[2] - xm[2];
+ double tx = ry * fz - rz * fy;
+ double ty = rz * fx - rx * fz;
+ double tz = rx * fy - ry * fx;
+ torque[0] += tx;
+ torque[1] += ty;
+ torque[2] += tz;
+}
+
+/* ----------------------------------------------------------------------
+ Test if two points a and b are in opposite sides of the line that
+ connects two points x1 and x2
+------------------------------------------------------------------------- */
+
+int PairBodyRoundedPolygon::opposite_sides(double* x1, double* x2,
+ double* a, double* b)
+{
+ double m_a = (x1[1] - x2[1])*(a[0] - x1[0]) + (x2[0] - x1[0])*(a[1] - x1[1]);
+ double m_b = (x1[1] - x2[1])*(b[0] - x1[0]) + (x2[0] - x1[0])*(b[1] - x1[1]);
+ // equal to zero when either a or b is inline with the line x1-x2
+ if (m_a * m_b <= 0)
+ return 1;
+ else
+ return 0;
+}
+
+/* ----------------------------------------------------------------------
+ Calculate the total velocity of a point (vertex, a point on an edge):
+ vi = vcm + omega ^ (p - xcm)
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::total_velocity(double* p, double *xcm,
+ double* vcm, double *angmom, double *inertia,
+ double *quat, double* vi)
+{
+ double r[3],omega[3],ex_space[3],ey_space[3],ez_space[3];
+ r[0] = p[0] - xcm[0];
+ r[1] = p[1] - xcm[1];
+ r[2] = p[2] - xcm[2];
+ MathExtra::q_to_exyz(quat,ex_space,ey_space,ez_space);
+ MathExtra::angmom_to_omega(angmom,ex_space,ey_space,ez_space,
+ inertia,omega);
+ vi[0] = omega[1]*r[2] - omega[2]*r[1] + vcm[0];
+ vi[1] = omega[2]*r[0] - omega[0]*r[2] + vcm[1];
+ vi[2] = omega[0]*r[1] - omega[1]*r[0] + vcm[2];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBodyRoundedPolygon::distance(const double* x2, const double* x1,
+ double& r)
+{
+ r = sqrt((x2[0] - x1[0]) * (x2[0] - x1[0])
+ + (x2[1] - x1[1]) * (x2[1] - x1[1])
+ + (x2[2] - x1[2]) * (x2[2] - x1[2]));
+}
+
diff --git a/src/BODY/pair_body_rounded_polygon.h b/src/BODY/pair_body_rounded_polygon.h
new file mode 100644
index 0000000000000000000000000000000000000000..09c93b832e0d322594b55a7f78ddbb309e158cf2
--- /dev/null
+++ b/src/BODY/pair_body_rounded_polygon.h
@@ -0,0 +1,128 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(body/rounded/polygon,PairBodyRoundedPolygon)
+
+#else
+
+#ifndef LMP_PAIR_BODY_ROUNDED_POLYGON_H
+#define LMP_PAIR_BODY_ROUNDED_POLYGON_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairBodyRoundedPolygon : public Pair {
+ public:
+ PairBodyRoundedPolygon(class LAMMPS *);
+ ~PairBodyRoundedPolygon();
+ void compute(int, int);
+ void settings(int, char **);
+ void coeff(int, char **);
+ void init_style();
+ double init_one(int, int);
+
+ struct Contact {
+ int ibody, jbody; // body (i.e. atom) indices (not tags)
+ int vertex; // vertex of the first polygon
+ int edge; // edge of the second polygon
+ double xv[3]; // coordinates of the vertex
+ double xe[3]; // coordinates of the projection of the vertex on the edge
+ double separation;// separation at contact
+ };
+
+ protected:
+ double **k_n; // normal repulsion strength
+ double **k_na; // normal attraction strength
+ double c_n; // normal damping coefficient
+ double c_t; // tangential damping coefficient
+ double mu; // normal friction coefficient during gross sliding
+ double delta_ua; // contact line (area for 3D models) modification factor
+ double cut_inner; // cutoff for interaction between vertex-edge surfaces
+
+ class AtomVecBody *avec;
+ class BodyRoundedPolygon *bptr;
+
+ double **discrete; // list of all sub-particles for all bodies
+ int ndiscrete; // number of discretes in list
+ int dmax; // allocated size of discrete list
+ int *dnum; // number of discretes per line, 0 if uninit
+ int *dfirst; // index of first discrete per each line
+ int nmax; // allocated size of dnum,dfirst vectors
+
+ double **edge; // list of all edge for all bodies
+ int nedge; // number of edge in list
+ int edmax; // allocated size of edge list
+ int *ednum; // number of edges per line, 0 if uninit
+ int *edfirst; // index of first edge per each line
+ int ednummax; // allocated size of ednum,edfirst vectors
+
+ double *enclosing_radius; // enclosing radii for all bodies
+ double *rounded_radius; // rounded radii for all bodies
+ double *maxerad; // per-type maximum enclosing radius
+
+ void allocate();
+ void body2space(int);
+
+ int vertex_against_edge(int i, int j, double k_n, double k_na,
+ double** x, double** f, double** torque,
+ tagint* tag, Contact* contact_list,
+ int &num_contacts, double &evdwl, double* facc);
+ void sphere_against_sphere(int i, int j, double delx, double dely, double delz,
+ double rsq, double k_n, double k_na,
+ double** x, double** v, double** f, int evflag);
+ int compute_distance_to_vertex(int ibody, int edge_index, double* xmi,
+ double rounded_radius, double* x0,
+ double x0_rounded_radius, double cut_inner,
+ double &d, double hi[3], double &t,
+ int &contact);
+ double contact_separation(const Contact& c1, const Contact& c2);
+ void contact_forces(Contact& contact, double j_a, double** x,
+ double** v, double** f, double** angmom,
+ double** torque, double &evdwl, double* facc);
+ void sum_torque(double* xm, double *x, double fx,
+ double fy, double fz, double* torque);
+ int opposite_sides(double* x1, double* x2, double* a, double* b);
+ void total_velocity(double* p, double *xcm, double* vcm, double *angmom,
+ double *inertia, double *quat, double* vi);
+ inline void distance(const double* x2, const double* x1, double& r);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Pair body/rounded/polygon requires atom style body rounded/polygon
+
+Self-explanatory.
+
+E: Pair body requires body style rounded/polygon
+
+This pair style is specific to the rounded/polygon body style.
+
+*/
diff --git a/src/BODY/pair_body_rounded_polyhedron.cpp b/src/BODY/pair_body_rounded_polyhedron.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8d5f9ec72cb2d439c433b1f57c01d1536c6e4527
--- /dev/null
+++ b/src/BODY/pair_body_rounded_polyhedron.cpp
@@ -0,0 +1,2348 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)
+ Ref: Wang, Yu, Langston, Fraige, Particle shape effects in discrete
+ element modelling of cohesive angular particles, Granular Matter 2011,
+ 13:1-12.
+ Note: The current implementation has not taken into account
+ the contact history for friction forces.
+------------------------------------------------------------------------- */
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include "pair_body_rounded_polyhedron.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "atom_vec_body.h"
+#include "body_rounded_polyhedron.h"
+#include "comm.h"
+#include "force.h"
+#include "fix.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "memory.h"
+#include "error.h"
+#include "math_extra.h"
+#include "math_const.h"
+
+using namespace LAMMPS_NS;
+using namespace MathExtra;
+using namespace MathConst;
+
+#define DELTA 10000
+#define EPSILON 1e-3
+#define MAX_FACE_SIZE 4 // maximum number of vertices per face (same as BodyRoundedPolyhedron)
+#define MAX_CONTACTS 16 // for 3D models
+
+//#define _POLYHEDRON_DEBUG
+
+enum {EE_INVALID=0,EE_NONE,EE_INTERACT};
+enum {EF_INVALID=0,EF_NONE,EF_PARALLEL,EF_SAME_SIDE_OF_FACE,
+ EF_INTERSECT_INSIDE,EF_INTERSECT_OUTSIDE};
+
+/* ---------------------------------------------------------------------- */
+
+PairBodyRoundedPolyhedron::PairBodyRoundedPolyhedron(LAMMPS *lmp) : Pair(lmp)
+{
+ dmax = nmax = 0;
+ discrete = NULL;
+ dnum = dfirst = NULL;
+
+ edmax = ednummax = 0;
+ edge = NULL;
+ ednum = edfirst = NULL;
+
+ facmax = facnummax = 0;
+ face = NULL;
+ facnum = facfirst = NULL;
+
+ enclosing_radius = NULL;
+ rounded_radius = NULL;
+ maxerad = NULL;
+
+ single_enable = 0;
+ restartinfo = 0;
+
+ c_n = 0.1;
+ c_t = 0.2;
+ mu = 0.0;
+ A_ua = 1.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairBodyRoundedPolyhedron::~PairBodyRoundedPolyhedron()
+{
+ memory->destroy(discrete);
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+
+ memory->destroy(edge);
+ memory->destroy(ednum);
+ memory->destroy(edfirst);
+
+ memory->destroy(face);
+ memory->destroy(facnum);
+ memory->destroy(facfirst);
+
+ memory->destroy(enclosing_radius);
+ memory->destroy(rounded_radius);
+ memory->destroy(maxerad);
+
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+
+ memory->destroy(k_n);
+ memory->destroy(k_na);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ int ni,nj,npi,npj,ifirst,jfirst,nei,nej,iefirst,jefirst;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,facc[3];
+ double rsq,eradi,eradj,k_nij,k_naij;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **v = atom->v;
+ double **f = atom->f;
+ double **torque = atom->torque;
+ double **angmom = atom->angmom;
+ int *body = atom->body;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ int nall = nlocal + atom->nghost;
+ int newton_pair = force->newton_pair;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // grow the per-atom lists if necessary and initialize
+
+ if (atom->nmax > nmax) {
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+ memory->destroy(ednum);
+ memory->destroy(edfirst);
+ memory->destroy(facnum);
+ memory->destroy(facfirst);
+ memory->destroy(enclosing_radius);
+ memory->destroy(rounded_radius);
+ nmax = atom->nmax;
+ memory->create(dnum,nmax,"pair:dnum");
+ memory->create(dfirst,nmax,"pair:dfirst");
+ memory->create(ednum,nmax,"pair:ednum");
+ memory->create(edfirst,nmax,"pair:edfirst");
+ memory->create(facnum,nmax,"pair:facnum");
+ memory->create(facfirst,nmax,"pair:facfirst");
+ memory->create(enclosing_radius,nmax,"pair:enclosing_radius");
+ memory->create(rounded_radius,nmax,"pair:rounded_radius");
+ }
+
+ ndiscrete = nedge = nface = 0;
+ for (i = 0; i < nall; i++)
+ dnum[i] = ednum[i] = facnum[i] = 0;
+
+ // loop over neighbors of my atoms
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = type[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ if (body[i] >= 0) {
+ if (dnum[i] == 0) body2space(i);
+ npi = dnum[i];
+ ifirst = dfirst[i];
+ nei = ednum[i];
+ iefirst = edfirst[i];
+ eradi = enclosing_radius[i];
+ }
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+
+ // body/body interactions
+
+ evdwl = 0.0;
+ facc[0] = facc[1] = facc[2] = 0;
+
+ if (body[i] < 0 || body[j] < 0) continue;
+
+ if (dnum[j] == 0) body2space(j);
+ npj = dnum[j];
+ jfirst = dfirst[j];
+ nej = ednum[j];
+ jefirst = edfirst[j];
+ eradj = enclosing_radius[j];
+
+ k_nij = k_n[itype][jtype];
+ k_naij = k_na[itype][jtype];
+
+ // no interaction
+
+ double r = sqrt(rsq);
+ if (r > eradi + eradj + cut_inner) continue;
+
+ // sphere-sphere interaction
+
+ if (npi == 1 && npj == 1) {
+ sphere_against_sphere(i, j, delx, dely, delz, rsq,
+ k_nij, k_naij, v, f, evflag);
+ continue;
+ }
+
+ // reset vertex and edge forces
+
+ for (ni = 0; ni < npi; ni++) {
+ discrete[ifirst+ni][3] = 0;
+ discrete[ifirst+ni][4] = 0;
+ discrete[ifirst+ni][5] = 0;
+ discrete[ifirst+ni][6] = 0;
+ }
+
+ for (nj = 0; nj < npj; nj++) {
+ discrete[jfirst+nj][3] = 0;
+ discrete[jfirst+nj][4] = 0;
+ discrete[jfirst+nj][5] = 0;
+ discrete[jfirst+nj][6] = 0;
+ }
+
+ for (ni = 0; ni < nei; ni++) {
+ edge[iefirst+ni][2] = 0;
+ edge[iefirst+ni][3] = 0;
+ edge[iefirst+ni][4] = 0;
+ edge[iefirst+ni][5] = 0;
+ }
+
+ for (nj = 0; nj < nej; nj++) {
+ edge[jefirst+nj][2] = 0;
+ edge[jefirst+nj][3] = 0;
+ edge[jefirst+nj][4] = 0;
+ edge[jefirst+nj][5] = 0;
+ }
+
+ // one of the two bodies is a sphere
+
+ if (npj == 1) {
+ sphere_against_face(i, j, k_nij, k_naij, x, v, f, torque,
+ angmom, evflag);
+ sphere_against_edge(i, j, k_nij, k_naij, x, v, f, torque,
+ angmom, evflag);
+ continue;
+ } else if (npi == 1) {
+ sphere_against_face(j, i, k_nij, k_naij, x, v, f, torque,
+ angmom, evflag);
+ sphere_against_edge(j, i, k_nij, k_naij, x, v, f, torque,
+ angmom, evflag);
+ continue;
+ }
+
+ int interact, num_contacts;
+ Contact contact_list[MAX_CONTACTS];
+
+ num_contacts = 0;
+
+ // check interaction between i's edges and j' faces
+ #ifdef _POLYHEDRON_DEBUG
+ printf("INTERACTION between edges of %d vs. faces of %d:\n", i, j);
+ #endif
+ interact = edge_against_face(i, j, k_nij, k_naij, x, contact_list,
+ num_contacts, evdwl, facc);
+
+ // check interaction between j's edges and i' faces
+ #ifdef _POLYHEDRON_DEBUG
+ printf("\nINTERACTION between edges of %d vs. faces of %d:\n", j, i);
+ #endif
+ interact = edge_against_face(j, i, k_nij, k_naij, x, contact_list,
+ num_contacts, evdwl, facc);
+
+ // check interaction between i's edges and j' edges
+ #ifdef _POLYHEDRON_DEBUG
+ printf("INTERACTION between edges of %d vs. edges of %d:\n", i, j);
+ #endif
+ interact = edge_against_edge(i, j, k_nij, k_naij, x, contact_list,
+ num_contacts, evdwl, facc);
+
+ // estimate the contact area
+ // also consider point contacts and line contacts
+
+ if (num_contacts > 0) {
+ double contact_area;
+ if (num_contacts == 1) {
+ contact_area = 0;
+ } else if (num_contacts == 2) {
+ contact_area = num_contacts * A_ua;
+ } else {
+ int m;
+ double xc[3],dx,dy,dz;
+ xc[0] = xc[1] = xc[2] = 0;
+ for (m = 0; m < num_contacts; m++) {
+ xc[0] += contact_list[m].xi[0];
+ xc[1] += contact_list[m].xi[1];
+ xc[2] += contact_list[m].xi[2];
+ }
+
+ xc[0] /= (double)num_contacts;
+ xc[1] /= (double)num_contacts;
+ xc[2] /= (double)num_contacts;
+
+ contact_area = 0.0;
+ for (m = 0; m < num_contacts; m++) {
+ dx = contact_list[m].xi[0] - xc[0];
+ dy = contact_list[m].xi[1] - xc[1];
+ dz = contact_list[m].xi[2] - xc[2];
+ contact_area += (dx*dx + dy*dy + dz*dz);
+ }
+ contact_area *= (MY_PI/(double)num_contacts);
+ }
+ rescale_cohesive_forces(x, f, torque, contact_list, num_contacts,
+ contact_area, k_nij, k_naij, facc);
+ }
+
+ if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0.0,
+ facc[0],facc[1],facc[2],delx,dely,delz);
+
+ } // end for jj
+ }
+
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ memory->create(setflag,n+1,n+1,"pair:setflag");
+ for (int i = 1; i <= n; i++)
+ for (int j = i; j <= n; j++)
+ setflag[i][j] = 0;
+
+ memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+ memory->create(k_n,n+1,n+1,"pair:k_n");
+ memory->create(k_na,n+1,n+1,"pair:k_na");
+ memory->create(maxerad,n+1,"pair:maxerad");
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::settings(int narg, char **arg)
+{
+ if (narg < 5) error->all(FLERR,"Illegal pair_style command");
+
+ c_n = force->numeric(FLERR,arg[0]);
+ c_t = force->numeric(FLERR,arg[1]);
+ mu = force->numeric(FLERR,arg[2]);
+ A_ua = force->numeric(FLERR,arg[3]);
+ cut_inner = force->numeric(FLERR,arg[4]);
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::coeff(int narg, char **arg)
+{
+ if (narg < 4 || narg > 5)
+ error->all(FLERR,"Incorrect args for pair coefficients");
+ if (!allocated) allocate();
+
+ int ilo,ihi,jlo,jhi;
+ force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
+ force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
+
+ double k_n_one = force->numeric(FLERR,arg[2]);
+ double k_na_one = force->numeric(FLERR,arg[3]);
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ k_n[i][j] = k_n_one;
+ k_na[i][j] = k_na_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::init_style()
+{
+ avec = (AtomVecBody *) atom->style_match("body");
+ if (!avec) error->all(FLERR,"Pair body/rounded/polyhedron requires atom style body");
+ if (strcmp(avec->bptr->style,"rounded/polyhedron") != 0)
+ error->all(FLERR,"Pair body/rounded/polyhedron requires body style rounded/polyhedron");
+ bptr = (BodyRoundedPolyhedron *) avec->bptr;
+
+ if (comm->ghost_velocity == 0)
+ error->all(FLERR,"Pair body/rounded/polyhedron requires ghost atoms store velocity");
+
+ neighbor->request(this);
+
+ // find the maximum enclosing radius for each atom type
+
+ int i, itype;
+ double eradi;
+ int* body = atom->body;
+ int* type = atom->type;
+ int ntypes = atom->ntypes;
+ int nlocal = atom->nlocal;
+
+ if (atom->nmax > nmax) {
+ memory->destroy(dnum);
+ memory->destroy(dfirst);
+ memory->destroy(ednum);
+ memory->destroy(edfirst);
+ memory->destroy(facnum);
+ memory->destroy(facfirst);
+ memory->destroy(enclosing_radius);
+ memory->destroy(rounded_radius);
+ nmax = atom->nmax;
+ memory->create(dnum,nmax,"pair:dnum");
+ memory->create(dfirst,nmax,"pair:dfirst");
+ memory->create(ednum,nmax,"pair:ednum");
+ memory->create(edfirst,nmax,"pair:edfirst");
+ memory->create(facnum,nmax,"pair:facnum");
+ memory->create(facfirst,nmax,"pair:facfirst");
+ memory->create(enclosing_radius,nmax,"pair:enclosing_radius");
+ memory->create(rounded_radius,nmax,"pair:rounded_radius");
+ }
+
+ ndiscrete = nedge = nface = 0;
+ for (i = 0; i < nlocal; i++)
+ dnum[i] = ednum[i] = facnum[i] = 0;
+
+ double *merad = NULL;
+ memory->create(merad,ntypes+1,"pair:merad");
+ for (i = 1; i <= ntypes; i++)
+ maxerad[i] = merad[i] = 0;
+
+ int ipour;
+ for (ipour = 0; ipour < modify->nfix; ipour++)
+ if (strcmp(modify->fix[ipour]->style,"pour") == 0) break;
+ if (ipour == modify->nfix) ipour = -1;
+
+ int idep;
+ for (idep = 0; idep < modify->nfix; idep++)
+ if (strcmp(modify->fix[idep]->style,"deposit") == 0) break;
+ if (idep == modify->nfix) idep = -1;
+
+ for (i = 1; i <= ntypes; i++) {
+ merad[i] = 0.0;
+ if (ipour >= 0) {
+ itype = i;
+ merad[i] =
+ *((double *) modify->fix[ipour]->extract("radius",itype));
+ }
+ if (idep >= 0) {
+ itype = i;
+ merad[i] =
+ *((double *) modify->fix[idep]->extract("radius",itype));
+ }
+ }
+
+ for (i = 0; i < nlocal; i++) {
+ itype = type[i];
+ if (body[i] >= 0) {
+ if (dnum[i] == 0) body2space(i);
+ eradi = enclosing_radius[i];
+ if (eradi > merad[itype]) merad[itype] = eradi;
+ } else
+ merad[itype] = 0;
+ }
+
+ MPI_Allreduce(&merad[1],&maxerad[1],ntypes,MPI_DOUBLE,MPI_MAX,world);
+
+ memory->destroy(merad);
+
+ sanity_check();
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairBodyRoundedPolyhedron::init_one(int i, int j)
+{
+ k_n[j][i] = k_n[i][j];
+ k_na[j][i] = k_na[i][j];
+
+ return (maxerad[i]+maxerad[j]);
+}
+
+/* ----------------------------------------------------------------------
+ convert N sub-particles in body I to space frame using current quaternion
+ store sub-particle space-frame displacements from COM in discrete list
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::body2space(int i)
+{
+ int ibonus = atom->body[i];
+ AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
+ int nsub = bptr->nsub(bonus);
+ double *coords = bptr->coords(bonus);
+ int body_num_edges = bptr->nedges(bonus);
+ double* edge_ends = bptr->edges(bonus);
+ int body_num_faces = bptr->nfaces(bonus);
+ double* face_pts = bptr->faces(bonus);
+ double eradius = bptr->enclosing_radius(bonus);
+ double rradius = bptr->rounded_radius(bonus);
+
+ // get the number of sub-particles (vertices)
+ // and the index of the first vertex of my body in the list
+
+ dnum[i] = nsub;
+ dfirst[i] = ndiscrete;
+
+ // grow the vertex list if necessary
+ // the first 3 columns are for coords, the last 3 for forces
+
+ if (ndiscrete + nsub > dmax) {
+ dmax += DELTA;
+ memory->grow(discrete,dmax,7,"pair:discrete");
+ }
+
+ double p[3][3];
+ MathExtra::quat_to_mat(bonus->quat,p);
+
+ for (int m = 0; m < nsub; m++) {
+ MathExtra::matvec(p,&coords[3*m],discrete[ndiscrete]);
+ discrete[ndiscrete][3] = 0;
+ discrete[ndiscrete][4] = 0;
+ discrete[ndiscrete][5] = 0;
+ discrete[ndiscrete][6] = 0;
+ ndiscrete++;
+ }
+
+ // get the number of edges (vertices)
+ // and the index of the first edge of my body in the list
+
+ ednum[i] = body_num_edges;
+ edfirst[i] = nedge;
+
+ // grow the edge list if necessary
+ // the first 2 columns are for vertex indices within body, the last 3 for forces
+
+ if (nedge + body_num_edges > edmax) {
+ edmax += DELTA;
+ memory->grow(edge,edmax,6,"pair:edge");
+ }
+
+ for (int m = 0; m < body_num_edges; m++) {
+ edge[nedge][0] = static_cast<int>(edge_ends[2*m+0]);
+ edge[nedge][1] = static_cast<int>(edge_ends[2*m+1]);
+ edge[nedge][2] = 0;
+ edge[nedge][3] = 0;
+ edge[nedge][4] = 0;
+ edge[nedge][5] = 0;
+ nedge++;
+ }
+
+ // get the number of faces and the index of the first face
+
+ facnum[i] = body_num_faces;
+ facfirst[i] = nface;
+
+ // grow the face list if necessary
+ // the first 3 columns are for vertex indices within body, the last 3 for forces
+
+ if (nface + body_num_faces > facmax) {
+ facmax += DELTA;
+ memory->grow(face,facmax,MAX_FACE_SIZE,"pair:face");
+ }
+
+ for (int m = 0; m < body_num_faces; m++) {
+ for (int k = 0; k < MAX_FACE_SIZE; k++)
+ face[nface][k] = static_cast<int>(face_pts[MAX_FACE_SIZE*m+k]);
+ nface++;
+ }
+
+ enclosing_radius[i] = eradius;
+ rounded_radius[i] = rradius;
+}
+
+/* ----------------------------------------------------------------------
+ Interaction between two spheres with different radii
+ according to the 2D model from Fraige et al.
+---------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::sphere_against_sphere(int ibody, int jbody,
+ double delx, double dely, double delz, double rsq,
+ double k_n, double k_na, double** v, double** f,
+ int evflag)
+{
+ double rradi,rradj,contact_dist;
+ double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
+ double rij,rsqinv,R,fx,fy,fz,fn[3],ft[3],fpair,shift,energy;
+ int nlocal = atom->nlocal;
+ int newton_pair = force->newton_pair;
+
+ rradi = rounded_radius[ibody];
+ rradj = rounded_radius[jbody];
+ contact_dist = rradi + rradj;
+
+ rij = sqrt(rsq);
+ R = rij - contact_dist;
+ shift = k_na * cut_inner;
+
+ energy = 0;
+
+ if (R <= 0) { // deformation occurs
+ fpair = -k_n * R - shift;
+ energy = (0.5 * k_n * R + shift) * R;
+ } else if (R <= cut_inner) { // not deforming but cohesive ranges overlap
+ fpair = k_na * R - shift;
+ energy = (-0.5 * k_na * R + shift) * R;
+ } else fpair = 0.0;
+
+ fx = delx*fpair/rij;
+ fy = dely*fpair/rij;
+ fz = delz*fpair/rij;
+
+ if (R <= 0) { // in contact
+
+ // relative translational velocity
+
+ vr1 = v[ibody][0] - v[jbody][0];
+ vr2 = v[ibody][1] - v[jbody][1];
+ vr3 = v[ibody][2] - v[jbody][2];
+
+ // normal component
+
+ rsqinv = 1.0/rsq;
+ vnnr = vr1*delx + vr2*dely + vr3*delz;
+ vn1 = delx*vnnr * rsqinv;
+ vn2 = dely*vnnr * rsqinv;
+ vn3 = delz*vnnr * rsqinv;
+
+ // tangential component
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // normal friction term at contact
+
+ fn[0] = -c_n * vn1;
+ fn[1] = -c_n * vn2;
+ fn[2] = -c_n * vn3;
+
+ // tangential friction term at contact,
+ // excluding the tangential deformation term for now
+
+ ft[0] = -c_t * vt1;
+ ft[1] = -c_t * vt2;
+ ft[2] = -c_t * vt3;
+
+ fx += fn[0] + ft[0];
+ fy += fn[1] + ft[1];
+ fz += fn[2] + ft[2];
+ }
+
+ f[ibody][0] += fx;
+ f[ibody][1] += fy;
+ f[ibody][2] += fz;
+
+ if (newton_pair || jbody < nlocal) {
+ f[jbody][0] -= fx;
+ f[jbody][1] -= fy;
+ f[jbody][2] -= fz;
+ }
+
+ if (evflag) ev_tally_xyz(ibody,jbody,nlocal,newton_pair,
+ energy,0.0,fx,fy,fz,delx,dely,delz);
+}
+
+/* ----------------------------------------------------------------------
+ Interaction bt the faces of a polyhedron (ibody) and a sphere (jbody)
+---------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::sphere_against_face(int ibody, int jbody,
+ double k_n, double k_na, double** x, double** v,
+ double** f, double** torque, double** angmom,
+ int evflag)
+{
+ int ni,nfi,inside,ifirst,iffirst,npi1,npi2,npi3,ibonus,tmp;
+ double xi1[3],xi2[3],xi3[3],ui[3],vi[3],vti[3],n[3],h[3],fn[3],ft[3],d;
+ double delx,dely,delz,rsq,rij,rsqinv,R,fx,fy,fz,fpair,shift,energy;
+ double rradi,rradj,contact_dist;
+ double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
+ double *quat, *inertia;
+ AtomVecBody::Bonus *bonus;
+
+ int nlocal = atom->nlocal;
+ int newton_pair = force->newton_pair;
+
+ ifirst = dfirst[ibody];
+ iffirst = facfirst[ibody];
+ nfi = facnum[ibody];
+
+ rradi = rounded_radius[ibody];
+ rradj = rounded_radius[jbody];
+ contact_dist = rradi + rradj;
+
+ for (ni = 0; ni < nfi; ni++) {
+
+ npi1 = static_cast<int>(face[iffirst+ni][0]);
+ npi2 = static_cast<int>(face[iffirst+ni][1]);
+ npi3 = static_cast<int>(face[iffirst+ni][2]);
+
+ // compute the space-fixed coordinates for the vertices of the face
+
+ xi1[0] = x[ibody][0] + discrete[ifirst+npi1][0];
+ xi1[1] = x[ibody][1] + discrete[ifirst+npi1][1];
+ xi1[2] = x[ibody][2] + discrete[ifirst+npi1][2];
+
+ xi2[0] = x[ibody][0] + discrete[ifirst+npi2][0];
+ xi2[1] = x[ibody][1] + discrete[ifirst+npi2][1];
+ xi2[2] = x[ibody][2] + discrete[ifirst+npi2][2];
+
+ xi3[0] = x[ibody][0] + discrete[ifirst+npi3][0];
+ xi3[1] = x[ibody][1] + discrete[ifirst+npi3][1];
+ xi3[2] = x[ibody][2] + discrete[ifirst+npi3][2];
+
+ // find the normal unit vector of the face
+
+ MathExtra::sub3(xi2, xi1, ui);
+ MathExtra::sub3(xi3, xi1, vi);
+ MathExtra::cross3(ui, vi, n);
+ MathExtra::norm3(n);
+
+ // skip if the COM of the two bodies are in the same side of the face
+
+ if (opposite_sides(n, xi1, x[ibody], x[jbody]) == 0) continue;
+
+ // find the projection of the sphere on the face
+
+ project_pt_plane(x[jbody], xi1, xi2, xi3, h, d, inside);
+
+ inside_polygon(ibody, ni, x[ibody], h, NULL, inside, tmp);
+ if (inside == 0) continue;
+
+ delx = h[0] - x[jbody][0];
+ dely = h[1] - x[jbody][1];
+ delz = h[2] - x[jbody][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ rij = sqrt(rsq);
+ R = rij - contact_dist;
+ shift = k_na * cut_inner;
+
+ energy = 0;
+
+ if (R <= 0) { // deformation occurs
+ fpair = -k_n * R - shift;
+ energy = (0.5 * k_n * R + shift) * R;
+ } else if (R <= cut_inner) { // not deforming but cohesive ranges overlap
+ fpair = k_na * R - shift;
+ energy = (-0.5 * k_na * R + shift) * R;
+ } else fpair = 0.0;
+
+ fx = delx*fpair/rij;
+ fy = dely*fpair/rij;
+ fz = delz*fpair/rij;
+
+ if (R <= 0) { // in contact
+
+ // compute the velocity of the vertex in the space-fixed frame
+
+ ibonus = atom->body[ibody];
+ bonus = &avec->bonus[ibonus];
+ quat = bonus->quat;
+ inertia = bonus->inertia;
+ total_velocity(h, x[ibody], v[ibody], angmom[ibody],
+ inertia, quat, vti);
+
+ // relative translational velocity
+
+ vr1 = vti[0] - v[jbody][0];
+ vr2 = vti[1] - v[jbody][1];
+ vr3 = vti[2] - v[jbody][2];
+
+ // normal component
+
+ rsqinv = 1.0/rsq;
+ vnnr = vr1*delx + vr2*dely + vr3*delz;
+ vn1 = delx*vnnr * rsqinv;
+ vn2 = dely*vnnr * rsqinv;
+ vn3 = delz*vnnr * rsqinv;
+
+ // tangential component
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // normal friction term at contact
+
+ fn[0] = -c_n * vn1;
+ fn[1] = -c_n * vn2;
+ fn[2] = -c_n * vn3;
+
+ // tangential friction term at contact,
+ // excluding the tangential deformation term for now
+
+ ft[0] = -c_t * vt1;
+ ft[1] = -c_t * vt2;
+ ft[2] = -c_t * vt3;
+
+ fx += fn[0] + ft[0];
+ fy += fn[1] + ft[1];
+ fz += fn[2] + ft[2];
+ }
+
+ f[ibody][0] += fx;
+ f[ibody][1] += fy;
+ f[ibody][2] += fz;
+ sum_torque(x[ibody], h, fx, fy, fz, torque[ibody]);
+
+ if (newton_pair || jbody < nlocal) {
+ f[jbody][0] -= fx;
+ f[jbody][1] -= fy;
+ f[jbody][2] -= fz;
+ }
+
+ if (evflag) ev_tally_xyz(ibody,jbody,nlocal,newton_pair,
+ energy,0.0,fx,fy,fz,delx,dely,delz);
+ }
+}
+
+/* ----------------------------------------------------------------------
+ Interaction bt the edges of a polyhedron (ibody) and a sphere (jbody)
+---------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::sphere_against_edge(int ibody, int jbody,
+ double k_n, double k_na, double** x, double** v,
+ double** f, double** torque, double** angmom,
+ int evflag)
+{
+ int ni,nei,ifirst,iefirst,npi1,npi2,ibonus;
+ double xi1[3],xi2[3],vti[3],h[3],fn[3],ft[3],d,t;
+ double delx,dely,delz,rij,rsqinv,R,fx,fy,fz,fpair,shift,energy;
+ double rradi,rradj,contact_dist;
+ double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
+ double *quat, *inertia;
+ AtomVecBody::Bonus *bonus;
+
+ int nlocal = atom->nlocal;
+ int newton_pair = force->newton_pair;
+
+ ifirst = dfirst[ibody];
+ iefirst = edfirst[ibody];
+ nei = ednum[ibody];
+
+ rradi = rounded_radius[ibody];
+ rradj = rounded_radius[jbody];
+ contact_dist = rradi + rradj;
+
+ for (ni = 0; ni < nei; ni++) {
+
+ npi1 = static_cast<int>(edge[iefirst+ni][0]);
+ npi2 = static_cast<int>(edge[iefirst+ni][1]);
+
+ // compute the space-fixed coordinates for the vertices of the face
+
+ xi1[0] = x[ibody][0] + discrete[ifirst+npi1][0];
+ xi1[1] = x[ibody][1] + discrete[ifirst+npi1][1];
+ xi1[2] = x[ibody][2] + discrete[ifirst+npi1][2];
+
+ xi2[0] = x[ibody][0] + discrete[ifirst+npi2][0];
+ xi2[1] = x[ibody][1] + discrete[ifirst+npi2][1];
+ xi2[2] = x[ibody][2] + discrete[ifirst+npi2][2];
+
+ // find the projection of the jbody's COM on the edge
+
+ project_pt_line(x[jbody], xi1, xi2, h, d, t);
+
+ if (d > contact_dist + cut_inner) continue;
+ if (t < 0 || t > 1) continue;
+
+ if (fabs(t) < EPSILON) {
+ if (static_cast<int>(discrete[ifirst+npi1][6]) == 1)
+ continue;
+ else {
+ h[0] = xi1[0];
+ h[1] = xi1[1];
+ h[2] = xi1[2];
+ discrete[ifirst+npi1][6] = 1;
+ }
+ }
+
+ if (fabs(t-1) < EPSILON) {
+ if (static_cast<int>(discrete[ifirst+npi2][6]) == 1)
+ continue;
+ else {
+ h[0] = xi2[0];
+ h[1] = xi2[1];
+ h[2] = xi2[2];
+ discrete[ifirst+npi2][6] = 1;
+ }
+ }
+
+ delx = h[0] - x[jbody][0];
+ dely = h[1] - x[jbody][1];
+ delz = h[2] - x[jbody][2];
+ rij = sqrt(delx*delx + dely*dely + delz*delz);
+ R = rij - contact_dist;
+ shift = k_na * cut_inner;
+
+ energy = 0;
+
+ if (R <= 0) { // deformation occurs
+ fpair = -k_n * R - shift;
+ energy = (0.5 * k_n * R + shift) * R;
+ } else if (R <= cut_inner) { // not deforming but cohesive ranges overlap
+ fpair = k_na * R - shift;
+ energy = (-0.5 * k_na * R + shift) * R;
+ } else fpair = 0.0;
+
+ fx = delx*fpair/rij;
+ fy = dely*fpair/rij;
+ fz = delz*fpair/rij;
+
+ if (R <= 0) { // in contact
+
+ // compute the velocity of the vertex in the space-fixed frame
+
+ ibonus = atom->body[ibody];
+ bonus = &avec->bonus[ibonus];
+ quat = bonus->quat;
+ inertia = bonus->inertia;
+ total_velocity(h, x[ibody], v[ibody], angmom[ibody],
+ inertia, quat, vti);
+
+ // relative translational velocity
+
+ vr1 = vti[0] - v[jbody][0];
+ vr2 = vti[1] - v[jbody][1];
+ vr3 = vti[2] - v[jbody][2];
+
+ // normal component
+
+ vnnr = vr1*delx + vr2*dely + vr3*delz;
+ vn1 = delx*vnnr * rsqinv;
+ vn2 = dely*vnnr * rsqinv;
+ vn3 = delz*vnnr * rsqinv;
+
+ // tangential component
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // normal friction term at contact
+
+ fn[0] = -c_n * vn1;
+ fn[1] = -c_n * vn2;
+ fn[2] = -c_n * vn3;
+
+ // tangential friction term at contact, excluding the tangential deformation term
+
+ ft[0] = -c_t * vt1;
+ ft[1] = -c_t * vt2;
+ ft[2] = -c_t * vt3;
+
+ fx += fn[0] + ft[0];
+ fy += fn[1] + ft[1];
+ fz += fn[2] + ft[2];
+ }
+
+ f[ibody][0] += fx;
+ f[ibody][1] += fy;
+ f[ibody][2] += fz;
+ sum_torque(x[ibody], h, fx, fy, fz, torque[ibody]);
+
+ if (newton_pair || jbody < nlocal) {
+ f[jbody][0] -= fx;
+ f[jbody][1] -= fy;
+ f[jbody][2] -= fz;
+ }
+
+ if (evflag) ev_tally_xyz(ibody,jbody,nlocal,newton_pair,
+ energy,0.0,fx,fy,fz,delx,dely,delz);
+ }
+
+}
+
+/* ----------------------------------------------------------------------
+ Determine the interaction mode between i's edges against j's faces
+
+ i = atom i (body i)
+ j = atom j (body j)
+ x = atoms' coordinates
+ f = atoms' forces
+ torque = atoms' torques
+ tag = atoms' tags
+ contact_list = list of contacts
+ num_contacts = number of contacts between i's edges and j's faces
+ Return:
+
+---------------------------------------------------------------------- */
+
+int PairBodyRoundedPolyhedron::edge_against_face(int ibody, int jbody,
+ double k_n, double k_na, double** x, Contact* contact_list,
+ int &num_contacts, double &evdwl, double* facc)
+{
+ int ni,nei,nj,nfj,contact,interact;
+ double rradi,rradj,energy;
+
+ nei = ednum[ibody];
+ rradi = rounded_radius[ibody];
+ nfj = facnum[jbody];
+ rradj = rounded_radius[jbody];
+
+ energy = 0;
+ interact = EF_NONE;
+
+ // loop through body i's edges
+
+ for (ni = 0; ni < nei; ni++) {
+
+ // loop through body j's faces
+
+ for (nj = 0; nj < nfj; nj++) {
+
+ // compute the distance between the face nj to the edge ni
+ #ifdef _POLYHEDRON_DEBUG
+ printf("Compute interaction between face %d of body %d with edge %d of body %d:\n",
+ nj, jbody, ni, ibody);
+ #endif
+
+ interact = interaction_face_to_edge(jbody, nj, x[jbody], rradj,
+ ibody, ni, x[ibody], rradi,
+ k_n, k_na, cut_inner,
+ contact_list, num_contacts,
+ energy, facc);
+ }
+
+ } // end for looping through the edges of body i
+
+ evdwl += energy;
+
+ return interact;
+}
+
+/* ----------------------------------------------------------------------
+ Determine the interaction mode between i's edges against j's edges
+
+ i = atom i (body i)
+ j = atom j (body j)
+ x = atoms' coordinates
+ f = atoms' forces
+ torque = atoms' torques
+ tag = atoms' tags
+ contact_list = list of contacts
+ num_contacts = number of contacts between i's edges and j's edges
+ Return:
+
+---------------------------------------------------------------------- */
+
+int PairBodyRoundedPolyhedron::edge_against_edge(int ibody, int jbody,
+ double k_n, double k_na, double** x, Contact* contact_list,
+ int &num_contacts, double &evdwl, double* facc)
+{
+ int ni,nei,nj,nej,contact,interact;
+ double rradi,rradj,energy;
+
+ nei = ednum[ibody];
+ rradi = rounded_radius[ibody];
+ nej = ednum[jbody];
+ rradj = rounded_radius[jbody];
+
+ energy = 0;
+ interact = EE_NONE;
+
+ // loop through body i's edges
+
+ for (ni = 0; ni < nei; ni++) {
+
+ for (nj = 0; nj < nej; nj++) {
+
+ // compute the distance between the edge nj to the edge ni
+ #ifdef _POLYHEDRON_DEBUG
+ printf("Compute interaction between edge %d of body %d with edge %d of body %d:\n",
+ nj, jbody, ni, ibody);
+ #endif
+
+ interact = interaction_edge_to_edge(ibody, ni, x[ibody], rradi,
+ jbody, nj, x[jbody], rradj,
+ k_n, k_na, cut_inner,
+ contact_list, num_contacts,
+ energy, facc);
+ }
+
+ } // end for looping through the edges of body i
+
+ evdwl += energy;
+
+ return interact;
+}
+
+/* -------------------------------------------------------------------------
+ Compute the interaction between a face of body i and an edge from
+ another body
+ Input:
+ ibody = body i (i.e. atom i)
+ face_index = face index of body i
+ xmi = atom i's coordinates (body i's center of mass)
+ rounded_radius_i = rounded radius of the body i
+ jbody = body i (i.e. atom j)
+ edge_index = coordinate of the tested edge from another body
+ xmj = atom j's coordinates (body j's center of mass)
+ rounded_radius_j = rounded radius of the body j
+ cut_inner = cutoff for vertex-vertex and vertex-edge interaction
+ Output:
+ d = Distance from a point x0 to an edge
+ hi = coordinates of the projection of x0 on the edge
+
+ contact = 0 no contact between the queried edge and the face
+ 1 contact detected
+ return
+ INVALID if the face index is invalid
+ NONE if there is no interaction
+------------------------------------------------------------------------- */
+
+int PairBodyRoundedPolyhedron::interaction_face_to_edge(int ibody,
+ int face_index,
+ double *xmi,
+ double rounded_radius_i,
+ int jbody,
+ int edge_index,
+ double *xmj,
+ double rounded_radius_j,
+ double k_n,
+ double k_na,
+ double cut_inner,
+ Contact* contact_list,
+ int &num_contacts,
+ double &energy,
+ double* facc)
+{
+ if (face_index >= facnum[ibody]) return EF_INVALID;
+
+ int ifirst,iffirst,jfirst,npi1,npi2,npi3;
+ int jefirst,npj1,npj2;
+ double xi1[3],xi2[3],xi3[3],xpj1[3],xpj2[3],ui[3],vi[3],n[3];
+
+ double** x = atom->x;
+ double** v = atom->v;
+ double** f = atom->f;
+ double** torque = atom->torque;
+ double** angmom = atom->angmom;
+
+ ifirst = dfirst[ibody];
+ iffirst = facfirst[ibody];
+ npi1 = static_cast<int>(face[iffirst+face_index][0]);
+ npi2 = static_cast<int>(face[iffirst+face_index][1]);
+ npi3 = static_cast<int>(face[iffirst+face_index][2]);
+
+ // compute the space-fixed coordinates for the vertices of the face
+
+ xi1[0] = xmi[0] + discrete[ifirst+npi1][0];
+ xi1[1] = xmi[1] + discrete[ifirst+npi1][1];
+ xi1[2] = xmi[2] + discrete[ifirst+npi1][2];
+
+ xi2[0] = xmi[0] + discrete[ifirst+npi2][0];
+ xi2[1] = xmi[1] + discrete[ifirst+npi2][1];
+ xi2[2] = xmi[2] + discrete[ifirst+npi2][2];
+
+ xi3[0] = xmi[0] + discrete[ifirst+npi3][0];
+ xi3[1] = xmi[1] + discrete[ifirst+npi3][1];
+ xi3[2] = xmi[2] + discrete[ifirst+npi3][2];
+
+ // find the normal unit vector of the face, ensure it point outward of the body
+
+ MathExtra::sub3(xi2, xi1, ui);
+ MathExtra::sub3(xi3, xi1, vi);
+ MathExtra::cross3(ui, vi, n);
+ MathExtra::norm3(n);
+
+ double xc[3], dot, ans[3];
+ xc[0] = (xi1[0] + xi2[0] + xi3[0])/3.0;
+ xc[1] = (xi1[1] + xi2[1] + xi3[1])/3.0;
+ xc[2] = (xi1[2] + xi2[2] + xi3[2])/3.0;
+ MathExtra::sub3(xc, xmi, ans);
+ dot = MathExtra::dot3(ans, n);
+ if (dot < 0) MathExtra::negate3(n);
+
+ // two ends of the edge from body j
+
+ jfirst = dfirst[jbody];
+ jefirst = edfirst[jbody];
+ npj1 = static_cast<int>(edge[jefirst+edge_index][0]);
+ npj2 = static_cast<int>(edge[jefirst+edge_index][1]);
+
+ xpj1[0] = xmj[0] + discrete[jfirst+npj1][0];
+ xpj1[1] = xmj[1] + discrete[jfirst+npj1][1];
+ xpj1[2] = xmj[2] + discrete[jfirst+npj1][2];
+
+ xpj2[0] = xmj[0] + discrete[jfirst+npj2][0];
+ xpj2[1] = xmj[1] + discrete[jfirst+npj2][1];
+ xpj2[2] = xmj[2] + discrete[jfirst+npj2][2];
+
+ // no interaction if two ends of the edge are on the same side with the COM wrt the face
+
+ if (opposite_sides(n, xi1, xmi, xpj1) == 0 &&
+ opposite_sides(n, xi1, xmi, xpj2) == 0)
+ return EF_NONE;
+
+ // determine the intersection of the edge to the face
+
+ double hi1[3], hi2[3], d1, d2, contact_dist, shift;
+ int inside1 = 0;
+ int inside2 = 0;
+
+ // enum {EF_PARALLEL=0,EF_SAME_SIDE_OF_FACE,EF_INTERSECT_INSIDE,EF_INTERSECT_OUTSIDE};
+ int interact = edge_face_intersect(xi1, xi2, xi3, xpj1, xpj2,
+ hi1, hi2, d1, d2, inside1, inside2);
+
+ inside_polygon(ibody, face_index, xmi, hi1, hi2, inside1, inside2);
+
+ contact_dist = rounded_radius_i + rounded_radius_j;
+ shift = k_na * cut_inner;
+
+ // both endpoints are on the same side of, or parallel to, the face
+ // and both are out of the interaction zone
+
+ if (interact == EF_SAME_SIDE_OF_FACE || interact == EF_PARALLEL) {
+
+ if (d1 > contact_dist + cut_inner && d2 > contact_dist + cut_inner)
+ return EF_NONE;
+
+ int num_outside = 0;
+ int jflag = 1;
+
+ #ifdef _POLYHEDRON_DEBUG
+ if (interact == EF_SAME_SIDE_OF_FACE)
+ printf(" - same side of face\n");
+ else if (interact == EF_PARALLEL)
+ printf(" - parallel\n");
+ printf(" face: xi1 (%f %f %f) xi2 (%f %f %f) xi3 (%f %f %f)\n",
+ xi1[0], xi1[1], xi1[2], xi2[0], xi2[1], xi2[2], xi3[0], xi3[1], xi3[2]);
+ printf(" edge: xpj1 (%f %f %f) xpj2 (%f %f %f)\n",
+ xpj1[0], xpj1[1], xpj1[2], xpj2[0], xpj2[1], xpj2[2]);
+ #endif
+
+ // xpj1 is in the interaction zone
+ // and its projection on the face is inside the triangle
+ // compute vertex-face interaction and accumulate force/torque to both bodies
+
+ if (d1 <= contact_dist + cut_inner) {
+ if (inside1) {
+ if (static_cast<int>(discrete[jfirst+npj1][6]) == 0) {
+ pair_force_and_torque(jbody, ibody, xpj1, hi1, d1, contact_dist,
+ k_n, k_na, shift, x, v, f, torque, angmom,
+ jflag, energy, facc);
+ #ifdef _POLYHEDRON_DEBUG
+ printf(" - compute pair force between vertex %d from edge %d of body %d "
+ "with face %d of body %d: d1 = %f\n",
+ npj1, edge_index, jbody, face_index, ibody, d1);
+ #endif
+
+ if (d1 <= contact_dist) {
+ // store the contact info
+ contact_list[num_contacts].ibody = ibody;
+ contact_list[num_contacts].jbody = jbody;
+ contact_list[num_contacts].xi[0] = hi1[0];
+ contact_list[num_contacts].xi[1] = hi1[1];
+ contact_list[num_contacts].xi[2] = hi1[2];
+ contact_list[num_contacts].xj[0] = xpj1[0];
+ contact_list[num_contacts].xj[1] = xpj1[1];
+ contact_list[num_contacts].xj[2] = xpj1[2];
+ contact_list[num_contacts].type = 0;
+ contact_list[num_contacts].separation = d1 - contact_dist;
+ num_contacts++;
+ }
+
+ discrete[jfirst+npj1][6] = 1;
+ }
+ } else {
+ num_outside++;
+ }
+ }
+
+ // xpj2 is in the interaction zone
+ // and its projection on the face is inside the triangle
+ // compute vertex-face interaction and accumulate force/torque to both bodies
+
+ if (d2 <= contact_dist + cut_inner) {
+ if (inside2) {
+ if (static_cast<int>(discrete[jfirst+npj2][6]) == 0) {
+ pair_force_and_torque(jbody, ibody, xpj2, hi2, d2, contact_dist,
+ k_n, k_na, shift, x, v, f, torque, angmom,
+ jflag, energy, facc);
+ #ifdef _POLYHEDRON_DEBUG
+ printf(" - compute pair force between vertex %d from edge %d of body %d "
+ "with face %d of body %d: d2 = %f\n",
+ npj2, edge_index, jbody, face_index, ibody, d2);
+ #endif
+
+ if (d2 <= contact_dist) {
+ // store the contact info
+ contact_list[num_contacts].ibody = ibody;
+ contact_list[num_contacts].jbody = jbody;
+ contact_list[num_contacts].xi[0] = hi2[0];
+ contact_list[num_contacts].xi[1] = hi2[1];
+ contact_list[num_contacts].xi[2] = hi2[2];
+ contact_list[num_contacts].xj[0] = xpj2[0];
+ contact_list[num_contacts].xj[1] = xpj2[1];
+ contact_list[num_contacts].xj[2] = xpj2[2];
+ contact_list[num_contacts].type = 0;
+ contact_list[num_contacts].separation = d2 - contact_dist;
+ num_contacts++;
+ }
+ discrete[jfirst+npj2][6] = 1;
+ }
+ } else {
+ num_outside++;
+ }
+ }
+
+ // both ends have projection outside of the face
+ // compute interaction between the edge with the three edges of the face
+
+ if (num_outside == 2) {
+
+ #ifdef _POLYHEDRON_DEBUG
+ printf(" - outside = 2\n");
+ printf(" - compute pair force between edge %d of body %d "
+ "with 3 edges of face %d of body %d\n",
+ edge_index, jbody, face_index, ibody);
+ #endif
+
+ interact = EF_INTERSECT_OUTSIDE;
+
+ }
+
+ } else if (interact == EF_INTERSECT_OUTSIDE) {
+
+ // compute interaction between the edge with the three edges of the face
+
+ #ifdef _POLYHEDRON_DEBUG
+ printf(" - intersect outside triangle\n");
+ printf(" - compute pair force between edge %d of body %d "
+ "with face %d of body %d\n", edge_index, jbody, face_index, ibody);
+ printf(" face: xi1 (%f %f %f) xi2 (%f %f %f) xi3 (%f %f %f)\n",
+ xi1[0], xi1[1], xi1[2], xi2[0], xi2[1], xi2[2], xi3[0], xi3[1], xi3[2]);
+ printf(" edge: xpj1 (%f %f %f) xpj2 (%f %f %f)\n",
+ xpj1[0], xpj1[1], xpj1[2], xpj2[0], xpj2[1], xpj2[2]);
+
+ #endif
+ } else if (interact == EF_INTERSECT_INSIDE) {
+
+ }
+
+ return interact;
+}
+
+/* -------------------------------------------------------------------------
+ Compute the distance between an edge of body i and an edge from
+ another body
+ Input:
+ ibody = body i (i.e. atom i)
+ face_index = face index of body i
+ xmi = atom i's coordinates (body i's center of mass)
+ rounded_radius_i = rounded radius of the body i
+ jbody = body i (i.e. atom j)
+ edge_index = coordinate of the tested edge from another body
+ xmj = atom j's coordinates (body j's center of mass)
+ rounded_radius_j = rounded radius of the body j
+ cut_inner = cutoff for vertex-vertex and vertex-edge interaction
+ Output:
+ d = Distance from a point x0 to an edge
+ hi = coordinates of the projection of x0 on the edge
+
+ contact = 0 no contact between the queried edge and the face
+ 1 contact detected
+ return
+ INVALID if the face index is invalid
+ NONE if there is no interaction
+------------------------------------------------------------------------- */
+
+int PairBodyRoundedPolyhedron::interaction_edge_to_edge(int ibody,
+ int edge_index_i,
+ double *xmi,
+ double rounded_radius_i,
+ int jbody,
+ int edge_index_j,
+ double *xmj,
+ double rounded_radius_j,
+ double k_n,
+ double k_na,
+ double cut_inner,
+ Contact* contact_list,
+ int &num_contacts,
+ double &energy,
+ double* facc)
+{
+ int ifirst,iefirst,jfirst,jefirst,npi1,npi2,npj1,npj2,interact;
+ double xi1[3],xi2[3],xpj1[3],xpj2[3];
+ double r,t1,t2,h1[3],h2[3];
+ double contact_dist, shift;
+
+ double** x = atom->x;
+ double** v = atom->v;
+ double** f = atom->f;
+ double** torque = atom->torque;
+ double** angmom = atom->angmom;
+
+ ifirst = dfirst[ibody];
+ iefirst = edfirst[ibody];
+ npi1 = static_cast<int>(edge[iefirst+edge_index_i][0]);
+ npi2 = static_cast<int>(edge[iefirst+edge_index_i][1]);
+
+ // compute the space-fixed coordinates for the edge ends
+
+ xi1[0] = xmi[0] + discrete[ifirst+npi1][0];
+ xi1[1] = xmi[1] + discrete[ifirst+npi1][1];
+ xi1[2] = xmi[2] + discrete[ifirst+npi1][2];
+
+ xi2[0] = xmi[0] + discrete[ifirst+npi2][0];
+ xi2[1] = xmi[1] + discrete[ifirst+npi2][1];
+ xi2[2] = xmi[2] + discrete[ifirst+npi2][2];
+
+ // two ends of the edge from body j
+
+ jfirst = dfirst[jbody];
+ jefirst = edfirst[jbody];
+ npj1 = static_cast<int>(edge[jefirst+edge_index_j][0]);
+ npj2 = static_cast<int>(edge[jefirst+edge_index_j][1]);
+
+ xpj1[0] = xmj[0] + discrete[jfirst+npj1][0];
+ xpj1[1] = xmj[1] + discrete[jfirst+npj1][1];
+ xpj1[2] = xmj[2] + discrete[jfirst+npj1][2];
+
+ xpj2[0] = xmj[0] + discrete[jfirst+npj2][0];
+ xpj2[1] = xmj[1] + discrete[jfirst+npj2][1];
+ xpj2[2] = xmj[2] + discrete[jfirst+npj2][2];
+
+ contact_dist = rounded_radius_i + rounded_radius_j;
+ shift = k_na * cut_inner;
+
+ int jflag = 1;
+ distance_bt_edges(xpj1, xpj2, xi1, xi2, h1, h2, t1, t2, r);
+
+ #ifdef _POLYHEDRON_DEBUG
+ double ui[3],uj[3];
+ MathExtra::sub3(xi1,xi2,ui);
+ MathExtra::norm3(ui);
+ MathExtra::sub3(xpj1,xpj2,uj);
+ MathExtra::norm3(uj);
+ double dot = MathExtra::dot3(ui, uj);
+ printf(" edge npi1 = %d (%f %f %f); npi2 = %d (%f %f %f) vs."
+ " edge npj1 = %d (%f %f %f); npj2 = %d (%f %f %f): "
+ "t1 = %f; t2 = %f; r = %f; dot = %f\n",
+ npi1, xi1[0], xi1[1], xi1[2], npi2, xi2[0], xi2[1], xi2[2],
+ npj1, xpj1[0], xpj1[1], xpj1[2], npj2, xpj2[0], xpj2[1], xpj2[2],
+ t1, t2, r, dot);
+ #endif
+
+ interact = EE_NONE;
+
+ // singularity case, ignore interactions
+
+ if (r < EPSILON) return interact;
+
+ // include the vertices for interactions
+
+ if (t1 >= 0 && t1 <= 1 && t2 >= 0 && t2 <= 1 &&
+ r < contact_dist + cut_inner) {
+ pair_force_and_torque(jbody, ibody, h1, h2, r, contact_dist,
+ k_n, k_na, shift, x, v, f, torque, angmom,
+ jflag, energy, facc);
+
+ interact = EE_INTERACT;
+ if (r <= contact_dist) {
+ // store the contact info
+ contact_list[num_contacts].ibody = ibody;
+ contact_list[num_contacts].jbody = jbody;
+ contact_list[num_contacts].xi[0] = h2[0];
+ contact_list[num_contacts].xi[1] = h2[1];
+ contact_list[num_contacts].xi[2] = h2[2];
+ contact_list[num_contacts].xj[0] = h1[0];
+ contact_list[num_contacts].xj[1] = h1[1];
+ contact_list[num_contacts].xj[2] = h1[2];
+ contact_list[num_contacts].type = 1;
+ contact_list[num_contacts].separation = r - contact_dist;
+ num_contacts++;
+ }
+ } else {
+
+ }
+
+ return interact;
+}
+
+/* ----------------------------------------------------------------------
+ Compute forces and torques between two bodies caused by the interaction
+ between a pair of points on either bodies (similar to sphere-sphere)
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::pair_force_and_torque(int ibody, int jbody,
+ double* pi, double* pj, double r, double contact_dist,
+ double k_n, double k_na, double shift, double** x,
+ double** v, double** f, double** torque, double** angmom,
+ int jflag, double& energy, double* facc)
+{
+ double delx,dely,delz,R,fx,fy,fz,fpair;
+
+ delx = pi[0] - pj[0];
+ dely = pi[1] - pj[1];
+ delz = pi[2] - pj[2];
+ R = r - contact_dist;
+ if (R <= 0) { // deformation occurs
+ fpair = -k_n * R - shift;
+ energy += (0.5 * k_n * R + shift) * R;
+ } else if (R <= cut_inner) { // not deforming but cohesive ranges overlap
+ fpair = k_na * R - shift;
+ energy += (-0.5 * k_na * R + shift) * R;
+ } else fpair = 0.0;
+
+ fx = delx*fpair/r;
+ fy = dely*fpair/r;
+ fz = delz*fpair/r;
+
+ #ifdef _POLYHEDRON_DEBUG
+ printf(" - R = %f; r = %f; k_na = %f; shift = %f; fpair = %f;"
+ " energy = %f; jflag = %d\n", R, r, k_na, shift, fpair,
+ energy, jflag);
+ #endif
+
+ if (R <= 0) {
+
+ // contact: accumulate normal and tangential contact force components
+
+ contact_forces(ibody, jbody, pi, pj, delx, dely, delz, fx, fy, fz,
+ x, v, angmom, f, torque, facc);
+ } else {
+
+ // accumulate force and torque to both bodies directly
+
+ f[ibody][0] += fx;
+ f[ibody][1] += fy;
+ f[ibody][2] += fz;
+ sum_torque(x[ibody], pi, fx, fy, fz, torque[ibody]);
+
+ facc[0] += fx; facc[1] += fy; facc[2] += fz;
+
+ if (jflag) {
+ f[jbody][0] -= fx;
+ f[jbody][1] -= fy;
+ f[jbody][2] -= fz;
+ sum_torque(x[jbody], pj, -fx, -fy, -fz, torque[jbody]);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ Compute contact forces between two bodies
+ modify the force stored at the vertex and edge in contact by j_a
+ sum forces and torque to the corresponding bodies
+ fx,fy,fz = unscaled cohesive forces
+ fn = normal friction component
+ ft = tangential friction component (-c_t * v_t)
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::contact_forces(int ibody, int jbody,
+ double *xi, double *xj, double delx, double dely, double delz,
+ double fx, double fy, double fz, double** x, double** v, double** angmom,
+ double** f, double** torque, double* facc)
+{
+ int ibonus,jbonus;
+ double rsq,rsqinv,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
+ double fn[3],ft[3],vi[3],vj[3];
+ double *quat, *inertia;
+ AtomVecBody::Bonus *bonus;
+
+ // compute the velocity of the vertex in the space-fixed frame
+
+ ibonus = atom->body[ibody];
+ bonus = &avec->bonus[ibonus];
+ quat = bonus->quat;
+ inertia = bonus->inertia;
+ total_velocity(xi, x[ibody], v[ibody], angmom[ibody],
+ inertia, quat, vi);
+
+ // compute the velocity of the point on the edge in the space-fixed frame
+
+ jbonus = atom->body[jbody];
+ bonus = &avec->bonus[jbonus];
+ quat = bonus->quat;
+ inertia = bonus->inertia;
+ total_velocity(xj, x[jbody], v[jbody], angmom[jbody],
+ inertia, quat, vj);
+
+ // vector pointing from the contact point on ibody to that on jbody
+
+ rsq = delx*delx + dely*dely + delz*delz;
+ rsqinv = 1.0/rsq;
+
+ // relative translational velocity
+
+ vr1 = vi[0] - vj[0];
+ vr2 = vi[1] - vj[1];
+ vr3 = vi[2] - vj[2];
+
+ // normal component
+
+ vnnr = vr1*delx + vr2*dely + vr3*delz;
+ vn1 = delx*vnnr * rsqinv;
+ vn2 = dely*vnnr * rsqinv;
+ vn3 = delz*vnnr * rsqinv;
+
+ // tangential component
+
+ vt1 = vr1 - vn1;
+ vt2 = vr2 - vn2;
+ vt3 = vr3 - vn3;
+
+ // normal friction term at contact
+
+ fn[0] = -c_n * vn1;
+ fn[1] = -c_n * vn2;
+ fn[2] = -c_n * vn3;
+
+ // tangential friction term at contact
+ // excluding the tangential deformation term for now
+
+ ft[0] = -c_t * vt1;
+ ft[1] = -c_t * vt2;
+ ft[2] = -c_t * vt3;
+
+ // cohesive forces will be scaled by j_a after contact area is computed
+ // mu * fne = tangential friction deformation during gross sliding
+ // see Eq. 4, Fraige et al.
+
+ fx = fn[0] + ft[0] + mu * fx;
+ fy = fn[1] + ft[1] + mu * fy;
+ fz = fn[2] + ft[2] + mu * fz;
+
+ f[ibody][0] += fx;
+ f[ibody][1] += fy;
+ f[ibody][2] += fz;
+ sum_torque(x[ibody], xi, fx, fy, fz, torque[ibody]);
+
+ f[jbody][0] -= fx;
+ f[jbody][1] -= fy;
+ f[jbody][2] -= fz;
+ sum_torque(x[jbody], xj, -fx, -fy, -fz, torque[jbody]);
+
+ facc[0] += fx; facc[1] += fy; facc[2] += fz;
+
+ #ifdef _POLYHEDRON_DEBUG
+ printf("contact ibody = %d: f = %f %f %f; torque = %f %f %f\n", ibody,
+ f[ibody][0], f[ibody][1], f[ibody][2],
+ torque[ibody][0], torque[ibody][1], torque[ibody][2]);
+ printf("contact jbody = %d: f = %f %f %f; torque = %f %f %f\n", jbody,
+ f[jbody][0], f[jbody][1], f[jbody][2],
+ torque[jbody][0], torque[jbody][1], torque[jbody][2]);
+ #endif
+}
+
+/* ----------------------------------------------------------------------
+ Rescale the forces and torques for all the contacts
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::rescale_cohesive_forces(double** x,
+ double** f, double** torque, Contact* contact_list, int &num_contacts,
+ double contact_area, double k_n, double k_na, double* facc)
+{
+ int m,ibody,jbody;
+ double delx,dely,delz,fx,fy,fz,R,fpair,r,shift;
+ double j_a = contact_area / (num_contacts * A_ua);
+ if (j_a < 1.0) j_a = 1.0;
+
+ shift = k_na * cut_inner;
+
+ for (m = 0; m < num_contacts; m++) {
+ ibody = contact_list[m].ibody;
+ jbody = contact_list[m].jbody;
+
+ delx = contact_list[m].xi[0] - contact_list[m].xj[0];
+ dely = contact_list[m].xi[1] - contact_list[m].xj[1];
+ delz = contact_list[m].xi[2] - contact_list[m].xj[2];
+ r = sqrt(delx*delx + dely*dely + delz*delz);
+ R = contact_list[m].separation;
+ if (R <= 0) { // deformation occurs
+ fpair = -k_n * R - shift;
+ } else if (R <= cut_inner) { // not deforming but cohesive ranges overlap
+ fpair = k_na * R - shift;
+ } else fpair = 0.0;
+
+ fpair *= j_a;
+ fx = delx*fpair/r;
+ fy = dely*fpair/r;
+ fz = delz*fpair/r;
+
+ f[ibody][0] += fx;
+ f[ibody][1] += fy;
+ f[ibody][2] += fz;
+ sum_torque(x[ibody], contact_list[m].xi, fx, fy, fz, torque[ibody]);
+
+ f[jbody][0] -= fx;
+ f[jbody][1] -= fy;
+ f[jbody][2] -= fz;
+ sum_torque(x[jbody], contact_list[m].xj, -fx, -fy, -fz, torque[jbody]);
+
+ facc[0] += fx; facc[1] += fy; facc[2] += fz;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ Accumulate torque to body from the force f=(fx,fy,fz) acting at point x
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::sum_torque(double* xm, double *x, double fx,
+ double fy, double fz, double* torque)
+{
+ double rx = x[0] - xm[0];
+ double ry = x[1] - xm[1];
+ double rz = x[2] - xm[2];
+ double tx = ry * fz - rz * fy;
+ double ty = rz * fx - rx * fz;
+ double tz = rx * fy - ry * fx;
+ torque[0] += tx;
+ torque[1] += ty;
+ torque[2] += tz;
+}
+
+/* ----------------------------------------------------------------------
+ Test if two points a and b are in opposite sides of a plane defined by
+ a normal vector n and a point x0
+------------------------------------------------------------------------- */
+
+int PairBodyRoundedPolyhedron::opposite_sides(double* n, double* x0,
+ double* a, double* b)
+{
+ double m_a = n[0]*(a[0] - x0[0])+n[1]*(a[1] - x0[1])+n[2]*(a[2] - x0[2]);
+ double m_b = n[0]*(b[0] - x0[0])+n[1]*(b[1] - x0[1])+n[2]*(b[2] - x0[2]);
+ // equal to zero when either a or b is on the plane
+ if (m_a * m_b <= 0)
+ return 1;
+ else
+ return 0;
+}
+
+/* ----------------------------------------------------------------------
+ Test if a line segment defined by two points a and b intersects with
+ a triangle defined by three points x1, x2 and x3
+------------------------------------------------------------------------- */
+
+int PairBodyRoundedPolyhedron::edge_face_intersect(double* x1, double* x2,
+ double* x3, double* a, double* b, double* h_a, double* h_b,
+ double& d_a, double& d_b, int& inside_a, int& inside_b)
+{
+ double s[3], u[3], v[3], n[3];
+
+ // line director
+
+ MathExtra::sub3(b, a, s);
+
+ // plane normal vector
+
+ MathExtra::sub3(x2, x1, u);
+ MathExtra::sub3(x3, x1, v);
+ MathExtra::cross3(u, v, n);
+ MathExtra::norm3(n);
+
+ // find the projection of a and b to the plane and the corresponding distances
+
+ project_pt_plane(a, x1, x2, x3, h_a, d_a, inside_a);
+
+ project_pt_plane(b, x1, x2, x3, h_b, d_b, inside_b);
+
+ // check if the line segment is parallel to the plane
+
+ double dot = MathExtra::dot3(s, n);
+ if (fabs(dot) < EPSILON) return EF_PARALLEL;
+
+ // solve for the intersection between the line and the plane
+
+ double m[3][3], invm[3][3], p[3], ans[3];
+ m[0][0] = -s[0];
+ m[0][1] = u[0];
+ m[0][2] = v[0];
+
+ m[1][0] = -s[1];
+ m[1][1] = u[1];
+ m[1][2] = v[1];
+
+ m[2][0] = -s[2];
+ m[2][1] = u[2];
+ m[2][2] = v[2];
+
+ MathExtra::sub3(a, x1, p);
+ MathExtra::invert3(m, invm);
+ MathExtra::matvec(invm, p, ans);
+
+ // p is reused for the intersection point
+ // s = b - a
+
+ double t = ans[0];
+ p[0] = a[0] + s[0] * t;
+ p[1] = a[1] + s[1] * t;
+ p[2] = a[2] + s[2] * t;
+
+ // check if p is inside the triangle, excluding the edges and vertices
+ // the edge-edge and edge-vertices are handled separately
+
+ int inside = 0;
+ if (ans[1] > 0 && ans[2] > 0 && ans[1] + ans[2] < 1)
+ inside = 1;
+
+ int interact;
+ if (t < 0 || t > 1) {
+ interact = EF_SAME_SIDE_OF_FACE;
+ } else {
+ if (inside == 1)
+ interact = EF_INTERSECT_INSIDE;
+ else
+ interact = EF_INTERSECT_OUTSIDE;
+ }
+
+ return interact;
+}
+
+/* ----------------------------------------------------------------------
+ Find the projection of q on the plane defined by point p and the normal
+ unit vector n: q_proj = q - dot(q - p, n) * n
+ and the distance d from q to the plane
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::project_pt_plane(const double* q,
+ const double* p, const double* n,
+ double* q_proj, double &d)
+{
+ double dot, ans[3], n_p[3];
+ n_p[0] = n[0]; n_p[1] = n[1]; n_p[2] = n[2];
+ MathExtra::sub3(q, p, ans);
+ dot = MathExtra::dot3(ans, n_p);
+ MathExtra::scale3(dot, n_p);
+ MathExtra::sub3(q, n_p, q_proj);
+ MathExtra::sub3(q, q_proj, ans);
+ d = MathExtra::len3(ans);
+}
+
+/* ----------------------------------------------------------------------
+ Check if points q1 and q2 are inside a convex polygon, i.e. a face of
+ a polyhedron
+ ibody = atom i's index
+ face_index = face index of the body
+ xmi = atom i's coordinates
+ q1 = tested point on the face (e.g. the projection of a point)
+ q2 = another point (can be NULL) for face-edge intersection
+ Output:
+ inside1 = 1 if q1 is inside the polygon, 0 otherwise
+ inside2 = 1 if q2 is inside the polygon, 0 otherwise
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::inside_polygon(int ibody, int face_index,
+ double* xmi, const double* q1, const double* q2,
+ int& inside1, int& inside2)
+
+{
+ int i,n,ifirst,iffirst,npi1,npi2;
+ double xi1[3],xi2[3],u[3],v[3],costheta,anglesum1,anglesum2,magu,magv;
+
+ ifirst = dfirst[ibody];
+ iffirst = facfirst[ibody];
+ anglesum1 = anglesum2 = 0;;
+ for (i = 0; i < MAX_FACE_SIZE; i++) {
+ npi1 = static_cast<int>(face[iffirst+face_index][i]);
+ if (npi1 < 0) break;
+ n = i + 1;
+ if (n <= MAX_FACE_SIZE - 1) {
+ npi2 = static_cast<int>(face[iffirst+face_index][n]);
+ if (npi2 < 0) npi2 = static_cast<int>(face[iffirst+face_index][0]);
+ } else {
+ npi2 = static_cast<int>(face[iffirst+face_index][0]);
+ }
+
+ xi1[0] = xmi[0] + discrete[ifirst+npi1][0];
+ xi1[1] = xmi[1] + discrete[ifirst+npi1][1];
+ xi1[2] = xmi[2] + discrete[ifirst+npi1][2];
+
+ xi2[0] = xmi[0] + discrete[ifirst+npi2][0];
+ xi2[1] = xmi[1] + discrete[ifirst+npi2][1];
+ xi2[2] = xmi[2] + discrete[ifirst+npi2][2];
+
+ MathExtra::sub3(xi1,q1,u);
+ MathExtra::sub3(xi2,q1,v);
+ magu = MathExtra::len3(u);
+ magv = MathExtra::len3(v);
+
+ // the point is at either vertices
+
+ if (magu * magv < EPSILON) inside1 = 1;
+ else {
+ costheta = MathExtra::dot3(u,v)/(magu*magv);
+ anglesum1 += acos(costheta);
+ }
+
+ if (q2 != NULL) {
+ MathExtra::sub3(xi1,q2,u);
+ MathExtra::sub3(xi2,q2,v);
+ magu = MathExtra::len3(u);
+ magv = MathExtra::len3(v);
+ if (magu * magv < EPSILON) inside2 = 1;
+ else {
+ costheta = MathExtra::dot3(u,v)/(magu*magv);
+ anglesum2 += acos(costheta);
+ }
+ }
+ }
+
+ if (fabs(anglesum1 - MY_2PI) < EPSILON) inside1 = 1;
+ else inside1 = 0;
+
+ if (q2 != NULL) {
+ if (fabs(anglesum2 - MY_2PI) < EPSILON) inside2 = 1;
+ else inside2 = 0;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ Find the projection of q on the plane defined by 3 points x1, x2 and x3
+ returns the distance d from q to the plane and whether the projected
+ point is inside the triangle defined by (x1, x2, x3)
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::project_pt_plane(const double* q,
+ const double* x1, const double* x2, const double* x3, double* q_proj,
+ double &d, int& inside)
+{
+ double u[3],v[3],n[3];
+
+ // plane normal vector
+
+ MathExtra::sub3(x2, x1, u);
+ MathExtra::sub3(x3, x1, v);
+ MathExtra::cross3(u, v, n);
+ MathExtra::norm3(n);
+
+ // solve for the intersection between the line and the plane
+
+ double m[3][3], invm[3][3], p[3], ans[3];
+ m[0][0] = -n[0];
+ m[0][1] = u[0];
+ m[0][2] = v[0];
+
+ m[1][0] = -n[1];
+ m[1][1] = u[1];
+ m[1][2] = v[1];
+
+ m[2][0] = -n[2];
+ m[2][1] = u[2];
+ m[2][2] = v[2];
+
+ MathExtra::sub3(q, x1, p);
+ MathExtra::invert3(m, invm);
+ MathExtra::matvec(invm, p, ans);
+
+ double t = ans[0];
+ q_proj[0] = q[0] + n[0] * t;
+ q_proj[1] = q[1] + n[1] * t;
+ q_proj[2] = q[2] + n[2] * t;
+
+ // check if the projection point is inside the triangle
+ // exclude the edges and vertices
+ // edge-sphere and sphere-sphere interactions are handled separately
+
+ inside = 0;
+ if (ans[1] > 0 && ans[2] > 0 && ans[1] + ans[2] < 1) {
+ inside = 1;
+ }
+
+ // distance from q to q_proj
+
+ MathExtra::sub3(q, q_proj, ans);
+ d = MathExtra::len3(ans);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::project_pt_line(const double* q,
+ const double* xi1, const double* xi2, double* h, double& d, double& t)
+{
+ double u[3],v[3],r[3],s;
+
+ MathExtra::sub3(xi2, xi1, u);
+ MathExtra::norm3(u);
+ MathExtra::sub3(q, xi1, v);
+
+ s = MathExtra::dot3(u, v);
+ h[0] = xi1[0] + s * u[0];
+ h[1] = xi1[1] + s * u[1];
+ h[2] = xi1[2] + s * u[2];
+
+ MathExtra::sub3(q, h, r);
+ d = MathExtra::len3(r);
+
+ if (fabs(xi2[0] - xi1[0]) > 0)
+ t = (h[0] - xi1[0])/(xi2[0] - xi1[0]);
+ else if (fabs(xi2[1] - xi1[1]) > 0)
+ t = (h[1] - xi1[1])/(xi2[1] - xi1[1]);
+ else if (fabs(xi2[2] - xi1[2]) > 0)
+ t = (h[2] - xi1[2])/(xi2[2] - xi1[2]);
+}
+
+/* ----------------------------------------------------------------------
+ compute the shortest distance between two edges (line segments)
+ x1, x2: two endpoints of the first edge
+ x3, x4: two endpoints of the second edge
+ h1: the end point of the shortest segment perpendicular to both edges
+ on the line (x1;x2)
+ h2: the end point of the shortest segment perpendicular to both edges
+ on the line (x3;x4)
+ t1: fraction of h1 in the segment (x1,x2)
+ t2: fraction of h2 in the segment (x3,x4)
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::distance_bt_edges(const double* x1,
+ const double* x2, const double* x3, const double* x4,
+ double* h1, double* h2, double& t1, double& t2, double& r)
+{
+ double u[3],v[3],n[3],dot;
+
+ // set the default returned values
+
+ t1 = -2;
+ t2 = 2;
+ r = 0;
+
+ // find the edge unit directors and their dot product
+
+ MathExtra::sub3(x2, x1, u);
+ MathExtra::norm3(u);
+ MathExtra::sub3(x4, x3, v);
+ MathExtra::norm3(v);
+ dot = MathExtra::dot3(u,v);
+ dot = fabs(dot);
+
+ // check if two edges are parallel
+ // find the two ends of the overlapping segment, if any
+
+ if (fabs(dot - 1.0) < EPSILON) {
+
+ double s1,s2,x13[3],x23[3],x13h[3];
+ double t13,t23,t31,t41,x31[3],x41[3];
+
+ MathExtra::sub3(x1,x3,x13); // x13 = x1 - x3
+ MathExtra::sub3(x2,x3,x23); // x23 = x2 - x3
+
+ s1 = MathExtra::dot3(x13,v);
+ x13h[0] = x13[0] - s1*v[0];
+ x13h[1] = x13[1] - s1*v[1];
+ x13h[2] = x13[2] - s1*v[2];
+ r = MathExtra::len3(x13h);
+
+ // x13 is the projection of x1 on x3-x4
+
+ x13[0] = x3[0] + s1*v[0];
+ x13[1] = x3[1] + s1*v[1];
+ x13[2] = x3[2] + s1*v[2];
+
+ // x23 is the projection of x2 on x3-x4
+
+ s2 = MathExtra::dot3(x23,v);
+ x23[0] = x3[0] + s2*v[0];
+ x23[1] = x3[1] + s2*v[1];
+ x23[2] = x3[2] + s2*v[2];
+
+ // find the fraction of the projection points on the edges
+
+ if (fabs(x4[0] - x3[0]) > 0)
+ t13 = (x13[0] - x3[0])/(x4[0] - x3[0]);
+ else if (fabs(x4[1] - x3[1]) > 0)
+ t13 = (x13[1] - x3[1])/(x4[1] - x3[1]);
+ else if (fabs(x4[2] - x3[2]) > 0)
+ t13 = (x13[2] - x3[2])/(x4[2] - x3[2]);
+
+ if (fabs(x4[0] - x3[0]) > 0)
+ t23 = (x23[0] - x3[0])/(x4[0] - x3[0]);
+ else if (fabs(x4[1] - x3[1]) > 0)
+ t23 = (x23[1] - x3[1])/(x4[1] - x3[1]);
+ else if (fabs(x4[2] - x3[2]) > 0)
+ t23 = (x23[2] - x3[2])/(x4[2] - x3[2]);
+
+ if (fabs(x23[0] - x13[0]) > 0)
+ t31 = (x3[0] - x13[0])/(x23[0] - x13[0]);
+ else if (fabs(x23[1] - x13[1]) > 0)
+ t31 = (x3[1] - x13[1])/(x23[1] - x13[1]);
+ else if (fabs(x23[2] - x13[2]) > 0)
+ t31 = (x3[2] - x13[2])/(x23[2] - x13[2]);
+
+ // x31 is the projection of x3 on x1-x2
+
+ x31[0] = x1[0] + t31*(x2[0] - x1[0]);
+ x31[1] = x1[1] + t31*(x2[1] - x1[1]);
+ x31[2] = x1[2] + t31*(x2[2] - x1[2]);
+
+ if (fabs(x23[0] - x13[0]) > 0)
+ t41 = (x4[0] - x13[0])/(x23[0] - x13[0]);
+ else if (fabs(x23[1] - x13[1]) > 0)
+ t41 = (x4[1] - x13[1])/(x23[1] - x13[1]);
+ else if (fabs(x23[2] - x13[2]) > 0)
+ t41 = (x4[2] - x13[2])/(x23[2] - x13[2]);
+
+ // x41 is the projection of x4 on x1-x2
+
+ x41[0] = x1[0] + t41*(x2[0] - x1[0]);
+ x41[1] = x1[1] + t41*(x2[1] - x1[1]);
+ x41[2] = x1[2] + t41*(x2[2] - x1[2]);
+
+ // determine two ends from the overlapping segments
+
+ int n1 = 0;
+ int n2 = 0;
+ if (t13 >= 0 && t13 <= 1) {
+ h1[0] = x1[0];
+ h1[1] = x1[1];
+ h1[2] = x1[2];
+ h2[0] = x13[0];
+ h2[1] = x13[1];
+ h2[2] = x13[2];
+ t1 = 0;
+ t2 = t13;
+ n1++;
+ n2++;
+ }
+ if (t23 >= 0 && t23 <= 1) {
+ if (n1 == 0) {
+ h1[0] = x2[0];
+ h1[1] = x2[1];
+ h1[2] = x2[2];
+ h2[0] = x23[0];
+ h2[1] = x23[1];
+ h2[2] = x23[2];
+ t1 = 1;
+ t2 = t23;
+ n1++;
+ n2++;
+ } else {
+ h1[0] = (x1[0]+x2[0])/2;
+ h1[1] = (x1[1]+x2[1])/2;
+ h1[2] = (x1[2]+x2[2])/2;
+ h2[0] = (x13[0]+x23[0])/2;
+ h2[1] = (x13[1]+x23[1])/2;
+ h2[2] = (x13[2]+x23[2])/2;
+ t1 = 0.5;
+ t2 = (t13+t23)/2;
+ n1++;
+ n2++;
+ }
+ }
+
+ if (n1 == 0 && n2 == 0) {
+ if (t31 >= 0 && t31 <= 1) {
+ h1[0] = x31[0];
+ h1[1] = x31[1];
+ h1[2] = x31[2];
+ h2[0] = x3[0];
+ h2[1] = x3[1];
+ h2[2] = x3[2];
+ t1 = t31;
+ t2 = 0;
+ n1++;
+ n2++;
+ }
+ if (t41 >= 0 && t41 <= 1) {
+ if (n1 == 0) {
+ h1[0] = x41[0];
+ h1[1] = x41[1];
+ h1[2] = x41[2];
+ h2[0] = x4[0];
+ h2[1] = x4[1];
+ h2[2] = x4[2];
+ t1 = t41;
+ t2 = 1;
+ n1++;
+ n2++;
+ } else {
+ h1[0] = (x31[0]+x41[0])/2;
+ h1[1] = (x31[1]+x41[1])/2;
+ h1[2] = (x31[2]+x41[2])/2;
+ h2[0] = (x3[0]+x4[0])/2;
+ h2[1] = (x3[1]+x4[1])/2;
+ h2[2] = (x3[2]+x4[2])/2;
+ t1 = (t31+t41)/2;
+ t2 = 0.5;
+ n1++;
+ n2++;
+ }
+ }
+ }
+
+ // if n1 == 0 and n2 == 0 at this point,
+ // which means no overlapping segments bt two parallel edges,
+ // return the default values of t1 and t2
+
+ return;
+
+ }
+
+ // find the vector n perpendicular to both edges
+
+ MathExtra::cross3(u, v, n);
+ MathExtra::norm3(n);
+
+ // find the intersection of the line (x3,x4) and the plane (x1,x2,n)
+ // s = director of the line (x3,x4)
+ // n_p = plane normal vector of the plane (x1,x2,n)
+
+ double s[3], n_p[3];
+ MathExtra::sub3(x4, x3, s);
+ MathExtra::sub3(x2, x1, u);
+ MathExtra::cross3(u, n, n_p);
+ MathExtra::norm3(n_p);
+
+ // solve for the intersection between the line and the plane
+
+ double m[3][3], invm[3][3], p[3], ans[3];
+ m[0][0] = -s[0];
+ m[0][1] = u[0];
+ m[0][2] = n[0];
+
+ m[1][0] = -s[1];
+ m[1][1] = u[1];
+ m[1][2] = n[1];
+
+ m[2][0] = -s[2];
+ m[2][1] = u[2];
+ m[2][2] = n[2];
+
+ MathExtra::sub3(x3, x1, p);
+ MathExtra::invert3(m, invm);
+ MathExtra::matvec(invm, p, ans);
+
+ t2 = ans[0];
+ h2[0] = x3[0] + s[0] * t2;
+ h2[1] = x3[1] + s[1] * t2;
+ h2[2] = x3[2] + s[2] * t2;
+
+ project_pt_plane(h2, x1, n, h1, r);
+
+ if (fabs(x2[0] - x1[0]) > 0)
+ t1 = (h1[0] - x1[0])/(x2[0] - x1[0]);
+ else if (fabs(x2[1] - x1[1]) > 0)
+ t1 = (h1[1] - x1[1])/(x2[1] - x1[1]);
+ else if (fabs(x2[2] - x1[2]) > 0)
+ t1 = (h1[2] - x1[2])/(x2[2] - x1[2]);
+}
+
+/* ----------------------------------------------------------------------
+ Calculate the total velocity of a point (vertex, a point on an edge):
+ vi = vcm + omega ^ (p - xcm)
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::total_velocity(double* p, double *xcm,
+ double* vcm, double *angmom, double *inertia, double *quat, double* vi)
+{
+ double r[3],omega[3],ex_space[3],ey_space[3],ez_space[3];
+ r[0] = p[0] - xcm[0];
+ r[1] = p[1] - xcm[1];
+ r[2] = p[2] - xcm[2];
+ MathExtra::q_to_exyz(quat,ex_space,ey_space,ez_space);
+ MathExtra::angmom_to_omega(angmom,ex_space,ey_space,ez_space,
+ inertia,omega);
+ vi[0] = omega[1]*r[2] - omega[2]*r[1] + vcm[0];
+ vi[1] = omega[2]*r[0] - omega[0]*r[2] + vcm[1];
+ vi[2] = omega[0]*r[1] - omega[1]*r[0] + vcm[2];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::sanity_check()
+{
+
+ double x1[3],x2[3],x3[3],x4[3],h_a[3],h_b[3],d_a,d_b,u[3],v[3],n[3];
+ double a[3],b[3],t_a,t_b;
+ int inside_a, inside_b;
+
+ x1[0] = 0; x1[1] = 3; x1[2] = 0;
+ x2[0] = 3; x2[1] = 0; x2[2] = 0;
+ x3[0] = 4; x3[1] = 3; x3[2] = 0;
+ x4[0] = 5; x4[1] = 3; x4[2] = 0;
+
+ a[0] = 0; a[1] = 0; a[2] = 0;
+ b[0] = 4; b[1] = 0; b[2] = 0;
+
+ project_pt_line(a, x1, x2, h_a, d_a, t_a);
+ project_pt_line(b, x1, x2, h_b, d_b, t_b);
+/*
+ printf("h_a: %f %f %f; h_b: %f %f %f; t_a = %f; t_b = %f; d = %f; d_b = %f\n",
+ h_a[0], h_a[1], h_a[2], h_b[0], h_b[1], h_b[2], t_a, t_b, d_a, d_b);
+*/
+/*
+ int mode = edge_face_intersect(x1, x2, x3, a, b, h_a, h_b, d_a, d_b,
+ inside_a, inside_b);
+
+ MathExtra::sub3(x2, x1, u);
+ MathExtra::sub3(x3, x1, v);
+ MathExtra::cross3(u, v, n);
+ MathExtra::norm3(n);
+*/
+/*
+ project_pt_plane(a, x1, x2, x3, h_a, d_a, inside_a);
+ printf("h_a: %f %f %f; d = %f: inside %d\n",
+ h_a[0], h_a[1], h_a[2], d_a, inside_a);
+ project_pt_plane(b, x1, x2, x3, h_b, d_b, inside_b);
+ printf("h_b: %f %f %f; d = %f: inside %d\n",
+ h_b[0], h_b[1], h_b[2], d_b, inside_b);
+*/
+/*
+ distance_bt_edges(x1, x2, x3, x4, h_a, h_b, t_a, t_b, d_a);
+ printf("h_a: %f %f %f; h_b: %f %f %f; t_a = %f; t_b = %f; d = %f\n",
+ h_a[0], h_a[1], h_a[2], h_b[0], h_b[1], h_b[2], t_a, t_b, d_a);
+*/
+}
+
diff --git a/src/BODY/pair_body_rounded_polyhedron.h b/src/BODY/pair_body_rounded_polyhedron.h
new file mode 100644
index 0000000000000000000000000000000000000000..d7ee1f013e66a837948bdd6bddfb5464630e50fe
--- /dev/null
+++ b/src/BODY/pair_body_rounded_polyhedron.h
@@ -0,0 +1,176 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(body/rounded/polyhedron,PairBodyRoundedPolyhedron)
+
+#else
+
+#ifndef LMP_PAIR_BODY_ROUNDED_POLYHEDRON_H
+#define LMP_PAIR_BODY_ROUNDED_POLYHEDRON_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairBodyRoundedPolyhedron : public Pair {
+ public:
+ PairBodyRoundedPolyhedron(class LAMMPS *);
+ ~PairBodyRoundedPolyhedron();
+ void compute(int, int);
+ void settings(int, char **);
+ void coeff(int, char **);
+ void init_style();
+ double init_one(int, int);
+
+ struct Contact {
+ int ibody, jbody; // body (i.e. atom) indices (not tags)
+ int type; // 0 = VERTEX-FACE; 1 = EDGE-EDGE
+ double fx,fy,fz; // unscaled cohesive forces at contact
+ double xi[3]; // coordinates of the contact point on ibody
+ double xj[3]; // coordinates of the contact point on jbody
+ double separation; // contact surface separation
+ };
+
+ protected:
+ double **k_n; // normal repulsion strength
+ double **k_na; // normal attraction strength
+ double c_n; // normal damping coefficient
+ double c_t; // tangential damping coefficient
+ double mu; // normal friction coefficient during gross sliding
+ double A_ua; // characteristic contact area
+ double cut_inner; // cutoff for interaction between vertex-edge surfaces
+
+ class AtomVecBody *avec;
+ class BodyRoundedPolyhedron *bptr;
+
+ double **discrete; // list of all sub-particles for all bodies
+ int ndiscrete; // number of discretes in list
+ int dmax; // allocated size of discrete list
+ int *dnum; // number of discretes per line, 0 if uninit
+ int *dfirst; // index of first discrete per each line
+ int nmax; // allocated size of dnum,dfirst vectors
+
+ double **edge; // list of all edge for all bodies
+ int nedge; // number of edge in list
+ int edmax; // allocated size of edge list
+ int *ednum; // number of edges per line, 0 if uninit
+ int *edfirst; // index of first edge per each line
+ int ednummax; // allocated size of ednum,edfirst vectors
+
+ double **face; // list of all edge for all bodies
+ int nface; // number of faces in list
+ int facmax; // allocated size of face list
+ int *facnum; // number of faces per line, 0 if uninit
+ int *facfirst; // index of first face per each line
+ int facnummax; // allocated size of facnum,facfirst vectors
+
+ double *enclosing_radius; // enclosing radii for all bodies
+ double *rounded_radius; // rounded radii for all bodies
+ double *maxerad; // per-type maximum enclosing radius
+
+ void allocate();
+ void body2space(int);
+
+ int edge_against_face(int ibody, int jbody, double k_n, double k_na,
+ double** x, Contact* contact_list, int &num_contacts,
+ double &evdwl, double* facc);
+ int edge_against_edge(int ibody, int jbody, double k_n, double k_na,
+ double** x,Contact* contact_list, int &num_contacts,
+ double &evdwl, double* facc);
+ void sphere_against_sphere(int ibody, int jbody, double delx, double dely, double delz,
+ double rsq, double k_n, double k_na,
+ double** v, double** f, int evflag);
+ void sphere_against_face(int ibody, int jbody,
+ double k_n, double k_na, double** x, double** v,
+ double** f, double** torque, double** angmom, int evflag);
+ void sphere_against_edge(int ibody, int jbody,
+ double k_n, double k_na, double** x, double** v,
+ double** f, double** torque, double** angmom, int evflag);
+
+ int interaction_face_to_edge(int ibody, int face_index, double* xmi,
+ double rounded_radius_i, int jbody, int edge_index,
+ double* xmj, double rounded_radius_j,
+ double k_n, double k_na, double cut_inner,
+ Contact* contact_list, int &num_contacts,
+ double& energy, double* facc);
+ int interaction_edge_to_edge(int ibody, int edge_index_i, double* xmi,
+ double rounded_radius_i, int jbody, int edge_index_j,
+ double* xmj, double rounded_radius_j,
+ double k_n, double k_na, double cut_inner,
+ Contact* contact_list, int &num_contacts,
+ double& energy, double* facc);
+
+ void contact_forces(int ibody, int jbody, double *xi, double *xj,
+ double delx, double dely, double delz, double fx, double fy, double fz,
+ double** x, double** v, double** angmom, double** f, double** torque,
+ double* facc);
+
+ void pair_force_and_torque(int ibody, int jbody, double* pi, double* pj,
+ double r, double contact_dist, double k_n,
+ double k_na, double shift, double** x, double** v,
+ double** f, double** torque, double** angmom,
+ int jflag, double& energy, double* facc);
+ void rescale_cohesive_forces(double** x, double** f, double** torque,
+ Contact* contact_list, int &num_contacts,
+ double contact_area, double k_n, double k_na, double* facc);
+
+ void sum_torque(double* xm, double *x, double fx, double fy, double fz, double* torque);
+ int opposite_sides(double* n, double* x0, double* a, double* b);
+ int edge_face_intersect(double* x1, double* x2, double* x3, double* a, double* b,
+ double* hi1, double* hi2, double& d1, double& d2,
+ int& inside_a, int& inside_b);
+ void project_pt_plane(const double* q, const double* p,
+ const double* n, double* q_proj, double &d);
+ void project_pt_plane(const double* q, const double* x1, const double* x2,
+ const double* x3, double* q_proj, double &d, int& inside);
+ void project_pt_line(const double* q, const double* xi1, const double* xi2,
+ double* h, double& d, double& t);
+ void inside_polygon(int ibody, int face_index, double* xmi,
+ const double* q1, const double* q2, int& inside1, int& inside2);
+
+ void distance_bt_edges(const double* x1, const double* x2,
+ const double* x3, const double* x4,
+ double* h1, double* h2, double& t1, double& t2, double& r);
+ void total_velocity(double* p, double *xcm, double* vcm, double *angmom,
+ double *inertia, double *quat, double* vi);
+ void sanity_check();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory. Check the input script or data file.
+
+E: Pair body/rounded/polyhedron requires atom style body rounded/polyhedron
+
+Self-explanatory.
+
+E: Pair body requires body style rounded/polyhedron
+
+This pair style is specific to the rounded/polyhedron body style.
+
+*/