From d951cfdf0bb4835b19f73a87496c6781cfb14b40 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Tue, 23 Nov 2021 18:55:07 -0800 Subject: [PATCH] plot geometry for dispersive materials without initializing structure object (#1827) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * plot geometry without initializing structure class * update docstrings * rotate epsilon grid by 90 degrees to properly orient axes * add support for dispersive ε * update markdown pages from docstrings * make frequency and resolution parameters of plot2D dictionary keys of eps_parameters * reinstate frequency parameter of plot2D and add warning that it is deprecated * fix order of frequency check --- doc/docs/Python_User_Interface.md | 45 +++-- doc/docs/Python_User_Interface.md.in | 2 +- python/meep.i | 12 +- python/simulation.py | 73 ++++--- python/tests/test_visualization.py | 2 +- python/visualization.py | 291 ++++++++++++++++----------- src/meepgeom.cpp | 27 ++- src/meepgeom.hpp | 3 +- 8 files changed, 270 insertions(+), 185 deletions(-) diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index 1904b9fc7..581dd3382 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -676,16 +676,21 @@ frequency-independent part of $\mu$ (the $\omega\to\infty$ limit).
```python -def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None): +def get_epsilon_grid(self, + xtics=None, + ytics=None, + ztics=None, + frequency=0): ```
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 +grid anywhere within the cell volume, compute the trace of the $\varepsilon(f)$ tensor at frequency +$f$ (in Meep units) from the `geometry` exactly at each grid point. `frequency` defaults to 0 which is +the instantaneous $\varepsilon$. (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. @@ -2067,7 +2072,7 @@ including outside the cell and a `near2far` object, returns the computed of fields $(E_x^1,E_y^1,E_z^1,H_x^1,H_y^1,H_z^1,E_x^2,E_y^2,E_z^2,H_x^2,H_y^2,H_z^2,...)$ in Cartesian coordinates and $(E_r^1,E_\phi^1,E_z^1,H_r^1,H_\phi^1,H_z^1,E_r^2,E_\phi^2,E_z^2,H_r^2,H_\phi^2,H_z^2,...)$ -in cylindrical coordinates for the frequencies 1,2,…,`nfreq`. +in cylindrical coordinates for the frequencies 1,2,...,`nfreq`.
@@ -2609,7 +2614,7 @@ fr = mp.FluxRegion(volume=mp.GDSII_vol(fname, layer, zmin, zmax)) ### Data Visualization -This module provides basic visualization functionality for the simulation domain. The spirit of the module is to provide functions that can be called with *no customization options whatsoever* and will do useful relevant things by default, but which can also be customized in cases where you *do* want to take the time to spruce up the output. The `Simulation` class provides the following methods: +This module provides basic visualization functionality for the simulation domain. The intent of the module is to provide functions that can be called with *no customization options whatsoever* and will do useful relevant things by default, but which can also be customized in cases where you *do* want to take the time to spruce up the output. The `Simulation` class provides the following methods: @@ -2648,14 +2653,18 @@ sim.run(...) field_func = lambda x: 20*np.log10(np.abs(x)) import matplotlib.pyplot as plt sim.plot2D(fields=mp.Ez, - field_parameters={'alpha':0.8, 'cmap':'RdBu', 'interpolation':'none', 'post_process':field_func}, - boundary_parameters={'hatch':'o', 'linewidth':1.5, 'facecolor':'y', 'edgecolor':'b', 'alpha':0.3}) + field_parameters={'alpha':0.8, 'cmap':'RdBu', 'interpolation':'none', 'post_process':field_func}, + boundary_parameters={'hatch':'o', 'linewidth':1.5, 'facecolor':'y', 'edgecolor':'b', 'alpha':0.3}) plt.show() plt.savefig('sim_domain.png') ``` +If you just want to quickly visualize the simulation domain without the fields (i.e., when +setting up your simulation), there is no need to invoke the `run` function prior to calling +`plot2D`. Just define the `Simulation` object followed by any DFT monitors and then +invoke `plot2D`. -Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects to be called -on all processes, but only generates a plot on the master process. +Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects +to be called on all processes, but only generates a plot on the master process. **Parameters:** @@ -2677,6 +2686,12 @@ on all processes, but only generates a plot on the master process. - `alpha=1.0`: transparency of geometry - `contour=False`: if `True`, plot a contour of the geometry rather than its image - `contour_linewidth=1`: line width of the contour lines if `contour=True` + - `frequency=None`: for materials with a [frequency-dependent + permittivity](Materials.md#material-dispersion) $\varepsilon(f)$, specifies the + frequency $f$ (in Meep units) of the real part of the permittivity to use in the + plot. Defaults to the `frequency` parameter of the [Source](#source) object. + - `resolution=None`: the resolution of the $\varepsilon$ grid. Defaults to the + `resolution` of the `Simulation` object. * `boundary_parameters`: a `dict` of optional plotting parameters that override the default parameters for the boundary layers. - `alpha=1.0`: transparency of boundary layers @@ -2713,10 +2728,6 @@ on all processes, but only generates a plot on the master process. - `alpha=0.6`: transparency of fields - `post_process=np.real`: post processing function to apply to fields (must be a function object) -* `frequency`: for materials with a [frequency-dependent - permittivity](Materials.md#material-dispersion) $\varepsilon(f)$, specifies the - frequency $f$ (in Meep units) of the real part of the permittivity to use in the - plot. Defaults to the `frequency` parameter of the [Source](#source) object.
@@ -4401,7 +4412,7 @@ def __init__(self, medium2, weights=None, grid_type='U_DEFAULT', - do_averaging=False, + do_averaging=True, beta=0, eta=0.5, damping=0): @@ -6893,7 +6904,7 @@ A class used to record the fields during timestepping (i.e., a [`run`](#run-func function). The object is initialized prior to timestepping by specifying the simulation object and the field component. The object can then be passed to any [step-function modifier](#step-function-modifiers). For example, one can record the -Ez fields at every one time unit using: +$E_z$ fields at every one time unit using: ```py animate = mp.Animate2D(sim, diff --git a/doc/docs/Python_User_Interface.md.in b/doc/docs/Python_User_Interface.md.in index 362f652ec..2ed09ca3d 100644 --- a/doc/docs/Python_User_Interface.md.in +++ b/doc/docs/Python_User_Interface.md.in @@ -480,7 +480,7 @@ This feature is only available if Meep is built with [libGDSII](Build_From_Sourc ### Data Visualization -This module provides basic visualization functionality for the simulation domain. The spirit of the module is to provide functions that can be called with *no customization options whatsoever* and will do useful relevant things by default, but which can also be customized in cases where you *do* want to take the time to spruce up the output. The `Simulation` class provides the following methods: +This module provides basic visualization functionality for the simulation domain. The intent of the module is to provide functions that can be called with *no customization options whatsoever* and will do useful relevant things by default, but which can also be customized in cases where you *do* want to take the time to spruce up the output. The `Simulation` class provides the following methods: @@ Simulation.plot2D @@ @@ Simulation.plot3D @@ diff --git a/python/meep.i b/python/meep.i index 948d6cc78..f1dc8012b 100644 --- a/python/meep.i +++ b/python/meep.i @@ -1095,12 +1095,12 @@ void _get_gradient(PyObject *grad, double scalegrad, PyObject *fields_a, PyObjec $1 = (double *)array_data($input); } -%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* grid_vals { +%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") std::complex* grid_vals { $1 = is_array($input); } -%typemap(in, fragment="NumPy_Macros") double* grid_vals { - $1 = (double *)array_data($input); +%typemap(in, fragment="NumPy_Macros") std::complex* grid_vals { + $1 = (std::complex *)array_data($input); } // typemap for solve_cw: @@ -2040,7 +2040,8 @@ void _get_epsilon_grid(geometric_object_list gobj_list, int nx, double *xtics, int ny, double *ytics, int nz, double *ztics, - double *grid_vals) { + std::complex *grid_vals, + double frequency) { meep_geom::get_epsilon_grid(gobj_list, mlist, _default_material, @@ -2051,7 +2052,8 @@ void _get_epsilon_grid(geometric_object_list gobj_list, nx, xtics, ny, ytics, nz, ztics, - grid_vals); + grid_vals, + frequency); } %} diff --git a/python/simulation.py b/python/simulation.py index 83742fc71..5fbf66155 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -2225,18 +2225,19 @@ 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): + def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None, frequency=0): """ 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 + grid anywhere within the cell volume, compute the trace of the $\\varepsilon(f)$ tensor at frequency + $f$ (in Meep units) from the `geometry` exactly at each grid point. `frequency` defaults to 0 which is + the instantaneous $\\varepsilon$. (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)))) + grid_vals = np.squeeze(np.empty((len(xtics), len(ytics), len(ztics)), dtype=np.complex128)) gv = self._create_grid_volume(False) mp._get_epsilon_grid(self.geometry, self.extra_materials, @@ -2247,7 +2248,8 @@ def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None): len(xtics), xtics, len(ytics), ytics, len(ztics), ztics, - grid_vals) + grid_vals, + frequency) return grid_vals def get_filename_prefix(self): @@ -2755,8 +2757,8 @@ def get_farfield(self, near2far, x): (Fourier-transformed) "far" fields at `x` as list of length 6`nfreq`, consisting of fields $(E_x^1,E_y^1,E_z^1,H_x^1,H_y^1,H_z^1,E_x^2,E_y^2,E_z^2,H_x^2,H_y^2,H_z^2,...)$ in Cartesian coordinates and - $(E_r^1,E_\phi^1,E_z^1,H_r^1,H_\phi^1,H_z^1,E_r^2,E_\phi^2,E_z^2,H_r^2,H_\phi^2,H_z^2,...)$ - in cylindrical coordinates for the frequencies 1,2,…,`nfreq`. + $(E_r^1,E_\\phi^1,E_z^1,H_r^1,H_\\phi^1,H_z^1,E_r^2,E_\\phi^2,E_z^2,H_r^2,H_\\phi^2,H_z^2,...)$ + in cylindrical coordinates for the frequencies 1,2,...,`nfreq`. """ return mp._get_farfield(near2far.swigobj, py_v3_to_vec(self.dimensions, x, is_cylindrical=self.is_cylindrical)) @@ -4084,14 +4086,18 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False, field_func = lambda x: 20*np.log10(np.abs(x)) import matplotlib.pyplot as plt sim.plot2D(fields=mp.Ez, - field_parameters={'alpha':0.8, 'cmap':'RdBu', 'interpolation':'none', 'post_process':field_func}, - boundary_parameters={'hatch':'o', 'linewidth':1.5, 'facecolor':'y', 'edgecolor':'b', 'alpha':0.3}) + field_parameters={'alpha':0.8, 'cmap':'RdBu', 'interpolation':'none', 'post_process':field_func}, + boundary_parameters={'hatch':'o', 'linewidth':1.5, 'facecolor':'y', 'edgecolor':'b', 'alpha':0.3}) plt.show() plt.savefig('sim_domain.png') ``` + If you just want to quickly visualize the simulation domain without the fields (i.e., when + setting up your simulation), there is no need to invoke the `run` function prior to calling + `plot2D`. Just define the `Simulation` object followed by any DFT monitors and then + invoke `plot2D`. - Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects to be called - on all processes, but only generates a plot on the master process. + Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects + to be called on all processes, but only generates a plot on the master process. **Parameters:** @@ -4113,6 +4119,12 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False, - `alpha=1.0`: transparency of geometry - `contour=False`: if `True`, plot a contour of the geometry rather than its image - `contour_linewidth=1`: line width of the contour lines if `contour=True` + - `frequency=None`: for materials with a [frequency-dependent + permittivity](Materials.md#material-dispersion) $\\varepsilon(f)$, specifies the + frequency $f$ (in Meep units) of the real part of the permittivity to use in the + plot. Defaults to the `frequency` parameter of the [Source](#source) object. + - `resolution=None`: the resolution of the $\\varepsilon$ grid. Defaults to the + `resolution` of the `Simulation` object. * `boundary_parameters`: a `dict` of optional plotting parameters that override the default parameters for the boundary layers. - `alpha=1.0`: transparency of boundary layers @@ -4149,23 +4161,26 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False, - `alpha=0.6`: transparency of fields - `post_process=np.real`: post processing function to apply to fields (must be a function object) - * `frequency`: for materials with a [frequency-dependent - permittivity](Materials.md#material-dispersion) $\\varepsilon(f)$, specifies the - frequency $f$ (in Meep units) of the real part of the permittivity to use in the - plot. Defaults to the `frequency` parameter of the [Source](#source) object. - """ - return vis.plot2D(self, ax=ax, output_plane=output_plane, fields=fields, labels=labels, - eps_parameters=eps_parameters, boundary_parameters=boundary_parameters, - source_parameters=source_parameters,monitor_parameters=monitor_parameters, - field_parameters=field_parameters, frequency=frequency, - plot_eps_flag=plot_eps_flag, plot_sources_flag=plot_sources_flag, - plot_monitors_flag=plot_monitors_flag, plot_boundaries_flag=plot_boundaries_flag, - **kwargs) - - - def plot_fields(self,**kwargs): - return vis.plot_fields(self,**kwargs) + return vis.plot2D(self, + ax=ax, + output_plane=output_plane, + fields=fields, + labels=labels, + eps_parameters=eps_parameters, + boundary_parameters=boundary_parameters, + source_parameters=source_parameters, + monitor_parameters=monitor_parameters, + field_parameters=field_parameters, + frequency=frequency, + plot_eps_flag=plot_eps_flag, + plot_sources_flag=plot_sources_flag, + plot_monitors_flag=plot_monitors_flag, + plot_boundaries_flag=plot_boundaries_flag, + **kwargs) + + def plot_fields(self, **kwargs): + return vis.plot_fields(self, **kwargs) def plot3D(self): """ diff --git a/python/tests/test_visualization.py b/python/tests/test_visualization.py index e2b9d2846..1b030716e 100644 --- a/python/tests/test_visualization.py +++ b/python/tests/test_visualization.py @@ -161,7 +161,7 @@ def test_animation_output(self): sim = setup_sim() # generate 2D simulation Animate = mp.Animate2D(sim,fields=mp.Ez, realtime=False, normalize=False) # Check without normalization - Animate_norm = mp.Animate2D(sim,mp.Ez,realtime=False,normalize=True) # Check with normalization + Animate_norm = mp.Animate2D(sim, mp.Ez, realtime=False, normalize=True) # Check with normalization # test both animation objects during same run sim.run( diff --git a/python/visualization.py b/python/visualization.py index b42a3dcf7..1e17f6a08 100644 --- a/python/visualization.py +++ b/python/visualization.py @@ -47,7 +47,9 @@ 'cmap':'binary', 'alpha':1.0, 'contour':False, - 'contour_linewidth':1 + 'contour_linewidth':1, + 'frequency':None, + 'resolution':None } default_boundary_parameters = { @@ -95,7 +97,7 @@ def filter_dict(dict_to_filter, func_with_kwargs): # ------------------------------------------------------- # # Routines to add legends to plot -def place_label(ax,label_text,x,y,centerx,centery,label_parameters=None): +def place_label(ax, label_text, x, y, centerx, centery, label_parameters=None): label_parameters = default_label_parameters if label_parameters is None else dict(default_label_parameters, **label_parameters) @@ -112,29 +114,37 @@ def place_label(ax,label_text,x,y,centerx,centery,label_parameters=None): else: ytext = offset - ax.annotate(label_text, xy=(x, y), xytext=(xtext,ytext), - textcoords='offset points', ha='center', va='bottom', - bbox=dict(boxstyle='round,pad=0.2', fc=color, alpha=alpha), - arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0.5', - color=color)) + ax.annotate(label_text, xy=(x, y), xytext=(xtext, ytext), + textcoords='offset points', ha='center', va='bottom', + bbox=dict(boxstyle='round,pad=0.2', fc=color, alpha=alpha), + arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0.5', + color=color)) return ax # ------------------------------------------------------- # # Helper functions used to plot volumes on a 2D plane -# Returns the intersection points of 2 Volumes. +# Returns the intersection points of two Volumes. # Volumes must be a line, plane, or rectangular prism # (since they are volume objects) -def intersect_volume_volume(volume1,volume2): +def intersect_volume_volume(volume1, volume2): # volume1 ............... [volume] # volume2 ............... [volume] # Represent the volumes by an "upper" and "lower" coordinate - U1 = [volume1.center.x+volume1.size.x/2,volume1.center.y+volume1.size.y/2,volume1.center.z+volume1.size.z/2] - L1 = [volume1.center.x-volume1.size.x/2,volume1.center.y-volume1.size.y/2,volume1.center.z-volume1.size.z/2] - - U2 = [volume2.center.x+volume2.size.x/2,volume2.center.y+volume2.size.y/2,volume2.center.z+volume2.size.z/2] - L2 = [volume2.center.x-volume2.size.x/2,volume2.center.y-volume2.size.y/2,volume2.center.z-volume2.size.z/2] + U1 = [volume1.center.x+volume1.size.x/2, + volume1.center.y+volume1.size.y/2, + volume1.center.z+volume1.size.z/2] + L1 = [volume1.center.x-volume1.size.x/2, + volume1.center.y-volume1.size.y/2, + volume1.center.z-volume1.size.z/2] + + U2 = [volume2.center.x+volume2.size.x/2, + volume2.center.y+volume2.size.y/2, + volume2.center.z+volume2.size.z/2] + L2 = [volume2.center.x-volume2.size.x/2, + volume2.center.y-volume2.size.y/2, + volume2.center.z-volume2.size.z/2] # Evaluate intersection U = np.min([U1,U2],axis=0) @@ -164,13 +174,13 @@ def intersect_volume_volume(volume1,volume2): # All of the 2D plotting routines need an output plane over which to plot. # The user has many options to specify this output plane. They can pass # the output_plane parameter, which is a 2D volume object. They can specify -# a volume using in_volume, which stores the volume as a C volume, not a python +# a volume using in_volume, which stores the volume as a C volume, not a Python # volume. They can also do nothing and plot the XY plane through Z=0. # # Not only do we need to check for all of these possibilities, but we also need # to check if the user accidentally specifies a plane that stretches beyond the # simulation domain. -def get_2D_dimensions(sim,output_plane): +def get_2D_dimensions(sim, output_plane): from meep.simulation import Volume # Pull correct plane from user @@ -182,12 +192,12 @@ def get_2D_dimensions(sim,output_plane): plane_center, plane_size = (sim.geometry_center, sim.cell_size) plane_volume = Volume(center=plane_center,size=plane_size) - if plane_size.x!=0 and plane_size.y!=0 and plane_size.z!=0: + if plane_size.x != 0 and plane_size.y != 0 and plane_size.z != 0: raise ValueError("Plane volume must be 2D (a plane).") - check_volume = Volume(center=sim.geometry_center,size=sim.cell_size) + check_volume = Volume(center=sim.geometry_center, size=sim.cell_size) - vertices = intersect_volume_volume(check_volume,plane_volume) + vertices = intersect_volume_volume(check_volume, plane_volume) if len(vertices) == 0: raise ValueError("The specified user volume is completely outside of the simulation domain.") @@ -203,10 +213,7 @@ def get_2D_dimensions(sim,output_plane): # ------------------------------------------------------- # # actual plotting routines -def plot_volume(sim,ax,volume,output_plane=None,plotting_parameters=None,label=None): - if not sim._is_initialized: - sim.init_sim() - +def plot_volume(sim, ax, volume, output_plane=None, plotting_parameters=None, label=None): import matplotlib.patches as patches from matplotlib import pyplot as plt from meep.simulation import Volume @@ -215,39 +222,57 @@ def plot_volume(sim,ax,volume,output_plane=None,plotting_parameters=None,label=N plotting_parameters = default_volume_parameters if plotting_parameters is None else dict(default_volume_parameters, **plotting_parameters) # Get domain measurements - sim_center, sim_size = get_2D_dimensions(sim,output_plane) + sim_center, sim_size = get_2D_dimensions(sim, output_plane) - plane = Volume(center=sim_center,size=sim_size) + plane = Volume(center=sim_center, size=sim_size) # Pull volume parameters size = volume.size center = volume.center - xmax = center.x+size.x/2 - xmin = center.x-size.x/2 - ymax = center.y+size.y/2 - ymin = center.y-size.y/2 - zmax = center.z+size.z/2 - zmin = center.z-size.z/2 + xmax = center.x + size.x/2 + xmin = center.x - size.x/2 + ymax = center.y + size.y/2 + ymin = center.y - size.y/2 + zmax = center.z + size.z/2 + zmin = center.z - size.z/2 # Add labels if requested if label is not None and mp.am_master(): if sim_size.x == 0: - ax = place_label(ax,label,center.y,center.z,sim_center.y,sim_center.z,label_parameters=plotting_parameters) + ax = place_label(ax, + label, + center.y, + center.z, + sim_center.y, + sim_center.z, + label_parameters=plotting_parameters) elif sim_size.y == 0: - ax = place_label(ax,label,center.x,center.z,sim_center.x,sim_center.z,label_parameters=plotting_parameters) + ax = place_label(ax, + label, + center.x, + center.z, + sim_center.x, + sim_center.z, + label_parameters=plotting_parameters) elif sim_size.z == 0: - ax = place_label(ax,label,center.x,center.y,sim_center.x,sim_center.y,label_parameters=plotting_parameters) + ax = place_label(ax, + label, + center.x, + center.y, + sim_center.x, + sim_center.y, + label_parameters=plotting_parameters) # Intersect plane with volume - intersection = intersect_volume_volume(volume,plane) + intersection = intersect_volume_volume(volume, plane) # Sort the points in a counter clockwise manner to ensure convex polygon is formed def sort_points(xy): xy = np.squeeze(xy) - xy_mean = np.mean(xy,axis=0) - theta = np.arctan2(xy[:,1]-xy_mean[1],xy[:,0]-xy_mean[0]) - return xy[np.argsort(theta,axis=0),:] + xy_mean = np.mean(xy, axis=0) + theta = np.arctan2(xy[:,1] - xy_mean[1], xy[:,0] - xy_mean[0]) + return xy[np.argsort(theta, axis=0), :] if mp.am_master(): # Point volume @@ -269,16 +294,16 @@ def sort_points(xy): elif len(intersection) == 2: line_args = {key:value for key, value in plotting_parameters.items() if key in ['color','linestyle','linewidth','alpha']} # Plot YZ - if sim_size.x==0: - ax.plot([a.y for a in intersection],[a.z for a in intersection], **line_args) + if sim_size.x == 0: + ax.plot([a.y for a in intersection], [a.z for a in intersection], **line_args) return ax - #Plot XZ - elif sim_size.y==0: - ax.plot([a.x for a in intersection],[a.z for a in intersection], **line_args) + # Plot XZ + elif sim_size.y == 0: + ax.plot([a.x for a in intersection], [a.z for a in intersection], **line_args) return ax # Plot XY - elif sim_size.z==0: - ax.plot([a.x for a in intersection],[a.y for a in intersection], **line_args) + elif sim_size.z == 0: + ax.plot([a.x for a in intersection], [a.y for a in intersection], **line_args) return ax else: return ax @@ -287,15 +312,15 @@ def sort_points(xy): elif len(intersection) > 2: planar_args = {key:value for key, value in plotting_parameters.items() if key in ['edgecolor','linewidth','facecolor','hatch','alpha']} # Plot YZ - if sim_size.x==0: + if sim_size.x == 0: ax.add_patch(patches.Polygon(sort_points([[a.y,a.z] for a in intersection]), **planar_args)) return ax - #Plot XZ + # Plot XZ elif sim_size.y==0: ax.add_patch(patches.Polygon(sort_points([[a.x,a.z] for a in intersection]), **planar_args)) return ax # Plot XY - elif sim_size.z==0: + elif sim_size.z == 0: ax.add_patch(patches.Polygon(sort_points([[a.x,a.y] for a in intersection]), **planar_args)) return ax else: @@ -304,16 +329,32 @@ def sort_points(xy): return ax return ax -def plot_eps(sim,ax,output_plane=None,eps_parameters=None,frequency=0): - if sim.structure is None: - sim.init_sim() - - +def plot_eps(sim, ax, output_plane=None, eps_parameters=None, frequency=None): # consolidate plotting parameters eps_parameters = default_eps_parameters if eps_parameters is None else dict(default_eps_parameters, **eps_parameters) + # Determine a frequency to plot all epsilon + if frequency is not None: + warnings.warn('The frequency parameter of plot2D has been deprecated. Use the frequency key of the eps_parameters dictionary instead.') + eps_parameters['frequency'] = frequency + if eps_parameters['frequency'] is None: + try: + # Continuous sources + eps_parameters['frequency'] = sim.sources[0].frequency + except: + try: + # Gaussian sources + eps_parameters['frequency'] = sim.sources[0].src.frequency + except: + try: + # Custom sources + eps_parameters['frequency'] = sim.sources[0].src.center_frequency + except: + # No sources + eps_parameters['frequency'] = 0 + # Get domain measurements - sim_center, sim_size = get_2D_dimensions(sim,output_plane) + sim_center, sim_size = get_2D_dimensions(sim, output_plane) xmin = sim_center.x - sim_size.x/2 xmax = sim_center.x + sim_size.x/2 @@ -322,28 +363,43 @@ def plot_eps(sim,ax,output_plane=None,eps_parameters=None,frequency=0): zmin = sim_center.z - sim_size.z/2 zmax = sim_center.z + sim_size.z/2 - center = Vector3(sim_center.x,sim_center.y,sim_center.z) - cell_size = Vector3(sim_size.x,sim_size.y,sim_size.z) + center = Vector3(sim_center.x, sim_center.y, sim_center.z) + cell_size = Vector3(sim_size.x, sim_size.y, sim_size.z) + + grid_resolution = eps_parameters['resolution'] if eps_parameters['resolution'] else sim.resolution + Nx = int((xmax - xmin) * grid_resolution + 1) + Ny = int((ymax - ymin) * grid_resolution + 1) + Nz = int((zmax - zmin) * grid_resolution + 1) if sim_size.x == 0: # Plot y on x axis, z on y axis (YZ plane) - extent = [ymin,ymax,zmin,zmax] + extent = [ymin, ymax, zmin, zmax] xlabel = 'Y' ylabel = 'Z' + xtics = np.array([sim_center.x]) + ytics = np.linspace(ymin, ymax, Ny) + ztics = np.linspace(zmin, zmax, Nz) elif sim_size.y == 0: # Plot x on x axis, z on y axis (XZ plane) - extent = [xmin,xmax,zmin,zmax] + extent = [xmin, xmax, zmin, zmax] xlabel = 'X' ylabel = 'Z' + xtics = np.linspace(xmin, xmax, Nx) + ytics = np.array([sim_center.y]) + ztics = np.linspace(zmin, zmax, Nz) elif sim_size.z == 0: # Plot x on x axis, y on y axis (XY plane) - extent = [xmin,xmax,ymin,ymax] + extent = [xmin, xmax, ymin, ymax] xlabel = 'X' ylabel = 'Y' + xtics = np.linspace(xmin, xmax, Nx) + ytics = np.linspace(ymin, ymax, Ny) + ztics = np.array([sim_center.z]) else: raise ValueError("A 2D plane has not been specified...") - eps_data = np.rot90(np.real(sim.get_array(center=center, size=cell_size, component=mp.Dielectric, frequency=frequency))) + eps_data = np.rot90(np.real(sim.get_epsilon_grid(xtics, ytics, ztics, eps_parameters['frequency']))) + if mp.am_master(): if eps_parameters['contour']: ax.contour(eps_data, 0, colors='black', origin='upper', extent=extent, linewidths=eps_parameters['contour_linewidth']) @@ -354,14 +410,11 @@ def plot_eps(sim,ax,output_plane=None,eps_parameters=None,frequency=0): return ax -def plot_boundaries(sim,ax,output_plane=None,boundary_parameters=None): - if not sim._is_initialized: - sim.init_sim() - +def plot_boundaries(sim, ax, output_plane=None, boundary_parameters=None): # consolidate plotting parameters boundary_parameters = default_boundary_parameters if boundary_parameters is None else dict(default_boundary_parameters, **boundary_parameters) - def get_boundary_volumes(thickness,direction,side): + def get_boundary_volumes(thickness, direction, side): from meep.simulation import Volume thickness = boundary.thickness @@ -403,20 +456,20 @@ def get_boundary_volumes(thickness,direction,side): import itertools for boundary in sim.boundary_layers: - # All 4 side are the same + # All four sides are the same if boundary.direction == mp.ALL and boundary.side == mp.ALL: if sim.dimensions == 1: dims = [mp.X] elif sim.dimensions == 2: - dims = [mp.X,mp.Y] + dims = [mp.X, mp.Y] elif sim.dimensions == 3: - dims = [mp.X,mp.Y,mp.Z] + dims = [mp.X, mp.Y, mp.Z] else: raise ValueError("Invalid simulation dimensions") for permutation in itertools.product(dims, [mp.Low, mp.High]): vol = get_boundary_volumes(boundary.thickness,*permutation) ax = plot_volume(sim,ax,vol,output_plane,plotting_parameters=boundary_parameters) - # 2 sides are the same + # two sides are the same elif boundary.side == mp.ALL: for side in [mp.Low, mp.High]: vol = get_boundary_volumes(boundary.thickness,boundary.direction,side) @@ -427,10 +480,7 @@ def get_boundary_volumes(thickness,direction,side): ax = plot_volume(sim,ax,vol,output_plane,plotting_parameters=boundary_parameters) return ax -def plot_sources(sim,ax,output_plane=None,labels=False,source_parameters=None): - if not sim._is_initialized: - sim.init_sim() - +def plot_sources(sim, ax, output_plane=None, labels=False, source_parameters=None): from meep.simulation import Volume # consolidate plotting parameters @@ -443,13 +493,10 @@ def plot_sources(sim,ax,output_plane=None,labels=False,source_parameters=None): ax = plot_volume(sim,ax,vol,output_plane,plotting_parameters=source_parameters,label=label) return ax -def plot_monitors(sim,ax,output_plane=None,labels=False,monitor_parameters=None): - if not sim._is_initialized: - sim.init_sim() - +def plot_monitors(sim, ax, output_plane=None, labels=False, monitor_parameters=None): from meep.simulation import Volume - # consolidate plotting parameters + # consolidate plotting parameters monitor_parameters = default_monitor_parameters if monitor_parameters is None else dict(default_monitor_parameters, **monitor_parameters) label = 'monitor' if labels else None @@ -460,7 +507,7 @@ def plot_monitors(sim,ax,output_plane=None,labels=False,monitor_parameters=None) ax = plot_volume(sim,ax,vol,output_plane,plotting_parameters=monitor_parameters,label=label) return ax -def plot_fields(sim,ax=None,fields=None,output_plane=None,field_parameters=None): +def plot_fields(sim, ax=None, fields=None, output_plane=None, field_parameters=None): if not sim._is_initialized: sim.init_sim() @@ -472,7 +519,7 @@ def plot_fields(sim,ax=None,fields=None,output_plane=None,field_parameters=None) # user specifies a field component if fields in [mp.Ex, mp.Ey, mp.Ez, mp.Hx, mp.Hy, mp.Hz]: # Get domain measurements - sim_center, sim_size = get_2D_dimensions(sim,output_plane) + sim_center, sim_size = get_2D_dimensions(sim, output_plane) xmin = sim_center.x - sim_size.x/2 xmax = sim_center.x + sim_size.x/2 @@ -481,22 +528,22 @@ def plot_fields(sim,ax=None,fields=None,output_plane=None,field_parameters=None) zmin = sim_center.z - sim_size.z/2 zmax = sim_center.z + sim_size.z/2 - center = Vector3(sim_center.x,sim_center.y,sim_center.z) - cell_size = Vector3(sim_size.x,sim_size.y,sim_size.z) + center = Vector3(sim_center.x, sim_center.y, sim_center.z) + cell_size = Vector3(sim_size.x, sim_size.y, sim_size.z) if sim_size.x == 0: # Plot y on x axis, z on y axis (YZ plane) - extent = [ymin,ymax,zmin,zmax] + extent = [ymin, ymax, zmin, zmax] xlabel = 'Y' ylabel = 'Z' elif sim_size.y == 0: # Plot x on x axis, z on y axis (XZ plane) - extent = [xmin,xmax,zmin,zmax] + extent = [xmin, xmax, zmin, zmax] xlabel = 'X' ylabel = 'Z' elif sim_size.z == 0: # Plot x on x axis, y on y axis (XY plane) - extent = [xmin,xmax,ymin,ymax] + extent = [xmin, xmax, ymin, ymax] xlabel = 'X' ylabel = 'Y' fields = sim.get_array(center=center, size=cell_size, component=fields) @@ -515,73 +562,72 @@ def plot_fields(sim,ax=None,fields=None,output_plane=None,field_parameters=None) return np.rot90(fields) return ax -def plot2D(sim,ax=None, output_plane=None, fields=None, labels=False, - eps_parameters=None,boundary_parameters=None, - source_parameters=None,monitor_parameters=None, +def plot2D(sim, ax=None, output_plane=None, fields=None, labels=False, + eps_parameters=None, boundary_parameters=None, + source_parameters=None, monitor_parameters=None, field_parameters=None, frequency=None, plot_eps_flag=True, plot_sources_flag=True, plot_monitors_flag=True, plot_boundaries_flag=True): - # Initialize the simulation - if sim.structure is None: - sim.init_sim() # Ensure a figure axis exists if ax is None and mp.am_master(): from matplotlib import pyplot as plt ax = plt.gca() - # Determine a frequency to plot all epsilon - if frequency is None: - try: - # Continuous sources - frequency = sim.sources[0].frequency - except: - try: - # Gaussian sources - frequency = sim.sources[0].src.frequency - except: - try: - # Custom sources - frequency = sim.sources[0].src.center_frequency - except: - # No sources - frequency = 0 # validate the output plane to ensure proper 2D coordinates from meep.simulation import Volume - sim_center, sim_size = get_2D_dimensions(sim,output_plane) - output_plane = Volume(center=sim_center,size=sim_size) + sim_center, sim_size = get_2D_dimensions(sim, output_plane) + output_plane = Volume(center=sim_center, size=sim_size) # Plot geometry if plot_eps_flag: - ax = plot_eps(sim,ax,output_plane=output_plane,eps_parameters=eps_parameters,frequency=frequency) + ax = plot_eps(sim, ax, output_plane=output_plane, + eps_parameters=eps_parameters, frequency=frequency) # Plot boundaries if plot_boundaries_flag: - ax = plot_boundaries(sim,ax,output_plane=output_plane,boundary_parameters=boundary_parameters) + ax = plot_boundaries(sim, ax, output_plane=output_plane, + boundary_parameters=boundary_parameters) # Plot sources if plot_sources_flag: - ax = plot_sources(sim,ax,output_plane=output_plane,labels=labels,source_parameters=source_parameters) + ax = plot_sources(sim, ax, output_plane=output_plane, + labels=labels, source_parameters=source_parameters) # Plot monitors if plot_monitors_flag: - ax = plot_monitors(sim,ax,output_plane=output_plane,labels=labels,monitor_parameters=monitor_parameters) + ax = plot_monitors(sim, ax, output_plane=output_plane, + labels=labels, monitor_parameters=monitor_parameters) # Plot fields - ax = plot_fields(sim,ax,fields,output_plane=output_plane,field_parameters=field_parameters) + if fields: + ax = plot_fields(sim, ax, fields, output_plane=output_plane, + field_parameters=field_parameters) return ax def plot3D(sim): from mayavi import mlab - if not sim._is_initialized: - sim.init_sim() - if sim.dimensions < 3: raise ValueError("Simulation must have 3 dimensions to visualize 3D") - eps_data = sim.get_epsilon() + xmin = sim.geometry_center.x - 0.5*sim.cell_size.x + xmax = sim.geometry_center.x + 0.5*sim.cell_size.x + ymin = sim.geometry_center.y - 0.5*sim.cell_size.y + ymax = sim.geometry_center.y + 0.5*sim.cell_size.y + zmin = sim.geometry_center.z - 0.5*sim.cell_size.z + zmax = sim.geometry_center.z + 0.5*sim.cell_size.z + + Nx = int(sim.cell_size.x * sim.resolution) + 1 + Ny = int(sim.cell_size.y * sim.resolution) + 1 + Nz = int(sim.cell_size.z * sim.resolution) + 1 + + xtics = np.linspace(xmin, xmax, Nx) + ytics = np.linspace(ymin, ymax, Ny) + ztics = np.linspace(zmin, zmax, Nz) + + eps_data = sim.get_epsilon_grid(xtics, ytics, ztics) s = mlab.contour3d(eps_data, colormap="YlGnBu") return s @@ -711,7 +757,7 @@ class Animate2D(object): function). The object is initialized prior to timestepping by specifying the simulation object and the field component. The object can then be passed to any [step-function modifier](#step-function-modifiers). For example, one can record the - Ez fields at every one time unit using: + $E_z$ fields at every one time unit using: ```py animate = mp.Animate2D(sim, @@ -792,7 +838,6 @@ def mod1(ax): self.plot_modifiers = plot_modifiers self.customization_args = customization_args - self.cumulative_fields = [] self._saved_frames = [] @@ -813,7 +858,7 @@ def __call__(self,sim,todo): # Initialize the plot if not self.init: filtered_plot2D = filter_dict(self.customization_args, plot2D) - ax = sim.plot2D(ax=self.ax,fields=self.fields,**filtered_plot2D) + ax = sim.plot2D(ax=self.ax, fields=self.fields, **filtered_plot2D) # Run the plot modifier functions if self.plot_modifiers: for k in range(len(self.plot_modifiers)): @@ -827,7 +872,7 @@ def __call__(self,sim,todo): else: # Update the plot filtered_plot_fields= filter_dict(self.customization_args, plot_fields) - fields = sim.plot_fields(fields=self.fields,**filtered_plot_fields) + fields = sim.plot_fields(fields=self.fields, **filtered_plot_fields) if mp.am_master(): self.ax.images[-1].set_data(fields) self.ax.images[-1].set_clim(vmin=0.8*np.min(fields), vmax=0.8*np.max(fields)) @@ -851,7 +896,7 @@ def __call__(self,sim,todo): if self.normalize and mp.am_master(): if mp.verbosity.meep > 0: print("Normalizing field data...") - fields = np.array(self.cumulative_fields) / np.max(np.abs(self.cumulative_fields),axis=(0,1,2)) + fields = np.array(self.cumulative_fields) / np.max(np.abs(self.cumulative_fields), axis=(0,1,2)) for k in range(len(self.cumulative_fields)): self.ax.images[-1].set_data(fields[k,:,:]) self.ax.images[-1].set_clim(vmin=-0.8, vmax=0.8) diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index dda2b5782..a5fdfb56f 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -2557,13 +2557,11 @@ void invert_tensor(std::complex t_inv[9], std::complex t[9]) { #undef minv(x,y) } -void eff_chi1inv_row_disp(meep::component c, std::complex chi1inv_row[3], - const meep::vec &r, double freq, geom_epsilon *geps) { +void get_chi1_tensor_disp(std::complex tensor[9], const meep::vec &r, double freq, geom_epsilon *geps) { // locate the proper material material_type md; geps->get_material_pt(md, r); const medium_struct *mm = &(md->medium); - std::complex tensor[9], tensor_inv[9]; // loop over all the tensor components for (int i = 0; i < 9; i++) { @@ -2574,7 +2572,7 @@ void eff_chi1inv_row_disp(meep::component c, std::complex chi1inv_row[3] double conductivityCur = vec_to_value(mm->D_conductivity_diag, dummy, i); a = std::complex(1.0, conductivityCur / (2*meep::pi*freq)); - // compute lorentzian component + // compute lorentzian component including the instantaneous ε b = cvec_to_value(mm->epsilon_diag, mm->epsilon_offdiag, i); for (const auto &mm_susc: mm->E_susceptibilities) { meep::lorentzian_susceptibility sus = meep::lorentzian_susceptibility( @@ -2586,7 +2584,12 @@ void eff_chi1inv_row_disp(meep::component c, std::complex chi1inv_row[3] // elementwise multiply tensor[i] = a * b; } +} +void eff_chi1inv_row_disp(meep::component c, std::complex chi1inv_row[3], + const meep::vec &r, double freq, geom_epsilon *geps) { + std::complex tensor[9], tensor_inv[9]; + get_chi1_tensor_disp(tensor, r, freq, geps); // invert the matrix invert_tensor(tensor_inv, tensor); @@ -3050,7 +3053,8 @@ void get_epsilon_grid(geometric_object_list gobj_list, int nx, const double *x, int ny, const double *y, int nz, const double *z, - double *grid_vals) { + std::complex *grid_vals, + double frequency) { double min_val[3], max_val[3]; for (int n = 0; n < 3; ++n) { int ndir = (n == 0) ? nx : ((n == 1) ? ny : nz); @@ -3066,10 +3070,17 @@ void get_epsilon_grid(geometric_object_list gobj_list, 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 + for (int k = 0; k < nz; ++k) { + /* obtain the trace of the ε tensor (dispersive or non) 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])); + if (frequency == 0) + grid_vals[k + nz*(j + ny*i)] = geps.chi1p1(meep::E_stuff, meep::vec(x[i],y[j],z[k])); + else { + std::complex tensor[9]; + get_chi1_tensor_disp(tensor, meep::vec(x[i],y[j],z[k]), frequency, &geps); + grid_vals[k + nz*(j + ny*i)] = (tensor[0] + tensor[4] + tensor[8]) / 3.0; + } + } } } // namespace meep_geom diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index d052d5b4a..f66c33e3b 100644 --- a/src/meepgeom.hpp +++ b/src/meepgeom.hpp @@ -281,7 +281,8 @@ void get_epsilon_grid(geometric_object_list gobj_list, int nx, const double *x, int ny, const double *y, int nz, const double *z, - double *grid_vals); + std::complex *grid_vals, + double frequency = 0); 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);