diff --git a/Project.toml b/Project.toml
index 3967801c..806d44d1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FinEtools"
 uuid = "91bb5406-6c9a-523d-811d-0644c4229550"
 authors = ["Petr Krysl <pkrysl@ucsd.edu>"]
-version = "6.0.13"
+version = "6.0.14"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
diff --git a/README.md b/README.md
index 985e4255..53b3c153 100644
--- a/README.md
+++ b/README.md
@@ -21,6 +21,7 @@ The package supports application packages, for instance:
 
 ## News
 
+- 04/18/2023: Implement resizing of assembly buffers. Makematrix! resets pointers.
 - 04/16/2023: Enabled creation of finite element sets from arbitrary arrays.
 - 03/15/2023: Changed strategy when assembling into the COO format.
 - 03/10/2023: Genericity of arguments enabled in the assembly module.
diff --git a/src/AssemblyModule.jl b/src/AssemblyModule.jl
index ec990df2..92e52458 100644
--- a/src/AssemblyModule.jl
+++ b/src/AssemblyModule.jl
@@ -38,6 +38,7 @@ mutable struct SysmatAssemblerSparse{IT, MBT, IBT} <: AbstractSysmatAssembler
     ndofs_row::IT
     ndofs_col::IT
     nomatrixresult::Bool
+    force_init::Bool
 end
 
 """
@@ -97,8 +98,8 @@ At this point all the buffers of the assembler have been cleared, and
 `makematrix!(a) ` is no longer possible.
 
 """
-function SysmatAssemblerSparse(z= zero(FFlt), nomatrixresult = false)
-    return SysmatAssemblerSparse(0, FFlt[z], FInt[0], FInt[0], 0, 0, 0, nomatrixresult)
+function SysmatAssemblerSparse(z = zero(FFlt), nomatrixresult = false)
+    return SysmatAssemblerSparse(0, FFlt[z], FInt[0], FInt[0], 0, 0, 0, nomatrixresult, false)
 end
 
 """
@@ -156,6 +157,7 @@ function startassembly!(
     end
     # Leave the buffers uninitialized, unless the user requests otherwise
     if force_init
+        self.force_init = force_init
         self.rowbuffer .= 1
         self.colbuffer .= 1
         self.matbuffer .= 0
@@ -186,7 +188,18 @@ function assemble!(
     nrows = length(dofnums_row)
     ncolumns = length(dofnums_col)
     p = self.buffer_pointer
-    @assert p + ncolumns * nrows <= self.buffer_length + 1
+    if p + ncolumns * nrows >= self.buffer_length
+        new_buffer_length = self.buffer_length + ncolumns * nrows * 1000
+        resize!(self.rowbuffer, new_buffer_length)
+        resize!(self.colbuffer, new_buffer_length)
+        resize!(self.matbuffer, new_buffer_length)
+        if self.force_init
+            self.rowbuffer[self.buffer_length+1:end] .= 1
+            self.colbuffer[self.buffer_length+1:end] .= 1
+            self.matbuffer[self.buffer_length+1:end] .= 0
+        end
+        self.buffer_length = new_buffer_length
+    end
     @assert size(mat) == (nrows, ncolumns)
     @inbounds for j in 1:ncolumns
         if 0 < dofnums_col[j] <= self.ndofs_col
@@ -243,7 +256,8 @@ function makematrix!(self::SysmatAssemblerSparse)
         self.ndofs_row,
         self.ndofs_col,
         )
-    self = SysmatAssemblerSparse(zero(eltype(self.matbuffer)))# get rid of the buffers
+    # Get ready for more assembling
+    self.buffer_pointer = 1
     return S
 end
 
@@ -266,6 +280,7 @@ mutable struct SysmatAssemblerSparseSymm{IT, MBT, IBT} <: AbstractSysmatAssemble
     buffer_pointer::IT
     ndofs::IT
     nomatrixresult::Bool
+    force_init::Bool
 end
 
 """
@@ -296,7 +311,7 @@ This is how a symmetric sparse matrix is assembled from two square dense matrice
 SysmatAssemblerSparse
 """
 function SysmatAssemblerSparseSymm(z= zero(FFlt), nomatrixresult = false)
-    return SysmatAssemblerSparseSymm(0, FFlt[z], FInt[0], FInt[0], 0, 0, nomatrixresult)
+    return SysmatAssemblerSparseSymm(0, FFlt[z], FInt[0], FInt[0], 0, 0, nomatrixresult, false)
 end
 
 """
@@ -356,6 +371,7 @@ function startassembly!(
         self.rowbuffer .= 1
         self.colbuffer .= 1
         self.matbuffer .= 0
+        self.force_init = force_init
     end
     return self
 end
@@ -385,7 +401,18 @@ function assemble!(
     nrows = length(dofnums)
     ncolumns = nrows
     p = self.buffer_pointer
-    @assert p + ncolumns * nrows <= self.buffer_length + 1
+    if p + ncolumns * nrows >= self.buffer_length
+        new_buffer_length = self.buffer_length + ncolumns * nrows * 1000
+        resize!(self.rowbuffer, new_buffer_length)
+        resize!(self.colbuffer, new_buffer_length)
+        resize!(self.matbuffer, new_buffer_length)
+        if self.force_init
+            self.rowbuffer[self.buffer_length+1:end] .= 1
+            self.colbuffer[self.buffer_length+1:end] .= 1
+            self.matbuffer[self.buffer_length+1:end] .= 0
+        end
+        self.buffer_length = new_buffer_length
+    end
     @assert size(mat) == (nrows, ncolumns)
     @inbounds for j in 1:ncolumns
         if 0 < dofnums[j] <= self.ndofs
@@ -448,7 +475,8 @@ function makematrix!(self::SysmatAssemblerSparseSymm)
     @inbounds for j in axes(S, 1)
         S[j, j] *= 0.5      # the diagonal is there twice; fix it;
     end
-    self = SysmatAssemblerSparseSymm(zero(eltype(self.matbuffer))) # get rid of the buffers
+    # Get ready for more assembling
+    self.buffer_pointer = 1
     return S
 end
 
@@ -615,7 +643,8 @@ function makematrix!(self::SysmatAssemblerSparseDiag)
         self.ndofs,
         self.ndofs,
     )
-    self = SysmatAssemblerSparseDiag(zero(eltype(self.matbuffer))) # get rid of the buffers
+    # Get ready for more assembling
+    self.buffer_pointer = 1
     return S
 end
 
@@ -945,7 +974,8 @@ function makematrix!(self::SysmatAssemblerSparseHRZLumpingSymm)
     @inbounds for j in axes(S, 1)
         S[j, j] /= 2.0      # the diagonal is there twice; fix it;
     end
-    self = SysmatAssemblerSparseHRZLumpingSymm(zero(eltype(self.matbuffer)))# get rid of the buffers
+    # Get ready for more assembling
+        self.buffer_pointer = 1
     return S
 end
 
diff --git a/test/test_basics.jl b/test/test_basics.jl
index ad2db669..19801de0 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -40,14 +40,14 @@ function test()
     assemble!(a, m1, i1, i1)
     assemble!(a, m2, i2, i2)
     Au = makematrix!(a)
-    # @show norm(Au-testA)
+    @test maximum(abs.(testA - Au)) < 1.0e-5
 
     a = SysmatAssemblerSparseSymm(0.0)
     startassembly!(a, 5, 5, 3, 7, 7)
     assemble!(a, m1, i1, i1)
     assemble!(a, m2, i2, i2)
     A = makematrix!(a)
-    # @show norm(A-testA)
+    @test maximum(abs.(A - Au)) < 1.0e-5
 
     @test maximum(abs.(testA - A)) < 1.0e-5
     @test abs.(maximum(A - transpose(A))) < 1.0e-5
@@ -1611,7 +1611,7 @@ function _test()
         colbuffer = view(a.colbuffer, istart:iend)
         # @show length(colbuffer), istart, iend, buffer_length
         buffer_pointer = 1
-        a1 = SysmatAssemblerSparse(buffer_length, matbuffer, rowbuffer, colbuffer, buffer_pointer, ndofs_row, ndofs_col, true)
+        a1 = SysmatAssemblerSparse(buffer_length, matbuffer, rowbuffer, colbuffer, buffer_pointer, ndofs_row, ndofs_col, true, false)
         for i in ch[1]
             assemble!(a1, assembly_line[i]...)
         end
@@ -1639,7 +1639,7 @@ function _test()
             rowbuffer .= 1
             colbuffer .= 1
             buffer_pointer = 1
-            a1 = SysmatAssemblerSparse(buffer_length, matbuffer, rowbuffer, colbuffer, buffer_pointer, ndofs_row, ndofs_col, true)
+            a1 = SysmatAssemblerSparse(buffer_length, matbuffer, rowbuffer, colbuffer, buffer_pointer, ndofs_row, ndofs_col, true, false)
             Threads.@spawn let r =  $ch[1]
                 for i in r
                     assemble!(a1, assembly_line[i]...)
@@ -1745,7 +1745,7 @@ function _test()
         rowbuffer = view(a.rowbuffer, istart:iend)
         colbuffer = view(a.colbuffer, istart:iend)
         buffer_pointer = 1
-        a1 = SysmatAssemblerSparse(buffer_length, matbuffer, rowbuffer, colbuffer, buffer_pointer, ndofs_row, ndofs_col, true)
+        a1 = SysmatAssemblerSparse(buffer_length, matbuffer, rowbuffer, colbuffer, buffer_pointer, ndofs_row, ndofs_col, true, false)
         for i in ch[1]
             assemble!(a1, assembly_line[i]...)
         end
@@ -1775,7 +1775,7 @@ function _test()
             colbuffer = view(a.colbuffer, istart:iend)
             buffer_pointer = 1
             Threads.@spawn let r =  $ch[1]
-                a1 = SysmatAssemblerSparse(buffer_length, $matbuffer, $rowbuffer, $colbuffer, $buffer_pointer, ndofs_row, ndofs_col, true)
+                a1 = SysmatAssemblerSparse(buffer_length, $matbuffer, $rowbuffer, $colbuffer, $buffer_pointer, ndofs_row, ndofs_col, true, false)
                 # @show ch[2], r
                 for i in r
                     assemble!(a1, assembly_line[i]...)
@@ -1886,7 +1886,7 @@ function _test()
         rowbuffer .= 1
         colbuffer .= 1
         buffer_pointer = 1
-        a1 = SysmatAssemblerSparse(buffer_length, matbuffer, rowbuffer, colbuffer, buffer_pointer, ndofs_row, ndofs_col, true)
+        a1 = SysmatAssemblerSparse(buffer_length, matbuffer, rowbuffer, colbuffer, buffer_pointer, ndofs_row, ndofs_col, true, false)
         for i in ch[1]
             assemble!(a1, assembly_line[i]...)
         end
@@ -1915,7 +1915,7 @@ function _test()
         matbuffer .= 0.0
         rowbuffer .= 1
         colbuffer .= 1
-        a1 = SysmatAssemblerSparse(buffer_length, matbuffer, rowbuffer, colbuffer, buffer_pointer, ndofs_row, ndofs_col, true)
+        a1 = SysmatAssemblerSparse(buffer_length, matbuffer, rowbuffer, colbuffer, buffer_pointer, ndofs_row, ndofs_col, true, false)
     end
 
     # Parallel execution
@@ -1948,3 +1948,318 @@ end
 _test()
 
 end # module
+
+module mbasicass001
+
+using Test
+using LinearAlgebra
+using FinEtools
+using FinEtools.AssemblyModule
+using SparseArrays
+using Random
+
+function _test()
+    Random.seed!(1234);
+        # @show Base.Threads.nthreads()
+    m = [
+    0.24406   0.599773    0.833404  0.0420141
+    0.786024  0.00206713  0.995379  0.780298
+    0.845816  0.198459    0.355149  0.224996
+    ]
+    m1 = m'*m;
+    m5 = m*m';
+    m = [
+    0.146618  -0.53471   0.614342    0.737833
+    0.479719  0.41354   -0.00760941  0.836455
+    0.254868  0.476189  0.460794    0.00919633
+    0.159064  0.261821  0.317078    0.77646
+    -0.643538  0.429817  0.59788     0.958909
+    ]
+    m2 = m'*m;
+    m = [
+    -0.146618  0.53471   0.614342    0.737833
+    0.479719  0.41354   0.00760941  -0.836455
+    0.479719  0.41354   0.00760941  -0.836455
+    0.254868  -0.476189  0.460794    0.00919633
+    0.159064  0.261821  0.317078    0.77646
+    0.159064  0.261821  0.317078    0.77646
+    ]
+    m3 = m'*m;
+    m4 = Matrix(m*m');
+    ms = [m1, m2, m3, m4, m5, ]
+    # Testing symmetric matrices
+    for i in eachindex(ms)
+        m = ms[i][1]
+        @assert norm(m - m') == 0
+    end
+    elem_mat_nrows = maximum(size.(ms, 1))
+    elem_mat_ncols = maximum(size.(ms, 2))
+    N = 133
+    p = randperm(N)
+    assembly_line = []
+    for i in 1:N
+        for k in eachindex(ms)
+            m = ms[k]
+            if length(p) <= sum(size(m))
+                p = randperm(N)
+            end
+            r = [popfirst!(p) for _ in 1:size(m, 1)]
+            c = [popfirst!(p) for _ in 1:size(m, 2)]
+            # Symmetric matrices: the rows are the same as columns
+            push!(assembly_line, (m, r, r))
+            push!(assembly_line, (m, c, c))
+        end
+    end
+
+    elem_mat_nmatrices = length(assembly_line)
+    ndofs_row = N
+    ndofs_col = N
+
+    refA = spzeros(ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        m, r, c = assembly_line[i]
+        @assert length(r) == length(c)
+        refA[r, c] += m
+    end
+
+    a = SysmatAssemblerSparse(0.0)
+    startassembly!(a, elem_mat_nrows, elem_mat_ncols, elem_mat_nmatrices, ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        assemble!(a, assembly_line[i]...)
+    end
+    A = makematrix!(a)
+    # @show time() - start
+    @test norm(A - refA) / norm(refA) < 1.0e-9
+
+
+    a = SysmatAssemblerSparseSymm(0.0)
+    startassembly!(a, elem_mat_nrows, elem_mat_ncols, elem_mat_nmatrices, ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        assemble!(a, assembly_line[i]...)
+    end
+    A = makematrix!(a)
+
+    # @show time() - start
+    @test norm(A - refA) / norm(refA) < 1.0e-9
+
+    return A, refA
+end
+
+_test()
+
+end # module
+
+module mbasicass002
+using Test
+using LinearAlgebra
+using FinEtools
+using FinEtools.AssemblyModule
+using SparseArrays
+using Random
+
+function _test()
+    Random.seed!(1234);
+        # @show Base.Threads.nthreads()
+    m = [
+    0.24406   0.599773    0.833404  0.0420141
+    0.786024  0.00206713  0.995379  0.780298
+    0.845816  0.198459    0.355149  0.224996
+    ]
+    m1 = m'*m;
+    m5 = m*m';
+    m = [
+    0.146618  -0.53471   0.614342    0.737833
+    0.479719  0.41354   -0.00760941  0.836455
+    0.254868  0.476189  0.460794    0.00919633
+    0.159064  0.261821  0.317078    0.77646
+    -0.643538  0.429817  0.59788     0.958909
+    ]
+    m2 = m'*m;
+    m = [
+    -0.146618  0.53471   0.614342    0.737833
+    0.479719  0.41354   0.00760941  -0.836455
+    0.479719  0.41354   0.00760941  -0.836455
+    0.254868  -0.476189  0.460794    0.00919633
+    0.159064  0.261821  0.317078    0.77646
+    0.159064  0.261821  0.317078    0.77646
+    ]
+    m3 = m'*m;
+    m4 = Matrix(m*m');
+    ms = [m1, m2, m3, m4, m5, ]
+    # Testing symmetric matrices
+    for i in eachindex(ms)
+        m = ms[i][1]
+        @assert norm(m - m') == 0
+    end
+    elem_mat_nrows = maximum(size.(ms, 1))
+    elem_mat_ncols = maximum(size.(ms, 2))
+    N = 1333
+    p = randperm(N)
+    assembly_line = []
+    for i in 1:N
+        for k in eachindex(ms)
+            m = ms[k]
+            if length(p) <= sum(size(m))
+                p = randperm(N)
+            end
+            r = [popfirst!(p) for _ in 1:size(m, 1)]
+            c = [popfirst!(p) for _ in 1:size(m, 2)]
+            # Symmetric matrices: the rows are the same as columns
+            push!(assembly_line, (m, r, r))
+            push!(assembly_line, (m, c, c))
+        end
+    end
+
+    elem_mat_nmatrices = length(assembly_line)
+    ndofs_row = N
+    ndofs_col = N
+
+    refA = spzeros(ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        m, r, c = assembly_line[i]
+        @assert length(r) == length(c)
+        refA[r, c] += m
+    end
+
+    a = SysmatAssemblerSparse(0.0)
+    # We are testing resizing of buffers by under sizing initially
+    startassembly!(a, elem_mat_nrows, elem_mat_ncols, 1, ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        assemble!(a, assembly_line[i]...)
+    end
+    A = makematrix!(a)
+    # @show time() - start
+    @test norm(A - refA) / norm(refA) < 1.0e-9
+
+
+    a = SysmatAssemblerSparseSymm(0.0)
+    # We are testing resizing of buffers by under sizing initially
+    startassembly!(a, elem_mat_nrows, elem_mat_ncols, 1, ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        assemble!(a, assembly_line[i]...)
+    end
+    A = makematrix!(a)
+
+    # @show time() - start
+    @test norm(A - refA) / norm(refA) < 1.0e-9
+
+    return A, refA
+end
+
+_test()
+
+end # module
+
+module mbasicass003
+
+using Test
+using LinearAlgebra
+using FinEtools
+using FinEtools.AssemblyModule
+using SparseArrays
+using Random
+
+function _test()
+    Random.seed!(1234);
+        # @show Base.Threads.nthreads()
+    m = [
+    0.24406   0.599773    0.833404  0.0420141
+    0.786024  0.00206713  0.995379  0.780298
+    0.845816  0.198459    0.355149  0.224996
+    ]
+    m1 = m'*m;
+    m5 = m*m';
+    m = [
+    0.146618  -0.53471   0.614342    0.737833
+    0.479719  0.41354   -0.00760941  0.836455
+    0.254868  0.476189  0.460794    0.00919633
+    0.159064  0.261821  0.317078    0.77646
+    -0.643538  0.429817  0.59788     0.958909
+    ]
+    m2 = m'*m;
+    m = [
+    -0.146618  0.53471   0.614342    0.737833
+    0.479719  0.41354   0.00760941  -0.836455
+    0.479719  0.41354   0.00760941  -0.836455
+    0.254868  -0.476189  0.460794    0.00919633
+    0.159064  0.261821  0.317078    0.77646
+    0.159064  0.261821  0.317078    0.77646
+    ]
+    m3 = m'*m;
+    m4 = Matrix(m*m');
+    ms = [m1, m2, m3, m4, m5, ]
+    # Testing symmetric matrices
+    for i in eachindex(ms)
+        m = ms[i][1]
+        @assert norm(m - m') == 0
+    end
+    elem_mat_nrows = maximum(size.(ms, 1))
+    elem_mat_ncols = maximum(size.(ms, 2))
+    N = 1333
+    p = randperm(N)
+    assembly_line = []
+    for i in 1:N
+        for k in eachindex(ms)
+            m = ms[k]
+            if length(p) <= sum(size(m))
+                p = randperm(N)
+            end
+            r = [popfirst!(p) for _ in 1:size(m, 1)]
+            c = [popfirst!(p) for _ in 1:size(m, 2)]
+            # Symmetric matrices: the rows are the same as columns
+            push!(assembly_line, (m, r, r))
+            push!(assembly_line, (m, c, c))
+        end
+    end
+
+    elem_mat_nmatrices = length(assembly_line)
+    ndofs_row = N
+    ndofs_col = N
+
+    refA = spzeros(ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        m, r, c = assembly_line[i]
+        @assert length(r) == length(c)
+        refA[r, c] += m
+    end
+
+    a = SysmatAssemblerSparse(0.0)
+    startassembly!(a, elem_mat_nrows, elem_mat_ncols, 110, ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        assemble!(a, assembly_line[i]...)
+    end
+    A = makematrix!(a)
+    # @show time() - start
+    @test norm(A - refA) / norm(refA) < 1.0e-9
+    # Here we test that we can start assembly again
+    startassembly!(a, elem_mat_nrows, elem_mat_ncols, 110, ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        assemble!(a, assembly_line[i]...)
+    end
+    A = makematrix!(a)
+    # @show time() - start
+    @test norm(A - refA) / norm(refA) < 1.0e-9
+
+    a = SysmatAssemblerSparseSymm(0.0)
+    startassembly!(a, elem_mat_nrows, elem_mat_ncols, 1, ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        assemble!(a, assembly_line[i]...)
+    end
+    A = makematrix!(a)
+    # @show time() - start
+    @test norm(A - refA) / norm(refA) < 1.0e-9
+    # Here we test that we can start assembly again
+    startassembly!(a, elem_mat_nrows, elem_mat_ncols, 1, ndofs_row, ndofs_col)
+    for i in eachindex(assembly_line)
+        assemble!(a, assembly_line[i]...)
+    end
+    A = makematrix!(a)
+    # @show time() - start
+    @test norm(A - refA) / norm(refA) < 1.0e-9
+
+    return A, refA
+end
+
+_test()
+
+end # module