From 32121873438497ed272e43ce06a38da491cca995 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Sun, 26 Jul 2015 08:30:56 -0400
Subject: [PATCH] common/TwoPhase.h: when averaging a system with external BC,
 do not compute the contribution from the inlet. tests/ComponenLabel.cpp: use
 the .LBM files instead of writing raw binary so that things can be visualized
 in VisIT

---
 common/TwoPhase.h        | 17 ++++++++++++++---
 tests/ComponentLabel.cpp | 38 +++++++++++++++++++++++++++++++++++++-
 2 files changed, 51 insertions(+), 4 deletions(-)

diff --git a/common/TwoPhase.h b/common/TwoPhase.h
index 1c0f87b3..04a8560a 100644
--- a/common/TwoPhase.h
+++ b/common/TwoPhase.h
@@ -452,10 +452,15 @@ void TwoPhase::UpdateMeshValues(){
 
 }
 void TwoPhase::ComputeLocal(){
-	int i,j,k,n;
+	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}};
 
-	for (k=1; k<Nz-1; k++){
+	// 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++){
 				//...........................................................................
@@ -558,6 +563,7 @@ void TwoPhase::ComputeLocal(){
 
 void TwoPhase::ComponentAverages(){
     int i,j,k,n;
+    int kmin,kmax;
 	int LabelWP,LabelNWP;
 	double TempLocal;
 
@@ -577,7 +583,12 @@ void TwoPhase::ComponentAverages(){
 		printf("Number of non-wetting phase components is %i \n",NumberComponents_NWP);
 	}
 
-	for (k=1; k<Nz-1; k++){
+	// 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++){
 
diff --git a/tests/ComponentLabel.cpp b/tests/ComponentLabel.cpp
index 04224156..0bbccdfe 100644
--- a/tests/ComponentLabel.cpp
+++ b/tests/ComponentLabel.cpp
@@ -434,6 +434,42 @@ int main(int argc, char **argv)
 	Averages.SortBlobs();
 	Averages.PrintComponents(timestep);
 
+    // Create the MeshDataStruct
+    fillHalo<double> fillData(Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
+    std::vector<IO::MeshDataStruct> meshData(1);
+    meshData[0].meshName = "domain";
+    meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
+    std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
+    std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
+    std::shared_ptr<IO::Variable> LabelWPVar( new IO::Variable() );
+    std::shared_ptr<IO::Variable> LabelNWPVar( new IO::Variable() );
+    PhaseVar->name = "phase";
+    PhaseVar->type = IO::VolumeVariable;
+    PhaseVar->dim = 1;
+    PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(PhaseVar);
+    SignDistVar->name = "SignDist";
+    SignDistVar->type = IO::VolumeVariable;
+    SignDistVar->dim = 1;
+    SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(SignDistVar);
+    LabelNWPVar->name = "LabelNWP";
+    LabelNWPVar->type = IO::VolumeVariable;
+    LabelNWPVar->dim = 1;
+    LabelNWPVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(LabelNWPVar);
+    LabelWPVar->name = "LabelWP";
+    LabelWPVar->type = IO::VolumeVariable;
+    LabelWPVar->dim = 1;
+    LabelWPVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(LabelWPVar);
+
+    fillData.copy(Averages.SDn,PhaseVar->data);
+    fillData.copy(Averages.SDs,SignDistVar->data);
+    fillData.copy(Averages.Label_WP,LabelWPVar->data);
+    fillData.copy(Averages.Label_NWP,LabelNWPVar->data);
+    IO::writeData( 0, meshData, 2 );
+/*
 	FILE *NWP_FILE;
 	NWP_FILE = fopen("NWP.dat","wb");
 	fwrite(Averages.Label_NWP.get(),4,Nx*Ny*Nz,NWP_FILE);
@@ -448,7 +484,7 @@ int main(int argc, char **argv)
 	DISTANCE = fopen("SignDist.dat","wb");
 	fwrite(Averages.SDs.get(),8,Nx*Ny*Nz,DISTANCE);
 	fclose(DISTANCE);
-	
+	*/
 	// ****************************************************
 	MPI_Barrier(MPI_COMM_WORLD);
 	MPI_Finalize();