Starting to draft primal HHO

Project.toml

@@ -16,4 +16,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 DrWatson = "1"
-Gridap = "0.17"
+Gridap = "0.18"

src/AddNaiveInnerMostBlockLevelMap.jl

@@ -9,18 +9,18 @@ function Gridap.Arrays.return_cache(
   ai=testitem(a.array)
 
   # 2. Generate cache of k applied to ai
-  ci=return_cache(k,ai)
+  ci=Gridap.Arrays.return_cache(k,ai)
 
   # 3. Generate array of caches
   cache_array=Array{typeof(ci),length(size(a.array))}(undef,size(a.array))
   for i in eachindex(a.array)
     if a.touched[i]
-      cache_array[i]=return_cache(k,a.array[i])
+      cache_array[i]=Gridap.Arrays.return_cache(k,a.array[i])
     end
   end
 
   # 4. Generate ArrayBlock with the output entries
-  bi       = evaluate!(ci,k,ai)
+  bi       = Gridap.Arrays.evaluate!(ci,k,ai)
   barray   = Array{typeof(bi),length(size(a.array))}(undef,size(a.array))
   btouched = copy(a.touched)
   b=Gridap.Fields.ArrayBlock(barray,btouched)
@@ -31,14 +31,14 @@ end
 function Gridap.Arrays.evaluate!(
   cache,
   k::AddNaiveInnerMostBlockLevelMap,
-  a::Gridap.Fields.ArrayBlock{<:ArrayBlock}) where{A}
+  a::Gridap.Fields.ArrayBlock{<:ArrayBlock})
 
   b,cache_array=cache
   Gridap.Helpers.@check size(cache_array) == size(a.array)
   Gridap.Helpers.@check b.touched == a.touched
   for i in eachindex(a.array)
     if a.touched[i]
-      b.array[i]=evaluate!(cache_array[i],k,a.array[i])
+      b.array[i]=Gridap.Arrays.evaluate!(cache_array[i],k,a.array[i])
     end
   end
   b
@@ -48,17 +48,29 @@ function Gridap.Arrays.return_cache(
   k::AddNaiveInnerMostBlockLevelMap,
   a::Gridap.Fields.ArrayBlock{T}) where T
 
+  n = length(size(a.array))
+  Gridap.Helpers.@check n==1 || n==2
+
   # Generate an ArrayBlock with the same rank/shape as "a"
   # but 1x1 blocks instead of scalar entries
-  ba_array=Array{VectorBlock{T}}(undef,size(a.array))
+  if n==1
+    ba_array=Array{VectorBlock{T}}(undef,size(a.array))
+  else
+    ba_array=Array{MatrixBlock{T}}(undef,size(a.array))
+  end
   ba_touched=copy(a.touched)
   ba=Gridap.Fields.ArrayBlock(ba_array,ba_touched)
 
   # Generate 1x1 Naive ArrayBlocks
   for i in eachindex(ba_array)
     if ba_touched[i]
-      sb_array=Vector{T}(undef,1)
-      sb_touched=Vector{Bool}(undef,1)
+      if n==1
+        sb_array=Vector{T}(undef,1)
+        sb_touched=Vector{Bool}(undef,1)
+      else
+        sb_array=Matrix{T}(undef,1,1)
+        sb_touched=Matrix{Bool}(undef,1,1)
+      end
       sb_touched[1]=true
       sb=Gridap.Fields.ArrayBlock(sb_array,sb_touched)
       ba_array[i]=sb

src/BackwardStaticCondensationMap.jl

@@ -29,9 +29,9 @@ function Gridap.Arrays.return_cache(k::BackwardStaticCondensationMap{IFT,BFT},
   b::Gridap.Fields.VectorBlock{<:AbstractVector{T}},
   x::Gridap.Fields.VectorBlock{<:AbstractVector{T}}) where {IFT,BFT,T}
   m=DensifyInnerMostBlockLevelMap()
-  cm=return_cache(m,x)
+  cm=Gridap.Arrays.return_cache(m,x)
   mx=evaluate!(cm,m,x)
-  ck=return_cache(k,A,b,mx)
+  ck=Gridap.Arrays.return_cache(k,A,b,mx)
   (m,cm,ck)
 end
 

src/GridapAPIExtensions.jl

@@ -76,15 +76,25 @@ function Gridap.Arrays.lazy_map(
   lazy_map(Gridap.Fields.LinearCombinationMap(:),i_to_values,i_to_basis_x)
 end
 
+function Gridap.Arrays.lazy_map(
+  ::typeof(evaluate), a::Gridap.Arrays.LazyArray{<:Fill{typeof(Gridap.Fields.linear_combination)}},
+                      x::Fill{<:Gridap.Fields.VectorBlock})
+
+  i_to_values = a.args[1]
+  i_to_basis = a.args[2]
+  i_to_basis_x = lazy_map(evaluate,i_to_basis,x)
+  lazy_map(Gridap.Fields.LinearCombinationMap(:),i_to_values,i_to_basis_x)
+end
+
 # The following three function overloads are required in order to be able to
 # have a FE function in place of the trial operand in an integral over the
 # cell boundaries
-function Gridap.Fields.linear_combination(u::Vector,
+function Gridap.Fields.linear_combination(u::AbstractArray,
                                           f::Gridap.Fields.ArrayBlock)
-  i::Int = findfirst(f.touched)
+  i = findfirst(f.touched)
   fi = f.array[i]
   ufi = Gridap.Fields.linear_combination(u,fi)
-  g = Vector{typeof(ufi)}(undef,length(f.touched))
+  g = Array{typeof(ufi)}(undef,size(f.touched))
   for i in eachindex(f.touched)
     if f.touched[i]
       g[i] = Gridap.Fields.linear_combination(u,f.array[i])
@@ -94,25 +104,25 @@ function Gridap.Fields.linear_combination(u::Vector,
 end
 
 function Gridap.Arrays.return_cache(k::Gridap.Fields.LinearCombinationMap,
-                                    u::Vector,
+                                    u::AbstractArray,
                                     fx::Gridap.Fields.ArrayBlock)
-  i::Int = findfirst(fx.touched)
+  i = findfirst(fx.touched)
   fxi = fx.array[i]
-  li = return_cache(k,u,fxi)
+  li = Gridap.Arrays.return_cache(k,u,fxi)
   ufxi = evaluate!(li,k,u,fxi)
-  l = Vector{typeof(li)}(undef,size(fx.array))
-  g = Vector{typeof(ufxi)}(undef,size(fx.array))
+  l = Array{typeof(li)}(undef,size(fx.array))
+  g = Array{typeof(ufxi)}(undef,size(fx.array))
   for i in eachindex(fx.array)
     if fx.touched[i]
-      l[i] = return_cache(k,u,fx.array[i])
+      l[i] = Gridap.Arrays.return_cache(k,u,fx.array[i])
     end
   end
   ArrayBlock(g,fx.touched),l
 end
 
 function Gridap.Arrays.evaluate!(cache,
                                  k::Gridap.Fields.LinearCombinationMap,
-                                 u::Vector,
+                                 u::AbstractArray,
                                  fx::Gridap.Fields.ArrayBlock)
   g,l = cache
   Gridap.Helpers.@check g.touched == fx.touched
@@ -144,13 +154,13 @@ function Gridap.Arrays.return_cache(k::Gridap.Fields.LinearCombinationMap,
                                     fx::Gridap.Fields.ArrayBlock)
   i::Int = findfirst(fx.touched)
   fxi = fx.array[i]
-  li = return_cache(k,u[i],fxi)
+  li = Gridap.Arrays.return_cache(k,u[i],fxi)
   ufxi = evaluate!(li,k,u[i],fxi)
   l = Vector{typeof(li)}(undef,size(fx.array))
   g = Vector{typeof(ufxi)}(undef,size(fx.array))
   for i in eachindex(fx.array)
     if fx.touched[i]
-      l[i] = return_cache(k,u[i],fx.array[i])
+      l[i] = Gridap.Arrays.return_cache(k,u[i],fx.array[i])
     end
   end
   ArrayBlock(g,fx.touched),l
@@ -252,7 +262,7 @@ end
 function Gridap.Arrays.return_cache(
   k::Broadcasting{typeof(∘)},
   f::Gridap.Fields.VectorBlock,
-  h::Gridap.Fields.VectorBlock) where {A,N}
+  h::Gridap.Fields.VectorBlock)
   Gridap.Helpers.@check length(f)==length(h)
   Gridap.Helpers.@check f.touched==h.touched
   fi = Gridap.Arrays.testitem(f)
@@ -510,7 +520,7 @@ function Gridap.Arrays.return_cache(
   k::typeof(setup_bulk_to_skeleton_l2_projected_fields),
   A::Array,
   B::AbstractArray{<:Gridap.Fields.Field})
-  return_cache(Gridap.Fields.linear_combination,A,B)
+  Gridap.Arrays.return_cache(Gridap.Fields.linear_combination,A,B)
 end
 
 function Gridap.Arrays.evaluate!(
@@ -527,7 +537,7 @@ function Gridap.Arrays.return_cache(
   B::Transpose{T,<:AbstractArray{<:Gridap.Fields.Field}}) where T
   c=return_cache(k,A,B.parent)
   v=evaluate!(c,k,A,B.parent)
-  return_cache(transpose,v), c
+  Gridap.Arrays.return_cache(transpose,v), c
 end
 
 function Gridap.Arrays.evaluate!(
@@ -616,7 +626,7 @@ function Gridap.Arrays.return_cache(
   A::Gridap.Fields.MatrixBlock,
   B::Gridap.Fields.MatrixBlock)
 
-  Gridap.Helpers.@check size(A.array) == size(B.array)
+  Gridap.Helpers.@check size(A.array)==size(B.array)
 
   nA=findall(A.touched)
   nB=findall(B.touched)
@@ -698,7 +708,7 @@ function Gridap.Arrays.return_cache(
   B::Gridap.Fields.VectorBlock{<:Matrix})
 
   Gridap.Helpers.@check size(A.array)[1] == size(A.array)[2]
-  Gridap.Helpers.@check size(A.array)[1] == size(A.array)[2]
+  Gridap.Helpers.@check size(A.array)[1] == size(B.array)[1]
 
   nA=findall(A.touched)
   nB=findall(B.touched)
@@ -725,7 +735,53 @@ function Gridap.Arrays.evaluate!(
   B::Gridap.Fields.VectorBlock{<:Matrix})
   r,c=cache
   Gridap.Helpers.@check size(A.array)[1] == size(A.array)[2]
+  Gridap.Helpers.@check size(A.array)[1] == size(B.array)[1]
+  Gridap.Helpers.@check length(findall(A.touched)) == 1
+  Gridap.Helpers.@check length(findall(B.touched)) == 1
+  Gridap.Helpers.@check r.touched[1]
+  ai=testitem(A)
+  bi=testitem(B)
+  r.array[1]=evaluate!(c,k,ai,bi)
+  r
+end
+
+# A [f,f]      , e.g., [2,2] nonzero block (matrix to be inverted)
+# B [f,1]      , e.g., [2,1] nonzero block (several RHS)
+# Return [1]
+function Gridap.Arrays.return_cache(
+  k::typeof(compute_bulk_to_skeleton_l2_projection_dofs),
+  A::Gridap.Fields.MatrixBlock{<:Matrix},
+  B::Gridap.Fields.MatrixBlock{<:Matrix})
+
+  Gridap.Helpers.@check size(A.array)[1] == size(A.array)[2]
+  Gridap.Helpers.@check size(A.array)[1] == size(B.array)[1]
+
+  nA=findall(A.touched)
+  nB=findall(B.touched)
+
+  Gridap.Helpers.@check length(nA)==length(nB)
+  Gridap.Helpers.@check length(nA)==1
+  Gridap.Helpers.@check nA[1][1]==nB[1][1]
+  Gridap.Helpers.@check nA[1][2]==nB[1][1]
+
+  ai=testitem(A)
+  bi=testitem(B)
+  T=Gridap.Arrays.return_type(k,ai,bi)
+  touched  = Vector{Bool}(undef,1)
+  touched .= true
+  r = Vector{T}(undef,1)
+  c = Gridap.Arrays.return_cache(k,ai,bi)
+  Gridap.Fields.ArrayBlock(r,touched), c
+end
+
+function Gridap.Arrays.evaluate!(
+  cache,
+  k::typeof(compute_bulk_to_skeleton_l2_projection_dofs),
+  A::Gridap.Fields.MatrixBlock{<:Matrix},
+  B::Gridap.Fields.MatrixBlock{<:Matrix})
+  r,c=cache
   Gridap.Helpers.@check size(A.array)[1] == size(A.array)[2]
+  Gridap.Helpers.@check size(A.array)[1] == size(B.array)[1]
   Gridap.Helpers.@check length(findall(A.touched)) == 1
   Gridap.Helpers.@check length(findall(B.touched)) == 1
   Gridap.Helpers.@check r.touched[1]
@@ -826,3 +882,53 @@ end
 function Gridap.FESpaces._compute_cell_ids(uh,ttrian::SkeletonTriangulation)
   collect(1:num_cells(ttrian))
 end
+
+# Required by LocalFEProjector
+function Gridap.Arrays.return_cache(::typeof(\),a::AbstractArray{<:Number},b::AbstractArray{<:Number})
+  c = a\b
+  Gridap.Arrays.CachedArray(c)
+end
+
+# Required by LocalFEProjector
+function Gridap.Arrays.evaluate!(cache,::typeof(\),a::AbstractMatrix{<:Number},b::AbstractVector{<:Number})
+  m = size(a,1)
+  Gridap.Arrays.setsize!(cache,(m,))
+  c = cache.array
+  c .= a\b
+  c
+end
+
+# Required by LocalFEProjector
+function Gridap.Arrays.evaluate!(cache,::typeof(\),a::AbstractMatrix{<:Number},b::AbstractMatrix{<:Number})
+  m = size(a,1)
+  n = size(b,2)
+  Gridap.Arrays.setsize!(cache,(m,n))
+  c = cache.array
+  c .= a\b
+  c
+end
+
+function Gridap.Arrays.evaluate!(cache,
+                   k::Gridap.Fields.DensifyInnerMostBlockLevelMap,
+                   a::Array{T}) where {T<:Number}
+   a
+end
+
+function Gridap.Fields.linear_combination(a::AbstractArray{<:Number},b::Transpose)
+  c=Gridap.Fields.linear_combination(a,b.parent)
+  Transpose(c)
+end
+
+function Gridap.Arrays.return_cache(k::Gridap.Fields.LinearCombinationMap,
+                                    v::AbstractMatrix,
+                                    fx::Gridap.Fields.TransposeFieldIndices)
+  return_cache(k,v,fx.matrix)
+end
+
+function Gridap.Arrays.evaluate!(cache,
+  k::Gridap.Fields.LinearCombinationMap,
+  v::AbstractMatrix,
+  fx::Gridap.Fields.TransposeFieldIndices)
+  gx=evaluate!(cache,k,v,fx.matrix)
+  Gridap.Fields.TransposeFieldIndices(gx)
+end

src/GridapHybrid.jl

@@ -33,8 +33,26 @@ export HybridFEOperator
 
 include("HybridLinearSolvers.jl")
 
-
 include("GridapAPIExtensions.jl")
 include("GridapTmpModifications.jl")
 
+# HHO-methods specific components
+include("OrthogonalBasisRefFEs.jl")
+include("OrthogonalBasisFESpaces.jl")
+export OrthogonalBasisRefFE
+export orthogonal_basis
+
+include("MonomialBasisRefFEs.jl")
+export MonomialBasisRefFE
+export monomial_basis
+
+include("LocalFEOperatorsTools.jl")
+include("ReconstructionFEOperators.jl")
+include("ProjectionFEOperators.jl")
+export LocalFEOperator
+export ReconstructionFEOperator
+export ProjectionFEOperator
+export SingleValued
+export MultiValued
+
 end # module