Error when evaluating fine-grid DoFs on coarse FEFunction or FEBasis

[amartinhuertas]

The following MWE:

using Gridap
using Gridap.ReferenceFEs
reffe=LagrangianRefFE(Float64,QUAD,1)
modelH=CartesianDiscreteModel((0,1,0,1),(1,1))
modelh=refine(modelH,2)
XH = TestFESpace(modelH,reffe)
Xh = TestFESpace(modelh,reffe)
xH = get_fe_basis(XH)
uH = zero(XH)
σh = Gridap.FESpaces.get_fe_dof_basis(Xh)
σh(xH) # Fails 
σh(uH) # Fails

produces the error below, and I think it should not. I did not take a careful look at the core of Gridap, but it seems there might be a missing change_domain ?

ERROR: 

A CellDof can only be evaluated on a CellField defined on the same Triangulation.

Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] macro expansion
   @ ~/git-repos/Gridap.jl/src/Helpers/Macros.jl:47 [inlined]
 [3] evaluate!(cache::Nothing, s::Gridap.CellData.CellDof{ReferenceDomain}, f::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}})
   @ Gridap.CellData ~/git-repos/Gridap.jl/src/CellData/CellDofs.jl:108
 [4] evaluate
   @ ~/git-repos/Gridap.jl/src/Arrays/Maps.jl:87 [inlined]
 [5] (::Gridap.CellData.CellDof{ReferenceDomain})(f::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}})
   @ Gridap.CellData ~/git-repos/Gridap.jl/src/CellData/CellDofs.jl:100
 [6] top-level scope
   @ ~/git-repos/Gridap.jl/mwe.jl:12

[JordiManyer]

Hey @amartinhuertas ! After having a look, don't you think this is a feature and not a bug?
This is the code that is run when evaluating a CellField on a CellDof:

function evaluate!(cache,s::CellDof,f::CellField)

  trian_f = get_triangulation(f)
  trian_s = get_triangulation(s)

  if trian_f !== trian_s
    @unreachable """\n
    A CellDof can only be evaluated on a CellField defined on the same Triangulation.
    """
  end

  b = change_domain(f,s.domain_style)
  lazy_map(evaluate,get_data(s),get_data(b))
end

From what I see, there is no change_domain whatsoever invoked here, which means this is a design choice. The same routine would fail if you tried evaluating a CellField defined on the body-fitted triangulation on a CellDof defined on the boundary (without any refinement).

I think the domain change is meant to be done earlier during the interpolation methods (i.e here), but not here.

What do you think?

[amartinhuertas]

What do you think?

I agree. It seems a design choice. However, I am not totally sure which is the final motivation behind such design choice.
Is this because one does not want to mix triangulations of different Dc values ? (just guessing ... I do not have an answer). @fverdugo can perhaps help us and say some words here?

Anyway, whichever the requirements/design choices were at the time this function was developed, these are not set in stone and they may involve in time. After we have introduced adaptivity, we, can, e.g. do:

uHh=interpolate(uH,Uh) 

which is actually the same as $\sigma_h(u_H)$. Thus, I see something inconsistent here, why do we allow uHh=interpolate(uH,Uh) and we dont allow $\sigma_h(u_H)$?