git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10722 f3b2605a-c512-4ea7-a41b-209d697bcdaa

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

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:

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

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();
 

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