diff --git a/src/solver/operator.jl b/src/solver/operator.jl
index 770fe26a..0972934d 100644
--- a/src/solver/operator.jl
+++ b/src/solver/operator.jl
@@ -29,7 +29,7 @@ end
 
 # TODO constructor which checks for axis compat
 struct BlockOperator{OPS <: Tuple}
-    # TODO maybe SMatrix?
+    # TODO custom "square matrix tuple"
     operators::OPS # stored row by row as in [1 2; 3 4]
 end
 
@@ -324,6 +324,7 @@ end
         end
         # Evaluate f
         fx = f.f(x,time)
+        # TODO replace with evaluate_coefficient
         # Evaluate all basis functions
         @inbounds for j ∈ 1:getnbasefunctions(cv)
             δu = shape_value(cv, qp, j)
diff --git a/src/utils.jl b/src/utils.jl
index 3dfd9a2f..7c0e571d 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -146,6 +146,8 @@ function mul(A::ThreadedSparseMatrixCSR, x::AbstractVector)
     return mul!(y, A, x)
 end
 *(A::ThreadedSparseMatrixCSR, v::AbstractVector) = mul(A,v)
+*(A::ThreadedSparseMatrixCSR, v::BlockArrays.FillArrays.AbstractZeros{<:Any, 1}) = mul(A,v)
+*(A::ThreadedSparseMatrixCSR, v::BlockArrays.ArrayLayouts.LayoutVector) = mul(A,v)
 
 Base.eltype(A::ThreadedSparseMatrixCSR) = Base.eltype(A.A)
 getrowptr(A::ThreadedSparseMatrixCSR) = getrowptr(A.A)