From fff80f27983d922d26da5c05d06b45b611350aeb Mon Sep 17 00:00:00 2001 From: mdecea Date: Fri, 10 Mar 2023 16:56:11 -0800 Subject: [PATCH 1/6] Adding functionality to select modes based on position or overlap --- femwell/mode_solver.py | 72 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 71 insertions(+), 1 deletion(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 48376810..8c6d1e5f 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -326,6 +326,50 @@ def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", directi return fig, axs +def select_mode_by_overlap(mode_basis, E_modes, H_modes, + bbox = None, + elements_list = None, + selection_basis = None): + """ + Selects the mode in the "modes" list that has the highest power contained + within the given bounding box, element or basis. + + If bbox is not None, we select the overlap with the given bounding box + [x_min, x_max, y_min, y_max] + + If elements_list is not None, then the overlap is calculated between the + given modes and the elements in the list + + If elements_list is None, then we calculate the overlap with the given + selection_basis. + """ + + overlaps = list() + + if bbox is not None: + + def sel_fun(x): + return (x[0] < bbox[1]) * (x[0] > bbox[0]) * (x[1] > bbox[2]) * (x[1] < bbox[3]) + + selection_basis = Basis(mode_basis.mesh, mode_basis.elem, + elements = lambda x: sel_fun(x) + ) + + elif elements_list is not None: + + selection_basis = Basis(mode_basis.mesh, mode_basis.elem, elements=elements_list) + + for E_f, H_f in zip(E_modes, H_modes): + overlaps.append(calculate_overlap(selection_basis, E_f, H_f, selection_basis, E_f, H_f)) + + ind_max = np.argmax(np.abs(overlaps)) + + print(overlaps) + + return ind_max + + + if __name__ == "__main__": from collections import OrderedDict @@ -384,7 +428,8 @@ def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", directi # basis0.plot(epsilon, colorbar=True).show() lams, basis, xs = compute_modes( - basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=6, order=2, radius=3 + basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=6, order=2, radius=3, + solver='scipy' ) print(lams) @@ -403,6 +448,8 @@ def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", directi plt.show() integrals = np.zeros((len(lams),) * 2, dtype=complex) + H_modes = list() + for i in range(len(lams)): for j in range(len(lams)): E_i = xs[i] @@ -420,7 +467,30 @@ def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", directi omega=2 * np.pi / 1.55 * scipy.constants.speed_of_light, ) integrals[i, j] = calculate_overlap(basis, E_i, H_i, basis, E_j, H_j) + H_modes.append(H_i) plt.imshow(np.real(integrals)) plt.colorbar() plt.show() + + # Create basis to select a certain simulation extent + def sel_fun(x): + print(x) + return (x[0] < 0) * (x[0] > -1) * (x[1] > 0) * (x[1] < 0.5) + selection_basis = Basis(basis.mesh, basis.elem, + #elements = lambda x: x[0] < 0 and x[0] > -1 and x[1] > 0 and x[1] < 0.5 + elements = lambda x: sel_fun(x) + ) + + print(select_mode_by_overlap(mode_basis = basis, + E_modes = xs, + H_modes = H_modes, + elements_list = ["core"], + selection_basis = None)) + + print(select_mode_by_overlap(mode_basis = basis, + E_modes = xs, + H_modes = H_modes, + elements_list = None, + selection_basis = selection_basis)) + From b676334f7f23ca9cfa54d141077f2bbcfad719e6 Mon Sep 17 00:00:00 2001 From: mdecea Date: Mon, 13 Mar 2023 09:16:34 -0700 Subject: [PATCH 2/6] Adding docs for select_mode --- docs/photonics/examples/selecting_modes.py | 275 +++++++++++++++++++++ femwell/mode_solver.py | 2 +- 2 files changed, 276 insertions(+), 1 deletion(-) create mode 100644 docs/photonics/examples/selecting_modes.py diff --git a/docs/photonics/examples/selecting_modes.py b/docs/photonics/examples/selecting_modes.py new file mode 100644 index 00000000..ca9fb7ec --- /dev/null +++ b/docs/photonics/examples/selecting_modes.py @@ -0,0 +1,275 @@ +# --- +# jupyter: +# jupytext: +# custom_cell_magics: kql +# formats: py:percent,ipynb +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.11.2 +# kernelspec: +# display_name: env_3.11 +# language: python +# name: python3 +# --- + +# %% [markdown] +# # Selecting modes in the mode solver + +# Sometimes we have structures where the mode of interest is +# not the mode with the highest effective index. There are a few +# ways to select modes of interest in femwell + +# %% tags=["hide-input"] +from collections import OrderedDict + +import matplotlib.pyplot as plt +import numpy as np +from scipy.constants import epsilon_0, speed_of_light +from scipy.integrate import solve_ivp +from shapely.geometry import Polygon +from skfem import Basis, ElementTriP0, Mesh +from skfem.io import from_meshio +import scipy.constants + +from femwell.mesh import mesh_from_OrderedDict +from femwell.mode_solver import ( + calculate_coupling_coefficient, + calculate_hfield, + calculate_overlap, + compute_modes, + plot_mode, + select_mode_by_overlap +) + +# %% [markdown] + +# We will use as an example a system with a Si and a SiN sections. +# This could happen, for example, in a system where we are trying +# to heat a SiN waveguide with a Si resistor + + +# %% tags=["remove-stderr", "hide-input"] + +w_sim = 6 +h_clad = 2 +h_box = 2 +w_sin = 1 +w_si = 0.4 +gap = 1.0 +h_sin = 0.4 +h_si = 0.22 + +wavelength = 1.55 +k0 = 2 * np.pi / wavelength + +polygons = OrderedDict( + sin=Polygon( + [ + (-w_sin - gap / 2, 0), + (-w_sin - gap / 2, h_sin), + (-gap / 2, h_sin), + (-gap / 2, 0), + ] + ), + si=Polygon( + [ + (w_si + gap / 2, 0), + (w_si + gap / 2, h_si), + (gap / 2, h_si), + (gap / 2, 0), + ] + ), + clad=Polygon( + [ + (-w_sim / 2, 0), + (-w_sim / 2, h_clad), + (w_sim / 2, h_clad), + (w_sim / 2, 0), + ] + ), + box=Polygon( + [ + (-w_sim / 2, 0), + (-w_sim / 2, -h_box), + (w_sim / 2, -h_box), + (w_sim / 2, 0), + ] + ), +) + +resolutions = dict( + sin={"resolution": 0.03, "distance": 1}, + si={"resolution": 0.03, "distance": 1}, +) + +mesh = from_meshio( + mesh_from_OrderedDict(polygons, resolutions, filename="mesh.msh", default_resolution_max=0.2) +) +mesh.draw().show() + +# %% + +basis0 = Basis(mesh, ElementTriP0(), intorder=4) + +epsilon = basis0.zeros() + 1.444**2 +epsilon[basis0.get_dofs(elements=("si"))] = 3.4777**2 +epsilon[basis0.get_dofs(elements=("sin"))] = 1.973**2 + + +# %% [markdown] + +# ## 0. Directly using femwell + +# If we use `find_modes`, these are the modes we get: + +# %% + +# basis0.plot(epsilon, colorbar=True).show() +lams, basis, xs = compute_modes( + basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=4 +) + + +plot_mode(basis, np.real(xs[0]), direction="x") +plt.show() +plot_mode(basis, np.real(xs[1]), direction="x") +plt.show() +plot_mode(basis, np.real(xs[2]), direction="x") +plt.show() +plot_mode(basis, np.real(xs[3]), direction="x") +plt.show() + +print(f'The effective index of the SiN mode is {np.real(lams[2])}') + +# %% [markdown] + +# We can see how to get the SiN mode (which is the mode of +# interest for us) we need to go to the third mode found by femwell. + +# Are there easier ways to get the SiN modes? Yes! + +# %% [markdown] + +# ## 1. Hack (not 100% accurate): Erasing the Si waveguide + +# One thing we can do to find the SiN mode is to "erase" the Si +# waveguide, or in other words assign the refractive index of SiO2 +# to the Si waveguide. + +# Of course, this is in general not desired, because this way we are +# missing the effect of the presence of the Si waveguide. + +# thi smight not be an issue in this example but there's many +# examples where this is not an acceptable option. + +# %% + +epsilon = basis0.zeros() + 1.444**2 +epsilon[basis0.get_dofs(elements=("si"))] = 1.444**2 +epsilon[basis0.get_dofs(elements=("sin"))] = 1.973**2 + +lams, basis, xs = compute_modes( + basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2 +) + + +plot_mode(basis, np.real(xs[0]), direction="x") +plt.show() +plot_mode(basis, np.real(xs[1]), direction="x") +plt.show() + +print(f'The effective index of the SiN mode is {np.real(lams[0])}') + +# %% [markdown] +# ## 2. Giving a guess effective index + +# We can use the `n_guess` parameter to `compute_modes` to +# select modes close to that effective index. + +# This is great, but of course we need to know what's that guess +# effective index. The way to do that would be to use option 1 above +# and then use that as the n_guess. + +# %% +epsilon = basis0.zeros() + 1.444**2 +epsilon[basis0.get_dofs(elements=("si"))] = 3.4777**2 +epsilon[basis0.get_dofs(elements=("sin"))] = 1.973**2 + +lams, basis, xs = compute_modes( + basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2, + n_guess = 1.62 +) + + +plot_mode(basis, np.real(xs[0]), direction="x") +plt.show() +plot_mode(basis, np.real(xs[1]), direction="x") +plt.show() + +print(f'The effective index of the SiN mode is {np.real(lams[1])}') + +# %% [markdown] + +# You can see how using `n_guess` can still give the wrong mode! + +# %% [markdown] + +# ## 3. Using `select_mode_by_overlap` + +# This allows to choose a mode that has the biggest overlap with +# a given structure. + +# There are two main ways to specify the structure: +# 1. Using the name of the polygon of interest +# 2. Giving a square bounding box of coordinates + +# You can also give it directly the selection_basis of the +# are of interest. + +# A requirement for using `select_mode_by_overlap` is to +# calculate the H field of the found modes. + +# %% + +lams, basis, xs = compute_modes( + basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=4 +) + +# Calculate H field +H_modes = list() +for i, E in enumerate(xs): + H_modes.append( + calculate_hfield( + basis, + E, + lams[i] * (2 * np.pi / wavelength), + omega=2 * np.pi / 1.55 * + scipy.constants.speed_of_light, + )) + +# Option 1: using an element name + +ind_mode = select_mode_by_overlap(basis, xs, H_modes, + bbox = None, + elements_list = "sin", + selection_basis = None) + +plot_mode(basis, np.real(xs[ind_mode]), direction="x") +plt.show() + +# Option 1: using bounding box + +# Format: [xmin, xmax, ymin, ymax] +bbox = [-2, 0, 0, 0.4] + +ind_mode = select_mode_by_overlap(basis, xs, H_modes, + bbox = bbox, + elements_list = None, + selection_basis = None) + +plot_mode(basis, np.real(xs[ind_mode]), direction="x") +plt.show() + +# %% diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 8c6d1e5f..fd223756 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -364,7 +364,7 @@ def sel_fun(x): ind_max = np.argmax(np.abs(overlaps)) - print(overlaps) + #print(overlaps) return ind_max From 57f9a8f4e038859cd77f2d0a87f1742e4340d5fe Mon Sep 17 00:00:00 2001 From: mdecea Date: Mon, 13 Mar 2023 11:01:53 -0700 Subject: [PATCH 3/6] adding bbox creation to utils --- docs/photonics/examples/selecting_modes.py | 2 -- femwell/mode_solver.py | 9 ++------- femwell/utils.py | 20 ++++++++++++++++++++ 3 files changed, 22 insertions(+), 9 deletions(-) diff --git a/docs/photonics/examples/selecting_modes.py b/docs/photonics/examples/selecting_modes.py index ca9fb7ec..69a8c7cd 100644 --- a/docs/photonics/examples/selecting_modes.py +++ b/docs/photonics/examples/selecting_modes.py @@ -271,5 +271,3 @@ plot_mode(basis, np.real(xs[ind_mode]), direction="x") plt.show() - -# %% diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index fd223756..2249f5bd 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -22,6 +22,7 @@ ) from skfem.helpers import cross, curl, dot, grad, inner from skfem.utils import solver_eigen_scipy +from utils import create_bbox_basis def compute_modes( @@ -347,13 +348,7 @@ def select_mode_by_overlap(mode_basis, E_modes, H_modes, overlaps = list() if bbox is not None: - - def sel_fun(x): - return (x[0] < bbox[1]) * (x[0] > bbox[0]) * (x[1] > bbox[2]) * (x[1] < bbox[3]) - - selection_basis = Basis(mode_basis.mesh, mode_basis.elem, - elements = lambda x: sel_fun(x) - ) + selection_basis = create_bbox_basis(mode_basis, bbox) elif elements_list is not None: diff --git a/femwell/utils.py b/femwell/utils.py index 0bd40a33..09c0f415 100644 --- a/femwell/utils.py +++ b/femwell/utils.py @@ -11,6 +11,7 @@ from scipy.sparse import spmatrix from skfem.utils import CondensedSystem, bmat +from skfem import Basis def mpc_symmetric( @@ -92,3 +93,22 @@ def mpc_symmetric( lambda x: np.concatenate((x, T @ x[len(U) :] + g)), ), ) + +def create_bbox_basis(basis, bbox): + """ + Creates a basis that only contains elements in the + given bounding box + + Args: + basis: original basis with all the elements + bbox: 4 element list with [x_min, x_max, y_min, y_max] + """ + + def sel_fun(x): + return (x[0] < bbox[1]) * (x[0] > bbox[0]) * (x[1] > bbox[2]) * (x[1] < bbox[3]) + + selection_basis = Basis(basis.mesh, basis.elem, + elements = lambda x: sel_fun(x) + ) + + return selection_basis From 0b6fdf11ed928ef7b2fad3a150196eb65e3a10aa Mon Sep 17 00:00:00 2001 From: mdecea Date: Tue, 14 Mar 2023 13:35:41 -0700 Subject: [PATCH 4/6] selection mode improvement --- docs/photonics/examples/selecting_modes.py | 75 ++++++++---------- femwell/mode_solver.py | 90 ++++++++++------------ femwell/utils.py | 22 +++--- 3 files changed, 80 insertions(+), 107 deletions(-) diff --git a/docs/photonics/examples/selecting_modes.py b/docs/photonics/examples/selecting_modes.py index 69a8c7cd..95ac41e7 100644 --- a/docs/photonics/examples/selecting_modes.py +++ b/docs/photonics/examples/selecting_modes.py @@ -17,8 +17,8 @@ # %% [markdown] # # Selecting modes in the mode solver -# Sometimes we have structures where the mode of interest is -# not the mode with the highest effective index. There are a few +# Sometimes we have structures where the mode of interest is +# not the mode with the highest effective index. There are a few # ways to select modes of interest in femwell # %% tags=["hide-input"] @@ -26,22 +26,23 @@ import matplotlib.pyplot as plt import numpy as np +import scipy.constants from scipy.constants import epsilon_0, speed_of_light from scipy.integrate import solve_ivp from shapely.geometry import Polygon from skfem import Basis, ElementTriP0, Mesh from skfem.io import from_meshio -import scipy.constants from femwell.mesh import mesh_from_OrderedDict from femwell.mode_solver import ( + argsort_modes_by_power_in_elements, calculate_coupling_coefficient, calculate_hfield, calculate_overlap, compute_modes, plot_mode, - select_mode_by_overlap ) +from femwell.utils import inside_bbox # %% [markdown] @@ -68,7 +69,7 @@ sin=Polygon( [ (-w_sin - gap / 2, 0), - (-w_sin - gap / 2, h_sin), + (-w_sin - gap / 2, h_sin), (-gap / 2, h_sin), (-gap / 2, 0), ] @@ -124,12 +125,10 @@ # If we use `find_modes`, these are the modes we get: -# %% +# %% # basis0.plot(epsilon, colorbar=True).show() -lams, basis, xs = compute_modes( - basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=4 -) +lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4, solver="scipy") plot_mode(basis, np.real(xs[0]), direction="x") @@ -141,7 +140,7 @@ plot_mode(basis, np.real(xs[3]), direction="x") plt.show() -print(f'The effective index of the SiN mode is {np.real(lams[2])}') +print(f"The effective index of the SiN mode is {np.real(lams[2])}") # %% [markdown] @@ -164,15 +163,13 @@ # thi smight not be an issue in this example but there's many # examples where this is not an acceptable option. -# %% +# %% epsilon = basis0.zeros() + 1.444**2 epsilon[basis0.get_dofs(elements=("si"))] = 1.444**2 epsilon[basis0.get_dofs(elements=("sin"))] = 1.973**2 -lams, basis, xs = compute_modes( - basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2 -) +lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2) plot_mode(basis, np.real(xs[0]), direction="x") @@ -180,27 +177,24 @@ plot_mode(basis, np.real(xs[1]), direction="x") plt.show() -print(f'The effective index of the SiN mode is {np.real(lams[0])}') +print(f"The effective index of the SiN mode is {np.real(lams[0])}") # %% [markdown] # ## 2. Giving a guess effective index # We can use the `n_guess` parameter to `compute_modes` to -# select modes close to that effective index. +# select modes close to that effective index. # This is great, but of course we need to know what's that guess # effective index. The way to do that would be to use option 1 above # and then use that as the n_guess. -# %% +# %% epsilon = basis0.zeros() + 1.444**2 epsilon[basis0.get_dofs(elements=("si"))] = 3.4777**2 epsilon[basis0.get_dofs(elements=("sin"))] = 1.973**2 -lams, basis, xs = compute_modes( - basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2, - n_guess = 1.62 -) +lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, n_guess=1.62) plot_mode(basis, np.real(xs[0]), direction="x") @@ -208,7 +202,7 @@ plot_mode(basis, np.real(xs[1]), direction="x") plt.show() -print(f'The effective index of the SiN mode is {np.real(lams[1])}') +print(f"The effective index of the SiN mode is {np.real(lams[1])}") # %% [markdown] @@ -216,9 +210,9 @@ # %% [markdown] -# ## 3. Using `select_mode_by_overlap` +# ## 3. Using `argsort_modes_by_power_in_elements` -# This allows to choose a mode that has the biggest overlap with +# This allows to choose a mode that has the biggest overlap with # a given structure. # There are two main ways to specify the structure: @@ -228,14 +222,12 @@ # You can also give it directly the selection_basis of the # are of interest. -# A requirement for using `select_mode_by_overlap` is to -# calculate the H field of the found modes. +# A requirement for using `argsort_modes_by_power_in_elements` is to +# calculate the H field of the found modes. -# %% +# %% -lams, basis, xs = compute_modes( - basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=4 -) +lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4) # Calculate H field H_modes = list() @@ -245,29 +237,26 @@ basis, E, lams[i] * (2 * np.pi / wavelength), - omega=2 * np.pi / 1.55 * - scipy.constants.speed_of_light, - )) + omega=2 * np.pi / 1.55 * scipy.constants.speed_of_light, + ) + ) # Option 1: using an element name -ind_mode = select_mode_by_overlap(basis, xs, H_modes, - bbox = None, - elements_list = "sin", - selection_basis = None) +ind_mode = argsort_modes_by_power_in_elements(basis, xs, H_modes, elements="sin")[0] plot_mode(basis, np.real(xs[ind_mode]), direction="x") plt.show() -# Option 1: using bounding box +# Option 2: using bounding box -# Format: [xmin, xmax, ymin, ymax] +# Format: [xmin, ymin, xmax, ymax] bbox = [-2, 0, 0, 0.4] -ind_mode = select_mode_by_overlap(basis, xs, H_modes, - bbox = bbox, - elements_list = None, - selection_basis = None) +elements = inside_bbox(bbox) +ind_mode = argsort_modes_by_power_in_elements(basis, xs, H_modes, elements)[0] plot_mode(basis, np.real(xs[ind_mode]), direction="x") plt.show() + +# %% diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 2249f5bd..4e2e74ed 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -22,7 +22,6 @@ ) from skfem.helpers import cross, curl, dot, grad, inner from skfem.utils import solver_eigen_scipy -from utils import create_bbox_basis def compute_modes( @@ -327,43 +326,23 @@ def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", directi return fig, axs -def select_mode_by_overlap(mode_basis, E_modes, H_modes, - bbox = None, - elements_list = None, - selection_basis = None): - """ - Selects the mode in the "modes" list that has the highest power contained - within the given bounding box, element or basis. - - If bbox is not None, we select the overlap with the given bounding box - [x_min, x_max, y_min, y_max] - - If elements_list is not None, then the overlap is calculated between the - given modes and the elements in the list +def argsort_modes_by_power_in_elements(mode_basis, E_modes, H_modes, elements): + """Sorts the modes in the "modes" list by the power contained + within the given elements. - If elements_list is None, then we calculate the overlap with the given - selection_basis. + Returns: + the indices sorted from highest to lowest power. """ - overlaps = list() + selection_basis = Basis(mode_basis.mesh, mode_basis.elem, elements=elements) - if bbox is not None: - selection_basis = create_bbox_basis(mode_basis, bbox) - - elif elements_list is not None: - - selection_basis = Basis(mode_basis.mesh, mode_basis.elem, elements=elements_list) + overlaps = [ + calculate_overlap(selection_basis, E_f, H_f, selection_basis, E_f, H_f) + for E_f, H_f in zip(E_modes, H_modes) + ] - for E_f, H_f in zip(E_modes, H_modes): - overlaps.append(calculate_overlap(selection_basis, E_f, H_f, selection_basis, E_f, H_f)) - - ind_max = np.argmax(np.abs(overlaps)) + return np.argsort(np.abs(overlaps))[::-1] - #print(overlaps) - - return ind_max - - if __name__ == "__main__": from collections import OrderedDict @@ -423,8 +402,7 @@ def select_mode_by_overlap(mode_basis, E_modes, H_modes, # basis0.plot(epsilon, colorbar=True).show() lams, basis, xs = compute_modes( - basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=6, order=2, radius=3, - solver='scipy' + basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=6, order=2, radius=3, solver="scipy" ) print(lams) @@ -467,25 +445,35 @@ def select_mode_by_overlap(mode_basis, E_modes, H_modes, plt.imshow(np.real(integrals)) plt.colorbar() plt.show() - + # Create basis to select a certain simulation extent def sel_fun(x): print(x) return (x[0] < 0) * (x[0] > -1) * (x[1] > 0) * (x[1] < 0.5) - selection_basis = Basis(basis.mesh, basis.elem, - #elements = lambda x: x[0] < 0 and x[0] > -1 and x[1] > 0 and x[1] < 0.5 - elements = lambda x: sel_fun(x) - ) - - print(select_mode_by_overlap(mode_basis = basis, - E_modes = xs, - H_modes = H_modes, - elements_list = ["core"], - selection_basis = None)) - - print(select_mode_by_overlap(mode_basis = basis, - E_modes = xs, - H_modes = H_modes, - elements_list = None, - selection_basis = selection_basis)) + selection_basis = Basis( + basis.mesh, + basis.elem, + # elements = lambda x: x[0] < 0 and x[0] > -1 and x[1] > 0 and x[1] < 0.5 + elements=lambda x: sel_fun(x), + ) + + print( + select_mode_by_overlap( + mode_basis=basis, + E_modes=xs, + H_modes=H_modes, + elements_list=["core"], + selection_basis=None, + ) + ) + + print( + select_mode_by_overlap( + mode_basis=basis, + E_modes=xs, + H_modes=H_modes, + elements_list=None, + selection_basis=selection_basis, + ) + ) diff --git a/femwell/utils.py b/femwell/utils.py index 09c0f415..56458c36 100644 --- a/femwell/utils.py +++ b/femwell/utils.py @@ -10,8 +10,8 @@ else: from scipy.sparse import spmatrix -from skfem.utils import CondensedSystem, bmat from skfem import Basis +from skfem.utils import CondensedSystem, bmat def mpc_symmetric( @@ -94,21 +94,17 @@ def mpc_symmetric( ), ) -def create_bbox_basis(basis, bbox): + +def inside_bbox(bbox): """ - Creates a basis that only contains elements in the - given bounding box + Creates a selection function that is True for elements in the + given bounding box coordinates. Args: - basis: original basis with all the elements - bbox: 4 element list with [x_min, x_max, y_min, y_max] + bbox: 4 element list with [xmin, ymin, xmax, ymax] """ def sel_fun(x): - return (x[0] < bbox[1]) * (x[0] > bbox[0]) * (x[1] > bbox[2]) * (x[1] < bbox[3]) - - selection_basis = Basis(basis.mesh, basis.elem, - elements = lambda x: sel_fun(x) - ) - - return selection_basis + return (x[0] < bbox[2]) * (x[0] > bbox[0]) * (x[1] > bbox[1]) * (x[1] < bbox[3]) + + return sel_fun From 0339aecf1bc05bf047685da3f4bc0d37ef3a1fa1 Mon Sep 17 00:00:00 2001 From: mdecea Date: Tue, 14 Mar 2023 13:39:33 -0700 Subject: [PATCH 5/6] Updating table of contents --- docs/_toc.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/_toc.yml b/docs/_toc.yml index 6d7602f7..9ecc6a37 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -18,6 +18,7 @@ parts: - file: photonics/examples/bent_waveguide.py - file: photonics/examples/vary_width.py - file: photonics/examples/vary_wavelength.py + - file: photonics/examples/selecting_modes.py - file: photonics/examples/fiber_overlap.py - file: photonics/examples/metal_heater_phase_shifter.py - file: photonics/examples/metal_heater_phase_shifter_transient.py From d77d83780c0bbae1a3182337a385633509661d12 Mon Sep 17 00:00:00 2001 From: mdecea Date: Tue, 14 Mar 2023 13:56:27 -0700 Subject: [PATCH 6/6] removing scipy solver --- docs/photonics/examples/selecting_modes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/photonics/examples/selecting_modes.py b/docs/photonics/examples/selecting_modes.py index 95ac41e7..c9a9ef53 100644 --- a/docs/photonics/examples/selecting_modes.py +++ b/docs/photonics/examples/selecting_modes.py @@ -128,7 +128,7 @@ # %% # basis0.plot(epsilon, colorbar=True).show() -lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4, solver="scipy") +lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4) plot_mode(basis, np.real(xs[0]), direction="x")