I am working on implementing a PDE model that simulates a certain physical phenomenon on the surface of a 3D mesh. The model involves calculating mixed partial derivatives of a scalar function defined on the vertices of the mesh. What I tried so far (which is not giving good results), is this: 1. First, I calculate the gradient of the scalar function over the mesh. This gives me the first derivative in the x, y, and z directions (the gradient vector field over the mesh). 2. Since the gradient field is defined per face (as opposed to per vertex), I interpolate the gradient field to the vertices by performing a weighted averaging over the triangles that share a vertex. I am using the areas of the triangles as the weights. 3. After this interpolation, the gradient (first order derivative is now define "interpolated" over the vertices). I use this interpolated vector field as input to the gradient operator again and get the gradient of each component in the vector. I know that this is incorrect (or at least inaccurate) since the scalar function originally defined over the mesh goes out the range of values it should be restricted to if the model was implemented correctly. So, how can I calculate an accurate descritized approximation of mixed partial derivatives over a 3D mesh?