diff --git a/extensions/GalerkinToolkitExamples/src/example001.jl b/extensions/GalerkinToolkitExamples/src/example001.jl
index 957a5256..03119a85 100644
--- a/extensions/GalerkinToolkitExamples/src/example001.jl
+++ b/extensions/GalerkinToolkitExamples/src/example001.jl
@@ -6,6 +6,7 @@ using StaticArrays
 using LinearAlgebra
 using SparseArrays
 using WriteVTK
+using PartitionedArrays: JaggedArray, dense_vector
 
 # This one implements a vanilla sequential iso-parametric Poisson solver by only
 # using the mesh interface.
@@ -118,7 +119,7 @@ function setup_isomap(params,integration)
     end
     rid_to_shape_vals = map(first,shape_funs)
     rid_to_shape_grads = map(last,shape_funs)
-    face_to_nodes = gk.face_nodes(mesh,d)
+    face_to_nodes = JaggedArray(gk.face_nodes(mesh,d))
     node_to_coords = gk.node_coordinates(mesh)
     isomap = (;face_to_nodes,node_to_coords,face_to_rid,rid_to_shape_vals,rid_to_shape_grads,d)
 end
@@ -142,6 +143,28 @@ function setup_user_funs(params)
     user_funs = (;u,f,g)
 end
 
+function setup_dofs(params,dirichlet_bcs,cell_isomap,face_isomap)
+    free_and_dirichlet_nodes = dirichlet_bcs.free_and_dirichlet_nodes
+    node_to_free_node = gk.permutation(free_and_dirichlet_nodes)
+    n_free = length(first(free_and_dirichlet_nodes))
+    n_dofs = n_free
+    cell_to_nodes = cell_isomap.face_to_nodes
+    face_to_nodes = face_isomap.face_to_nodes
+    function node_to_dof(node)
+        free_node = node_to_free_node[node]
+        if free_node <= n_free
+            return Int(free_node)
+        end
+        Int(-(free_node-n_free))
+    end
+    cell_to_dofs_data = node_to_dof.(cell_to_nodes.data)
+    face_to_dofs_data = node_to_dof.(face_to_nodes.data)
+    cell_to_dofs = JaggedArray(cell_to_dofs_data,cell_to_nodes.ptrs)
+    face_to_dofs = JaggedArray(face_to_dofs_data,face_to_nodes.ptrs)
+    dofs = (;cell_to_dofs,face_to_dofs,n_dofs)
+    dofs
+end
+
 function setup(params)
     dirichlet_bcs = setup_dirichlet_bcs(params)
     neumann_bcs = setup_neumann_bcs(params)
@@ -149,9 +172,10 @@ function setup(params)
     face_integration = setup_integration(params,:faces)
     cell_isomap = setup_isomap(params,cell_integration)
     face_isomap = setup_isomap(params,face_integration)
+    dofs = setup_dofs(params,dirichlet_bcs,cell_isomap,face_isomap)
     user_funs = setup_user_funs(params)
     solver = params[:solver]
-    state = (;solver,dirichlet_bcs,neumann_bcs,cell_integration,cell_isomap,face_integration,face_isomap,user_funs)
+    state = (;solver,dirichlet_bcs,neumann_bcs,cell_integration,cell_isomap,face_integration,face_isomap,user_funs,dofs)
 end
 
 function add_basic_info(results,params,state)
@@ -167,7 +191,79 @@ function add_basic_info(results,params,state)
     results
 end
 
-function assemble_system_loop(state)
+function assemble_system(state)
+    args = assemble_sytem_coo(state)
+    I,J,V,II,VV = args
+    n_dofs = state.dofs.n_dofs
+    A = sparse(I,J,V,n_dofs,n_dofs)
+    b = dense_vector(II,VV,n_dofs)
+    A,b
+end
+
+function assemble_sytem_coo(state)
+    args = assemble_sybmolic(state)
+    ii_coo = assemble_numeric_cells!(args...,state)
+    assemble_numeric_faces!(args...,state,ii_coo)
+    set_dirichlet_bcs!(args...,state)
+    filter_dirichlet!(args...)
+    args
+end
+
+function assemble_sybmolic(state)
+    cell_to_dofs = state.dofs.cell_to_dofs
+    face_to_dofs = state.dofs.face_to_dofs
+    neum_face_to_face = state.neumann_bcs.neum_face_to_face
+
+    n_coo = 0
+    nn_coo = 0
+    ncells = length(cell_to_dofs)
+    nfaces = length(face_to_dofs)
+    for cell in 1:ncells
+        dofs = cell_to_dofs[cell]
+        ndofs = length(dofs)
+        n_coo += ndofs*ndofs
+        nn_coo += ndofs
+    end
+    for face in neum_face_to_face
+        dofs = face_to_dofs[face]
+        ndofs = length(dofs)
+        nn_coo += ndofs
+    end
+
+    I_coo = zeros(Int,n_coo)
+    J_coo = zeros(Int,n_coo)
+    V_coo = zeros(Float64,n_coo)
+    II_coo = zeros(Int,nn_coo)
+    VV_coo = zeros(Float64,nn_coo)
+
+    n_coo = 0
+    nn_coo = 0
+    for cell in 1:ncells
+        dofs = cell_to_dofs[cell]
+        ndofs = length(dofs)
+        for i in 1:ndofs
+           nn_coo += 1
+           II_coo[nn_coo] = dofs[i]
+            for j in 1:ndofs
+                n_coo += 1
+                I_coo[n_coo] = dofs[i]
+                J_coo[n_coo] = dofs[j]
+            end
+        end
+    end
+    for face in neum_face_to_face
+        dofs = face_to_dofs[face]
+        ndofs = length(dofs)
+        for i in 1:ndofs
+           nn_coo += 1
+           II_coo[nn_coo] = dofs[i]
+       end
+    end
+
+    I_coo,J_coo,V_coo,II_coo,VV_coo
+end
+
+function assemble_numeric_cells!(I_coo,J_coo,V_coo,II_coo,VV_coo,state)
 
     cell_to_nodes = state.cell_isomap.face_to_nodes
     cell_to_rid = state.cell_isomap.face_to_rid
@@ -179,49 +275,22 @@ function assemble_system_loop(state)
     f = state.user_funs.f
     w = state.cell_integration.rid_to_weights
     d = state.cell_integration.d
-
-    # Count coo entries
-    n_coo = 0
     ncells = length(cell_to_nodes)
-    node_to_free_node = gk.permutation(free_and_dirichlet_nodes)
-    nfree = length(first(free_and_dirichlet_nodes))
-    for cell in 1:ncells
-        nodes = cell_to_nodes[cell]
-        for nj in nodes
-            if node_to_free_node[nj] > nfree
-                continue
-            end
-            for ni in nodes
-                if node_to_free_node[ni] > nfree
-                    continue
-                end
-                n_coo += 1
-            end
-        end
-    end
-
-    # Allocate coo values
-    I_coo = zeros(Int32,n_coo)
-    J_coo = zeros(Int32,n_coo)
-    V_coo = zeros(Float64,n_coo)
-    b = zeros(nfree)
 
     # Allocate auxiliary buffers
     ∇x = map(i->similar(i,size(i,2)),∇s)
     Aes = map(i->zeros(size(i,2),size(i,2)),∇s)
     fes = map(i->zeros(size(i,2)),∇s)
-    ues = map(i->zeros(size(i,2)),∇s)
 
     Tx = SVector{d,Float64}
     TJ = typeof(zero(Tx)*zero(Tx)')
 
-    # Fill coo values
     i_coo = 0
+    ii_coo = 0
     for cell in 1:ncells
         rid = cell_to_rid[cell]
         nodes = cell_to_nodes[cell]
         Ae = Aes[rid]
-        ue = ues[rid]
         fe = fes[rid]
         nl = length(nodes)
         ∇se = ∇s[rid]
@@ -231,15 +300,7 @@ function assemble_system_loop(state)
         nq = length(we)
         fill!(Ae,zero(eltype(Ae)))
         fill!(fe,zero(eltype(Ae)))
-        fill!(ue,zero(eltype(Ae)))
-        for k in 1:nl
-            nk = nodes[k]
-            gk = node_to_free_node[nk]
-            if  gk <= nfree
-                continue
-            end
-            ue[k] = u_dirichlet[gk-nfree]
-        end
+        # Integrate
         for iq in 1:nq
             Jt = zero(TJ) 
             xint = zero(Tx)
@@ -266,54 +327,34 @@ function assemble_system_loop(state)
                 fe[k] +=  fx*se[iq,k]*dV
             end
         end
+        # Set the result in the output array
         for i in 1:nl
-            for j in 1:nl
-                fe[i] -= Ae[i,j]*ue[j]
-            end
-        end
-        for i in 1:nl
-            ni = nodes[i]
-            gi = node_to_free_node[ni]
-            if gi > nfree
-                continue
-            end
-            b[gi] += fe[i]
+            ii_coo += 1
+            VV_coo[ii_coo] = fe[i]
         end
         for j in 1:nl
-            nj = nodes[j]
-            gj = node_to_free_node[nj]
-            if  gj > nfree
-                continue
-            end
             for i in 1:nl
-                ni = nodes[i]
-                gi = node_to_free_node[ni]
-                if gi > nfree
-                    continue
-                end
                 i_coo += 1
-                I_coo[i_coo] = gi
-                J_coo[i_coo] = gj
                 V_coo[i_coo] = Ae[i,j]
             end
         end
     end
-    add_neumann_bcs!(b,state)
-    I_coo,J_coo,V_coo,b
+
+    ii_coo
 end
 
-function add_neumann_bcs!(b,state)
+function assemble_numeric_faces!(I_coo,J_coo,V_coo,II_coo,VV_coo,state,ii_coo)
 
     neum_face_to_face = state.neumann_bcs.neum_face_to_face
     if length(neum_face_to_face) == 0
         return nothing
     end
+
     face_to_rid = state.face_isomap.face_to_rid
     face_to_nodes = state.face_isomap.face_to_nodes
     ∇s_f = state.face_isomap.rid_to_shape_grads
     s_f = state.face_isomap.rid_to_shape_vals
     node_to_x = state.face_isomap.node_to_coords
-    free_and_dirichlet_nodes = state.dirichlet_bcs.free_and_dirichlet_nodes
     g = state.user_funs.g
     w_f = state.face_integration.rid_to_weights
     d = state.cell_integration.d
@@ -321,8 +362,6 @@ function add_neumann_bcs!(b,state)
     fes = map(i->zeros(size(i,2)),s_f)
     Tx = SVector{d,Float64}
     Ts = typeof(zero(SVector{d-1,Float64})*zero(Tx)')
-    node_to_free_node = gk.permutation(free_and_dirichlet_nodes)
-    nfree = length(first(free_and_dirichlet_nodes))
 
     for face in neum_face_to_face
         rid = face_to_rid[face]
@@ -352,23 +391,63 @@ function add_neumann_bcs!(b,state)
             end
         end
         for i in 1:nl
-            ni = nodes[i]
-            gi = node_to_free_node[ni]
-            if gi > nfree
+            ii_coo += 1
+            VV_coo[ii_coo] = fe[i]
+        end
+    end
+end
+
+function set_dirichlet_bcs!(I_coo,J_coo,V_coo,II_coo,VV_coo,state)
+
+    cell_to_dofs = state.dofs.cell_to_dofs
+    ncells = length(cell_to_dofs)
+    u_dirichlet = state.dirichlet_bcs.u_dirichlet
+    i_coo = 0
+    ii_coo = 0
+    for cell in 1:ncells
+        dofs = cell_to_dofs[cell]
+        ndofs = length(dofs)
+        for i in 1:ndofs
+            if !(dofs[i]>0)
                 continue
             end
-            b[gi] += fe[i]
+            for j in 1:ndofs
+                if dofs[j]>0
+                    continue
+                end
+                uj = u_dirichlet[-dofs[j]]
+                ij = (j-1)*ndofs + i
+                Aij = V_coo[ij+i_coo]
+                VV_coo[i+ii_coo] -= Aij*uj
+            end
         end
+        i_coo += ndofs*ndofs
+        ii_coo += ndofs
     end
 
-    nothing
 end
 
-function assemble_system(state)
-    I,J,V,b = assemble_system_loop(state)
-    nfree = length(b)
-    A = sparse(I,J,V,nfree,nfree)
-    A,b
+function filter_dirichlet!(I_coo,J_coo,V_coo,II_coo,VV_coo)
+    # TODO this is only needed when the functions in charge of
+    # compressing the matrix/vector are not able to skip negative
+    # values
+    # NB. It will not work for empty problems
+    for p in 1:length(I_coo)
+        i = I_coo[p]
+        j = J_coo[p]
+        if i<=0 || j<=0
+            I_coo[p] = 1
+            J_coo[p] = 1
+            V_coo[p] = 0
+        end
+    end
+    for p in 1:length(II_coo)
+        i = II_coo[p]
+        if i<=0
+            II_coo[p] = 1
+            VV_coo[p] = 0
+        end
+    end
 end
 
 function solve_system(A,b,params)
@@ -481,4 +560,211 @@ function export_results(uh,params,state)
     end
 end
 
+### From here, it is old stuff to to be deleted
+
+function assemble_system_loop(state)
+
+    cell_to_nodes = state.cell_isomap.face_to_nodes
+    cell_to_rid = state.cell_isomap.face_to_rid
+    free_and_dirichlet_nodes = state.dirichlet_bcs.free_and_dirichlet_nodes
+    node_to_x = state.cell_isomap.node_to_coords
+    ∇s = state.cell_isomap.rid_to_shape_grads
+    s = state.cell_isomap.rid_to_shape_vals
+    u_dirichlet = state.dirichlet_bcs.u_dirichlet
+    f = state.user_funs.f
+    w = state.cell_integration.rid_to_weights
+    d = state.cell_integration.d
+
+    # Count coo entries
+    n_coo = 0
+    ncells = length(cell_to_nodes)
+    node_to_free_node = gk.permutation(free_and_dirichlet_nodes)
+    nfree = length(first(free_and_dirichlet_nodes))
+    for cell in 1:ncells
+        nodes = cell_to_nodes[cell]
+        for nj in nodes
+            if node_to_free_node[nj] > nfree
+                continue
+            end
+            for ni in nodes
+                if node_to_free_node[ni] > nfree
+                    continue
+                end
+                n_coo += 1
+            end
+        end
+    end
+
+    # Allocate coo values
+    I_coo = zeros(Int32,n_coo)
+    J_coo = zeros(Int32,n_coo)
+    V_coo = zeros(Float64,n_coo)
+    b = zeros(nfree)
+
+    # Allocate auxiliary buffers
+    ∇x = map(i->similar(i,size(i,2)),∇s)
+    Aes = map(i->zeros(size(i,2),size(i,2)),∇s)
+    fes = map(i->zeros(size(i,2)),∇s)
+    ues = map(i->zeros(size(i,2)),∇s)
+
+    Tx = SVector{d,Float64}
+    TJ = typeof(zero(Tx)*zero(Tx)')
+
+    # Fill coo values
+    i_coo = 0
+    for cell in 1:ncells
+        rid = cell_to_rid[cell]
+        nodes = cell_to_nodes[cell]
+        Ae = Aes[rid]
+        ue = ues[rid]
+        fe = fes[rid]
+        nl = length(nodes)
+        ∇se = ∇s[rid]
+        ∇xe = ∇x[rid]
+        se = s[rid]
+        we = w[rid]
+        nq = length(we)
+        fill!(Ae,zero(eltype(Ae)))
+        fill!(fe,zero(eltype(Ae)))
+        fill!(ue,zero(eltype(Ae)))
+        for k in 1:nl
+            nk = nodes[k]
+            gk = node_to_free_node[nk]
+            if  gk <= nfree
+                continue
+            end
+            ue[k] = u_dirichlet[gk-nfree]
+        end
+        for iq in 1:nq
+            Jt = zero(TJ) 
+            xint = zero(Tx)
+            for k in 1:nl
+                x = node_to_x[nodes[k]]
+                ∇sqx = ∇se[iq,k]
+                sqx = se[iq,k]
+                Jt += ∇sqx*x'
+                xint += sqx*x
+            end
+            detJt = det(Jt)
+            invJt = inv(Jt)
+            dV = abs(detJt)*we[iq]
+            for k in 1:nl
+                ∇xe[k] = invJt*∇se[iq,k]
+            end
+            for j in 1:nl
+                for i in 1:nl
+                    Ae[i,j] += ∇xe[i]⋅∇xe[j]*dV
+                end
+            end
+            fx = f(xint)
+            for k in 1:nl
+                fe[k] +=  fx*se[iq,k]*dV
+            end
+        end
+        for i in 1:nl
+            for j in 1:nl
+                fe[i] -= Ae[i,j]*ue[j]
+            end
+        end
+        for i in 1:nl
+            ni = nodes[i]
+            gi = node_to_free_node[ni]
+            if gi > nfree
+                continue
+            end
+            b[gi] += fe[i]
+        end
+        for j in 1:nl
+            nj = nodes[j]
+            gj = node_to_free_node[nj]
+            if  gj > nfree
+                continue
+            end
+            for i in 1:nl
+                ni = nodes[i]
+                gi = node_to_free_node[ni]
+                if gi > nfree
+                    continue
+                end
+                i_coo += 1
+                I_coo[i_coo] = gi
+                J_coo[i_coo] = gj
+                V_coo[i_coo] = Ae[i,j]
+            end
+        end
+    end
+    add_neumann_bcs!(b,state)
+    I_coo,J_coo,V_coo,b
+end
+
+function add_neumann_bcs!(b,state)
+
+    neum_face_to_face = state.neumann_bcs.neum_face_to_face
+    if length(neum_face_to_face) == 0
+        return nothing
+    end
+    face_to_rid = state.face_isomap.face_to_rid
+    face_to_nodes = state.face_isomap.face_to_nodes
+    ∇s_f = state.face_isomap.rid_to_shape_grads
+    s_f = state.face_isomap.rid_to_shape_vals
+    node_to_x = state.face_isomap.node_to_coords
+    free_and_dirichlet_nodes = state.dirichlet_bcs.free_and_dirichlet_nodes
+    g = state.user_funs.g
+    w_f = state.face_integration.rid_to_weights
+    d = state.cell_integration.d
+
+    fes = map(i->zeros(size(i,2)),s_f)
+    Tx = SVector{d,Float64}
+    Ts = typeof(zero(SVector{d-1,Float64})*zero(Tx)')
+    node_to_free_node = gk.permutation(free_and_dirichlet_nodes)
+    nfree = length(first(free_and_dirichlet_nodes))
+
+    for face in neum_face_to_face
+        rid = face_to_rid[face]
+        nodes = face_to_nodes[face]
+        ∇se = ∇s_f[rid]
+        se = s_f[rid]
+        we = w_f[rid]
+        nq = length(we)
+        fe = fes[rid]
+        fill!(fe,zero(eltype(fe)))
+        nl = length(nodes)
+        for iq in 1:nq
+            Jt = zero(Ts) 
+            xint = zero(Tx)
+            for k in 1:nl
+                x = node_to_x[nodes[k]]
+                ∇sqx = ∇se[iq,k]
+                sqx = se[iq,k]
+                Jt += ∇sqx*x'
+                xint += sqx*x
+            end
+            J = transpose(Jt)
+            dS = sqrt(det(Jt*J))*we[iq]
+            gx = g(xint)
+            for k in 1:nl
+                fe[k] +=  gx*se[iq,k]*dS
+            end
+        end
+        for i in 1:nl
+            ni = nodes[i]
+            gi = node_to_free_node[ni]
+            if gi > nfree
+                continue
+            end
+            b[gi] += fe[i]
+        end
+    end
+
+    nothing
+end
+
+function assemble_system_old(state)
+    I,J,V,b = assemble_system_loop(state)
+    nfree = length(b)
+    A = sparse(I,J,V,nfree,nfree)
+    A,b
+end
+
+
 end # module
diff --git a/extensions/GalerkinToolkitExamples/test/example001_tests.jl b/extensions/GalerkinToolkitExamples/test/example001_tests.jl
index ceda93e7..5751fb7b 100644
--- a/extensions/GalerkinToolkitExamples/test/example001_tests.jl
+++ b/extensions/GalerkinToolkitExamples/test/example001_tests.jl
@@ -5,8 +5,15 @@ using GalerkinToolkitExamples: Example001
 using Test
 
 tol = 1.0e-10
+
+params = Dict{Symbol,Any}()
+params[:mesh] = gk.cartesian_mesh((0,1,0,1),(3,3))
+results = Example001.main(params)
+@test results[:eh1] < tol
+@test results[:el2] < tol
+
 params = Dict{Symbol,Any}()
-params[:mesh] = gk.cartesian_mesh((0,3,0,2),(10,5),complexify=false)
+params[:mesh] = gk.cartesian_mesh((0,3,0,2),(5,10),complexify=false)
 params[:dirichlet_tags] = ["1-face-1","1-face-3","1-face-4"]
 params[:neumann_tags] = ["1-face-2"]
 params[:u] = (x) -> sum(x)