Merge pull request #126 from bmergl/feature/lsvie_4

Add Lippmann-Schwinger volume integral equation (HH3D VIE)

examples/hh_lsvie.jl

@@ -0,0 +1,123 @@
+using BEAST
+using CompScienceMeshes
+using StaticArrays
+using LinearAlgebra
+using Test
+using Plots
+using SphericalScattering
+
+# Layered Dielectic Sphere
+# MoM (Lippmann Schwinger VIE) vs. Analytical Solution (SphericalScattering)
+
+
+
+# Environment
+ε1 = 1.0*ε0
+μ1 = μ0
+
+# Layered Dielectic Sphere
+ε2 = 5.0*ε0 # outer shell
+μ2 = μ0
+
+ε3 = 20.0*ε0
+μ3 = μ0
+
+ε4 = 7.0*ε0
+μ4 = μ0
+
+ε5 = 30.0*ε0 # inner core
+μ5 = μ0
+
+r = 1.0
+radii = SVector(0.2*r, 0.6*r, 0.8*r, r) # inner core stops at first radius
+filling = SVector(Medium(ε5, μ5), Medium(ε4, μ4), Medium(ε3, μ3), Medium(ε2, μ2))
+
+
+
+
+# Mesh, Basis
+h = 0.1
+mesh = CompScienceMeshes.tetmeshsphere(r,h)
+X = lagrangec0d1(mesh; dirichlet = false)
+@show numfunctions(X)
+
+bnd = boundary(mesh)
+strc = X -> strace(X, bnd)
+
+
+# VIE Operators
+function generate_tau(ε_env, radii, filling)
+
+    function tau(x::SVector{U,T}) where {U,T}
+        for (i, radius) in enumerate(radii)
+            norm(x) <= radius && (return T(1.0 - filling[i].ε/ε_env))
+        end
+        #return 0.0
+        error("Evaluated contrast outside the sphere!")
+    end
+
+    return tau
+end
+τ = generate_tau(ε1, radii, filling) #contrast function
+
+I = Identity()
+V = VIE.hhvolume(tau = τ, wavenumber = 0.0)
+B = VIE.hhboundary(tau = τ, wavenumber = 0.0)
+Y = VIE.hhvolumegradG(tau = τ, wavenumber = 0.0)
+
+
+# Exitation
+dirE = SVector(1.0, 0.0, 0.0)
+dirgradΦ = - dirE
+amp = 1.0
+Φ_inc = VIE.linearpotential(direction = dirgradΦ, amplitude = amp)
+
+
+# Assembly
+b = assemble(Φ_inc, X)
+
+Z_I = assemble(I, X, X)
+
+Z_V = assemble(V, X, X)
+Z_B = assemble(B, strc(X), X)
+
+Z_Y = assemble(Y, X, X)
+
+
+# System matrix
+Z_version1 = Z_I + Z_V + Z_B
+Z_version2 = Z_I + Z_Y
+
+
+# Solution of the linear system of equations
+u_version1 = Z_version1 \ Vector(b)
+u_version2 = Z_version2 \ Vector(b)
+
+
+
+## Observe the potential on a x-line in dirgradΦ direction
+
+# x-line at y0, z0 - Potential only inside the sphere mesh valid!
+y0 = 0.0
+z0 = 0.0
+x_range = range(0.0, stop = 1.0*r, length = 200)
+points_x = [point(x, y0, z0) for x in x_range]
+x = collect(x_range)
+
+# SphericalScattering solution on x-line
+sp = LayeredSphere(radii = radii, filling = filling)
+ex = UniformField(direction = dirE, amplitude = amp)
+Φ_x = real.(field(sp, ex, ScalarPotential(points_x)))
+
+# MoM solution on x-line
+Φ_MoM_version1_x = BEAST.grideval(points_x, u_version1, X, type = Float64)
+Φ_MoM_version2_x = BEAST.grideval(points_x, u_version2, X, type = Float64)
+
+# Plot
+plot(x, Φ_x, label = "SphericalScattering")
+plot!(x, Φ_MoM_version1_x, label = "MoM: S=I+B+V, boundary+volume")
+plot!(x, Φ_MoM_version2_x, label = "MoM: S=I+Y, gradgreen")
+xlims!(0.0, 1.0) # sphere center to radius r
+ylims!(-0.3, 0.0)
+title!("Potential Φ(x, y0, z0)")
+xlabel!("x")

src/bases/local/laglocal.jl

@@ -137,6 +137,23 @@ function strace(x::LagrangeRefSpace, cell, localid, face)
     Q
 end
 
+function strace(x::LagrangeRefSpace{T, 1, 4, 4}, cell, localid, face) where {T}
+
+    #T = scalartype(x)
+    t = zeros(T, 3, 4)
+    for (k,fvert) in enumerate(face.vertices)
+        for (l,cvert) in enumerate(cell.vertices)
+            nrm = norm(fvert - cvert)
+            if isapprox(nrm, 0, atol=sqrt(eps(T)))
+                t[k,l] = T(1.0)
+                break
+            end
+        end
+    end
+
+    return t
+end
+
 
 function restrict(refs::LagrangeRefSpace{T,0}, dom1, dom2) where T
     #Q = eye(T, numfunctions(refs))
@@ -287,4 +304,29 @@ const _dof_lag0perm_matrix = [
 
     @SMatrix[1]         # 6. {3,2,1}
 
-]
+]
+
+function (ϕ::LagrangeRefSpace{T, 1, 4, 4})(lag) where T
+
+    u, v, w = parametric(lag)
+
+    tu = tangents(lag, 1)
+    tv = tangents(lag, 2)
+    tw = tangents(lag, 3)
+
+    B = [tu tv tw]
+    A = inv(transpose(B))
+
+    # gradient in u,v,w (unit tetrahedron)
+    gr1=SVector{3, T}(1.0, 0.0, 0.0)
+    gr2=SVector{3, T}(0.0, 1.0, 0.0)
+    gr3=SVector{3, T}(0.0, 0.0, 1.0)
+    gr4=SVector{3, T}(-1.0, -1.0, -1.0)
+
+    return SVector((
+        (value = u, gradient = A*gr1),
+        (value = v, gradient = A*gr2),
+        (value = w, gradient = A*gr3),
+        (value = T(1.0)-u-v-w, gradient = A*gr4)
+    ))
+end

src/bases/trace.jl

@@ -139,6 +139,60 @@ function strace(X::Space, γ, dim1::Type{Val{2}})
     strace(X, Σ, fns)
 end
 
+function strace(X::Space, γ, dim1::Type{Val{3}})
+
+    x = refspace(X)
+    E, ad, P = assemblydata(X)
+    igeo = geometry(X)
+
+    @assert dimension(γ) == dimension(igeo)-1
+
+    ogeo = boundary(igeo)
+    on_target = overlap_gpredicate(γ)
+    ogeo = submesh(ogeo) do m,f
+        ch = chart(m,f)
+        on_target(ch)
+    end
+
+    D = connectivity(igeo, ogeo, abs)
+    rows, vals = rowvals(D), nonzeros(D)
+
+    T = scalartype(X)
+    S = Shape{T}
+    fns = [Vector{S}() for i in 1:numfunctions(X)]
+
+    for (p,el) in enumerate(E)
+
+        for (q,fc) in enumerate(faces(el))
+            on_target(fc) || continue
+
+            r = 0
+            for k in nzrange(D,P[p])
+                vals[k] == q && (r = rows[k]; break)
+            end
+            @assert r != 0
+
+            fc1 = chart(ogeo, r)
+            Q = strace(x, el, q, fc1)
+
+            for i in 1:size(Q,1)
+                for j in 1:size(Q,2)
+                    for (m,a) in ad[p,j]
+                        v = a*Q[i,j]
+                        isapprox(v,0,atol=sqrt(eps(T))) && continue
+                        push!(fns[m], Shape(r, i, v))
+                    end
+                end
+            end
+
+        end
+
+    end
+
+    strace(X, ogeo, fns)
+
+end
+
 """
 currently not working!
 """

src/volumeintegral/sauterschwab_ints.jl

@@ -71,6 +71,35 @@ function reorder_dof(space::LagrangeRefSpace{Float64,0,3,1},I)
     return SVector{1,Int64}(1),SVector{1,Int64}(1)
 end
 
+function reorder_dof(space::LagrangeRefSpace{T,1,3,3},I) where T
+    n = length(I)
+    K = zeros(MVector{3,Int64})
+    for i in 1:n
+        for j in 1:n
+            if I[j] == i
+                K[i] = j
+                break
+            end
+        end
+    end
+
+    return SVector(K),SVector{3,Int64}(1,1,1)
+end
+
+function reorder_dof(space::LagrangeRefSpace{T,1,4,4},I) where T
+    n = length(I)
+    K = zeros(MVector{4,Int64})
+    for i in 1:n
+        for j in 1:n
+            if I[j] == i
+                K[i] = j
+                break
+            end
+        end
+    end
+
+    return SVector(K),SVector{4,Int64}(1,1,1,1)
+end
 
 function momintegrals!(op::VIEOperator,
     test_local_space::RefSpace, trial_local_space::RefSpace,