Plotting and animation

.travis.yml

@@ -17,6 +17,7 @@ addons:
     - libopenmpi-dev
     - mpi-default-bin
     - openmpi-bin
+    - ffmpeg
 
 ##################################################
 # environment variables applied to all build cases
@@ -64,7 +65,7 @@ matrix:
           - *common_deps
           - libhdf5-serial-dev
     install:
-      - pip install numpy mpi4py scipy h5py
+      - pip install numpy mpi4py scipy h5py matplotlib
   - python: "2.7"
     env:
       - MPICONF="--with-mpi"
@@ -79,7 +80,7 @@ matrix:
           - *common_deps
           - libhdf5-serial-dev
     install:
-      - pip install numpy mpi4py scipy h5py
+      - pip install numpy mpi4py scipy h5py matplotlib
   - python: "2.7"
     env:
       - MPICONF="--with-mpi"
@@ -95,7 +96,7 @@ matrix:
           - *common_deps
           - libhdf5-openmpi-dev
     install:
-      - pip install numpy mpi4py scipy
+      - pip install numpy mpi4py scipy matplotlib
       - pip install --no-binary=h5py h5py
   - python: "3.7"
     env:
@@ -110,7 +111,7 @@ matrix:
           - *common_deps
           - libhdf5-serial-dev
     install:
-      - pip install numpy mpi4py scipy h5py coverage coveralls
+      - pip install numpy mpi4py scipy h5py coverage coveralls matplotlib
   - python: "3.7"
     env:
       - MPICONF="--with-mpi"
@@ -125,7 +126,7 @@ matrix:
           - *common_deps
           - libhdf5-serial-dev
     install:
-      - pip install numpy mpi4py scipy h5py
+      - pip install numpy mpi4py scipy h5py matplotlib
   - python: "3.7"
     env:
       - MPICONF="--with-mpi"
@@ -141,7 +142,7 @@ matrix:
           - *common_deps
           - libhdf5-openmpi-dev
     install:
-      - pip install numpy mpi4py scipy
+      - pip install numpy mpi4py scipy matplotlib
       - pip install --no-binary=h5py h5py
 
 

doc/docs/Python_User_Interface.md

@@ -1525,10 +1525,112 @@ fr = mp.FluxRegion(volume=mp.GDSII_vol(fname, layer, zmin, zmax))
 
 ### Data Visualization
 
+This module attempts to provide some simple visualization functions that should just work in most cases to provide some graphical intuition for what you're working with. 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.
+
+**`Simulation.plot2D(ax=None, output_plane=None fields=None, labels=False, eps_parameters=None,boundary_parameters=None, source_parameters=None,monitor_parameters=None, field_parameters=None)`**
+—
+Plots a 2D cross section of the simulation domain using `matplotlib`. The plot includes the simulation geometry, boundary layers, sources, and monitors. Can also superimpose fields on top of 2D slice. Requires [matplotlib](https://matplotlib.org).
+
+* `ax`: A `matplotlib` axis object. `plot2D()` will add plot objects, like lines, patches, and scatter plots, to this object. If no `ax` is supplied, then the routine will create a new figure and grab its axis.
+* `output_plane`: A `Volume` object that specifies the plane over which to plot. Must be a 2D volume and a subset of the entire simulation volume (i.e. it should not span outside the domain at all).
+* `fields`: The field component (`mp.Ex`, `mp.Ey`, `mp.Ez`, `mp.Hx`, `mp.Hy`, `mp.Hz`) to superimpose over the simulation geometry. Default is `None`, where no fields are superimposed.
+* `labels`: If `True`, then labels will appear over each of the simulation elements.
+* `eps_parameters`: A `dict` of optional plotting parameters that override the default plotting parameters for the geometry.
+   * `interpolation='spline36'`: interpolation algorithm used to upsample the pixels.
+   * `cmap='binary'`: the color map of the geometry
+   * `alpha=1.0`: transparency of geometry  
+* `boundary_parameters`:A `dict` of optional plotting parameters that override the default plotting parameters for the boundary layers. Available parameters include:
+   * `alpha=1.0`: transparency of boundary layers
+   * `facecolor='g'`: color of polygon face
+   * `edgecolor='g'`: color of outline stroke
+   * `linewidth=1`: line width of outline stroke
+   * `hatch='\'`: hatching pattern
+* `source_parameters`: A `dict` of optional plotting parameters that override the default plotting parameters for the sources.
+   * `color='r`: color of line and pt sources
+   * `alpha=1.0`: transparency of source
+   * `facecolor='none'`: color of polygon face for planar sources
+   * `edgecolor='r'`: color of outline stroke for planar sources
+   * `linewidth=1`: line width of outline stroke
+   * `hatch='\'`: hatching pattern
+   * `label_color='r'`: color of source labels
+   * `label_alpha=0.3`: transparency of source label box
+   * `offset=20`: distance from source center and label box
+* `monitor_parameters`: A `dict` of optional plotting parameters that override the default plotting parameters for the monitors.
+   * `color='g`: color of line and point monitors
+   * `alpha=1.0`: transparency of monitors
+   * `facecolor='none'`: color of polygon face for planar monitors
+   * `edgecolor='r'`: color of outline stroke for planar monitors
+   * `linewidth=1`: line width of outline stroke
+   * `hatch='\'`: hatching pattern
+   * `label_color='g'`: color of source labels
+   * `label_alpha=0.3`: transparency of monitor label box
+   * `offset=20`: distance from monitor center and label box
+* `field_parameters`: A `dict` of optional plotting parameters that override the default plotting parameters for the fields.
+   * `interpolation='spline36'`: interpolation function used to upsample field pixels
+   * `cmap='RdBu'`: color map for field pixels
+   * `alpha=0.6`: transparency of fields
+
+**`Simulation.plot3D()`**
+— Uses Mayavi to render a 3D simulation domain. The simulation object must be 3D. Can also be embedded in jupyter notebooks.
+
 **`Simulation.visualize_chunks()`**
 —
 Displays an interactive image of how the cell is divided into chunks. Each rectangular region is a chunk, and each color represents a different processor. Requires [matplotlib](https://matplotlib.org).
 
+#### Animate2D
+A class used to record the fields while the simulation is timestepping (as a result of a `run` function). The user first intializes the object by specifying the simulation object and the field component of interest. The object can then be passed to any arbitrary step function modifier. For example, one can record the `E_z` component every 1 MEEP time units for 25 units with:
+```python
+animate = mp.Animate2D(sim,mp.Ez)
+sim.run(mp.at_every(1,animate),until=25)
+```
+
+By default, the object saves each frame as a png image into memory (not disk). This is typically more memory efficient than storing the actual fields. If the user sets the `normalize` argument, then the object will save the actual field information as a `numpy` array to be normalized for post processing. The user can also choose to update the fields of a figure in realtime by setting the `realtime` flag. This does not work for IPython or Jupyter notebooks, however.
+
+Once the simulation is run, the user can output the animation as an interactive JSHTML object, an mp4, or a GIF.
+
+Multiple Animate2D objects can be initialized and passed to the run function to track different volume locations (using `mp.in_volume`) or field components.
+
+Properties:
+
+**`sim`**
+— simulation object.
+
+**`fields`**
+— field component to record at each time step.
+
+**`f=None`**
+— optional `matplotlib` figure object that the routine will update on each call. If not supplied, then a new one will be created upon initialization.
+
+**`realtime=True`**
+— whether or not to update a figure window in realtime as the simulation progresses. Disabled by default. Not compatible with IPython/Jupyter notebooks.
+
+**`normalize=False`**
+— records fields at each time step in memory in a numpy array and then normalizes the result for future post-processing.
+
+**`plot_modifiers=None`**
+— A list of functions that can modify the figure's `axis` object. Each function modifier accepts a single argument, an `axis` object, and must return that same axis object. The following modifier changes the `xlabel`:
+```python
+def mod1(ax):
+    ax.set_xlabel('Testing')
+    return ax
+
+plot_modifiers = [mod1]
+```
+
+**`**customization_args`**
+— customization keyword arguments pass to `plot2D()` (i.e. `labels`, `eps_parameters`, `boundary_parameters`, etc.)
+
+Methods:
+
+**`Animate2D.to_jshtml(fps)`**
+— Outputs an interactable JSHTML animation object that is embeddable in Jupyter notebooks. The object is packaged with controls to manipulate the video's playback. User must specify a frame rate `fps` in frames per second.
+
+**`Animate2D.to_mp4(fps,filename)`**
+— Generates and outputs an mp4 video file of the animation with the filename, `filename`, and the frame rate, `fps`. Default encoding is h264 with yuv420p format. Requires `ffmpeg`.
+
+**`Animate2D.to_mp4(fps,filename)`**
+— Generates and outputs a GIF file of the animation with the filename, `filename` and the frame rate, `fps`. Requires `ffmpeg`.
+
 Run and Step Functions
 ----------------------
 

python/Makefile.am

@@ -68,6 +68,7 @@ TESTS =                                   \
     $(TEST_DIR)/simulation.py             \
     $(TEST_DIR)/source.py                 \
     $(TEST_DIR)/user_defined_material.py  \
+    $(TEST_DIR)/visualization.py          \
     $(WVG_SRC_TEST)
 
 if WITH_COVERAGE
@@ -187,6 +188,7 @@ HL_IFACE =                  \
     $(srcdir)/geom.py       \
     $(srcdir)/simulation.py \
     $(srcdir)/source.py     \
+    $(srcdir)/visualization.py  \
     $(srcdir)/materials.py
 
 pkgpython_PYTHON = __init__.py $(HL_IFACE)
@@ -224,3 +226,5 @@ clean-local:
 
 distclean-local:
 	rm -f *.h5
+	rm -f *.mp4
+	rm -f *.gif

python/examples/README.md

@@ -0,0 +1,42 @@
+# Tutorials and Examples
+
+Meep simulations are Python scripts which involve specifying the device geometry, materials, current sources, monitor fields, and everything else necessary to set up a calculation. A Python script provides the flexibility to customize the simulation for practically any application particularly those involving parameter sweeps and optimization.
+
+Python libraries such as NumPy, SciPy, and Matplotlib can be used to augment the simulation functionality and will also be demonstrated. Much of the functionality of the low-level C++ interface has been abstracted in Python which means that you don't need to be an experienced programmer to set up simulations. Reasonable defaults are available where necessary.
+
+Several tutorials and examples are found here. These tutorials are meant to illustrate Meep's various features in an interactive and application-oriented manner. 
+
+## iPython/Jupyter Notebooks
+
+Jupyter notebooks are interactive, browser based framework for displaying text, running python code, and visualizing results. There are several ways to read these notebooks online:
+
+### Local
+
+The recommended method to run the tutorial notebooks is by 1.) installing `meep` via `conda`, 2.) cloning this repo to your local drive, and 3.) launch a notebook server in this directory using `jupyter notebook`.
+
+### nbviewer
+
+`nbviewer` is a web platform that can render interactive features found inside Jupyter notebooks openly stored on web-servers. While `nbviewer` can't run python code, it can execute stored javascript code used to animate the simulations. 
+
+### GitHub
+
+GitHub is able to render some of the smaller notebooks as plain text. However, they are not interactive and are often too large for GitHub.
+
+### Tutorials
+
+Below are summaries for each tutorial, along with the features the tutorials highlight. While there is no particular order to the tutorials, they progressively incorporate more complicated features.
+
+1. __`straight-waveguide.ipynb`__ -
+A simple 2D straight waveguide tutorial that explores basic meep features like `geometry`, `sources`, and `PML` layers. The tutorial also explores basic visualization and animation features.
+
+2. __`bent-waveguide.ipynb`__-
+A followup to the 2D straight waveguide tutorial by adding a bend.
+
+3. __`bend-flux.ipynb`__-
+Using the previous bent waveguide example, this tutorial calculates the loss, transmission, and reflection that the bent waveguide undergoes.
+
+4. __`ring.ipynb`__ -
+Computes the resonant mode frequencies of a 2D ring resonator using `harminv`.
+
+5. __`visualization.ipynb`__ -
+Demonstrates various visualization and animation features.

python/meep.i

@@ -1514,6 +1514,12 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where
         SourceTime,
         check_positive,
     )
+    from .visualization import (
+        plot2D,
+        plot3D,
+        plot_fields,
+        Animate2D
+    )
 
     if with_mpi():
         try:

python/simulation.py

@@ -18,6 +18,7 @@
 import meep as mp
 from meep.geom import Vector3, init_do_averaging
 from meep.source import EigenModeSource, check_positive
+import meep.visualization as vis
 
 
 try:
@@ -158,6 +159,48 @@ def __init__(self, center, size=Vector3(), dims=2, is_cylindrical=False):
         vec2 = py_v3_to_vec(self.dims, v2, is_cylindrical)
 
         self.swigobj = mp.volume(vec1, vec2)
+
+    def get_vertices(self):
+        xmin = self.center.x - self.size.x/2
+        xmax = self.center.x + self.size.x/2
+        ymin = self.center.y - self.size.y/2
+        ymax = self.center.y + self.size.y/2
+        zmin = self.center.z - self.size.z/2
+        zmax = self.center.z + self.size.z/2
+
+        # Iterate over and remove duplicates for collapsed dimensions (i.e. min=max))
+        return [Vector3(x,y,z) for x in list(set([xmin,xmax])) for y in list(set([ymin,ymax])) for z in list(set([zmin,zmax]))]
+
+
+    def get_edges(self):
+        vertices = self.get_vertices()
+        edges = []
+
+        # Useful for importing weird geometries and the sizes are slightly off
+        def nearly_equal(a,b,sig_fig=10):
+            return a==b or (abs(a-b) < 10**(-sig_fig))
+
+        for iter1 in range(len(vertices)):
+            for iter2 in range(iter1+1,len(vertices)):
+                if ((iter1 != iter2) and 
+                nearly_equal((vertices[iter1]-vertices[iter2]).norm(),self.size.x) or 
+                nearly_equal((vertices[iter1]-vertices[iter2]).norm(),self.size.y) or
+                nearly_equal((vertices[iter1]-vertices[iter2]).norm(),self.size.z)):
+                    edges.append([vertices[iter1],vertices[iter2]])
+        return edges
+
+    def pt_in_volume(self,pt):
+        xmin = self.center.x - self.size.x/2
+        xmax = self.center.x + self.size.x/2
+        ymin = self.center.y - self.size.y/2
+        ymax = self.center.y + self.size.y/2
+        zmin = self.center.z - self.size.z/2
+        zmax = self.center.z + self.size.z/2
+
+        if (pt.x >= xmin and pt.x <= xmax and pt.y >= ymin and pt.y <= ymax and pt.z >= zmin and pt.z <= zmax):
+            return True
+        else:
+            return False
 
 
 class FluxRegion(object):
@@ -2142,65 +2185,17 @@ def get_sfield_r(self):
     def get_sfield_p(self):
         return self.get_array(mp.Sp, cmplx=True)
 
-    def visualize_chunks(self):
-        if self.structure is None:
-            self.init_sim()
-        try:
-            import matplotlib.pyplot as plt
-            import matplotlib.cm
-            import matplotlib.colors
-            from mpl_toolkits.mplot3d import Axes3D
-            from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
-        except ImportError:
-            warnings.warn("matplotlib is required for visualize_chunks", ImportWarning)
-            return
-
-        vols = self.structure.get_chunk_volumes()
-        owners = self.structure.get_chunk_owners()
-
-        fig = plt.figure()
-        ax = fig.add_subplot(111, projection='3d')
-        chunk_colors = matplotlib.cm.rainbow(np.linspace(0, 1, mp.count_processors()))
-
-        def plot_box(box, proc):
-            low = mp.Vector3(box.low.x, box.low.y, box.low.z)
-            high = mp.Vector3(box.high.x, box.high.y, box.high.z)
-            points = [low, high]
-
-            x_len = mp.Vector3(high.x) - mp.Vector3(low.x)
-            y_len = mp.Vector3(y=high.y) - mp.Vector3(y=low.y)
-            xy_len = mp.Vector3(high.x, high.y) - mp.Vector3(low.x, low.y)
-
-            points += [low + x_len]
-            points += [low + y_len]
-            points += [low + xy_len]
-            points += [high - x_len]
-            points += [high - y_len]
-            points += [high - xy_len]
-            points = np.array([np.array(v) for v in points])
-
-            edges = [
-                [points[0], points[2], points[4], points[3]],
-                [points[1], points[5], points[7], points[6]],
-                [points[0], points[3], points[5], points[7]],
-                [points[1], points[4], points[2], points[6]],
-                [points[3], points[4], points[1], points[5]],
-                [points[0], points[7], points[6], points[2]]
-            ]
+    def plot2D(self,**kwargs):
+        return vis.plot2D(self,**kwargs)
+
+    def plot_fields(self,**kwargs):
+        return vis.plot_fields(self,**kwargs)
+
+    def plot3D(self):
+        return vis.plot3D(self)
 
-            faces = Poly3DCollection(edges, linewidths=1, edgecolors='k')
-            color_with_alpha = matplotlib.colors.to_rgba(chunk_colors[proc], alpha=0.2)
-            faces.set_facecolor(color_with_alpha)
-            ax.add_collection3d(faces)
-
-            # Plot the points themselves to force the scaling of the axes
-            ax.scatter(points[:, 0], points[:, 1], points[:, 2], s=0)
-
-        if mp.am_master():
-            for i, v in enumerate(vols):
-                plot_box(mp.gv2box(v.surroundings()), owners[i])
-            ax.set_aspect('auto')
-            plt.show()
+    def visualize_chunks(self):
+        vis.visualize_chunks(self)
 
 
 def _create_boundary_region_from_boundary_layers(boundary_layers, gv):