Add option to use fallback integration for user material functions (#771

)

* Add option to use fallback integration for user material functions

* Handle divide by zero

python/geom.py

@@ -22,6 +22,10 @@ def check_nonnegative(prop, val):
     else:
         raise ValueError("{} cannot be negative. Got {}".format(prop, val))
 
+def init_do_averaging(mat_func):
+    if not hasattr(mat_func, 'do_averaging'):
+        mat_func.do_averaging = False
+
 
 class Vector3(object):
 
@@ -303,8 +307,10 @@ class GeometricObject(object):
 
     def __init__(self, material=Medium(), center=Vector3(), epsilon_func=None):
         if type(material) is not Medium and callable(material):
+            init_do_averaging(material)
             material.eps = False
         elif epsilon_func:
+            init_do_averaging(epsilon_func)
             epsilon_func.eps = True
             material = epsilon_func
 

python/simulation.py

@@ -15,7 +15,7 @@
 import numpy as np
 
 import meep as mp
-from meep.geom import Vector3
+from meep.geom import Vector3, init_do_averaging
 from meep.source import EigenModeSource, check_positive
 
 
@@ -945,9 +945,11 @@ def _init_structure(self, k=False):
         absorbers = [bl for bl in self.boundary_layers if type(bl) is Absorber]
 
         if self.material_function:
+            init_do_averaging(self.material_function)
             self.material_function.eps = False
             self.default_material = self.material_function
         elif self.epsilon_func:
+            init_do_averaging(self.epsilon_func)
             self.epsilon_func.eps = True
             self.default_material = self.epsilon_func
         elif self.epsilon_input_file:

python/solver.py

@@ -12,8 +12,8 @@
 import numpy as np
 import meep as mp
 from . import mode_solver, with_hermitian_epsilon
+from meep.geom import init_do_averaging
 from meep.simulation import get_num_args
-
 try:
     basestring
 except NameError:
@@ -123,6 +123,7 @@ def __init__(self,
         grid_size = self._adjust_grid_size()
 
         if type(self.default_material) is not mp.Medium and callable(self.default_material):
+            init_do_averaging(self.default_material)
             self.default_material.eps = False
 
         self.mode_solver = mode_solver(

python/typemap_utils.cpp

@@ -405,13 +405,18 @@ static int pymaterial_to_material(PyObject *po, material_type *mt) {
     if (!pymedium_to_medium(po, &md->medium)) { return 0; }
   } else if (PyFunction_Check(po)) {
     PyObject *eps = PyObject_GetAttrString(po, "eps");
-    if (!eps) { return 0; }
-    if (eps == Py_True) {
-      md = make_user_material(py_epsilon_func_wrap, po);
+    PyObject *py_do_averaging = PyObject_GetAttrString(po, "do_averaging");
+    bool do_averaging = false;
+    if (py_do_averaging) {
+      do_averaging = PyObject_IsTrue(py_do_averaging);
+    }
+    if (eps && eps == Py_True) {
+      md = make_user_material(py_epsilon_func_wrap, po, do_averaging);
     } else {
-      md = make_user_material(py_user_material_func_wrap, po);
+      md = make_user_material(py_user_material_func_wrap, po, do_averaging);
     }
-    Py_DECREF(eps);
+    Py_XDECREF(eps);
+    Py_XDECREF(py_do_averaging);
   } else if (IsPyString(po)) {
     const char *eps_input_file = PyObject_ToCharPtr(po);
     md = make_file_material(eps_input_file);

src/material_data.hpp

@@ -171,6 +171,7 @@ struct material_data {
   // these fields used only if which_subclass==MATERIAL_USER
   user_material_func user_func;
   void *user_data;
+  bool do_averaging;
 
   // these fields used only if which_subclass==MATERIAL_FILE
   meep::realnum *epsilon_data;

src/meepgeom.cpp

@@ -346,7 +346,7 @@ class geom_epsilon : public meep::material_function {
                                double tol, int maxeval);
 
   void eff_chi1inv_matrix(meep::component c, symmetric_matrix *chi1inv_matrix,
-                          const meep::volume &v, double tol, int maxeval);
+                          const meep::volume &v, double tol, int maxeval, bool &fallback);
 
   void fallback_chi1inv_row(meep::component c, double chi1inv_row[3], const meep::volume &v,
                             double tol, int maxeval);
@@ -685,38 +685,46 @@ static bool get_front_object(const meep::volume &v, geom_box_tree geometry_tree,
 void geom_epsilon::eff_chi1inv_row(meep::component c, double chi1inv_row[3], const meep::volume &v,
                                    double tol, int maxeval) {
   symmetric_matrix meps_inv;
-  eff_chi1inv_matrix(c, &meps_inv, v, tol, maxeval);
-
-  switch (component_direction(c)) {
-    case meep::X:
-    case meep::R:
-      chi1inv_row[0] = meps_inv.m00;
-      chi1inv_row[1] = meps_inv.m01;
-      chi1inv_row[2] = meps_inv.m02;
-      break;
-    case meep::Y:
-    case meep::P:
-      chi1inv_row[0] = meps_inv.m01;
-      chi1inv_row[1] = meps_inv.m11;
-      chi1inv_row[2] = meps_inv.m12;
-      break;
-    case meep::Z:
-      chi1inv_row[0] = meps_inv.m02;
-      chi1inv_row[1] = meps_inv.m12;
-      chi1inv_row[2] = meps_inv.m22;
-      break;
-    case meep::NO_DIRECTION: chi1inv_row[0] = chi1inv_row[1] = chi1inv_row[2] = 0; break;
+  bool fallback;
+  eff_chi1inv_matrix(c, &meps_inv, v, tol, maxeval, fallback);;
+
+  if (fallback) {
+    fallback_chi1inv_row(c, chi1inv_row, v, tol, maxeval);
+  }
+  else {
+    switch (component_direction(c)) {
+      case meep::X:
+      case meep::R:
+        chi1inv_row[0] = meps_inv.m00;
+        chi1inv_row[1] = meps_inv.m01;
+        chi1inv_row[2] = meps_inv.m02;
+        break;
+      case meep::Y:
+      case meep::P:
+        chi1inv_row[0] = meps_inv.m01;
+        chi1inv_row[1] = meps_inv.m11;
+        chi1inv_row[2] = meps_inv.m12;
+        break;
+      case meep::Z:
+        chi1inv_row[0] = meps_inv.m02;
+        chi1inv_row[1] = meps_inv.m12;
+        chi1inv_row[2] = meps_inv.m22;
+        break;
+      case meep::NO_DIRECTION: chi1inv_row[0] = chi1inv_row[1] = chi1inv_row[2] = 0; break;
+    }
   }
 }
 
 void geom_epsilon::eff_chi1inv_matrix(meep::component c, symmetric_matrix *chi1inv_matrix,
-                                      const meep::volume &v, double tol, int maxeval) {
+                                      const meep::volume &v, double tol, int maxeval,
+                                      bool &fallback) {
   const geometric_object *o;
   material_type mat, mat_behind;
   symmetric_matrix meps;
   vector3 p, shiftby, normal;
+  fallback = false;
 
-  if (maxeval == 0 || !get_front_object(v, geometry_tree, p, &o, shiftby, mat, mat_behind)) {
+  if (maxeval == 0) {
   noavg:
     get_material_pt(mat, v.center());
   trivial:
@@ -725,13 +733,15 @@ void geom_epsilon::eff_chi1inv_matrix(meep::component c, symmetric_matrix *chi1i
     return;
   }
 
-  // FIXME: reimplement support for fallback integration, without
-  //        messing up anisotropic support
-  //  if (!get_front_object(v, geometry_tree,
-  //                        p, &o, shiftby, mat, mat_behind)) {
-  //     fallback_chi1inv_row(c, chi1inv_row, v, tol, maxeval);
-  //     return;
-  //  }
+  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) {
+      fallback = true;
+      return;
+    } else {
+      goto trivial;
+    }
+  }
 
   /* check for trivial case of only one object/material */
   if (material_type_equal(mat, mat_behind)) goto trivial;
@@ -903,11 +913,12 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c, double chi1inv_row[3]
 
   symmetric_matrix chi1p1, chi1p1_inv;
   material_type material;
+  meep::vec gradient(normal_vector(meep::type(c), v));
   get_material_pt(material, v.center());
   material_epsmu(meep::type(c), material, &chi1p1, &chi1p1_inv);
   material_gc(material);
   if (chi1p1.m01 != 0 || chi1p1.m02 != 0 || chi1p1.m12 != 0 || chi1p1.m00 != chi1p1.m11 ||
-      chi1p1.m11 != chi1p1.m22 || chi1p1.m00 != chi1p1.m22) {
+      chi1p1.m11 != chi1p1.m22 || chi1p1.m00 != chi1p1.m22 || meep::abs(gradient) == 0) {
     int rownum = meep::component_direction(c) % 3;
     if (rownum == 0) {
       chi1inv_row[0] = chi1p1.m00;
@@ -972,7 +983,6 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c, double chi1inv_row[3]
     minveps = 1.0 / (meps = eps(v.center()));
 
   {
-    meep::vec gradient(normal_vector(meep::type(c), v));
     double n[3] = {0, 0, 0};
     double nabsinv = 1.0 / meep::abs(gradient);
     LOOP_OVER_DIRECTIONS(gradient.dim, k) { n[k % 3] = gradient.in_direction(k) * nabsinv; }
@@ -1525,11 +1535,12 @@ material_type make_dielectric(double epsilon) {
   return md;
 }
 
-material_type make_user_material(user_material_func user_func, void *user_data) {
+material_type make_user_material(user_material_func user_func, void *user_data, bool do_averaging) {
   material_data *md = new material_data();
   md->which_subclass = material_data::MATERIAL_USER;
   md->user_func = user_func;
   md->user_data = user_data;
+  md->do_averaging = do_averaging;
   return md;
 }
 

src/meepgeom.hpp

@@ -166,7 +166,7 @@ void set_materials_from_geometry(meep::structure *s, geometric_object_list g,
                                  material_type_list extra_materials = material_type_list());
 
 material_type make_dielectric(double epsilon);
-material_type make_user_material(user_material_func user_func, void *user_data);
+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);
 
 vector3 vec_to_vector3(const meep::vec &pt);

tests/user-defined-material.cpp

@@ -189,7 +189,7 @@ int main(int argc, char *argv[]) {
   data.rOuter = R2;
 
   meep_geom::material_type my_material =
-      meep_geom::make_user_material(my_material_func, (void *)&data);
+      meep_geom::make_user_material(my_material_func, (void *)&data, false);
 
   geometric_object_list g = {0, 0};
   vector3 center = {0, 0, 0};