From 80e3aba9a64782ec8e1d8cbb1530fe812ef2a900 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sun, 22 Nov 2015 14:39:19 -0500 Subject: [PATCH] Refactoring Euler characteristic --- common/TwoPhase.cpp | 33 +++++--------------------------- common/pmmc.h | 46 +++++++++++++++++++++++++++++++++++++++++++++ tests/TestTorus.cpp | 2 +- 3 files changed, 52 insertions(+), 29 deletions(-) diff --git a/common/TwoPhase.cpp b/common/TwoPhase.cpp index 7c9e6839..4263cd9f 100644 --- a/common/TwoPhase.cpp +++ b/common/TwoPhase.cpp @@ -437,33 +437,7 @@ void TwoPhase::ComputeLocal() } } } -/* // Compute local contrubition to Euler characteristic based on 6 adjacency - for (int p=0;p<8;p++){ - // binary id for the wetting phase - double binid= PhaseID(i+cube[p][0],j+cube[p][1],k+cube[p][2]; - if (binid != 1) binid=0; // non-wetting phase - CubeValues(cube[p][0],cube[p][1],cube[p][2])=binid; - } - // update faces edges and cubes for NWP - epc_ncubes += CubeValues(0,0,0)*CubeValues(1,0,0)*CubeValues(0,1,0)*CubeValues(0,0,1)* - CubeValues(1,1,0)*CubeValues(1,0,1)*CubeValues(0,1,1)*CubeValues(1,1,1); - // three faces (others shared by other cubes) - epc_nface += CubeValues(1,0,0)*CubeValues(1,1,0)*CubeValues(1,0,1)*CubeValues(1,1,1); - epc_nface += CubeValues(0,1,0)*CubeValues(1,1,0)*CubeValues(0,1,1)*CubeValues(1,1,1); - epc_nface += CubeValues(0,0,1)*CubeValues(0,1,1)*CubeValues(1,0,1)*CubeValues(1,1,1); - // six of twelve edges (others shared by other cubes) - epc_nedge += CubeValues(1,1,1)*CubeValues(1,1,0); - epc_nedge += CubeValues(1,1,1)*CubeValues(1,0,1); - epc_nedge += CubeValues(1,1,1)*CubeValues(0,1,1); - epc_nedge += CubeValues(1,0,1)*CubeValues(1,0,0); - epc_nedge += CubeValues(1,0,1)*CubeValues(0,0,1); - epc_nedge += CubeValues(0,1,1)*CubeValues(0,0,1); - // four of eight vertices - epc_nvert += CubeValues(1,1,0); - epc_nvert += CubeValues(1,0,1); - epc_nvert += CubeValues(0,1,1); - epc_nvert += CubeValues(1,1,1); -*/ + //........................................................................... // Construct the interfaces and common curve pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue, @@ -535,7 +509,10 @@ void TwoPhase::ComputeLocal() Kn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values, i,j,k,n_nw_pts,n_nw_tris); - 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); + + // Compute the euler characteristic from phase count + euler += mink_phase_epc6(PhaseID,CubeValues,1,i,j,k); } } diff --git a/common/pmmc.h b/common/pmmc.h index 291da3df..61258241 100644 --- a/common/pmmc.h +++ b/common/pmmc.h @@ -4427,6 +4427,52 @@ inline double geomavg_EulerCharacteristic(DTMutableList &Points, IntArray return EulerChar; } +inline double mink_phase_epc6(IntArray &PhaseID, DoubleArray &CubeValues, int PhaseLabel, + int &i, int &j, int &k){ + + /* + * Compute the Euler Poincare characteristic for a phase based on 6 adjacency + * compute the local contribution within the specified cube must sum to get total + * PhaseID -- label of phases on a mesh + * PhaseLabel -- label assigned to phase that we are computing + * i,j,k -- identifies the local cube + */ + + double epc_ncubes, epc_nface, epc_nedge, epc_nvert; + epc_ncubes=epc_nface=epc_nedge=epc_nvert = 0; + + // Assign CubeValues to be phase indicator for + for (int kk=0; kk<2; kk++){ + for (int jj=0; jj<2; jj++){ + for (int ii=0; ii<2; ii++){ + if ( PhaseID(i+ii,j+jj,k+kk) == PhaseLabel) CubeValues(ii,jj,kk)=1; + else CubeValues(ii,jj,kk)=0; + } + } + } + // update faces edges and cubes for NWP + epc_ncubes += CubeValues(0,0,0)*CubeValues(1,0,0)*CubeValues(0,1,0)*CubeValues(0,0,1)* + CubeValues(1,1,0)*CubeValues(1,0,1)*CubeValues(0,1,1)*CubeValues(1,1,1); + // three faces (others shared by other cubes) + epc_nface += CubeValues(1,0,0)*CubeValues(1,1,0)*CubeValues(1,0,1)*CubeValues(1,1,1); + epc_nface += CubeValues(0,1,0)*CubeValues(1,1,0)*CubeValues(0,1,1)*CubeValues(1,1,1); + epc_nface += CubeValues(0,0,1)*CubeValues(0,1,1)*CubeValues(1,0,1)*CubeValues(1,1,1); + // six of twelve edges (others shared by other cubes) + epc_nedge += CubeValues(1,1,1)*CubeValues(1,1,0); + epc_nedge += CubeValues(1,1,1)*CubeValues(1,0,1); + epc_nedge += CubeValues(1,1,1)*CubeValues(0,1,1); + epc_nedge += CubeValues(1,0,1)*CubeValues(1,0,0); + epc_nedge += CubeValues(1,0,1)*CubeValues(0,0,1); + epc_nedge += CubeValues(0,1,1)*CubeValues(0,0,1); + // four of eight vertices + epc_nvert += CubeValues(1,1,0); + epc_nvert += CubeValues(1,0,1); + epc_nvert += CubeValues(0,1,1); + epc_nvert += CubeValues(1,1,1); + + double chi= epc_nvert - epc_nedge + epc_nface - epc_ncubes; + return chi; +} inline double mink_EulerCharacteristic(DTMutableList &Points, IntArray &Triangles, DoubleArray &CubeValues, int &npts, int &ntris, int &i, int &j, int &k){ diff --git a/tests/TestTorus.cpp b/tests/TestTorus.cpp index e936fbe6..220ae377 100644 --- a/tests/TestTorus.cpp +++ b/tests/TestTorus.cpp @@ -27,7 +27,7 @@ int main(int argc, char **argv) if ( rank==0 ) { printf("-----------------------------------------------------------\n"); - printf("Labeling Blobs from Two-Phase Lattice Boltzmann Simulation \n"); + printf("Unit test for torus (Euler-Poincarie characteristic) \n"); printf("-----------------------------------------------------------\n"); }