plot geometry for dispersive materials without initializing structure object

doc/docs/Python_User_Interface.md

@@ -676,16 +676,21 @@ frequency-independent part of $\mu$ (the $\omega\to\infty$ limit).
 <div class="class_members" markdown="1">
 
 ```python
-def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None):
+def get_epsilon_grid(self,
+                     xtics=None,
+                     ytics=None,
+                     ztics=None,
+                     frequency=0):
 ```
 
 <div class="method_docstring" markdown="1">
 
 Given three 1d NumPy arrays (`xtics`,`ytics`,`ztics`) which define the coordinates of a Cartesian
-grid anywhere within the cell volume, compute the trace of the $\varepsilon$ tensor from the `geometry`
-exactly at each grid point. (For [`MaterialGrid`](#materialgrid)s, the $\varepsilon$ at each grid
-point is computed using bilinear interpolation from the nearest `MaterialGrid` points and possibly also
-projected to form a level set.) Note that this is different from `get_epsilon_point` which computes
+grid anywhere within the cell volume, compute the trace of the $\varepsilon(f)$ tensor at frequency
+$f$ (in Meep units) from the `geometry` exactly at each grid point. `frequency` defaults to 0 which is
+the instantaneous $\varepsilon$. (For [`MaterialGrid`](#materialgrid)s, the $\varepsilon$ at each
+grid point is computed using bilinear interpolation from the nearest `MaterialGrid` points and possibly
+also projected to form a level set.) Note that this is different from `get_epsilon_point` which computes
 $\varepsilon$ by bilinearly interpolating from the nearest Yee grid points. This function is useful for
 sampling the material geometry to any arbitrary resolution. The return value is a NumPy array with shape
 equivalent to `numpy.meshgrid(xtics,ytics,ztics)`. Empty dimensions are collapsed.
@@ -2067,7 +2072,7 @@ including outside the cell and a `near2far` object, returns the computed
 of fields $(E_x^1,E_y^1,E_z^1,H_x^1,H_y^1,H_z^1,E_x^2,E_y^2,E_z^2,H_x^2,H_y^2,H_z^2,...)$
 in Cartesian coordinates and
 $(E_r^1,E_\phi^1,E_z^1,H_r^1,H_\phi^1,H_z^1,E_r^2,E_\phi^2,E_z^2,H_r^2,H_\phi^2,H_z^2,...)$
-in cylindrical coordinates for the frequencies 1,2,…,`nfreq`.
+in cylindrical coordinates for the frequencies 1,2,...,`nfreq`.
 
 </div>
 
@@ -2609,7 +2614,7 @@ fr = mp.FluxRegion(volume=mp.GDSII_vol(fname, layer, zmin, zmax))
 
 ### Data Visualization
 
-This module provides basic visualization functionality for the simulation domain. The spirit of the module is to provide functions that can be called with *no customization options whatsoever* and will do useful relevant things by default, but which can also be customized in cases where you *do* want to take the time to spruce up the output. The `Simulation` class provides the following methods:
+This module provides basic visualization functionality for the simulation domain. The intent of the module is to provide functions that can be called with *no customization options whatsoever* and will do useful relevant things by default, but which can also be customized in cases where you *do* want to take the time to spruce up the output. The `Simulation` class provides the following methods:
 
 
 <a id="Simulation.plot2D"></a>
@@ -2648,14 +2653,18 @@ sim.run(...)
 field_func = lambda x: 20*np.log10(np.abs(x))
 import matplotlib.pyplot as plt
 sim.plot2D(fields=mp.Ez,
-        field_parameters={'alpha':0.8, 'cmap':'RdBu', 'interpolation':'none', 'post_process':field_func},
-        boundary_parameters={'hatch':'o', 'linewidth':1.5, 'facecolor':'y', 'edgecolor':'b', 'alpha':0.3})
+           field_parameters={'alpha':0.8, 'cmap':'RdBu', 'interpolation':'none', 'post_process':field_func},
+           boundary_parameters={'hatch':'o', 'linewidth':1.5, 'facecolor':'y', 'edgecolor':'b', 'alpha':0.3})
 plt.show()
 plt.savefig('sim_domain.png')
 ```
+If you just want to quickly visualize the simulation domain without the fields (i.e., when
+setting up your simulation), there is no need to invoke the `run` function prior to calling
+`plot2D`. Just define the `Simulation` object followed by any DFT monitors and then
+invoke `plot2D`.
 
-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.
+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:**
 
@@ -2677,6 +2686,12 @@ on all processes, but only generates a plot on the master process.
     - `alpha=1.0`: transparency of geometry
     - `contour=False`: if `True`, plot a contour of the geometry rather than its image
     - `contour_linewidth=1`: line width of the contour lines if `contour=True`
+    - `frequency=None`: for materials with a [frequency-dependent
+      permittivity](Materials.md#material-dispersion) $\varepsilon(f)$, specifies the
+      frequency $f$ (in Meep units) of the real part of the permittivity to use in the
+      plot. Defaults to the `frequency` parameter of the [Source](#source) object.
+    - `resolution=None`: the resolution of the $\varepsilon$ grid. Defaults to the
+      `resolution` of the `Simulation` object.
 * `boundary_parameters`: a `dict` of optional plotting parameters that override
   the default parameters for the boundary layers.
     - `alpha=1.0`: transparency of boundary layers
@@ -2713,10 +2728,6 @@ on all processes, but only generates a plot on the master process.
     - `alpha=0.6`: transparency of fields
     - `post_process=np.real`: post processing function to apply to fields (must be
       a function object)
-* `frequency`: for materials with a [frequency-dependent
-  permittivity](Materials.md#material-dispersion) $\varepsilon(f)$, specifies the
-  frequency $f$ (in Meep units) of the real part of the permittivity to use in the
-  plot. Defaults to the `frequency` parameter of the [Source](#source) object.
 
 </div>
 
@@ -4401,7 +4412,7 @@ def __init__(self,
              medium2,
              weights=None,
              grid_type='U_DEFAULT',
-             do_averaging=False,
+             do_averaging=True,
              beta=0,
              eta=0.5,
              damping=0):
@@ -6893,7 +6904,7 @@ A class used to record the fields during timestepping (i.e., a [`run`](#run-func
 function). The object is initialized prior to timestepping by specifying the
 simulation object and the field component. The object can then be passed to any
 [step-function modifier](#step-function-modifiers). For example, one can record the
-E<sub>z</sub> fields at every one time unit using:
+$E_z$ fields at every one time unit using:
 
 ```py
 animate = mp.Animate2D(sim,

doc/docs/Python_User_Interface.md.in

@@ -480,7 +480,7 @@ This feature is only available if Meep is built with [libGDSII](Build_From_Sourc
 
 ### Data Visualization
 
-This module provides basic visualization functionality for the simulation domain. The spirit of the module is to provide functions that can be called with *no customization options whatsoever* and will do useful relevant things by default, but which can also be customized in cases where you *do* want to take the time to spruce up the output. The `Simulation` class provides the following methods:
+This module provides basic visualization functionality for the simulation domain. The intent of the module is to provide functions that can be called with *no customization options whatsoever* and will do useful relevant things by default, but which can also be customized in cases where you *do* want to take the time to spruce up the output. The `Simulation` class provides the following methods:
 
 @@ Simulation.plot2D @@
 @@ Simulation.plot3D @@

python/meep.i

@@ -1097,12 +1097,12 @@ void _get_gradient(PyObject *grad, double scalegrad, PyObject *fields_a, PyObjec
     $1 = (double *)array_data($input);
 }
 
-%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* grid_vals {
+%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") std::complex<double>* grid_vals {
     $1 = is_array($input);
 }
 
-%typemap(in, fragment="NumPy_Macros") double* grid_vals {
-    $1 = (double *)array_data($input);
+%typemap(in, fragment="NumPy_Macros") std::complex<double>* grid_vals {
+    $1 = (std::complex<double> *)array_data($input);
 }
 
 // typemap for solve_cw:
@@ -2041,7 +2041,8 @@ void _get_epsilon_grid(geometric_object_list gobj_list,
                        int nx, double *xtics,
                        int ny, double *ytics,
                        int nz, double *ztics,
-                       double *grid_vals) {
+                       std::complex<double> *grid_vals,
+                       double frequency) {
      meep_geom::get_epsilon_grid(gobj_list,
                                  mlist,
                                  _default_material,
@@ -2052,7 +2053,8 @@ void _get_epsilon_grid(geometric_object_list gobj_list,
                                  nx, xtics,
                                  ny, ytics,
                                  nz, ztics,
-                                 grid_vals);
+                                 grid_vals,
+                                 frequency);
 }
 
 %}

python/simulation.py

@@ -2225,18 +2225,19 @@ def get_mu_point(self, pt, frequency=0):
         v3 = py_v3_to_vec(self.dimensions, pt, self.is_cylindrical)
         return self.fields.get_mu(v3,frequency)
 
-    def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None):
+    def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None, frequency=0):
         """
         Given three 1d NumPy arrays (`xtics`,`ytics`,`ztics`) which define the coordinates of a Cartesian
-        grid anywhere within the cell volume, compute the trace of the $\\varepsilon$ tensor from the `geometry`
-        exactly at each grid point. (For [`MaterialGrid`](#materialgrid)s, the $\\varepsilon$ at each grid
-        point is computed using bilinear interpolation from the nearest `MaterialGrid` points and possibly also
-        projected to form a level set.) Note that this is different from `get_epsilon_point` which computes
+        grid anywhere within the cell volume, compute the trace of the $\\varepsilon(f)$ tensor at frequency
+        $f$ (in Meep units) from the `geometry` exactly at each grid point. `frequency` defaults to 0 which is
+        the instantaneous $\\varepsilon$. (For [`MaterialGrid`](#materialgrid)s, the $\\varepsilon$ at each
+        grid point is computed using bilinear interpolation from the nearest `MaterialGrid` points and possibly
+        also projected to form a level set.) Note that this is different from `get_epsilon_point` which computes
         $\\varepsilon$ by bilinearly interpolating from the nearest Yee grid points. This function is useful for
         sampling the material geometry to any arbitrary resolution. The return value is a NumPy array with shape
         equivalent to `numpy.meshgrid(xtics,ytics,ztics)`. Empty dimensions are collapsed.
         """
-        grid_vals = np.squeeze(np.empty((len(xtics), len(ytics), len(ztics))))
+        grid_vals = np.squeeze(np.empty((len(xtics), len(ytics), len(ztics)), dtype=np.complex128))
         gv = self._create_grid_volume(False)
         mp._get_epsilon_grid(self.geometry,
                              self.extra_materials,
@@ -2247,7 +2248,8 @@ def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None):
                              len(xtics), xtics,
                              len(ytics), ytics,
                              len(ztics), ztics,
-                             grid_vals)
+                             grid_vals,
+                             frequency)
         return grid_vals
 
     def get_filename_prefix(self):
@@ -2755,8 +2757,8 @@ def get_farfield(self, near2far, x):
         (Fourier-transformed) "far" fields at `x` as list of length 6`nfreq`, consisting
         of fields $(E_x^1,E_y^1,E_z^1,H_x^1,H_y^1,H_z^1,E_x^2,E_y^2,E_z^2,H_x^2,H_y^2,H_z^2,...)$
         in Cartesian coordinates and
-        $(E_r^1,E_\phi^1,E_z^1,H_r^1,H_\phi^1,H_z^1,E_r^2,E_\phi^2,E_z^2,H_r^2,H_\phi^2,H_z^2,...)$
-        in cylindrical coordinates for the frequencies 1,2,…,`nfreq`.
+        $(E_r^1,E_\\phi^1,E_z^1,H_r^1,H_\\phi^1,H_z^1,E_r^2,E_\\phi^2,E_z^2,H_r^2,H_\\phi^2,H_z^2,...)$
+        in cylindrical coordinates for the frequencies 1,2,...,`nfreq`.
         """
         return mp._get_farfield(near2far.swigobj, py_v3_to_vec(self.dimensions, x, is_cylindrical=self.is_cylindrical))
 
@@ -4084,14 +4086,18 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
         field_func = lambda x: 20*np.log10(np.abs(x))
         import matplotlib.pyplot as plt
         sim.plot2D(fields=mp.Ez,
-                field_parameters={'alpha':0.8, 'cmap':'RdBu', 'interpolation':'none', 'post_process':field_func},
-                boundary_parameters={'hatch':'o', 'linewidth':1.5, 'facecolor':'y', 'edgecolor':'b', 'alpha':0.3})
+                   field_parameters={'alpha':0.8, 'cmap':'RdBu', 'interpolation':'none', 'post_process':field_func},
+                   boundary_parameters={'hatch':'o', 'linewidth':1.5, 'facecolor':'y', 'edgecolor':'b', 'alpha':0.3})
         plt.show()
         plt.savefig('sim_domain.png')
         ```
+        If you just want to quickly visualize the simulation domain without the fields (i.e., when
+        setting up your simulation), there is no need to invoke the `run` function prior to calling
+        `plot2D`. Just define the `Simulation` object followed by any DFT monitors and then
+        invoke `plot2D`.
 
-        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.
+        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:**
 
@@ -4113,6 +4119,12 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
             - `alpha=1.0`: transparency of geometry
             - `contour=False`: if `True`, plot a contour of the geometry rather than its image
             - `contour_linewidth=1`: line width of the contour lines if `contour=True`
+            - `frequency=None`: for materials with a [frequency-dependent
+              permittivity](Materials.md#material-dispersion) $\\varepsilon(f)$, specifies the
+              frequency $f$ (in Meep units) of the real part of the permittivity to use in the
+              plot. Defaults to the `frequency` parameter of the [Source](#source) object.
+            - `resolution=None`: the resolution of the $\\varepsilon$ grid. Defaults to the
+              `resolution` of the `Simulation` object.
         * `boundary_parameters`: a `dict` of optional plotting parameters that override
           the default parameters for the boundary layers.
             - `alpha=1.0`: transparency of boundary layers
@@ -4149,23 +4161,26 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
             - `alpha=0.6`: transparency of fields
             - `post_process=np.real`: post processing function to apply to fields (must be
               a function object)
-        * `frequency`: for materials with a [frequency-dependent
-          permittivity](Materials.md#material-dispersion) $\\varepsilon(f)$, specifies the
-          frequency $f$ (in Meep units) of the real part of the permittivity to use in the
-          plot. Defaults to the `frequency` parameter of the [Source](#source) object.
-
         """
-        return vis.plot2D(self, ax=ax, output_plane=output_plane, fields=fields, labels=labels,
-            eps_parameters=eps_parameters, boundary_parameters=boundary_parameters,
-            source_parameters=source_parameters,monitor_parameters=monitor_parameters,
-            field_parameters=field_parameters, frequency=frequency,
-            plot_eps_flag=plot_eps_flag, plot_sources_flag=plot_sources_flag,
-            plot_monitors_flag=plot_monitors_flag, plot_boundaries_flag=plot_boundaries_flag,
-            **kwargs)
-
-
-    def plot_fields(self,**kwargs):
-        return vis.plot_fields(self,**kwargs)
+        return vis.plot2D(self,
+                          ax=ax,
+                          output_plane=output_plane,
+                          fields=fields,
+                          labels=labels,
+                          eps_parameters=eps_parameters,
+                          boundary_parameters=boundary_parameters,
+                          source_parameters=source_parameters,
+                          monitor_parameters=monitor_parameters,
+                          field_parameters=field_parameters,
+                          frequency=frequency,
+                          plot_eps_flag=plot_eps_flag,
+                          plot_sources_flag=plot_sources_flag,
+                          plot_monitors_flag=plot_monitors_flag,
+                          plot_boundaries_flag=plot_boundaries_flag,
+                          **kwargs)
+
+    def plot_fields(self, **kwargs):
+        return vis.plot_fields(self, **kwargs)
 
     def plot3D(self):
         """

python/tests/test_visualization.py

@@ -161,7 +161,7 @@ def test_animation_output(self):
         sim = setup_sim() # generate 2D simulation
 
         Animate = mp.Animate2D(sim,fields=mp.Ez, realtime=False, normalize=False) # Check without normalization
-        Animate_norm = mp.Animate2D(sim,mp.Ez,realtime=False,normalize=True) # Check with normalization
+        Animate_norm = mp.Animate2D(sim, mp.Ez, realtime=False, normalize=True) # Check with normalization
 
         # test both animation objects during same run
         sim.run(