collapse empty dimensions of time-domain grid slices by default and remove snap argument from get_array_metadata

doc/docs/Python_User_Interface.md

@@ -2457,8 +2457,8 @@ plt.show()
 plt.savefig('sim_domain.png')
 ```
 
-Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function must be called
-from all parallel processes, but only produces a plot on the `meep.am_master()` process.
+Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects to be called
+on all processes, but only generates a plot on the master process.
 
 **Parameters:**
 
@@ -3349,8 +3349,7 @@ def get_array(self,
               size=None,
               cmplx=None,
               arr=None,
-              frequency=0,
-              omega=0):
+              frequency=0):
 ```
 
 <div class="method_docstring" markdown="1">
@@ -3370,7 +3369,7 @@ current simulation time.
   supplying a `Volume`.
 
 + `component`: field/material component (i.e., `mp.Ex`, `mp.Hy`, `mp.Sz`,
-  `mp.Dielectric`, etc). Defaults to `mp.Ez`.
+  `mp.Dielectric`, etc). Defaults to `None`.
 
 + `cmplx`: `boolean`; if `True`, return complex-valued data otherwise return
   real-valued data (default).
@@ -3425,9 +3424,9 @@ Returns the Fourier-transformed fields as a NumPy array.
 + `dft_obj`: a `dft_flux`, `dft_force`, `dft_fields`, or `dft_near2far` object
   obtained from calling the appropriate `add` function (e.g., `mp.add_flux`).
 
-+ `component`: a field component (e.g., `mp.Ez`)
++ `component`: a field component (e.g., `mp.Ez`).
 
-+ `num_freq`: the index of the frequency: an integer in the range `0...nfreq-1`,
++ `num_freq`: the index of the frequency. An integer in the range `0...nfreq-1`,
   where `nfreq` is the number of frequencies stored in `dft_obj` as set by the
   `nfreq` parameter to `add_dft_fields`, `add_dft_flux`, etc.
 
@@ -3449,8 +3448,6 @@ def get_array_metadata(self,
                        center=None,
                        size=None,
                        dft_cell=None,
-                       collapse=False,
-                       snap=False,
                        return_pw=False):
 ```
 
@@ -3460,9 +3457,9 @@ This routine provides geometric information useful for interpreting the arrays
 returned by `get_array` or `get_dft_array` for the spatial region defined by `vol`
 or `center`/`size`. In both cases, the return value is a tuple `(x,y,z,w)`, where:
 
-+ `x,y,z` are 1d NumPy arrays storing the $x,y,z$ coordinates of the points in the
++ `x,y,z` are tuples storing the $x,y,z$ coordinates of the points in the
   grid slice
-+ `w` is an array of the same dimensions as the array returned by
++ `w` is a NumPy array of the same dimensions as the array returned by
   `get_array`/`get_dft_array`, whose entries are the weights in a cubature rule
   for integrating over the spatial region (with the points in the cubature rule
   being just the grid points contained in the region). Thus, if $Q(\mathbf{x})$ is
@@ -3472,15 +3469,23 @@ or `center`/`size`. In both cases, the return value is a tuple `(x,y,z,w)`, wher
 $$ \int_{\mathcal V} Q(\mathbf{x})d\mathbf{x} \approx \sum_{n} w_n Q_n.$$
 
 This is a 1-, 2-, or 3-dimensional integral depending on the number of dimensions
-in which $\mathcal{V}$ has zero extent. If the $\{Q_n\}$ samples are stored in an
+in which $\mathcal{V}$ has zero extent. If the $Q_n$ samples are stored in an
 array `Q` of the same dimensions as `w`, then evaluating the sum on the RHS is
 just one line: `np.sum(w*Q).`
 
 A convenience parameter `dft_cell` is provided as an alternative to `vol` or
-`center/size`; set `dft_cell` to a `dft_flux` or `dft_fields` object to define the
-region covered by the array. If the `dft` argument is provided then all other
-arguments (`vol`, `center`, and `size`) are ignored. If no arguments are provided,
-then the entire cell is used.
+`center`/`size`. Set `dft_cell` to a `dft_flux` or `dft_fields` object to define the
+region covered by the array. If the `dft_cell` argument is provided then all other
+arguments related to the spatial region (`vol`, `center`, and `size`) are ignored.
+If no arguments are provided, then the entire cell is used.
+
+For empty dimensions of the grid slice `get_array_metadata` will collapse
+the *two* elements corresponding to the nearest Yee grid points into a *single*
+element using linear interpolation.
+
+If `return_pw=True`, the return value is a 2-tuple `(p,w)` where `p` (points) is a
+list of `mp.Vector3`s with the same dimensions as `w` (weights). Otherwise, by
+default the return value is a 4-tuple `(x,y,z,w)`.
 
 </div>
 
@@ -3559,6 +3564,7 @@ def get_source(self,
 Return an array of complex values of the [source](#source) amplitude for
 `component` over the given `vol` or `center`/`size`. The array has the same
 dimensions as that returned by [`get_array`](#array-slices).
+Not supported for [cylindrical coordinates](Python_Tutorials/Cylindrical_Coordinates.md).
 
 </div>
 

python/meep.i

@@ -545,7 +545,7 @@ 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 snap_empty_dimensions,
+                                      bool collapse_empty_dimensions,
                                       meep::component cgrid = Centered, PyObject *min_max_loc = NULL) {
     // Return value: New reference
     meep::direction dirs[3] = {meep::X, meep::X, meep::X};
@@ -554,7 +554,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where
     meep::vec* min_max_loc_vec_ptr = min_max_loc_vec;
     if (!min_max_loc) min_max_loc_vec_ptr = NULL;
 
-    int rank = f->get_array_slice_dimensions(where, dims, dirs, collapse_empty_dimensions, snap_empty_dimensions, min_max_loc_vec_ptr, 0, cgrid);
+    int rank = f->get_array_slice_dimensions(where, dims, dirs, collapse_empty_dimensions, min_max_loc_vec_ptr, 0, cgrid);
 
     PyObject *py_dirs = PyList_New(3);
     for (Py_ssize_t i = 0; i < 3; ++i) {
@@ -1558,7 +1558,7 @@ kpoint_list get_eigenmode_coefficients_and_kpoints(meep::fields *f, meep::dft_fl
                                                    meep::kpoint_func user_kpoint_func, void *user_kpoint_data,
                                                    double *cscale, meep::direction d);
 PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where, size_t dims[3],
-                                      bool collapse_empty_dimensions, bool snap_empty_dimensions,
+                                      bool collapse_empty_dimensions,
                                       meep::component cgrid = Centered, PyObject *min_max_loc = NULL);
 
 %ignore eps_func;

python/simulation.py

@@ -3094,7 +3094,7 @@ def h5topng(self, rm_h5, option, *step_funcs):
         cmd = re.sub(r'\$EPS', self.last_eps_filename, opts)
         return convert_h5(rm_h5, cmd, *step_funcs)
 
-    def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None, arr=None, frequency=0, omega=0):
+    def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None, arr=None, frequency=0):
         """
         Takes as input a subregion of the cell and the field/material component. The
         method returns a NumPy array containing values of the field/material at the
@@ -3111,7 +3111,7 @@ def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None
           supplying a `Volume`.
 
         + `component`: field/material component (i.e., `mp.Ex`, `mp.Hy`, `mp.Sz`,
-          `mp.Dielectric`, etc). Defaults to `mp.Ez`.
+          `mp.Dielectric`, etc). Defaults to `None`.
 
         + `cmplx`: `boolean`; if `True`, return complex-valued data otherwise return
           real-valued data (default).
@@ -3157,14 +3157,10 @@ 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, True)
+        _, dirs = mp._get_array_slice_dimensions(self.fields, v, dim_sizes, True)
 
         dims = [s for s in dim_sizes if s != 0]
 
-        if omega != 0:
-            frequency = omega
-            warnings.warn("get_array: omega has been deprecated; use frequency instead", RuntimeWarning)
-
         if cmplx is None:
             cmplx = frequency != 0 or (component < mp.Dielectric and not self.fields.is_real)
 
@@ -3198,9 +3194,9 @@ def get_dft_array(self, dft_obj, component, num_freq):
         + `dft_obj`: a `dft_flux`, `dft_force`, `dft_fields`, or `dft_near2far` object
           obtained from calling the appropriate `add` function (e.g., `mp.add_flux`).
 
-        + `component`: a field component (e.g., `mp.Ez`)
+        + `component`: a field component (e.g., `mp.Ez`).
 
-        + `num_freq`: the index of the frequency: an integer in the range `0...nfreq-1`,
+        + `num_freq`: the index of the frequency. An integer in the range `0...nfreq-1`,
           where `nfreq` is the number of frequencies stored in `dft_obj` as set by the
           `nfreq` parameter to `add_dft_fields`, `add_dft_flux`, etc.
         """
@@ -3225,31 +3221,29 @@ def get_source(self, component, vol=None, center=None, size=None):
         Return an array of complex values of the [source](#source) amplitude for
         `component` over the given `vol` or `center`/`size`. The array has the same
         dimensions as that returned by [`get_array`](#array-slices).
+        Not supported for [cylindrical coordinates](Python_Tutorials/Cylindrical_Coordinates.md).
         """
         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)
         dim_sizes = np.zeros(3, dtype=np.uintp)
-        mp._get_array_slice_dimensions(self.fields, v, dim_sizes, False, True)
+        mp._get_array_slice_dimensions(self.fields, v, dim_sizes, 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)
+        self.fields.get_source_slice(v, component, arr)
         return arr
 
-    # 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):
+                           return_pw=False):
         """
         This routine provides geometric information useful for interpreting the arrays
         returned by `get_array` or `get_dft_array` for the spatial region defined by `vol`
         or `center`/`size`. In both cases, the return value is a tuple `(x,y,z,w)`, where:
 
-        + `x,y,z` are 1d NumPy arrays storing the $x,y,z$ coordinates of the points in the
+        + `x,y,z` are tuples storing the $x,y,z$ coordinates of the points in the
           grid slice
-        + `w` is an array of the same dimensions as the array returned by
+        + `w` is a NumPy array of the same dimensions as the array returned by
           `get_array`/`get_dft_array`, whose entries are the weights in a cubature rule
           for integrating over the spatial region (with the points in the cubature rule
           being just the grid points contained in the region). Thus, if $Q(\\mathbf{x})$ is
@@ -3259,41 +3253,43 @@ def get_array_metadata(self, vol=None, center=None, size=None, dft_cell=None,
         $$ \\int_{\\mathcal V} Q(\\mathbf{x})d\\mathbf{x} \\approx \\sum_{n} w_n Q_n.$$
 
         This is a 1-, 2-, or 3-dimensional integral depending on the number of dimensions
-        in which $\\mathcal{V}$ has zero extent. If the $\\{Q_n\\}$ samples are stored in an
+        in which $\\mathcal{V}$ has zero extent. If the $Q_n$ samples are stored in an
         array `Q` of the same dimensions as `w`, then evaluating the sum on the RHS is
         just one line: `np.sum(w*Q).`
 
         A convenience parameter `dft_cell` is provided as an alternative to `vol` or
-        `center/size`; set `dft_cell` to a `dft_flux` or `dft_fields` object to define the
-        region covered by the array. If the `dft` argument is provided then all other
-        arguments (`vol`, `center`, and `size`) are ignored. If no arguments are provided,
-        then the entire cell is used.
+        `center`/`size`. Set `dft_cell` to a `dft_flux` or `dft_fields` object to define the
+        region covered by the array. If the `dft_cell` argument is provided then all other
+        arguments related to the spatial region (`vol`, `center`, and `size`) are ignored.
+        If no arguments are provided, then the entire cell is used.
+
+        For empty dimensions of the grid slice `get_array_metadata` will collapse
+        the *two* elements corresponding to the nearest Yee grid points into a *single*
+        element using linear interpolation.
+
+        If `return_pw=True`, the return value is a 2-tuple `(p,w)` where `p` (points) is a
+        list of `mp.Vector3`s with the same dimensions as `w` (weights). Otherwise, by
+        default the return value is a 4-tuple `(x,y,z,w)`.
         """
         if dft_cell:
-            vol, collapse = dft_cell.where, True
+            vol = dft_cell.where
         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)
-        xyzw_vector = self.fields.get_array_metadata(v, collapse, snap)
+        xyzw_vector = self.fields.get_array_metadata(v)
         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)
+        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, collapse=True)
-
     def get_array_slice_dimensions(self, component, vol=None, center=None, size=None):
         """
         Computes the dimensions of a dft array for a particular `component` (`mp.Ez`, `mp.Ey`, etc.).
@@ -3310,7 +3306,7 @@ def get_array_slice_dimensions(self, component, vol=None, center=None, size=None
             v = self._volume_from_kwargs(vol, center, size)
         dim_sizes = np.zeros(3, dtype=np.uintp)
         corners = []
-        _,_ = mp._get_array_slice_dimensions(self.fields, v, dim_sizes, False, False, component,corners)
+        _,_ = mp._get_array_slice_dimensions(self.fields, v, dim_sizes, True, component, corners)
         dim_sizes[dim_sizes==0] = 1
         min_corner = corners[0]
         max_corner = corners[1]

python/tests/array_metadata.py

@@ -70,16 +70,13 @@ def vec_func(r):
                             boundary_layers=pml_layers)
 
         dft_obj = sim.add_dft_fields([mp.Ez], fcen, 0, 1, where=nonpml_vol)
-        dft_obj_yee = sim.add_dft_fields([mp.Hy], fcen, 0, 1, where=nonpml_vol,yee_grid=True)
         sim.run(until_after_sources=100)
 
-        Hy_yee_dims,_,_ = sim.get_array_slice_dimensions(mp.Hy, dft_obj.where)
-
         Ez = sim.get_dft_array(dft_obj, mp.Ez, 0)
         (X,Y,Z,W) = sim.get_array_metadata(dft_cell=dft_obj)
-        Eps = sim.get_array(vol=nonpml_vol,component=mp.Dielectric)
+        Eps = sim.get_array(vol=nonpml_vol, component=mp.Dielectric)
         EpsE2 = np.real(Eps*np.conj(Ez)*Ez)
-        xm, ym = np.meshgrid(X,Y)
+        xm, ym = np.meshgrid(X, Y)
         vec_func_sum = np.sum(W*(xm**2 + 2*ym**2))
         pulse_modal_volume = np.sum(W*EpsE2)/np.max(EpsE2) * vec_func_sum
 

python/tests/faraday_rotation.py

@@ -41,12 +41,11 @@ def check_rotation(self, mat, L, fsrc, zsrc, resolution, tmax, zout, kpred, tol=
                                  boundary_layers=pml_layers,
                                  default_material=mat, resolution=resolution)
 
-        record_vol = mp.Volume(center=mp.Vector3(0, 0, zout))
         record_Ex, record_Ey, record_t = [], [], []
 
         def record_ex_ey(sim):
-            record_Ex.append(sim.get_array(vol=record_vol, component=mp.Ex))
-            record_Ey.append(sim.get_array(vol=record_vol, component=mp.Ey))
+            record_Ex.append(sim.get_field_point(mp.Ex, mp.Vector3(0, 0, zout)))
+            record_Ey.append(sim.get_field_point(mp.Ey, mp.Vector3(0, 0, zout)))
             record_t.append(sim.meep_time())
 
         self.sim.run(mp.after_time(0.5*tmax, mp.at_every(1e-6, record_ex_ey)), until=tmax)