Enable Eigenmode Features with Dispersive Materials (NanoComp#919)

* first stab

* begin with test files

* bug fixes

* update test

* found one bug

* bug fixes, h5 output, and array-slice output

* fix serial bugs

* more serial fixes

* fix mu issue

* disable h5 support

* big fix

* generalize matrix inverse and fix parallel bug

* fix memory issue

* add tutorial

* minor fixes

* refactor

* check1

* case2

* add docs

* cleanup

* check3

* fix test

* 1 more fix!

* updates

* det fix

doc/docs/Python_User_Interface.md

@@ -1017,9 +1017,9 @@ Sets the condition of the boundary on the specified side in the specified direct
 —
 Given a `component` or `derived_component` constant `c` and a `Vector3` `pt`, returns the value of that component at that point.
 
-**`get_epsilon_point(pt)`**
+**`get_epsilon_point(pt, omega=0)`**
 —
-Equivalent to `get_field_point(mp.Dielectric, pt)`.
+Given a frequency `omega` and a `Vector3` `pt`, returns the average eigenvalue of the permittivity tensor at that location and frequency.
 
 **`initialize_field(c, func)`**
 —
@@ -1800,7 +1800,7 @@ See also [Field Functions](Field_Functions.md), and [Synchronizing the Magnetic
 The output functions described above write the data for the fields and materials for the entire cell to an HDF5 file. This is useful for post-processing as you can later read in the HDF5 file to obtain field/material data as a NumPy array. However, in some cases it is convenient to bypass the disk altogether to obtain the data *directly* in the form of a NumPy array without writing/reading HDF5 files. Additionally, you may want the field/material data on just a subregion (or slice) of the entire volume. This functionality is provided by the `get_array` method which takes as input a subregion of the cell and the field/material component. The method returns a NumPy array containing values of the field/material at the current simulation time.
 
 ```python
- get_array(vol=None, center=None, size=None, component=mp.Ez, cmplx=False, arr=None)
+ get_array(vol=None, center=None, size=None, component=mp.Ez, cmplx=False, arr=None, omega=0)
 ```
 
 with the following input parameters:
@@ -1815,7 +1815,9 @@ with the following input parameters:
 
 + `arr`: optional field to pass a pre-allocated NumPy array of the correct size, which will be overwritten with the field/material data instead of allocating a new array.  Normally, this will be the array returned from a previous call to `get_array` for a similar slice, allowing one to re-use `arr` (e.g., when fetching the same slice repeatedly at different times).
 
-For convenience, the following wrappers for `get_array` over the entire cell are available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, `get_tot_pwr()`, `get_Xfield()`, `get_Xfield_x()`, `get_Xfield_y()`, `get_Xfield_z()`, `get_Xfield_r()`, `get_Xfield_p()` where `X` is one of `h`, `b`, `e`, `d`, or `s`. The routines `get_Xfield_*` all return an array type consistent with the fields (real or complex).
++ `omega`: optional frequency point over which the average eigenvalue of the dielectric and permeability tensors are evaluated (defaults to 0).
+
+For convenience, the following wrappers for `get_array` over the entire cell are available: `get_epsilon()`, `get_mu()`, `get_hpwr()`, `get_dpwr()`, `get_tot_pwr()`, `get_Xfield()`, `get_Xfield_x()`, `get_Xfield_y()`, `get_Xfield_z()`, `get_Xfield_r()`, `get_Xfield_p()` where `X` is one of `h`, `b`, `e`, `d`, or `s`. The routines `get_Xfield_*` all return an array type consistent with the fields (real or complex). The routines `get_epsilon()` and `get_mu()` accept the optional omega parameter (defaults to 0).
 
 **Note on array-slice dimensions:** The routines `get_epsilon`, `get_Xfield_z`, etc. use as default `size=meep.Simulation.fields.total_volume()` which for simulations involving Bloch-periodic boundaries (via `k_point`) will result in arrays that have slightly *different* dimensions than e.g. `get_array(center=meep.Vector3(), size=cell_size, component=meep.Dielectric`, etc. (i.e., the slice spans the entire cell volume `cell_size`). Neither of these approaches is "wrong", they are just slightly different methods of fetching the boundaries. The key point is that if you pass the same value for the `size` parameter, or use the default, the slicing routines always give you the same-size array for all components. You should *not* try to predict the exact size of these arrays; rather, you should simply rely on Meep's output.
 

python/Makefile.am

@@ -18,12 +18,14 @@ endif # WITH_MPI
 
 if WITH_MPB
   BINARY_GRATING_TEST = $(TEST_DIR)/binary_grating.py
+  DISPERSIVE_EIGENMODE_TEST = $(TEST_DIR)/dispersive_eigenmode.py
   KDOM_TEST = $(TEST_DIR)/kdom.py
   MODE_COEFFS_TEST = $(TEST_DIR)/mode_coeffs.py
   MODE_DECOMPOSITION_TEST = $(TEST_DIR)/mode_decomposition.py
   WVG_SRC_TEST = $(TEST_DIR)/wvg_src.py
 else
   BINARY_GRATING_TEST =
+  DISPERSIVE_EIGENMODE_TEST =
   KDOM_TEST =
   MODE_COEFFS_TEST =
   MODE_DECOMPOSITION_TEST =
@@ -41,6 +43,7 @@ TESTS =                                   \
     $(TEST_DIR)/cavity_farfield.py        \
     $(TEST_DIR)/chunks.py                 \
     $(TEST_DIR)/cyl_ellipsoid.py          \
+    ${DISPERSIVE_EIGENMODE_TEST}          \
     $(TEST_DIR)/dft_energy.py             \
     $(TEST_DIR)/dft_fields.py             \
     $(TEST_DIR)/faraday_rotation.py       \