From 30c3edb198406eb9c817c205c408f25fa5373327 Mon Sep 17 00:00:00 2001
From: Chris Hogan <c096h800@gmail.com>
Date: Fri, 15 Mar 2019 11:17:09 -0500
Subject: [PATCH 1/2] Add option to use fallback integration for user material
functions
---
python/geom.py | 6 +++
python/simulation.py | 4 +-
python/solver.py | 3 +-
python/typemap_utils.cpp | 15 ++++---
src/material_data.hpp | 1 +
src/meepgeom.cpp | 75 +++++++++++++++++++--------------
src/meepgeom.hpp | 2 +-
tests/user-defined-material.cpp | 2 +-
8 files changed, 67 insertions(+), 41 deletions(-)
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..686cabc05 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;
@@ -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};
From 0844b7f26be4f2ae27e4d7e2de7b6e49c8063e7c Mon Sep 17 00:00:00 2001
From: Chris Hogan <c096h800@gmail.com>
Date: Thu, 21 Mar 2019 12:37:13 -0500
Subject: [PATCH 2/2] Handle divide by zero
---
src/meepgeom.cpp | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index 686cabc05..b7611c82c 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -913,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;
@@ -982,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; }