yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3959
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Wed 2016-11-02 11:40:55 +0100
message:
  Add cell label for phase cluster. NW-res label 0, W-res label 1, disconnected w-clusters start from 2. label can be saved in savePhaseVtk.
modified:
  pkg/pfv/TwoPhaseFlowEngine.cpp
  pkg/pfv/TwoPhaseFlowEngine.hpp
  pkg/pfv/UnsaturatedEngine.cpp


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'pkg/pfv/TwoPhaseFlowEngine.cpp'
--- pkg/pfv/TwoPhaseFlowEngine.cpp	2015-09-24 07:25:27 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.cpp	2016-11-02 10:40:55 +0000
@@ -39,6 +39,7 @@
 		computePoreBodyVolume();//save capillary volume of all cells, for fast calculating saturation
 		computeSolidLine();//save cell->info().solidLine[j][y]
 		initializeReservoirs();//initial pressure, reservoir flags and local pore saturation
+		if(isCellLabelActivated) updateReservoirLabel();
 		solver->noCache = true;
 }
 
@@ -410,6 +411,13 @@
 		if (isDrawable){vtkfile.write_data(cell->info().saturation);}
 	}
 	vtkfile.end_data();
+
+	vtkfile.begin_data("Label",CELL_DATA,SCALARS,FLOAT);
+	for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
+		bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal()  && cell->vertex(3)->info().isReal();
+		if (isDrawable){vtkfile.write_data(cell->info().label);}
+	}
+	vtkfile.end_data();
 }
 void TwoPhaseFlowEngine::computePoreThroatRadiusTrickyMethod1()
 {
@@ -654,5 +662,56 @@
   cell->info().p() = Pg - Pc; */
 }
 
+void TwoPhaseFlowEngine:: updateReservoirLabel()
+{
+    RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+    FiniteCellsIterator cellEnd = tri.finite_cells_end();
+    for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+      if (cell->info().isNWRes) cell->info().label=0;
+      else if (cell->info().isWRes) cell->info().label=1;
+      else if (cell->info().label>1) continue;
+      else cell->info().label=-1;
+    }
+}
+
+void TwoPhaseFlowEngine:: updateCellLabel()
+{
+    int currentLabel = getMaxCellLabel();
+    updateReservoirLabel();
+    RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+    FiniteCellsIterator cellEnd = tri.finite_cells_end();
+    for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+        if (cell->info().label==-1) {
+            updateSingleCellLabelRecursion(cell,currentLabel+1);
+            currentLabel++;
+        }
+    }
+}
+
+int TwoPhaseFlowEngine:: getMaxCellLabel()
+{
+    int maxLabel=-1;
+    RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+    FiniteCellsIterator cellEnd = tri.finite_cells_end();
+    for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+        if (cell->info().label>maxLabel) maxLabel=cell->info().label;
+    }
+    return maxLabel;
+}
+
+void TwoPhaseFlowEngine:: updateSingleCellLabelRecursion(CellHandle cell, int label)
+{
+    cell->info().label=label;
+    for (int facet = 0; facet < 4; facet ++) {
+        CellHandle nCell = cell->neighbor(facet);
+        if (solver->T[solver->currentTes].Triangulation().is_infinite(nCell)) continue;
+//         if (nCell->info().Pcondition) continue;
+//         if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+        //TODO:the following condition may relax to relate to nCell->info().hasInterface
+        if ( (nCell->info().saturation==cell->info().saturation) && (nCell->info().label!=cell->info().label) )
+            updateSingleCellLabelRecursion(nCell,label);
+    }
+}
+
 #endif //TwoPhaseFLOW
  

=== modified file 'pkg/pfv/TwoPhaseFlowEngine.hpp'
--- pkg/pfv/TwoPhaseFlowEngine.hpp	2016-10-24 11:41:05 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.hpp	2016-11-02 10:40:55 +0000
@@ -34,6 +34,7 @@
 	int windowsID;//a temp cell info for experiment comparison(used by chao)
 	double solidLine [4][4];//the length of intersecting line between sphere and facet. [i][j] is for facet "i" and sphere (facetVertices)"[i][j]". Last component [i][3] for 1/sumLines in the facet "i" (used by chao).
 	
+	int label;//for marking disconnected clusters. initally all set to -1; first update -> connect to NW-res: 0; connect to W-res: 1; then label disconnected W-clusters by 2,3,4...
 	TwoPhaseCellInfo (void)
 	{
 		isWRes = true; isNWRes = false; isTrapW = false; isTrapNW = false;
@@ -45,6 +46,7 @@
 		poreBodyVolume = 0;
 		windowsID = 0;
 		for (int k=0; k<4;k++) for (int l=0; l<4;l++) solidLine[k][l]=0;
+		label=-1;
 	}
 	
 };
@@ -91,7 +93,11 @@
 	void computePoreThroatRadiusMethod3();
 	void savePoreNetwork();
 	
-	
+	void updateReservoirLabel();
+	void updateCellLabel();
+	void updateSingleCellLabelRecursion(CellHandle cell, int label);
+	int getMaxCellLabel();
+
 	
 	boost::python::list cellporeThroatRadius(unsigned int id){ // Temporary function to allow for simulations in Python, can be easily accessed in c++
 	  boost::python::list ids;
@@ -124,6 +130,7 @@
 	CELL_SCALAR_GETTER(Real,.poreBodyRadius,cellInSphereRadius) //Temporary function to allow for simulations in Python	
 	CELL_SCALAR_GETTER(Real,.poreBodyVolume,cellVoidVolume) //Temporary function to allow for simulations in Python	
 	CELL_SCALAR_SETTER(bool,.hasInterface,setCellHasInterface) //Temporary function to allow for simulations in Python
+	CELL_SCALAR_GETTER(int,.label,cellLabel)
 
 	YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(TwoPhaseFlowEngine,TwoPhaseFlowEngineT,"documentation here",
 	((double,surfaceTension,0.0728,,"Water Surface Tension in contact with air at 20 Degrees Celsius is: 0.0728(N/m)"))
@@ -135,6 +142,7 @@
 	((double,dtDynTPF,0.0,,"Parameter which stores the smallest time step, based on the residence time"))
 	((int,entryPressureMethod,1,,"integer to define the method used to determine the pore throat radii and the according entry pressures. 1)radius of entry pore throat based on MS-P method; 2) radius of the inscribed circle; 3) radius of the circle with equivalent surface area of the pore throat."))
 	((double,partiallySaturatedPores,false,,"Include partially saturated pores or not?"))
+	((bool, isCellLabelActivated, false,, "Activate cell labels for marking disconnected wetting clusters. NW-reservoir label 0; W-reservoir label 1; disconnected W-clusters label from 2. "))
 
 	,/*TwoPhaseFlowEngineT()*/,
 	,
@@ -157,6 +165,7 @@
 	.def("getCellVoidVolume",&TwoPhaseFlowEngine::cellVoidVolume,"get the volume of pore space in each pore unit")
 	.def("setCellHasInterface",&TwoPhaseFlowEngine::setCellHasInterface,"change wheter a cell has a NW-W interface")
 	.def("savePoreNetwork",&TwoPhaseFlowEngine::savePoreNetwork,"Extract the pore network of the granular material")
+	.def("getCellLabel",&TwoPhaseFlowEngine::cellLabel,"get cell label. 0 for NW-reservoir; 1 for W-reservoir; others for disconnected W-clusters.")
 	
 	)
 	DECLARE_LOGGER;

=== modified file 'pkg/pfv/UnsaturatedEngine.cpp'
--- pkg/pfv/UnsaturatedEngine.cpp	2016-10-24 11:41:05 +0000
+++ pkg/pfv/UnsaturatedEngine.cpp	2016-11-02 10:40:55 +0000
@@ -281,6 +281,9 @@
 	if ( cell->info().isTrapNW) {cell->info().p()=bndCondValue[2]+cell->info().trapCapP;}
    }
     if(solver->debugOut) {cout<<"----invasion1.update trapped W-phase/NW-phase Pressure----"<<endl;}
+
+    if(isCellLabelActivated) updateCellLabel();
+    if(solver->debugOut) {cout<<"----update cell labels----"<<endl;}
 }
 
 ///search trapped W-phase or NW-phase, define trapCapP=Pn-Pw. assign isTrapW/isTrapNW info.