# Approximating the gradient of a Finite Element solution on nodes

I have been working on a Finite Element implementation that approximates the solution to the following PDE in 3D using tetrahedral elements and piecewise linear basis functions.

$$$$\Delta \phi = \nabla\cdot \mathbf{M}$$$$

where $$\mathbf{H} = -\nabla\phi$$, with a nonlinear relation between $$\mathbf{H}$$ and $$\mathbf{M}$$ (Hysteresis). After applying Dirichlet boundary conditions, I believe I succeeded in deriving a mesh-dependent divergence matrix, such that the discretized right hand side of the above equation can be written as $$F = D_M M$$where $$M$$ is now a vector of length $$3N_n$$, where $$N_n$$ is the number of mesh nodes, and $$D_M$$ is a $$N_n\times 3N_n$$ sparse matrix. This can all be done 'relatively rigorously' using the weak formulation of the above PDE.

However, I am now attempting to derive a similar matrix in order to be able to write $$\mathbf{H} = -G\phi$$ where now $$\phi$$ is the FEM solution to the above equation, and $$H$$ is a vector of size $$3N_n$$ containing the three components of the magnetic field at all nodal points. Now, the weak formulation of the above PDE cannot be used (I think), which requires another approach.

So far, I have come up with this. Neglecting boundary condition considerations for now, the solution of the PDE above is approximated by weighted sum of the chosen standard piecewise linear basis functions: $$\phi(\mathbf{x}) \approx \sum_{j=1}^{N_n}\phi_j v_j(\mathbf{x})$$ Consider node $$\mathbf{x}_i$$. Then \begin{align} -\mathbf{H}(\mathbf{x}_j) &= \nabla\phi(x_j)\\ &= \sum_{j=1}^{N_n}\phi_j \nabla v_j(\mathbf{x}_i)\\ &= \sum_{j\in n_i} \phi_j \nabla v_j(\mathbf{x}_i)\\ &= \phi_i \nabla v_i(\mathbf{x}_i) + \sum_{j\in n_i\setminus i} \phi_j \nabla v_j(\mathbf{x}_i)\\ \end{align} where $$n_i := i\cup \{\text{node indices, attached to node x_i}\}$$, and where the last two steps follow because the support of each basis function is the convex set with as boundary points, the nodes of set $$n_i$$.

I know that the gradient of the FEM solution is piecewise constant on elements and formally does not exist on nodes. One can try a (weighted) average of the surrounding element gradients perhaps, to obtain an approximation of the gradient on nodes? Does anyone have any ideas on how to find the best approach to compute the 'discretized gradient' matrix $$G$$ ?