diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index f7ef5c1b7..663c88076 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/doc/docs/Python_User_Interface.md
@@ -4156,6 +4156,89 @@ media](FAQ.md#why-does-my-simulation-diverge-if-0).
+---
+
+
+### MaterialGrid
+
+```python
+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).
+
+
+
+
+
+
+
+
+
+```python
+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).
+
+
+
+
+
+
+
+
+
+```python
+def update_parameters(self, x):
+```
+
+
+
+Reset the design grid `design_parameters` to `x`.
+
+
+
+
+
---
diff --git a/doc/docs/Python_User_Interface.md.in b/doc/docs/Python_User_Interface.md.in
index 236117d35..f77f50ea9 100644
--- a/doc/docs/Python_User_Interface.md.in
+++ b/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] @@
diff --git a/python/geom.py b/python/geom.py
index cc820a5b3..9a68965ce 100755
--- a/python/geom.py
+++ b/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)
diff --git a/python/tests/material_grid.py b/python/tests/material_grid.py
index 94caebf90..8a7b3fd03 100644
--- a/python/tests/material_grid.py
+++ b/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()