Merge pull request #247 from gridap/add_fe_space_with_linear_constraints

Added FESpaceWithLinearConstraints

src/Arrays/IdentityVectors.jl

@@ -29,3 +29,7 @@ end
 function reindex(a::AppliedArray,b::IdentityVector)
   a
 end
+
+function reindex(a::Fill,b::IdentityVector)
+  a
+end

src/Arrays/Tables.jl

@@ -442,3 +442,8 @@ function from_dict(::Type{Table{T,P}}, dict::Dict{Symbol,Any}) where {T,P}
   ptrs::Vector{P} = dict[:ptrs]
   Table(data,ptrs)
 end
+
+function Base.copy(a::Table)
+  Table(copy(a.data),copy(a.ptrs))
+end
+

src/FESpaces/CellDofBases.jl

@@ -53,3 +53,12 @@ end
 RefStyle(::Type{<:GenericCellDofBasis{R}}) where R = Val{R}()
 
 get_array(a::GenericCellDofBasis) = a.array
+
+function reindex(cf::CellDofBasis,a::AbstractVector)
+  similar_object(cf,reindex(get_array(cf),a))
+end
+
+function similar_object(cf::CellDofBasis,a::AbstractArray)
+  GenericCellDofBasis(RefStyle(cf),a)
+end
+

src/FESpaces/FESpaces.jl

@@ -89,6 +89,8 @@ export get_cell_basis
 export zero_free_values
 export constraint_style
 export has_constraints
+export get_cell_isconstrained
+export get_cell_constraints
 export get_constraint_kernel_matrix_cols
 export get_constraint_kernel_matrix_rows
 export get_constraint_kernel_vector
@@ -202,6 +204,7 @@ export CLagrangianFESpace
 export DivConformingFESpace
 export CurlConformingFESpace
 export DirichletFESpace
+export FESpaceWithLinearConstraints
 export ExtendedFESpace
 
 export @law
@@ -269,6 +272,8 @@ include("DirichletFESpaces.jl")
 
 include("ExtendedFESpaces.jl")
 
+include("FESpacesWithLinearConstraints.jl")
+
 include("FESpaceFactories.jl")
 
 include("StateLaws.jl")

src/FESpaces/FESpacesInterfaces.jl

@@ -67,21 +67,39 @@ end
 
 constraint_style(f::T) where T<:FESpace = constraint_style(T)
 
+function get_cell_isconstrained(f::FESpace)
+  @abstractmethod
+end
+
+function get_cell_constraints(f::FESpace)
+  @abstractmethod
+end
+
 """
 """
 function get_constraint_kernel_matrix_cols(f::FESpace)
-  @abstractmethod
+  _default_constraint_kernel(f,constraint_style(f))
 end
 
 """
 """
 function get_constraint_kernel_matrix_rows(f::FESpace)
-  @abstractmethod
+  _default_constraint_kernel(f,constraint_style(f))
 end
 
 """
 """
 function get_constraint_kernel_vector(f::FESpace)
+  _default_constraint_kernel(f,constraint_style(f))
+end
+
+function _default_constraint_kernel(f,::Val{false})
+  function default_kernel(m,isconstr,constr)
+    m
+  end
+end
+
+function _default_constraint_kernel(f,::Val{true})
   @abstractmethod
 end
 
@@ -150,8 +168,10 @@ function _apply_constraints_matrix_cols(::Val{false},f,cellmat,cellids)
 end
 
 function _apply_constraints_matrix_cols(::Val{true},f,cellmat,cellids)
+  cell_to_isconstr = reindex(get_cell_isconstrained(f),cellids)
+  cell_to_constr = reindex(get_cell_constraints(f),cellids)
   k = get_constraint_kernel_matrix_cols(f)
-  apply(k,cellmat,cellids)
+  apply(k,cellmat,cell_to_isconstr,cell_to_constr)
 end
 
 """
@@ -165,8 +185,10 @@ function _apply_constraints_matrix_rows(::Val{false},f,cellmat,cellids)
 end
 
 function _apply_constraints_matrix_rows(::Val{true},f,cellmat,cellids)
+  cell_to_isconstr = reindex(get_cell_isconstrained(f),cellids)
+  cell_to_constr = reindex(get_cell_constraints(f),cellids)
   k = get_constraint_kernel_matrix_rows(f)
-  apply(k,cellmat,cellids)
+  apply(k,cellmat,cell_to_isconstr,cell_to_constr)
 end
 
 """
@@ -180,8 +202,10 @@ function _apply_constraints_vector(::Val{false},f,cellvec,cellids)
 end
 
 function _apply_constraints_vector(::Val{true},f,cellvec,cellids)
+  cell_to_isconstr = reindex(get_cell_isconstrained(f),cellids)
+  cell_to_constr = reindex(get_cell_constraints(f),cellids)
   k = get_constraint_kernel_vector(f)
-  apply(k,cellvec,cellids)
+  apply(k,cellvec,cell_to_isconstr,cell_to_constr)
 end
 
 """
@@ -195,9 +219,11 @@ function _apply_constraints_matrix_and_vector_cols(::Val{false},f,cellmatvec,cel
 end
 
 function _apply_constraints_matrix_and_vector_cols(::Val{true},f,cellmatvec,cellids)
+  cell_to_isconstr = reindex(get_cell_isconstrained(f),cellids)
+  cell_to_constr = reindex(get_cell_constraints(f),cellids)
   kmat = get_constraint_kernel_matrix_cols(f)
   k = MatKernel(kmat)
-  apply(k,cellmatvec,cellids)
+  apply(k,cellmatvec,cell_to_isconstr,cell_to_constr)
 end
 
 """
@@ -211,10 +237,12 @@ function _apply_constraints_matrix_and_vector_rows(::Val{false},f,cellmatvec,cel
 end
 
 function _apply_constraints_matrix_and_vector_rows(::Val{true},f,cellmatvec,cellids)
+  cell_to_isconstr = reindex(get_cell_isconstrained(f),cellids)
+  cell_to_constr = reindex(get_cell_constraints(f),cellids)
   kmat = get_constraint_kernel_matrix_rows(f)
   kvec = get_constraint_kernel_vector(f)
   k = MatVecKernel(kmat,kvec)
-  apply(k,cellmatvec,cellids)
+  apply(k,cellmatvec,cell_to_isconstr,cell_to_constr)
 end
 
 # Helpers
@@ -224,47 +252,47 @@ struct MatVecKernel{A<:Kernel,B<:Kernel} <: Kernel
   kvec::B
 end
 
-function kernel_cache(k::MatVecKernel,matvec,cellid)
+function kernel_cache(k::MatVecKernel,matvec,isconstr,constr)
   mat, vec = matvec
-  cmat = kernel_cache(k.kmat,mat,cellid)
-  cvec = kernel_cache(k.kvec,vec,cellid)
+  cmat = kernel_cache(k.kmat,mat,isconstr,constr)
+  cvec = kernel_cache(k.kvec,vec,isconstr,constr)
   (cmat, cvec)
 end
 
-function kernel_return_type(k::MatVecKernel,matvec,cellid)
+function kernel_return_type(k::MatVecKernel,matvec,isconstr,constr)
   mat, vec = matvec
-  A = kernel_return_type(k.kmat,mat,cellid)
-  B = kernel_return_type(k.kvec,vec,cellid)
+  A = kernel_return_type(k.kmat,mat,isconstr,constr)
+  B = kernel_return_type(k.kvec,vec,isconstr,constr)
   Tuple{A,B}
 end
 
-@inline function apply_kernel!(cache,k::MatVecKernel,matvec,cellid)
+@inline function apply_kernel!(cache,k::MatVecKernel,matvec,isconstr,constr)
   cmat, cvec = cache
   mat, vec = matvec
-  a = apply_kernel!(cmat,k.kmat,mat,cellid)
-  b = apply_kernel!(cvec,k.kvec,vec,cellid)
+  a = apply_kernel!(cmat,k.kmat,mat,isconstr,constr)
+  b = apply_kernel!(cvec,k.kvec,vec,isconstr,constr)
   (a,b)
 end
 
 struct MatKernel{A<:Kernel} <: Kernel
   kmat::A
 end
 
-function kernel_cache(k::MatKernel,matvec,cellid)
+function kernel_cache(k::MatKernel,matvec,isconstr,constr)
   mat, vec = matvec
-  cmat = kernel_cache(k.kmat,mat,cellid)
+  cmat = kernel_cache(k.kmat,mat,isconstr,constr)
   cmat
 end
 
-function kernel_return_type(k::MatKernel,matvec,cellid)
+function kernel_return_type(k::MatKernel,matvec,isconstr,constr)
   mat, vec = matvec
-  A = kernel_return_type(k.kmat,mat,cellid)
+  A = kernel_return_type(k.kmat,mat,isconstr,constr)
   B = typeof(vec)
   Tuple{A,B}
 end
 
-@inline function apply_kernel!(cmat,k::MatKernel,matvec,cellid)
+@inline function apply_kernel!(cmat,k::MatKernel,matvec,isconstr,constr)
   mat, vec = matvec
-  a = apply_kernel!(cmat,k.kmat,mat,cellid)
+  a = apply_kernel!(cmat,k.kmat,mat,isconstr,constr)
   (a,vec)
 end