Add interpolation support for Helmholtz RHS

This commit enables support of primal piecewise linear and primal
piecewise constant RHSs

TODO: Probably, either an extension or a rewrite of DofInterpolate is
advisable for the LagrangeBasis (for example, also dealing with higher-order
functions and dual functions)

src/helmholtz3d/hh3dexc.jl

@@ -11,6 +11,10 @@ function (f::HH3DPlaneWave)(r)
     return a * exp(-γ*dot(d,r))
 end
 
+function (f::HH3DPlaneWave)(r::CompScienceMeshes.MeshPointNM)
+    return f(cartesian(r))
+end
+
 scalartype(f::HH3DPlaneWave{P,K,T}) where {P,K,T} = promote_type(eltype(P), K, T)
 
 """
@@ -32,6 +36,9 @@ function (f::HH3DLinearPotential)(r)
     return a * dot(d, r)
 end
 
+function (f::HH3DLinearPotential)(r::CompScienceMeshes.MeshPointNM)
+    return f(cartesian(r))
+end
 struct gradHH3DLinearPotential{T,P}
     direction::P
     amplitude::T
@@ -44,6 +51,11 @@ function (f::gradHH3DLinearPotential)(r)
     return a * d
 end
 
+function (f::gradHH3DLinearPotential)(r::CompScienceMeshes.MeshPointNM)
+    r = cartesian(mp)
+    return dot(normal(mp), f(r))
+end
+
 function grad(m::HH3DLinearPotential)
     return gradHH3DLinearPotential(m.direction, m.amplitude)
 end
@@ -66,6 +78,10 @@ end
 
 scalartype(x::HH3DMonopole{P,K,T}) where {P,K,T} = promote_type(eltype(P), K, T)
 
+function (f::HH3DMonopole)(r::CompScienceMeshes.MeshPointNM)
+    return f(cartesian(r))
+end
+
 function (f::HH3DMonopole)(r)
     γ = f.gamma
     p = f.position
@@ -92,6 +108,11 @@ function (f::gradHH3DMonopole)(r)
     return -a * vecR * exp(-γ * R) / R^2 * (γ + 1 / R)
 end
 
+function (f::gradHH3DMonopole)(mp::CompScienceMeshes.MeshPointNM)
+    r = cartesian(mp)
+    return dot(normal(mp), f(r))
+end
+
 function grad(m::HH3DMonopole)
     return gradHH3DMonopole(m.position, m.gamma, m.amplitude)
 end

src/interpolation.jl

@@ -1,9 +1,11 @@
 
 function DofInterpolate(basis::LagrangeBasis, field) 
 
+    T = promote_type(scalartype(basis), scalartype(field))
+
     num_bfs = numfunctions(basis)
 
-    res = Vector(undef, num_bfs)
+    res = Vector{T}(undef, num_bfs)
 
     for  b in 1 : num_bfs
         bfs = basis.fns[b]
@@ -33,11 +35,43 @@ function DofInterpolate(basis::LagrangeBasis, field)
     return res
 end
 
+### Piecewise constant elements require separate treatment
+### TODO: Probably a rewrite is advisable to also take into account
+### dual elements properly.
+function DofInterpolate(basis::LagrangeBasis{0,-1,M,T,NF,P}, field) where {M, T, NF, P}
+
+    TT = promote_type(scalartype(basis), scalartype(field))
+
+    num_bfs = numfunctions(basis)
+
+    res = Vector{TT}(undef, num_bfs)
+
+    for  b in 1 : num_bfs
+        bfs = basis.fns[b]
+
+        basis.pos[b]
+
+        shape = bfs[1]
+
+        cellid = shape.cellid
+
+        tria = chart(basis.geo, cellid)
+
+        v = neighborhood(tria, [1/3, 1/3])
+
+        res[b] = field(v)
+    end
+
+    return res
+end
+
 function DofInterpolate(basis::RTBasis, field)
 
+    T = promote_type(scalartype(basis), scalartype(field))
+
     num_bfs = numfunctions(basis)
 
-    res = Vector(undef, num_bfs)
+    res = Vector{T}(undef, num_bfs)
 
     for b in 1 : num_bfs
 
@@ -76,9 +110,11 @@ end
 
 function DofInterpolate(basis::BDMBasis, field)
 
+    T = promote_type(scalartype(basis), scalartype(field))
+
     num_bfs = numfunctions(basis)
 
-    res = Vector(undef, num_bfs)
+    res = Vector{T}(undef, num_bfs)
 
     for b in 1 : num_bfs
 
@@ -129,9 +165,11 @@ end
 
 function DofInterpolate(basis::NDLCDBasis, field)
 
+    T = promote_type(scalartype(basis), scalartype(field))
+
     num_bfs = numfunctions(basis)
 
-    res = Vector(undef, num_bfs)
+    res = Vector{T}(undef, num_bfs)
 
     for b in 1 : num_bfs
 
@@ -173,9 +211,11 @@ end
 
 function DofInterpolate(basis::BEAST.BDM3DBasis, field)
 
+    T = promote_type(scalartype(basis), scalartype(field))
+
     num_bfs = numfunctions(basis)
 
-    res = Vector(undef, num_bfs)
+    res = Vector{T}(undef, num_bfs)
 
     for b in 1 : num_bfs
 

test/test_interpolation.jl

@@ -0,0 +1,28 @@
+using BEAST
+using CompScienceMeshes
+using Test
+
+Γ =  meshrectangle(1.0, 1.0, 0.4)
+
+XN = lagrangecxd0(Γ) # lagrangec0d1(Γₑ₁; dirichlet=true) # Dirichlet=true -> no boundary function
+XD = lagrangec0d1(Γ)
+
+mp = Helmholtz3D.monopole(position=SVector(0.0, 0.0, 3.0))
+gradmp = grad(mp)
+
+φ = DofInterpolate(XD, mp)
+
+@test eltype(φ) == Float64
+
+for i in eachindex(φ)
+    @test mp(XD.pos[i]) == φ[i] 
+end
+
+σ = DofInterpolate(XN, gradmp)
+
+@test eltype(σ) == Float64
+
+for i in eachindex(σ)
+    # Orientation of meshrectangle is in -z direction
+    @test dot(SVector(0.0, 0.0, -1.0), gradmp(XN.pos[i])) == σ[i]
+end