From 06bb2fda8269a6e9ddff817c34b4920f82dfdd46 Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Sun, 15 Sep 2019 14:54:57 -0500 Subject: [PATCH 1/5] WIP: geom_box_tree for faster fragment_stats --- python/meep.i | 3 +-- src/meepgeom.cpp | 14 ++++++++++++++ src/meepgeom.hpp | 2 ++ 3 files changed, 17 insertions(+), 2 deletions(-) diff --git a/python/meep.i b/python/meep.i index 081b6dfe8..4596202cd 100644 --- a/python/meep.i +++ b/python/meep.i @@ -1660,8 +1660,7 @@ meep::structure *create_structure_and_set_materials(vector3 cell_size, } // Return params to default state - meep_geom::fragment_stats::resolution = 0; - meep_geom::fragment_stats::split_chunks_evenly = false; + meep_geom::fragment_stats::reset(); return s; } diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index ca6699d56..0d6fa85a2 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -1639,6 +1639,7 @@ int fragment_stats::maxeval = 0; int fragment_stats::resolution = 0; meep::ndim fragment_stats::dims = meep::D1; geometric_object_list fragment_stats::geom = {}; +geom_box_tree fragment_stats::geom_tree = NULL; std::vector fragment_stats::dft_data_list; std::vector fragment_stats::pml_1d_vols; std::vector fragment_stats::pml_2d_vols; @@ -1706,6 +1707,7 @@ fragment_stats compute_fragment_stats( geom_box box = make_box_from_cell(cell_size); fragment_stats stats = init_stats(box, tol, maxeval, gv); stats.compute(); + fragment_stats::reset(); return stats; } @@ -1717,6 +1719,13 @@ fragment_stats::fragment_stats(geom_box &bx) num_pixels_in_box = get_pixels_in_box(&bx); } +void fragment_stats::reset() { + resolution = 0; + split_chunks_evenly = false; + destroy_geom_box_tree(geom_tree); + geom_tree = NULL; +} + void fragment_stats::init_libctl(material_type default_mat, bool ensure_per, meep::grid_volume *gv, vector3 cell_size, vector3 cell_center, geometric_object_list *geom_) { @@ -1727,6 +1736,8 @@ void fragment_stats::init_libctl(material_type default_mat, bool ensure_per, mee dimensions = meep::number_of_directions(gv->dim); geometry_lattice.size = cell_size; geom_fix_object_list(*geom_); + geom_box cell_box = make_box_from_cell(cell_size); + geom_tree = create_geom_box_tree0(*geom_, cell_box); } bool fragment_stats::has_non_medium_material() { @@ -1764,6 +1775,9 @@ void fragment_stats::compute_stats() { update_stats_from_material((material_type)default_material, num_pixels_in_box); } + geom_box_tree current_tree = geom_tree; + std::vector potential_objs_in_box; + for (int i = 0; i < geom.num_items; ++i) { geometric_object *go = &geom.items[i]; double overlap = box_overlap_with_object(box, *go, tol, maxeval); diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index 46d4c02a4..9c7b1b8aa 100644 --- a/src/meepgeom.hpp +++ b/src/meepgeom.hpp @@ -67,6 +67,7 @@ struct fragment_stats { static int resolution; static meep::ndim dims; static geometric_object_list geom; + static geom_box_tree geom_tree; static std::vector dft_data_list; static std::vector pml_1d_vols; static std::vector pml_2d_vols; @@ -79,6 +80,7 @@ struct fragment_stats { static void init_libctl(meep_geom::material_type default_mat, bool ensure_per, meep::grid_volume *gv, vector3 cell_size, vector3 cell_center, geometric_object_list *geom_list); + static void reset(); size_t num_anisotropic_eps_pixels; size_t num_anisotropic_mu_pixels; From 10420b5196bab1f0dc1ff85ecf33e7dca876cc0d Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Sun, 22 Sep 2019 15:33:37 -0500 Subject: [PATCH 2/5] WIP: speed up choose_chunkdivision with geom_box_tree --- src/meepgeom.cpp | 101 +++++++++++++++++++++++++++++++---------------- src/meepgeom.hpp | 2 + 2 files changed, 68 insertions(+), 35 deletions(-) diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index 0d6fa85a2..0f4968fa9 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -1731,13 +1731,15 @@ void fragment_stats::init_libctl(material_type default_mat, bool ensure_per, mee geometric_object_list *geom_) { geom_initialize(); default_material = default_mat; - ensure_periodicity = ensure_per; geometry_center = cell_center; dimensions = meep::number_of_directions(gv->dim); geometry_lattice.size = cell_size; geom_fix_object_list(*geom_); geom_box cell_box = make_box_from_cell(cell_size); + // Don't let libctl create duplicate objects in the tree + ensure_periodicity = 0; geom_tree = create_geom_box_tree0(*geom_, cell_box); + // ensure_periodicity = ensure_per; } bool fragment_stats::has_non_medium_material() { @@ -1768,50 +1770,79 @@ void fragment_stats::update_stats_from_material(material_type mat, size_t pixels } } +void fragment_stats::tree_search(geom_box_tree t, std::vector &overlaps) { + if (!t || !geom_boxes_intersect(&t->b, &box)) + return; + + if (t->nobjects > 0) { + geom_box intersection; + geom_box_intersection(&intersection, &t->b, &box); + + for (int i = 0; i < t->nobjects; ++i) { + if (geom_boxes_intersect(&intersection, &t->objects[i].box)) { + geom_box shifted_box = {vector3_minus(intersection.low, t->objects[i].shiftby), + vector3_minus(intersection.high, t->objects[i].shiftby)}; + double overlap = box_overlap_with_object(shifted_box, *t->objects[i].o, tol, maxeval); + if (overlap > 0.0) { + for (int j = 0; j < geom.num_items; ++j) { + if (t->objects[i].o == geom.items + j) { + overlaps[j] += overlap; + } + } + } + } + } + } + + tree_search(t->t1, overlaps); + tree_search(t->t2, overlaps); +} + void fragment_stats::compute_stats() { if (geom.num_items == 0) { // If there is no geometry, count the default material for the whole fragment update_stats_from_material((material_type)default_material, num_pixels_in_box); } - - geom_box_tree current_tree = geom_tree; - std::vector potential_objs_in_box; - - for (int i = 0; i < geom.num_items; ++i) { - geometric_object *go = &geom.items[i]; - double overlap = box_overlap_with_object(box, *go, tol, maxeval); - - bool anisotropic_pixels_already_added = false; - if (eps_averaging) { - // If the object doesn't overlap the entire box, that implies that - // an object interface intercepts the box, which means we treat - // the entire box as anisotropic. This method could give some false - // positives if there is another object with the same material behind - // the current object, but in practice it is probably reasonable to - // assume that there is a material interface somwhere in the box so - // we won't worry about fancier edge-detection methods for now. - if (overlap != 1.0) { - anisotropic_pixels_already_added = true; - num_anisotropic_eps_pixels += num_pixels_in_box; - if (mu_not_1(go->material)) { - num_anisotropic_mu_pixels += num_pixels_in_box; + else { + std::vector overlaps(geom.num_items, 0.0); + tree_search(geom_tree, overlaps); + + for (int i = 0; i < geom.num_items; ++i) { + geometric_object *go = &geom.items[i]; + double overlap = overlaps[i]; + + bool anisotropic_pixels_already_added = false; + if (eps_averaging) { + // If the object doesn't overlap the entire box, that implies that + // an object interface intercepts the box, which means we treat + // the entire box as anisotropic. This method could give some false + // positives if there is another object with the same material behind + // the current object, but in practice it is probably reasonable to + // assume that there is a material interface somwhere in the box so + // we won't worry about fancier edge-detection methods for now. + if (overlap != 1.0) { + anisotropic_pixels_already_added = true; + num_anisotropic_eps_pixels += num_pixels_in_box; + if (mu_not_1(go->material)) { + num_anisotropic_mu_pixels += num_pixels_in_box; + } } } - } - // Count contributions from material of object - size_t pixels = (size_t)ceil(overlap * num_pixels_in_box); - if (pixels > 0) { - material_type mat = (material_type)go->material; - update_stats_from_material(mat, pixels, anisotropic_pixels_already_added); - } + // Count contributions from material of object + size_t pixels = (size_t)ceil(overlap * num_pixels_in_box); + if (pixels > 0) { + material_type mat = (material_type)go->material; + update_stats_from_material(mat, pixels, anisotropic_pixels_already_added); + } - // Count contributions from default_material - size_t default_material_pixels = num_pixels_in_box - pixels; - if (default_material_pixels > 0) { - update_stats_from_material((material_type)default_material, default_material_pixels, - anisotropic_pixels_already_added); + // Count contributions from default_material + size_t default_material_pixels = num_pixels_in_box - pixels; + if (default_material_pixels > 0) { + update_stats_from_material((material_type)default_material, default_material_pixels, + anisotropic_pixels_already_added); + } } } } diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index 9c7b1b8aa..cb0c84ab6 100644 --- a/src/meepgeom.hpp +++ b/src/meepgeom.hpp @@ -118,6 +118,8 @@ struct fragment_stats { void compute_dft_stats(); void compute_pml_stats(); void compute_absorber_stats(); + // TODO rename + void tree_search(geom_box_tree t, std::vector &overlaps); }; fragment_stats From db9fe3584d3bbab11762a5a785846f5312ee7beb Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Sun, 29 Sep 2019 11:59:18 -0500 Subject: [PATCH 3/5] WIP: Speed up fragment_stats with geom_box_tree --- python/tests/fragment_stats.py | 14 ++++++++++++ src/meepgeom.cpp | 41 +++++++++++++++++++++------------- src/meepgeom.hpp | 3 +-- 3 files changed, 41 insertions(+), 17 deletions(-) diff --git a/python/tests/fragment_stats.py b/python/tests/fragment_stats.py index 029f34a91..6d569a8d9 100644 --- a/python/tests/fragment_stats.py +++ b/python/tests/fragment_stats.py @@ -307,6 +307,20 @@ def test_3d_cell_smaller_than_minimum_fragment_size(self): self.assertEqual(fs.box.high.z, 0.5) self.assertEqual(fs.num_pixels_in_box, 1000) + def test_geom_box_tree(self): + + block_size = mp.Vector3(3, 3) + mat = mp.Medium(index=13) + geometry = [ + mp.Block(size=block_size, center=mp.Vector3(-3, 3), material=mat), + mp.Block(size=block_size, center=mp.Vector3(3, 3), material=mat), + mp.Block(size=block_size, center=mp.Vector3(3, -3), material=mat), + mp.Block(size=block_size, center=mp.Vector3(-3, -3), material=mat) + ] + fs = self.get_fragment_stats(mp.Vector3(), mp.Vector3(9, 9), 2, geom=geometry) + self.check_stats(fs, 32400, 0, 0, 0, 0) + self.assertEqual(fs.num_pixels_in_box, 8100) + class TestPMLToVolList(unittest.TestCase): diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index 0f4968fa9..3081859b2 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -1736,10 +1736,17 @@ void fragment_stats::init_libctl(material_type default_mat, bool ensure_per, mee geometry_lattice.size = cell_size; geom_fix_object_list(*geom_); geom_box cell_box = make_box_from_cell(cell_size); - // Don't let libctl create duplicate objects in the tree + int temp_dims = dimensions; + int temp_periodicity = ensure_periodicity; + dimensions = 3; ensure_periodicity = 0; geom_tree = create_geom_box_tree0(*geom_, cell_box); - // ensure_periodicity = ensure_per; + // if (meep::am_master()) display_geom_box_tree(2, geom_tree); + // int depth, nobjects; + // geom_box_tree_stats(geom_tree, &depth, &nobjects); + // master_printf("depth: %d, nobjects: %d\n", depth, nobjects); + dimensions = temp_dims; + ensure_periodicity = temp_periodicity; } bool fragment_stats::has_non_medium_material() { @@ -1770,7 +1777,7 @@ void fragment_stats::update_stats_from_material(material_type mat, size_t pixels } } -void fragment_stats::tree_search(geom_box_tree t, std::vector &overlaps) { +void fragment_stats::compute_overlaps_from_tree(geom_box_tree t, std::vector &overlaps) { if (!t || !geom_boxes_intersect(&t->b, &box)) return; @@ -1787,6 +1794,7 @@ void fragment_stats::tree_search(geom_box_tree t, std::vector &overlaps) for (int j = 0; j < geom.num_items; ++j) { if (t->objects[i].o == geom.items + j) { overlaps[j] += overlap; + break; } } } @@ -1794,8 +1802,8 @@ void fragment_stats::tree_search(geom_box_tree t, std::vector &overlaps) } } - tree_search(t->t1, overlaps); - tree_search(t->t2, overlaps); + compute_overlaps_from_tree(t->t1, overlaps); + compute_overlaps_from_tree(t->t2, overlaps); } void fragment_stats::compute_stats() { @@ -1806,14 +1814,16 @@ void fragment_stats::compute_stats() { } else { std::vector overlaps(geom.num_items, 0.0); - tree_search(geom_tree, overlaps); + compute_overlaps_from_tree(geom_tree, overlaps); + + size_t total_obj_pixels = 0; + bool anisotropic_pixels_already_added = false; for (int i = 0; i < geom.num_items; ++i) { geometric_object *go = &geom.items[i]; double overlap = overlaps[i]; - bool anisotropic_pixels_already_added = false; - if (eps_averaging) { + if (eps_averaging && !anisotropic_pixels_already_added) { // If the object doesn't overlap the entire box, that implies that // an object interface intercepts the box, which means we treat // the entire box as anisotropic. This method could give some false @@ -1833,16 +1843,17 @@ void fragment_stats::compute_stats() { // Count contributions from material of object size_t pixels = (size_t)ceil(overlap * num_pixels_in_box); if (pixels > 0) { + total_obj_pixels += pixels; material_type mat = (material_type)go->material; update_stats_from_material(mat, pixels, anisotropic_pixels_already_added); } - - // Count contributions from default_material - size_t default_material_pixels = num_pixels_in_box - pixels; - if (default_material_pixels > 0) { - update_stats_from_material((material_type)default_material, default_material_pixels, - anisotropic_pixels_already_added); - } + } + // Count contributions from default_material + // TODO: This is wrong if objects overlap each other. + size_t default_material_pixels = num_pixels_in_box - total_obj_pixels; + if (default_material_pixels > 0) { + update_stats_from_material((material_type)default_material, default_material_pixels, + anisotropic_pixels_already_added); } } } diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index cb0c84ab6..88d93e4db 100644 --- a/src/meepgeom.hpp +++ b/src/meepgeom.hpp @@ -118,8 +118,7 @@ struct fragment_stats { void compute_dft_stats(); void compute_pml_stats(); void compute_absorber_stats(); - // TODO rename - void tree_search(geom_box_tree t, std::vector &overlaps); + void compute_overlaps_from_tree(geom_box_tree t, std::vector &overlaps); }; fragment_stats From 404d6104d7fc4ff232960756068a6d085c15fb9f Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Tue, 8 Oct 2019 20:38:47 -0500 Subject: [PATCH 4/5] Working much better now --- python/tests/fragment_stats.py | 10 ++-- src/meepgeom.cpp | 86 ++++++++++++++-------------------- src/meepgeom.hpp | 2 +- 3 files changed, 40 insertions(+), 58 deletions(-) diff --git a/python/tests/fragment_stats.py b/python/tests/fragment_stats.py index 6d569a8d9..9dd3a4921 100644 --- a/python/tests/fragment_stats.py +++ b/python/tests/fragment_stats.py @@ -58,7 +58,7 @@ def get_fragment_stats(self, block_size, cell_size, dims, box_center=mp.Vector3( gv = sim._create_grid_volume(False) stats = sim._compute_fragment_stats(gv) - + stats.print_stats() return stats def _test_1d(self, sym, pml=[]): @@ -308,17 +308,17 @@ def test_3d_cell_smaller_than_minimum_fragment_size(self): self.assertEqual(fs.num_pixels_in_box, 1000) def test_geom_box_tree(self): - block_size = mp.Vector3(3, 3) - mat = mp.Medium(index=13) + mat = mp.Medium(index=13, E_chi2_diag=mp.Vector3(1, 1)) geometry = [ mp.Block(size=block_size, center=mp.Vector3(-3, 3), material=mat), mp.Block(size=block_size, center=mp.Vector3(3, 3), material=mat), mp.Block(size=block_size, center=mp.Vector3(3, -3), material=mat), - mp.Block(size=block_size, center=mp.Vector3(-3, -3), material=mat) + mp.Block(size=block_size, center=mp.Vector3(-3, -3), material=mat), + mp.Block(size=mp.Vector3(9, 1.5), center=mp.Vector3(0, 0), material=mat) ] fs = self.get_fragment_stats(mp.Vector3(), mp.Vector3(9, 9), 2, geom=geometry) - self.check_stats(fs, 32400, 0, 0, 0, 0) + self.check_stats(fs, 8100, 0, 9900, 0, 0) self.assertEqual(fs.num_pixels_in_box, 8100) diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index 3081859b2..6bdaa6bcf 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -1741,10 +1741,10 @@ void fragment_stats::init_libctl(material_type default_mat, bool ensure_per, mee dimensions = 3; ensure_periodicity = 0; geom_tree = create_geom_box_tree0(*geom_, cell_box); - // if (meep::am_master()) display_geom_box_tree(2, geom_tree); - // int depth, nobjects; - // geom_box_tree_stats(geom_tree, &depth, &nobjects); - // master_printf("depth: %d, nobjects: %d\n", depth, nobjects); + if (meep::am_master()) display_geom_box_tree(2, geom_tree); + int depth, nobjects; + geom_box_tree_stats(geom_tree, &depth, &nobjects); + master_printf("depth: %d, nobjects: %d\n", depth, nobjects); dimensions = temp_dims; ensure_periodicity = temp_periodicity; } @@ -1777,7 +1777,7 @@ void fragment_stats::update_stats_from_material(material_type mat, size_t pixels } } -void fragment_stats::compute_overlaps_from_tree(geom_box_tree t, std::vector &overlaps) { +void fragment_stats::compute_overlaps_from_tree(geom_box_tree t, bool &anisotropic_pixels_already_added) { if (!t || !geom_boxes_intersect(&t->b, &box)) return; @@ -1791,10 +1791,32 @@ void fragment_stats::compute_overlaps_from_tree(geom_box_tree t, std::vectorobjects[i].shiftby)}; double overlap = box_overlap_with_object(shifted_box, *t->objects[i].o, tol, maxeval); if (overlap > 0.0) { - for (int j = 0; j < geom.num_items; ++j) { - if (t->objects[i].o == geom.items + j) { - overlaps[j] += overlap; - break; + // Count contributions from material of object + size_t num_pixels_in_intersection_box = get_pixels_in_box(&intersection); + size_t pixels = (size_t)ceil(overlap * num_pixels_in_intersection_box); + if (pixels > 0) { + material_type mat = (material_type)t->objects[i].o->material; + update_stats_from_material(mat, pixels, false); + } + size_t default_material_pixels = num_pixels_in_intersection_box - pixels; + if (default_material_pixels > 0) { + update_stats_from_material((material_type)default_material, default_material_pixels, false); + } + + if (eps_averaging && !anisotropic_pixels_already_added) { + // If the object doesn't overlap the entire box, that implies that + // an object interface intercepts the box, which means we treat + // the entire box as anisotropic. This method could give some false + // positives if there is another object with the same material behind + // the current object, but in practice it is probably reasonable to + // assume that there is a material interface somwhere in the box so + // we won't worry about fancier edge-detection methods for now. + if (overlap != 1.0) { + anisotropic_pixels_already_added = true; + num_anisotropic_eps_pixels += num_pixels_in_box; + if (mu_not_1(t->objects[i].o->material)) { + num_anisotropic_mu_pixels += num_pixels_in_box; + } } } } @@ -1802,8 +1824,8 @@ void fragment_stats::compute_overlaps_from_tree(geom_box_tree t, std::vectort1, overlaps); - compute_overlaps_from_tree(t->t2, overlaps); + compute_overlaps_from_tree(t->t1, anisotropic_pixels_already_added); + compute_overlaps_from_tree(t->t2, anisotropic_pixels_already_added); } void fragment_stats::compute_stats() { @@ -1813,48 +1835,8 @@ void fragment_stats::compute_stats() { update_stats_from_material((material_type)default_material, num_pixels_in_box); } else { - std::vector overlaps(geom.num_items, 0.0); - compute_overlaps_from_tree(geom_tree, overlaps); - - size_t total_obj_pixels = 0; bool anisotropic_pixels_already_added = false; - - for (int i = 0; i < geom.num_items; ++i) { - geometric_object *go = &geom.items[i]; - double overlap = overlaps[i]; - - if (eps_averaging && !anisotropic_pixels_already_added) { - // If the object doesn't overlap the entire box, that implies that - // an object interface intercepts the box, which means we treat - // the entire box as anisotropic. This method could give some false - // positives if there is another object with the same material behind - // the current object, but in practice it is probably reasonable to - // assume that there is a material interface somwhere in the box so - // we won't worry about fancier edge-detection methods for now. - if (overlap != 1.0) { - anisotropic_pixels_already_added = true; - num_anisotropic_eps_pixels += num_pixels_in_box; - if (mu_not_1(go->material)) { - num_anisotropic_mu_pixels += num_pixels_in_box; - } - } - } - - // Count contributions from material of object - size_t pixels = (size_t)ceil(overlap * num_pixels_in_box); - if (pixels > 0) { - total_obj_pixels += pixels; - material_type mat = (material_type)go->material; - update_stats_from_material(mat, pixels, anisotropic_pixels_already_added); - } - } - // Count contributions from default_material - // TODO: This is wrong if objects overlap each other. - size_t default_material_pixels = num_pixels_in_box - total_obj_pixels; - if (default_material_pixels > 0) { - update_stats_from_material((material_type)default_material, default_material_pixels, - anisotropic_pixels_already_added); - } + compute_overlaps_from_tree(geom_tree, anisotropic_pixels_already_added); } } diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index 88d93e4db..61b4e78e2 100644 --- a/src/meepgeom.hpp +++ b/src/meepgeom.hpp @@ -118,7 +118,7 @@ struct fragment_stats { void compute_dft_stats(); void compute_pml_stats(); void compute_absorber_stats(); - void compute_overlaps_from_tree(geom_box_tree t, std::vector &overlaps); + void compute_overlaps_from_tree(geom_box_tree t, bool &anisotropic_pixels_already_added); }; fragment_stats From 6778a696cff07c0ae4b8df065ce2a75fb84dad5b Mon Sep 17 00:00:00 2001 From: Chris Hogan Date: Wed, 9 Oct 2019 20:20:33 -0500 Subject: [PATCH 5/5] Remove debug code and fix anisotropic bug --- python/tests/fragment_stats.py | 3 +- src/meepgeom.cpp | 51 +++++++++++++++------------------- 2 files changed, 24 insertions(+), 30 deletions(-) diff --git a/python/tests/fragment_stats.py b/python/tests/fragment_stats.py index 9dd3a4921..f10f38dd2 100644 --- a/python/tests/fragment_stats.py +++ b/python/tests/fragment_stats.py @@ -58,7 +58,6 @@ def get_fragment_stats(self, block_size, cell_size, dims, box_center=mp.Vector3( gv = sim._create_grid_volume(False) stats = sim._compute_fragment_stats(gv) - stats.print_stats() return stats def _test_1d(self, sym, pml=[]): @@ -318,7 +317,7 @@ def test_geom_box_tree(self): mp.Block(size=mp.Vector3(9, 1.5), center=mp.Vector3(0, 0), material=mat) ] fs = self.get_fragment_stats(mp.Vector3(), mp.Vector3(9, 9), 2, geom=geometry) - self.check_stats(fs, 8100, 0, 9900, 0, 0) + self.check_stats(fs, 8100, 0, 9910, 0, 0) self.assertEqual(fs.num_pixels_in_box, 8100) diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index 6bdaa6bcf..e8506bb84 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -1732,21 +1732,14 @@ void fragment_stats::init_libctl(material_type default_mat, bool ensure_per, mee geom_initialize(); default_material = default_mat; geometry_center = cell_center; - dimensions = meep::number_of_directions(gv->dim); geometry_lattice.size = cell_size; geom_fix_object_list(*geom_); geom_box cell_box = make_box_from_cell(cell_size); - int temp_dims = dimensions; - int temp_periodicity = ensure_periodicity; dimensions = 3; ensure_periodicity = 0; geom_tree = create_geom_box_tree0(*geom_, cell_box); - if (meep::am_master()) display_geom_box_tree(2, geom_tree); - int depth, nobjects; - geom_box_tree_stats(geom_tree, &depth, &nobjects); - master_printf("depth: %d, nobjects: %d\n", depth, nobjects); - dimensions = temp_dims; - ensure_periodicity = temp_periodicity; + dimensions = meep::number_of_directions(gv->dim); + ensure_periodicity = ensure_per; } bool fragment_stats::has_non_medium_material() { @@ -1790,34 +1783,36 @@ void fragment_stats::compute_overlaps_from_tree(geom_box_tree t, bool &anisotrop geom_box shifted_box = {vector3_minus(intersection.low, t->objects[i].shiftby), vector3_minus(intersection.high, t->objects[i].shiftby)}; double overlap = box_overlap_with_object(shifted_box, *t->objects[i].o, tol, maxeval); + + if (eps_averaging && !anisotropic_pixels_already_added) { + // If the object doesn't overlap the entire box, that implies that + // an object interface intercepts the box, which means we treat + // the entire box as anisotropic. This method could give some false + // positives if there is another object with the same material behind + // the current object, but in practice it is probably reasonable to + // assume that there is a material interface somwhere in the box so + // we won't worry about fancier edge-detection methods for now. + if (overlap != 1.0) { + anisotropic_pixels_already_added = true; + num_anisotropic_eps_pixels += num_pixels_in_box; + if (mu_not_1(t->objects[i].o->material)) { + num_anisotropic_mu_pixels += num_pixels_in_box; + } + } + } + if (overlap > 0.0) { // Count contributions from material of object size_t num_pixels_in_intersection_box = get_pixels_in_box(&intersection); size_t pixels = (size_t)ceil(overlap * num_pixels_in_intersection_box); if (pixels > 0) { material_type mat = (material_type)t->objects[i].o->material; - update_stats_from_material(mat, pixels, false); + update_stats_from_material(mat, pixels, anisotropic_pixels_already_added); } size_t default_material_pixels = num_pixels_in_intersection_box - pixels; if (default_material_pixels > 0) { - update_stats_from_material((material_type)default_material, default_material_pixels, false); - } - - if (eps_averaging && !anisotropic_pixels_already_added) { - // If the object doesn't overlap the entire box, that implies that - // an object interface intercepts the box, which means we treat - // the entire box as anisotropic. This method could give some false - // positives if there is another object with the same material behind - // the current object, but in practice it is probably reasonable to - // assume that there is a material interface somwhere in the box so - // we won't worry about fancier edge-detection methods for now. - if (overlap != 1.0) { - anisotropic_pixels_already_added = true; - num_anisotropic_eps_pixels += num_pixels_in_box; - if (mu_not_1(t->objects[i].o->material)) { - num_anisotropic_mu_pixels += num_pixels_in_box; - } - } + update_stats_from_material((material_type)default_material, default_material_pixels, + anisotropic_pixels_already_added); } } }