diff --git a/examples/composedoperator.jl b/examples/composedoperator.jl
new file mode 100644
index 00000000..5cf9f764
--- /dev/null
+++ b/examples/composedoperator.jl
@@ -0,0 +1,21 @@
+##### user manual composed operator and potential
+
+using BEAST
+using CompScienceMeshes
+
+Γ = meshcuboid(1.0,1.0,1.0,0.2);
+X = lagrangecxd0(Γ)
+
+S = BEAST.PotentialIntegralOperator{2}(BEAST.HH3DGreen(1im),*,B->B)
+St = BEAST.trace(S,1.0)
+assemble(St,X,X)
+
+SB = Helmholtz3D.singlelayer(wavenumber=1.0)
+assemble(SB,X,X)
+
+Y = raviartthomas(Γ)
+K = BEAST.PotentialIntegralOperator{2}(BEAST.HH3DGradGreen(0*1im),×,B->B)
+Kt = BEAST.ttrace(K,1.0;testfunction_tangential=true)
+assemble(Kt,Y,Y)
+
+KB = Maxwell3D.doublelayer(wavenumber=0.0)
\ No newline at end of file
diff --git a/examples/vp_global_multi_trace.jl b/examples/vp_global_multi_trace.jl
new file mode 100644
index 00000000..553f7a2d
--- /dev/null
+++ b/examples/vp_global_multi_trace.jl
@@ -0,0 +1,168 @@
+using BEAST
+using CompScienceMeshes
+using LinearAlgebra
+using SparseArrays
+using PlotlyJS
+using StaticArrays
+
+ω = 10.0^8*2*pi
+ϵ0 = 8.854e-12
+μ0 = 4π*1e-7
+κ0 = ω*√(μ0*ϵ0)
+
+ϵr = [1,2,2]
+μr = [2,1,2]
+
+h = 0.5
+fn = joinpath(@__DIR__, "assets/rectangular_torus.geo")
+Γ11 = CompScienceMeshes.meshgeo(fn; dim=2, h)
+Γ11 = Mesh([point(-z,y,x) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11)))
+function onlys(a)
+    c = sort.(a)
+    b = [i for i in a if length(findall(==(sort(i)),c))==1]
+    return b
+end
+
+Γ12 = -Mesh([point(x,y,-z-4) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11)))
+Γ1 = weld(Γ11,Γ12;glueop=onlys)
+Γ2 = Mesh([point(x,-y-4,-z-4) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11)))
+Γ3 = -Mesh([point(x,-y-4,z) for (x,y,z) in vertices(Γ11)], deepcopy(cells(Γ11)))
+
+
+Γ = [Γ1,Γ2,Γ3] # the meshes
+
+κ = [ω*√(μ0*ϵ0*ϵri*μri) for (ϵri,μri) in zip(ϵr,μr)]
+
+########RHS definition (lorentz gauge:)
+x0 = @SVector [0.0,0.0,0.0] # used to check the type of the incidnet fields
+A(x) = x[1]*sqrt(ϵ0*μ0)exp(-1im*κ0*x[3])* (@SVector [0.0,0.0,1.0])
+curlA(x) = -(@SVector [0.0,1.0,0.0])*sqrt(ϵ0*μ0)*exp(-1.0im*κ0 *x[3])
+phi(x) =  x[1]*exp(-1.0im*κ0 *x[3])
+gradphi(x) = (@SVector [1.0,0.0,-1im*κ0*x[1]])*exp(-1.0im*κ0 *x[3])
+
+####################################################
+#
+# From here on do nothing
+#
+####################################################
+@hilbertspace a1 a2
+@hilbertspace b1 b2
+
+#### potentials + operators
+G0 = BEAST.HH3DGreen(1im*κ0)
+G = [BEAST.HH3DGreen(1im*k) for k in κ]
+dG0 = BEAST.HH3DGradGreen(1im*κ0)
+dG = [BEAST.HH3DGradGreen(1im*k) for k in κ]
+
+    U11 = -BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0)
+    U12 = BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->n*b),-1.0)
+    U22 = BEAST.trace(BEAST.PotentialIntegralOperator{2}(dG0,(x,y)->transpose(x)*y,b->n*b),-1.0)
+
+DDL0 = U11[a1,b1] + U12[a1,b2] + U22[a2,b2]
+
+    U11 = [-BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG]
+    U12 = [BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,*,b->n*b),1.0) for g in G]
+    U22 = [BEAST.trace(BEAST.PotentialIntegralOperator{2}(g,(x,y)->transpose(x)*y,b->n*b),1.0) for g in dG]
+
+DDL = [U11i[a1,b1] + er*mr*U12i[a1,b2] + U22i[a2,b2] for (U11i,U12i,U22i,er,mr) in zip(U11,U12,U22,ϵr,μr)]
+
+    U11 = -BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0)
+    U21 = -BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0)
+    U22 = BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(dG0,*,b->b),-1.0)
+
+NDL0 = U11[a1,b1] + U21[a2,b1] + U22[a2,b2]
+
+    U11 = [-BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG]
+    U21 = [-BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in G]
+    U22 = [BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in dG]
+
+NDL = [U11i[a1,b1] + er*mr*U21i[a2,b1] + U22i[a2,b2] for (U11i,U21i,U22i,er,mr) in zip(U11,U21,U22,ϵr,μr)]
+
+    U11 = -BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0)
+    U12 = BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(dG0,*,b->b),-1.0)
+    U21 = -BEAST.trace(BEAST.PotentialIntegralOperator{2}(dG0,(x,y)->transpose(x)*y,b->b),-1.0)
+    U22 = -κ0^2* BEAST.trace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b),-1.0)
+
+NDSL0 = U11[a1,b1] + U12[a1,b2] + U21[a2,b1] + U22[a2,b2]
+
+    U11 = [-BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in G]
+    U12 = [BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for g in dG]
+    U21 = [-BEAST.trace(BEAST.PotentialIntegralOperator{2}(g,(x,y)->transpose(x)*y,b->b),1.0) for g in dG]
+    U22 = [-k^2*BEAST.trace(BEAST.PotentialIntegralOperator{2}(g,*,b->b),1.0) for (g,k) in zip(G,κ)]
+
+NDSL = [mr*U11i[a1,b1] + 1/er*U12i[a1,b2] + 1/er*U21i[a2,b1] + 1/(er^2*mr)*U22i[a2,b2] for (U11i,U12i,U21i,U22i,er,mr) in zip(U11,U12,U21,U22,ϵr,μr)]
+
+    U11 = -κ0^2 * BEAST.pvrtrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->b))-BEAST.CompDoubleInt(B->divergence(n×B),*,G0,*,B->divergence(B))
+    U12 = BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->n*b),-1.0)
+    U21 = -BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(dG0,×,b->b),-1.0)
+    U22 =  BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(G0,*,b->n*b),-1.0)
+
+DNSL0 = U11[a1,b1] + U12[a1,b2] + U21[a2,b1] + U22[a2,b2]
+
+    U11 = [-k^2 * BEAST.pvrtrace(BEAST.PotentialIntegralOperator{2}(g,*,b->b))-BEAST.CompDoubleInt(B->divergence(n×B),*,g,*,B->divergence(B)) for (g,k) in zip(G,κ)]
+    U12 = [BEAST.rtrace(BEAST.PotentialIntegralOperator{2}(g,×,b->n*b),1.0) for g in dG]
+    U21 = [-BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,×,b->b),1.0) for g in dG]
+    U22 = [BEAST.ntrace(BEAST.PotentialIntegralOperator{2}(g,*,b->n*b),1.0) for (g,k) in zip(G,κ)]
+
+DNSL = [1/mr*U11i[a1,b1] + er*U12i[a1,b2] + er*U21i[a2,b1] + er^2*mr*U22i[a2,b2] for (U11i,U12i,U21i,U22i,er,mr) in zip(U11,U12,U21,U22,ϵr,μr)]
+
+### definition per domain bilinear forms
+
+A0 = DNSL0[a1,b1] + NDSL0[a2,b2] + DDL0[a2,b1] + NDL0[a1,b2]
+Ad = [DNSLi[a1,b1] + NDSLi[a2,b2] + DDLi[a2,b1] + NDLi[a1,b2] for (DNSLi,NDSLi,DDLi,NDLi) in zip(DNSL,NDSL,DDL,NDL)]
+
+p = BEAST.hilbertspace(:p, length(κ))
+q = BEAST.hilbertspace(:q, length(κ))
+
+LHSA = A0[p,q] + BEAST.Variational.DirectProductKernel(Ad)[p,q]
+
+####################################################
+#
+# RHS definition
+#
+####################################################
+
+### wrapping incident field
+
+w_A = BEAST.FunctionWrapper{typeof(A(x0))}(A)
+w_curlA = BEAST.FunctionWrapper{typeof(curlA(x0))}(curlA)
+w_phi = BEAST.FunctionWrapper{typeof(phi(x0))}(phi)
+w_gradphi = BEAST.FunctionWrapper{typeof(gradphi(x0))}(gradphi)
+
+### traces RHS
+
+t_nxA = BEAST.CrossTraceMW(w_A)
+t_nxcurlA = BEAST.CrossTraceMW(w_curlA)
+t_ndA = BEAST.NDotTrace(w_A)
+t_div_rescA =  BEAST.Trace(w_phi)
+
+t_d = t_nxA[a1] + -1im*κ0^2/ω*t_div_rescA[a2] 
+t_n = t_nxcurlA[a1] + t_ndA[a2]
+RHSAi = t_d[a2] + t_n[a1]
+
+t_phi = BEAST.Trace(w_phi)
+t_ngradphi = BEAST.NDotTrace(w_gradphi)
+
+RHSA = RHSAi[p]
+
+### Spaces
+RT = raviartthomas.(Γ)
+ND = Ref(n) .× RT
+MND = Ref(n) .× (Ref(n) .× (Ref(n) .× RT))
+L0 = lagrangecxd0.(Γ)
+L1 = lagrangec0d1.(Γ)
+
+Xd = [RTi×L1i for (RTi,L1i) in zip(RT,L1)]
+Xn = [RTi×L0i for (RTi,L0i) in zip(RT,L0)]
+
+Yn = [NDi×L1i for (NDi,L1i) in zip(MND,L1)]
+Yd = [NDi×L0i for (NDi,L0i) in zip(ND,L0)]
+
+Qₕ = [BEAST.DirectProductSpace([Xdi,Xni]) for (Xdi,Xni) in zip(Xd,Xn)]
+Pₕ = [BEAST.DirectProductSpace([Yni,Ydi]) for (Ydi,Yni) in zip(Yd,Yn)]
+
+deq = BEAST.discretise(LHSA==RHSA, (p .∈ Pₕ)..., (q .∈ Qₕ)...)
+# Z = assemble(LHSA,BEAST.DirectProductSpace(Pₕ),BEAST.DirectProductSpace(Qₕ))
+# B = assemble(RHSA,BEAST.DirectProductSpace(Pₕ))
+u = solve(deq)
+
diff --git a/src/BEAST.jl b/src/BEAST.jl
index 879a8e94..5425ad7d 100644
--- a/src/BEAST.jl
+++ b/src/BEAST.jl
@@ -190,6 +190,9 @@ include("bases/stagedtimestep.jl")
 include("bases/timebasis.jl")
 include("bases/tensorbasis.jl")
 
+include("bases/composedbasis.jl")
+include("bases/local/localcomposedbasis.jl")
+
 include("operator.jl")
 
 include("quadrature/quadstrats.jl")
@@ -286,8 +289,11 @@ include("solvers/itsolver.jl")
 
 include("utils/plotlyglue.jl")
 
-
-
+include("composedoperators/composedoperator.jl")
+include("composedoperators/displacementmesh.jl")
+include("composedoperators/potentials.jl")
+include("composedoperators/trace.jl")
+include("composedoperators/analytic_excitation.jl")
 
 const x̂ = point(1, 0, 0)
 const ŷ = point(0, 1, 0)
diff --git a/src/bases/composedbasis.jl b/src/bases/composedbasis.jl
new file mode 100644
index 00000000..de2e135c
--- /dev/null
+++ b/src/bases/composedbasis.jl
@@ -0,0 +1,76 @@
+struct FunctionWrapper{T} 
+    func::Function
+end
+
+function (f::FunctionWrapper{T})(x)::T where {T}
+    return f.func(x)
+end
+
+function (f::FunctionWrapper{T})(x::CompScienceMeshes.MeshPointNM)::T where {T}
+    return f.func(cartesian(x))
+end
+scalartype(F::FunctionWrapper{T}) where {T} = eltype(T)
+
+
+abstract type _BasisOperations{T} <: Space{T} end
+
+struct _BasisTimes{T} <: _BasisOperations{T}
+    el1
+    el2
+end
+_BasisTimes(el1::NormalVector,el2::Space) = _BasisTimes{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el1,el2)
+_BasisTimes(el2::Space,el1::NormalVector) = _BasisTimes{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el2,el1)
+_BasisTimes(el1::Space,el2::FunctionWrapper) = _BasisTimes{promote_type(scalartype(el1),scalartype(el2))}(el1,el2)
+_BasisTimes(el1::FunctionWrapper,el2::Space) = _BasisTimes{promote_type(scalartype(el1),scalartype(el2))}(el1,el2)
+
+struct _BasisCross{T} <: _BasisOperations{T}
+    el1
+    el2
+end
+_BasisCross(el1::NormalVector,el2::Space) = _BasisCross{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el1,el2)
+_BasisCross(el2::Space,el1::NormalVector) = _BasisCross{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el2,el1)
+_BasisCross(el1::FunctionWrapper,el2::Space) = _BasisCross{promote_type(scalartype(el1),scalartype(el2))}(el1,el2)
+_BasisCross(el1::Space,el2::FunctionWrapper) = _BasisCross{promote_type(scalartype(el1),scalartype(el2))}(el1,el2)
+
+struct _BasisDot{T} <: _BasisOperations{T}
+    el1
+    el2
+end
+_BasisDot(el1::NormalVector,el2::Space) = _BasisDot{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el1,el2)
+_BasisDot(el2::Space,el1::NormalVector) = _BasisDot{promote_type(eltype(vertextype(geometry(el2))),scalartype(el2))}(el2,el1)
+_BasisDot(el1::Space,el2::FunctionWrapper) = _BasisDot{promote_type(scalartype(el1),scalartype(el2))}(el1,el2)
+_BasisDot(el1::FunctionWrapper,el2::Space) = _BasisDot{promote_type(scalartype(el1),scalartype(el2))}(el1,el2)
+
+refspace(a::_BasisTimes{T}) where {T} = _LocalBasisTimes(T,refspace(a.el1),refspace(a.el2))
+refspace(a::_BasisCross{T}) where {T} = _LocalBasisCross(T,refspace(a.el1),refspace(a.el2))
+refspace(a::_BasisDot{T}) where {T} = _LocalBasisDot(T,refspace(a.el1),refspace(a.el2))
+
+refspace(a::FunctionWrapper) = a 
+refspace(a::NormalVector) = a
+
+numfunctions(a::Union{NormalVector,FunctionWrapper}) = missing
+numfunctions(a::Union{NormalVector,FunctionWrapper},simp) = missing
+numfunctions(a::_BasisOperations) = coalesce(numfunctions(a.el1) , numfunctions(a.el2))
+
+geometry(a::_BasisOperations) = coalesce(geometry(a.el1),geometry(a.el2))
+basisfunction(a::_BasisOperations,i) = coalesce(basisfunction(a.el1,i),basisfunction(a.el2,i))
+positions(a::_BasisOperations) = coalesce(positions(a.el1),positions(a.el2))
+
+geometry(a::Union{NormalVector,FunctionWrapper}) = missing
+basisfunction(a::Union{NormalVector,FunctionWrapper},i) = missing
+positions(a::Union{NormalVector,FunctionWrapper}) = missing
+
+
+subset(a::T,I) where {T <: _BasisOperations} = T(subset(a.el1,I),subset(a.el2,I))
+subset(a::Union{NormalVector,FunctionWrapper},I) = a
+
+#### mathematical support will be called if there are no dedicated routines
+×(n::NormalVector,s::Space) = _BasisCross(n,s)
+×(s::Space,n::NormalVector) = _BasisCross(s,n)
+
+⋅(n::NormalVector,s::Space) = _BasisDot(n,s)
+⋅(s::Space,n::NormalVector) = _BasisDot(s,n)
+
+*(n::NormalVector,s::Space) = _BasisTimes(n,s)
+*(s::Space,n::NormalVector) = _BasisTimes(s,n)
+
diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl
index cbc3a617..552f21cb 100644
--- a/src/bases/lagrange.jl
+++ b/src/bases/lagrange.jl
@@ -296,10 +296,15 @@ function singleduallagd0(fine, F, v; interpolatory=false)
 
     T = coordtype(fine)
     fn = Shape{T}[]
+    vol = T(0.0)
+    for cellid in F
+        ptch = chart(fine, cellid)
+        vol += volume(ptch)
+    end
     for cellid in F
         # cell = cells(fine)[cellid]
         ptch = chart(fine, cellid)
-        coeff = interpolatory ? T(1.0) : 1 / volume(ptch) / length(F)
+        coeff = interpolatory ? T(1.0) : 1 / vol
         refid = 1
         push!(fn, Shape(cellid, refid, coeff))
     end
@@ -498,7 +503,6 @@ end
 
 
 function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}})
-
     T = coordtype(mesh)
     num_faces = dimension(mesh)+1
 
@@ -559,7 +563,9 @@ function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}})
                 fine_idcs = cells(refined)[c]
                 local_id = something(findfirst(isequal(v), fine_idcs),0)
                 @assert local_id != 0
-                shape = Shape(c, local_id, 1/n/2)
+                #shape = Shape(c, local_id, 1/n/2)
+                shape = Shape(c, local_id, 1/2)
+                
                 push!(fns[i], shape)
             end
         end
@@ -572,7 +578,7 @@ function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}})
                 fine_idcs = cells(refined)[c]
                 local_id = something(findfirst(isequal(v), fine_idcs),0)
                 @assert local_id != 0
-                shape = Shape(c, local_id, 1/n/2)
+                shape = Shape(c, local_id, 1/(n/2))
                 push!(fns[i], shape)
             end
         end
@@ -1045,4 +1051,22 @@ function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order)
     for f in mesh pos[nV + nE + (f-1)*nf + 1: nV + nE + f*nf] .= Ref(cartesian(center(chart(mesh, f)))) end
 
     return LagrangeBasis{order,0,localdim}(mesh, fns, pos)
+end
+
+function _surface(X)
+    C = unitfunctioncxd0(X.geo)
+    G = assemble(Identity(),X,X)
+    u = Vector(assemble(Identity(),X,C)[1:end,1])
+    dot(u,G\u)
+end
+@testitem "duallagrangec0d1" begin
+    using CompScienceMeshes
+    Γ = meshcuboid(1.0,1.0,1.0,0.5)
+    @test BEAST._surface(duallagrangec0d1(Γ)) ≈ 6.0
+end
+
+@testitem "duallagrangecxd0" begin
+    using CompScienceMeshes
+    Γ = meshcuboid(1.0,1.0,1.0,0.5)
+    @test BEAST._surface(duallagrangecxd0(Γ)) ≈ 6.0
 end
\ No newline at end of file
diff --git a/src/bases/local/localcomposedbasis.jl b/src/bases/local/localcomposedbasis.jl
new file mode 100644
index 00000000..447539b7
--- /dev/null
+++ b/src/bases/local/localcomposedbasis.jl
@@ -0,0 +1,89 @@
+abstract type _LocalBasisOperations{T} <: RefSpace{T} end
+numfunctions(a::_LocalBasisOperations) = coalesce(numfunctions(a.el1) , numfunctions(a.el2))
+numfunctions(a::_LocalBasisOperations,simp) = coalesce(numfunctions(a.el1,simp) , numfunctions(a.el2,simp))
+struct TriangleSupport end
+struct TetraSupport end
+
+
+support(a::_LocalBasisOperations) = coalesce(support(a.el1),support(a.el2))
+support(a::NormalVector) = missing
+support(a::FunctionWrapper) = missing
+support(a::LagrangeRefSpace{T,N,3}) where {T,N} = TriangleSupport()
+support(a::LagrangeRefSpace{T,N,4}) where {T,N} = TetraSupport()
+restrict(a::_LocalBasisOperations,bcell,cell) = restrict(_refspace(a),bcell,cell)
+
+_refspace(a::_LocalBasisOperations) = coalesce(_refspace(a.el1) , _refspace(a.el2))
+_refspace(a::NormalVector) = missing
+_refspace(a::FunctionWrapper) = missing
+_refspace(a::RefSpace) = a
+struct _LocalBasisTimes{T,U,V} <: _LocalBasisOperations{T}
+    el1::U
+    el2::V
+    function _LocalBasisTimes(::Type{T},el1::U,el2::V) where {T,U,V}
+        new{T,U,V}(el1,el2)
+    end
+end
+
+struct _LocalBasisCross{T,U,V} <: _LocalBasisOperations{T}
+    el1::U
+    el2::V
+    function _LocalBasisCross(::Type{T},el1::U,el2::V) where {T,U,V}
+        new{T,U,V}(el1,el2)
+    end
+end
+
+struct _LocalBasisDot{T,U,V} <: _LocalBasisOperations{T}
+    el1::U
+    el2::V
+    function _LocalBasisDot(::Type{T},el1::U,el2::V) where {T,U,V}
+        new{T,U,V}(el1,el2)
+    end
+end
+
+function (op::U where {U <: _LocalBasisOperations})(p)
+    l = op.el1(p)
+    r = op.el2(p)
+
+    _modbasis(l,r) do a,b
+        operation(op)(a,b)
+    end
+end
+
+function (op::NormalVector)(x::CompScienceMeshes.MeshPointNM)
+    return normal(x)
+end
+
+function _modbasis_genleft(::Type{U}, ::Type{V}) where {U<:SVector{N}, V} where {N}
+    ex = :(SVector{N}(()))
+    for n in 1:N
+
+        push!(ex.args[2].args, :(value = f(a[$n].value, b),))
+    end
+
+    return ex
+end
+function _modbasis_genright(::Type{U}, ::Type{V}) where {V<:SVector{N}, U} where {N}
+    ex = :(SVector{N}(()))
+
+    for n in 1:N
+
+        push!(ex.args[2].args, :(value = f(a, b[$n].value),))
+    end
+
+    return ex
+end
+@generated function _modbasis(f, a::SVector{N,T}, b::U) where {N, T <:NamedTuple,U}
+    ex = _modbasis_genleft(a,b)
+
+    return ex
+end
+@generated function _modbasis(f, a::U, b::SVector{N,T}) where {N,T<:NamedTuple,U}
+    ex = _modbasis_genright(a,b)
+
+    return ex
+end
+
+operation(a::_LocalBasisTimes) = *
+operation(a::_LocalBasisCross) = ×
+operation(a::_LocalBasisDot) = (x,y) --> transpose(x)*y
+
diff --git a/src/bases/ndspace.jl b/src/bases/ndspace.jl
index 0393e5bb..10042d8d 100644
--- a/src/bases/ndspace.jl
+++ b/src/bases/ndspace.jl
@@ -40,12 +40,27 @@ function nedelec(surface, edges=skeleton(surface,1))
     NDBasis(surface, fns, pos)
 end
 
+# function LinearAlgebra.cross(::NormalVector, s::NDBasis)
+#     # @assert CompScienceMeshes.isoriented(s.geo)
+#     RTBasis(s.geo, s.fns, s.pos)
+# end
+
 function LinearAlgebra.cross(::NormalVector, s::NDBasis)
-    # @assert CompScienceMeshes.isoriented(s.geo)
-    RTBasis(s.geo, s.fns, s.pos)
+    @assert CompScienceMeshes.isoriented(s.geo)
+    fns = similar(s.fns)
+    for (i,fn) in pairs(s.fns)
+        fns[i] = [Shape(sh.cellid, sh.refid, -sh.coeff) for sh in fn]
+    end
+    RTBasis(s.geo, fns, s.pos)
+end
+function LinearAlgebra.cross(s::NDBasis, ::NormalVector)
+    @assert CompScienceMeshes.isoriented(s.geo)
+    fns = similar(s.fns)
+    for (i,fn) in pairs(s.fns)
+        fns[i] = [Shape(sh.cellid, sh.refid, sh.coeff) for sh in fn]
+    end
+    RTBasis(s.geo, fns, s.pos)
 end
-
-
 function curl(space::NDBasis)
-    divergence(n × space)
+    divergence(n × (n × (n × space)))
 end
diff --git a/src/bases/rtspace.jl b/src/bases/rtspace.jl
index b6a280a6..b4e4d42f 100644
--- a/src/bases/rtspace.jl
+++ b/src/bases/rtspace.jl
@@ -264,7 +264,15 @@ function LinearAlgebra.cross(::NormalVector, s::RTBasis)
     @assert CompScienceMeshes.isoriented(s.geo)
     fns = similar(s.fns)
     for (i,fn) in pairs(s.fns)
-        fns[i] = [Shape(sh.cellid, sh.refid, -sh.coeff) for sh in fn]
+        fns[i] = [Shape(sh.cellid, sh.refid, sh.coeff) for sh in fn]
     end
     NDBasis(s.geo, fns, s.pos)
 end
+function LinearAlgebra.cross(s::RTBasis,::NormalVector)
+    @assert CompScienceMeshes.isoriented(s.geo)
+    fns = similar(s.fns)
+    for (i,fn) in pairs(s.fns)
+        fns[i] = [Shape(sh.cellid, sh.refid, -sh.coeff) for sh in fn]
+    end
+    NDBasis(s.geo, fns, s.pos)
+end
\ No newline at end of file
diff --git a/src/composedoperators/analytic_excitation.jl b/src/composedoperators/analytic_excitation.jl
new file mode 100644
index 00000000..48163f17
--- /dev/null
+++ b/src/composedoperators/analytic_excitation.jl
@@ -0,0 +1,10 @@
+struct Trace{T,F} <: Functional
+    field::F
+end
+  
+Trace(f::F) where {F} = Trace{scalartype(f), F}(f)
+Trace{T}(f::F) where {T,F} = Trace{T,F}(f)
+scalartype(s::Trace{T}) where {T} = T
+
+(ϕ::Trace)(p) = ϕ.field(cartesian(p))
+integrand(::Trace, g, ϕ) = *(g.value, ϕ)
diff --git a/src/composedoperators/composedoperator.jl b/src/composedoperators/composedoperator.jl
new file mode 100644
index 00000000..13a5e9f3
--- /dev/null
+++ b/src/composedoperators/composedoperator.jl
@@ -0,0 +1,249 @@
+
+abstract type AbstractCompInt <: AbstractOperator end
+abstract type Kernel{T} end
+
+#Make these structures, those are not yet pulled trough the smartmath function
+struct CompDoubleInt{T} <: AbstractCompInt
+    tfunc
+    op1
+    kernel::Kernel{T}
+    op2
+    bfunc
+end
+
+struct CompSingleInt{T} <: AbstractCompInt
+    tfunc
+    op1
+    kernel::Kernel{T}
+    op2
+    bfunc
+end
+scalartype(op::CompDoubleInt{T}) where {T} = T 
+scalartype(op::CompSingleInt{T}) where {T} = T 
+
+
+#the functionallity is pased to the basis at this point.
+struct CompDoubleKern{T,U,V,W} <: IntegralOperator
+    op1::U
+    kernel::W
+    op2::V
+end
+function CompDoubleKern(p1,k::Kernel{T},p2) where {T}
+    CompDoubleKern{T,typeof(p1),typeof(p2),typeof(k)}(p1,k,p2)
+end
+struct CompSingleKern{T,U,V,W} <: LocalOperator
+    op1::U
+    kernel::W
+    op2::V
+end
+function CompSingleKern(p1,k::Kernel{T},p2) where {T}
+    CompSingleKern{T,typeof(p1),typeof(p2),typeof(k)}(p1,k,p2)
+end
+scalartype(op::CompDoubleKern{T}) where {T} = T 
+scalartype(op::CompSingleKern{T}) where {T} = T 
+
+integralop(a::CompDoubleInt) = CompDoubleKern(a.op1,a.kernel,a.op2)
+integralop(a::CompSingleInt) = CompSingleKern(a.op1,a.kernel,a.op2)
+
+
+function assemble!(op::AbstractCompInt, tfs::AbstractSpace, bfs::AbstractSpace,
+    store, threading = Threading{:multi};
+    quadstrat=defaultquadstrat(op, tfs, bfs))
+
+    assemble!(integralop(op),op.tfunc(tfs),op.bfunc(bfs),store,threading;quadstrat)
+
+end
+
+
+
+####### kernel definition (also posibilities to do symbolic actions)
+
+# struct NullKernel <: Kernel{Int} end
+# struct NxDotNy{T} <: Kernel{T} 
+#     mult::Kernel{T}
+# end
+# NxDotNy() = NxDotNy(IdentityKernel())
+# NxDotNy(a::NxDotNy) = @error "dont do nx⋅ny*nx⋅ny (yet)"
+
+struct HH3DGreen{T} <: Kernel{T}
+    gamma::T
+end
+struct HH3DGradGreen{T} <: Kernel{T}
+    gamma::T
+end
+
+function (op::HH3DGreen)(x,y)
+    gamma = op.gamma
+
+    r = cartesian(x) - cartesian(y)
+    R = norm(r)
+    iR = 1/R
+    green = exp(-gamma*R)*(i4pi*iR)
+    green
+end
+function (op::HH3DGreen)(x::Union{SVector,Vector},y)
+    gamma = op.gamma
+
+    r = x - cartesian(y)
+    R = norm(r)
+    iR = 1/R
+    green = exp(-gamma*R)*(i4pi*iR)
+    green
+end
+
+function (op::HH3DGradGreen)(x,y)
+    gamma = op.gamma
+
+    r = cartesian(x) - cartesian(y)
+    R = norm(r)
+    iR = 1/R
+    green = exp(-gamma*R)*(iR*i4pi)
+    gradgreen = -(gamma + iR) * green * (iR * r)
+
+    return gradgreen
+end
+function (op::HH3DGradGreen)(x::Union{SVector,Vector},y)
+    gamma = op.gamma
+
+    r = x - cartesian(y)
+    R = norm(r)
+    iR = 1/R
+    green = exp(-gamma*R)*(iR*i4pi)
+    gradgreen = -(gamma + iR) * green * (iR * r)
+    tt = gradgreen
+
+    return tt
+end
+
+##### integrand evaluation
+function integrand(op::CompSingleKern,kerneldata,x,f,g,disp)
+    kernel = op.kernel(x,disp)
+    op1 = op.op1
+    op2 = op.op2
+
+    op1(f.value,op2(kernel,g.value))
+
+end
+
+function (op::Integrand{<:CompDoubleKern})(x,y,f,g)
+    kernel = op.operator.kernel(x,y)
+    op1 = op.operator.op1
+    op2 = op.operator.op2
+    _integrands(f,g) do fi, gi
+        op1(fi.value,op2(kernel,gi.value))
+    end
+end
+
+defaultquadstrat(op::CompDoubleInt,tfs::Space,bfs::Space) = DoubleNumSauterQstrat(5,5,5,5,5,5)
+defaultquadstrat(op::CompSingleInt,tfs::Space,bfs::Space) = SingleNumQStrat(6)
+##### assembling the single integral 
+function assemble!(biop::CompSingleKern, tfs::Space, bfs::Space, store,
+    threading::Type{Threading{:multi}};
+    quadstrat=defaultquadstrat(biop, tfs, bfs))
+
+    return assemble_local_mixed!(biop, tfs, bfs, store; quadstrat)
+end
+
+function assemble_local_mixed!(biop::CompSingleKern, tfs::Space{T}, bfs::Space{T}, store;
+    quadstrat=defaultquadstrat(biop, tfs, bfs)) where {T}
+
+    tol = sqrt(eps(T))
+
+    trefs = refspace(tfs)
+    brefs = refspace(bfs)
+
+    tels, tad = assemblydata(tfs)
+    bels, bad = assemblydata(bfs)
+
+    tgeo = geometry(tfs)
+    bgeo = geometry(bfs)
+
+    tdom = domain(chart(tgeo, first(tgeo)))
+    bdom = domain(chart(bgeo, first(bgeo)))
+
+    num_trefs = numfunctions(trefs, tdom)
+    num_brefs = numfunctions(brefs, bdom)
+
+    qd = quaddata(biop, trefs, brefs, tels, bels, quadstrat)
+
+    # store the bcells in an octree
+    tree = elementstree(bels)
+
+    print("dots out of 10: ")
+    todo, done, pctg = length(tels), 0, 0
+    for (p,tcell) in enumerate(tels)
+
+        tc, ts = boundingbox(tcell.vertices)
+        pred = (c,s) -> boxesoverlap(c,s,tc,ts)
+
+        for box in boxes(tree, pred)
+            for q in box
+                bcell = bels[q]
+
+                if overlap(tcell, bcell)
+
+                    isct = intersection(tcell, bcell)
+                    disp = displacement(displacementchart(geometry(tfs),p),displacementchart(geometry(bfs),q)) #number to multiply with test normal, if zero use the orientation of displacementvector
+                    for cell in isct
+                        volume(cell) < tol && continue
+
+                        P = restrict(brefs, bcell, cell)
+                        Q = restrict(trefs, tcell, cell)
+
+                        qr = quadrule(biop, trefs, brefs, cell, qd, quadstrat)
+                        zlocal = cellinteractions(biop, trefs, brefs, cell, qr, disp)
+                        zlocal = Q * zlocal * P'
+
+                        for i in 1 : num_trefs
+                            for j in 1 : num_brefs
+                                for (m,a) in tad[p,i]
+                                    for (n,b) in bad[q,j]
+                                        store(a * zlocal[i,j] * b, m, n)
+                                    end # next basis function this cell supports
+                                end # next test function this cell supports
+                            end # next refshape on basis side
+                        end # next refshape on test side
+
+                    end # next cell in intersection
+                end # if overlap
+            end # next cell in the basis geometry
+        end # next box in the octree
+
+        done += 1
+        new_pctg = round(Int, done / todo * 100)
+        if new_pctg > pctg + 9
+            print(".")
+            pctg = new_pctg
+        end
+    end # next cell in the test geometry
+
+    println("")
+end
+kernelvals(::CompSingleKern,x) = nothing
+function cellinteractions(biop::CompSingleKern, trefs::U, brefs::V, cell, qr, disp) where {U<:RefSpace{T},V<:RefSpace{T}} where {T}
+
+    num_tshs = length(qr[1][3])
+    num_bshs = length(qr[1][4])
+
+    zlocal = zeros(T, num_tshs, num_bshs)
+    for q in qr
+
+        w, mp, tvals, bvals = q[1], q[2], q[3], q[4]
+        j = w * jacobian(mp)
+        kernel = kernelvals(biop, mp)
+
+        for m in 1 : num_tshs
+            tval = tvals[m]
+
+            for n in 1 : num_bshs
+                bval = bvals[n]
+
+                igd = integrand(biop, kernel, mp, tval, bval,disp)
+                zlocal[m,n] += j * igd
+
+            end
+        end
+    end
+
+    return zlocal
+end
diff --git a/src/composedoperators/displacementmesh.jl b/src/composedoperators/displacementmesh.jl
new file mode 100644
index 00000000..e6079ce4
--- /dev/null
+++ b/src/composedoperators/displacementmesh.jl
@@ -0,0 +1,40 @@
+##### displacement meshes
+abstract type DisplacementMesh{T} <: CompScienceMeshes.AbstractMesh{3,3,T} end
+
+"""
+GlobalDisplacementMesh(mesh,epsilon)
+    every chart is displaced with epsilon allong the normal
+"""
+struct GlobalDisplacementMesh{T} <: DisplacementMesh{T}
+    mesh::Mesh{<:Any,<:Any,T}
+    epsilon 
+end
+struct DisplacementChart
+    chart 
+    epsilon
+end
+CompScienceMeshes.mesh(m::GlobalDisplacementMesh) = m.mesh
+displacementchart(m::GlobalDisplacementMesh,i) = DisplacementChart(chart(m.mesh,i),m.epsilon)
+displacementchart(m::Mesh,i) = DisplacementChart(chart(m,i),-1.0)
+CompScienceMeshes.chart(m::DisplacementMesh,i) = chart(mesh(m),i)
+CompScienceMeshes.cells(m::DisplacementMesh) = cells(mesh(m))
+CompScienceMeshes.celltype(m::DisplacementMesh) = CompScienceMeshes.celltype(mesh(m))
+CompScienceMeshes.indextype(m::DisplacementMesh,a) = CompScienceMeshes.indextype(mesh(m),a)
+CompScienceMeshes.celltype(m::DisplacementMesh,a) = CompScienceMeshes.celltype(mesh(m),a)
+CompScienceMeshes.indices(m::DisplacementMesh,n) = CompScienceMeshes.indices(mesh(m),n)
+CompScienceMeshes.vertices(m::DisplacementMesh) = CompScienceMeshes.vertices(mesh(m))
+CompScienceMeshes.numvertices(m::DisplacementMesh) = CompScienceMeshes.numvertices(mesh(m))
+CompScienceMeshes.numcells(m::DisplacementMesh) = CompScienceMeshes.numcells(mesh(m))
+CompScienceMeshes.vertextype(m::DisplacementMesh) = CompScienceMeshes.vertextype(mesh(m))
+
+
+function displacement(a::DisplacementChart,b::DisplacementChart)
+    rel_orient = sign(dot(normal(a.chart),normal(b.chart)))
+    if a.epsilon*rel_orient ≈ b.epsilon
+        return 0
+    end
+    testdisp = a.epsilon*normal(a.chart)
+    trialdisp = b.epsilon*normal(b.chart)
+    
+    return sign(dot(trialdisp-testdisp,normal(a.chart)))
+end
\ No newline at end of file
diff --git a/src/composedoperators/potentials.jl b/src/composedoperators/potentials.jl
new file mode 100644
index 00000000..e7752b7c
--- /dev/null
+++ b/src/composedoperators/potentials.jl
@@ -0,0 +1,26 @@
+#D: dimension of manifold over which integration is performed
+struct PotentialIntegralOperator{D}
+    kernel
+    op2
+    bfunc
+end
+
+struct PotentialIntegralOperatorKern{U,V}
+    kernel::U
+    op2::V
+end
+
+function integrand(a::PotentialIntegralOperatorKern,krn,y,f,x)
+    return a.op2(a.kernel(y,x),f.value)
+end
+function potential(op::PotentialIntegralOperator, points, coeffs, basis; 
+    type=SVector{3,ComplexF64},
+	quadstrat=defaultquadstrat(op, basis))
+
+    return potential(PotentialIntegralOperatorKern(op.kernel,op.op2),points,coeffs,op.bfunc(basis);type,quadstrat)
+end
+
+defaultquadstrat(op::PotentialIntegralOperator, basis) = SingleNumQStrat(6)
+quaddata(op::PotentialIntegralOperatorKern,rs,els,qs::SingleNumQStrat) = quadpoints(rs,els,(qs.quad_rule,))
+quadrule(op::PotentialIntegralOperatorKern,refspace,p,y,q,el,qdata,qs::SingleNumQStrat) = qdata[1,q]
+kernelvals(op::PotentialIntegralOperatorKern,y,p) = nothing
\ No newline at end of file
diff --git a/src/composedoperators/trace.jl b/src/composedoperators/trace.jl
new file mode 100644
index 00000000..ad5b1736
--- /dev/null
+++ b/src/composedoperators/trace.jl
@@ -0,0 +1,179 @@
+struct DisplacementVector{T} <: Kernel{T} 
+    interior::Bool
+end
+
+
+function (u::DisplacementVector)(x,disp)
+    if disp == 0
+        return (u.interior - !u.interior) * normal(x)
+    else
+        return sign(disp)*normal(x)
+    end
+end
+
+
+"""
+    _trace(kernel,interior,Val{N})
+
+    computes the jump contribution of a kernel, where N is the dimension of the potential integral.
+"""
+_trace(::Kernel,interior,_) = 0, 0
+function _trace(::HH3DGradGreen{T},interior::Bool,::Val{2}) where {T}
+    return DisplacementVector{real(T)}(interior) , real(T)(0.5)
+end
+
+"""
+    ttrace(Potential,interior::Bool;testfunction_tangential=false)
+
+Compute the tangential trace, -n×n× Potential, of a potential operator mapping it to a boundary operator.
+
+This function assumes the normalvector on the mesh to point outwards. 
+Global multi-trace structures are supported when this function is called.
+A small acceleration can be gained by setting testfunction_tangential to true.
+"""
+function ttrace(a::PotentialIntegralOperator{D},interior::Bool;testfunction_tangential=false) where {D}
+    t,f = _trace(a.kernel,interior,Val{D}())
+    if testfunction_tangential
+        doubleint = CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc)
+        t!=0 && (singleint = f*CompSingleInt(B->B,(x,y)->transpose(x)*y,t,a.op2,a.bfunc))
+    else
+
+        doubleint = -1*CompDoubleInt(B-> n×(n×B),(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc)
+        t!=0 && (singleint = (-f)*CompSingleInt(B -> n×(n×B),(x,y)->transpose(x)*y,t,a.op2,a.bfunc))
+    end
+    if t==0
+     
+        return doubleint
+    else
+       
+        #eventualy a scan to be non zero can be added here.
+        return doubleint + singleint
+    end
+end
+"""
+    ntrace(Potential,interior::Bool)
+
+Compute the normal trace, n⋅ Potential, of a potential operator mapping it to a boundary operator.
+
+This function assumes the normalvector on the mesh to point outwards. 
+Global multi-trace structures are supported when this function is called.
+"""
+function ntrace(a::PotentialIntegralOperator{D},interior::Bool) where {D}
+    t,f = _trace(a.kernel,interior,Val{D}())
+    doubleint = CompDoubleInt(B-> B*n,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc)
+    t!=0 && (singleint = f*CompSingleInt(B-> B*n,(x,y)->transpose(x)*y,t,a.op2,a.bfunc))
+    
+    if t==0
+        return doubleint
+    else
+        #eventualy a scan to be non zero can be added here.
+        return doubleint + singleint
+    end
+end
+"""
+    trace(Potential,interior::Bool)
+
+Compute the trace of a potential operator mapping it to a boundary operator.
+
+This function assumes the normalvector on the mesh to point outwards. 
+Global multi-trace structures are supported when this function is called.
+"""
+function trace(a::PotentialIntegralOperator{D},interior::Bool) where {D}
+    t,f = _trace(a.kernel,interior,Val{D}())
+    doubleint = CompDoubleInt(B->B,*,a.kernel,a.op2,a.bfunc)
+    t!=0 && (singleint = f*CompSingleInt(B->B,*,t,a.op2,a.bfunc))
+    
+    if t==0
+        return doubleint
+    else
+        #eventualy a scan to be non zero can be added here.
+        return doubleint + singleint
+    end
+end
+"""
+    rtrace(Potential,interior::Bool)
+
+Compute the rotated trace, n × Potential, of a potential operator mapping it to a boundary operator.
+
+This function assumes the normalvector on the mesh to point outwards. 
+Global multi-trace structures are supported when this function is called.
+"""
+function rtrace(a::PotentialIntegralOperator{D},interior::Bool) where {D}
+    t,f = _trace(a.kernel,interior,Val{D}())
+    doubleint = CompDoubleInt(B-> B × n,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc)
+    t!=0 && (singleint = f*CompSingleInt(B-> B × n,(x,y)->transpose(x)*y,t,a.op2,a.bfunc))
+    
+    if t==0
+        return doubleint
+    else
+        #eventualy a scan to be non zero can be added here.
+        return doubleint + singleint
+    end
+end
+"""
+    pvttrace(Potential,interior::Bool;testfunction_tangential=false)
+
+Compute the principal value of the tangential trace, -n×n× Potential.
+
+A small acceleration can be gained by setting testfunction_tangential to true.
+"""
+# the principal value of the trace.
+function pvttrace(a::PotentialIntegralOperator{D};testfunction_tangential=false) where {D}
+    if testfunction_tangential
+        doubleint = CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc)
+    else
+        doubleint = -1*CompDoubleInt(B-> n×(n×B),(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc)
+    end
+    return doubleint
+end
+
+"""
+    pvntrace(Potential,interior::Bool)
+
+Compute the principal value of the normal trace, n⋅ Potential.
+
+This function assumes the normalvector on the mesh to point outwards. 
+"""
+function pvntrace(a::PotentialIntegralOperator{D}) where {D}
+
+    doubleint = CompDoubleInt(B-> B*n,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc)
+
+    return doubleint
+end
+
+"""
+    pvtrace(Potential,interior::Bool)
+
+Compute the principal value of the  trace.
+"""
+function pvtrace(a::PotentialIntegralOperator{D}) where {D}
+    doubleint = CompDoubleInt(B-> B,*,a.kernel,a.op2,a.bfunc)
+
+    return doubleint
+end
+
+"""
+    pvrtrace(Potential,interior::Bool)
+
+Compute the principal value of the rotated trace, n× Potential.
+
+This function assumes the normalvector on the mesh to point outwards. 
+"""
+function pvrtrace(a::PotentialIntegralOperator{D}) where {D}
+    doubleint = CompDoubleInt(B-> B × n,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc)
+    return doubleint
+end
+
+"""
+    volumetrace(Potential)
+
+Mapping the potential to a volume operator.
+"""
+function volumetrace(a::PotentialIntegralOperator{D}) where {D}
+    CompDoubleInt(B->B,(x,y)->transpose(x)*y,a.kernel,a.op2,a.bfunc)
+end
+
+
+
+
+
diff --git a/src/identityop.jl b/src/identityop.jl
index c55d1dce..0209bae7 100644
--- a/src/identityop.jl
+++ b/src/identityop.jl
@@ -22,15 +22,15 @@ end
 
 const LinearRefSpaceTriangle = Union{RTRefSpace, RT2RefSpace, NDRefSpace, ND2RefSpace, BDMRefSpace, NCrossBDMRefSpace}
 defaultquadstrat(::LocalOperator, ::LinearRefSpaceTriangle, ::LinearRefSpaceTriangle) = SingleNumQStrat(6)
-function quaddata(op::LocalOperator, g::LinearRefSpaceTriangle, f::LinearRefSpaceTriangle, tels, bels,
-        qs::SingleNumQStrat)
+# function quaddata(op::LocalOperator, g::LinearRefSpaceTriangle, f::LinearRefSpaceTriangle, tels, bels,
+#         qs::SingleNumQStrat)
 
-    u, w = trgauss(qs.quad_rule)
-    qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
-    A = _alloc_workspace(qd, g, f, tels, bels)
+#     u, w = trgauss(qs.quad_rule)
+#     qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
+#     A = _alloc_workspace(qd, g, f, tels, bels)
 
-    return qd, A
-end
+#     return qd, A
+# end
 
 defaultquadstrat(::LocalOperator, ::subReferenceSpace, ::subReferenceSpace) = SingleNumQStrat(6)
 function quaddata(op::LocalOperator, g::subReferenceSpace, f::subReferenceSpace, tels, bels,
@@ -65,14 +65,14 @@ function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,2} where {T,Deg},
 end
 
 defaultquadstrat(::LocalOperator, ::LagrangeRefSpace{T,D1,3}, ::LagrangeRefSpace{T,D2,3}) where {T,D1,D2} = SingleNumQStrat(6)
-function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,3} where {T,Deg},
-    f::LagrangeRefSpace, tels::Vector, bels::Vector, qs::SingleNumQStrat)
+# function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,3} where {T,Deg},
+#     f::LagrangeRefSpace, tels::Vector, bels::Vector, qs::SingleNumQStrat)
 
-    u, w = trgauss(qs.quad_rule)
-    qd = [(w[i], SVector(u[1,i], u[2,i])) for i in 1:length(w)]
-    A = _alloc_workspace(qd, g, f, tels, bels)
-    return qd, A
-end
+#     u, w = trgauss(qs.quad_rule)
+#     qd = [(w[i], SVector(u[1,i], u[2,i])) for i in 1:length(w)]
+#     A = _alloc_workspace(qd, g, f, tels, bels)
+#     return qd, A
+# end
 
 defaultquadstrat(::LocalOperator, ::LagrangeRefSpace{T,D1,4}, ::LagrangeRefSpace{T,D2,4}) where {T,D1,D2} = SingleNumQStrat(6)
 function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,4} where {T,Deg},
@@ -86,20 +86,22 @@ function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,4} where {T,Deg},
      return qd, A
 end
 
+support(::LinearRefSpaceTriangle) = TriangleSupport()
+support(::LinearRefSpaceTetr) = TetraSupport()
 
 defaultquadstrat(::LocalOperator, ::GWPDivRefSpace{<:Real,D1},
     ::GWPDivRefSpace{<:Real,D2}) where {D1,D2} = SingleNumQStrat(7)
 
 
-function quaddata(op::LocalOperator, g::GWPDivRefSpace, f::GWPDivRefSpace,
-    tels::Vector, bels::Vector,
-    qs::SingleNumQStrat)
+# function quaddata(op::LocalOperator, g::GWPDivRefSpace, f::GWPDivRefSpace,
+#     tels::Vector, bels::Vector,
+#     qs::SingleNumQStrat)
 
-    u, w = trgauss(qs.quad_rule)
-    qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
-    A = _alloc_workspace(qd, g, f, tels, bels)
-    return qd, A
-end
+#     u, w = trgauss(qs.quad_rule)
+#     qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
+#     A = _alloc_workspace(qd, g, f, tels, bels)
+#     return qd, A
+# end
 
 
 function quadrule(op::LocalOperator, ψ::RefSpace, ϕ::RefSpace, τ, (qd,A), qs::SingleNumQStrat)
@@ -110,3 +112,38 @@ function quadrule(op::LocalOperator, ψ::RefSpace, ϕ::RefSpace, τ, (qd,A), qs:
     end
     return A
 end
+
+
+defaultquadstrat(::LocalOperator, _, _) = SingleNumQStrat(6)
+# defaultquadstrat(::LocalOperator, ::TetraSupport, ::TetraSupport) = SingleNumQStrat(3)
+# defaultquadstrat(op::LocalOperator, a::_LocalBasisOperations, b) = defaultquadstrat(op,support(a),support(b))
+# defaultquadstrat(op::LocalOperator, a, b::_LocalBasisOperations) = defaultquadstrat(op,support(a),support(b))
+# defaultquadstrat(op::LocalOperator, a::_LocalBasisOperations, b::_LocalBasisOperations) = defaultquadstrat(op,support(a),support(b))
+
+# function quaddata(op::LocalOperator, g, f, tels, bels,
+#     qs::SingleNumQStrat)
+#     gs = support(g)
+#     fs = support(f)
+#     quaddata(op,g,f,gs,fs,tels,bels,qs)
+# end
+
+# function quaddata(op::LocalOperator, g, f,supg::TriangleSupport,supf::TriangleSupport, tels, bels,
+#         qs::SingleNumQStrat)
+
+#     u, w = trgauss(qs.quad_rule)
+#     qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
+#     A = _alloc_workspace(qd, g, f, tels, bels)
+
+#     return qd, A
+# end
+
+function quaddata(op::LocalOperator, g, f, 
+    tels::Vector{<:CompScienceMeshes.Simplex{U,2}}, bels::Vector{<:CompScienceMeshes.Simplex{U,2}},
+    qs::SingleNumQStrat) where {U}
+
+    u, w = trgauss(qs.quad_rule)
+    qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
+    A = _alloc_workspace(qd, g, f, tels, bels)
+
+    return qd, A
+end
\ No newline at end of file
diff --git a/src/quadrature/nonconformingoverlapqrule.jl b/src/quadrature/nonconformingoverlapqrule.jl
index 12c563b0..ec8fc8ff 100644
--- a/src/quadrature/nonconformingoverlapqrule.jl
+++ b/src/quadrature/nonconformingoverlapqrule.jl
@@ -1,7 +1,10 @@
 struct NonConformingOverlapQRule{S}
     conforming_qstrat::S
 end
-
+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 
+end
 
 function momintegrals!(op,
     test_local_space, basis_local_space,
@@ -23,6 +26,9 @@ function momintegrals!(op,
     test_charts = [ch for ch in test_charts if volume(ch) .> 1e6 * eps(T)]
     bsis_charts = [ch for ch in bsis_charts if volume(ch) .> 1e6 * eps(T)]
 
+    test_charts = [ch for ch in test_charts if tangent_rank(ch)]
+    bsis_charts = [ch for ch in bsis_charts if tangent_rank(ch)]
+
     isempty(test_charts) && return
     isempty(bsis_charts) && return
 
diff --git a/src/solvers/solver.jl b/src/solvers/solver.jl
index 92542222..6bbfb93a 100644
--- a/src/solvers/solver.jl
+++ b/src/solvers/solver.jl
@@ -257,8 +257,8 @@ function assemble(bf::BilForm, X::DirectProductSpace, Y::DirectProductSpace;
             y = op[end](op[1:end-1]..., y)
         end
         
-        a = term.coeff * term.kernel
-        z = materialize(a, x, y)
+        a = term.kernel
+        z = term.coeff * materialize(a, x, y)
 
         Smap = lift(z, Block(term.test_id), Block(term.trial_id), U, V)
         T = promote_type(T, eltype(Smap))
diff --git a/test/runtests.jl b/test/runtests.jl
index bbcee67f..9d544756 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -79,6 +79,10 @@ include("test_variational.jl")
 include("test_handlers.jl")
 include("test_gridfunction.jl")
 
+include("test_composed_basis.jl")
+include("test_composed_operator.jl")
+include("test_analytic_excitation.jl")
+
 using TestItemRunner
 @run_package_tests
 
diff --git a/test/test_analytic_excitation.jl b/test/test_analytic_excitation.jl
new file mode 100644
index 00000000..b3d5e3de
--- /dev/null
+++ b/test/test_analytic_excitation.jl
@@ -0,0 +1,41 @@
+using BEAST
+using CompScienceMeshes
+using LinearAlgebra
+using Test
+using StaticArrays
+
+# define the plane wave excitation analytically
+Γ = meshcuboid(1.0,1.0,1.0,0.2)
+Y = unitfunctionc0d1(Γ)
+X = raviartthomas(Γ)
+
+### test crosstracemw
+a = (@SVector [1.0,0.0,0.0])
+eloc(x) = a*exp(-1.0im*dot((@SVector [1.0,0.0,0.0]),x))# not physical but to test the traces.
+ex = BEAST.CrossTraceMW(BEAST.Trace(BEAST.FunctionWrapper{SVector{3,ComplexF64}}(eloc)))
+# probeer testen te verzinnen die integraal nul geven.
+
+rhs = assemble(ex, n× X)
+
+# loop projector
+∂Γ = boundary(Γ)
+edges = setminus(skeleton(Γ,1), ∂Γ)
+verts = setminus(skeleton(Γ,0), skeleton(∂Γ,0))
+Λ = Matrix(connectivity(verts, edges, sign))
+
+PΛ = Λ * pinv(Λ'*Λ) * Λ'
+
+z = PΛ  * rhs
+@test isapprox(norm(z), 0.0;atol = sqrt(eps()))
+
+### test NDotTraceMW
+
+a = (@SVector [1.0,0.0,0.0])
+eloc(x) = a
+ex = BEAST.NDotTrace(BEAST.Trace(BEAST.FunctionWrapper{SVector{3,Float64}}(eloc)))
+# probeer testen te verzinnen die integraal nul geven.
+
+rhs = assemble(ex, Y)
+@test isapprox(norm(rhs), 0.0;atol = sqrt(eps()))
+
+
diff --git a/test/test_basis.jl b/test/test_basis.jl
index 64a7d0ef..8e83f4b7 100644
--- a/test/test_basis.jl
+++ b/test/test_basis.jl
@@ -251,8 +251,8 @@ for i in eachindex(crl.fns)
     crli = sort(crl.fns[i], by=sh->(sh.cellid, sh.refid))
     nxgradi = sort(nxgrad.fns[i], by=sh->(sh.cellid, sh.refid))
     for j in eachindex(crl.fns[i])
-        # @test crl.fns[i][j] == nxgrad.fns[i][j]
-        @test crli[j].coeff == -nxgradi[j].coeff
+        @test crl.fns[i][j] == nxgrad.fns[i][j]
+        #@test crli[j].coeff == -nxgradi[j].coeff
     end
 end
 
diff --git a/test/test_composed_basis.jl b/test/test_composed_basis.jl
new file mode 100644
index 00000000..527e2f3b
--- /dev/null
+++ b/test/test_composed_basis.jl
@@ -0,0 +1,26 @@
+######## test multiplied basis
+using CompScienceMeshes
+using LinearAlgebra
+using BEAST
+using Test
+using StaticArrays
+
+Γ = meshcuboid(1.0,1.0,1.0,1.0)
+X = raviartthomas(Γ)
+func = BEAST.FunctionWrapper{Float64}(x->2.0)
+Y = BEAST._BasisTimes(X,func)
+K = Maxwell3D.doublelayer(wavenumber=1.0)
+m1 = assemble(K,Y,Y)
+m2 = assemble(K,X,X)
+@test m1 ≈ 4*m2
+
+
+L = lagrangecxd0(Γ)
+Z = BEAST._BasisTimes(L,n)
+m3 = assemble(K,Z,Z)
+@test norm(m3) ≈ 0.08420116178577139
+
+m1 = assemble(NCross(),X,X)
+m2 = assemble(Identity(),X×n,X)
+
+@test m1 ≈ m2
diff --git a/test/test_composed_operator.jl b/test/test_composed_operator.jl
new file mode 100644
index 00000000..71f2cec2
--- /dev/null
+++ b/test/test_composed_operator.jl
@@ -0,0 +1,53 @@
+######## test multiplied basis
+using CompScienceMeshes
+using LinearAlgebra
+using BEAST
+using Test
+using StaticArrays
+
+
+# different options to compute selfpatch trace
+Γ = meshcuboid(1.0,1.0,1.0,0.4)
+Γd = BEAST.GlobalDisplacementMesh(Γ,1.0)
+X = raviartthomas(Γ)
+Xd = raviartthomas(Γd)
+
+∇G = BEAST.HH3DGradGreen(0.0)
+Z = BEAST.PotentialIntegralOperator{2}(∇G,×,b->b)
+K = BEAST.ttrace(Z,false;testfunction_tangential=false) #-1.0 trace taken allong the test normal (from outside to inside)
+M = assemble(K,X,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6),BEAST.SingleNumQStrat(6)])
+
+
+
+K_compare = Maxwell3D.doublelayer(wavenumber = 0.0) - 0.5*NCross() #trace from outside to inside
+M_compare = assemble(K_compare,X,X;quadstrat = [BEAST.SingleNumQStrat(6),BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6)])
+
+M_displaced = assemble(K,Xd,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6),BEAST.SingleNumQStrat(6)])
+
+@test M_displaced ≈ M
+@test M_displaced ≈ M_compare 
+
+####################################################
+#
+# Material interactions (global multi-trace context)
+#
+####################################################
+
+Γ_mirror = -mirrormesh(Γ,(@SVector [1.0,0.0,0.0]),(@SVector [0.0,0.0,0.0]))
+Γd_mirror = BEAST.GlobalDisplacementMesh(Γ_mirror,-1.0)
+Xm = raviartthomas(Γ_mirror)
+Xdm = raviartthomas(Γd_mirror)
+
+K = ttrace(Z,true) #number does not matter here, trace part should just be spawned, this number only matters fr selfpatch or if user is not carefull with displacement meshes
+M = assemble(K,Xm,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6),BEAST.SingleNumQStrat(6)])
+
+K_compare = Maxwell3D.doublelayer(wavenumber = 0.0) + 0.5*NCross() #trace from outside to inside
+M_compare = assemble(K_compare,Xm,X;quadstrat = [BEAST.SingleNumQStrat(6),BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6)])
+
+M_displaced = assemble(K,Xdm,X;quadstrat = [BEAST.DoubleNumSauterQstrat(6,6,6,6,6,6),BEAST.SingleNumQStrat(6)])
+
+@test M_displaced ≈ M
+@test M_displaced ≈ M_compare 
+
+
+
diff --git a/test/test_dipole.jl b/test/test_dipole.jl
index f11e277f..f96f287b 100644
--- a/test/test_dipole.jl
+++ b/test/test_dipole.jl
@@ -57,7 +57,7 @@ for U in [Float32,Float64]
     nf_H_EFIE = potential(BEAST.MWDoubleLayerField3D(𝓚), pts, j_EFIE, X)
     ff_E_EFIE = potential(MWFarField3D(𝓣), pts, j_EFIE, X)
     ff_H_EFIE = potential(BEAST.MWDoubleLayerFarField3D(𝓚), pts, j_EFIE, X)
-    ff_H_EFIE_rotated = potential(n × BEAST.MWDoubleLayerFarField3D(𝓚), pts, -j_EFIE, n × X)
+    ff_H_EFIE_rotated = -potential(n × BEAST.MWDoubleLayerFarField3D(𝓚), pts, -j_EFIE, n × X)
     ff_H_EFIE_doublerotated = potential(n × BEAST.MWDoubleLayerRotatedFarField3D(n × 𝓚), pts, -j_EFIE, X)