diff --git a/src/FESpaces/ZeroMeanFESpaces.jl b/src/FESpaces/ZeroMeanFESpaces.jl
index 2549adead..12e46172b 100644
--- a/src/FESpaces/ZeroMeanFESpaces.jl
+++ b/src/FESpaces/ZeroMeanFESpaces.jl
@@ -35,7 +35,8 @@ function FEFunction(
     dirichlet_values,
     f.vol_i,
     f.vol,
-    f.space.dof_to_fix)
+    f.space.dof_to_fix
+  )
   fv = lazy_map(+,free_values,Fill(c,length(free_values)))
   dv = dirichlet_values .+ c
   FEFunction(f.space,fv,dv)
@@ -57,12 +58,12 @@ function _compute_new_fixedval(fv,dv,vol_i,vol,fixed_dof)
     c += fv[i-1]*vol_i[i]
   end
   c = -c/vol
-  c
+  return c
 end
 
 # This is required, otherwise we end up calling `FEFunction` with a fixed value of zero, 
 # which does not properly interpolate the function provided. 
-# With this change, we modify are interpolating in the unconstrained space and then
+# With this change, we are interpolating in the unconstrained space and then
 # substracting the mean.
 function interpolate!(object,free_values,fs::ZeroMeanFESpace)
   dirichlet_values = zero_dirichlet_values(fs)
diff --git a/test/FESpacesTests/ZeroMeanFESpacesTests.jl b/test/FESpacesTests/ZeroMeanFESpacesTests.jl
index 4f98ecc5f..e37a3266f 100644
--- a/test/FESpacesTests/ZeroMeanFESpacesTests.jl
+++ b/test/FESpacesTests/ZeroMeanFESpacesTests.jl
@@ -1,6 +1,7 @@
 module ZeroMeanFESpacesTests
 
 using Test
+using Gridap
 using Gridap.Arrays
 using Gridap.Geometry
 using Gridap.Fields
@@ -45,7 +46,44 @@ eh = vh - g
 @test abs(sum(∫(vh)*dΩ)) < 1.0e-10
 @test sum(∫(eh*eh)*dΩ) < 1.0e-10
 
-V = FESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:L2,constraint=:zeromean)
-@test isa(V,ZeroMeanFESpace)
+# Check Stokes problem properties
+
+u_ex(x) = VectorValue(x[2],-x[1])
+p_ex(x) = x[1] + 2*x[2]
+p_mean = sum(∫(p_ex)*dΩ)
+
+order = 2
+reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
+reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P)
+
+V  = FESpace(model,reffe_u;dirichlet_tags="boundary")
+U  = TrialFESpace(V,u_ex)
+Q  = FESpace(model,reffe_p;conformity=:L2)
+Q0 = FESpace(model,reffe_p;conformity=:L2,constraint=:zeromean)
+@test isa(Q0,ZeroMeanFESpace)
+
+X = MultiFieldFESpace([U,Q0])
+Y = MultiFieldFESpace([V,Q0])
+
+ph_i = interpolate(p_ex, Q0)
+@test abs(sum(∫(ph_i)*dΩ)) < 1.0e-10
+@test abs(sum(∫(ph_i - p_ex + p_mean)*dΩ)) < 1.0e-10
+
+a((u,p),(v,q)) = ∫(∇(u)⊙∇(v) - (∇⋅v)*p - q*(∇⋅u))dΩ
+l((v,q)) = a((u_ex,p_ex),(v,q))
+
+op = AffineFEOperator(a,l,X,Y)
+uh, ph = solve(op)
+
+l2_error(u,v) = sqrt(sum(∫((u-v)⋅(u-v))*dΩ))
+@test l2_error(uh,u_ex) < 1.0e-10
+@test abs(sum(∫(ph)*dΩ)) < 1.0e-10
+@test l2_error(ph,ph_i) < 1.0e-10
+
+b(u,q) = ∫(q*(∇⋅u))dΩ
+B = assemble_vector(q -> b(uh,q),Q)
+B0 = assemble_vector(q -> b(u_ex,q),Q0)
+@test abs(sum(B)) < 1.0e-10
+@test abs(sum(B0)) < 1.0e-10
 
 end # module