diff --git a/examples/disabled/tdhh3d_neumann.jl b/examples/disabled/tdhh3d_neumann.jl
index 0c7f6508..adf753ab 100644
--- a/examples/disabled/tdhh3d_neumann.jl
+++ b/examples/disabled/tdhh3d_neumann.jl
@@ -2,8 +2,8 @@ using CompScienceMeshes
 using BEAST
 using LinearAlgebra
 
-# G = meshsphere(1.0, 0.25)
-G = meshsphere(1.0, 0.30)
+#G = meshsphere(1.0, 0.30)
+G = CompScienceMeshes.meshmobius(h=0.2)
 
 c = 1.0
 S = BEAST.HH3DSingleLayerTDBIO(c)
@@ -20,9 +20,9 @@ h = BEAST.gradient(e)
 
 # X = lagrangecxd0(G)
 X = lagrangecxd0(G)
-Y = duallagrangec0d1(G)
-# Y = lagrangecxd0(G)
-
+#Y = duallagrangec0d1(G)
+Y = lagrangecxd0(G)
+numfunctions(X)
 # X = duallagrangecxd0(G, boundary(G))
 # Y = lagrangec0d1(G)
 
@@ -31,21 +31,29 @@ Y = duallagrangec0d1(G)
 Δt, Nt = 0.16, 300
 # T = timebasisc0d1(Δt, Nt)
 P = timebasiscxd0(Δt, Nt)
-H = timebasisc0d1(Δt, Nt)
+#H = timebasisc0d1(Δt, Nt)
 δ = timebasisdelta(Δt, Nt)
 
 # assemble the right hand side
 
 bd  = assemble(n⋅h,     X ⊗ P)
-Z1d = assemble(Id ⊗ Id, X ⊗ P, X ⊗ P, Val{:bandedstorage})
-Z0d = assemble(D,       X ⊗ P, X ⊗ P, Val{:bandedstorage})
+Z1d = assemble(Id ⊗ Id, X ⊗ P, X ⊗ P)
+Z0d = assemble(D,       X ⊗ P, X ⊗ P)
 Zd = Z0d + (-0.5)*Z1d
 u = marchonintime(inv(Zd[:,:,1]), Zd, bd, Nt)
 
 bs = assemble(e, X ⊗ δ)
-Zs = assemble(S, X ⊗ δ, X ⊗ P, Val{:bandedstorage})
+Zs = assemble(S, X ⊗ δ, X ⊗ P)
 v = marchonintime(inv(Zs[:,:,1]), Zs, -bs, Nt)
 
+
+tdacusticsl = @discretise D[j′,j] == 1.0(n⋅h)[j′]   j∈ (X ⊗ P)  j′∈ (X ⊗ δ)
+xacusticsl = solve(tdacusticsl)
+
+
+import Plots
+Plots.plot(xacusticsl[1000,2900:3000])
+
 U, Δω, ω0 = fouriertransform(u, Δt, 0.0, 2)
 V, Δω, ω0 = fouriertransform(v, Δt, 0.0, 2)
 
diff --git a/src/bases/local/bdmlocal.jl b/src/bases/local/bdmlocal.jl
index 1c83f5f7..605963d9 100644
--- a/src/bases/local/bdmlocal.jl
+++ b/src/bases/local/bdmlocal.jl
@@ -21,6 +21,8 @@ end
 
 divergence(ref::BDMRefSpace, sh, el) = [Shape(sh.cellid, 1, sh.coeff/(2*volume(el)))]
 
+numfunctions(x::BDMRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 6
+
 const _vert_perms_bdm = [
     (1,2,3),
     (2,3,1),
diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl
index 21355e0f..0adea811 100644
--- a/src/bases/local/laglocal.jl
+++ b/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)),))
diff --git a/src/bases/local/ncrossbdmlocal.jl b/src/bases/local/ncrossbdmlocal.jl
index 5f7fc211..cb332499 100644
--- a/src/bases/local/ncrossbdmlocal.jl
+++ b/src/bases/local/ncrossbdmlocal.jl
@@ -18,3 +18,5 @@ function (f::NCrossBDMRefSpace{T})(p) where T
         (value= n × (u*tu)/j,         curl=d),
         (value= n × (v*tv)/j,         curl=d),]
 end
+
+numfunctions(x::NCrossBDMRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 6
diff --git a/src/helmholtz3d/timedomain/tdhh3dexc.jl b/src/helmholtz3d/timedomain/tdhh3dexc.jl
index 91287032..5f2d328f 100644
--- a/src/helmholtz3d/timedomain/tdhh3dexc.jl
+++ b/src/helmholtz3d/timedomain/tdhh3dexc.jl
@@ -10,6 +10,7 @@ function planewave(direction, speedoflight::Number, signature, amplitude=one(spe
     PlaneWaveHH3DTD(direction, speedoflight, signature, amplitude)
 end
 
+scalartype(p::PlaneWaveHH3DTD) = scalartype(p.amplitude)
 
 *(a, f::PlaneWaveHH3DTD) = PlaneWaveHH3DTD(f.direction, f.speed_of_light, f.signature, a*f.amplitude)
 
diff --git a/src/postproc/segcurrents.jl b/src/postproc/segcurrents.jl
index a02e37c9..844f2eec 100644
--- a/src/postproc/segcurrents.jl
+++ b/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
diff --git a/src/volumeintegral/sauterschwab_ints.jl b/src/volumeintegral/sauterschwab_ints.jl
index 68385e7a..6b8fc46d 100644
--- a/src/volumeintegral/sauterschwab_ints.jl
+++ b/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
diff --git a/test/runtests.jl b/test/runtests.jl
index bbcee67f..15214f95 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -77,8 +77,12 @@ include("test_td_tensoroperator.jl")
 include("test_variational.jl")
 
 include("test_handlers.jl")
+include("test_ncrossbdm.jl")
+include("test_curl_lagc0d1_lagc0d2.jl")
 include("test_gridfunction.jl")
 
+include("test_hh_lsvie.jl")
+
 using TestItemRunner
 @run_package_tests
 
diff --git a/test/test_hh_lsvie.jl b/test/test_hh_lsvie.jl
index 899c55b1..cc16cb35 100644
--- a/test/test_hh_lsvie.jl
+++ b/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)