Skip to content

Example implementation of 1D curvilinear elements #180

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
sbadrian wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Example implementation of 1D curvilinear elements #180

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 10 additions & 3 deletions src/bases/lagrange.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1194,10 +1194,17 @@ function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,2}; order)
return LagrangeBasis{order,-1,NF}(mesh, fns, pos)
end

function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order)
function lagrangecx(mesh::CompScienceMeshes.AbstractMesh; order)

T = coordtype(mesh)
NF = binomial(2+order, 2)

if dimension(mesh) == 1
NF = 1 + order
elseif dimension(mesh) == 2
NF = binomial(2+order, 2)
else
error("We do not yet support tetrahedron")
end
P = vertextype(mesh)
S = Shape{T}

Expand Down Expand Up @@ -1320,7 +1327,7 @@ end



function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order)
function lagrangec0(mesh::CompScienceMeshes.Mesh{<:Any,3}; order)

T = coordtype(mesh)
atol = sqrt(eps(T))
Expand Down
4 changes: 2 additions & 2 deletions src/helmholtz2d/hh2dops.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -230,8 +230,8 @@ function quaddata(op::HelmholtzOperator2D,
end

function quadrule(op::HelmholtzOperator2D, g::LagrangeRefSpace, f::LagrangeRefSpace,
i, τ::CompScienceMeshes.Simplex{<:Any, 1},
j, σ::CompScienceMeshes.Simplex{<:Any, 1},
i, τ::CompScienceMeshes.AbstractSimplex{<:Any, 1},
j, σ::CompScienceMeshes.AbstractSimplex{<:Any, 1},
qd, qs::DoubleNumSauterQstrat)

hits = _numhits(τ, σ)
Expand Down
2 changes: 1 addition & 1 deletion src/identityop.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ end
defaultquadstrat(::LocalOperator, ::LagrangeRefSpace{T,D1,2}, ::LagrangeRefSpace{T,D2,2}) where {T,D1,D2} = SingleNumQStrat(6)
function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,2} where {T,Deg},
f::LagrangeRefSpace, tels::Vector, bels::Vector, qs::SingleNumQStrat)
U = typeof(tels[1].volume)
U = typeof(volume(tels[1])) # works for straight and curvilinear charts
u, w = legendre(qs.quad_rule, 0.0, 1.0)
qd = [(w[i],u[i]) for i in eachindex(w)]
A = _alloc_workspace(Tuple{U,U}.(qd), g, f, tels, bels)
Expand Down
4 changes: 2 additions & 2 deletions src/quadrature/doublenumsauterqstrat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -228,8 +228,8 @@ function _numhits(τ, σ)
hits = 0
dtol = 1.0e3 * eps(T)
dmin2 = floatmax(T)
for t in vertices(τ)
for s in vertices(σ)
for t in nodes(τ)
for s in nodes(σ)
d2 = LinearAlgebra.norm_sqr(t-s)
d = norm(t-s)
dmin2 = min(dmin2, d2)
Expand Down
109 changes: 35 additions & 74 deletions test/test_hh2d_mie_higher_order.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,8 @@ using LinearAlgebra
using Test
using SpecialFunctions


# We check only
for order = [2; 3]
# We check only order 5 (implies four internal control vertices)
for order = [5]

ε0 = 8.854187821e-12
μ0 = 4π*1e-7
Expand All @@ -21,7 +20,7 @@ for order = [2; 3]
h = λ/10
a = 1.0 # radius of the scatterer

circle = CompScienceMeshes.meshcircle(a, h)
circle = CompScienceMeshes.gmshcircle(a, h; order=order)
X0 = lagrangecx(circle, order=order)
X1 = lagrangec0(circle, order=order)

Expand All @@ -31,20 +30,15 @@ for order = [2; 3]

## Analytical solutions

dbesselj(n,x) = (besselj(n-1, x) - besselj(n+1, x)) / 2
dhankelh2(n,x) = (hankelh2(n-1, x) - hankelh2(n+1, x)) / 2
dbesselj(n,x) = (besselj(n-1, x) - besselj(n+1, x)) / 2
dhankelh2(n,x) = (hankelh2(n-1, x) - hankelh2(n+1, x)) / 2

cart2polar(x,y) = SVector(sqrt(x^2 + y^2), atan(y, x))




## Analytical Solutions

# TM planewave solution (based on Sec 6.4.2, Jin, Theory and Computation of Electromagnetic Fields])
function TM_pec_planewave(E0, k, a, ρ, φ)
n = 0

# Jin (6.4.11)
a_n(n) = -(1.0im)^(-n)*(besselj(n, k*a) / hankelh2(n, k*a))
valEz(n) = a_n(n) * hankelh2(n, k*ρ)*exp(im*n*φ)
Expand All @@ -55,33 +49,22 @@ for order = [2; 3]
retHφ = valHφ(0)
retHρ = valHρ(0)

while true
tol = eps(eltype(real(retEz))) * 1e1
n = 0
while n < 100
n += 1
bufEz = valEz(n) + valEz(-n)
if abs(bufEz) >= abs(retEz) * eps(eltype(real(retEz))) * 1e3
retEz += bufEz
else
break
end
end

n=0
while true
n += 1
bufEz = valEz(n) + valEz(-n)
bufHφ = valHφ(n) + valHφ(-n)
if abs(bufHφ) >= abs(retHφ) * eps(eltype(real(retHφ))) * 1e3
retHφ += bufHφ
else
break
end
end

n=0
while true
n += 1
bufHρ = valHρ(n) + valHρ(-n)
if abs(bufHρ) >= abs(retHρ) * eps(eltype(real(retHρ))) * 1e3
retHρ += bufHρ

if abs(bufEz) >= abs(retEz) * tol ||
abs(bufHφ) >= abs(retHφ) * tol ||
abs(bufHρ) >= abs(retHρ) * tol

bufEz != NaN ? retEz += bufEz : ~
bufHφ != NaN ? retHφ += bufHφ : ~
bufHρ != NaN ? retHρ += bufHρ : ~
else
break
end
Expand All @@ -98,10 +81,8 @@ for order = [2; 3]
return [TM_pec_planewave(E0, k, a, cart2polar(p[1], p[2])...)[2] for p in pts]
end


# TE planewave solution
function TE_pec_planewave(H0, k, a, ρ, φ)
n = 0

# Jin (6.4.19)
b(n) = -(1.0im)^(-n)*dbesselj(n,k*a)/dhankelh2(n,k*a)
Expand All @@ -115,33 +96,22 @@ for order = [2; 3]
retEφ = valEφ(0)
retEρ = valEρ(0)

while true
n += 1
bufHz = valHz(n) + valHz(-n)
if abs(bufHz) >= abs(retHz) * eps(eltype(real(retHz))) * 1e3
retHz += bufHz
else
break
end
end

tol = eps(eltype(real(retHz))) * 1e1
n = 0
while true
n += 1
bufEφ = valEφ(n) + valEφ(-n)
if abs(bufEφ) >= abs(retEφ) * eps(eltype(real(retEφ))) * 1e3
retEφ += bufEφ
else
break
end
end

n = 0
while true
while n < 100
n += 1
bufHz = valHz(n) + valHz(-n)
bufEφ = valEφ(n) + valEφ(-n)
bufEρ = valEρ(n) + valEρ(-n)
if abs(bufEρ) >= abs(retEρ) * eps(eltype(real(retEρ))) * 1e3
retEρ += bufEρ

if abs(bufHz) >= abs(retHz) * tol ||
abs(bufEφ) >= abs(retEφ) * tol ||
abs(bufEρ) >= abs(retEρ) * tol

bufHz != NaN ? retHz += bufHz : ~
bufEφ != NaN ? retEφ += bufEφ : ~
bufEρ != NaN ? retEρ += bufEρ : ~
else
break
end
Expand All @@ -158,9 +128,6 @@ for order = [2; 3]
return [TE_pec_planewave(H0, k, a, cart2polar(p[1],p[2])...)[2] for p in pts]
end




## TM Scattering

# TM-EFIE system matrix
Expand All @@ -185,17 +152,13 @@ for order = [2; 3]
Ez_pw_sca_num = -potential(HH2DSingleLayerNear(𝒮), pts, j_TMEFIE_pw, X0; type=ComplexF64)
Ez_pw_sca_ana = TM_pec_planewave_E(E0, k, a, pts)

# Will not improve until we have curvilinear elements
# But must not become worse either
# TODO: Once curvilinear elements are support: refine order of functions and
# mesh in a lock step and update these values
Ez_pw_sca_bounds = [0.002; 0.002; 0.0014]
Ez_pw_sca_bounds = [0.002; 1.6e-8; 2.3e-9; 5.0e-13; 5.1e-13]
@test norm(Ez_pw_sca_ana - Ez_pw_sca_num) ./ norm(Ez_pw_sca_ana) < Ez_pw_sca_bounds[order]

Ht_pw_sca_num = -potential(HH2DDoubleLayerTransposedNear(𝒟ᵀ), pts, j_TMEFIE_pw, X0; type=SVector{2,ComplexF64})
Ht_pw_sca_ana = TM_pec_planewave_H(E0, k, a, pts)

Ht_pw_sca_bounds = [0.002; 0.002; 0.0015]
Ht_pw_sca_bounds = [0.002; 1.6e-8; 2.3e-9; 5e-13; 5.1e-13]
@test norm(Ref(ẑ) .× Ht_pw_sca_num - Ht_pw_sca_ana) ./ norm(Ht_pw_sca_ana) < Ht_pw_sca_bounds[order]

# TM-MFIE
Expand All @@ -204,11 +167,9 @@ for order = [2; 3]

j_TMMFIE_pw = M_TMMFIE \ ht_pw_inc

j_pw_sca_bounds = [0.05; 0.001; 0.0009]
j_pw_sca_bounds = [0.05; 0.0002; 5e-6; 4e-7; 6.3e-9]
@test norm(j_TMMFIE_pw - j_TMEFIE_pw) / norm(j_TMEFIE_pw) < j_pw_sca_bounds[order]



## TE Scattering

I = 1.0
Expand Down Expand Up @@ -238,7 +199,7 @@ for order = [2; 3]
Hz_pw_sca_num = potential(HH2DDoubleLayerNear(𝒟), pts, j_TEMFIE_pw, X1; type=ComplexF64)
Hz_pw_sca_ana = TE_pec_planewave_H(H0, k, a, pts)

Hz_pw_sca_bounds = [0.0015; 0.0015; 0.0015]
Hz_pw_sca_bounds = [0.0015; 7.2e-7; 2.4e-9; 1.8e-10; 3.2e-10]
@test norm(Hz_pw_sca_num - Hz_pw_sca_ana) /norm(Hz_pw_sca_ana) < Hz_pw_sca_bounds[order]

Et_pw_inc = 1 / (im * ω * ε0) * curl(Hz_pw_inc)
Expand All @@ -249,14 +210,14 @@ for order = [2; 3]

j_TEEFIE_pw = M_TEEFIE \ et_pw_inc

j_pw_sca_bounds = [0.007; 0.0002; 1.19e-5]
j_pw_sca_bounds = [0.007; 0.0002; 1e-5; 2.3e-7; 1.2e-8]
@test norm(j_TEEFIE_pw - j_TEMFIE_pw)/norm(j_TEMFIE_pw) < j_pw_sca_bounds[order]

Et_pw_sca_num = -potential(HH2DHyperSingularNear(𝒩), pts, j_TEEFIE_pw, X1; type=SVector{2, ComplexF64})
Et_pw_sca_ana = TE_pec_planewave_E(H0, k, a, pts)

# We compute the scattered Ez component (scalar)
Et_pw_sca_bounds = [0.003; 0.0016; 0.0016]
Et_pw_sca_bounds = [0.003; 1.6e-5; 1.9e-8; 3.3e-11; 5.3e-13]
@test norm(Ref(ẑ) .× Et_pw_sca_num - Et_pw_sca_ana) / norm(Et_pw_sca_ana) < Et_pw_sca_bounds[order]

end
Loading
Loading