From 9a19777d4f0e3cd878dca237a06e47d452d8be3f Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 5 Sep 2013 22:42:44 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10722 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/RIGID/fix_rigid.cpp | 46 ++++++++++++++++++++++++++- src/RIGID/fix_rigid.h | 4 ++- src/RIGID/fix_rigid_small.cpp | 58 +++++++++++++++++++++++++++++++++++ src/RIGID/fix_rigid_small.h | 4 ++- src/dump_local.cpp | 4 +++ 5 files changed, 113 insertions(+), 3 deletions(-) diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index 3ec7225549..c4051d495f 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -2547,7 +2547,6 @@ double FixRigid::compute_scalar() double wbody[3],rot[3][3]; double t = 0.0; - for (int i = 0; i < nbody; i++) { t += masstotal[i] * (fflag[i][0]*vcm[i][0]*vcm[i][0] + fflag[i][1]*vcm[i][1]*vcm[i][1] + @@ -2593,6 +2592,51 @@ void *FixRigid::extract(const char *str, int &dim) return NULL; } +/* ---------------------------------------------------------------------- + return translational KE for all rigid bodies + KE = 1/2 M Vcm^2 +------------------------------------------------------------------------- */ + +double FixRigid::extract_ke() +{ + double ke = 0.0; + for (int i = 0; i < nbody; i++) + ke += masstotal[i] * + (vcm[i][0]*vcm[i][0] + vcm[i][1]*vcm[i][1] + vcm[i][2]*vcm[i][2]); + + return 0.5*ke; +} + +/* ---------------------------------------------------------------------- + return rotational KE for all rigid bodies + Erotational = 1/2 I wbody^2 +------------------------------------------------------------------------- */ + +double FixRigid::extract_erotational() +{ + double wbody[3],rot[3][3]; + + double erotate = 0.0; + for (int i = 0; i < nbody; i++) { + + // wbody = angular velocity in body frame + + MathExtra::quat_to_mat(quat[i],rot); + MathExtra::transpose_matvec(rot,angmom[i],wbody); + if (inertia[i][0] == 0.0) wbody[0] = 0.0; + else wbody[0] /= inertia[i][0]; + if (inertia[i][1] == 0.0) wbody[1] = 0.0; + else wbody[1] /= inertia[i][1]; + if (inertia[i][2] == 0.0) wbody[2] = 0.0; + else wbody[2] /= inertia[i][2]; + + erotate += inertia[i][0]*wbody[0]*wbody[0] + + inertia[i][1]*wbody[1]*wbody[1] + inertia[i][2]*wbody[2]*wbody[2]; + } + + return 0.5*erotate; +} + /* ---------------------------------------------------------------------- return attributes of a rigid body 15 values per body diff --git a/src/RIGID/fix_rigid.h b/src/RIGID/fix_rigid.h index e30c1b9681..a08d9c87d5 100644 --- a/src/RIGID/fix_rigid.h +++ b/src/RIGID/fix_rigid.h @@ -53,7 +53,9 @@ class FixRigid : public Fix { void reset_dt(); void zero_momentum(int); void zero_rotation(int); - virtual void *extract(const char*,int &); + virtual void *extract(const char*, int &); + double extract_ke(); + double extract_erotational(); double compute_array(int, int); protected: diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 986a2660f8..98a6224d30 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -3061,6 +3061,64 @@ void *FixRigidSmall::extract(const char *str, int &dim) return NULL; } +/* ---------------------------------------------------------------------- + return translational KE for all rigid bodies + KE = 1/2 M Vcm^2 + sum local body results across procs +------------------------------------------------------------------------- */ + +double FixRigidSmall::extract_ke() +{ + double *vcm; + + double ke = 0.0; + for (int i = 0; i < nlocal_body; i++) { + vcm = body[i].vcm; + ke += body[i].mass * (vcm[0]*vcm[0] + vcm[1]*vcm[1] + vcm[2]*vcm[2]); + } + + double keall; + MPI_Allreduce(&ke,&keall,1,MPI_DOUBLE,MPI_SUM,world); + + return 0.5*keall; +} + +/* ---------------------------------------------------------------------- + return rotational KE for all rigid bodies + Erotational = 1/2 I wbody^2 +------------------------------------------------------------------------- */ + +double FixRigidSmall::extract_erotational() +{ + double wbody[3],rot[3][3]; + double *inertia; + + double erotate = 0.0; + for (int i = 0; i < nlocal_body; i++) { + + // for Iw^2 rotational term, need wbody = angular velocity in body frame + // not omega = angular velocity in space frame + + inertia = body[i].inertia; + MathExtra::quat_to_mat(body[i].quat,rot); + MathExtra::transpose_matvec(rot,body[i].angmom,wbody); + if (inertia[0] == 0.0) wbody[0] = 0.0; + else wbody[0] /= inertia[0]; + if (inertia[1] == 0.0) wbody[1] = 0.0; + else wbody[1] /= inertia[1]; + if (inertia[2] == 0.0) wbody[2] = 0.0; + else wbody[2] /= inertia[2]; + + erotate += inertia[0]*wbody[0]*wbody[0] + inertia[1]*wbody[1]*wbody[1] + + inertia[2]*wbody[2]*wbody[2]; + } + + double erotateall; + MPI_Allreduce(&erotate,&erotateall,1,MPI_DOUBLE,MPI_SUM,world); + + return 0.5*erotateall; +} + /* ---------------------------------------------------------------------- return temperature of collection of rigid bodies non-active DOF are removed by fflag/tflag and in tfactor diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index 0249901765..650eb34306 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -62,7 +62,9 @@ class FixRigidSmall : public Fix { void reset_dt(); void zero_momentum(int); void zero_rotation(int); - void *extract(const char*,int &); + void *extract(const char*, int &); + double extract_ke(); + double extract_erotational(); double compute_scalar(); double memory_usage(); diff --git a/src/dump_local.cpp b/src/dump_local.cpp index 2812fb9a5c..93bf4982d1 100644 --- a/src/dump_local.cpp +++ b/src/dump_local.cpp @@ -147,6 +147,10 @@ void DumpLocal::init_style() vformat[i] = strcat(vformat[i]," "); } + // setup boundary string + + domain->boundary_string(boundstr); + // find current ptr for each compute,fix,variable // check that fix frequency is acceptable