Understanding the difference between Eulerian and Lagrangian view point w.r.t to the Material Derivative: Consider a situation where we have a particle system <-- Lagrangian View Point Also, consider that we are trying to find out the temperature at a point in the system, where the particle system lives <- Eulerian View point Now, the we say that the material derivative of a particle in the system is zero, what does it mean ? Consider an example, in which we are trying to find out the temperature a point (say 0,0,0) The question is, what is the temperature at this point ? Now if we say that the temperature at this point (0,0,0 -> Eulerian) is given by the temperature of the particles (<- Lagrangian) at that moment of time, at that point (0,0,0), then if the Material Derivative is zero then, the temperature of the particle (at 0,0,0 or any other point in space) does not change, but if the particles are moving, then the temperature at a point (0,0,0) will change, as a new particles will replace the old particles to give us a new temperature at that point ! Trilinear Interpolation Another important thing to keep in my while making your fluid solver is Trilinear Interpolation Following Wikipedia, I would get the following pseudo-code for trilinear interpolation. Say Xg, Yg, Zg is a the points on the Grid, while Xp, Yp, Zp is the point whose value (say density) is to be found. Then following this approach, we get the following: voxFaceX = (int) Xp; voxFaceY = (int )Yp; voxFaceZ = (int) Zp; Xd = Xp - voxFaceX Yd = Yp - voxFaceY Zd = Zp - voxFaceZ First we interpolate along the z- axis i1 = ( v(voxFaceX, voxFaceY, voxFaceZ) * ( 1 - voxFaceZ) ) + (v(voxFaceX, voxFaceY, voxFaceZ +1) * voxFaceZ ) i2 = ( v(voxFaceX, voxFaceY + 1, voxFaceZ) * ( 1 - voxFaceZ) ) + (v(voxFaceX, voxFaceY +1 , voxFaceZ +1) * voxFaceZ ) j1 = ( v(voxFaceX + 1, voxFaceY, voxFaceZ) * ( 1 - voxFaceZ) ) + (v(voxFaceX + 1, voxFaceY, voxFaceZ +1) * voxFaceZ ) j2 = ( v(voxFaceX + 1, voxFaceY, voxFaceZ) * ( 1 - voxFaceZ) ) + (v(voxFaceX + 1 , voxFaceY + 1, voxFaceZ +1) * voxFaceZ ) w1 = (i1 * ( 1 - Yd)) + (i2*Yd) w2 = (j1 * ( 1 - Yd)) + (j2*Yd) d = w1 * ( 1 - Xd) + w2*Xd Here is the result of my implementation: |