diff --git a/NEWS.md b/NEWS.md
index 72f923b96..6ecdcc3d4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -9,6 +9,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixed
 - Inheritance relationship for DiscreteModelPortion. Since PR [#645](https://github.com/gridap/Gridap.jl/pull/645).
 
+### Added
+- Additional high level API assembly functions. Since PR [#652](https://github.com/gridap/Gridap.jl/pull/652).
+
 ## [0.16.4] - 2021-08-17
 
 ### Added
diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl
index 87b48504b..5fc45f834 100644
--- a/src/FESpaces/Assemblers.jl
+++ b/src/FESpaces/Assemblers.jl
@@ -270,44 +270,88 @@ function assemble_matrix(f::Function,a::Assembler,U::FESpace,V::FESpace)
   assemble_matrix(a,collect_cell_matrix(U,V,f(u,v)))
 end
 
+function assemble_matrix!(f::Function,A::AbstractMatrix,a::Assembler,U::FESpace,V::FESpace)
+  v = get_fe_basis(V)
+  u = get_trial_fe_basis(U)
+  assemble_matrix!(A,a,collect_cell_matrix(U,V,f(u,v)))
+end
+
 function assemble_vector(f::Function,a::Assembler,V::FESpace)
   v = get_fe_basis(V)
   assemble_vector(a,collect_cell_vector(V,f(v)))
 end
 
+function assemble_vector!(f::Function,b::AbstractVector,a::Assembler,V::FESpace)
+  v = get_fe_basis(V)
+  assemble_vector!(b,a,collect_cell_vector(V,f(v)))
+end
+
 function assemble_matrix_and_vector(f::Function,b::Function,a::Assembler,U::FESpace,V::FESpace)
   v = get_fe_basis(V)
   u = get_trial_fe_basis(U)
   assemble_matrix_and_vector(a,collect_cell_matrix_and_vector(U,V,f(u,v),b(v)))
 end
 
+function assemble_matrix_and_vector!(f::Function,b::Function,M::AbstractMatrix,r::AbstractVector,a::Assembler,U::FESpace,V::FESpace)
+  v = get_fe_basis(V)
+  u = get_trial_fe_basis(U)
+  assemble_matrix_and_vector!(M,r,a,collect_cell_matrix_and_vector(U,V,f(u,v),b(v)))
+end
+
 function assemble_matrix(f,a::Assembler,U::FESpace,V::FESpace)
   assemble_matrix(a,collect_cell_matrix(U,V,f))
 end
 
+function assemble_matrix!(f,A::AbstractMatrix,a::Assembler,U::FESpace,V::FESpace)
+  assemble_matrix!(A,a,collect_cell_matrix(U,V,f))
+end
+
 function assemble_vector(f,a::Assembler,V::FESpace)
   assemble_vector(a,collect_cell_vector(V,f))
 end
 
+function assemble_vector!(f,b::AbstractVector,a::Assembler,V::FESpace)
+  assemble_vector!(b,a,collect_cell_vector(V,f))
+end
+
 function assemble_matrix_and_vector(f,b,a::Assembler,U::FESpace,V::FESpace)
   assemble_matrix_and_vector(a,collect_cell_matrix_and_vector(U,V,f,b))
 end
 
+function assemble_matrix_and_vector!(f,b,M::AbstractMatrix,r::AbstractVector,a::Assembler,U::FESpace,V::FESpace)
+  assemble_matrix_and_vector!(M,r,a,collect_cell_matrix_and_vector(U,V,f,b))
+end
+
 function assemble_matrix(f,U::FESpace,V::FESpace)
   a = SparseMatrixAssembler(U,V)
   assemble_matrix(f,a,U,V)
 end
 
+function assemble_matrix!(f,A::AbstractMatrix,U::FESpace,V::FESpace)
+  a = SparseMatrixAssembler(U,V)
+  assemble_matrix!(f,A,a,U,V)
+end
+
 function assemble_vector(f,V::FESpace)
   a = SparseMatrixAssembler(V,V)
   assemble_vector(f,a,V)
 end
 
+function assemble_vector!(f,b::AbstractVector,V::FESpace)
+  a = SparseMatrixAssembler(V,V)
+  assemble_vector!(f,b,a,V)
+end
+
 function assemble_matrix_and_vector(f,b,U::FESpace,V::FESpace)
   a = SparseMatrixAssembler(U,V)
   assemble_matrix_and_vector(f,b,a,U,V)
 end
 
+function assemble_matrix_and_vector!(f,b,M::AbstractMatrix,r::AbstractVector,U::FESpace,V::FESpace)
+  a = SparseMatrixAssembler(U,V)
+  assemble_matrix_and_vector!(f,b,M,r,a,U,V)
+end
+
 # Abstract interface for computing the data to be sent to the assembler
 
 function collect_cell_matrix(trial::FESpace,test::FESpace,mat_contributions)
diff --git a/test/FESpacesTests/AssemblersTests.jl b/test/FESpacesTests/AssemblersTests.jl
index 70696182a..1e2bee7f8 100644
--- a/test/FESpacesTests/AssemblersTests.jl
+++ b/test/FESpacesTests/AssemblersTests.jl
@@ -65,13 +65,40 @@ A,b = assemble_matrix_and_vector(assem,data)
 x = A\b
 uh = FEFunction(U,x)
 
-A = assemble_matrix(a,U,V)
-b = assemble_vector(ℓ,V)
-A,b = assemble_matrix_and_vector(a,ℓ,U,V)
-
-A = assemble_matrix(a(du,dv),U,V)
-b = assemble_vector(ℓ(dv),V)
-A,b = assemble_matrix_and_vector(a(du,dv),ℓ(dv),U,V)
+A1 = assemble_matrix(a,U,V)
+b1 = assemble_vector(ℓ,V)
+A2,b2 = assemble_matrix_and_vector(a,ℓ,U,V)
+
+A12 = copy(A1); A12[1,1]=rand()
+b12 = copy(b1); b12[1]=rand()
+A22 = copy(A2); A22[1,1]=rand()
+b22 = copy(b2); b22[1]=rand()
+
+tol = 1.e-14
+assemble_matrix!(a,A12,U,V)
+assemble_vector!(ℓ,b12,V)
+assemble_matrix_and_vector!(a,ℓ,A22,b22,U,V)
+@test norm(A12-A1) < tol
+@test norm(b12-b1) < tol
+@test norm(A22-A2) < tol
+@test norm(b22-b2) < tol
+
+A1 = assemble_matrix(a(du,dv),U,V)
+b1 = assemble_vector(ℓ(dv),V)
+A2,b2 = assemble_matrix_and_vector(a(du,dv),ℓ(dv),U,V)
+
+A12 = copy(A1); A12[1,1]=rand()
+b12 = copy(b1); b12[1]=rand()
+A22 = copy(A2); A22[1,1]=rand()
+b22 = copy(b2); b22[1]=rand()
+
+assemble_matrix!(a(du,dv),A12,U,V)
+assemble_vector!(ℓ(dv),b12,V)
+assemble_matrix_and_vector!(a(du,dv),ℓ(dv),A22,b22,U,V)
+@test norm(A12-A1) < tol
+@test norm(b12-b1) < tol
+@test norm(A22-A2) < tol
+@test norm(b22-b2) < tol
 
 V = TestFESpace(
   model,