Treat fragments with object interfaces as anisotropic (#975)

* Treat fragments with object interfaces as anisotropic

* Only count full box as anisotropic if eps_averaging=True

python/simulation.py

@@ -990,6 +990,7 @@ def _compute_fragment_stats(self, gv):
             self.subpixel_tol,
             self.subpixel_maxeval,
             self.ensure_periodicity,
+            self.eps_averaging
         )
 
         mirror_symmetries = [sym for sym in self.symmetries if isinstance(sym, Mirror)]

python/tests/fragment_stats.py

@@ -75,8 +75,8 @@ def _test_1d(self, sym, pml=[]):
 
         sym_factor = 2 if sym else 1
         self.check_stats(fs,
-                         a_eps=100 / sym_factor,
-                         a_mu=100 / sym_factor,
+                         a_eps=300 / sym_factor,
+                         a_mu=300 / sym_factor,
                          nonlin=300 / sym_factor,
                          susc=300 / sym_factor,
                          cond=300 / sym_factor)
@@ -101,7 +101,7 @@ def test_1d_with_overlap(self):
         # A z=30 cell, with a block covering the middle 20 units.
         mat = mp.Medium(H_susceptibilities=[mp.DrudeSusceptibility()])
         fs = self.get_fragment_stats(mp.Vector3(z=20), mp.Vector3(z=30), 1, def_mat=mat)
-        self.check_stats(fs, a_eps=200, a_mu=200, nonlin=600, susc=700, cond=600)
+        self.check_stats(fs, a_eps=300, a_mu=300, nonlin=600, susc=700, cond=600)
 
     def test_1d_with_partial_fragment(self):
         # A cell with z=26, with a 16 unit block in the center
@@ -114,7 +114,7 @@ def test_1d_with_partial_fragment(self):
         ])
         fs = self.get_fragment_stats(mp.Vector3(z=16), mp.Vector3(z=26), 1, dft_vecs=dft_vecs)
 
-        self.check_stats(fs, a_eps=160, a_mu=160, nonlin=480, susc=480, cond=480)
+        self.check_stats(fs, a_eps=260, a_mu=260, nonlin=480, susc=480, cond=480)
         # Check dft stats
         self.assertEqual(fs.num_dft_pixels, 10400)
 
@@ -152,8 +152,8 @@ def _test_2d(self, sym, pml=[]):
         # Middle fragment contains entire block
         sym_factor = 4 if sym else 1
         self.check_stats(fs,
-                         a_eps=10000 / sym_factor,
-                         a_mu=10000 / sym_factor,
+                         a_eps=90000 / sym_factor,
+                         a_mu=90000 / sym_factor,
                          nonlin=30000 / sym_factor,
                          susc=30000 / sym_factor,
                          cond=30000 / sym_factor)
@@ -218,8 +218,8 @@ def _test_3d(self, sym, pml=[]):
 
         sym_factor = 8 if sym else 1
         self.check_stats(fs,
-                         a_eps=1000000 / sym_factor,
-                         a_mu=1000000 / sym_factor,
+                         a_eps=27000000 / sym_factor,
+                         a_mu=27000000 / sym_factor,
                          nonlin=3000000 / sym_factor,
                          susc=3000000 / sym_factor,
                          cond=3000000 / sym_factor)
@@ -265,7 +265,7 @@ def test_cyl(self):
         self.assertEqual(fs.box.low.z, -15)
         self.assertEqual(fs.box.high.x, 15)
         self.assertEqual(fs.box.high.z, 15)
-        self.check_stats(fs, a_eps=10000, a_mu=10000, nonlin=30000, susc=30000, cond=30000)
+        self.check_stats(fs, a_eps=90000, a_mu=90000, nonlin=30000, susc=30000, cond=30000)
         self.assertEqual(fs.num_dft_pixels, 2040000)
 
     def test_no_geometry(self):

src/meepgeom.cpp

@@ -1647,6 +1647,7 @@ std::vector<meep::volume> fragment_stats::pml_2d_vols;
 std::vector<meep::volume> fragment_stats::pml_3d_vols;
 std::vector<meep::volume> fragment_stats::absorber_vols;
 bool fragment_stats::split_chunks_evenly = false;
+bool fragment_stats::eps_averaging = false;
 
 static geom_box make_box_from_cell(vector3 cell_size) {
   double edgex = cell_size.x / 2;
@@ -1693,14 +1694,15 @@ fragment_stats compute_fragment_stats(
     material_type default_mat, std::vector<dft_data> dft_data_list_,
     std::vector<meep::volume> pml_1d_vols_, std::vector<meep::volume> pml_2d_vols_,
     std::vector<meep::volume> pml_3d_vols_, std::vector<meep::volume> absorber_vols_, double tol,
-    int maxeval, bool ensure_per) {
+    int maxeval, bool ensure_per, bool eps_averaging) {
 
   fragment_stats::geom = geom_;
   fragment_stats::dft_data_list = dft_data_list_;
   fragment_stats::pml_1d_vols = pml_1d_vols_;
   fragment_stats::pml_2d_vols = pml_2d_vols_;
   fragment_stats::pml_3d_vols = pml_3d_vols_;
   fragment_stats::absorber_vols = absorber_vols_;
+  fragment_stats::eps_averaging = eps_averaging;
 
   fragment_stats::init_libctl(default_mat, ensure_per, gv, cell_size, cell_center, &geom_);
   geom_box box = make_box_from_cell(cell_size);
@@ -1738,11 +1740,14 @@ bool fragment_stats::has_non_medium_material() {
   return false;
 }
 
-void fragment_stats::update_stats_from_material(material_type mat, size_t pixels) {
+void fragment_stats::update_stats_from_material(material_type mat, size_t pixels,
+                                                bool anisotropic_pixels_already_added) {
   switch (mat->which_subclass) {
     case material_data::MEDIUM: {
       medium_struct *med = &mat->medium;
-      count_anisotropic_pixels(med, pixels);
+      if (!anisotropic_pixels_already_added) {
+        count_anisotropic_pixels(med, pixels);
+      }
       count_nonlinear_pixels(med, pixels);
       count_susceptibility_pixels(med, pixels);
       count_nonzero_conductivity_pixels(med, pixels);
@@ -1765,17 +1770,36 @@ void fragment_stats::compute_stats() {
     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;
+        }
+      }
+    }
+
     // 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);
+      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);
+      update_stats_from_material((material_type)default_material, default_material_pixels,
+                                 anisotropic_pixels_already_added);
     }
   }
 }

src/meepgeom.hpp

@@ -73,6 +73,7 @@ struct fragment_stats {
   static std::vector<meep::volume> pml_3d_vols;
   static std::vector<meep::volume> absorber_vols;
   static bool split_chunks_evenly;
+  static bool eps_averaging;
 
   static bool has_non_medium_material();
   static void init_libctl(meep_geom::material_type default_mat, bool ensure_per,
@@ -105,7 +106,8 @@ struct fragment_stats {
   void print_stats() const;
 
 private:
-  void update_stats_from_material(material_type mat, size_t pixels);
+  void update_stats_from_material(material_type mat, size_t pixels,
+                                  bool anisotropic_pixels_already_added=false);
   void compute_stats();
   void count_anisotropic_pixels(medium_struct *med, size_t pixels);
   void count_nonlinear_pixels(medium_struct *med, size_t pixels);
@@ -122,7 +124,7 @@ compute_fragment_stats(geometric_object_list geom, meep::grid_volume *gv, vector
                        std::vector<dft_data> dft_data_list, std::vector<meep::volume> pml_1d_vols,
                        std::vector<meep::volume> pml_2d_vols, std::vector<meep::volume> pml_3d_vols,
                        std::vector<meep::volume> absorber_vols, double tol, int maxeval,
-                       bool ensure_per);
+                       bool ensure_per, bool eps_averaging);
 
 /***************************************************************/
 /* these routines create and append absorbing layers to an     */