How to find electric current density from electric potential in 3d finite element

Hello,
The value of electric potential(∅) is known at every node in a 3d finite element mesh. The relation between electric current density(i) and electric potential(∅) is i=k.∇∅, I am writing a code in c, I want to know how to find the gradient of electric potential(∇∅) at every node so as to get current density(i). Please help.

The OP said a finite element mesh. The previous answer looks more appropriate for a regular finite differnce mesh, not an FE mesh with arbitrary geometry.

You can calculate the gradient at each node of each element using the element shape functions. (This will be similar, but simpler, than calculating strains from displacements in a structural FE program).

The problem is that the only output that is consistent with the FE method is gradient averaged over the volume of each element. When several elements meet at a node, you will get a different gradient value at the node for each of the elements.

If you just want to plot the data you can do something fairly simple, for example average of the different nodal values (possibly weighted by the volume of the elements). Or you can use the difference between the element values at each node as a measure of the accuracy of the solution.

If you want to do some mathematical post processing using the gradients, it would be better to do it based on the value within the volume of each element, rather than some more or less arbitrary nodal averaging method.

Method suggested by uart is correct, but it suits a regular structured rectangular mesh. Because I am dealing with unstructured mesh, I would prefer to use the method proposed by AlphaZero. Thanks both for their inputs.
Neverthless in strict sense ∇∅ is an elemental solution and we need to use some averaging method to calculate its value at nodes.

Neverthless in strict sense ∇∅ is an elemental solution and we need to use some averaging method to calculate its value at nodes.

I would rephrase that as "the FE approximation for ∇∅ is discontinuous across the element boundaries". Whether you "need" to get rid of the discontinuities depends what you want to use the solution for.

If you really don't want discontinuities, you could formulate an element with a different variational principle so the gradients are nodal variables. That has been done for structural analysis (but few people actually use those element formulations). I don't have any hands-on experiemnce of FE methods in electromagnetism and I don't know what the state of the art is in that field.