enables body style rounded/polygon to write to data files

doc/src/Howto_body.rst

@@ -219,21 +219,19 @@ in the *Bodies* section of the data file:
    x1 y1 z1
    ...
    xN yN zN
-   i j j k k ...
    diameter
 
-where M = 6 + 3\*N + 2\*N + 1, and N is the number of vertices in the
-body particle.
+where M = 6 + 3\*N + 1, and N is the number of vertices in the body
+particle.
 
 The integer line has a single value N.  The floating point line(s)
-list 6 moments of inertia followed by the coordinates of the N
+list 6 moments of inertia, followed by the coordinates of the N
 vertices (x1 to zN) as 3N values (with z = 0.0 for each), followed by
-2N vertex indices corresponding to the end points of the N edges,
-followed by a single diameter value = the rounded diameter of the
-circle that surrounds each vertex. The diameter value can be different
-for each body particle. These floating-point values can be listed on
-as many lines as you wish; see the :doc:`read_data <read_data>` command
-for more details.
+a diameter value = the rounded diameter of the circle that surrounds
+each vertex. The diameter value can be different for each body
+particle. These floating-point values can be listed on as many lines
+as you wish; see the :doc:`read_data <read_data>` command for more
+details.
 
 The 6 moments of inertia (ixx,iyy,izz,ixy,ixz,iyz) should be the
 values consistent with the current orientation of the rigid body
@@ -260,10 +258,6 @@ is consistent with the 6 moments of inertia: ixx iyy izz ixy ixz iyz =
    -0.7071 0.7071 0
    0.7071 0.7071 0
    0.7071 -0.7071 0
-   0 1
-   1 2
-   2 3
-   3 0
    1.0
 
 A rod in 2D, whose length is 4.0, mass 1.0, rounded at two ends

src/BODY/body_rounded_polygon.cpp

@@ -26,6 +26,7 @@
 #include "math_extra.h"
 #include "memory.h"
 #include "error.h"
+#include "fmt/format.h"
 
 using namespace LAMMPS_NS;
 
@@ -181,7 +182,8 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
   bonus->ninteger = 1;
   bonus->ivalue = icp->get(bonus->iindex);
   bonus->ivalue[0] = nsub;
-  bonus->ndouble = 3*nsub + 2*nsub + 1 + 1;
+  if (nsub < 3) bonus->ndouble = 3*nsub + 2 + 1 + 1;
+  else bonus->ndouble = 3*nsub + 2*nsub + 1 + 1;
   bonus->dvalue = dcp->get(bonus->ndouble,bonus->dindex);
 
   // diagonalize inertia tensor
@@ -223,7 +225,7 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
   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
 
@@ -239,10 +241,13 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
   // find the enclosing radius of the body from the maximum displacement
 
   int i,m;
-  double delta[3], rsq, erad, rrad;
-  double erad2 = 0;
+  double rsq,erad,rrad;
+  double delta[3];
+    
+  double erad2 = 0.0;
   int j = 6;
   int k = 0;
+  
   for (i = 0; i < nsub; i++) {
     delta[0] = dfile[j];
     delta[1] = dfile[j+1];
@@ -256,17 +261,17 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
     k += 3;
   }
 
-  // .. the next 2*nsub elements are edge ends
-
-  int nedges;
-  if (nsub == 1) { // spheres
-    nedges = 0;
+  // the next 2 or 2*nsub elements are edge ends
+  // the final two values are the enclosing radius and rounded radius
+  // set atom->radius = enclosing + rounded radii (except for spheres)
+  
+  // spheres have just 1 edge
+  
+  if (nsub == 1) {
     bonus->dvalue[k] = 0;
-    *(&bonus->dvalue[k]+1) = 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;
@@ -276,42 +281,36 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
 
     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;
-    }
+  // rods have just 1 edge
+    
+  } else if (nsub == 2) {
+    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++) {
+  // polygons have Nsub edges
+    
+  } else {
+    for (i = 0; i < nsub; i++) {
       bonus->dvalue[k] = i;
       m = i+1;
-      if (m == nedges) m = 0;
-      *(&bonus->dvalue[k]+1) = m;
+      if (m == nsub) 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;
@@ -327,8 +326,64 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
 
 int BodyRoundedPolygon::pack_data_body(tagint atomID, int ibonus, double *buf)
 {
+  int m,ilast;
+  double values[3],p[3][3],pdiag[3][3],ispace[3][3];
+
   AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
-  return 0;
+
+  double *quat = bonus->quat;
+  double *inertia = bonus->inertia;
+  int *ivalue = bonus->ivalue;
+  double *dvalue = bonus->dvalue;
+
+  int nsub = ivalue[0];
+
+  if (buf) {
+
+    // ID ninteger ndouble
+
+    m = 0;
+    buf[m++] = ubuf(atomID).d;
+    buf[m++] = ubuf(1).d;
+    buf[m++] = ubuf(6 + 3*nsub + 1).d;
+
+    // single integer nsub
+
+    buf[m++] = ubuf(nsub).d;
+
+    // 6 moments of inertia
+
+    MathExtra::quat_to_mat(quat,p);
+    MathExtra::times3_diag(p,inertia,pdiag);
+    MathExtra::times3_transpose(pdiag,p,ispace);
+
+    buf[m++] = ispace[0][0];
+    buf[m++] = ispace[1][1];
+    buf[m++] = ispace[2][2];
+    buf[m++] = ispace[0][1];
+    buf[m++] = ispace[0][2];
+    buf[m++] = ispace[1][2];
+
+    // 3*nsub particle coords = displacement from COM in box frame
+
+    for (int i = 0; i < nsub; i++) {
+      MathExtra::matvec(p,&dvalue[3*i],values);
+      buf[m++] = values[0];
+      buf[m++] = values[1];
+      buf[m++] = values[2];
+    }
+
+    // rounded diameter = 2 * last dvalue = rounded radius
+    // for nsub = 1,2: skip one edge and enclosing radius
+    // for nsub > 2: skip Nsub edges and enclosing radius
+
+    if (nsub < 3) ilast = 3*nsub + 2 + 1;
+    else ilast = 3*nsub + 2*nsub + 1;
+    buf[m++] = 2.0 * dvalue[ilast];
+    
+  } else m = 3 + 1 + 6 + 3*nsub + 1;
+
+  return m;
 }
 
 /* ----------------------------------------------------------------------
@@ -337,7 +392,33 @@ int BodyRoundedPolygon::pack_data_body(tagint atomID, int ibonus, double *buf)
 
 int BodyRoundedPolygon::write_data_body(FILE *fp, double *buf)
 {
-  return 0;
+  int m = 0;
+
+  // atomID ninteger ndouble
+
+  fmt::print(fp,"{} {} {}\n",ubuf(buf[m]).i,ubuf(buf[m+1]).i,ubuf(buf[m+2]).i);
+  m += 3;
+
+  const int nsub = (int) ubuf(buf[m++]).i;
+  fmt::print(fp,"{}\n",nsub);
+
+  // inertia
+
+  fmt::print(fp,"{} {} {} {} {} {}\n",
+             buf[m+0],buf[m+1],buf[m+2],buf[m+3],buf[m+4],buf[m+5]);
+  m += 6;
+
+  // nsub vertices
+  
+  for (int i = 0; i < nsub; i++, m+=3)
+    fmt::print(fp,"{} {} {}\n",buf[m],buf[m+1],buf[m+2]);
+
+  // rounded diameter
+  
+  double diameter = buf[m++];
+  fmt::print(fp,"{}\n",diameter);
+
+  return m;
 }
 
 /* ----------------------------------------------------------------------