write_data support for body style rounded/polyhedron

bond/non-body (point particle) interactions.  More details of each

style are described below.

More styles may be added in the future.  See the :doc:`Modify body <Modify_body>` doc page for details on how to add a new body

style to the code.

More styles may be added in the future.  See the :doc:`Modify body

<Modify_body>` doc page for details on how to add a new body style to

the code.

----------

The :doc:`pair_style body/nparticle <pair_body_nparticle>` command can be used

with this body style to compute body/body and body/non-body interactions.

For output purposes via the :doc:`compute body/local <compute_body_local>` and :doc:`dump local <dump>`

commands, this body style produces one datum for each of the N

sub-particles in a body particle.  The datum has 3 values:

.. parsed-literal::

   1 = x position of sub-particle

   2 = y position of sub-particle

   3 = z position of sub-particle

These values are the current position of the sub-particle within the

simulation domain, not a displacement from the center-of-mass (COM) of

the body particle itself.  These values are calculated using the

current COM and orientation of the body particle.

For images created by the :doc:`dump image <dump_image>` command, if the

*body* keyword is set, then each body particle is drawn as a

collection of spheres, one for each sub-particle.  The size of each

sphere is determined by the *bflag1* parameter for the *body* keyword.

The *bflag2* argument is ignored.

----------

**Specifics of body style rounded/polygon:**

structures used internally by each particle.

When the :doc:`read_data <read_data>` command reads a data file for this

body style, the following information must be provided for each entry

body style, the following information must be provided for each body

in the *Bodies* section of the data file:

.. parsed-literal::

   1 2 3 4

   diameter

where M = 6 + 3\*N + 2\*E + 4\*F + 1, and N is the number of vertices in

the body particle, E = number of edges, F = number of faces.

where M = 6 + 3\*N + 2\*E + 4\*F + 1, and N is the number of vertices

in the body particle, E = number of edges, F = number of faces.  For N

= 1 or 2, the format is simpler.  E and F are ignored and no edges or

faces are listed, so that M = 6 + 3\*N + 1.

The integer line has three values: number of vertices (N), number of

edges (E) and number of faces (F). The floating point line(s) list 6

faces.  The last value is the diameter value = the rounded diameter of

the sphere 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.  Because the

maximum number of vertices per face is hard-coded to be 4

(i.e. quadrilaterals), faces with more than 4 vertices need to be

split into triangles or quadrilaterals.  For triangular faces, the

last vertex index should be set to -1.

The ordering of the 4 vertices within a face should follow

the right-hand rule so that the normal vector of the face points

outwards from the center of mass.

listed on as many lines as you wish; see the :doc:`read_data

<read_data>` command for more details.

Note that vertices are numbered from 0 to N-1 inclusive.  The 2

vertices in each edge can be in any order.  Faces can be triangles or

quadrilaterals.  In both cases 4 vertices must be specified.  For a

triangle the 4th vertex is -1.  The 4 vertices within each triangle or

quadrilateral face should be ordered by the right-hand rule so that

the normal vector of the face points outwards from the center of mass.

For polyhedron with faces with more than 4 vertices, you should

split the complex face into multiple simple faces, each of

which is a triangle or quadrilateral.

The 6 moments of inertia (ixx,iyy,izz,ixy,ixz,iyz) should be the

values consistent with the current orientation of the rigid body

.. parsed-literal::

   1 1 13

   2

   1 3 13

   2 0 0 

   0 1.33333 1.33333 0 0 0

   -2 0 0

   2 0 0

.. parsed-literal::

   1 1 10

   1

   1 3 10

   1 0 0

   0.9 0.9 0.9 0 0 0

   0 0 0

   3.0

The :doc:`pair_style body/rounded/polhedron <pair_body_rounded_polyhedron>` command can

be used with this body style to compute body/body interactions.  The

:doc:`fix wall/body/polyhedron <fix_wall_body_polygon>` command can be

used with this body style to compute the interaction of body particles

with a wall.

The :doc:`pair_style body/rounded/polhedron

<pair_body_rounded_polyhedron>` command can be used with this body

style to compute body/body interactions.  The :doc:`fix

wall/body/polyhedron <fix_wall_body_polygon>` command can be used with

this body style to compute the interaction of body particles with a

wall.

----------

For output purposes via the :doc:`compute body/local <compute_body_local>` and :doc:`dump local <dump>`

commands, this body style produces one datum for each of the N

sub-particles in a body particle.  The datum has 3 values:

**Output specifics for all body styles:**

For the :doc:`compute body/local <compute_body_local>` and :doc:`dump

local <dump>` commands, all 3 of the body styles described on his page

produces one datum for each of the N vertices (of sub-particles) in a

body particle.  The datum has 3 values:

.. parsed-literal::

   1 = x position of vertex

   1 = x position of vertex (or sub-particle)

   2 = y position of vertex

   3 = z position of vertex

the body particle itself.  These values are calculated using the

current COM and orientation of the body particle.

For images created by the :doc:`dump image <dump_image>` command, if the

*body* keyword is set, then each body particle is drawn as a polygon

consisting of N line segments.  Note that the line segments are drawn

between the N vertices, which does not correspond exactly to the

physical extent of the body (because the :doc:`pair_style rounded/polygon <pair_body_rounded_polygon>` defines finite-size

spheres at those point and the line segments between the spheres are

tangent to the spheres).  The drawn diameter of each line segment is

determined by the *bflag1* parameter for the *body* keyword.  The

*bflag2* argument is ignored.

The :doc:`dump image <dump_image>` command and its *body* keyword can

be used to render body particles.

For the *nparticle* body style, each body is drawn as a

collection of spheres, one for each sub-particle.  The size of each

sphere is determined by the *bflag1* parameter for the *body* keyword.

The *bflag2* argument is ignored.

For the *rounded/polygon* body style, each body is drawn as a polygon

with N line segments.  For the *rounded/polyhedron* body style, each

face of each body is drawn as a polygon with N line segments.  The

drawn diameter of each line segment is determined by the *bflag1*

parameter for the *body* keyword.  The *bflag2* argument is ignored.

Note that for both the *rounded/polygon* and *rounded/polyhedron*

styles, line segments are drawn between the pairs of vertices.

Depending on the diameters of the line segments this may be slightly

different than the physical extent of the body as calculated by the

:doc:`pair_style rounded/polygon <pair_body_rounded_polygon>` or

:doc:`pair_style rounded/polyhedron <pair_body_rounded_polyhedron>`

commands.  Conceptually, the pair styles define the surface of a 2d or

3d body by lines or planes that are tangent to the finite-size spheres

of specified diameter which are placed on each vertex position.

----------

    buf[m++] = ubuf(1).d;

    buf[m++] = ubuf(6 + 3*nsub).d;

    // single integer Nsub

    // single integer nsub

    buf[m++] = ubuf(nsub).d;

    buf[m++] = ispace[0][2];

    buf[m++] = ispace[1][2];

    // 3*Nsub particle coords = displacement from COM in box frame

    // 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+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]);

#include "math_extra.h"

#include "memory.h"

#include "error.h"

#include "fmt/format.h"

using namespace LAMMPS_NS;

  size_forward = 0;

  // 3 integers: 1 for no. of vertices, 1 for no. of edges, 1 for no. of faces

  // 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

  // 3 integers: nvertices, nedges, nfaces

  // doubles = 3*nvertices + 2*nedge + MAX_FACE_SIZE*nfaces +

  //           1 double for enclosing radius + 1 double for rounded radius

  size_border = 3 + 3*nmax + 2*nmax + MAX_FACE_SIZE*nmax + 1 + 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

  return nedges; //(nvertices+nfaces-2); // Euler formula: V-E+F=2

}

/* ---------------------------------------------------------------------- */

/* ---------------------------------------------------------------------- */

int BodyRoundedPolyhedron::unpack_border_body(AtomVecBody::Bonus *bonus,

                                           double *buf)

int BodyRoundedPolyhedron::unpack_border_body(AtomVecBody::Bonus *bonus, double *buf)

{

  int nsub = static_cast<int> (buf[0]);

  int ned = static_cast<int> (buf[1]);

               "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

  // nsub == 1,2:

  //   6 for inertia + 3*nsub + 1 for rounded radius

  // nsub > 2:

  //   6 for inertia + 3*nsub + 2*nedges + MAX_FACE_SIZE*nfaces + 1 for rounded radius

  int nedges,nentries;

  if (nsub == 1 || nsub == 2) {

  if (nsub < 3) {

    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->ivalue[0] = nsub;

  bonus->ivalue[1] = ned;

  bonus->ivalue[2] = nfac;

  if (nsub == 1 || nsub == 2) bonus->ndouble = 3*nsub + 2*nsub + 1 + 1;

  if (nsub < 3) bonus->ndouble = 3*nsub + 2 + 1 + 1;

  else bonus->ndouble = 3*nsub + 2*nedges + MAX_FACE_SIZE*nfac + 1 + 1;

  bonus->dvalue = dcp->get(bonus->ndouble,bonus->dindex);

    k += 3;

  }

  // .. the next 2*nsub elements are edge ends

  // the next 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) { // spheres

    nedges = 0;

  if (nsub == 1) {

    bonus->dvalue[k] = 0;

    *(&bonus->dvalue[k]+1) = 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

    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++) {

  // rods have just 1 edge

  } else if (nsub == 2) {

    bonus->dvalue[k] = 0;

      *(&bonus->dvalue[k]+1) = 1;

    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

  // polyhedra have Nedges and Nfaces

  } else {

    // edges

    for (i = 0; i < nedges; i++) {

      bonus->dvalue[k] = dfile[j];

      *(&bonus->dvalue[k]+1) = dfile[j+1];

      bonus->dvalue[k+1] = dfile[j+1];

      k += 2;

      j += 2;

    }

    for (i = 0; i < nfac; i++) {

      for (m = 0; m < MAX_FACE_SIZE; m++)

        *(&bonus->dvalue[k]+m) = dfile[j+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;

int BodyRoundedPolyhedron::pack_data_body(tagint atomID, int ibonus, double *buf)

{

  int i,j,m;

  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];

  int nedge = ivalue[1];

  int nface = ivalue[2];

  if (buf) {

    // ID ninteger ndouble

    m = 0;

    buf[m++] = ubuf(atomID).d;

    buf[m++] = ubuf(3).d;

    if (nsub < 3) buf[m++] = ubuf(6 + 3*nsub + 1).d;

    else buf[m++] = ubuf(6 + 3*nsub + 2*nedge + MAX_FACE_SIZE*nface + 1).d;

    // 3 integers nsub,nedge,nface

    buf[m++] = ubuf(nsub).d;

    buf[m++] = ubuf(nedge).d;

    buf[m++] = ubuf(nface).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 (i = 0; i < nsub; i++) {

      MathExtra::matvec(p,&dvalue[3*i],values);

      buf[m++] = values[0];

      buf[m++] = values[1];

      buf[m++] = values[2];

    }

    // 2*nedge edge indices

    // 4*nface face indices

    j = 3*nsub;

    if (nsub < 3) j += 2;

    else {

      for (i = 0; i < nedge; i++) {

	buf[m++] = static_cast<int> (dvalue[j++]);

	buf[m++] = static_cast<int> (dvalue[j++]);

      }

      for (i = 0; i < nface; i++) {

	buf[m++] = static_cast<int> (dvalue[j++]);

	buf[m++] = static_cast<int> (dvalue[j++]);

	buf[m++] = static_cast<int> (dvalue[j++]);

	buf[m++] = static_cast<int> (dvalue[j++]);

      }

    }

    // rounded diameter = 2 * last dvalue = rounded radius

    // j+1 to skip enclosing radius

    buf[m++] = 2.0 * dvalue[j+1];

  } else {

    m = 3 + 3 + 6 + 3*nsub + 1;

    if (nsub > 2) m += 2*nedge + MAX_FACE_SIZE*nface;

  }

  return m;

}

/* ----------------------------------------------------------------------

int BodyRoundedPolyhedron::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;

  // nvert, nedge, nface

  const int nsub = (int) ubuf(buf[m++]).i;

  const int nedge = (int) ubuf(buf[m++]).i;

  const int nface = (int) ubuf(buf[m++]).i;

  fmt::print(fp,"{} {} {}\n",nsub,nedge,nface);

  // 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]);

  // nedge 2-tuples and nface 4-tuples

  // unless nsub = 1 or 2

  if (nsub > 2) {

    for (int i = 0; i < nedge; i++, m+=2)

      fmt::print(fp,"{} {}\n",static_cast<int> (buf[m]),static_cast<int> (buf[m+1]));

    for (int i = 0; i < nface; i++, m+=4)

      fmt::print(fp,"{} {} {} {}\n",

		 static_cast<int> (buf[m]),static_cast<int> (buf[m+1]),

		 static_cast<int> (buf[m+2]),static_cast<int> (buf[m+3]));

  }

  // rounded diameter

  double diameter = buf[m++];

  fmt::print(fp,"{}\n",diameter);

  return m;

}

/* ----------------------------------------------------------------------