Moment based reffes

src/Adaptivity/MacroFEs.jl

@@ -229,26 +229,32 @@ function Arrays.return_cache(a::MacroFEBasis,xc::AbstractArray{<:Point})
   k = CoarseToFinePointMap()
   geo_cache = return_cache(k,rr,xc)
   xf, ids = evaluate!(geo_cache,k,rr,xc)
+  xf_cache = array_cache(xf)
 
-  eval_caches = map(return_cache,a.fine_data,xf)
+  # NOTE: xf may be empty for some subcells, so it's safer to use testvalue
+  xt = testvalue(xc)
+  eval_caches = map(ffields -> return_cache(ffields,xt),a.fine_data)
 
-  T = eltype(evaluate!(first(eval_caches),first(a.fine_data),first(xf)))
+  T = eltype(evaluate!(first(eval_caches),first(a.fine_data),xt))
   res_cache = CachedArray(zeros(T,length(xc),length(a)))
-  return res_cache, k, geo_cache, eval_caches
+  return res_cache, k, geo_cache, eval_caches, xf_cache
 end
 
 function Arrays.evaluate!(caches, a::MacroFEBasis,xc::AbstractArray{<:Point})
-  res_cache, k, geo_cache, eval_caches = caches
+  res_cache, k, geo_cache, eval_caches, xf_cache = caches
   setsize!(res_cache,(length(xc),length(a)))
   res = res_cache.array
 
   fill!(res,zero(eltype(res)))
   xf, ids = evaluate!(geo_cache,k,a.rrule,xc)
 
-  for fcell in 1:num_subcells(a.rrule)
-    r = xf.ptrs[fcell]:xf.ptrs[fcell+1]-1
-    vals = evaluate!(eval_caches[fcell],a.fine_data[fcell],view(xf.data,r))
-    I  = view(ids.data,r)
+  for fcell in 1:num_subcells(a.rrule) 
+    xf_k = getindex!(xf_cache,xf,fcell)
+    if isempty(xf_k)
+      continue
+    end
+    vals = evaluate!(eval_caches[fcell],a.fine_data[fcell],xf_k)
+    I = view(ids,fcell)
     J = a.ids.fcell_to_cids[fcell]
     res[I,J] .= vals
   end

src/FESpaces/CurlConformingFESpaces.jl

@@ -2,8 +2,8 @@
 function get_cell_dof_basis(
   model::DiscreteModel,
   cell_reffe::AbstractArray{<:GenericRefFE{Nedelec}},
-  ::CurlConformity)
-
+  ::CurlConformity
+)
   cell_dofs = lazy_map(get_dof_basis,cell_reffe)
   cell_ownids = lazy_map(get_face_own_dofs,cell_reffe)
   cell_map = get_cell_map(Triangulation(model))

src/FESpaces/DivConformingFESpaces.jl

@@ -24,132 +24,129 @@
 #   shape function corresponding to the global DoF.
 # * We do NOT have to use the signed determinant, but its absolute value, in the Piola Map.
 
-struct TransformRTDofBasis{Dc,Dp} <: Map end ;
-
-function get_cell_dof_basis(model::DiscreteModel,
-                            cell_reffe::AbstractArray{<:GenericRefFE{<:DivConforming}},
-                            ::DivConformity,
-                            sign_flip=get_sign_flip(model, cell_reffe))
-    cell_map  = get_cell_map(Triangulation(model))
-    phi       = cell_map[1]
-    Jt        = lazy_map(Broadcasting(∇),cell_map)
-    x         = lazy_map(get_nodes,lazy_map(get_dof_basis,cell_reffe))
-    Jtx       = lazy_map(evaluate,Jt,x)
-    reffe     = cell_reffe[1]
-    Dc        = num_dims(reffe)
-    # @santiagobadia: A hack here, for RT returns Float64 and for BDM VectorValue{Float64}
-    et        = eltype(return_type(get_prebasis(reffe)))
-    pt        = Point{Dc,et}
-    Dp        = first(size(return_type(phi,zero(pt))))
-    k         = TransformRTDofBasis{Dc,Dp}()
-    lazy_map(k,cell_reffe,Jtx,sign_flip)
+struct TransformRTDofBasis{Dc,Dp} <: Map end
+
+function get_cell_dof_basis(
+  model::DiscreteModel{Dc,Dp},
+  cell_reffe::AbstractArray{<:GenericRefFE{<:DivConforming}},
+  ::DivConformity,
+  sign_flip = get_sign_flip(model, cell_reffe)
+) where {Dc,Dp}
+  cell_map  = get_cell_map(get_grid(model))
+  Jt  = lazy_map(Broadcasting(∇),cell_map)
+  x   = lazy_map(get_nodes,lazy_map(get_dof_basis,cell_reffe))
+  Jtx = lazy_map(evaluate,Jt,x)
+  k = TransformRTDofBasis{Dc,Dp}()
+  lazy_map(k,cell_reffe,Jtx,sign_flip)
 end
 
-function get_cell_shapefuns(model::DiscreteModel,
-                            cell_reffe::AbstractArray{<:GenericRefFE{<:DivConforming}},
-                            ::DivConformity,
-                            sign_flip=get_sign_flip(model, cell_reffe))
-    cell_reffe_shapefuns=lazy_map(get_shapefuns,cell_reffe)
-    k=ContraVariantPiolaMap()
-    lazy_map(k,
-             cell_reffe_shapefuns,
-             get_cell_map(Triangulation(model)),
-             lazy_map(Broadcasting(constant_field), sign_flip))
+function get_cell_shapefuns(
+  model::DiscreteModel,
+  cell_reffe::AbstractArray{<:GenericRefFE{<:DivConforming}},
+  ::DivConformity,
+  sign_flip = get_sign_flip(model, cell_reffe)
+)
+  cell_map = get_cell_map(get_grid(model))
+  cell_shapefuns = lazy_map(get_shapefuns,cell_reffe)
+  k = ContraVariantPiolaMap()
+  lazy_map(k,cell_shapefuns,cell_map,lazy_map(Broadcasting(constant_field),sign_flip))
 end
 
 struct SignFlipMap{T} <: Map
   model::T
+  facet_owners::Vector{Int32}
 end
 
-function return_cache(k::SignFlipMap,reffe,cell_id)
-  model = k.model
-  D = num_cell_dims(model)
-  gtopo = get_grid_topology(model)
-
-  # Extract composition among cells and facets
-  cell_wise_facets_ids = get_faces(gtopo, D, D - 1)
-  cache_cell_wise_facets_ids = array_cache(cell_wise_facets_ids)
+function SignFlipMap(model)
+  facet_owners = compute_facet_owners(model)
+  SignFlipMap(model,facet_owners)
+end
 
-  # Extract cells around facets
-  cells_around_facets = get_faces(gtopo, D - 1, D)
-  cache_cells_around_facets = array_cache(cells_around_facets)
+function return_cache(k::SignFlipMap,reffe,facet_own_dofs,cell)
+  model = k.model
+  Dc = num_cell_dims(model)
+  topo = get_grid_topology(model)
 
-  (cell_wise_facets_ids,
-   cache_cell_wise_facets_ids,
-   cells_around_facets,
-   cache_cells_around_facets,
-   CachedVector(Bool))
+  cell_facets = get_faces(topo, Dc, Dc-1)
+  cell_facets_cache = array_cache(cell_facets)
 
+  return cell_facets,cell_facets_cache,CachedVector(Bool)
 end
 
-function evaluate!(cache,k::SignFlipMap,reffe,cell_id)
-  model = k.model
-
-  cell_wise_facets_ids,
-  cache_cell_wise_facets_ids,
-  cells_around_facets,
-  cache_cells_around_facets,
-  sign_flip_cached = cache
+function evaluate!(cache,k::SignFlipMap,reffe,facet_own_dofs,cell)
+  cell_facets,cell_facets_cache,sign_flip_cache = cache
+  facet_owners = k.facet_owners
 
-  setsize!(sign_flip_cached, (num_dofs(reffe),))
-  sign_flip = sign_flip_cached.array
+  setsize!(sign_flip_cache, (num_dofs(reffe),))
+  sign_flip = sign_flip_cache.array
   sign_flip .= false
 
-  D = num_dims(reffe)
-  face_own_dofs = get_face_own_dofs(reffe)
-  facet_lid = get_offsets(get_polytope(reffe))[D] + 1
-  cell_facets_ids = getindex!(cache_cell_wise_facets_ids,
-                              cell_wise_facets_ids,
-                              cell_id)
-  for facet_gid in cell_facets_ids
-      facet_cells_around = getindex!(cache_cells_around_facets,
-           cells_around_facets,
-           facet_gid)
-      is_slave = (findfirst((x) -> (x == cell_id), facet_cells_around) == 2)
-      if is_slave
-          for dof in face_own_dofs[facet_lid]
-              sign_flip[dof] = true
-          end
+  facets = getindex!(cell_facets_cache,cell_facets,cell)
+  for (lfacet,facet) in enumerate(facets)
+    owner = facet_owners[facet]
+    if owner != cell
+      for dof in facet_own_dofs[lfacet]
+        sign_flip[dof] = true
       end
-      facet_lid = facet_lid + 1
+    end
   end
-  sign_flip
+
+  return sign_flip
 end
 
-function get_sign_flip(model::DiscreteModel,
-                       cell_reffe::AbstractArray{<:GenericRefFE{<:DivConforming}})
-    lazy_map(SignFlipMap(model),
-            cell_reffe,
-            IdentityVector(Int32(num_cells(model))))
+function get_sign_flip(
+  model::DiscreteModel{Dc},
+  cell_reffe::AbstractArray{<:GenericRefFE{<:DivConforming}}
+) where Dc
+  # Comment: lazy_maps on cell_reffes are very optimised, since they are CompressedArray/FillArray
+  get_facet_own_dofs(reffe) = view(get_face_own_dofs(reffe),get_dimrange(get_polytope(reffe),Dc-1))
+  cell_facet_own_dofs = lazy_map(get_facet_own_dofs,cell_reffe)
+  cell_ids = IdentityVector(Int32(num_cells(model)))
+  lazy_map(SignFlipMap(model),cell_reffe,cell_facet_own_dofs,cell_ids)
 end
 
-function return_cache(::TransformRTDofBasis{Dc,Dp},
-                      reffe::GenericRefFE{<:DivConforming},
-                      Jtx,
-                      ::AbstractVector{Bool}) where {Dc,Dp}
-  p = get_polytope(reffe)
-  prebasis = get_prebasis(reffe)
-  order = get_order(prebasis)
+function compute_facet_owners(model::DiscreteModel{Dc,Dp}) where {Dc,Dp}
+  topo = get_grid_topology(model)
+  facet_to_cell = get_faces(topo, Dc-1, Dc)
+
+  nfacets = num_faces(topo, Dc-1)
+  owners = Vector{Int32}(undef, nfacets)
+  for facet in 1:nfacets
+    facet_cells = view(facet_to_cell, facet)
+    owners[facet] = first(facet_cells)
+  end
+
+  return owners
+end
+
+function return_cache(
+  ::TransformRTDofBasis{Dc,Dp},
+  reffe::GenericRefFE{<:DivConforming},
+  Jtx,
+  ::AbstractVector{Bool}
+) where {Dc,Dp}
   # @santiagobadia: Hack as above
-  et = eltype(return_type(prebasis))
+  et = eltype(return_type(get_prebasis(reffe)))
   dofs = get_dof_basis(reffe)
-  nodes, nf_nodes, nf_moments =  get_nodes(dofs),
-                                 get_face_nodes_dofs(dofs),
-                                 get_face_moments(dofs)
+
+  nodes = get_nodes(dofs)
+  nf_nodes = get_face_nodes_dofs(dofs)
+  nf_moments = get_face_moments(dofs)
   db = MomentBasedDofBasis(nodes,nf_moments,nf_nodes)
   face_moments = [ similar(i,VectorValue{Dp,et})  for i in nf_moments ]
 
-  cache = (db.nodes, db.face_nodes, nf_moments, face_moments)
-  cache
+  return db.nodes, db.face_nodes, nf_moments, face_moments
 end
 
-function evaluate!(cache,
-                   ::TransformRTDofBasis,
-                   reffe::GenericRefFE{<:DivConforming},
-                   Jt_q,
-                   sign_flip::AbstractVector{Bool})
+function evaluate!(
+  cache,
+  ::TransformRTDofBasis,
+  reffe::GenericRefFE{<:DivConforming},
+  Jt_q,
+  sign_flip::AbstractVector{Bool}
+)
   nodes, nf_nodes, nf_moments, face_moments = cache
-  face_own_dofs=get_face_own_dofs(reffe)
+  face_own_dofs = get_face_own_dofs(reffe)
   for face in 1:length(face_moments)
     nf_moments_face   = nf_moments[face]
     face_moments_face = face_moments[face]
@@ -170,20 +167,23 @@ end
 
 
 # Support for DIV operator
-function DIV(f::LazyArray{<:Fill{T}}) where T
-  df=DIV(f.args[1])
-  k=f.maps.value
+
+function DIV(f::LazyArray{<:Fill})
+  df = DIV(f.args[1])
+  k  = f.maps.value
   lazy_map(k,df)
 end
+
 function DIV(f::LazyArray{<:Fill{Broadcasting{Operation{ContraVariantPiolaMap}}}})
   ϕrgₖ       = f.args[1]
-  fsign_flip = f.args[4]
-  div_ϕrgₖ = lazy_map(Broadcasting(divergence),ϕrgₖ)
-  fsign_flip=lazy_map(Broadcasting(Operation(x->(-1)^x)), fsign_flip)
+  fsign_flip = f.args[3]
+  div_ϕrgₖ   = lazy_map(Broadcasting(divergence),ϕrgₖ)
+  fsign_flip = lazy_map(Broadcasting(Operation(x->(-1)^x)), fsign_flip)
   lazy_map(Broadcasting(Operation(*)),fsign_flip,div_ϕrgₖ)
 end
+
 function DIV(a::LazyArray{<:Fill{typeof(linear_combination)}})
-  i_to_basis = DIV(a.args[2])
+  i_to_basis  = DIV(a.args[2])
   i_to_values = a.args[1]
   lazy_map(linear_combination,i_to_values,i_to_basis)
 end

src/Fields/FieldArrays.jl

@@ -509,6 +509,22 @@ for T in (:(Point),:(AbstractArray{<:Point}))
       evaluate!(r,bm,rs...)
     end
 
+    function return_cache(k::BroadcastOpFieldArray{typeof(∘)},x::$T)
+      f, g = k.args
+      cg = return_cache(g,x)
+      gx = evaluate!(cg,g,x)
+      cf = return_cache(f,gx)
+      return cg, cf
+    end
+
+    function evaluate!(cache, k::BroadcastOpFieldArray{typeof(∘)},x::$T)
+      cg, cf = cache
+      f, g = k.args
+      gx = evaluate!(cg,g,x)
+      fgx = evaluate!(cf,f,gx)
+      return fgx
+    end
+
   end
 end
 
@@ -694,3 +710,43 @@ for op in (:*,:⋅,:⊙,:⊗)
     end
   end
 end
+
+# Optimisations to 
+# lazy_map(Broadcasting(constant_field),a::AbstractArray{<:AbstractArray{<:Number}})
+
+struct ConstantFieldArray{T,N,A} <: AbstractArray{ConstantField{T},N}
+  values::A
+  function ConstantFieldArray(values::AbstractArray{T,N}) where {T,N}
+    A = typeof(values)
+    new{T,N,A}(values)
+  end
+end
+
+Base.size(a::ConstantFieldArray) = size(a.values)
+Base.axes(a::ConstantFieldArray) = axes(a.values)
+Base.getindex(a::ConstantFieldArray,i::Integer) = ConstantField(a.values[i])
+
+function return_value(::Broadcasting{typeof(constant_field)},values::AbstractArray{<:Number})
+  ConstantFieldArray(values)
+end
+
+function evaluate!(cache,::Broadcasting{typeof(constant_field)},values::AbstractArray{<:Number})
+  ConstantFieldArray(values)
+end
+
+function evaluate!(c,f::ConstantFieldArray,x::Point)
+  return f.values
+end
+
+function return_cache(f::ConstantFieldArray{T},x::AbstractArray{<:Point}) where T
+  return CachedArray(zeros(T,(size(x)...,size(f)...)))
+end
+
+function evaluate!(c,f::ConstantFieldArray{T},x::AbstractArray{<:Point}) where T
+  setsize!(c,(size(x)...,size(f)...))
+  r = c.array
+  for i in eachindex(x)
+    r[i,:] .= f.values
+  end
+  return r
+end