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

src/compute_reduce.cpp

@@ -16,6 +16,7 @@
 #include "compute_reduce.h"
 #include "atom.h"
 #include "update.h"
+#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "force.h"
@@ -40,12 +41,22 @@ enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
 ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
-  if (narg < 5) error->all("Illegal compute reduce command");
+  int iarg;
+  if (strcmp(style,"reduce") == 0) {
+    if (narg < 5) error->all("Illegal compute reduce command");
+    iarg = 3;
+  } else if (strcmp(style,"reduce/region") == 0) {
+    if (narg < 6) error->all("Illegal compute reduce/region command");
+    iregion = domain->find_region(arg[3]);
+    if (iregion == -1) error->all("Compute reduce region ID does not exist");
+    iarg = 4;
+  }
 
-  if (strcmp(arg[3],"sum") == 0) mode = SUM;
-  else if (strcmp(arg[3],"min") == 0) mode = MINN;
-  else if (strcmp(arg[3],"max") == 0) mode = MAXX;
+  if (strcmp(arg[iarg],"sum") == 0) mode = SUM;
+  else if (strcmp(arg[iarg],"min") == 0) mode = MINN;
+  else if (strcmp(arg[iarg],"max") == 0) mode = MAXX;
   else error->all("Illegal compute reduce command");
+  iarg++;
 
   // parse remaining values
 
@@ -55,7 +66,6 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
   value2index = new int[narg-4];
   nvalues = 0;
 
-  int iarg = 4;
   while (iarg < narg) {
     ids[nvalues] = NULL;
 

src/compute_reduce.h

@@ -21,14 +21,14 @@ namespace LAMMPS_NS {
 class ComputeReduce : public Compute {
  public:
   ComputeReduce(class LAMMPS *, int, char **);
-  ~ComputeReduce();
+  virtual ~ComputeReduce();
   void init();
   double compute_scalar();
   void compute_vector();
   double memory_usage();
 
- private:
-  int mode,nvalues;
+ protected:
+  int mode,nvalues,iregion;
   int *which,*argindex,*value2index;
   char **ids;
   double *onevec;
@@ -36,7 +36,7 @@ class ComputeReduce : public Compute {
   int maxatom;
   double *varatom;
 
-  double compute_one(int);
+  virtual double compute_one(int);
   void combine(double &, double);
 };
 

src/compute_reduce_region.cpp

@@ -0,0 +1,138 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "string.h"
+#include "stdlib.h"
+#include "compute_reduce_region.h"
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "domain.h"
+#include "region.h"
+#include "fix.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{SUM,MINN,MAXX};
+enum{X,V,F,COMPUTE,FIX,VARIABLE};
+enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+
+#define BIG 1.0e20
+
+/* ---------------------------------------------------------------------- */
+
+ComputeReduceRegion::ComputeReduceRegion(LAMMPS *lmp, int narg, char **arg) :
+  ComputeReduce(lmp, narg, arg) {}
+
+/* ---------------------------------------------------------------------- */
+
+double ComputeReduceRegion::compute_one(int m)
+{
+  int i;
+
+  Region *region = domain->regions[iregion];
+
+  // invoke the appropriate attribute,compute,fix,variable
+  // compute scalar quantity by summing over atom scalars
+  // only include atoms in group
+
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  int n = value2index[m];
+  int j = argindex[m];
+
+  double one;
+  if (mode == SUM) one = 0.0;
+  else if (mode == MINN) one = BIG;
+  else if (mode == MAXX) one = -BIG;
+
+  if (which[m] == X) {
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+	combine(one,x[i][j]);
+  } else if (which[m] == V) {
+    double **v = atom->v;
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+	combine(one,v[i][j]);
+  } else if (which[m] == F) {
+    double **f = atom->f;
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+	combine(one,f[i][j]);
+    
+  // invoke compute if not previously invoked
+
+  } else if (which[m] == COMPUTE) {
+    Compute *compute = modify->compute[n];
+    if (!(compute->invoked_flag & INVOKED_PERATOM)) {
+      compute->compute_peratom();
+      compute->invoked_flag |= INVOKED_PERATOM;
+    }
+
+    if (j == 0) {
+      double *compute_scalar = compute->scalar_atom;
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+	  combine(one,compute_scalar[i]);
+    } else {
+      double **compute_vector = compute->vector_atom;
+      int jm1 = j - 1;
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+	  combine(one,compute_vector[i][jm1]);
+    }
+
+  // access fix fields, check if frequency is a match
+
+  } else if (which[m] == FIX) {
+    if (update->ntimestep % modify->fix[n]->peratom_freq)
+      error->all("Fix used in compute reduce not computed at compatible time");
+
+    if (j == 0) {
+      double *fix_scalar = modify->fix[n]->scalar_atom;
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+	  combine(one,fix_scalar[i]);
+    } else {
+      double **fix_vector = modify->fix[n]->vector_atom;
+      int jm1 = j - 1;
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+	  combine(one,fix_vector[i][jm1]);
+    }
+    
+  // evaluate atom-style variable
+
+  } else if (which[m] == VARIABLE) {
+    if (nlocal > maxatom) {
+      maxatom = atom->nmax;
+      memory->sfree(varatom);
+      varatom =	(double *) 
+	memory->smalloc(maxatom*sizeof(double),"compute/reduce:varatom");
+    }
+
+    input->variable->compute_atom(n,igroup,varatom,1,0);
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+	combine(one,varatom[i]);
+  }
+
+  return one;
+}

src/compute_reduce_region.h

@@ -0,0 +1,32 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef COMPUTE_REDUCE_REGION_H
+#define COMPUTE_REDUCE_REGION_H
+
+#include "compute_reduce.h"
+
+namespace LAMMPS_NS {
+
+class ComputeReduceRegion : public ComputeReduce {
+ public:
+  ComputeReduceRegion(class LAMMPS *, int, char **);
+  ~ComputeReduceRegion() {}
+
+ private:
+  double compute_one(int);
+};
+
+}
+
+#endif

src/style.h

@@ -87,6 +87,7 @@ CommandStyle(write_restart,WriteRestart)
 #include "compute_pe_atom.h"
 #include "compute_pressure.h"
 #include "compute_reduce.h"
+#include "compute_reduce_region.h"
 #include "compute_erotate_sphere.h"
 #include "compute_stress_atom.h"
 #include "compute_temp.h"
@@ -110,6 +111,7 @@ ComputeStyle(pe,ComputePE)
 ComputeStyle(pe/atom,ComputePEAtom)
 ComputeStyle(pressure,ComputePressure)
 ComputeStyle(reduce,ComputeReduce)
+ComputeStyle(reduce/region,ComputeReduceRegion)
 ComputeStyle(erotate/sphere,ComputeERotateSphere)
 ComputeStyle(stress/atom,ComputeStressAtom)
 ComputeStyle(temp,ComputeTemp)