Python test for subpixel smoothing of prisms (#1135)

* Python test for subpixel smoothing of prisms

* add test involving marching squares algorithm

* whoops, wrong link to conda pack

doc/docs/FAQ.md

@@ -109,7 +109,7 @@ Note also that, as a consequence of the above analysis, ε must go to a positive
 
 ### Why are there strange peaks in my reflectance/transmittance spectrum when modeling planar or periodic structures?
 
-There are two possible explanations: (1) the simulation run time may be too short and your results have not sufficiently [converged](#checking-convergence), or (2) you may be using a higher-dimensional cell with multiple periods (a supercell) which introduces unwanted additional modes due to band folding. Modeling flat/planar structures typically requires a 1d cell and periodic structures a single unit cell in 2d/3d. For more details, see Section 4.6 ("Sources in Supercells") in [Chapter 4](http://arxiv.org/abs/arXiv:1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707). Note that a 1d cell must be along the $z$ direction with only the E<sub>x</sub> and H<sub>y</sub> field components permitted.
+There are two possible explanations: (1) the simulation run time may be too short or the resolution may be too low and thus your results have not sufficiently [converged](#checking-convergence), or (2) you may be using a higher-dimensional cell with multiple periods (a supercell) which introduces unwanted additional modes due to band folding. One indication that band folding is present is that small changes in the resolution produce large and unexpected changes in the frequency spectra. Modeling flat/planar structures typically requires a 1d cell and periodic structures a single unit cell in 2d/3d. For more details, see Section 4.6 ("Sources in Supercells") in [Chapter 4](http://arxiv.org/abs/arXiv:1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707). Note that a 1d cell must be along the $z$ direction with only the E<sub>x</sub> and H<sub>y</sub> field components permitted.
 
 ### How do I model the solar radiation spectrum?
 

doc/docs/Installation.md

@@ -123,6 +123,10 @@ print(mp.__version__)
 
 This will show something like `1.11.0-1-g415bc8eb` where the first three digits (`1.11.0`) refer to a stable tarball release, the following digit is the number of commits after this stable release, and the eight characters following the `g` in the final string refer to the commit hash.
 
+### Non-Networked Systems
+
+To install the PyMeep Conda package on a [non-networked system](https://docs.anaconda.com/anaconda/user-guide/tasks/install-packages/#installing-packages-on-a-non-networked-air-gapped-computer), using the bz2 tarball of the [official release](https://anaconda.org/conda-forge/pymeep/files) or [nightly build](https://anaconda.org/simpetus/pymeep/files) will *not* work without the dependencies. A possible workaround is [Conda-Pack](https://github.com/conda/conda-pack).
+
 Installation on Linux
 -------------------------
 

doc/docs/Python_User_Interface.md

@@ -1852,7 +1852,7 @@ with the following input parameters:
 
 + `omega`: optional frequency point over which the average eigenvalue of the dielectric and permeability tensors are evaluated (defaults to 0).
 
-For convenience, the following wrappers for `get_array` over the entire cell are available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, `get_tot_pwr()`, `get_Xfield()`, `get_Xfield_x()`, `get_Xfield_y()`, `get_Xfield_z()`, `get_Xfield_r()`, `get_Xfield_p()` where `X` is one of `h`, `b`, `e`, `d`, or `s`. The routines `get_Xfield_*` all return an array type consistent with the fields (real or complex). The routines `get_epsilon()` and `get_mu()` accept the optional omega parameter (defaults to 0).
+For convenience, the following wrappers for `get_array` over the entire cell are available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, `get_tot_pwr()`, `get_Xfield()`, `get_Xfield_x()`, `get_Xfield_y()`, `get_Xfield_z()`, `get_Xfield_r()`, `get_Xfield_p()` where `X` is one of `h`, `b`, `e`, `d`, or `s`. The routines `get_Xfield_*` all return an array type consistent with the fields (real or complex). The routines `get_epsilon()` and `get_mu()` accept the optional `omega` parameter (defaults to 0).
 
 **Note on array-slice dimensions:** The routines `get_epsilon`, `get_Xfield_z`, etc. use as default `size=meep.Simulation.fields.total_volume()` which for simulations involving Bloch-periodic boundaries (via `k_point`) will result in arrays that have slightly *different* dimensions than e.g. `get_array(center=meep.Vector3(), size=cell_size, component=meep.Dielectric`, etc. (i.e., the slice spans the entire cell volume `cell_size`). Neither of these approaches is "wrong", they are just slightly different methods of fetching the boundaries. The key point is that if you pass the same value for the `size` parameter, or use the default, the slicing routines always give you the same-size array for all components. You should *not* try to predict the exact size of these arrays; rather, you should simply rely on Meep's output.
 

python/Makefile.am

@@ -65,6 +65,7 @@ TESTS =                                   \
     $(TEST_DIR)/n2f_periodic.py           \
     $(TEST_DIR)/oblique_source.py         \
     $(TEST_DIR)/physical.py               \
+    $(TEST_DIR)/prism.py                  \
     $(TEST_DIR)/pw_source.py              \
     $(TEST_DIR)/refl_angular.py           \
     $(TEST_DIR)/ring.py                   \

python/tests/prism.py

@@ -0,0 +1,125 @@
+from __future__ import division
+
+import os
+import unittest
+import numpy as np
+import meep as mp
+
+class TestPrism(unittest.TestCase):
+
+  def prism_marching_squares(self,npts):
+    resolution = 50
+
+    cell = mp.Vector3(3,3)
+
+    data_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), 'data'))
+    vertices_file = os.path.join(data_dir, 'prism_vertices.npz')
+    vertices_obj = np.load(vertices_file)
+
+    ## prism vertices precomputed for a circle of radius 1.0 using
+    ## marching squares algorithm of skimage.measure.find_contours
+    ## ref: https://github.com/NanoComp/meep/issues/1060
+    vertices_data = vertices_obj["N{}".format(npts)]
+    vertices = [mp.Vector3(v[0],v[1],0) for v in vertices_data]
+
+    geometry = [mp.Prism(vertices,
+                         height=mp.inf,
+                         material=mp.Medium(epsilon=12))]
+
+    sim = mp.Simulation(cell_size=cell,
+                        geometry=geometry,
+                        resolution=resolution)
+
+    sim.init_sim()
+
+    prism_eps = sim.integrate_field_function([mp.Dielectric], lambda r,eps: eps)
+
+    sim.reset_meep()
+
+    geometry = [mp.Cylinder(radius=1.0,
+                            center=mp.Vector3(),
+                            height=mp.inf,
+                            material=mp.Medium(epsilon=12))]
+
+    sim = mp.Simulation(cell_size=cell,
+                        geometry=geometry,
+                        resolution=resolution)
+
+    sim.init_sim()
+
+    cyl_eps = sim.integrate_field_function([mp.Dielectric], lambda r,eps: eps)
+
+    print("epsilon-sum:, {} (prism-msq), {} (cylinder), {} (relative error)".format(abs(prism_eps),abs(cyl_eps),abs((prism_eps-cyl_eps)/cyl_eps)))
+
+    return abs((prism_eps-cyl_eps)/cyl_eps)
+
+
+  def prism_circle(self,npts,r):
+    resolution = 50
+
+    cell = mp.Vector3(3,3)
+
+    ### prism vertices computed as npts equally-spaced points
+    ### along the circumference of a circle with radius r
+    angles = 2*np.pi/npts * np.arange(npts)
+    vertices = [mp.Vector3(r*np.cos(ang),r*np.sin(ang)) for ang in angles]
+    geometry = [mp.Prism(vertices,
+                         height=mp.inf,
+                         material=mp.Medium(epsilon=12))]
+
+    sim = mp.Simulation(cell_size=cell,
+                        geometry=geometry,
+                        resolution=resolution)
+
+    sim.init_sim()
+
+    prism_eps = sim.integrate_field_function([mp.Dielectric], lambda r,eps: eps)
+
+    sim.reset_meep()
+
+    geometry = [mp.Cylinder(radius=r,
+                            center=mp.Vector3(),
+                            height=mp.inf,
+                            material=mp.Medium(epsilon=12))]
+
+    sim = mp.Simulation(cell_size=cell,
+                        geometry=geometry,
+                        resolution=resolution)
+
+    sim.init_sim()
+
+    cyl_eps = sim.integrate_field_function([mp.Dielectric], lambda r,eps: eps)
+
+    print("epsilon-sum:, {} (prism-cyl), {} (cylinder), {} (relative error)".format(abs(prism_eps),abs(cyl_eps),abs((prism_eps-cyl_eps)/cyl_eps)))
+
+    return abs((prism_eps-cyl_eps)/cyl_eps)
+
+
+  def test_prism(self):
+    print("Testing Prism object using marching squares algorithm...")
+    d = [self.prism_marching_squares(92),
+         self.prism_marching_squares(192),
+         self.prism_marching_squares(392)]
+
+    self.assertLess(d[1],d[0])
+    self.assertLess(d[2],d[1])
+
+    print("Testing Prism object using circle formula...")
+    r = 1.0458710786934182
+    d = [self.prism_circle(51,r),
+         self.prism_circle(101,r),
+         self.prism_circle(201,r)]
+
+    self.assertLess(d[1],d[0])
+    self.assertLess(d[2],d[1])
+
+    r = 1.2896871096581341
+    d = [self.prism_circle(31,r),
+         self.prism_circle(61,r),
+         self.prism_circle(121,r)]
+
+    self.assertLess(d[1],d[0])
+    self.assertLess(d[2],d[1])
+
+if __name__ == '__main__':
+  unittest.main()