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).
@@ -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.
-
```python
-def update_parameters(self, x):
+def update_weights(self, x):
```
-Reset the design grid `design_parameters` to `x`.
+Reset the `weights` to `x`.
diff --git a/libpympb/pympb.hpp b/libpympb/pympb.hpp
index 76499cc3c..e57ece0ce 100644
--- a/libpympb/pympb.hpp
+++ b/libpympb/pympb.hpp
@@ -5,7 +5,7 @@
#include "ctlgeom.h"
#include "mpb.h"
-#include "meepgeom.hpp"
+#include "../src/meepgeom.hpp"
namespace py_mpb {
diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py
index ec77ee761..dcb45d09b 100644
--- a/python/adjoint/optimization_problem.py
+++ b/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
# -------------------------------------------- #
diff --git a/python/examples/adjoint_optimization/01-Introduction.ipynb b/python/examples/adjoint_optimization/01-Introduction.ipynb
index 3aa445414..48668958e 100644
--- a/python/examples/adjoint_optimization/01-Introduction.ipynb
+++ b/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",
diff --git a/python/examples/adjoint_optimization/02-Waveguide_Bend.ipynb b/python/examples/adjoint_optimization/02-Waveguide_Bend.ipynb
index 2bfdcd15a..595106841 100644
--- a/python/examples/adjoint_optimization/02-Waveguide_Bend.ipynb
+++ b/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",
diff --git a/python/examples/adjoint_optimization/03-Filtered_Waveguide_Bend.ipynb b/python/examples/adjoint_optimization/03-Filtered_Waveguide_Bend.ipynb
index a1f1e2731..6a68a85a4 100644
--- a/python/examples/adjoint_optimization/03-Filtered_Waveguide_Bend.ipynb
+++ b/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."
]
},
{
diff --git a/python/examples/adjoint_optimization/04-Splitter.ipynb b/python/examples/adjoint_optimization/04-Splitter.ipynb
index f2fc13612..5bcf768ae 100644
--- a/python/examples/adjoint_optimization/04-Splitter.ipynb
+++ b/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)))"
]
},
diff --git a/python/examples/adjoint_optimization/05-Minimax.ipynb b/python/examples/adjoint_optimization/05-Minimax.ipynb
index d32494fdf..4331dffd3 100644
--- a/python/examples/adjoint_optimization/05-Minimax.ipynb
+++ b/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",
diff --git a/python/examples/adjoint_optimization/Near2Far-Optimization-with-Epigraph-Formulation.ipynb b/python/examples/adjoint_optimization/Near2Far-Optimization-with-Epigraph-Formulation.ipynb
index c42bc8af3..ce5b485f8 100644
--- a/python/examples/adjoint_optimization/Near2Far-Optimization-with-Epigraph-Formulation.ipynb
+++ b/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",
diff --git a/python/geom.py b/python/geom.py
index 9a68965ce..72682c727 100755
--- a/python/geom.py
+++ b/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):
"""
diff --git a/python/meep.i b/python/meep.i
index 3bf36a9dd..fcece4087 100644
--- a/python/meep.i
+++ b/python/meep.i
@@ -744,7 +744,7 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c,
if (((material_data *)$1.material)->medium.H_susceptibilities.items) {
delete[] ((material_data *)$1.material)->medium.H_susceptibilities.items;
}
- delete[] ((material_data *)$1.material)->design_parameters;
+ delete[] ((material_data *)$1.material)->weights;
delete[] ((material_data *)$1.material)->epsilon_data;
delete (material_data *)$1.material;
geometric_object_destroy($1);
@@ -786,7 +786,7 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c,
delete[] ((material_data *)$1.items[i].material)->medium.H_susceptibilities.items;
}
delete[] ((material_data *)$1.items[i].material)->epsilon_data;
- delete[] ((material_data *)$1.items[i].material)->design_parameters;
+ delete[] ((material_data *)$1.items[i].material)->weights;
delete (material_data *)$1.items[i].material;
geometric_object_destroy($1.items[i]);
}
@@ -978,7 +978,7 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
if ($1->medium.H_susceptibilities.items) {
delete[] $1->medium.H_susceptibilities.items;
}
- delete[] $1->design_parameters;
+ delete[] $1->weights;
delete[] $1->epsilon_data;
delete $1;
}
@@ -1351,7 +1351,7 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
if ($1.items[i]->medium.H_susceptibilities.items) {
delete[] $1.items[i]->medium.H_susceptibilities.items;
}
- delete[] $1.items[i]->design_parameters;
+ delete[] $1.items[i]->weights;
delete[] $1.items[i]->epsilon_data;
}
delete[] $1.items;
diff --git a/python/solver.py b/python/solver.py
index 38f31bf5f..acd967ede 100644
--- a/python/solver.py
+++ b/python/solver.py
@@ -22,7 +22,7 @@
U_MIN = 0
U_PROD = 1
-U_SUM = 2
+U_MEAN = 2
class MPBArray(np.ndarray):
diff --git a/python/tests/adjoint_solver.py b/python/tests/adjoint_solver.py
index ff61f915b..8d0df40ea 100644
--- a/python/tests/adjoint_solver.py
+++ b/python/tests/adjoint_solver.py
@@ -56,8 +56,8 @@ def forward_simulation(design_params,mon_type):
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
mp.air,
silicon,
- design_parameters=design_params.reshape(Nx,Ny),
- grid_type='U_SUM')
+ weights=design_params.reshape(Nx,Ny),
+ grid_type='U_MEAN')
matgrid_geometry = [mp.Block(center=mp.Vector3(),
size=mp.Vector3(design_shape.x,design_shape.y,0),
@@ -109,7 +109,7 @@ def adjoint_solver(design_params, mon_type):
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
mp.air,
silicon,
- design_parameters=np.ones((Nx,Ny)))
+ weights=np.ones((Nx,Ny)))
matgrid_region = mpa.DesignRegion(matgrid,
volume=mp.Volume(center=mp.Vector3(),
diff --git a/python/tests/material_grid.py b/python/tests/material_grid.py
index 8a7b3fd03..3b0fbdcb3 100644
--- a/python/tests/material_grid.py
+++ b/python/tests/material_grid.py
@@ -31,7 +31,7 @@ def compute_resonant_mode(res):
matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny),
mp.air,
mp.Medium(index=3.5),
- design_parameters=filtered_design_params,
+ weights=filtered_design_params,
do_averaging=True,
beta=1000,
eta=0.5)
@@ -75,8 +75,8 @@ def test_material_grid(self):
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)
+ ## and is better than linear convergence because of subpixel smoothing
+ self.assertLess(abs(freq_matgrid[1]-freq_ref)*(res[1]**2)/2,abs(freq_matgrid[0]-freq_ref)*(res[0]**2))
if __name__ == '__main__':
unittest.main()
diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp
index 04154ecd7..f0606d2d8 100644
--- a/python/typemap_utils.cpp
+++ b/python/typemap_utils.cpp
@@ -454,7 +454,7 @@ static int pymaterial_grid_to_material_grid(PyObject *po, material_data *md) {
switch (gt_enum) {
case 0: md->material_grid_kinds = material_data::U_MIN; break;
case 1: md->material_grid_kinds = material_data::U_PROD; break;
- case 2: md->material_grid_kinds = material_data::U_SUM; break;
+ case 2: md->material_grid_kinds = material_data::U_MEAN; break;
case 3: md->material_grid_kinds = material_data::U_DEFAULT; break;
default: meep::abort("Invalid material grid enumeration code: %d.\n", gt_enum);
}
@@ -474,15 +474,15 @@ static int pymaterial_grid_to_material_grid(PyObject *po, material_data *md) {
meep::abort("MaterialGrid medium2 failed to init.");
}
- // Initialize design parameters
- PyObject *po_dp = PyObject_GetAttrString(po, "design_parameters");
+ // Initialize weights
+ PyObject *po_dp = PyObject_GetAttrString(po, "weights");
PyArrayObject *pao = (PyArrayObject *)po_dp;
- if (!PyArray_Check(pao)) { meep::abort("MaterialGrid design_parameters failed to init."); }
+ if (!PyArray_Check(pao)) { meep::abort("MaterialGrid weights failed to init."); }
if (!PyArray_ISCARRAY(pao)) {
- meep::abort("Numpy array design_parameters must be C-style contiguous.");
+ meep::abort("Numpy array weights must be C-style contiguous.");
}
- md->design_parameters = new realnum[PyArray_SIZE(pao)];
- memcpy(md->design_parameters, (realnum *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(realnum));
+ md->weights = new realnum[PyArray_SIZE(pao)];
+ memcpy(md->weights, (realnum *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(realnum));
// if needed, combine sus structs to main object
PyObject *py_e_sus_m1 = PyObject_GetAttrString(po_medium1, "E_susceptibilities");
diff --git a/src/material_data.hpp b/src/material_data.hpp
index f6ddb36d3..104246b42 100644
--- a/src/material_data.hpp
+++ b/src/material_data.hpp
@@ -187,7 +187,7 @@ struct material_data {
// these fields used only if which_subclass==MATERIAL_GRID
vector3 grid_size;
- meep::realnum *design_parameters;
+ meep::realnum *weights;
medium_struct medium_1;
medium_struct medium_2;
meep::realnum beta;
@@ -206,7 +206,7 @@ struct material_data {
stalling convergence, although we try to avoid this by making the
minimum u = 1e-4 instead of 0.
- For U_SUM: The gradient is divided by the number of overlapping grids.
+ For U_MEAN: The gradient is divided by the number of overlapping grids.
This doesn't have the property that u=0 in one grid makes the total
u=0, unfortunately, which is desirable if u=0 indicates "drilled holes".
@@ -214,11 +214,11 @@ struct material_data {
the object on top always wins and everything underneath is ignored.
Specifically, that means that u = the top material grid value at that point.
*/
- enum { U_MIN = 0, U_PROD = 1, U_SUM = 2, U_DEFAULT = 3 } material_grid_kinds;
+ enum { U_MIN = 0, U_PROD = 1, U_MEAN = 2, U_DEFAULT = 3 } material_grid_kinds;
material_data()
: which_subclass(MEDIUM), medium(), user_func(NULL), user_data(NULL), epsilon_data(NULL),
- design_parameters(NULL), medium_1(), medium_2() {
+ weights(NULL), medium_1(), medium_2() {
epsilon_dims[0] = 0;
epsilon_dims[1] = 0;
epsilon_dims[2] = 0;
@@ -247,7 +247,7 @@ material_type make_user_material(user_material_func user_func, void *user_data);
material_type make_file_material(char *epsilon_input_file);
material_type make_material_grid(bool do_averaging, double beta, double eta);
void read_epsilon_file(const char *eps_input_file);
-void update_design_parameters(material_type matgrid, double *design_parameters);
+void update_weights(material_type matgrid, double *weights);
}; // namespace meep_geom
diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index fc2430553..d0690729a 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -312,7 +312,7 @@ meep::realnum material_grid_val(vector3 p, material_data *md) {
if (!is_material_grid(md)) { meep::abort("Invalid material grid detected.\n"); }
meep::realnum u;
- u = meep::linear_interpolate(p.x, p.y, p.z, md->design_parameters, md->grid_size.x,
+ u = meep::linear_interpolate(p.x, p.y, p.z, md->weights, md->grid_size.x,
md->grid_size.y, md->grid_size.z, 1);
return u;
}
@@ -354,7 +354,7 @@ meep::realnum matgrid_val(vector3 p, geom_box_tree tp, int oi, material_data *md
? umin
: (md->material_grid_kinds == material_data::U_PROD
? uprod
- : (md->material_grid_kinds == material_data::U_SUM ? usum / matgrid_val_count
+ : (md->material_grid_kinds == material_data::U_MEAN ? usum / matgrid_val_count
: udefault)));
// project interpolated grid point
@@ -397,7 +397,7 @@ static void interp_tensors(vector3 diag_in_1, vector3 offdiag_in_1, vector3 diag
void epsilon_material_grid(material_data *md, meep::realnum u) {
// NOTE: assume p lies on normalized grid within (0,1)
- if (!(md->design_parameters)) meep::abort("material params were not initialized!");
+ if (!(md->weights)) meep::abort("material params were not initialized!");
medium_struct *mm = &(md->medium);
medium_struct *m1 = &(md->medium_1);
@@ -1819,9 +1819,9 @@ material_type make_material_grid(bool do_averaging, double beta, double eta) {
return md;
}
-void update_design_parameters(material_type matgrid, double *design_parameters) {
+void update_weights(material_type matgrid, double *weights) {
int N = matgrid->grid_size.x * matgrid->grid_size.y * matgrid->grid_size.z;
- memcpy(matgrid->design_parameters, design_parameters, N * sizeof(meep::realnum));
+ memcpy(matgrid->weights, weights, N * sizeof(meep::realnum));
}
/******************************************************************************/
@@ -2461,7 +2461,7 @@ void material_grids_addgradient_point(meep::realnum *v, std::complex
fie
return; /* no material grids at this point */
// Calculate the number of material grids if we are averaging values
- if ((tp) && ((kind = mg->material_grid_kinds) == material_data::U_SUM)) {
+ if ((tp) && ((kind = mg->material_grid_kinds) == material_data::U_MEAN)) {
int matgrid_val_count = 0;
geom_box_tree tp_sum;
tp_sum = geom_tree_search(p, geometry_tree, &ois);
@@ -2489,7 +2489,7 @@ void material_grids_addgradient_point(meep::realnum *v, std::complex fie
it's up to the user to only have one unique design grid in this volume.*/
vector3 sz = mg->grid_size;
meep::realnum *vcur = v, *ucur;
- ucur = mg->design_parameters;
+ ucur = mg->weights;
uval = matgrid_val(p, tp, oi, mg);
do {
vector3 pb = to_geom_box_coords(p, &tp->objects[oi]);
@@ -2506,7 +2506,7 @@ void material_grids_addgradient_point(meep::realnum *v, std::complex fie
vector3 pb = to_geom_box_coords(p, &tp->objects[oi]);
vector3 sz = mg->grid_size;
meep::realnum *vcur = v, *ucur;
- ucur = mg->design_parameters;
+ ucur = mg->weights;
uval = matgrid_val(p, tp, oi, mg);
add_interpolate_weights(pb.x, pb.y, pb.z, vcur, sz.x, sz.y, sz.z, 1,
get_material_gradient(uval, fields_a, fields_f, freq, mg, field_dir) *
diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp
index e25f50846..cfb0bc3a7 100644
--- a/src/meepgeom.hpp
+++ b/src/meepgeom.hpp
@@ -201,7 +201,7 @@ geom_box gv2box(const meep::volume &v);
/***************************************************************/
// material grid functions
/***************************************************************/
-void update_design_parameters(material_type matgrid, double *design_parameters);
+void update_weights(material_type matgrid, double *weights);
meep::realnum matgrid_val(vector3 p, geom_box_tree tp, int oi, material_data *md);
meep::realnum material_grid_val(vector3 p, material_data *md);
geom_box_tree calculate_tree(const meep::volume &v, geometric_object_list g);