Revised array metadata

python/meep.i

@@ -408,7 +408,7 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c,
 template<typename dft_type>
 PyObject *_get_dft_array(meep::fields *f, dft_type dft, meep::component c, int num_freq) {
     int rank;
-    int dims[3];
+    size_t dims[3];
     std::complex<double> *dft_arr = f->get_dft_array(dft, c, num_freq, &rank, dims);
 
     if (rank==0 || dft_arr==NULL) // this can happen e.g. if component c vanishes by symmetry
@@ -417,7 +417,7 @@ PyObject *_get_dft_array(meep::fields *f, dft_type dft, meep::component c, int n
     size_t length = 1;
     npy_intp *arr_dims = new npy_intp[rank];
     for (int i = 0; i < rank; ++i) {
-        arr_dims[i] = dims[i];
+        arr_dims[i] = dims[i];       // implicit size_t -> int cast, presumed safe for individual array dimensions
         length *= dims[i];
     }
 
@@ -496,9 +496,9 @@ kpoint_list get_eigenmode_coefficients_and_kpoints(meep::fields *f, meep::dft_fl
 }
 
 PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where, size_t dims[3],
-                                      bool collapse_empty_dimensions) {
+                                      bool collapse_empty_dimensions, bool snap_empty_dimensions) {
     meep::direction dirs[3] = {meep::X, meep::X, meep::X};
-    int rank = f->get_array_slice_dimensions(where, dims, dirs, collapse_empty_dimensions);
+    int rank = f->get_array_slice_dimensions(where, dims, dirs, collapse_empty_dimensions, snap_empty_dimensions);
 
     PyObject *py_dirs = PyList_New(3);
     for (Py_ssize_t i = 0; i < 3; ++i) {
@@ -1343,7 +1343,7 @@ kpoint_list get_eigenmode_coefficients_and_kpoints(meep::fields *f, meep::dft_fl
                                                    meep::kpoint_func user_kpoint_func, void *user_kpoint_data,
                                                    bool verbose);
 PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where, size_t dims[3],
-                                      bool collapse_empty_dimensions);
+                                      bool collapse_empty_dimensions, , bool snap_empty_dimensions);
 
 %ignore eps_func;
 %ignore inveps_func;

python/simulation.py

@@ -1706,7 +1706,7 @@ def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None
         else:
             v = self._volume_from_kwargs(vol, center, size)
 
-        _, dirs = mp._get_array_slice_dimensions(self.fields, v, dim_sizes, False)
+        _, dirs = mp._get_array_slice_dimensions(self.fields, v, dim_sizes, False, True)
 
         dims = [s for s in dim_sizes if s != 0]
 
@@ -1757,42 +1757,41 @@ def get_source(self, component, vol=None, center=None, size=None):
         else:
             v = self._volume_from_kwargs(vol, center, size)
         dim_sizes = np.zeros(3, dtype=np.uintp)
-        mp._get_array_slice_dimensions(self.fields, v, dim_sizes, False)
+        mp._get_array_slice_dimensions(self.fields, v, dim_sizes, False, True)
         dims = [s for s in dim_sizes if s != 0]
         arr = np.zeros(dims, dtype=np.complex128)
         self.fields.get_source_slice(v, component ,arr)
         return arr
 
-    def get_source_slice(self, component, vol=None, center=None, size=None):
-        warnings.warn('get_source_slice is deprecated. Please use get_source instead', DeprecationWarning)
-        return self.get_source(component, vol=vol, center=center, size=size)
-
-    def get_array_metadata(self, vol=None, center=None, size=None, dft=None):
-        if dft:
-            vol = dft.where
+    # if return_pw, the return value is (points, weights) where points is a
+    # mp.Vector3-valued array of the same dimensions as weights.
+    # otherwise return value is 4-tuple (xtics, ytics, ztics, weights).
+    def get_array_metadata(self, vol=None, center=None, size=None, dft_cell=None,
+                           collapse=False, snap=False, return_pw=False):
+        if dft_cell:
+            vol, collapse = dft_cell.where, True
         if vol is None and center is None and size is None:
             v = self.fields.total_volume()
         else:
             v = self._volume_from_kwargs(vol, center, size)
-        dims = np.zeros(3, dtype=np.uintp)
-        rank, dirs = mp._get_array_slice_dimensions(self.fields, v, dims, bool(dft))
-
-        nxyz = np.ones(3, dtype=np.intp)
-        for r in range(0,rank):
-            nxyz[dirs[r]]=dims[r]
-        nw=nxyz[0]*nxyz[1]*nxyz[2]
-        xtics=np.zeros(nxyz[0],dtype=np.float64)
-        ytics=np.zeros(nxyz[1],dtype=np.float64)
-        ztics=np.zeros(nxyz[2],dtype=np.float64)
-        weights=np.zeros(nw,dtype=np.float64)
-        self.fields.get_array_metadata(v, xtics, ytics, ztics, weights, bool(dft))
-        return (xtics,ytics,ztics,np.reshape(weights,dims[np.nonzero(dims)]))
+        xyzw_vector = self.fields.get_array_metadata(v, collapse, snap)
+        offset, tics = 0, []
+        for n in range(3):
+            N = int(xyzw_vector[offset])
+            tics.append( xyzw_vector[offset+1:offset+1+N] )
+            offset += 1+N
+        wshape=[len(t) for t in tics if len(t)>1]
+        weights=np.reshape(xyzw_vector[offset:],wshape)
+        if return_pw:
+            points=[ mp.Vector3(x,y,z) for x in tics[0] for y in tics[1] for z in tics[2] ]
+            return points,weights
+        return tics + [weights]
 
     def get_dft_array_metadata(self, dft_cell=None, vol=None, center=None, size=None):
         warnings.warn('get_dft_array_metadata is deprecated. Please use get_array_metadata instead',
                       DeprecationWarning)
         return self.get_array_metadata(vol=dft_cell.where if dft_cell is not None else vol,
-                                       center=center, size=size)
+                                       center=center, size=size, collapse=True)
 
     def get_eigenmode_coefficients(self, flux, bands, eig_parity=mp.NO_PARITY, eig_vol=None,
                                    eig_resolution=0, eig_tolerance=1e-12, kpoint_func=None, verbose=False):

python/tests/simulation.py

@@ -599,7 +599,7 @@ def print_field(sim):
     def test_source_slice(self):
         sim = self.init_simple_simulation()
         sim.run(until=5)
-        slice = sim.get_source_slice(mp.Ez)
+        slice = sim.get_source(mp.Ez)
         print(slice)
 
 if __name__ == '__main__':

src/array_slice.cpp

@@ -46,6 +46,14 @@ typedef struct {
   direction ds[3];
   size_t slice_size;
 
+  // if non-null, min_max_loc[0,1] are filled in by get_array_slice_dimensions_chunkloop
+  // with the (coordinate-wise) minimum and maximum grid points encountered
+  // in looping over the slice region.
+  vec *min_max_loc;
+
+  // if the component of the field array is NO_COMPONENT,
+  // the array is populated with the grid weights
+
   // the function to output and related info (offsets for averaging, etc.)
   // note: either fun *or* rfun should be non-NULL (not both)
   field_function fun;
@@ -106,13 +114,23 @@ static void get_array_slice_dimensions_chunkloop(fields_chunk *fc, int ichnk, co
   UNUSED(dV0);
   UNUSED(dV1);
   UNUSED(shift_phase);
-  UNUSED(fc);
+
   array_slice_data *data = (array_slice_data *)data_;
   ivec isS = S.transform(is, sn) + shift;
   ivec ieS = S.transform(ie, sn) + shift;
   data->min_corner = min(data->min_corner, min(isS, ieS));
   data->max_corner = max(data->max_corner, max(isS, ieS));
   data->num_chunks++;
+
+  if (data->min_max_loc) {
+     vec rshift(shift * (0.5 * fc->gv.inva));
+     LOOP_OVER_IVECS(fc->gv, is, ie, idx) {
+       IVEC_LOOP_LOC(fc->gv, loc);
+       loc = S.transform(loc, sn) + rshift;
+       data->min_max_loc[0] = min(loc, data->min_max_loc[0]);
+       data->min_max_loc[1] = max(loc, data->min_max_loc[1]);
+      }
+   }
 }
 
 /*****************************************************************/
@@ -186,6 +204,7 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
 
   //-----------------------------------------------------------------------//
   // Find output chunk dimensions and strides, etc.
+  //-----------------------------------------------------------------------//
 
   int start[3] = {0, 0, 0}, count[3] = {1, 1, 1};
   ptrdiff_t offset[3] = {0, 0, 0};
@@ -283,7 +302,11 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
     // interpolate fields at the four nearest grid points
     // to get the value of the field component for this point
     for (int i = 0; i < num_components; ++i) {
-      if (cS[i] == Dielectric) {
+      // special case for fetching grid point coordinates and weights
+      if (cS[i] == NO_COMPONENT) {
+          fields[i] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2);
+       }
+      else if (cS[i] == Dielectric) {
         double tr = 0.0;
         for (int k = 0; k < data->ninveps; ++k) {
           const realnum *ie = fc->s->chi1inv[iecs[k]][ieds[k]];
@@ -344,7 +367,9 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
 /* get_array_slice.                                            */
 /***************************************************************/
 int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], direction dirs[3],
-                                       bool collapse_empty_dimensions, void *caller_data) {
+                                       bool collapse_empty_dimensions,
+                                       bool snap_empty_dimensions,
+                                       vec *min_max_loc, void *caller_data) {
   am_now_working_on(FieldOutput);
 
   // use a local data structure if the caller didn't provide one
@@ -356,10 +381,28 @@ int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], dire
   data->max_corner = gv.round_vec(where.get_min_corner()) - one_ivec(gv.dim);
   data->num_chunks = 0;
 
-  loop_in_chunks(get_array_slice_dimensions_chunkloop, (void *)data, where, Centered, true, true);
+  data->min_max_loc = min_max_loc;
+  vec *min_loc = 0, *max_loc = 0;
+  if (min_max_loc)
+   { min_loc = min_max_loc + 0;
+     max_loc = min_max_loc + 1;
+     LOOP_OVER_DIRECTIONS(gv.dim,d)  {
+       min_loc->set_direction(d,+infinity);
+       max_loc->set_direction(d,-infinity);
+     }
+    }
+
+  bool use_symmetry = true;
+  loop_in_chunks(get_array_slice_dimensions_chunkloop, (void *)data, where, Centered,
+                 use_symmetry, snap_empty_dimensions);
 
-  data->max_corner = max_to_all(data->max_corner);
   data->min_corner = -max_to_all(-data->min_corner); // i.e., min_to_all
+  data->max_corner = max_to_all(data->max_corner);
+  if (min_max_loc)
+   LOOP_OVER_DIRECTIONS(gv.dim,d) {
+       min_loc->set_direction(d, -1.0*max_to_all(-1.0*min_loc->in_direction(d)));
+       max_loc->set_direction(d,      max_to_all(     max_loc->in_direction(d)));
+    }
   data->num_chunks = sum_to_all(data->num_chunks);
   if (data->num_chunks == 0 || !(data->min_corner <= data->max_corner))
     return 0; // no data to write;
@@ -398,13 +441,13 @@ void *fields::do_get_array_slice(const volume &where, std::vector<component> com
   /* call get_array_slice_dimensions to get slice dimensions and */
   /* partially initialze an array_slice_data struct              */
   /***************************************************************/
-  // by tradition, empty dimensions in time-domain field arrays are *not* collapsed;
+  // by tradition, empty dimensions in time-domain field arrays are *not* collapsed
   // TODO make this a caller-specifiable parameter to get_array_slice()?
-  bool collapse = false;
+  bool collapse = false, snap = true;
   size_t dims[3];
   direction dirs[3];
   array_slice_data data;
-  int rank = get_array_slice_dimensions(where, dims, dirs, collapse, &data);
+  int rank = get_array_slice_dimensions(where, dims, dirs, collapse, snap, 0, &data);
   size_t slice_size = data.slice_size;
   if (rank == 0 || slice_size == 0) return 0; // no data to write
 
@@ -551,4 +594,140 @@ cdouble *fields::get_source_slice(const volume &where, component source_slice_co
                                        true);
 }
 
+/**********************************************************************/
+/* increment a multi-index <n_1, n_2, ..., n_r> in which n_i runs over*/
+/* the range 0 <= n_i < nMax_i. return true if this is the increment  */
+/* that brings the multiindex to all zeros.                           */
+/**********************************************************************/
+bool increment(size_t *n, size_t *nMax, int rank) {
+  for (rank--; rank >= 0; rank--)
+    if (++n[rank] < nMax[rank])
+      return false;
+    else
+      n[rank] = 0;
+  return true;
+}
+
+// data_size = 1,2 for real,complex-valued array
+double *collapse_array(double *array, int *rank, size_t dims[3], direction dirs[3], volume where, int data_size=1) {
+
+  /*--------------------------------------------------------------*/
+  /*- detect empty dimensions and compute rank and strides for    */
+  /*- collapsed array                                             */
+  /*--------------------------------------------------------------*/
+  int full_rank = *rank;
+  if (full_rank == 0) return array;
+
+  int reduced_rank = 0;
+  size_t reduced_dims[3], reduced_stride[3] = {1, 1, 1};
+  direction reduced_dirs[3];
+  for(int r=0; r<full_rank; r++) {
+    if (where.in_direction(dirs[r]) == 0.0)
+      reduced_stride[r] = 0; // degenerate dimension, to be collapsed
+    else
+      { reduced_dirs[reduced_rank]   = dirs[r];
+        reduced_dims[reduced_rank++] = dims[r];
+      }
+  }
+  if (reduced_rank == full_rank) return array; // nothing to collapse
+
+  /*--------------------------------------------------------------*/
+  /*- set up strides into full and reduced arrays                 */
+  /*--------------------------------------------------------------*/
+  size_t stride[3] = {1, 1, 1}; // non-reduced array strides
+  if (full_rank == 2)
+    stride[0] = dims[1]; // rstride is already all set in this case
+  else if (full_rank == 3) {
+    stride[0] = dims[1] * dims[2];
+    stride[1] = dims[2];
+    if (reduced_stride[0] != 0)
+      reduced_stride[0] = reduced_dims[1];
+    else if (reduced_stride[1] != 0)
+      reduced_stride[1] = reduced_dims[1];
+    // else: two degenerate dimensions->reduced array is 1-diml, no strides needed
+  }
+
+  /*--------------------------------------------------------------*/
+  /*- allocate reduced array and compress full array into it     -*/
+  /*--------------------------------------------------------------*/
+  size_t reduced_grid_size  = reduced_dims[0] * (reduced_rank == 2 ? reduced_dims[1] : 1);
+  size_t reduced_array_size = data_size * reduced_grid_size;
+  double *reduced_array     = new double[ reduced_array_size ];
+  if (!reduced_array) abort("%s:%i: out of memory (%zu)", __FILE__, __LINE__, reduced_array_size);
+  memset(reduced_array, 0, reduced_array_size * sizeof(double));
+
+  size_t n[3] = {0, 0, 0};
+  do {
+    size_t index = n[0] * stride[0] + n[1] * stride[1] + n[2] * stride[2];
+    size_t rindex = n[0] * reduced_stride[0] + n[1] * reduced_stride[1] + n[2] * reduced_stride[2];
+    for(int i=0; i<data_size; i++)
+     reduced_array[ data_size*rindex + i ] += array[ data_size*index + i ];
+  } while (!increment(n, dims, full_rank));
+
+  *rank = reduced_rank;
+  for(int r=0; r<reduced_rank; r++)
+   { dims[r] = reduced_dims[r];
+     dirs[r] = reduced_dirs[r];
+   }
+  delete[] array;
+  return reduced_array;
+}
+
+cdouble *collapse_array(cdouble *array, int *rank, size_t dims[3], direction dirs[3], volume where)
+{ return (cdouble *)collapse_array((double *)array, rank, dims, dirs, where, 2); }
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+std::vector<double> fields::get_array_metadata(const volume &where,
+                                               bool collapse_empty_dimensions,
+                                               bool snap_empty_dimensions) {
+
+  /* get extremal corners of subgrid and array of weights, collapsed if necessary */
+  size_t dims[3];
+  direction dirs[3];
+  vec min_max_loc[2];   // extremal points in subgrid
+  int rank = get_array_slice_dimensions(where, dims, dirs, false /*collapse_empty_dimensions*/,
+                                        snap_empty_dimensions, min_max_loc);
+  int full_rank=rank;
+  direction full_dirs[3];
+  for(int fr=0; fr<rank; fr++)
+   full_dirs[fr]=dirs[fr];
+
+  double *weights = get_array_slice(where, NO_COMPONENT);
+  if (collapse_empty_dimensions)
+   weights=collapse_array(weights, &rank, dims, dirs, where);
+
+  /* get length and endpoints of x,y,z tics arrays */
+  size_t nxyz[3]={1,1,1};
+  double xyzmin[3]={0.0,0.0,0.0}, xyzmax[3]={0.0,0.0,0.0};
+  for(int fr=0, rr=0; fr<full_rank; fr++)
+   { direction d = full_dirs[fr];
+     int nd = d-X;
+     if ( where.in_direction(d)==0.0 && collapse_empty_dimensions )
+      { xyzmin[nd]=xyzmax[nd]=where.in_direction_min(d);
+        nxyz[nd]=1;
+      }
+     else
+      { nxyz[nd]   = dims[rr++];
+        xyzmin[nd] = min_max_loc[0].in_direction(d);
+        xyzmax[nd] = min_max_loc[1].in_direction(d);
+      }
+   }
+
+  /* pack all data into a single vector with each tics array preceded by its */
+  /* length: [ NX, xtics[:], NY, ytics[:], NZ, ztics[:], weights[:] ]        */
+  std::vector<double> xyzw;
+  for(int nd=0; nd<3; nd++)
+   { xyzw.push_back( (double)nxyz[nd] );
+     for(size_t n=0; n<nxyz[nd]; n++)
+      xyzw.push_back( xyzmin[nd] + n*gv.inva );
+   }
+  for(unsigned nw=0; nw<(nxyz[0]*nxyz[1]*nxyz[2]); nw++)
+   xyzw.push_back(weights[nw]);
+
+  return xyzw;
+
+} // get_array_metadata
+
 } // namespace meep