refactor DCEL

analysis/Minkowski.cpp

@@ -15,33 +15,12 @@
 
 // Constructor
 Minkowski::Minkowski(std::shared_ptr <Domain> dm):
-	n_obj_pts(0), n_obj_tris(0), kstart(0), kfinish(0), isovalue(0), Volume(0),
+	kstart(0), kfinish(0), isovalue(0), Volume(0),
-    LOGFILE(NULL), Dm(dm), vol_n(0), vol_n_global(0)
+    LOGFILE(NULL), Dm(dm), Vi(0), Vi_global(0)
 {
 	Nx=dm->Nx; Ny=dm->Ny; Nz=dm->Nz;
 	Volume=double((Nx-2)*(Ny-2)*(Nz-2))*double(Dm->nprocx()*Dm->nprocy()*Dm->nprocz());
 	
-	TempID = new char[Nx*Ny*Nz];
-
-	// Global arrays
-	PhaseID.resize(Nx,Ny,Nz);       PhaseID.fill(0);
-	SDn.resize(Nx,Ny,Nz);           SDn.fill(0);
-	MeanCurvature.resize(Nx,Ny,Nz); MeanCurvature.fill(0);
-	GaussCurvature.resize(Nx,Ny,Nz); GaussCurvature.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);
-	//.........................................
-	// Allocate cube storage space
-	CubeValues.resize(2,2,2);
-	obj_tris.resize(3,20);
-	obj_pts=DTMutableList<Point>(20);
-	tmp=DTMutableList<Point>(20);
-	//.........................................
-	Values.resize(20);
-	//DistanceValues.resize(20);
-	NormalVector.resize(60);
-	
 	if (Dm->rank()==0){
 		LOGFILE = fopen("minkowski.csv","a+");
 		if (fseek(LOGFILE,0,SEEK_SET) == fseek(LOGFILE,0,SEEK_CUR))
@@ -60,117 +39,20 @@ Minkowski::~Minkowski()
     if ( LOGFILE!=NULL ) { fclose(LOGFILE); }
 }
 
-
-void Minkowski::Initialize()
-{
-	isovalue=0.0;
-	vol_n = euler = Jn = An = Kn = 0.0;
-}
-
 double Minkowski::V(){
-	return vol_n_global;
+	return Vi_global;
 }
 
 double Minkowski::A(){
-	return An_global;
+	return Ai_global;
 }
 
 double Minkowski::J(){
-	return Jn_global;
+	return Ji_global;
 }
 
 double Minkowski::X(){
-	return euler_global;
+	return Xi_global;
-}
-
-
-void Minkowski::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);
-	//...........................................................................
-	//...........................................................................
-	// Compute the mesh curvature of the phase indicator field
-	pmmc_MeshCurvature(SDn, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
-	//...........................................................................
-	//...........................................................................
-	Dm->CommunicateMeshHalo(MeanCurvature);
-	//...........................................................................
-	Dm->CommunicateMeshHalo(GaussCurvature);
-	//...........................................................................
-	// 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 {
-					// non-wetting phase
-					PhaseID(i,j,k) = 1;
-				}
-			}
-		}
-	}
-
-}
-void Minkowski::ComputeLocal()
-{
-	int i,j,k,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;
-
-	vol_n = euler = Jn = An = Kn = 0.0;
-	for (k=kmin; k<kmax; k++){
-		for (j=1; j<Ny-1; j++){
-			for (i=1; i<Nx-1; i++){
-				//...........................................................................
-				n_obj_pts=0;
-				n_obj_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;
-					// 1-D index for this cube corner
-					if ( SDn(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < 0 ){
-						vol_n += 0.125;
-					}
-				}
-
-				n_obj_pts=n_obj_tris=0;
-				// Compute the non-wetting phase surface and associated area
-				An += geomavg_MarchingCubes(SDn,isovalue,i,j,k,obj_pts,n_obj_pts,obj_tris,n_obj_tris);
-
-				Jn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,obj_pts,obj_tris,Values,
-										i,j,k,n_obj_pts,n_obj_tris);
-				// Compute Euler characteristic from integral of gaussian curvature
-				Kn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,obj_pts,obj_tris,Values,
-						i,j,k,n_obj_pts,n_obj_tris);
-
-				euler += geomavg_EulerCharacteristic(obj_pts,obj_tris,n_obj_pts,n_obj_tris,i,j,k);
-
-			}
-		}
-	}
-	
 }
 
 void Minkowski::ComputeScalar(const DoubleArray Field, const double isovalue)
@@ -243,64 +125,10 @@ void Minkowski::ComputeScalar(const DoubleArray Field, const double isovalue)
 
 }
 
-/*
-void Minkowski::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,isovalue,Label_NWP,Dm->Comm);
-	
-}
-*/
-void Minkowski::Reduce()
-{
-	//...........................................................................
-	MPI_Barrier(Dm->Comm);
-	// Phase averages
-	MPI_Allreduce(&vol_n,&vol_n_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_Barrier(Dm->Comm);
-
-	// normalize to per unit volume basis
-	vol_n_global /=Volume;
-	An_global /=Volume;
-	Jn_global /=Volume;
-	euler_global /=Volume;
-	euler_global /= (2*PI);
-
-}
-
-void Minkowski::NonDimensionalize(double D)
-{
-	An_global *= D;
-	Jn_global *= D*D;
-	euler_global *= D*D*D;
-}
-
 void Minkowski::PrintAll()
 {
 	if (Dm->rank()==0){
-		fprintf(LOGFILE,"%.5g %.5g %.5g %.5g\n",vol_n_global, An_global, Jn_global, euler_global);			// minkowski measures
+		fprintf(LOGFILE,"%.5g %.5g %.5g %.5g\n",Vi_global, Ai_global, Ji_global, Xi_global);			// minkowski measures
 		fflush(LOGFILE);
 	}
 }

analysis/Minkowski.h

@@ -4,7 +4,7 @@
 
 #include <vector>
 
-#include "analysis/decl.h"
+#include "analysis/dcel.h"
 #include "common/Domain.h"
 #include "common/Communication.h"
 #include "analysis/analysis.h"
@@ -18,31 +18,11 @@
 
 
 class Minkowski{
-
 	//...........................................................................
-	int n_obj_pts;
-	int n_obj_tris;
-	//...........................................................................
-	int nc;
 	int kstart,kfinish;
 
 	double isovalue;
 	double Volume;
-	// initialize lists for vertices for surfaces, common line
-	DTMutableList<Point> obj_pts;
-	DTMutableList<Point> tmp;
-
-	// initialize triangle lists for surfaces
-	IntArray obj_tris;
-
-	// Temporary storage arrays
-	DoubleArray CubeValues;
-	DoubleArray Values;
-	DoubleArray NormalVector;
-
-	DoubleArray RecvBuffer;
-
-	char *TempID;
 
 	// CSV / text file where time history of averages is saved
 	FILE *LOGFILE;
@@ -54,24 +34,12 @@ public:
 	// Averaging variables
 	//...........................................................................
 	// local averages (to each MPI process)
-	double vol_n;						// volumes the exclude the interfacial region
-	// Global averages (all processes)
-	double vol_n_global;			// volumes the exclude the interfacial region
-	double euler,Kn,Jn,An;
-	double euler_global,Kn_global,Jn_global,An_global;
 	double Ai,Ji,Xi,Vi;
+	// Global averages (all processes)
 	double Ai_global,Ji_global,Xi_global,Vi_global;
 
 	//...........................................................................
 	int Nx,Ny,Nz;
-	IntArray PhaseID;	// Phase ID array (solid=0, non-wetting=1, wetting=2)
-	DoubleArray SDn;
-	DoubleArray MeanCurvature;
-	DoubleArray GaussCurvature;
-	DoubleArray SDn_x;		// Gradient of the signed distance
-	DoubleArray SDn_y;
-	DoubleArray SDn_z;
-
 	double V();
 	double A();
 	double J();
@@ -80,14 +48,9 @@ public:
 	//...........................................................................
 	Minkowski(std::shared_ptr <Domain> Dm);
 	~Minkowski();
-	void Initialize();
-	void UpdateMeshValues();
-	void ComputeLocal();
-	void Reduce();
-	void NonDimensionalize(double D);
-	void PrintAll();
-	int GetCubeLabel(int i, int j, int k, IntArray &BlobLabel);
 	void ComputeScalar(const DoubleArray Field, const double isovalue);
+	void PrintAll();
+
 };
 
 #endif

tests/CMakeLists.txt

@@ -54,8 +54,8 @@ ADD_LBPM_TEST( TestForceMoments  ../example/Piston/input.db)
 ADD_LBPM_TEST( TestForceD3Q19 )
 ADD_LBPM_TEST( TestMomentsD3Q19 )
 ADD_LBPM_TEST( TestInterfaceSpeed  ../example/Bubble/input.db)
-ADD_LBPM_TEST( TestSphereCurvature )
+ADD_LBPM_TEST( Test_dcel_Minkowski )
-ADD_LBPM_TEST( TestTriNormal )
+ADD_LBPM_TEST( Test_dcel_TriNormal )
 #ADD_LBPM_TEST_1_2_4( TestTwoPhase )
 ADD_LBPM_TEST_1_2_4( TestBlobIdentify )
 #ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 )