add body support, change bonus buf to 1d so can be variable length

This commit is contained in:
Steve Plimpton
2020-07-07 11:34:30 -06:00
parent 7918919d30
commit 955bff7c4e
20 changed files with 234 additions and 138 deletions

View File

@ -21,6 +21,7 @@
#include "force.h"
#include "memory.h"
#include "error.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
@ -192,6 +193,92 @@ void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
}
}
/* ----------------------------------------------------------------------
pack data struct for one body into buf for writing to data file
if buf is NULL, just return buffer size
------------------------------------------------------------------------- */
int BodyNparticle::pack_data_body(tagint atomID, int ibonus, double *buf)
{
int m;
double values[3],p[3][3],pdiag[3][3],ispace[3][3];
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
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).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];
}
} else m = 3 + 1 + 6 + 3*nsub;
return m;
}
/* ----------------------------------------------------------------------
write info for one body to data file
------------------------------------------------------------------------- */
int BodyNparticle::write_data_body(FILE *fp, double *buf)
{
int m = 0;
tagint atomID = (tagint) ubuf(buf[m++]).i;
int ninteger = (int) ubuf(buf[m++]).i;
int ndouble = (int) ubuf(buf[m++]).i;
fmt::print(fp,"{} {} {}\n",atomID,ninteger,ndouble);
int nsub = (int) ubuf(buf[m++]).i;
fmt::print(fp,"{}\n",nsub);
double *inertia = &buf[m];
fmt::print(fp,"{} {} {} {} {} {}\n",
inertia[0],inertia[1],inertia[2],inertia[3],inertia[4],inertia[5]);
m += 6;
for (int i = 0; i < nsub; i++)
fmt::print(fp,"{} {} {}\n",buf[m++],buf[m++],buf[m++]);
return m;
}
/* ----------------------------------------------------------------------
return radius of body particle defined by ifile/dfile params
params are ordered as in data file

View File

@ -35,6 +35,8 @@ class BodyNparticle : public Body {
int pack_border_body(struct AtomVecBody::Bonus *, double *);
int unpack_border_body(struct AtomVecBody::Bonus *, double *);
void data_body(int, int, int, int *, double *);
int pack_data_body(tagint, int, double *);
int write_data_body(FILE *, double *);
double radius_body(int, int, int *, double *);
int noutrow(int);

View File

@ -320,6 +320,26 @@ void BodyRoundedPolygon::data_body(int ibonus, int ninteger, int ndouble,
}
}
/* ----------------------------------------------------------------------
pack data struct for one body into buf for writing to data file
if buf is NULL, just return buffer size
------------------------------------------------------------------------- */
int BodyRoundedPolygon::pack_data_body(tagint atomID, int ibonus, double *buf)
{
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
return 0;
}
/* ----------------------------------------------------------------------
write info for one body to data file
------------------------------------------------------------------------- */
int BodyRoundedPolygon::write_data_body(FILE *fp, double *buf)
{
return 0;
}
/* ----------------------------------------------------------------------
return radius of body particle defined by ifile/dfile params
params are ordered as in data file

View File

@ -39,6 +39,8 @@ class BodyRoundedPolygon : public Body {
int pack_border_body(struct AtomVecBody::Bonus *, double *);
int unpack_border_body(struct AtomVecBody::Bonus *, double *);
void data_body(int, int, int, int *, double *);
int pack_data_body(tagint, int, double *);
int write_data_body(FILE *, double *);
double radius_body(int, int, int *, double *);
int noutrow(int);

View File

@ -378,6 +378,26 @@ void BodyRoundedPolyhedron::data_body(int ibonus, int ninteger, int ndouble,
}
}
/* ----------------------------------------------------------------------
pack data struct for one body into buf for writing to data file
if buf is NULL, just return buffer size
------------------------------------------------------------------------- */
int BodyRoundedPolyhedron::pack_data_body(tagint atomID, int ibonus, double *buf)
{
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
return 0;
}
/* ----------------------------------------------------------------------
write info for one body to data file
------------------------------------------------------------------------- */
int BodyRoundedPolyhedron::write_data_body(FILE *fp, double *buf)
{
return 0;
}
/* ----------------------------------------------------------------------
return radius of body particle defined by ifile/dfile params
params are ordered as in data file

View File

@ -41,6 +41,8 @@ class BodyRoundedPolyhedron : public Body {
int pack_border_body(struct AtomVecBody::Bonus *, double *);
int unpack_border_body(struct AtomVecBody::Bonus *, double *);
void data_body(int, int, int, int *, double *);
int pack_data_body(tagint, int, double *);
int write_data_body(FILE *, double *);
double radius_body(int, int, int *, double *);
int noutrow(int);

View File

@ -2304,16 +2304,6 @@ void AtomVec::write_improper(FILE *fp, int n, tagint **buf, int index)
}
}
/* ----------------------------------------------------------------------
return size_data_bonus
only AtomVecHybrid overrides this, so it can select which bonus data via flag
------------------------------------------------------------------------- */
int AtomVec::size_data_bonus_query(int /*flag*/)
{
return size_data_bonus;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
------------------------------------------------------------------------- */

View File

@ -145,9 +145,8 @@ class AtomVec : protected Pointers {
virtual int pack_improper(tagint **);
virtual void write_improper(FILE *, int, tagint **, int);
virtual int size_data_bonus_query(int);
virtual int pack_data_bonus(double **, int) {return 0;}
virtual void write_data_bonus(FILE *, int, double **, int) {}
virtual int pack_data_bonus(double *, int) {return 0;}
virtual void write_data_bonus(FILE *, int, double *, int) {}
virtual int property_atom(char *) {return -1;}
virtual void pack_property_atom(int, double *, int, int) {}

View File

@ -581,26 +581,20 @@ void AtomVecBody::pack_data_pre(int ilocal)
/* ----------------------------------------------------------------------
pack bonus body info for writing to data file
if buf is NULL, just return count of bodies
if buf is NULL, just return buffer size
------------------------------------------------------------------------- */
int AtomVecBody::pack_data_bonus(double **buf, int /*flag*/)
int AtomVecBody::pack_data_bonus(double *buf, int /*flag*/)
{
int i,j;
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
// NOTE: needs to call Body sub-class to fill buffer
int m = 0;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
if (buf) {
buf[m][0] = ubuf(tag[i]).d;
j = body[i];
}
m++;
m += bptr->pack_data_body(tag[i],body[i],buf);
}
return m;
@ -610,11 +604,11 @@ int AtomVecBody::pack_data_bonus(double **buf, int /*flag*/)
write bonus body info to data file
------------------------------------------------------------------------- */
void AtomVecBody::write_data_bonus(FILE *fp, int n, double **buf, int /*flag*/)
void AtomVecBody::write_data_bonus(FILE *fp, int n, double *buf, int /*flag*/)
{
// NOTE: needs to call Body sub-class to do the write
for (int i = 0; i < n; i++) {
int i = 0;
while (i < n) {
i += bptr->write_data_body(fp,&buf[i]);
}
}

View File

@ -63,8 +63,8 @@ class AtomVecBody : public AtomVec {
void pack_data_pre(int);
void pack_data_post(int);
int pack_data_bonus(double **, int);
void write_data_bonus(FILE *, int, double **, int);
int pack_data_bonus(double *, int);
void write_data_bonus(FILE *, int, double *, int);
// methods used by other classes to query/set body info

View File

@ -484,10 +484,10 @@ void AtomVecEllipsoid::pack_data_post(int ilocal)
/* ----------------------------------------------------------------------
pack bonus ellipsoid info for writing to data file
if buf is NULL, just return count of ellipsoids
if buf is NULL, just return buffer size
------------------------------------------------------------------------- */
int AtomVecEllipsoid::pack_data_bonus(double **buf, int /*flag*/)
int AtomVecEllipsoid::pack_data_bonus(double *buf, int /*flag*/)
{
int i,j;
@ -498,17 +498,16 @@ int AtomVecEllipsoid::pack_data_bonus(double **buf, int /*flag*/)
for (i = 0; i < nlocal; i++) {
if (ellipsoid[i] < 0) continue;
if (buf) {
buf[m][0] = ubuf(tag[i]).d;
buf[m++] = ubuf(tag[i]).d;
j = ellipsoid[i];
buf[m][1] = bonus[j].shape[0];
buf[m][2] = bonus[j].shape[1];
buf[m][3] = bonus[j].shape[2];
buf[m][4] = bonus[j].quat[0];
buf[m][5] = bonus[j].quat[1];
buf[m][6] = bonus[j].quat[2];
buf[m][7] = bonus[j].quat[3];
}
m++;
buf[m++] = bonus[j].shape[0];
buf[m++] = bonus[j].shape[1];
buf[m++] = bonus[j].shape[2];
buf[m++] = bonus[j].quat[0];
buf[m++] = bonus[j].quat[1];
buf[m++] = bonus[j].quat[2];
buf[m++] = bonus[j].quat[3];
} else m += size_data_bonus;
}
return m;
@ -518,12 +517,14 @@ int AtomVecEllipsoid::pack_data_bonus(double **buf, int /*flag*/)
write bonus ellipsoid info to data file
------------------------------------------------------------------------- */
void AtomVecEllipsoid::write_data_bonus(FILE *fp, int n, double **buf, int /*flag*/)
void AtomVecEllipsoid::write_data_bonus(FILE *fp, int n, double *buf, int /*flag*/)
{
for (int i = 0; i < n; i++) {
fmt::print(fp,"{} {} {} {} {} {} {} {}",
(tagint) ubuf(buf[i][0]).i,buf[i][1],buf[i][2],buf[i][3],
buf[i][4],buf[i][5],buf[i][6],buf[i][7]);
int i = 0;
while (i < n) {
fmt::print(fp,"{} {} {} {} {} {} {} {}\n",
(tagint) ubuf(buf[i]).i,
buf[i+1],buf[i+2],buf[i+3],buf[i+4],buf[i+5],buf[i+6],buf[i+7]);
i += size_data_bonus;
}
}

View File

@ -56,8 +56,8 @@ class AtomVecEllipsoid : public AtomVec {
void pack_data_pre(int);
void pack_data_post(int);
int pack_data_bonus(double **, int);
void write_data_bonus(FILE *, int, double **, int);
int pack_data_bonus(double *, int);
void write_data_bonus(FILE *, int, double *, int);
// unique to AtomVecEllipsoid

View File

@ -459,54 +459,33 @@ void AtomVecHybrid::pack_data_post(int ilocal)
}
/* ----------------------------------------------------------------------
return size_data_bonus
match flag to sub-style
pack bonus info for writing to data file, match flag to sub-style
------------------------------------------------------------------------- */
int AtomVecHybrid::size_data_bonus_query(int flag)
int AtomVecHybrid::pack_data_bonus(double *buf, int flag)
{
for (int k = 0; k < nstyles; k++) {
if (flag == ELLIPSOID && strcmp(keywords[k],"ellipsoid") == 0)
return styles[k]->size_data_bonus;
if (flag == LINE && strcmp(keywords[k],"line") == 0)
return styles[k]->size_data_bonus;
if (flag == TRIANGLE && strcmp(keywords[k],"tri") == 0)
return styles[k]->size_data_bonus;
// this will not work, body style does not set size_data_bonus
// if (flag == BODY && strcmp(keywords[k],"body") == 0)
// return styles[k]->size_data_bonus;
}
}
if (flag == ELLIPSOID && strcmp(keywords[k],"ellipsoid") != 0) continue;
if (flag == LINE && strcmp(keywords[k],"line") != 0) continue;
if (flag == TRIANGLE && strcmp(keywords[k],"tri") != 0) continue;
if (flag == BODY && strcmp(keywords[k],"body") != 0) continue;
/* ----------------------------------------------------------------------
pack bonus info for writing to data file
match flag to sub-style
------------------------------------------------------------------------- */
int AtomVecHybrid::pack_data_bonus(double **buf, int flag)
{
for (int k = 0; k < nstyles; k++) {
if (flag == ELLIPSOID && strcmp(keywords[k],"ellipsoid") == 0)
return styles[k]->pack_data_bonus(buf,flag);
if (flag == LINE && strcmp(keywords[k],"line") == 0)
return styles[k]->pack_data_bonus(buf,flag);
if (flag == TRIANGLE && strcmp(keywords[k],"tri") == 0)
return styles[k]->pack_data_bonus(buf,flag);
}
}
/* ----------------------------------------------------------------------
write bonus info to data file
write bonus info to data file, match flag to sub-style
------------------------------------------------------------------------- */
void AtomVecHybrid::write_data_bonus(FILE *fp, int n, double **buf, int flag)
void AtomVecHybrid::write_data_bonus(FILE *fp, int n, double *buf, int flag)
{
for (int k = 0; k < nstyles; k++) {
if (flag == ELLIPSOID && strcmp(keywords[k],"ellipsoid") == 0)
styles[k]->write_data_bonus(fp,n,buf,flag);
if (flag == LINE && strcmp(keywords[k],"line") == 0)
styles[k]->write_data_bonus(fp,n,buf,flag);
if (flag == TRIANGLE && strcmp(keywords[k],"tri") == 0)
if (flag == ELLIPSOID && strcmp(keywords[k],"ellipsoid") != 0) continue;
if (flag == LINE && strcmp(keywords[k],"line") != 0) continue;
if (flag == TRIANGLE && strcmp(keywords[k],"tri") != 0) continue;
if (flag == BODY && strcmp(keywords[k],"body") != 0) continue;
styles[k]->write_data_bonus(fp,n,buf,flag);
}
}

View File

@ -58,9 +58,8 @@ class AtomVecHybrid : public AtomVec {
void pack_data_pre(int);
void pack_data_post(int);
int size_data_bonus_query(int);
int pack_data_bonus(double **, int);
void write_data_bonus(FILE *, int, double **, int);
int pack_data_bonus(double *, int);
void write_data_bonus(FILE *, int, double *, int);
int property_atom(char *);
void pack_property_atom(int, double *, int, int);

View File

@ -457,10 +457,10 @@ void AtomVecLine::pack_data_post(int ilocal)
/* ----------------------------------------------------------------------
pack bonus line info for writing to data file
if buf is NULL, just return count of lines
if buf is NULL, just return buffer size
------------------------------------------------------------------------- */
int AtomVecLine::pack_data_bonus(double **buf, int /*flag*/)
int AtomVecLine::pack_data_bonus(double *buf, int /*flag*/)
{
int i,j;
double length,theta;
@ -474,7 +474,7 @@ int AtomVecLine::pack_data_bonus(double **buf, int /*flag*/)
for (i = 0; i < nlocal; i++) {
if (line[i] < 0) continue;
if (buf) {
buf[m][0] = ubuf(tag[i]).d;
buf[m++] = ubuf(tag[i]).d;
j = line[i];
length = bonus[j].length;
theta = bonus[j].theta;
@ -484,12 +484,11 @@ int AtomVecLine::pack_data_bonus(double **buf, int /*flag*/)
y1 = yc - 0.5*sin(theta)*length;
x2 = xc + 0.5*cos(theta)*length;
y2 = yc + 0.5*sin(theta)*length;
buf[m][1] = x1;
buf[m][2] = y1;
buf[m][3] = x2;
buf[m][4] = y2;
}
m++;
buf[m++] = x1;
buf[m++] = y1;
buf[m++] = x2;
buf[m++] = y2;
} else m += size_data_bonus;
}
return m;
@ -499,12 +498,13 @@ int AtomVecLine::pack_data_bonus(double **buf, int /*flag*/)
write bonus line info to data file
------------------------------------------------------------------------- */
void AtomVecLine::write_data_bonus(FILE *fp, int n, double **buf, int /*flag*/)
void AtomVecLine::write_data_bonus(FILE *fp, int n, double *buf, int /*flag*/)
{
for (int i = 0; i < n; i++) {
fmt::print(fp,"{} {} {} {} {}",
(tagint) ubuf(buf[i][0]).i,
buf[i][1],buf[i][2],buf[i][3],buf[i][4]);
int i = 0;
while (i < n) {
fmt::print(fp,"{} {} {} {} {}\n",
(tagint) ubuf(buf[i]).i,buf[i+1],buf[i+2],buf[i+3],buf[i+4]);
i += size_data_bonus;
}
}

View File

@ -56,8 +56,8 @@ class AtomVecLine : public AtomVec {
void pack_data_pre(int);
void pack_data_post(int);
int pack_data_bonus(double **, int);
void write_data_bonus(FILE *, int, double **, int);
int pack_data_bonus(double *, int);
void write_data_bonus(FILE *, int, double *, int);
// unique to AtomVecLine

View File

@ -688,10 +688,10 @@ void AtomVecTri::pack_data_post(int ilocal)
/* ----------------------------------------------------------------------
pack bonus tri info for writing to data file
if buf is NULL, just return count of lines
if buf is NULL, just return buffer size
------------------------------------------------------------------------- */
int AtomVecTri::pack_data_bonus(double **buf, int /*flag*/)
int AtomVecTri::pack_data_bonus(double *buf, int /*flag*/)
{
int i,j;
double xc,yc,zc;
@ -706,7 +706,7 @@ int AtomVecTri::pack_data_bonus(double **buf, int /*flag*/)
for (i = 0; i < nlocal; i++) {
if (tri[i] < 0) continue;
if (buf) {
buf[m][0] = ubuf(tag[i]).d;
buf[m++] = ubuf(tag[i]).d;
j = tri[i];
MathExtra::quat_to_mat(bonus[j].quat,p);
MathExtra::matvec(p,bonus[j].c1,dc1);
@ -715,15 +715,15 @@ int AtomVecTri::pack_data_bonus(double **buf, int /*flag*/)
xc = x[i][0];
yc = x[i][1];
zc = x[i][2];
buf[m][1] = xc + dc1[0];
buf[m][2] = yc + dc1[1];
buf[m][3] = zc + dc1[2];
buf[m][4] = xc + dc2[0];
buf[m][5] = yc + dc2[1];
buf[m][6] = zc + dc2[2];
buf[m][7] = xc + dc3[0];
buf[m][8] = yc + dc3[1];
buf[m][9] = zc + dc3[2];
buf[m++] = xc + dc1[0];
buf[m++] = yc + dc1[1];
buf[m++] = zc + dc1[2];
buf[m++] = xc + dc2[0];
buf[m++] = yc + dc2[1];
buf[m++] = zc + dc2[2];
buf[m++] = xc + dc3[0];
buf[m++] = yc + dc3[1];
buf[m++] = zc + dc3[2];
}
m++;
}
@ -735,14 +735,14 @@ int AtomVecTri::pack_data_bonus(double **buf, int /*flag*/)
write bonus tri info to data file
------------------------------------------------------------------------- */
void AtomVecTri::write_data_bonus(FILE *fp, int n, double **buf, int /*flag*/)
void AtomVecTri::write_data_bonus(FILE *fp, int n, double *buf, int /*flag*/)
{
for (int i = 0; i < n; i++) {
fmt::print(fp,"{} {} {} {} {} {} {} {} {} {} {}",
(tagint) ubuf(buf[i][0]).i,
buf[i][1],buf[i][2],buf[i][3],
buf[i][4],buf[i][5],buf[i][6],
buf[i][7],buf[i][8],buf[i][9]);
int i = 0;
while (i < n) {
fmt::print(fp,"{} {} {} {} {} {} {} {} {} {} {}\n",
(tagint) ubuf(buf[i]).i,buf[i+1],buf[i+2],buf[i+3],
buf[i+4],buf[i+5],buf[i+6],buf[i+7],buf[i+8],buf[i+9]);
i += size_data_bonus;
}
}

View File

@ -58,8 +58,8 @@ class AtomVecTri : public AtomVec {
void pack_data_pre(int);
void pack_data_post(int);
int pack_data_bonus(double **, int);
void write_data_bonus(FILE *, int, double **, int);
int pack_data_bonus(double *, int);
void write_data_bonus(FILE *, int, double *, int);
// unique to AtomVecTri

View File

@ -43,6 +43,9 @@ class Body : protected Pointers {
double *) {return 0;}
virtual void data_body(int, int, int, int*, double *) = 0;
virtual int pack_data_body(tagint, int, double *) = 0;
virtual int write_data_body(FILE *, double *) = 0;
virtual int noutrow(int) = 0;
virtual int noutcol() = 0;
virtual void output(int, int, double *) = 0;

View File

@ -658,16 +658,15 @@ void WriteData::impropers()
void WriteData::bonus(int flag)
{
// communication buffer for all my Bonus info
// maxrow X ncol = largest buffer needed by any proc
// maxvalues = largest buffer needed by any proc
int ncol = atom->avec->size_data_bonus_query(flag);
int sendrow = atom->avec->pack_data_bonus(NULL,flag);
int maxrow;
MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);
int nvalues = atom->avec->pack_data_bonus(NULL,flag);
int maxvalues;
MPI_Allreduce(&nvalues,&maxvalues,1,MPI_INT,MPI_MAX,world);
double **buf;
if (me == 0) memory->create(buf,MAX(1,maxrow),ncol,"write_data:buf");
else memory->create(buf,MAX(1,sendrow),ncol,"write_data:buf");
double *buf;
if (me == 0) memory->create(buf,MAX(1,maxvalues),"write_data:buf");
else memory->create(buf,MAX(1,nvalues),"write_data:buf");
// pack my bonus data into buf
@ -690,19 +689,18 @@ void WriteData::bonus(int flag)
for (int iproc = 0; iproc < nprocs; iproc++) {
if (iproc) {
MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_DOUBLE,iproc,0,world,&request);
MPI_Irecv(buf,maxvalues,MPI_DOUBLE,iproc,0,world,&request);
MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_DOUBLE,&recvrow);
recvrow /= ncol;
} else recvrow = sendrow;
MPI_Get_count(&status,MPI_DOUBLE,&nvalues);
}
atom->avec->write_data_bonus(fp,recvrow,buf,flag);
atom->avec->write_data_bonus(fp,nvalues,buf,flag);
}
} else {
MPI_Recv(&tmp,0,MPI_INT,0,0,world,MPI_STATUS_IGNORE);
MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_DOUBLE,0,0,world);
MPI_Rsend(buf,nvalues,MPI_DOUBLE,0,0,world);
}
memory->destroy(buf);