Dear community. I would like to know if there is way to compute the nodal forces on the nodes of one finite element that mechanically balance the mechanical deformation of such an element, i.e., such that the virtual power of such nodal forces equals the strain energy of the element. In other words, I wonder how I could implement the following integral,

\mathbf{f}^{n}=\int_{\mathcal{B}e}\sigma^t \cdot\nabla N^n dV

where \nabla N^n dV Is the gradient of the shape function of the node n; and \mathcal{B}e the volume domain of the single element.

To this end, I guess that I need a way to locate an individual finite element within a larger mesh, to obtain the derivative of the nodal shape functions, and to compute the integral on that element. Then, do the same for the rest of nodes of the element.

Thank you, dokken.
How can I get the displacement vector uh?
I guess it must contain only the solved nodal degrees of freedom. I have tried u.vector().array() using FEniCSx, but I get the error “‘petsc4py.PETSc.Vec’ object is not callable”.
Also, which Function from dolfinx.fem.assemble should I use?

No, I work with finite strains and the deformation gradient (\mathbf{F}). More precisely, it is not the Cauchy stress tensor the one I want to use, but the Eshelby one, but this does not affect the implementation. I define u as a Function from a Function Space, and then, \mathbf{F} = \mathbf{I} + \text{ufl.grad}(\mathbf{u}). With this implementation, I aim to compute the material (configurational) nodal forces.
Hence, is this u solution function what I should use in the assemble you proposed?

Since you wrote uh, I am not sure if that is simply my Function \mathbf{u} already solved, or just the nodal d.o.f.s ordered in a different manner.

That for the assemble works, thanks. However, with print(f.get_local()) I get the following error: AttributeError: 'petsc4py.PETSc.Vec' object has no attribute 'get_local'

As you didn’t post a u in your initial code, I used uh to illustrate the numerical displacement function (stemming from one process or another). uh should therefore be replaced with your u

This seems to work. Thank you for the explanations.

Then, the only thing I still need to wrap up this implementation would be to get a vector with the coordinates of the nodes for which the nodal forces in the vector f.array correspond to. I would appreciate your help again here.

One additional concern related to the nodal forces:
How could I write the forces on the nodes, i.e., the vector \mathbf{f}, to my .xdmf file to visualize them on Paraview?