Computing element-wise spatial derivatives of the solution

General
4

Projecting a linear function into a quadratic function space will still render the function linear in time.

Imagine that you have uh = x
Then up = x, and thus the double derivative is still zero.
This can be illustrated with the following:

import dolfinx
from mpi4py import MPI
import numpy as np
import ufl

cell = ufl.Cell("triangle", geometric_dimension=2)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, 1))
x = np.array([[0., 0.], [1., 0.], [0., 1.]])
cells = np.array([[0, 1, 2]], dtype=np.int32)
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, x, domain)

V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
W = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 2))

u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0]+3*x[1])

w = dolfinx.fem.Function(W)

uh = ufl.TrialFunction(W)
vh = ufl.TestFunction(W)
a = ufl.inner(uh, vh) * ufl.dx
L = ufl.inner(u, vh) * ufl.dx
problem = dolfinx.fem.petsc.LinearProblem(a, L, u=w)
problem.solve()

print(w.x.array)

when we get

[-2.22478286e-14  1.00000000e+00  3.00000000e+00  2.00000000e+00
  1.50000000e+00  5.00000000e-01]

which matches the ordering of

i.e. the function can be expressed as:

\begin{align} u &= \sum_{i=0}^{5} c_i \phi_i\\ & = 0 \cdot (2x^2+4xy-3x+2y^2-3y+1) \\ &+ 1 \cdot (x (2x-1)) + 3 \cdot ( y (2y-1))\\ &+ 2\cdot (4xy)\\ &+ 1.5 \cdot (4y (-x-y+1))\\ &+ 0.5 \cdot ( 4x(-x-y+1))\\ &= 0\\ &+2x^2-x + 6y^2 - 3y\\ &+8xy\\ &-6xy- 6y^2 + 6y\\ & - 2x^2 - 2xy +2x\\ &= x+3y \end{align}

and thus also the double derivative of this function will be zero.