Python API for GDSII regions (NanoComp#518)

* Python API for GDSII regions * Update docs * Add docs for GDSII_vol * FieldsRegion uses full volume by default
bencbartlett · Sep 20, 2018 · 85346ef · 85346ef
1 parent 7df7bcc
commit 85346ef
Show file tree

Hide file tree

Showing 9 changed files with 115 additions and 96 deletions.
diff --git a/doc/docs/Python_Tutorials/GDSII_Import.md b/doc/docs/Python_Tutorials/GDSII_Import.md
@@ -34,12 +34,14 @@ As an alternative, the volume regions of the source and flux monitors could have
 The simulation script is [`coupler.py`](https://github.com/stevengj/meep/blob/master/python/examples/coupler.py).
 
 ```python
-import meep as mp
 import argparse
+import os
+import meep as mp
 
-resolution = 25 # pixels/um
 
-gdsII_file = 'coupler.gds'
+resolution = 25 # pixels/um
+examples_dir = os.path.realpath(os.path.dirname(__file__))
+gdsII_file = os.path.join(examples_dir, 'coupler.gds')
 CELL_LAYER = 0
 PORT1_LAYER = 1
 PORT2_LAYER = 2
@@ -65,19 +67,6 @@ lcen = 1.55
 fcen = 1/lcen
 df = 0.2*fcen
 
-# extract center and size of meep::volume
-def get_center_and_size(v):
-    rmin = v.get_min_corner()
-    rmax = v.get_max_corner()
-    v3rmin = mp.Vector3(rmin.x(),rmin.y(),rmin.z())
-    v3rmax = mp.Vector3(rmax.x(),rmax.y(),rmax.z())
-    if v.dim<mp.D3:
-      v3rmin.z = v3rmax.z=0
-    if v.dim<mp.D2:
-      v3rmin.y = v3rmax.y=0
-    center = 0.5*(v3rmin+v3rmax)
-    size = v3rmax-v3rmin
-    return (center,size)
 
 def main(args):
     d  = args.d
@@ -93,23 +82,23 @@ def main(args):
     upper_branch = mp.get_GDSII_prisms(silicon, gdsII_file, UPPER_BRANCH_LAYER, si_zmin, si_zmax)
     lower_branch = mp.get_GDSII_prisms(silicon, gdsII_file, LOWER_BRANCH_LAYER, si_zmin, si_zmax)
 
-    (cell_center,cell_size) = get_center_and_size(mp.get_GDSII_volume(gdsII_file,CELL_LAYER, cell_zmin, cell_zmax))
-    (p1_center,p1_size) = get_center_and_size(mp.get_GDSII_volume(gdsII_file,PORT1_LAYER, si_zmin, si_zmax))
-    (p2_center,p2_size) = get_center_and_size(mp.get_GDSII_volume(gdsII_file,PORT2_LAYER, si_zmin, si_zmax))
-    (p3_center,p3_size) = get_center_and_size(mp.get_GDSII_volume(gdsII_file,PORT3_LAYER, si_zmin, si_zmax))
-    (p4_center,p4_size) = get_center_and_size(mp.get_GDSII_volume(gdsII_file,PORT4_LAYER, si_zmin, si_zmax))
-    (src_center,src_size) = get_center_and_size(mp.get_GDSII_volume(gdsII_file,SOURCE_LAYER, si_zmin, si_zmax))
+    cell = mp.GDSII_vol(gdsII_file, CELL_LAYER, cell_zmin, cell_zmax)
+    p1 = mp.GDSII_vol(gdsII_file, PORT1_LAYER, si_zmin, si_zmax)
+    p2 = mp.GDSII_vol(gdsII_file, PORT2_LAYER, si_zmin, si_zmax)
+    p3 = mp.GDSII_vol(gdsII_file, PORT3_LAYER, si_zmin, si_zmax)
+    p4 = mp.GDSII_vol(gdsII_file, PORT4_LAYER, si_zmin, si_zmax)
+    src_vol = mp.GDSII_vol(gdsII_file, SOURCE_LAYER, si_zmin, si_zmax)
 
     # displace upper and lower branches of coupler (as well as source and flux regions)
     if d != default_d:
         delta_y = 0.5*(d-default_d)
         delta = mp.Vector3(y=delta_y)
-        p1_center += delta
-        p2_center -= delta
-        p3_center += delta
-        p4_center -= delta
-        src_center += delta
-        cell_size += 2*delta
+        p1.center += delta
+        p2.center -= delta
+        p3.center += delta
+        p4.center -= delta
+        src_vol.center += delta
+        cell.size += 2*delta
         for np in range(len(lower_branch)):
             lower_branch[np].center -= delta
             for nv in range(len(lower_branch[np].vertices)):
@@ -123,39 +112,39 @@ def main(args):
 
     if args.three_d:
         oxide_center = mp.Vector3(z=-0.5*t_oxide)
-        oxide_size = mp.Vector3(cell_size.x,cell_size.y,t_oxide)
+        oxide_size = mp.Vector3(cell.size.x,cell.size.y,t_oxide)
         oxide_layer = [mp.Block(material=oxide, center=oxide_center, size=oxide_size)]
         geometry = geometry+oxide_layer
 
     sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df),
-                                  size=src_size,
-                                  center=src_center,
+                                  size=src_vol.size,
+                                  center=src_vol.center,
                                   eig_match_freq=True)]
 
     sim = mp.Simulation(resolution=resolution,
-                        cell_size=cell_size,
+                        cell_size=cell.size,
                         boundary_layers=[mp.PML(dpml)],
                         sources=sources,
                         geometry=geometry)
 
-    p1_region = mp.FluxRegion(center=p1_center,size=p1_size)
+    p1_region = mp.FluxRegion(volume=p1)
     flux1 = sim.add_flux(fcen,0,1,p1_region)
-    p2_region = mp.FluxRegion(center=p2_center,size=p2_size)
+    p2_region = mp.FluxRegion(volume=p2)
     flux2 = sim.add_flux(fcen,0,1,p2_region)
-    p3_region = mp.FluxRegion(center=p3_center,size=p3_size)
+    p3_region = mp.FluxRegion(volume=p3)
     flux3 = sim.add_flux(fcen,0,1,p3_region)
-    p4_region = mp.FluxRegion(center=p4_center,size=p4_size)
+    p4_region = mp.FluxRegion(volume=p4)
     flux4 = sim.add_flux(fcen,0,1,p4_region)
 
-    sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ez,p3_center,1e-8))
+    sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ez,p3.center,1e-8))
 
     p1_flux = mp.get_fluxes(flux1)
     p2_flux = mp.get_fluxes(flux2)
     p3_flux = mp.get_fluxes(flux3)
     p4_flux = mp.get_fluxes(flux4)
 
     mp.master_printf("data:, {}, {}, {}, {}".format(d,-p2_flux[0]/p1_flux[0],p3_flux[0]/p1_flux[0],p4_flux[0]/p1_flux[0]))
-            
+
 if __name__ == '__main__':
     parser = argparse.ArgumentParser()
     parser.add_argument('-d',type=float,default=0.1,help='branch separation (default: 0.1 um)')
@@ -184,4 +173,4 @@ The results are plotted below. When the two waveguide branches are spaced suffic
 
 <center>
 ![](../images/coupler3D.png)
-</center>
+</center>
diff --git a/doc/docs/Python_Tutorials/Optical_Forces.md b/doc/docs/Python_Tutorials/Optical_Forces.md
@@ -114,8 +114,8 @@ sim.change_sources(new_sources)
 flx_reg = mp.FluxRegion(direction=mp.Z, center=mp.Vector3(), size=mp.Vector3(1.2 * (2 * a + s), 1.2 * a))
 wvg_pwr = sim.add_flux(f, 0, 1, flx_reg)
 
-frc_reg1 = mp.ForceRegion(mp.Vector3(0.5 * s), mp.X, weight=1.0, size=mp.Vector3(y=a))
-frc_reg2 = mp.ForceRegion(mp.Vector3(0.5 * s + a), mp.X, weight=-1.0, size=mp.Vector3(y=a))
+frc_reg1 = mp.ForceRegion(mp.Vector3(0.5 * s), direction=mp.X, weight=1.0, size=mp.Vector3(y=a))
+frc_reg2 = mp.ForceRegion(mp.Vector3(0.5 * s + a), direction=mp.X, weight=-1.0, size=mp.Vector3(y=a))
 wvg_force = sim.add_force(f, 0, 1, frc_reg1, frc_reg2)
 
 runtime = 5000
@@ -137,4 +137,4 @@ We run this simulation over a range of non-zero separation distances and compare
 
 <center>
 ![](../images/Waveguide_forces.png)
-</center>
+</center>
diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
@@ -769,6 +769,9 @@ Properties:
 **`weight` [`complex`]**
 —A weight factor to multiply the flux by when it is computed. Default is 1.0.
 
+**`volume` [`Volume`]**
+—A `meep.Volume` can be used to specify the flux region instead of a center and a size.
+
 Note that the flux is always computed in the *positive* coordinate direction, although this can effectively be flipped by using a `weight` of -1.0. This is useful, for example, if you want to compute the outward flux through a box, so that the sides of the box add instead of subtract.
 
 ### Volume
@@ -1031,6 +1034,10 @@ The direction of the force that you wish to compute (e.g. `X`, `Y`, etcetera). U
 —
 A weight factor to multiply the force by when it is computed. Default is 1.0.
 
+**`volume` [`Volume`]**
+—
+A `meep.Volume` can be used to specify the force region instead of a center and a size.
+
 In most circumstances, you should define a set of `ForceRegion`s whose union is a closed surface lying in vacuum and enclosing the object that is experiencing the force.
 
 **`Simulation.add_force(fcen, df, nfreq, ForceRegions...)`**
@@ -1209,6 +1216,14 @@ This feature is only available if Meep is built with [libGDSII](Build_From_Sourc
 —
 Returns a list of `GeometricObject`s with `material` (`mp.Medium`) on layer `layer` of a GDSII file `gdsii_filename`.
 
+**`mp.GDSII_vol(fname, layer, zmin, zmax)`**
+—
+Returns a `mp.Volume` read from a GDSII file `fname` on `layer` with `zmin` and `zmax`. This function is useful for creating a `FluxRegion` from a GDSII file as follows
+
+```python
+fr = mp.FluxRegion(volume=mp.GDSII_vol(fname, layer, zmin, zmax))
+```
+
 Run and Step Functions
 ----------------------
 

diff --git a/python/examples/coupler.py b/python/examples/coupler.py
@@ -1,9 +1,11 @@
-import meep as mp
 import argparse
+import os
+import meep as mp
 
-resolution = 25 # pixels/um
 
-gdsII_file = 'coupler.gds'
+resolution = 25 # pixels/um
+examples_dir = os.path.realpath(os.path.dirname(__file__))
+gdsII_file = os.path.join(examples_dir, 'coupler.gds')
 CELL_LAYER = 0
 PORT1_LAYER = 1
 PORT2_LAYER = 2
@@ -44,23 +46,23 @@ def main(args):
     upper_branch = mp.get_GDSII_prisms(silicon, gdsII_file, UPPER_BRANCH_LAYER, si_zmin, si_zmax)
     lower_branch = mp.get_GDSII_prisms(silicon, gdsII_file, LOWER_BRANCH_LAYER, si_zmin, si_zmax)
 
-    (cell_center,cell_size) = mp.get_center_and_size(mp.get_GDSII_volume(gdsII_file,CELL_LAYER, cell_zmin, cell_zmax))
-    (p1_center,p1_size) = mp.get_center_and_size(mp.get_GDSII_volume(gdsII_file,PORT1_LAYER, si_zmin, si_zmax))
-    (p2_center,p2_size) = mp.get_center_and_size(mp.get_GDSII_volume(gdsII_file,PORT2_LAYER, si_zmin, si_zmax))
-    (p3_center,p3_size) = mp.get_center_and_size(mp.get_GDSII_volume(gdsII_file,PORT3_LAYER, si_zmin, si_zmax))
-    (p4_center,p4_size) = mp.get_center_and_size(mp.get_GDSII_volume(gdsII_file,PORT4_LAYER, si_zmin, si_zmax))
-    (src_center,src_size) = mp.get_center_and_size(mp.get_GDSII_volume(gdsII_file,SOURCE_LAYER, si_zmin, si_zmax))
+    cell = mp.GDSII_vol(gdsII_file, CELL_LAYER, cell_zmin, cell_zmax)
+    p1 = mp.GDSII_vol(gdsII_file, PORT1_LAYER, si_zmin, si_zmax)
+    p2 = mp.GDSII_vol(gdsII_file, PORT2_LAYER, si_zmin, si_zmax)
+    p3 = mp.GDSII_vol(gdsII_file, PORT3_LAYER, si_zmin, si_zmax)
+    p4 = mp.GDSII_vol(gdsII_file, PORT4_LAYER, si_zmin, si_zmax)
+    src_vol = mp.GDSII_vol(gdsII_file, SOURCE_LAYER, si_zmin, si_zmax)
 
     # displace upper and lower branches of coupler (as well as source and flux regions)
     if d != default_d:
         delta_y = 0.5*(d-default_d)
         delta = mp.Vector3(y=delta_y)
-        p1_center += delta
-        p2_center -= delta
-        p3_center += delta
-        p4_center -= delta
-        src_center += delta
-        cell_size += 2*delta
+        p1.center += delta
+        p2.center -= delta
+        p3.center += delta
+        p4.center -= delta
+        src_vol.center += delta
+        cell.size += 2*delta
         for np in range(len(lower_branch)):
             lower_branch[np].center -= delta
             for nv in range(len(lower_branch[np].vertices)):
@@ -74,39 +76,39 @@ def main(args):
 
     if args.three_d:
         oxide_center = mp.Vector3(z=-0.5*t_oxide)
-        oxide_size = mp.Vector3(cell_size.x,cell_size.y,t_oxide)
+        oxide_size = mp.Vector3(cell.size.x,cell.size.y,t_oxide)
         oxide_layer = [mp.Block(material=oxide, center=oxide_center, size=oxide_size)]
         geometry = geometry+oxide_layer
 
     sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df),
-                                  size=src_size,
-                                  center=src_center,
+                                  size=src_vol.size,
+                                  center=src_vol.center,
                                   eig_match_freq=True)]
 
     sim = mp.Simulation(resolution=resolution,
-                        cell_size=cell_size,
+                        cell_size=cell.size,
                         boundary_layers=[mp.PML(dpml)],
                         sources=sources,
                         geometry=geometry)
 
-    p1_region = mp.FluxRegion(center=p1_center,size=p1_size)
+    p1_region = mp.FluxRegion(volume=p1)
     flux1 = sim.add_flux(fcen,0,1,p1_region)
-    p2_region = mp.FluxRegion(center=p2_center,size=p2_size)
+    p2_region = mp.FluxRegion(volume=p2)
     flux2 = sim.add_flux(fcen,0,1,p2_region)
-    p3_region = mp.FluxRegion(center=p3_center,size=p3_size)
+    p3_region = mp.FluxRegion(volume=p3)
     flux3 = sim.add_flux(fcen,0,1,p3_region)
-    p4_region = mp.FluxRegion(center=p4_center,size=p4_size)
+    p4_region = mp.FluxRegion(volume=p4)
     flux4 = sim.add_flux(fcen,0,1,p4_region)
 
-    sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ez,p3_center,1e-8))
+    sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ez,p3.center,1e-8))
 
     p1_flux = mp.get_fluxes(flux1)
     p2_flux = mp.get_fluxes(flux2)
     p3_flux = mp.get_fluxes(flux3)
     p4_flux = mp.get_fluxes(flux4)
 
     mp.master_printf("data:, {}, {}, {}, {}".format(d,-p2_flux[0]/p1_flux[0],p3_flux[0]/p1_flux[0],p4_flux[0]/p1_flux[0]))
-            
+
 if __name__ == '__main__':
     parser = argparse.ArgumentParser()
     parser.add_argument('-d',type=float,default=0.1,help='branch separation (default: 0.1 um)')

diff --git a/python/examples/parallel-wvgs-force.py b/python/examples/parallel-wvgs-force.py
@@ -70,8 +70,8 @@ def main(args):
     flx_reg = mp.FluxRegion(direction=mp.Z, center=mp.Vector3(), size=mp.Vector3(1.2 * (2 * a + s), 1.2 * a))
     wvg_pwr = sim.add_flux(f, 0, 1, flx_reg)
 
-    frc_reg1 = mp.ForceRegion(mp.Vector3(0.5 * s), mp.X, weight=1.0, size=mp.Vector3(y=a))
-    frc_reg2 = mp.ForceRegion(mp.Vector3(0.5 * s + a), mp.X, weight=-1.0, size=mp.Vector3(y=a))
+    frc_reg1 = mp.ForceRegion(mp.Vector3(0.5 * s), direction=mp.X, weight=1.0, size=mp.Vector3(y=a))
+    frc_reg2 = mp.ForceRegion(mp.Vector3(0.5 * s + a), direction=mp.X, weight=-1.0, size=mp.Vector3(y=a))
     wvg_force = sim.add_force(f, 0, 1, frc_reg1, frc_reg2)
 
     runtime = 5000

diff --git a/python/meep.i b/python/meep.i
@@ -1246,6 +1246,7 @@ py_eigenmode_data _get_eigenmode(meep::fields *f, double omega_src, meep::direct
         dft_ldos,
         display_progress,
         during_sources,
+        GDSII_vol,
         get_center_and_size,
         get_flux_freqs,
         get_fluxes,

diff --git a/python/simulation.py b/python/simulation.py
@@ -156,40 +156,37 @@ def __init__(self, center, size=Vector3(), dims=2, is_cylindrical=False):
 
 class FluxRegion(object):
 
-    def __init__(self, center, size=Vector3(), direction=mp.AUTOMATIC, weight=1.0):
-        self.center = center
-        self.size = size
-        self.direction = direction
-        self.weight = complex(weight)
-
-
-ModeRegion = FluxRegion
+    def __init__(self, center=None, size=Vector3(), direction=mp.AUTOMATIC, weight=1.0, volume=None):
+        if center is None and volume is None:
+            raise ValueError("Either center or volume required")
 
+        if volume:
+            self.center = volume.center
+            self.size = volume.size
+        else:
+            self.center = center
+            self.size = size
 
-class ForceRegion(object):
-
-    def __init__(self, center, direction, size=mp.Vector3(), weight=1.0):
-        self.center = center
         self.direction = direction
-        self.size = size
         self.weight = complex(weight)
 
 
-class Near2FarRegion(object):
-
-    def __init__(self, center, size=mp.Vector3(), direction=mp.AUTOMATIC, weight=1.0):
-        self.center = center
-        self.size = size
-        self.direction = direction
-        self.weight = complex(weight)
+ModeRegion = FluxRegion
+Near2FarRegion = FluxRegion
+ForceRegion = FluxRegion
 
 
 class FieldsRegion(object):
 
     def __init__(self, where=None, center=None, size=None):
+        if where:
+            self.center = where.center
+            self.size = where.size
+        else:
+            self.center = center
+            self.size = size
+
         self.where = where
-        self.center = center
-        self.size = size
 
 
 class DftObj(object):
@@ -2230,3 +2227,18 @@ def get_center_and_size(v):
     center = 0.5 * (v3rmin + v3rmax)
     size = v3rmax - v3rmin
     return center, size
+
+
+def GDSII_vol(fname, layer, zmin, zmax):
+    meep_vol = mp.get_GDSII_volume(fname, layer, zmin, zmax)
+    dims = meep_vol.dim + 1
+    is_cyl = False
+
+    if dims == 4:
+        # cylindrical
+        dims = 2
+        is_cyl = True
+
+    center, size = get_center_and_size(meep_vol)
+
+    return Volume(center, size, dims, is_cyl)
diff --git a/python/tests/force.py b/python/tests/force.py
@@ -20,7 +20,7 @@ def setUp(self):
                                  boundary_layers=[pml_layers],
                                  sources=[sources])
 
-        fr = mp.ForceRegion(mp.Vector3(y=1.27), mp.Y, size=mp.Vector3(4.38))
+        fr = mp.ForceRegion(mp.Vector3(y=1.27), direction=mp.Y, size=mp.Vector3(4.38))
         self.myforce = self.sim.add_force(fcen, 0, 1, fr)
 
     def test_force(self):