Gridap.jl
Gridap.jl copied to clipboard
Laplacian of product of `CellField`s doesn't work.
The Laplacian of product of CellField
s (formed using the same FESpace
and Triangulation
) fails. I came across the following,
https://github.com/gridap/Gridap.jl/blob/672e3d9027246bfa677c07b4e9fabb9e1fde9eb4/src/Fields/FieldsInterfaces.jl#L94-L102
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)
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
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.