diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index fe0b08edcb..143fbbf09c 100644 --- a/src/compute_reduce.cpp +++ b/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; diff --git a/src/compute_reduce.h b/src/compute_reduce.h index 265f85605f..82fe48ffb6 100644 --- a/src/compute_reduce.h +++ b/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); }; diff --git a/src/compute_reduce_region.cpp b/src/compute_reduce_region.cpp new file mode 100644 index 0000000000..8c80d75392 --- /dev/null +++ b/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; +} diff --git a/src/compute_reduce_region.h b/src/compute_reduce_region.h new file mode 100644 index 0000000000..6aaf1741e5 --- /dev/null +++ b/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 diff --git a/src/style.h b/src/style.h index 8cd3a0ffcc..4af24afe8b 100644 --- a/src/style.h +++ b/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)