more work on dump image

src/dump_grid.cpp

@@ -26,9 +26,6 @@
 #include "region.h"
 #include "update.h"
 
-// DEBUG
-#include "comm.h"
-
 #include <cstring>
 
 using namespace LAMMPS_NS;
@@ -652,6 +649,37 @@ int DumpGrid::parse_fields(int narg, char **arg)
 
   for (int iarg = 0; iarg < narg; iarg++) {
 
+    char *id;
+    int igrid,idata,index;
+    
+    int iflag =
+      utils::check_grid_reference(igrid,idata,index,lmp);
+    //utils::check_grid_reference((char *) "Dump grid",
+    //                              arg[iarg],igrid,idata,index,lmp);
+
+    // arg is not a Grid reference
+    
+    if (iflag < 0) return iarg;
+    
+    // grid reference is to a compute
+    
+    if (iflag == ArgInfo::COMPUTE) {
+
+      
+    // grid reference is to a fix
+
+    } else if (iflag == ArgInfo::FIX) {
+
+    }
+  }
+
+  return narg;
+}
+
+
+
+
+/*   
     ArgInfo argi(arg[iarg], ArgInfo::COMPUTE | ArgInfo::FIX);
     argindex[iarg] = argi.get_index1();
     auto name = argi.get_name();
@@ -762,6 +790,7 @@ int DumpGrid::parse_fields(int narg, char **arg)
 
   return narg;
 }
+*/
 
 /* ----------------------------------------------------------------------
    add Compute to list of Compute objects used by dump

src/dump_image.cpp

@@ -21,10 +21,13 @@
 #include "atom_vec_tri.h"
 #include "body.h"
 #include "comm.h"
+#include "compute.h"
 #include "domain.h"
 #include "error.h"
 #include "fix.h"
 #include "force.h"
+#include "grid2d.h"
+#include "grid3d.h"
 #include "image.h"
 #include "input.h"
 #include "math_const.h"
@@ -33,6 +36,7 @@
 #include "modify.h"
 #include "molecule.h"
 #include "tokenizer.h"
+#include "update.h"
 #include "variable.h"
 
 #include <cmath>
@@ -179,10 +183,14 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg],"grid") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
       gridflag = YES;
-      if (strcmp(arg[iarg+1],"type") == 0) lcolor = TYPE;
+
-      else error->all(FLERR,"Illegal dump image command");
+      char *id;
-      ldiam = NUMERIC;
+      int igrid,idata,index;
-      ldiamvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp);
+      int iflag =
+        utils::check_grid_reference(igrid,idata,index,lmp);
+      //        utils::check_grid_reference((char *) "Dump image",
+      //                             arg[iarg+1],igrid,idata,index,lmp);
+      
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"line") == 0) {
@@ -453,7 +461,23 @@ void DumpImage::init_style()
 
   DumpCustom::init_style();
 
-  // check variables
+  // for grid output, find current ptr for compute or fix
+  // check that fix frequency is acceptable
+
+  if (gridflag) {
+    if (id_grid_compute) {
+      grid_compute = modify->get_compute_by_id(id_grid_compute);
+      if (!grid_compute)
+        error->all(FLERR,"Could not find dump image grid compute ID {}",id_grid_compute);
+    } else if (id_grid_fix) {
+      grid_fix = modify->get_fix_by_id(id_grid_fix);
+      if (!grid_fix) error->all(FLERR,"Could not find dump image fix ID {}",id_grid_fix);
+      if (nevery % grid_fix->peratom_freq)
+        error->all(FLERR,"Dump image and grid fix not computed at compatible times");
+    }
+  }
+  
+  // check image variables
 
   if (thetastr) {
     thetavar = input->variable->find(thetastr);
@@ -562,7 +586,7 @@ void DumpImage::write()
     memory->create(buf,maxbuf*size_one,"dump:buf");
   }
 
-  // pack buf with color & diameter
+  // pack atom buf with color & diameter
 
   pack(nullptr);
 
@@ -582,26 +606,95 @@ void DumpImage::write()
     two[1] = hi;
     MPI_Allreduce(two,twoall,2,MPI_DOUBLE,MPI_MAX,world);
     int flag = image->map_minmax(0,-twoall[0],twoall[1]);
-    if (flag) error->all(FLERR,"Invalid color map min/max values");
+    if (flag) error->all(FLERR,"Invalid atom color map min/max values");
+  }
+
+  // pack grid gbuf with grid cell values
+  // ngrid_mine = # of grid cells this proc owns
+
+  if (gridflag) {
+    if (domain->dimension == 2) {
+      if (grid_compute)
+        grid2d = (Grid2d *) grid_compute->get_grid_by_index(grid_index);
+      else if (grid_fix)
+        grid2d = (Grid2d *) grid_fix->get_grid_by_index(grid_index);
+      grid2d->get_bounds_owned(nxlo_in,nxhi_in,nylo_in,nyhi_in);
+    } else {
+      if (grid_compute)
+        grid3d = (Grid3d *) grid_compute->get_grid_by_index(grid_index);
+      else if (grid_fix)
+        grid3d = (Grid3d *) grid_fix->get_grid_by_index(grid_index);
+      grid3d->get_bounds_owned(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in);
+    }
+
+    // invoke Compute for per-grid quantities
+    // only if within a run or minimize
+    // else require the compute is current
+    // this prevents the compute from being invoked by the WriteDump class
+
+    if (grid_compute) {
+      if (update->whichflag == 0) {
+        if (grid_compute->invoked_pergrid != update->ntimestep)
+          error->all(FLERR,"Grid compute {} used in dump image between runs is not current",
+                     grid_compute->id);
+      } else {
+        if (!(grid_compute->invoked_flag & Compute::INVOKED_PERGRID)) {
+          grid_compute->compute_pergrid();
+          grid_compute->invoked_flag |= Compute::INVOKED_PERGRID;
+        }
+      }
+    }
+
+    // access grid data and load gbuf
+
+    /*
+    if (index == 0) {
+      double **vec2d;
+      if (field2source[n] == COMPUTE)
+        vec2d = (double **)
+          compute[field2index[n]]->get_griddata_by_index(field2data[n]);
+      else if (field2source[n] == FIX)
+        vec2d = (double **)
+          fix[field2index[n]]->get_griddata_by_index(field2data[n]);
+      for (int iy = nylo_in; iy <= nyhi_in; iy++)
+        for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
+          buf[n] = vec2d[iy][ix];
+          n += size_one;
+        }
+    } else {
+      double ***array2d;
+      if (field2source[n] == COMPUTE)
+        array2d = (double ***)
+          compute[field2index[n]]->get_griddata_by_index(field2data[n]);
+      else if (field2source[n] == FIX)
+        array2d = (double ***)
+          fix[field2index[n]]->get_griddata_by_index(field2data[n]);
+      index--;
+      for (int iy = nylo_in; iy <= nyhi_in; iy++)
+        for (int ix = nxlo_in; ix <= nxhi_in; ix++) {
+          buf[n] = array2d[iy][ix][index];
+          n += size_one;
+        }
+    }
+    */
+
   }
 
   // set minmax color range if using dynamic grid color map
 
-  if (acolor == ATTRIBUTE && image->map_dynamic(1)) {
+  if (gridflag && image->map_dynamic(1)) {
     double two[2],twoall[2];
     double lo = BIG;
     double hi = -BIG;
-    int m = 0;
+    for (int i = 0; i < ngrid_owned; i++) {
-    for (int i = 0; i < nchoose; i++) {
+      lo = MIN(lo,gbuf[i]);
-      lo = MIN(lo,buf[m]);
+      hi = MAX(hi,gbuf[i]);
-      hi = MAX(hi,buf[m]);
-      m += size_one;
     }
     two[0] = -lo;
     two[1] = hi;
     MPI_Allreduce(two,twoall,2,MPI_DOUBLE,MPI_MAX,world);
-    int flag = image->map_minmax(0,-twoall[0],twoall[1]);
+    int flag = image->map_minmax(1,-twoall[0],twoall[1]);
-    if (flag) error->all(FLERR,"Invalid color map min/max values");
+    if (flag) error->all(FLERR,"Invalid grid color map min/max values");
   }
 
   // create image on each proc, then merge them
@@ -702,8 +795,9 @@ void DumpImage::view_params()
 }
 
 /* ----------------------------------------------------------------------
-   create image for atoms on this proc
+   create image for all data this proc owns
-   every pixel has depth
+   all procs draw simulation box edges if requested
+   every drawn pixel has depth so merge can decide which to keep
 ------------------------------------------------------------------------- */
 
 void DumpImage::create_image()
@@ -767,6 +861,16 @@ void DumpImage::create_image()
     }
   }
 
+  // render my grid cells
+
+  if (gridflag) {
+
+
+    // draw 2 or 12 triangles
+    //image->draw_triangle(x,y,z,color);
+    
+  }
+
   // render atoms that are lines
 
   if (lineflag) {

src/dump_image.h

@@ -41,8 +41,6 @@ class DumpImage : public DumpCustom {
   int acolor, adiam;    // what determines color/diam of atoms
   double adiamvalue;    // atom diameter value
 
-  int gridflag;         // 0/1 for draw grid cells
-  
   int lineflag;                   // 0/1 for draw atoms as lines
   int lcolor, ldiam;              // what determines color/diam of lines
   double ldiamvalue;              // line diameter value
@@ -81,16 +79,27 @@ class DumpImage : public DumpCustom {
   double *diamtype, *diamelement, *bdiamtype;          // per-type diameters
   double **colortype, **colorelement, **bcolortype;    // per-type colors
 
+  int gridflag;         // 0/1 for draw grid cells
+  class Grid2d *grid2d;
+  class Grid3d *grid3d;
+  char *id_grid_compute,*id_grid_fix;
+  class Compute *grid_compute;
+  class Fix *grid_fix;
+  int grid_index;
+  double *gbuf;
+  int nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in;
+  int ngrid_owned;
+  
   class AtomVecLine *avec_line;    // ptrs to atom style (sub)classes
   class AtomVecTri *avec_tri;
   class AtomVecBody *avec_body;
 
   class Fix *fixptr;    // ptr to Fix that provides image data
 
-  class Image *image;    // class that renders each image
+  class Image *image;   // class that renders each image
 
-  int *chooseghost;    // extended choose array for comm
+  int *chooseghost;     // extended choose array for comm
-  double **bufcopy;    // buffer for communicating bond/atom info
+  double **bufcopy;     // buffer for communicating bond/atom info
   int maxbufcopy;
 
   void init_style() override;

src/fix_ave_grid.cpp

@@ -312,7 +312,7 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) :
     for (int i = 0; i < nvalues; i++) {
       if (which[i] == ArgInfo::COMPUTE) {
 
-        auto words = utils::parse_gridid(FLERR,ids[i],error);
+        auto words = utils::parse_grid_id(FLERR,ids[i],error);
         const auto &idcompute = words[0];
         const auto &gname = words[1];
         const auto &dname = words[2];
@@ -350,7 +350,7 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) :
 
       } else if (which[i] == ArgInfo::FIX) {
 
-        auto words = utils::parse_gridid(FLERR,ids[i],error);
+        auto words = utils::parse_grid_id(FLERR,ids[i],error);
         const auto &idfix = words[0];
         const auto &gname = words[1];
         const auto &dname = words[2];

src/utils.cpp

@@ -13,6 +13,7 @@
 
 #include "utils.h"
 
+#include "arg_info.h"
 #include "atom.h"
 #include "comm.h"
 #include "compute.h"
@@ -632,7 +633,7 @@ int utils::expand_args(const char *file, int line, int narg, char **arg, int mod
     // match grids
 
     if (strmatch(word, "^[cf]_\\w+:\\w+:\\w+\\[\\d*\\*\\d*\\]")) {
-      auto gridid = utils::parse_gridid(FLERR, word, lmp->error);
+      auto gridid = utils::parse_grid_id(FLERR, word, lmp->error);
 
       size_t first = gridid[2].find('[');
       size_t second = gridid[2].find(']', first + 1);
@@ -863,13 +864,123 @@ char *utils::expand_type(const char *file, int line, const std::string &str, int
     return nullptr;
 }
 
+/* ----------------------------------------------------------------------
+   Check grid reference for valid Compute or Fix which produces per-grid data
+   errstr = name of calling command used if error is generated
+   ref = grid reference as it appears in an input script
+   return arguments:
+     id = ID of compute or fix
+     igrid = index of which grid in compute/fix
+     idata = index of which data field in igrid
+     index = index into data field (0 for vector, 1-N for column of array)
+   method return = ArgInfo::COMPUTE or ArgInfo::FIX or -1 for neither
+     caller decides what to do if not COMPUTE or FIX
+------------------------------------------------------------------------- */
+
+//int check_grid_reference(char *errstr, char *ref, int &igrid, int &idata, int &index,
+//                         LAMMPS *lmp)
+int check_grid_reference(int &igrid, int &idata, int &index, LAMMPS *lmp)
+{
+  char *ref;
+  char *errstr;
+  
+  ArgInfo argi(ref, ArgInfo::COMPUTE | ArgInfo::FIX);
+  index = argi.get_index1();
+  auto name = argi.get_name();
+
+  switch (argi.get_type()) {
+    
+    case ArgInfo::UNKNOWN: {
+      lmp->error->all(FLERR,"%s grid reference %s is invalid",errstr,ref);
+    } break;
+      
+    // compute value = c_ID
+
+    case ArgInfo::COMPUTE: {
+
+      // split name = idcompute:gname:dname into 3 strings
+  
+      auto words = utils::parse_grid_id(FLERR,name,lmp->error);
+      const auto &idcompute = words[0];
+      const auto &gname = words[1];
+      const auto &dname = words[2];
+      
+      auto icompute = lmp->modify->get_compute_by_id(idcompute);
+      if (!icompute) lmp->error->all(FLERR,"%s compute ID {} not found",errstr,idcompute);
+      if (icompute->pergrid_flag == 0)
+        lmp->error->all(FLERR,"%s compute {} does not compute per-grid info",errstr,idcompute);
+      
+      int dim;
+      igrid = icompute->get_grid_by_name(gname,dim);
+      if (igrid < 0)
+        lmp->error->all(FLERR,"%s compute {} does not recognize grid name {}",errstr,idcompute,gname);
+
+      int ncol;
+      idata = icompute->get_griddata_by_name(igrid,dname,ncol);
+      if (idata < 0)
+        lmp->error->all(FLERR,"%s compute {} does not recognize data name {}",errstr,idcompute,dname);
+
+      if (argi.get_dim() == 0 && ncol)
+        lmp->error->all(FLERR,"%s compute {} data {} is not per-grid vector",errstr,idcompute,dname);
+      if (argi.get_dim() && ncol == 0)
+        lmp->error->all(FLERR,"%s compute {} data {} is not per-grid array",errstr,idcompute,dname);
+      if (argi.get_dim() && argi.get_index1() > ncol)
+        lmp->error->all(FLERR,"%s compute {} array {} is accessed out-of-range",errstr,idcompute,dname);
+
+      //id = utils::strdup(idcompute);
+      return ArgInfo::COMPUTE;
+    } break;
+
+    // fix value = f_ID
+
+    case ArgInfo::FIX: {
+
+      // split name = idfix:gname:dname into 3 strings
+
+      auto words = utils::parse_grid_id(FLERR,name,lmp->error);
+      const auto &idfix = words[0];
+      const auto &gname = words[1];
+      const auto &dname = words[2];
+
+      auto ifix = lmp->modify->get_fix_by_id(idfix);
+      if (!ifix) lmp->error->all(FLERR,"%s fix ID {} not found",errstr,idfix);
+      if (ifix->pergrid_flag == 0)
+        lmp->error->all(FLERR,"%s fix {} does not compute per-grid info",errstr,idfix);
+      //if (nevery % ifix->pergrid_freq)
+      //  lmp->error->all(FLERR,"%s fix {} not computed at compatible time",errstr,if);
+
+      int dim;
+      int igrid = ifix->get_grid_by_name(gname,dim);
+      if (igrid < 0)
+        lmp->error->all(FLERR,"%s fix {} does not recognize grid name {}",errstr,idfix,gname);
+
+      int ncol;
+      int idata = ifix->get_griddata_by_name(igrid,dname,ncol);
+      if (idata < 0)
+        lmp->error->all(FLERR,"%s fix {} does not recognize data name {}",errstr,idfix,dname);
+
+      if (argi.get_dim() == 0 && ncol)
+        lmp->error->all(FLERR,"%s fix {} data {} is not per-grid vector",errstr,idfix,dname);
+      if (argi.get_dim() > 0 && ncol == 0)
+        lmp->error->all(FLERR,"%s fix {} data {} is not per-grid array",errstr,idfix,dname);
+      if (argi.get_dim() > 0 && argi.get_index1() > ncol)
+        lmp->error->all(FLERR,"%s fix {} array {} is accessed out-of-range",errstr,idfix,dname);
+
+      //id = utils::strdup(idfix);
+      return ArgInfo::FIX;
+    } break;
+  }
+
+  return -1;
+}
+
 /* ----------------------------------------------------------------------
    Parse grid reference into id:gridname:dataname
    return vector of 3 substrings
 ------------------------------------------------------------------------- */
 
-std::vector<std::string> utils::parse_gridid(const char *file, int line, const std::string &name,
+std::vector<std::string> utils::parse_grid_id(const char *file, int line, const std::string &name,
-                                             Error *error)
+                                              Error *error)
 {
   auto words = Tokenizer(name, ":").as_vector();
   if (words.size() != 3) {

src/utils.h

@@ -379,6 +379,15 @@ namespace utils {
 
   char *expand_type(const char *file, int line, const std::string &str, int mode, LAMMPS *lmp);
 
+
+
+
+
+  //int check_grid_reference(char *errstr, char *ref, int &igrid, int &idata, int &index,
+  //                         LAMMPS *lmp);
+  int check_grid_reference(int &igrid, int &idata, int &index, LAMMPS *lmp);
+
+  
   /*! Parse grid reference into 3 sub-strings
    *
    * Format of grid ID reference = id:gname:dname
@@ -387,8 +396,8 @@ namespace utils {
    * \param name = complete grid ID
    * \return std::vector<std::string> containing the 3 sub-strings  */
 
-  std::vector<std::string> parse_gridid(const char *file, int line, const std::string &name,
+  std::vector<std::string> parse_grid_id(const char *file, int line, const std::string &name,
-                                        Error *error);
+                                         Error *error);
 
   /*! Make C-style copy of string in new storage
    *