Calculating gradient and how to access x, y components of it

General
1

I am calculating gradient of the following expression whose grad is (2x,2y) but I am getting a vector of ones as shown below. I know this is very basic but I am confused what is wrong, I am also confused about dim of x component of grad calculated on 1 by 1 grid. After calculating gradient I want to access both x and y components of the gradient. Here is the minimal example

 # define mesh and function space
msh = create_unit_square(MPI.COMM_WORLD,1,1)
V = functionspace(msh, ("Lagrange", 1))
x = SpatialCoordinate(msh)
u = Function(V)
u.interpolate(lambda x: x[0]**2 + x[1]**2)
gdim = msh.geometry.dim
dg = dolfinx.fem.functionspace(msh,('DG',0,(gdim,)))
grad_int = dolfinx.fem.Function(dg)
grad_uh = dolfinx.fem.Expression(grad(u),dg.element.interpolation_points())
grad_int.interpolate(grad_uh)
grad_int.x.array

this is the output
image

2

Consider that the spaces you’re employing are not rich enough to appropriately represent your exact functions on such a coarse mesh. See the following adjustment to your code such that you can compare your approximate gradient with the exact gradient.

In this example I’ve filled out the spaces such that each function is exactly represented; however, you can play with these spaces to see the implications of inexact representation of functions.

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

msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD,1,1)
x = ufl.SpatialCoordinate(msh)

# Exact function
V = dolfinx.fem.functionspace(msh, ("Lagrange", 2))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0]**2 + x[1]**2)

# Approximate gradient
dg = dolfinx.fem.functionspace(msh,('DG',1,(msh.geometry.dim,)))
grad_uh = dolfinx.fem.Function(dg)
grad_uh.interpolate(dolfinx.fem.Expression(ufl.grad(u),dg.element.interpolation_points()))

# Exact gradient
grad_u = dolfinx.fem.Function(dolfinx.fem.functionspace(msh, ("Lagrange", 1, (msh.geometry.dim,))))
grad_u.interpolate(lambda x: np.stack((2 * x[0], 2 * x[1])))

# Error comparison of exact vs approximate gradient
l2_error = msh.comm.allreduce(dolfinx.fem.assemble_scalar(
	dolfinx.fem.form((grad_uh - grad_u) ** 2 * ufl.dx))) ** 0.5
print(f"L2 error: {l2_error}")

which for me yields:

L2 error: 6.409875621278547e-17
1 Like
3

Thanks for your reply, I see. In addition, I want to access both x, y components of grad as an array I see grad_uh.x.array gives x component how to get y component as well.

4

grad_uh.x.array yields the underlying array of finite element function coefficients wrapped by a numpy array. The x does not mean geometric component.

Adapted from Problem solving Harmonic Oscillator - #7 by nate, consider appending the following to the code to evaluate the finite element functions at positions in space:

from dolfinx import geometry
# Points at which to evaluate gradient
x = np.array([
	[0.0, 0.0, 0.0],
	[1.0, 0.0, 0.0],
	[1.0, 1.0, 0.0],
	[0.0, 1.0, 0.0]])
bbox = geometry.bb_tree(msh, msh.topology.dim)
candidate_cells = geometry.compute_collisions_points(bbox, x)
collided_cells = geometry.compute_colliding_cells(
    msh, candidate_cells, x)
collided_cells = collided_cells.array[collided_cells.offsets[:-1]]

grad_uh_values = grad_uh.eval(x, collided_cells)
print(grad_uh_values)

Which yields the evaluations of the gradient:

[[4.30845614e-17 4.30845614e-17]
 [2.00000000e+00 6.79377410e-17]
 [2.00000000e+00 2.00000000e+00]
 [6.79377410e-17 2.00000000e+00]]
1 Like
5

Thank you Dr. @nate I appreciate your response.