Primal HHO remarks, issues found, to-discuss, etc.

[amartinhuertas]

• We need div(M grad(u)) with M being a polytopal-wise constant tensor, and u being a FEBasis. At present, we can "only" compute laplacian(u), i.e., M=I. See here and here. [NOTE] This need can be actually by-passed using integration by parts, at the expense of having two terms instead of one.
• I have the feeling that LocalFEOperator is becoming a monster struct with many responsibilities (reconstruction operator, reduction operator, L2-projection cell-wise and facet-wise operator, etc.) Perhaps we may consider splitting it into several structs. Just an idea. This requires some thinking.
• Question: The gradient reconstruction space is a higher order space than the one where we seek the solution in. Exact integration of these higher order polynomials requires a higher order quadrature, right? Should we use mixed quadratures for integration? (solution space versus reconstruction space?)
• At present, we are writing the reconstructed gradient bi-linear form term by blocks as:

  function r(u,v)
    uK,u∂K=R(u)
    vK,v∂K=R(v)
    ∫(∇(vK)⋅∇(uK))dΩ + ∫(∇(vK)⋅∇(u∂K))dΩ + ∫(∇(v∂K)⋅∇(uK))dΩ + ∫(∇(v∂K)⋅∇(u∂K))dΩ
  end

where P is the reconstruction operator. This is justified, among others, because we want to do static condensation, but also because in the stabilization bi-linear form term we need the result of the reconstruction operator with a block layout. Is this reasonable or should we think how to do it scalar-wise and then re-block after that?

• Discuss the current design of LocalFEOperator based on the field_type_at_common_faces trait. This is currently being used to accommodate in the same type the reduction operator and the difference operator (the latter required for the stabilization).
• These terms ∫(δv∂K_K*δu∂K_K)d∂K+∫(δv∂K_K*δu∂K_∂K)d∂K+∫(δv∂K_∂K*δu∂K_K)d∂K in the bilinear form stabilization are zero numerically. Are they mathematically zero?
• Discuss the need for the keyword argument vector_type in

  GridapHybrid.jl/test/PoissonHHOTests.jl

  Line 18 in cb54ef2

  VKR_C = TestFESpace(Ω, reffe_c ; conformity=:L2, vector_type=Vector{Float64})