I think I found a way to fix the code. On the line where Gauss' law is used:
vol = vol + abs(P(1) * N(1) * a);
The absolute value there is not correct and should be removed - when the vector field and normal vector are antiparallel the addition to the volume needs to be negative. With the current implementation you can change the volume by shifting the X coordinate. If you get rid of the absolute value everything is fine, so long as you can ensure that the normals on the surface point OUT. Although the addition of the absolute value seems intuitive, it is not how Gauss's law works.
The normals can be computed with another code on the filie exchange: computeNormalVectorTriangulation by David Gingras. David's code will give you normals that point outwards so long as the triangles are defined in a certain way... Refer to his code.
This update seems to give the right volume and remove the dependence on translation of vertices, at least in my testing.
This is an excellent file that does exactly what it says it does. It worked the first time. I don't know why someone else gave it a bad rating - I think it's perfect as is.
