From 24bb33b94d4defc274860cdb65431d401ff63ef2 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Thu, 14 Feb 2013 23:03:03 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9499
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/RIGID/fix_rigid_small.cpp | 96 ++++++++++++++++-------------------
src/RIGID/fix_rigid_small.h | 10 ++--
2 files changed, 51 insertions(+), 55 deletions(-)
diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index 53568571d9..908ea6859e 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -1307,7 +1307,7 @@ void FixRigidSmall::create_bodies()
// when done, have full bbox for every rigid body my atoms are part of
frsptr = this;
- comm->ring(m,sizeof(double),buf,1,bbox_update,NULL);
+ comm->ring(m,sizeof(double),buf,1,ring_bbox,NULL);
// ctr = center pt of each rigid body my atoms are part of
@@ -1347,7 +1347,7 @@ void FixRigidSmall::create_bodies()
// when done, have idclose for every rigid body my atoms are part of
frsptr = this;
- comm->ring(m,sizeof(double),buf,2,nearest_update,NULL);
+ comm->ring(m,sizeof(double),buf,2,ring_nearest,NULL);
// set bodytag of all owned atoms, based on idclose
@@ -1373,7 +1373,7 @@ void FixRigidSmall::create_bodies()
update bounding box for rigid bodies my atoms are part of
------------------------------------------------------------------------- */
-void FixRigidSmall::bbox_update(int n, char *cbuf)
+void FixRigidSmall::ring_bbox(int n, char *cbuf)
{
std::map<int,int> *hash = frsptr->hash;
double **bbox = frsptr->bbox;
@@ -1405,7 +1405,7 @@ void FixRigidSmall::bbox_update(int n, char *cbuf)
update nearest atom to body center for rigid bodies my atoms are part of
------------------------------------------------------------------------- */
-void FixRigidSmall::nearest_update(int n, char *cbuf)
+void FixRigidSmall::ring_nearest(int n, char *cbuf)
{
std::map<int,int> *hash = frsptr->hash;
double **ctr = frsptr->ctr;
@@ -1593,12 +1593,13 @@ void FixRigidSmall::setup_bodies()
// dx,dy,dz = coords relative to center-of-mass
// symmetric 3x3 inertia tensor stored in Voigt notation as 6-vector
+ memory->create(itensor,nlocal_body+nghost_body,6,"rigid/small:itensor");
+ for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++)
+ for (i = 0; i < 6; i++) itensor[ibody][i] = 0.0;
+
double dx,dy,dz,rad;
double *inertia;
- for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++)
- for (i = 0; i < 6; i++) body[ibody].itensor[i] = 0.0;
-
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
Body *b = &body[atom2body[i]];
@@ -1612,7 +1613,7 @@ void FixRigidSmall::setup_bodies()
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
- inertia = b->itensor;
+ inertia = itensor[atom2body[i]];
inertia[0] += massone * (dy*dy + dz*dz);
inertia[1] += massone * (dx*dx + dz*dz);
inertia[2] += massone * (dx*dx + dy*dy);
@@ -1630,7 +1631,7 @@ void FixRigidSmall::setup_bodies()
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
- inertia = body[atom2body[i]].itensor;
+ inertia = itensor[atom2body[i]];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
@@ -1673,10 +1674,6 @@ void FixRigidSmall::setup_bodies()
}
}
- for (ibody = 0; ibody < nlocal_body; ibody++) {
- inertia = body[atom2body[i]].itensor;
- }
-
// reverse communicate inertia tensor of all bodies
commflag = ITENSOR;
@@ -1689,17 +1686,15 @@ void FixRigidSmall::setup_bodies()
int ierror;
double cross[3];
double tensor[3][3],evectors[3][3];
- double *itensor,*ex,*ey,*ez;
+ double *ex,*ey,*ez;
for (ibody = 0; ibody < nlocal_body; ibody++) {
- itensor = body[ibody].itensor;
-
- tensor[0][0] = itensor[0];
- tensor[1][1] = itensor[1];
- tensor[2][2] = itensor[2];
- tensor[1][2] = tensor[2][1] = itensor[3];
- tensor[0][2] = tensor[2][0] = itensor[4];
- tensor[0][1] = tensor[1][0] = itensor[5];
+ tensor[0][0] = itensor[ibody][0];
+ tensor[1][1] = itensor[ibody][1];
+ tensor[2][2] = itensor[ibody][2];
+ tensor[1][2] = tensor[2][1] = itensor[ibody][3];
+ tensor[0][2] = tensor[2][0] = itensor[ibody][4];
+ tensor[0][1] = tensor[1][0] = itensor[ibody][5];
inertia = body[ibody].inertia;
ierror = MathExtra::jacobi(tensor,inertia,evectors);
@@ -1808,11 +1803,11 @@ void FixRigidSmall::setup_bodies()
// extended particles may contribute extra terms to moments of inertia
for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++)
- for (i = 0; i < 6; i++) body[ibody].itensor[i] = 0.0;
+ for (i = 0; i < 6; i++) itensor[ibody][i] = 0.0;
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
- inertia = body[atom2body[i]].itensor;
+ inertia = itensor[atom2body[i]];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
@@ -1835,7 +1830,7 @@ void FixRigidSmall::setup_bodies()
for (i = 0; i < nlocal; i++) {
if (atom2body[i] < 0) continue;
- inertia = body[atom2body[i]].itensor;
+ inertia = itensor[atom2body[i]];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
@@ -1887,33 +1882,32 @@ void FixRigidSmall::setup_bodies()
double norm;
for (ibody = 0; ibody < nlocal_body; ibody++) {
inertia = body[ibody].inertia;
- itensor = body[ibody].itensor;
if (inertia[0] == 0.0) {
- if (fabs(itensor[0]) > TOLERANCE)
+ if (fabs(itensor[ibody][0]) > TOLERANCE)
error->all(FLERR,"Fix rigid: Bad principal moments");
} else {
- if (fabs((itensor[0]-inertia[0])/inertia[0]) >
+ if (fabs((itensor[ibody][0]-inertia[0])/inertia[0]) >
TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
}
if (inertia[1] == 0.0) {
- if (fabs(itensor[1]) > TOLERANCE)
+ if (fabs(itensor[ibody][1]) > TOLERANCE)
error->all(FLERR,"Fix rigid: Bad principal moments");
} else {
- if (fabs((itensor[1]-inertia[1])/inertia[1]) >
+ if (fabs((itensor[ibody][1]-inertia[1])/inertia[1]) >
TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
}
if (inertia[2] == 0.0) {
- if (fabs(itensor[2]) > TOLERANCE)
+ if (fabs(itensor[ibody][2]) > TOLERANCE)
error->all(FLERR,"Fix rigid: Bad principal moments");
} else {
- if (fabs((itensor[2]-inertia[2])/inertia[2]) >
+ if (fabs((itensor[ibody][2]-inertia[2])/inertia[2]) >
TOLERANCE) error->all(FLERR,"Fix rigid: Bad principal moments");
}
norm = (inertia[0] + inertia[1] + inertia[2]) / 3.0;
- if (fabs(itensor[3]/norm) > TOLERANCE ||
- fabs(itensor[4]/norm) > TOLERANCE ||
- fabs(itensor[5]/norm) > TOLERANCE)
+ if (fabs(itensor[ibody][3]/norm) > TOLERANCE ||
+ fabs(itensor[ibody][4]/norm) > TOLERANCE ||
+ fabs(itensor[ibody][5]/norm) > TOLERANCE)
error->all(FLERR,"Fix rigid: Bad principal moments");
}
}
@@ -2265,7 +2259,7 @@ void FixRigidSmall::unpack_comm(int n, int first, double *buf)
int FixRigidSmall::pack_reverse_comm(int n, int first, double *buf)
{
int i,j,m,last;
- double *fcm,*torque,*vcm,*angmom,*xcm,*itensor;
+ double *fcm,*torque,*vcm,*angmom,*xcm;
m = 0;
last = first + n;
@@ -2309,13 +2303,13 @@ int FixRigidSmall::pack_reverse_comm(int n, int first, double *buf)
} else if (commflag == ITENSOR) {
for (i = first; i < last; i++) {
if (bodyown[i] < 0) continue;
- itensor = body[bodyown[i]].itensor;
- buf[m++] = itensor[0];
- buf[m++] = itensor[1];
- buf[m++] = itensor[2];
- buf[m++] = itensor[3];
- buf[m++] = itensor[4];
- buf[m++] = itensor[5];
+ j = bodyown[i];
+ buf[m++] = itensor[j][0];
+ buf[m++] = itensor[j][1];
+ buf[m++] = itensor[j][2];
+ buf[m++] = itensor[j][3];
+ buf[m++] = itensor[j][4];
+ buf[m++] = itensor[j][5];
}
} else if (commflag == DOF) {
@@ -2338,7 +2332,7 @@ int FixRigidSmall::pack_reverse_comm(int n, int first, double *buf)
void FixRigidSmall::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,k;
- double *fcm,*torque,*vcm,*angmom,*xcm,*itensor;
+ double *fcm,*torque,*vcm,*angmom,*xcm;
int *tag = atom->tag;
int nlocal = atom->nlocal;
@@ -2388,13 +2382,13 @@ void FixRigidSmall::unpack_reverse_comm(int n, int *list, double *buf)
for (i = 0; i < n; i++) {
j = list[i];
if (bodyown[j] < 0) continue;
- itensor = body[bodyown[j]].itensor;
- itensor[0] += buf[m++];
- itensor[1] += buf[m++];
- itensor[2] += buf[m++];
- itensor[3] += buf[m++];
- itensor[4] += buf[m++];
- itensor[5] += buf[m++];
+ k = bodyown[j];
+ itensor[k][0] += buf[m++];
+ itensor[k][1] += buf[m++];
+ itensor[k][2] += buf[m++];
+ itensor[k][3] += buf[m++];
+ itensor[k][4] += buf[m++];
+ itensor[k][5] += buf[m++];
}
} else if (commflag == DOF) {
diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h
index 8e84144a62..237b3f173b 100644
--- a/src/RIGID/fix_rigid_small.h
+++ b/src/RIGID/fix_rigid_small.h
@@ -80,7 +80,6 @@ class FixRigidSmall : public Fix {
double torque[3]; // torque on COM
double quat[4]; // quaternion for orientation of body
double inertia[3]; // 3 principal components of inertia
- double itensor[6]; // 6 space-frame components of inertia tensor
double ex_space[3]; // principal axes in space coords
double ey_space[3];
double ez_space[3];
@@ -122,7 +121,10 @@ class FixRigidSmall : public Fix {
class AtomVecLine *avec_line;
class AtomVecTri *avec_tri;
- int **counts;
+ // temporary per-body storage
+
+ int **counts; // counts of atom types in bodies
+ double **itensor; // 6 space-frame components of inertia tensor
// Langevin thermostatting
@@ -149,8 +151,8 @@ class FixRigidSmall : public Fix {
// callback functions for ring communication
- static void bbox_update(int, char *);
- static void nearest_update(int, char *);
+ static void ring_bbox(int, char *);
+ static void ring_nearest(int, char *);
// debug
--
GitLab