diff --git a/src/BEAST.jl b/src/BEAST.jl
index 7c729532..ce280bd4 100644
--- a/src/BEAST.jl
+++ b/src/BEAST.jl
@@ -185,15 +185,6 @@ include("bases/tensorbasis.jl")
 include("operator.jl")
 
 include("quadrature/quadstrats.jl")
-
-include("excitation.jl")
-include("localop.jl")
-include("multiplicativeop.jl")
-include("identityop.jl")
-include("integralop.jl")
-include("dyadicop.jl")
-include("interpolation.jl")
-
 include("quadrature/doublenumqstrat.jl")
 include("quadrature/doublenumsauterqstrat.jl")
 include("quadrature/doublenumwiltonsauterqstrat.jl")
@@ -202,12 +193,23 @@ include("quadrature/selfsauterdnumotherwiseqstrat.jl")
 include("quadrature/nonconformingintegralopqstrat.jl")
 include("quadrature/commonfaceoverlappingedgeqstrat.jl")
 include("quadrature/strategies/cfcvsautercewiltonpdnumqstrat.jl")
+include("quadrature/strategies/testrefinestrialqstrat.jl")
+
+include("excitation.jl")
+include("localop.jl")
+include("multiplicativeop.jl")
+include("identityop.jl")
+include("integralop.jl")
+include("dyadicop.jl")
+include("interpolation.jl")
 
+include("quadrature/rules/momintegrals.jl")
 include("quadrature/doublenumints.jl")
 include("quadrature/singularityextractionints.jl")
 include("quadrature/sauterschwabints.jl")
 include("quadrature/nonconformingoverlapqrule.jl")
 include("quadrature/nonconformingtouchqrule.jl")
+include("quadrature/rules/testrefinestrialqrule.jl")
 
 include("postproc.jl")
 include("postproc/segcurrents.jl")
diff --git a/src/decoupled/dpops.jl b/src/decoupled/dpops.jl
index cfae115d..502b3e6c 100644
--- a/src/decoupled/dpops.jl
+++ b/src/decoupled/dpops.jl
@@ -17,13 +17,14 @@ function integrand(op::CurlSingleLayerDP3D, kernel_vals,
     return -α * dot(nx × gx, ∇G * fy)
 end
 
-function momintegrals!(op::CurlSingleLayerDP3D,
-    test_local_space::RTRefSpace, trial_local_space::LagrangeRefSpace,
-    test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)
+function momintegrals!(out, op::CurlSingleLayerDP3D,
+    test_local_space::RTRefSpace, tptr, test_triangular_element,
+    trial_local_space::LagrangeRefSpace, bptr, trial_triangular_element,
+    qrule::SauterSchwabStrategy)
 
     I, J, K, L = SauterSchwabQuadrature.reorder(
         test_triangular_element.vertices,
-        trial_triangular_element.vertices, strat)
+        trial_triangular_element.vertices, qrule)
 
     test_triangular_element  = simplex(test_triangular_element.vertices[I]...)
     trial_triangular_element = simplex(trial_triangular_element.vertices[J]...)
@@ -55,7 +56,7 @@ function momintegrals!(op::CurlSingleLayerDP3D,
         return @SMatrix[dot(f[i].value × nx, R[j]) for i in 1:3, j in 1:3]
     end
 
-    Q = sauterschwab_parameterized(igd, strat)
+    Q = sauterschwab_parameterized(igd, qrule)
     for j ∈ 1:3
         for i ∈ 1:3
             out[i,j] += Q[K[i],L[j]]
diff --git a/src/helmholtz2d/helmholtzop.jl b/src/helmholtz2d/helmholtzop.jl
index 3e10ef0e..ee3806f1 100644
--- a/src/helmholtz2d/helmholtzop.jl
+++ b/src/helmholtz2d/helmholtzop.jl
@@ -29,7 +29,9 @@ function testfunc1()
     print("test function!")
 end
 
-defaultquadstrat(op::HelmholtzOperator2D, tfs, bfs) = DoubleNumWiltonSauterQStrat(4,3,4,3,4,4,4,4)
+# defaultquadstrat(op::HelmholtzOperator2D, tfs, bfs) = DoubleNumWiltonSauterQStrat(4,3,4,3,4,4,4,4)
+defaultquadstrat(op::HelmholtzOperator2D, tfs, bfs) = DoubleNumQStrat(4,3)
+
 
 function quaddata(op::HelmholtzOperator2D, g::LagrangeRefSpace, f::LagrangeRefSpace, tels, bels,
         qs::DoubleNumWiltonSauterQStrat)
diff --git a/src/helmholtz3d/hh3d_sauterschwabqr.jl b/src/helmholtz3d/hh3d_sauterschwabqr.jl
index 19e4ce52..8b137891 100644
--- a/src/helmholtz3d/hh3d_sauterschwabqr.jl
+++ b/src/helmholtz3d/hh3d_sauterschwabqr.jl
@@ -1,265 +1 @@
-# function pulled_back_integrand(op::HH3DSingleLayerFDBIO,
-#     test_local_space::LagrangeRefSpace,
-#     trial_local_space::LagrangeRefSpace,
-#     test_chart, trial_chart)
 
-#     (u,v) -> begin
-
-#         x = neighborhood(test_chart,u)
-#         y = neighborhood(trial_chart,v)
-
-#         f = test_local_space(x)
-#         g = trial_local_space(y)
-
-#         j = jacobian(x) * jacobian(y)
-
-#         α = op.alpha
-#         γ = gamma(op)
-#         R = norm(cartesian(x)-cartesian(y))
-#         G = exp(-γ*R)/(4*π*R)
-
-#         αjG = α*G*j
-
-#         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,
-#     trial_local_space::LagrangeRefSpace,
-#     test_chart, trial_chart)
-
-#     (u,v) -> begin
-
-#         x = neighborhood(test_chart,u)
-#         y = neighborhood(trial_chart,v)
-
-#         nx = normal(x)
-#         ny = normal(y)
-
-#         f = test_local_space(x)
-#         g = trial_local_space(y)
-
-#         j = jacobian(x) * jacobian(y)
-
-#         α = op.alpha
-#         β = op.beta
-#         γ = gamma(op)
-#         R = norm(cartesian(x)-cartesian(y))
-#         G = exp(-γ*R)/(4*π*R)
-
-#         αjG = ny*α*G*j
-#         βjG = β*G*j
-
-#         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,
-#     trial_local_space::LagrangeRefSpace,
-#     test_chart, trial_chart)
-
-#     (u,v) -> begin
-
-#         x = neighborhood(test_chart,u)
-#         y = neighborhood(trial_chart,v)
-
-#         ny = normal(y)
-#         f = test_local_space(x)
-#         g = trial_local_space(y)
-
-#         j = jacobian(x) * jacobian(y)
-
-#         α = op.alpha
-#         γ = gamma(op)
-
-#         r = cartesian(x) - cartesian(y)
-#         R = norm(r)
-#         G = exp(-γ*R)/(4*π*R)
-#         inv_R = 1/R
-#         ∇G = -(γ + inv_R) * G * inv_R * r
-#         αnyj∇G = dot(ny,-α*∇G*j)
-
-#         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,
-#     trial_local_space::LagrangeRefSpace,
-#     test_chart, trial_chart)
-
-#     (u,v) -> begin
-
-#         x = neighborhood(test_chart,u)
-#         y = neighborhood(trial_chart,v)
-
-#         nx = normal(x)
-#         f = test_local_space(x)
-#         g = trial_local_space(y)
-
-#         j = jacobian(x) * jacobian(y)
-
-#         α = op.alpha
-#         γ = gamma(op)
-
-#         r = cartesian(x) - cartesian(y)
-#         R = norm(r)
-#         G = exp(-γ*R)/(4*π*R)
-#         inv_R = 1/R
-#         ∇G = -(γ + inv_R) * G * inv_R * r
-#         αnxj∇G = dot(nx,α*∇G*j)
-
-#         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::Helmholtz3DOp,
-#     test_local_space::LagrangeRefSpace{<:Any,0}, trial_local_space::LagrangeRefSpace{<:Any,0},
-#     test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)
-
-#     I, J, K, L = SauterSchwabQuadrature.reorder(
-#         test_triangular_element.vertices,
-#         trial_triangular_element.vertices, strat)
-
-#     test_triangular_element  = simplex(
-#         test_triangular_element.vertices[I[1]],
-#         test_triangular_element.vertices[I[2]],
-#         test_triangular_element.vertices[I[3]])
-
-#     trial_triangular_element = simplex(
-#         trial_triangular_element.vertices[J[1]],
-#         trial_triangular_element.vertices[J[2]],
-#         trial_triangular_element.vertices[J[3]])
-
-#     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] * σ
-
-#     nothing
-# end
-
-
-# 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(
-#         test_triangular_element.vertices,
-#         trial_triangular_element.vertices, strat)
-
-#     test_triangular_element  = simplex(
-#         test_triangular_element.vertices[I[1]],
-#         test_triangular_element.vertices[I[2]],
-#         test_triangular_element.vertices[I[3]])
-
-#     trial_triangular_element = simplex(
-#         trial_triangular_element.vertices[J[1]],
-#         trial_triangular_element.vertices[J[2]],
-#         trial_triangular_element.vertices[J[3]])
-
-#     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)
-#     for j ∈ 1:3, i ∈ 1:3
-#         out[i,j] += G[K[i],L[j]] * σ
-#     end
-
-#     nothing
-# end
-
-# function momintegrals!(op::Helmholtz3DOp,
-#     test_local_space::LagrangeRefSpace{<:Any,0}, trial_local_space::LagrangeRefSpace{<:Any,1},
-#     test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)
-
-#     I, J, K, L = SauterSchwabQuadrature.reorder(
-#         test_triangular_element.vertices,
-#         trial_triangular_element.vertices, strat)
-
-#     test_triangular_element  = simplex(
-#         test_triangular_element.vertices[I[1]],
-#         test_triangular_element.vertices[I[2]],
-#         test_triangular_element.vertices[I[3]])
-
-#     trial_triangular_element = simplex(
-#         trial_triangular_element.vertices[J[1]],
-#         trial_triangular_element.vertices[J[2]],
-#         trial_triangular_element.vertices[J[3]])
-
-#     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)
-
-#     for i ∈ 1:3
-#         out[1,i] += G[L[i]] * σ
-#     end
-
-#     nothing
-# end
-
-# function momintegrals!(op::Helmholtz3DOp,
-#     test_local_space::LagrangeRefSpace{<:Any,1}, trial_local_space::LagrangeRefSpace{<:Any,0},
-#     test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)
-
-#     I, J, K, L = SauterSchwabQuadrature.reorder(
-#         test_triangular_element.vertices,
-#         trial_triangular_element.vertices, strat)
-
-#     test_triangular_element  = simplex(
-#         test_triangular_element.vertices[I[1]],
-#         test_triangular_element.vertices[I[2]],
-#         test_triangular_element.vertices[I[3]])
-
-#     trial_triangular_element = simplex(
-#         trial_triangular_element.vertices[J[1]],
-#         trial_triangular_element.vertices[J[2]],
-#         trial_triangular_element.vertices[J[3]])
-
-#         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)
-
-#     for i ∈ 1:3
-#         out[i,1] += G[K[i]] * σ
-#     end
-
-#     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/integralop.jl b/src/integralop.jl
index 289d02b3..0610b475 100644
--- a/src/integralop.jl
+++ b/src/integralop.jl
@@ -1,4 +1,4 @@
-abstract type IntegralOperator <: Operator end
+
 
 defaultquadstrat(op::IntegralOperator, tfs::RefSpace, bfs::RefSpace) =
     DoubleNumSauterQstrat(2,3,5,5,4,3)
@@ -69,9 +69,6 @@ Computes the matrix of operator biop wrt the finite element spaces tfs and bfs
 function assemblechunk!(biop::IntegralOperator, tfs::Space, bfs::Space, store;
         quadstrat=defaultquadstrat(biop, tfs, bfs))
 
-    # test_elements, tad, tcells = assemblydata(tfs)
-    # bsis_elements, bad, bcells = assemblydata(bfs)
-
     tr = assemblydata(tfs); tr == nothing && return
     br = assemblydata(bfs); br == nothing && return
 
@@ -88,10 +85,15 @@ function assemblechunk!(biop::IntegralOperator, tfs::Space, bfs::Space, store;
     zlocal = zeros(scalartype(biop, tfs, bfs), 2num_tshapes, 2num_bshapes)
 
     if CompScienceMeshes.refines(tgeo, bgeo)
-        assemblechunk_body_test_refines_trial!(biop,
+        # assemblechunk_body_test_refines_trial!(biop,
+        #     tfs, test_elements, tad, tcells,
+        #     bfs, bsis_elements, bad, bcells,
+        #     qd, zlocal, store; quadstrat)
+        qs = TestRefinesTrialQStrat(quadstrat)
+        assemblechunk_body!(biop,
             tfs, test_elements, tad, tcells,
             bfs, bsis_elements, bad, bcells,
-            qd, zlocal, store; quadstrat)
+            qd, zlocal, store; quadstrat=qs)
     elseif CompScienceMeshes.refines(bgeo, tgeo)
         assemblechunk_body_trial_refines_test!(biop,
             tfs, test_elements, tad, tcells,
@@ -99,28 +101,32 @@ function assemblechunk!(biop::IntegralOperator, tfs::Space, bfs::Space, store;
             qd, zlocal, store; quadstrat)
     else
         assemblechunk_body!(biop,
-            tshapes, test_elements, tad,
-            bshapes, bsis_elements, bad,
+            tfs, test_elements, tad, tcells,
+            bfs, bsis_elements, bad, bcells,
             qd, zlocal, store; quadstrat)
     end
 end
 
 
 function assemblechunk_body!(biop,
-        test_shapes, test_elements, test_assembly_data,
-        trial_shapes, trial_elements, trial_assembly_data,
+        test_space, test_elements, test_assembly_data, test_cell_ptrs,
+        trial_space, trial_elements, trial_assembly_data, trial_cell_ptrs,
         qd, zlocal, store; quadstrat)
 
+    test_shapes = refspace(test_space)
+    trial_shapes = refspace(trial_space)
+
     myid = Threads.threadid()
     myid == 1 && print("dots out of 10: ")
     todo, done, pctg = length(test_elements), 0, 0
-    for (p,tcell) in enumerate(test_elements)
-        for (q,bcell) in enumerate(trial_elements)
+    for (p,(tcell,tptr)) in enumerate(zip(test_elements, test_cell_ptrs))
+        for (q,(bcell,bptr)) in enumerate(zip(trial_elements, trial_cell_ptrs))
 
         fill!(zlocal, 0)
         qrule = quadrule(biop, test_shapes, trial_shapes, p, tcell, q, bcell, qd, quadstrat)
-        # @show "integralop" qrule
-        momintegrals!(biop, test_shapes, trial_shapes, tcell, bcell, zlocal, qrule)
+        momintegrals!(zlocal, biop,
+            test_space,  tptr, tcell,
+            trial_space, bptr, bcell, qrule)
         I = length(test_assembly_data[p])
         J = length(trial_assembly_data[q])
         for j in 1 : J, i in 1 : I
@@ -342,7 +348,9 @@ function assembleblock_body!(biop::IntegralOperator,
 
             fill!(zlocals[Threads.threadid()], 0)
             qrule = quadrule(biop, test_shapes, trial_shapes, p, tcell, q, bcell, quadrature_data, quadstrat)
-            momintegrals!(biop, test_shapes, trial_shapes, tcell, bcell, zlocals[Threads.threadid()], qrule)
+            momintegrals!(zlocals[Threads.threadid()], biop,
+                tfs,  nothing, tcell,
+                bfs, nothing, bcell, qrule)
 
             for j in 1 : size(zlocals[Threads.threadid()],2)
                 for i in 1 : size(zlocals[Threads.threadid()],1)
@@ -531,14 +539,14 @@ function assemblerow!(biop::IntegralOperator, test_functions::Space, trial_funct
 
     assemblerow_body!(biop,
         test_functions, test_elements, test_shapes,
-        trial_assembly_data, trial_elements, trial_shapes,
+        trial_assembly_data, trial_functions, trial_elements, trial_shapes,
         zlocal, quadrature_data, store; quadstrat)
 end
 
 
 function assemblerow_body!(biop,
     test_functions, test_elements, test_shapes,
-    trial_assembly_data, trial_elements, trial_shapes,
+    trial_assembly_data, trial_functions, trial_elements, trial_shapes,
     zlocal, quadrature_data, store; quadstrat)
 
     test_function = test_functions.fns[1]
@@ -551,7 +559,10 @@ function assemblerow_body!(biop,
 
             fill!(zlocal, 0)
             qrule = quadrule(biop, test_shapes, trial_shapes, p, tcell, q, bcell, quadrature_data, quadstrat)
-            momintegrals!(biop, test_shapes, trial_shapes, tcell, bcell, zlocal, qrule)
+            momintegrals!(zlocal, biop,
+                test_functions, nothing, tcell,
+                trial_functions, nothing, bcell,
+                qrule)
 
             for j in 1:size(zlocal,2)
                 for (n,b) in trial_assembly_data[q,j]
@@ -580,14 +591,14 @@ function assemblecol!(biop::IntegralOperator, test_functions::Space, trial_funct
     @assert numfunctions(trial_functions) == 1
 
     assemblecol_body!(biop,
-        test_assembly_data, test_elements,  test_shapes,
+        test_assembly_data, test_functions, test_elements,  test_shapes,
         trial_functions,   trial_elements, trial_shapes,
         zlocal, quadrature_data, store; quadstrat)
 end
 
 
 function assemblecol_body!(biop,
-    test_assembly_data, test_elements, test_shapes,
+    test_assembly_data, test_functions, test_elements, test_shapes,
     trial_functions, trial_elements, trial_shapes,
     zlocal, quadrature_data, store; quadstrat)
 
@@ -602,7 +613,9 @@ function assemblecol_body!(biop,
 
             fill!(zlocal, 0)
             qrule = quadrule(biop, test_shapes, trial_shapes, p, tcell, q, bcell, quadrature_data, quadstrat)
-            momintegrals!(biop, test_shapes, trial_shapes, tcell, bcell, zlocal, qrule)
+            momintegrals!(zlocal, biop,
+                test_functions, nothing, tcell,
+                trial_functions, nothing, bcell, qrule)
 
             for i in 1:size(zlocal,1)
                 for (m,a) in test_assembly_data[p,i]
diff --git a/src/maxwell/bogaertints.jl b/src/maxwell/bogaertints.jl
index d5f08e49..9307bb25 100644
--- a/src/maxwell/bogaertints.jl
+++ b/src/maxwell/bogaertints.jl
@@ -1,4 +1,7 @@
-function momintegrals!(op::MWSingleLayer3D, g::RTRefSpace, f::RTRefSpace, t, s, z, strat::BogaertStrategy)
+function momintegrals!(z, op::MWSingleLayer3D,
+    g::RTRefSpace, tptr, t,
+    f::RTRefSpace, bptr, s,
+    strat::BogaertStrategy)
 
     T, GG  = GetIntegrals(t, s, op.gamma, strat)
 
@@ -52,8 +55,10 @@ end
 
 
 
-function momintegrals!(op::MWDoubleLayer3D, g::RTRefSpace, f::RTRefSpace,
-    τ, σ, z, strat::BogaertStrategy)
+function momintegrals!(z, op::MWDoubleLayer3D,
+    g::RTRefSpace, tptr, τ,
+    f::RTRefSpace, bptr, σ,
+    strat::BogaertStrategy)
 
     # Get the primitives
     r = τ.vertices
diff --git a/src/operator.jl b/src/operator.jl
index e1958954..b54b2a74 100644
--- a/src/operator.jl
+++ b/src/operator.jl
@@ -6,13 +6,11 @@ struct Threading{T} end
 import Base: transpose, +, -, *
 
 abstract type AbstractOperator end
-
-#@linearspace AbstractOperator{T} T
-
 """
 *Atomic operator*: one that assemblechunk can deal with
 """
 abstract type Operator <: AbstractOperator end
+abstract type IntegralOperator <: Operator end
 
 mutable struct TransposedOperator <: AbstractOperator
     op::AbstractOperator
diff --git a/src/quadrature/doublenumints.jl b/src/quadrature/doublenumints.jl
index 1efdf398..0efc2ccb 100644
--- a/src/quadrature/doublenumints.jl
+++ b/src/quadrature/doublenumints.jl
@@ -9,7 +9,8 @@ momintegrals!(biop, tshs, bshs, tcell, bcell, interactions, strat)
 
 Function for the computation of moment integrals using simple double quadrature.
 """
-function momintegrals!(biop, tshs, bshs, tcell, bcell, z, strat::DoubleQuadRule)
+function momintegrals!(biop,
+    tshs, bshs, tcell, bcell, z, strat::DoubleQuadRule)
 
     igd = Integrand(biop, tshs, bshs, tcell, bcell)
 
diff --git a/src/quadrature/doublenumqstrat.jl b/src/quadrature/doublenumqstrat.jl
index 2c940e31..90a6d766 100644
--- a/src/quadrature/doublenumqstrat.jl
+++ b/src/quadrature/doublenumqstrat.jl
@@ -1,7 +1,10 @@
 function quaddata(operator::IntegralOperator,
-    local_test_basis::RefSpace, local_trial_basis::RefSpace,
+    local_test_basis, local_trial_basis,
     test_elements, trial_elements, qs::DoubleNumQStrat)
 
+    # local_test_basis = refspace(test_basis)
+    # local_trial_basis = refspace(trial_basis)
+
     test_quad_data  = quadpoints(local_test_basis,  test_elements,  (qs.outer_rule,))
     trial_quad_data = quadpoints(local_trial_basis, trial_elements, (qs.inner_rule,))
 
@@ -10,7 +13,7 @@ end
 
 
 function quadrule(operator::IntegralOperator,
-    local_test_basis::RefSpace, local_trial_basis::RefSpace,
+    local_test_basis, local_trial_basis,
     test_id, test_element, trial_id, trial_element,
     quad_data, qs::DoubleNumQStrat)
 
diff --git a/src/quadrature/doublenumwiltonsauterqstrat.jl b/src/quadrature/doublenumwiltonsauterqstrat.jl
index 2cdde450..70918398 100644
--- a/src/quadrature/doublenumwiltonsauterqstrat.jl
+++ b/src/quadrature/doublenumwiltonsauterqstrat.jl
@@ -1,8 +1,10 @@
 function quaddata(op::IntegralOperator,
-    test_local_space::RefSpace, trial_local_space::RefSpace,
+    test_local_space, trial_local_space,
     test_charts, trial_charts, qs::DoubleNumWiltonSauterQStrat)
 
     T = coordtype(test_charts[1])
+    # test_local_space = refspace(test_space)
+    # trial_local_space = refspace(trial_space)
 
     tqd = quadpoints(test_local_space,  test_charts,  (qs.outer_rule_far,qs.outer_rule_near))
     bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule_far,qs.inner_rule_near))
@@ -16,7 +18,7 @@ function quaddata(op::IntegralOperator,
 end
 
 
-function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace,  i, τ, j, σ, qd,
+function quadrule(op::IntegralOperator, g, f,  i, τ, j, σ, qd,
     qs::DoubleNumWiltonSauterQStrat)
 
     T = eltype(eltype(τ.vertices))
diff --git a/src/quadrature/nonconformingoverlapqrule.jl b/src/quadrature/nonconformingoverlapqrule.jl
index 3e2e3a26..9276b7ce 100644
--- a/src/quadrature/nonconformingoverlapqrule.jl
+++ b/src/quadrature/nonconformingoverlapqrule.jl
@@ -3,7 +3,8 @@ struct NonConformingOverlapQRule{S}
 end
 
 
-function momintegrals!(op, test_local_space, basis_local_space,
+function momintegrals!(op,
+    test_local_space, basis_local_space,
     test_chart::CompScienceMeshes.Simplex, basis_chart::CompScienceMeshes.Simplex,
     out, qrule::NonConformingOverlapQRule)
 
@@ -22,12 +23,17 @@ function momintegrals!(op, test_local_space, basis_local_space,
     test_charts = [ch for ch in test_charts if volume(ch) .> 1e6 * eps(T)]
     bsis_charts = [ch for ch in bsis_charts if volume(ch) .> 1e6 * eps(T)]
 
+    isempty(test_charts) && return
+    isempty(bsis_charts) && return
+
     # @assert volume(test_chart) ≈ sum(volume.(test_charts))
     # if volume(basis_chart) ≈ sum(volume.(bsis_charts)) else
     #     @show volume(basis_chart)
     #     @show sum(volume.(bsis_charts))
     #     error()
     # end
+    # test_local_space = refspace(test_functions)
+    # basis_local_space = refspace(basis_functions)
 
     qstrat = CommonFaceOverlappingEdgeQStrat(qrule.conforming_qstrat)
     qdata = quaddata(op, test_local_space, basis_local_space,
@@ -43,8 +49,9 @@ function momintegrals!(op, test_local_space, basis_local_space,
             Q = restrict(basis_local_space, basis_chart, bchart)
             zlocal = zero(out)
 
-            momintegrals!(op, test_local_space, basis_local_space,
-                tchart, bchart, zlocal, qrule)
+            momintegrals!(zlocal, op,
+                test_local_space, nothing, tchart,
+                basis_local_space, nothing, bchart, qrule)
 
             for i in axes(P,1)
                 for j in axes(Q,1)
diff --git a/src/quadrature/nonconformingtouchqrule.jl b/src/quadrature/nonconformingtouchqrule.jl
index d0c13696..3c31236d 100644
--- a/src/quadrature/nonconformingtouchqrule.jl
+++ b/src/quadrature/nonconformingtouchqrule.jl
@@ -4,9 +4,13 @@ struct NonConformingTouchQRule{S}
     bsis_overlapping_edge_index::Int
 end
 
-function momintegrals!(op, test_locspace, bsis_locspace,
-        τ::CompScienceMeshes.Simplex, σ::CompScienceMeshes.Simplex,
-        out, qrule::NonConformingTouchQRule)
+function momintegrals!(op,
+    test_locspace, bsis_locspace,
+    τ::CompScienceMeshes.Simplex, σ::CompScienceMeshes.Simplex,
+    out, qrule::NonConformingTouchQRule)
+
+    # test_locspace = refspace(test_functions)
+    # bsis_locspace = refspace(bsis_functions)
 
     T = coordtype(τ)
     P = eltype(τ.vertices)
@@ -58,8 +62,9 @@ function momintegrals!(op, test_locspace, bsis_locspace,
             Q = restrict(bsis_locspace, σ, bchart)
             zlocal = zero(out)
 
-            momintegrals!(op, test_locspace, bsis_locspace,
-                tchart, bchart, zlocal, qrule)
+            momintegrals!(zlocal, op,
+                test_locspace, nothing, tchart,
+                bsis_locspace, nothing, bchart, qrule)
 
             # out .+= P * zlocal * Q'
             for i in axes(P,1)
diff --git a/src/quadrature/rules/momintegrals.jl b/src/quadrature/rules/momintegrals.jl
new file mode 100644
index 00000000..6b749033
--- /dev/null
+++ b/src/quadrature/rules/momintegrals.jl
@@ -0,0 +1,13 @@
+function momintegrals!(out, op,
+    test_functions::Space, test_cellptr, test_chart,
+    trial_functions::Space, trial_cellptr, trial_chart,
+    quadrule)
+
+    local_test_space = refspace(test_functions)
+    local_trial_space = refspace(trial_functions)
+
+    momintegrals!(op,
+        local_test_space, local_trial_space,
+        test_chart, trial_chart,
+        out, quadrule)
+end
\ No newline at end of file
diff --git a/src/quadrature/rules/testrefinestrialqrule.jl b/src/quadrature/rules/testrefinestrialqrule.jl
new file mode 100644
index 00000000..08bcec5a
--- /dev/null
+++ b/src/quadrature/rules/testrefinestrialqrule.jl
@@ -0,0 +1,38 @@
+struct TestRefinesTrialQRule{S}
+    conforming_qstrat::S
+end
+
+function momintegrals!(out, op,
+    test_functions::Space, test_cell, test_chart,
+    trial_functions::Space, trial_cell, trial_chart,
+    qr::TestRefinesTrialQRule)
+
+    test_local_space = refspace(test_functions)
+    trial_local_space = refspace(trial_functions)
+
+    test_mesh = geometry(test_functions)
+    trial_mesh = geometry(trial_functions)
+
+    parent_mesh = CompScienceMeshes.parent(test_mesh)
+    trial_charts = [chart(test_mesh, p) for p in CompScienceMeshes.children(parent_mesh, trial_cell)]
+
+    quadstrat = qr.conforming_qstrat
+    qd = quaddata(op, test_local_space, trial_local_space,
+        [test_chart], trial_charts, quadstrat)
+
+    for (q,chart) in enumerate(trial_charts)
+        qr = quadrule(op, test_local_space, trial_local_space,
+            1, test_chart, q ,chart, qd, quadstrat)
+            
+        Q = restrict(trial_local_space, trial_chart, chart)
+        zlocal = zero(out)
+
+        momintegrals!(zlocal, op,
+            test_functions, nothing, test_chart,
+            trial_functions, nothing, chart, qr)
+
+        for j in 1:numfunctions(trial_local_space)
+            for i in 1:numfunctions(test_local_space)
+                for k in 1:size(Q, 2)
+                    out[i,j] += zlocal[i,k] * Q[j,k]
+end end end end end
\ No newline at end of file
diff --git a/src/quadrature/sauterschwabints.jl b/src/quadrature/sauterschwabints.jl
index c2c97fd2..ab3da50b 100644
--- a/src/quadrature/sauterschwabints.jl
+++ b/src/quadrature/sauterschwabints.jl
@@ -111,8 +111,12 @@ end
 
 
 function momintegrals!(op::Operator,
-    test_local_space::RefSpace, trial_local_space::RefSpace,
-    test_chart, trial_chart, out, rule::SauterSchwabStrategy)
+    test_local_space, trial_local_space,
+    test_chart, trial_chart,
+    out, rule::SauterSchwabStrategy)
+
+    # test_local_space = refspace(test_space)
+    # trial_local_space = refspace(trial_space)
 
     I, J, _, _ = SauterSchwabQuadrature.reorder(
         vertices(test_chart),
@@ -153,11 +157,12 @@ function momintegrals_test_refines_trial!(out, op,
     trial_functions, trial_cell, trial_chart,
     quadrule, quadstrat)
 
-    test_local_space = refspace(test_functions)
-    trial_local_space = refspace(trial_functions)
+    # test_local_space = refspace(test_functions)
+    # trial_local_space = refspace(trial_functions)
 
-    momintegrals!(op, test_local_space, trial_local_space,
-        test_chart, trial_chart, out, quadrule)
+    momintegrals!(out, op,
+        test_functions, test_cell, test_chart,
+        trial_functions, trial_cell, trial_chart, quadrule)
 end
 
 # const MWOperator3D = Union{MWSingleLayer3D, MWDoubleLayer3D}
@@ -170,7 +175,7 @@ function momintegrals_test_refines_trial!(out, op,
     trial_local_space = refspace(trial_functions)
 
     test_mesh = geometry(test_functions)
-    trial_mesh = geometry(trial_functions)
+    # trial_mesh = geometry(trial_functions)
 
     parent_mesh = CompScienceMeshes.parent(test_mesh)
     trial_charts = [chart(test_mesh, p) for p in CompScienceMeshes.children(parent_mesh, trial_cell)]
@@ -186,8 +191,9 @@ function momintegrals_test_refines_trial!(out, op,
         Q = restrict(trial_local_space, trial_chart, chart)
         zlocal = zero(out)
 
-        momintegrals!(op, test_local_space, trial_local_space,
-            test_chart, chart, zlocal, qr)
+        momintegrals!(zlocal, op,
+            test_functions, test_cell, test_chart,
+            trial_functions, nothing, chart, qr)
 
         for j in 1:numfunctions(trial_local_space)
             for i in 1:numfunctions(test_local_space)
@@ -205,8 +211,9 @@ function momintegrals_trial_refines_test!(out, op,
     test_local_space = refspace(test_functions)
     trial_local_space = refspace(trial_functions)
 
-    momintegrals!(op, test_local_space, trial_local_space,
-        test_chart, trial_chart, out, quadrule)
+    momintegrals!(out, op,
+        test_functions, test_cell, test_chart,
+        trial_functions, trial_cell, trial_chart, quadrule)
 end
 
 
@@ -233,8 +240,9 @@ function momintegrals_trial_refines_test!(out, op,
 
         Q = restrict(test_local_space, test_chart, chart)
         zlocal = zero(out)
-        momintegrals!(op, test_local_space, trial_local_space,
-            chart, trial_chart, zlocal, qr)
+        momintegrals!(zlocal, op,
+            test_functions, nothing, chart,
+            trial_functions, trial_cell, trial_chart, qr)
 
         for j in 1:numfunctions(trial_local_space)
             for i in 1:numfunctions(test_local_space)
diff --git a/src/quadrature/singularityextractionints.jl b/src/quadrature/singularityextractionints.jl
index 10c4685f..6cb3377b 100644
--- a/src/quadrature/singularityextractionints.jl
+++ b/src/quadrature/singularityextractionints.jl
@@ -1,7 +1,9 @@
 abstract type SingularityExtractionRule end
 regularpart_quadrule(qr::SingularityExtractionRule) = qr.regularpart_quadrule
 
-function momintegrals!(op, g, f, t, s, z, qrule::SingularityExtractionRule)
+function momintegrals!(op,
+    g, f,t, s,
+    z, qrule::SingularityExtractionRule)
 
     womps = qrule.outer_quad_points
 
diff --git a/src/quadrature/strategies/testrefinestrialqstrat.jl b/src/quadrature/strategies/testrefinestrialqstrat.jl
new file mode 100644
index 00000000..3e4e34e9
--- /dev/null
+++ b/src/quadrature/strategies/testrefinestrialqstrat.jl
@@ -0,0 +1,28 @@
+struct TestRefinesTrialQStrat{S}
+    conforming_qstrat::S
+end
+
+function quaddata(a, X, Y, tels, bels, qs::TestRefinesTrialQStrat)
+    return quaddata(a, X, Y, tels, bels, qs.conforming_qstrat)
+end
+
+function quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd,
+    qs::TestRefinesTrialQStrat)
+
+    hits = _numhits(τ, σ)
+    if hits > 0
+        return TestRefinesTrialQRule(qs.conforming_qstrat)
+    end
+
+    # if CompScienceMeshes.overlap(τ, σ)
+    #     return NonConformingOverlapQRule(qs.conforming_qstrat)
+    # end
+
+    # for (i,λ) in pairs(faces(τ))
+    #     for (j,μ) in pairs(faces(σ))
+    #         if CompScienceMeshes.overlap(λ, μ)
+    #             return NonConformingOverlapQRule(qs.conforming_qstrat)
+    # end end end
+
+    return quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)    
+end
\ No newline at end of file
diff --git a/src/volumeintegral/sauterschwab_ints.jl b/src/volumeintegral/sauterschwab_ints.jl
index 31164385..68385e7a 100644
--- a/src/volumeintegral/sauterschwab_ints.jl
+++ b/src/volumeintegral/sauterschwab_ints.jl
@@ -101,9 +101,10 @@ function reorder_dof(space::LagrangeRefSpace{T,1,4,4},I) where T
     return SVector(K),SVector{4,Int64}(1,1,1,1)
 end
 
-function momintegrals!(op::VIEOperator,
-    test_local_space::RefSpace, trial_local_space::RefSpace,
-    test_tetrahedron_element, trial_tetrahedron_element, out, strat::SauterSchwab3DStrategy)
+function momintegrals!(out, op::VIEOperator,
+    test_local_space::RefSpace, test_ptr, test_tetrahedron_element,
+    trial_local_space::RefSpace, trial_ptr, trial_tetrahedron_element,
+    strat::SauterSchwab3DStrategy)
 
     #Find permutation of vertices to match location of singularity to SauterSchwab
     J, I= SauterSchwab3D.reorder(strat.sing)
@@ -159,7 +160,10 @@ function momintegrals!(op::VIEOperator,
     nothing
 end
 
-function momintegrals!(biop::VIEOperator, tshs, bshs, tcell, bcell, z, strat::DoubleQuadRule)
+function momintegrals!(z, biop::VIEOperator,
+    tshs, tptr, tcell,
+    bshs, bptr, bcell,
+    strat::DoubleQuadRule)
 
     # memory allocation here is a result from the type instability on strat
     # which is on purpose, i.e. the momintegrals! method is chosen based
diff --git a/test/test_assemblerow.jl b/test/test_assemblerow.jl
index 8e82630a..948e1b93 100644
--- a/test/test_assemblerow.jl
+++ b/test/test_assemblerow.jl
@@ -12,7 +12,7 @@ for T in [Float32, Float64]
 	numfunctions(X)
 
 	##
-	X1 = subset(X,1:1)
+	local X1 = subset(X,1:1)
 	numfunctions(X1)
 
 	T1 = assemble(t,X1,X)
diff --git a/test/test_basis.jl b/test/test_basis.jl
index 38e153b7..c6f36017 100644
--- a/test/test_basis.jl
+++ b/test/test_basis.jl
@@ -17,6 +17,8 @@ for T in [Float32, Float64]
     identityop    = Identity()
     doublelayer   = DoubleLayer(κ)
 
+    @show BEAST.defaultquadstrat(hypersingular, X, X)
+    # @show @which BEAST.defaultquadstrat(hypersingular, X, X)
     @time N = assemble(hypersingular, X, X)
     @time I = assemble(identityop, X, X)
 
diff --git a/test/test_wiltonints.jl b/test/test_wiltonints.jl
index aba1fd27..ed0acfbe 100644
--- a/test/test_wiltonints.jl
+++ b/test/test_wiltonints.jl
@@ -343,7 +343,9 @@ tqd = BE.quadpoints(x, [t], (12,))
 bqd = BE.quadpoints(x, [s], (13,))
 
 DQ_strategy = BE.DoubleQuadRule(tqd[1,1], bqd[1,1])
-BEAST.momintegrals!(op, x, x, t, s, z1, DQ_strategy)
+BEAST.momintegrals!(z1, op,
+    X, nothing, t,
+    X, nothing, s, DQ_strategy)
 
 SE_strategy = BE.WiltonSERule(
   tqd[1,1],
@@ -352,6 +354,8 @@ SE_strategy = BE.WiltonSERule(
 	bqd[1,1],
   ),
 )
-BEAST.momintegrals!(op, x, x, t, s, z2, SE_strategy)
+BEAST.momintegrals!(z2, op,
+    X, nothing, t,
+    X, nothing, s, SE_strategy)
 
 @test norm(z1-z2) < 1.0e-7