OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
27898d890ff83b1341c3d1cc4654a89d96720523
LBPM/analysis/TwoPhase.cpp
Mark Berrill 27898d890f Fixing bug in Domain
2018-05-16 14:00:35 -04:00

1500 lines
56 KiB
C++
Raw Blame History

#include "analysis/TwoPhase.h"
#include "analysis/pmmc.h"
#include "common/Domain.h"
#include "common/Communication.h"
#include "analysis/analysis.h"
#include "shared_ptr.h"
#include "common/Utilities.h"
#include "common/MPI_Helpers.h"
#include "IO/MeshDatabase.h"
#include "IO/Reader.h"
#include "IO/Writer.h"
#define BLOB_AVG_COUNT 35
// Array access for averages defined by the following
#define VOL 0
#define TRIMVOL 1
#define PRS 2
#define AWN 3
#define AWS 4
#define ANS 5
#define LWNS 6
#define JWN 7
#define KWN 8
#define CWNS 9
#define KNWNS 10
#define KGWNS 11
#define VX 12
#define VY 13
#define VZ 14
#define VSQ 15
#define VWNX 16
#define VWNY 17
#define VWNZ 18
#define VWNSX 19
#define VWNSY 20
#define VWNSZ 21
#define GWNXX 22
#define GWNYY 23
#define GWNZZ 24
#define GWNXY 25
#define GWNXZ 26
#define GWNYZ 27
#define TRAWN 28
#define TRJWN 29
#define CMX 30
#define CMY 31
#define CMZ 32
#define EULER 33
#define INTCURV 34
#define PI 3.14159265359
// Constructor
TwoPhase::TwoPhase(Domain &dm):
n_nw_pts(0), n_ns_pts(0), n_ws_pts(0), n_nws_pts(0), n_local_sol_pts(0), n_local_nws_pts(0),
n_nw_tris(0), n_ns_tris(0), n_ws_tris(0), n_nws_seg(0), n_local_sol_tris(0),
nc(0), kstart(0), kfinish(0), fluid_isovalue(0), solid_isovalue(0), Volume(0),
TIMELOG(NULL), NWPLOG(NULL), WPLOG(NULL),
Dm(dm), NumberComponents_WP(0), NumberComponents_NWP(0), trimdist(0),
porosity(0), poreVol(0), awn(0), ans(0), aws(0), lwns(0), wp_volume(0), nwp_volume(0),
As(0), dummy(0), vol_w(0), vol_n(0), sat_w(0), sat_w_previous(0),
pan(0), paw(0), pan_global(0), paw_global(0), vol_w_global(0), vol_n_global(0),
awn_global(0), ans_global(0),aws_global(0), lwns_global(0), efawns(0), efawns_global(0),
Jwn(0), Jwn_global(0), Kwn(0), Kwn_global(0), KNwns(0), KNwns_global(0),
KGwns(0), KGwns_global(0), trawn(0), trawn_global(0), trJwn(0), trJwn_global(0),
trRwn(0), trRwn_global(0), nwp_volume_global(0), wp_volume_global(0),
As_global(0), wwndnw_global(0), wwnsdnwn_global(0), Jwnwwndnw_global(0), dEs(0), dAwn(0), dAns(0)
{
Nx=dm.Nx; Ny=dm.Ny; Nz=dm.Nz;
Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm.nprocx()*Dm.nprocy()*Dm.nprocz()*1.0;
TempID = new char[Nx*Ny*Nz];
// Global arrays
PhaseID.resize(Nx,Ny,Nz); PhaseID.fill(0);
Label_WP.resize(Nx,Ny,Nz); Label_WP.fill(0);
Label_NWP.resize(Nx,Ny,Nz); Label_NWP.fill(0);
SDn.resize(Nx,Ny,Nz); SDn.fill(0);
SDs.resize(Nx,Ny,Nz); SDs.fill(0);
Phase.resize(Nx,Ny,Nz); Phase.fill(0);
Press.resize(Nx,Ny,Nz); Press.fill(0);
dPdt.resize(Nx,Ny,Nz); dPdt.fill(0);
MeanCurvature.resize(Nx,Ny,Nz); MeanCurvature.fill(0);
GaussCurvature.resize(Nx,Ny,Nz); GaussCurvature.fill(0);
SDs_x.resize(Nx,Ny,Nz); SDs_x.fill(0); // Gradient of the signed distance
SDs_y.resize(Nx,Ny,Nz); SDs_y.fill(0);
SDs_z.resize(Nx,Ny,Nz); SDs_z.fill(0);
SDn_x.resize(Nx,Ny,Nz); SDn_x.fill(0); // Gradient of the signed distance
SDn_y.resize(Nx,Ny,Nz); SDn_y.fill(0);
SDn_z.resize(Nx,Ny,Nz); SDn_z.fill(0);
DelPhi.resize(Nx,Ny,Nz); DelPhi.fill(0);
Phase_tplus.resize(Nx,Ny,Nz); Phase_tplus.fill(0);
Phase_tminus.resize(Nx,Ny,Nz); Phase_tminus.fill(0);
Vel_x.resize(Nx,Ny,Nz); Vel_x.fill(0); // Gradient of the phase indicator field
Vel_y.resize(Nx,Ny,Nz); Vel_y.fill(0);
Vel_z.resize(Nx,Ny,Nz); Vel_z.fill(0);
//.........................................
// Allocate cube storage space
CubeValues.resize(2,2,2);
nw_tris.resize(3,20);
ns_tris.resize(3,20);
ws_tris.resize(3,20);
nws_seg.resize(2,20);
local_sol_tris.resize(3,18);
nw_pts=DTMutableList<Point>(20);
ns_pts=DTMutableList<Point>(20);
ws_pts=DTMutableList<Point>(20);
nws_pts=DTMutableList<Point>(20);
local_nws_pts=DTMutableList<Point>(20);
local_sol_pts=DTMutableList<Point>(20);
tmp=DTMutableList<Point>(20);
//.........................................
Values.resize(20);
DistanceValues.resize(20);
KGwns_values.resize(20);
KNwns_values.resize(20);
InterfaceSpeed.resize(20);
NormalVector.resize(60);
//.........................................
van.resize(3);
vaw.resize(3);
vawn.resize(3);
vawns.resize(3);
Gwn.resize(6);
Gns.resize(6);
Gws.resize(6);
van_global.resize(3);
vaw_global.resize(3);
vawn_global.resize(3);
vawns_global.resize(3);
Gwn_global.resize(6);
Gns_global.resize(6);
Gws_global.resize(6);
//.........................................
if (Dm.rank()==0){
TIMELOG = fopen("timelog.tcat","a+");
if (fseek(TIMELOG,0,SEEK_SET) == fseek(TIMELOG,0,SEEK_CUR))
{
// If timelog is empty, write a short header to list the averages
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"time dEs "); // Timestep, Change in Surface Energy
fprintf(TIMELOG,"sw pw pn awn ans aws Jwn Kwn lwns cwns KNwns KGwns "); // Scalar averages
fprintf(TIMELOG,"vawx vawy vawz vanx vany vanz "); // Velocity averages
fprintf(TIMELOG,"vawnx vawny vawnz vawnsx vawnsy vawnsz ");
fprintf(TIMELOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors
fprintf(TIMELOG,"Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz ");
fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
fprintf(TIMELOG,"trawn trJwn trRwn "); //trimmed curvature,
fprintf(TIMELOG,"wwndnw wwnsdnwn Jwnwwndnw "); //kinematic quantities,
fprintf(TIMELOG,"Euler Kn Jn An\n"); //miknowski measures,
}
NWPLOG = fopen("components.NWP.tcat","a+");
fprintf(NWPLOG,"time label vol pn awn ans Jwn Kwn lwns cwns ");
fprintf(NWPLOG,"vx vy vz vwnx vwny vwnz vwnsx vwnsy vwnsz vsq ");
fprintf(NWPLOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz Cx Cy Cz trawn trJwn Kn Euler\n");
WPLOG = fopen("components.WP.tcat","a+");
fprintf(WPLOG,"time label vol pw awn ans Jwn Kwn lwns cwns ");
fprintf(WPLOG,"vx vy vz vwnx vwny vwnz vwnsx vwnsy vwnsz vsq ");
fprintf(WPLOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz Cx Cy Cz trawn trJwn\n");
}
else{
char LocalRankString[8];
sprintf(LocalRankString,"%05d",Dm.rank());
char LocalRankFilename[40];
sprintf(LocalRankFilename,"%s%s","timelog.tcat.",LocalRankString);
TIMELOG = fopen(LocalRankFilename,"a+");
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"time "); // Timestep, Change in Surface Energy
fprintf(TIMELOG,"sw pw pn awn ans aws Jwn Kwn lwns cwns KNwns KGwns "); // Scalar averages
fprintf(TIMELOG,"vawx vawy vawz vanx vany vanz "); // Velocity averages
fprintf(TIMELOG,"vawnx vawny vawnz vawnsx vawnsy vawnsz ");
fprintf(TIMELOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors
fprintf(TIMELOG,"Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz ");
fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
fprintf(TIMELOG,"trawn trJwn trRwn "); //trimmed curvature,
fprintf(TIMELOG,"wwndnw wwnsdnwn Jwnwwndnw "); //kinematic quantities,
fprintf(TIMELOG,"Euler Kn Jn An\n"); //miknowski measures,
}
}
// Destructor
TwoPhase::~TwoPhase()
{
delete [] TempID;
if ( TIMELOG!=NULL ) { fclose(TIMELOG); }
if ( NWPLOG!=NULL ) { fclose(NWPLOG); }
if ( WPLOG!=NULL ) { fclose(WPLOG); }
}
void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, DoubleArray &DistData)
{
double factor,temp,value;
/* factor=0.5/Beta;
// Initialize to -1,1 (segmentation)
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
value = ColorData(i,j,k);
// Set phase ID field for non-wetting phase
// here NWP is tagged with 1
// all other phases tagged with 0
int n = k*Nx*Ny+j*Nx+i;
if (value > 0) TempID[n] = 1;
else TempID[n] = 0;
// Distance threshhold
// temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm.
// distance should be negative outside the NWP
// distance should be positive inside of the NWP
temp = factor*log((1.0+value)/(1.0-value));
if (value > 0.8) DistData(i,j,k) = 2.94*factor;
else if (value < -0.8) DistData(i,j,k) = -2.94*factor;
else DistData(i,j,k) = temp;
// Basic threshold
// distance to the NWP
// negative inside NWP, positive outside
//if (value > 0) DistData(i,j,k) = -0.5;
//else DistData(i,j,k) = 0.5;
// Initialize directly
// DistData(i,j,k) = value;
}
}
}
Eikonal(DistData,TempID,Dm,50);
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
DistData(i,j,k) += 1.0;
}
}
}
*/
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
DistData(i,j,k) = ColorData(i,j,k);
}
}
}
}
void TwoPhase::ComputeDelPhi()
{
int i,j,k;
double fx,fy,fz;
Dm.CommunicateMeshHalo(Phase);
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
// Compute all of the derivatives using finite differences
fx = 0.5*(Phase(i+1,j,k) - Phase(i-1,j,k));
fy = 0.5*(Phase(i,j+1,k) - Phase(i,j-1,k));
fz = 0.5*(Phase(i,j,k+1) - Phase(i,j,k-1));
DelPhi(i,j,k) = sqrt(fx*fx+fy*fy+fz*fz);
}
}
}
}
void TwoPhase::Initialize()
{
trimdist=1.0;
fluid_isovalue=solid_isovalue=0.0;
// Initialize the averaged quantities
awn = aws = ans = lwns = 0.0;
nwp_volume = wp_volume = 0.0;
As = 0.0;
pan = paw = 0.0;
vaw(0) = vaw(1) = vaw(2) = 0.0;
van(0) = van(1) = van(2) = 0.0;
vawn(0) = vawn(1) = vawn(2) = 0.0;
vawns(0) = vawns(1) = vawns(2) = 0.0;
Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
Gws(0) = Gws(1) = Gws(2) = 0.0;
Gws(3) = Gws(4) = Gws(5) = 0.0;
Gns(0) = Gns(1) = Gns(2) = 0.0;
Gns(3) = Gns(4) = Gns(5) = 0.0;
vol_w = vol_n =0.0;
KGwns = KNwns = 0.0;
Jwn = Kwn = efawns = 0.0;
trJwn = trawn = trRwn = 0.0;
euler = Jn = An = Kn = 0.0;
wwndnw = 0.0; wwnsdnwn = 0.0; Jwnwwndnw=0.0;
}
/*
void TwoPhase::SetupCubes(Domain &Dm)
{
int i,j,k;
kstart = 1;
kfinish = Nz-1;
if (Dm.BoundaryCondition !=0 && Dm.kproc==0) kstart = 4;
if (Dm.BoundaryCondition !=0 && Dm.kproc==Dm.nprocz-1) kfinish = Nz-4;
nc=0;
for (k=kstart; k<kfinish; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
cubeList(0,nc) = i;
cubeList(1,nc) = j;
cubeList(2,nc) = k;
nc++;
}
}
}
ncubes = nc;
}
*/
void TwoPhase::UpdateSolid()
{
Dm.CommunicateMeshHalo(SDs);
//...........................................................................
// Gradient of the Signed Distance function
//...........................................................................
pmmc_MeshGradient(SDs,SDs_x,SDs_y,SDs_z,Nx,Ny,Nz);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_x);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_y);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_z);
//...........................................................................
}
void TwoPhase::UpdateMeshValues()
{
int i,j,k,n;
//...........................................................................
Dm.CommunicateMeshHalo(SDn);
//...........................................................................
// Compute the gradients of the phase indicator and signed distance fields
pmmc_MeshGradient(SDn,SDn_x,SDn_y,SDn_z,Nx,Ny,Nz);
//...........................................................................
// Gradient of the phase indicator field
//...........................................................................
Dm.CommunicateMeshHalo(SDn_x);
//...........................................................................
Dm.CommunicateMeshHalo(SDn_y);
//...........................................................................
Dm.CommunicateMeshHalo(SDn_z);
//...........................................................................
Dm.CommunicateMeshHalo(SDs);
pmmc_MeshGradient(SDs,SDs_x,SDs_y,SDs_z,Nx,Ny,Nz);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_x);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_y);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_z);
//...........................................................................
// Compute the mesh curvature of the phase indicator field
pmmc_MeshCurvature(SDn, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
//...........................................................................
// Update the time derivative of non-dimensional density field
// Map Phase_tplus and Phase_tminus
for (int n=0; n<Nx*Ny*Nz; n++) dPdt(n) = 0.125*(Phase_tplus(n) - Phase_tminus(n));
//...........................................................................
Dm.CommunicateMeshHalo(Press);
//...........................................................................
Dm.CommunicateMeshHalo(Vel_x);
//...........................................................................
Dm.CommunicateMeshHalo(Vel_y);
//...........................................................................
Dm.CommunicateMeshHalo(Vel_z);
//...........................................................................
Dm.CommunicateMeshHalo(MeanCurvature);
//...........................................................................
Dm.CommunicateMeshHalo(GaussCurvature);
//...........................................................................
Dm.CommunicateMeshHalo(DelPhi);
//...........................................................................
// Initializing the blob ID
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
n = k*Nx*Ny+j*Nx+i;
if (Dm.id[n] == 0){
// Solid phase
PhaseID(i,j,k) = 0;
}
else if (SDn(i,j,k) < 0.0){
// wetting phase
PhaseID(i,j,k) = 2;
}
else {
// non-wetting phase
PhaseID(i,j,k) = 1;
}
}
}
}
}
void TwoPhase::ComputeLocal()
{
int i,j,k,n,kmin,kmax;
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}};
// If external boundary conditions are set, do not average over the inlet
kmin=1; kmax=Nz-1;
if (Dm.BoundaryCondition > 0 && Dm.kproc() == 0) kmin=4;
if (Dm.BoundaryCondition > 0 && Dm.kproc() == Dm.nprocz()-1) kmax=Nz-4;
for (k=kmin; k<kmax; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
//...........................................................................
n_nw_pts=n_ns_pts=n_ws_pts=n_nws_pts=n_local_sol_pts=n_local_nws_pts=0;
n_nw_tris=n_ns_tris=n_ws_tris=n_nws_seg=n_local_sol_tris=0;
//...........................................................................
// Compute volume averages
for (int p=0;p<8;p++){
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
if ( Dm.id[n] != 0 ){
// 1-D index for this cube corner
// compute the norm of the gradient of the phase indicator field
// 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;
// velocity
van(0) += 0.125*Vel_x(n);
van(1) += 0.125*Vel_y(n);
van(2) += 0.125*Vel_z(n);
// volume the excludes the interfacial region
if (DelPhi(n) < 1e-4){
vol_n += 0.125;
// pressure
pan += 0.125*Press(n);
}
}
else{
wp_volume += 0.125;
// velocity
vaw(0) += 0.125*Vel_x(n);
vaw(1) += 0.125*Vel_y(n);
vaw(2) += 0.125*Vel_z(n);
if (DelPhi(n) < 1e-4){
// volume the excludes the interfacial region
vol_w += 0.125;
// pressure
paw += 0.125*Press(n);
}
}
}
}
//...........................................................................
// Construct the interfaces and common curve
pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue,
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);
// wn interface averages
if (n_nw_pts > 0){
awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris);
Kwn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
// Compute the normal speed of the interface
wwndnw += pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris,
NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
//for (int p=0; p <n_nw_tris; p++) wwndnw += InterfaceSpeed(p);
for (int p=0; p <n_nw_tris; p++) Jwnwwndnw += InterfaceSpeed(p)*Values(p);
// Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels)
pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues,
i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trJwn);
pmmc_CubeTrimSurfaceInterpInverseValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues,
i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trRwn);
}
// wns common curve averages
if (n_local_nws_pts > 0){
efawns += pmmc_CubeContactAngle(CubeValues,Values,SDn_x,SDn_y,SDn_z,SDs_x,SDs_y,SDs_z,
local_nws_pts,i,j,k,n_local_nws_pts);
wwnsdnwn += pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z,
local_nws_pts,i,j,k,n_local_nws_pts);
pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
SDs_z, KNwns_values, KGwns_values, KNwns, KGwns,
nws_pts, n_nws_pts, i, j, k);
lwns += pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
}
// Solid interface averagees
if (n_local_sol_tris > 0){
As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);
// Compute the surface orientation and the interfacial area
ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris);
aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris);
}
//...........................................................................
// Compute the integral curvature of the non-wetting phase
n_nw_pts=n_nw_tris=0;
// Compute the non-wetting phase surface and associated area
An += geomavg_MarchingCubes(SDn,fluid_isovalue,i,j,k,nw_pts,n_nw_pts,nw_tris,n_nw_tris);
// Compute the integral of mean curvature
if (n_nw_pts>0){
pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues,
i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,dummy);
}
Jn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,
i,j,k,n_nw_pts,n_nw_tris);
// Compute Euler characteristic from integral of gaussian curvature
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);
}
}
}
}
void TwoPhase::AssignComponentLabels()
{
int LabelNWP=1;
int LabelWP=2;
// NOTE: labeling the wetting phase components is tricky! One sandstone media had over 800,000 components
// NumberComponents_WP = ComputeGlobalPhaseComponent(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,PhaseID,LabelWP,Label_WP);
// treat all wetting phase is connected
NumberComponents_WP=1;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
Label_WP(i,j,k) = 0;
//if (SDs(i,j,k) > 0.0) PhaseID(i,j,k) = 0;
//else if (Phase(i,j,k) > 0.0) PhaseID(i,j,k) = LabelNWP;
//else PhaseID(i,j,k) = LabelWP;
}
}
}
// Fewer non-wetting phase features are present
//NumberComponents_NWP = ComputeGlobalPhaseComponent(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,PhaseID,LabelNWP,Label_NWP);
NumberComponents_NWP = ComputeGlobalBlobIDs(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,SDs,SDn,solid_isovalue,fluid_isovalue,Label_NWP,Dm.Comm);
}
void TwoPhase::ComponentAverages()
{
int i,j,k,n;
int kmin,kmax;
int LabelWP,LabelNWP;
double TempLocal;
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}};
ComponentAverages_WP.resize(BLOB_AVG_COUNT,NumberComponents_WP);
ComponentAverages_NWP.resize(BLOB_AVG_COUNT,NumberComponents_NWP);
ComponentAverages_WP.fill(0.0);
ComponentAverages_NWP.fill(0.0);
if (Dm.rank()==0){
printf("Number of wetting phase components is %i \n",NumberComponents_WP);
printf("Number of non-wetting phase components is %i \n",NumberComponents_NWP);
}
// If external boundary conditions are set, do not average over the inlet
kmin=1; kmax=Nz-1;
if (Dm.BoundaryCondition > 0 && Dm.kproc() == 0) kmin=4;
if (Dm.BoundaryCondition > 0 && Dm.kproc() == Dm.nprocz()-1) kmax=Nz-4;
for (k=kmin; k<kmax; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
LabelWP=GetCubeLabel(i,j,k,Label_WP);
LabelNWP=GetCubeLabel(i,j,k,Label_NWP);
n_nw_pts=n_ns_pts=n_ws_pts=n_nws_pts=n_local_sol_pts=n_local_nws_pts=0;
n_nw_tris=n_ns_tris=n_ws_tris=n_nws_seg=n_local_sol_tris=0;
// Initialize the averaged quantities
awn = aws = ans = lwns = 0.0;
vawn(0) = vawn(1) = vawn(2) = 0.0;
vawns(0) = vawns(1) = vawns(2) = 0.0;
Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
Gws(0) = Gws(1) = Gws(2) = 0.0;
Gws(3) = Gws(4) = Gws(5) = 0.0;
Gns(0) = Gns(1) = Gns(2) = 0.0;
Gns(3) = Gns(4) = Gns(5) = 0.0;
KGwns = KNwns = 0.0;
Jwn = Kwn = efawns = 0.0;
trawn=trJwn=0.0;
euler=0.0;
//...........................................................................
//...........................................................................
// Compute volume averages
for (int p=0;p<8;p++){
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
if ( Dm.id[n] != 0 ){
// 1-D index for this cube corner
// compute the norm of the gradient of the phase indicator field
// Compute the non-wetting phase volume contribution
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0.0 && !(LabelNWP < 0) ){
// volume
ComponentAverages_NWP(VOL,LabelNWP) += 0.125;
// velocity
ComponentAverages_NWP(VX,LabelNWP) += 0.125*Vel_x(n);
ComponentAverages_NWP(VY,LabelNWP) += 0.125*Vel_y(n);
ComponentAverages_NWP(VZ,LabelNWP) += 0.125*Vel_z(n);
// center of mass
ComponentAverages_NWP(CMX,LabelNWP) += 0.125*(i+cube[p][0]+Dm.iproc()*Nx);
ComponentAverages_NWP(CMY,LabelNWP) += 0.125*(j+cube[p][1]+Dm.jproc()*Ny);
ComponentAverages_NWP(CMZ,LabelNWP) += 0.125*(k+cube[p][2]+Dm.kproc()*Nz);
// twice the kinetic energy
ComponentAverages_NWP(VSQ,LabelNWP) += 0.125*(Vel_x(n)*Vel_x(n)+Vel_y(n)*Vel_y(n)+Vel_z(n)*Vel_z(n));
// volume the for pressure averaging excludes the interfacial region
if (DelPhi(n) < 1e-4 ){
ComponentAverages_NWP(TRIMVOL,LabelNWP) += 0.125;
ComponentAverages_NWP(PRS,LabelNWP) += 0.125*Press(n);
}
}
else if (!(LabelWP < 0)){
ComponentAverages_WP(VOL,LabelWP) += 0.125;
// velocity
ComponentAverages_WP(VX,LabelWP) += 0.125*Vel_x(n);
ComponentAverages_WP(VY,LabelWP)+= 0.125*Vel_y(n);
ComponentAverages_WP(VZ,LabelWP) += 0.125*Vel_z(n);
// Center of mass
ComponentAverages_WP(CMX,LabelWP) += 0.125*(i+cube[p][0]+Dm.iproc()*Nx);
ComponentAverages_WP(CMY,LabelWP) += 0.125*(j+cube[p][1]+Dm.jproc()*Ny);
ComponentAverages_WP(CMZ,LabelWP) += 0.125*(k+cube[p][2]+Dm.kproc()*Nz);
// twice the kinetic energy
ComponentAverages_WP(VSQ,LabelWP) += 0.125*(Vel_x(n)*Vel_x(n)+Vel_y(n)*Vel_y(n)+Vel_z(n)*Vel_z(n));
// volume the for pressure averaging excludes the interfacial region
if (DelPhi(n) < 1e-4){
ComponentAverages_WP(TRIMVOL,LabelWP) += 0.125;
ComponentAverages_WP(PRS,LabelWP) += 0.125*Press(n);
}
}
}
}
//...........................................................................
// Construct the interfaces and common curve
pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue,
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);
//...........................................................................
// wn interface averages
if (n_nw_pts>0 && LabelNWP >=0 && LabelWP >=0 ){
// Mean curvature
TempLocal = pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
ComponentAverages_WP(JWN,LabelWP) += TempLocal;
ComponentAverages_NWP(JWN,LabelNWP) += TempLocal;
// Trimmed Mean curvature
pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues,
i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,trJwn);
ComponentAverages_WP(TRAWN,LabelWP) += trawn;
ComponentAverages_WP(TRJWN,LabelWP) += trJwn;
ComponentAverages_NWP(TRAWN,LabelNWP) += trawn;
ComponentAverages_NWP(TRJWN,LabelNWP) += trJwn;
// Gaussian curvature
TempLocal = pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
ComponentAverages_WP(KWN,LabelWP) += TempLocal;
ComponentAverages_NWP(KWN,LabelNWP) += TempLocal;
// Compute the normal speed of the interface
pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris,
NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
ComponentAverages_WP(VWNX,LabelWP) += vawn(0);
ComponentAverages_WP(VWNY,LabelWP) += vawn(1);
ComponentAverages_WP(VWNZ,LabelWP) += vawn(2);
ComponentAverages_NWP(VWNX,LabelNWP) += vawn(0);
ComponentAverages_NWP(VWNY,LabelNWP) += vawn(1);
ComponentAverages_NWP(VWNZ,LabelNWP) += vawn(2);
// Interfacial Area
TempLocal = pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris);
ComponentAverages_WP(AWN,LabelWP) += TempLocal;
ComponentAverages_NWP(AWN,LabelNWP) += TempLocal;
ComponentAverages_WP(GWNXX,LabelWP) += Gwn(0);
ComponentAverages_WP(GWNYY,LabelWP) += Gwn(1);
ComponentAverages_WP(GWNZZ,LabelWP) += Gwn(2);
ComponentAverages_WP(GWNXY,LabelWP) += Gwn(3);
ComponentAverages_WP(GWNXZ,LabelWP) += Gwn(4);
ComponentAverages_WP(GWNYZ,LabelWP) += Gwn(5);
ComponentAverages_NWP(GWNXX,LabelNWP) += Gwn(0);
ComponentAverages_NWP(GWNYY,LabelNWP) += Gwn(1);
ComponentAverages_NWP(GWNZZ,LabelNWP) += Gwn(2);
ComponentAverages_NWP(GWNXY,LabelNWP) += Gwn(3);
ComponentAverages_NWP(GWNXZ,LabelNWP) += Gwn(4);
ComponentAverages_NWP(GWNYZ,LabelNWP) += Gwn(5);
}
//...........................................................................
// Common curve averages
if (n_local_nws_pts > 0 && LabelNWP >=0 && LabelWP >=0){
// Contact angle
TempLocal = pmmc_CubeContactAngle(CubeValues,Values,SDn_x,SDn_y,SDn_z,SDs_x,SDs_y,SDs_z,
local_nws_pts,i,j,k,n_local_nws_pts);
ComponentAverages_WP(CWNS,LabelWP) += TempLocal;
ComponentAverages_NWP(CWNS,LabelNWP) += TempLocal;
// Kinematic velocity of the common curve
pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z,
local_nws_pts,i,j,k,n_local_nws_pts);
ComponentAverages_WP(VWNSX,LabelWP) += vawns(0);
ComponentAverages_WP(VWNSY,LabelWP) += vawns(1);
ComponentAverages_WP(VWNSZ,LabelWP) += vawns(2);
ComponentAverages_NWP(VWNSX,LabelNWP) += vawns(0);
ComponentAverages_NWP(VWNSY,LabelNWP) += vawns(1);
ComponentAverages_NWP(VWNSZ,LabelNWP) += vawns(2);
// Curvature of the common curve
pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
SDs_z, KNwns_values, KGwns_values, KNwns, KGwns,
nws_pts, n_nws_pts, i, j, k);
ComponentAverages_WP(KNWNS,LabelWP) += KNwns;
ComponentAverages_WP(KGWNS,LabelWP) += KGwns;
ComponentAverages_NWP(KNWNS,LabelNWP) += KNwns;
ComponentAverages_NWP(KGWNS,LabelNWP) += KGwns;
// Length of the common curve
TempLocal = pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
ComponentAverages_NWP(LWNS,LabelNWP) += TempLocal;
}
//...........................................................................
// Solid interface averages
if (n_local_sol_pts > 0 && LabelWP >=0){
As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);
// Compute the surface orientation and the interfacial area
TempLocal = pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris);
ComponentAverages_WP(AWS,LabelWP) += TempLocal;
}
if (n_ns_pts > 0 && LabelNWP >=0 ){
TempLocal = pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris);
ComponentAverages_NWP(ANS,LabelNWP) += TempLocal;
}
//...........................................................................
/* Compute the Euler characteristic
* count all vertices, edges and faces (triangles)
* Euler Number = vertices - edges + faces
* double geomavg_EulerCharacteristic(PointList, PointCount, TriList, TriCount);
*/
if (LabelNWP >= 0 ){
// find non-wetting phase interface
// 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);
// Compute Euler characteristic from integral of gaussian curvature
euler = pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,
i,j,k,n_nw_pts,n_nw_tris);
ComponentAverages_NWP(INTCURV,LabelNWP) += euler;
// Compute the Euler characteristic from vertices - faces + edges
euler = geomavg_EulerCharacteristic(nw_pts,nw_tris,n_nw_pts,n_nw_tris,i,j,k);
ComponentAverages_NWP(EULER,LabelNWP) += euler;
}
}
}
}
MPI_Barrier(Dm.Comm);
if (Dm.rank()==0){
printf("Component averages computed locally -- reducing result... \n");
}
// Globally reduce the non-wetting phase averages
RecvBuffer.resize(BLOB_AVG_COUNT,NumberComponents_NWP);
/* for (int b=0; b<NumberComponents_NWP; b++){
MPI_Barrier(Dm.Comm);
MPI_Allreduce(&ComponentAverages_NWP(0,b),&RecvBuffer(0),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,Dm.Comm);
for (int idx=0; idx<BLOB_AVG_COUNT; idx++) ComponentAverages_NWP(idx,b)=RecvBuffer(idx);
}
*/
MPI_Barrier(Dm.Comm);
MPI_Allreduce(ComponentAverages_NWP.data(),RecvBuffer.data(),BLOB_AVG_COUNT*NumberComponents_NWP, MPI_DOUBLE,MPI_SUM,Dm.Comm);
// MPI_Reduce(ComponentAverages_NWP.data(),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
if (Dm.rank()==0){
printf("rescaling... \n");
}
for (int b=0; b<NumberComponents_NWP; b++){
for (int idx=0; idx<BLOB_AVG_COUNT; idx++) ComponentAverages_NWP(idx,b)=RecvBuffer(idx,b);
}
for (int b=0; b<NumberComponents_NWP; b++){
if (ComponentAverages_NWP(VOL,b) > 0.0){
double Vn,pn,awn,ans,Jwn,Kwn,lwns,cwns,vsq;
Vn = ComponentAverages_NWP(VOL,b);
awn = ComponentAverages_NWP(AWN,b);
ans = ComponentAverages_NWP(ANS,b);
van(0) = ComponentAverages_NWP(VX,b)/Vn;
van(1) = ComponentAverages_NWP(VY,b)/Vn;
van(2) = ComponentAverages_NWP(VZ,b)/Vn;
vsq = ComponentAverages_NWP(VSQ,b)/Vn;
NULL_USE(ans);
if (ComponentAverages_NWP(TRIMVOL,b) > 0.0){
pn = ComponentAverages_NWP(PRS,b)/ComponentAverages_NWP(TRIMVOL,b);
}
else pn = 0.0;
if (awn != 0.0){
Jwn = ComponentAverages_NWP(JWN,b)/awn;
Kwn = ComponentAverages_NWP(KWN,b)/awn;
vawn(0) = ComponentAverages_NWP(VWNSX,b)/awn;
vawn(1) = ComponentAverages_NWP(VWNSY,b)/awn;
vawn(2) = ComponentAverages_NWP(VWNSZ,b)/awn;
Gwn(0) = ComponentAverages_NWP(GWNXX,b)/awn;
Gwn(1) = ComponentAverages_NWP(GWNYY,b)/awn;
Gwn(2) = ComponentAverages_NWP(GWNZZ,b)/awn;
Gwn(3) = ComponentAverages_NWP(GWNXY,b)/awn;
Gwn(4) = ComponentAverages_NWP(GWNXZ,b)/awn;
Gwn(5) = ComponentAverages_NWP(GWNYZ,b)/awn;
}
else Jwn=Kwn=0.0;
trawn = ComponentAverages_NWP(TRAWN,b);
trJwn = ComponentAverages_NWP(TRJWN,b);
if (trawn > 0.0) trJwn /= trawn;
lwns = ComponentAverages_NWP(LWNS,b);
if (lwns != 0.0){
cwns = ComponentAverages_NWP(CWNS,b)/lwns;
vawns(0) = ComponentAverages_NWP(VWNSX,b)/lwns;
vawns(1) = ComponentAverages_NWP(VWNSY,b)/lwns;
vawns(2) = ComponentAverages_NWP(VWNSZ,b)/lwns;
}
else cwns=0.0;
ComponentAverages_NWP(PRS,b) = pn;
ComponentAverages_NWP(VX,b) = van(0);
ComponentAverages_NWP(VY,b) = van(1);
ComponentAverages_NWP(VZ,b) = van(2);
ComponentAverages_NWP(VSQ,b) = vsq;
ComponentAverages_NWP(JWN,b) = Jwn;
ComponentAverages_NWP(KWN,b) = Kwn;
ComponentAverages_NWP(VWNX,b) = vawn(0);
ComponentAverages_NWP(VWNY,b) = vawn(1);
ComponentAverages_NWP(VWNZ,b) = vawn(2);
ComponentAverages_NWP(GWNXX,b) = Gwn(0);
ComponentAverages_NWP(GWNYY,b) = Gwn(1);
ComponentAverages_NWP(GWNZZ,b) = Gwn(2);
ComponentAverages_NWP(GWNXY,b) = Gwn(3);
ComponentAverages_NWP(GWNXZ,b) = Gwn(4);
ComponentAverages_NWP(GWNYZ,b) = Gwn(5);
ComponentAverages_NWP(CWNS,b) = cwns;
ComponentAverages_NWP(VWNSX,b) = vawns(0);
ComponentAverages_NWP(VWNSY,b) = vawns(1);
ComponentAverages_NWP(VWNSZ,b) = vawns(2);
ComponentAverages_NWP(CMX,b) /= Vn;
ComponentAverages_NWP(CMY,b) /= Vn;
ComponentAverages_NWP(CMZ,b) /= Vn;
ComponentAverages_NWP(TRJWN,b) = trJwn;
ComponentAverages_NWP(EULER,b) /= (2*PI);
}
}
if (Dm.rank()==0){
printf("reduce WP averages... \n");
}
// reduce the wetting phase averages
for (int b=0; b<NumberComponents_WP; b++){
MPI_Barrier(Dm.Comm);
// MPI_Allreduce(&ComponentAverages_WP(0,b),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Reduce(&ComponentAverages_WP(0,b),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
for (int idx=0; idx<BLOB_AVG_COUNT; idx++) ComponentAverages_WP(idx,b)=RecvBuffer(idx);
}
for (int b=0; b<NumberComponents_WP; b++){
if (ComponentAverages_WP(VOL,b) > 0.0){
double Vw,pw,awn,ans,Jwn,Kwn,lwns,cwns,vsq;
Vw = ComponentAverages_WP(VOL,b);
awn = ComponentAverages_WP(AWN,b);
ans = ComponentAverages_WP(ANS,b);
vaw(0) = ComponentAverages_WP(VX,b)/Vw;
vaw(1) = ComponentAverages_WP(VY,b)/Vw;
vaw(2) = ComponentAverages_WP(VZ,b)/Vw;
vsq = ComponentAverages_WP(VSQ,b)/Vw;
NULL_USE(ans);
if (ComponentAverages_WP(TRIMVOL,b) > 0.0){
pw = ComponentAverages_WP(PRS,b)/ComponentAverages_WP(TRIMVOL,b);
}
else pw = 0.0;
if (awn != 0.0){
Jwn = ComponentAverages_WP(JWN,b)/awn;
Kwn = ComponentAverages_WP(KWN,b)/awn;
vawn(0) = ComponentAverages_WP(VWNSX,b)/awn;
vawn(1) = ComponentAverages_WP(VWNSY,b)/awn;
vawn(2) = ComponentAverages_WP(VWNSZ,b)/awn;
Gwn(0) = ComponentAverages_WP(GWNXX,b)/awn;
Gwn(1) = ComponentAverages_WP(GWNYY,b)/awn;
Gwn(2) = ComponentAverages_WP(GWNZZ,b)/awn;
Gwn(3) = ComponentAverages_WP(GWNXY,b)/awn;
Gwn(4) = ComponentAverages_WP(GWNXZ,b)/awn;
Gwn(5) = ComponentAverages_WP(GWNYZ,b)/awn;
}
else Jwn=Kwn=0.0;
trawn = ComponentAverages_WP(TRAWN,b);
trJwn = ComponentAverages_WP(TRJWN,b);
if (trawn > 0.0) trJwn /= trawn;
lwns = ComponentAverages_WP(LWNS,b);
if (lwns != 0.0){
cwns = ComponentAverages_WP(CWNS,b)/lwns;
vawns(0) = ComponentAverages_WP(VWNSX,b)/lwns;
vawns(1) = ComponentAverages_WP(VWNSY,b)/lwns;
vawns(2) = ComponentAverages_WP(VWNSZ,b)/lwns;
}
else cwns=0.0;
ComponentAverages_WP(PRS,b) = pw;
ComponentAverages_WP(VX,b) = vaw(0);
ComponentAverages_WP(VY,b) = vaw(1);
ComponentAverages_WP(VZ,b) = vaw(2);
ComponentAverages_WP(VSQ,b) = vsq;
ComponentAverages_WP(JWN,b) = Jwn;
ComponentAverages_WP(KWN,b) = Kwn;
ComponentAverages_WP(VWNX,b) = vawn(0);
ComponentAverages_WP(VWNY,b) = vawn(1);
ComponentAverages_WP(VWNZ,b) = vawn(2);
ComponentAverages_WP(GWNXX,b) = Gwn(0);
ComponentAverages_WP(GWNYY,b) = Gwn(1);
ComponentAverages_WP(GWNZZ,b) = Gwn(2);
ComponentAverages_WP(GWNXY,b) = Gwn(3);
ComponentAverages_WP(GWNXZ,b) = Gwn(4);
ComponentAverages_WP(GWNYZ,b) = Gwn(5);
ComponentAverages_WP(CWNS,b) = cwns;
ComponentAverages_WP(VWNSX,b) = vawns(0);
ComponentAverages_WP(VWNSY,b) = vawns(1);
ComponentAverages_WP(VWNSZ,b) = vawns(2);
ComponentAverages_WP(TRJWN,b) = trJwn;
}
}
}
void TwoPhase::WriteSurfaces(int logcount)
{
int i,j,k;
int ncubes=(Nx-1)*(Ny-1)*(Nz-1);
Point P,A,B,C;
std::shared_ptr<IO::TriList> wn_mesh( new IO::TriList() );
wn_mesh->A.reserve(8*ncubes);
wn_mesh->B.reserve(8*ncubes);
wn_mesh->C.reserve(8*ncubes);
std::shared_ptr<IO::TriList> ns_mesh( new IO::TriList() );
ns_mesh->A.reserve(8*ncubes);
ns_mesh->B.reserve(8*ncubes);
ns_mesh->C.reserve(8*ncubes);
std::shared_ptr<IO::TriList> ws_mesh( new IO::TriList() );
ws_mesh->A.reserve(8*ncubes);
ws_mesh->B.reserve(8*ncubes);
ws_mesh->C.reserve(8*ncubes);
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){ // Get cube from the list
//...........................................................................
// Construct the interfaces and common curve
pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue,
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);
//.......................................................................................
// Write the triangle lists to text file
for (int r=0;r<n_nw_tris;r++){
A = nw_pts(nw_tris(0,r));
B = nw_pts(nw_tris(1,r));
C = nw_pts(nw_tris(2,r));
// compare the trianlge orientation against the color gradient
// Orientation of the triangle
double tri_normal_x = (A.y-B.y)*(B.z-C.z) - (A.z-B.z)*(B.y-C.y);
double tri_normal_y = (A.z-B.z)*(B.x-C.x) - (A.x-B.x)*(B.z-C.z);
double tri_normal_z = (A.x-B.x)*(B.y-C.y) - (A.y-B.y)*(B.x-C.x);
double normal_x = SDn_x(i,j,k);
double normal_y = SDn_y(i,j,k);
double normal_z = SDn_z(i,j,k);
// If the normals don't point in the same direction, flip the orientation of the triangle
// Right hand rule for triangle orientation is used to determine rendering for most software
if (normal_x*tri_normal_x + normal_y*tri_normal_y + normal_z*tri_normal_z < 0.0){
P = A;
A = C;
C = P;
}
// Remap the points
A.x += 1.0*Dm.iproc()*(Nx-2);
A.y += 1.0*Dm.jproc()*(Nx-2);
A.z += 1.0*Dm.kproc()*(Nx-2);
B.x += 1.0*Dm.iproc()*(Nx-2);
B.y += 1.0*Dm.jproc()*(Nx-2);
B.z += 1.0*Dm.kproc()*(Nx-2);
C.x += 1.0*Dm.iproc()*(Nx-2);
C.y += 1.0*Dm.jproc()*(Nx-2);
C.z += 1.0*Dm.kproc()*(Nx-2);
wn_mesh->A.push_back(A);
wn_mesh->B.push_back(B);
wn_mesh->C.push_back(C);
}
for (int r=0;r<n_ws_tris;r++){
A = ws_pts(ws_tris(0,r));
B = ws_pts(ws_tris(1,r));
C = ws_pts(ws_tris(2,r));
// Remap the points
A.x += 1.0*Dm.iproc()*(Nx-2);
A.y += 1.0*Dm.jproc()*(Nx-2);
A.z += 1.0*Dm.kproc()*(Nx-2);
B.x += 1.0*Dm.iproc()*(Nx-2);
B.y += 1.0*Dm.jproc()*(Nx-2);
B.z += 1.0*Dm.kproc()*(Nx-2);
C.x += 1.0*Dm.iproc()*(Nx-2);
C.y += 1.0*Dm.jproc()*(Nx-2);
C.z += 1.0*Dm.kproc()*(Nx-2);
ws_mesh->A.push_back(A);
ws_mesh->B.push_back(B);
ws_mesh->C.push_back(C);
}
for (int r=0;r<n_ns_tris;r++){
A = ns_pts(ns_tris(0,r));
B = ns_pts(ns_tris(1,r));
C = ns_pts(ns_tris(2,r));
// Remap the points
A.x += 1.0*Dm.iproc()*(Nx-2);
A.y += 1.0*Dm.jproc()*(Nx-2);
A.z += 1.0*Dm.kproc()*(Nx-2);
B.x += 1.0*Dm.iproc()*(Nx-2);
B.y += 1.0*Dm.jproc()*(Nx-2);
B.z += 1.0*Dm.kproc()*(Nx-2);
C.x += 1.0*Dm.iproc()*(Nx-2);
C.y += 1.0*Dm.jproc()*(Nx-2);
C.z += 1.0*Dm.kproc()*(Nx-2);
ns_mesh->A.push_back(A);
ns_mesh->B.push_back(B);
ns_mesh->C.push_back(C);
}
}
}
}
std::vector<IO::MeshDataStruct> meshData(4);
meshData[0].meshName = "wn-tris";
meshData[0].mesh = wn_mesh;
meshData[1].meshName = "ws-tris";
meshData[1].mesh = ws_mesh;
meshData[2].meshName = "ns-tris";
meshData[2].mesh = ns_mesh;
IO::writeData( logcount, meshData, Dm.Comm );
}
void TwoPhase::Reduce()
{
int i;
double iVol_global=1.0/Volume;
//...........................................................................
MPI_Barrier(Dm.Comm);
MPI_Allreduce(&nwp_volume,&nwp_volume_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&wp_volume,&wp_volume_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&awn,&awn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&ans,&ans_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&aws,&aws_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&lwns,&lwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&As,&As_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Jwn,&Jwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Kwn,&Kwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&KGwns,&KGwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&KNwns,&KNwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&wwndnw,&wwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&wwnsdnwn,&wwnsdnwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Jwnwwndnw,&Jwnwwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
// Phase averages
MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&paw,&paw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&pan,&pan_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vaw(0),&vaw_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&van(0),&van_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vawn(0),&vawn_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vawns(0),&vawns_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Gwn(0),&Gwn_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Gns(0),&Gns_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Gws(0),&Gws_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&trawn,&trawn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&trJwn,&trJwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&trRwn,&trRwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&euler,&euler_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&An,&An_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Jn,&Jn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Kn,&Kn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Barrier(Dm.Comm);
// Normalize the phase averages
// (density of both components = 1.0)
if (vol_w_global > 0.0){
paw_global = paw_global / vol_w_global;
}
if (wp_volume_global > 0.0){
vaw_global(0) = vaw_global(0) / wp_volume_global;
vaw_global(1) = vaw_global(1) / wp_volume_global;
vaw_global(2) = vaw_global(2) / wp_volume_global;
}
if (vol_n_global > 0.0){
pan_global = pan_global / vol_n_global;
}
if (nwp_volume_global > 0.0){
van_global(0) = van_global(0) / nwp_volume_global;
van_global(1) = van_global(1) / nwp_volume_global;
van_global(2) = van_global(2) / nwp_volume_global;
}
// Normalize surface averages by the interfacial area
if (awn_global > 0.0){
Jwn_global /= awn_global;
Kwn_global /= awn_global;
wwndnw_global /= awn_global;
for (i=0; i<3; i++) vawn_global(i) /= awn_global;
for (i=0; i<6; i++) Gwn_global(i) /= awn_global;
}
if (lwns_global > 0.0){
efawns_global /= lwns_global;
KNwns_global /= lwns_global;
KGwns_global /= lwns_global;
for (i=0; i<3; i++) vawns_global(i) /= lwns_global;
}
if (trawn_global > 0.0){
trJwn_global /= trawn_global;
trRwn_global /= trawn_global;
trRwn_global = 2.0*fabs(trRwn_global);
trJwn_global = fabs(trJwn_global);
}
if (ans_global > 0.0) for (i=0; i<6; i++) Gns_global(i) /= ans_global;
if (aws_global > 0.0) for (i=0; i<6; i++) Gws_global(i) /= aws_global;
euler_global /= (2*PI);
//sat_w = 1.0 - nwp_volume_global*iVol_global/porosity;
sat_w = 1.0 - nwp_volume_global/(nwp_volume_global+wp_volume_global);
// Compute the specific interfacial areas and common line length (dimensionless per unit volume)
awn_global = awn_global*iVol_global;
ans_global = ans_global*iVol_global;
aws_global = aws_global*iVol_global;
dEs = dEs*iVol_global;
lwns_global = lwns_global*iVol_global;
}
void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT)
{
awn_global *= D;
ans_global *= D;
ans_global *= D;
lwns_global *= D*D;
}
void TwoPhase::PrintAll(int timestep)
{
if (Dm.rank()==0){
fprintf(TIMELOG,"%i %.5g ",timestep,dEs); // change in surface energy
fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure
fprintf(TIMELOG,"%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas
fprintf(TIMELOG,"%.5g %.5g ",Jwn_global, Kwn_global); // curvature of wn interface
fprintf(TIMELOG,"%.5g ",lwns_global); // common curve length
fprintf(TIMELOG,"%.5g ",efawns_global); // average contact angle
fprintf(TIMELOG,"%.5g %.5g ",KNwns_global, KGwns_global); // curvature of wn interface
fprintf(TIMELOG,"%.5g %.5g %.5g ",vaw_global(0),vaw_global(1),vaw_global(2)); // average velocity of w phase
fprintf(TIMELOG,"%.5g %.5g %.5g ",van_global(0),van_global(1),van_global(2)); // average velocity of n phase
fprintf(TIMELOG,"%.5g %.5g %.5g ",vawn_global(0),vawn_global(1),vawn_global(2)); // velocity of wn interface
fprintf(TIMELOG,"%.5g %.5g %.5g ",vawns_global(0),vawns_global(1),vawns_global(2)); // velocity of wn interface
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",
Gwn_global(0),Gwn_global(1),Gwn_global(2),Gwn_global(3),Gwn_global(4),Gwn_global(5)); // orientation of wn interface
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",
Gns_global(0),Gns_global(1),Gns_global(2),Gns_global(3),Gns_global(4),Gns_global(5)); // orientation of ns interface
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",
Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface
fprintf(TIMELOG,"%.5g %.5g %.5g ",trawn_global, trJwn_global, trRwn_global); // Trimmed curvature
fprintf(TIMELOG,"%.5g %.5g %.5g ",wwndnw_global, wwnsdnwn_global, Jwnwwndnw_global); // kinematic quantities
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",euler_global, Kn_global, Jn_global, An_global); // minkowski measures
fflush(TIMELOG);
}
else{
sat_w = 1.0 - nwp_volume/(nwp_volume+wp_volume);
fprintf(TIMELOG,"%i ",timestep); // change in surface energy
fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw,pan); // saturation and pressure
fprintf(TIMELOG,"%.5g %.5g %.5g ",awn,ans,aws); // interfacial areas
fprintf(TIMELOG,"%.5g %.5g ",Jwn, Kwn); // curvature of wn interface
fprintf(TIMELOG,"%.5g ",lwns); // common curve length
fprintf(TIMELOG,"%.5g ",efawns); // average contact angle
fprintf(TIMELOG,"%.5g %.5g ",KNwns, KGwns); // curvature of wn interface
fprintf(TIMELOG,"%.5g %.5g %.5g ",vaw(0),vaw(1),vaw(2)); // average velocity of w phase
fprintf(TIMELOG,"%.5g %.5g %.5g ",van(0),van(1),van(2)); // average velocity of n phase
fprintf(TIMELOG,"%.5g %.5g %.5g ",vawn(0),vawn(1),vawn(2)); // velocity of wn interface
fprintf(TIMELOG,"%.5g %.5g %.5g ",vawns(0),vawns(1),vawns(2)); // velocity of wn interface
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",
Gwn(0),Gwn(1),Gwn(2),Gwn(3),Gwn(4),Gwn(5)); // orientation of wn interface
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",
Gns(0),Gns(1),Gns(2),Gns(3),Gns(4),Gns(5)); // orientation of ns interface
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",
Gws(0),Gws(1),Gws(2),Gws(3),Gws(4),Gws(5)); // orientation of ws interface
fprintf(TIMELOG,"%.5g %.5g %.5g ",trawn, trJwn, trRwn); // Trimmed curvature
fprintf(TIMELOG,"%.5g %.5g %.5g ",wwndnw, wwnsdnwn, Jwnwwndnw); // kinematic quantities
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",euler, Kn, Jn, An); // minkowski measures
fflush(TIMELOG);
}
}
void TwoPhase::PrintComponents(int timestep)
{
if (Dm.rank()==0){
printf("PRINT %i COMPONENT AVEREAGES: time = %i \n",(int)ComponentAverages_NWP.size(1),timestep);
for (int b=0; b<NumberComponents_NWP; b++){
//if (ComponentAverages_NWP(TRIMVOL,b) > 0.0){
fprintf(NWPLOG,"%i ",timestep-5);
if ( Label_NWP_map.empty() ) {
// print index
fprintf(NWPLOG,"%i ",b);
} else {
// print final global id
fprintf(NWPLOG,"%i ",Label_NWP_map[b]);
}
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VOL,b));
// fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(TRIMVOL,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(PRS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(AWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(ANS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(JWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(KWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(LWNS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CWNS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNSX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNSY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNSZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VSQ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNXX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNYY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNZZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNXY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNXZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNYZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CMX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CMY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CMZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(TRAWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(TRJWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(INTCURV,b));
fprintf(NWPLOG,"%.5g\n",ComponentAverages_NWP(EULER,b));
// fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(NVERT,b));
// fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(NSIDE,b));
// fprintf(NWPLOG,"%.5g\n",ComponentAverages_NWP(NFACE,b));
//}
}
fflush(NWPLOG);
for (int b=0; b<NumberComponents_WP; b++){
if (ComponentAverages_WP(TRIMVOL,b) > 0.0){
fprintf(WPLOG,"%i ",timestep-5);
fprintf(WPLOG,"%i ",b);
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VOL,b));
// fprintf(WPLOG,"%.5g ",ComponentAverages_WP(TRIMVOL,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(PRS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(AWN,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(AWS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(JWN,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(KWN,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(LWNS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CWNS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNSX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNSY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNSZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VSQ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNXX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNYY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNZZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNXY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNXZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNYZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CMX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CMY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CMZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(TRAWN,b));
fprintf(WPLOG,"%.5g\n",ComponentAverages_WP(TRJWN,b));
}
}
fflush(WPLOG);
}
}
inline int TwoPhase::GetCubeLabel(int i, int j, int k, IntArray &BlobLabel)
{
int label;
label=BlobLabel(i,j,k);
label=max(label,BlobLabel(i+1,j,k));
label=max(label,BlobLabel(i,j+1,k));
label=max(label,BlobLabel(i+1,j+1,k));
label=max(label,BlobLabel(i,j,k+1));
label=max(label,BlobLabel(i+1,j,k+1));
label=max(label,BlobLabel(i,j+1,k+1));
label=max(label,BlobLabel(i+1,j+1,k+1));
return label;
}
void TwoPhase::SortBlobs()
{
//printf("Sorting the blobs based on volume \n");
//printf("-----------------------------------------------\n");
int TempLabel,a,aa,bb,i,j,k,idx;
double TempValue;
//.......................................................................
// Sort NWP components by volume
//.......................................................................
IntArray OldLabel(NumberComponents_NWP);
for (a=0; a<NumberComponents_NWP; a++) OldLabel(a) = a;
// Sort the blob averages based on volume
for (aa=0; aa<NumberComponents_NWP-1; aa++)
{
for ( bb=aa+1; bb<NumberComponents_NWP; bb++)
{
if (ComponentAverages_NWP(VOL,aa) < ComponentAverages_NWP(VOL,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<BLOB_AVG_COUNT; idx++)
{
TempValue = ComponentAverages_NWP(idx,bb);
ComponentAverages_NWP(idx,bb) = ComponentAverages_NWP(idx,aa);
ComponentAverages_NWP(idx,aa) = TempValue;
}
}
}
}
IntArray NewLabel(NumberComponents_NWP);
for (aa=0; aa<NumberComponents_NWP; 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 (Label_NWP(i,j,k) > -1)
{
TempLabel = NewLabel(Label_NWP(i,j,k));
Label_NWP(i,j,k) = TempLabel;
}
}
}
}
//.......................................................................
// Sort WP components by volume
//.......................................................................
OldLabel.resize(NumberComponents_WP);
for (a=0; a<NumberComponents_WP; a++) OldLabel(a) = a;
// Sort the blob averages based on volume
for (aa=0; aa<NumberComponents_WP-1; aa++)
{
for ( bb=aa+1; bb<NumberComponents_WP; bb++)
{
if (ComponentAverages_WP(VOL,aa) < ComponentAverages_WP(VOL,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<BLOB_AVG_COUNT; idx++)
{
TempValue = ComponentAverages_WP(idx,bb);
ComponentAverages_WP(idx,bb) = ComponentAverages_WP(idx,aa);
ComponentAverages_WP(idx,aa) = TempValue;
}
}
}
}
NewLabel.resize(NumberComponents_WP);
for (aa=0; aa<NumberComponents_WP; 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 (Label_WP(i,j,k) > -1)
{
TempLabel = NewLabel(Label_WP(i,j,k));
Label_WP(i,j,k) = TempLabel;
}
}
}
}
//.......................................................................
}
Reference in New Issue View Git Blame Copy Permalink