Laplacian of product of CellFields doesn't work.

[kishore-nori]

The Laplacian of product of CellFields (formed using the same FESpace and Triangulation) fails. I came across the following,

Gridap.jl/src/Fields/FieldsInterfaces.jl

Lines 94 to 102 in 672e3d9

	function push_∇∇(∇∇a::Field,ϕ::Field)
	@notimplemented """\n
	Second order derivatives of quantities defined in the reference domain not implemented yet.
	This is a feature that we want to have at some point in Gridap.
	If you are ready to help with this implementation, please contact the
	Gridap administrators.
	"""
	end

But the error doesn't come from here, the evaluation doesn't hit the above, at least from preliminary debugging observation. The following is the MWE:

using Gridap
using Test

domain = (0.,1.,0.,1.)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

g(x) = exp(-norm(x))
reffe = ReferenceFE(lagrangian,Float64,2)
V = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary")
U = TrialFESpace(V,g)

n = num_free_dofs(U)
uh = FEFunction(U,rand(n))
qh = FEFunction(U,rand(n))

p = Point(0.1,0.2)

cf1 = ∇(uh + qh)
@show @test cf1(p) ≈ ∇(uh)(p) + ∇(qh)(p) # works 

cf2 = ∇(uh*qh)
@show @test cf2(p) ≈ (qh*∇(uh))(p) + (uh*∇(qh))(p) # works 

cf3 = Δ(uh + qh) # equiv to tr(∇∇(⋅)), making it an OperationCellField
@show @test cf3(p) ≈ Δ(uh)(p) + Δ(qh)(p) # works 

cf4 = ∇∇(uh + qh)
@show @test cf4(p) ≈ ∇∇(uh)(p) + ∇∇(qh)(p) # works 

# Although the CellField is built, when evaluated it results in error
cf5 = ∇∇(uh*qh) # works 
cf5(p) # fails 

cf6 = Δ(uh*qh) # even the construction fails 

I would be happy to fix the related issue and implement the second order derivatives based on your guidance and suggestions, as the @notimplemented message directs.

cc @amartinhuertas, @oriolcg (I think latest commits along these lines where your contributions)

[carlodev]

I don't know if this is completely correct. I looked also as push_∇(∇a::Field,ϕ::Field) = pinvJt(∇(ϕ))⋅∇a:

function push_∇∇(∇∇a::Field,ϕ::Field)
  s = pinvJt(∇(ϕ))
  tr(s⋅s⋅∇∇a)
end

[kishore-nori]

Hi @carlodev, thank you for the possible solution. I'll work it out on paper and get back. At first sight it seems to be alright to me, although I am not sure on whether the above is only for a ϕ which is affine linear. An other thing is, like I mentioned above, the error doesn't seem to be from the unavailability of above, at least from my very preliminary debugging experience, rather, it seems to me that the higher order chain rules are unavailable as well. I am in middle of something else right now, will start working on it after 10 days and get back here, sorry. For now, until the above is fully resolved, we will have to rely on manual product rule expansion of such operators.