Ve = VectorElement("CG", mesh.ufl_cell(), 1)
Re = VectorElement("R", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10))
V = FunctionSpace(mesh, Ve)
v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W)
w = Function(W)
dx = Measure('dx')(subdomain_data=subdomains)
Eps = Constant(((1, 0), (0, 0)))
F = sum([inner(sigma(dv, i, Eps), eps(v_))*dx(i) for i in range(nphases)])
a, L = lhs(F), rhs(F)
a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx
solve(a == L, w, [])

I wanted shape function derivatives w.r.t spatial coordinates (x,y) at each cell (triangle) of the mesh.
This means dN1/dx, dN2/dx; dN1/dy, dN2/dy data for each cell in the mesh.

Can I get a list of shape function derivatives of each cell ? /
is there any way to use function w for computing shape function derivatives like project function ?.

Note that shape derivatives are not computer per cell, but for each mesh node (usually in fenics a vector function space of first order Lagrange elements).

Thanks @dokken for your response and sharing the reference. I have doubts regarding your statement for above case.

“Note that shape derivatives are not computer per cell, but for each mesh node (usually in fenics a vector function space of first order Lagrange elements).”

Does this mean that, at each node we get dN1/dx, dN2/dx, dN3/dx; dn1/dy, dN2/dy , dN3/dy (6X1) evaluated at each node

For entire mesh, would it be mesh.num_vertices() X 6 scalar matrix ?

would shape function derivative be a scalar of 6X1 dimension for each node (after using w.compute_vertex_values )?

I would suggest you look at the references I have provided. The first one:

goes into detail about the mathematics, and the expected input/output, while the second reference goes more into the usage of this derivative for optimization problems

The resources you suggested shows Automatic differentiation algorithm to find change in functional w.r.t. some variable and looks more of time dependent application. I wanted to compute the derivatives of shape functions w.r.t x, y (for an element, evaluated through guass points). I don’t find way I can interpret the shape functions as it is formed through argument (test and trial functions in backend) and therefore, couldn’t find way to minimize the basis or shape functions.

I wanted to use the shape function derivatives w.r.t spatial coordinates for each element like
dN1/dx_element= (Summation dN1/dx_i * weight_i) , i=1… no. of guass points.

I wanted the direction related to my application. Kindly help me.

I dont understand what you mean by this.
A shape derivative for me is the change of an integral when you perturb the integration domain
This can be achieved in FEniCS by.

u = TrialFunction(V)
v = TestFunction(V)
a = inner(u,v)*dx
x= SpatialCoordinate(mesh)
dadx =derivative(a,x)

Based on the requirement, shape derivative means derivative of lagrangian basis functions used in general finite element formulation. \:u=\:\sum _{i=1}^Nu_i\:\phi i\: , N is the number of nodes. \phi -lagrangian shape/basis function.

This u is used as Function (Coefficient) u=Function(W) and \phi as argument in UFL language. dv= TrialFunction(W)

Here, \phi is the lagrangian shape functions and the derivative which I want is \frac{d}{dx}\left(\phi \right)

Can you tell me how this derivative is elated to perturbation of integral which you described? I apologize for any mistakes I committed during explanation.