enable support for subpixel smoothing of MaterialGrid (#1503)

* enable support for subpixel smoothing of MaterialGrid

* apply hyperbolic tangent projection/thresholding after bilinear interpolation

* fixes

* minor fix to if statement

python/geom.py

@@ -523,7 +523,7 @@ def _get_epsmu(self, diag, offdiag, susceptibilities, conductivity_diag, conduct
         return np.squeeze(epsmu)
 
 class MaterialGrid(object):
-    def __init__(self,grid_size,medium1,medium2,design_parameters=None,grid_type="U_DEFAULT"):
+    def __init__(self,grid_size,medium1,medium2,design_parameters=None,grid_type="U_DEFAULT",do_averaging=False,beta=0,eta=0.5):
         self.grid_size = mp.Vector3(*grid_size)
         self.medium1 = medium1
         self.medium2 = medium2
@@ -536,7 +536,9 @@ def isclose(a, b, rel_tol=1e-09, abs_tol=0.0):
         if isclose(self.grid_size.z,0):
             self.grid_size.z = 1
         self.num_params=int(self.grid_size.x*self.grid_size.y*self.grid_size.z)
-
+        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:

python/typemap_utils.cpp

@@ -522,8 +522,20 @@ static int pymaterial_to_material(PyObject *po, material_type *mt) {
     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();
+    PyObject *py_do_averaging = PyObject_GetAttrString(po, "do_averaging");
+    bool do_averaging = false;
+    if (py_do_averaging) { do_averaging = PyObject_IsTrue(py_do_averaging); }
+    PyObject *py_beta = PyObject_GetAttrString(po, "beta");
+    double beta = 0;
+    if (py_beta) { beta = PyFloat_AsDouble(py_beta); }
+    PyObject *py_eta = PyObject_GetAttrString(po, "eta");
+    double eta = 0;
+    if (py_eta) { eta = PyFloat_AsDouble(py_eta); }
+    md = make_material_grid(do_averaging,beta,eta);
     if (!pymaterial_grid_to_material_grid(po, md)) { return 0; }
+    Py_XDECREF(py_do_averaging);
+    Py_XDECREF(py_beta);
+    Py_XDECREF(py_eta);
   }
   else if (PyFunction_Check(po)) {
     PyObject *eps = PyObject_GetAttrString(po, "eps");

src/material_data.hpp

@@ -190,6 +190,8 @@ struct material_data {
   meep::realnum *design_parameters;
   medium_struct medium_1;
   medium_struct medium_2;
+  meep::realnum beta;
+  meep::realnum eta;
   /*
   There are several possible scenarios when material grids overlap -- these
   different scenarios enable different applications.
@@ -243,7 +245,7 @@ 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();
+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);
 

src/meepgeom.cpp

@@ -350,12 +350,22 @@ meep::realnum matgrid_val(vector3 p, geom_box_tree tp, int oi, material_data *md
     ++matgrid_val_count;
   }
 
-  return (md->material_grid_kinds == material_data::U_MIN
-              ? umin
-              : (md->material_grid_kinds == material_data::U_PROD
-                     ? uprod
-                     : (md->material_grid_kinds == material_data::U_SUM ? usum / matgrid_val_count
-                                                                        : udefault)));
+  meep::realnum u_interp = (md->material_grid_kinds == material_data::U_MIN
+                            ? umin
+                            : (md->material_grid_kinds == material_data::U_PROD
+                               ? uprod
+                               : (md->material_grid_kinds == material_data::U_SUM ? usum / matgrid_val_count
+                                  : udefault)));
+
+  // project interpolated grid point
+  meep::realnum u_proj;
+  if (md->beta != 0) {
+    double tanh_beta_eta = tanh(md->beta*md->eta);
+    u_proj = (tanh_beta_eta + tanh(md->beta*(u_interp-md->eta))) /
+      (tanh_beta_eta + tanh(md->beta*(1-md->eta)));
+  }
+
+  return (md->beta != 0) ? u_proj : u_interp;
 }
 static void cinterp_tensors(vector3 diag_in_1, cvector3 offdiag_in_1, vector3 diag_in_2,
                             cvector3 offdiag_in_2, vector3 *diag_out, cvector3 *offdiag_out,
@@ -697,8 +707,10 @@ void geom_epsilon::get_material_pt(material_type &material, const meep::vec &r)
       geom_box_tree tp;
 
       tp = geom_tree_search(p, restricted_tree, &oi);
-      u = matgrid_val(p, tp, oi, md); // interpolate onto material grid
-      epsilon_material_grid(md, u);   // interpolate material from material grid point
+      // interpolate and (possibly) project onto material grid
+      u = matgrid_val(p, tp, oi, md);
+      // interpolate material from material grid point
+      epsilon_material_grid(md, u);
 
       return;
     // material read from file: interpolate to get properties at r
@@ -869,7 +881,6 @@ void geom_epsilon::eff_chi1inv_row(meep::component c, double chi1inv_row[3], con
   symmetric_matrix meps_inv;
   bool fallback;
   eff_chi1inv_matrix(c, &meps_inv, v, tol, maxeval, fallback);
-  ;
 
   if (fallback) { fallback_chi1inv_row(c, chi1inv_row, v, tol, maxeval); }
   else {
@@ -916,7 +927,9 @@ void geom_epsilon::eff_chi1inv_matrix(meep::component c, symmetric_matrix *chi1i
 
   if (!get_front_object(v, geometry_tree, p, &o, shiftby, mat, mat_behind)) {
     get_material_pt(mat, v.center());
-    if (mat && mat->which_subclass == material_data::MATERIAL_USER && mat->do_averaging) {
+    if (mat && (mat->which_subclass == material_data::MATERIAL_USER ||
+                mat->which_subclass == material_data::MATERIAL_GRID)
+        && mat->do_averaging) {
       fallback = true;
       return;
     }
@@ -1797,10 +1810,12 @@ material_type make_file_material(const char *eps_input_file) {
 /******************************************************************************/
 /* Material grid functions                                                    */
 /******************************************************************************/
-material_type make_material_grid() {
+material_type make_material_grid(bool do_averaging, double beta, double eta) {
   material_data *md = new material_data();
   md->which_subclass = material_data::MATERIAL_GRID;
-  md->do_averaging = false;
+  md->do_averaging = do_averaging;
+  md->beta = beta;
+  md->eta = eta;
   return md;
 }
 

src/meepgeom.hpp

@@ -174,7 +174,7 @@ void set_materials_from_geometry(meep::structure *s, geometric_object_list g,
 material_type make_dielectric(double epsilon);
 material_type make_user_material(user_material_func user_func, void *user_data, bool do_averaging);
 material_type make_file_material(const char *eps_input_file);
-material_type make_material_grid();
+material_type make_material_grid(bool do_averaging, double beta, double eta);
 
 vector3 vec_to_vector3(const meep::vec &pt);
 meep::vec vector3_to_vec(const vector3 v3);