diff --git a/.github/workflows/format_check.yml b/.github/workflows/format_check.yml
index 935035142..bb667589a 100644
--- a/.github/workflows/format_check.yml
+++ b/.github/workflows/format_check.yml
@@ -26,7 +26,7 @@ jobs:
     if: contains(github.event.pull_request.labels.*.name, 'draft') == false
     strategy:
       matrix:
-        julia-version: [1.10]
+        julia-version: ['1.10']
         julia-arch: [x86]
         os: [ubuntu-latest]
     steps:
diff --git a/src/LazyOperators/lazy_operator.jl b/src/LazyOperators/lazy_operator.jl
index 9a62b388e..30631ed1b 100644
--- a/src/LazyOperators/lazy_operator.jl
+++ b/src/LazyOperators/lazy_operator.jl
@@ -287,6 +287,11 @@ get_operator(a::NullOperator) = nothing
 get_args(a::NullOperator) = (nothing,)
 materialize(a::NullOperator, x) = a
 
+# Base.length(::NullOperator) = 1
+# Base.iterate(a::NullOperator) = (a, nothing)
+# Base.iterate(a::NullOperator, state) = nothing
+Base.map(f, a::NullOperator) = f(a)
+
 function show_lazy_operator(op::NullOperator; level = 1, indent = 4, islast = (true,))
     level == 1 && println("\n---------------")
     print_tree_prefix(level, indent, islast)
diff --git a/src/assembler/assembler.jl b/src/assembler/assembler.jl
index 08ad64a90..f1ebab3ea 100644
--- a/src/assembler/assembler.jl
+++ b/src/assembler/assembler.jl
@@ -283,15 +283,11 @@ function _append_contribution!(X, I, J, U, V, values, elementInfo::CellInfo, dom
     icell = cellindex(elementInfo)
     nU = Val(get_ndofs(U, shape(celltype(elementInfo))))
     nV = Val(get_ndofs(V, shape(celltype(elementInfo))))
-    jdofs = get_dofs(U, icell, nU) # columns correspond to the TrialFunction
-    idofs = get_dofs(V, icell, nV) # lines correspond to the TestFunction
-    _idofs, _jdofs = _cartesian_product(idofs, jdofs)
+    Udofs = get_dofs(U, icell, nU) # columns correspond to the TrialFunction
+    Vdofs = get_dofs(V, icell, nV) # lines correspond to the TestFunction
     unwrapValues = _unwrap_cell_integrate(V, values)
-    append!(I, _idofs)
-    append!(J, _jdofs)
-    for _v in unwrapValues
-        append!(X, _v)
-    end
+    matrixvalues = _pack_bilinear_cell_contribution(unwrapValues, Udofs, Vdofs)
+    _append_bilinear!(I, J, X, Vdofs, Udofs, matrixvalues)
     return nothing
 end
 
@@ -348,6 +344,14 @@ function _append_bilinear!(
     nothing
 end
 
+function _pack_bilinear_cell_contribution(
+    values,
+    col_dofs_U::SVector{NU},
+    row_dofs_V::SVector{NV},
+) where {NU, NV}
+    return SMatrix{NV, NU}(values[j][i] for i in 1:NV, j in 1:NU)
+end
+
 function _pack_bilinear_face_contribution(
     values,
     col_dofs_U_n::SVector{NUn},
diff --git a/src/integration/integration.jl b/src/integration/integration.jl
index a08276aff..a6e56d307 100644
--- a/src/integration/integration.jl
+++ b/src/integration/integration.jl
@@ -101,6 +101,7 @@ end
 end
 _apply_quadrature2(wu) = reduce(_mapsum, wu)
 _mapsum(a, b) = map(+, a, b)
+_mapsum(a::NullOperator, b::NullOperator) = a
 _mapquad(ω, u) = map(Base.Fix1(*, ω), u)
 
 """
diff --git a/test/issues/issue_130.jl b/test/issues/issue_130.jl
new file mode 100644
index 000000000..32df16504
--- /dev/null
+++ b/test/issues/issue_130.jl
@@ -0,0 +1,55 @@
+@testset "issue #130" begin
+    function run()
+        # Settings
+        Zin = 4
+        Zout = 8
+        nr = 3
+        nz = 3
+        degree = 1
+        qdeg = 2 * degree + 1 # quad degree
+
+        # Mesh and measures
+        mesh = rectangle_mesh(
+            nr,
+            nz;
+            ymin = Zin,
+            ymax = Zout,
+            bnd_names = ("AXIS", "WALL", "INLET", "OUTLET"),
+        )
+        Ω = CellDomain(mesh)
+        dΩ = Measure(Ω, qdeg)
+
+        # FESpace
+        fs = FunctionSpace(:Lagrange, degree)
+        U_up = TrialFESpace(fs, mesh)
+        U_ur = TrialFESpace(fs, mesh)
+        U_uz = TrialFESpace(fs, mesh)
+
+        V_up = TestFESpace(U_up)
+        V_ur = TestFESpace(U_ur)
+        V_uz = TestFESpace(U_uz)
+
+        U = MultiFESpace(U_up, U_ur, U_uz)
+        V = MultiFESpace(V_up, V_ur, V_uz)
+
+        # Bilinear forms
+        function test(u, ∇u, v, ∇v)
+            up, ur, uz = u
+            vp, vr, vz = v
+            ∇up, ∇ur, ∇uz = ∇u
+
+            ∂u∂r = ∇ur ⋅ SA[1, 0]
+            return ∂u∂r * vp
+        end
+
+        a((up, ur, uz), (vp, vr, vz)) = ∫(∇(ur) ⋅ SA[1, 0] * vp)dΩ
+        A = assemble_bilinear(a, U, V)
+
+        b(u, v) = ∫(test ∘ (u, map(∇, u), v, map(∇, v)))dΩ
+        B = assemble_bilinear(b, U, V)
+
+        @test A == B
+    end
+
+    run()
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index acec2b2cf..2e0f3da30 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -163,4 +163,8 @@ fname2sum = Dict(r[2] => r[1] for r in eachrow(f))
     custom_include("./operator/test_algebra.jl")
     custom_include("./dof/test_meshdata.jl")
     custom_include("./writers/test_vtk.jl")
+
+    @testset "Issues" begin
+        custom_include("./issues/issue_130.jl")
+    end
 end