Inhomogeneous case
Problem
This problem is the second example in [3] page .
- $q(x_{1}, x_{2}) = 9$ in a disk with the center $(0, 0.5)$ and radius $0.3$ and $q(x_{1}, x_{2}) = 1$ outside the disk.
- Neumann boundary condition on $\partial \Omega$.
Code
We need the same packages as the homogeneous case
using Ferrite: Lagrange, RefTriangle
using ClosedWaveguideDispersion
Here we need to define a new refractive index.
function n(x)
r = sqrt(x[1]^2 + (x[2] - 0.5)^2)
if r <= 0.3
return 9.0
else
return 1.0
end
end;
p = 1.0;
h = 1.0;
N = 100;
Then we generate the mesh and set up variables for finite element method.
# Set up the grid
grid = setup_grid(lc=0.05, period=p, height=h)
# Define the interpolation
ip = Lagrange{RefTriangle, 1}()
# Set up the FE values: CellValues
cv = setup_fevs(ip)
# Set up the DofHandler
dh = setup_dofs(grid, ip)
Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}
Fields:
:u, Ferrite.Lagrange{Ferrite.RefTriangle, 1}()
Dofs per cell: 3
Total dofs: 513
For the boundary condition, we only need to use setup_bdcs
to impose the periodic boundary condition in $x_{1}$ direction.
# Set the boundary conditions
cst = setup_bdcs(dh, period=p)
ConstraintHandler:
BCs:
We assemble the matrices for the generalized linear eigenvalue problems.
# Allocate the matrices
A = allocate_matries(dh, cst)
B = allocate_matries(dh, cst)
513×513 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 3509 stored entries:
⎡⠛⣤⣤⢃⣚⣃⢂⢙⢀⡀⠀⠀⢀⠐⠃⠠⡐⠀⠀⠀⠀⡐⠠⠐⠠⢤⡰⠉⠠⠀⠰⠘⠁⠤⢀⠘⢢⠆⢥⡃⎤
⎢⠤⢛⣱⣾⡤⡁⢀⡀⠄⢪⠀⡀⠀⠀⡈⢈⠀⠈⠁⠀⠰⡀⠀⡄⠰⠊⣉⠀⠨⠰⢈⠀⢕⠱⠰⠫⠀⠁⠨⠄⎥
⎢⠾⢸⠄⠫⡻⢎⢜⣎⢵⡈⠀⠀⠈⠀⠀⢉⠐⠀⠡⠈⠂⠦⠀⠀⠀⠖⠊⡉⠁⠈⠨⠀⠄⠀⠠⠄⢢⠀⠾⡪⎥
⎢⣌⢐⠀⠰⡲⢵⣻⣾⢽⡆⠀⠀⠐⢠⢍⡠⠀⠄⡁⡃⠆⠄⡀⠀⠀⠀⢀⠰⠀⠀⠀⠲⠀⠀⠨⢐⡐⠀⡀⣋⎥
⎢⠀⠰⡠⣁⡑⠳⠳⠷⠿⣧⢀⣀⠀⣄⢁⠃⡒⠀⡇⠈⢠⠀⠈⠂⠂⠀⠂⠂⠀⠈⠀⠤⠀⠀⠐⠐⠱⣤⠁⠓⎥
⎢⠀⠀⠀⠠⠀⠀⠀⠀⠀⢰⣿⣿⣟⣽⠐⡴⠣⢐⢨⠀⣜⢃⠐⢄⢑⠩⠀⠀⠀⠈⠆⠀⠀⠀⢂⠄⠀⢒⠀⠀⎥
⎢⢀⠐⠀⠀⠂⠀⠐⣀⠀⢤⣟⣽⣟⣽⡼⣾⠂⠠⠑⠂⢴⠈⡐⡄⢁⢀⢄⣀⠀⠲⡂⠅⠈⠐⡦⠁⠘⠨⠁⠀⎥
⎢⠉⡀⡂⢈⡄⢀⠃⡱⠥⠐⢐⡤⣲⣯⣿⢟⡨⣨⠃⡒⡁⡀⠆⣒⡀⠠⢄⠀⠀⠠⡄⠂⠀⠀⠆⠆⢀⢀⠁⠐⎥
⎢⠐⠈⡀⠀⠐⠀⠀⠄⠘⠈⢉⢂⠈⡀⡂⣪⣟⣽⡞⠬⣊⠥⠘⡐⠃⢀⡀⢀⣤⣦⠂⡤⠀⠀⡀⢂⠀⢀⠂⠂⎥
⎢⠀⠀⠁⠀⡁⠂⠥⠨⡉⠉⠂⠒⠱⠀⢩⠠⡚⡍⠟⣥⡗⣘⣥⣈⠊⠀⡀⠀⠀⡊⠄⠄⠈⠥⠁⠀⠀⠐⢋⢉⎥
⎢⢀⠠⠐⠢⠨⡄⠈⠅⠀⠒⠶⢙⡐⠓⠁⠨⠎⡜⣙⢩⠛⢄⡁⣱⢄⠆⠦⡘⣀⡀⠅⠂⣠⠩⠞⠐⠀⠈⠈⠄⎥
⎢⢀⠂⠀⠤⠀⠀⠀⠈⠢⠀⠐⢄⠐⠬⢨⢡⢒⠠⡁⢻⢅⣨⡵⣯⣿⡎⣄⡥⠄⠀⡀⢤⠁⠈⡅⠍⠐⡶⠐⡋⎥
⎢⠀⣆⡰⠂⢠⠄⠀⠀⠈⠀⡕⡐⠁⢐⠀⡈⠉⢀⠊⠀⠠⠕⡻⠿⡱⣮⣥⠂⡤⠡⡥⢨⠉⠀⠄⡍⠈⠷⢐⠀⎥
⎢⡔⠊⠃⠘⡎⠠⢀⡐⠨⠀⠀⠀⠀⢱⠀⠑⠀⢈⠀⠈⣈⠣⠄⡽⠡⠛⠟⣥⣤⠍⡄⣥⠄⣤⣁⠬⠉⠀⠐⠀⎥
⎢⠀⠂⢂⡂⡁⠀⠀⠀⡀⠀⡀⠀⢠⡀⠀⡀⠠⣿⡠⠠⠀⠸⠀⠁⠄⡋⡄⠟⡿⣯⣷⣯⢢⠐⠂⠲⣣⠥⠀⠀⎥
⎢⣐⠂⠂⠐⠂⠂⢠⡀⠀⡄⠈⠁⠌⠌⠠⠉⠈⡤⠀⠅⠡⠁⠀⣌⡁⣋⠄⣭⡽⣿⢟⣵⣶⠙⣥⡠⣥⡱⠐⠀⎥
⎢⠁⡄⢕⡑⠀⠁⠀⠀⠀⠀⠀⠀⢂⠀⠀⠀⠀⠀⠆⡄⡄⡚⡁⠀⠃⠀⠀⣥⢈⠒⣜⠛⣻⣾⡬⣫⠖⠩⠌⠀⎥
⎢⣀⠐⡴⡂⠀⠆⢂⢂⢐⠀⠈⠔⠌⠋⠨⠅⠠⢈⠁⠀⢚⠁⡅⠍⡄⠥⡁⡜⢨⡀⠁⡻⡦⣫⣟⣽⣫⣒⠐⠀⎥
⎢⠨⠖⠄⠀⠈⠒⠐⠈⠑⣦⢠⢀⡒⡀⠀⢐⠀⢀⢀⠀⡀⠀⢰⡤⢦⡄⠃⠀⠍⡞⢅⡻⡜⡁⢫⢺⢿⣷⡃⠛⎥
⎣⠥⠳⠂⠆⡺⡣⡤⢨⢥⠀⠀⠀⠁⠀⢁⠀⠨⠀⡏⢐⠂⠄⡴⠠⠐⠐⠐⠀⠀⠀⠐⠀⠂⠁⠐⠀⣭⠈⢕⣵⎦
Then we solve the generalized eigenvalue problems at each $\alpha$ in $bz$,
# Discretize the Brillouin zone
bz = collect(range(-π/p, π/p, N))
# Calculate the dispersion diagram
μ = calc_diagram(cv, dh, cst, A, B, n, bz, nevs=7)
100×7 Matrix{Float64}:
1.93232 4.45248 6.3289 13.1637 15.8762 19.9986 22.1028
1.92551 4.46018 6.32281 13.1704 15.8692 19.9133 22.183
1.91114 4.47794 6.31192 13.1848 15.8586 19.802 22.2934
1.88955 4.50541 6.29632 13.2068 15.8445 19.6717 22.4269
1.86125 4.54209 6.27615 13.2363 15.8268 19.5281 22.5779
1.82684 4.58736 6.25158 13.2731 15.8056 19.3754 22.7422
1.78699 4.64056 6.2228 13.3169 15.7808 19.2168 22.9166
1.7424 4.70096 6.19007 13.3676 15.7525 19.0546 23.0988
1.69377 4.76789 6.15363 13.4248 15.7206 18.8904 23.2873
1.64174 4.84067 6.11378 13.4883 15.6853 18.7256 23.4807
⋮ ⋮
1.69377 4.76789 6.15363 13.4248 15.7206 18.8904 23.2873
1.7424 4.70096 6.19007 13.3676 15.7525 19.0546 23.0988
1.78699 4.64056 6.2228 13.3169 15.7808 19.2168 22.9166
1.82684 4.58736 6.25158 13.2731 15.8056 19.3754 22.7422
1.86125 4.54209 6.27615 13.2363 15.8268 19.5281 22.5779
1.88955 4.50541 6.29632 13.2068 15.8445 19.6717 22.4269
1.91114 4.47794 6.31192 13.1848 15.8586 19.802 22.2934
1.92551 4.46018 6.32281 13.1704 15.8692 19.9133 22.183
1.93232 4.45248 6.3289 13.1637 15.8762 19.9986 22.1028
Finally, we can get our dispersion diagram.
# Plot the dispersion diagram
plot_diagram(bz, μ, period=p)

This page was generated using Literate.jl.