tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

clean whitespace

Browse Source
This commit is contained in:
Steve Plimpton
2024-07-10 16:13:24 -06:00
parent 04490971e4
commit 0812e77d27
12 changed files with 211 additions and 211 deletions
Download Patch File Download Diff File Expand all files Collapse all files

56
src/GRANSURF/fix_surface_global.cpp
View File

@ -1231,7 +1231,7 @@ int FixSurfaceGlobal::endpt_neigh_check(int i, int j, int jflag)
ncheck = connect2d[j].np2;
neighs = connect2d[j].neigh_p2;
}
// check overlap with each neighbor line
// if any line has interior overlap, another line computes
// if all lines have endpt overlap, line with lowest index computes
@ -1498,7 +1498,7 @@ int FixSurfaceGlobal::edge_neigh_check(int i, int j, int jflag)
int ncheck;
int *neighs;
if (jflag == -1) {
if (connect3d[j].ne1 == 1) return 0;
ncheck = connect3d[j].ne1;
@ -1568,7 +1568,7 @@ int FixSurfaceGlobal::corner_neigh_check(int i, int j, int jflag)
ncheck = connect3d[j].nc3;
neighs = connect3d[j].neigh_c3;
}
// check overlap with each neighbor tri
// if any tri has interior or edge overlap, another tri computes
// if all tris have corner pt overlap, tri with lowest index computes
@ -1581,11 +1581,11 @@ int FixSurfaceGlobal::corner_neigh_check(int i, int j, int jflag)
int k,kflag;
double rsq;
double dr[3],contact[3];
int trimin = j;
for (int m = 0; m < ncheck; m++) {
k = neighs[m];
k = neighs[m];
if (k == j) continue; // skip self tri
kflag = overlap_sphere_tri(i,k,contact,dr,rsq);
if (kflag > 0) return 1;
@ -1620,13 +1620,13 @@ void FixSurfaceGlobal::extract_from_molecules(char *molID)
tris = nullptr;
npoints = nlines = ntris = 0;
int maxpoints = 0;
int imol = atom->find_molecule(molID);
if (imol == -1)
error->all(FLERR,"Molecule template ID for fix surface/global does not exist");
// loop over one or more molecules in molID
Molecule **onemols = &atom->molecules[imol];
int nmol = onemols[0]->nset;
@ -1651,7 +1651,7 @@ void FixSurfaceGlobal::extract_from_molecules(char *molID)
// create a map
// key = xyz coords of a point
// value = index into unique points vector
std::map<std::tuple<double,double,double>,int> hash;
// offset line/tri index lists by previous npoints
@ -1770,7 +1770,7 @@ void FixSurfaceGlobal::connectivity2d_global()
{
connect2d = (Connect2d *) memory->smalloc(nlines*sizeof(Connect2d),
"surface/global:connect2d");
// setup end point connectivity lists
// count # of lines containing each end point
// create ragged 2d array to contain all point indices, then fill it
@ -1779,7 +1779,7 @@ void FixSurfaceGlobal::connectivity2d_global()
int *counts;
memory->create(counts,npoints,"surface/global:count");
for (int i = 0; i < npoints; i++) counts[i] = 0;
for (int i = 0; i < nlines; i++) {
@ -1800,12 +1800,12 @@ void FixSurfaceGlobal::connectivity2d_global()
connect2d[i].np1 = counts[lines[i].p1];
if (connect2d[i].np1 == 1) connect2d[i].neigh_p1 = nullptr;
else connect2d[i].neigh_p1 = plist[lines[i].p1];
connect2d[i].np2 = counts[lines[i].p2];
if (connect2d[i].np2 == 1) connect2d[i].neigh_p2 = nullptr;
else connect2d[i].neigh_p2 = plist[lines[i].p2];
}
memory->destroy(counts);
}
@ -1816,27 +1816,27 @@ void FixSurfaceGlobal::connectivity2d_global()
void FixSurfaceGlobal::connectivity3d_global()
{
int p1,p2,p3;
connect3d = (Connect3d *) memory->smalloc(ntris*sizeof(Connect3d),
"surface/global:connect3d");
int **tri2edge;
memory->create(tri2edge,ntris,3,"surfface/global::tri2edge");
// create a map
// key = <p1,p2> indices of 2 points, in either order
// value = index into count of unique edges
std::map<std::tuple<int,int>,int> hash;
int nedges = 0;
for (int i = 0; i < ntris; i++) {
p1 = tris[i].p1;
p2 = tris[i].p2;
p3 = tris[i].p3;
auto key1 = std::make_tuple(p1,p2);
auto key2 = std::make_tuple(p2,p1);
if (hash.find(key1) == hash.end() && hash.find(key2) == hash.end()) {
hash[key1] = nedges;
tri2edge[i][0] = nedges;
@ -1847,7 +1847,7 @@ void FixSurfaceGlobal::connectivity3d_global()
key1 = std::make_tuple(p2,p3);
key2 = std::make_tuple(p3,p2);
if (hash.find(key1) == hash.end() && hash.find(key2) == hash.end()) {
hash[key1] = nedges;
tri2edge[i][1] = nedges;
@ -1858,7 +1858,7 @@ void FixSurfaceGlobal::connectivity3d_global()
key1 = std::make_tuple(p3,p1);
key2 = std::make_tuple(p1,p3);
if (hash.find(key1) == hash.end() && hash.find(key2) == hash.end()) {
hash[key1] = nedges;
tri2edge[i][2] = nedges;
@ -1888,7 +1888,7 @@ void FixSurfaceGlobal::connectivity3d_global()
memory->create_ragged(elist,nedges,counts,"surface/global:elist");
for (int i = 0; i < nedges; i++) counts[i] = 0;
for (int i = 0; i < ntris; i++) {
elist[tri2edge[i][0]][counts[tri2edge[i][0]]++] = i;
elist[tri2edge[i][1]][counts[tri2edge[i][1]]++] = i;
@ -1899,11 +1899,11 @@ void FixSurfaceGlobal::connectivity3d_global()
connect3d[i].ne1 = counts[tri2edge[i][0]];
if (connect3d[i].ne1 == 1) connect3d[i].neigh_e1 = nullptr;
else connect3d[i].neigh_e1 = elist[tri2edge[i][0]];
connect3d[i].ne2 = counts[tri2edge[i][1]];
if (connect3d[i].ne2 == 1) connect3d[i].neigh_e2 = nullptr;
else connect3d[i].neigh_e2 = elist[tri2edge[i][1]];
connect3d[i].ne3 = counts[tri2edge[i][2]];
if (connect3d[i].ne3 == 1) connect3d[i].neigh_e3 = nullptr;
else connect3d[i].neigh_e3 = elist[tri2edge[i][2]];
@ -1911,7 +1911,7 @@ void FixSurfaceGlobal::connectivity3d_global()
memory->destroy(counts);
memory->destroy(tri2edge);
// setup corner point connectivity lists
// count # of tris containing each point
// create ragged 2d array to contain all tri indices, then fill it
@ -1919,7 +1919,7 @@ void FixSurfaceGlobal::connectivity3d_global()
// nc123 counts and vectors include self tri
memory->create(counts,npoints,"surface/global:count");
for (int i = 0; i < npoints; i++) counts[i] = 0;
for (int i = 0; i < ntris; i++) {
@ -1942,11 +1942,11 @@ void FixSurfaceGlobal::connectivity3d_global()
connect3d[i].nc1 = counts[tris[i].p1];
if (connect3d[i].nc1 == 1) connect3d[i].neigh_c1 = nullptr;
else connect3d[i].neigh_c1 = clist[tris[i].p1];
connect3d[i].nc2 = counts[tris[i].p2];
if (connect3d[i].nc2 == 1) connect3d[i].neigh_c2 = nullptr;
else connect3d[i].neigh_c2 = clist[tris[i].p2];
connect3d[i].nc3 = counts[tris[i].p3];
if (connect3d[i].nc3 == 1) connect3d[i].neigh_c3 = nullptr;
else connect3d[i].neigh_c3 = clist[tris[i].p3];

4
src/GRANSURF/fix_surface_global.h
View File

@ -27,7 +27,7 @@ namespace LAMMPS_NS {
class FixSurfaceGlobal : public Fix {
public:
// neighbor lists for spheres with surfs and shear history
// accessed by fix shear/history
@ -104,7 +104,7 @@ class FixSurfaceGlobal : public Fix {
Tri *tris; // global list of tris
int npoints,nlines,ntris; // count of each
int nsurf; // count of lines or tris for 2d/3d
int **plist; // ragged 2d array for global end pt lists
int **elist; // ragged 2d array for global edge lists
int **clist; // ragged 2d array for global corner pt lists

154
src/GRANSURF/fix_surface_local.cpp
View File

@ -65,7 +65,7 @@ FixSurfaceLocal::FixSurfaceLocal(LAMMPS *lmp, int narg, char **arg) :
nlocal_connect = nghost_connect = nmax_connect = 0;
connect2d = nullptr;
connect3d = nullptr;
// NOTE: what is optimal way to initialize MyPoolChunk?
tcp = nullptr;
@ -233,18 +233,18 @@ void FixSurfaceLocal::pre_neighbor()
grow connect array by DELTA
------------------------------------------------------------------------- */
void FixSurfaceLocal::grow_connect()
void FixSurfaceLocal::grow_connect()
{
nmax_connect = nmax_connect/DELTA_BONUS * DELTA_BONUS;
nmax_connect += DELTA_BONUS;
if (nmax_connect < 0) error->one(FLERR,"Per-processor system is too big");
if (dimension == 2)
connect2d = (Connect2d *)
if (dimension == 2)
connect2d = (Connect2d *)
memory->srealloc(connect2d,nmax_connect*sizeof(Connect2d),
"surface/local:connect2d");
else
connect3d = (Connect3d *)
connect3d = (Connect3d *)
memory->srealloc(connect3d,nmax_connect*sizeof(Connect3d),
"surface/local:connect3d");
}
@ -256,7 +256,7 @@ void FixSurfaceLocal::grow_connect()
void FixSurfaceLocal::copy_arrays(int i, int j, int delflag)
{
// if deleting atom J via delflag and J has connect data, then delete it
if (dimension == 2) {
if (delflag && index[j] >= 0) {
int k = index[j];
@ -275,7 +275,7 @@ void FixSurfaceLocal::copy_arrays(int i, int j, int delflag)
}
// if atom I has connect data, reset I's connect.ilocal to loc J
// do NOT do this if self-copy (I=J) since I's connection data
// do NOT do this if self-copy (I=J) since I's connection data
// is already deleted
if (index[i] >= 0 && i != j) {
@ -367,9 +367,9 @@ int FixSurfaceLocal::pack_border(int n, int *list, double *buf)
buf[m++] = ubuf(connect3d[ic].nc2).d;
buf[m++] = ubuf(connect3d[ic].nc3).d;
buf[m++] = ubuf(connect3d[ic].flags).d;
// NOTE: which pack method to use
nc = connect3d[ic].nc1;
for (k = 0; k < nc; k++)
buf[m++] = ubuf(connect3d[ic].neigh_c1[k]).d;
@ -408,7 +408,7 @@ int FixSurfaceLocal::unpack_border(int n, int first, double *buf)
connect2d[j].neigh_p1 = (tagint) ubuf(buf[m++]).i;
connect2d[j].neigh_p2 = (tagint) ubuf(buf[m++]).i;
connect2d[j].flags = (int) ubuf(buf[m++]).i;
connect2d[j].ilocal = i;
index[i] = j;
nghost_connect++;
@ -423,7 +423,7 @@ int FixSurfaceLocal::unpack_border(int n, int first, double *buf)
else {
j = nlocal_connect + nghost_connect;
if (j == nmax_connect) grow_connect();
connect3d[j].neigh_e1 = (tagint) ubuf(buf[m++]).i;
connect3d[j].neigh_e2 = (tagint) ubuf(buf[m++]).i;
connect3d[j].neigh_e3 = (tagint) ubuf(buf[m++]).i;
@ -431,13 +431,13 @@ int FixSurfaceLocal::unpack_border(int n, int first, double *buf)
connect3d[j].nc2 = (int) ubuf(buf[m++]).i;
connect3d[j].nc3 = (int) ubuf(buf[m++]).i;
connect3d[j].flags = (int) ubuf(buf[m++]).i;
connect3d[j].neigh_c1 = tcp->get(connect3d[j].nc1,connect3d[j].indexc1);
connect3d[j].neigh_c2 = tcp->get(connect3d[j].nc2,connect3d[j].indexc2);
connect3d[j].neigh_c3 = tcp->get(connect3d[j].nc3,connect3d[j].indexc3);
// NOTE: which unpack method to use
nc = connect3d[j].nc1;
for (k = 0; k < nc; k++)
connect3d[j].neigh_c1[k] = (tagint) ubuf(buf[m++]).i;
@ -447,7 +447,7 @@ int FixSurfaceLocal::unpack_border(int n, int first, double *buf)
nc = connect3d[j].nc3;
for (k = 0; k < nc; k++)
connect3d[j].neigh_c3[k] = (tagint) ubuf(buf[m++]).i;
connect3d[j].ilocal = i;
index[i] = j;
nghost_connect++;
@ -491,7 +491,7 @@ int FixSurfaceLocal::pack_exchange(int i, double *buf)
buf[m++] = ubuf(connect3d[j].nc2).d;
buf[m++] = ubuf(connect3d[j].nc3).d;
buf[m++] = ubuf(connect3d[j].flags).d;
// NOTE: which method to use
n = connect3d[j].nc1;
@ -503,7 +503,7 @@ int FixSurfaceLocal::pack_exchange(int i, double *buf)
n = connect3d[j].nc3;
for (k = 0; k < n; k++)
buf[m++] = ubuf(connect3d[j].neigh_c3[k]).d;
/*
memcpy(&buf[m],connect3d[i].neigh_c1,connect3d[i].nc1*sizeof(tagint));
if (tagintdoubleratio == 1) m += connect3d[i].nc1;
@ -559,13 +559,13 @@ int FixSurfaceLocal::unpack_exchange(int nlocal, double *buf)
connect3d[nlocal_connect].flags = (int) ubuf(buf[m++]).i;
connect3d[nlocal_connect].ilocal = nlocal;
connect3d[nlocal_connect].neigh_c1 =
connect3d[nlocal_connect].neigh_c1 =
tcp->get(connect3d[nlocal_connect].nc1,
connect3d[nlocal_connect].indexc1);
connect3d[nlocal_connect].neigh_c2 =
connect3d[nlocal_connect].neigh_c2 =
tcp->get(connect3d[nlocal_connect].nc2,
connect3d[nlocal_connect].indexc2);
connect3d[nlocal_connect].neigh_c3 =
connect3d[nlocal_connect].neigh_c3 =
tcp->get(connect3d[nlocal_connect].nc3,
connect3d[nlocal_connect].indexc3);
@ -633,7 +633,7 @@ void FixSurfaceLocal::connectivity2d_local()
// error check
avec_line = (AtomVecLine *) atom->style_match("line");
if (!avec_line)
if (!avec_line)
error->all(FLERR,"Fix surface/local requires atom style line");
// calculate current endpts of owned lines
@ -792,14 +792,14 @@ void FixSurfaceLocal::connectivity2d_local()
fptr = this;
nmatch = errormatch = 0;
/*
if (ebuf)
comm->ring(nline,5*sizeof(double),(void *) &ebuf[0][0],0,linematch,nullptr);
else
else
comm->ring(nline,5*sizeof(double),nullptr,0,linematch,nullptr);
*/
// check for errors = matches with more than 2 points
// print stats on # of matches
@ -817,7 +817,7 @@ void FixSurfaceLocal::connectivity2d_local()
if (comm->me == 0) {
if (screen)
fprintf(screen," matched %g line end point pairs\n",0.5*all);
if (logfile)
if (logfile)
fprintf(logfile," matched %g line end point pairs\n",0.5*all);
}
@ -957,7 +957,7 @@ void FixSurfaceLocal::calculate_endpts(int n)
int FixSurfaceLocal::pt2bin2d(double *pt)
{
if (pt[0] < binlo[0] || pt[0] >= binhi[0] ||
if (pt[0] < binlo[0] || pt[0] >= binhi[0] ||
pt[1] < binlo[1] || pt[1] >= binhi[1])
return -1;
@ -995,7 +995,7 @@ int FixSurfaceLocal::overlap2bin2d(double *pt, double eps, int *indices)
int i,j;
int n = 0;
for (j = jlo; j <= jhi; j++)
for (j = jlo; j <= jhi; j++)
for (i = ilo; i <= ihi; i++)
indices[n++] = j*nbinx + i;
return n;
@ -1018,7 +1018,7 @@ void FixSurfaceLocal::connectivity3d_local()
// error check
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri)
if (!avec_tri)
error->all(FLERR,"Fix surface/local requires atom style tri");
// calculate current corner pts of owned triangles
@ -1197,25 +1197,25 @@ void FixSurfaceLocal::connectivity3d_local()
fptr = this;
nmatch1 = nmatch2 = errormatch1 = errormatch2 = 0;
vecflag = 0;
/*
if (tbuf)
if (tbuf)
comm->ring(ntri,10*sizeof(double),(void *) &tbuf[0][0],0,trimatch,nullptr);
else
else
comm->ring(ntri,10*sizeof(double),nullptr,0,trimatch,nullptr);
*/
for (i = 0; i < nlocal; i++) {
if (tri[i] < 0) continue;
j = tri[i];
connect3d[j].indexc1 = connect3d[j].indexc2 = connect3d[j].indexc3 = -1;
if (connect3d[j].nc1)
if (connect3d[j].nc1)
connect3d[j].neigh_c1 = tcp->get(connect3d[j].nc1,connect3d[j].indexc1);
else connect3d[j].neigh_c1 = nullptr;
if (connect3d[j].nc2)
if (connect3d[j].nc2)
connect3d[j].neigh_c2 = tcp->get(connect3d[j].nc2,connect3d[j].indexc2);
else connect3d[j].neigh_c2 = nullptr;
if (connect3d[j].nc3)
if (connect3d[j].nc3)
connect3d[j].neigh_c3 = tcp->get(connect3d[j].nc3,connect3d[j].indexc3);
else connect3d[j].neigh_c3 = nullptr;
connect3d[j].neigh_e1 = connect3d[j].neigh_e2 = connect3d[j].neigh_e3 = 0;
@ -1227,14 +1227,14 @@ void FixSurfaceLocal::connectivity3d_local()
fptr = this;
nmatch1 = nmatch2 = errormatch1 = errormatch2 = 0;
vecflag = 1;
/*
if (tbuf)
if (tbuf)
comm->ring(ntri,10*sizeof(double),(void *) &tbuf[0][0],0,trimatch,nullptr);
else
else
comm->ring(ntri,10*sizeof(double),nullptr,0,trimatch,nullptr);
*/
// check for errors = matches with misaligned or more than 2 edges
// print stats on # of matches
@ -1258,20 +1258,20 @@ void FixSurfaceLocal::connectivity3d_local()
MPI_Allreduce(&nmatch1,&all,1,MPI_INT,MPI_SUM,world);
if (comm->me == 0) {
if (screen)
if (screen)
fprintf(screen,
" matched %g tri edge pairs in fix surface/local\n",0.5*all);
if (logfile)
if (logfile)
fprintf(logfile,
" matched %g tri edge pairs in fix surface/local\n",0.5*all);
}
MPI_Allreduce(&nmatch2,&all,1,MPI_INT,MPI_SUM,world);
if (comm->me == 0) {
if (screen)
if (screen)
fprintf(screen,
" matched %g tri corners in fix surface/local\n",0.5*all);
if (logfile)
if (logfile)
fprintf(logfile,
" matched %g tri corners in fix surface/local\n",0.5*all);
}
@ -1355,7 +1355,7 @@ void FixSurfaceLocal::trimatch(int n, char *cbuf)
// loop over all tri corner pts in bins[jfirst:jlast], owned by me
// skip bin pt if from same tri ID as buf pt is from
for (j = jfirst; j < jlast; j++) {
iatom = pts[j].iatom;
iconnect = pts[j].iconnect;
@ -1419,9 +1419,9 @@ void FixSurfaceLocal::trimatch(int n, char *cbuf)
else {
connect3d[iconnect].neigh_e1 = tagother;
if (vecflag) {
connect3d[iconnect].neigh_c1[connect3d[iconnect].nc1] =
connect3d[iconnect].neigh_c1[connect3d[iconnect].nc1] =
tagother;
connect3d[iconnect].neigh_c2[connect3d[iconnect].nc2] =
connect3d[iconnect].neigh_c2[connect3d[iconnect].nc2] =
tagother;
}
connect3d[iconnect].nc1++;
@ -1432,9 +1432,9 @@ void FixSurfaceLocal::trimatch(int n, char *cbuf)
else {
connect3d[iconnect].neigh_e2 = tagother;
if (vecflag) {
connect3d[iconnect].neigh_c2[connect3d[iconnect].nc2] =
connect3d[iconnect].neigh_c2[connect3d[iconnect].nc2] =
tagother;
connect3d[iconnect].neigh_c3[connect3d[iconnect].nc3] =
connect3d[iconnect].neigh_c3[connect3d[iconnect].nc3] =
tagother;
}
connect3d[iconnect].nc2++;
@ -1445,9 +1445,9 @@ void FixSurfaceLocal::trimatch(int n, char *cbuf)
else {
connect3d[iconnect].neigh_e3 = tagother;
if (vecflag) {
connect3d[iconnect].neigh_c3[connect3d[iconnect].nc3] =
connect3d[iconnect].neigh_c3[connect3d[iconnect].nc3] =
tagother;
connect3d[iconnect].neigh_c1[connect3d[iconnect].nc1] =
connect3d[iconnect].neigh_c1[connect3d[iconnect].nc1] =
tagother;
}
connect3d[iconnect].nc3++;
@ -1466,19 +1466,19 @@ void FixSurfaceLocal::trimatch(int n, char *cbuf)
(*nmatch2)++;
if (ptwhich == 1) {
if (vecflag) {
connect3d[iconnect].neigh_c1[connect3d[iconnect].nc1] =
connect3d[iconnect].neigh_c1[connect3d[iconnect].nc1] =
tagother;
}
connect3d[iconnect].nc1++;
} else if (ptwhich == 2) {
if (vecflag) {
connect3d[iconnect].neigh_c2[connect3d[iconnect].nc2] =
connect3d[iconnect].neigh_c2[connect3d[iconnect].nc2] =
tagother;
}
connect3d[iconnect].nc2++;
} else if (ptwhich == 3) {
if (vecflag) {
connect3d[iconnect].neigh_c3[connect3d[iconnect].nc3] =
connect3d[iconnect].neigh_c3[connect3d[iconnect].nc3] =
tagother;
}
connect3d[iconnect].nc3++;
@ -1530,7 +1530,7 @@ void FixSurfaceLocal::calculate_corners(int n)
int FixSurfaceLocal::pt2bin3d(double *pt)
{
if (pt[0] < binlo[0] || pt[0] >= binhi[0] ||
if (pt[0] < binlo[0] || pt[0] >= binhi[0] ||
pt[1] < binlo[1] || pt[1] >= binhi[1] ||
pt[2] < binlo[2] || pt[2] >= binhi[2])
return -1;
@ -1560,7 +1560,7 @@ int FixSurfaceLocal::overlap2bin3d(double *pt, double eps, int *indices)
int jhi = static_cast<int> ((pt[1]+eps-binlo[1]) * invbiny);
int klo = static_cast<int> ((pt[2]-eps-binlo[2]) * invbinz);
int khi = static_cast<int> ((pt[2]+eps-binlo[2]) * invbinz);
ilo = MAX(ilo,0);
ihi = MIN(ihi,nbinx-1);
jlo = MAX(jlo,0);
@ -1575,10 +1575,10 @@ int FixSurfaceLocal::overlap2bin3d(double *pt, double eps, int *indices)
int i,j,k;
int n = 0;
for (k = klo; k <= khi; k++)
for (j = jlo; j <= jhi; j++)
for (k = klo; k <= khi; k++)
for (j = jlo; j <= jhi; j++)
for (i = ilo; i <= ihi; i++)
indices[n++] = k*nbiny*nbinx + j*nbinx + i;
indices[n++] = k*nbiny*nbinx + j*nbinx + i;
return n;
}
@ -1607,7 +1607,7 @@ void FixSurfaceLocal::extract_from_molecules(char *molID)
if (surf[i] >= 0) flag = 1;
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
if (flagall)
if (flagall)
error->all(FLERR,"Fix surface/local molecule file when surfaces already exist");
// populate global point/line/tri data structs
@ -1617,16 +1617,16 @@ void FixSurfaceLocal::extract_from_molecules(char *molID)
tris = nullptr;
npoints = nlines = ntris = 0;
int maxpoints = 0;
int imol = atom->find_molecule(molID);
if (imol == -1)
error->all(FLERR,"Molecule template ID for fix surface/local does not exist");
// loop over one or more molecules in molID
Molecule **onemols = &atom->molecules[imol];
int nmol = onemols[0]->nset;
for (int m = 0; m < nmol; m++) {
if (dimension == 2)
if (onemols[m]->lineflag == 0)
@ -1648,7 +1648,7 @@ void FixSurfaceLocal::extract_from_molecules(char *molID)
// create a map
// key = xyz coords of a point
// value = index into unique points vector
std::map<std::tuple<double,double,double>,int> hash;
// offset line/tri index lists by previous npoints
@ -1800,7 +1800,7 @@ void FixSurfaceLocal::connectivity2d_global()
p1 = lines[i].p1;
p2 = lines[i].p2;
if (ptflag[p1][0] < 0)
if (ptflag[p1][0] < 0)
error->all(FLERR,"Fix surface/local point part of more than 2 lines");
else if (ptflag[p1][0] == 1)
error->all(FLERR,"Fix surface/local pair of lines are misoriented");
@ -1814,7 +1814,7 @@ void FixSurfaceLocal::connectivity2d_global()
connect2dall[j].neigh_p2 = i+1;
}
if (ptflag[p2][0] < 0)
if (ptflag[p2][0] < 0)
error->all(FLERR,"Fix surface/local point part of more than 2 lines");
else if (ptflag[p2][0] == 2)
error->all(FLERR,"Fix surface/local pair of lines are misoriented");
@ -1851,7 +1851,7 @@ void FixSurfaceLocal::connectivity3d_global()
"surface/local:connect3dall");
for (int i = 0; i < ntris; i++) {
connect3dall[i].neigh_e1 = connect3dall[i].neigh_e2 =
connect3dall[i].neigh_e1 = connect3dall[i].neigh_e2 =
connect3dall[i].neigh_e3 = 0;
connect3dall[i].flags = 0;
}
@ -1976,7 +1976,7 @@ void FixSurfaceLocal::assign2d()
AtomVec *avec = atom->avec;
avec_line = (AtomVecLine *) atom->style_match("line");
if (!avec_line)
if (!avec_line)
error->all(FLERR,"Fix surface/local requires atom style line");
// set bounds for my proc
@ -2055,7 +2055,7 @@ void FixSurfaceLocal::assign2d()
char *values[4];
for (i = 0; i < 4; i++) values[i] = &pts[i][0];
std::vector<std::string> svalues(5);
for (i = 0; i < nlines; i++) {
// midpt of line
@ -2084,7 +2084,7 @@ void FixSurfaceLocal::assign2d()
// create a default particle of correct style
avec->create_atom(lines[i].type,xmid);
// change it to be a line
int n = atom->nlocal - 1;
@ -2104,12 +2104,12 @@ void FixSurfaceLocal::assign2d()
// svalues[0] is not used here
// used when caller of data_atom_bonus() is via read_data command
svalues[1] = (const char *) values[0];
svalues[2] = (const char *) values[1];
svalues[3] = (const char *) values[2];
svalues[4] = (const char *) values[3];
avec_line->data_atom_bonus(n,svalues);
// copy global connectivity enty for line I to local list
@ -2146,7 +2146,7 @@ void FixSurfaceLocal::assign3d()
AtomVec *avec = atom->avec;
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri)
if (!avec_tri)
error->all(FLERR,"Fix surface/local requires atom style tri");
// set bounds for my proc
@ -2258,7 +2258,7 @@ void FixSurfaceLocal::assign3d()
// create a default particle of correct style
avec->create_atom(tris[i].type,xmid);
// change it to be a triangle
int n = atom->nlocal - 1;
@ -2314,7 +2314,7 @@ void FixSurfaceLocal::assign3d()
self = connect3d[nlocal_connect].neigh_c1;
for (j = 0; j < nc; j++) self[j] += idmaxall;
if (nc) {
neigh = connect3d[nlocal_connect].neigh_c1 =
neigh = connect3d[nlocal_connect].neigh_c1 =
tcp->get(nc,connect3d[nlocal_connect].indexc1);
m = 0;
for (j = 0; j <= nc; j++)
@ -2325,7 +2325,7 @@ void FixSurfaceLocal::assign3d()
self = connect3d[nlocal_connect].neigh_c2;
for (j = 0; j < nc; j++) self[j] += idmaxall;
if (nc) {
neigh = connect3d[nlocal_connect].neigh_c2 =
neigh = connect3d[nlocal_connect].neigh_c2 =
tcp->get(nc,connect3d[nlocal_connect].indexc2);
m = 0;
for (j = 0; j <= nc; j++)
@ -2336,7 +2336,7 @@ void FixSurfaceLocal::assign3d()
self = connect3d[nlocal_connect].neigh_c3;
for (j = 0; j < nc; j++) self[j] += idmaxall;
if (nc) {
neigh = connect3d[nlocal_connect].neigh_c3 =
neigh = connect3d[nlocal_connect].neigh_c3 =
tcp->get(nc,connect3d[nlocal_connect].indexc3);
m = 0;
for (j = 0; j <= nc; j++)

6
src/GRANSURF/fix_surface_local.h
View File

@ -57,7 +57,7 @@ class FixSurfaceLocal : public Fix {
// size of local/ghost connection info vectors
int nlocal_connect,nghost_connect,nmax_connect;
int nlocal_connect,nghost_connect,nmax_connect;
FixSurfaceLocal(class LAMMPS *, int, char **);
virtual ~FixSurfaceLocal();
@ -137,13 +137,13 @@ class FixSurfaceLocal : public Fix {
double epssq; // distance tolerance for end pts
// from different lines to be connected
int nmatch; // # of line connections
int nmatch1,nmatch2; // # of tri connections
int errormatch; // # of errors with line connectivity
int errormatch1,errormatch2; // # of errors with tri connectivity
int vecflag; // 0/1 whether tri matching should also
int vecflag; // 0/1 whether tri matching should also
// store variable-length vecs of corner connections
// static variable for ring communication callback to access class data

26
src/GRANSURF/pair_line_gran_hooke.cpp
View File

@ -25,7 +25,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairLineGranHooke::PairLineGranHooke(LAMMPS *lmp) :
PairLineGranHooke::PairLineGranHooke(LAMMPS *lmp) :
PairLineGranHookeHistory(lmp)
{
history = 0;
@ -131,10 +131,10 @@ void PairLineGranHooke::compute(int eflag, int vflag)
if (line[i] >= 0 || line[j] < 0)
error->one(FLERR,"Pair line/gran iteraction is invalid");
// check for overlap of sphere and line segment
// check for overlap of sphere and line segment
// jflag = 0 for no overlap, 1 for interior line pt, -1/-2 for end pts
// if no overlap, just continue
// for overlap, also return:
// for overlap, also return:
// contact = nearest point on line to sphere center
// dr = vector from contact pt to sphere center
// rsq = squared length of dr
@ -169,13 +169,13 @@ void PairLineGranHooke::compute(int eflag, int vflag)
ds[1] = contact[1] - x[j][1];
ds[2] = contact[2] - x[j][2];
// vline = velocity of contact pt on line, translation + rotation
// vline = velocity of contact pt on line, translation + rotation
vline[0] = v[j][0] + (omega[j][1]*ds[2] - omega[j][2]*ds[1]);
vline[1] = v[j][1] + (omega[j][2]*ds[0] - omega[j][0]*ds[2]);
vline[2] = v[j][2] + (omega[j][0]*ds[1] - omega[j][1]*ds[0]);
// relative translational velocity
// relative translational velocity
vr1 = v[i][0] - vline[0];
vr2 = v[i][1] - vline[1];
@ -194,7 +194,7 @@ void PairLineGranHooke::compute(int eflag, int vflag)
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// relative rotational velocity
// relative rotational velocity
wr1 = (radi*omega[i][0]) * rinv;
wr2 = (radi*omega[i][1]) * rinv;
@ -240,11 +240,11 @@ void PairLineGranHooke::compute(int eflag, int vflag)
fs2 = -ft*vtr2 * factor_couple;
fs3 = -ft*vtr3 * factor_couple;
// total force on sphere
// total force on sphere
fx = dr[0]*ccel + fs1;
fy = dr[1]*ccel + fs2;
fz = dr[2]*ccel + fs3;
fx = dr[0]*ccel + fs1;
fy = dr[1]*ccel + fs2;
fz = dr[2]*ccel + fs3;
// sphere force & torque
@ -260,9 +260,9 @@ void PairLineGranHooke::compute(int eflag, int vflag)
torque[i][2] -= radi*tor3;
// line force & torque
// torque applied at contact pt
// use total force for torque
// since opposite force is not norm/tang to line at its end pt
// torque applied at contact pt
// use total force for torque
// since opposite force is not norm/tang to line at its end pt
f[j][0] -= fx;
f[j][1] -= fy;

48
src/GRANSURF/pair_line_gran_hooke_history.cpp
View File

@ -38,7 +38,7 @@ PairLineGranHookeHistory::PairLineGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
single_enable = 0;
history = 1;
size_history = 3;
nmax = 0;
mass_rigid = nullptr;
@ -114,7 +114,7 @@ void PairLineGranHookeHistory::compute(int eflag, int vflag)
// mass_body = mass of each rigid body
// forward comm mass_rigid so have it for ghost lines
// also grab current line connectivity info from FixSurfaceLocal
if (neighbor->ago == 0) {
if (fix_rigid) {
int tmp;
@ -128,9 +128,9 @@ void PairLineGranHookeHistory::compute(int eflag, int vflag)
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (body[i] >= 0)
mass_rigid[i] = mass_body[body[i]];
mass_rigid[i] = mass_body[body[i]];
else
mass_rigid[i] = 0.0;
mass_rigid[i] = 0.0;
comm->forward_comm(this);
}
@ -213,10 +213,10 @@ void PairLineGranHookeHistory::compute(int eflag, int vflag)
if (line[i] >= 0 || line[j] < 0)
error->one(FLERR,"Pair line/gran iteraction is invalid");
// check for overlap of sphere and line segment
// check for overlap of sphere and line segment
// jflag = 0 for no overlap, 1 for interior line pt, -1/-2 for end pts
// if no overlap, just continue
// for overlap, also return:
// for overlap, also return:
// contact = nearest point on line to sphere center
// dr = vector from contact pt to sphere center
// rsq = squared length of dr
@ -257,13 +257,13 @@ void PairLineGranHookeHistory::compute(int eflag, int vflag)
ds[1] = contact[1] - x[j][1];
ds[2] = contact[2] - x[j][2];
// vline = velocity of contact pt on line, translation + rotation
// vline = velocity of contact pt on line, translation + rotation
vline[0] = v[j][0] + (omega[j][1]*ds[2] - omega[j][2]*ds[1]);
vline[1] = v[j][1] + (omega[j][2]*ds[0] - omega[j][0]*ds[2]);
vline[2] = v[j][2] + (omega[j][0]*ds[1] - omega[j][1]*ds[0]);
// relative translational velocity
// relative translational velocity
vr1 = v[i][0] - vline[0];
vr2 = v[i][1] - vline[1];
@ -282,7 +282,7 @@ void PairLineGranHookeHistory::compute(int eflag, int vflag)
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// relative rotational velocity
// relative rotational velocity
wr1 = (radi*omega[i][0]) * rinv;
wr2 = (radi*omega[i][1]) * rinv;
@ -363,11 +363,11 @@ void PairLineGranHookeHistory::compute(int eflag, int vflag)
} else fs1 = fs2 = fs3 = 0.0;
}
// total force on sphere
// total force on sphere
fx = dr[0]*ccel + fs1;
fy = dr[1]*ccel + fs2;
fz = dr[2]*ccel + fs3;
fx = dr[0]*ccel + fs1;
fy = dr[1]*ccel + fs2;
fz = dr[2]*ccel + fs3;
// sphere force & torque
@ -383,9 +383,9 @@ void PairLineGranHookeHistory::compute(int eflag, int vflag)
torque[i][2] -= radi*tor3;
// line force & torque
// torque applied at contact pt
// use total force for torque
// since opposite force is not norm/tang to line at its end pt
// torque applied at contact pt
// use total force for torque
// since opposite force is not norm/tang to line at its end pt
f[j][0] -= fx;
f[j][1] -= fy;
@ -544,7 +544,7 @@ void PairLineGranHookeHistory::init_style()
fsl = nullptr;
for (int m = 0; m < modify->nfix; m++) {
if (strcmp(modify->fix[m]->style,"surface/local") == 0) {
if (fsl)
if (fsl)
error->all(FLERR,"Pair line/gran requires single fix surface/local");
fsl = (FixSurfaceLocal *) modify->fix[m];
}
@ -564,7 +564,7 @@ void PairLineGranHookeHistory::init_style()
int nlocal = atom->nlocal;
int flag = 0;
for (int i = 0; i < nlocal; i++) {
if (line[i] < 0) continue;
if (line[i] < 0) continue;
if (mask[i] & groupbit) flag = 1;
}
int any;
@ -626,7 +626,7 @@ void PairLineGranHookeHistory::init_style()
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (line[i] >= 0)
if (line[i] >= 0)
onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]);
else {
if (mask[i] & freeze_group_bit)
@ -666,7 +666,7 @@ double PairLineGranHookeHistory::init_one(int i, int j)
cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]);
return cutoff;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
@ -761,7 +761,7 @@ int PairLineGranHookeHistory::pack_forward_comm(int n, int *list, double *buf,
/* ---------------------------------------------------------------------- */
void PairLineGranHookeHistory::unpack_forward_comm(int n, int first,
void PairLineGranHookeHistory::unpack_forward_comm(int n, int first,
double *buf)
{
int i,m,last;
@ -840,7 +840,7 @@ void PairLineGranHookeHistory::calculate_endpts()
based on geometry.cpp::distsq_point_line() in SPARTA
------------------------------------------------------------------------- */
int PairLineGranHookeHistory::overlap_sphere_line(int i, int j, double *pt,
int PairLineGranHookeHistory::overlap_sphere_line(int i, int j, double *pt,
double *r, double &rsq)
{
double p1[3],p2[3];
@ -867,7 +867,7 @@ int PairLineGranHookeHistory::overlap_sphere_line(int i, int j, double *pt,
// alpha = fraction of distance from P1 to P2 that P is located at
// P = projected point on infinite line that is nearest to Xsphere center
// alpha can be any value
double alpha = MathExtra::dot3(a,b) / MathExtra::lensq3(b);
// pt = point on line segment that is nearest to Xsphere center
@ -939,7 +939,7 @@ int PairLineGranHookeHistory::endpt_neigh_check(int i, int j, int jflag)
int kflag = overlap_sphere_line(i,k,contact,dr,rsq);
if (kflag > 0) return 1;
if (kflag == 0)
if (kflag == 0)
error->one(FLERR,"Pair line/gran neighbor line overlap is invalid");
if (tag[j] < tag[k]) return 0;

32
src/GRANSURF/pair_tri_gran_hooke.cpp
View File

@ -130,11 +130,11 @@ void PairTriGranHooke::compute(int eflag, int vflag)
if (tri[i] >= 0 || tri[j] < 0)
error->one(FLERR,"Pair tri/gran iteraction is invalid");
// check for overlap of sphere and triangle
// check for overlap of sphere and triangle
// jflag = 0 for no overlap, 1 for interior line pt,
// -1/-2/-3 for 3 edges, -4/-5/-6 for 3 corner pts
// -1/-2/-3 for 3 edges, -4/-5/-6 for 3 corner pts
// if no overlap, just continue
// for overlap, also returns:
// for overlap, also returns:
// contact = nearest point on tri to sphere center
// dr = vector from contact pt to sphere center
// rsq = squared length of dr
@ -170,14 +170,14 @@ void PairTriGranHooke::compute(int eflag, int vflag)
ds[1] = contact[1] - x[j][1];
ds[2] = contact[2] - x[j][2];
// vtri = velocity of contact pt on tri, translation + rotation
// omega for tri was set from angmom by calculate_corners()
// vtri = velocity of contact pt on tri, translation + rotation
// omega for tri was set from angmom by calculate_corners()
vtri[0] = v[j][0] + (omega[j][1]*ds[2] - omega[j][2]*ds[1]);
vtri[1] = v[j][1] + (omega[j][2]*ds[0] - omega[j][0]*ds[2]);
vtri[2] = v[j][2] + (omega[j][0]*ds[1] - omega[j][1]*ds[0]);
// relative translational velocity
// relative translational velocity
vr1 = v[i][0] - vtri[0];
vr2 = v[i][1] - vtri[1];
@ -196,7 +196,7 @@ void PairTriGranHooke::compute(int eflag, int vflag)
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// relative rotational velocity
// relative rotational velocity
wr1 = (radi*omega[i][0]) * rinv;
wr2 = (radi*omega[i][1]) * rinv;
@ -242,18 +242,18 @@ void PairTriGranHooke::compute(int eflag, int vflag)
fs2 = -ft*vtr2 * factor_couple;
fs3 = -ft*vtr3 * factor_couple;
// total force on sphere
// total force on sphere
fx = dr[0]*ccel + fs1;
fy = dr[1]*ccel + fs2;
fz = dr[2]*ccel + fs3;
fx = dr[0]*ccel + fs1;
fy = dr[1]*ccel + fs2;
fz = dr[2]*ccel + fs3;
// sphere force & torque
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
tor1 = rinv * (dr[1]*fs3 - dr[2]*fs2);
tor2 = rinv * (dr[2]*fs1 - dr[0]*fs3);
tor3 = rinv * (dr[0]*fs2 - dr[1]*fs1);
@ -262,14 +262,14 @@ void PairTriGranHooke::compute(int eflag, int vflag)
torque[i][2] -= radi*tor3;
// tri force & torque
// torque applied at contact pt
// use total force for torque
// since opposite force is not norm/tang to tri at its contact pt
// torque applied at contact pt
// use total force for torque
// since opposite force is not norm/tang to tri at its contact pt
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
torque[j][0] -= ds[1]*fz - ds[2]*fy;
torque[j][1] -= ds[2]*fx - ds[0]*fz;
torque[j][2] -= ds[0]*fy - ds[1]*fx;

54
src/GRANSURF/pair_tri_gran_hooke_history.cpp
View File

@ -128,9 +128,9 @@ void PairTriGranHookeHistory::compute(int eflag, int vflag)
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (body[i] >= 0)
mass_rigid[i] = mass_body[body[i]];
mass_rigid[i] = mass_body[body[i]];
else
mass_rigid[i] = 0.0;
mass_rigid[i] = 0.0;
comm->forward_comm(this);
}
@ -205,11 +205,11 @@ void PairTriGranHookeHistory::compute(int eflag, int vflag)
if (tri[i] >= 0 || tri[j] < 0)
error->one(FLERR,"Pair tri/gran iteraction is invalid");
// check for overlap of sphere and triangle
// check for overlap of sphere and triangle
// jflag = 0 for no overlap, 1 for interior line pt,
// -1/-2/-3 for 3 edges, -4/-5/-6 for 3 corner pts
// -1/-2/-3 for 3 edges, -4/-5/-6 for 3 corner pts
// if no overlap, just continue
// for overlap, also returns:
// for overlap, also returns:
// contact = nearest point on tri to sphere center
// dr = vector from contact pt to sphere center
// rsq = squared length of dr
@ -255,14 +255,14 @@ void PairTriGranHookeHistory::compute(int eflag, int vflag)
ds[1] = contact[1] - x[j][1];
ds[2] = contact[2] - x[j][2];
// vtri = velocity of contact pt on tri, translation + rotation
// omega for tri was set from angmom by calculate_corners()
// vtri = velocity of contact pt on tri, translation + rotation
// omega for tri was set from angmom by calculate_corners()
vtri[0] = v[j][0] + (omega[j][1]*ds[2] - omega[j][2]*ds[1]);
vtri[1] = v[j][1] + (omega[j][2]*ds[0] - omega[j][0]*ds[2]);
vtri[2] = v[j][2] + (omega[j][0]*ds[1] - omega[j][1]*ds[0]);
// relative translational velocity
// relative translational velocity
vr1 = v[i][0] - vtri[0];
vr2 = v[i][1] - vtri[1];
@ -281,7 +281,7 @@ void PairTriGranHookeHistory::compute(int eflag, int vflag)
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// relative rotational velocity
// relative rotational velocity
wr1 = (radi*omega[i][0]) * rinv;
wr2 = (radi*omega[i][1]) * rinv;
@ -362,18 +362,18 @@ void PairTriGranHookeHistory::compute(int eflag, int vflag)
} else fs1 = fs2 = fs3 = 0.0;
}
// total force on sphere
// total force on sphere
fx = dr[0]*ccel + fs1;
fy = dr[1]*ccel + fs2;
fz = dr[2]*ccel + fs3;
fx = dr[0]*ccel + fs1;
fy = dr[1]*ccel + fs2;
fz = dr[2]*ccel + fs3;
// sphere force & torque
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
tor1 = rinv * (dr[1]*fs3 - dr[2]*fs2);
tor2 = rinv * (dr[2]*fs1 - dr[0]*fs3);
tor3 = rinv * (dr[0]*fs2 - dr[1]*fs1);
@ -382,14 +382,14 @@ void PairTriGranHookeHistory::compute(int eflag, int vflag)
torque[i][2] -= radi*tor3;
// tri force & torque
// torque applied at contact pt
// use total force for torque
// since opposite force is not norm/tang to tri at its contact pt
// torque applied at contact pt
// use total force for torque
// since opposite force is not norm/tang to tri at its contact pt
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
torque[j][0] -= ds[1]*fz - ds[2]*fy;
torque[j][1] -= ds[2]*fx - ds[0]*fz;
torque[j][2] -= ds[0]*fy - ds[1]*fx;
@ -527,13 +527,13 @@ void PairTriGranHookeHistory::init_style()
modify->replace_fix("NEIGH_HISTORY_HHTRI_DUMMY" + std::to_string(instance_me), cmd, 1));
fix_history->pair = this;
}
// set ptr to FixSurfaceLocal for surf connectivity info
fsl = nullptr;
for (int m = 0; m < modify->nfix; m++) {
if (strcmp(modify->fix[m]->style,"surface/local") == 0) {
if (fsl)
if (fsl)
error->all(FLERR,"Pair tri/gran requires single fix surface/local");
fsl = (FixSurfaceLocal *) modify->fix[m];
}
@ -616,7 +616,7 @@ void PairTriGranHookeHistory::init_style()
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (tri[i] >= 0)
if (tri[i] >= 0)
onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]);
else {
if (mask[i] & freeze_group_bit)
@ -656,7 +656,7 @@ double PairTriGranHookeHistory::init_one(int i, int j)
cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]);
return cutoff;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
@ -857,7 +857,7 @@ void PairTriGranHookeHistory::corners2norm(double *corners, double *norm)
based on geometry.cpp::distsq_point_tri() in SPARTA
------------------------------------------------------------------------- */
int PairTriGranHookeHistory::overlap_sphere_tri(int i, int j, double *pt,
int PairTriGranHookeHistory::overlap_sphere_tri(int i, int j, double *pt,
double *r, double &rsq)
{
int e12flag,e23flag,e31flag,o12flag,o23flag,o31flag;
@ -1024,8 +1024,8 @@ int PairTriGranHookeHistory::overlap_sphere_tri(int i, int j, double *pt,
based on geometry.cpp::distsq_point_line() in SPARTA
------------------------------------------------------------------------- */
int PairTriGranHookeHistory::nearest_point_line(double *x,
double *p1, double *p2,
int PairTriGranHookeHistory::nearest_point_line(double *x,
double *p1, double *p2,
double *pt)
{
double a[3],b[3];
@ -1039,7 +1039,7 @@ int PairTriGranHookeHistory::nearest_point_line(double *x,
// alpha = fraction of distance from P1 to P2 that P is located at
// P = projected point on infinite line that is nearest to X
// alpha can be any value
double alpha = MathExtra::dot3(a,b) / MathExtra::lensq3(b);
// pt = point on line segment that is nearest to X
@ -1099,7 +1099,7 @@ int PairTriGranHookeHistory::edge_neigh_check(int i, int j, int jflag)
int kflag = overlap_sphere_tri(i,k,contact,dr,rsq);
if (kflag > 0) return 1;
if (kflag == 0)
if (kflag == 0)
error->one(FLERR,"Pair tri/gran neighbor tri overlap is invalid");
tagint *tag = atom->tag;

2
src/GRANSURF/pair_tri_gran_hooke_history.h
View File

@ -50,7 +50,7 @@ class PairTriGranHookeHistory : public Pair {
double dt;
int freeze_group_bit;
int history;
int neighprev;
double *onerad_dynamic,*onerad_frozen;
double *maxrad_dynamic,*maxrad_frozen;

8
src/molecule.cpp
View File

@ -448,7 +448,7 @@ void Molecule::read(int flag)
if (text.empty()) continue;
// search line for header keywords and set corresponding variable
try {
ValueTokenizer values(text);
@ -546,7 +546,7 @@ void Molecule::read(int flag)
if (ndihedrals < 0) error->all(FLERR, "Invalid dihedral count in molecule file");
if (nimpropers < 0) error->all(FLERR, "Invalid improper count in molecule file");
if (natoms == 0 && nlines == 0 && ntris == 0)
if (natoms == 0 && nlines == 0 && ntris == 0)
error->all(FLERR,"Molecule file must define either atoms or lines or triangles");
if (nlines && domain->dimension != 2)
@ -627,7 +627,7 @@ void Molecule::read(int flag)
if (flag)
line_segments(line);
else
skip_lines(nlines, line, keyword);
skip_lines(nlines, line, keyword);
} else if (keyword == "Tris") {
triflag = 1;
if (flag)
@ -1109,7 +1109,7 @@ void Molecule::line_segments(char *line)
} catch (TokenizerException &e) {
error->all(FLERR, "Invalid line in Lines section of molecule file: {}\n{}", e.what(), line);
}
// check all line molecule-IDs and types
// add toffset to line type

2
src/neighbor.cpp
View File

@ -2111,7 +2111,7 @@ int Neighbor::choose_pair(NeighRequest *rq)
for (int i = 0; i < npclass; i++) {
mask = pairmasks[i];
//printf(" PAIR NAMES i %d %d name %s mask %d\n",i,nrequest,
// pairnames[i],pairmasks[i]);

30
src/stl_reader.cpp
View File

@ -47,34 +47,34 @@ int STLReader::read_file(const char *filename, double **caller_tris)
if (me == 0) {
// open file (text or binary)
FILE *fp = fopen(filename, "rb");
if (fp == nullptr) error->one(FLERR, "Cannot open file {}: {}", filename, utils::getsyserror());
// first try reading the file in ASCII format
TextFileReader reader(fp, "STL mesh");
try {
char *line = reader.next_line();
if (!line || !utils::strmatch(line, "^solid"))
throw TokenizerException("Invalid STL mesh file format", "");
line += 6;
if (utils::strmatch(line, "^binary"))
throw TokenizerException("Invalid STL mesh file format", "");
utils::logmesg(lmp, "Reading STL object {} from text file {}\n", utils::trim(line), filename);
while ((line = reader.next_line())) {
// next line is facet line with 5 words
auto values = utils::split_words(line);
// otherwise next line should be endsolid and are done reading file
if ((values.size() != 5) || !utils::strmatch(values[0], "^facet")) {
if (!utils::strmatch(values[0], "^endsolid"))
throw TokenizerException("Error reading endsolid", "");
@ -93,18 +93,18 @@ int STLReader::read_file(const char *filename, double **caller_tris)
maxtris += DELTA;
memory->grow(tris,maxtris,9,"STLReader:tris");
}
for (int k = 0; k < 3; ++k) {
line = reader.next_line(4);
values = utils::split_words(line);
if ((values.size() != 4) || !utils::strmatch(values[0], "^vertex"))
throw TokenizerException("Error reading vertex", "");
tris[ntris][3*k+0] = utils::numeric(FLERR, values[1], false, lmp);
tris[ntris][3*k+1] = utils::numeric(FLERR, values[2], false, lmp);
tris[ntris][3*k+2] = utils::numeric(FLERR, values[3], false, lmp);
}
line = reader.next_line(1);
if (!line || !utils::strmatch(line, "^ *endloop"))
throw TokenizerException("Error reading endloop", "");
@ -114,11 +114,11 @@ int STLReader::read_file(const char *filename, double **caller_tris)
ntris++;
}
} catch (std::exception &e) {
// if read of text file failed for the first line, try reading as binary
if (utils::strmatch(e.what(), "^Invalid STL mesh file format")) {
char title[80];
float triangle[12];
@ -138,7 +138,7 @@ int STLReader::read_file(const char *filename, double **caller_tris)
}
// NOTE: worry about unsigned int versus signed int
memory->create(tris,ntris,9,"STLReader:tris");
for (uint32_t i = 0U; i < ntri; ++i) {
@ -152,7 +152,7 @@ int STLReader::read_file(const char *filename, double **caller_tris)
for (int k = 0; k < 3; ++k)
tris[i][m++] = triangle[3 * j + 3 + k];
}
} else {
error->all(FLERR, "Error reading triangles from file {}: {}", filename, e.what());
}