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

Refactoring Euler characteristic

Browse Source
This commit is contained in:
James E McClure 2015-11-22 14:39:19 -05:00
parent 2e874fc41a
commit 80e3aba9a6
3 changed files with 52 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
common/TwoPhase.cpp
View File

@ -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 // Construct the interfaces and common curve
pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue, 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, Kn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,
i,j,k,n_nw_pts,n_nw_tris); 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);
} }
} }

46
common/pmmc.h
View File

@ -4427,6 +4427,52 @@ inline double geomavg_EulerCharacteristic(DTMutableList<Point> &Points, IntArray
return EulerChar; 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<Point> &Points, IntArray &Triangles, inline double mink_EulerCharacteristic(DTMutableList<Point> &Points, IntArray &Triangles,
DoubleArray &CubeValues, int &npts, int &ntris, int &i, int &j, int &k){ DoubleArray &CubeValues, int &npts, int &ntris, int &i, int &j, int &k){

2
tests/TestTorus.cpp
View File

@ -27,7 +27,7 @@ int main(int argc, char **argv)
if ( rank==0 ) { if ( rank==0 ) {
printf("-----------------------------------------------------------\n"); printf("-----------------------------------------------------------\n");
printf("Labeling Blobs from Two-Phase Lattice Boltzmann Simulation \n"); printf("Unit test for torus (Euler-Poincarie characteristic) \n");
printf("-----------------------------------------------------------\n"); printf("-----------------------------------------------------------\n");
} }