diff --git a/Manifest.toml b/Manifest.toml
index bb44ef8e..11df4bd3 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -10,38 +10,38 @@ uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 
 [[ArrayInterface]]
 deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
-git-tree-sha1 = "b8d49c34c3da35f220e7295659cd0bab8e739fed"
+git-tree-sha1 = "1d6835607e9f214cb4210310868f8cf07eb0facc"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "3.1.33"
+version = "3.1.34"
 
 [[ArrayLayouts]]
 deps = ["FillArrays", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "a745df3fdce8aff885e6f4a3b7427cc27ce5f86c"
+git-tree-sha1 = "7a92ea1dd16472d18ca1ffcbb7b3cc67d7e78a3f"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
-version = "0.7.6"
+version = "0.7.7"
 
 [[Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
 [[BSON]]
-git-tree-sha1 = "92b8a8479128367aaab2620b8e73dff632f5ae69"
+git-tree-sha1 = "ebcd6e22d69f21249b7b8668351ebf42d6dc87a1"
 uuid = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
-version = "0.3.3"
+version = "0.3.4"
 
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
 [[BlockArrays]]
 deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "782362509cf50a51092f513e92783c18ab0b6d51"
+git-tree-sha1 = "5524e27323cf4c4505699c3fb008c3f772269945"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-version = "0.16.8"
+version = "0.16.9"
 
 [[ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "a325370b9dd0e6bf5656a6f1a7ae80755f8ccc46"
+git-tree-sha1 = "d9e40e3e370ee56c5b57e0db651d8f92bce98fea"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.7.2"
+version = "1.10.1"
 
 [[CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
@@ -130,9 +130,9 @@ version = "1.11.1"
 
 [[FillArrays]]
 deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
-git-tree-sha1 = "29890dfbc427afa59598b8cfcc10034719bd7744"
+git-tree-sha1 = "8756f9935b7ccc9064c6eef0bff0ad643df733a3"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "0.12.6"
+version = "0.12.7"
 
 [[FiniteDiff]]
 deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
@@ -141,16 +141,16 @@ uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
 version = "2.8.1"
 
 [[ForwardDiff]]
-deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
-git-tree-sha1 = "c4203b60d37059462af370c4f3108fb5d155ff13"
+deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
+git-tree-sha1 = "63777916efbcb0ab6173d09a658fb7f2783de485"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.20"
+version = "0.10.21"
 
 [[Gridap]]
 deps = ["AbstractTrees", "BSON", "BlockArrays", "Combinatorics", "DocStringExtensions", "FastGaussQuadrature", "FileIO", "FillArrays", "ForwardDiff", "JLD2", "JSON", "LineSearches", "LinearAlgebra", "NLsolve", "NearestNeighbors", "QuadGK", "Random", "SparseArrays", "SparseMatricesCSR", "StaticArrays", "Test", "WriteVTK"]
 git-tree-sha1 = "d3bfeb769259b08f7385bdad91681f38c50fc519"
 repo-rev = "gridap_distributed"
-repo-url = "https://github.com/gridap/Gridap.jl.git"
+repo-url = "https://github.com/gridap/Gridap.jl"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 version = "0.17.0"
 
@@ -163,10 +163,16 @@ version = "0.1.0"
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[InverseFunctions]]
+deps = ["Test"]
+git-tree-sha1 = "f0c6489b12d28fb4c2103073ec7452f3423bd308"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.1"
+
 [[IrrationalConstants]]
-git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
+git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
 uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
-version = "0.1.0"
+version = "0.1.1"
 
 [[IterativeSolvers]]
 deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"]
@@ -176,9 +182,9 @@ version = "0.9.1"
 
 [[JLD2]]
 deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"]
-git-tree-sha1 = "192934b3e2a94e897ce177423fd6cf7bdf464bce"
+git-tree-sha1 = "46b7834ec8165c541b0b5d1c8ba63ec940723ffb"
 uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-version = "0.4.14"
+version = "0.4.15"
 
 [[JLLWrappers]]
 deps = ["Preferences"]
@@ -234,10 +240,10 @@ deps = ["Libdl"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[LogExpFunctions]]
-deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a"
+deps = ["ChainRulesCore", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "6193c3815f13ba1b78a51ce391db8be016ae9214"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.3"
+version = "0.3.4"
 
 [[Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
@@ -335,13 +341,13 @@ version = "0.12.3"
 
 [[Parsers]]
 deps = ["Dates"]
-git-tree-sha1 = "a8709b968a1ea6abc2dc1967cb1db6ac9a00dfb6"
+git-tree-sha1 = "98f59ff3639b3d9485a03a72f3ab35bab9465720"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.0.5"
+version = "2.0.6"
 
 [[PartitionedArrays]]
 deps = ["Distances", "IterativeSolvers", "LinearAlgebra", "MPI", "Printf", "SparseArrays", "SparseMatricesCSR"]
-git-tree-sha1 = "763146e768b343fec44c521f5c7822116adc24ad"
+git-tree-sha1 = "4b79fc86fd30f5063d97d9c40291976ede6fd4d5"
 repo-rev = "gridap_distributed"
 repo-url = "https://github.com/fverdugo/PartitionedArrays.jl.git"
 uuid = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
@@ -409,16 +415,16 @@ deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [[SparseMatricesCSR]]
-deps = ["LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "dd6dda1484d4031a47d27614ceeb89f511cd9ca4"
+deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "ad906b39ce5e05ec509495dfc7b38d11ce9ab40b"
 uuid = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
-version = "0.6.1"
+version = "0.6.5"
 
 [[SpecialFunctions]]
 deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "793793f1df98e3d7d554b65a107e9c9a6399a6ed"
+git-tree-sha1 = "2d57e14cd614083f132b6224874296287bfa3979"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "1.7.0"
+version = "1.8.0"
 
 [[Static]]
 deps = ["IfElse"]
@@ -441,6 +447,10 @@ git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
 uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 version = "1.0.0"
 
+[[SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
 [[TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
diff --git a/Project.toml b/Project.toml
index 27dddced..e175e57e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,7 +13,9 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
+FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12"
 Gridap = "0.17"
+MPI = "0.16, 0.17, 0.18, 0.19"
 PartitionedArrays = "0.2"
 julia = "1.3"
 
diff --git a/src/Algebra.jl b/src/Algebra.jl
index 4e4a2b06..69cee515 100644
--- a/src/Algebra.jl
+++ b/src/Algebra.jl
@@ -380,37 +380,47 @@ function local_views(a::PVectorCounter,rows)
   a.counters
 end
 
-function Arrays.nz_allocation(a::PVectorCounter)
+function Arrays.nz_allocation(a::PVectorCounter{<:FullyAssembledRows})
   dofs = a.rows
   values = map_parts(nz_allocation,a.counters)
-  PVectorAllocation(a.par_strategy,values,dofs)
+  PVectorAllocationTrackOnlyValues(a.par_strategy,values,dofs)
 end
 
-struct PVectorAllocation{A,B,C}
+struct PVectorAllocationTrackOnlyValues{A,B,C}
   par_strategy::A
   values::B
   rows::C
 end
 
-function local_views(a::PVectorAllocation)
+function local_views(a::PVectorAllocationTrackOnlyValues)
   a.values
 end
 
-function local_views(a::PVectorAllocation,rows)
+function local_views(a::PVectorAllocationTrackOnlyValues,rows)
   @check rows === a.rows
   a.values
 end
 
-function Algebra.create_from_nz(a::PVectorAllocation{<:FullyAssembledRows})
-  @notimplemented
+function Algebra.create_from_nz(a::PVectorAllocationTrackOnlyValues{<:FullyAssembledRows})
+  # Create PRange for the rows of the linear system
+  parts = get_part_ids(a.values)
+  ngdofs = length(a.rows)
+  nodofs = map_parts(num_oids,a.rows.partition)
+  first_grdof = map_parts(first_gdof_from_ids,a.rows.partition)
+
+  # This one has no ghost rows
+  rows = PRange(parts,ngdofs,nodofs,first_grdof)
+
+  _rhs_callback(a,rows)
 end
 
-function Algebra.create_from_nz(a::PVectorAllocation{<:SubAssembledRows})
-  @notimplemented
+function Algebra.create_from_nz(a::PVectorAllocationTrackOnlyValues{<:SubAssembledRows})
+  # This point MUST NEVER be reached. If reached there is an inconsistency
+  # in the parallel code in charge of vector assembly
+  @assert false
 end
 
 function _rhs_callback(c_fespace,rows)
-
   # The ghost values in b_fespace are aligned with the FESpace
   # but not with the ghost values in the rows of A
   b_fespace = PVector(c_fespace.values,c_fespace.rows)
@@ -442,14 +452,14 @@ end
 
 function Algebra.create_from_nz(a::PVector)
   # For FullyAssembledRows the underlying Exchanger should
-  # not have ghost layer making asseble! do nothing (TODO check)
+  # not have ghost layer making assemble! do nothing (TODO check)
   assemble!(a)
   a
 end
 
 function Algebra.create_from_nz(
   a::DistributedAllocationCOO{<:FullyAssembledRows},
-  c_fespace::PVectorAllocation{<:FullyAssembledRows})
+  c_fespace::PVectorAllocationTrackOnlyValues{<:FullyAssembledRows})
 
   function callback(rows)
     _rhs_callback(c_fespace,rows)
@@ -459,9 +469,15 @@ function Algebra.create_from_nz(
   A,b
 end
 
+struct PVectorAllocationTrackTouchedAndValues{A,B,C}
+  allocations::A
+  values::B
+  rows::C
+end
+
 function Algebra.create_from_nz(
   a::DistributedAllocationCOO{<:SubAssembledRows},
-  c_fespace::PVectorAllocation{<:SubAssembledRows})
+  c_fespace::PVectorAllocationTrackOnlyValues{<:SubAssembledRows})
 
   function callback(rows)
     _rhs_callback(c_fespace,rows)
@@ -476,3 +492,119 @@ function Algebra.create_from_nz(
   A,b
 end
 
+struct ArrayAllocationTrackTouchedAndValues{A}
+  touched::Vector{Bool}
+  values::A
+end
+
+Gridap.Algebra.LoopStyle(::Type{<:ArrayAllocationTrackTouchedAndValues}) = Gridap.Algebra.Loop()
+
+
+function local_views(a::PVectorAllocationTrackTouchedAndValues,rows)
+  @check rows === a.rows
+  a.allocations
+end
+
+@inline function Arrays.add_entry!(c::Function,a::ArrayAllocationTrackTouchedAndValues,v,i,j)
+  @notimplemented
+end
+@inline function Arrays.add_entry!(c::Function,a::ArrayAllocationTrackTouchedAndValues,v,i)
+  if i>0
+    if !(a.touched[i])
+      a.touched[i]=true
+    end
+    if v!=nothing
+      a.values[i]=c(v,a.values[i])
+    end
+  end
+  nothing
+end
+@inline function Arrays.add_entries!(c::Function,a::ArrayAllocationTrackTouchedAndValues,v,i,j)
+  @notimplemented
+end
+@inline function Arrays.add_entries!(c::Function,a::ArrayAllocationTrackTouchedAndValues,v,i)
+  if v != nothing
+    for (ve,ie) in zip(v,i)
+      Arrays.add_entry!(c,a,ve,ie)
+    end
+  else
+    for ie in i
+      Arrays.add_entry!(c,a,nothing,ie)
+    end
+  end
+  nothing
+end
+
+function Arrays.nz_allocation(a::DistributedCounterCOO{<:SubAssembledRows},
+                              b::PVectorCounter{<:SubAssembledRows})
+  A = nz_allocation(a)
+  dofs = b.rows
+  values = map_parts(nz_allocation,b.counters)
+  B=PVectorAllocationTrackOnlyValues(b.par_strategy,values,dofs)
+  A,B
+end
+
+function Arrays.nz_allocation(a::PVectorCounter{<:SubAssembledRows})
+  dofs = a.rows
+  values = map_parts(nz_allocation,a.counters)
+  touched = map_parts(values) do values
+     fill!(Vector{Bool}(undef,length(values)),false)
+  end
+  allocations=map_parts(values,touched) do values,touched
+    ArrayAllocationTrackTouchedAndValues(touched,values)
+  end
+  PVectorAllocationTrackTouchedAndValues(allocations,values,dofs)
+end
+
+function local_views(a::PVectorAllocationTrackTouchedAndValues)
+  a.allocations
+end
+
+function Algebra.create_from_nz(a::PVectorAllocationTrackTouchedAndValues)
+   parts = get_part_ids(a.values)
+   rdofs = a.rows # dof ids of the test space
+   ngrdofs = length(rdofs)
+   nordofs = map_parts(num_oids,rdofs.partition)
+   first_grdof = map_parts(first_gdof_from_ids,rdofs.partition)
+   rneigs_snd = rdofs.exchanger.parts_snd
+   rneigs_rcv = rdofs.exchanger.parts_rcv
+
+   # Find the ghost rows
+   hrow_to_hrdof=map_parts(local_views(a.allocations),rdofs.partition) do allocation, indices
+    lids_touched=findall(allocation.touched)
+    nhlids = count((x)->indices.lid_to_ohid[x]<0,lids_touched)
+    hlids = Vector{Int32}(undef,nhlids)
+    cur=1
+    for lid in lids_touched
+      hlid=indices.lid_to_ohid[lid]
+      if hlid<0
+        hlids[cur]=-hlid
+        cur=cur+1
+      end
+    end
+    hlids
+   end
+   hrow_to_gid, hrow_to_part = map_parts(
+       find_gid_and_part,hrow_to_hrdof,rdofs.partition)
+
+   # Create the range for rows
+   rows = PRange(
+           parts,
+           ngrdofs,
+           nordofs,
+           first_grdof,
+           hrow_to_gid,
+           hrow_to_part,
+           rneigs_snd,
+           rneigs_rcv)
+
+   b = _rhs_callback(a,rows)
+   t2 = async_assemble!(b)
+
+   # Wait the transfer to finish
+   if t2 !== nothing
+     map_parts(schedule,t2)
+     map_parts(wait,t2)
+   end
+   b
+end
diff --git a/src/FESpaces.jl b/src/FESpaces.jl
index d62f72a4..a0af7f2a 100644
--- a/src/FESpaces.jl
+++ b/src/FESpaces.jl
@@ -365,12 +365,13 @@ function FESpaces.collect_cell_matrix_and_vector(
   end
 end
 
-struct DistributedSparseMatrixAssembler{A,B,C,D,E} <: SparseMatrixAssembler
-  assems::A
-  matrix_builder::B
-  vector_builder::C
-  rows::D
-  cols::E
+struct DistributedSparseMatrixAssembler{A,B,C,D,E,F} <: SparseMatrixAssembler
+  strategy::A
+  assems::B
+  matrix_builder::C
+  vector_builder::D
+  rows::E
+  cols::F
 end
 
 local_views(a::DistributedSparseMatrixAssembler) = a.assems
@@ -442,6 +443,5 @@ function FESpaces.SparseMatrixAssembler(
   vector_builder = PVectorBuilder(Tv,par_strategy)
   rows = get_free_dof_ids(test)
   cols = get_free_dof_ids(trial)
-  DistributedSparseMatrixAssembler(assems,matrix_builder,vector_builder,rows,cols)
+  DistributedSparseMatrixAssembler(par_strategy,assems,matrix_builder,vector_builder,rows,cols)
 end
-
diff --git a/src/Geometry.jl b/src/Geometry.jl
index a645475d..a66500d1 100644
--- a/src/Geometry.jl
+++ b/src/Geometry.jl
@@ -101,7 +101,7 @@ function Geometry.CartesianDiscreteModel(
   nc = desc.partition
   msg = """
   A CartesianDiscreteModel needs a Cartesian subdomain partition
-  of the rigth dimensions.
+  of the right dimensions.
   """
   @assert length(size(parts)) == length(nc) msg
   gcids = PCartesianIndices(parts,nc,PArrays.with_ghost)
@@ -241,7 +241,7 @@ end
 # For the moment remove_ghost_cells
 # refers to the triangulation faces
 # pointing into the ghost cells
-# in the underlying background discrete 
+# in the underlying background discrete
 # model. This might change when solving
 # multi-field PDEs with one of the fields
 # defined on the boundary (e.g. a Lagrange multiplier)
diff --git a/test/FESpacesTests.jl b/test/FESpacesTests.jl
index f052ab4b..787c2465 100644
--- a/test/FESpacesTests.jl
+++ b/test/FESpacesTests.jl
@@ -8,7 +8,11 @@ using PartitionedArrays
 using Test
 
 function main(parts)
+  main(parts,SubAssembledRows())
+  main(parts,FullyAssembledRows())
+end
 
+function main(parts,das)
   output = mkpath(joinpath(@__DIR__,"output"))
 
   domain = (0,4,0,4)
@@ -31,22 +35,25 @@ function main(parts)
   eh = u - uh
 
   dΩ = Measure(Ω,3)
-  cont = ∫( abs2(eh) )dΩ
+  cont  = ∫( abs2(eh) )dΩ
   @test sqrt(sum(cont)) < 1.0e-9
 
   # Assembly
+  Ωass  = Triangulation(das,model)
+  dΩass = Measure(Ωass,3)
   dv = get_fe_basis(V)
   du = get_trial_fe_basis(U)
-  a(u,v) = ∫( ∇(v)⋅∇(u) )dΩ
-  l(v) = ∫( 0*dv )dΩ
-  assem = SparseMatrixAssembler(U,V)
+  a(u,v) = ∫( ∇(v)⋅∇(u) )dΩass
+  l(v) = ∫( 0*dv )dΩass
+  assem = SparseMatrixAssembler(U,V,das)
 
   data = collect_cell_matrix_and_vector(U,V,a(du,dv),l(dv),zh)
-  A,b = assemble_matrix_and_vector(assem,data)
-  x = A\b
-  r = A*x -b
-  uh = FEFunction(U,x)
-  eh = u - uh
+  A1,b1 = assemble_matrix_and_vector(assem,data)
+  x1 = A1\b1
+  r1 = A1*x1 -b1
+  uh1 = FEFunction(U,x1)
+  eh1 = u - uh1
+  @test sqrt(sum(∫( abs2(eh1) )dΩ)) < 1.0e-9
 
   writevtk(Ω,joinpath(output,"Ω"), nsubcells=10,
     celldata=["err"=>cont[Ω]],
@@ -54,35 +61,44 @@ function main(parts)
 
   writevtk(Γ,joinpath(output,"Γ"),cellfields=["uh"=>uh])
 
-  @test sqrt(sum(∫( abs2(eh) )dΩ)) < 1.0e-9
+  A2,b2 = allocate_matrix_and_vector(assem,data)
+  assemble_matrix_and_vector!(A2,b2,assem,data)
+  x2 = A2\b2
+  r2 = A2*x2 -b2
+  uh = FEFunction(U,x2)
+  eh2 = u - uh
+  sqrt(sum(∫( abs2(eh2) )dΩ)) < 1.0e-9
 
-  op = AffineFEOperator(a,l,U,V)
+  op = AffineFEOperator(a,l,U,V,das)
   solver = LinearFESolver(BackslashSolver())
   uh = solve(solver,op)
   eh = u - uh
   @test sqrt(sum(∫( abs2(eh) )dΩ)) < 1.0e-9
 
   data = collect_cell_matrix(U,V,a(du,dv))
-  A2 = assemble_matrix(assem,data)
-
-  Ωl = Triangulation(with_ghost,model)
-  dΩl = Measure(Ωl,3)
-  al(u,v) = ∫( ∇(v)⋅∇(u) )dΩl
-  ll(v) = ∫( 0*v )dΩl
-
-  assem = SparseMatrixAssembler(U,V,FullyAssembledRows())
-  data = collect_cell_matrix_and_vector(U,V,al(du,dv),ll(dv),zh)
-  A,b = assemble_matrix_and_vector(assem,data)
-  x = A\b
-  r = A*x -b
-  uh = FEFunction(U,x)
-  eh = u - uh
-  @test sqrt(sum(∫( abs2(eh) )dΩ)) < 1.0e-9
+  A3 = assemble_matrix(assem,data)
+  x3 = A3\op.op.vector
+  uh = FEFunction(U,x3)
+  eh3 = u - uh
+  sqrt(sum(∫( abs2(eh3) )dΩ)) < 1.0e-9
+
+  A4 = allocate_matrix(assem,data)
+  assemble_matrix!(A4,assem,data)
+  x4 = A4\op.op.vector
+  uh = FEFunction(U,x4)
+  eh4 = u - uh
+  sqrt(sum(∫( abs2(eh4) )dΩ)) < 1.0e-9
 
-  op = AffineFEOperator(al,ll,U,V,FullyAssembledRows())
-  uh = solve(solver,op)
-  eh = u - uh
-  @test sqrt(sum(∫( abs2(eh) )dΩ)) < 1.0e-9
+  dv = get_fe_basis(V)
+  l=∫(1*dv)dΩass
+  vecdata=collect_cell_vector(V,l)
+  assem = SparseMatrixAssembler(U,V,das)
+  b1=assemble_vector(assem,vecdata)
+  @test abs(sum(b1)-length(b1)) < 1.0e-12
+
+  b2=allocate_vector(assem,vecdata)
+  assemble_vector!(b2,assem,vecdata)
+  @test abs(sum(b2)-length(b2)) < 1.0e-12
 
 end