get_epsilon_grid function for evaluating ε on user-specified grid (Na…

…noComp#1522)

* get_epsilon_grid function for evaluating ε on user-specified grid

* refactor fragment_stats::init_libctl out of class fragment_stats

* clean up

* initialize coordinates of user-specified volume with empty dimensions to zero

* remove if statment from get_epsilon_grid

* modify find_array_min_max to accept arguments by reference rather than return complex double

* fixes

doc/docs/Python_User_Interface.md

@@ -666,6 +666,29 @@ frequency-independent part of $\mu$ (the $\omega\to\infty$ limit).
 
 </div>
 
+<a id="Simulation.get_epsilon_grid"></a>
+
+<div class="class_members" markdown="1">
+
+```python
+def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None):
+```
+
+<div class="method_docstring" markdown="1">
+
+Given three 1d NumPy arrays (`xtics`,`ytics`,`ztics`) which define the coordinates of a Cartesian
+grid anywhere within the cell volume, compute the trace of the $\varepsilon$ tensor from the `geometry`
+exactly at each grid point. (For [`MaterialGrid`](#materialgrid)s, the $\varepsilon$ at each grid
+point is computed using bilinear interpolation from the nearest `MaterialGrid` points and possibly also
+projected to form a level set.) Note that this is different from `get_epsilon_point` which computes
+$\varepsilon$ by bilinearly interpolating from the nearest Yee grid points. This function is useful for
+sampling the material geometry to any arbitrary resolution. The return value is a NumPy array with shape
+equivalent to `numpy.meshgrid(xtics,ytics,ztics)`. Empty dimensions are collapsed.
+
+</div>
+
+</div>
+
 <a id="Simulation.initialize_field"></a>
 
 <div class="class_members" markdown="1">

doc/docs/Python_User_Interface.md.in

@@ -99,6 +99,7 @@ Meep supports a large number of functions to perform computations on the fields.
 @@ Simulation.get_field_point @@
 @@ Simulation.get_epsilon_point @@
 @@ Simulation.get_mu_point @@
+@@ Simulation.get_epsilon_grid @@
 @@ Simulation.initialize_field @@
 @@ Simulation.add_dft_fields @@
 @@ Simulation.flux_in_box @@

python/meep.i

@@ -1087,6 +1087,40 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
 %apply int INPLACE_ARRAY1[ANY] { int [3] };
 %apply double INPLACE_ARRAY1[ANY] { double [3] };
 
+// typemaps for meep_geom::get_epsilon_grid
+
+%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* xtics {
+    $1 = is_array($input);
+}
+
+%typemap(in, fragment="NumPy_Macros") double* xtics {
+    $1 = (double *)array_data($input);
+}
+
+%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* ytics {
+    $1 = is_array($input);
+}
+
+%typemap(in, fragment="NumPy_Macros") double* ytics {
+    $1 = (double *)array_data($input);
+}
+
+%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* ztics {
+    $1 = is_array($input);
+}
+
+%typemap(in, fragment="NumPy_Macros") double* ztics {
+    $1 = (double *)array_data($input);
+}
+
+%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* grid_vals {
+    $1 = is_array($input);
+}
+
+%typemap(in, fragment="NumPy_Macros") double* grid_vals {
+    $1 = (double *)array_data($input);
+}
+
 // typemap for solve_cw:
 
 %typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") std::complex<double>* eigfreq {
@@ -1872,8 +1906,8 @@ meep::structure *create_structure_and_set_materials(vector3 cell_size,
     meep_geom::fragment_stats::resolution = gv.a;
     meep_geom::fragment_stats::dims = gv.dim;
     meep_geom::fragment_stats::split_chunks_evenly = split_chunks_evenly;
-    meep_geom::fragment_stats::init_libctl(_default_material, _ensure_periodicity,
-                                           &gv, cell_size, center, &gobj_list);
+    meep_geom::init_libctl(_default_material, _ensure_periodicity,
+                           &gv, cell_size, center, &gobj_list);
 
     if (output_chunk_costs) {
          meep::volume thev = gv.surroundings();
@@ -1927,4 +1961,29 @@ meep::structure *create_structure_and_set_materials(vector3 cell_size,
 
     return s;
 }
+
+void _get_epsilon_grid(geometric_object_list gobj_list,
+                       meep_geom::material_type_list mlist,
+                       meep_geom::material_type _default_material,
+                       bool _ensure_periodicity,
+                       meep::grid_volume gv,
+                       vector3 cell_size,
+                       vector3 cell_center,
+                       int nx, double *xtics,
+                       int ny, double *ytics,
+                       int nz, double *ztics,
+                       double *grid_vals) {
+     meep_geom::get_epsilon_grid(gobj_list,
+                                 mlist,
+                                 _default_material,
+                                 _ensure_periodicity,
+                                 gv,
+                                 cell_size,
+                                 cell_center,
+                                 nx, xtics,
+                                 ny, ytics,
+                                 nz, ztics,
+                                 grid_vals);
+}
+
 %}

python/simulation.py

@@ -2080,6 +2080,31 @@ def get_mu_point(self, pt, frequency=0):
         v3 = py_v3_to_vec(self.dimensions, pt, self.is_cylindrical)
         return self.fields.get_mu(v3,frequency)
 
+    def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None):
+        """
+        Given three 1d NumPy arrays (`xtics`,`ytics`,`ztics`) which define the coordinates of a Cartesian
+        grid anywhere within the cell volume, compute the trace of the $\\varepsilon$ tensor from the `geometry`
+        exactly at each grid point. (For [`MaterialGrid`](#materialgrid)s, the $\\varepsilon$ at each grid
+        point is computed using bilinear interpolation from the nearest `MaterialGrid` points and possibly also
+        projected to form a level set.) Note that this is different from `get_epsilon_point` which computes
+        $\\varepsilon$ by bilinearly interpolating from the nearest Yee grid points. This function is useful for
+        sampling the material geometry to any arbitrary resolution. The return value is a NumPy array with shape
+        equivalent to `numpy.meshgrid(xtics,ytics,ztics)`. Empty dimensions are collapsed.
+        """
+        grid_vals = np.squeeze(np.empty((len(xtics), len(ytics), len(ztics))))
+        gv = self._create_grid_volume(False)
+        mp._get_epsilon_grid(self.geometry,
+                             self.extra_materials,
+                             self.default_material,
+                             self.ensure_periodicity and not not self.k_point,
+                             gv,
+                             self.cell_size, self.geometry_center,
+                             len(xtics), xtics,
+                             len(ytics), ytics,
+                             len(ztics), ztics,
+                             grid_vals)
+        return grid_vals
+
     def get_filename_prefix(self):
         """
         Return the current prefix string that is prepended, by default, to all file names.

src/meepgeom.cpp

@@ -565,7 +565,6 @@ geom_epsilon::geom_epsilon(geometric_object_list g, material_type_list mlist,
     if (num_print < geometry.num_items && meep::verbosity > 0)
       master_printf("%*s...(+ %d objects not shown)...\n", 5, "", geometry.num_items - num_print);
   }
-
   geom_fix_object_list(geometry);
   geom_box box = gv2box(v);
   geometry_tree = create_geom_box_tree0(geometry, box);
@@ -1926,7 +1925,7 @@ fragment_stats compute_fragment_stats(
   fragment_stats::extra_materials = extra_materials_;
   fragment_stats::eps_averaging = eps_averaging;
 
-  fragment_stats::init_libctl(default_mat, ensure_per, gv, cell_size, cell_center, &geom_);
+  init_libctl(default_mat, ensure_per, gv, cell_size, cell_center, &geom_);
   geom_box box = make_box_from_cell(cell_size);
   fragment_stats stats = init_stats(box, tol, maxeval, gv);
   stats.compute();
@@ -1941,9 +1940,9 @@ fragment_stats::fragment_stats(geom_box &bx)
   num_pixels_in_box = get_pixels_in_box(&bx);
 }
 
-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_) {
+void init_libctl(material_type default_mat, bool ensure_per, meep::grid_volume *gv,
+                 vector3 cell_size, vector3 cell_center,
+                 geometric_object_list *geom_) {
   geom_initialize();
   default_material = default_mat;
   ensure_periodicity = ensure_per;
@@ -2567,4 +2566,48 @@ void material_grids_addgradient(meep::realnum *v, size_t ng, std::complex<double
   }
 }
 
+static void find_array_min_max(int n, const double *data, double &min_val, double &max_val) {
+  min_val = data[0];
+  max_val = data[0];
+  for (int i = 1; i < n; ++i) {
+    if (data[i] < min_val)
+      min_val = data[i];
+    if (data[i] > max_val)
+      max_val = data[i];
+  }
+  return;
+}
+
+void get_epsilon_grid(geometric_object_list gobj_list,
+                      material_type_list mlist,
+                      material_type _default_material,
+                      bool _ensure_periodicity,
+                      meep::grid_volume gv,
+                      vector3 cell_size,
+                      vector3 cell_center,
+                      int nx, const double *x,
+                      int ny, const double *y,
+                      int nz, const double *z,
+                      double *grid_vals) {
+  double min_val[3], max_val[3];
+  for (int n = 0; n < 3; ++n) {
+    int ndir = (n == 0) ? nx : ((n == 1) ? ny : nz);
+    if (ndir < 1) meep::abort("get_epsilon_grid: ndir < 1.");
+    const double *adir = (n == 0) ? x : ((n == 1) ? y : z);
+    find_array_min_max(ndir, adir, min_val[n], max_val[n]);
+  }
+  const meep::volume vol(meep::vec(min_val[0],min_val[1],min_val[2]),
+                         meep::vec(max_val[0],max_val[1],max_val[2]));
+  init_libctl(_default_material, _ensure_periodicity, &gv,
+              cell_size, cell_center, &gobj_list);
+  dim = gv.dim;
+  geom_epsilon geps(gobj_list, mlist, vol);
+  for (int i = 0; i < nx; ++i)
+    for (int j = 0; j < ny; ++j)
+      for (int k = 0; k < nz; ++k)
+        /* obtain the trace of the \varepsilon tensor for each
+           grid point in row-major order (the order used by NumPy) */
+        grid_vals[k + nz*(j + ny*i)] = geps.chi1p1(meep::E_stuff, meep::vec(x[i],y[j],z[k]));
+}
+
 } // namespace meep_geom

src/meepgeom.hpp

@@ -77,9 +77,6 @@ struct fragment_stats {
   static bool eps_averaging;
 
   static bool has_non_medium_material();
-  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);
 
   size_t num_anisotropic_eps_pixels;
   size_t num_anisotropic_mu_pixels;
@@ -197,6 +194,20 @@ bool is_medium(void *md, medium_struct **m);
 bool is_metal(meep::field_type ft, const material_type *material);
 void check_offdiag(medium_struct *m);
 geom_box gv2box(const meep::volume &v);
+void get_epsilon_grid(geometric_object_list gobj_list,
+                      material_type_list mlist,
+                      material_type _default_material,
+                      bool _ensure_periodicity,
+                      meep::grid_volume gv,
+                      vector3 cell_size,
+                      vector3 cell_center,
+                      int nx, const double *x,
+                      int ny, const double *y,
+                      int nz, const double *z,
+                      double *grid_vals);
+void init_libctl(material_type default_mat, bool ensure_per,
+                 meep::grid_volume *gv, vector3 cell_size, vector3 cell_center,
+                 geometric_object_list *geom_list);
 
 /***************************************************************/
 // material grid functions