additional improvements for MaterialGrid (NanoComp#1512)

doc/docs/Python_User_Interface.md

@@ -4168,8 +4168,8 @@ 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).
+`material` argument of a [`Block`](#block) geometric object or the `default_material` argument of the
+[`Simulation`](#Simulation) constructor (similar to a material function).
 
 </div>
 
@@ -4184,7 +4184,7 @@ def __init__(self,
              grid_size,
              medium1,
              medium2,
-             design_parameters=None,
+             weights=None,
              grid_type='U_DEFAULT',
              do_averaging=False,
              beta=0,
@@ -4195,45 +4195,46 @@ def __init__(self,
 
 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.
+The input are two materials `medium1` and `medium2` along with a weight function $u(x)$ which is defined on the Cartesian grid points
+by the NumPy array `weights` of size `grid_size` (a 3-tuple or `Vector3` of integers). Elements of the `weights` array must be in the
+range [0,1] where 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. The
+materials are bilinearly interpolated from the Cartesian grid points to Meep's [Yee grid](Yee_Lattice.md).
 
-[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
+For improving accuracy, [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 weights $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$ where $\beta=\infty$ gives an unsmoothed, discontinuous interface. The projection operator is $(\tanh(\beta\times\eta)
++\tanh(\beta\times(u-\eta)))/(\tanh(\beta\times\eta)+\tanh(\beta\times(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.
+defining a *discontinuous* function from 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.
+can be specified using the `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).
+objects (with no intervening objects), any overlapping points are computed using the method `grid_type` which is one of
+`"U_MIN"` (minimum of the overlapping grid values), `"U_PROD"` (product), `"U_MEAN"` (mean), `"U_DEFAULT"`
+(topmost material grid). In general, these `"U_*"` options allow you to combine any material grids that overlap
+in space with no intervening objects.
 
 </div>
 
 </div>
 
-<a id="MaterialGrid.update_parameters"></a>
+<a id="MaterialGrid.update_weights"></a>
 
 <div class="class_members" markdown="1">
 
 ```python
-def update_parameters(self, x):
+def update_weights(self, x):
 ```
 
 <div class="method_docstring" markdown="1">
 
-Reset the design grid `design_parameters` to `x`.
+Reset the `weights` to `x`.
 
 </div>
 

libpympb/pympb.hpp

@@ -5,7 +5,7 @@
 
 #include "ctlgeom.h"
 #include "mpb.h"
-#include "meepgeom.hpp"
+#include "../src/meepgeom.hpp"
 
 namespace py_mpb {
 

python/adjoint/optimization_problem.py

@@ -16,7 +16,7 @@ def __init__(self,design_parameters,volume=None, size=None, center=mp.Vector3(),
         self.num_design_params=design_parameters.num_params
         self.MaterialGrid=MaterialGrid
     def update_design_parameters(self,design_parameters):
-        self.design_parameters.update_parameters(design_parameters)
+        self.design_parameters.update_weights(design_parameters)
     def get_gradient(self,sim,fields_a,fields_f,frequencies):
         for c in range(3):
             fields_a[c] = fields_a[c].flatten(order='C')
@@ -310,7 +310,7 @@ def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,fi
         for k in fd_gradient_idx:
 
             b0 = np.ones((self.num_design_params[design_variables_idx],))
-            b0[:] = (self.design_regions[design_variables_idx].design_parameters.design_parameters)
+            b0[:] = (self.design_regions[design_variables_idx].design_parameters.weights)
             # -------------------------------------------- #
             # left function evaluation
             # -------------------------------------------- #

python/examples/adjoint_optimization/01-Introduction.ipynb

@@ -110,7 +110,7 @@
     "Nx = design_region_resolution\n",
     "Ny = design_region_resolution\n",
     "\n",
-    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_SUM')\n",
+    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')\n",
     "design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(1, 1, 0)))\n",
     "\n",
     "\n",

python/examples/adjoint_optimization/02-Waveguide_Bend.ipynb

@@ -75,7 +75,7 @@
     "Nx = design_region_resolution\n",
     "Ny = design_region_resolution\n",
     "\n",
-    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_SUM')\n",
+    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')\n",
     "design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(1, 1, 0)))\n",
     "\n",
     "\n",

python/examples/adjoint_optimization/03-Filtered_Waveguide_Bend.ipynb

@@ -139,7 +139,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Let's set up our design region. The defauly grid type behaves just like any other geometry object -- the last object in the tree overides any objects underneath it. In some cases, it's important to *average* overlapping design regions, i.e. to enforce particular symmetries. To enable this capability, we'll set our design variables to use the `U_SUM` grid type."
+    "Let's set up our design region. The defauly grid type behaves just like any other geometry object -- the last object in the tree overides any objects underneath it. In some cases, it's important to *average* overlapping design regions, i.e. to enforce particular symmetries. To enable this capability, we'll set our design variables to use the `U_MEAN` grid type."
    ]
   },
   {
@@ -151,7 +151,7 @@
     "Nx = int(design_region_resolution*design_region_width)\n",
     "Ny = int(design_region_resolution*design_region_height)\n",
     "\n",
-    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_SUM')\n",
+    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')\n",
     "design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(design_region_width, design_region_height, 0)))"
    ]
   },
@@ -215,7 +215,7 @@
    "source": [
     "Let's build the geometry and simulation objects. We can force a rotational symmetry by adding extra blocks with the same design variables, but different basis vectors (`e1`, `e2`, etc.) In our case, we need one additional block object rotated by 90 degrees.\n",
     "\n",
-    "Because our design variables are using the `U_SUM` scheme, the adjoint solver will search for all of the material grids at each point in our design region and average the overlapping design variables. This way, we can enforce our symmetry constraint."
+    "Because our design variables are using the `U_MEAN` scheme, the adjoint solver will search for all of the material grids at each point in our design region and average the overlapping design variables. This way, we can enforce our symmetry constraint."
    ]
   },
   {

python/examples/adjoint_optimization/04-Splitter.ipynb

@@ -122,7 +122,7 @@
     "Nx = int(design_region_resolution*design_region_width)\n",
     "Ny = int(design_region_resolution*design_region_height)\n",
     "\n",
-    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_SUM')\n",
+    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')\n",
     "design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(design_region_width, design_region_height, 0)))"
    ]
   },

python/examples/adjoint_optimization/05-Minimax.ipynb

@@ -108,7 +108,7 @@
     "Nx = int(design_region_resolution*design_region_width)\n",
     "Ny = int(design_region_resolution*design_region_height)\n",
     "\n",
-    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_SUM')\n",
+    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')\n",
     "design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(design_region_width, design_region_height, 0)))\n",
     "\n",
     "x_g = np.linspace(-design_region_width/2,design_region_width/2,Nx)\n",

python/examples/adjoint_optimization/Near2Far-Optimization-with-Epigraph-Formulation.ipynb

@@ -79,7 +79,7 @@
     "Nx = int(design_region_resolution*design_region_width)\n",
     "Ny = int(design_region_resolution*design_region_height)\n",
     "\n",
-    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_SUM')\n",
+    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')\n",
     "design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(design_region_width, design_region_height, 0)))\n",
     "\n",
     "def mapping(x,eta,beta):\n",

python/geom.py

@@ -525,36 +525,37 @@ def _get_epsmu(self, diag, offdiag, susceptibilities, conductivity_diag, conduct
 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).
+    `material` argument of a [`Block`](#block) geometric object or the `default_material` argument of the
+    [`Simulation`](#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):
+    def __init__(self,grid_size,medium1,medium2,weights=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.
+        The input are two materials `medium1` and `medium2` along with a weight function $u(x)$ which is defined on the Cartesian grid points
+        by the NumPy array `weights` of size `grid_size` (a 3-tuple or `Vector3` of integers). Elements of the `weights` array must be in the
+        range [0,1] where 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. The
+        materials are bilinearly interpolated from the Cartesian grid points to Meep's [Yee grid](Yee_Lattice.md).
 
-        [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
+        For improving accuracy, [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 weights $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$ where $\\beta=\\infty$ gives an unsmoothed, discontinuous interface. The projection operator is $(\\tanh(\\beta\\times\\eta)
+        +\\tanh(\\beta\\times(u-\\eta)))/(\\tanh(\\beta\\times\\eta)+\\tanh(\\beta\\times(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.
+        defining a *discontinuous* function from 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.
+        can be specified using the `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).
+        objects (with no intervening objects), any overlapping points are computed using the method `grid_type` which is one of
+        `"U_MIN"` (minimum of the overlapping grid values), `"U_PROD"` (product), `"U_MEAN"` (mean), `"U_DEFAULT"`
+        (topmost material grid). In general, these `"U_*"` options allow you to combine any material grids that overlap
+        in space with no intervening objects.
         """
         self.grid_size = mp.Vector3(*grid_size)
         self.medium1 = medium1
@@ -571,31 +572,31 @@ def isclose(a, b, rel_tol=1e-09, abs_tol=0.0):
         self.do_averaging = do_averaging
         self.beta = beta
         self.eta = eta
-        if design_parameters is None:
-            self.design_parameters = np.zeros((self.num_params,))
-        elif design_parameters.size != self.num_params:
-            raise ValueError("design_parameters of shape {} do not match user specified grid dimension: {}".format(design_parameters.size,self.grid_size))
+        if weights is None:
+            self.weights = np.zeros((self.num_params,))
+        elif weights.size != self.num_params:
+            raise ValueError("weights of shape {} do not match user specified grid dimension: {}".format(weights.size,self.grid_size))
         else:
-            self.design_parameters = design_parameters.flatten().astype(np.float64)
+            self.weights = weights.flatten().astype(np.float64)
 
         grid_type_dict = {
             "U_MIN":0,
             "U_PROD":1,
-            "U_SUM":2,
+            "U_MEAN":2,
             "U_DEFAULT":3
         }
         if grid_type not in grid_type_dict:
-            raise ValueError("Invalid grid_type: {}. Must be either U_MIN, U_PROD, U_SUM, or U_DEFAULT".format(grid_type_dict))
+            raise ValueError("Invalid grid_type: {}. Must be either U_MIN, U_PROD, U_MEAN, or U_DEFAULT".format(grid_type_dict))
         self.grid_type = grid_type_dict[grid_type]
 
         self.swigobj = None
-    def update_parameters(self,x):
+    def update_weights(self,x):
         """
-        Reset the design grid `design_parameters` to `x`.
+        Reset the `weights` 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)
+            raise ValueError("weights of shape {} do not match user specified grid dimension: {}".format(self.weights.size,self.grid_size))
+        self.weights[:]=x.flatten().astype(np.float64)
 
 class Susceptibility(object):
     """