Fix IntegrationMap on empty subds with physical quads

src/FESpaces/Assemblers.jl

@@ -434,6 +434,7 @@ function collect_cell_vector(test::FESpace,a::DomainContribution)
   w = []
   r = []
   for strian in get_domains(a)
+    num_cells(strian) == 0 && continue
     scell_vec = get_contribution(a,strian)
     cell_vec, trian = move_contributions(scell_vec,strian)
     @assert ndims(eltype(cell_vec)) == 1
@@ -497,13 +498,15 @@ function _pair_contribution_when_possible(biform,liform)
   mat = DomainContribution()
   vec = DomainContribution()
   for (trian,t) in biform.dict
+    num_cells(trian) == 0 && continue
     if haskey(liform.dict,trian)
       matvec.dict[trian] = pair_arrays(t,liform.dict[trian])
     else
       mat.dict[trian] = t
     end
   end
   for (trian,t) in liform.dict
+    num_cells(trian) == 0 && continue
     if ! haskey(biform.dict,trian)
       vec.dict[trian] = t
     end
@@ -516,11 +519,13 @@ function _pair_contribution_when_possible(biform,liform,uhd)
   matvec = DomainContribution()
   mat = DomainContribution()
   for (trian,t) in _matvec.dict
+    num_cells(trian) == 0 && continue
     cellvals = get_cell_dof_values(uhd,trian)
     cellmask = get_cell_is_dirichlet(uhd,trian)
     matvec.dict[trian] = attach_dirichlet(t,cellvals,cellmask)
   end
   for (trian,t) in _mat.dict
+    num_cells(trian) == 0 && continue
     cellvals = get_cell_dof_values(uhd,trian)
     cellmask = get_cell_is_dirichlet(uhd,trian)
     matvec.dict[trian] = attach_dirichlet(t,cellvals,cellmask)

src/Fields/FieldArrays.jl

@@ -277,7 +277,7 @@ end
 
 function evaluate!(cache,k::LinearCombinationMap{<:Integer},v::AbstractArray,fx::AbstractVector)
   z = zero(return_type(outer,testitem(fx),testitem(v)))
-  @check length(fx) == size(v,1)
+  length(fx) != size(v,1) && return z
   @inbounds for i in eachindex(fx)
     # We need to do the product in this way
     # so that the gradient also works
@@ -579,6 +579,7 @@ function evaluate!(
   nj = size(b,3)
   setsize!(cache,(np,ni,nj))
   r = cache.array
+  size(b,1) == 0 && return r
   for p in 1:np
     for j in 1:nj
       bpj = b[p,1,j]

src/Fields/FieldsInterfaces.jl

@@ -496,6 +496,15 @@
 
 struct IntegrationMap <: Map end
 
+function return_value(k::IntegrationMap,ax::AbstractVector,w)
+  if size(ax,1) == length(w)
+    evaluate(k,ax,w)
+  else
+    T = typeof( testitem(ax)*testitem(w) + testitem(ax)*testitem(w) )
+    zero(T)
+  end
+end
+
 function evaluate!(cache,k::IntegrationMap,ax::AbstractVector,w)
   T = typeof( testitem(ax)*testitem(w) + testitem(ax)*testitem(w) )
   z = zero(T)
@@ -518,6 +527,15 @@
   z
 end
 
+function return_value(k::IntegrationMap,ax::AbstractArray,w)
+  if size(ax,1) == length(w)
+    evaluate(k,ax,w)
+  else
+    c = return_cache(k,ax,w)
+    c.array
+  end
+end
+
 function return_cache(k::IntegrationMap,ax::AbstractArray,w)
   T = typeof( testitem(ax)*testitem(w) + testitem(ax)*testitem(w) )
   r = zeros(T,size(ax)[2:end])
@@ -680,7 +698,11 @@
 function return_cache(f::VoidField,x::AbstractVector{<:Point})
   c = return_cache(f.field,x)
   fx = evaluate!(c,f.field,x)
-  z = similar(fx)
+  if(eltype(fx) != Any)
+    z = similar(fx)
+  else
+    z = Float64[]
+  end
   c, CachedArray(z)
 end