Is anyone familiar with the PETACAL subroutine of NASA-VOF2D? In the codes, NF values are claculated based on two variables: PFX and PFY. They are said to be the partial derivatives of F in X and Y directions. But I don't understand how they derived those formulas for surface slopes. If the codes are extended to 3D, how to express PFX, PFY and PFZ? Thank you!

In fact you need to calculate the slope (normal of surface) of interface, then based on principal component of normal select one NF, e.g. if normal vector is: (3., 1, 2) principal direction is positive x-direction and NF is 1 (i don't mention this in nasa code, probably for this nf is 1).

so you must calculate normal: normal is derivative of F (volume fraction) field: in simple way you can sum all 9 neighbour of each side of cell and e.g. for x-component calculate difference of this summation of left and right face over DELT_x. better way is using waighted averaging.

e.g. (in fortran, )

c normal vector: mx,my,mz; c loca(i-1,j-1,k-1)=(i,j,k) in 3d array

mm1 = F(loca(i-1,j-1,k-1))+F(loca(i-1,j-1,k+1))

% +F(loca(i-1,j+1,k-1))

% +F(loca(i-1,j+1,k+1))+2.0d0*(F(loca(i-1,j-1,k))

% +F(loca(i-1,j+1,k))

% +F(loca(i-1,j,k-1))+F(loca(i-1,j,k+1)))

% +4.0d0*F(loca(i-1,j,k))

mm2 = F(loca(i+1,j-1,k-1))+F(loca(i+1,j-1,k+1))

% +F(loca(i+1,j+1,k-1))

% +F(loca(i+1,j+1,k+1))+2.0d0*(F(loca(i+1,j-1,k))

% +F(loca(i+1,j+1,k))

% +F(loca(i+1,j,k-1))+F(loca(i+1,j,k+1)))

% +4.0d0*F(loca(i+1,j,k))

mx = mm1 - mm2

mm1 = F(loca(i-1,j-1,k-1))+F(loca(i-1,j-1,k+1))

% +F(loca(i+1,j-1,k-1))

% +F(loca(i+1,j-1,k+1))+2.0d0*(F(loca(i-1,j-1,k))

% +F(loca(i+1,j-1,k))

% +F(loca(i,j-1,k-1))+F(loca(i,j-1,k+1)))

% +4.0d0*F(loca(i,j-1,k))

mm2 = F(loca(i-1,j+1,k-1))+F(loca(i-1,j+1,k+1))

% +F(loca(i+1,j+1,k-1))

% +F(loca(i+1,j+1,k+1))+2.0d0*(F(loca(i-1,j+1,k))

% +F(loca(i+1,j+1,k))

% +F(loca(i,j+1,k-1))+F(loca(i,j+1,k+1)))

% +4.0d0*F(loca(i,j+1,k))

my = mm1 - mm2

mm1 = F(loca(i-1,j-1,k-1))+F(loca(i-1,j+1,k-1))

% +F(loca(i+1,j-1,k-1))

% +F(loca(i+1,j+1,k-1))+2.0d0*(F(loca(i-1,j,k-1))

% +F(loca(i+1,j,k-1))

% +F(loca(i,j-1,k-1))+F(loca(i,j+1,k-1)))

% +4.0d0*F(loca(i,j,k-1))

mm2 = F(loca(i-1,j-1,k+1))+F(loca(i-1,j+1,k+1))

% +F(loca(i+1,j-1,k+1))

% +F(loca(i+1,j+1,k+1))+2.0d0*(F(loca(i-1,j,k+1))

% +F(loca(i+1,j,k+1))

% +F(loca(i,j-1,k+1))+F(loca(i,j+1,k+1)))

% +4.0d0*F(loca(i,j,k+1))

mz = mm1 - mm2

after calculating (mx, my, mz) compare component of select appropriate direction then NF, if three component is nearly equal you can set NF as 0 or use weighted averaging instead using one interpolating cell.

rt: Thank you very much for your help! I still have some questions with your algorithm. 1. Why are the weights of F(loca(i-1,j-1,k)etc. greater than those of other neighboring cells? I understand the weight of the center cell is greater. But I don't understand why F(loca(i-1,j-1,k)etc are special. 2. What is the meaning of "local(x,y,z)"? Does F(local(x,y,z)) means F(x,y,z)?

Thank you!

consider 9 neighbor of one face (e.g. left face) i divide them to 3 types:

1) first degree weight=4 shared with one face (e.g. i-1,j,k) 2) second degree weight=2 shared with one edge (e.g. i-1,j-1,k) 3) third degree weight=1 shared with one point (e.g. i-1,j-1,k-1)

>>Does F(local(x,y,z)) means F(x,y,z)? yes for better performance i vectorize array (convert 3d array to 1d) and Loca find virtual position (i,j,k) in 1D array.

I see. Thank you very much!