diff --git a/src/dump_grid.cpp b/src/dump_grid.cpp index 71ddcee5f2..a11b042de4 100644 --- a/src/dump_grid.cpp +++ b/src/dump_grid.cpp @@ -26,9 +26,6 @@ #include "region.h" #include "update.h" -// DEBUG -#include "comm.h" - #include 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 diff --git a/src/dump_image.cpp b/src/dump_image.cpp index 4247ebfa95..f76368d0c4 100644 --- a/src/dump_image.cpp +++ b/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 @@ -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"); - ldiam = NUMERIC; - ldiamvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); + + char *id; + int igrid,idata,index; + 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 < nchoose; i++) { - lo = MIN(lo,buf[m]); - hi = MAX(hi,buf[m]); - m += size_one; + for (int i = 0; i < ngrid_owned; i++) { + lo = MIN(lo,gbuf[i]); + hi = MAX(hi,gbuf[i]); } 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]); - if (flag) error->all(FLERR,"Invalid color map min/max values"); + int flag = image->map_minmax(1,-twoall[0],twoall[1]); + 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 - every pixel has depth + create image for all data this proc owns + 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) { diff --git a/src/dump_image.h b/src/dump_image.h index a455e5df19..c4ce2224f4 100644 --- a/src/dump_image.h +++ b/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 - double **bufcopy; // buffer for communicating bond/atom info + int *chooseghost; // extended choose array for comm + double **bufcopy; // buffer for communicating bond/atom info int maxbufcopy; void init_style() override; diff --git a/src/fix_ave_grid.cpp b/src/fix_ave_grid.cpp index 06c0175503..fe28b9a8e2 100644 --- a/src/fix_ave_grid.cpp +++ b/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]; diff --git a/src/utils.cpp b/src/utils.cpp index c4d60a799b..77d3699f18 100644 --- a/src/utils.cpp +++ b/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 utils::parse_gridid(const char *file, int line, const std::string &name, - Error *error) +std::vector utils::parse_grid_id(const char *file, int line, const std::string &name, + Error *error) { auto words = Tokenizer(name, ":").as_vector(); if (words.size() != 3) { diff --git a/src/utils.h b/src/utils.h index 8bf2c542ee..417884442a 100644 --- a/src/utils.h +++ b/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 containing the 3 sub-strings */ - std::vector parse_gridid(const char *file, int line, const std::string &name, - Error *error); + std::vector parse_grid_id(const char *file, int line, const std::string &name, + Error *error); /*! Make C-style copy of string in new storage *