diff --git a/Project.toml b/Project.toml
index 67f335a0..fe353052 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "BEAST"
 uuid = "bb4162c7-ba94-5a20-af32-d8ec4428bdd1"
-version = "2.10.0"
+version = "2.11.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -67,7 +67,7 @@ ProgressMeter = "1.11.0"
 PlotlyBase = "0.8.21"
 Requires = "1"
 SauterSchwab3D = "0.2"
-SauterSchwabQuadrature = "2.4.0"
+SauterSchwabQuadrature = "2.5.0"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1, 2"
 StaticArrays = "0.8.3, 0.9, 0.10, 0.11, 0.12, 1"
 Suppressor = "0.2.8"
diff --git a/src/BEAST.jl b/src/BEAST.jl
index bde48f83..e7440983 100644
--- a/src/BEAST.jl
+++ b/src/BEAST.jl
@@ -230,6 +230,8 @@ include("quadrature/strategies/testrefinestrialqstrat.jl")
 include("quadrature/strategies/trialrefinestestqstrat.jl")
 include("quadrature/strategies/nonconftestbaryrefoftrialqstrat.jl")
 include("quadrature/strategies/timedomain/nothingqstrat.jl")
+include("quadrature/quadraturebuffer.jl")
+include("quadrature/quadrulecallbacks.jl")
 
 
 include("excitation.jl")
diff --git a/src/coloring.jl b/src/coloring.jl
index 2916c4f5..d019890c 100644
--- a/src/coloring.jl
+++ b/src/coloring.jl
@@ -41,7 +41,7 @@ function GraphsColoring.conflicts(
     end
 
     for i in eachindex(conflictindices)
-        conflictindices[i] = unique(conflictindices[i])
+        unique!(conflictindices[i])
     end
 
     return eachindex(elements),
diff --git a/src/decoupled/dpops.jl b/src/decoupled/dpops.jl
index 502b3e6c..bac57ed7 100644
--- a/src/decoupled/dpops.jl
+++ b/src/decoupled/dpops.jl
@@ -20,9 +20,11 @@ end
 function momintegrals!(out, op::CurlSingleLayerDP3D,
     test_local_space::RTRefSpace, tptr, test_triangular_element,
     trial_local_space::LagrangeRefSpace, bptr, trial_triangular_element,
-    qrule::SauterSchwabStrategy)
+    qrule::SauterSchwabStrategy, qbuffer)
 
-    I, J, K, L = SauterSchwabQuadrature.reorder(
+    sbuffer = sauterschwab_buffer(qbuffer, qrule)
+    I, J, K, L = sbuffer.I, sbuffer.J, sbuffer.K, sbuffer.L
+    sauterschwab_reorder!(I, J, K, L,
         test_triangular_element.vertices,
         trial_triangular_element.vertices, qrule)
 
diff --git a/src/integralop.jl b/src/integralop.jl
index 8e9fba70..a2a2e975 100644
--- a/src/integralop.jl
+++ b/src/integralop.jl
@@ -97,7 +97,7 @@ function assemblechunk!(biop::IntegralOperator, tfs::Space, bfs::Space, store;
     else
         quadstrat
     end
-    
+
     qd = quaddata(biop, tshapes, bshapes, test_elements, bsis_elements, qs)
     zlocal = zeros(scalartype(biop, tfs, bfs), 2num_tshapes, 2num_bshapes)
     # @show "after" qs
@@ -142,7 +142,6 @@ end
     @test BlockArrays.blocksizes(M) == [(n1,n1) (n1,n2); (n2,n1) (n2,n2)]
 end
 
-
 function assemblechunk_body!(biop, test_space, trial_space,
     test_elements, test_element_ptrs, test_assembly_data, active_test_els,
     trial_elements, trial_element_ptrs, trial_assembly_data, active_trial_els,
@@ -156,6 +155,7 @@ function assemblechunk_body!(biop, test_space, trial_space,
         @local begin
             zlocal = zeros(scalartype(biop, test_space, trial_space), num_tshapes, num_bshapes)
             tadjq = Vector{eltype(trial_assembly_data.data)}(undef, size(trial_assembly_data.data,1))
+            qbuffer = quadraturebuffer(quadstrat)
         end
         P = active_test_els[p]
         tcell = test_elements[P]
@@ -166,10 +166,11 @@ function assemblechunk_body!(biop, test_space, trial_space,
             bptr = trial_element_ptrs[Q]
 
             fill!(zlocal, 0)
-            qrule = quadrule(biop, refspace(test_space), refspace(trial_space),
+            apply = ApplyMomintegrals(zlocal, biop,
+                test_space, tptr, tcell, trial_space, bptr, bcell,
+                qbuffer)
+            quadrule(apply, biop, refspace(test_space), refspace(trial_space),
                 P, tcell, Q, bcell, qd, quadstrat)
-            momintegrals!(zlocal, biop,
-                test_space,  tptr, tcell, trial_space, bptr, bcell, qrule)
             for j in 1 : num_bshapes
                 tadjq .= @view trial_assembly_data.data[:,j,q]
                 for i in 1 : num_tshapes
@@ -184,6 +185,7 @@ function assemblechunk_body!(biop, test_space, trial_space,
 end end end end end end end
 
 
+
 # function assemblechunk_body_colored!(biop,
 #         test_space, testad, testelementids,
 #         trial_space, trialad, trialelementids,
@@ -202,7 +204,7 @@ end end end end end end end
 #         @set scheduler = scheduler
 #         @local zlocal = zeros(scalartype(biop, test_space, trial_space), num_tshapes, num_bshapes)
 #         tcell, tptr = test_elements[p], test_cell_ptrs[p]
-    
+
 #         for q in trialelementids
 #             bcell, bptr = trial_elements[q], trial_cell_ptrs[q]
 #             fill!(zlocal, 0)
@@ -432,7 +434,7 @@ end
 #                             m′ = get(test_id_in_blk, m, 0)
 #                             m′ == 0 && continue
 #                             store(a*zlocal[i,j]*b, m′, n′)
-#     end end end end end end 
+#     end end end end end end
 #     # put!(zlocals, zlocal)
 # end
 
@@ -629,6 +631,7 @@ function assemblerow_body!(biop,
     zlocal, quadrature_data, store; quadstrat)
 
     test_function = test_functions.fns[1]
+    qbuffer = quadraturebuffer(quadstrat)
     for shape in test_function
         p = shape.cellid
         i = shape.refid
@@ -637,11 +640,12 @@ function assemblerow_body!(biop,
         for (q,bcell) in enumerate(trial_elements)
 
             fill!(zlocal, 0)
-            qrule = quadrule(biop, test_shapes, trial_shapes, p, tcell, q, bcell, quadrature_data, quadstrat)
-            momintegrals!(zlocal, biop,
+            apply = ApplyMomintegrals(zlocal, biop,
                 test_functions, nothing, tcell,
                 trial_functions, nothing, bcell,
-                qrule)
+                qbuffer)
+            quadrule(apply, biop, test_shapes, trial_shapes,
+                p, tcell, q, bcell, quadrature_data, quadstrat)
 
             for j in 1:size(zlocal,2)
                 for (n,b) in trial_assembly_data[q,j]
@@ -682,6 +686,7 @@ function assemblecol_body!(biop,
     zlocal, quadrature_data, store; quadstrat)
 
     trial_function = trial_functions.fns[1]
+    qbuffer = quadraturebuffer(quadstrat)
     for shape in trial_function
         q = shape.cellid
         j = shape.refid
@@ -691,15 +696,14 @@ function assemblecol_body!(biop,
         for (p,tcell) in enumerate(test_elements)
 
             fill!(zlocal, 0)
-            qrule = quadrule(biop, test_shapes, trial_shapes, p, tcell, q, bcell, quadrature_data, quadstrat)
-            momintegrals!(zlocal, biop,
+            apply = ApplyMomintegrals(zlocal, biop,
                 test_functions, nothing, tcell,
-                trial_functions, nothing, bcell, qrule)
+                trial_functions, nothing, bcell,
+                qbuffer)
+            quadrule(apply, biop, test_shapes, trial_shapes,
+                p, tcell, q, bcell, quadrature_data, quadstrat)
 
             for i in 1:size(zlocal,1)
                 for (m,a) in test_assembly_data[p,i]
                     store(a*zlocal[i,j]*b, m, 1)
 end end end end end
-
-
-
diff --git a/src/maxwell/bogaertints.jl b/src/maxwell/bogaertints.jl
index 9307bb25..dcee5824 100644
--- a/src/maxwell/bogaertints.jl
+++ b/src/maxwell/bogaertints.jl
@@ -1,7 +1,7 @@
 function momintegrals!(z, op::MWSingleLayer3D,
     g::RTRefSpace, tptr, t,
     f::RTRefSpace, bptr, s,
-    strat::BogaertStrategy)
+    strat::BogaertStrategy, qbuffer)
 
     T, GG  = GetIntegrals(t, s, op.gamma, strat)
 
@@ -58,7 +58,7 @@ end
 function momintegrals!(z, op::MWDoubleLayer3D,
     g::RTRefSpace, tptr, τ,
     f::RTRefSpace, bptr, σ,
-    strat::BogaertStrategy)
+    strat::BogaertStrategy, qbuffer)
 
     # Get the primitives
     r = τ.vertices
diff --git a/src/maxwell/sauterschwabints_bdm_rt.jl b/src/maxwell/sauterschwabints_bdm_rt.jl
index b03d6ee8..8aba3321 100644
--- a/src/maxwell/sauterschwabints_bdm_rt.jl
+++ b/src/maxwell/sauterschwabints_bdm_rt.jl
@@ -60,9 +60,11 @@ end
 
 function momintegrals!(op::MWDoubleLayer3D,
     test_local_space::BDMRefSpace, trial_local_space::RTRefSpace,
-    test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)
+    test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy, qbuffer)
 
-    I, J, K, L = SauterSchwabQuadrature.reorder(
+    sbuffer = sauterschwab_buffer(qbuffer, strat)
+    I, J, K, L = sbuffer.I, sbuffer.J, sbuffer.K, sbuffer.L
+    sauterschwab_reorder!(I, J, K, L,
         test_triangular_element.vertices,
         trial_triangular_element.vertices, strat)
 
diff --git a/src/maxwell/wiltonints.jl b/src/maxwell/wiltonints.jl
index 9b715446..7f6bb59d 100644
--- a/src/maxwell/wiltonints.jl
+++ b/src/maxwell/wiltonints.jl
@@ -5,6 +5,8 @@ mutable struct WiltonSERule{P,Q} <: SingularityExtractionRule
     regularpart_quadrule::Q
 end
 
+quadraturebuffer(qrule::WiltonSERule) = quadraturebuffer(qrule.regularpart_quadrule)
+
 function innerintegrals!(op::MWSingleLayer3DSng, p, g, f, t, s, z,
         strat::WiltonSERule, dx)
 
diff --git a/src/operator.jl b/src/operator.jl
index f26dd0fe..657f4e56 100644
--- a/src/operator.jl
+++ b/src/operator.jl
@@ -246,12 +246,13 @@ function assemble!(operator::Operator, testfunctions::Space, trialfunctions::Spa
 
     testelements, testad, trialelements, trialad, qdata, _ = assembleblock_primer(
         operator, testfunctions, trialfunctions; quadstrat=quadstrat)
-    
+
     # testad = (testelements, testad, 1:length(testelements))
-    # trialad = (trialelements, trialad, 1:length(trialelements)) 
+    # trialad = (trialelements, trialad, 1:length(trialelements))
 
     testelementcolors = color(testfunctions, scheduler; addata=(testelements, testad, 1:length(testelements)))
     trialelementcolors = [1:length(trialelements)]
+    # trialelementcolors = color(trialfunctions, scheduler; addata=(trialelements, trialad, 1:length(trialelements)))
 
     qs = if CompScienceMeshes.refines(geometry(testfunctions), geometry(trialfunctions))
         TestRefinesTrialQStrat(quadstrat)
@@ -265,7 +266,7 @@ function assemble!(operator::Operator, testfunctions::Space, trialfunctions::Spa
 
     for (i, coloredtestelements) in enumerate(testelementcolors)
         testad1 = AssemblyData(testad.data[:,:,coloredtestelements])
-        for (j, coloredtrialelements) in enumerate(trialelementcolors)            
+        for (j, coloredtrialelements) in enumerate(trialelementcolors)
             # assemblechunk_body_colored!(operator,
             #     testfunctions, testad, coloredtestelements,
             #     trialfunctions, trialad, coloredtrialelements,
@@ -276,7 +277,7 @@ function assemble!(operator::Operator, testfunctions::Space, trialfunctions::Spa
                 trialelements, eachindex(trialelements), trialad, coloredtrialelements,
                 qdata, nothing, store; quadstrat=qs, scheduler)
             next!(pbar; step = length(testelementcolors[i]) * length(trialelementcolors[j]))
-    end end 
+    end end
     finish!(pbar)
 end
 
@@ -347,7 +348,7 @@ function assemble!(op::AbstractOperator, tfs::DirectProductSpace, bfs::DirectPro
     store, threading=Threading{:multi};
     quadstrat=defaultquadstrat(op, tfs[1], bfs[1]),
     kwargs...)
-    
+
     I = Int[0]
     for s in tfs.factors push!(I, last(I) + numfunctions(s)) end
     for (i,s) in enumerate(tfs.factors)
@@ -377,7 +378,7 @@ function assemble!(op::BlockDiagonalOperator, U::DirectProductSpace, V::DirectPr
     store, threading=Threading{:multi};
     quadstrat = defaultquadstrat(op, U, V),
     kwargs...)
-    
+
     @assert length(U.factors) == length(V.factors)
     I = Int[0]; for u in U.factors push!(I, last(I) + numfunctions(u)) end
     J = Int[0]; for v in V.factors push!(J, last(J) + numfunctions(v)) end
@@ -404,7 +405,7 @@ function assemble!(op::BlockFullOperators, U::DirectProductSpace, V::DirectProdu
     store, threading;
     quadstrat = defaultquadstrat(op, U, V),
     kwargs...)
-    
+
     I = Int[0]; for u in U.factors push!(I, last(I) + numfunctions(u)) end
     J = Int[0]; for v in V.factors push!(J, last(J) + numfunctions(v)) end
 
diff --git a/src/operators/quasilocalops.jl b/src/operators/quasilocalops.jl
index b4d1b6ea..5963f422 100644
--- a/src/operators/quasilocalops.jl
+++ b/src/operators/quasilocalops.jl
@@ -70,6 +70,7 @@ function assemblechunk!(op::QuasiLocalOperator, tfs::Space, bfs::Space, store321
 
     qd = quaddata(op, trefs, brefs, tels, bels, quadstrat)
     zlocal = zeros(T, num_trefs, num_brefs)
+    qbuffer = quadraturebuffer(quadstrat)
     tree = elementstree(bels, 1.1)
 
     δ = oprange(op)
@@ -93,10 +94,12 @@ function assemblechunk!(op::QuasiLocalOperator, tfs::Space, bfs::Space, store321
                 @assert q <= size(bad.data, 3)
 
                 fill!(zlocal, 0)
-                qrule = quadrule(op, trefs, brefs, p, tcell, q, bcell, qd, quadstrat)
-                momintegrals!(zlocal, op,
+                apply = ApplyMomintegrals(zlocal, op,
                     tfs, tptr, tcell,
-                    bfs, bptr, bcell, qrule)
+                    bfs, bptr, bcell,
+                    qbuffer)
+                quadrule(apply, op, trefs, brefs,
+                    p, tcell, q, bcell, qd, quadstrat)
 
                 for j in 1 : length(bad[q])
                     for i in 1 : length(tad[p])
diff --git a/src/quadrature/SauterSchwabQuadrature1D.jl b/src/quadrature/SauterSchwabQuadrature1D.jl
index 6193eae0..6d973d26 100644
--- a/src/quadrature/SauterSchwabQuadrature1D.jl
+++ b/src/quadrature/SauterSchwabQuadrature1D.jl
@@ -6,9 +6,9 @@ export SauterSchwabStrategy1D
 export CommonEdge, CommonVertex
 
 # functions
-export sauterschwab_parameterized1D, _MRWrules, reorder
+export sauterschwab_parameterized1D, _MRWrules, reorder, reorder!
 
 # -------- included files
 include("doublesauterschwabint.jl")
 
-end
\ No newline at end of file
+end
diff --git a/src/quadrature/commonfaceoverlappingedgeqstrat.jl b/src/quadrature/commonfaceoverlappingedgeqstrat.jl
index 183223ba..2cc06ad0 100644
--- a/src/quadrature/commonfaceoverlappingedgeqstrat.jl
+++ b/src/quadrature/commonfaceoverlappingedgeqstrat.jl
@@ -10,17 +10,23 @@ end
 function quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd,
     qs::CommonFaceOverlappingEdgeQStrat)
 
+    return quadrule(ReturnQuadrule(), a, 𝒳, 𝒴, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, a, 𝒳, 𝒴, i, τ, j, σ, qd,
+    qs::CommonFaceOverlappingEdgeQStrat)
+
     if CompScienceMeshes.overlap(τ, σ)
-        return quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)
+        return quadrule(f, a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)
     end
 
     for (i,λ) in pairs(faces(τ))
         for (j,μ) in pairs(faces(σ))
             if CompScienceMeshes.overlap(λ, μ)
-                return NonConformingTouchQRule(qs.conforming_qstrat, i, j)
+                return f(NonConformingTouchQRule(qs.conforming_qstrat, i, j))
     end end end
 
     # Either positive distance, common face, or common vertex, which can
     # be handled directly by the parent quadrature strategy
-    return quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)    
-end
\ No newline at end of file
+    return quadrule(f, a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)
+end
diff --git a/src/quadrature/doublenumints.jl b/src/quadrature/doublenumints.jl
index db383410..a8bd54b5 100644
--- a/src/quadrature/doublenumints.jl
+++ b/src/quadrature/doublenumints.jl
@@ -10,19 +10,19 @@ 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)
+    tshs, bshs, tcell, bcell, z, strat::DoubleQuadRule, qbuffer)
 
     igd = Integrand(biop, tshs, bshs, tcell, bcell)
 
     womps = strat.outer_quad_points
     wimps = strat.inner_quad_points
-    
+
     for womp in womps
         tgeo = womp.point
         tvals = womp.value
         M = length(tvals)
         jx = womp.weight
-        
+
         for wimp in wimps
             bgeo = wimp.point
             bvals = wimp.value
@@ -43,4 +43,4 @@ function momintegrals!(biop,
 end
 
 
-_TransposedStrat(a::DoubleQuadRule) = a
\ No newline at end of file
+_TransposedStrat(a::DoubleQuadRule) = a
diff --git a/src/quadrature/doublenumsauterqstrat.jl b/src/quadrature/doublenumsauterqstrat.jl
index af965e4b..1fc7392e 100644
--- a/src/quadrature/doublenumsauterqstrat.jl
+++ b/src/quadrature/doublenumsauterqstrat.jl
@@ -15,7 +15,7 @@ function quaddata(op::IntegralOperator,
 
     tqd = quadpoints(test_local_space,  test_charts,  (qs.outer_rule,))
     bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,))
-     
+
     leg = (
       convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_vert,0,1)),
       convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_edge,0,1)),
@@ -36,7 +36,7 @@ function quaddata(op::IntegralOperator,
     t_qp = quadpoints(test_local_space,  test_charts,  (qs.outer_rule,))
     b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,))
 
-    sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_common_vert,0,1), 
+    sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_common_vert,0,1),
                SauterSchwab3D._shunnham2D(qs.sauter_schwab_common_edge),
                SauterSchwab3D._shunnham3D(qs.sauter_schwab_common_face),
                SauterSchwab3D._shunnham4D(qs.sauter_schwab_common_tetr),)
@@ -53,7 +53,7 @@ function quaddata(op::IntegralOperator,
     t_qp = quadpoints(test_local_space,  test_charts,  (qs.outer_rule,))
     b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,))
 
-    sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_common_vert,0,1), 
+    sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_common_vert,0,1),
                SauterSchwab3D._shunnham2D(qs.sauter_schwab_common_edge),
                SauterSchwab3D._shunnham3D(qs.sauter_schwab_common_face),
                SauterSchwab3D._shunnham4D(qs.sauter_schwab_common_tetr),)
@@ -71,7 +71,7 @@ function quaddata(op::IntegralOperator,
     t_qp = quadpoints(test_local_space,  test_charts,  (qs.outer_rule,))
     b_qp = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,))
 
-    sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_common_vert,0,1), 
+    sing_qp = (SauterSchwab3D._legendre(qs.sauter_schwab_common_vert,0,1),
                SauterSchwab3D._shunnham2D(qs.sauter_schwab_common_edge),
                SauterSchwab3D._shunnham3D(qs.sauter_schwab_common_face),
                SauterSchwab3D._shunnham4D(qs.sauter_schwab_common_tetr),)
@@ -79,7 +79,15 @@ function quaddata(op::IntegralOperator,
 end
 
 
-function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace,
+function quadrule(op::IntegralOperator, g::RefSpace, h::RefSpace,
+    i, τ::CompScienceMeshes.Simplex{<:Any, 2},
+    j, σ::CompScienceMeshes.Simplex{<:Any, 2},
+    qd, qs::DoubleNumSauterQstrat)
+
+    return quadrule(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, op::IntegralOperator, g::RefSpace, h::RefSpace,
     i, τ::CompScienceMeshes.Simplex{<:Any, 2},
     j, σ::CompScienceMeshes.Simplex{<:Any, 2},
     qd, qs::DoubleNumSauterQstrat)
@@ -87,40 +95,77 @@ function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace,
     hits = _numhits(τ, σ)
     @assert hits <= 3
 
-    hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3])
-    hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
-    hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])
+    hits == 3 && return f(SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3]))
+    hits == 2 && return f(SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2]))
+    hits == 1 && return f(SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1]))
 
-    return DoubleQuadRule(
+    return f(DoubleQuadRule(
         qd.tpoints[1,i],
-        qd.bpoints[1,j],)
+        qd.bpoints[1,j],))
 end
 
 struct _TransposedStrat{A}
     strat::A
-end 
+end
+
+struct _TransposedQuadruleCallback{F <: QuadruleCallback} <: QuadruleCallback
+    f::F
+end
+
+function (f::_TransposedQuadruleCallback)(qrule)
+    return f.f(_TransposedStrat(qrule))
+end
+
+function quadrule(op::IntegralOperator, g::RefSpace, h::RefSpace,
+    i, τ::CompScienceMeshes.Simplex{<:Any, 3},
+    j, σ::CompScienceMeshes.Simplex{<:Any, 3},
+    qd, qs::DoubleNumSauterQstrat)
+    return quadrule(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
+end
+function quadrule(op::IntegralOperator, g::RefSpace, h::RefSpace,
+    i, τ::CompScienceMeshes.Simplex{<:Any, 3},
+    j, σ::CompScienceMeshes.Simplex{<:Any, 2},
+    qd, qs::DoubleNumSauterQstrat)
+    return quadrule(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(op::IntegralOperator, g::RefSpace, h::RefSpace,
+    i, τ::CompScienceMeshes.Simplex{<:Any, 2},
+    j, σ::CompScienceMeshes.Simplex{<:Any, 3},
+    qd, qs::DoubleNumSauterQstrat)
+    return quadrule(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
+end
 
-function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, 
+function quadrule(f::QuadruleCallback, op::IntegralOperator, g::RefSpace, h::RefSpace,
     i, τ::CompScienceMeshes.Simplex{<:Any, 3},
-    j, σ::CompScienceMeshes.Simplex{<:Any, 3}, 
-    qd, qs::DoubleNumSauterQstrat) 
-    qr_volume(op, g, f, i, τ, j, σ, qd, qs)
+    j, σ::CompScienceMeshes.Simplex{<:Any, 3},
+    qd, qs::DoubleNumSauterQstrat)
+
+    return qr_volume(f, op, g, h, i, τ, j, σ, qd, qs)
 end
-function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, 
+
+function quadrule(f::QuadruleCallback, op::IntegralOperator, g::RefSpace, h::RefSpace,
     i, τ::CompScienceMeshes.Simplex{<:Any, 3},
-    j, σ::CompScienceMeshes.Simplex{<:Any, 2}, 
-    qd, qs::DoubleNumSauterQstrat) 
-    qr_boundary(op, g, f, i, τ, j, σ, qd, qs)
+    j, σ::CompScienceMeshes.Simplex{<:Any, 2},
+    qd, qs::DoubleNumSauterQstrat)
+
+    return qr_boundary(f, op, g, h, i, τ, j, σ, qd, qs)
 end
 
-function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace, 
+function quadrule(f::QuadruleCallback, op::IntegralOperator, g::RefSpace, h::RefSpace,
     i, τ::CompScienceMeshes.Simplex{<:Any, 2},
-    j, σ::CompScienceMeshes.Simplex{<:Any, 3}, 
-    qd, qs::DoubleNumSauterQstrat) 
-    _TransposedStrat(qr_boundary(op, g, f, i, τ, j, σ, qd, qs))
+    j, σ::CompScienceMeshes.Simplex{<:Any, 3},
+    qd, qs::DoubleNumSauterQstrat)
+
+    return qr_boundary(_TransposedQuadruleCallback(f), op, g, h, i, τ, j, σ, qd, qs)
+end
+
+function qr_volume(op::IntegralOperator, g::RefSpace, h::RefSpace, i, τ, j, σ, qd, qs)
+
+    return qr_volume(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
 end
 
-function qr_volume(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ, j, σ, qd, qs)
+function qr_volume(f::QuadruleCallback, op::IntegralOperator, g::RefSpace, h::RefSpace, i, τ, j, σ, qd, qs)
 
     dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))
 
@@ -149,21 +194,26 @@ function qr_volume(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ, j, σ,
     #singData = SauterSchwab3D.Singularity{D,hits}(idx_t, idx_s )
    @assert hits <= 4
 
-    hits == 4 && return SauterSchwab3D.CommonVolume6D_S(SauterSchwab3D.Singularity6DVolume(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[4]))
-    hits == 3 && return SauterSchwab3D.CommonFace6D_S(SauterSchwab3D.Singularity6DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3]))
-    hits == 2 && return SauterSchwab3D.CommonEdge6D_S(SauterSchwab3D.Singularity6DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3],qd.sing_qp[4]))
-    hits == 1 && return SauterSchwab3D.CommonVertex6D_S(SauterSchwab3D.Singularity6DPoint(idx_t,idx_s),qd.sing_qp[3])
+    hits == 4 && return f(SauterSchwab3D.CommonVolume6D_S(SauterSchwab3D.Singularity6DVolume(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[4])))
+    hits == 3 && return f(SauterSchwab3D.CommonFace6D_S(SauterSchwab3D.Singularity6DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])))
+    hits == 2 && return f(SauterSchwab3D.CommonEdge6D_S(SauterSchwab3D.Singularity6DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3],qd.sing_qp[4])))
+    hits == 1 && return f(SauterSchwab3D.CommonVertex6D_S(SauterSchwab3D.Singularity6DPoint(idx_t,idx_s),qd.sing_qp[3]))
 
 
 
-    return DoubleQuadRule(
+    return f(DoubleQuadRule(
         qd[1][1,i],
-        qd[2][1,j])
+        qd[2][1,j]))
 
 end
 
 
-function qr_boundary(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ, j,  σ, qd, qs)
+function qr_boundary(op::IntegralOperator, g::RefSpace, h::RefSpace, i, τ, j, σ, qd, qs)
+
+    return qr_boundary(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
+end
+
+function qr_boundary(f::QuadruleCallback, op::IntegralOperator, g::RefSpace, h::RefSpace, i, τ, j,  σ, qd, qs)
 
     dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))
 
@@ -191,20 +241,28 @@ function qr_boundary(op::IntegralOperator, g::RefSpace, f::RefSpace, i, τ, j,
 
     @assert hits <= 3
     #singData = SauterSchwab3D.Singularity{D,hits}(idx_t, idx_s )
-   
 
-    hits == 3 && return SauterSchwab3D.CommonFace5D_S(SauterSchwab3D.Singularity5DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3]))
-    hits == 2 && return SauterSchwab3D.CommonEdge5D_S(SauterSchwab3D.Singularity5DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3]))
-    hits == 1 && return SauterSchwab3D.CommonVertex5D_S(SauterSchwab3D.Singularity5DPoint(idx_t,idx_s),(qd.sing_qp[3],qd.sing_qp[2]))
 
+    hits == 3 && return f(SauterSchwab3D.CommonFace5D_S(SauterSchwab3D.Singularity5DFace(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])))
+    hits == 2 && return f(SauterSchwab3D.CommonEdge5D_S(SauterSchwab3D.Singularity5DEdge(idx_t,idx_s),(qd.sing_qp[1],qd.sing_qp[2],qd.sing_qp[3])))
+    hits == 1 && return f(SauterSchwab3D.CommonVertex5D_S(SauterSchwab3D.Singularity5DPoint(idx_t,idx_s),(qd.sing_qp[3],qd.sing_qp[2])))
 
-    return DoubleQuadRule(
+
+    return f(DoubleQuadRule(
         qd[1][1,i],
-        qd[2][1,j])
+        qd[2][1,j]))
 
 end
 
-function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace,
+function quadrule(op::IntegralOperator, g::RefSpace, h::RefSpace,
+    i, τ::CompScienceMeshes.Quadrilateral,
+    j, σ::CompScienceMeshes.Quadrilateral,
+    qd, qs::DoubleNumSauterQstrat)
+
+    return quadrule(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, op::IntegralOperator, g::RefSpace, h::RefSpace,
     i, τ::CompScienceMeshes.Quadrilateral,
     j, σ::CompScienceMeshes.Quadrilateral,
     qd, qs::DoubleNumSauterQstrat)
@@ -213,13 +271,13 @@ function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace,
     @assert hits != 3
     @assert hits <= 4
 
-    hits == 4 && return SauterSchwabQuadrature.CommonFaceQuad(qd.gausslegendre[3])
-    hits == 2 && return SauterSchwabQuadrature.CommonEdgeQuad(qd.gausslegendre[2])
-    hits == 1 && return SauterSchwabQuadrature.CommonVertexQuad(qd.gausslegendre[1])
+    hits == 4 && return f(SauterSchwabQuadrature.CommonFaceQuad(qd.gausslegendre[3]))
+    hits == 2 && return f(SauterSchwabQuadrature.CommonEdgeQuad(qd.gausslegendre[2]))
+    hits == 1 && return f(SauterSchwabQuadrature.CommonVertexQuad(qd.gausslegendre[1]))
 
-    return DoubleQuadRule(
+    return f(DoubleQuadRule(
         qd.tpoints[1,i],
-        qd.bpoints[1,j],)
+        qd.bpoints[1,j],))
 end
 
 
@@ -237,4 +295,4 @@ function _numhits(τ, σ)
         end
     end
     return hits
-end
\ No newline at end of file
+end
diff --git a/src/quadrature/doublenumwiltonbogaertqstrat.jl b/src/quadrature/doublenumwiltonbogaertqstrat.jl
index f1aff891..e9eb0977 100644
--- a/src/quadrature/doublenumwiltonbogaertqstrat.jl
+++ b/src/quadrature/doublenumwiltonbogaertqstrat.jl
@@ -17,14 +17,20 @@ function quaddata(op::IntegralOperator,
     return (tpoints=tqd, bpoints=bqd)
 end
 
-function quadrule(op::IntegralOperator, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd,
+function quadrule(op::IntegralOperator, g::RTRefSpace, h::RTRefSpace, i, τ, j, σ, qd,
+    qs::DoubleNumWiltonBogaertQStrat)
+
+    return quadrule(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, op::IntegralOperator, g::RTRefSpace, h::RTRefSpace, i, τ, j, σ, qd,
     qs::DoubleNumWiltonBogaertQStrat)
 
     dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))
     xtol = 0.2
-  
+
     k = norm(gamma(op))
-  
+
     hits = 0
     xmin = xtol
     for t in τ.vertices
@@ -37,21 +43,21 @@ function quadrule(op::IntegralOperator, g::RTRefSpace, f::RTRefSpace, i, τ, j,
         end
       end
     end
-  
-    hits == 3   && return BogaertSelfPatchStrategy(5)
-    hits == 2   && return BogaertEdgePatchStrategy(8, 4)
-    hits == 1   && return BogaertPointPatchStrategy(2, 3)
+
+    hits == 3   && return f(BogaertSelfPatchStrategy(5))
+    hits == 2   && return f(BogaertEdgePatchStrategy(8, 4))
+    hits == 1   && return f(BogaertPointPatchStrategy(2, 3))
     rmin = xmin/k
-    xmin < xtol && return WiltonSERule(
+    xmin < xtol && return f(WiltonSERule(
       qd.tpoints[1,i],
       DoubleQuadRule(
         qd.tpoints[2,i],
         qd.bpoints[2,j],
       ),
-    )
-    return DoubleQuadRule(
+    ))
+    return f(DoubleQuadRule(
       qd.tpoints[1,i],
       qd.bpoints[1,j],
-    )
-  
-  end
\ No newline at end of file
+    ))
+
+  end
diff --git a/src/quadrature/doublenumwiltonsauterqstrat.jl b/src/quadrature/doublenumwiltonsauterqstrat.jl
index 6853ce55..8e0e6bb2 100644
--- a/src/quadrature/doublenumwiltonsauterqstrat.jl
+++ b/src/quadrature/doublenumwiltonsauterqstrat.jl
@@ -29,7 +29,17 @@ function quaddata(op::IntegralOperator,
 end
 
 
-function quadrule(op::IntegralOperator, g, f,  i, τ, j, σ, qd,
+function quadrule(op::IntegralOperator, g, h,
+    i, τ::CompScienceMeshes.Simplex{<:Any,2},
+    j, σ::CompScienceMeshes.Simplex{<:Any,2}, qd,
+    qs::DoubleNumWiltonSauterQStrat)
+
+    return quadrule(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, op::IntegralOperator, g, h,
+    i, τ::CompScienceMeshes.Simplex{<:Any,2},
+    j, σ::CompScienceMeshes.Simplex{<:Any,2}, qd,
     qs::DoubleNumWiltonSauterQStrat)
 
     T = eltype(eltype(τ.vertices))
@@ -48,22 +58,22 @@ function quadrule(op::IntegralOperator, g, f,  i, τ, j, σ, qd,
 
     @assert hits <= 3
 
-    hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3])
-    hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
-    hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])
+    hits == 3 && return f(SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3]))
+    hits == 2 && return f(SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2]))
+    hits == 1 && return f(SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1]))
 
     h2 = volume(σ)
     xtol2 = 0.2 * 0.2
     k2 = abs2(gamma(op))
     if max(dmin2*k2, dmin2/16h2) < xtol2
-        return WiltonSERule(
+        return f(WiltonSERule(
             qd.tpoints[2,i],
             DoubleQuadRule(
                 qd.tpoints[2,i],
-                qd.bpoints[2,j],),)
+                qd.bpoints[2,j],),))
     end
 
-    return DoubleQuadRule(
+    return f(DoubleQuadRule(
         qd.tpoints[1,i],
-        qd.bpoints[1,j],)
-end
\ No newline at end of file
+        qd.bpoints[1,j],))
+end
diff --git a/src/quadrature/doublesauterschwabint.jl b/src/quadrature/doublesauterschwabint.jl
index ae3a5cc6..a0b4d058 100644
--- a/src/quadrature/doublesauterschwabint.jl
+++ b/src/quadrature/doublesauterschwabint.jl
@@ -77,8 +77,18 @@ function _MRWrules(n,a,b)
     return collect(zip(x,w))
 end
 
-function sauterschwab_parameterized1D(integrand, strategy::SauterSchwabStrategy1D)
-    return sum(w1 * w2 * strategy(integrand, η, ξ) for (η, w1) in strategy.qpsi, (ξ, w2) in strategy.qpso)
+function sauterschwab_parameterized1D(integrand, strategy::S) where {S <: SauterSchwabStrategy1D}
+    η, w1 = strategy.qpsi[1]
+    ξ, w2 = strategy.qpso[1]
+    acc = zero(w1 * w2 * strategy(integrand, η, ξ))
+
+    for (η, w1) in strategy.qpsi
+        for (ξ, w2) in strategy.qpso
+            acc += w1 * w2 * strategy(integrand, η, ξ)
+        end
+    end
+
+    return acc
 end
 
 # In the reference domain [0, 1] x [0,1], we assume
@@ -88,7 +98,7 @@ end
 # ==> Therefore, the vertices must be mapped that the segments
 # with vertices [vt1, vt2] and [vs1, vs2] meet at vertices
 # vt2 and vs2 (since the barycentric coordinate ξ = 0 is at vt2/vs2)
-function reorder(t, s, strat::CommonVertex)
+function reorder!(I, J, K, L, t, s, strat::CommonVertex)
 
     T = eltype(t[1])
     tol = 1e3 * eps(T)
@@ -98,8 +108,6 @@ function reorder(t, s, strat::CommonVertex)
     # Find the permutation P of t and s that make
     # Pt = [P, A1, A2]
     # Ps = [P, B1, B2]
-    I = zeros(Int, 1)
-    J = zeros(Int, 1)
     e = 1
     for i in 1:2
         v = t[i]
@@ -115,10 +123,23 @@ function reorder(t, s, strat::CommonVertex)
         e == 2 && break
     end
 
-    prepend!(I, setdiff([1, 2], I))
-    prepend!(J, setdiff([1, 2], J))
+    i2 = I[1]
+    j2 = J[1]
+    for i in 1:2
+        if i != i2
+            I[1] = i
+            break
+        end
+    end
+    for j in 1:2
+        if j != j2
+            J[1] = j
+            break
+        end
+    end
+    I[2] = i2
+    J[2] = j2
 
-    K = zeros(Int, 2)
     for i in 1:2
         for j in 1:2
             if I[j] == i
@@ -128,7 +149,6 @@ function reorder(t, s, strat::CommonVertex)
         end
     end
 
-    L = zeros(Int, 2)
     for i in 1:2
         for j in 1:2
             if J[j] == i
@@ -141,25 +161,33 @@ function reorder(t, s, strat::CommonVertex)
     return I, J, K, L
 end
 
-function reorder(t, s, strat::CommonEdge)
+function reorder(t, s, strat::Union{CommonVertex,CommonEdge})
+    I = Vector{Int}(undef, 2)
+    J = Vector{Int}(undef, 2)
+    K = Vector{Int}(undef, 2)
+    L = Vector{Int}(undef, 2)
+    return reorder!(I, J, K, L, t, s, strat)
+end
+
+function reorder!(I, J, K, L, t, s, strat::CommonEdge)
 
     T = eltype(t[1])
     tol = 1e3 * eps(T)
     # tol = 1e5 * eps(T)
     # tol = sqrt(eps(T))
 
-    I = [1, 2]
-    J = zeros(Int, 2)
+    I[1] = 1
+    I[2] = 2
     v = t[1]
     w = s[1]
     if norm(w - v) < tol
-        J[:] = I[:]
+        J[1] = I[1]
+        J[2] = I[2]
     else # If first vertices do not coincide -> swap
         J[1] = 2
         J[2] = 1
     end
 
-    K = zeros(Int, 2)
     for i in 1:2
         for j in 1:2
             if I[j] == i
@@ -169,7 +197,6 @@ function reorder(t, s, strat::CommonEdge)
         end
     end
 
-    L = zeros(Int, 2)
     for i in 1:2
         for j in 1:2
             if J[j] == i
@@ -180,4 +207,4 @@ function reorder(t, s, strat::CommonEdge)
     end
 
     return I, J, K, L
-end
\ No newline at end of file
+end
diff --git a/src/quadrature/nonconformingintegralopqstrat.jl b/src/quadrature/nonconformingintegralopqstrat.jl
index abd2c394..55f21929 100644
--- a/src/quadrature/nonconformingintegralopqstrat.jl
+++ b/src/quadrature/nonconformingintegralopqstrat.jl
@@ -9,17 +9,23 @@ end
 function quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd,
     qs::NonConformingIntegralOpQStrat)
 
+    return quadrule(ReturnQuadrule(), a, 𝒳, 𝒴, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, a, 𝒳, 𝒴, i, τ, j, σ, qd,
+    qs::NonConformingIntegralOpQStrat)
+
     if CompScienceMeshes.overlap(τ, σ)
-        return NonConformingOverlapQRule(qs.conforming_qstrat)
+        return f(NonConformingOverlapQRule(qs.conforming_qstrat))
     end
 
     for (i,λ) in pairs(faces(τ))
         for (j,μ) in pairs(faces(σ))
             if CompScienceMeshes.overlap(λ, μ)
-                return NonConformingTouchQRule(qs.conforming_qstrat, i, j)
+                return f(NonConformingTouchQRule(qs.conforming_qstrat, i, j))
     end end end
 
     # Either positive distance or common vertex, both can
     # be handled directly by the parent quadrature strategy
-    return quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)    
-end
\ No newline at end of file
+    return quadrule(f, a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)
+end
diff --git a/src/quadrature/nonconformingoverlapqrule.jl b/src/quadrature/nonconformingoverlapqrule.jl
index 94186092..4a64fc04 100644
--- a/src/quadrature/nonconformingoverlapqrule.jl
+++ b/src/quadrature/nonconformingoverlapqrule.jl
@@ -1,6 +1,8 @@
 struct NonConformingOverlapQRule{S}
     conforming_qstrat::S
 end
+
+quadraturebuffer(qrule::NonConformingOverlapQRule) = quadraturebuffer(qrule.conforming_qstrat)
 # function tangent_rank(p::CompScienceMeshes.Simplex{U,D}) where {U,D}
 #     G = [dot(p.tangents[i], p.tangents[j]) for i in 1:D, j in 1:D]
 #     return rank(G) == D 
@@ -9,7 +11,7 @@ end
 function momintegrals!(op,
     test_local_space, basis_local_space,
     test_chart::CompScienceMeshes.Simplex, basis_chart::CompScienceMeshes.Simplex,
-    out, qrule::NonConformingOverlapQRule)
+    out, qrule::NonConformingOverlapQRule, qbuffer)
 
     num_tshapes = numfunctions(test_local_space, domain(test_chart))
     num_bshapes = numfunctions(basis_local_space, domain(basis_chart))
@@ -61,12 +63,12 @@ function momintegrals!(op,
         for (q,bchart) in enumerate(bsis_charts)
             restrict!(Q, basis_local_space, basis_chart, bchart, trial_overlaps[q])
 
-            qrule = quadrule(op, test_local_space, basis_local_space,
-                p, tchart, q, bchart, qdata, qstrat)
-
             fill!(zlocal, 0)
-            momintegrals!(op, test_local_space, basis_local_space,
-                tchart, bchart, zlocal, qrule)
+            apply = ApplyLocalMomintegrals(op,
+                test_local_space, basis_local_space,
+                tchart, bchart, zlocal, qbuffer)
+            quadrule(apply, op, test_local_space, basis_local_space,
+                p, tchart, q, bchart, qdata, qstrat)
 
             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 e2b63866..dfa6aa09 100644
--- a/src/quadrature/nonconformingtouchqrule.jl
+++ b/src/quadrature/nonconformingtouchqrule.jl
@@ -4,10 +4,12 @@ struct NonConformingTouchQRule{S}
     bsis_overlapping_edge_index::Int
 end
 
+quadraturebuffer(qrule::NonConformingTouchQRule) = quadraturebuffer(qrule.conforming_qstrat)
+
 function momintegrals!(op,
     test_locspace, bsis_locspace,
     τ::CompScienceMeshes.Simplex, σ::CompScienceMeshes.Simplex,
-    out, qrule::NonConformingTouchQRule)
+    out, qrule::NonConformingTouchQRule, qbuffer)
 
     num_tshapes = numfunctions(test_locspace, domain(τ))
     num_bshapes = numfunctions(bsis_locspace, domain(σ))
@@ -56,12 +58,12 @@ function momintegrals!(op,
         for (q,bchart) in enumerate(σs)
             restrict!(Q, bsis_locspace, σ, bchart, trial_overlaps[q])
             
-            qrule = quadrule(op, test_locspace, bsis_locspace,
-                p, tchart, q, bchart, qdata, qstrat)
-
             fill!(zlocal, 0)
-            momintegrals!(op, test_locspace, bsis_locspace,
-                tchart, bchart, zlocal, qrule)
+            apply = ApplyLocalMomintegrals(op,
+                test_locspace, bsis_locspace,
+                tchart, bchart, zlocal, qbuffer)
+            quadrule(apply, op, test_locspace, bsis_locspace,
+                p, tchart, q, bchart, qdata, qstrat)
 
             for i in axes(P,1)
                 for j in axes(Q,1)
diff --git a/src/quadrature/quadraturebuffer.jl b/src/quadrature/quadraturebuffer.jl
new file mode 100644
index 00000000..0f7d78ee
--- /dev/null
+++ b/src/quadrature/quadraturebuffer.jl
@@ -0,0 +1,36 @@
+quadraturebuffer(quadstrat) = (;)
+
+function quadraturebuffer(qs::Union{TestRefinesTrialQStrat,TrialRefinesTestQStrat,NonConformingIntegralOpQStrat,CommonFaceOverlappingEdgeQStrat})
+    return quadraturebuffer(qs.conforming_qstrat)
+end
+
+function quadraturebuffer(qs::NonConfTestBaryRefOfTrialQStrat)
+    return quadraturebuffer(qs.conforming_qstrat)
+end
+
+function _sauterschwab_buffer(n)
+    return (;
+        I = Vector{Int}(undef, n),
+        J = Vector{Int}(undef, n),
+        K = Vector{Int}(undef, n),
+        L = Vector{Int}(undef, n),
+    )
+end
+
+function quadraturebuffer(::Union{DoubleNumSauterQstrat,DoubleNumWiltonSauterQStrat,SelfSauterOtherwiseDNumQStrat,CommonFaceVertexSauterCommonEdgeWiltonPostitiveDistanceNumQStrat})
+    return (;
+        edge = _sauterschwab_buffer(2),
+        triangle = _sauterschwab_buffer(3),
+        quadrilateral = _sauterschwab_buffer(4),
+    )
+end
+
+
+
+function quadraturebuffer(::Union{SauterSchwabQuadrature.CommonVertex,SauterSchwabQuadrature.CommonEdge,SauterSchwabQuadrature.CommonFace})
+    return _sauterschwab_buffer(3)
+end
+
+function quadraturebuffer(::Union{SauterSchwabQuadrature.CommonVertexQuad,SauterSchwabQuadrature.CommonEdgeQuad,SauterSchwabQuadrature.CommonFaceQuad})
+    return _sauterschwab_buffer(4)
+end
diff --git a/src/quadrature/quadrulecallbacks.jl b/src/quadrature/quadrulecallbacks.jl
new file mode 100644
index 00000000..28b3531b
--- /dev/null
+++ b/src/quadrature/quadrulecallbacks.jl
@@ -0,0 +1,35 @@
+struct ApplyMomintegrals{Z,Op,TS,TP,TC,BS,BP,BC,Buf} <: QuadruleCallback
+    zlocal::Z
+    biop::Op
+    test_space::TS
+    tptr::TP
+    tcell::TC
+    trial_space::BS
+    bptr::BP
+    bcell::BC
+    qbuffer::Buf
+end
+
+function (f::ApplyMomintegrals)(qrule)
+    return momintegrals!(f.zlocal, f.biop,
+        f.test_space, f.tptr, f.tcell,
+        f.trial_space, f.bptr, f.bcell,
+        qrule, f.qbuffer)
+end
+
+struct ApplyLocalMomintegrals{Op,TS,BS,TC,BC,Z,Buf} <: QuadruleCallback
+    op::Op
+    test_local_space::TS
+    trial_local_space::BS
+    test_chart::TC
+    trial_chart::BC
+    out::Z
+    qbuffer::Buf
+end
+
+function (f::ApplyLocalMomintegrals)(qrule)
+    return momintegrals!(f.op,
+        f.test_local_space, f.trial_local_space,
+        f.test_chart, f.trial_chart,
+        f.out, qrule, f.qbuffer)
+end
diff --git a/src/quadrature/quadstrats.jl b/src/quadrature/quadstrats.jl
index ec9f2e58..a003098a 100644
--- a/src/quadrature/quadstrats.jl
+++ b/src/quadrature/quadstrats.jl
@@ -49,6 +49,17 @@ struct SingleNumQStrat{R} <: AbstractQuadStrat
     quad_rule::R
 end
 
+abstract type QuadruleCallback end
+
+struct ReturnQuadrule <: QuadruleCallback end
+
+(::ReturnQuadrule)(qrule) = qrule
+
+function quadrule(f::QuadruleCallback, args...)
+    return f(quadrule(args...))
+end
+
+
 function quadinfo(op, tfs, bfs; quadstrat=defaultquadstrat(op, tfs, bfs))
 
     tels, tad = assemblydata(tfs)
diff --git a/src/quadrature/rules/momintegrals.jl b/src/quadrature/rules/momintegrals.jl
index a599c539..b65b1adc 100644
--- a/src/quadrature/rules/momintegrals.jl
+++ b/src/quadrature/rules/momintegrals.jl
@@ -1,7 +1,7 @@
 function momintegrals!(out, op,
     test_functions, test_cellptr, test_chart,
     trial_functions, trial_cellptr, trial_chart,
-    quadrule)
+    quadrule, qbuffer)
 
     local_test_space = refspace(test_functions)
     local_trial_space = refspace(trial_functions)
@@ -9,5 +9,5 @@ function momintegrals!(out, op,
     momintegrals!(op,
         local_test_space, local_trial_space,
         test_chart, trial_chart,
-        out, quadrule)
-end
\ No newline at end of file
+        out, quadrule, qbuffer)
+end
diff --git a/src/quadrature/rules/testinbaryrefoftrialqrule.jl b/src/quadrature/rules/testinbaryrefoftrialqrule.jl
index 6f5239b7..6e487d53 100644
--- a/src/quadrature/rules/testinbaryrefoftrialqrule.jl
+++ b/src/quadrature/rules/testinbaryrefoftrialqrule.jl
@@ -2,10 +2,12 @@ struct TestInBaryRefOfTrialQRule{S}
     conforming_qstrat::S
 end
 
+quadraturebuffer(qrule::TestInBaryRefOfTrialQRule) = quadraturebuffer(qrule.conforming_qstrat)
+
 function BEAST.momintegrals!(out, op,
     test_functions, test_cell, test_chart,
     trial_functions, trial_cell, trial_chart,
-    qr::TestInBaryRefOfTrialQRule)
+    qr::TestInBaryRefOfTrialQRule, qbuffer)
 
     test_local_space = refspace(test_functions)
     trial_local_space = refspace(trial_functions)
@@ -59,15 +61,15 @@ function BEAST.momintegrals!(out, op,
     Q = zeros(T, num_tshapes, num_tshapes)
     out1 = zero(out)
     for (q,chart) in enumerate(trial_charts)
-        qr1 = BEAST.quadrule(op, test_local_space, trial_local_space,
-            1, test_chart, q ,chart, qd, quadstrat)
-            
         BEAST.restrict!(Q, trial_local_space, trial_chart, chart, X[q])
 
         fill!(out1, 0)
-        BEAST.momintegrals!(out1, op,
+        apply = BEAST.ApplyMomintegrals(out1, op,
             test_functions, nothing, test_chart,
-            trial_functions, nothing, chart, qr1)
+            trial_functions, nothing, chart,
+            qbuffer)
+        BEAST.quadrule(apply, op, test_local_space, trial_local_space,
+            1, test_chart, q, chart, qd, quadstrat)
 
         for j in 1:num_bshapes
             for i in 1:num_tshapes
diff --git a/src/quadrature/rules/testrefinestrialqrule.jl b/src/quadrature/rules/testrefinestrialqrule.jl
index baf388f7..7188d452 100644
--- a/src/quadrature/rules/testrefinestrialqrule.jl
+++ b/src/quadrature/rules/testrefinestrialqrule.jl
@@ -2,10 +2,12 @@ struct TestRefinesTrialQRule{S}
     conforming_qstrat::S
 end
 
+quadraturebuffer(qrule::TestRefinesTrialQRule) = quadraturebuffer(qrule.conforming_qstrat)
+
 function momintegrals!(out, op,
     test_functions::Space, test_cell, test_chart,
     trial_functions::Space, trial_cell, trial_chart,
-    qr::TestRefinesTrialQRule)
+    qr::TestRefinesTrialQRule, qbuffer)
 
     test_local_space = refspace(test_functions)
     trial_local_space = refspace(trial_functions)
@@ -38,13 +40,13 @@ function momintegrals!(out, op,
     for (q,chart) in enumerate(trial_charts)
         restrict!(Q, trial_local_space, trial_chart, chart, trial_overlaps[q])
 
-        qr = quadrule(op, test_local_space, trial_local_space,
-            1, test_chart, q ,chart, qd, quadstrat)
-
         fill!(zlocal, 0)
-        momintegrals!(zlocal, op,
+        apply = ApplyMomintegrals(zlocal, op,
             test_functions, nothing, test_chart,
-            trial_functions, nothing, chart, qr)
+            trial_functions, nothing, chart,
+            qbuffer)
+        quadrule(apply, op, test_local_space, trial_local_space,
+            1, test_chart, q, chart, qd, quadstrat)
 
         for j in 1:num_bshapes
             for i in 1:num_tshapes
diff --git a/src/quadrature/rules/trialrefinestestqrule.jl b/src/quadrature/rules/trialrefinestestqrule.jl
index f0c74f54..30cf41bc 100644
--- a/src/quadrature/rules/trialrefinestestqrule.jl
+++ b/src/quadrature/rules/trialrefinestestqrule.jl
@@ -2,10 +2,12 @@ struct TrialRefinesTestQRule{S}
     conforming_qstrat::S
 end
 
+quadraturebuffer(qrule::TrialRefinesTestQRule) = quadraturebuffer(qrule.conforming_qstrat)
+
 function momintegrals!(out, op,
     test_functions::Space, test_cell, test_chart,
     trial_functions::Space, trial_cell, trial_chart,
-    qr::TrialRefinesTestQRule)
+    qr::TrialRefinesTestQRule, qbuffer)
 
     test_local_space = refspace(test_functions)
     trial_local_space = refspace(trial_functions)
@@ -27,14 +29,14 @@ function momintegrals!(out, op,
         test_charts, [trial_chart], quadstrat)
 
     for (p,chart) in enumerate(test_charts)
-        qr = quadrule(op, test_local_space, trial_local_space,
-            p, chart, 1, trial_chart, qd, quadstrat)
-
         Q = restrict(test_local_space, test_chart, chart)
         zlocal = zero(out)
-        momintegrals!(zlocal, op,
+        apply = ApplyMomintegrals(zlocal, op,
             test_functions, nothing, chart,
-            trial_functions, trial_cell, trial_chart, qr)
+            trial_functions, trial_cell, trial_chart,
+            qbuffer)
+        quadrule(apply, op, test_local_space, trial_local_space,
+            p, chart, 1, trial_chart, qd, quadstrat)
 
         for j in 1:num_bshapes
             for i in 1:num_tshapes
diff --git a/src/quadrature/sauterschwabints.jl b/src/quadrature/sauterschwabints.jl
index f1d0ab7e..78fab1bf 100644
--- a/src/quadrature/sauterschwabints.jl
+++ b/src/quadrature/sauterschwabints.jl
@@ -8,10 +8,10 @@ end
 
 
 function (igd::Integrand)(u,v)
-    
+
     x = neighborhood(igd.test_chart,u)
     y = neighborhood(igd.trial_chart,v)
-    
+
     f = igd.local_test_space(x)
     g = igd.local_trial_space(y)
 
@@ -27,10 +27,10 @@ function (igd::Integrand{<:IntegralOperator,<:DivRefSpace,<:DivRefSpace})(u,v)
 
     x = CompScienceMeshes.neighborhood_lazy(igd.test_chart,u)
     y = CompScienceMeshes.neighborhood_lazy(igd.trial_chart,v)
-    
+
     p = neighborhood(test_domain, u)
     q = neighborhood(bsis_domain, v)
-    
+
     f̂ = igd.local_test_space(p)
     ĝ = igd.local_trial_space(q)
 
@@ -131,7 +131,11 @@ function pulledback_integrand(igd,
     ichart2 = CompScienceMeshes.permute_vertices(dom2, J)
 
     PulledBackIntegrand(igd, ichart1, ichart2)
-end 
+end
+
+function quadraturebuffer(::SauterSchwabQuadrature1D.SauterSchwabStrategy1D)
+    return _sauterschwab_buffer(2)
+end
 
 function sauterschwab_parameterized(igdp, rule::SauterSchwabStrategy)
     return SauterSchwabQuadrature.sauterschwab_parameterized(igdp, rule)
@@ -141,28 +145,32 @@ function sauterschwab_parameterized(igdp, rule::SauterSchwabQuadrature1D.SauterS
     return SauterSchwabQuadrature1D.sauterschwab_parameterized1D(igdp, rule)
 end
 
-function sauterschwab_reorder(test_vertices, trial_vertices, rule::SauterSchwabStrategy)
-    I, J, _, _ = SauterSchwabQuadrature.reorder(test_vertices, trial_vertices, rule)
+function sauterschwab_reorder!(I, J, K, L, test_vertices, trial_vertices, rule::SauterSchwabStrategy)
+    I, J, _, _ = SauterSchwabQuadrature.reorder!(I, J, K, L, test_vertices, trial_vertices, rule)
 
     return I, J
 end
 
-function sauterschwab_reorder(test_vertices, trial_vertices, rule::SauterSchwabQuadrature1D.SauterSchwabStrategy1D)
-    I, J, _, _ = SauterSchwabQuadrature1D.reorder(test_vertices, trial_vertices, rule)
+function sauterschwab_reorder!(I, J, K, L, test_vertices, trial_vertices, rule::SauterSchwabQuadrature1D.SauterSchwabStrategy1D)
+    I, J, _, _ = SauterSchwabQuadrature1D.reorder!(I, J, K, L, test_vertices, trial_vertices, rule)
 
     return I, J
 end
 
+
+sauterschwab_buffer(qbuffer, rule) = qbuffer
+sauterschwab_buffer(qbuffer, rule::SauterSchwabQuadrature1D.SauterSchwabStrategy1D) = hasproperty(qbuffer, :edge) ? qbuffer.edge : qbuffer
+sauterschwab_buffer(qbuffer, rule::Union{SauterSchwabQuadrature.CommonVertex,SauterSchwabQuadrature.CommonEdge,SauterSchwabQuadrature.CommonFace}) = hasproperty(qbuffer, :triangle) ? qbuffer.triangle : qbuffer
+sauterschwab_buffer(qbuffer, rule::Union{SauterSchwabQuadrature.CommonVertexQuad,SauterSchwabQuadrature.CommonEdgeQuad,SauterSchwabQuadrature.CommonFaceQuad}) = hasproperty(qbuffer, :quadrilateral) ? qbuffer.quadrilateral : qbuffer
+
 function momintegrals!(op::Operator,
     test_local_space, trial_local_space,
     test_chart, trial_chart,
-    out, rule::Union{SauterSchwabStrategy,SauterSchwabQuadrature1D.SauterSchwabStrategy1D})
+    out, rule::Union{SauterSchwabStrategy,SauterSchwabQuadrature1D.SauterSchwabStrategy1D}, qbuffer)
 
-    I, J = sauterschwab_reorder(
-        vertices(test_chart),
-        vertices(trial_chart),
-        rule
-    )
+    sbuffer = sauterschwab_buffer(qbuffer, rule)
+    I, J, K, L = sbuffer.I, sbuffer.J, sbuffer.K, sbuffer.L
+    sauterschwab_reorder!(I, J, K, L, vertices(test_chart), vertices(trial_chart), rule)
 
     num_tshapes = numfunctions(test_local_space, domain(test_chart))
     num_bshapes = numfunctions(trial_local_space, domain(trial_chart))
@@ -184,7 +192,7 @@ end
 function momintegrals!(op::Operator,
     test_local_space, trial_local_space,
     test_chart, trial_chart,
-    out, rule::SauterSchwab3DStrategy)
+    out, rule::SauterSchwab3DStrategy, qbuffer)
 
     I, J = SauterSchwab3D.reorder(rule.sing)
 
@@ -236,7 +244,7 @@ reversestrat(a::T) where {T <: SauterSchwab3D.SauterSchwab3DStrategy} = T(revers
 function momintegrals!(op::Operator,
     test_local_space, trial_local_space,
     test_chart, trial_chart,
-    out, rule::_TransposedStrat{<:SauterSchwab3DStrategy})
+    out, rule::_TransposedStrat{<:SauterSchwab3DStrategy}, qbuffer)
     rule2 = reversestrat(rule.strat)
     J, I = SauterSchwab3D.reorder(rule2.sing)
     #I2,J2 = SauterSchwab3D.reorder(rule.strat.sing)
@@ -249,4 +257,4 @@ function momintegrals!(op::Operator,
     G = SauterSchwab3D.sauterschwab_parameterized(igdp, rule2)
     out[1:num_tshapes, 1:num_bshapes] .+= G
     nothing
-end
\ No newline at end of file
+end
diff --git a/src/quadrature/selfsauterdnumotherwiseqstrat.jl b/src/quadrature/selfsauterdnumotherwiseqstrat.jl
index 1cbca3a4..752e1341 100644
--- a/src/quadrature/selfsauterdnumotherwiseqstrat.jl
+++ b/src/quadrature/selfsauterdnumotherwiseqstrat.jl
@@ -12,7 +12,7 @@ function quaddata(op::IntegralOperator,
 
     tqd = quadpoints(test_local_space,  test_charts,  (qs.outer_rule,))
     bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,))
-     
+
     leg = (
       convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common,0,1)),
       )
@@ -21,7 +21,13 @@ function quaddata(op::IntegralOperator,
 end
 
 
-function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace,  i, τ, j, σ, qd,
+function quadrule(op::IntegralOperator, g::RefSpace, h::RefSpace, i, τ, j, σ, qd,
+    qs::SelfSauterOtherwiseDNumQStrat)
+
+    return quadrule(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, op::IntegralOperator, g::RefSpace, h::RefSpace, i, τ, j, σ, qd,
     qs::SelfSauterOtherwiseDNumQStrat)
 
     T = eltype(eltype(τ.vertices))
@@ -40,9 +46,9 @@ function quadrule(op::IntegralOperator, g::RefSpace, f::RefSpace,  i, τ, j, σ,
 
     @assert hits <= 3
 
-    hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[1])
+    hits == 3 && return f(SauterSchwabQuadrature.CommonFace(qd.gausslegendre[1]))
 
-    return DoubleQuadRule(
+    return f(DoubleQuadRule(
         qd.tpoints[1,i],
-        qd.bpoints[1,j],)
-end
\ No newline at end of file
+        qd.bpoints[1,j],))
+end
diff --git a/src/quadrature/singularityextractionints.jl b/src/quadrature/singularityextractionints.jl
index 6cb3377b..896faa06 100644
--- a/src/quadrature/singularityextractionints.jl
+++ b/src/quadrature/singularityextractionints.jl
@@ -3,7 +3,7 @@ regularpart_quadrule(qr::SingularityExtractionRule) = qr.regularpart_quadrule
 
 function momintegrals!(op,
     g, f,t, s,
-    z, qrule::SingularityExtractionRule)
+    z, qrule::SingularityExtractionRule, qbuffer)
 
     womps = qrule.outer_quad_points
 
@@ -11,7 +11,7 @@ function momintegrals!(op,
     rop = regularpart(op)
 
     regqrule = regularpart_quadrule(qrule)
-    momintegrals!(rop, g, f, t, s, z, regqrule)
+    momintegrals!(rop, g, f, t, s, z, regqrule, qbuffer)
 
     for p in 1 : length(womps)
         x = womps[p].point
diff --git a/src/quadrature/strategies/cfcvsautercewiltonpdnumqstrat.jl b/src/quadrature/strategies/cfcvsautercewiltonpdnumqstrat.jl
index 5e86b67c..686028b7 100644
--- a/src/quadrature/strategies/cfcvsautercewiltonpdnumqstrat.jl
+++ b/src/quadrature/strategies/cfcvsautercewiltonpdnumqstrat.jl
@@ -24,7 +24,13 @@ function quaddata(op::IntegralOperator, test_local_space, bsis_local_space,
             qs.sauter_schwab_common_vert,))
 end
 
-function quadrule(op::IntegralOperator, g, f,  i, τ, j, σ,
+function quadrule(op::IntegralOperator, g, h, i, τ, j, σ,
+    qd, qs::CommonFaceVertexSauterCommonEdgeWiltonPostitiveDistanceNumQStrat)
+
+    return quadrule(ReturnQuadrule(), op, g, h, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, op::IntegralOperator, g, h, i, τ, j, σ,
     qd, qs::CommonFaceVertexSauterCommonEdgeWiltonPostitiveDistanceNumQStrat)
 
     T = eltype(eltype(τ.vertices))
@@ -43,27 +49,27 @@ function quadrule(op::IntegralOperator, g, f,  i, τ, j, σ,
 
     @assert hits <= 3
 
-    hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3])
-    # hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
+    hits == 3 && return f(SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3]))
+    # hits == 2 && return f(SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2]))
     if hits == 2
-        return WiltonSERule(
+        return f(WiltonSERule(
             qd.tpoints[2,i],
-            SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2]),)
+            SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2]),))
     end
-    hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])
+    hits == 1 && return f(SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1]))
 
     h2 = volume(σ)
     xtol2 = 0.2 * 0.2
     k2 = abs2(gamma(op))
     if max(dmin2*k2, dmin2/16h2) < xtol2
-        return WiltonSERule(
+        return f(WiltonSERule(
             qd.tpoints[2,i],
             DoubleQuadRule(
                 qd.tpoints[2,i],
-                qd.bpoints[2,j],),)
+                qd.bpoints[2,j],),))
     end
 
-    return DoubleQuadRule(
+    return f(DoubleQuadRule(
         qd.tpoints[1,i],
-        qd.bpoints[1,j],)
-end
\ No newline at end of file
+        qd.bpoints[1,j],))
+end
diff --git a/src/quadrature/strategies/nonconftestbaryrefoftrialqstrat.jl b/src/quadrature/strategies/nonconftestbaryrefoftrialqstrat.jl
index 51c1394d..70e84ac1 100644
--- a/src/quadrature/strategies/nonconftestbaryrefoftrialqstrat.jl
+++ b/src/quadrature/strategies/nonconftestbaryrefoftrialqstrat.jl
@@ -9,11 +9,7 @@ end
 function BEAST.quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd,
     quadstrat::NonConfTestBaryRefOfTrialQStrat)
 
-    # return TestInBaryRefOfTrialQRule(quadstrat.conforming_qstrat)
-    nh = BEAST._numhits(τ, σ)
-    nh > 0 && return TestInBaryRefOfTrialQRule(quadstrat.conforming_qstrat)
-    return BEAST.quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd,
-        quadstrat.conforming_qstrat)
+    return BEAST.quadrule(BEAST.ReturnQuadrule(), a, 𝒳, 𝒴, i, τ, j, σ, qd, quadstrat)
 end
 
 @testitem "NonConfTestBaryRefOfTrialQStrat" begin
@@ -45,4 +41,14 @@ end
     @test norm(Kyx2 - Kyx3) < 0.002
     @test norm(Kyx2 - Kyx4) < 0.002
     @test norm(Kyx3 - Kyx4) < 1e-12
-end
\ No newline at end of file
+end
+
+function BEAST.quadrule(f::BEAST.QuadruleCallback, a, 𝒳, 𝒴, i, τ, j, σ, qd,
+    quadstrat::NonConfTestBaryRefOfTrialQStrat)
+
+    # return f(TestInBaryRefOfTrialQRule(quadstrat.conforming_qstrat))
+    nh = BEAST._numhits(τ, σ)
+    nh > 0 && return f(TestInBaryRefOfTrialQRule(quadstrat.conforming_qstrat))
+    return BEAST.quadrule(f, a, 𝒳, 𝒴, i, τ, j, σ, qd,
+        quadstrat.conforming_qstrat)
+end
diff --git a/src/quadrature/strategies/testrefinestrialqstrat.jl b/src/quadrature/strategies/testrefinestrialqstrat.jl
index e295c6e5..3a9e54f5 100644
--- a/src/quadrature/strategies/testrefinestrialqstrat.jl
+++ b/src/quadrature/strategies/testrefinestrialqstrat.jl
@@ -9,10 +9,16 @@ end
 function quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd,
     qs::TestRefinesTrialQStrat)
 
+    return quadrule(ReturnQuadrule(), a, 𝒳, 𝒴, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, a, 𝒳, 𝒴, i, τ, j, σ, qd,
+    qs::TestRefinesTrialQStrat)
+
     hits = _numhits(τ, σ)
     if hits > 0
-        return TestRefinesTrialQRule(qs.conforming_qstrat)
+        return f(TestRefinesTrialQRule(qs.conforming_qstrat))
     end
 
-    return quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)    
-end
\ No newline at end of file
+    return quadrule(f, a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)
+end
diff --git a/src/quadrature/strategies/trialrefinestestqstrat.jl b/src/quadrature/strategies/trialrefinestestqstrat.jl
index 7cd70e03..d9571df3 100644
--- a/src/quadrature/strategies/trialrefinestestqstrat.jl
+++ b/src/quadrature/strategies/trialrefinestestqstrat.jl
@@ -9,10 +9,16 @@ end
 function quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd,
     qs::TrialRefinesTestQStrat)
 
+    return quadrule(ReturnQuadrule(), a, 𝒳, 𝒴, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, a, 𝒳, 𝒴, i, τ, j, σ, qd,
+    qs::TrialRefinesTestQStrat)
+
     hits = _numhits(τ, σ)
     if hits > 0
-        return TrialRefinesTestQRule(qs.conforming_qstrat)
+        return f(TrialRefinesTestQRule(qs.conforming_qstrat))
     end
 
-    return quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)    
-end
\ No newline at end of file
+    return quadrule(f, a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)
+end
diff --git a/src/volumeintegral/sauterschwab_ints.jl b/src/volumeintegral/sauterschwab_ints.jl
index bb0050e4..843c6c58 100644
--- a/src/volumeintegral/sauterschwab_ints.jl
+++ b/src/volumeintegral/sauterschwab_ints.jl
@@ -22,8 +22,14 @@ function reorder(test_chart::CompScienceMeshes.Simplex{<:Any,2,<:Any,3}, trial_c
     return I, J
 end
 
-function reorder(test_chart::CompScienceMeshes.Simplex, trial_chart::CompScienceMeshes.Simplex, 
-    strat::SauterSchwab3DStrategy)
+function reorder(test_chart::CompScienceMeshes.Simplex{<:Any,2,<:Any,3}, trial_chart::CompScienceMeshes.Simplex{<:Any,3,<:Any,4},
+    strat::SauterSchwab3DStrategy, qbuffer)
+    J, I = SauterSchwab3D.reorder(strat.sing) # 5D integral: ∫∫_Γ ∫∫∫_Ω
+    return I, J
+end
+
+function reorder(test_chart::CompScienceMeshes.Simplex, trial_chart::CompScienceMeshes.Simplex,
+    strat::SauterSchwab3DStrategy, qbuffer)
     I, J = SauterSchwab3D.reorder(strat.sing) # 6D integral: ∫∫∫_Ω ∫∫∫_Ω, 5D integral: ∫∫∫_Ω ∫∫_Γ
     return I, J
 end
@@ -35,7 +41,7 @@ end
 function momintegrals!(out, op::VIEOperator,
     test_functions::Space, test_ptr, test_chart,
     trial_functions::Space, trial_ptr, trial_chart,
-    strat::SauterSchwab3DStrategy)
+    strat::SauterSchwab3DStrategy, qbuffer)
 
     test_local_space = refspace(test_functions)
     trial_local_space = refspace(trial_functions)
@@ -43,7 +49,7 @@ function momintegrals!(out, op::VIEOperator,
     num_tshapes = numfunctions(test_local_space, domain(test_chart))
     num_bshapes = numfunctions(trial_local_space, domain(trial_chart))
 
-    I, J = reorder(test_chart, trial_chart, strat)
+    I, J = reorder(test_chart, trial_chart, strat, qbuffer)
         
     igd = Integrand(op, test_local_space, trial_local_space, test_chart, trial_chart)
     igdp = pulledback_integrand(igd, I, test_chart, J, trial_chart)
diff --git a/test/test_bogaertints.jl b/test/test_bogaertints.jl
index e364025a..84cc03e2 100644
--- a/test/test_bogaertints.jl
+++ b/test/test_bogaertints.jl
@@ -30,10 +30,10 @@ SE_strategy = BE.WiltonSERule(
 	bqd[1,1],
   ),
 )
-BEAST.momintegrals!(op, x, x, t, t, z1, SE_strategy)
+BEAST.momintegrals!(op, x, x, t, t, z1, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
 
 EE_strategy = BEAST.BogaertSelfPatchStrategy(20)
-BEAST.momintegrals!(op, x, x, t, t, z2, EE_strategy)
+BEAST.momintegrals!(op, x, x, t, t, z2, EE_strategy, BEAST.quadraturebuffer(EE_strategy))
 
 Γ = meshrectangle(1.0, 1.0, 1.0)
 t = chart(Γ, Γ.faces[1])
@@ -42,7 +42,7 @@ s = chart(Γ, Γ.faces[2])
 
 z3 = zeros(ComplexF64, n, n)
 EE_strategy = BEAST.BogaertEdgePatchStrategy(13,30)
-BEAST.momintegrals!(op, x, x, t, s, z3, EE_strategy)
+BEAST.momintegrals!(op, x, x, t, s, z3, EE_strategy, BEAST.quadraturebuffer(EE_strategy))
 
 z4 = zeros(ComplexF64, n, n)
 tqd = BE.quadpoints(x, [t], (12,13))
@@ -54,7 +54,7 @@ SE_strategy = BE.WiltonSERule(
 	bqd[1,1],
   ),
 )
-BEAST.momintegrals!(op, x, x, t, s, z4, SE_strategy)
+BEAST.momintegrals!(op, x, x, t, s, z4, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
 
 p = [
     point(0.0, 0.0, 0.0),
@@ -68,7 +68,7 @@ s = simplex(p[1], p[4], p[5])
 
 z5 = zeros(ComplexF64, n, n)
 EE_strategy = BEAST.BogaertPointPatchStrategy(9,10)
-BEAST.momintegrals!(op, x, x, t, s, z5, EE_strategy)
+BEAST.momintegrals!(op, x, x, t, s, z5, EE_strategy, BEAST.quadraturebuffer(EE_strategy))
 
 z6 = zeros(ComplexF64, n, n)
 tqd = BE.quadpoints(x, [t], (12,13))
@@ -80,7 +80,7 @@ SE_strategy = BE.WiltonSERule(
 	bqd[1,1],
   ),
 )
-BEAST.momintegrals!(op, x, x, t, s, z6, SE_strategy)
+BEAST.momintegrals!(op, x, x, t, s, z6, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
 
 @test norm(z1-z2)/norm(z1) < 1.0e-6
 @test norm(z3-z4)/norm(z3) < 3.0e-6
@@ -127,9 +127,9 @@ s3 = BE.DoubleQuadRule(
   tqd[1,1],
 )
 
-BEAST.momintegrals!(op, x, x, t1, t2, z1, s1)
-BEAST.momintegrals!(op, x, x, t1, t2, z2, s2)
-BEAST.momintegrals!(op, x, x, t1, t2, z3, s3)
+BEAST.momintegrals!(op, x, x, t1, t2, z1, s1, BEAST.quadraturebuffer(s1))
+BEAST.momintegrals!(op, x, x, t1, t2, z2, s2, BEAST.quadraturebuffer(s2))
+BEAST.momintegrals!(op, x, x, t1, t2, z3, s3, BEAST.quadraturebuffer(s3))
 
 @test norm(z2-z1)/norm(z2) < 1.0e-3
 
@@ -148,9 +148,9 @@ s3 = BE.DoubleQuadRule(
   tqd[1,1],
 )
 
-fill!(z1, 0); BEAST.momintegrals!(op, x, x, t2, t3, z1, s1)
-fill!(z2, 0); BEAST.momintegrals!(op, x, x, t2, t3, z2, s2)
-fill!(z3, 0); BEAST.momintegrals!(op, x, x, t2, t3, z3, s3)
+fill!(z1, 0); BEAST.momintegrals!(op, x, x, t2, t3, z1, s1, BEAST.quadraturebuffer(s1))
+fill!(z2, 0); BEAST.momintegrals!(op, x, x, t2, t3, z2, s2, BEAST.quadraturebuffer(s2))
+fill!(z3, 0); BEAST.momintegrals!(op, x, x, t2, t3, z3, s3, BEAST.quadraturebuffer(s3))
 
 @test norm(z2-z1)/norm(z2) < 2.0e-6
 
@@ -169,9 +169,9 @@ s3 = BE.DoubleQuadRule(
   tqd[1,1],
 )
 
-fill!(z1, 0); BEAST.momintegrals!(op, x, x, t1, t1, z1, s1)
-fill!(z2, 0); BEAST.momintegrals!(op, x, x, t1, t1, z2, s2)
-fill!(z3, 0); BEAST.momintegrals!(op, x, x, t1, t1, z3, s3)
+fill!(z1, 0); BEAST.momintegrals!(op, x, x, t1, t1, z1, s1, BEAST.quadraturebuffer(s1))
+fill!(z2, 0); BEAST.momintegrals!(op, x, x, t1, t1, z2, s2, BEAST.quadraturebuffer(s2))
+fill!(z3, 0); BEAST.momintegrals!(op, x, x, t1, t1, z3, s3, BEAST.quadraturebuffer(s3))
 
 
 
diff --git a/test/test_hh3dints.jl b/test/test_hh3dints.jl
index 008ad7c8..95fec37a 100644
--- a/test/test_hh3dints.jl
+++ b/test/test_hh3dints.jl
@@ -47,9 +47,9 @@ t1 = simplex(
   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)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t1, z_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t1, z_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t1, z_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_se ≈ z_ss rtol=1e-5
   @test z_double ≈ z_ss rtol = 1e-1
@@ -70,9 +70,9 @@ t1 = simplex(
   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)
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t1, z_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t1, z_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t1, z_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ss ≈ z_se rtol = 1e-5
   @test z_ss ≈ z_double rtol = 1e-1
@@ -93,9 +93,9 @@ t1 = simplex(
   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)
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t1, z_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t1, z_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t1, z_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ss ≈ z_se rtol = 1e-5
   @test z_ss ≈ z_double rtol = 1e-1
@@ -116,9 +116,9 @@ t1 = simplex(
   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)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t1, z_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t1, z_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t1, z_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ss ≈ z_se rtol = 1e-5
   @test z_ss ≈ z_double rtol = 1e-1
@@ -168,9 +168,9 @@ end
   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)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_ss rtol=1e-4
@@ -180,9 +180,9 @@ end
   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)
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_double rtol=1e-4
@@ -193,9 +193,9 @@ end
   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)
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_ss rtol=1e-4
@@ -216,9 +216,9 @@ end
   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)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_double rtol=1e-4
@@ -229,9 +229,9 @@ end
   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)
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_ss rtol=1e-4
@@ -242,9 +242,9 @@ end
   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)
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-3
   @test z_cv_double ≈ z_cv_se rtol = 1e-5
@@ -255,9 +255,9 @@ end
   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)
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   # @show z_cv_double
   # @show z_cv_se
@@ -285,9 +285,9 @@ end
   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)
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_ss rtol=1e-4
@@ -298,9 +298,9 @@ end
   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)
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_ss rtol=1e-4
@@ -311,9 +311,9 @@ end
   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)
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_ss rtol=1e-4
@@ -337,9 +337,9 @@ end
   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)
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_ss rtol=1e-4
@@ -350,9 +350,9 @@ end
   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)
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_ss rtol=1e-4
@@ -363,9 +363,9 @@ end
   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)
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_cv_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_cv_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_cv_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_cv_double ≈ z_cv_ss rtol = 1e-4
   @test z_cv_se ≈ z_cv_ss rtol=1e-4
@@ -414,9 +414,9 @@ end
   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)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-3
   @test z_ce_se ≈ z_ce_ss rtol=1e-4
@@ -427,9 +427,9 @@ end
   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)
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op2, lag0, lag0, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_ss ≈ z_ce_double rtol = 1e-1
   @test z_ce_se ≈ z_ce_ss rtol = 1e-2
@@ -440,9 +440,9 @@ end
   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)
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op3, lag0, lag0, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-1
   @test z_ce_se ≈ z_ce_ss rtol=1e-2
@@ -463,9 +463,9 @@ end
   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)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-3
   @test z_ce_se ≈ z_ce_double rtol=1e-3
@@ -476,9 +476,9 @@ end
   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)
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-1
   @test z_ce_se ≈ z_ce_ss rtol=1e-4
@@ -489,9 +489,9 @@ end
   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)
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-1
   @test z_ce_ss ≈ z_ce_se rtol = 1e-2
@@ -502,9 +502,9 @@ end
   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)
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-3
   @test z_ce_se ≈ z_ce_ss rtol=1e-4
@@ -528,9 +528,9 @@ end
   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)
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag1, lag0, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-3
   @test z_ce_se ≈ z_ce_ss rtol=1e-4
@@ -541,9 +541,9 @@ end
   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)
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op2, lag1, lag0, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-1
   @test z_ce_se ≈ z_ce_ss rtol=1e-4
@@ -554,9 +554,9 @@ end
   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)
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-1
   @test z_ce_se ≈ z_ce_ss rtol=1e-2
@@ -580,9 +580,9 @@ end
   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)
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op1, lag0, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-3
   @test z_ce_se ≈ z_ce_ss rtol=1e-4
@@ -593,9 +593,9 @@ end
   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)
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-1
   @test z_ce_se ≈ z_ce_ss rtol=1e-4
@@ -606,9 +606,9 @@ end
   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)
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_ce_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
+  BEAST.momintegrals!(op3, lag0, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_double ≈ z_ce_ss rtol = 1e-1
   @test z_ce_se ≈ z_ce_ss rtol=1e-1
@@ -648,8 +648,8 @@ end
   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)
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag0, lag0, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
   
   @test z_ce_se ≈ z_ce_double rtol=1e-6
 
@@ -684,8 +684,8 @@ end
   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)
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op1, lag1, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_se ≈ z_ce_double rtol=1e-4
 
@@ -701,8 +701,8 @@ end
     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)
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op2, lag0, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
   
   @test z_ce_se ≈ z_ce_double rtol = 1e-4
 
@@ -717,8 +717,8 @@ end
     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)
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op2, lag1, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_se ≈ z_ce_double rtol = 1e-4
   # doublelayer transposed
@@ -733,8 +733,8 @@ end
   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)
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op3, lag1, lag0, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_se ≈ z_ce_double rtol=1e-4
 
@@ -749,8 +749,8 @@ end
   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)
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op3, lag1, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
 
   @test z_ce_se ≈ z_ce_double rtol=1e-4
   # hypersingular
@@ -766,8 +766,8 @@ end
   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)
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_ce_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+  BEAST.momintegrals!(op4, lag1, lag1, t1, t2, z_ce_double, Double_strategy, BEAST.quadraturebuffer(Double_strategy))
   @test z_ce_se ≈ z_ce_double rtol = 1e-5
 
 end
diff --git a/test/test_sauterschwabints.jl b/test/test_sauterschwabints.jl
index ef5d190a..85f0da0a 100644
--- a/test/test_sauterschwabints.jl
+++ b/test/test_sauterschwabints.jl
@@ -45,16 +45,16 @@ SS_strategy = SauterSchwabQuadrature.CommonFace(BEAST._legendre(8,T(0.0),T(1.0))
 z_se = zeros(3,3)
 z_ss = zeros(3,3)
 
-BEAST.momintegrals!(op1, rt, rt, t1, t1, z_se, SE_strategy)
-BEAST.momintegrals!(op1, rt, rt, t1, t1, z_ss, SS_strategy)
+BEAST.momintegrals!(op1, rt, rt, t1, t1, z_se, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+BEAST.momintegrals!(op1, rt, rt, t1, t1, z_ss, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
 
 @test z_se ≈ z_ss atol=1e-4
 
 z_se2 = zeros(3,3)
 z_ss2 = zeros(3,3)
 
-BEAST.momintegrals!(op2, rt, rt, t1, t1, z_se2, SE_strategy)
-BEAST.momintegrals!(op2, rt, rt, t1, t1, z_ss2, SS_strategy)
+BEAST.momintegrals!(op2, rt, rt, t1, t1, z_se2, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+BEAST.momintegrals!(op2, rt, rt, t1, t1, z_ss2, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
 
 @test z_se2 ≈ z_ss2 atol=1e-4
 
@@ -73,18 +73,18 @@ z_cv_ss_2 = zeros(3,3)
 z_cv_se_3 = zeros(3,3)
 z_cv_ss_3 = zeros(3,3)
 
-BEAST.momintegrals!(op1, rt, rt, t1, t2, z_cv_se_1, SE_strategy)
-BEAST.momintegrals!(op1, rt, rt, t1, t2, z_cv_ss_1, SS_strategy)
+BEAST.momintegrals!(op1, rt, rt, t1, t2, z_cv_se_1, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+BEAST.momintegrals!(op1, rt, rt, t1, t2, z_cv_ss_1, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
 
 @test z_cv_se_1 ≈ z_cv_ss_1 atol=1e-7
 
-BEAST.momintegrals!(op2, rt, rt, t1, t2, z_cv_se_2, SE_strategy)
-BEAST.momintegrals!(op2, rt, rt, t1, t2, z_cv_ss_2, SS_strategy)
+BEAST.momintegrals!(op2, rt, rt, t1, t2, z_cv_se_2, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+BEAST.momintegrals!(op2, rt, rt, t1, t2, z_cv_ss_2, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
 
 @test z_cv_se_2 ≈ z_cv_ss_2 atol=1e-7
 
-BEAST.momintegrals!(op3, rt, rt, t1, t2, z_cv_se_3, SE_strategy)
-BEAST.momintegrals!(op3, rt, rt, t1, t2, z_cv_ss_3, SS_strategy)
+BEAST.momintegrals!(op3, rt, rt, t1, t2, z_cv_se_3, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+BEAST.momintegrals!(op3, rt, rt, t1, t2, z_cv_ss_3, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
 
 @show z_cv_se_3
 @show z_cv_ss_3
@@ -119,26 +119,26 @@ z_ce_ss_2 = zeros(3,3)
 z_ce_se_3 = zeros(3,3)
 z_ce_ss_3 = zeros(3,3)
 
-BEAST.momintegrals!(op1, rt, rt, t1, t2, z_ce_se_1, SE_strategy)
-BEAST.momintegrals!(op1, rt, rt, t1, t2, z_ce_ss_1, SS_strategy)
+BEAST.momintegrals!(op1, rt, rt, t1, t2, z_ce_se_1, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+BEAST.momintegrals!(op1, rt, rt, t1, t2, z_ce_ss_1, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
 
 @show norm(z_ce_se_1 - z_ce_ss_1)
 @test z_ce_se_1 ≈ z_ce_ss_1 atol=1e-5
 
-BEAST.momintegrals!(op2, rt, rt, t1, t2, z_ce_se_2, SE_strategy)
-BEAST.momintegrals!(op2, rt, rt, t1, t2, z_ce_ss_2, SS_strategy)
+BEAST.momintegrals!(op2, rt, rt, t1, t2, z_ce_se_2, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+BEAST.momintegrals!(op2, rt, rt, t1, t2, z_ce_ss_2, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
 
 @show norm(z_ce_se_2 - z_ce_ss_2)
 @test z_ce_se_2 ≈ z_ce_ss_2 atol=1e-5
 
-BEAST.momintegrals!(op3, rt, rt, t1, t2, z_ce_se_3, SE_strategy)
-BEAST.momintegrals!(op3, rt, rt, t1, t2, z_ce_ss_3, SS_strategy)
+BEAST.momintegrals!(op3, rt, rt, t1, t2, z_ce_se_3, SE_strategy, BEAST.quadraturebuffer(SE_strategy))
+BEAST.momintegrals!(op3, rt, rt, t1, t2, z_ce_ss_3, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
 @show norm(z_ce_se_3 - z_ce_ss_3)
 @test z_ce_se_3 ≈ z_ce_ss_3 atol=1e-5
 
 SS_strategy = SauterSchwabQuadrature.CommonEdge(BEAST._legendre(18,T(0.0),T(1.0)))
 z_ce_ss_3_18 = zeros(3,3)
-BEAST.momintegrals!(op3, rt, rt, t1, t2, z_ce_ss_3_18, SS_strategy)
+BEAST.momintegrals!(op3, rt, rt, t1, t2, z_ce_ss_3_18, SS_strategy, BEAST.quadraturebuffer(SS_strategy))
 @show norm(z_ce_ss_3 - z_ce_ss_3_18)
 @test z_ce_ss_3 ≈ z_ce_ss_3_18 atol=1e-14
-#end
\ No newline at end of file
+#end
diff --git a/test/test_ss_nested_meshes.jl b/test/test_ss_nested_meshes.jl
index 9e71a826..8b5160ea 100644
--- a/test/test_ss_nested_meshes.jl
+++ b/test/test_ss_nested_meshes.jl
@@ -58,7 +58,8 @@ out_ss = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, dom
 BEAST.momintegrals!(out_ss, 𝒜,
     Y, p, test_chart,
     X, 1, trial_chart,
-    BEAST.TestRefinesTrialQRule(qs_strat))
+    BEAST.TestRefinesTrialQRule(qs_strat),
+    BEAST.quadraturebuffer(BEAST.TestRefinesTrialQRule(qs_strat)))
 
 
 wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints))
@@ -66,7 +67,8 @@ out_dw = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, dom
 BEAST.momintegrals!(out_dw, 𝒜,
     Y, p, test_chart,
     X, 1, trial_chart,
-    wiltonsingext)
+    wiltonsingext,
+    BEAST.quadraturebuffer(wiltonsingext))
 
 @show norm(out_ss-out_dw) / norm(out_dw)
 
@@ -78,11 +80,11 @@ trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]
 
 sauterschwab = BEAST.SauterSchwabQuadrature.CommonFace(BEAST._legendre(10,0.0,1.0))
 out_ss1 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss1,sauterschwab)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss1,sauterschwab, BEAST.quadraturebuffer(sauterschwab))
 
 wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints))
 out_dw1 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw1,wiltonsingext)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw1,wiltonsingext, BEAST.quadraturebuffer(wiltonsingext))
 
 @test out_ss1 ≈ out_dw1 rtol=3e-6
 
@@ -95,11 +97,11 @@ trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]
 
 sauterschwab = BEAST.SauterSchwabQuadrature.CommonEdge(BEAST._legendre(10,0.0,1.0))
 out_ss2 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss2,sauterschwab)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss2,sauterschwab, BEAST.quadraturebuffer(sauterschwab))
 
 wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints))
 out_dw2 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw2,wiltonsingext)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw2,wiltonsingext, BEAST.quadraturebuffer(wiltonsingext))
 
 @test out_ss2 ≈ out_dw2 rtol=4e-6
 
@@ -112,11 +114,11 @@ trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]
 
 sauterschwab = BEAST.SauterSchwabQuadrature.CommonEdge(BEAST._legendre(10,0.0,1.0))
 out_ss3 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss3,sauterschwab)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss3,sauterschwab, BEAST.quadraturebuffer(sauterschwab))
 
 wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints))
 out_dw3 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw3,wiltonsingext)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw3,wiltonsingext, BEAST.quadraturebuffer(wiltonsingext))
 
 @test out_ss3 ≈ out_dw3 rtol=3e-6
 
@@ -130,11 +132,11 @@ trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]
 
 sauterschwab = BEAST.SauterSchwabQuadrature.CommonVertex(BEAST._legendre(10,0.0,1.0))
 out_ss4 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss4,sauterschwab)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss4,sauterschwab, BEAST.quadraturebuffer(sauterschwab))
 
 wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints))
 out_dw4 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw4,wiltonsingext)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw4,wiltonsingext, BEAST.quadraturebuffer(wiltonsingext))
 
 @test out_ss4 ≈ out_dw4 rtol=2e-9
 
@@ -147,11 +149,11 @@ trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]
 
 sauterschwab = BEAST.SauterSchwabQuadrature.CommonVertex(BEAST._legendre(10,0.0,1.0))
 out_ss5 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss5,sauterschwab)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss5,sauterschwab, BEAST.quadraturebuffer(sauterschwab))
 
 wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints))
 out_dw5 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw5,wiltonsingext)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw5,wiltonsingext, BEAST.quadraturebuffer(wiltonsingext))
 
 @test out_ss4 ≈ out_dw4 rtol=3e-8
 
@@ -164,11 +166,11 @@ trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1]
 
 sauterschwab = BEAST.SauterSchwabQuadrature.CommonVertex(BEAST._legendre(10,0.0,1.0))
 out_ss6 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss6,sauterschwab)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss6,sauterschwab, BEAST.quadraturebuffer(sauterschwab))
 
 wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints))
 out_dw6 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart)))
-BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw6,wiltonsingext)
+BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw6,wiltonsingext, BEAST.quadraturebuffer(wiltonsingext))
 
 @test out_ss4 ≈ out_dw4 rtol=6e-8
 
diff --git a/test/test_wiltonints.jl b/test/test_wiltonints.jl
index 1f1469b7..9bb8d457 100644
--- a/test/test_wiltonints.jl
+++ b/test/test_wiltonints.jl
@@ -345,7 +345,8 @@ bqd = BE.quadpoints(x, [s], (13,))
 DQ_strategy = BE.DoubleQuadRule(tqd[1,1], bqd[1,1])
 BEAST.momintegrals!(z1, op,
     X, nothing, t,
-    X, nothing, s, DQ_strategy)
+    X, nothing, s, DQ_strategy,
+    BEAST.quadraturebuffer(DQ_strategy))
 
 SE_strategy = BE.WiltonSERule(
   tqd[1,1],
@@ -356,6 +357,7 @@ SE_strategy = BE.WiltonSERule(
 )
 BEAST.momintegrals!(z2, op,
     X, nothing, t,
-    X, nothing, s, SE_strategy)
+    X, nothing, s, SE_strategy,
+    BEAST.quadraturebuffer(SE_strategy))
 
 @test norm(z1-z2) < 1.0e-7
diff --git a/test/untested/test_mwops.jl b/test/untested/test_mwops.jl
index b2ec8586..65361d80 100644
--- a/test/untested/test_mwops.jl
+++ b/test/untested/test_mwops.jl
@@ -63,5 +63,5 @@ qd = quaddata(K, rt, rt, elements(τ), elements(σ))
 zlocal = zeros(promote_type(scalartype(K), scalartype(S), scalartype(T)), 3, 3)
 strat = BEAST.quadrule(K, rt, rt, 1, t, 1, s, qd)
 strat =  BEAST.DoubleQuadRule(qd.tpoints[1,1],  qd.bpoints[1,1])
-BEAST.momintegrals!(K, rt, rt, t, s, zlocal, strat)
+BEAST.momintegrals!(K, rt, rt, t, s, zlocal, strat, BEAST.quadraturebuffer(strat))
 z = zlocal[3,3]