Working on component labeling

This commit is contained in:
James E McClure 2015-07-11 13:06:57 -04:00
parent 7c9807b7ac
commit 83c0a8f0e7

View File

@ -493,7 +493,7 @@ int main(int argc, char **argv)
int number_NWP_components = ComputeLocalPhaseComponent(PhaseLabel,1,NWP,false);
int number_WP_components = ComputeLocalPhaseComponent(PhaseLabel,2,WP,false);
DoubleArray BlobAverages(NUM_AVERAGES,nblobs);
DoubleArray BlobAverages(NUM_AVERAGES,number_NWP_components);
// Map the signed distance for the analysis
for (i=0; i<Nx*Ny*Nz; i++) SignDist(i) -= (1.0);
@ -513,433 +513,6 @@ int main(int argc, char **argv)
printf("Media porosity is %f \n",porosity);
/* ****************************************************************
RUN TCAT AVERAGING ON EACH BLOB
****************************************************************** */
int n_nw_tris_beg, n_ns_tris_beg, n_ws_tris_beg, n_nws_seg_beg;
int start=0,finish;
int a,c;
printf("-----------------------------------------------\n");
printf("Computing TCAT averages based on connectivity \n");
printf("The number of non-wetting phase features is %i \n",nblobs-1);
printf("-----------------------------------------------\n");
// Wetting phase averages assume global connectivity
As = 0.0;
vol_w = 0.0;
paw = 0.0;
aws = 0.0;
vaw(0) = vaw(1) = vaw(2) = 0.0;
Gws(0) = Gws(1) = Gws(2) = 0.0;
Gws(3) = Gws(4) = Gws(5) = 0.0;
// Don't compute the last blob unless specified
// the last blob is the entire wetting phase
// nblobs -=1;
#ifdef WP
nblobs+=1;
#endif
for (a=0;a<nblobs;a++){
finish = start+b(a);
/* ****************************************************************
RUN PMMC ON EACH BLOB
****************************************************************** */
// Store beginning points for surfaces for blob p
n_nw_tris_beg = n_nw_tris;
n_ns_tris_beg = n_ns_tris;
n_ws_tris_beg = n_ws_tris;
n_nws_seg_beg = n_nws_seg;
// Loop over all cubes
nwp_volume = 0.0;
// Compute phase averages
vol_n =0.0;
pan = 0.0;
awn = ans = lwns = 0.0;
van(0) = van(1) = van(2) = 0.0;
vawn(0) = vawn(1) = vawn(2) = 0.0;
Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
Gns(0) = Gns(1) = Gns(2) = 0.0;
Gns(3) = Gns(4) = Gns(5) = 0.0;
Jwn = Kwn = efawns = 0.0;
trJwn = trawn = trRwn = 0.0;
for (c=start;c<finish;c++){
// Get cube from the list
i = blobs(0,c);
j = blobs(1,c);
k = blobs(2,c);
// Use the cube to compute volume averages
for (p=0;p<8;p++){
if ( SignDist(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
n = i+cube[p][0] + Nx*(j+cube[p][1]) + Nx*Ny*(k+cube[p][2]);
// Compute the non-wetting phase volume contribution
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 )
nwp_volume += 0.125;
// volume averages over the non-wetting phase
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0.99 ){
// volume the excludes the interfacial region
vol_n += 0.125;
// pressure
pan += 0.125*Press(n);
// velocity
van(0) += 0.125*Vel_x(n);
van(1) += 0.125*Vel_y(n);
van(2) += 0.125*Vel_z(n);
}
// volume averages over the wetting phase
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < -0.99 ){
// volume the excludes the interfacial region
vol_w += 0.125;
// pressure
paw += 0.125*Press(n);
// velocity
vaw(0) += 0.125*Vel_x(n);
vaw(1) += 0.125*Vel_y(n);
vaw(2) += 0.125*Vel_z(n);
}
}
}
// Interface and common curve averages
n_local_sol_tris = 0;
n_local_sol_pts = 0;
n_local_nws_pts = 0;
//...........................................................................
// Construct the interfaces and common curve
pmmc_ConstructLocalCube(SignDist, Phase, vS, vF,
nw_pts, nw_tris, values, ns_pts, ns_tris, ws_pts, ws_tris,
local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
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);
// Integrate the contact angle
efawns += pmmc_CubeContactAngle(CubeValues,Values,Phase_x,Phase_y,Phase_z,SignDist_x,SignDist_y,SignDist_z,
local_nws_pts,i,j,k,n_local_nws_pts);
// Integrate the mean curvature
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
Kwn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
// Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels)
pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SignDist,nw_pts,nw_tris,Values,DistValues,
i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,trJwn);
pmmc_CubeTrimSurfaceInterpInverseValues(CubeValues,MeanCurvature,SignDist,nw_pts,nw_tris,Values,DistValues,
i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trRwn);
// Compute the normal speed of the interface
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);
As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);
// Compute the surface orientation and the interfacial area
awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris);
ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris);
aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris);
lwns += pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
//...........................................................................
//*******************************************************************
// Reset the triangle counts to zero
n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0;
n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
n_nw_tris_beg = n_nw_tris;
n_ns_tris_beg = n_ns_tris;
n_ws_tris_beg = n_ws_tris;
n_nws_seg_beg = n_nws_seg;
//*******************************************************************
}
start = finish;
if (a < nblobs-1){
if (vol_n > 0.0){
pan /= vol_n;
for (i=0;i<3;i++) van(i) /= vol_n;
}
if (awn > 0.0){
Jwn /= awn;
Kwn /= awn;
for (i=0;i<3;i++) vawn(i) /= awn;
for (i=0;i<6;i++) Gwn(i) /= awn;
}
if (lwns > 0.0){
efawns /= lwns;
}
if (ans > 0.0){
for (i=0;i<6;i++) Gns(i) /= ans;
}
if (trawn > 0.0){
trJwn /= trawn;
}
BlobAverages(0,a) = nwp_volume;
BlobAverages(1,a) = pan;
BlobAverages(2,a) = awn;
BlobAverages(3,a) = ans;
BlobAverages(4,a) = Jwn;
BlobAverages(5,a) = Kwn;
BlobAverages(6,a) = lwns;
BlobAverages(7,a) = efawns;
BlobAverages(8,a) = van(0);
BlobAverages(9,a) = van(1);
BlobAverages(10,a) = van(2);
BlobAverages(11,a) = vawn(0);
BlobAverages(12,a) = vawn(1);
BlobAverages(13,a) = vawn(2);
BlobAverages(14,a) = Gwn(0);
BlobAverages(15,a) = Gwn(1);
BlobAverages(16,a) = Gwn(2);
BlobAverages(17,a) = Gwn(3);
BlobAverages(18,a) = Gwn(4);
BlobAverages(19,a) = Gwn(5);
BlobAverages(20,a) = Gns(0);
BlobAverages(21,a) = Gns(1);
BlobAverages(22,a) = Gns(2);
BlobAverages(23,a) = Gns(3);
BlobAverages(24,a) = Gns(4);
BlobAverages(25,a) = Gns(5);
BlobAverages(26,a) = trawn;
BlobAverages(27,a) = trJwn;
BlobAverages(28,a) = vol_n;
BlobAverages(29,a) = trRwn;
printf("Computed TCAT averages for feature = %i \n", a);
}
} // End of the blob loop
NULL_USE(n_nw_tris_beg);
NULL_USE(n_ns_tris_beg);
NULL_USE(n_ws_tris_beg);
NULL_USE(n_nws_seg_beg);
nblobs -= 1;
printf("-----------------------------------------------\n");
printf("Sorting the blobs based on volume \n");
printf("-----------------------------------------------\n");
int TempLabel,aa,bb;
double TempValue;
IntArray OldLabel(nblobs);
for (a=0; a<nblobs; a++) OldLabel(a) = a;
// Sort the blob averages based on volume
for (aa=0; aa<nblobs-1; aa++){
for ( bb=aa+1; bb<nblobs; bb++){
if (BlobAverages(0,aa) < BlobAverages(0,bb)){
// Exchange location of blobs aa and bb
//printf("Switch blob %i with %i \n", OldLabel(aa),OldLabel(bb));
// switch the label
TempLabel = OldLabel(bb);
OldLabel(bb) = OldLabel(aa);
OldLabel(aa) = TempLabel;
// switch the averages
for (idx=0; idx<NUM_AVERAGES; idx++){
TempValue = BlobAverages(idx,bb);
BlobAverages(idx,bb) = BlobAverages(idx,aa);
BlobAverages(idx,aa) = TempValue;
}
}
}
}
IntArray NewLabel(nblobs);
for (aa=0; aa<nblobs; aa++){
// Match the new label for original blob aa
bb=0;
while (OldLabel(bb) != aa) bb++;
NewLabel(aa) = bb;
}
// Re-label the blob ID
printf("Re-labeling the blobs, now indexed by volume \n");
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
if (LocalBlobID(i,j,k) > -1){
TempLabel = NewLabel(LocalBlobID(i,j,k));
LocalBlobID(i,j,k) = TempLabel;
}
}
}
}
FILE *BLOBLOG= fopen("blobs.tcat","a");
for (a=0; a<nblobs; a++){
//printf("Blob id =%i \n",a);
//printf("Original Blob id = %i \n",OldLabel(a));
//printf("Blob volume (voxels) = %f \n", BlobAverages(0,a));
for (idx=0; idx<28; idx++){
fprintf(BLOBLOG,"%.8g ",BlobAverages(idx,a));
}
fprintf(BLOBLOG,"\n");
}
fclose(BLOBLOG);
double iVol = 1.0/Nx/Ny/Nz;
sw = 1.0;
// Compute the Sauter mean grain diamter
double D = 6.0*Nx*Ny*Nz*(1.0-porosity) / As;
double pw,pn,pc,awnD,ansD,awsD,JwnD,trJwnD,lwnsDD,cwns;
pw = paw/vol_w;
printf("paw = %f \n", paw/vol_w);
printf("vol_w = %f \n", vol_w);
printf("-----------------------------------------------\n");
double pwn=0.0;
vol_n = nwp_volume = 0.0;
pan = 0.0;
awn = ans = lwns = 0.0;
van(0) = van(1) = van(2) = 0.0;
vawn(0) = vawn(1) = vawn(2) = 0.0;
Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
Gns(0) = Gns(1) = Gns(2) = 0.0;
Gns(3) = Gns(4) = Gns(5) = 0.0;
Jwn = Kwn = efawns = 0.0;
trJwn = trawn = trRwn = 0.0;
// Write out the "equilibrium" state with a 0.5 % change in saturation"
// Always write the largest blob
// sw, pw, pn, pc*D/gamma, awn*D, aws*D, ans*D,lwns*D*D, Jwn*D, trJwn*D, cwns, nblobs
printf("Computing equilibria from blobs \n");
printf("Sauter mean diamter = %f \n",D);
printf("WARNING: lazy coder hard-coded the surface tension as 0.058 \n");
FILE *BLOBSTATES= fopen("blobstates.tcat","a");
fprintf(BLOBSTATES,"%.5g %.5g %.5g\n",vol_w,pw,aws);
// Compute the averages over the entire non-wetting phsae
for (a=0; a<nblobs; a++){
nwp_volume += BlobAverages(0,a);
pwn += (BlobAverages(1,a)-pw)*BlobAverages(2,a);
pan += BlobAverages(1,a)*BlobAverages(28,a);
awn += BlobAverages(2,a);
ans += BlobAverages(3,a);
Jwn += BlobAverages(4,a)*BlobAverages(2,a);
Kwn += BlobAverages(5,a)*BlobAverages(2,a);
lwns += BlobAverages(6,a);
efawns += BlobAverages(7,a)*BlobAverages(6,a);
van(0) += BlobAverages(8,a)*BlobAverages(28,a);
van(1) += BlobAverages(9,a)*BlobAverages(28,a);
van(2) += BlobAverages(10,a)*BlobAverages(28,a);
vawn(0) += BlobAverages(11,a)*BlobAverages(2,a);
vawn(1) += BlobAverages(12,a)*BlobAverages(2,a);
vawn(2) += BlobAverages(13,a)*BlobAverages(2,a);
Gwn(0) += BlobAverages(14,a)*BlobAverages(2,a);
Gwn(1) += BlobAverages(15,a)*BlobAverages(2,a);
Gwn(2) += BlobAverages(16,a)*BlobAverages(2,a);
Gwn(3) += BlobAverages(17,a)*BlobAverages(2,a);
Gwn(4) += BlobAverages(18,a)*BlobAverages(2,a);
Gwn(5) += BlobAverages(19,a)*BlobAverages(2,a);
Gns(0) += BlobAverages(20,a)*BlobAverages(3,a);
Gns(1) += BlobAverages(21,a)*BlobAverages(3,a);
Gns(2) += BlobAverages(22,a)*BlobAverages(3,a);
Gns(3) += BlobAverages(23,a)*BlobAverages(3,a);
Gns(4) += BlobAverages(24,a)*BlobAverages(3,a);
Gns(5) += BlobAverages(25,a)*BlobAverages(3,a);
trawn += BlobAverages(26,a);
trJwn += BlobAverages(27,a)*BlobAverages(26,a);
vol_n += BlobAverages(28,a);
}
// Subtract off portions of non-wetting phase in order of size
for (a=nblobs-1; a>0; a--){
// Subtract the features one-by-one
nwp_volume -= BlobAverages(0,a);
pan -= BlobAverages(1,a)*BlobAverages(28,a);
pwn -= (BlobAverages(1,a)-pw)*BlobAverages(2,a);
awn -= BlobAverages(2,a);
ans -= BlobAverages(3,a);
Jwn -= BlobAverages(4,a)*BlobAverages(2,a);
Kwn -= BlobAverages(5,a)*BlobAverages(2,a);
lwns -= BlobAverages(6,a);
efawns -= BlobAverages(7,a)*BlobAverages(6,a);
van(0) -= BlobAverages(8,a)*BlobAverages(28,a);
van(1) -= BlobAverages(9,a)*BlobAverages(28,a);
van(2) -= BlobAverages(10,a)*BlobAverages(28,a);
vawn(0) -= BlobAverages(11,a)*BlobAverages(2,a);
vawn(1) -= BlobAverages(12,a)*BlobAverages(2,a);
vawn(2) -= BlobAverages(13,a)*BlobAverages(2,a);
Gwn(0) -= BlobAverages(14,a)*BlobAverages(2,a);
Gwn(1) -= BlobAverages(15,a)*BlobAverages(2,a);
Gwn(2) -= BlobAverages(16,a)*BlobAverages(2,a);
Gwn(3) -= BlobAverages(17,a)*BlobAverages(2,a);
Gwn(4) -= BlobAverages(18,a)*BlobAverages(2,a);
Gwn(5) -= BlobAverages(19,a)*BlobAverages(2,a);
Gns(0) -= BlobAverages(20,a)*BlobAverages(3,a);
Gns(1) -= BlobAverages(21,a)*BlobAverages(3,a);
Gns(2) -= BlobAverages(22,a)*BlobAverages(3,a);
Gns(3) -= BlobAverages(23,a)*BlobAverages(3,a);
Gns(4) -= BlobAverages(24,a)*BlobAverages(3,a);
Gns(5) -= BlobAverages(25,a)*BlobAverages(3,a);
trawn -= BlobAverages(26,a);
trJwn -= BlobAverages(27,a)*BlobAverages(26,a);
vol_n -= BlobAverages(28,a);
// Update wetting phase averages
aws += BlobAverages(3,a);
Gws(0) += BlobAverages(20,a)*BlobAverages(3,a);
Gws(1) += BlobAverages(21,a)*BlobAverages(3,a);
Gws(2) += BlobAverages(22,a)*BlobAverages(3,a);
Gws(3) += BlobAverages(23,a)*BlobAverages(3,a);
Gws(4) += BlobAverages(24,a)*BlobAverages(3,a);
Gws(5) += BlobAverages(25,a)*BlobAverages(3,a);
if (fabs(1.0 - nwp_volume*iVol/porosity - sw) > 0.0025 || a == 1){
sw = 1.0 - nwp_volume*iVol/porosity;
JwnD = -Jwn*D/awn;
pn = pan/vol_n;
trJwnD = -trJwn*D/trawn;
cwns = -efawns / lwns;
awnD = awn*D*iVol;
awsD = aws*D*iVol;
ansD = ans*D*iVol;
lwnsDD = lwns*D*D*iVol;
pc = pwn*D/0.058/awn; // hard-coded surface tension due to being lazy
fprintf(BLOBSTATES,"%.5g %.5g %.5g ",sw,pn,pw);
fprintf(BLOBSTATES,"%.5g %.5g %.5g %.5g ",awnD,awsD,ansD,lwnsDD);
fprintf(BLOBSTATES,"%.5g %.5g %.5g %.5g %i\n",pc,JwnD,trJwnD,cwns,a);
}
}
fclose(BLOBSTATES);
start = 0;
for (a=0;a<nblobs;a++){
finish = start+b(a);
for (c=start;c<finish;c++){
// Get cube from the list
i = blobs(0,c);
j = blobs(1,c);
k = blobs(2,c);
// Label the entire cube so that interfaces can be re-labled easily
LocalBlobID(i,j,k) = NewLabel(a);
LocalBlobID(i+1,j,k) = NewLabel(a);
LocalBlobID(i,j+1,k) = NewLabel(a);
LocalBlobID(i+1,j+1,k) = NewLabel(a);
LocalBlobID(i,j,k+1) = NewLabel(a);
LocalBlobID(i+1,j,k+1) = NewLabel(a);
LocalBlobID(i,j+1,k+1) = NewLabel(a);
LocalBlobID(i+1,j+1,k+1) = NewLabel(a);
}
start=finish;
}
FILE *BLOBS;
BLOBS = fopen("Blobs.dat","wb");