From 9254a1601e753b57a69bb57da5acdfc8d04ae198 Mon Sep 17 00:00:00 2001 From: Christopher Hogan Date: Thu, 7 Jun 2018 11:50:28 -0500 Subject: [PATCH] Add Python support for Prisms (#345) * Add python Prism wrapper * Increase tolerance and reuse compare_arrays from tests/mpb.py * Simplify vector assignment * Get new scheme data with fixed size blocks * Revert "Get new scheme data with fixed size blocks" This reverts commit 0203ab4aa1fc164e9f3b747a84e0b97280bbbd70. * height is double * Extend prism out of cell and increase height. --- python/geom.py | 9 ++++++++ python/meep.i | 3 ++- python/tests/bend_flux.py | 27 +++++++++++++++-------- python/tests/geom.py | 9 ++++++++ python/tests/mpb.py | 34 ++++++++++------------------- python/tests/utils.py | 15 +++++++++++++ python/typemap_utils.cpp | 46 +++++++++++++++++++++++++++++++++++++++ 7 files changed, 110 insertions(+), 33 deletions(-) create mode 100644 python/tests/utils.py diff --git a/python/geom.py b/python/geom.py index 5a207bc5a..01ff7ff06 100644 --- a/python/geom.py +++ b/python/geom.py @@ -341,6 +341,15 @@ def __init__(self, **kwargs): super(Ellipsoid, self).__init__(**kwargs) +class Prism(GeometricObject): + + def __init__(self, vertices, height, axis=Vector3(z=1), **kwargs): + self.vertices = vertices + self.height = height + self.axis = axis + super(Prism, self).__init__(**kwargs) + + class Matrix(object): def __init__(self, c1=Vector3(), c2=Vector3(), c3=Vector3()): diff --git a/python/meep.i b/python/meep.i index f0898f0e7..5ee7ff4ed 100644 --- a/python/meep.i +++ b/python/meep.i @@ -67,7 +67,7 @@ static int get_attr_int(PyObject *py_obj, int *result, const char *name) { } *result = PyInteger_AsLong(py_attr); - Py_DECREF(py_attr); + Py_XDECREF(py_attr); return 1; } @@ -1044,6 +1044,7 @@ void display_geometric_object_info(int indentby, GEOMETRIC_OBJECT o); Medium, NoisyDrudeSusceptibility, NoisyLorentzianSusceptibility, + Prism, Sphere, Susceptibility, Vector3, diff --git a/python/tests/bend_flux.py b/python/tests/bend_flux.py index e1c019409..f351edf46 100644 --- a/python/tests/bend_flux.py +++ b/python/tests/bend_flux.py @@ -3,6 +3,7 @@ import unittest import numpy as np import meep as mp +from utils import compare_arrays class TestBendFlux(unittest.TestCase): @@ -15,15 +16,24 @@ def init(self, no_bend=False): w = 1 wvg_ycen = -0.5 * (sy - w - (2 * pad)) wvg_xcen = 0.5 * (sx - w - (2 * pad)) + height = 100 if no_bend: - geometry = [mp.Block(mp.Vector3(mp.inf, w, mp.inf), center=mp.Vector3(0, wvg_ycen), - material=mp.Medium(epsilon=12))] + no_bend_vertices = [mp.Vector3(-0.5 * sx - 5, wvg_ycen - 0.5 * w), + mp.Vector3(+0.5 * sx + 5, wvg_ycen - 0.5 * w), + mp.Vector3(+0.5 * sx + 5, wvg_ycen + 0.5 * w), + mp.Vector3(-0.5 * sx - 5, wvg_ycen + 0.5 * w)] + + geometry = [mp.Prism(no_bend_vertices, height, material=mp.Medium(epsilon=12))] else: - geometry = [mp.Block(mp.Vector3(sx - pad, w, mp.inf), center=mp.Vector3(-0.5 * pad, wvg_ycen), - material=mp.Medium(epsilon=12)), - mp.Block(mp.Vector3(w, sy - pad, mp.inf), center=mp.Vector3(wvg_xcen, 0.5 * pad), - material=mp.Medium(epsilon=12))] + bend_vertices = [mp.Vector3(-0.5 * sx, wvg_ycen - 0.5 * w), + mp.Vector3(wvg_xcen + 0.5 * w, wvg_ycen - 0.5 * w), + mp.Vector3(wvg_xcen + 0.5 * w, 0.5 * sy), + mp.Vector3(wvg_xcen - 0.5 * w, 0.5 * sy), + mp.Vector3(wvg_xcen - 0.5 * w, wvg_ycen + 0.5 * w), + mp.Vector3(-0.5 * sx, wvg_ycen + 0.5 * w)] + + geometry = [mp.Prism(bend_vertices, height, material=mp.Medium(epsilon=12))] fcen = 0.15 df = 0.1 @@ -84,12 +94,11 @@ def test_bend_flux(self): (0.117171717172, 0.0147547920552, 0.0146151488236), (0.118181818182, 0.0194782085272, 0.0192840042241), (0.119191919192, 0.0254987474079, 0.0252348211592), - ] res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - np.testing.assert_allclose(expected, res[:20]) + compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-7) # Real run self.sim = None @@ -124,7 +133,7 @@ def test_bend_flux(self): res = list(zip(mp.get_flux_freqs(self.trans), mp.get_fluxes(self.trans), mp.get_fluxes(self.refl))) - np.testing.assert_allclose(expected, res[:20]) + compare_arrays(self, np.array(expected), np.array(res[:20]), tol=1e-3) if __name__ == '__main__': diff --git a/python/tests/geom.py b/python/tests/geom.py index e23623bcd..86d881839 100644 --- a/python/tests/geom.py +++ b/python/tests/geom.py @@ -245,6 +245,15 @@ def test_contains_point(self): self.assertIn(zeros(), e) +class TestPrism(unittest.TestCase): + + def test_contains_point(self): + vertices = [gm.Vector3(-1, 1), gm.Vector3(1, 1), gm.Vector3(1, -1), gm.Vector3(-1, -1)] + p = gm.Prism(vertices, height=1) + self.assertIn(zeros(), p) + self.assertNotIn(gm.Vector3(2, 2), p) + + class TestMedium(unittest.TestCase): def test_D_conductivity(self): diff --git a/python/tests/mpb.py b/python/tests/mpb.py index 368a36096..34a5541d7 100644 --- a/python/tests/mpb.py +++ b/python/tests/mpb.py @@ -14,6 +14,7 @@ from scipy.optimize import ridder import meep as mp from meep import mpb +from utils import compare_arrays class TestModeSolver(unittest.TestCase): @@ -178,20 +179,7 @@ def check_fields_against_h5(self, ref_path, field, suffix=''): ref_arr[1::3] = ref_y.ravel() ref_arr[2::3] = ref_z.ravel() - self.compare_arrays(ref_arr, field) - - def compare_arrays(self, exp, res, tol=1e-3): - exp_1d = exp.ravel() - res_1d = res.ravel() - - norm_exp = np.linalg.norm(exp_1d) - norm_res = np.linalg.norm(res_1d) - - if norm_exp == 0: - self.assertEqual(norm_res, 0) - else: - diff = np.linalg.norm(res_1d - exp_1d) / norm_exp - self.assertLess(diff, tol) + compare_arrays(self, ref_arr, field) def compare_h5_files(self, ref_path, res_path, tol=1e-3): with h5py.File(ref_path) as ref: @@ -200,7 +188,7 @@ def compare_h5_files(self, ref_path, res_path, tol=1e-3): if k == 'description': self.assertEqual(ref[k].value, res[k].value) else: - self.compare_arrays(ref[k].value, res[k].value, tol=tol) + compare_arrays(self, ref[k].value, res[k].value, tol=tol) def test_update_band_range_data(self): brd = [] @@ -407,7 +395,7 @@ def test_compute_field_energy(self): expected_energy_in_dielectric = 0.6990769686037558 - self.compare_arrays(np.array(expected_energy), np.array(energy)) + compare_arrays(self, np.array(expected_energy), np.array(energy)) self.assertAlmostEqual(expected_energy_in_dielectric, energy_in_dielectric, places=3) def test_compute_group_velocity(self): @@ -693,7 +681,7 @@ def get_dpwr(ms, band): with h5py.File(ref_path, 'r') as f: expected = f['data-new'].value - self.compare_arrays(expected, converted_dpwr[-1]) + compare_arrays(self, expected, converted_dpwr[-1]) def test_hole_slab(self): from mpb_hole_slab import ms @@ -958,7 +946,7 @@ def test_tri_rods(self): expected_re = f['z.r-new'].value expected_im = f['z.i-new'].value expected = np.vectorize(complex)(expected_re, expected_im) - self.compare_arrays(expected, new_efield) + compare_arrays(self, expected, new_efield) ms.run_te() @@ -992,7 +980,7 @@ def test_tri_rods(self): with h5py.File(ref_path, 'r') as f: ref = f['data-new'].value - self.compare_arrays(ref, new_eps, tol=1e-3) + compare_arrays(self, ref, new_eps, tol=1e-3) def test_subpixel_averaging(self): ms = self.init_solver() @@ -1065,7 +1053,7 @@ def test_run_te_with_mu_material(self): mu = ms.get_mu() with h5py.File(data_path, 'r') as f: data = f['data'].value - self.compare_arrays(data, mu) + compare_arrays(self, data, mu) def test_output_tot_pwr(self): ms = self.init_solver() @@ -1083,7 +1071,7 @@ def test_output_tot_pwr(self): with h5py.File(ref_path, 'r') as f: expected = f['data'].value - self.compare_arrays(expected, arr) + compare_arrays(self, expected, arr) def test_get_eigenvectors(self): ms = self.init_solver() @@ -1196,7 +1184,7 @@ def test_epsilon_input_file(self): ] self.check_band_range_data(expected_brd, ms.band_range_data) - self.compare_arrays(expected_freqs, ms.all_freqs[-1]) + compare_arrays(self, expected_freqs, ms.all_freqs[-1]) self.check_gap_list(expected_gap_list, ms.gap_list) @@ -1340,7 +1328,7 @@ def test_multiply_bloch_in_map_data(self): md = mpb.MPBData(rectify=True, resolution=32, periods=3) result2 = md.convert(efield_no_bloch) - self.compare_arrays(result1, result2, tol=1e-7) + compare_arrays(self, result1, result2, tol=1e-7) def test_poynting(self): ms = self.init_solver() diff --git a/python/tests/utils.py b/python/tests/utils.py new file mode 100644 index 000000000..2e011d1f7 --- /dev/null +++ b/python/tests/utils.py @@ -0,0 +1,15 @@ +import numpy as np + + +def compare_arrays(test_instance, exp, res, tol=1e-3): + exp_1d = exp.ravel() + res_1d = res.ravel() + + norm_exp = np.linalg.norm(exp_1d) + norm_res = np.linalg.norm(res_1d) + + if norm_exp == 0: + test_instance.assertEqual(norm_res, 0) + else: + diff = np.linalg.norm(res_1d - exp_1d) / norm_exp + test_instance.assertLess(diff, tol) diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp index 646b098a0..69f1e94ef 100644 --- a/python/typemap_utils.cpp +++ b/python/typemap_utils.cpp @@ -541,6 +541,49 @@ static int pyellipsoid_to_ellipsoid(PyObject *py_ell, geometric_object *e) { return 1; } +static int pyprism_to_prism(PyObject *py_prism, geometric_object *p) { + material_type material; + double height; + vector3 axis; + + if (!get_attr_material(py_prism, &material) || + !get_attr_dbl(py_prism, &height, "height") || + !get_attr_v3(py_prism, &axis, "axis")) { + + return 0; + } + + PyObject *py_vert_list = PyObject_GetAttrString(py_prism, "vertices"); + + if (!py_vert_list) { + PyErr_PrintEx(0); + return 0; + } + + if (!PyList_Check(py_vert_list)) { + PyErr_SetString(PyExc_TypeError, "Expected Prism.vertices to be a list\n"); + return 0; + } + + int num_vertices = PyList_Size(py_vert_list); + vector3 *vertices = new vector3[num_vertices]; + + for (Py_ssize_t i = 0; i < num_vertices; ++i) { + vector3 v3; + if (!pyv3_to_v3(PyList_GetItem(py_vert_list, i), &v3)) { + return 0; + } + vertices[i] = v3; + } + + *p = make_prism(material, vertices, num_vertices, height, axis); + + delete [] vertices; + Py_DECREF(py_vert_list); + + return 1; +} + static int py_gobj_to_gobj(PyObject *po, geometric_object *o) { int success = 0; std::string go_type = py_class_name_as_string(po); @@ -563,6 +606,9 @@ static int py_gobj_to_gobj(PyObject *po, geometric_object *o) { else if (go_type == "Ellipsoid") { success = pyellipsoid_to_ellipsoid(po, o); } + else if (go_type == "Prism") { + success = pyprism_to_prism(po, o); + } else { PyErr_Format(PyExc_TypeError, "Error: %s is not a valid GeometricObject type\n", go_type.c_str()); return 0;