Added first HDG test

Author: JordiManyer

src/FESpaces/PatchAssemblers.jl

@@ -109,10 +109,11 @@ function get_patch_dofs(space::FESpace,ptopo::PatchTopology)
   patch_to_faces = get_patch_faces(ptopo,Df)
 
   patch_to_dofs = map(patch_to_faces) do pfaces
-    tfaces = filter(!iszero,face_to_tface[pfaces])
+    tfaces = filter(x->x>0,face_to_tface[pfaces])
     dofs = SortedSet{Int}()
     for tface in tfaces
-      push!(dofs,view(tface_to_dofs,tface)...)
+      tface_dofs = filter(x->x>0,view(tface_to_dofs,tface))
+      !isempty(tface_dofs) && push!(dofs,tface_dofs...)
     end
     collect(dofs)
   end |> Table
@@ -235,17 +236,52 @@ struct PatchAssemblyMap{T} <: Map
   cell_data
 end
 
+function _alloc_cache(::PatchAssemblyStrategy,s::NTuple{N,Int}) where N
+  zeros(s...)
+end
+
+function _alloc_cache(::PatchAssemblyStrategy,s::NTuple{N,<:BlockedOneTo}) where N
+  bs = map(blocklength,s)
+  ss = map(blocklengths,s)
+  array = [ zeros(ntuple(i -> ss[i][I[i]], Val(N))) for I in CartesianIndices(bs) ]
+  ArrayBlock(array,fill(true,bs))
+end
+
+function _alloc_cache(b::ArrayBlockView,s::NTuple{N,<:BlockedOneTo}) where N
+  bs = map(blocklength,s)
+  ss = map(blocklengths,s)
+  array = [ zeros(ntuple(i -> ss[i][I[i]], Val(N))) for I in CartesianIndices(bs) ]
+  bmap = ifelse(N == 2, b.block_map, map(idx -> CartesianIndex(idx[1]), diag(b.block_map)))
+  ArrayBlockView(ArrayBlock(array,fill(true,bs)),bmap)
+end
+
+function _resize_cache!(a,::PatchAssemblyStrategy,s::NTuple{N,Int}) where N
+  setsize!(a,s)
+end
+
+function _resize_cache!(a,::PatchAssemblyStrategy,s::NTuple{N,<:BlockedOneTo}) where N
+  bs = map(blocklength,s)
+  ss = map(blocklengths,s)
+  for I in CartesianIndices(bs)
+    setsize!(a.array[I],ntuple(i -> ss[i][I[i]], Val(N)))
+  end
+end
+
+function _resize_cache!(a,b::ArrayBlockView,s::NTuple{N,<:BlockedOneTo}) where N
+  bs = map(blocklength,s)
+  ss = map(blocklengths,s)
+  for I in CartesianIndices(bs)
+    setsize!(a.array.array[I],ntuple(i -> ss[i][I[i]], Val(N)))
+  end
+end
+
+_unview(a) = a
+_unview(a::ArrayBlockView) = a.array
+
 # Mat & Vec assembly
 
 function Arrays.return_cache(k::PatchAssemblyMap,patch)
-  _zeros(s::NTuple{N,Int}) where N = zeros(s...)
-  function _zeros(s::NTuple{N,BlockedOneTo}) where N
-    bs = map(blocklength,s)
-    ss = map(blocklengths,s)
-    array = [ zeros(ntuple(i -> ss[i][I[i]], Val(N))) for I in CartesianIndices(bs) ]
-    ArrayBlock(array,fill(true,bs))
-  end
-  res = _zeros(k.patch_sizes[patch])
+  res = _alloc_cache(k.assem.strategy,k.patch_sizes[patch])
   caches = patch_assembly_cache(res,k.cell_data)
 
   c_res = CachedArray(res)
@@ -255,33 +291,19 @@ function Arrays.return_cache(k::PatchAssemblyMap,patch)
 end
 
 function Arrays.evaluate!(cache,k::PatchAssemblyMap,patch)
-  _setsize!(a,s::NTuple{N,Int}) where N = setsize!(a,s)
-  function _setsize!(a,s::NTuple{N,BlockedOneTo}) where N
-    bs = map(blocklength,s)
-    ss = map(blocklengths,s)
-    for I in CartesianIndices(bs)
-      setsize!(a.array[I],ntuple(i -> ss[i][I[i]], Val(N)))
-    end
-  end
   c_res, uwr_cache, caches = cache
-  _setsize!(c_res,k.patch_sizes[patch])
+  _resize_cache!(c_res,k.assem.strategy,k.patch_sizes[patch])
   res = evaluate!(uwr_cache,Fields.unwrap_cached_array,c_res)
   Fields._zero_entries!(res)
   patch_assembly!(caches,res,k.cell_data,patch)
-  return res
+  return _unview(res)
 end
 
 # Mat-Vec assembly
 
 function Arrays.return_cache(k::PatchAssemblyMap{<:Tuple{<:Tuple,<:Tuple}},patch)
-  _zeros(s::NTuple{N,Int}) where N = zeros(s...)
-  function _zeros(s::NTuple{N,BlockedOneTo}) where N
-    bs = map(blocklength,s)
-    ss = map(blocklengths,s)
-    array = [ zeros(ntuple(i -> ss[i][I[i]], Val(N))) for I in CartesianIndices(bs) ]
-    ArrayBlock(array,fill(true,bs))
-  end
-  mat, vec = _zeros(k.patch_sizes[patch][1]), _zeros(k.patch_sizes[patch][2])
+  mat = _alloc_cache(k.assem.strategy,k.patch_sizes[patch][1])
+  vec = _alloc_cache(k.assem.strategy,k.patch_sizes[patch][2])
 
   matvecdata, matdata, vecdata = k.cell_data
   mat_caches = patch_assembly_cache(mat,matdata)
@@ -296,19 +318,11 @@ function Arrays.return_cache(k::PatchAssemblyMap{<:Tuple{<:Tuple,<:Tuple}},patch
 end
 
 function Arrays.evaluate!(cache,k::PatchAssemblyMap{<:Tuple{<:Tuple,<:Tuple}},patch)
-  _setsize!(a,s::NTuple{N,Int}) where N = setsize!(a,s)
-  function _setsize!(a,s::NTuple{N,BlockedOneTo}) where N
-    bs = map(blocklength,s)
-    ss = map(blocklengths,s)
-    for I in CartesianIndices(bs)
-      setsize!(a.array[I],ntuple(i -> ss[i][I[i]], Val(N)))
-    end
-  end
   c_mat, c_vec, uwm_cache, uwv_cache, matvec_caches, mat_caches, vec_caches = cache
   matvecdata, matdata, vecdata = k.cell_data
 
-  _setsize!(c_mat,k.patch_sizes[patch][1])
-  _setsize!(c_vec,k.patch_sizes[patch][2])
+  _resize_cache!(c_mat,k.assem.strategy,k.patch_sizes[patch][1])
+  _resize_cache!(c_vec,k.assem.strategy,k.patch_sizes[patch][2])
   mat = evaluate!(uwm_cache,Fields.unwrap_cached_array,c_mat)
   vec = evaluate!(uwv_cache,Fields.unwrap_cached_array,c_vec)
 
@@ -319,11 +333,11 @@ function Arrays.evaluate!(cache,k::PatchAssemblyMap{<:Tuple{<:Tuple,<:Tuple}},pa
   patch_assembly!(mat_caches,mat,matdata,patch)
   patch_assembly!(vec_caches,vec,vecdata,patch)
 
-  return mat, vec
+  return _unview(mat), _unview(vec)
 end
 
-const MatOrMatBlock = Union{AbstractMatrix,MatrixBlock}
-const VecOrVecBlock = Union{AbstractVector,VectorBlock}
+const MatOrMatBlock = Union{AbstractMatrix,MatrixBlock,MatrixBlockView}
+const VecOrVecBlock = Union{AbstractVector,VectorBlock,VectorBlockView}
 
 function patch_assembly_cache(mat::MatOrMatBlock,cell_matdata)
   caches = ()
@@ -451,7 +465,7 @@ function Arrays.evaluate!(cache,k::StaticCondensationMap,matvec)
   @check size(mat.array) == (2,2)
   @check size(vec.array) == (2,)
 
-  Kii, Kib, Kbi, Kbb = mat.array
+  Kii, Kbi, Kib, Kbb = mat.array
   bi, bb = vec.array
 
   f = lu!(Kii,k.pivot;check=false)

src/Fields/ArrayBlocks.jl

@@ -1519,3 +1519,26 @@ end
 
 LinearAlgebra.diag(a::MatrixBlockView) = view(a.array.array, diag(a.block_map))
 LinearAlgebra.diag(a::MatrixBlock) = view(a.array,diag(CartesianIndices(a.array)))
+
+_zero_entries!(a::ArrayBlockView) = _zero_entries!(a.array)
+
+function Arrays.CachedArray(a::ArrayBlockView)
+  ArrayBlockView(CachedArray(a.array),a.block_map)
+end
+
+function unwrap_cached_array(a::ArrayBlockView)
+  cache = return_cache(unwrap_cached_array,a)
+  evaluate!(cache,unwrap_cached_array,a)
+end
+
+function return_cache(::typeof(unwrap_cached_array),a::ArrayBlockView)
+  cache = return_cache(unwrap_cached_array,a.array)
+  array = evaluate!(cache,unwrap_cached_array,a.array)
+  return ArrayBlockView(array,a.block_map), cache
+end
+
+function evaluate!(cache,::typeof(unwrap_cached_array),a::ArrayBlockView)
+  r, c = cache
+  evaluate!(c,unwrap_cached_array,a.array)
+  return r
+end

test/GeometryTests/PatchTriangulationTests.jl

@@ -54,7 +54,8 @@ x1 = lazy_map(FESpaces.LocalSolveMap(),A,b)
 x2 = lazy_map(FESpaces.LocalSolveMap(),Ab)
 
 
-# MultiField
+# MultiField + Static condensation
+
 using Gridap.MultiField
 using Gridap.CellData
 

test/GridapTests/HDGTests.jl

@@ -0,0 +1,93 @@
+
+using Gridap
+using Gridap.Geometry, Gridap.FESpaces, Gridap.MultiField
+using Gridap.CellData, Gridap.Fields
+
+function get_abs_normal_vector(trian)
+  function normal(c)
+    t = c[2] - c[1]
+    n = VectorValue(-t[2],t[1])
+    n = n/norm(n)
+    return n
+  end
+  face_coords = get_cell_coordinates(trian)
+  face_normals = lazy_map(constant_field,lazy_map(normal,face_coords))
+  return CellData.GenericCellField(face_normals,trian,ReferenceDomain())
+end
+
+function statically_condensed_assembly(ptopo,X,Y,M,M_test,a,l)
+  # Lazily assemble the global patch-systems and 
+  # perform the static condensation
+  assem = FESpaces.PatchAssembler(ptopo,X,Y)
+  full_matvecs = assemble_matrix_and_vector(a,l,assem,X,Y)
+  sc_matvecs = lazy_map(FESpaces.StaticCondensationMap(),full_matvecs)
+
+  # Regular assembly of the statically-assembled systems
+  patch_rows = assem.strategy.array.array[2,2].patch_rows
+  patch_cols = assem.strategy.array.array[2,2].patch_cols
+  matvecdata = ([sc_matvecs,],[patch_rows,],[patch_cols,])
+  matdata = ([],[],[]) # dummy matdata
+  vecdata = ([],[],[]) # dummy vecdata
+  data = (matvecdata, matdata, vecdata)
+  A, b = assemble_matrix_and_vector(SparseMatrixAssembler(M,M_test),data)
+  return A, b
+end
+
+u(x) = sin(2*π*x[1])*sin(2*π*x[2])*(1-x[1])*x[2]*(1-x[2])
+q(x) = -∇(u)(x)
+f(x) = (∇ ⋅ q)(x)
+
+nc = (4,4)
+model = UnstructuredDiscreteModel(CartesianDiscreteModel((0,1,0,1),nc))
+D = num_cell_dims(model)
+Ω = Triangulation(ReferenceFE{D}, model)
+Γ = Triangulation(ReferenceFE{D-1}, model)
+
+ptopo = Geometry.PatchTopology(model)
+Ωp = Geometry.PatchTriangulation(model,ptopo)
+Γp = Geometry.PatchBoundaryTriangulation(model,ptopo)
+
+# Reference FEs
+order = 1
+reffeV = ReferenceFE(lagrangian, VectorValue{D, Float64}, order; space=:P)
+reffeQ = ReferenceFE(lagrangian, Float64, order; space=:P)
+reffeM = ReferenceFE(lagrangian, Float64, order; space=:P)
+
+# HDG test FE Spaces
+V_test = TestFESpace(Ω, reffeV; conformity=:L2)
+Q_test = TestFESpace(Ω, reffeQ; conformity=:L2)
+M_test = TestFESpace(Γ, reffeM; conformity=:L2, dirichlet_tags="boundary")
+
+# HDG trial FE Spaces
+V = TrialFESpace(V_test)
+Q = TrialFESpace(Q_test)
+M = TrialFESpace(M_test, u)
+
+mfs = MultiField.BlockMultiFieldStyle(2,(2,1))
+Y = MultiFieldFESpace([V_test, Q_test, M_test];style=mfs)
+X = MultiFieldFESpace([V, Q, M];style=mfs)
+
+τ = 1.0 # HDG stab parameter
+
+degree = 2*(order+1)
+dΩp = Measure(Ωp,degree)
+dΓp = Measure(Γp,degree)
+nrel = get_normal_vector(Γp)
+nabs = get_abs_normal_vector(Γp)
+n = (nrel⋅nabs)⋅nabs
+
+Πn(u) = u⋅n
+Π(u) = change_domain(u,Γp,DomainStyle(u))
+a((qh,uh,sh),(vh,wh,lh)) = ∫( qh⋅vh - uh*(∇⋅vh) - qh⋅∇(wh) )dΩp + ∫(sh*Πn(vh))dΓp +
+                           ∫((Πn(qh) + τ*(Π(uh) - sh))*(Π(wh) + lh))dΓp
+l((vh,wh,hatmh)) = ∫( f*wh )*dΩp
+
+A, b = statically_condensed_assembly(ptopo,X,Y,M,M_test,a,l)
+
+solver = LUSolver()
+ns = numerical_setup(symbolic_setup(solver,A),A)
+x = zeros(size(b))
+solve!(x,ns,b)
+sh = FEFunction(M_test,x)
+
+