VIE bugfix

src/bases/local/laglocal.jl

@@ -9,8 +9,9 @@ numfunctions(s::LagrangeRefSpace{T,1,3}, ch::CompScienceMeshes.ReferenceSimplex{
 numfunctions(s::LagrangeRefSpace{T,2,3}, ch::CompScienceMeshes.ReferenceSimplex{2}) where {T} = 6
 numfunctions(s::LagrangeRefSpace{T,Dg}, ch::CompScienceMeshes.ReferenceSimplex{D}) where {T,Dg,D} = binomial(D+Dg,Dg)
 
-valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} =
-        SVector{numfunctions(ref), Tuple{T,T}}
+# valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} =
+#         SVector{numfunctions(ref), Tuple{T,T}}
+valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} = T
 
 # Evaluate constant lagrange elements on anything
 (ϕ::LagrangeRefSpace{T,0})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),))

src/postproc/segcurrents.jl

@@ -29,21 +29,25 @@ function grideval(points, coeffs, basis; type=nothing)
 
     V = valuetype(refs, eltype(charts))
     T = promote_type(eltype(coeffs), eltype(V))
-    P = similar_type(V, T)
+    if !(V <: SVector)
+        P = T
+    else
+        P = similar_type(V, T)
+    end
 
-    type != nothing && (P = type)
+    type !== nothing && (P = type)
 
     values = zeros(P, size(points))
 
     chart_tree = BEAST.octree(charts)
     for (j,point) in enumerate(points)
         i = CompScienceMeshes.findchart(charts, chart_tree, point)
-        if i != nothing
+        if i !== nothing
             # @show i
             chart = charts[i]
             u = carttobary(chart, point)
             vals = refs(neighborhood(chart,u))
-            for r in 1 : numfunctions(refs)
+            for r in 1 : numfunctions(refs, domain(chart))
                 for (m,w) in ad[i, r]
                     values[j] += w * coeffs[m] * vals[r][1]
                 end

src/volumeintegral/sauterschwab_ints.jl

@@ -66,7 +66,7 @@ function reorder_dof(space::RTRefSpace,I)
     return SVector(K),SVector{3,Int64}(1,1,1)
 end
 
-function reorder_dof(space::LagrangeRefSpace{Float64,0,3,1},I)
+function reorder_dof(space::LagrangeRefSpace{T,0,3,1},I) where T
 
     return SVector{1,Int64}(1),SVector{1,Int64}(1)
 end
@@ -102,16 +102,19 @@ function reorder_dof(space::LagrangeRefSpace{T,1,4,4},I) where T
 end
 
 function momintegrals!(out, op::VIEOperator,
-    test_local_space::RefSpace, test_ptr, test_tetrahedron_element,
-    trial_local_space::RefSpace, trial_ptr, trial_tetrahedron_element,
+    test_functions::Space, test_ptr, test_tetrahedron_element,
+    trial_functions::Space, trial_ptr, trial_tetrahedron_element,
     strat::SauterSchwab3DStrategy)
 
+    local_test_space = refspace(test_functions)
+    local_trial_space = refspace(trial_functions)
+
     #Find permutation of vertices to match location of singularity to SauterSchwab
     J, I= SauterSchwab3D.reorder(strat.sing)
 
     #Get permutation and rel. orientatio of DoFs 
-    K,O1 = reorder_dof(test_local_space, I)
-    L,O2 = reorder_dof(trial_local_space, J)
+    K,O1 = reorder_dof(local_test_space, I)
+    L,O2 = reorder_dof(local_trial_space, J)
     #Apply permuation to elements
 
     if length(I) == 4
@@ -146,7 +149,7 @@ function momintegrals!(out, op::VIEOperator,
 
     #Define integral (returns a function that only needs barycentric coordinates)
     igd = VIEIntegrand(test_tetrahedron_element, trial_tetrahedron_element,
-        op, test_local_space, trial_local_space)
+        op, local_test_space, local_trial_space)
 
     #Evaluate integral
     Q = SauterSchwab3D.sauterschwab_parameterized(igd, strat)
@@ -161,8 +164,8 @@ function momintegrals!(out, op::VIEOperator,
 end
 
 function momintegrals!(z, biop::VIEOperator,
-    tshs, tptr, tcell,
-    bshs, bptr, bcell,
+    test_functions::Space, tptr, tcell,
+    trial_functions::Space, bptr, bcell,
     strat::DoubleQuadRule)
 
     # memory allocation here is a result from the type instability on strat

test/runtests.jl

@@ -79,6 +79,8 @@ include("test_variational.jl")
 include("test_handlers.jl")
 include("test_gridfunction.jl")
 
+include("test_hh_lsvie.jl")
+
 using TestItemRunner
 @run_package_tests
 

test/test_hh_lsvie.jl

@@ -11,12 +11,12 @@ using Test
 @testset "Lippmann Schwinger Volume Integral Equation" begin
 
     # Environment
-    ε1 = 1.0*ε0
-    μ1 = μ0
+    ε1 = 1.0*SphericalScattering.ε0
+    μ1 = SphericalScattering.μ0
 
     # Dielectic Sphere
-    ε2 = 5.0*ε0
-    μ2 = μ0
+    ε2 = 5.0*SphericalScattering.ε0
+    μ2 = SphericalScattering.μ0
     r = 1.0
     sp = DielectricSphere(; radius = r, filling = Medium(ε2, μ2))
 
@@ -56,7 +56,7 @@ using Test
 
 
     # Assembly
-    b = assemble(Φ_inc, X)
+    b = real.(assemble(Φ_inc, X))
 
     Z_I = assemble(I, X, X)
 
@@ -69,8 +69,8 @@ using Test
     Z_version2 = Z_I + Z_Y
 
     # MoM solution
-    u_version1 = Z_version1 \ Vector(b)
-    u_version2 = Z_version2 \ Vector(b)
+    u_version1 = Z_version1 \ b
+    u_version2 = Z_version2 \ b
 
     # Observation points
     range_ = range(-1.0*r,stop=1.0*r,length=14)
@@ -84,8 +84,8 @@ using Test
     Φ = field(sp, ex, ScalarPotential(points_sp))
 
     # MoM solution inside the dielectric sphere
-    Φ_MoM_version1 = BEAST.grideval(points_sp, u_version1, X, type=Float64)
-    Φ_MoM_version2 = BEAST.grideval(points_sp, u_version2, X, type=Float64)
+    Φ_MoM_version1 = BEAST.grideval(points_sp, u_version1, X)
+    Φ_MoM_version2 = BEAST.grideval(points_sp, u_version2, X)