diff --git a/src/bases/divergence.jl b/src/bases/divergence.jl
index a4b32add..f66d4b43 100644
--- a/src/bases/divergence.jl
+++ b/src/bases/divergence.jl
@@ -48,6 +48,11 @@ function curl(x::Space)
     curl(x, geo, crl)
 end
 
+function gradient(X::DirectProductSpace{T}) where T
+    x = Space{T}[gradient(s) for s in X.factors]
+    DirectProductSpace(x)
+end
+
 
 function gradient(x::Space)
     ref = refspace(x)
diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl
index ef26ac56..729a8ab1 100644
--- a/src/bases/lagrange.jl
+++ b/src/bases/lagrange.jl
@@ -744,15 +744,15 @@ function duallagrangec0d1(mesh, mesh2, pred, ::Type{Val{2}})
   geometry = mesh2
   cellids2, ncells2 = vertextocellmap(mesh2)
 
-  fns = Vector{Vector{Shape{T}}}(undef,num_cells1)
+  fns = Vector{Vector{Shape{T}}}()
   pos = Vector{vertextype(mesh)}()
   # We will iterate over the coarse mesh segments to assign all the functions to it.
   for segment_coarse in 1 : num_cells1
     # For the dual Lagrange there is a 6 shapes per segment
-    numshapes = (ncells1[segment_coarse]*4) -2
+    numshapes = (ncells1[mesh.faces[segment_coarse].indices[1]]*4) -2
     shapes = Vector{Shape{T}}(undef,numshapes)
     # Now we will get all the smaller faces within the coarse segment
-    #i.e The coose segment will have two points, and these tow points are connected to two segmesnts in the finer mesh
+    #i.e The coose segment will have two points, and these tow points are connected to two segments in the finer mesh
     # This will give us a 4 smaller faces per Dual lagrange basis, we store them first in all_faces
     # all_faces= Array{SVector{2,Int}}(undef,4)                      # faces in the original segment (4)
     CT = CompScienceMeshes.SimplexGraph{2}
@@ -799,7 +799,7 @@ function duallagrangec0d1(mesh, mesh2, pred, ::Type{Val{2}})
         shapes[6]= Shape(cellids2[mesh.faces[segment_coarse][2],1],2,0.5)
     end
     # Now assign all of these shapes to the relevent segment in the coarse mesh
-    fns[segment_coarse]=shapes
+    push!(fns, shapes)
     push!(pos, cartesian(center(chart(mesh, segment_coarse))))
   end
 
diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl
index 2b140c88..7a6f4e5a 100644
--- a/src/bases/local/laglocal.jl
+++ b/src/bases/local/laglocal.jl
@@ -211,7 +211,19 @@ function (f::LagrangeRefSpace{T,10,2})(mp) where T
 end
 
 # Evaluete linear lagrange elements on a triangle
-function (f::LagrangeRefSpace{T,1,3})(t) where T
+function (f::LagrangeRefSpace{T,1,3})(t::CompScienceMeshes.MeshPointNM{T,C,D,2}) where {T,C,D}
+    u,v,w, = barycentric(t)
+
+    j = jacobian(t)
+    p = t.patch
+    SVector(
+        (value=u, curl=(p[3]-p[2])/j),
+        (value=v, curl=(p[1]-p[3])/j),
+        (value=w, curl=(p[2]-p[1])/j))
+end
+
+# Evaluete linear lagrange elements on a triangle
+function (f::LagrangeRefSpace{T,1,3})(t::CompScienceMeshes.MeshPointNM{T,C,D,3}) where {T,C,D}
     u,v,w, = barycentric(t)
 
     j = jacobian(t)
diff --git a/src/bases/local/ndlocal.jl b/src/bases/local/ndlocal.jl
index 3b3919f0..9cd96afc 100644
--- a/src/bases/local/ndlocal.jl
+++ b/src/bases/local/ndlocal.jl
@@ -12,8 +12,27 @@ function valuetype(ref::NDRefSpace{T}, charttype::Type) where {T}
     SVector{universedimension(charttype),T}
 end
 
+function (ϕ::NDRefSpace)(nbd::CompScienceMeshes.MeshPointNM{T,C,D,2}) where {T,C,D}
 
-function (ϕ::NDRefSpace)(nbd)
+    u, v = parametric(nbd)
+    j = jacobian(nbd)
+
+    tu = tangents(nbd,1)
+    tv = tangents(nbd,2)
+
+    d = 2/j
+
+    rotate_left(t) = StaticArrays.SVector{2, Float64}(t[2], -t[1])
+
+    return SVector((
+        (value = rotate_left(tu*(u-1) + tv*v    ) / j, curl = d),
+        (value = rotate_left(tu*u     + tv*(v-1)) / j, curl = d),
+        (value = rotate_left(tu*u     + tv*v    ) / j, curl = d)
+    ))
+
+end
+
+function (ϕ::NDRefSpace)(nbd::CompScienceMeshes.MeshPointNM{T,C,D,3}) where {T,C,D}
 
     u, v = parametric(nbd)
     n = normal(nbd)
diff --git a/src/excitation.jl b/src/excitation.jl
index 3c7d2c52..f3fc8982 100644
--- a/src/excitation.jl
+++ b/src/excitation.jl
@@ -39,7 +39,9 @@ function assemble!(field::Functional, tfs::Space, store;
     quadstrat=defaultquadstrat(field, tfs))
 
     qs = quadstrat(field, tfs)
-    tels, tad = assemblydata(tfs)
+    
+    tr = assemblydata(tfs); tr == nothing && return
+    tels, tad = tr
 
     trefs = refspace(tfs)
     qd = quaddata(field, trefs, tels, qs)
diff --git a/src/localop.jl b/src/localop.jl
index 8b3ad69a..ada78f1a 100644
--- a/src/localop.jl
+++ b/src/localop.jl
@@ -300,8 +300,11 @@ function assemble_local_mixed!(biop::LocalOperator, tfs::Space{T}, bfs::Space{T}
     num_trefs = numfunctions(trefs, tdom)
     num_brefs = numfunctions(brefs, bdom)
 
-    tels, tad = assemblydata(tfs)
-    bels, bad = assemblydata(bfs)
+    tr = assemblydata(tfs); tr == nothing && return
+    br = assemblydata(bfs); br == nothing && return
+
+    tels, tad = tr
+    bels, bad = br
 
     qd = quaddata(biop, trefs, brefs, tels, bels, quadstrat)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 03887f44..2c723d0b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -41,6 +41,7 @@ include("test_interpolate_and_restrict.jl")
 include("test_rt3d.jl")
 include("test_gradient.jl")
 include("test_mult.jl")
+include("test_assemble_2D_surface.jl")
 
 include("test_gram.jl")
 include("test_vector_gram.jl")
diff --git a/test/test_assemble_2D_surface.jl b/test/test_assemble_2D_surface.jl
new file mode 100644
index 00000000..94be3bc4
--- /dev/null
+++ b/test/test_assemble_2D_surface.jl
@@ -0,0 +1,110 @@
+using Test
+
+import Permutations
+import StaticArrays
+using SparseArrays
+using CompScienceMeshes
+using BEAST
+
+"""
+Test if the assemble methode for a 2D mesh with udim=2 works. 
+Test are constructing a rectangle with udim=2 and udim=3 with langrange and nedelec interpolation. 
+And assemble these two end checking that matrix from rectangle of udim=2 is equal to the matrix from the ractangle of udim=3.
+"""
+
+
+
+# overwrite the existing method of CompScienceMeshes for the positive value of the volume.
+function CompScienceMeshes._normals(tangents::StaticArrays.SVector{2,StaticArrays.SVector{2,T}}, ::Type{Val{0}}) where {T}
+
+    t = tangents[1]
+    s = tangents[2]
+    v = abs(t[1]*s[2] - t[2]*s[1])/2
+    # n[3] = tangents[1] × tangents[2]
+    # l = norm(n)
+
+    P = StaticArrays.SVector{2,T}
+    StaticArrays.SVector{0,P}(), v
+end
+
+
+function permutate_vector(X, Y)
+    tol = sqrt(eps(eltype(X[1])))
+
+    permut = Vector{Int32}()  
+    temp = collect(1:length(X))
+
+    for p in Y
+        index = findfirst(isapprox(p;atol = tol), X)
+        @assert !isnothing(index)   # vertex from Y exist in X
+        @assert temp[index] != 0    # vertex from Y unique in X
+
+        temp[index] = 0
+        append!(permut, index)
+    end
+
+    for i in temp
+        if i !=0
+            push!(permut, i)
+    end end
+
+    return permut
+end
+
+
+
+
+
+
+@testset "assemble 2D langrange" begin
+    h = 1/5
+    Ω2D = meshrectangle(1.0, 1.0, h, 2) # udim = 2
+    Ω3D = meshrectangle(1.0, 1.0, h, 3) # udim = 3
+
+    #permutate_mesh(Ω2D, Ω3D)
+
+    # lagrange spaces
+    X2D = lagrangec0d1(Ω2D, skeleton(Ω2D,0))
+    X3D = lagrangec0d1(Ω3D, skeleton(Ω3D,0))
+
+    σ =  permutate_vector(X2D.pos, map(vec -> StaticArrays.SVector(vec[1], vec[2]), X3D.pos))
+
+    @test map(vec -> StaticArrays.SVector(vec[1], vec[2]), X3D.pos) == X2D.pos[σ]
+
+    #@test X2D.geo.vertices == map(x-> x[1:2], X3D.geo.vertices)
+
+    I = BEAST.Identity()
+
+    A2D = assemble(I, X2D, X2D)
+    A3D = assemble(I, X3D, X3D)
+
+    @test sum(abs.(permute(sparse(A2D), σ, σ) - sparse(A3D))) ≈ 0 atol=1e-8
+
+end
+
+@testset "assemble 2D nedelec" begin
+    h = 1/5
+    Ω2D = meshrectangle(1.0, 1.0, h, 2) # udim = 2
+    Ω3D = meshrectangle(1.0, 1.0, h, 3) # udim = 3
+
+    #permutate_mesh(Ω2D, Ω3D)
+
+    # lagrange spaces
+    X2D = lagrangec0d1(Ω2D, skeleton(Ω2D,0))
+    X3D = lagrangec0d1(Ω3D, skeleton(Ω3D,0))
+
+    σ =  permutate_vector(X2D.pos, map(vec -> StaticArrays.SVector(vec[1], vec[2]), X3D.pos))
+
+    @test map(vec -> StaticArrays.SVector(vec[1], vec[2]), X3D.pos) == X2D.pos[σ]
+
+    @hilbertspace u
+    @hilbertspace v
+
+    I = BEAST.Identity()
+    IGG = @varform I[gradient(u), gradient(v)]
+
+    A2D = assemble(IGG, X2D, X2D)
+    A3D = assemble(IGG, X3D, X3D)
+
+    @test sum(abs.(permute(sparse(A2D), σ, σ) - sparse(A3D))) ≈ 0 atol=1e-8
+end
\ No newline at end of file