Merge pull request #43 from gridap/extend_assembler_to_multiple_fields

This PR extends the FE assembly routines to support weak forms with several terms, which are integrated in different domains. E.g., problems with Neumann BCs.

This PR addresses issue #42

NEWS.md

@@ -11,8 +11,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `NonIterableCellMap`, a cell map that has iteration intentionally disabled. Available since commit [956a537](https://github.com/gridap/Gridap.jl/commit/956a5374db6c3b9546e85e0d4d49ae0560057565).
 - `BoundaryTriangulation` an integration mesh used to integrate `CellField` and `CellBasis` objects restricted on a surface. Available since commit [e981f3c](https://github.com/gridap/Gridap.jl/commit/e981f3c221f3624cfc6764efa47f22652fc22b4f)
 - Function `restrict` for restricting `CellField` and `CellBasis` objects to surfaces. Available since commit [e981f3c](https://github.com/gridap/Gridap.jl/commit/e981f3c221f3624cfc6764efa47f22652fc22b4f)
+- `IdentityCellNumber`, an indexable cell number that simply returns the given index. Also efficient implementation of `reindex` for this type (i.e. do nothing). Available since commit [b6b4c32](https://github.com/gridap/Gridap.jl/commit/b6b4c32c8c4b826a41ba64c770ac8a1c394e16f0)
 
 ### Changed
+- Changed the signature of `assemble`, `apply_constraints`, `apply_constraints_rows`, and `apply_constraints_cols` to support FE assembly of several terms, which  are integrated in different domains. The old API of `asseble` is still functional, but not for the `apply_constraints` et al.
+
+
 ### Removed
 ### Deprecated
 ### Fixed

src/CellValues/IdentityCellNumbers.jl

@@ -0,0 +1,32 @@
+module IdentityCellNumbers
+
+using Gridap
+
+export IdentityCellNumber
+import Gridap: reindex
+import Base: size
+import Base: getindex
+
+struct IdentityCellNumber{T} <: IndexCellValue{T,1}
+  length::Int
+end
+
+function IdentityCellNumber(::Type{T},l::Integer) where T <: Integer
+  IdentityCellNumber{T}(l)
+end
+
+function getindex(c::IdentityCellNumber{T},i::Integer) where T
+  @assert i > 0
+  @assert i <= c.length
+  j::T = i
+  j
+end
+
+size(c::IdentityCellNumber) = (c.length,)
+
+function reindex(values::IndexCellValue, indices::IdentityCellNumber)
+  @assert length(values) == length(indices)
+  values
+end
+
+end # module

src/CellValues/files.jl

@@ -53,4 +53,7 @@ include("CompressedCellValues.jl")
 include("NonIterableCellMaps.jl")
 @reexport using Gridap.NonIterableCellMaps
 
+include("IdentityCellNumbers.jl")
+@reexport using Gridap.IdentityCellNumbers
+
 

src/FESpaces/Assemblers.jl

@@ -11,10 +11,6 @@ export SparseMatrixAssembler
 export assemble
 export assemble!
 
-#@fverdugo for the moment the abstract interface of Assembler
-# (and therefore its concrete implementations)
-# assumes single field, and single term
-
 """
 Abstract assembly operator
 Parametrized by the type of returned matrix and vector
@@ -24,31 +20,53 @@ abstract type Assembler{M<:AbstractMatrix,V<:AbstractVector} end
 """
 Assembly of a vector allocating output
 """
-function assemble(::Assembler{M,V},::CellVector)::V where {M,V}
+function assemble(
+  ::Assembler{M,V},
+  ::Vararg{Tuple{<:CellVector,<:CellNumber}})::V where {M,V}
   @abstractmethod
 end
 
 """
 Assembly of a matrix allocating output
 """
-function assemble(::Assembler{M,V},::CellMatrix)::M where {M,V}
+function assemble(
+  ::Assembler{M,V},
+  ::Vararg{Tuple{<:CellMatrix,<:CellNumber}})::M where {M,V}
   @abstractmethod
 end
 
 """
 In-place assembly of a vector (allows optimizations)
 """
-function assemble!(::V,::Assembler{M,V},::CellVector)::V where {M,V}
+function assemble!(
+  ::V,
+  ::Assembler{M,V},
+  ::Vararg{Tuple{<:CellVector,<:CellNumber}})::V where {M,V}
   @abstractmethod
 end
 
 """
 In-place assembly of a matrix (allows a LOT of optimizations)
 """
-function assemble!(::M,::Assembler{M,V},::CellMatrix)::M where {M,V}
+function assemble!(
+  ::M,
+  ::Assembler{M,V},
+  ::Vararg{Tuple{<:CellMatrix,<:CellNumber}})::M where {M,V}
   @abstractmethod
 end
 
+function assemble(a::Assembler,cv::CellArray)
+  l = length(cv)
+  ide = IdentityCellNumber(Int,l)
+  assemble(a,(cv,ide))
+end
+
+function assemble!(r,a::Assembler,cv::CellArray)
+  l = length(cv)
+  ide = IdentityCellNumber(Int,l)
+  assemble!(r,a,(cv,ide))
+end
+
 """
 Assembler that produces SparseMatrices from the SparseArrays package
 """
@@ -62,19 +80,28 @@ function SparseMatrixAssembler(test::FESpace{D,Z,T}, trial::FESpace{D,Z,T}) wher
   SparseMatrixAssembler{E}(trial,test)
 end
 
-function assemble(this::SparseMatrixAssembler{E}, vals::CellVector) where E
+function assemble(
+  this::SparseMatrixAssembler{E},
+  vals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where E
+
   n = num_free_dofs(this.testfesp)
   vec = zeros(E,n)
-  assemble!(vec,this,vals)
+  assemble!(vec,this,vals...)
   vec
 end
 
 function assemble!(
-  vec::Vector{E},this::SparseMatrixAssembler{E}, vals::CellVector) where E
-  _vals = apply_constraints(this.testfesp, vals)
-  _rows = celldofids(this.testfesp)
+  vec::Vector{E},
+  this::SparseMatrixAssembler{E},
+  allvals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where E
+
   vec .= zero(E)
-  _assemble_vector!(vec, _vals, _rows)
+  rows = celldofids(this.testfesp)
+  for (vals, cellids) in allvals
+    _vals = apply_constraints(this.testfesp, vals, cellids)
+    _rows = reindex(rows,cellids)
+    _assemble_vector!(vec, _vals, _rows)
+  end
 end
 
 function _assemble_vector!(vec,vals,rows)
@@ -87,17 +114,28 @@ function _assemble_vector!(vec,vals,rows)
   end
 end
 
-function assemble(this::SparseMatrixAssembler{E}, vals::CellMatrix) where E
-  _vals = apply_constraints_rows(this.testfesp, vals)
-  rows_m = celldofids(this.testfesp)
-  _vals = apply_constraints_cols(this.trialfesp, _vals)
-  cols_m = celldofids(this.trialfesp)
-  args = _assemble_sparse_matrix_values(_vals,rows_m,cols_m,Int,E)
-  sparse(args...)
-end
+function assemble(
+  this::SparseMatrixAssembler{E},
+  allvals::Vararg{Tuple{<:CellMatrix,<:CellNumber}}) where E
 
-function _assemble_sparse_matrix_values(vals,rows,cols,I,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)
+    _assemble_sparse_matrix_values!(
+      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)
   for (rows_c, cols_c, vals_c) in zip(rows,cols,vals)
      for (j,gidcol) in enumerate(cols_c)
        if gidcol > 0
@@ -111,44 +149,16 @@ function _assemble_sparse_matrix_values(vals,rows,cols,I,E)
       end
     end
   end
-  (aux_row, aux_col, aux_val)
 end
 
 function assemble!(
-  mat::SparseMatrixCSC{E}, this::SparseMatrixAssembler{E}, vals::CellMatrix) where E
+  mat::SparseMatrixCSC{E},
+  this::SparseMatrixAssembler{E},
+  vals::Vararg{Tuple{<:CellMatrix,<: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)
+  m = assemble(this,vals...)
   mat.nzval .= m.nzval
 end
 
-# Draft of multi field assembler
-#function _assemble_sparse_matrix_values(mf_vals,mf_rows,mf_cols,I,E)
-#  aux_row = I[]; aux_col = I[]; aux_val = E[]
-#  for (mf_rows_c, mf_cols_c, mf_vals_c) in zip(mf_rows,mf_cols,mf_vals)
-#    for (vals_c, (ifield, jfield)) in eachblock(mf_vals_c)
-#      rows_c = mf_rows_c[ifield]
-#      cols_c = mf_cols_c[jfield]
-#      row_offset = row_offsets[ifield]
-#      col_offset = col_offsets[jfield]
-#      _asseble_cell_values!(aux_row,aux_col,aux_val,rows_c,cols_c,vals_c,col_offset,row_offset)
-#    end
-#  end
-#  (aux_row, aux_col, aux_val)
-#end
-#
-#function _asseble_cell_values!(aux_row,aux_col,aux_val,rows_c,cols_c,vals_c,col_offset,row_offset)
-# for (j,gidcol) in enumerate(cols_c)
-#   if gidcol > 0
-#     for (i,gidrow) in enumerate(rows_c)
-#       if gidrow > 0
-#         push!(aux_row, gidrow+row_offset)
-#         push!(aux_col, gidcol+col_offset)
-#         push!(aux_val, vals_c[i,j])
-#       end
-#     end
-#   end
-# end
-#end
-
-end # module Assemblers
+end # module

src/FESpaces/FESpaces.jl

@@ -51,15 +51,15 @@ num_diri_dofs(::FESpace)::Int = @abstractmethod
 
 diri_tags(::FESpace)::Vector{Int} = @abstractmethod
 
-function apply_constraints(::FESpace, cellvec::CellVector)::CellVector
+function apply_constraints(::FESpace, cellvec::CellVector, cellids::CellNumber)::CellVector
   @abstractmethod
 end
 
-function apply_constraints_rows(::FESpace, cellmat::CellMatrix)::CellMatrix
+function apply_constraints_rows(::FESpace, cellmat::CellMatrix, cellids::CellNumber)::CellMatrix
   @abstractmethod
 end
 
-function apply_constraints_cols(::FESpace, cellmat::CellMatrix)::CellMatrix
+function apply_constraints_cols(::FESpace, cellmat::CellMatrix, cellids::CellNumber)::CellMatrix
   @abstractmethod
 end
 
@@ -279,18 +279,18 @@ num_diri_dofs(f::FESpaceWithDirichletData) = num_diri_dofs(f.fespace)
 diri_tags(f::FESpaceWithDirichletData) = diri_tags(f.fespace)
 
 function apply_constraints(
-  f::FESpaceWithDirichletData, cellvec::CellVector)
-  apply_constraints(f.fespace,cellvec)
+  f::FESpaceWithDirichletData, cellvec::CellVector, cellids::CellNumber)
+  apply_constraints(f.fespace,cellvec,cellids)
 end
 
 function apply_constraints_rows(
-  f::FESpaceWithDirichletData, cellmat::CellMatrix)
-  apply_constraints_rows(f.fespace,cellmat)
+  f::FESpaceWithDirichletData, cellmat::CellMatrix, cellids::CellNumber)
+  apply_constraints_rows(f.fespace,cellmat,cellids)
 end
 
 function apply_constraints_cols(
-  f::FESpaceWithDirichletData, cellmat::CellMatrix)
-  apply_constraints_cols(f.fespace,cellmat)
+  f::FESpaceWithDirichletData, cellmat::CellMatrix, cellids::CellNumber)
+  apply_constraints_cols(f.fespace,cellmat,cellids)
 end
 
 function celldofids(f::FESpaceWithDirichletData)
@@ -377,17 +377,17 @@ num_diri_dofs(this::ConformingFESpace) = this.num_diri_dofs
 diri_tags(f::ConformingFESpace) = f.diri_tags
 
 function apply_constraints(
-  this::ConformingFESpace, cellvec::CellVector)
+  this::ConformingFESpace, cellvec::CellVector, cellids::CellNumber)
   cellvec
 end
 
 function apply_constraints_rows(
-  this::ConformingFESpace, cellmat::CellMatrix)
+  this::ConformingFESpace, cellmat::CellMatrix, cellids::CellNumber)
   cellmat
 end
 
 function apply_constraints_cols(
-  this::ConformingFESpace, cellmat::CellMatrix)
+  this::ConformingFESpace, cellmat::CellMatrix, cellids::CellNumber)
   cellmat
 end