GitHub

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)
Example block output

This page was generated using Literate.jl.