Gradient of basic functions evaluated in each point of the mesh

For the \phi_i's, the basis functions obtained by the approximation of a P1 (linear) finite element space related to the discretization of a cylinder \Omega,

The setup is as follows:

r=15
large=70
begin=[-15, 0, 0]
ending=[45, 0, 0]
cylinder=Cylinder(Point(begin),Point(ending), r ,r)
mesh_cilinder = generate_mesh(cylinder, 10)
Space = FunctionSpace(mesh_cilinder,'P', 1)
basis_function=TrialFunction(Space)

I am trying to compute the gradients of the basis functions evaluated at each vertex of the mesh. Specifically, I want to compute the following matrices: \big[\nabla_x \phi_i(x_j)\big]_{i,j}, \big[\nabla_y \phi_i(x_j)\big]_{i,j}, \big[\nabla_z \phi_i(x_j)\big]_{i,j}, where x_j are the coordinates of the mesh vertices.

My approach was to create an expression for each basis function and then interpolate it over the entire function space:

coordinates_nodes=mesh_cilinder.coordinates()
cells_nodes=mesh_cilinder.cells()
Space = FunctionSpace(mesh_cilinder,'P', 1)
basis_function=TrialFunction(Space)

basis_function_list=[ ]
integral_basis=[ ]
n=Space.dim()


for i in range(n):
    # Create a function that represents the i-th basis function   
    x0=coordinates_nodes[i][0]
    x1=coordinates_nodes[i][1]
    x2=coordinates_nodes[i][2]
    phi_node = interpolate(Expression('(x[0] == x0 && x[1] == x1 && x[2] == x2) ? 1 : 0', degree=3, x0=x0,x1=x1,x2=x2),Space)
    basis_function_list.append(phi_node)

and after that compute the gradient and evaluate in the vertex of the mesh, for example for phi_0

u0=interpolate(basis_function_list[0],Space)
gradf   = grad(u0)
u,v,w   = gradf
project(u).vector().get_local()
project(v).vector().get_local()
project(w).vector().get_local()

However, the values I obtained for the gradients are not what I expected, since they are not constant, however the derivatives of the basis functions should be piecewise constant since I am using P1.

Could you please help me understand the correct way to compute the gradients of the basis functions and evaluate them at the mesh vertices for computing the matrices?

Please read https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21: on top of all written there, it’s very difficult to read your code, because you haven’t formatted it properly with “```”.

From what I can gather, the space in which you want to apply interpolate is a DG 0 space (and not the P 1 space you define), because, as you say, the gradient of a P 1 function is piecewise constant function.

1 Like

Thanks, for your answer it worked, however it takes a lot of time when the size of discretization is small. might you have another idea for computing that matrices efficiently?

Switch to DOLFINx and use dolfinx.fem.Expression to compute it, as you can interpolate the gradient onto a compatible DG-0 or DG-1 space and easily extract the values.