From 52f0a23611665d53f580148deebde981f35b394a Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Fri, 25 Oct 2019 15:08:40 +0200 Subject: [PATCH 1/7] Add SparseMatricesCSR SparseMatrixAssembler produces SparseMatrices also from SparseMatricesCSR packages --- Manifest.toml | 6 +++++ Project.toml | 1 + src/FESpaces/Assemblers.jl | 55 +++++++++++++++++++++++++------------- 3 files changed, 44 insertions(+), 18 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index dec049782..7e5179bbc 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -245,6 +245,12 @@ uuid = "6462fe0b-24de-5631-8697-dd941f90decc" deps = ["LinearAlgebra", "Random"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +[[SparseMatricesCSR]] +deps = ["LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "7896f9c00e5413cbaecd758ade9fe67bf50ad32e" +uuid = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" +version = "0.3.1" + [[SpecialFunctions]] deps = ["BinDeps", "BinaryProvider", "Libdl"] git-tree-sha1 = "3bdd374b6fd78faf0119b8c5d538788dbf910c6e" diff --git a/Project.toml b/Project.toml index 3340c86cd..fcad92900 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TensorPolynomialBases = "f77b161d-e03a-5d68-b365-d5306aabeaa6" TensorValues = "31c64edf-cdeb-50e4-845d-ed2334360c20" diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index bceea139a..a1dcd2e53 100644 --- a/src/FESpaces/Assemblers.jl +++ b/src/FESpaces/Assemblers.jl @@ -5,6 +5,7 @@ using Gridap using Gridap.Helpers using SparseArrays +using SparseMatricesCSR export Assembler export SparseMatrixAssembler @@ -80,21 +81,28 @@ function assemble!(r,a::Assembler,cv::CellMatrix) end """ -Assembler that produces SparseMatrices from the SparseArrays package +Assembler that produces SparseMatrices from the SparseArrays and SparseMatricesCSR packages """ -struct SparseMatrixAssembler{E} <: Assembler{SparseMatrixCSC{E,Int},Vector{E}} +struct SparseMatrixAssembler{E,M} <: Assembler{AbstractSparseMatrix{E,Int},Vector{E}} testfesp::FESpace trialfesp::FESpace + + function SparseMatrixAssembler( + ::Type{M}, + testfesp::FESpace{D,Z,T}, + trialfesp::FESpace{D,Z,T}) where {M<:AbstractSparseMatrix,D,Z,T} + E = eltype(T) + new{E,M}(testfesp,trialfesp) + end end function SparseMatrixAssembler(test::FESpace{D,Z,T}, trial::FESpace{D,Z,T}) where {D,Z,T} - E = eltype(T) - SparseMatrixAssembler{E}(trial,test) + SparseMatrixAssembler(SparseMatrixCSC, trial,test) end function assemble( - this::SparseMatrixAssembler{E}, - vals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where E + this::SparseMatrixAssembler{E,M}, + vals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where {E,M} n = num_free_dofs(this.testfesp) vec = zeros(E,n) @@ -104,8 +112,8 @@ end function assemble!( vec::Vector{E}, - this::SparseMatrixAssembler{E}, - allvals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where E + this::SparseMatrixAssembler{E,M}, + allvals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where {E,M} vec .= zero(E) rows = celldofids(this.testfesp) @@ -127,8 +135,8 @@ function _assemble_vector!(vec,vals,rows) end function assemble( - this::SparseMatrixAssembler{E}, - allvals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where E + this::SparseMatrixAssembler{E,M}, + allvals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where {E,M} I = Int aux_row = I[]; aux_col = I[]; aux_val = E[] @@ -141,21 +149,19 @@ function assemble( rows_m = reindex(_rows_m, cellids_row) vals_m = apply_constraints_cols(this.trialfesp, _vals, cellids_col) cols_m = reindex(_cols_m, cellids_col) - _assemble_sparse_matrix_values!( + _assemble_sparse_matrix_values!(M, aux_row,aux_col,aux_val,vals_m,rows_m,cols_m) end sparse(aux_row,aux_col,aux_val) end -function _assemble_sparse_matrix_values!(aux_row,aux_col,aux_val,vals,rows,cols) +function _assemble_sparse_matrix_values!(::Type{M},aux_row,aux_col,aux_val,vals,rows,cols) where {M} for (rows_c, cols_c, vals_c) in zip(rows,cols,vals) for (j,gidcol) in enumerate(cols_c) if gidcol > 0 for (i,gidrow) in enumerate(rows_c) if gidrow > 0 - push!(aux_row, gidrow) - push!(aux_col, gidcol) - push!(aux_val, vals_c[i,j]) + push_coo!(M,aux_row, aux_col, aux_val, gidrow, gidcol, vals_c[i,j]) end end end @@ -165,12 +171,25 @@ end function assemble!( mat::SparseMatrixCSC{E}, - this::SparseMatrixAssembler{E}, - vals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where E + this::SparseMatrixAssembler{E,M}, + vals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where {E,M} # This routine can be optimized a lot taking into a count the sparsity graph of mat # For the moment we create an intermediate matrix and then transfer the nz values m = assemble(this,vals...) - mat.nzval .= m.nzval + mat.nzval .= nonzeros(m) end +function sparse_from_coo(::Type{<:SparseMatrixCSC}, args...) + sparse(args...) +end + +function sparse_from_coo(::Type{<:SparseMatrixCSR}, args...) + sparsecsr(args...) +end + +function sparse_from_coo(::Type{<:SymSparseMatrixCSR}, args...) + symsparsecsr(args...) +end + + end # module From f10a40870ae4ee10e95dad0d770ee69e301218c9 Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Mon, 28 Oct 2019 10:50:56 +0100 Subject: [PATCH 2/7] Finish SparseMatricesCSR integration into Assemblers.jl Add SparseMAtrixCSR and SymSparseMatrixCSR to AssemblersTests.jl --- src/FESpaces/Assemblers.jl | 8 +-- test/FESpacesTests/AssemblersTests.jl | 86 ++++++++++++++------------- 2 files changed, 49 insertions(+), 45 deletions(-) diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index a1dcd2e53..f6ac8697a 100644 --- a/src/FESpaces/Assemblers.jl +++ b/src/FESpaces/Assemblers.jl @@ -143,7 +143,7 @@ function assemble( _rows_m = celldofids(this.testfesp) _cols_m = celldofids(this.trialfesp) - +@show allvals for (vals, cellids_row, cellids_col) in allvals _vals = apply_constraints_rows(this.testfesp, vals, cellids_row) rows_m = reindex(_rows_m, cellids_row) @@ -152,7 +152,7 @@ function assemble( _assemble_sparse_matrix_values!(M, aux_row,aux_col,aux_val,vals_m,rows_m,cols_m) end - sparse(aux_row,aux_col,aux_val) + sparse_from_coo(M,aux_row,aux_col,aux_val) end function _assemble_sparse_matrix_values!(::Type{M},aux_row,aux_col,aux_val,vals,rows,cols) where {M} @@ -170,13 +170,13 @@ function _assemble_sparse_matrix_values!(::Type{M},aux_row,aux_col,aux_val,vals, end function assemble!( - mat::SparseMatrixCSC{E}, + mat::AbstractSparseMatrix{E,Int}, this::SparseMatrixAssembler{E,M}, vals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where {E,M} # This routine can be optimized a lot taking into a count the sparsity graph of mat # For the moment we create an intermediate matrix and then transfer the nz values m = assemble(this,vals...) - mat.nzval .= nonzeros(m) + nonzeros(mat) .= nonzeros(m) end function sparse_from_coo(::Type{<:SparseMatrixCSC}, args...) diff --git a/test/FESpacesTests/AssemblersTests.jl b/test/FESpacesTests/AssemblersTests.jl index 60ba9eb5c..3017624ed 100644 --- a/test/FESpacesTests/AssemblersTests.jl +++ b/test/FESpacesTests/AssemblersTests.jl @@ -2,76 +2,80 @@ module AssemblersTests using Test using Gridap +using SparseArrays +using SparseMatricesCSR -model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,2)) +for T in [SparseMatrixCSC, SparseMatrixCSR, SymSparseMatrixCSR] + model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,2)) -tags = [1,2,3,4,6,5] + tags = [1,2,3,4,6,5] -order = 1 -fespace = H1ConformingFESpace(Float64,model,order,tags) + order = 1 + fespace = H1ConformingFESpace(Float64,model,order,tags) -assem = SparseMatrixAssembler(fespace,fespace) + assem = SparseMatrixAssembler(T,fespace,fespace) -trian = Triangulation(model) + trian = Triangulation(model) -quad = CellQuadrature(trian,degree=2) + quad = CellQuadrature(trian,degree=2) -basis = CellBasis(fespace) + basis = CellBasis(fespace) -a(v,u) = varinner(∇(v), ∇(u)) + a(v,u) = varinner(∇(v), ∇(u)) -bfun(x) = x[2] + bfun(x) = x[2] -b(v) = varinner(v,CellField(trian,bfun)) + b(v) = varinner(v,CellField(trian,bfun)) -mmat = integrate(a(basis,basis),trian,quad) + mmat = integrate(a(basis,basis),trian,quad) -bvec = integrate(b(basis),trian,quad) + bvec = integrate(b(basis),trian,quad) -vec = assemble(assem, bvec) + vec = assemble(assem, bvec) -mat = assemble(assem, mmat) + mat = assemble(assem, mmat) -x = mat \ vec + x = mat \ vec -assemble!(vec,assem, bvec) + assemble!(vec,assem, bvec) -assemble!(mat,assem, mmat) + assemble!(mat,assem, mmat) -x2 = mat \ vec + x2 = mat \ vec -@test x ≈ x2 + @test x ≈ x2 -@test vec ≈ [0.0625, 0.125, 0.0625] + @test vec ≈ [0.0625, 0.125, 0.0625] -@test mat[1, 1] ≈ 1.333333333333333 -@test mat[2, 1] ≈ -0.33333333333333 -@test mat[1, 2] ≈ -0.33333333333333 -@test mat[2, 2] ≈ 2.666666666666666 -@test mat[3, 2] ≈ -0.33333333333333 -@test mat[2, 3] ≈ -0.33333333333333 -@test mat[3, 3] ≈ 1.333333333333333 + @test mat[1, 1] ≈ 1.333333333333333 + @test mat[2, 1] ≈ -0.33333333333333 + @test mat[1, 2] ≈ -0.33333333333333 + @test mat[2, 2] ≈ 2.666666666666666 + @test mat[3, 2] ≈ -0.33333333333333 + @test mat[2, 3] ≈ -0.33333333333333 + @test mat[3, 3] ≈ 1.333333333333333 -cellids = IdentityCellNumber(Int,length(bvec)) + cellids = IdentityCellNumber(Int,length(bvec)) -vec2 = assemble(assem, (bvec,cellids)) + vec2 = assemble(assem, (bvec,cellids)) -@test vec2 == vec + @test vec2 == vec -vec3 = assemble(assem, (bvec,cellids), (bvec,cellids)) -@test vec3 ≈ 2*vec + vec3 = assemble(assem, (bvec,cellids), (bvec,cellids)) + @test vec3 ≈ 2*vec -assemble!(vec3, assem, (bvec,cellids), (bvec,cellids)) -@test vec3 ≈ 2*vec + assemble!(vec3, assem, (bvec,cellids), (bvec,cellids)) + @test vec3 ≈ 2*vec -mat2 = assemble(assem, (mmat,cellids,cellids) ) + mat2 = assemble(assem, (mmat,cellids,cellids) ) -@test mat2 == mat + @test mat2 == mat -mat3 = assemble(assem, (mmat,cellids,cellids), (mmat,cellids,cellids)) -@test mat3 ≈ 2*mat + mat3 = assemble(assem, (mmat,cellids,cellids), (mmat,cellids,cellids)) + @test mat3 ≈ 2*mat -assemble!(mat3, assem, (mmat,cellids,cellids), (mmat,cellids,cellids)) -@test mat3 ≈ 2*mat + assemble!(mat3, assem, (mmat,cellids,cellids), (mmat,cellids,cellids)) + @test mat3 ≈ 2*mat +end end # module AssemblersTests From 0c48dd0833f3cc929d95c16bf06bce42333b7740 Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Mon, 28 Oct 2019 15:19:00 +0100 Subject: [PATCH 3/7] MultiSparseMatrixAssembler produces SparseMatrices also from SparseMatricesCSR packages Finish SparseMatricesCSR integration into MultiAssemblers.jl Add SparseMAtrixCSR and SymSparseMatrixCSR to MultiAssemblersTests.jl --- src/FESpaces/Assemblers.jl | 3 +- src/MultiField/MultiAssemblers.jl | 80 +++++++++++----- test/MultiFieldTests/MultiAssemblersTests.jl | 96 ++++++++++---------- 3 files changed, 111 insertions(+), 68 deletions(-) diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index f6ac8697a..17ffa85b4 100644 --- a/src/FESpaces/Assemblers.jl +++ b/src/FESpaces/Assemblers.jl @@ -11,6 +11,8 @@ export Assembler export SparseMatrixAssembler export assemble export assemble! +export push_coo! +export sparse_from_coo """ Abstract assembly operator @@ -143,7 +145,6 @@ function assemble( _rows_m = celldofids(this.testfesp) _cols_m = celldofids(this.trialfesp) -@show allvals for (vals, cellids_row, cellids_col) in allvals _vals = apply_constraints_rows(this.testfesp, vals, cellids_row) rows_m = reindex(_rows_m, cellids_row) diff --git a/src/MultiField/MultiAssemblers.jl b/src/MultiField/MultiAssemblers.jl index f3103863a..76c87139d 100644 --- a/src/MultiField/MultiAssemblers.jl +++ b/src/MultiField/MultiAssemblers.jl @@ -12,6 +12,9 @@ export MultiSparseMatrixAssembler import Gridap.Assemblers: assemble import Gridap.Assemblers: assemble! import Gridap.Assemblers: SparseMatrixAssembler +import Gridap.Assemblers: SparseMatrixAssembler +import Gridap.Assemblers: push_coo! +import Gridap.Assemblers: sparse_from_coo abstract type MultiAssembler{M<:AbstractMatrix,V<:AbstractVector} end @@ -70,34 +73,70 @@ end """ MultiAssembler that produces SparseMatrices from the SparseArrays package """ -struct MultiSparseMatrixAssembler{E} <: MultiAssembler{SparseMatrixCSC{E,Int},Vector{E}} +struct MultiSparseMatrixAssembler{E,M} <: MultiAssembler{AbstractSparseMatrix{E,Int},Vector{E}} testfesps::MultiFESpace{E} trialfesps::MultiFESpace{E} + + function MultiSparseMatrixAssembler( + ::Type{M}, + testfesp::MultiFESpace{E}, + trialfesp::MultiFESpace{E}) where {M<:AbstractSparseMatrix,E} + new{E,M}(testfesp,trialfesp) + end +end + +function MultiSparseMatrixAssembler( + testfesp::MultiFESpace{E}, + trialfesp::MultiFESpace{E}) where {E} + MultiSparseMatrixAssembler(SparseMatrixCSC, testfesps, trialfesp) end function MultiSparseMatrixAssembler( + ::Type{M}, testfesps::Vector{<:FESpaceWithDirichletData}, - trialfesps::Vector{<:FESpaceWithDirichletData}) + trialfesps::Vector{<:FESpaceWithDirichletData}) where {M<:AbstractSparseMatrix} V = MultiFESpace(testfesps) U = MultiFESpace(trialfesps) - MultiSparseMatrixAssembler(V,U) + MultiSparseMatrixAssembler(M,V,U) +end + +function MultiSparseMatrixAssembler( + testfesps::Vector{<:FESpaceWithDirichletData}, + trialfesps::Vector{<:FESpaceWithDirichletData}) + MultiSparseMatrixAssembler(SparseMatrixCSC,V,U) +end + +function SparseMatrixAssembler( + ::Type{M}, + testfesps::Vector{<:FESpaceWithDirichletData}, + trialfesps::Vector{<:FESpaceWithDirichletData}) where {M<:AbstractSparseMatrix} + MultiSparseMatrixAssembler(M,testfesps,trialfesps) end function SparseMatrixAssembler( testfesps::Vector{<:FESpaceWithDirichletData}, trialfesps::Vector{<:FESpaceWithDirichletData}) - MultiSparseMatrixAssembler(testfesps,trialfesps) + MultiSparseMatrixAssembler(SparseMatrixCSC,testfesps,trialfesps) end function MultiSparseMatrixAssembler( - testfesps::MultiFESpace{E,S}, trialfesps::MultiFESpace{E,S}) where {E,S} + ::Type{M}, + testfesps::MultiFESpace{E,S}, + trialfesps::MultiFESpace{E,S}) where {M<:AbstractSparseMatrix,E,S} @notimplementedif S != ConsequtiveMultiFieldStyle - MultiSparseMatrixAssembler{E}(testfesps,trialfesps) + MultiSparseMatrixAssembler{E,M}(M,testfesps,trialfesps) +end + +function MultiSparseMatrixAssembler( + testfesps::MultiFESpace{E,S}, + trialfesps::MultiFESpace{E,S}) where {E,S} + @notimplementedif S != ConsequtiveMultiFieldStyle + MultiSparseMatrixAssembler(SparseMatrixCSC,testfesps,trialfesps) end function assemble( - this::MultiSparseMatrixAssembler{E}, - allvals::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where E + this::MultiSparseMatrixAssembler{E,M}, + allvals::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {E,M} n = num_free_dofs(this.testfesps) vec = zeros(E,n) @@ -107,8 +146,8 @@ end function assemble!( vec::Vector{E}, - this::MultiSparseMatrixAssembler{E}, - allmcv::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where E + this::MultiSparseMatrixAssembler{E,M}, + allmcv::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {E,M} vec .= zero(E) V = this.testfesps @@ -123,8 +162,8 @@ function assemble!( end function assemble( - this::MultiSparseMatrixAssembler{E}, - allmcm::Vararg{Tuple{<:MultiCellMatrix,<:CellNumber,<:CellNumber}}) where {E} + this::MultiSparseMatrixAssembler{E,M}, + allmcm::Vararg{Tuple{<:MultiCellMatrix,<:CellNumber,<:CellNumber}}) where {E,M} V = this.testfesps U = this.trialfesps @@ -140,21 +179,21 @@ function assemble( mf_vals = apply_constraints_cols(U,mf_vals,cellids_col) mf_cols = reindex(_mf_cols, cellids_col) i_to_fieldid = mf_vals.fieldids - _assemble_mat!( + _assemble_mat!(M, aux_row,aux_col,aux_val,mf_rows,mf_cols,mf_vals, i_to_fieldid,offsets_row,offsets_col) end - sparse(aux_row, aux_col, aux_val) + sparse_from_coo(M,aux_row, aux_col, aux_val) end function assemble!( - mat::SparseMatrixCSC{E}, + mat::AbstractSparseMatrix{E}, this::MultiSparseMatrixAssembler{E}, allvals::Vararg{Tuple{<:MultiCellMatrix,<:CellNumber,<:CellNumber}}) where E # This routine can be optimized a lot taking into a count the sparsity graph of mat # For the moment we create an intermediate matrix and then transfer the nz values m = assemble(this,allvals...) - mat.nzval .= m.nzval + nonzeros(mat) .= nonzeros(m) end function _assemble_vec!(vec,mf_rows,mf_vals,i_to_fieldid,offsets) @@ -172,9 +211,9 @@ function _assemble_vec!(vec,mf_rows,mf_vals,i_to_fieldid,offsets) end end -function _assemble_mat!( +function _assemble_mat!(::Type{M}, aux_row,aux_col,aux_val,mf_rows,mf_cols,mf_vals, - i_to_fieldid,offsets_row,offsets_col) + i_to_fieldid,offsets_row,offsets_col) where {M<:AbstractSparseMatrix} for (mf_rows_c,mf_cols_c,mf_vals_c) in zip(mf_rows,mf_cols,mf_vals) for (i,vals_c) in enumerate(mf_vals_c) @@ -187,9 +226,8 @@ function _assemble_mat!( if gidcol > 0 for (lidrow,gidrow) in enumerate(rows_c) if gidrow > 0 - push!(aux_row, gidrow+offset_row) - push!(aux_col, gidcol+offset_col) - push!(aux_val, vals_c[lidrow,lidcol]) + push_coo!(M,aux_row, aux_col, aux_val, + gidrow+offset_row, gidcol+offset_col, vals_c[lidrow,lidcol]) end end end diff --git a/test/MultiFieldTests/MultiAssemblersTests.jl b/test/MultiFieldTests/MultiAssemblersTests.jl index 15039daf2..f27446a60 100644 --- a/test/MultiFieldTests/MultiAssemblersTests.jl +++ b/test/MultiFieldTests/MultiAssemblersTests.jl @@ -2,77 +2,81 @@ module MultiAssemblersTests using Test using Gridap +using SparseArrays +using SparseMatricesCSR -order = 1 -diritag = "boundary" -model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,3)) -fespace = H1ConformingFESpace(Float64,model,order,diritag) +for T in [SparseMatrixCSC, SparseMatrixCSR, SymSparseMatrixCSR] + order = 1 + diritag = "boundary" + model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,3)) + fespace = H1ConformingFESpace(Float64,model,order,diritag) -ufun1(x) = x[1] + x[2] -U1 = TrialFESpace(fespace,ufun1) + ufun1(x) = x[1] + x[2] + U1 = TrialFESpace(fespace,ufun1) -ufun2(x) = x[1] + x[2] -U2 = TrialFESpace(fespace,ufun2) + ufun2(x) = x[1] + x[2] + U2 = TrialFESpace(fespace,ufun2) -V1 = TestFESpace(fespace) -V2 = V1 + V1 = TestFESpace(fespace) + V2 = V1 -V = [V1,V2] -U = [U1,U2] + V = [V1,V2] + U = [U1,U2] -assem = SparseMatrixAssembler(V,U) + assem = SparseMatrixAssembler(T,V,U) -@test isa(assem,MultiSparseMatrixAssembler) + @test isa(assem,MultiSparseMatrixAssembler) -a(v,u) = varinner(∇(v),∇(u)) + a(v,u) = varinner(∇(v),∇(u)) -m(v,u) = varinner(v,u) + m(v,u) = varinner(v,u) -bfun1(x) = 1.0 -b1(v) = varinner(v,CellField(trian,bfun1)) + bfun1(x) = 1.0 + b1(v) = varinner(v,CellField(trian,bfun1)) -bfun2(x) = 2.0 -b2(v) = varinner(v,CellField(trian,bfun2)) + bfun2(x) = 2.0 + b2(v) = varinner(v,CellField(trian,bfun2)) -trian = Triangulation(model) -quad = CellQuadrature(trian,degree=2) + trian = Triangulation(model) + quad = CellQuadrature(trian,degree=2) -v1 = CellBasis(V1) -v2 = CellBasis(V2) + v1 = CellBasis(V1) + v2 = CellBasis(V2) -u1 = CellBasis(U1) -u2 = CellBasis(U2) + u1 = CellBasis(U1) + u2 = CellBasis(U2) -mat11 = integrate(a(v1,u1),trian,quad) -mat12 = integrate(m(v1,u2),trian,quad) -mat22 = integrate(m(v2,u2),trian,quad) + mat11 = integrate(a(v1,u1),trian,quad) + mat12 = integrate(m(v1,u2),trian,quad) + mat22 = integrate(m(v2,u2),trian,quad) -mat = MultiCellMatrix([mat11,mat12,mat22],[(1,1),(1,2),(2,2)]) + mat = MultiCellMatrix([mat11,mat12,mat22],[(1,1),(1,2),(2,2)]) -vec1 = integrate(b1(v1),trian,quad) -vec2 = integrate(b2(v2),trian,quad) + vec1 = integrate(b1(v1),trian,quad) + vec2 = integrate(b2(v2),trian,quad) -vec = MultiCellVector([vec1,vec2],[(1,),(2,)]) + vec = MultiCellVector([vec1,vec2],[(1,),(2,)]) -v = assemble(assem,vec) -@test v ≈ [0.166666666666, 0.166666666666, 0.3333333333333, 0.333333333333] + v = assemble(assem,vec) + @test v ≈ [0.166666666666, 0.166666666666, 0.3333333333333, 0.333333333333] -mm = assemble(assem,mat) + mm = assemble(assem,mat) -assemble!(mm,assem,mat) + assemble!(mm,assem,mat) -cellids = IdentityCellNumber(Int,length(vec)) + cellids = IdentityCellNumber(Int,length(vec)) -v2 = assemble(assem,(vec,cellids),(vec,cellids)) -@test v2 ≈ 2*v + v2 = assemble(assem,(vec,cellids),(vec,cellids)) + @test v2 ≈ 2*v -assemble!(v2,assem,(vec,cellids),(vec,cellids)) -@test v2 ≈ 2*v + assemble!(v2,assem,(vec,cellids),(vec,cellids)) + @test v2 ≈ 2*v -mm2 = assemble(assem,(mat,cellids,cellids),(mat,cellids,cellids)) -@test mm2 ≈ 2*mm + mm2 = assemble(assem,(mat,cellids,cellids),(mat,cellids,cellids)) + @test mm2 ≈ 2*mm -assemble!(mm2,assem,(mat,cellids,cellids),(mat,cellids,cellids)) -@test mm2 ≈ 2*mm + assemble!(mm2,assem,(mat,cellids,cellids),(mat,cellids,cellids)) + @test mm2 ≈ 2*mm +end end # module MultiAssemblersTests From 53d940ed44ebc072add8c1ab444aefa88a970fa7 Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Mon, 28 Oct 2019 16:20:30 +0100 Subject: [PATCH 4/7] Solve review comments in #118 --- src/FESpaces/Assemblers.jl | 11 +++++------ src/MultiField/MultiAssemblers.jl | 14 +++++++------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index 17ffa85b4..c295b82db 100644 --- a/src/FESpaces/Assemblers.jl +++ b/src/FESpaces/Assemblers.jl @@ -11,7 +11,6 @@ export Assembler export SparseMatrixAssembler export assemble export assemble! -export push_coo! export sparse_from_coo """ @@ -85,7 +84,7 @@ end """ Assembler that produces SparseMatrices from the SparseArrays and SparseMatricesCSR packages """ -struct SparseMatrixAssembler{E,M} <: Assembler{AbstractSparseMatrix{E,Int},Vector{E}} +struct SparseMatrixAssembler{E,M} <: Assembler{M,Vector{E}} testfesp::FESpace trialfesp::FESpace @@ -103,8 +102,8 @@ function SparseMatrixAssembler(test::FESpace{D,Z,T}, trial::FESpace{D,Z,T}) wher end function assemble( - this::SparseMatrixAssembler{E,M}, - vals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where {E,M} + this::SparseMatrixAssembler{E}, + vals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where {E} n = num_free_dofs(this.testfesp) vec = zeros(E,n) @@ -114,8 +113,8 @@ end function assemble!( vec::Vector{E}, - this::SparseMatrixAssembler{E,M}, - allvals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where {E,M} + this::SparseMatrixAssembler{E}, + allvals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where {E} vec .= zero(E) rows = celldofids(this.testfesp) diff --git a/src/MultiField/MultiAssemblers.jl b/src/MultiField/MultiAssemblers.jl index 76c87139d..0f46b211f 100644 --- a/src/MultiField/MultiAssemblers.jl +++ b/src/MultiField/MultiAssemblers.jl @@ -5,6 +5,7 @@ using Gridap.Helpers using Gridap.MultiFESpaces: _compute_offsets using SparseArrays +using SparseMatricesCSR export MultiAssembler export MultiSparseMatrixAssembler @@ -13,8 +14,7 @@ import Gridap.Assemblers: assemble import Gridap.Assemblers: assemble! import Gridap.Assemblers: SparseMatrixAssembler import Gridap.Assemblers: SparseMatrixAssembler -import Gridap.Assemblers: push_coo! -import Gridap.Assemblers: sparse_from_coo +using Gridap.Assemblers: sparse_from_coo abstract type MultiAssembler{M<:AbstractMatrix,V<:AbstractVector} end @@ -73,7 +73,7 @@ end """ MultiAssembler that produces SparseMatrices from the SparseArrays package """ -struct MultiSparseMatrixAssembler{E,M} <: MultiAssembler{AbstractSparseMatrix{E,Int},Vector{E}} +struct MultiSparseMatrixAssembler{E,M} <: MultiAssembler{M,Vector{E}} testfesps::MultiFESpace{E} trialfesps::MultiFESpace{E} @@ -135,8 +135,8 @@ function MultiSparseMatrixAssembler( end function assemble( - this::MultiSparseMatrixAssembler{E,M}, - allvals::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {E,M} + this::MultiSparseMatrixAssembler{E}, + allvals::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {E} n = num_free_dofs(this.testfesps) vec = zeros(E,n) @@ -146,8 +146,8 @@ end function assemble!( vec::Vector{E}, - this::MultiSparseMatrixAssembler{E,M}, - allmcv::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {E,M} + this::MultiSparseMatrixAssembler{E}, + allmcv::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {E} vec .= zero(E) V = this.testfesps From ee467fec6fd37ccd6e76a6b55ba7d9e00c86e16b Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Mon, 28 Oct 2019 16:23:33 +0100 Subject: [PATCH 5/7] Solve review comments in #118 --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index fcad92900..1952980f7 100644 --- a/Project.toml +++ b/Project.toml @@ -29,6 +29,7 @@ LineSearches = "7.0.1" NLsolve = "4.1.0" QuadGK = "2.1.0" Reexport = "0.2.0" +SparseMatricesCSR = "0.3.1" StaticArrays = "0.10.3" TensorPolynomialBases = "0.1.4" TensorValues = "0.3.5" From 6ef127f127180c0b4ae205ee9b92e1558e8fc83f Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Mon, 28 Oct 2019 17:44:14 +0100 Subject: [PATCH 6/7] Add finalize_coo! call in assemble funcion for Assemblers and MultiAssemblers --- src/FESpaces/Assemblers.jl | 3 +++ src/MultiField/MultiAssemblers.jl | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/FESpaces/Assemblers.jl b/src/FESpaces/Assemblers.jl index c295b82db..605a5611a 100644 --- a/src/FESpaces/Assemblers.jl +++ b/src/FESpaces/Assemblers.jl @@ -152,6 +152,9 @@ function assemble( _assemble_sparse_matrix_values!(M, aux_row,aux_col,aux_val,vals_m,rows_m,cols_m) end + num_rows = num_free_dofs(this.testfesp) + num_cols = num_free_dofs(this.trialfesp) + finalize_coo!(M,aux_row,aux_col,aux_val,num_rows,num_cols) sparse_from_coo(M,aux_row,aux_col,aux_val) end diff --git a/src/MultiField/MultiAssemblers.jl b/src/MultiField/MultiAssemblers.jl index 0f46b211f..d3f774139 100644 --- a/src/MultiField/MultiAssemblers.jl +++ b/src/MultiField/MultiAssemblers.jl @@ -183,6 +183,9 @@ function assemble( aux_row,aux_col,aux_val,mf_rows,mf_cols,mf_vals, i_to_fieldid,offsets_row,offsets_col) end + num_rows = num_free_dofs(this.testfesps) + num_cols = num_free_dofs(this.trialfesps) + finalize_coo!(M,aux_row,aux_col,aux_val,num_rows,num_cols) sparse_from_coo(M,aux_row, aux_col, aux_val) end From 08edf7ed10d49cc5a140af7d5b71f0c0453f23f7 Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Tue, 29 Oct 2019 11:57:45 +0100 Subject: [PATCH 7/7] Add new signatures for LinearFEOperator and NonLinearFEOperator functions accepting the the matrix template type --- src/FESpaces/FEOperators.jl | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/FESpaces/FEOperators.jl b/src/FESpaces/FEOperators.jl index 568e6075a..129b6c5fd 100644 --- a/src/FESpaces/FEOperators.jl +++ b/src/FESpaces/FEOperators.jl @@ -4,6 +4,7 @@ using Gridap using Gridap.Helpers using LinearAlgebra +using SparseArrays import Gridap: apply export apply! @@ -187,6 +188,16 @@ function LinearFEOperator( LinearFEOperator(testfesp,trialfesp,assem,terms...) end +function LinearFEOperator( + ::Type{M}, + testfesp::FESpaceLike, + trialfesp::FESpaceLike, + terms::Vararg{<:AffineFETerm}) where {M} + + assem = SparseMatrixAssembler(M,testfesp,trialfesp) + LinearFEOperator(testfesp,trialfesp,assem,terms...) +end + function LinearFEOperator( biform::Function, liform::Function, @@ -290,6 +301,17 @@ function NonLinearFEOperator( NonLinearFEOperator(testfesp,trialfesp,assem,terms) end +function NonLinearFEOperator( + ::Type{M}, + testfesp::FESpaceLike, + trialfesp::FESpaceLike, + terms::Vararg{<:FETerm}) where {M} + + assem = SparseMatrixAssembler(M,testfesp,trialfesp) + NonLinearFEOperator(testfesp,trialfesp,assem,terms...) +end + + function NonLinearFEOperator( res::Function, jac::Function,