WIP:MaterialGrid

python/geom.py

@@ -278,6 +278,43 @@ def _get_epsmu(self, diag, offdiag, susceptibilities, conductivity_diag, conduct
         # Convert list matrix to 3D numpy array size [freqs,3,3]
         return np.squeeze(epsmu)
 
+class MaterialGrid(object):
+    def __init__(self,grid_size,medium1,medium2,design_parameters=None,grid_type="U_DEFAULT"):
+        self.grid_size = grid_size
+        self.medium1 = medium1
+        self.medium2 = medium2
+        if self.grid_size.x == 0:
+            self.grid_size.x = 1
+        elif self.grid_size.y == 0:
+            self.grid_size.y = 1
+        elif self.grid_size.z == 0:
+            self.grid_size.z = 1
+        self.num_params=self.grid_size.x*self.grid_size.y*self.grid_size.z
+
+        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,grid_size))
+        else:
+            self.design_parameters = design_parameters.flatten().astype(np.float64)
+
+        grid_type_dict = {
+            "U_MIN":0,
+            "U_PROD":1,
+            "U_SUM":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))
+        self.grid_type = grid_type_dict[grid_type]
+
+        self.swigobj = None
+    def update_parameters(self,x):
+        if x.size != self.num_params:
+            raise ValueError("design_parameters of shape {} do not match user specified grid dimension: {}".format(design_parameters.size,grid_size))
+        self.design_parameters=x.flatten().astype(np.float64)
+    def get_gradient(self,fields,grid):
+        return
 
 class Susceptibility(object):
 

python/meep.i

@@ -679,6 +679,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)->epsilon_data;
         delete (material_data *)$1.material;
         geometric_object_destroy($1);
@@ -720,6 +721,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;
         geometric_object_destroy($1.items[i]);
     }
@@ -841,6 +843,7 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c,
     if ($1->medium.H_susceptibilities.items) {
         delete[] $1->medium.H_susceptibilities.items;
     }
+    delete[] $1->design_parameters;
     delete[] $1->epsilon_data;
     delete $1;
 }
@@ -1190,6 +1193,7 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c,
             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]->epsilon_data;
         }
         delete[] $1.items;
@@ -1404,6 +1408,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where
         GyrotropicSaturatedSusceptibility,
         Lattice,
         LorentzianSusceptibility,
+        MaterialGrid,
         Matrix,
         Medium,
         MultilevelAtom,

python/tests/material_grid.py

@@ -0,0 +1,75 @@
+import meep as mp
+import numpy as np
+import matplotlib.pyplot as plt
+from meep.materials import Si, SiO2
+
+resolution = 50 # pixels/μm
+
+cell_size = mp.Vector3(14,14)
+
+pml_layers = [mp.PML(thickness=2)]
+
+# rotation angle (in degrees) of waveguide, counter clockwise (CCW) around z-axis
+#rot_angle = np.radians(20)
+
+w = 3.0 # width of waveguide
+
+m1 = SiO2#mp.Medium(epsilon=2)
+m2 = Si#mp.Medium(epsilon=12)
+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 = []
+
+for r in np.linspace(0,360,4,endpoint=False):
+    rot_angle = np.radians(r)
+    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)]
+
+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,
+                    sources=sources,
+                    extra_materials = [m1,m2],
+                    #eps_averaging=False,
+                    geometry=geometry
+                    #symmetries=[mp.Mirror(mp.Y)] if rot_angle == 0 else []
+                    )
+
+sim.plot2D(frequency=1/1.55)
+plt.show()
+sim.reset_meep()
+mg.update_parameters(dp*0)
+sim.init_sim()
+sim.plot2D(frequency=1/1.55)
+plt.show()

python/typemap_utils.cpp

@@ -73,6 +73,15 @@ static PyObject *py_material_object() {
   return material_object;
 }
 
+static PyObject *py_material_grid_object() {
+  static PyObject *material_object = NULL;
+  if (material_object == NULL) {
+    PyObject *geom_mod = get_geom_mod();
+    material_object = PyObject_GetAttrString(geom_mod, "MaterialGrid");
+  }
+  return material_object;
+}
+
 static PyObject *py_vector3_object() {
   static PyObject *vector3_object = NULL;
   if (vector3_object == NULL) {
@@ -421,13 +430,69 @@ static int py_list_to_susceptibility_list(PyObject *po, susceptibility_list *sl)
   return 1;
 }
 
+static int pymaterial_grid_to_material_grid(PyObject *po, material_data *md) {
+  //po must be a python MaterialGrid object
+
+  // specify the type of material grid
+  long gt_enum = PyInt_AsLong(PyObject_GetAttrString(po, "grid_type"));
+  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 3: md->material_grid_kinds=material_data::U_DEFAULT; break;
+    default: meep::abort("Invalid material grid enumeration code: %d.\n",gt_enum);
+  }
+
+  // initialize grid size
+  if (!get_attr_v3(po, &md->grid_size, "grid_size")) { meep::abort("MaterialGrid grid_size failed to init.");}
+
+  // initialize user specified materials
+  PyObject *po_medium1 = PyObject_GetAttrString(po, "medium1");
+  if (!pymedium_to_medium(po_medium1, &md->medium_1)) { meep::abort("MaterialGrid medium1 failed to init."); }
+  PyObject *po_medium2 = PyObject_GetAttrString(po, "medium2");
+  if (!pymedium_to_medium(po_medium2, &md->medium_2)) { meep::abort("MaterialGrid medium2 failed to init."); }
+
+  // Initialize design parameters
+  PyObject *po_dp = PyObject_GetAttrString(po, "design_parameters");
+  PyArrayObject *pao = (PyArrayObject *)po_dp;
+  if (!PyArray_Check(pao)) { meep::abort("MaterialGrid design_parameters failed to init.");}
+  if (!PyArray_ISCARRAY(pao)) { meep::abort("Numpy array design_parameters 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));
+
+  // if needed, combine sus structs to main object
+  PyObject *py_e_sus_m1 = PyObject_GetAttrString(po_medium1, "E_susceptibilities");
+  PyObject *py_e_sus_m2 = PyObject_GetAttrString(po_medium2, "E_susceptibilities");
+  PyObject *py_h_sus_m1 = PyObject_GetAttrString(po_medium1, "H_susceptibilities");
+  PyObject *py_h_sus_m2 = PyObject_GetAttrString(po_medium2, "H_susceptibilities");
+  if (!py_e_sus_m1 || !py_e_sus_m2 || !py_h_sus_m1 || !py_h_sus_m2) { return 0;}
+
+  for (int i=0; i<PyList_Size(py_e_sus_m2);i++){
+    if (PyList_Append(py_e_sus_m1, PyList_GetItem(py_e_sus_m2, i)) != 0) {meep::abort("unable to merge e sus lists.\n");}
+  }
+  for (int i=0; i<PyList_Size(py_h_sus_m2);i++){
+    if (PyList_Append(py_h_sus_m1, PyList_GetItem(py_h_sus_m2, i)) != 0) {meep::abort("unable to merge h sus lists.\n");}
+  }
+
+  if (!py_list_to_susceptibility_list(py_e_sus_m1, &md->medium.E_susceptibilities) ||
+      !py_list_to_susceptibility_list(py_h_sus_m1, &md->medium.H_susceptibilities)) {
+    return 0;
+  }
+
+  return 1;
+}
+
 static int pymaterial_to_material(PyObject *po, material_type *mt) {
   material_data *md;
 
   if (PyObject_IsInstance(po, py_material_object())) {
     md = make_dielectric(1);
     if (!pymedium_to_medium(po, &md->medium)) { return 0; }
   }
+  else if (PyObject_IsInstance(po, py_material_grid_object())) { // Material grid subclass
+    md = make_material_grid();
+    if (!pymaterial_grid_to_material_grid(po, md)) { return 0; }
+  }
   else if (PyFunction_Check(po)) {
     PyObject *eps = PyObject_GetAttrString(po, "eps");
     PyObject *py_do_averaging = PyObject_GetAttrString(po, "do_averaging");
@@ -462,7 +527,7 @@ static int pymaterial_to_material(PyObject *po, material_type *mt) {
                   md->epsilon_dims[1], md->epsilon_dims[2]);
   }
   else {
-    meep::abort("Expected a Medium, a function, or a filename");
+    meep::abort("Expected a Medium, a Material Grid, a function, or a filename");
   }
 
   *mt = md;

src/fields.cpp

@@ -678,9 +678,9 @@ realnum linear_interpolate(realnum rx, realnum ry, realnum rz, realnum *data, in
     rz = 1.0 - rz;
 
   /* get the point corresponding to r in the epsilon array grid: */
-  x = pmod(int(rx * nx), nx);
-  y = pmod(int(ry * ny), ny);
-  z = pmod(int(rz * nz), nz);
+  x = rx == 1.0 ? nx-1 : pmod(int(rx * nx), nx);
+  y = ry == 1.0 ? ny-1 : pmod(int(ry * ny), ny);
+  z = rz == 1.0 ? nz-1 : pmod(int(rz * nz), nz);
 
   /* get the difference between (x,y,z) and the actual point
      ... we shift by 0.5 to center the data points in the pixels */

src/material_data.hpp

@@ -159,16 +159,21 @@ typedef void (*user_material_func)(vector3 x, void *user_data, medium_struct *me
 //                 'medium' field is filled in appropriately at
 //                 each evaluation point by calling the user's
 //                 routine.
+//  MATERIAL_GRID: material properties position-dependent, described
+//                 by user-supplied array of grid points. In this case 
+//                 the 'medium' field is filled in appropriately at
+//                 each evaluation point by interpolating the array.
 //  PERFECT_METAL: the 'medium' field is never referenced in this case.
 struct material_data {
   enum {
     MEDIUM,
     MATERIAL_FILE, // formerly MATERIAL_TYPE_SELF
     MATERIAL_USER, // formerly MATERIAL_FUNCTION
+    MATERIAL_GRID,
     PERFECT_METAL
   } which_subclass;
 
-  // this field is used for all material types except PERFECT_METAL
+  // this field is used for all material types except PERFECT_METAL and MATERIAL_GRID
   medium_struct medium;
 
   // these fields used only if which_subclass==MATERIAL_USER
@@ -180,11 +185,44 @@ struct material_data {
   meep::realnum *epsilon_data;
   size_t epsilon_dims[3];
 
+  // these fields used only if which_subclass==MATERIAL_GRID
+  vector3 grid_size;
+  meep::realnum* design_parameters;
+  medium_struct medium_1;
+  medium_struct medium_2;
+  /*
+  There are several possible scenarios when material grids overlap -- these
+  different scenarios enable different applications.
+
+  For U_MIN: Where multiple grids overlap, only those grids that contribute
+  the minimum u contribute, and for other grids the gradient is zero.
+  This unfortunately makes the gradient only piecewise continuous.
+
+  For U_PROD: The gradient is multiplied by the product of u's from
+  overlapping grids, divided by the u from the current grid.  This
+  unfortunately makes the gradient zero when two or more u's are zero,
+  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.
+  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".
+
+  For U_DEFAULT: Expect the default behavior with libctl objects; that is
+  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;
+
   material_data()
-      : which_subclass(MEDIUM), medium(), user_func(NULL), user_data(NULL), epsilon_data(NULL) {
+      : which_subclass(MEDIUM), medium(), user_func(NULL), user_data(NULL), epsilon_data(NULL), design_parameters(NULL), medium_1(), medium_2() {
     epsilon_dims[0] = 0;
     epsilon_dims[1] = 0;
     epsilon_dims[2] = 0;
+    grid_size.x = 0;
+    grid_size.y = 0;
+    grid_size.z = 0;
+    material_grid_kinds = U_DEFAULT;
   }
 };
 
@@ -204,7 +242,9 @@ extern material_type vacuum;
 material_type make_dielectric(double epsilon);
 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();
 void read_epsilon_file(const char *eps_input_file);
+void update_design_parameters(material_type matgrid, double* design_parameters);
 
 }; // namespace meep_geom