diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index d284c3901..ece3b7eda 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -3349,7 +3349,8 @@ def get_array(self, size=None, cmplx=None, arr=None, - frequency=0): + frequency=0, + snap=False): ```
@@ -3381,7 +3382,15 @@ current simulation time. fetching the same slice repeatedly at different times). + `frequency`: optional frequency point over which the average eigenvalue of the - dielectric and permeability tensors are evaluated (defaults to 0). + $\varepsilon$ and $\mu$ tensors are evaluated (defaults to 0). + ++ `snap`: By default, the elements of the grid slice are obtained using a bilinear + interpolation of the nearest Yee grid points. Empty dimensions of the grid slice + are "collapsed" into a single element. However, if `snap` is set to `True`, this + interpolation behavior is disabled and the grid slice is instead "snapped" + everywhere to the nearest grid point. (Empty slice dimensions are still of size + one.) This feature is mainly useful for comparing results with the + [`output_` routines](#output-functions) (e.g., `output_epsilon`, `output_efield_z`, etc.). For convenience, the following wrappers for `get_array` over the entire cell are available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, @@ -3389,7 +3398,8 @@ available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, `get_Xfield_z()`, `get_Xfield_r()`, `get_Xfield_p()` where `X` is one of `h`, `b`, `e`, `d`, or `s`. The routines `get_Xfield_*` all return an array type consistent with the fields (real or complex). The routines `get_epsilon()` and `get_mu()` -accept the optional `frequency` parameter (defaults to 0). +accept the optional argument `frequency` (defaults to 0) and all routines accept +`snap` (defaults to `False`). **Note on array-slice dimensions:** The routines `get_epsilon`, `get_Xfield_z`, etc. use as default `size=meep.Simulation.fields.total_volume()` which for @@ -3457,7 +3467,7 @@ 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 tuples storing the $x,y,z$ coordinates of the points in the ++ `x,y,z` are 1d NumPy arrays storing the $x,y,z$ coordinates of the points in the grid slice + `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 diff --git a/python/simulation.py b/python/simulation.py index ad7f5f857..9a18ca246 100644 --- a/python/simulation.py +++ b/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): + def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None, arr=None, frequency=0, snap=False): """ 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 @@ -3123,7 +3123,15 @@ def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None fetching the same slice repeatedly at different times). + `frequency`: optional frequency point over which the average eigenvalue of the - dielectric and permeability tensors are evaluated (defaults to 0). + $\\varepsilon$ and $\\mu$ tensors are evaluated (defaults to 0). + + + `snap`: By default, the elements of the grid slice are obtained using a bilinear + interpolation of the nearest Yee grid points. Empty dimensions of the grid slice + are "collapsed" into a single element. However, if `snap` is set to `True`, this + interpolation behavior is disabled and the grid slice is instead "snapped" + everywhere to the nearest grid point. (Empty slice dimensions are still of size + one.) This feature is mainly useful for comparing results with the + [`output_` routines](#output-functions) (e.g., `output_epsilon`, `output_efield_z`, etc.). For convenience, the following wrappers for `get_array` over the entire cell are available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, @@ -3131,7 +3139,8 @@ def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None `get_Xfield_z()`, `get_Xfield_r()`, `get_Xfield_p()` where `X` is one of `h`, `b`, `e`, `d`, or `s`. The routines `get_Xfield_*` all return an array type consistent with the fields (real or complex). The routines `get_epsilon()` and `get_mu()` - accept the optional `frequency` parameter (defaults to 0). + accept the optional argument `frequency` (defaults to 0) and all routines accept + `snap` (defaults to `False`). **Note on array-slice dimensions:** The routines `get_epsilon`, `get_Xfield_z`, etc. use as default `size=meep.Simulation.fields.total_volume()` which for @@ -3179,9 +3188,9 @@ def get_array(self, component=None, vol=None, center=None, size=None, cmplx=None arr = np.zeros(dims, dtype=np.complex128 if cmplx else np.float64) if np.iscomplexobj(arr): - self.fields.get_complex_array_slice(v, component, arr, frequency) + self.fields.get_complex_array_slice(v, component, arr, frequency, snap) else: - self.fields.get_array_slice(v, component, arr, frequency) + self.fields.get_array_slice(v, component, arr, frequency, snap) return arr @@ -3241,7 +3250,7 @@ def get_array_metadata(self, vol=None, center=None, size=None, dft_cell=None, 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 tuples storing the $x,y,z$ coordinates of the points in the + + `x,y,z` are 1d NumPy arrays storing the $x,y,z$ coordinates of the points in the grid slice + `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 @@ -3288,7 +3297,7 @@ def get_array_metadata(self, vol=None, center=None, size=None, dft_cell=None, 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] + return tuple(tics) + (weights,) def get_array_slice_dimensions(self, component, vol=None, center=None, size=None): """ @@ -3671,153 +3680,150 @@ def output_times(self, fname): fname += '.csv' self.fields.output_times(fname) - def get_epsilon(self,frequency=0,omega=0): - if omega != 0: - frequency = omega - warnings.warn("get_epsilon: omega has been deprecated; use frequency instead", RuntimeWarning) - return self.get_array(component=mp.Dielectric,frequency=frequency) + def get_epsilon(self,frequency=0,snap=False): + return self.get_array(component=mp.Dielectric,frequency=frequency,snap=snap) - def get_mu(self): - return self.get_array(component=mp.Permeability) + def get_mu(self,frequency=0,snap=False): + return self.get_array(component=mp.Permeability,frequency=frequency,snap=snap) - def get_hpwr(self): - return self.get_array(component=mp.H_EnergyDensity) + def get_hpwr(self,snap=False): + return self.get_array(component=mp.H_EnergyDensity,snap=snap) - def get_dpwr(self): - return self.get_array(component=mp.D_EnergyDensity) + def get_dpwr(self,snap=False): + return self.get_array(component=mp.D_EnergyDensity,snap=snap) - def get_tot_pwr(self): - return self.get_array(component=mp.EnergyDensity) + def get_tot_pwr(self,snap=False): + return self.get_array(component=mp.EnergyDensity,snap=snap) - def get_hfield(self): + def get_hfield(self,snap=False): if self.is_cylindrical: - r = self.get_array(mp.Hr, cmplx=not self.fields.is_real) - p = self.get_array(mp.Hp, cmplx=not self.fields.is_real) + r = self.get_array(mp.Hr, cmplx=not self.fields.is_real, snap=snap) + p = self.get_array(mp.Hp, cmplx=not self.fields.is_real, snap=snap) return np.stack([r, p], axis=-1) else: - x = self.get_array(mp.Hx, cmplx=not self.fields.is_real) - y = self.get_array(mp.Hy, cmplx=not self.fields.is_real) - z = self.get_array(mp.Hz, cmplx=not self.fields.is_real) + x = self.get_array(mp.Hx, cmplx=not self.fields.is_real, snap=snap) + y = self.get_array(mp.Hy, cmplx=not self.fields.is_real, snap=snap) + z = self.get_array(mp.Hz, cmplx=not self.fields.is_real, snap=snap) return np.stack([x, y, z], axis=-1) - def get_hfield_x(self): - return self.get_array(mp.Hx, cmplx=not self.fields.is_real) + def get_hfield_x(self,snap=False): + return self.get_array(mp.Hx, cmplx=not self.fields.is_real, snap=snap) - def get_hfield_y(self): - return self.get_array(mp.Hy, cmplx=not self.fields.is_real) + def get_hfield_y(self,snap=False): + return self.get_array(mp.Hy, cmplx=not self.fields.is_real, snap=snap) - def get_hfield_z(self): - return self.get_array(mp.Hz, cmplx=not self.fields.is_real) + def get_hfield_z(self,snap=False): + return self.get_array(mp.Hz, cmplx=not self.fields.is_real, snap=snap) - def get_hfield_r(self): - return self.get_array(mp.Hr, cmplx=not self.fields.is_real) + def get_hfield_r(self,snap=False): + return self.get_array(mp.Hr, cmplx=not self.fields.is_real, snap=snap) - def get_hfield_p(self): - return self.get_array(mp.Hp, cmplx=not self.fields.is_real) + def get_hfield_p(self,snap=False): + return self.get_array(mp.Hp, cmplx=not self.fields.is_real, snap=snap) - def get_bfield(self): + def get_bfield(self,snap=False): if self.is_cylindrical: - r = self.get_array(mp.Br, cmplx=not self.fields.is_real) - p = self.get_array(mp.Bp, cmplx=not self.fields.is_real) + r = self.get_array(mp.Br, cmplx=not self.fields.is_real, snap=snap) + p = self.get_array(mp.Bp, cmplx=not self.fields.is_real, snap=snap) return np.stack([r, p], axis=-1) else: - x = self.get_array(mp.Bx, cmplx=not self.fields.is_real) - y = self.get_array(mp.By, cmplx=not self.fields.is_real) - z = self.get_array(mp.Bz, cmplx=not self.fields.is_real) + x = self.get_array(mp.Bx, cmplx=not self.fields.is_real, snap=snap) + y = self.get_array(mp.By, cmplx=not self.fields.is_real, snap=snap) + z = self.get_array(mp.Bz, cmplx=not self.fields.is_real, snap=snap) return np.stack([x, y, z], axis=-1) - def get_bfield_x(self): - return self.get_array(mp.Bx, cmplx=not self.fields.is_real) + def get_bfield_x(self,snap=False): + return self.get_array(mp.Bx, cmplx=not self.fields.is_real, snap=snap) - def get_bfield_y(self): - return self.get_array(mp.By, cmplx=not self.fields.is_real) + def get_bfield_y(self,snap=False): + return self.get_array(mp.By, cmplx=not self.fields.is_real, snap=snap) - def get_bfield_z(self): - return self.get_array(mp.Bz, cmplx=not self.fields.is_real) + def get_bfield_z(self,snap=False): + return self.get_array(mp.Bz, cmplx=not self.fields.is_real, snap=snap) - def get_bfield_r(self): - return self.get_array(mp.Br, cmplx=not self.fields.is_real) + def get_bfield_r(self,snap=False): + return self.get_array(mp.Br, cmplx=not self.fields.is_real, snap=snap) - def get_bfield_p(self): - return self.get_array(mp.Bp, cmplx=not self.fields.is_real) + def get_bfield_p(self,snap=False): + return self.get_array(mp.Bp, cmplx=not self.fields.is_real, snap=snap) - def get_efield(self): + def get_efield(self,snap=False): if self.is_cylindrical: - r = self.get_array(mp.Er, cmplx=not self.fields.is_real) - p = self.get_array(mp.Ep, cmplx=not self.fields.is_real) + r = self.get_array(mp.Er, cmplx=not self.fields.is_real, snap=snap) + p = self.get_array(mp.Ep, cmplx=not self.fields.is_real, snap=snap) return np.stack([r, p], axis=-1) else: - x = self.get_array(mp.Ex, cmplx=not self.fields.is_real) - y = self.get_array(mp.Ey, cmplx=not self.fields.is_real) - z = self.get_array(mp.Ez, cmplx=not self.fields.is_real) + x = self.get_array(mp.Ex, cmplx=not self.fields.is_real, snap=snap) + y = self.get_array(mp.Ey, cmplx=not self.fields.is_real, snap=snap) + z = self.get_array(mp.Ez, cmplx=not self.fields.is_real, snap=snap) return np.stack([x, y, z], axis=-1) - def get_efield_x(self): - return self.get_array(mp.Ex, cmplx=not self.fields.is_real) + def get_efield_x(self,snap=False): + return self.get_array(mp.Ex, cmplx=not self.fields.is_real, snap=snap) - def get_efield_y(self): - return self.get_array(mp.Ey, cmplx=not self.fields.is_real) + def get_efield_y(self,snap=False): + return self.get_array(mp.Ey, cmplx=not self.fields.is_real, snap=snap) - def get_efield_z(self): - return self.get_array(mp.Ez, cmplx=not self.fields.is_real) + def get_efield_z(self,snap=False): + return self.get_array(mp.Ez, cmplx=not self.fields.is_real, snap=snap) - def get_efield_r(self): - return self.get_array(mp.Er, cmplx=not self.fields.is_real) + def get_efield_r(self,snap=False): + return self.get_array(mp.Er, cmplx=not self.fields.is_real, snap=snap) - def get_efield_p(self): - return self.get_array(mp.Ep, cmplx=not self.fields.is_real) + def get_efield_p(self,snap=False): + return self.get_array(mp.Ep, cmplx=not self.fields.is_real, snap=snap) - def get_dfield(self): + def get_dfield(self,snap=False): if self.is_cylindrical: - r = self.get_array(mp.Dr, cmplx=not self.fields.is_real) - p = self.get_array(mp.Dp, cmplx=not self.fields.is_real) + r = self.get_array(mp.Dr, cmplx=not self.fields.is_real, snap=snap) + p = self.get_array(mp.Dp, cmplx=not self.fields.is_real, snap=snap) return np.stack([r, p], axis=-1) else: - x = self.get_array(mp.Dx, cmplx=not self.fields.is_real) - y = self.get_array(mp.Dy, cmplx=not self.fields.is_real) - z = self.get_array(mp.Dz, cmplx=not self.fields.is_real) + x = self.get_array(mp.Dx, cmplx=not self.fields.is_real, snap=snap) + y = self.get_array(mp.Dy, cmplx=not self.fields.is_real, snap=snap) + z = self.get_array(mp.Dz, cmplx=not self.fields.is_real, snap=snap) return np.stack([x, y, z], axis=-1) - def get_dfield_x(self): - return self.get_array(mp.Dx, cmplx=not self.fields.is_real) + def get_dfield_x(self,snap=False): + return self.get_array(mp.Dx, cmplx=not self.fields.is_real, snap=snap) - def get_dfield_y(self): - return self.get_array(mp.Dy, cmplx=not self.fields.is_real) + def get_dfield_y(self,snap=False): + return self.get_array(mp.Dy, cmplx=not self.fields.is_real, snap=snap) - def get_dfield_z(self): - return self.get_array(mp.Dz, cmplx=not self.fields.is_real) + def get_dfield_z(self,snap=False): + return self.get_array(mp.Dz, cmplx=not self.fields.is_real, snap=snap) - def get_dfield_r(self): - return self.get_array(mp.Dr, cmplx=not self.fields.is_real) + def get_dfield_r(self,snap=False): + return self.get_array(mp.Dr, cmplx=not self.fields.is_real, snap=snap) - def get_dfield_p(self): - return self.get_array(mp.Dp, cmplx=not self.fields.is_real) + def get_dfield_p(self,snap=False): + return self.get_array(mp.Dp, cmplx=not self.fields.is_real, snap=snap) - def get_sfield(self): + def get_sfield(self,snap=False): if self.is_cylindrical: - r = self.get_array(mp.Sr, cmplx=not self.fields.is_real) - p = self.get_array(mp.Sp, cmplx=not self.fields.is_real) + r = self.get_array(mp.Sr, cmplx=not self.fields.is_real, snap=snap) + p = self.get_array(mp.Sp, cmplx=not self.fields.is_real, snap=snap) return np.stack([r, p], axis=-1) else: - x = self.get_array(mp.Sx, cmplx=not self.fields.is_real) - y = self.get_array(mp.Sy, cmplx=not self.fields.is_real) - z = self.get_array(mp.Sz, cmplx=not self.fields.is_real) + x = self.get_array(mp.Sx, cmplx=not self.fields.is_real, snap=snap) + y = self.get_array(mp.Sy, cmplx=not self.fields.is_real, snap=snap) + z = self.get_array(mp.Sz, cmplx=not self.fields.is_real, snap=snap) return np.stack([x, y, z], axis=-1) - def get_sfield_x(self): - return self.get_array(mp.Sx, cmplx=not self.fields.is_real) + def get_sfield_x(self,snap=False): + return self.get_array(mp.Sx, cmplx=not self.fields.is_real, snap=snap) - def get_sfield_y(self): - return self.get_array(mp.Sy, cmplx=not self.fields.is_real) + def get_sfield_y(self,snap=False): + return self.get_array(mp.Sy, cmplx=not self.fields.is_real, snap=snap) - def get_sfield_z(self): - return self.get_array(mp.Sz, cmplx=not self.fields.is_real) + def get_sfield_z(self,snap=False): + return self.get_array(mp.Sz, cmplx=not self.fields.is_real, snap=snap) - def get_sfield_r(self): - return self.get_array(mp.Sr, cmplx=not self.fields.is_real) + def get_sfield_r(self,snap=False): + return self.get_array(mp.Sr, cmplx=not self.fields.is_real, snap=snap) - def get_sfield_p(self): - return self.get_array(mp.Sp, cmplx=not self.fields.is_real) + def get_sfield_p(self,snap=False): + return self.get_array(mp.Sp, cmplx=not self.fields.is_real, snap=snap) def plot2D(self, ax=None, output_plane=None, fields=None, labels=False, diff --git a/python/tests/dispersive_eigenmode.py b/python/tests/dispersive_eigenmode.py index 60521b756..afe20b3f6 100644 --- a/python/tests/dispersive_eigenmode.py +++ b/python/tests/dispersive_eigenmode.py @@ -63,7 +63,7 @@ def verify_output_and_slice(self,material,frequency): sim.init_sim() # Check to make sure the get_slice routine is working with frequency - n_slice = np.sqrt(np.max(sim.get_epsilon(frequency))) + n_slice = np.sqrt(np.max(sim.get_epsilon(frequency=frequency,snap=True))) self.assertAlmostEqual(n,n_slice, places=4) # Check to make sure h5 output is working with frequency diff --git a/python/tests/simulation.py b/python/tests/simulation.py index 883ce050c..b874b3eae 100644 --- a/python/tests/simulation.py +++ b/python/tests/simulation.py @@ -373,10 +373,10 @@ def test_get_array_output(self): mp.output_tot_pwr(sim) mp.output_efield(sim) - eps_arr = sim.get_epsilon() - efield_z_arr = sim.get_efield_z() - energy_arr = sim.get_tot_pwr() - efield_arr = sim.get_efield() + eps_arr = sim.get_epsilon(snap=True) + efield_z_arr = sim.get_efield_z(snap=True) + energy_arr = sim.get_tot_pwr(snap=True) + efield_arr = sim.get_efield(snap=True) fname_fmt = os.path.join(self.temp_dir, 'test_get_array_output-{}-000020.00.h5') diff --git a/src/array_slice.cpp b/src/array_slice.cpp index 14d085be8..74f836fa0 100644 --- a/src/array_slice.cpp +++ b/src/array_slice.cpp @@ -76,6 +76,7 @@ typedef struct { component invmu_cs[3]; direction invmu_ds[3]; + bool snap; } array_slice_data; #define UNUSED(x) (void)x // silence compiler warnings @@ -332,7 +333,7 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr frequency)); if (abs(tr) == 0.0) tr += 4.0; // default inveps == 1 } - fields[i] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1.0) * (4.0 * data->ninveps) / tr; + fields[i] = (data->snap ? 1.0 : IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1.0)) * (4.0 * data->ninveps) / tr; } else if (cS[i] == Permeability) { complex tr(0.0, 0.0); @@ -344,7 +345,7 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr frequency)); if (abs(tr) == 0.0) tr += 4.0; // default invmu == 1 } - fields[i] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1.0) * (4.0 * data->ninvmu) / tr; + fields[i] = (data->snap ? 1.0 : IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1.0)) * (4.0 * data->ninvmu) / tr; } else { double f[2]; @@ -355,7 +356,7 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr fc->f[cS[i]][k][idx + off[2 * i] + off[2 * i + 1]]); else f[k] = 0; - fields[i] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1.0) * complex(f[0], f[1]) * ph[i]; + fields[i] = (data->snap ? 1.0 : IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1.0)) * complex(f[0], f[1]) * ph[i]; } } @@ -561,7 +562,7 @@ cdouble *collapse_array(cdouble *array, int *rank, size_t dims[3], direction dir /**********************************************************************/ void *fields::do_get_array_slice(const volume &where, std::vector components, field_function fun, field_rfunction rfun, void *fun_data, - void *vslice, double frequency) { + void *vslice, double frequency, bool snap) { am_now_working_on(FieldOutput); /***************************************************************/ @@ -573,7 +574,6 @@ void *fields::do_get_array_slice(const volume &where, std::vector com array_slice_data data; int rank = get_array_slice_dimensions(where, dims, dirs, false, 0, &data); size_t slice_size = data.slice_size; - bool complex_data = (rfun == 0); cdouble *zslice; double *slice; @@ -591,6 +591,7 @@ void *fields::do_get_array_slice(const volume &where, std::vector com } data.vslice = vslice_uncollapsed; + data.snap = snap; data.fun = fun; data.rfun = rfun; data.fun_data = fun_data; @@ -674,40 +675,42 @@ void *fields::do_get_array_slice(const volume &where, std::vector com /***************************************************************/ double *fields::get_array_slice(const volume &where, std::vector components, field_rfunction rfun, void *fun_data, double *slice, - double frequency) { + double frequency, bool snap) { return (double *)do_get_array_slice(where, components, 0, rfun, fun_data, (void *)slice, - frequency); + frequency, snap); } cdouble *fields::get_complex_array_slice(const volume &where, std::vector components, field_function fun, void *fun_data, cdouble *slice, - double frequency) { + double frequency, bool snap) { return (cdouble *)do_get_array_slice(where, components, fun, 0, fun_data, (void *)slice, - frequency); + frequency, snap); } -double *fields::get_array_slice(const volume &where, component c, double *slice, double frequency) { +double *fields::get_array_slice(const volume &where, component c, double *slice, double frequency, + bool snap) { std::vector components(1); components[0] = c; return (double *)do_get_array_slice(where, components, 0, default_field_rfunc, 0, (void *)slice, - frequency); + frequency, snap); } double *fields::get_array_slice(const volume &where, derived_component c, double *slice, - double frequency) { + double frequency, bool snap) { int nfields; component carray[12]; field_rfunction rfun = derived_component_func(c, gv, nfields, carray); std::vector cs(carray, carray + nfields); - return (double *)do_get_array_slice(where, cs, 0, rfun, &nfields, (void *)slice, frequency); + return (double *)do_get_array_slice(where, cs, 0, rfun, &nfields, (void *)slice, + frequency, snap); } cdouble *fields::get_complex_array_slice(const volume &where, component c, cdouble *slice, - double frequency) { + double frequency, bool snap) { std::vector components(1); components[0] = c; return (cdouble *)do_get_array_slice(where, components, default_field_func, 0, 0, (void *)slice, - frequency); + frequency, snap); } cdouble *fields::get_source_slice(const volume &where, component source_slice_component, diff --git a/src/meep.hpp b/src/meep.hpp index f51bad148..cf9e54349 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1658,24 +1658,27 @@ class fields { // must eventually be caller-deallocated via delete[]. double *get_array_slice(const volume &where, std::vector components, field_rfunction rfun, void *fun_data, double *slice = 0, - double frequency = 0); + double frequency = 0, + bool snap = false); std::complex *get_complex_array_slice(const volume &where, std::vector components, field_function fun, void *fun_data, std::complex *slice = 0, - double frequency = 0); + double frequency = 0, + bool snap = false); // alternative entry points for when you have no field // function, i.e. you want just a single component or // derived component.) double *get_array_slice(const volume &where, component c, double *slice = 0, - double frequency = 0); + double frequency = 0, bool snap = false); double *get_array_slice(const volume &where, derived_component c, double *slice = 0, - double frequency = 0); + double frequency = 0, bool snap = false); std::complex *get_complex_array_slice(const volume &where, component c, std::complex *slice = 0, - double frequency = 0); + double frequency = 0, + bool snap = false); // like get_array_slice, but for *sources* instead of fields std::complex *get_source_slice(const volume &where, component source_slice_component, @@ -1684,7 +1687,7 @@ class fields { // master routine for all above entry points void *do_get_array_slice(const volume &where, std::vector components, field_function fun, field_rfunction rfun, void *fun_data, void *vslice, - double frequency = 0); + double frequency = 0, bool snap = false); /* fetch and return coordinates and integration weights of grid points covered by an array slice, */ /* packed into a vector with format [NX, xtics[:], NY, ytics[:], NZ, ztics[:], weights[:] ] */ diff --git a/tests/Makefile.am b/tests/Makefile.am index 6a9964d60..53eae02ae 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -112,7 +112,7 @@ TESTS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrica if WITH_MPI LOG_COMPILER = $(RUNCODE) else - TESTS += cyl-ellipsoid-ll + TESTS += cyl-ellipsoid-ll array-slice-ll endif # Note: this requires GNU make diff --git a/tests/array-slice-ll-ref.h5 b/tests/array-slice-ll-ref.h5 index 9d359db83..c54e9c945 100644 Binary files a/tests/array-slice-ll-ref.h5 and b/tests/array-slice-ll-ref.h5 differ diff --git a/tests/array-slice-ll.cpp b/tests/array-slice-ll.cpp index 782cb2b1c..ba9ce857b 100644 --- a/tests/array-slice-ll.cpp +++ b/tests/array-slice-ll.cpp @@ -224,13 +224,13 @@ int main(int argc, char *argv[]) { // rank = f.get_array_slice_dimensions(v1d, dims1D, dirs1D, true); if (rank != 1 || dims1D[0] != NX) abort("incorrect dimensions for 1D slice"); - cdouble *slice1d = f.get_complex_array_slice(v1d, Hz); + cdouble *slice1d = f.get_complex_array_slice(v1d, Hz, 0, 0, true); double RelErr1D = Compare(slice1d, file_slice1d, NX, "Hz_1d"); master_printf("1D: rel error %e\n", RelErr1D); rank = f.get_array_slice_dimensions(v2d, dims2D, dirs2D, true); if (rank != 2 || dims2D[0] != NX || dims2D[1] != NY) abort("incorrect dimensions for 2D slice"); - double *slice2d = f.get_array_slice(v2d, Sy); + double *slice2d = f.get_array_slice(v2d, Sy, 0, 0, true); double RelErr2D = Compare(slice2d, file_slice2d, NX * NY, "Sy_2d"); master_printf("2D: rel error %e\n", RelErr2D);