overlap fails when bases are different

[AmitRotem]

overlap integral of 2 waveguides with different widths;

#
from collections import OrderedDict
import numpy as np
import shapely
from shapely.ops import clip_by_rect
from skfem import Basis, ElementTriP0
from skfem.io.meshio import from_meshio

from femwell.maxwell.waveguide import compute_modes
from femwell.mesh import mesh_from_OrderedDict
from femwell.visualization import plot_domains


def trapezoid(wg_width, wg_thickness, angle, x_shift=0.0, y_shift=0.0):
    lower_shift = wg_thickness*np.tan((90-angle)*np.pi/180)
    bottom_left  = (-wg_width/2+lower_shift, 0.0)
    bottom_right = (+wg_width/2-lower_shift, 0.0)
    top_left     = (-wg_width/2            , wg_thickness)
    top_right    = (+wg_width/2            , wg_thickness)
    core = shapely.geometry.Polygon((bottom_right, top_right, top_left, bottom_left))
    return shapely.affinity.translate(core, x_shift, y_shift)


def single_rail_setup(wg_width, wg_thickness=0.22, angle=80, clad_thickness=2.2, core_index=3.5, clad_index=1.444,
                      box_index=1.444, air_index=1.0003, max_height=None, env_shape_res=8, env_xfact=1,
                      wg_res=0.02, max_res=1, grading_dist=3):
    if max_height is None:
        max_height = min(wg_width, wg_thickness)*7
    # geometry
    core = trapezoid(wg_width, wg_thickness, angle)
    env = shapely.affinity.scale(core.buffer(max_height, resolution=env_shape_res), xfact=env_xfact)

    polygons = OrderedDict(
        core=core,
        sub=clip_by_rect(env, -np.inf, -np.inf, np.inf, 0),
        clad=clip_by_rect(env, -np.inf, 0, np.inf, clad_thickness),
    )
    if max_height>clad_thickness-wg_thickness:
        polygons['air'] = clip_by_rect(env, -np.inf, clad_thickness, np.inf, np.inf)

    # mesh
    resolutions = dict(core={"resolution": wg_res, "distance": grading_dist})
    mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=max_res))

    # refraction index
    basis0 = Basis(mesh, ElementTriP0())
    epsilon = basis0.zeros()
    index_dict = {"core": core_index, "sub": box_index, "clad": clad_index}
    if max_height>clad_thickness-wg_thickness:
        index_dict["air"] = air_index
    for subdomain, n in index_dict.items():
        epsilon[basis0.get_dofs(elements=subdomain)] = n**2

    return basis0, epsilon, mesh


def single_rail(wavelength, num_modes=1, mode_order=2, **kwargs):
    # setup
    basis0, epsilon, mesh = single_rail_setup(**kwargs)
    # compute modes
    modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=num_modes, order=mode_order)
    return modes

mode1 = single_rail(1.56, wg_width=0.5)
mode2 = single_rail(1.56, wg_width=0.51)
mode1[0].calculate_overlap(mode2[0])

Seems to get pass this line (last changed in #140)

femwell/femwell/maxwell/waveguide.py

Lines 680 to 682 in b925773

	if basis_i == basis_j or (
	np.isclose(basis_i.X, basis_j.X).all() and np.isclose(basis_i.W, basis_j.W).all()
	):

and fail with

return A[0] * B[1] - A[1] * B[0]
           ~~~~~^~~~~~
ValueError: operands could not be broadcast together with shapes (2123,3) (2161,3)

[HelgeGehring]

You'd need to have both simulations on the same grid.
Otherwise the quadrature points won't match - of course we could interpolate but that'll lead to quite some inaccuracies.

Thus, what I'd recommend:

Make one mesh in which you include both geometries, I.e. core1 and core2 (whereas core2 is the larger one and is only core2-core1 in the final mesh as you put core1 on higher priority in the meshing)

then you do two simulations - in one you set core2 to background (to simulate your smaller waveguide) and in the other one you set it also to core_index (to simulate the larger waveguide)

That way you have the exact same mesh for both simulations (saved one meshing also!) and thus the quadrature points are also at the exact same location which leads to great accuracy.

Hope this helps!
Let me know how this works for you! (and feel free to make an example out of this if you want :) )

[AmitRotem]

That's definitely the smart way of doing it!
It seems very simple to adopt if core_j covers core_(j-1).
Is there a simple way to extend this to multiple geometries with non-trivial overlap?

In any case, I would still expect to have the interpolation as fallback when doing this the naive way.
In my example; basis_i == basis_j returns False but (np.isclose(basis_i.X, basis_j.X).all() and np.isclose(basis_i.W, basis_j.W).all() returns True, and the overlap is calculated as if the meshs are identical, which is not the case.

[HelgeGehring]

Hmm, having too many geometries in one mesh will lead to a lot of very small mesh elements - how many are you planning / what geometries? and what are you planning to do with that?

The test (np.isclose(basis_i.X, basis_j.X).all() and np.isclose(basis_i.W, basis_j.W).all()) seems not enough to me 🤔
It just tests the locations and weights of the quadrature element within the reference element I think....
We'd need to check there more things of being equal, probably we'd just want to check if basis.mesh is the exact same (in principal we should just check if all properties of basis.mesh are the same.... but in most cases I think the mesh would just be reused) and also if basis.element is an object of the same class?
That way you'd go straight to the interpolation 🤔

[AmitRotem]

Eigen mode propagation, e.g., tapers, bends, splitters. usually through 20-30 slices.
I could probably hand craft a solution to each individual scenario, but it would be nice to have something generalized.

For the simple case; `core_j` covers `core_(j-1)`

# parameters
## loop over single rail widths
params = [dict(
    wg_width=wg_width,
    wg_thickness=0.22,
    angle=80,
) for wg_width in np.linspace(0.4,0.6,3)]

## other parameters
clad_thickness=2.2
core_index=1.997
clad_index=1.444
box_index=1.444
air_index=1.0003

## mesh parameters
env_shape_res=8
env_xfact=1
wg_res=0.1
max_res=1
grading_dist=3
max_height=None
if max_height is None:
    min_wg_width = np.min([param['wg_width'] for param in params])
    min_wg_thickness = np.min([param['wg_thickness'] for param in params])
    max_height = min(min_wg_width, min_wg_thickness)*7

# geometry
core_list = [trapezoid(param['wg_width'], param['wg_thickness'], param['angle'], param['x_shift']) for param in params]
for core_current,core_next in zip(core_list[:-1], core_list[1:]):
    #* if this assert fails mesh_from_OrderedDict will fail
    assert core_current.covered_by(core_next)
# simulation environment
env = shapely.affinity.scale(shapely.ops.unary_union(core_list).buffer(max_height, resolution=env_shape_res), xfact=env_xfact)
# substrate
sub = clip_by_rect(env, -np.inf, -np.inf, np.inf, 0)
# clad - `clad_thickness` above the substrate+core
clad = shapely.affinity.scale(shapely.ops.unary_union(core_list+[sub]).buffer(clad_thickness, resolution=env_shape_res), xfact=env_xfact)
clad = clad.intersection(env)

## build polygons dict
polygons = OrderedDict()
for i, core in enumerate(core_list):
    polygons[f'core_{i}'] = core
polygons['sub'] = sub
polygons['clad'] = clad
### if clad covers the environment, we do not need to add it (also, mesh_from_OrderedDict will fail if we add it)
if not env.covered_by(clad):
    polygons['env'] = env

# mesh
resolutions = {f'core_{i}': {"resolution": wg_res, "distance": grading_dist} for i in range(len(core_list))}
mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=max_res))
## plot mesh
mesh.draw()
plt.title("Mesh")
plt.show(block=False)
## plot domains
plot_domains(mesh)
plt.title("Regions")
plt.grid()
plt.show(block=False)

# refraction index
## defualt index dict - all cores have the clad index
base_index_dict = {**{f"core_{i}": clad_index for i in range(len(core_list))}, **{"sub": box_index, "clad": clad_index}}
if not env.covered_by(clad):
    base_index_dict['env'] = air_index
## list of index dicts - each one corresponds to a core
index_dict_list = [base_index_dict.copy() for _ in range(len(core_list))]
for i in range(len(core_list)):
    for j in range(i+1):
        # core_j is covered by core_i
        # and they are ordered in priority, so we need to set the index of all core_i
        index_dict_list[i][f"core_{j}"] = core_index
## create basis and epsilon
basis0 = Basis(mesh, ElementTriP0())
epsilon_list = [basis0.zeros() for _ in range(len(core_list))]
for index_dict, epsilon in zip(index_dict_list, epsilon_list):
    # set refractive index for each subdomain
    for subdomain, n in index_dict.items():
        epsilon[basis0.get_dofs(elements=subdomain)] = n**2

## plot refraction index
for i,epsilon in enumerate(epsilon_list):
    basis0.plot(np.sqrt(epsilon), colorbar=True)
    plt.title(f"Refractive index for core {i}")
plt.show(block=False)

# compute modes
modes_list = [compute_modes(basis0, epsilon, wavelength=1.55, num_modes=2, order=2) for epsilon in tqdm(epsilon_list)]

# calculate overlaps
np.array([[mi[0].calculate_overlap(mj[0]) for mi in modes_list] for mj in modes_list])

[HelgeGehring]

Okay, for eigenmode expansion having all in the same mesh would be too much (I might have tried :D )

Key for that will be to use super meshing, see https://github.com/kinnala/scikit-fem/blob/master/skfem/supermeshing.py
That creates a mesh which includes all nodes from j/j-1. From what I understand the results should be precise if we calculate the overlap integrals based on the quadrature points of the supermesh.