From fa7543b714e068a8a5bba7058343574fa0a3e26e Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 27 Jan 2016 20:33:55 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14506 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GRANULAR/fix_pour.cpp | 24 +++++++- src/atom_vec_tri.cpp | 112 ++++++++++++++++++++++++++++++-------- src/atom_vec_tri.h | 2 +- src/displace_atoms.cpp | 14 +++-- src/dump_image.cpp | 20 ++++++- src/dump_image.h | 3 +- src/fix_move.cpp | 25 ++++----- 7 files changed, 150 insertions(+), 50 deletions(-) diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 3e1ea6b89c..662d46667e 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -156,7 +156,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : memory->create(coords,natom_max,4,"pour:coords"); memory->create(imageflags,natom_max,"pour:imageflags"); - // find current max atom and molecule IDs if necessary + // find max atom and molecule IDs just once if (idnext) find_maxid(); @@ -385,7 +385,7 @@ void FixPour::pre_exchange() atom->nghost = 0; atom->avec->clear_bonus(); - // find current max atom and molecule IDs if necessary + // find current max atom and molecule IDs on every insertion step if (!idnext) find_maxid(); @@ -737,12 +737,19 @@ void FixPour::find_maxid() return 1 if yes, 0 if no for ATOM mode, use delta with maximum size for inserted atoms for MOLECULE mode, use delta with max radius of inserted molecules + if ignore line/tri set, ignore line or tri particles account for PBC in overlap decision via outside() and minimum_image() ------------------------------------------------------------------------- */ int FixPour::overlap(int i) { double delta; + + if (ignoreflag) { + if (ignoreline && atom->line[i] >= 0) return 0; + if (ignoretri && atom->tri[i] >= 0) return 0; + } + if (mode == ATOM) delta = atom->radius[i] + radius_max; else delta = atom->radius[i] + molradius_max; @@ -860,6 +867,7 @@ void FixPour::options(int narg, char **arg) shakeflag = 0; idshake = NULL; idnext = 0; + ignoreflag = ignoreline = ignoretri = 0; dstyle = ONE; radius_max = radius_one = 0.5; radius_poly = frac_poly = NULL; @@ -925,6 +933,13 @@ void FixPour::options(int narg, char **arg) else if (strcmp(arg[iarg+1],"next") == 0) idnext = 1; else error->all(FLERR,"Illegal fix pour command"); iarg += 2; + + } else if (strcmp(arg[iarg],"ignore") == 0) { + if (atom->line_flag) ignoreline = 1; + if (atom->tri_flag) ignoretri = 1; + if (ignoreline || ignoretri) ignoreflag = 1; + iarg += 1; + } else if (strcmp(arg[iarg],"diam") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command"); if (strcmp(arg[iarg+1],"one") == 0) { @@ -938,6 +953,7 @@ void FixPour::options(int narg, char **arg) dstyle = RANGE; radius_lo = 0.5 * force->numeric(FLERR,arg[iarg+2]); radius_hi = 0.5 * force->numeric(FLERR,arg[iarg+3]); + if (radius_lo > radius_hi) error->all(FLERR,"Illegal fix pour command"); radius_max = radius_hi; iarg += 4; } else if (strcmp(arg[iarg+1],"poly") == 0) { @@ -968,6 +984,7 @@ void FixPour::options(int narg, char **arg) if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command"); density_lo = force->numeric(FLERR,arg[iarg+1]); density_hi = force->numeric(FLERR,arg[iarg+2]); + if (density_lo > density_hi) error->all(FLERR,"Illegal fix pour command"); iarg += 3; } else if (strcmp(arg[iarg],"vol") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command"); @@ -985,6 +1002,8 @@ void FixPour::options(int narg, char **arg) vxhi = force->numeric(FLERR,arg[iarg+2]); vylo = force->numeric(FLERR,arg[iarg+3]); vyhi = force->numeric(FLERR,arg[iarg+4]); + if (vxlo > vxhi || vylo > vyhi) + error->all(FLERR,"Illegal fix pour command"); vz = force->numeric(FLERR,arg[iarg+5]); iarg += 6; } else { @@ -993,6 +1012,7 @@ void FixPour::options(int narg, char **arg) vxhi = force->numeric(FLERR,arg[iarg+2]); vy = force->numeric(FLERR,arg[iarg+3]); vz = 0.0; + if (vxlo > vxhi) error->all(FLERR,"Illegal fix pour command"); iarg += 4; } } else error->all(FLERR,"Illegal fix pour command"); diff --git a/src/atom_vec_tri.cpp b/src/atom_vec_tri.cpp index cbfb5af8e2..7dc65d5f0f 100644 --- a/src/atom_vec_tri.cpp +++ b/src/atom_vec_tri.cpp @@ -41,7 +41,7 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp) size_forward = 7; size_reverse = 6; size_border = 26; - size_velocity = 6; + size_velocity = 9; size_data_atom = 8; size_data_vel = 7; size_data_bonus = 10; @@ -49,7 +49,9 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp) atom->tri_flag = 1; atom->molecule_flag = atom->rmass_flag = 1; - atom->radius_flag = atom->angmom_flag = atom->torque_flag = 1; + atom->radius_flag = atom->omega_flag = atom->angmom_flag = 1; + atom->torque_flag = 1; + atom->sphere_flag = 1; nlocal_bonus = nghost_bonus = nmax_bonus = 0; bonus = NULL; @@ -100,6 +102,7 @@ void AtomVecTri::grow(int n) molecule = memory->grow(atom->molecule,nmax,"atom:molecule"); rmass = memory->grow(atom->rmass,nmax,"atom:rmass"); radius = memory->grow(atom->radius,nmax,"atom:radius"); + omega = memory->grow(atom->omega,nmax,3,"atom:omega"); angmom = memory->grow(atom->angmom,nmax,3,"atom:angmom"); torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque"); tri = memory->grow(atom->tri,nmax,"atom:tri"); @@ -119,7 +122,8 @@ void AtomVecTri::grow_reset() mask = atom->mask; image = atom->image; x = atom->x; v = atom->v; f = atom->f; molecule = atom->molecule; rmass = atom->rmass; - radius = atom->radius; angmom = atom->angmom; torque = atom->torque; + radius = atom->radius; omega = atom->omega; + angmom = atom->angmom; torque = atom->torque; tri = atom->tri; } @@ -158,6 +162,9 @@ void AtomVecTri::copy(int i, int j, int delflag) molecule[j] = molecule[i]; rmass[j] = rmass[i]; radius[j] = radius[i]; + omega[j][0] = omega[i][0]; + omega[j][1] = omega[i][1]; + omega[j][2] = omega[i][2]; angmom[j][0] = angmom[i][0]; angmom[j][1] = angmom[i][1]; angmom[j][2] = angmom[i][2]; @@ -343,6 +350,9 @@ int AtomVecTri::pack_comm_vel(int n, int *list, double *buf, buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; buf[m++] = angmom[j][0]; buf[m++] = angmom[j][1]; buf[m++] = angmom[j][2]; @@ -373,6 +383,9 @@ int AtomVecTri::pack_comm_vel(int n, int *list, double *buf, buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; buf[m++] = angmom[j][0]; buf[m++] = angmom[j][1]; buf[m++] = angmom[j][2]; @@ -402,6 +415,9 @@ int AtomVecTri::pack_comm_vel(int n, int *list, double *buf, buf[m++] = v[j][1]; buf[m++] = v[j][2]; } + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; buf[m++] = angmom[j][0]; buf[m++] = angmom[j][1]; buf[m++] = angmom[j][2]; @@ -479,6 +495,9 @@ void AtomVecTri::unpack_comm_vel(int n, int first, double *buf) v[i][0] = buf[m++]; v[i][1] = buf[m++]; v[i][2] = buf[m++]; + omega[i][0] = buf[m++]; + omega[i][1] = buf[m++]; + omega[i][2] = buf[m++]; angmom[i][0] = buf[m++]; angmom[i][1] = buf[m++]; angmom[i][2] = buf[m++]; @@ -728,6 +747,9 @@ int AtomVecTri::pack_border_vel(int n, int *list, double *buf, buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; buf[m++] = angmom[j][0]; buf[m++] = angmom[j][1]; buf[m++] = angmom[j][2]; @@ -782,6 +804,9 @@ int AtomVecTri::pack_border_vel(int n, int *list, double *buf, buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; buf[m++] = angmom[j][0]; buf[m++] = angmom[j][1]; buf[m++] = angmom[j][2]; @@ -835,6 +860,9 @@ int AtomVecTri::pack_border_vel(int n, int *list, double *buf, buf[m++] = v[j][1]; buf[m++] = v[j][2]; } + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; buf[m++] = angmom[j][0]; buf[m++] = angmom[j][1]; buf[m++] = angmom[j][2]; @@ -1002,6 +1030,9 @@ void AtomVecTri::unpack_border_vel(int n, int first, double *buf) v[i][0] = buf[m++]; v[i][1] = buf[m++]; v[i][2] = buf[m++]; + omega[i][0] = buf[m++]; + omega[i][1] = buf[m++]; + omega[i][2] = buf[m++]; angmom[i][0] = buf[m++]; angmom[i][1] = buf[m++]; angmom[i][2] = buf[m++]; @@ -1082,6 +1113,9 @@ int AtomVecTri::pack_exchange(int i, double *buf) buf[m++] = ubuf(molecule[i]).d; buf[m++] = rmass[i]; buf[m++] = radius[i]; + buf[m++] = omega[i][0]; + buf[m++] = omega[i][1]; + buf[m++] = omega[i][2]; buf[m++] = angmom[i][0]; buf[m++] = angmom[i][1]; buf[m++] = angmom[i][2]; @@ -1143,6 +1177,9 @@ int AtomVecTri::unpack_exchange(double *buf) molecule[nlocal] = (tagint) ubuf(buf[m++]).i; rmass[nlocal] = buf[m++]; radius[nlocal] = buf[m++]; + omega[nlocal][0] = buf[m++]; + omega[nlocal][1] = buf[m++]; + omega[nlocal][2] = buf[m++]; angmom[nlocal][0] = buf[m++]; angmom[nlocal][1] = buf[m++]; angmom[nlocal][2] = buf[m++]; @@ -1197,8 +1234,8 @@ int AtomVecTri::size_restart() int n = 0; int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) - if (tri[i] >= 0) n += 34; - else n += 18; + if (tri[i] >= 0) n += 37; + else n += 21; if (atom->nextra_restart) for (int iextra = 0; iextra < atom->nextra_restart; iextra++) @@ -1231,6 +1268,9 @@ int AtomVecTri::pack_restart(int i, double *buf) buf[m++] = ubuf(molecule[i]).d; buf[m++] = rmass[i]; buf[m++] = radius[i]; + buf[m++] = omega[i][0]; + buf[m++] = omega[i][1]; + buf[m++] = omega[i][2]; buf[m++] = angmom[i][0]; buf[m++] = angmom[i][1]; buf[m++] = angmom[i][2]; @@ -1298,6 +1338,9 @@ int AtomVecTri::unpack_restart(double *buf) molecule[nlocal] = (tagint) ubuf(buf[m++]).i; rmass[nlocal] = buf[m++]; radius[nlocal] = buf[m++]; + omega[nlocal][0] = buf[m++]; + omega[nlocal][1] = buf[m++]; + omega[nlocal][2] = buf[m++]; angmom[nlocal][0] = buf[m++]; angmom[nlocal][1] = buf[m++]; angmom[nlocal][2] = buf[m++]; @@ -1366,6 +1409,9 @@ void AtomVecTri::create_atom(int itype, double *coord) molecule[nlocal] = 0; radius[nlocal] = 0.5; rmass[nlocal] = 4.0*MY_PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal]; + omega[nlocal][0] = 0.0; + omega[nlocal][1] = 0.0; + omega[nlocal][2] = 0.0; angmom[nlocal][0] = 0.0; angmom[nlocal][1] = 0.0; angmom[nlocal][2] = 0.0; @@ -1415,6 +1461,9 @@ void AtomVecTri::data_atom(double *coord, imageint imagetmp, char **values) v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; v[nlocal][2] = 0.0; + omega[nlocal][0] = 0.0; + omega[nlocal][1] = 0.0; + omega[nlocal][2] = 0.0; angmom[nlocal][0] = 0.0; angmom[nlocal][1] = 0.0; angmom[nlocal][2] = 0.0; @@ -1591,21 +1640,27 @@ void AtomVecTri::data_vel(int m, char **values) v[m][0] = atof(values[0]); v[m][1] = atof(values[1]); v[m][2] = atof(values[2]); - angmom[m][0] = atof(values[3]); - angmom[m][1] = atof(values[4]); - angmom[m][2] = atof(values[5]); + omega[m][0] = atof(values[3]); + omega[m][1] = atof(values[4]); + omega[m][2] = atof(values[5]); + angmom[m][0] = atof(values[6]); + angmom[m][1] = atof(values[7]); + angmom[m][2] = atof(values[8]); } /* ---------------------------------------------------------------------- - unpack hybrid quantities from one tri in Velocities section of data file + unpack hybrid quantities from one line in Velocities section of data file ------------------------------------------------------------------------- */ int AtomVecTri::data_vel_hybrid(int m, char **values) { - angmom[m][0] = atof(values[0]); - angmom[m][1] = atof(values[1]); - angmom[m][2] = atof(values[2]); - return 3; + omega[m][0] = atof(values[0]); + omega[m][1] = atof(values[1]); + omega[m][2] = atof(values[2]); + angmom[m][0] = atof(values[3]); + angmom[m][1] = atof(values[4]); + angmom[m][2] = atof(values[5]); + return 6; } /* ---------------------------------------------------------------------- @@ -1703,9 +1758,12 @@ void AtomVecTri::pack_vel(double **buf) buf[i][1] = v[i][0]; buf[i][2] = v[i][1]; buf[i][3] = v[i][2]; - buf[i][4] = angmom[i][0]; - buf[i][5] = angmom[i][1]; - buf[i][6] = angmom[i][2]; + buf[i][4] = omega[i][0]; + buf[i][5] = omega[i][1]; + buf[i][6] = omega[i][2]; + buf[i][7] = angmom[i][0]; + buf[i][8] = angmom[i][1]; + buf[i][9] = angmom[i][2]; } } @@ -1715,10 +1773,13 @@ void AtomVecTri::pack_vel(double **buf) int AtomVecTri::pack_vel_hybrid(int i, double *buf) { - buf[0] = angmom[i][0]; - buf[1] = angmom[i][1]; - buf[2] = angmom[i][2]; - return 3; + buf[0] = omega[i][0]; + buf[1] = omega[i][1]; + buf[2] = omega[i][2]; + buf[3] = angmom[i][0]; + buf[4] = angmom[i][1]; + buf[5] = angmom[i][2]; + return 6; } /* ---------------------------------------------------------------------- @@ -1729,9 +1790,10 @@ void AtomVecTri::write_vel(FILE *fp, int n, double **buf) { for (int i = 0; i < n; i++) fprintf(fp,TAGINT_FORMAT - " %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e\n", + " %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e " + "%-1.16e %-1.16e %-1.16e\n", (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][4],buf[i][5],buf[i][6],buf[i][7],buf[i][8],buf[i][9]); } /* ---------------------------------------------------------------------- @@ -1740,8 +1802,9 @@ void AtomVecTri::write_vel(FILE *fp, int n, double **buf) int AtomVecTri::write_vel_hybrid(FILE *fp, double *buf) { - fprintf(fp," %-1.16e %-1.16e %-1.16e",buf[0],buf[1],buf[2]); - return 3; + fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e", + buf[0],buf[1],buf[2],buf[3],buf[4],buf[5]); + return 6; } /* ---------------------------------------------------------------------- @@ -1763,6 +1826,7 @@ bigint AtomVecTri::memory_usage() if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax); if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax); if (atom->memcheck("radius")) bytes += memory->usage(radius,nmax); + if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3); if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3); if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax*comm->nthreads,3); diff --git a/src/atom_vec_tri.h b/src/atom_vec_tri.h index ce08ecc2b0..b5dc2cfde0 100644 --- a/src/atom_vec_tri.h +++ b/src/atom_vec_tri.h @@ -92,7 +92,7 @@ class AtomVecTri : public AtomVec { double **x,**v,**f; tagint *molecule; double *rmass,*radius; - double **angmom,**torque; + double **omega,**angmom,**torque; int *tri; int nlocal_bonus,nghost_bonus,nmax_bonus; diff --git a/src/displace_atoms.cpp b/src/displace_atoms.cpp index 8d2e501641..46f697f099 100644 --- a/src/displace_atoms.cpp +++ b/src/displace_atoms.cpp @@ -232,8 +232,16 @@ void DisplaceAtoms::command(int narg, char **arg) runit[2] = axis[2]/len; double angle = MY_PI*theta/180.0; - double sine = sin(angle); double cosine = cos(angle); + double sine = sin(angle); + + double qcosine = cos(0.5*angle); + double qsine = sin(0.5*angle); + qrotate[0] = qcosine; + qrotate[1] = runit[0]*qsine; + qrotate[2] = runit[1]*qsine; + qrotate[3] = runit[2]*qsine; + double ddotr; // flags for additional orientation info stored by some atom styles @@ -304,10 +312,6 @@ void DisplaceAtoms::command(int narg, char **arg) else if (body_flag && body[i] >= 0) quat = avec_body->bonus[body[i]].quat; if (quat) { - qrotate[0] = cosine; - qrotate[1] = runit[0]*sine; - qrotate[2] = runit[1]*sine; - qrotate[3] = runit[2]*sine; MathExtra::quatquat(qrotate,quat,qnew); quat[0] = qnew[0]; quat[1] = qnew[1]; diff --git a/src/dump_image.cpp b/src/dump_image.cpp index e587e661cb..de0db88777 100644 --- a/src/dump_image.cpp +++ b/src/dump_image.cpp @@ -178,11 +178,13 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : iarg += 3; } else if (strcmp(arg[iarg],"tri") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command"); + if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command"); triflag = YES; if (strcmp(arg[iarg+1],"type") == 0) tcolor = TYPE; else error->all(FLERR,"Illegal dump image command"); - iarg += 2; + tstyle = force->inumeric(FLERR,arg[iarg+2]); + tdiamvalue = force->numeric(FLERR,arg[iarg+3]); + iarg += 4; } else if (strcmp(arg[iarg],"body") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command"); @@ -806,12 +808,19 @@ void DumpImage::create_image() } // render atoms that are triangles + // tstyle = 1 for tri only, 2 for edges only, 3 for both if (triflag) { + int tridraw = 1; + if (tstyle == 2) tridraw = 0; + int edgedraw = 1; + if (tstyle == 1) edgedraw = 0; + double **x = atom->x; int *tri = atom->tri; int *type = atom->type; + for (i = 0; i < nchoose; i++) { j = clist[i]; if (tri[j] < 0) continue; @@ -828,7 +837,12 @@ void DumpImage::create_image() MathExtra::add3(pt2,x[i],pt2); MathExtra::add3(pt3,x[i],pt3); - image->draw_triangle(pt1,pt2,pt3,color); + if (tridraw) image->draw_triangle(pt1,pt2,pt3,color); + if (edgedraw) { + image->draw_cylinder(pt1,pt2,color,tdiamvalue,3); + image->draw_cylinder(pt2,pt3,color,tdiamvalue,3); + image->draw_cylinder(pt3,pt1,color,tdiamvalue,3); + } } } diff --git a/src/dump_image.h b/src/dump_image.h index af6652c3ea..879710838f 100644 --- a/src/dump_image.h +++ b/src/dump_image.h @@ -45,7 +45,8 @@ class DumpImage : public DumpCustom { int lcolor,ldiam; // what determines color/diam of lines double ldiamvalue; // line diameter value int triflag; // 0/1 for draw atoms as triangles - int tcolor; // what determines color of tris + int tcolor,tstyle; // what determines color/style of tris + double tdiamvalue; // tri edge diameter value int bodyflag; // 0/1 for draw atoms as bodies int bodycolor; // what determines color of bodies double bodyflag1,bodyflag2; // user-specified params for drawing bodies diff --git a/src/fix_move.cpp b/src/fix_move.cpp index a0530f2206..ca3727520c 100644 --- a/src/fix_move.cpp +++ b/src/fix_move.cpp @@ -193,11 +193,6 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : "Fix move cannot define z or vz variable for 2d problem"); } - if (atom->angmom_flag && comm->me == 0) - error->warning(FLERR,"Fix move does not update angular momentum"); - if (atom->ellipsoid_flag && comm->me == 0) - error->warning(FLERR,"Fix move does not update quaternions"); - // setup scaling and apply scaling factors to velocity & amplitude if ((mstyle == LINEAR || mstyle == WIGGLE || mstyle == ROTATE) && @@ -609,8 +604,15 @@ void FixMove::initial_integrate(int vflag) } else if (mstyle == ROTATE) { double arg = omega_rotate * delta; - double sine = sin(arg); double cosine = cos(arg); + double sine = sin(arg); + + double qcosine = cos(0.5*arg); + double qsine = sin(0.5*arg); + qrotate[0] = qcosine; + qrotate[1] = runit[0]*qsine; + qrotate[2] = runit[1]*qsine; + qrotate[3] = runit[2]*qsine; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -646,12 +648,13 @@ void FixMove::initial_integrate(int vflag) if (extra_flag) { - // omega for spheres and lines + // omega for spheres, lines, tris if (omega_flag) { flag = 0; if (radius_flag && radius[i] > 0.0) flag = 1; if (line_flag && line[i] >= 0.0) flag = 1; + if (tri_flag && tri[i] >= 0.0) flag = 1; if (flag) { omega[i][0] = omega_rotate*runit[0]; omega[i][1] = omega_rotate*runit[1]; @@ -706,13 +709,7 @@ void FixMove::initial_integrate(int vflag) quat = avec_tri->bonus[tri[i]].quat; else if (body_flag && body[i] >= 0) quat = avec_body->bonus[body[i]].quat; - if (quat) { - qrotate[0] = cosine; - qrotate[1] = runit[0]*sine; - qrotate[2] = runit[1]*sine; - qrotate[3] = runit[2]*sine; - MathExtra::quatquat(qrotate,qoriginal[i],quat); - } + if (quat) MathExtra::quatquat(qrotate,qoriginal[i],quat); } }