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

Added common curve orientation tensor

Browse Source
This commit is contained in:
James McClure
2014-12-30 09:36:05 -05:00
parent bbbd0ada1e
commit 06b993166d
2 changed files with 55 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

31
common/pmmc.h
View File

@@ -4260,6 +4260,37 @@ inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, Do
}
}
}
inline void pmmc_CurveOrientation(DoubleArray &Orientation, DTMutableList<Point> &Points, int npts, int ic, int jc, int kc){
for (p=0; p<npts-1; p++){
// Extract the line segment
A = Points(p);
B = Points(p+1);
P.x = 0.5*(A.x+B.x) - 1.0*i;
P.y = 0.5*(A.y+B.y) - 1.0*j;
P.z = 0.5*(A.z+B.z) - 1.0*k;
A.x -= 1.0*i;
A.y -= 1.0*j;
A.z -= 1.0*k;
B.x -= 1.0*i;
B.y -= 1.0*j;
B.z -= 1.0*k;
norm = 1.0/((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));// tangent vector
twnsx = norm*(B.x - A.x);
twnsy = norm*(B.y - A.y);
twnsz = norm*(B.z - A.z);
Orientation(0) += 1.0 - twnsx*twnsx; // Gxx
Orientation(1) += 1.0 - twnsy*twnsy; // Gyy
Orientation(2) += 1.0 - twnsz*twnsz; // Gzz
Orientation(3) += 1.0 - twnsx*twnsy; // Gxy
Orientation(4) += 1.0 - twnsx*twnsz; // Gxz
Orientation(5) += 1.0 - twnsy*twnsz; // Gyz
}
}
inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, DoubleArray &KN, DoubleArray &KG,
double &KNavg, double &KGavg, DTMutableList<Point> &Points, int npts, int ic, int jc, int kc){

25
tests/TestContactAngle.cpp
View File

@@ -48,6 +48,7 @@ int main (int argc, char *argv[])
double As;
double efawns,Jwn;
double KNwns,KGwns;
DoubleArray Gwns(6);
// bool add=1; // Set to false if any corners contain nw-phase ( F > fluid_isovalue)
int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners
@@ -132,6 +133,7 @@ int main (int argc, char *argv[])
awn = aws = ans = lwns = 0.0;
nwp_volume = 0.0;
As = 0.0;
for (i=0; i<6; i++) Gwns(i) = 0.0;
for (c=0;c<ncubes;c++){
// Get cube from the list
@@ -167,6 +169,8 @@ int main (int argc, char *argv[])
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues, MeanCurvature, nw_pts, nw_tris,
wn_curvature, i, j, k, n_nw_pts, n_nw_tris);
pmmc_CurveOrientation(Gwns, 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));
@@ -184,6 +188,8 @@ int main (int argc, char *argv[])
KNwns /= lwns;
Jwn /= awn;
efawns /= lwns;
for (i=0; i<6; i++) Gwns(i) /= lwns;
printf("Analysis complete. \n");
double CAPHEIGHT = CAPRAD-sqrt(CAPRAD*CAPRAD-RADIUS*RADIUS); // height of the sphereical cap
@@ -200,6 +206,8 @@ int main (int argc, char *argv[])
printf("Normal curvature (wns) = %f, Analytical = %f \n", KNwns, 1.0/RADIUS);
printf("-------------------------------- \n");
//.........................................................................
printf("Gwns=%i,%i,%i,%i,%i,%i",Gwns(0),Gwns(1),Gwns(2),Gwns(3),Gwns(4),Gwns(5));
int toReturn = 0;
if (fabs(efawns - 1.0*RADIUS/CAPRAD)/(1.0*RADIUS/CAPRAD) > 0.01){
@@ -207,8 +215,23 @@ int main (int argc, char *argv[])
printf("tests/TestContactAngle.cpp: exceeded error tolerance for the contact angle \n");
}
else{
printf("Contact angle: passed test");
printf("Passed test: contact angle");
}
if (fabs(ans - 2*PI*RADIUS*(N-2)-4*PI*RADIUS*(CAPRAD-CAPHEIGHT))/(4*PI*RADIUS*(CAPRAD-CAPHEIGHT)) > 0.01 ){
printf("tests/TestContactAngle.cpp: exceeded error tolerance for ns area \n");
}
else{
printf("Passed test: ans");
}
if (fabs(awn-2*PI*(CAPHEIGHT*CAPHEIGHT+RADIUS*RADIUS))/(2*PI*(CAPHEIGHT*CAPHEIGHT+RADIUS*RADIUS)) > 0.01){
printf("tests/TestContactAngle.cpp: exceeded error tolerance for wn area \n");
}
else{
printf("Passed test: awn");
}
return toReturn;
}