Finite element method with transparent boundary condition

Introduction

Consider a flat surface illuminated by a plane wave. The incident field is given by

\[u^{\text{in}} = e^{i\alpha x_{1} - i\beta x_{2}},\]

where $\alpha = k\sin\theta, \beta = k\cos\theta$; $k > 0$ is the wave number, and $\theta \in (-\pi/2, \pi/2)$ is the incident angle. We assume that the total field $u = u^{\text{in}} + u^{\text{sc}}$, composed of the incident field $u^{\text{in}}$ and the scattered field $u^{\text{sc}}$, satisfies the Dirichelt boundary condition on the bottom $x_{2} = 0$.

Due to the periodic structure and incident plane wave, we can work on periodic function space, that is, $v = e^{-i\alpha x_{1}}u$. So the exact solution is $v_{\text{exact}} = -2i\sin(\beta x_{2})$.

Numerical examples

First, we load Nmrc and Ferrite (Ferrite is used to define the interpolation).

using Nmrc, Ferrite

Parameters

We assume that the wavenumber is $1.0$ and the incident angle is $\pi/3$.

# incident plane wave
k = 1.0
θ = π/3
inc = Incident(k, θ)

# Number of the truncated terms
N = 10;

# Height of the artificial boundary
height = 2.0;

Generating a mesh by Gmsh

We generate a mesh on a periodic cell which is a rectangle $(0.0, 2\pi) \times (0.0, 2.0)$.

# Generate the mesh in a periodic cell
grid = periodic_cell(lc=0.1, period=2π, height=height)

Ferrite.Grid{2, Ferrite.Triangle, Float64} with 2980 Ferrite.Triangle cells and 1574 nodes

Setting up finite element space and boundary conditions

We choose $P1$ element here. Then CellValues and FacetValues can be constructed by the function setup_vals.

ip = Lagrange{RefTriangle, 1}()
cv, fv = setup_vals(ip)

(Ferrite.CellValues{Ferrite.FunctionValues{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Vec{2, Float64}}, Nothing, Nothing}, Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}, Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}, Vector{Float64}}(Ferrite.FunctionValues{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Vec{2, Float64}}, Nothing, Nothing}(Ferrite.Lagrange{Ferrite.RefTriangle, 1}(), [0.16666666666667 0.16666666666667 0.66666666666667; 0.16666666666667 0.66666666666667 0.16666666666667; 0.6666666666666601 0.16666666666666008 0.16666666666666005], [0.16666666666667 0.16666666666667 0.66666666666667; 0.16666666666667 0.66666666666667 0.16666666666667; 0.6666666666666601 0.16666666666666008 0.16666666666666005], Tensors.Vec{2, Float64}[[NaN, NaN] [NaN, NaN] [NaN, NaN]; [NaN, NaN] [NaN, NaN] [NaN, NaN]; [NaN, NaN] [NaN, NaN] [NaN, NaN]], Tensors.Vec{2, Float64}[[1.0, 0.0] [1.0, 0.0] [1.0, 0.0]; [0.0, 1.0] [0.0, 1.0] [0.0, 1.0]; [-1.0, -1.0] [-1.0, -1.0] [-1.0, -1.0]], nothing, nothing), Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}(Ferrite.Lagrange{Ferrite.RefTriangle, 1}(), [0.16666666666667 0.16666666666667 0.66666666666667; 0.16666666666667 0.66666666666667 0.16666666666667; 0.6666666666666601 0.16666666666666008 0.16666666666666005], Tensors.Vec{2, Float64}[[1.0, 0.0] [1.0, 0.0] [1.0, 0.0]; [0.0, 1.0] [0.0, 1.0] [0.0, 1.0]; [-1.0, -1.0] [-1.0, -1.0] [-1.0, -1.0]], nothing), Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}([0.166666666666665, 0.166666666666665, 0.166666666666665], Tensors.Vec{2, Float64}[[0.16666666666667, 0.16666666666667], [0.16666666666667, 0.66666666666667], [0.66666666666667, 0.16666666666667]]), [NaN, NaN, NaN]), Ferrite.FacetValues{Ferrite.FunctionValues{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Vec{2, Float64}}, Nothing, Nothing}, Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}, Ferrite.FacetQuadratureRule{Ferrite.RefTriangle, Vector{Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}}}, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}, Vector{Ferrite.FunctionValues{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Vec{2, Float64}}, Nothing, Nothing}}, Vector{Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}}}(Ferrite.FunctionValues{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Vec{2, Float64}}, Nothing, Nothing}[Ferrite.FunctionValues{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Vec{2, Float64}}, Nothing, Nothing}(Ferrite.Lagrange{Ferrite.RefTriangle, 1}(), [0.7886751345948129 0.21132486540518713; 0.21132486540518713 0.7886751345948129; 0.0 0.0], [0.7886751345948129 0.21132486540518713; 0.21132486540518713 0.7886751345948129; 0.0 0.0], Tensors.Vec{2, Float64}[[NaN, NaN] [NaN, NaN]; [NaN, NaN] [NaN, NaN]; [NaN, NaN] [NaN, NaN]], Tensors.Vec{2, Float64}[[1.0, 0.0] [1.0, 0.0]; [0.0, 1.0] [0.0, 1.0]; [-1.0, -1.0] [-1.0, -1.0]], nothing, nothing), Ferrite.FunctionValues{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Vec{2, Float64}}, Nothing, Nothing}(Ferrite.Lagrange{Ferrite.RefTriangle, 1}(), [0.0 0.0; 0.7886751345948129 0.21132486540518713; 0.21132486540518713 0.7886751345948129], [0.0 0.0; 0.7886751345948129 0.21132486540518713; 0.21132486540518713 0.7886751345948129], Tensors.Vec{2, Float64}[[NaN, NaN] [NaN, NaN]; [NaN, NaN] [NaN, NaN]; [NaN, NaN] [NaN, NaN]], Tensors.Vec{2, Float64}[[1.0, 0.0] [1.0, 0.0]; [0.0, 1.0] [0.0, 1.0]; [-1.0, -1.0] [-1.0, -1.0]], nothing, nothing), Ferrite.FunctionValues{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Vec{2, Float64}}, Nothing, Nothing}(Ferrite.Lagrange{Ferrite.RefTriangle, 1}(), [0.21132486540518713 0.7886751345948129; 0.0 0.0; 0.7886751345948129 0.21132486540518713], [0.21132486540518713 0.7886751345948129; 0.0 0.0; 0.7886751345948129 0.21132486540518713], Tensors.Vec{2, Float64}[[NaN, NaN] [NaN, NaN]; [NaN, NaN] [NaN, NaN]; [NaN, NaN] [NaN, NaN]], Tensors.Vec{2, Float64}[[1.0, 0.0] [1.0, 0.0]; [0.0, 1.0] [0.0, 1.0]; [-1.0, -1.0] [-1.0, -1.0]], nothing, nothing)], Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}[Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}(Ferrite.Lagrange{Ferrite.RefTriangle, 1}(), [0.7886751345948129 0.21132486540518713; 0.21132486540518713 0.7886751345948129; 0.0 0.0], Tensors.Vec{2, Float64}[[1.0, 0.0] [1.0, 0.0]; [0.0, 1.0] [0.0, 1.0]; [-1.0, -1.0] [-1.0, -1.0]], nothing), Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}(Ferrite.Lagrange{Ferrite.RefTriangle, 1}(), [0.0 0.0; 0.7886751345948129 0.21132486540518713; 0.21132486540518713 0.7886751345948129], Tensors.Vec{2, Float64}[[1.0, 0.0] [1.0, 0.0]; [0.0, 1.0] [0.0, 1.0]; [-1.0, -1.0] [-1.0, -1.0]], nothing), Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}(Ferrite.Lagrange{Ferrite.RefTriangle, 1}(), [0.21132486540518713 0.7886751345948129; 0.0 0.0; 0.7886751345948129 0.21132486540518713], Tensors.Vec{2, Float64}[[1.0, 0.0] [1.0, 0.0]; [0.0, 1.0] [0.0, 1.0]; [-1.0, -1.0] [-1.0, -1.0]], nothing)], Ferrite.FacetQuadratureRule{Ferrite.RefTriangle, Vector{Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}}}(Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}[Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}([0.4999999999999999, 0.4999999999999999], Tensors.Vec{2, Float64}[[0.7886751345948129, 0.21132486540518713], [0.21132486540518713, 0.7886751345948129]]), Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}([0.4999999999999999, 0.4999999999999999], Tensors.Vec{2, Float64}[[0.0, 0.7886751345948129], [0.0, 0.21132486540518713]]), Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}([0.4999999999999999, 0.4999999999999999], Tensors.Vec{2, Float64}[[0.21132486540518713, 0.0], [0.7886751345948129, 0.0]])]), [NaN, NaN], Tensors.Vec{2, Float64}[[NaN, NaN], [NaN, NaN]], 1))

DofHandler handles numbering and distribution of degrees of freedom (dofs) of our solution.

dh = setup_dh(grid, ip)

Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}
  Fields:
    :u, Ferrite.Lagrange{Ferrite.RefTriangle, 1}()
  Dofs per cell: 3
  Total dofs: 1574

We use setup_bcs to construct our ConstraintHandler which contains the information of the boundary conditions. In this tutorial, we impose the Dirichelt boundary condition on the lower boundary and the periodic boundary condition on the left and right boundaries.

cst = setup_bcs(dh; period=2π)

ConstraintHandler:
  BCs:
    Field: u, Components: [1]

Assembling the stiffness matrix

We allocate the stiffness matrix $A$, the TBC matrix $F$ and the load vector $f$. Due to the DtN map, we need to extract the dofs associated to the artificial boundary.

# Extract dofs on the "top" boundary
top = getfacetset(grid, "top")
dofsDtN = dofs_on_dtn(dh, :u, top)

64-element Vector{Int64}:
   10
  109
  111
  122
  143
  144
  166
  191
  217
  233
    ⋮
 1494
 1497
 1501
 1508
 1538
 1542
 1543
 1565
 1566

Then we create the sparse pattern of $A$ and $F$. For more details, please refer to allocate_stiff_matrix.

# Allocate the stiffness matrix and the load vector
A = allocate_stiff_matrix(dh, cst, dofsDtN)
F = allocate_stiff_matrix(dh, cst, dofsDtN)
f = zeros(ComplexF64, ndofs(dh));

# Assemble the load vector f
f = assemble_load(fv, dh, top, f, inc, height)

# Assemble the matrix A
A = assemble_A(cv, dh, A, inc)

# Assemble the TBC matrix
F = assemble_tbc(fv, dh, inc, top, F, N, dofsDtN)

1574×1574 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 14666 stored entries:
⎡⠿⣧⣿⣿⣣⣩⣟⣜⠈⠀⠄⠔⠠⠐⠴⠴⠈⣰⠤⠀⣅⡂⠀⠀⢠⢘⢐⢣⢂⡒⡆⣟⣿⠁⢊⢩⢭⠸⢹⡊⎤
⎢⣿⣿⣿⣿⣿⣿⣧⢷⠶⠁⠄⠁⠈⡠⢓⠁⠨⢃⠁⠠⣤⡓⠁⠈⢠⡀⢰⡆⠡⢋⡀⣿⣿⡇⣸⢼⣶⣺⣾⣾⎥
⎢⡍⣺⣿⣿⣿⣿⡿⣿⢟⡁⠂⡡⠀⡤⡄⡉⡶⢈⢅⠈⣄⡢⢸⠀⠉⠩⡸⠐⠑⢀⠐⣭⣯⡅⢬⢨⣽⢨⣿⣯⎥
⎢⣛⢽⢭⣟⣿⣯⢻⣶⣆⠇⡀⠀⠈⠬⡍⡻⢦⢧⠏⠘⢲⠄⠐⠈⢐⢀⡀⢈⠀⣖⡓⣭⢻⠁⠈⠬⢯⠬⢽⣫⎥
⎢⠂⠀⠜⠃⠟⠱⠬⠝⠿⣧⣤⣝⡯⢘⣐⣎⢵⢿⣱⠽⡯⠷⠤⠬⢐⠄⠀⡃⠃⢁⣍⠤⢑⠀⠐⠀⠘⠯⠩⠢⎥
⎢⢀⠅⠄⠁⠌⡠⠀⠈⣄⢿⠿⣧⣿⣸⣦⠶⠚⣢⣦⣞⡌⠁⠈⢠⡧⠇⢠⠽⡆⢈⣨⠀⠈⠋⡄⠶⠤⠂⠸⢠⎥
⎢⢀⠂⠂⡠⠀⡤⡂⡄⣋⢋⣛⣻⣿⣿⣿⣽⡞⣫⣯⣾⡴⢒⡐⣀⠢⠽⠀⠀⠁⡑⢋⠆⠀⠞⠹⣔⢑⡰⠁⢀⎥
⎢⢐⡇⠝⠐⡄⠩⣧⡩⡰⢼⢨⡟⣟⣿⢿⣷⣿⣷⡿⣇⠖⠰⠂⠀⡨⢵⠀⣢⠨⣟⣥⠗⢸⠁⠰⣍⢣⠀⠀⠹⎥
⎢⢂⣠⠦⢂⡘⢋⠬⣗⣵⣗⠺⣠⡾⣩⢿⣿⣿⣿⣿⣿⣦⡮⢄⠌⣰⣨⠨⡊⢀⡟⠘⠀⠐⢬⡫⡋⣀⣃⠒⠀⎥
⎢⠀⠃⠁⡀⡁⠑⣋⠁⣕⡞⣨⢿⣫⣿⠿⢯⣿⣿⣿⣿⣧⡾⣷⡛⠲⣀⠂⠲⠨⡌⠌⠄⠀⢰⡭⣿⠐⡤⣿⠈⎥
⎢⠡⠹⢤⠻⠠⡹⠘⠖⢯⡏⠆⠉⢰⢋⢘⡁⡨⡿⣩⡿⢿⣷⣿⣧⣾⣸⣎⢪⠐⣼⣯⠀⠰⡌⣽⡍⢉⡟⠉⠊⎥
⎢⠀⠀⡁⠀⠒⠒⡐⠀⡀⡇⠂⣀⠐⢨⠈⠀⡀⠕⣽⠻⠿⣿⣿⣿⣿⣻⣟⢹⡀⠣⣯⠀⠨⣟⣉⠁⠉⡁⠗⠂⎥
⎢⣀⢒⠀⠲⡇⡀⠐⢐⠐⠔⠭⠏⣌⡆⢆⣎⡐⣺⠘⢢⣚⣻⣿⣻⣿⣿⣧⣿⣶⣽⡟⠀⠠⠿⢟⢉⣖⣁⡁⠆⎥
⎢⠴⣐⠰⠶⢒⠊⡀⢈⠤⠠⣄⡖⠀⠀⠠⣠⡢⠢⢨⡀⡪⣙⣟⣙⣭⣿⣿⣿⣿⣿⣿⠀⠈⡋⡎⡔⢱⢎⠷⠁⎥
⎢⢨⠰⡥⢂⠑⢀⢠⢤⠍⢀⡈⢉⢅⠠⣦⢦⣤⠴⡂⠦⣐⣤⠤⡈⣜⣿⣿⣿⣿⣿⣿⡄⢰⣲⡿⡆⡀⠠⠆⠄⎥
⎢⣬⢭⣤⣬⡔⣤⡝⣬⠃⡝⠂⠚⠫⠔⢥⠟⠒⠀⠂⠅⠋⠛⠋⠛⠛⠉⠛⠛⠛⠿⠟⣥⣤⡿⢳⢰⣤⢠⣤⣤⎥
⎢⠟⠛⠿⠿⠏⠿⠟⠒⠑⠐⡦⠀⣠⠄⠖⠒⡐⣄⢀⣀⡐⠦⣦⢦⣤⡆⡦⠠⢰⣲⣤⡿⠿⣧⣼⣾⣿⢼⡿⠻⎥
⎢⡎⣐⣒⣞⡂⣓⡂⡄⠐⠀⢠⡍⢓⢦⡔⢦⡯⠪⣧⣯⡗⠿⠇⠘⡟⢑⢊⠭⠻⠯⢙⣒⣲⣿⣿⣿⣿⣿⣷⣐⎥
⎢⣃⡓⣸⣻⡓⣛⡋⡗⡶⡄⠠⠃⢑⡰⠉⠒⠤⢸⠐⡤⣧⠴⠇⠠⠜⢹⡱⢖⠀⡈⠀⣛⣛⣟⣿⣿⣿⣿⣿⣘⎥
⎣⡳⠲⣺⣿⡿⣿⡷⣳⠣⡂⠒⣂⠁⢀⣄⡀⠘⠀⡛⠛⡣⠀⠹⠁⠡⠌⠝⠃⠈⠅⠀⣿⣿⡋⢙⢻⣛⢻⣿⣿⎦

Be careful! Calculations between two sparse matrices may destory the structure of sparse matrices. Then imposing the boundary condition may fail. So we use sub_preserve_structure to subtract $F$ from $A$ with keeping the structure of $A$.

# Add the TBC matrix to A = A - F
A = sub_preserve_structure(A, F)

# Impose the boundary condition
apply!(A, f, cst)

# Solve the linear system
u = A\f

# Apply the boundary condition again
apply!(u, cst)

1574-element Vector{ComplexF64}:
   3.158683268810116e-5 - 0.15317244858344584im
  1.7158948105089617e-5 - 0.09651732153593437im
    3.79994723326924e-5 - 0.20674399138797583im
 -2.5584384023723487e-5 - 0.09624018671542418im
  -5.393305399234056e-5 - 0.20446036151544966im
 -4.0360981259559305e-5 - 0.15296601009255503im
 0.00010349732894413891 - 0.18883980316651755im
   4.906580304024992e-5 - 0.08636434114583684im
  4.5931041288458144e-5 - 0.08636330998468002im
   5.562684550090564e-5 - 1.6827193257988013im
                        ⋮
  -0.010124949146285431 - 1.682889105761293im
   0.001299890206179951 - 1.6419189991224463im
                    0.0 + 0.0im
                    0.0 + 0.0im
   5.812820521565083e-5 - 0.23554401819482904im
  -3.317216149964981e-5 - 0.14140319196319232im
 0.00016500884652190343 - 0.22843204821949303im
  0.0003688976164628844 - 0.49466242858564047im
  0.0002486683289360293 - 0.32632210400148437im

Export the VTK

Finally, we write the real part and imaginary part of the solution into vtk files which can be visualized by Paraview.

VTKGridFile("real_u", grid) do vtk
   write_solution(vtk, dh, real.(u))
end;

VTKGridFile("imag_u", grid) do vtk
    write_solution(vtk, dh, imag.(u))
end;

This page was generated using Literate.jl.