OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Updated common/pmmc.h and associated unit tests

Browse Source
This commit is contained in:
James E McClure 2014-10-23 16:58:38 -04:00
parent eb254d9b9c
commit 4db486ce7c
3 changed files with 73 additions and 227 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

283
common/pmmc.h
View File

@ -4153,7 +4153,7 @@ inline double pmmc_CubeSurfaceOrientation(DoubleArray &Orientation, DTMutableLis
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, DoubleArray &ReturnVector, inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, DoubleArray &ReturnVector,
DoubleArray &Fx, DoubleArray &Fy, DoubleArray &Fz, DoubleArray &Fx, DoubleArray &Fy, DoubleArray &Fz,
DoubleArray &Sx, DoubleArray &Sy, DoubleArray &Sz, DoubleArray &Sx, DoubleArray &Sy, DoubleArray &Sz,
DTMutableList<Point> &Points, int i, int j, int k, int npts) DTMutableList<Point> &Points, int i, int j, int k, int npts)
{ {
int p; int p;
@ -4164,43 +4164,37 @@ inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, Do
double ns_x, ns_y, ns_z; double ns_x, ns_y, ns_z;
double nwn_x, nwn_y, nwn_z; double nwn_x, nwn_y, nwn_z;
double nwns_x, nwns_y, nwns_z; double nwns_x, nwns_y, nwns_z;
Point A,B; Point P,A,B;
lwns = 0.0; lwns = 0.0;
//.............................................................................
// Compute the normal and the speed at points on the interface
CubeValues(0,0,0) = dPdt(i,j,k);
CubeValues(1,0,0) = dPdt(i+1,j,k);
CubeValues(0,1,0) = dPdt(i,j+1,k);
CubeValues(1,1,0) = dPdt(i+1,j+1,k);
CubeValues(0,0,1) = dPdt(i,j,k+1);
CubeValues(1,0,1) = dPdt(i+1,j,k+1);
CubeValues(0,1,1) = dPdt(i,j+1,k+1);
CubeValues(1,1,1) = dPdt(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz TriLinPoly Px,Py,Pz,SDx,SDy,SDz,Pt;
a = CubeValues(0,0,0); Px.assign(Fx,i,j,k);
b = CubeValues(1,0,0)-a; Py.assign(Fy,i,j,k);
c = CubeValues(0,1,0)-a; Pz.assign(Fz,i,j,k);
d = CubeValues(0,0,1)-a; SDx.assign(Sx,i,j,k);
e = CubeValues(1,1,0)-a-b-c; SDy.assign(Sy,i,j,k);
f = CubeValues(1,0,1)-a-b-d; SDz.assign(Sz,i,j,k);
g = CubeValues(0,1,1)-a-c-d; Pt.assign(dPdt,i,j,k);
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p < npts-1; p++){ for (p=0; p < npts-1; p++){
// Extract the line segment // Extract the line segment
A = Points(p); A = Points(p);
B = Points(p+1); B = Points(p+1);
// Midpoint of the line
P.x = 0.5*(A.x+B.x);
P.y = 0.5*(A.y+B.y);
P.z = 0.5*(A.z+B.z);
// Compute the curve tangent // Compute the curve tangent
tangent_x = A.x - B.x; tangent_x = A.x - B.x;
tangent_y = A.y - B.y; tangent_y = A.y - B.y;
tangent_z = A.z - B.z; tangent_z = A.z - B.z;
// Get the normal to the Color gradient // Get the normal to the Color gradient
nwn_x = Fx(i,j,k); nwn_x = Px.eval(P);
nwn_y = Fy(i,j,k); nwn_y = Py.eval(P);
nwn_z = Fz(i,j,k); nwn_z = Pz.eval(P);
norm = nwn_x*nwn_x + nwn_y*nwn_y + nwn_z*nwn_z; norm = nwn_x*nwn_x + nwn_y*nwn_y + nwn_z*nwn_z;
// Compute the interface speed
zeta = -Pt.eval(P) / norm;
if (norm > 0.0){ if (norm > 0.0){
nwn_x /= norm; nwn_x /= norm;
nwn_y /= norm; nwn_y /= norm;
@ -4208,9 +4202,9 @@ inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, Do
} }
// Get the normal to the solid surface // Get the normal to the solid surface
ns_x = Sx(i,j,k); ns_x = SDx.eval(P);
ns_y = Sy(i,j,k); ns_y = SDy.eval(P);;
ns_z = Sz(i,j,k); ns_z = SDz.eval(P);
norm = ns_x*ns_x + ns_y*ns_y + ns_z*ns_z; norm = ns_x*ns_x + ns_y*ns_y + ns_z*ns_z;
if (norm > 0.0){ if (norm > 0.0){
ns_x /= norm; ns_x /= norm;
@ -4235,24 +4229,16 @@ inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, Do
nwns_y /= norm; nwns_y /= norm;
nwns_z /= norm; nwns_z /= norm;
} }
// Compute the interface speed at the midpoint of the segment
x = 0.5*(A.x+B.x);
y = 0.5*(A.y+B.y);
z = 0.5*(A.z+B.z);
zeta = a + b*x + c*y + d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
//norm = ns_x*ns_x + ns_y*ns_y + ns_z*ns_z;
norm = nwn_x*nwn_x + nwn_y*nwn_y + nwn_z*nwn_z;
// Compute the common curve velocity
if (norm > 0.0){
ReturnVector(0) += -(zeta/norm)*nwns_x;
ReturnVector(1) += -(zeta/norm)*nwns_y;
ReturnVector(2) += -(zeta/norm)*nwns_z;
}
// Compute the length of the segment // Compute the length of the segment
s = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z)); s = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
// Add the length to the common line // Add the length to the common line
lwns += s; // lwns += s;
// Compute the common curve velocity
if (norm > 0.0){
ReturnVector(0) += zeta*(nwn_x*nwns_x+nwn_y*nwns_y+nwn_z*nwns_z)*nwns_x*s;
ReturnVector(1) += zeta*(nwn_x*nwns_x+nwn_y*nwns_y+nwn_z*nwns_z)*nwns_y*s;
ReturnVector(2) += zeta*(nwn_x*nwns_x+nwn_y*nwns_y+nwn_z*nwns_z)*nwns_z*s;
}
} }
} }
@ -4310,9 +4296,9 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, DoubleArray &KN,
Ayx = syy*fz + sy*fyz - syz*fy - sz*fyy; Ayx = syy*fz + sy*fyz - syz*fy - sz*fyy;
Ayy = syz*fx + sz*fxy - sxy*fz - sx*fyz; Ayy = syz*fx + sz*fxy - sxy*fz - sx*fyz;
Ayz = sxy*fy + sx*fyy - syy*fx - sy*fxy; Ayz = sxy*fy + sx*fyy - syy*fx - sy*fxy;
Axx = syz*fz + sy*fzz - szz*fy - sz*fyz; Azx = syz*fz + sy*fzz - szz*fy - sz*fyz;
Axy = szz*fx + sz*fxz - sxz*fz - sx*fzz; Azy = szz*fx + sz*fxz - sxz*fz - sx*fzz;
Axz = sxz*fy + sx*fyz - syz*fx - sy*fxz; Azz = sxz*fy + sx*fyz - syz*fx - sy*fxz;
// Normal to solid surface // Normal to solid surface
Sx.Corners(i-ic,j-jc,k-kc) = sx; Sx.Corners(i-ic,j-jc,k-kc) = sx;
@ -4352,13 +4338,14 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, DoubleArray &KN,
// tangent vector // tangent vector
twnsx = Tx.eval(P); twnsx = Tx.eval(P);
twnsy = Ty.eval(P); twnsy = Ty.eval(P);
twnsx = Tz.eval(P); twnsz = Tz.eval(P);
// normal vector and curvature to the wns curve // normal vector and curvature to the wns curve
nwnsx = Nx.eval(P); nwnsx = Nx.eval(P);
nwnsy = Ny.eval(P); nwnsy = Ny.eval(P);
nwnsx = Nz.eval(P); nwnsx = Nz.eval(P);
K = sqrt(nwnsx*nwnsx+nwnsy*nwnsy+nwnsz*nwnsz); K = sqrt(nwnsx*nwnsx+nwnsy*nwnsy+nwnsz*nwnsz);
if (K == 0.0) K=1.0;
nwnsx /= K; nwnsx /= K;
nwnsy /= K; nwnsy /= K;
nwnsz /= K; nwnsz /= K;
@ -4368,6 +4355,7 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, DoubleArray &KN,
nsy = Sy.eval(P); nsy = Sy.eval(P);
nsz = Sz.eval(P); nsz = Sz.eval(P);
norm = sqrt(nsx*nsx+nsy*nsy+nsz*nsz); norm = sqrt(nsx*nsx+nsy*nsy+nsz*nsz);
if (norm == 0.0) norm=1.0;
nsx /= norm; nsx /= norm;
nsy /= norm; nsy /= norm;
nsz /= norm; nsz /= norm;
@ -4377,6 +4365,7 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, DoubleArray &KN,
nwsy = twnsz*nsx-twnsx*nsz; nwsy = twnsz*nsx-twnsx*nsz;
nwsz = twnsx*nsy-twnsy*nsx; nwsz = twnsx*nsy-twnsy*nsx;
norm = sqrt(nwsx*nwsx+nwsy*nwsy+nwsz*nwsz); norm = sqrt(nwsx*nwsx+nwsy*nwsy+nwsz*nwsz);
if (norm == 0.0) norm=1.0;
nwsx /= norm; nwsx /= norm;
nwsy /= norm; nwsy /= norm;
nwsz /= norm; nwsz /= norm;
@ -4392,8 +4381,11 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, DoubleArray &KN,
A = Points(p); A = Points(p);
B = Points(p+1); B = Points(p+1);
length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z)); length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
KNavg += 0.5*length*(KN(p)+KN(p+1)); if (length > 0.0){
KGavg += 0.5*length*(KG(p)+KG(p+1)); KNavg += 0.5*length*(KN(p)+KN(p+1));
KGavg += 0.5*length*(KG(p)+KG(p+1));
// KGavg += length;
}
} }
} }
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
@ -4402,149 +4394,21 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
DoubleArray &SurfaceVector, DoubleArray &SurfaceValues, DoubleArray &AvgVel, DoubleArray &SurfaceVector, DoubleArray &SurfaceValues, DoubleArray &AvgVel,
int i, int j, int k, int npts, int ntris) int i, int j, int k, int npts, int ntris)
{ {
Point A,B,C; Point A,B,C,P;
int p; int p;
double vA,vB,vC; double vA,vB,vC;
double vAx,vBx,vCx,vAy,vBy,vCy,vAz,vBz,vCz; double vAx,vBx,vCx,vAy,vBy,vCy,vAz,vBz,vCz;
double x,y,z; double x,y,z;
double s,s1,s2,s3,temp; double s,s1,s2,s3,temp;
double a,b,c,d,e,f,g,h; double px,py,pz;
double norm, zeta; double norm, zeta;
// ................x component ............................. TriLinPoly Px,Py,Pz,Pt;
// Copy the x derivative values for the cube Px.assign(P_x,i,j,k);
CubeValues(0,0,0) = P_x(i,j,k); Py.assign(P_y,i,j,k);
CubeValues(1,0,0) = P_x(i+1,j,k); Pz.assign(P_z,i,j,k);
CubeValues(0,1,0) = P_x(i,j+1,k); Pt.assign(dPdt,i,j,k);
CubeValues(1,1,0) = P_x(i+1,j+1,k);
CubeValues(0,0,1) = P_x(i,j,k+1);
CubeValues(1,0,1) = P_x(i+1,j,k+1);
CubeValues(0,1,1) = P_x(i,j+1,k+1);
CubeValues(1,1,1) = P_x(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
SurfaceVector(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
// ................y component .............................
// Copy the y derivative values for the cube
CubeValues(0,0,0) = P_y(i,j,k);
CubeValues(1,0,0) = P_y(i+1,j,k);
CubeValues(0,1,0) = P_y(i,j+1,k);
CubeValues(1,1,0) = P_y(i+1,j+1,k);
CubeValues(0,0,1) = P_y(i,j,k+1);
CubeValues(1,0,1) = P_y(i+1,j,k+1);
CubeValues(0,1,1) = P_y(i,j+1,k+1);
CubeValues(1,1,1) = P_y(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
SurfaceVector(npts+p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
// ................z component .............................
// Copy the z derivative values for the cube
CubeValues(0,0,0) = P_z(i,j,k);
CubeValues(1,0,0) = P_z(i+1,j,k);
CubeValues(0,1,0) = P_z(i,j+1,k);
CubeValues(1,1,0) = P_z(i+1,j+1,k);
CubeValues(0,0,1) = P_z(i,j,k+1);
CubeValues(1,0,1) = P_z(i+1,j,k+1);
CubeValues(0,1,1) = P_z(i,j+1,k+1);
CubeValues(1,1,1) = P_z(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
SurfaceVector(2*npts+p) = a + b*x + c*y + d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
//.............................................................................
// Compute the normal and the speed at points on the interface
// Copy the curvature values for the cube
CubeValues(0,0,0) = dPdt(i,j,k);
CubeValues(1,0,0) = dPdt(i+1,j,k);
CubeValues(0,1,0) = dPdt(i,j+1,k);
CubeValues(1,1,0) = dPdt(i+1,j+1,k);
CubeValues(0,0,1) = dPdt(i,j,k+1);
CubeValues(1,0,1) = dPdt(i+1,j,k+1);
CubeValues(0,1,1) = dPdt(i,j+1,k+1);
CubeValues(1,1,1) = dPdt(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
// evaluate time derivative on the surface
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
zeta = a + b*x + c*y + d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
// compute the normal
x = SurfaceVector(p);
y = SurfaceVector(npts+p);
z = SurfaceVector(2*npts+p);
norm = sqrt(x*x+y*y+z*z);
// Save the surface values and normal vector
if (norm > 0.0){
SurfaceValues(p) = -zeta/norm;
SurfaceVector(p) = x/norm;
SurfaceVector(npts+p) = y/norm;
SurfaceVector(2*npts+p) = z/norm;
}
else{
SurfaceValues(p) = 0.0;
SurfaceVector(p) = 0.0;
SurfaceVector(npts+p) = 0.0;
SurfaceVector(2*npts+p) = 0.0;
}
}
//............................................................................. //.............................................................................
// Compute the average speed of the interface // Compute the average speed of the interface
for (int r=0; r<ntris; r++){ for (int r=0; r<ntris; r++){
@ -4557,44 +4421,21 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z)); s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3); s = 0.5*(s1+s2+s3);
temp = s*(s-s1)*(s-s2)*(s-s3); temp = s*(s-s1)*(s-s2)*(s-s3);
// Compute the centroid P
P.x = 0.33333333333333333*(A.x+B.x+C.x);
P.y = 0.33333333333333333*(A.y+B.y+C.y);
P.z = 0.33333333333333333*(A.z+B.z+C.z);
if (temp > 0.0){ if (temp > 0.0){
// Surface value (speed)
vA = SurfaceValues(Triangles(0,r));
vB = SurfaceValues(Triangles(1,r));
vC = SurfaceValues(Triangles(2,r));
// Increment the averaged values
// x component
vAx = SurfaceVector(Triangles(0,r))*vA;
vBx = SurfaceVector(Triangles(1,r))*vB;
vCx = SurfaceVector(Triangles(2,r))*vC;
// y component
vAy = SurfaceVector(npts+Triangles(0,r))*vA;
vBy = SurfaceVector(npts+Triangles(1,r))*vB;
vCy = SurfaceVector(npts+Triangles(2,r))*vC;
// z component
vAz = SurfaceVector(2*npts+Triangles(0,r))*vA;
vBz = SurfaceVector(2*npts+Triangles(1,r))*vB;
vCz = SurfaceVector(2*npts+Triangles(2,r))*vC;
AvgVel(0) += sqrt(temp)*0.33333333333333333*(vAx+vBx+vCx); x = Px.eval(P);
AvgVel(1) += sqrt(temp)*0.33333333333333333*(vAy+vBy+vCy); y = Py.eval(P);
AvgVel(2) += sqrt(temp)*0.33333333333333333*(vAz+vBz+vCz); z = Pz.eval(P);
norm = sqrt(x*x+y*y+z*z);
// Update the Averages. Differentiate between advancing (0,1,2) and receding (3,4,5) interfaces zeta = -Pt.eval(P) / norm;
// All points on a triangle have the same orientation in the color gradient
/* if (vA > 0.0){ AvgVel(0) += sqrt(temp)*zeta*x/norm;
// Advancing interface AvgVel(1) += sqrt(temp)*zeta*y/norm;
AvgVel(0) += sqrt(temp)*0.33333333333333333*(vAx+vBx+vCx); AvgVel(2) += sqrt(temp)*zeta*z/norm;
AvgVel(1) += sqrt(temp)*0.33333333333333333*(vAy+vBy+vCy);
AvgVel(2) += sqrt(temp)*0.33333333333333333*(vAz+vBz+vCz);
}
else{
// Receding interface
AvgVel(3) += sqrt(temp)*0.33333333333333333*(vAx+vBx+vCx);
AvgVel(4) += sqrt(temp)*0.33333333333333333*(vAy+vBy+vCy);
AvgVel(5) += sqrt(temp)*0.33333333333333333*(vAz+vBz+vCz);
}
*/
} }
} }
//............................................................................. //.............................................................................

4
tests/TestCylinderAreas.cpp
View File

@ -274,7 +274,7 @@ int main (int argc, char *argv[])
toReturn += 1; toReturn += 1;
printf("TestCylinderArea.cpp: error tolerance exceeded for wn area \n"); printf("TestCylinderArea.cpp: error tolerance exceeded for wn area \n");
} }
if (fabs(ans - (2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT)))/(2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT)) > 0.02 ){ if (fabs(ans - (2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT))/(2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT)> 0.02 ){
toReturn += 2; toReturn += 2;
printf("TestCylinderArea.cpp: error tolerance exceeded for ns area \n"); printf("TestCylinderArea.cpp: error tolerance exceeded for ns area \n");
} }
@ -290,6 +290,6 @@ int main (int argc, char *argv[])
toReturn += 5; toReturn += 5;
printf("TestCylinderArea.cpp: error tolerance exceeded for common curve length \n"); printf("TestCylinderArea.cpp: error tolerance exceeded for common curve length \n");
} }
return toReturn; return toReturn;
} }

13
tests/TestInterfaceSpeed.cpp
View File

@ -94,6 +94,7 @@ int main (int argc, char *argv[])
DoubleArray InterfaceSpeed(20); DoubleArray InterfaceSpeed(20);
DoubleArray NormalVector(60); DoubleArray NormalVector(60);
DoubleArray vawn(6); DoubleArray vawn(6);
DoubleArray vawns(3);
int c; int c;
//........................................................................... //...........................................................................
@ -202,6 +203,9 @@ int main (int argc, char *argv[])
pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris, pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris,
NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris); NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, Phase_x,Phase_y,Phase_z,Sx,Sy,Sz,
local_nws_pts,i,j,k,n_local_nws_pts);
pmmc_CurveCurvature(Phase, SignDist, KNwns_values, KGwns_values, KNwns, KGwns, nws_pts, n_nws_pts, i, j, k); pmmc_CurveCurvature(Phase, SignDist, KNwns_values, KGwns_values, KNwns, KGwns, nws_pts, n_nws_pts, i, j, k);
// if (n_nw_pts>0) printf("speed %f \n",InterfaceSpeed(0)); // if (n_nw_pts>0) printf("speed %f \n",InterfaceSpeed(0));
@ -219,6 +223,7 @@ int main (int argc, char *argv[])
Jwn /= awn; Jwn /= awn;
efawns /= lwns; efawns /= lwns;
for (i=0;i<6;i++) vawn(i) /= awn; for (i=0;i<6;i++) vawn(i) /= awn;
for (i=0;i<3;i++) vawns(i) /= lwns;
printf("-------------------------------- \n"); printf("-------------------------------- \n");
printf("NWP volume = %f \n", nwp_volume); printf("NWP volume = %f \n", nwp_volume);
@ -226,12 +231,12 @@ int main (int argc, char *argv[])
printf("Area ns = %f, Analytical = %f \n", ans, 2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT); printf("Area ns = %f, Analytical = %f \n", ans, 2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT);
printf("Area ws = %f, Analytical = %f \n", aws, 4*PI*RADIUS*HEIGHT); printf("Area ws = %f, Analytical = %f \n", aws, 4*PI*RADIUS*HEIGHT);
printf("Area s = %f, Analytical = %f \n", As, 2*PI*RADIUS*(N-2)); printf("Area s = %f, Analytical = %f \n", As, 2*PI*RADIUS*(N-2));
printf("Geodesic curvature (wns) = %f, Analytical = %f \n", KGwns, RADIUS); printf("Geodesic curvature (wns) = %f, Analytical = %f \n", KGwns, 0.0);
printf("Normal curvature (wns) = %f, Analytical = %f \n", KNwns, 0); printf("Normal curvature (wns) = %f, Analytical = %f \n", KNwns, 1.0/RADIUS);
printf("Length wns = %f, Analytical = %f \n", lwns, 4*PI*RADIUS); printf("Length wns = %f, Analytical = %f \n", lwns, 4*PI*RADIUS);
// printf("Cos(theta_wns) = %f, Analytical = %f \n",efawns/lwns,1.0*RADIUS/CAPRAD); // printf("Cos(theta_wns) = %f, Analytical = %f \n",efawns/lwns,1.0*RADIUS/CAPRAD);
printf("Advancing Interface Velocity = %f,%f,%f \n",vawn(0),vawn(1),vawn(2)); printf("Interface Velocity = %f,%f,%f \n",vawn(0),vawn(1),vawn(2));
printf("Receding Interface Velocity = %f,%f,%f \n",vawn(3),vawn(4),vawn(5)); printf("Common Curve Velocity = %f,%f,%f \n",vawns(0),vawns(1),vawns(2));
printf("-------------------------------- \n"); printf("-------------------------------- \n");
//......................................................................... //.........................................................................