diff --git a/python/geom.py b/python/geom.py index ead5b6389..a9b535c90 100644 --- a/python/geom.py +++ b/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 diff --git a/python/simulation.py b/python/simulation.py index f7b947104..83eb13053 100644 --- a/python/simulation.py +++ b/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: diff --git a/python/solver.py b/python/solver.py index 0fcf760c9..481796922 100644 --- a/python/solver.py +++ b/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( diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp index 94cf28f2b..af093874d 100644 --- a/python/typemap_utils.cpp +++ b/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); diff --git a/src/material_data.hpp b/src/material_data.hpp index bd9cce0a1..0610d7fdb 100644 --- a/src/material_data.hpp +++ b/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; diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index cdc90356e..b7611c82c 100644 --- a/src/meepgeom.cpp +++ b/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; } diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index 067f671e5..449f8dd49 100644 --- a/src/meepgeom.hpp +++ b/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); diff --git a/tests/user-defined-material.cpp b/tests/user-defined-material.cpp index 2ea7d444c..154e7be82 100644 --- a/tests/user-defined-material.cpp +++ b/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};