From 4187a055feda19595c49af9e10a95a5f4f2ca5c2 Mon Sep 17 00:00:00 2001
From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com>
Date: Thu, 2 Nov 2023 16:52:05 -0700
Subject: [PATCH] adjust for bent waveguide

---
 docs/julia/waveguide.jl            | 37 ++++++++++++++++++++----------
 src/Maxwell/Waveguide/Waveguide.jl | 17 ++++++++++++--
 2 files changed, 40 insertions(+), 14 deletions(-)

diff --git a/docs/julia/waveguide.jl b/docs/julia/waveguide.jl
index 4fc8293b..07234f1b 100644
--- a/docs/julia/waveguide.jl
+++ b/docs/julia/waveguide.jl
@@ -32,23 +32,36 @@ clip_by_rect = pyimport("shapely.ops").clip_by_rect
 OrderedDict = pyimport("collections").OrderedDict
 mesh_from_OrderedDict = pyimport("femwell.mesh").mesh_from_OrderedDict
 
-wg_width = 2.5
-wg_thickness = 0.3
-core = shapely.geometry.box(-wg_width / 2, 0, +wg_width / 2, wg_thickness)
-env = shapely.affinity.scale(core.buffer(5, resolution = 8), xfact = 0.5)
+radius = 5
+wg_width = 0.5
+wg_thickness = 0.22
+slab_width = 3
+slab_thickness = 0.11
+core = shapely.geometry.box(radius - wg_width / 2, 0, radius + wg_width / 2, wg_thickness)
+slab = shapely.geometry.box(
+    radius - slab_width / 2,
+    0,
+    radius + slab_width / 2,
+    slab_thickness,
+)
+env = shapely.affinity.scale(core.buffer(3, resolution = 8), xfact = 0.5)
 
 polygons = OrderedDict(
     core = core,
+    slab = slab,
     box = clip_by_rect(env, -np.inf, -np.inf, np.inf, 0),
     clad = clip_by_rect(env, -np.inf, 0, np.inf, np.inf),
 )
 
-resolutions = Dict("core" => Dict("resolution" => 0.03, "distance" => 0.5))
+resolutions = Dict(
+    "core" => Dict("resolution" => 0.03, "distance" => 0.5),
+    "slab" => Dict("resolution" => 0.03, "distance" => 0.5),
+)
 
 mesh = mesh_from_OrderedDict(
     polygons,
     resolutions,
-    default_resolution_max = 10,
+    default_resolution_max = 1,
     filename = "mesh.msh",
 )
 
@@ -74,19 +87,19 @@ model = GmshDiscreteModel("mesh.msh")
 
 labels = get_face_labeling(model)
 
-epsilons = ["core" => 1.9963^2, "box" => 1.444^2, "clad" => 1.0^2]
+epsilons = ["core" => 3.5^2 + 0.001im, "slab" => 3.5^2, "box" => 1.444^2, "clad" => 1.0^2]
 ε(tag) = Dict(get_tag_from_name(labels, u) => v for (u, v) in epsilons)[tag]
 
 
 #dΩ = Measure(Ω, 1)
 τ = CellField(get_face_tag(labels, num_cell_dims(model)), Ω)
 
-modes = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 2, order = 1)
+modes = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 1, order = 1, radius = 5)
 println(n_eff(modes[1]))
 write_mode_to_vtk("mode", modes[1])
 
 plot_mode(modes[1])
-plot_mode(modes[2])
+#plot_mode(modes[2])
 modes
 
 # %% [markdown]
@@ -94,7 +107,7 @@ modes
 # Let's add a minor perturbation and calculate the effective refractive index using perturbation theory:
 
 # %% tags=["remove-stderr"]
-epsilons_p = ["core" => -0.1im, "box" => -0.0im, "clad" => -0.0im]
+epsilons_p = ["core" => -0.1im, "box" => -0.0im, "slab" => -0.0im, "clad" => -0.0im]
 ε_p(tag) = Dict(get_tag_from_name(labels, u) => v for (u, v) in epsilons_p)[tag]
 
 println(perturbed_neff(modes[1], ε_p ∘ τ))
@@ -104,7 +117,7 @@ println(perturbed_neff(modes[1], ε_p ∘ τ))
 
 # %% tags=["remove-stderr"]
 
-modes_p = calculate_modes(model, ε ∘ τ + ε_p ∘ τ, λ = 1.55, num = 2, order = 1)
-println(n_eff(modes_p[1]))
+#modes_p = calculate_modes(model, ε ∘ τ + ε_p ∘ τ, λ = 1.55, num = 2, order = 1)
+#println(n_eff(modes_p[1]))
 
 # %%
diff --git a/src/Maxwell/Waveguide/Waveguide.jl b/src/Maxwell/Waveguide/Waveguide.jl
index e3c026cd..a36ccf8e 100644
--- a/src/Maxwell/Waveguide/Waveguide.jl
+++ b/src/Maxwell/Waveguide/Waveguide.jl
@@ -113,14 +113,27 @@ function calculate_modes(
         )dΩ
     rhs((u1, u2), (v1, v2)) = ∫(-1 / μ_r * u1 ⊙ v1 / k0^2)dΩ
 
+    radius_function(x) = x[1]
+    r = interpolate_everywhere(radius_function, V2)
+    lhs_radial((u1, u2), (v1, v2)) =
+        ∫(
+            -(r / radius * curl(u1) ⋅ curl(v1)) +
+            (radius / r * gradient(u2) ⋅ v1) +
+            (k0^2 * ε * r / radius * u1 ⋅ v1) - (radius / r * gradient(u2) ⋅ gradient(v2)) +
+            (k0^2 * ε * radius / r * u2 * v2),
+        )dΩ
+
+    rhs_radial((u1, u2), (v1, v2)) =
+        ∫((radius / r * u1 ⋅ v1) - (radius / r * u1 ⋅ gradient(v2)))dΩ
+
     epsilons = ε(get_cell_points(Measure(Ω, 1)))
     k0_guess =
         isnothing(k0_guess) ? k0^2 * maximum(maximum.(maximum.(real.(epsilons)))) * 1.1 :
         k0_guess
 
     assem = Gridap.FESpaces.SparseMatrixAssembler(U, V)
-    A = assemble_matrix(lhs, assem, U, V)
-    B = assemble_matrix(rhs, assem, U, V)
+    A = assemble_matrix(lhs_radial, assem, U, V)
+    B = assemble_matrix(rhs_radial, assem, U, V)
     if all(imag(A.nzval) .== 0) && all(imag(B.nzval) .== 0)
         A, B = real(A), real(B)
     end