PA-Two

Problem

PA-Two is an example in [4] and satisfies the following conditions:

• a rectangular hole of size $0.5 \times 0.5$ in the cell center.
• Neumann boundary condition on $\partial \Omega$.
• homogeneous Dirichlet boundary condition at the hole boundary.

Code

For this example, we also need to load Gmsh and FerriteGmsh to generate the mesh

using Gmsh
using Ferrite
using FerriteGmsh
using ClosedWaveguideDispersion

We need to customize our mesh and boundary conditions.

function my_grid(;lc=0.05, period=1.0, height=1.0, holewidth=0.5, holeheight=0.5)
    # Initialize gmsh
    gmsh.initialize()
    gmsh.option.setNumber("General.Verbosity", 2)

    # Add the points
    p1 = gmsh.model.geo.addPoint(-period/2, 0, 0, lc)
    p2 = gmsh.model.geo.addPoint(period/2, 0, 0, lc)
    p3 = gmsh.model.geo.addPoint(period/2, height, 0, lc)
    p4 = gmsh.model.geo.addPoint(-period/2, height, 0, lc)
    p5 = gmsh.model.geo.addPoint(-holewidth/2, (height - holeheight)/2, 0, lc)
    p6 = gmsh.model.geo.addPoint(holewidth/2, (height - holeheight)/2, 0, lc)
    p7 = gmsh.model.geo.addPoint(holewidth/2, (height + holeheight)/2, 0, lc)
    p8 = gmsh.model.geo.addPoint(-holewidth/2, (height + holeheight)/2, 0, lc)

    # Add the lines
    l1 = gmsh.model.geo.addLine(p1, p2)
    l2 = gmsh.model.geo.addLine(p2, p3)
    l3 = gmsh.model.geo.addLine(p3, p4)
    l4 = gmsh.model.geo.addLine(p4, p1)
    l5 = gmsh.model.geo.addLine(p5, p6)
    l6 = gmsh.model.geo.addLine(p6, p7)
    l7 = gmsh.model.geo.addLine(p7, p8)
    l8 = gmsh.model.geo.addLine(p8, p5)

    # Create loops and the domain
    outerLoop = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
    innerLoop = gmsh.model.geo.addCurveLoop([l5, l6, l7, l8])
    domain = gmsh.model.geo.addPlaneSurface([outerLoop, innerLoop])

    # Synchronize the model
    gmsh.model.geo.synchronize()

    # Create the physical domains
    gmsh.model.addPhysicalGroup(1, [l1], -1, "bottom")
    gmsh.model.addPhysicalGroup(1, [l2], -1, "right")
    gmsh.model.addPhysicalGroup(1, [l3], -1, "top")
    gmsh.model.addPhysicalGroup(1, [l4], -1, "left")
    gmsh.model.addPhysicalGroup(1, [l5, l6, l7, l8], -1, "Γ")
    gmsh.model.addPhysicalGroup(2, [domain], -1, "Ω")

    # Set Periodic boundary condition
    transform = [1, 0, 0, period, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1]
    gmsh.model.mesh.setPeriodic(1, [l2], [l4], transform)

    # Generate a 2D mesh
    gmsh.model.mesh.generate(2)

    grid = mktempdir() do dir
        path = joinpath(dir, "hole.msh")
        gmsh.write(path)
        togrid(path)
    end

    # Finalize the Gmsh library
    gmsh.finalize()

    return grid
end;

We still use homogeneous medium in the waveguide.

function n(x)
    return 1.0
end;

We need to implement a new function to impose the boundary condition.

function my_bdcs(dh::DofHandler; period=1.0)
    cst = ConstraintHandler(dh)

    # Periodic boundary condition on the left and right
    # side of the periodic cell
    pfacets = collect_periodic_facets(dh.grid, "right", "left", x -> x + Ferrite.Vec{2}((period, 0.0)))
    pbc = PeriodicDirichlet(:u, pfacets)
    add!(cst, pbc)

    # Set Dirichlet boundary condition on the inner boundary
    dfacets = getfacetset(dh.grid, "Γ")
    dbc = Dirichlet(:u, dfacets, x -> 0)
    add!(cst, dbc)

    close!(cst)

    return cst
end

my_bdcs (generic function with 1 method)

Set parameters

p = 1.0;
h = 1.0;
hw = 0.5;
hh = 0.5;
N = 150;

Set up the grid

grid = my_grid(lc=0.05, period=p, height=h, holewidth=hw, holeheight=hh)

Ferrite.Grid{2, Ferrite.Triangle, Float64} with 730 Ferrite.Triangle cells and 425 nodes

Basic steps need 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)

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

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

425×425 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2819 stored entries:
⎡⠻⢆⡄⢸⣦⠘⡠⠀⠈⠀⠁⠑⠀⠂⠀⠀⠼⠀⡀⠐⠀⣀⠁⢀⠀⠆⠀⠀⠀⡠⠀⢤⠆⠂⢀⠘⠃⠀⠀⠀⎤
⎢⣀⣉⣿⢟⠸⣱⡀⠠⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠧⡀⠁⣾⠰⠉⠀⠁⠀⠀⠀⠀⢈⠉⠇⡸⠈⣆⠢⠀⠠⠀⎥
⎢⣈⠛⢖⣢⡵⣯⡵⠀⣤⠀⡀⡨⠅⡀⠀⠀⠂⠀⢂⠀⠠⡢⢠⡤⠀⠂⢀⠀⠀⡍⠀⢀⠀⠀⠼⢨⣀⠀⠀⠐⎥
⎢⠀⠊⠀⡈⠑⠋⠿⣧⡈⠇⠀⠁⠀⠀⠀⠀⠢⠀⠀⠃⠀⡀⣨⠠⠀⠠⡸⠀⣠⡀⠐⡃⢂⣠⠸⡘⠁⠀⠀⡐⎥
⎢⠂⠀⠀⠐⠀⠛⠦⠌⡿⣯⢺⠂⣩⢀⣸⡀⠀⠀⠰⠂⠀⡄⠂⠐⠀⠐⠁⠀⠀⠉⠂⠁⠄⠀⡀⠈⠀⡀⢁⠀⎥
⎢⢅⠀⠀⠀⡀⡨⠄⠀⠺⠒⠻⣦⣰⡄⣀⡣⠤⡅⢤⡈⡁⠂⠀⠀⣀⡁⠘⠀⠀⠆⡀⠂⠀⠀⡂⠀⠎⠁⠂⠀⎥
⎢⠠⠀⠀⠀⠁⠡⠀⠀⠃⢚⠐⠾⠱⣦⡜⠦⠀⠆⡄⣠⠀⡄⠀⠀⡈⠈⠀⠁⠈⠆⠀⠂⠀⠀⠀⠀⠀⡘⢁⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠒⠺⠤⡸⠲⡍⣱⣾⢀⡄⠉⠱⢅⠁⠂⠀⠂⠀⠀⠀⠀⡀⣠⠀⡀⠀⢀⠀⠀⠲⣜⠀⎥
⎢⠒⠃⠀⠀⠈⠀⠈⠂⠀⠀⠄⠧⠠⠄⠀⠴⠟⣥⢠⡤⡶⠄⠒⡄⣀⠄⠈⠜⠅⠄⠀⠀⠆⠀⣵⠀⢠⢃⠴⣀⎥
⎢⢀⠈⠉⠣⠈⠐⠤⠀⠰⠂⡀⠳⠀⣩⢇⡀⠀⡶⠿⢇⡂⠊⠴⢘⡖⣂⡀⠀⠀⠀⡂⢀⠀⢀⠐⠠⡈⡠⢿⠨⎥
⎢⠀⢠⣡⣤⠠⡢⠀⠠⠀⠤⠡⠈⠀⠤⠅⠑⠘⠏⡨⠈⣻⣾⣱⣵⢐⢐⢖⡃⡐⡂⠁⠔⢍⠀⠈⢀⠠⠀⠠⠀⎥
⎢⠁⢀⡔⠂⠀⡶⠂⡚⢈⠀⠀⠀⠀⠀⠈⠀⠘⠤⣐⢃⢕⣾⢟⢕⠥⣀⣀⠂⢀⠀⠀⠀⠅⠄⠰⠘⠢⠁⠀⠧⎥
⎢⠠⠄⠄⠀⠠⠀⠀⡀⢀⠀⠄⠸⡂⠈⠈⠀⠀⠜⠸⢩⢐⢐⠁⢣⣻⢞⣗⡂⣶⠀⡂⠀⢄⢀⠾⠤⠠⠙⠌⠁⎥
⎢⠀⠀⠀⠀⠀⠐⠒⠊⠁⠀⠒⠀⠄⠀⠀⠀⣂⠄⠀⠈⠼⠱⠠⠘⠹⠹⠵⣧⡶⠊⢄⠐⠀⠈⠃⠠⠂⢀⣄⠕⎥
⎢⠀⡠⠀⠀⡄⠤⠀⠺⡄⠀⠠⠄⠢⠄⠀⠠⠁⠅⠀⠀⠰⠨⠀⠐⠘⠛⡸⠋⠻⢆⣈⡡⠈⠘⠄⢇⢈⢙⠀⠀⎥
⎢⠀⣄⡆⠐⠀⢀⠴⠠⠌⠀⠠⠈⠠⠀⠀⠚⠀⠀⠈⢈⢁⠄⠀⠀⠈⠈⢀⠑⠆⡸⠻⣦⠀⢀⠐⠊⡐⠘⠁⠀⎥
⎢⠨⠁⣉⡡⠀⠀⠈⣰⠀⠁⠀⠀⠀⠀⠀⠈⠈⠁⠀⢀⠃⠑⠁⠅⠀⢑⡀⠀⣂⠀⠀⢀⢿⢗⣁⡰⢀⣂⠈⠄⎥
⎢⣀⠐⠢⢤⡒⣃⣒⠢⡀⠈⠈⠈⠀⠀⠀⠐⠑⠛⠐⡀⠂⢀⣐⠂⠚⡇⠉⡀⠤⢅⡰⠀⢁⡸⠛⢄⠾⡑⡈⣉⎥
⎢⠉⠀⠈⠂⠀⠘⠁⠀⠀⠠⠎⠁⣀⠠⢠⡀⠤⢒⠂⡨⠀⠂⠌⠂⣄⠂⠈⢀⣆⢐⣐⠈⠠⢰⢞⠣⢱⣶⡰⠔⎥
⎣⠀⠀⠀⠂⢀⠀⢀⠠⠁⠐⠈⠀⠁⠐⠒⠙⠐⢣⡛⡓⠀⠂⠤⡄⠆⠁⢄⠝⠀⠀⠁⠀⠂⠄⡆⢨⢐⠎⠛⣤⎦

We discretize the Brillouin zone and calculate the dispersion diagram.

bz = collect(range(-π/p, 2π/p, N))
μ = calc_diagram(cv, dh, cst, A, B, n, bz, nevs=5)

150×5 Matrix{Float64}:
 27.0388  30.6954  48.0029  48.0246  67.2005
 27.0229  30.6755  48.0172  48.0494  67.1889
 27.0006  30.6455  48.0488  48.0989  67.1805
 26.9722  30.6057  48.0976  48.1731  67.1753
 26.9377  30.5564  48.1634  48.2715  67.1732
 26.8973  30.498   48.2457  48.3936  67.1745
 26.8514  30.4311  48.3444  48.5391  67.1792
 26.8001  30.3562  48.4589  48.7072  67.1873
 26.7438  30.2738  48.5888  48.8975  67.199
 26.6828  30.1848  48.7337  49.1092  67.2144
  ⋮                                  
 25.0369  27.659   57.2561  63.5987  72.456
 25.0159  27.6259  57.4625  64.1934  72.7955
 24.9992  27.5981  57.6546  64.7854  73.1449
 24.9867  27.5756  57.8319  65.3714  73.5019
 24.9785  27.5585  57.9939  65.9469  73.8631
 24.9748  27.5469  58.1402  66.506   74.2237
 24.9755  27.5409  58.2706  67.0402  74.5775
 24.9808  27.5404  58.3847  67.5383  74.9157
 24.9905  27.5455  58.4824  67.9844  75.2268

Now we can plot the dispersion diagram.

plot_diagram(bz, μ, period=p)

Plain code

using Gmsh
using Ferrite
using FerriteGmsh
using ClosedWaveguideDispersion

function my_grid(;lc=0.05, period=1.0, height=1.0, holewidth=0.5, holeheight=0.5)
    # Initialize gmsh
    gmsh.initialize()
    gmsh.option.setNumber("General.Verbosity", 2)

    # Add the points
    p1 = gmsh.model.geo.addPoint(-period/2, 0, 0, lc)
    p2 = gmsh.model.geo.addPoint(period/2, 0, 0, lc)
    p3 = gmsh.model.geo.addPoint(period/2, height, 0, lc)
    p4 = gmsh.model.geo.addPoint(-period/2, height, 0, lc)
    p5 = gmsh.model.geo.addPoint(-holewidth/2, (height - holeheight)/2, 0, lc)
    p6 = gmsh.model.geo.addPoint(holewidth/2, (height - holeheight)/2, 0, lc)
    p7 = gmsh.model.geo.addPoint(holewidth/2, (height + holeheight)/2, 0, lc)
    p8 = gmsh.model.geo.addPoint(-holewidth/2, (height + holeheight)/2, 0, lc)

    # Add the lines
    l1 = gmsh.model.geo.addLine(p1, p2)
    l2 = gmsh.model.geo.addLine(p2, p3)
    l3 = gmsh.model.geo.addLine(p3, p4)
    l4 = gmsh.model.geo.addLine(p4, p1)
    l5 = gmsh.model.geo.addLine(p5, p6)
    l6 = gmsh.model.geo.addLine(p6, p7)
    l7 = gmsh.model.geo.addLine(p7, p8)
    l8 = gmsh.model.geo.addLine(p8, p5)

    # Create loops and the domain
    outerLoop = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
    innerLoop = gmsh.model.geo.addCurveLoop([l5, l6, l7, l8])
    domain = gmsh.model.geo.addPlaneSurface([outerLoop, innerLoop])

    # Synchronize the model
    gmsh.model.geo.synchronize()

    # Create the physical domains
    gmsh.model.addPhysicalGroup(1, [l1], -1, "bottom")
    gmsh.model.addPhysicalGroup(1, [l2], -1, "right")
    gmsh.model.addPhysicalGroup(1, [l3], -1, "top")
    gmsh.model.addPhysicalGroup(1, [l4], -1, "left")
    gmsh.model.addPhysicalGroup(1, [l5, l6, l7, l8], -1, "Γ")
    gmsh.model.addPhysicalGroup(2, [domain], -1, "Ω")

    # Set Periodic boundary condition
    transform = [1, 0, 0, period, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1]
    gmsh.model.mesh.setPeriodic(1, [l2], [l4], transform)

    # Generate a 2D mesh
    gmsh.model.mesh.generate(2)

    grid = mktempdir() do dir
        path = joinpath(dir, "hole.msh")
        gmsh.write(path)
        togrid(path)
    end

    # Finalize the Gmsh library
    gmsh.finalize()

    return grid
end

function n(x)
    return 1.0
end

function my_bdcs(dh::DofHandler; period=1.0)
    cst = ConstraintHandler(dh)

    # Periodic boundary condition on the left and right
    # side of the periodic cell
    pfacets = collect_periodic_facets(dh.grid, "right", "left", x -> x + Ferrite.Vec{2}((period, 0.0)))
    pbc = PeriodicDirichlet(:u, pfacets)
    add!(cst, pbc)

    # Set Dirichlet boundary condition on the inner boundary
    dfacets = getfacetset(dh.grid, "Γ")
    dbc = Dirichlet(:u, dfacets, x -> 0)
    add!(cst, dbc)

    close!(cst)

    return cst
end

p = 1.0
h = 1.0
hw = 0.5
hh = 0.5
N = 150

grid = my_grid(lc=0.05, period=p, height=h, holewidth=hw, holeheight=hh)

ip = Lagrange{RefTriangle, 1}()
cv = setup_fevs(ip)
dh = setup_dofs(grid, ip)
cst = my_bdcs(dh, period=p)

A = allocate_matries(dh, cst)
B = allocate_matries(dh, cst)

bz = collect(range(-π/p, 2π/p, N))
μ = calc_diagram(cv, dh, cst, A, B, n, bz, nevs=5)
plot_diagram(bz, μ, period=p)

This page was generated using Literate.jl.