Updates for standalone VTK file writing

exputil/CMakeLists.txt

@@ -8,10 +8,10 @@ set(UTIL_SRC nrutil.cc elemfunc.cc euler.cc euler_slater.cc # Hankel.cc
   rotmatrix.cc wordSplit.cc FileUtils.cc BarrierWrapper.cc stack.cc
   localmpi.cc TableGrid.cc writePVD.cc libvars.cc TransformFFT.cc QDHT.cc
   YamlCheck.cc parseVersionString.cc EXPmath.cc laguerre_polynomial.cpp
-  YamlConfig.cc orthoTest.cc OrthoFunction.cc)
+  YamlConfig.cc orthoTest.cc OrthoFunction.cc VtkGrid.cc)
 
 if(HAVE_VTK)
-  list(APPEND UTIL_SRC VtkGrid.cc VtkPCA.cc)
+  list(APPEND UTIL_SRC VtkPCA.cc)
 endif()
 
 set(OPTIMIZATION_SRC SimAnn.cc)
@@ -51,13 +51,9 @@ set(common_INCLUDE_DIRS $<INSTALL_INTERFACE:include>
   ${DEP_INC} ${EIGEN3_INCLUDE_DIR} ${HDF5_INCLUDE_DIRS}
   ${FFTW_INCLUDE_DIRS})
 
-# set(common_LINKLIBS ${DEP_LIB} OpenMP::OpenMP_CXX MPI::MPI_CXX
-#   yaml-cpp ${VTK_LIBRARIES} ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES}
-#   ${HDF5_HL_LIBRARIES} ${FFTW_DOUBLE_LIB})
-
-set(common_LINKLIB ${DEP_LIB} OpenMP::OpenMP_CXX MPI::MPI_CXX
-  yaml-cpp ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}
-  ${FFTW_DOUBLE_LIB})
+set(common_LINKLIB ${DEP_LIB} OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp
+  ${VTK_LIBRARIES} ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES}
+  ${HDF5_HL_LIBRARIES} ${FFTW_DOUBLE_LIB})
 
 if(ENABLE_CUDA)
   list(APPEND common_LINKLIB CUDA::toolkit CUDA::cudart)

exputil/VtkGrid.cc

@@ -1,5 +1,8 @@
+#include <iomanip>
 #include <VtkGrid.H>
 
+#ifdef HAVE_VTK
+
 VtkGrid::VtkGrid(int nx, int ny, int nz,
 		 double xmin, double xmax,
 		 double ymin, double ymax,
@@ -126,3 +129,161 @@ void VtkGrid::Write(const std::string& name)
   // writer->SetDataModeToBinary();
   writer->Write();
 }
+
+#else
+
+VtkGrid::VtkGrid(int nx, int ny, int nz,
+		 double xmin, double xmax,
+		 double ymin, double ymax,
+		 double zmin, double zmax) :
+  ThreeDGrid(nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax)
+{
+  // Set knots
+  coord["X"].resize(nx);
+  coord["Y"].resize(ny);
+  coord["Z"].resize(nz);
+
+  for (int i=0; i<nx; i++) coord["X"][i] = xmin + (xmax - xmin)*i/(nx-1);
+  for (int i=0; i<ny; i++) coord["Y"][i] = ymin + (ymax - ymin)*i/(ny-1);
+
+  if (nz>1) {
+    for (int i=0; i<nz; i++) coord["Z"][i] = zmin + (zmax - zmin)*i/(nz-1);
+  } else {
+    coord["Z"][0] = 0.0;
+  } 
+}
+
+void VtkGrid::replaceAll(std::string& str,
+			 const std::string& from,
+			 const std::string& to)
+{
+  // Sanity check: nothing to replace
+  if (from.empty()) return;
+
+  // Look for strings to replace
+  size_t start_pos = 0;
+  while((start_pos = str.find(from, start_pos)) != std::string::npos) {
+    str.replace(start_pos, from.length(), to);
+    start_pos += to.length();	// In case 'to' contains 'from'
+  }
+}
+
+
+void VtkGrid::Add(const std::vector<double>& data, const std::string& name)
+{
+  // Remove XML sensitive characters; paraview bug?
+  std::string newName(name);
+  replaceAll(newName, ">", ".gt.");
+  replaceAll(newName, "<", ".lt.");
+
+  dataSet[newName].resize(nx*ny*nz);
+
+  // Insert grid data
+  //
+  int I = 0;
+
+  for (int i=0; i<nx; i++) {
+    for (int j=0; j<ny; j++) {
+      for (int k=0; k<nz; k++) {
+	dataSet[newName][I++] = static_cast<float>(data[(k*ny + j)*nx + i]);
+      }
+    }
+  }
+
+}
+
+
+void VtkGrid::writeBeg(std::ofstream & fout)
+{
+  int zero = 0;
+ fout << "<?xml version=\"1.0\"?>" << std::endl;
+  fout << "<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"LittleEndian\" header_type=\"UInt32\">" << std::endl;
+  fout << "  <RectilinearGrid WholeExtent=\""
+       << zero << " " << nx-1 << " "
+       << zero << " " << ny-1 << " "
+       << zero << " " << nz-1
+       << "\">" << std::endl;
+  fout << "  <Piece Extent=\""
+       << zero << " " << nx-1 << " "
+       << zero << " " << ny-1 << " "
+       << zero << " " << nz-1
+       << "\">" << std::endl;
+}
+
+void VtkGrid::writeFields(std::ofstream & fout)
+{
+ fout << "    <PointData>" << std::endl;
+  for (auto v : dataSet) {
+    // Get ranges
+    auto vmin = std::min_element(v.second.begin(), v.second.end());
+    auto vmax = std::max_element(v.second.begin(), v.second.end());
+    fout << "      <DataArray type=\"Float32\" Name=\"" << v.first << "\" "
+	 << " format=\"ascii\" RangeMin=\"" << *vmin << "\" RangeMax=\"" << *vmax
+	 << "\">" << std::endl;
+    int cnt = 0;
+    for (auto & d : v.second) {
+      if (cnt % 6 == 0) fout << "         ";
+      fout << std::scientific << std::setprecision(8) << d << " ";
+      if (++cnt % 6 == 0) fout << std::endl;
+    }
+    if (cnt % 6) fout << std::endl;
+    fout << "      </DataArray>" << std::endl;
+  }
+  fout << "    </PointData>" << std::endl;
+
+  // No cell data here so this stanza is empty, but it's part of the spec
+  fout << "    <CellData>" << std::endl;
+  fout << "    </CellData>" << std::endl;
+}
+
+void VtkGrid::writeCoords(std::ofstream & fout)
+{
+  fout << "    <Coordinates>" << std::endl;
+  for (auto & c : coord) {
+    fout << "      <DataArray type=\"Float32\" Name=\""
+	 << c.first << "\" format=\"ascii\" RangeMin=\""
+	 << *std::min_element(c.second.begin(), c.second.end())
+       << "\" RangeMax=\""
+	 << *std::max_element(c.second.begin(), c.second.end()) << "\">"
+	 << std::endl;
+
+    int cnt = 0;
+    for (auto & v : c.second) {
+      if (cnt % 6 == 0) fout << "        ";
+      fout << std::scientific << std::setprecision(8) << v << " ";
+      if (++cnt % 6 == 0) fout << std::endl;
+    }
+    if (cnt % 6) fout << std::endl;
+    fout << "      </DataArray>" << std::endl;
+  }
+  fout << "    </Coordinates>" << std::endl;
+}
+
+void VtkGrid::writeEnd(std::ofstream & fout)
+{
+  fout << "  </Piece>" << std::endl;
+  fout << "  </RectilinearGrid>" << std::endl;
+  fout << "</VTKFile>" << std::endl;
+}
+
+void VtkGrid::Write(const std::string& name)
+{
+  // Create the filename with the correct extension for Paraview
+  std::ostringstream filename;
+  filename << name << ".vtr";
+
+  std::ofstream fout(filename.str());
+  if (fout) {
+
+  } else {
+    throw std::runtime_error
+      (std::string("VtkGrid::Write: could not open file <") + filename.str() + ">");
+  }
+
+  writeBeg(fout);		// VTK XML preamble; open stanzas
+  writeFields(fout);		// Write each data field in the dataset map
+  writeCoords(fout); 		// Write the coordinate grids
+  writeEnd(fout);		// VTK XML finish; close open stanzas
+}
+
+#endif

include/DataGrid.H

@@ -8,10 +8,7 @@
 #include <config_exp.h>
 
 #include <TableGrid.H>
-#ifdef HAVE_VTK
 #include <VtkGrid.H>
-#endif
-
 
 /**
    This implementation of ThreeDGrid instantiates at VtkGrid if VTK is
@@ -32,19 +29,22 @@ public:
   DataGrid(int nx, int ny, int nz,
 	   double xmin, double xmax,
 	   double ymin, double ymax,
-	   double zmin, double zmax)
+	   double zmin, double zmax,
+	   std::string type="VTK")
   {
-#ifdef HAVE_VTK
-    ptr = std::make_shared<VtkGrid>(nx, ny, nz,
-				    xmin, xmax,
-				    ymin, ymax,
-				    zmin, zmax);
-#else    
-    ptr = std::make_shared<TableGrid>(nx, ny, nz,
+    std::transform(type.begin(), type.end(), type.begin(), 
+                   [](unsigned char c){ return std::tolower(c); });
+
+    if (type.compare("vtk") == 0)
+      ptr = std::make_shared<VtkGrid>(nx, ny, nz,
 				      xmin, xmax,
 				      ymin, ymax,
 				      zmin, zmax);
-#endif
+    else
+      ptr = std::make_shared<TableGrid>(nx, ny, nz,
+					xmin, xmax,
+					ymin, ymax,
+					zmin, zmax);
   }
 
   //! Add data

include/VtkGrid.H

@@ -1,10 +1,15 @@
 #ifndef _VtkGrid_H
 #define _VtkGrid_H
 
-#include <vector>
+#include <algorithm>
+#include <fstream>
 #include <sstream>
+#include <vector>
 #include <memory>
 #include <limits>
+#include <map>
+
+#ifdef HAVE_VTK
 
 //
 // VTK stuff
@@ -55,6 +60,56 @@ public:
   void Write(const std::string& name);
 };
 
+#else
+
+#include <ThreeDGrid.H>
+
+/**
+   A ThreeDGrid implementation of the rectangular grid database
+   implement in VTK and designed for consumption by PyVTK or ParaView.
+ */
+class VtkGrid : public ThreeDGrid
+{
+private:
+  // Grid coordinates
+  std::map<std::string, std::vector<float>> coord;
+
+  // Dataset
+  std::map<std::string, std::vector<float>> dataSet;
+
+  // Write header
+  void writeBeg(std::ofstream& file);
+
+  // Write footer
+  void writeEnd(std::ofstream& file);
+
+  // Write scalar fields
+  void writeFields(std::ofstream& file);
+
+  // Write coordinates
+  void writeCoords(std::ofstream& file);
+
+  // String replacement utility
+  void replaceAll(std::string& str,
+		  const std::string& from,
+		  const std::string& to);
+
+public:
+  //! Constructor
+  VtkGrid(int nx, int ny, int nz,
+	 double xmin, double xmax,
+	 double ymin, double ymax,
+	 double zmin, double zmax);
+
+  //! Add data for two-dimensional cylindrical basis
+  void Add(const std::vector<double>& data, const std::string& name);
+
+  //! Write output file
+  void Write(const std::string& name);
+};
+
+#endif
+
 typedef std::shared_ptr<VtkGrid> VtkGridPtr;
 
 #endif