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

Euler characteristic refactored and working for TestBlobAnalyze for a sphere. Using the entire NWP surface rather than trying to stitch the wn+ns surfaces together

Browse Source
This commit is contained in:
James E McClure
2015-08-30 14:25:00 -04:00
parent 4a8399e835
commit 355f093022
1 changed files with 2 additions and 96 deletions
Download Patch File Download Diff File Expand all files Collapse all files

98
common/TwoPhase.cpp
View File

@@ -584,99 +584,6 @@ void TwoPhase::ComponentAverages()
n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
i, j, k, Nx, Ny, Nz);
if (n_nw_pts+n_ns_pts > 0){
/* double euler;
euler = geomavg_EulerCharacteristic(nw_pts,nw_tris,n_nw_pts,n_nw_tris,i,j,k);
euler += geomavg_EulerCharacteristic(ns_pts,ns_tris,n_ns_pts,n_ns_tris,i,j,k);
// adjust for double-counted vertices and edges from the common curve
if (n_nws_pts > 0) euler += 1.0*n_nws_pts;
int nvert = n_nw_pts+n_ns_pts-n_nws_pts;
int nside = 2*nvert-3;
int nface = nvert-2;
//...........................................................
// Check that this point is not on a previously computed face
// Note direction that the marching cubes algorithm marches
// In parallel, other sub-domains fill in the lower boundary
for (int p=0; p<n_nw_pts; p++){
Point PT = nw_pts(p);
if (PT.x - double(i) < 1e-12) nvert-=1;
else if (PT.y - double(j) < 1e-12) nvert-=1;
else if (PT.z - double(k) < 1e-12) nvert-=1;
}
for (int p=0; p<n_ns_pts; p++){
Point PT = ns_pts(p);
if (PT.x - double(i) < 1e-12) nvert-=1;
else if (PT.y - double(j) < 1e-12) nvert-=1;
else if (PT.z - double(k) < 1e-12) nvert-=1;
}
// Remove previously computed edges
for (int p=0; p<n_nw_tris; p++){
Point A = nw_pts(nw_tris(0,p));
Point B = nw_pts(nw_tris(1,p));
Point C = nw_pts(nw_tris(2,p));
// Check side A-B
bool newside = true;
if (A.x - double(i) < 1e-12 && B.x - double(i) < 1e-12) newside=false;
if (A.y - double(j) < 1e-12 && B.y - double(j) < 1e-12) newside=false;
if (A.z - double(k) < 1e-12 && B.z - double(k) < 1e-12) newside=false;
if (!newside) nside-=1;
// Check side A-C
newside = true;
if (A.x - double(i)< 1e-12 && C.x - double(i) < 1e-12) newside=false;
if (A.y - double(j)< 1e-12 && C.y - double(j) < 1e-12) newside=false;
if (A.z - double(k)< 1e-12 && C.z - double(k) < 1e-12) newside=false;
if (!newside) nside-=1;
// Check side B-C
newside = true;
if (B.x - double(i) < 1e-12 && C.x - double(i) < 1e-12) newside=false;
if (B.y - double(j) < 1e-12 && C.y - double(j) < 1e-12) newside=false;
if (B.z - double(k) < 1e-12 && C.z - double(k) < 1e-12) newside=false;
if (!newside) nside-=1;
}
for (int p=0; p<n_ns_tris; p++){
Point A = ns_pts(ns_tris(0,p));
Point B = ns_pts(ns_tris(1,p));
Point C = ns_pts(ns_tris(2,p));
// Check side A-B
bool newside = true;
if (A.x - double(i) < 1e-12 && B.x - double(i) < 1e-12) newside=false;
if (A.y - double(j) < 1e-12 && B.y - double(j) < 1e-12) newside=false;
if (A.z - double(k) < 1e-12 && B.z - double(k) < 1e-12) newside=false;
if (!newside) nside-=1;
// Check side A-C
newside = true;
if (A.x - double(i) < 1e-12 && C.x - double(i) < 1e-12) newside=false;
if (A.y - double(j) < 1e-12 && C.y - double(j) < 1e-12) newside=false;
if (A.z - double(k) < 1e-12 && C.z - double(k) < 1e-12) newside=false;
if (!newside) nside-=1;
// Check side B-C
newside = true;
if (B.x - double(i) < 1e-12 && C.x - double(i) < 1e-12) newside=false;
if (B.y - double(j) < 1e-12 && C.y - double(j) < 1e-12) newside=false;
if (B.z - double(k) < 1e-12 && C.z - double(k) < 1e-12) newside=false;
if (!newside) nside-=1;
}
int euler=nvert-nside+nface; // euler characteristic for the cube
ComponentAverages_NWP(EULER,LabelNWP) += 1.0*euler;
// Counts for the vertices, sides and faces from the local cube
//ComponentAverages_NWP(NVERT,LabelNWP) += nvert;
//ComponentAverages_NWP(NSIDE,LabelNWP) += nside;
//ComponentAverages_NWP(NFACE,LabelNWP) += nface;
*
*/
}
//...........................................................................
// wn interface averages
if (n_nw_pts>0 && LabelNWP >=0 && LabelWP >=0 ){
@@ -782,11 +689,10 @@ void TwoPhase::ComponentAverages()
* double geomavg_EulerCharacteristic(PointList, PointCount, TriList, TriCount);
*/
if (n_nw_pts + n_ns_pts > 0 ){
printf("n_nw_pts=%i \n",n_nw_pts);
double euler;
n_nw_pts=n_nw_tris=0;
geomavg_MarchingCubes(SDn,fluid_isovalue,i,j,k,nw_pts,n_nw_pts,nw_tris,n_nw_tris);
printf("n_nw_pts=%i, n_nw_tris=%i \n",n_nw_pts,n_nw_tris);
double euler = geomavg_EulerCharacteristic(nw_pts,nw_tris,n_nw_pts,n_nw_tris,i,j,k);
euler = geomavg_EulerCharacteristic(nw_pts,nw_tris,n_nw_pts,n_nw_tris,i,j,k);
ComponentAverages_NWP(EULER,LabelNWP) += euler;
}
}