How compute the global stiffness matrix derivatives with respect to cell coeeficients

Hi all! How compute the global stiffness matrix derivatives with respect to cell coeeficients?
Consider a bilinear form:

a(u, v)=\int_{\Omega}d_h^p\mathbb{\sigma}:\mathbb{\epsilon}d\Omega
where d_h is a cell-wise function, h\in\mathcal{T} is the cell of the fem domain \mathcal{T}, p is a constant, \sigma is the stress tensor and \epsilon is the strain tensor.

Assemble the a(u,v), the stiffness matrix K is obtained.
Now my question is how to compute the derivative of K w.r.t. d_h using dolfinx and ufl.

You can use dadp = ufl.derivative(a, d) to get the derivative with respect to a coefficient d in your variational form.

However, please note that if a is bi-linear, then this becomes a tri-linear form, something that is extremely memory intensive to assemble.

It would help if you added some context regarding what you would like to achieve with this derivative,
and some accompanying code.