adding possibility to output 3D fieldmaps as VTK

src/Classic/Fields/Fieldmap.cpp

@@ -737,6 +737,76 @@ void Fieldmap::setFieldGap(double gap) {
 
 }
 
+void Fieldmap::write3DField(unsigned int nx,
+                            unsigned int ny,
+                            unsigned int nz,
+                            const std::pair<double, double> &xrange,
+                            const std::pair<double, double> &yrange,
+                            const std::pair<double, double> &zrange,
+                            const std::vector<Vector_t> &ef,
+                            const std::vector<Vector_t> &bf) {
+    const size_t numpoints = nx * ny * nz;
+    if (Ippl::myNode() != 0 ||
+        (ef.size() != numpoints && bf.size() != numpoints)) return;
+
+    size_t extensionStart = Filename_m.find_last_of('.');
+    std::ofstream of;
+    of.open (std::string ("data/" + Filename_m.substr(0,extensionStart) + ".vtk").c_str ());
+    assert (of.is_open ());
+    of.precision (6);
+
+    const double hx = (xrange.second - xrange.first) / (nx - 1);
+    const double hy = (yrange.second - yrange.first) / (ny - 1);
+    const double hz = (zrange.second - zrange.first) / (nz - 1);
+
+    of << "# vtk DataFile Version 2.0" << std::endl;
+    of << "generated by 3D fieldmaps" << std::endl;
+    of << "ASCII" << std::endl << std::endl;
+    of << "DATASET RECTILINEAR_GRID" << std::endl;
+    of << "DIMENSIONS " << nx << " " << ny << " " << nz << std::endl;
+
+    of << "X_COORDINATES " << nx << " float" << std::endl;
+    of << xrange.first;
+    for (unsigned int i = 1; i < nx - 1; ++ i) {
+        of << " " << xrange.first + i * hx;
+    }
+    of << " " << xrange.second << std::endl;
+
+    of << "Y_COORDINATES " << ny << " float" << std::endl;
+    of << yrange.first;
+    for (unsigned int i = 1; i < ny - 1; ++ i) {
+        of << " " << yrange.first + i * hy;
+    }
+    of << " " << yrange.second << std::endl;
+
+    of << "Z_COORDINATES " << nz << " float" << std::endl;
+    of << zrange.first;
+    for (unsigned int i = 1; i < nz - 1; ++ i) {
+        of << " " << zrange.first + i * hz;
+    }
+    of << " " << zrange.second << std::endl;
+
+    of << "POINT_DATA " << numpoints << std::endl;
+
+    if (ef.size() == numpoints) {
+        of << "VECTORS EField float" << std::endl;
+        // of << "LOOKUP_TABLE default" << std::endl;
+        for (size_t i = 0; i < numpoints; ++ i) {
+            of << ef[i](0) << " " << ef[i](1) << " " << ef[i](2) << std::endl;
+        }
+        // of << std::endl;
+    }
+
+    if (bf.size() == numpoints) {
+        of << "VECTORS BField float" << std::endl;
+        // of << "LOOKUP_TABLE default" << std::endl;
+        for (size_t i = 0; i < numpoints; ++ i) {
+            of << bf[i](0) << " " << bf[i](1) << " " << bf[i](2) << std::endl;
+        }
+        // of << std::endl;
+    }
+}
+
 REGISTER_PARSE_TYPE(int);
 REGISTER_PARSE_TYPE(unsigned int);
 REGISTER_PARSE_TYPE(double);

src/Classic/Fields/Fieldmap.h

@@ -148,6 +148,15 @@ protected:
                   gsl_spline *splineCoefficients,
                   gsl_interp_accel *splineAccelerator);
 
+    void write3DField(unsigned int nx,
+                      unsigned int ny,
+                      unsigned int nz,
+                      const std::pair<double, double> &xrange,
+                      const std::pair<double, double> &yrange,
+                      const std::pair<double, double> &zrange,
+                      const std::vector<Vector_t> &ef,
+                      const std::vector<Vector_t> &bf);
+
 public:
     virtual void readMap() = 0;
     virtual void freeMap() = 0;