unit test and documentation for MaterialGrid (#1508)

* unit test and documentation for MaterialGrid

* fixes and tweaks

* more fixes

* fix test

doc/docs/Python_User_Interface.md

@@ -4156,6 +4156,89 @@ media](FAQ.md#why-does-my-simulation-diverge-if-0).
 
 </div>
 
+---
+<a id="MaterialGrid"></a>
+
+### MaterialGrid
+
+```python
+class MaterialGrid(object):
+```
+
+<div class="class_docstring" markdown="1">
+
+This class is used to specify materials interpolated from a discrete Cartesian grid. A class object is passed as the
+`material` argument of a [`GeometricObject`](#GeometricObject) or the `default_material` argument of the `Simulation`
+constructor (similar to a material function).
+
+</div>
+
+
+
+<a id="MaterialGrid.__init__"></a>
+
+<div class="class_members" markdown="1">
+
+```python
+def __init__(self,
+             grid_size,
+             medium1,
+             medium2,
+             design_parameters=None,
+             grid_type='U_DEFAULT',
+             do_averaging=False,
+             beta=0,
+             eta=0.5):
+```
+
+<div class="method_docstring" markdown="1">
+
+Creates a `MaterialGrid` object.
+
+The input are two materials `medium1` and `medium2` which are linearly interpolated at each grid point using
+a NumPy array `design_parameters` of size `grid_size` (a 3-tuple or `Vector3` of integers) with floating-point values in
+the range [0,1] used to define a linear weight (i.e., 0 is `medium1` and 1 is `medium2`). Currently, only two material
+types are supported: (1) frequency-independent isotropic $\varepsilon$ or $\mu$ and (2) `LorentzianSusceptibility`.
+`medium1` and `medium2` must both be the same type.
+
+[Subpixel smoothing](Subpixel_Smoothing.md) can be enabled by specifying `do_averaging=True`. If you want to use a
+material grid to define a (nearly) discontinuous, piecewise-constant material that is either `medium1` or `medium2`
+almost everywhere, you can optionally enable a (smoothed) *projection* feature by setting the parameter `beta` to a
+positive value. When the projection feature is enabled, the design parameters $u(x)$ can be thought of as a [level-set
+function](https://en.wikipedia.org/wiki/Level-set_method) defining an interface at $u(x)=\eta$ with a smoothing factor
+$\beta$ ($\beta=\infty$ gives an unsmoothed, discontinuous interface). The projection operator is $\(tanh(\beta*\eta)
++ \tanh(\beta*(u-\eta))) / (\tanh(\beta*\eta) + \tanh(\beta*(1-\eta)))$ involving the parameters `beta`
+($\beta$: "smoothness" of the turn on) and `eta` ($\eta$: erosion/dilation). The level set provides a general approach for
+defining a *discontinuous* function of the otherwise continuously varying (via the bilinear interpolation) grid values.
+The subpixel smoothing is based on an adaptive quadrature scheme with properties `subpixel_maxeval` and `subpixel_tol` which
+can be specified using the [`Simulation`](#Simulation) constructor.
+
+Grids which are symmetric (e.g., mirror, rotation) must be explicitly defined. One way to implement this is by overlapping
+a given `MaterialGrid` object with a symmetrized copy of itself. In the case of spatially overlapping `MaterialGrid`
+objects (with no intervening objects), any overlapping points are combined using the method `grid_type` which is one of
+`"U_MIN"` (minimum of the overlapping grid values), `"U_PROD"` (product), `"U_SUM"` (mean), `"U_DEFAULT"`
+(topmost material at grid point).
+
+</div>
+
+</div>
+
+<a id="MaterialGrid.update_parameters"></a>
+
+<div class="class_members" markdown="1">
+
+```python
+def update_parameters(self, x):
+```
+
+<div class="method_docstring" markdown="1">
+
+Reset the design grid `design_parameters` to `x`.
+
+</div>
+
+</div>
+
 ---
 <a id="Susceptibility"></a>
 

doc/docs/Python_User_Interface.md.in

@@ -714,6 +714,8 @@ The following classes are available directly via the `meep` package.
 
 @@ Medium[methods-with-docstrings] @@
 
+@@ MaterialGrid[methods-with-docstrings] @@
+
 @@ Susceptibility[methods-with-docstrings] @@
 
 @@ LorentzianSusceptibility[methods-with-docstrings] @@

python/geom.py

@@ -523,7 +523,39 @@ def _get_epsmu(self, diag, offdiag, susceptibilities, conductivity_diag, conduct
         return np.squeeze(epsmu)
 
 class MaterialGrid(object):
+    """
+    This class is used to specify materials interpolated from a discrete Cartesian grid. A class object is passed as the
+    `material` argument of a [`GeometricObject`](#GeometricObject) or the `default_material` argument of the `Simulation`
+    constructor (similar to a material function).
+    """
     def __init__(self,grid_size,medium1,medium2,design_parameters=None,grid_type="U_DEFAULT",do_averaging=False,beta=0,eta=0.5):
+        """
+        Creates a `MaterialGrid` object.
+
+        The input are two materials `medium1` and `medium2` which are linearly interpolated at each grid point using
+        a NumPy array `design_parameters` of size `grid_size` (a 3-tuple or `Vector3` of integers) with floating-point values in
+        the range [0,1] used to define a linear weight (i.e., 0 is `medium1` and 1 is `medium2`). Currently, only two material
+        types are supported: (1) frequency-independent isotropic $\\varepsilon$ or $\\mu$ and (2) `LorentzianSusceptibility`.
+        `medium1` and `medium2` must both be the same type.
+
+        [Subpixel smoothing](Subpixel_Smoothing.md) can be enabled by specifying `do_averaging=True`. If you want to use a
+        material grid to define a (nearly) discontinuous, piecewise-constant material that is either `medium1` or `medium2`
+        almost everywhere, you can optionally enable a (smoothed) *projection* feature by setting the parameter `beta` to a
+        positive value. When the projection feature is enabled, the design parameters $u(x)$ can be thought of as a [level-set
+        function](https://en.wikipedia.org/wiki/Level-set_method) defining an interface at $u(x)=\\eta$ with a smoothing factor
+        $\\beta$ ($\\beta=\\infty$ gives an unsmoothed, discontinuous interface). The projection operator is $(\\tanh(\\beta*\\eta)
+        + \\tanh(\\beta*(u-\\eta))) / (\\tanh(\\beta*\\eta) + \\tanh(\\beta*(1-\\eta)))$ involving the parameters `beta`
+        ($\\beta$: "smoothness" of the turn on) and `eta` ($\\eta$: erosion/dilation). The level set provides a general approach for
+        defining a *discontinuous* function of the otherwise continuously varying (via the bilinear interpolation) grid values.
+        The subpixel smoothing is based on an adaptive quadrature scheme with properties `subpixel_maxeval` and `subpixel_tol` which
+        can be specified using the [`Simulation`](#Simulation) constructor.
+
+        Grids which are symmetric (e.g., mirror, rotation) must be explicitly defined. One way to implement this is by overlapping
+        a given `MaterialGrid` object with a symmetrized copy of itself. In the case of spatially overlapping `MaterialGrid`
+        objects (with no intervening objects), any overlapping points are combined using the method `grid_type` which is one of
+        `"U_MIN"` (minimum of the overlapping grid values), `"U_PROD"` (product), `"U_SUM"` (mean), `"U_DEFAULT"`
+        (topmost material at grid point).
+        """
         self.grid_size = mp.Vector3(*grid_size)
         self.medium1 = medium1
         self.medium2 = medium2
@@ -558,6 +590,9 @@ def isclose(a, b, rel_tol=1e-09, abs_tol=0.0):
 
         self.swigobj = None
     def update_parameters(self,x):
+        """
+        Reset the design grid `design_parameters` to `x`.
+        """
         if x.size != self.num_params:
             raise ValueError("design_parameters of shape {} do not match user specified grid dimension: {}".format(self.design_parameters.size,self.grid_size))
         self.design_parameters[:]=x.flatten().astype(np.float64)

python/tests/material_grid.py

@@ -1,63 +1,82 @@
 import meep as mp
 import numpy as np
-from meep.materials import Si, SiO2
+from scipy.ndimage import gaussian_filter
 import unittest
 
-class TestMaterialGrid(unittest.TestCase):
-    def gen_sim(self):
-        resolution = 10 # pixels/um
-        cell_size = mp.Vector3(14,14)
-        pml_layers = [mp.PML(thickness=2)]
-        w = 3.0 # width of waveguide
-        m1 = SiO2
-        m2 = Si
-        n = 10
-        gs = mp.Vector3(n,n)
-        np.random.seed(1)
-        dp = np.random.rand(n*n)
-        mg = mp.MaterialGrid(gs,m1,m2,dp,"U_SUM")
-        geometry = []
-        rot_angle = np.radians(0)
-        geometry += [mp.Block(center=mp.Vector3(),
-                            size=mp.Vector3(w,w,mp.inf),
-                            e1=mp.Vector3(x=1).rotate(mp.Vector3(z=1), rot_angle),
-                            e2=mp.Vector3(y=1).rotate(mp.Vector3(z=1), rot_angle),
-                            material=mg)]
-        geometry += [mp.Block(center=mp.Vector3(),
-                            size=mp.Vector3(w,w,mp.inf),
-                            e1=mp.Vector3(x=-1).rotate(mp.Vector3(z=1), np.pi/2),
-                            e2=mp.Vector3(y=1).rotate(mp.Vector3(z=1), np.pi/2),
-                            material=mg)]
-        fsrc = 1/1.55 # frequency of eigenmode or constant-amplitude source
-        bnum = 1    # band number of eigenmode
-        kpoint = mp.Vector3(x=1).rotate(mp.Vector3(z=1), rot_angle)
-        compute_flux = True # compute flux (True) or plot the field profile (False)
-        eig_src = True # eigenmode (True) or constant-amplitude (False) source
-        if eig_src:
-            sources = [mp.EigenModeSource(src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc) if compute_flux else mp.ContinuousSource(fsrc),
-                                        center=mp.Vector3(),
-                                        size=mp.Vector3(y=3*w),
-                                        direction=mp.NO_DIRECTION,
-                                        eig_kpoint=kpoint,
-                                        eig_band=bnum,
-                                        eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z,
-                                        eig_match_freq=True)]
-        else:
-            sources = [mp.Source(src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc) if compute_flux else mp.ContinuousSource(fsrc),
-                                center=mp.Vector3(),
-                                size=mp.Vector3(y=3*w),
-                                component=mp.Ez)]
-
-        sim = mp.Simulation(cell_size=cell_size,
-                            resolution=resolution,
-                            boundary_layers=pml_layers,
+def compute_resonant_mode(res):
+        cell_size = mp.Vector3(1,1,0)
+
+        rad = 0.301943
+
+        fcen = 0.3
+        df = 0.2*fcen
+        sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
+                             component=mp.Hz,
+                             center=mp.Vector3(-0.1057,0.2094,0))]
+
+        k_point = mp.Vector3(0.3892,0.1597,0)
+
+        design_shape = mp.Vector3(1,1,0)
+        design_region_resolution = 50
+        Nx = int(design_region_resolution*design_shape.x)
+        Ny = int(design_region_resolution*design_shape.y)
+        x = np.linspace(-0.5*design_shape.x,0.5*design_shape.x,Nx)
+        y = np.linspace(-0.5*design_shape.y,0.5*design_shape.y,Ny)
+        xv, yv = np.meshgrid(x,y)
+        design_params = np.sqrt(np.square(xv) + np.square(yv)) < rad
+        filtered_design_params = gaussian_filter(design_params,
+                                                 sigma=3.0,
+                                                 output=np.double)
+
+        matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
+                                  mp.air,
+                                  mp.Medium(index=3.5),
+                                  design_parameters=filtered_design_params,
+                                  do_averaging=True,
+                                  beta=1000,
+                                  eta=0.5)
+
+        geometry = [mp.Block(center=mp.Vector3(),
+                             size=mp.Vector3(design_shape.x,design_shape.y,0),
+                             material=matgrid)]
+
+        sim = mp.Simulation(resolution=res,
+                            cell_size=cell_size,
+                            geometry=geometry,
                             sources=sources,
-                            extra_materials = [m1,m2],
-                            geometry=geometry
-                            )
-        return sim
-
-    def test_eval(self):
-        sim = self.gen_sim()
+                            k_point=k_point)
+
+        h = mp.Harminv(mp.Hz, mp.Vector3(0.3718,-0.2076), fcen, df)
+        sim.run(mp.after_sources(h),
+                until_after_sources=200)
+
+        try:
+            for m in h.modes:
+                print("harminv:, {}, {}, {}".format(res,m.freq,m.Q))
+            freq = h.modes[0].freq
+        except:
+            print("No resonant modes found.")
+
+        sim.reset_meep()
+        return freq
+
+class TestMaterialGrid(unittest.TestCase):
+
+    def test_material_grid(self):
+        ## reference frequency computed using MaterialGrid at resolution = 300
+        freq_ref = 0.3068839373003908
+
+        res = [25, 50]
+        freq_matgrid = []
+        for r in res:
+            freq_matgrid.append(compute_resonant_mode(r))
+            ## verify that the resonant mode is approximately equivalent to
+            ## the reference value
+            self.assertAlmostEqual(freq_ref, freq_matgrid[-1], 2)
+
+        ## verify that the relative error is decreasing with increasing resolution
+        ## and is better than linear convergence
+        self.assertLess(abs(freq_matgrid[1]-freq_ref),abs(freq_matgrid[0]-freq_ref)/2)
+
 if __name__ == '__main__':
     unittest.main()