Point interpolation of shape function derivatives

[SouzaEM]

Is it possible to calculate the values of the shape function derivatives using vertex meshes? As suggested in #3782 and from the FWI example (https://www.firedrakeproject.org/demos/full_waveform_inversion.py.html), it is possible to obtain the shape function values as follows:

from firedrake import *

mesh = UnitSquareMesh(2, 2)

x_v = [np.r_[0.25, 0.25]]
mesh_v = VertexOnlyMesh(mesh, x_v)

V_v = FunctionSpace(mesh_v, "DG", 0)
phi_v = assemble(TestFunction(V_v)*dx)

V = FunctionSpace(mesh, "CG", 1)
phi = Cofunction(V.dual()).interpolate(phi_v)
print(phi.dat.data)

However, as the function spaces defined on vertex meshes may only be DG0, I was not able to assemble the gradient of the test functions.

[Ig-dolci]

Do you want a Cofunction that is a result of the assembly of shape function gradient in DG0? Then interpolate it into V space? This does not look like the right thing to do once DG0 has dimension 0 and a constant function in DG0 does not have a well-defined gradient. Or I am wrong and it is possible...

[SouzaEM]

Do you suggest something like this?

from firedrake import *

mesh = UnitSquareMesh(2, 2)

V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)

from firedrake.__future__ import interpolate

W = VectorFunctionSpace(mesh, "DG", 0)
w = assemble(interpolate(grad(v), W))

Then how do I integrate this code with point interpolation? I need the shape function derivatives in few points of the domain.

[Ig-dolci]

I think we cannot compute Cofunction gradients with respect to the coordinates. However, I have a workaround that may help:

phi_primal = phi.riesz_representation(riesz_map='l2')
grad_phi = grad(phi_primal)
grad_phi_v0 = assemble(interpolate(grad_phi[0], V_v) * dx)
grad_phi_v1 = assemble(interpolate(grad_phi[1], V_v) * dx)

[SouzaEM]

I tried running your suggestion with the following:

from firedrake import *

mesh = UnitSquareMesh(2, 2)

x_v = [np.r_[0, 0]]
mesh_v = VertexOnlyMesh(mesh, x_v)

V_v = FunctionSpace(mesh_v, "DG", 0)

V = FunctionSpace(mesh, "CG", 1)
phi = Cofunction(V.dual())

phi_primal = phi.riesz_representation(riesz_map='l2')
grad_phi = grad(phi_primal)
grad_phi_v0 = assemble(interpolate(grad_phi[0], V_v) * dx)
grad_phi_v1 = assemble(interpolate(grad_phi[1], V_v) * dx)

print(grad_phi_v0)
print(grad_phi_v1)

I changed the reference point to (0, 0) to make sure the derivative of the shape functions would not be zero at the reference point. However, the result is zero for grad_phi_v0 and grad_phi_v1. Also, the output is only one value instead of one value per shape function. Did I do something wrong?

[Ig-dolci]

Yes. See the following code:

from firedrake import *

mesh = UnitSquareMesh(2, 2)

x_v = [np.r_[0, 0]]
mesh_v = VertexOnlyMesh(mesh, x_v)

V_v = FunctionSpace(mesh_v, "DG", 0)
phi_v = assemble(TestFunction(V_v)*dx)

V = FunctionSpace(mesh, "CG", 1)
phi = Cofunction(V.dual()).interpolate(phi_v)
phi_primal = phi.riesz_representation(riesz_map='l2')
grad_phi = grad(phi_primal)
grad_phi_v0 = assemble(interpolate(grad_phi[0], V_v) * dx)
grad_phi_v1 = assemble(interpolate(grad_phi[1], V_v) * dx)
print(grad_phi_v0, grad_phi_v1)

[SouzaEM]

Thank you! This is really what I need. I will just add guidance for obtaining the derivatives of the other shape functions. Your code is returning the (-2, -2), which are the derivatives for the first shape function:

$\phi_1(x, y) = 1 - 2x - 2y$
$\frac{\partial\phi_1}{\partial x}(x, y) = -2$
$\frac{\partial\phi_1}{\partial y}(x, y) = -2$

For the second and third shape functions we have:
$\phi_2(x, y) = 2y$
$\frac{\partial\phi_2}{\partial x}(x, y) = 0$
$\frac{\partial\phi_2}{\partial y}(x, y) = 2$

and:
$\phi_3(x, y) = 2x$
$\frac{\partial\phi_3}{\partial x}(x, y) = 2$
$\frac{\partial\phi_3}{\partial y}(x, y) = 0$

Their derivatives may be obtained by appending the following lines to the end of your code:

# For Φ_2
phi.dat.data[0] = 0
phi.dat.data[1] = 1
grad_phi_v0 = assemble(interpolate(grad_phi[0], V_v) * dx)
grad_phi_v1 = assemble(interpolate(grad_phi[1], V_v) * dx)
print(grad_phi_v0, grad_phi_v1)

# For Φ_3
phi.dat.data[1] = 0
phi.dat.data[2] = 1
grad_phi_v0 = assemble(interpolate(grad_phi[0], V_v) * dx)
grad_phi_v1 = assemble(interpolate(grad_phi[1], V_v) * dx)
print(grad_phi_v0, grad_phi_v1)

[SouzaEM]

I will try to be more specific in what I'm trying to achieve with the point interpolation feature. I'm trying to obtain the following functionality in the global mesh:

from FIAT.lagrange import Lagrange as CG
from FIAT.reference_element import UFCTriangle

p_ref = [[0.25, 0.25]]  # Reference point
degree = 1              # Degree of the shape function
order = 1               # Order of the derivative

T = UFCTriangle()
element = CG(T, degree)
tab = element.tabulate(order, p_ref)

print(f"    \u03A6(0.25, 0.25) = {tab[(0,0)].transpose()}")
print(f"d\u03A6/dx(0.25, 0.25) = {tab[(1,0)].transpose()}")
print(f"d\u03A6/dy(0.25, 0.25) = {tab[(0,1)].transpose()}")

which gives the following for the reference triangle:

    Φ(0.25, 0.25) = [[0.5  0.25 0.25]]
dΦ/dx(0.25, 0.25) = [[-1.  1.  0.]]
dΦ/dy(0.25, 0.25) = [[-1.00000000e+00  5.55111512e-17  1.00000000e+00]]

The motivation for this is to implement moment tensor sources in elastic wave propagation with a displacement formulation. This requires the values of dΦ/dx and dΦ/dy. By studying the example of FWI, I learned to obtain Φ with the point interpolation feature, but I could not figure out how to calculate dΦ/dx and dΦ/dy.