diff --git a/src/BEAST.jl b/src/BEAST.jl
index bb9f8bbe..79e2fae0 100644
--- a/src/BEAST.jl
+++ b/src/BEAST.jl
@@ -214,7 +214,7 @@ include("maxwell/spotential.jl")
 include("maxwell/maxwell.jl")
 include("maxwell/sourcefield.jl")
 
-# Support for the Helmholtz equatio
+# Support for the Helmholtz equation
 include("helmholtz2d/helmholtzop.jl")
 
 include("helmholtz3d/hh3dexc.jl")
@@ -223,6 +223,7 @@ include("helmholtz3d/nitsche.jl")
 include("helmholtz3d/hh3dnear.jl")
 include("helmholtz3d/hh3dfar.jl")
 include("helmholtz3d/hh3d_sauterschwabqr.jl")
+include("helmholtz3d/wiltonints.jl")
 
 #suport for Volume Integral equation
 include("volumeintegral/vie.jl")
diff --git a/src/helmholtz3d/hh3d_sauterschwabqr.jl b/src/helmholtz3d/hh3d_sauterschwabqr.jl
index ea2298ec..49412cde 100644
--- a/src/helmholtz3d/hh3d_sauterschwabqr.jl
+++ b/src/helmholtz3d/hh3d_sauterschwabqr.jl
@@ -1,6 +1,6 @@
 function pulled_back_integrand(op::HH3DSingleLayerFDBIO,
-    test_local_space::LagrangeRefSpace{<:Any,0},
-    trial_local_space::LagrangeRefSpace{<:Any,0},
+    test_local_space::LagrangeRefSpace,
+    trial_local_space::LagrangeRefSpace,
     test_chart, trial_chart)
 
     (u,v) -> begin
@@ -20,14 +20,14 @@ function pulled_back_integrand(op::HH3DSingleLayerFDBIO,
 
         αjG = α*G*j
 
-        SA[f[1].value * αjG * g[1].value]
+        SMatrix{length(f),length(g)}((f[j].value * αjG * g[i].value for i in 1:length(g) for j in 1:length(f) )...)
     end
 end
 
 
 function pulled_back_integrand(op::HH3DHyperSingularFDBIO,
-    test_local_space::LagrangeRefSpace{<:Any,1},
-    trial_local_space::LagrangeRefSpace{<:Any,1},
+    test_local_space::LagrangeRefSpace,
+    trial_local_space::LagrangeRefSpace,
     test_chart, trial_chart)
 
     (u,v) -> begin
@@ -52,25 +52,16 @@ function pulled_back_integrand(op::HH3DHyperSingularFDBIO,
         αjG = ny*α*G*j
         βjG = β*G*j
 
-        A = SA[ αjG*g[1].value, αjG*g[2].value, αjG*g[3].value]
-        B = SA[ βjG*g[1].curl, βjG*g[2].curl, βjG*g[3].curl ] 
-
-        SMatrix{3,3}((
-            dot(nx*f[1].value, A[1]) + dot(f[1].curl, B[1]),
-            dot(nx*f[2].value, A[1]) + dot(f[2].curl, B[1]),
-            dot(nx*f[3].value, A[1]) + dot(f[3].curl, B[1]),
-            dot(nx*f[1].value, A[2]) + dot(f[1].curl, B[2]),
-            dot(nx*f[2].value, A[2]) + dot(f[2].curl, B[2]),
-            dot(nx*f[3].value, A[2]) + dot(f[3].curl, B[2]),
-            dot(nx*f[1].value, A[3]) + dot(f[1].curl, B[3]),
-            dot(nx*f[2].value, A[3]) + dot(f[2].curl, B[3]),
-            dot(nx*f[3].value, A[3]) + dot(f[3].curl, B[3])))
+        A = SA[(αjG*g[i].value for i in 1:length(g))...]
+        B = SA[(βjG*g[i].curl for i in 1:length(g))...]
+
+        SMatrix{length(f),length(g)}((((dot(nx*f[j].value,A[i])+dot(f[j].curl,B[i])) for i in 1:length(g) for j in 1:length(f))...))
     end
 end
 
 function pulled_back_integrand(op::HH3DDoubleLayerFDBIO,
-    test_local_space::LagrangeRefSpace{<:Any,0},
-    trial_local_space::LagrangeRefSpace{<:Any,1},
+    test_local_space::LagrangeRefSpace,
+    trial_local_space::LagrangeRefSpace,
     test_chart, trial_chart)
 
     (u,v) -> begin
@@ -94,16 +85,13 @@ function pulled_back_integrand(op::HH3DDoubleLayerFDBIO,
         ∇G = -(γ + inv_R) * G * inv_R * r
         αnyj∇G = dot(ny,-α*∇G*j)
 
-        SMatrix{1,3}(
-            f[1].value * αnyj∇G * g[1].value,
-            f[1].value * αnyj∇G * g[2].value,
-            f[1].value * αnyj∇G * g[3].value)
+        SMatrix{length(f),length(g)}((f[j].value * αnyj∇G * g[i].value  for i in 1:length(g) for j in 1:length(f))...)
     end
 end
 
 function pulled_back_integrand(op::HH3DDoubleLayerTransposedFDBIO,
-    test_local_space::LagrangeRefSpace{<:Any,1},
-    trial_local_space::LagrangeRefSpace{<:Any,0},
+    test_local_space::LagrangeRefSpace,
+    trial_local_space::LagrangeRefSpace,
     test_chart, trial_chart)
 
     (u,v) -> begin
@@ -127,14 +115,11 @@ function pulled_back_integrand(op::HH3DDoubleLayerTransposedFDBIO,
         ∇G = -(γ + inv_R) * G * inv_R * r
         αnxj∇G = dot(nx,α*∇G*j)
 
-        SMatrix{3,1}(
-            f[1].value * αnxj∇G * g[1].value,
-            f[2].value * αnxj∇G * g[1].value,
-            f[3].value * αnxj∇G * g[1].value)
+        SMatrix{length(f),length(g)}((f[j].value * αnxj∇G * g[i].value for i in 1:length(g) for j in 1:length(f))...)
     end
 end
 
-function momintegrals!(op::HH3DSingleLayerFDBIO,
+function momintegrals!(op::Helmholtz3DOp,
     test_local_space::LagrangeRefSpace{<:Any,0}, trial_local_space::LagrangeRefSpace{<:Any,0},
     test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)
 
@@ -152,19 +137,21 @@ function momintegrals!(op::HH3DSingleLayerFDBIO,
         trial_triangular_element.vertices[J[2]],
         trial_triangular_element.vertices[J[3]])
 
-    # igd = MWSL3DIntegrand(test_triangular_element, trial_triangular_element,
-    #     op, test_local_space, trial_local_space)
+    test_sign = Combinatorics.levicivita(I)
+    trial_sign = Combinatorics.levicivita(J)
+    σ = momintegrals_sign(op, test_sign, trial_sign)
+
     igd = pulled_back_integrand(op, test_local_space, trial_local_space,
         test_triangular_element, trial_triangular_element)
     G = SauterSchwabQuadrature.sauterschwab_parameterized(igd, strat)
-    out[1,1] += G[1,1]
+    out[1,1] += G[1,1] * σ
 
     nothing
 end
 
 
-function momintegrals!(op::HH3DHyperSingularFDBIO,
-    test_local_space::LagrangeRefSpace{<:Any,1}, trial_local_space::LagrangeRefSpace{<:Any,1},
+function momintegrals!(op::Helmholtz3DOp,
+    test_local_space::LagrangeRefSpace, trial_local_space::LagrangeRefSpace,
     test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)
 
     I, J, K, L = SauterSchwabQuadrature.reorder(
@@ -183,13 +170,14 @@ function momintegrals!(op::HH3DHyperSingularFDBIO,
 
     test_sign = Combinatorics.levicivita(I)
     trial_sign = Combinatorics.levicivita(J)
-    σ = test_sign * trial_sign
+
+    σ = momintegrals_sign(op, test_sign, trial_sign)
 
     igd = pulled_back_integrand(op, test_local_space, trial_local_space,
         test_triangular_element, trial_triangular_element)
     G = SauterSchwabQuadrature.sauterschwab_parameterized(igd, strat)
     for j ∈ 1:3, i ∈ 1:3
-        out[i,j] += G[K[i],L[j]] * σ 
+        out[i,j] += G[K[i],L[j]] * σ
     end
 
     nothing
@@ -213,9 +201,10 @@ function momintegrals!(op::Helmholtz3DOp,
         trial_triangular_element.vertices[J[2]],
         trial_triangular_element.vertices[J[3]])
 
+    test_sign = Combinatorics.levicivita(I)
     trial_sign = Combinatorics.levicivita(J)
 
-    σ = trial_sign
+    σ = momintegrals_sign(op, test_sign, trial_sign)
 
     igd = pulled_back_integrand(op, test_local_space, trial_local_space,
         test_triangular_element, trial_triangular_element)
@@ -247,7 +236,9 @@ function momintegrals!(op::Helmholtz3DOp,
         trial_triangular_element.vertices[J[3]])
 
         test_sign = Combinatorics.levicivita(I)
-        σ = test_sign
+        trial_sign = Combinatorics.levicivita(J)
+
+        σ = momintegrals_sign(op, test_sign, trial_sign)
 
     igd = pulled_back_integrand(op, test_local_space, trial_local_space,
         test_triangular_element, trial_triangular_element)
@@ -259,3 +250,16 @@ function momintegrals!(op::Helmholtz3DOp,
 
     nothing
 end
+
+function momintegrals_sign(op::HH3DSingleLayerFDBIO, test_sign, trial_sign)
+    return 1
+end
+function momintegrals_sign(op::HH3DDoubleLayerFDBIO, test_sign, trial_sign)
+    return trial_sign
+end
+function momintegrals_sign(op::HH3DDoubleLayerTransposedFDBIO, test_sign, trial_sign)
+    return test_sign
+end
+function momintegrals_sign(op::HH3DHyperSingularFDBIO, test_sign, trial_sign)
+    return test_sign * trial_sign
+end
diff --git a/src/helmholtz3d/hh3dops.jl b/src/helmholtz3d/hh3dops.jl
index 19e99e8f..47e2b39c 100644
--- a/src/helmholtz3d/hh3dops.jl
+++ b/src/helmholtz3d/hh3dops.jl
@@ -16,6 +16,23 @@ struct HH3DHyperSingularFDBIO{T,K} <: Helmholtz3DOp
     gamma::K
 end
 
+struct HH3DHyperSingularReg{T,K} <: Helmholtz3DOpReg
+    "coefficient of the weakly singular term"
+    alpha::T
+    "coefficient of the hyper singular term"
+    beta::T
+    "`im*κ` with `κ` the wave number"
+    gamma::K
+end
+
+struct HH3DHyperSingularSng{T,K} <: Helmholtz3DOp
+    "coefficient of the weakly singular term"
+    alpha::T
+    "coefficient of the hyper singular term"
+    beta::T
+    "`im*κ` with `κ` the wave number"
+    gamma::K
+end
 HH3DHyperSingularFDBIO(gamma) = HH3DHyperSingularFDBIO(gamma^2, one(gamma), gamma)
 scalartype(op::HH3DHyperSingularFDBIO) = promote_type(typeof(op.alpha), typeof(op.beta), typeof(op.gamma))
 
@@ -45,11 +62,31 @@ struct HH3DDoubleLayerFDBIO{T,K} <: Helmholtz3DOp
     gamma::K
 end
 
+struct HH3DDoubleLayerReg{T,K} <: Helmholtz3DOpReg
+    alpha::T
+    gamma::K
+end
+
+struct HH3DDoubleLayerSng{T,K} <: Helmholtz3DOp
+    alpha::T
+    gamma::K
+end
+
 struct HH3DDoubleLayerTransposedFDBIO{T,K} <: Helmholtz3DOp
     alpha::T
     gamma::K
 end
 
+struct HH3DDoubleLayerTransposedReg{T,K} <: Helmholtz3DOpReg
+    alpha::T
+    gamma::K
+end
+
+struct HH3DDoubleLayerTransposedSng{T,K} <: Helmholtz3DOp
+    alpha::T
+    gamma::K
+end
+
 defaultquadstrat(::Helmholtz3DOp, ::LagrangeRefSpace, ::LagrangeRefSpace) = DoubleNumWiltonSauterQStrat(2,3,2,3,4,4,4,4)
 
 function quaddata(op::Helmholtz3DOp, test_refspace::LagrangeRefSpace,
@@ -116,9 +153,9 @@ end
 
 
 
-function quadrule(op::HH3DSingleLayerFDBIO,
-        test_refspace::LagrangeRefSpace{T,0} where T,
-        trial_refspace::LagrangeRefSpace{T,0} where T,
+function quadrule(op::Helmholtz3DOp,
+        test_refspace::LagrangeRefSpace,
+        trial_refspace::LagrangeRefSpace,
         i, test_element, j, trial_element, qd,
         qs::DoubleNumWiltonSauterQStrat)
 
@@ -165,10 +202,12 @@ function quadrule(op::HH3DSingleLayerFDBIO,
         qd.bsis_qp[1,j])
 end
 
+regularpart(op::HH3DHyperSingularFDBIO) = HH3DHyperSingularReg(op.alpha, op.beta, op.gamma)
+singularpart(op::HH3DHyperSingularFDBIO) = HH3DHyperSingularSng(op.alpha, op.beta, op.gamma)
 
-function quadrule(op::HH3DHyperSingularFDBIO,
-    test_refspace::LagrangeRefSpace{T,1} where T,
-    trial_refspace::LagrangeRefSpace{T,1} where T,
+#= function quadrule(op::HH3DHyperSingularFDBIO,
+    test_refspace::LagrangeRefSpace,
+    trial_refspace::LagrangeRefSpace,
     i, test_element, j, trial_element, qd,
     qs::DoubleNumWiltonSauterQStrat)
 
@@ -182,19 +221,21 @@ function quadrule(op::HH3DHyperSingularFDBIO,
     hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
     hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])
 
-    # test_quadpoints  = qd.test_qp
-    # trial_quadpoints = qd.bsis_qp
-
-    # hits != 0 && return WiltonSERule(
-    #     test_quadpoints[1,i],
-    #     DoubleQuadRule(
-    #         test_quadpoints[1,i],
-    #         trial_quadpoints[1,j]))
+    test_quadpoints  = qd.test_qp
+    trial_quadpoints = qd.bsis_qp
+    h2 = volume(trial_element)
+    xtol2 = 0.2 * 0.2
+    k2 = abs2(op.gamma)
+    max(dmin2*k2, dmin2/16h2) < xtol2 && return WiltonSERule(
+        test_quadpoints[1,i],
+        DoubleQuadRule(
+            test_quadpoints[1,i],
+            trial_quadpoints[1,j]))
 
     return DoubleQuadRule(
         qd.test_qp[1,i],
         qd.bsis_qp[1,j])
-end
+end =#
 
 
 function quadrule(op::HH3DHyperSingularFDBIO,
@@ -271,7 +312,7 @@ function quadrule(op::Helmholtz3DOp,
         quadrature_data[2][1,j])
 end
 
-function quadrule(op::HH3DDoubleLayerTransposedFDBIO,
+#= function quadrule(op::HH3DDoubleLayerTransposedFDBIO,
     test_refspace::LagrangeRefSpace{T,1} where T,
     trial_refspace::LagrangeRefSpace{T,0} where T,
     i, test_element, j, trial_element, quadrature_data,
@@ -295,7 +336,7 @@ function quadrule(op::HH3DDoubleLayerTransposedFDBIO,
     trial_quadpoints = quadrature_data[2]
     test_quadpoints  = quadrature_data.test_qp
     trial_quadpoints = quadrature_data.bsis_qp
-#=     h2 = volume(trial_element)
+    h2 = volume(trial_element)
     xtol2 = 0.2 * 0.2
     k2 = abs2(op.gamma)
 
@@ -303,13 +344,14 @@ function quadrule(op::HH3DDoubleLayerTransposedFDBIO,
         test_quadpoints[1,i],
         DoubleQuadRule(
             test_quadpoints[1,i],
-            trial_quadpoints[1,j])) =#
+            trial_quadpoints[1,j]))
+
     return DoubleQuadRule(
         quadrature_data[1][1,i],
         quadrature_data[2][1,j])
-end
+end =#
 
-function quadrule(op::HH3DDoubleLayerFDBIO,
+#= function quadrule(op::HH3DDoubleLayerFDBIO,
     test_refspace::LagrangeRefSpace{T,0} where T,
     trial_refspace::LagrangeRefSpace{T,1} where T,
     i, test_element, j, trial_element, quadrature_data,
@@ -345,7 +387,7 @@ function quadrule(op::HH3DDoubleLayerFDBIO,
     return DoubleQuadRule(
         quadrature_data[1][1,i],
         quadrature_data[2][1,j])
-end
+end =#
 
 
 function quadrule(op::Helmholtz3DOp,
@@ -433,60 +475,23 @@ function (igd::Integrand{<:HH3DSingleLayerFDBIO})(x,y,f,g)
     end
 end
 
-function (igd::Integrand{<:HH3DSingleLayerReg})(x,y,f,g)
-    α = igd.operator.alpha
-    γ = igd.operator.gamma
-
-    r = cartesian(x) - cartesian(y)
-    R = norm(r)
-    iR = 1 / R
-    green = (expm1(-γ*R) + γ*R - 0.5*γ^2*R^2) / (4pi*R)
-    αG = α * green
-
-    _integrands(f,g) do fi, gi
-        dot(gi.value, αG*fi.value)
-    end
-end
-
-
 function integrand(op::Union{HH3DSingleLayerFDBIO,HH3DSingleLayerReg},
-        kernel, test_values, test_element, trial_values, trial_element)
-
-    α = op.alpha
-    G = kernel.green
-
-    g = test_values.value
-    f = trial_values.value
-
-    α*dot(g, G*f)
-end
-
+    kernel, test_values, test_element, trial_values, trial_element)
 
-function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood,
-        test_refspace::LagrangeRefSpace{T,0} where {T},
-        trial_refspace::LagrangeRefSpace{T,0} where {T},
-        test_elements, trial_element, zlocal, quadrature_rule::WiltonSERule, dx)
+α = op.alpha
+G = kernel.green
 
-    γ = op.gamma
-    α = op.alpha
+g = test_values.value
+f = trial_values.value
 
-    s1, s2, s3 = trial_element.vertices
-
-    x = cartesian(test_neighborhood)
-    n = normalize((s1-s3)×(s2-s3))
-    ρ = x - dot(x - s1, n) * n
-
-    scal, vec = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1})
-    ∫G = (scal[2] - γ*scal[3] + 0.5*γ^2*scal[4]) / (4π)
-
-    zlocal[1,1] += α * ∫G * dx
-    return nothing
+α*dot(g, G*f)
 end
 
 
-
 HH3DDoubleLayerFDBIO(gamma) = HH3DDoubleLayerFDBIO(one(gamma), gamma)
 
+regularpart(op::HH3DDoubleLayerFDBIO) = HH3DDoubleLayerReg(op.alpha, op.gamma)
+singularpart(op::HH3DDoubleLayerFDBIO) = HH3DDoubleLayerSng(op.alpha, op.gamma)
 
 function (igd::Integrand{<:HH3DDoubleLayerFDBIO})(x,y,f,g)
     γ = igd.operator.gamma
@@ -515,6 +520,9 @@ end
 
  HH3DDoubleLayerTransposedFDBIO(gamma) = HH3DDoubleLayerTransposedFDBIO(one(gamma), gamma)
 
+regularpart(op::HH3DDoubleLayerTransposedFDBIO) = HH3DDoubleLayerTransposedReg(op.alpha, op.gamma)
+singularpart(op::HH3DDoubleLayerTransposedFDBIO) = HH3DDoubleLayerTransposedSng(op.alpha, op.gamma)
+
 function (igd::Integrand{<:HH3DDoubleLayerTransposedFDBIO})(x,y,f,g)
     γ = igd.operator.gamma
 
diff --git a/src/helmholtz3d/wiltonints.jl b/src/helmholtz3d/wiltonints.jl
new file mode 100644
index 00000000..c1cf7386
--- /dev/null
+++ b/src/helmholtz3d/wiltonints.jl
@@ -0,0 +1,432 @@
+# single layer
+function (igd::Integrand{<:HH3DSingleLayerReg})(x,y,f,g)
+    α = igd.operator.alpha
+    γ = igd.operator.gamma
+
+    r = cartesian(x) - cartesian(y)
+    R = norm(r)
+    iR = 1 / R
+    green = (expm1(-γ*R) - 0.5*γ^2*R^2) / (4pi*R)
+    αG = α * green
+
+    _integrands(f,g) do fi, gi
+        dot(gi.value, αG*fi.value)
+    end
+end
+
+function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood,
+    test_refspace::LagrangeRefSpace{T,0} where {T},
+    trial_refspace::LagrangeRefSpace{T,0} where {T},
+    test_elements, trial_element, zlocal, quadrature_rule::WiltonSERule, dx)
+
+γ = op.gamma
+α = op.alpha
+
+s1, s2, s3 = trial_element.vertices
+
+x = cartesian(test_neighborhood)
+n = normalize((s1-s3)×(s2-s3))
+ρ = x - dot(x - s1, n) * n
+
+scal, vec = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1})
+∫G = (scal[2] + 0.5*γ^2*scal[4]) / (4π)
+
+zlocal[1,1] += α * ∫G * dx
+return nothing
+end
+
+# single layer with patch basis and pyramid testing
+function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood,
+    test_refspace::LagrangeRefSpace{T,1} where {T},
+    trial_refspace::LagrangeRefSpace{T,0} where {T},
+    test_elements, trial_element, zlocal, quadrature_rule::WiltonSERule, dx)
+
+γ = op.gamma
+α = op.alpha
+
+s1, s2, s3 = trial_element.vertices
+t1, t2, t3 = test_elements.vertices
+
+
+x = cartesian(test_neighborhood)
+n = normalize((s1-s3)×(s2-s3))
+ρ = x - dot(x - s1, n) * n
+
+scal, vec = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1})
+∫G = (scal[2] + 0.5*γ^2*scal[4]) / (4π)
+Atot = 1/2*norm(cross(t3-t1,t3-t2))
+for i in 1:numfunctions(test_refspace)
+    Ai = 1/2*norm(cross(test_elements.vertices[mod1(i-1,3)]-x,test_elements.vertices[mod1(i+1,3)]-x))
+    g = Ai/Atot
+    for j in 1:numfunctions(trial_refspace)
+        zlocal[i,j] += α * ∫G * g * dx
+    end
+end
+return nothing
+end
+
+# single layer with pyramid basis and patch testing
+function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood,
+    test_refspace::LagrangeRefSpace{T,0} where {T},
+    trial_refspace::LagrangeRefSpace{T,1} where {T},
+    test_elements, trial_element, zlocal, quadrature_rule::WiltonSERule, dx)
+
+γ = op.gamma
+α = op.alpha
+
+s1, s2, s3 = trial_element.vertices
+
+x = cartesian(test_neighborhood)
+n = normalize((s1-s3)×(s2-s3))
+ρ = x - dot(x - s1, n) * n
+
+∫Rⁿ, ∫RⁿN = WiltonInts84.higherorder(s1,s2,s3,x,3)
+∫G = (∫RⁿN[2] + 0.5*γ^2*∫RⁿN[3]) / (4π)
+
+for i in 1:numfunctions(trial_refspace)
+zlocal[1,i] += α * ∫G[i] * dx
+end
+return nothing
+end
+
+#single layer with pyramid basis and pyramid testing
+function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood,
+    test_refspace::LagrangeRefSpace{T,1} where {T},
+    trial_refspace::LagrangeRefSpace{T,1} where {T},
+    test_elements, trial_element, zlocal, quadrature_rule::WiltonSERule, dx)
+
+γ = op.gamma
+α = op.alpha
+
+s1, s2, s3 = trial_element.vertices
+t1, t2, t3 = test_elements.vertices
+
+x = cartesian(test_neighborhood)
+n = normalize((s1-s3)×(s2-s3))
+ρ = x - dot(x - s1, n) * n
+
+∫Rⁿ, ∫RⁿN = WiltonInts84.higherorder(s1,s2,s3,x,3)
+∫G = (∫RⁿN[2] + 0.5*γ^2*∫RⁿN[3]) / (4π)
+
+Atot = 1/2*norm(cross(t3-t1,t3-t2))
+for i in 1:numfunctions(test_refspace)
+    Ai = 1/2*norm(cross(test_elements.vertices[mod1(i-1,3)]-x,test_elements.vertices[mod1(i+1,3)]-x))
+    g = Ai/Atot
+    for j in 1:numfunctions(trial_refspace)
+        zlocal[i,j] += α * ∫G[j] * g * dx
+    end
+end
+return nothing
+end
+
+# double layer transposed
+function (igd::Integrand{<:HH3DDoubleLayerTransposedReg})(x,y,f,g)
+    γ = igd.operator.gamma
+    
+    r = cartesian(x) - cartesian(y)
+    R = norm(r)
+    iR = 1/R
+    γR = γ*R
+    expo = exp(-γR)
+
+    gradgreen = ( -(γR + 1)*expo + (1 - 0.5*γR^2) ) * (i4pi*iR^3) * r
+
+    n = normal(x)
+    fvalue = getvalue(f)
+    gvalue = getvalue(g)
+
+    return _krondot(fvalue,gvalue) * dot(n, gradgreen)
+end
+
+#double layer transposed with patch basis and pyramid testing
+function innerintegrals!(op::HH3DDoubleLayerTransposedSng, test_neighborhood,
+    test_refspace::LagrangeRefSpace{T,0,3,1} where {T},
+    trial_refspace::LagrangeRefSpace{T,0,3,1} where {T},
+    test_elements, trial_element, zlocal, quadrature_rule::WiltonSERule, dx)
+    γ = op.gamma
+    α = op.alpha
+
+    s1, s2, s3 = trial_element.vertices
+    t1, t2, t3 = test_elements.vertices
+    x = cartesian(test_neighborhood)
+    n = normalize((t1-t3)×(t2-t3))
+    ρ = x - dot(x - s1, n) * n
+
+    scal, vec, grad = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1})
+
+    ∫∇G = -(grad[1]+0.5*γ^2*grad[3])/(4π)
+    ∫n∇G = dot(n,∫∇G)
+    for i in 1:numfunctions(test_refspace)
+        for j in 1:numfunctions(trial_refspace)
+            zlocal[i,j] += α * ∫n∇G * dx
+        end
+    end
+
+    return nothing
+end
+#double layer transposed with patch basis and pyramid testing
+function innerintegrals!(op::HH3DDoubleLayerTransposedSng, test_neighborhood,
+    test_refspace::LagrangeRefSpace{T,1,3,3} where {T},
+    trial_refspace::LagrangeRefSpace{T,0,3,1} where {T},
+    test_elements, trial_element, zlocal, quadrature_rule::WiltonSERule, dx)
+    γ = op.gamma
+    α = op.alpha
+
+    s1, s2, s3 = trial_element.vertices
+    t1, t2, t3 = test_elements.vertices
+    x = cartesian(test_neighborhood)
+    n = normalize((t1-t3)×(t2-t3))
+    ρ = x - dot(x - s1, n) * n
+
+    scal, vec, grad = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1})
+
+    ∫∇G = -(grad[1]+0.5*γ^2*grad[3])/(4π)
+    ∫n∇G = dot(n,∫∇G)
+    Atot = 1/2*norm(cross(t3-t1,t3-t2))
+    for i in 1:numfunctions(test_refspace)
+        Ai = 1/2*norm(cross(test_elements.vertices[mod1(i-1,3)]-x,test_elements.vertices[mod1(i+1,3)]-x))
+        g = Ai/Atot
+        for j in 1:numfunctions(trial_refspace)
+            zlocal[i,j] += α * ∫n∇G * g * dx
+        end
+    end
+
+    return nothing
+end
+#double layer transposed with pyramid basis and patch testing
+function innerintegrals!(op::HH3DDoubleLayerTransposedSng, test_neighborhood,
+    test_refspace::LagrangeRefSpace{T,0,3,1} where {T},
+    trial_refspace::LagrangeRefSpace{T,1,3,3} where {T},
+    test_elements, trial_element, zlocal, quadrature_rule::WiltonSERule, dx)
+    γ = op.gamma
+    α = op.alpha
+
+    s1, s2, s3 = trial_element.vertices
+    t1, t2, t3 = test_elements.vertices
+    x = cartesian(test_neighborhood)
+    n = normalize((t1-t3)×(t2-t3))
+    ρ = x - dot(x - s1, n) * n
+
+    _, _, _, grad = WiltonInts84.higherorder(s1,s2,s3,x,3)
+
+    ∫∇G = -(grad[1] + 0.5*γ^2*grad[2]) / (4π)
+    for i in 1:numfunctions(test_refspace)
+        for j in 1:numfunctions(trial_refspace)
+            ∫n∇G = dot(n,∫∇G[j])
+            zlocal[i,j] += α * ∫n∇G  * dx
+        end
+    end
+
+    return nothing
+end
+#double layer transposed with pyramid basis and pyramid testing
+function innerintegrals!(op::HH3DDoubleLayerTransposedSng, test_neighborhood,
+    test_refspace::LagrangeRefSpace{T,1,3,3} where {T},
+    trial_refspace::LagrangeRefSpace{T,1,3,3} where {T},
+    test_elements, trial_element, zlocal, quadrature_rule::WiltonSERule, dx)
+    γ = op.gamma
+    α = op.alpha
+
+    s1, s2, s3 = trial_element.vertices
+    t1, t2, t3 = test_elements.vertices
+    x = cartesian(test_neighborhood)
+    n = normalize((t1-t3)×(t2-t3))
+    ρ = x - dot(x - s1, n) * n
+
+    _, _, _, grad = WiltonInts84.higherorder(s1,s2,s3,x,3)
+
+    ∫∇G = -(grad[1] + 0.5*γ^2*grad[2]) / (4π)
+
+    Atot = 1/2*norm(cross(t3-t1,t3-t2))
+    for i in 1:numfunctions(test_refspace)
+        Ai = 1/2*norm(cross(test_elements.vertices[mod1(i-1,3)]-x,test_elements.vertices[mod1(i+1,3)]-x))
+        g = Ai/Atot
+        for j in 1:numfunctions(trial_refspace)
+            ∫n∇G = dot(n,∫∇G[j])
+            zlocal[i,j] += α * ∫n∇G * g * dx
+        end
+    end
+
+    return nothing
+end
+# double layer
+function (igd::Integrand{<:HH3DDoubleLayerReg})(x,y,f,g)
+    γ = igd.operator.gamma
+
+    r = cartesian(x) - cartesian(y)
+    R = norm(r)
+    iR = 1/R
+    γR = γ*R
+    expo = exp(-γR)
+
+    gradgreen = ( -(γR + 1)*expo + (1 - 0.5*γR^2) ) * (i4pi*iR^3) * r
+
+    n = normal(y)
+    fvalue = getvalue(f)
+    gvalue = getvalue(g)
+
+    return _krondot(fvalue,gvalue) * dot(n, -gradgreen)
+
+end
+#double layer with pyramid basis and pyramid testing
+function innerintegrals!(op::HH3DDoubleLayerSng, p,
+    g::LagrangeRefSpace{T,1} where {T},
+    f::LagrangeRefSpace{T,1} where {T},
+    t, s, z, quadrature_rule::WiltonSERule, dx)
+    γ = op.gamma
+    α = op.alpha
+
+    s1, s2, s3 = s.vertices
+    t1, t2, t3 = t.vertices
+
+    x = cartesian(p)
+    n = normalize((s1-s3)×(s2-s3))
+    ρ = x - dot(x - s1, n) * n
+
+    _, _, _, grad = WiltonInts84.higherorder(s1,s2,s3,x,3)
+
+    ∫∇G = -(grad[1] + 0.5*γ^2*grad[2]) / (4π)
+
+    Atot = 1/2*norm(cross(t3-t1,t3-t2))
+    for i in 1:numfunctions(g)
+        Ai = 1/2*norm(cross(t.vertices[mod1(i-1,3)]-x,t.vertices[mod1(i+1,3)]-x))
+        g = Ai/Atot
+        for j in 1:numfunctions(f)     
+        z[i,j] += α * dot(n,-∫∇G[j]) * g * dx
+        end
+    end
+
+    return nothing
+end
+#double layer with patch basis and pyramid testing
+function innerintegrals!(op::HH3DDoubleLayerSng, p,
+    g::LagrangeRefSpace{T,0} where {T},
+    f::LagrangeRefSpace{T,0} where {T},
+    t, s, z, quadrature_rule::WiltonSERule, dx)
+    γ = op.gamma
+    α = op.alpha
+
+    s1, s2, s3 = s.vertices
+    t1, t2, t3 = t.vertices
+
+    x = cartesian(p)
+    n = normalize((s1-s3)×(s2-s3))
+    ρ = x - dot(x - s1, n) * n
+
+    scal, vec, grad = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1})
+
+    ∫∇G = -(grad[1]+0.5*γ^2*grad[3])/(4π)
+
+    for i in 1:numfunctions(g)
+        for j in 1:numfunctions(f)
+            z[i,j] += α * dot(n,-∫∇G) * dx #why the minus?
+        end
+    end
+    return nothing
+end
+#double layer with patch basis and pyramid testing
+function innerintegrals!(op::HH3DDoubleLayerSng, p,
+    g::LagrangeRefSpace{T,1} where {T},
+    f::LagrangeRefSpace{T,0} where {T},
+    t, s, z, quadrature_rule::WiltonSERule, dx)
+    γ = op.gamma
+    α = op.alpha
+
+    s1, s2, s3 = s.vertices
+    t1, t2, t3 = t.vertices
+
+    x = cartesian(p)
+    n = normalize((s1-s3)×(s2-s3))
+    ρ = x - dot(x - s1, n) * n
+
+    scal, vec, grad = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1})
+
+    ∫∇G = -(grad[1]+0.5*γ^2*grad[3])/(4π)
+
+    Atot = 1/2*norm(cross(t3-t1,t3-t2))
+for i in 1:numfunctions(g)
+    Ai = 1/2*norm(cross(t.vertices[mod1(i-1,3)]-x,t.vertices[mod1(i+1,3)]-x))
+    g = Ai/Atot
+    for j in 1:numfunctions(f)
+       z[i,j] += α * dot(n,-∫∇G) * g * dx
+    end
+end
+    return nothing
+end
+
+#double layer with pyramid basis and patch testing
+function innerintegrals!(op::HH3DDoubleLayerSng, p,
+    g::LagrangeRefSpace{T,0} where {T},
+    f::LagrangeRefSpace{T,1} where {T},
+    t, s, z, quadrature_rule::WiltonSERule, dx)
+    γ = op.gamma
+    α = op.alpha
+
+    s1, s2, s3 = s.vertices
+
+    x = cartesian(p)
+    n = normalize((s1-s3)×(s2-s3))
+    ρ = x - dot(x - s1, n) * n
+
+    _, _, _, grad = WiltonInts84.higherorder(s1,s2,s3,x,3)
+
+    ∫∇G = -(grad[1] + 0.5*γ^2*grad[2]) / (4π)
+
+  
+for i in 1:numfunctions(g)
+    for j in 1:numfunctions(f)
+       z[i,j] += α * dot(n,-∫∇G[j]) * dx
+    end
+end
+
+    return nothing
+end
+
+
+function (igd::Integrand{<:HH3DHyperSingularReg})(x,y,f,g)
+    α = igd.operator.alpha
+    β = igd.operator.beta
+    γ = igd.operator.gamma
+
+    r = cartesian(x) - cartesian(y)
+    R = norm(r)
+    iR = 1 / R
+    green = (expm1(-γ*R) - 0.5*γ^2*R^2) / (4pi*R)
+    nx = normal(x)
+    ny = normal(y)
+
+    _integrands(f,g) do fi, gi
+        α*dot(nx,ny)*gi.value*fi.value*green + β*dot(gi.curl,fi.curl)*green
+    end
+end
+
+function innerintegrals!(op::HH3DHyperSingularSng, p,
+    g::LagrangeRefSpace{T,1} where {T},
+    f::LagrangeRefSpace{T,1} where {T},
+    t, s, z, quadrature_rule::WiltonSERule, dx)
+    α = op.alpha
+    β = op.beta
+    γ = op.gamma
+    s1, s2, s3 = s.vertices
+    t1, t2, t3 = t.vertices
+    x = cartesian(p)
+    nx = normalize((s1-s3)×(s2-s3))
+    ny = normalize((t1-t3)×(t2-t3))
+
+    ∫Rⁿ, ∫RⁿN = WiltonInts84.higherorder(s1,s2,s3,x,3)
+    greenconst = (∫Rⁿ[2] + 0.5*γ^2*∫Rⁿ[3]) / (4π)
+    greenlinear = (∫RⁿN[2] + 0.5*γ^2*∫RⁿN[3] ) / (4π)
+
+    jt = volume(t) * factorial(dimension(t))
+    js = volume(s) * factorial(dimension(s))
+    curlt = [(t3-t2)/jt,(t1-t3)/jt,(t2-t1)/jt]
+    curls = [(s3-s2)/js,(s1-s3)/js,(s2-s1)/js]
+    Atot = 1/2*norm(cross(t3-t1,t3-t2))
+    for i in 1:numfunctions(g)
+        Ai = 1/2*norm(cross(t.vertices[mod1(i-1,3)]-x,t.vertices[mod1(i+1,3)]-x))
+        g = Ai/Atot
+        for j in 1:numfunctions(f)
+           z[i,j] += β * dot(curlt[i],curls[j])*greenconst*dx + α*dot(nx,ny) * greenlinear[j]*g*dx
+        end
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index f5a0e4c4..d2b96799 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -50,6 +50,7 @@ include("test_dipole.jl")
 
 include("test_wiltonints.jl")
 include("test_sauterschwabints.jl")
+include("test_hh3dints.jl")
 include("test_ss_nested_meshes.jl")
 include("test_nitsche.jl")
 include("test_nitschehh3d.jl")
diff --git a/test/test_hh3dints.jl b/test/test_hh3dints.jl
new file mode 100644
index 00000000..a441a647
--- /dev/null
+++ b/test/test_hh3dints.jl
@@ -0,0 +1,769 @@
+using Test
+using LinearAlgebra
+
+using BEAST, CompScienceMeshes, SauterSchwabQuadrature, StaticArrays
+
+T=Float64
+
+
+# operators
+op1 = Helmholtz3D.singlelayer(gamma=T(1.0)im+T(1.0))
+op2 = Helmholtz3D.doublelayer(gamma=T(1.0)im+T(1.0))
+op3 = Helmholtz3D.doublelayer_transposed(gamma=T(1.0)im+T(1.0))
+op4 = Helmholtz3D.hypersingular(gamma=T(1.0)im+T(1.0))
+## Common face case:
+@testset "Common Face" begin
+  # triangles
+t1 = simplex(
+  T.(@SVector [0.180878, -0.941848, -0.283207]),
+  T.(@SVector [0.0, -0.980785, -0.19509]),
+  T.(@SVector [0.0, -0.92388, -0.382683]))
+
+  lag0 = BEAST.LagrangeRefSpace{T,0,3,1}() # patch basis
+  lag1 = BEAST.LagrangeRefSpace{T,1,3,3}() #pyramid basis
+
+  tqd0 = BEAST.quadpoints(lag0, [t1], (13,)) #test quadrature data
+  tqd1 = BEAST.quadpoints(lag1, [t1], (13,))
+
+  bqd0 = BEAST.quadpoints(lag0, [t1], (12,)) #basis quadrature data
+  bqd1 = BEAST.quadpoints(lag1, [t1], (12,))
+
+  SE_strategy = BEAST.WiltonSERule( #wilton
+    tqd0[1,1],
+    BEAST.DoubleQuadRule(
+    tqd0[1,1],
+    bqd0[1,1],
+    ),
+  )
+
+  SS_strategy = SauterSchwabQuadrature.CommonFace(BEAST._legendre(12,T(0.0),T(1.0))) #sauter
+#single layer with patch-patch
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd0[1,1],
+    bqd0[1,1]
+  )
+
+  z_se = [0.0im]
+  z_ss = [0.0im]
+  z_double = [0.0im]
+
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t1, z_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t1, z_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t1, z_double, Double_strategy)
+
+  @test z_se ≈ z_ss rtol=1e-5
+  @test z_double ≈ z_ss rtol = 1e-1
+  @test norm(z_se-z_ss) < norm(z_double-z_ss)
+
+  #single layer with patch-pyramid
+  SE_strategy = BEAST.WiltonSERule( #wilton
+  tqd1[1,1],
+  BEAST.DoubleQuadRule(
+  tqd1[1,1],
+  bqd0[1,1],
+  ),)
+  Double_strategy = BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd0[1,1],
+    )
+
+  z_se = zeros(complex(T),3,1)
+  z_ss = zeros(complex(T),3,1)
+  z_double = zeros(complex(T),3,1)
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t1, z_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t1, z_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t1, z_double, Double_strategy)
+
+  @test z_ss ≈ z_se rtol = 1e-5
+  @test z_ss ≈ z_double rtol = 1e-1
+  @test norm(z_ss-z_se) < norm(z_ss-z_double)
+
+#single layer with pyramid-patch
+  SE_strategy = BEAST.WiltonSERule( #wilton
+  tqd0[1,1],
+  BEAST.DoubleQuadRule(
+  tqd0[1,1],
+  bqd1[1,1],
+  ),)
+  Double_strategy = BEAST.DoubleQuadRule(
+    tqd0[1,1],
+    bqd1[1,1],
+    )
+
+  z_se = zeros(complex(T),1,3)
+  z_ss = zeros(complex(T),1,3)
+  z_double = zeros(complex(T),1,3)
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t1, z_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t1, z_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t1, z_double, Double_strategy)
+
+  @test z_ss ≈ z_se rtol = 1e-5
+  @test z_ss ≈ z_double rtol = 1e-1
+  @test norm(z_ss-z_se) < norm(z_ss-z_double)
+
+  #single layer with pyramid pyramid
+  SE_strategy = BEAST.WiltonSERule( #wilton
+  tqd1[1,1],
+  BEAST.DoubleQuadRule(
+  tqd1[1,1],
+  bqd1[1,1],
+  ),)
+  Double_strategy = BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd1[1,1],
+    )
+
+  z_se = zeros(complex(T),3,3)
+  z_ss = zeros(complex(T),3,3)
+  z_double = zeros(complex(T),3,3)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t1, z_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t1, z_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t1, z_double, Double_strategy)
+
+  @test z_ss ≈ z_se rtol = 1e-5
+  @test z_ss ≈ z_double rtol = 1e-1
+  @test norm(z_ss-z_se) < norm(z_ss-z_double)
+
+end
+# common vertex case
+
+# In the common vertex case, for the single layer and hypersingular,
+# the double quadrature apparently is better than the wilton method.
+# It is not clear to me why this is. If someone knows why this happens
+# I would be grateful for an explanation.
+
+@testset "Common vertex" begin
+  t1 = simplex(
+    T.(@SVector [0.180878, -0.941848, -0.283207]),
+    T.(@SVector [0.0, -0.980785, -0.19509]),
+    T.(@SVector [0.0, -0.92388, -0.382683]))
+  t2 = simplex(
+    T.(@SVector [0.373086, -0.881524, -0.289348]),
+    T.(@SVector [0.180878, -0.941848, -0.283207]),
+    T.(@SVector [0.294908, -0.944921, -0.141962]))
+  lag0 = BEAST.LagrangeRefSpace{T,0,3,1}() # patch basis
+  lag1 = BEAST.LagrangeRefSpace{T,1,3,3}() #pyramid basis
+
+  tqd0 = BEAST.quadpoints(lag0, [t1,t2], (12,)) #test quadrature data
+  tqd1 = BEAST.quadpoints(lag1, [t1,t2], (12,))
+
+  bqd0 = BEAST.quadpoints(lag0, [t1,t2], (13,)) #basis quadrature data
+  bqd1 = BEAST.quadpoints(lag1, [t1,t2], (13,))
+#patch patch
+  SE_strategy = BEAST.WiltonSERule(
+    tqd0[1,1],
+    BEAST.DoubleQuadRule(
+    tqd0[1,1],
+    bqd0[1,2]))
+
+  SS_strategy = SauterSchwabQuadrature.CommonVertex(BEAST._legendre(8,T(0.0),T(1.0)))
+
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd0[1,1],
+    bqd0[1,2]
+  )
+  #single layer
+
+  z_cv_se = zeros(complex(T),1,1)
+  z_cv_ss = zeros(complex(T),1,1)
+  z_cv_double = zeros(complex(T),1,1)
+
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_ss rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se)/norm(z_cv_ss) < norm(z_cv_ss-z_cv_double)/norm(z_cv_ss) broken = true
+  #doublelayer
+  z_cv_se = zeros(complex(T),1,1)
+  z_cv_ss = zeros(complex(T),1,1)
+  z_cv_double = zeros(complex(T),1,1)
+
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_double rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double)
+
+  # doublelayer_transposed
+  z_cv_se = zeros(complex(T),1,1)
+  z_cv_ss = zeros(complex(T),1,1)
+  z_cv_double = zeros(complex(T),1,1)
+
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_ss rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double)
+
+#pyramid pyramid
+  # singlelayer
+  SE_strategy = BEAST.WiltonSERule(
+    tqd1[1,1],
+    BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd1[1,2]))
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd1[1,1],
+    bqd1[1,2]
+  )
+  z_cv_se = zeros(complex(T),3,3)
+  z_cv_ss = zeros(complex(T),3,3)
+  z_cv_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_double rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double) broken = true
+
+  #doublelayer 
+  z_cv_se = zeros(complex(T),3,3)
+  z_cv_ss = zeros(complex(T),3,3)
+  z_cv_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_ss rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double)
+
+  #doublelayer_transposed
+  z_cv_se = zeros(complex(T),3,3)
+  z_cv_ss = zeros(complex(T),3,3)
+  z_cv_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-3
+  @test z_cv_double ≈ z_cv_se rtol = 1e-5
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double)
+
+  #hypersingular
+  z_cv_se = zeros(complex(T),3,3)
+  z_cv_ss = zeros(complex(T),3,3)
+  z_cv_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_ss rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double) broken = true
+
+#pyramid test patch basis
+  SE_strategy = BEAST.WiltonSERule(
+    tqd1[1,1],
+    BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd0[1,2]))
+
+  SS_strategy = SauterSchwabQuadrature.CommonVertex(BEAST._legendre(8,T(0.0),T(1.0)))
+
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd1[1,1],
+    bqd0[1,2]
+  )
+  #singlelayer
+  z_cv_se = zeros(complex(T),3,1)
+  z_cv_ss = zeros(complex(T),3,1)
+  z_cv_double = zeros(complex(T),3,1)
+
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_ss rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double) broken = true
+
+  #doublelayer
+  z_cv_se = zeros(complex(T),3,1)
+  z_cv_ss = zeros(complex(T),3,1)
+  z_cv_double = zeros(complex(T),3,1)
+
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_ss rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se)/norm(z_cv_ss) < norm(z_cv_ss-z_cv_double)/norm(z_cv_ss)
+  #doublelayer_transposed
+
+  z_cv_se = zeros(complex(T),3,1)
+  z_cv_ss = zeros(complex(T),3,1)
+  z_cv_double = zeros(complex(T),3,1)
+
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_ss rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double)
+
+  #pyramid basis patch test
+  SE_strategy = BEAST.WiltonSERule(
+    tqd0[1,1],
+    BEAST.DoubleQuadRule(
+    tqd0[1,1],
+    bqd1[1,2]))
+
+  SS_strategy = SauterSchwabQuadrature.CommonVertex(BEAST._legendre(8,T(0.0),T(1.0)))
+
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd0[1,1],
+    bqd1[1,2]
+  )
+  #singlelayer
+  z_cv_se = zeros(complex(T),1,3)
+  z_cv_ss = zeros(complex(T),1,3)
+  z_cv_double = zeros(complex(T),1,3)
+
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_ss rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double) broken = true
+
+  #double layer
+  z_cv_se = zeros(complex(T),1,3)
+  z_cv_ss = zeros(complex(T),1,3)
+  z_cv_double = zeros(complex(T),1,3)
+
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_ss rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double)
+
+  #double layer transposed
+  z_cv_se = zeros(complex(T),1,3)
+  z_cv_ss = zeros(complex(T),1,3)
+  z_cv_double = zeros(complex(T),1,3)
+
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_cv_se, SE_strategy)
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_cv_ss, SS_strategy)
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_cv_double, Double_strategy)
+
+  @test z_cv_double ≈ z_cv_ss rtol = 1e-4
+  @test z_cv_se ≈ z_cv_ss rtol=1e-4
+  @test norm(z_cv_ss-z_cv_se) < norm(z_cv_ss-z_cv_double)
+
+end
+
+@testset "Common edge" begin
+  # common edge
+  t1 = simplex(
+      T.(@SVector [0.180878, -0.941848, -0.283207]),
+      T.(@SVector [0.0, -0.980785, -0.19509]),
+      T.(@SVector [0.0, -0.92388, -0.382683])
+      )
+  t2 = simplex(
+      T.(@SVector [0.158174, -0.881178, -0.44554]),
+      T.(@SVector [0.180878, -0.941848, -0.283207]),
+      T.(@SVector [0.0, -0.92388, -0.382683])
+      )
+
+  lag0 = BEAST.LagrangeRefSpace{T,0,3,1}() # patch basis
+  lag1 = BEAST.LagrangeRefSpace{T,1,3,3}() #pyramid basis
+
+  tqd0 = BEAST.quadpoints(lag0, [t1,t2], (12,)) #test quadrature data
+  tqd1 = BEAST.quadpoints(lag1, [t1,t2], (12,))
+
+  bqd0 = BEAST.quadpoints(lag0, [t1,t2], (13,)) #basis quadrature data
+  bqd1 = BEAST.quadpoints(lag1, [t1,t2], (13,))
+
+  #patch patch
+  SE_strategy = BEAST.WiltonSERule(
+    tqd0[1,1],
+    BEAST.DoubleQuadRule(
+    tqd0[1,1],
+    bqd0[1,2]))
+
+  SS_strategy = SauterSchwabQuadrature.CommonEdge(BEAST._legendre(8,T(0.0),T(1.0)))
+
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd0[1,1],
+    bqd0[1,2]
+  )
+  #single layer
+
+  z_ce_se = zeros(complex(T),1,1)
+  z_ce_ss = zeros(complex(T),1,1)
+  z_ce_double = zeros(complex(T),1,1)
+
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-3
+  @test z_ce_se ≈ z_ce_ss rtol=1e-4
+  @test norm(z_ce_ss-z_ce_se) < norm(z_ce_ss-z_ce_double)
+
+  #doublelayer
+  z_ce_se = zeros(complex(T),1,1)
+  z_ce_ss = zeros(complex(T),1,1)
+  z_ce_double = zeros(complex(T),1,1)
+
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_ss ≈ z_ce_double rtol = 1e-1
+  @test z_ce_se ≈ z_ce_ss rtol = 1e-2
+  @test norm(z_ce_ss-z_ce_se) < norm(z_ce_ss-z_ce_double)
+
+  # doublelayer_transposed
+  z_ce_se = zeros(complex(T),1,1)
+  z_ce_ss = zeros(complex(T),1,1)
+  z_ce_double = zeros(complex(T),1,1)
+
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-1
+  @test z_ce_se ≈ z_ce_ss rtol=1e-2
+  @test norm(z_ce_ss-z_ce_se) < norm(z_ce_ss-z_ce_double)
+
+#pyramid pyramid
+  # singlelayer
+  SE_strategy = BEAST.WiltonSERule(
+    tqd1[1,1],
+    BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd1[1,2]))
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd1[1,1],
+    bqd1[1,2]
+  )
+  z_ce_se = zeros(complex(T),3,3)
+  z_ce_ss = zeros(complex(T),3,3)
+  z_ce_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-3
+  @test z_ce_se ≈ z_ce_double rtol=1e-3
+  @test norm(z_ce_ss-z_ce_se) < norm(z_ce_ss-z_ce_double)
+
+  #doublelayer 
+  z_ce_se = zeros(complex(T),3,3)
+  z_ce_ss = zeros(complex(T),3,3)
+  z_ce_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-1
+  @test z_ce_se ≈ z_ce_ss rtol=1e-4
+  @test norm(z_ce_ss-z_ce_se)/norm(z_ce_ss) < norm(z_ce_ss-z_ce_double)/norm(z_ce_ss)
+
+  #doublelayer_transposed
+  z_ce_se = zeros(complex(T),3,3)
+  z_ce_ss = zeros(complex(T),3,3)
+  z_ce_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-1
+  @test z_ce_ss ≈ z_ce_se rtol = 1e-2
+  @test norm(z_ce_ss-z_ce_se) < norm(z_ce_ss-z_ce_double)
+
+  #hypersingular
+  z_ce_se = zeros(complex(T),3,3)
+  z_ce_ss = zeros(complex(T),3,3)
+  z_ce_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-3
+  @test z_ce_se ≈ z_ce_ss rtol=1e-4
+  @test norm(z_ce_ss-z_ce_se) < norm(z_ce_ss-z_ce_double)
+
+#pyramid test patch basis
+  SE_strategy = BEAST.WiltonSERule(
+    tqd1[1,1],
+    BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd0[1,2]))
+
+  SS_strategy = SauterSchwabQuadrature.CommonEdge(BEAST._legendre(8,T(0.0),T(1.0)))
+
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd1[1,1],
+    bqd0[1,2]
+  )
+  #singlelayer
+  z_ce_se = zeros(complex(T),3,1)
+  z_ce_ss = zeros(complex(T),3,1)
+  z_ce_double = zeros(complex(T),3,1)
+
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-3
+  @test z_ce_se ≈ z_ce_ss rtol=1e-4
+  @test norm(z_ce_ss-z_ce_se) < norm(z_ce_ss-z_ce_double)
+  
+  #doublelayer
+  z_ce_se = zeros(complex(T),3,1)
+  z_ce_ss = zeros(complex(T),3,1)
+  z_ce_double = zeros(complex(T),3,1)
+
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-1
+  @test z_ce_se ≈ z_ce_ss rtol=1e-4
+  @test norm(z_ce_ss-z_ce_se)/norm(z_ce_ss) < norm(z_ce_ss-z_ce_double)/norm(z_ce_ss)
+  #doublelayer_transposed
+  
+  z_ce_se = zeros(complex(T),3,1)
+  z_ce_ss = zeros(complex(T),3,1)
+  z_ce_double = zeros(complex(T),3,1)
+
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-1
+  @test z_ce_se ≈ z_ce_ss rtol=1e-2
+  @test norm(z_ce_ss-z_ce_se)/norm(z_ce_ss) < norm(z_ce_ss-z_ce_double)/norm(z_ce_ss)
+
+  #pyramid basis patch test
+  SE_strategy = BEAST.WiltonSERule(
+    tqd0[1,1],
+    BEAST.DoubleQuadRule(
+    tqd0[1,1],
+    bqd1[1,2]))
+
+  SS_strategy = SauterSchwabQuadrature.CommonEdge(BEAST._legendre(8,T(0.0),T(1.0)))
+
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd0[1,1],
+    bqd1[1,2]
+  )
+  #singlelayer
+  z_ce_se = zeros(complex(T),1,3)
+  z_ce_ss = zeros(complex(T),1,3)
+  z_ce_double = zeros(complex(T),1,3)
+
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-3
+  @test z_ce_se ≈ z_ce_ss rtol=1e-4
+  @test norm(z_ce_ss-z_ce_se)/norm(z_ce_ss) < norm(z_ce_ss-z_ce_double)/norm(z_ce_ss)
+
+  #double layer
+  z_ce_se = zeros(complex(T),1,3)
+  z_ce_ss = zeros(complex(T),1,3)
+  z_ce_double = zeros(complex(T),1,3)
+
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-1
+  @test z_ce_se ≈ z_ce_ss rtol=1e-4
+  @test norm(z_ce_ss-z_ce_se)/norm(z_ce_ss) < norm(z_ce_ss-z_ce_double)/norm(z_ce_ss)
+
+  #double layer transposed
+  z_ce_se = zeros(complex(T),1,3)
+  z_ce_ss = zeros(complex(T),1,3)
+  z_ce_double = zeros(complex(T),1,3)
+
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_ce_ss, SS_strategy)
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_double ≈ z_ce_ss rtol = 1e-1
+  @test z_ce_se ≈ z_ce_ss rtol=1e-1
+  @test norm(z_ce_ss-z_ce_se)/norm(z_ce_ss) < norm(z_ce_ss-z_ce_double)/norm(z_ce_ss)
+
+end
+
+@testset "Seperated triangles" begin
+  t1 = simplex(
+      T.(@SVector [0.0, -0.980785, -0.19509]),
+      T.(@SVector [0.180878, -0.941848, -0.283207]),
+      T.(@SVector [0.0, -0.92388, -0.382683]))
+  t2 = simplex(
+      T.(@SVector [0.0, -0.980785, -0.230102]),
+      T.(@SVector [0.180878, -0.941848, -0.383207]),
+      T.(@SVector [0.0, -0.92388, -0.411962]))
+
+  lag0 = BEAST.LagrangeRefSpace{T,0,3,1}() # patch basis
+  lag1 = BEAST.LagrangeRefSpace{T,1,3,3}() #pyramid basis
+  
+  tqd0 = BEAST.quadpoints(lag0, [t1,t2], (12,)) #test quadrature data
+  tqd1 = BEAST.quadpoints(lag1, [t1,t2], (12,))
+  
+  bqd0 = BEAST.quadpoints(lag0, [t1,t2], (13,)) #basis quadrature data
+  bqd1 = BEAST.quadpoints(lag1, [t1,t2], (13,))
+  # singlelayer
+  SE_strategy = BEAST.WiltonSERule(
+    tqd0[1,1],
+    BEAST.DoubleQuadRule(
+    tqd0[1,1],
+    bqd0[1,2]))
+  
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd0[1,1],
+    bqd0[1,2]
+  )
+  z_ce_se = zeros(complex(T),1,1)
+  z_ce_double = zeros(complex(T),1,1)
+  
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_ce_double, Double_strategy)
+  
+  @test z_ce_se ≈ z_ce_double rtol=1e-6
+
+  # calculate a reference value with hcubature to desired accuracy
+  #=using HCubature
+  r(u,v,t) = (1-u)*t[1]+u*((1-v)*t[2]+v*t[3])
+  gamma=1.0+1.0im
+  jac(u,v,t) = u * norm(cross(t[1],t[2])+cross(t[2],t[3])+cross(t[3],t[1]))
+  function outerf(x)
+    r2(x2) = r(x2[1],x2[2],t1) - r(x[1],x[2],t2)
+    innerf(x2) = exp(-gamma*norm(r2(x2)))/(4*pi*norm(r2(x2))) * jac(x2[1],x2[2],t1)
+    hcubature(innerf, (0.0,0.0), (1.0,1.0), rtol = 1e-8)[1] * jac(x[1],x[2],t2)
+  end
+
+  @show refval = hcubature(outerf, (0.0,0.0), (1.0,1.0), rtol = 1e-8)[1] =#
+  refval = 0.00033563892481993545 - 2.22592609372769e-5im # can be calculated with above code
+  
+  @test z_ce_se[1] ≈ refval rtol = 1e-7
+  @test z_ce_double[1] ≈ refval rtol = 1e-6
+
+
+  SE_strategy = BEAST.WiltonSERule(
+    tqd1[1,1],
+    BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd1[1,2]))
+
+  Double_strategy = BEAST.DoubleQuadRule( #doublequadstrat
+    tqd1[1,1],
+    bqd1[1,2]
+  )
+  z_ce_se = zeros(complex(T),3,3)
+  z_ce_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_se ≈ z_ce_double rtol=1e-4
+
+  # doublelayer
+  SE_strategy = BEAST.WiltonSERule(
+    tqd0[1,1],
+    BEAST.DoubleQuadRule(
+    tqd0[1,1],
+    bqd1[1,2]))
+
+  Double_strategy = BEAST.DoubleQuadRule(
+    tqd0[1,1],
+    bqd1[1,2])
+  z_ce_se = zeros(complex(T),1,3)
+  z_ce_double = zeros(complex(T),1,3)
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_ce_double, Double_strategy)
+  
+  @test z_ce_se ≈ z_ce_double rtol = 1e-4
+
+  SE_strategy = BEAST.WiltonSERule(
+    tqd1[1,1],
+    BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd1[1,2]))
+
+  Double_strategy = BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd1[1,2])
+  z_ce_se = zeros(complex(T),3,3)
+  z_ce_double = zeros(complex(T),3,3)
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_se ≈ z_ce_double rtol = 1e-4
+  # doublelayer transposed
+  SE_strategy = BEAST.WiltonSERule(
+      tqd1[1,1],
+      BEAST.DoubleQuadRule(
+          tqd1[1,1],
+          bqd0[1,2]))
+  Double_strategy = BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd0[1,2])
+  z_ce_se = zeros(complex(T),3,1)
+  z_ce_double = zeros(complex(T),3,1)
+
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_se ≈ z_ce_double rtol=1e-4
+
+  SE_strategy = BEAST.WiltonSERule(
+      tqd1[1,1],
+      BEAST.DoubleQuadRule(
+          tqd1[1,1],
+          bqd1[1,2]))
+  Double_strategy = BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd1[1,2])
+  z_ce_se = zeros(complex(T),3,3)
+  z_ce_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_ce_double, Double_strategy)
+
+  @test z_ce_se ≈ z_ce_double rtol=1e-4
+  # hypersingular
+  SE_strategy = BEAST.WiltonSERule(
+      tqd1[1,1],
+      BEAST.DoubleQuadRule(
+          tqd1[1,1],
+          bqd1[1,2]))
+  Double_strategy = BEAST.DoubleQuadRule(
+    tqd1[1,1],
+    bqd1[1,2])
+
+  z_ce_se = zeros(complex(T),3,3)
+  z_ce_double = zeros(complex(T),3,3)
+
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_ce_se, SE_strategy)
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_ce_double, Double_strategy)
+  @test z_ce_se ≈ z_ce_double rtol = 1e-5
+
+end