Compute Gradients over FE mesh using Nodal Solution

I have computed the solution to an electrostatics problem using the finite element method in a slightly non-standard workflow (I could not figure out how to accomplish this within the PDE modeling environment, see here):

1. Define finite element mesh outside the PDE modeling environment
2. Create an electrostatics model and import the mesh
3. Use FEM = assembleFEMatrices(model) to assemble the global stiffness matrix
4. Manually prescribe Dirichlet boundary conditions at mesh nodes by constructing an H matrix (not sure what the technical term here is - a condensation operator, a boundary condition applicator, a row eliminator?) and using the pdenullorth(H). By the way - it's strange that this function does not have any documentation that I can find!
5. Solve the FE problem using a backslash operation. Here's steps 3-5 (tetMesh struct contains the tetrahedral mesh I've imported, and catNodes & anNodes identify the nodes where I'd like to prescribe voltage boundary conditions):

%assemble matrices (right now F is always zero)

FEM = assembleFEMatrices(model,"KF");

%create a vector of prescribed voltages

Vd = sparse(tetMesh.nNodes,1);

Vd(tetMesh.catNodes) = 1000;

Vd(tetMesh.anNodes) = 0;

%form the consdensing matrix

H = sparse(1:tetMesh.nDirichlet,[tetMesh.catNodes tetMesh.anNodes],...

ones(tetMesh.nDirichlet,1),tetMesh.nDirichlet,tetMesh.nNodes);

[B,~] = pdenullorth(H);

% form the reduced matrix algebra problem

Kc = B'*FEM.K*B;

Fc = B'*(FEM.F-(FEM.K*Vd));

%compute the FE solution

u = B*(Kc\Fc) + Vd;

This works fine - I can now plot the nodal solutions (the electric potential) over the problem domain:

Now, I would like to compute the derived quantities associated with normal electrostatics solutions - namely the Electric Field, which is simply the gradient of the Electric Potential Field.

Ordinarily result = solvepde(model) would compute this quantity automatically - but I computed my solution simply using \. Is there a way to perform this using the FE mesh that the problem was solved on? I have a hunch that the shape functions and their derivatives embedded in the FE mesh will produce a more accurate gradient field than if I project the solution to a uniform grid, compute a gradient, and project back. For this reason I would like to avoid using scatteredInterpolant().