Draw perfect metals (ε = -∞) in dark red in plot2D

doc/docs/Python_User_Interface.md

@@ -3048,6 +3048,10 @@ to be called on all processes, but only generates a plot on the master process.
     - `resolution=None`: the resolution of the $\varepsilon$ grid. Defaults to the
       `resolution` of the `Simulation` object.
     - `colorbar=False`: whether to add a colorbar to the plot's parent Figure based on epsilon values.
+    - `pec_color='darkred'`: the color used to draw perfect metals / perfect
+      electric conductors, whose permittivity is $-\infty$. Drawing these
+      separately keeps them from saturating the color scale of the
+      finite-permittivity materials.
 * `boundary_parameters`: a `dict` of optional plotting parameters that override
   the default parameters for the boundary layers.
     - `alpha=1.0`: transparency of boundary layers

python/simulation.py

@@ -5052,6 +5052,10 @@ def plot2D(
             - `resolution=None`: the resolution of the $\\varepsilon$ grid. Defaults to the
               `resolution` of the `Simulation` object.
             - `colorbar=False`: whether to add a colorbar to the plot's parent Figure based on epsilon values.
+            - `pec_color='darkred'`: the color used to draw perfect metals / perfect
+              electric conductors, whose permittivity is $-\\infty$. Drawing these
+              separately keeps them from saturating the color scale of the
+              finite-permittivity materials.
         * `boundary_parameters`: a `dict` of optional plotting parameters that override
           the default parameters for the boundary layers.
             - `alpha=1.0`: transparency of boundary layers

python/tests/test_visualization.py

@@ -3,6 +3,7 @@
 # boundary conditions. Checks for subdomain plots.
 #
 # Also tests the animation run function, mp4 output, jshtml output, and git output.
+import importlib.util
 import os
 import unittest
 from subprocess import call
@@ -18,6 +19,7 @@
 import io
 
 from matplotlib import pyplot as plt
+from matplotlib.contour import ContourSet
 
 
 def hash_figure(fig):
@@ -327,6 +329,94 @@ def test_plot2D(self):
             hash_figure(f)
             # self.assertAlmostEqual(hash_figure(f),68926258)
 
+    @staticmethod
+    def setup_pec_sim():
+        return mp.Simulation(
+            cell_size=mp.Vector3(6, 6),
+            resolution=20,
+            geometry=[
+                mp.Block(
+                    size=mp.Vector3(2, 2),
+                    center=mp.Vector3(-1.2),
+                    material=mp.Medium(epsilon=12),
+                ),
+                mp.Block(
+                    size=mp.Vector3(2, 2),
+                    center=mp.Vector3(1.2),
+                    material=mp.metal,
+                ),
+            ],
+            sources=[
+                mp.Source(
+                    src=mp.ContinuousSource(frequency=0.15),
+                    component=mp.Ez,
+                    center=mp.Vector3(0, -2),
+                )
+            ],
+        )
+
+    def test_plot2D_pec(self):
+        # Metals have eps = -inf, which must not saturate the color scale of
+        # the finite-permittivity materials plotted alongside them.
+        sim = self.setup_pec_sim()
+
+        f = plt.figure()
+        ax = sim.plot2D(ax=f.gca(), eps_parameters={"colorbar": True})
+        if mp.am_master():
+            im = ax.get_images()[0]
+
+            # The metal pixels are masked out of the color scale...
+            self.assertTrue(np.all(np.isfinite(im.get_clim())))
+            self.assertAlmostEqual(im.get_clim()[0], 1)
+            self.assertAlmostEqual(im.get_clim()[1], 12)
+
+            # ...and drawn in the (opaque) `pec_color` instead.
+            self.assertTrue(
+                np.allclose(
+                    im.get_cmap().get_bad(), matplotlib.colors.to_rgba("darkred")
+                )
+            )
+            self.assertTrue(np.any(im.get_array().mask))
+
+            # The default pec_color must not coincide with any color of the
+            # default colormap, which is grayscale and therefore spans every
+            # gray; otherwise metals look like some finite permittivity.
+            ramp = im.get_cmap()(np.linspace(0, 1, 256))[:, :3]
+            pec_rgb = np.array(matplotlib.colors.to_rgb("darkred"))
+            self.assertGreater(np.min(np.linalg.norm(ramp - pec_rgb, axis=1)), 0.25)
+
+            # The module-level defaults must not have been mutated.
+            self.assertEqual(mp.visualization.default_eps_parameters["cmap"], "binary")
+
+        # A user-specified color is honored.
+        f = plt.figure()
+        ax = sim.plot2D(ax=f.gca(), eps_parameters={"pec_color": "red"})
+        if mp.am_master():
+            self.assertTrue(
+                np.allclose(
+                    ax.get_images()[0].get_cmap().get_bad(),
+                    matplotlib.colors.to_rgba("red"),
+                )
+            )
+
+    @unittest.skipIf(
+        importlib.util.find_spec("contourpy") is None, "contourpy is not installed"
+    )
+    def test_plot2D_pec_contour(self):
+        # The contour path must not emit the metal's -inf as a contour level:
+        # the geometry only has the two finite permittivities 1 and 12, so the
+        # smallest level is 1. Without masking the metal it would be -mp.inf.
+        sim = self.setup_pec_sim()
+        f = plt.figure()
+        ax = sim.plot2D(ax=f.gca(), eps_parameters={"contour": True})
+        if mp.am_master():
+            contours = [c for c in ax.collections if isinstance(c, ContourSet)]
+            self.assertEqual(len(contours), 1)
+            levels = contours[0].levels
+            self.assertTrue(np.all(np.isfinite(levels)))
+            self.assertGreaterEqual(np.min(levels), 1)
+            self.assertLessEqual(np.max(levels), 12)
+
     @unittest.skipIf(call(["which", "ffmpeg"]) != 0, "ffmpeg is not installed")
     def test_animation_output(self):
         # ------------------------- #

python/visualization.py

@@ -14,6 +14,7 @@
 
 ## Typing imports
 from matplotlib.axes import Axes
+from matplotlib.colors import Colormap
 from matplotlib.figure import Figure
 from typing import Callable, Union, Any, Tuple, List, Optional
 
@@ -61,6 +62,7 @@
     "frequency": None,
     "resolution": None,
     "colorbar": False,
+    "pec_color": "darkred",
 }
 
 default_colorbar_parameters = {
@@ -556,9 +558,21 @@ def sort_points(xy):
     return ax
 
 
+def _get_colormap(cmap: Union[str, Colormap]) -> Colormap:
+    """Return a Colormap object given a colormap name or Colormap object."""
+    import matplotlib as mpl
+
+    if isinstance(cmap, Colormap):
+        return cmap
+    try:
+        return mpl.colormaps[cmap]  # matplotlib >= 3.5
+    except AttributeError:
+        return mpl.cm.get_cmap(cmap)  # matplotlib < 3.5
+
+
 def _add_colorbar(
     ax: Axes,
-    cmap: str,
+    cmap: Union[str, Colormap],
     vmin: float,
     vmax: float,
     default_label: Optional[str] = None,
@@ -569,7 +583,7 @@ def _add_colorbar(
     from mpl_toolkits.axes_grid1 import make_axes_locatable
 
     if colorbar_parameters is None:
-        colorbar_parameters = copy.deepcopy(default_colorbar_parameters)
+        colorbar_parameters = copy.copy(default_colorbar_parameters)
     else:
         colorbar_parameters = dict(default_colorbar_parameters, **colorbar_parameters)
 
@@ -578,10 +592,9 @@ def _add_colorbar(
         colorbar_parameters["label"] = default_label
 
     # Create a map between field/eps values and colors in the colormap.
-    # Note: cm.get_cmap() is deprecated for matplotlib>=3.6, use mpl.colormaps[cmap] instead if necessary.
     sm = mpl.cm.ScalarMappable(
         norm=mpl.colors.Normalize(vmin, vmax),
-        cmap=mpl.cm.get_cmap(cmap),
+        cmap=_get_colormap(cmap),
     )
 
     # Pop specific values out of colorbar params so user can add any kwargs to plt.colorbar
@@ -607,7 +620,7 @@ def plot_eps(
 ) -> Union[Axes, Any]:
     # consolidate plotting parameters
     if eps_parameters is None:
-        eps_parameters = default_eps_parameters
+        eps_parameters = copy.copy(default_eps_parameters)
     else:
         eps_parameters = dict(default_eps_parameters, **eps_parameters)
 
@@ -676,11 +689,25 @@ def plot_eps(
         if not ax:
             return eps_data
 
+        # Metals / perfect electric conductors have eps = -inf (see python/meep.i,
+        # where mp.inf is the sentinel 1e20), which would otherwise saturate the
+        # color scale. Mask them out so that the scale is set by the finite
+        # permittivities, and draw them in `pec_color` instead. Note that the
+        # default `pec_color` is deliberately not a gray: the default `cmap`
+        # ('binary') spans every gray, so a gray would be indistinguishable from
+        # some finite permittivity.
+        eps_data = np.ma.masked_where(
+            ~np.isfinite(eps_data) | (eps_data <= -mp.inf), eps_data
+        )
+        cmap = copy.copy(_get_colormap(eps_parameters["cmap"]))
+        cmap.set_bad(color=eps_parameters["pec_color"])
+        eps_parameters["cmap"] = cmap
+
         if eps_parameters["contour"]:
             ax.contour(
                 eps_data,
                 0,
-                levels=np.unique(eps_data),
+                levels=np.unique(eps_data.compressed()),
                 colors="black",
                 origin="upper",
                 extent=extent,
@@ -690,11 +717,12 @@ def plot_eps(
             ax.imshow(eps_data, extent=extent, **filter_dict(eps_parameters, ax.imshow))
 
         if eps_parameters["colorbar"]:
+            finite_eps = eps_data.compressed()
             _add_colorbar(
                 ax=ax,
                 cmap=eps_parameters["cmap"],
-                vmin=np.amin(eps_data),
-                vmax=np.amax(eps_data),
+                vmin=np.amin(finite_eps) if finite_eps.size else 0.0,
+                vmax=np.amax(finite_eps) if finite_eps.size else 1.0,
                 default_label=r"$\epsilon_r$",
                 colorbar_parameters=colorbar_parameters,
             )
@@ -1072,7 +1100,7 @@ def plot_1d_index(
 ) -> Union[Axes, Any]:
     """Plots the refractive-index profile n(z) of a 1D simulation."""
     if index_parameters is None:
-        index_parameters = default_1d_index_parameters
+        index_parameters = copy.copy(default_1d_index_parameters)
     else:
         index_parameters = dict(default_1d_index_parameters, **index_parameters)