add chunkloop template

src/WeightingFieldCalculator.cxx

@@ -11,15 +11,88 @@ WeightingFieldCalculator::WeightingFieldCalculator(CylinderGeometry& geom, const
   antenna.AddToGeometry(*f, geom);
 }
 
+void saving_chunkloop(meep::fields_chunk *fc, int ichunk, meep::component cgrid, meep::ivec is, meep::ivec ie,
+		      meep::vec s0, meep::vec s1, meep::vec e0, meep::vec e1, double dV0, double dV1,
+		      meep::ivec shift, std::complex<double> shift_phase,
+		      const meep::symmetry &S, int sn, void *chunkloop_data)
+{
+
+  // make sure to give enough flexibility to store different fields / field combinations  
+  std::cout << "creating file" << std::endl;
+  std::string out_dir((char*)chunkloop_data);  
+  std::cout << "out_dir = " << out_dir << std::endl; 
+
+  // index vectors for start and end of chunk
+  meep::ivec isS = S.transform(is, sn) + shift;
+  meep::ivec ieS = S.transform(ie, sn) + shift;
+
+  size_t bufsz = 1;
+  size_t dims[2];
+  size_t stride[3] = {0, 1, 0};
+  int rank = 0;
+  LOOP_OVER_DIRECTIONS(fc->gv.dim, d) {
+    size_t cur_len = std::max(0, (ie.in_direction(d) - is.in_direction(d)) / 2 + 1);
+    bufsz *= cur_len;
+    dims[rank++] = cur_len;
+  }
+  stride[2] = dims[0];
+
+  std::cout << "bufsz = " << bufsz << std::endl;  
+  std::cout << "rank = " << rank << std::endl;
+  std::cout << "dims[0] = " << dims[0] << std::endl;
+  std::cout << "dims[1] = " << dims[1] << std::endl;
+
+  std::cout << "stride[0] = " << stride[0] << std::endl;
+  std::cout << "stride[1] = " << stride[1] << std::endl;
+  std::cout << "stride[2] = " << stride[2] << std::endl;
+
+  double buff[bufsz];
+
+  int metabuff[2];
+  metabuff[0] = is.in_direction(meep::R);
+  metabuff[1] = is.in_direction(meep::Z);
+
+  // some preliminary setup
+  meep::vec rshift(shift * (0.5*fc->gv.inva));  // shift into unit cell for PBC geometries
+
+  // prepare the list of field components to fetch at each grid point
+  meep::component components[] = {meep::Ex, meep::Ey, meep::Ez};
+  meep::chunkloop_field_components data(fc, cgrid, shift_phase, S, sn, 3, components);
+
+  // loop over all grid points in chunk
+  LOOP_OVER_IVECS(fc->gv, is, ie, idx) {
+    // get grid indices and coordinates of parent point
+    IVEC_LOOP_ILOC(fc->gv, iparent);  // grid indices
+    IVEC_LOOP_LOC(fc->gv, rparent);   // cartesian coordinates
+
+    // apply symmetry transform to get grid indices and coordinates of child point
+    meep::ivec ichild = S.transform(iparent, sn) + shift;
+    meep::vec rchild = S.transform(rparent, sn) + rshift;
+
+    // fetch field components at child point
+    data.update_values(idx);
+    std::complex<double> Ex_val = data.values[0];
+    std::complex<double> Ey_val = data.values[1];
+    std::complex<double> Ez_val = data.values[2];
+
+    ptrdiff_t idx2 = loop_i1 * stride[0] + loop_i2 * stride[2] + loop_i3 * stride[1];
+    buff[idx2] = Ez_val.real();
+  }
+
+  std::cout << "writing chunk ...";
+}
+
 void WeightingFieldCalculator::Calculate(double t_end) {
 
+  std::string data = "testdata";
   f -> output_hdf5(meep::Dielectric, gv -> surroundings());
 
   while (f -> time() < t_end) {
     if(meep::am_master()) {
       std::cout << "Simulation time: " << f -> time() << std::endl;
     }
     f -> step();
+    f -> loop_in_chunks(saving_chunkloop, (void*)data.c_str(), f -> total_volume());
     f -> output_hdf5(meep::Ez, gv -> surroundings());
   }
 }