Merge_beast

Project.toml

@@ -16,10 +16,12 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LiftedMaps = "d22a30c1-52ac-4762-a8c9-5838452405e0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 NestedUnitRanges = "032820ab-dc03-4b49-91f4-7d58d4da98b3"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SauterSchwab3D = "0a13313b-1c00-422e-8263-562364ed9544"
 SauterSchwabQuadrature = "535c7bfe-2023-5c1d-b712-654ef9d93a38"

examples/ex_globalmultitrace.jl

@@ -15,6 +15,7 @@ h = 0.125
 Γ2 =  -Mesh([point(-x,y,z) for (x,y,z) in vertices(Γ1)], deepcopy(cells(Γ1)))
 Γ = [Γ1, Γ2]
 
+
 # Beyond this point the problem description is completely independent
 # of the number of domains and their relative positioning
 

examples/pmchwt.jl

@@ -24,16 +24,18 @@ end
 μ0 = 4π*1e-7
 c = 1/√(ϵ0*μ0)
 
-λ = 2.9979563769321627
-ω = 2π*c/λ
+# λ = 2.9979563769321627
+# ω = 2π*c/λ
+ω = 10.0^8*2*pi
 
-Ω = CompScienceMeshes.tetmeshsphere(λ,0.1*λ)
+# Ω = CompScienceMeshes.tetmeshsphere(λ,0.1*λ)
+Ω = CompScienceMeshes.tetmeshsphere(1.0,0.2)
 Γ = boundary(Ω)
 X = raviartthomas(Γ)
 @show numfunctions(X)
 
 ϵr = 2.0
-μr = 1.0
+μr = 3.0
 
 κ, η = ω/c, √(μ0/ϵ0)
 κ′, η′ = κ*√(ϵr*μr), η*√(μr/ϵr)
@@ -61,7 +63,7 @@ pmchwt = @discretise(
          (K+K′)[l,j] + (α*T+α′*T′)[l,m] == -e[k] - h[l],
     j∈X, m∈X, k∈X, l∈X)
 
-u = solve(pmchwt)
+u = BEAST.solve(pmchwt)
 
 Θ, Φ = range(0.0,stop=2π,length=100), 0.0
 ffpoints = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]
@@ -97,10 +99,10 @@ Plotly.plot(patch(Γ, norm.(fcrm)))
 
 
 
-Z = range(-6,6,length=200)
-Y = range(-4,4,length=200)
-nfpoints = [point(0,y,z) for z in Z, y in Y]
-
+Zz = range(-1,1,length=130)
+Yy = range(-1,1,length=100)
+nfpoints = [point(z,0,y) for z in Zz, y in Yy]
+nfpoints = pts
 import Base.Threads: @spawn
 task1 = @spawn nearfield(u[m],u[j],X,X,κ,η,nfpoints,E,H)
 task2 = @spawn nearfield(-u[m],-u[j],X,X,κ′,η′,nfpoints)
@@ -116,11 +118,11 @@ Plots.contour(real.(getindex.(H_tot,2)))
 
 Plots.heatmap(Z, Y, clamp.(real.(getindex.(E_tot,1)),-1.5,1.5))
 Plots.heatmap(Z, Y, clamp.(imag.(getindex.(E_tot,1)),-1.5,1.5))
-Plots.heatmap(Z, Y, real.(getindex.(H_tot,2)))
-Plots.heatmap(Z, Y, imag.(getindex.(H_tot,2)))
+Plots.heatmap(Yy, Zz, norm.(getindex.(H_tot,2)))
+Plots.heatmap(Y, Z, norm.(getindex.(H_tot,2)))
 
 Plots.plot(real.(getindex.(E_tot[:,51],1)))
 Plots.plot(real.(getindex.(H_tot[:,51],2)))
 
-
-
+Plots.heatmap(Zz,Yy,norm.(getindex.(-Bbb[1].+Bbb[3].-Bbb[2]/2,2))/μ0,clim=(0,0.02))
+Plots.heatmap(Zz,Yy,real(getindex.(H.(pts),2)))