Extended assembly of matrices to skeleton terms

Author: fverdugo

src/FESpaces/Assemblers.jl

@@ -31,7 +31,7 @@ Assembly of a matrix allocating output
 """
 function assemble(
   ::Assembler{M,V},
-  ::Vararg{Tuple{<:CellMatrix,<:CellNumber}})::M where {M,V}
+  ::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}})::M where {M,V}
   @abstractmethod
 end
 
@@ -51,22 +51,34 @@ In-place assembly of a matrix (allows a LOT of optimizations)
 function assemble!(
   ::M,
   ::Assembler{M,V},
-  ::Vararg{Tuple{<:CellMatrix,<:CellNumber}})::M where {M,V}
+  ::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}})::M where {M,V}
   @abstractmethod
 end
 
-function assemble(a::Assembler,cv::CellArray)
+function assemble(a::Assembler,cv::CellVector)
   l = length(cv)
   ide = IdentityCellNumber(Int,l)
   assemble(a,(cv,ide))
 end
 
-function assemble!(r,a::Assembler,cv::CellArray)
+function assemble(a::Assembler,cv::CellMatrix)
+  l = length(cv)
+  ide = IdentityCellNumber(Int,l)
+  assemble(a,(cv,ide,ide))
+end
+
+function assemble!(r,a::Assembler,cv::CellVector)
   l = length(cv)
   ide = IdentityCellNumber(Int,l)
   assemble!(r,a,(cv,ide))
 end
 
+function assemble!(r,a::Assembler,cv::CellMatrix)
+  l = length(cv)
+  ide = IdentityCellNumber(Int,l)
+  assemble!(r,a,(cv,ide,ide))
+end
+
 """
 Assembler that produces SparseMatrices from the SparseArrays package
 """
@@ -116,19 +128,19 @@ end
 
 function assemble(
   this::SparseMatrixAssembler{E},
-  allvals::Vararg{Tuple{<:CellMatrix,<:CellNumber}}) where E
+  allvals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where E
 
   I = Int
   aux_row = I[]; aux_col = I[]; aux_val = E[]
 
   _rows_m = celldofids(this.testfesp)
   _cols_m = celldofids(this.trialfesp)
 
-  for (vals,cellids) in allvals
-    _vals = apply_constraints_rows(this.testfesp, vals, cellids)
-    rows_m = reindex(_rows_m, cellids)
-    vals_m = apply_constraints_cols(this.trialfesp, _vals, cellids)
-    cols_m = reindex(_cols_m, cellids)
+  for (vals, cellids_row, cellids_col) in allvals
+    _vals = apply_constraints_rows(this.testfesp, vals, cellids_row)
+    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!(
       aux_row,aux_col,aux_val,vals_m,rows_m,cols_m)
   end
@@ -154,7 +166,7 @@ end
 function assemble!(
   mat::SparseMatrixCSC{E},
   this::SparseMatrixAssembler{E},
-  vals::Vararg{Tuple{<:CellMatrix,<:CellNumber}}) where E
+  vals::Vararg{Tuple{<:CellMatrix,<: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,vals...)

src/FESpaces/FETerms.jl

@@ -222,19 +222,21 @@ function _append_vector_contribution!(m,r,c)
   push!(m,(r,c))
 end
 
-#function _append_vector_contribution!(m,r::SkeletonCellVector,c)
-#  push!(m,(r.cellvector1,c[1]))
-#  push!(m,(r.cellvector2,c[2]))
-#end
+function _append_vector_contribution!(m,r::SkeletonCellVector,c)
+  push!(m,(r.cellvector1,c[1]))
+  push!(m,(r.cellvector2,c[2]))
+end
 
 function _append_matrix_contribution!(m,r,c)
-  push!(m,(r,c))
+  push!(m,(r,c,c))
 end
 
-#function _append_matrix_contribution!(m,r::SkeletonCellMatrix,c)
-#  push!(m,(r.cellmatrix11,c[1]))
-#  push!(m,(r.cellmatrix12,c[2]))
-#end
+function _append_matrix_contribution!(m,r::SkeletonCellMatrix,c)
+  push!(m,(r.cellmatrix11,c[1],c[1]))
+  push!(m,(r.cellmatrix12,c[1],c[2]))
+  push!(m,(r.cellmatrix21,c[2],c[1]))
+  push!(m,(r.cellmatrix22,c[2],c[2]))
+end
 
 const _jac = setup_cell_jacobian
 

src/Integration/SkeletonCellFields.jl

@@ -2,6 +2,11 @@ module SkeletonCellFields
 
 using Gridap
 
+export SkeletonPair
+export SkeletonCellBasis
+export SkeletonCellVector
+export SkeletonCellMatrix
+
 export jump
 import Gridap: restrict
 import Gridap: gradient

test/FESpacesTests/AssemblersTests.jl

@@ -64,14 +64,14 @@ vec3 = assemble(assem, (bvec,cellids), (bvec,cellids))
 assemble!(vec3, assem, (bvec,cellids), (bvec,cellids))
 @test vec3 ≈ 2*vec
 
-mat2 = assemble(assem, (mmat,cellids) )
+mat2 = assemble(assem, (mmat,cellids,cellids) )
 
 @test mat2 == mat
 
-mat3 = assemble(assem, (mmat,cellids), (mmat,cellids))
+mat3 = assemble(assem, (mmat,cellids,cellids), (mmat,cellids,cellids))
 @test mat3 ≈ 2*mat
 
-assemble!(mat3, assem, (mmat,cellids), (mmat,cellids))
+assemble!(mat3, assem, (mmat,cellids,cellids), (mmat,cellids,cellids))
 @test mat3 ≈ 2*mat
 
 end # module AssemblersTests

test/FESpacesTests/FETermsTests.jl

@@ -152,19 +152,19 @@ assem = SparseMatrixAssembler(V,U)
 
 cms = setup_cell_jacobian(uh,v,du,t1)
 
-cms = setup_cell_jacobian(uh,v,du,t1,t2)
+cms = setup_cell_jacobian(uh,v,du,t1,t2,t6)
 
 mat = assemble(assem,cms...)
 
-cvs = setup_cell_residual(uh,v,t1,t2)
+cvs = setup_cell_residual(uh,v,t1,t2,t6)
 
 vec = assemble(assem,cvs...)
 
-cvs = setup_cell_vector(v,uhd,t1,t2)
+cvs = setup_cell_vector(v,uhd,t1,t2,t6)
 
 vec = assemble(assem,cvs...)
 
-cms = setup_cell_matrix(v,du,t1,t2)
+cms = setup_cell_matrix(v,du,t1,t2,t6)
 
 mat = assemble(assem,cms...)