85 lines
3.6 KiB
C++
85 lines
3.6 KiB
C++
//--------------------------------------------------------------------------------------------------------
|
|
inline double ContactAngle(DTPoint3D &pt, DTVectorField3D &gradF, DTVectorField3D &gradS)
|
|
{
|
|
double theta;
|
|
int m = gradS.Grid().m();
|
|
int n = gradS.Grid().n();
|
|
int o = gradS.Grid().o();
|
|
|
|
double dx = gradS.Grid().dx();
|
|
double dy = gradS.Grid().dy();
|
|
double dz = gradS.Grid().dz();
|
|
|
|
DTPoint3D origin = gradS.Grid().Origin();
|
|
|
|
// Get the gradient arrays
|
|
DTFloatArray Fx = gradF.X();
|
|
DTFloatArray Fy = gradF.Y();
|
|
DTFloatArray Fz = gradF.Z();
|
|
|
|
DTFloatArray Sx = gradS.X();
|
|
DTFloatArray Sy = gradS.Y();
|
|
DTFloatArray Sz = gradS.Z();
|
|
|
|
int i,j,k;
|
|
// Determine the cube containing this point
|
|
// i = int(floor((pt.x-origin.x)/dx)); // real space
|
|
// j = int(floor((pt.y-origin.y)/dy));
|
|
// k = int(floor((pt.z-origin.z)/dz));
|
|
i = int(floor(pt.x));
|
|
j = int(floor(pt.y));
|
|
k = int(floor(pt.z));
|
|
|
|
double x,y,z; // map to 0,1 cube
|
|
// x = (pt.x-origin.x-i*dx)/dx; // real space
|
|
// y = (pt.y-origin.y-j*dy)/dy;
|
|
// z = (pt.z-origin.z-k*dz)/dz;
|
|
x = pt.x - double(i);
|
|
y = pt.y - double(j);
|
|
z = pt.z - double(k);
|
|
|
|
// evaluate contact angle at corners of the cube
|
|
DTMutableDoubleArray cube(2,2,2);
|
|
// theta = acos ( -(gradF*gradS) / (|gradF| |gradS|) )
|
|
cube(0,0,0) = -( Fx(i,j,k)*Sx(i,j,k)+Fy(i,j,k)*Sy(i,j,k)+Fz(i,j,k)*Sz(i,j,k) )
|
|
/( sqrt(pow(Fx(i,j,k),2)+pow(Fy(i,j,k),2)+pow(Fz(i,j,k),2))
|
|
*sqrt(pow(Sx(i,j,k),2)+pow(Sy(i,j,k),2)+pow(Sz(i,j,k),2)) );
|
|
cube(1,0,0) = -( Fx(i+1,j,k)*Sx(i+1,j,k)+Fy(i+1,j,k)*Sy(i+1,j,k)+Fz(i+1,j,k)*Sz(i+1,j,k) )
|
|
/( sqrt(pow(Fx(i+1,j,k),2)+pow(Fy(i+1,j,k),2)+pow(Fz(i+1,j,k),2))
|
|
*sqrt(pow(Sx(i+1,j,k),2)+pow(Sy(i+1,j,k),2)+pow(Sz(i+1,j,k),2)) );
|
|
cube(0,1,0) = -( Fx(i,j+1,k)*Sx(i,j+1,k)+Fy(i,j+1,k)*Sy(i,j+1,k)+Fz(i,j+1,k)*Sz(i,j+1,k) )
|
|
/( sqrt(pow(Fx(i,j+1,k),2)+pow(Fy(i,j+1,k),2)+pow(Fz(i,j+1,k),2))
|
|
*sqrt(pow(Sx(i,j+1,k),2)+pow(Sy(i,j+1,k),2)+pow(Sz(i,j+1,k),2)) );
|
|
cube(0,0,1) = -( Fx(i,j,k+1)*Sx(i,j,k+1)+Fy(i,j,k+1)*Sy(i,j,k+1)+Fz(i,j,k+1)*Sz(i,j,k+1) )
|
|
/( sqrt(pow(Fx(i,j,k+1),2)+pow(Fy(i,j,k+1),2)+pow(Fz(i,j,k+1),2))
|
|
*sqrt(pow(Sx(i,j,k+1),2)+pow(Sy(i,j,k+1),2)+pow(Sz(i,j,k+1),2)) );
|
|
cube(1,1,0) = -( Fx(i+1,j+1,k)*Sx(i+1,j+1,k)+Fy(i+1,j+1,k)*Sy(i+1,j+1,k)+Fz(i+1,j+1,k)*Sz(i+1,j+1,k) )
|
|
/( sqrt(pow(Fx(i+1,j+1,k),2)+pow(Fy(i+1,j+1,k),2)+pow(Fz(i+1,j+1,k),2))
|
|
*sqrt(pow(Sx(i+1,j+1,k),2)+pow(Sy(i+1,j+1,k),2)+pow(Sz(i+1,j+1,k),2)) );
|
|
cube(1,0,1) = -( Fx(i+1,j,k+1)*Sx(i+1,j,k+1)+Fy(i+1,j,k+1)*Sy(i+1,j,k+1)+Fz(i+1,j,k+1)*Sz(i+1,j,k+1) )
|
|
/( sqrt(pow(Fx(i+1,j,k+1),2)+pow(Fy(i+1,j,k+1),2)+pow(Fz(i+1,j,k+1),2))
|
|
*sqrt(pow(Sx(i+1,j,k+1),2)+pow(Sy(i+1,j,k+1),2)+pow(Sz(i+1,j,k+1),2)) );
|
|
cube(0,1,1) = -( Fx(i,j+1,k+1)*Sx(i,j+1,k+1)+Fy(i,j+1,k+1)*Sy(i,j+1,k+1)+Fz(i,j+1,k+1)*Sz(i,j+1,k+1) )
|
|
/( sqrt(pow(Fx(i,j+1,k+1),2)+pow(Fy(i,j+1,k+1),2)+pow(Fz(i,j+1,k+1),2))
|
|
*sqrt(pow(Sx(i,j+1,k+1),2)+pow(Sy(i,j+1,k+1),2)+pow(Sz(i,j+1,k+1),2)) );
|
|
cube(1,1,1) = -( Fx(i+1,j+1,k+1)*Sx(i+1,j+1,k+1)+Fy(i+1,j+1,k+1)*Sy(i+1,j+1,k+1)+Fz(i+1,j+1,k+1)*Sz(i+1,j+1,k+1) )
|
|
/( sqrt(pow(Fx(i+1,j+1,k+1),2)+pow(Fy(i+1,j+1,k+1),2)+pow(Fz(i+1,j+1,k+1),2))
|
|
*sqrt(pow(Sx(i+1,j+1,k+1),2)+pow(Sy(i+1,j+1,k+1),2)+pow(Sz(i+1,j+1,k+1),2)) );
|
|
|
|
// construct polynomial approximation for contact angle
|
|
double a,b,c,d,e,f,g,h; // trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
|
|
a = cube(0,0,0);
|
|
b = cube(1,0,0)-a;
|
|
c = cube(0,1,0)-a;
|
|
d = cube(0,0,1)-a;
|
|
e = cube(1,1,0)-a-b-c;
|
|
f = cube(1,0,1)-a-b-d;
|
|
g = cube(0,1,1)-a-c-d;
|
|
h = cube(1,1,1)-a-b-c-d-e-f-g;
|
|
|
|
// evaluate at x,y,z
|
|
theta = acos (a+b*x+c*y+d*z+e*x*y+f*x*z+g*y*z+h*x*y*z );
|
|
|
|
return theta;
|
|
}
|