Hi Elizabeth,
the problem is that the higher order degrees of freedom for the H1(...order=2) are not related to your nodes in the mesh. H1 uses a hierarchical basis, where the shape functions are related to mesh edges, faces and cells. Their dof-numbers are numbers after total number of mesh-points.
Every second of your points (in each directions) are not used for placing unknowns, thus the dof-inconsistency warning.
I recommend to import your griddata as a VoxelCoefficientFunction, and then interpolate it to a GridFunction via
We cannot form the gradient directly of a voxelcf (at least not a the moment), but you can use the gradient of the GridFunction you interpolated to.
VoxelCF can use (tri-)linear interpolation by now, but higher order interpolation could be added with little coding.
Best,
Joachim