GitHub

Homogeneous case with Neumann boundary condition

Problem

This problem in given as an example in [3].

  • The function $q(x_{1}, x_{2}) = 1$ is a constant function.
  • Neumann boundary condition on $\partial \Omega$

Code

First we should load ClosedWaveguideDispersion.jl for implemented functions. We also need Ferrite.jl to define the interpolation.

using ClosedWaveguideDispersion
using Ferrite: Lagrange, RefTriangle

Since we consider the homogeneous waveguide, we need to define the refraction index.

# Refractive index
function n(x)
   return 1.0
end;

# Period of the waveguide
p = 1.0;
# Height of the waveguide
h = 1.0;
# Number of points in the discrete Brillouin zone
N = 100;
100

We use setup_grid to generate mesh for the periodic cell with period p and height h.

# Set up the grid
grid = setup_grid(lc=0.05, period=p, height=h)
Ferrite.Grid{2, Ferrite.Triangle, Float64} with 944 Ferrite.Triangle cells and 513 nodes

Then we need to define the interpolation, CellValues and DofHandler which are needed by Ferrite.jl.

# 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

In this example, the Neumann boundary condition is satisfied naturally. So we only need to impose the periodic boundary condition in setup_bdcs.

# Set the boundary conditions
cst = setup_bdcs(dh, period=p)
ConstraintHandler:
  BCs:

Now we should prepare for the assembly of the eigenvalue problem.

# Allocate the matrices
A = allocate_matries(dh, cst)
B = allocate_matries(dh, cst)

# Discretize the Brillouin zone
bz = collect(range(-π/p, π/p, N))
100-element Vector{Float64}:
 -3.141592653589793
 -3.0781261353354537
 -3.0146596170811146
 -2.951193098826775
 -2.887726580572436
 -2.8242600623180967
 -2.7607935440637577
 -2.697327025809418
 -2.633860507555079
 -2.5703939893007397
  ⋮
  2.633860507555079
  2.697327025809418
  2.7607935440637577
  2.8242600623180967
  2.887726580572436
  2.951193098826775
  3.0146596170811146
  3.0781261353354537
  3.141592653589793

In calc_diagram, B is generated only one time and A is generated with respect to $\alpha$ in bz. Then we solve the generalized linear eigenvalue problem $Ax = \mu Bx$ at each $\alpha$.

# Calculate the dispersion diagram
μ = calc_diagram(cv, dh, cst, A, B, n, bz, nevs=7)
100×7 Matrix{Float64}:
 9.86864  10.1161  19.7522  20.1265  49.5831  50.343   89.0705
 9.47394  10.5189  19.3575  20.5293  49.1884  50.7458  87.8782
 9.08729  10.9297  18.9709  20.9402  48.8018  51.1567  86.6941
 8.70869  11.3487  18.5923  21.3592  48.4233  51.5756  85.518
 8.33815  11.7757  18.2218  21.7862  48.0527  52.0025  84.3499
 7.97567  12.2107  17.8594  22.2212  47.6903  52.4376  83.1899
 7.62124  12.6538  17.505   22.6644  47.3359  52.8806  82.038
 7.27487  13.1049  17.1586  23.1155  46.9896  53.3318  80.8941
 6.93655  13.5641  16.8203  23.5748  46.6513  53.7909  79.7582
 6.60628  14.0314  16.4901  24.042   46.321   54.2582  78.6305
 ⋮                                             ⋮       
 6.93655  13.5641  16.8203  23.5748  46.6513  53.7909  79.7582
 7.27487  13.1049  17.1586  23.1155  46.9896  53.3318  80.8941
 7.62124  12.6538  17.505   22.6644  47.3359  52.8806  82.038
 7.97567  12.2107  17.8594  22.2212  47.6903  52.4376  83.1899
 8.33815  11.7757  18.2218  21.7862  48.0527  52.0025  84.3499
 8.70869  11.3487  18.5923  21.3592  48.4233  51.5756  85.518
 9.08729  10.9297  18.9709  20.9402  48.8018  51.1567  86.6941
 9.47394  10.5189  19.3575  20.5293  49.1884  50.7458  87.8782
 9.86864  10.1161  19.7522  20.1265  49.5831  50.343   89.0705

Finally, we can plot our dispersion diagram.:

# Plot the dispersion diagram
plot_diagram(bz, μ, period=p)
Example block output

Plain code

using ClosedWaveguideDispersion
using Ferrite: Lagrange, RefTriangle

# Parameters
function n(x)
   return 1.0 
end

p = 1.0
h = 1.0
N = 100

# 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)

# Set the boundary conditions
cst = setup_bdcs(dh, period=p)

# Allocate the matrices
A = allocate_matries(dh, cst)
B = allocate_matries(dh, cst)

# 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)

# Plot the dispersion diagram
plot_diagram(bz, μ, period=p)