yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3981
committer: bchareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Fri 2016-11-25 18:09:35 +0100
message:
  update phaseCluster attributes: volume, interfacial area, etc. + various smallfixes in the cluster labelling
modified:
  pkg/pfv/TwoPhaseFlowEngine.cpp
  pkg/pfv/TwoPhaseFlowEngine.hpp


--
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	2016-11-23 14:35:05 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.cpp	2016-11-25 17:09:35 +0000
@@ -41,7 +41,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();
+// 		if(isCellLabelActivated) updateReservoirLabel();
 		solver->noCache = true;
 }
 
@@ -619,12 +619,13 @@
 
 void TwoPhaseFlowEngine:: updateReservoirLabel()
 {
+    clusters[0]->reset(); clusters[0]->label=0;
+    clusters[1]->reset(); clusters[1]->label=1;
     RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
-    FiniteCellsIterator cellEnd = tri.finite_cells_end();
-    clusters[0]->reset(); clusters[1]->reset(); 
+    FiniteCellsIterator cellEnd = tri.finite_cells_end();     
     for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
-      if (cell->info().isNWRes) {cell->info().label=0; clusters[0]->pores.push_back(cell);}
-      else if (cell->info().isWRes) {cell->info().label=1; clusters[1]->pores.push_back(cell);
+      if (cell->info().isNWRes) clusterGetPore(clusters[0].get(),cell);
+      else if (cell->info().isWRes) { clusterGetPore(clusters[1].get(),cell);
 	      for (int facet = 0; facet < 4; facet ++) if (!cell->neighbor(facet)->info().isWRes) clusterGetFacet(clusters[1].get(),cell,facet);}
       else if (cell->info().label>1) continue;
       else cell->info().label=-1;
@@ -633,43 +634,52 @@
 
 void TwoPhaseFlowEngine::clusterGetFacet(PhaseCluster* cluster, CellHandle cell, int facet) {
 	cell->info().hasInterface = true;
-	cluster->interfacialArea += cell->info().poreThroatRadius[facet];//FIXME: define area correctly
+	cluster->interfacialArea += std::abs(cell->info().poreThroatRadius[facet]);//FIXME: define area correctly
 	if (cluster->entryRadius < cell->info().poreThroatRadius[facet]){
 		cluster->entryRadius = cell->info().poreThroatRadius[facet];
-		cluster->entryPore = cell->info().index;
-	}
-}
+		cluster->entryPore = cell->info().id;}
+}
+
+void TwoPhaseFlowEngine::clusterGetPore(PhaseCluster* cluster, CellHandle cell) {
+	cell->info().label=cluster->label;
+	cluster->volume+=cell->info().poreBodyVolume;
+	cluster->pores.push_back(cell);
+}
+// 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::updateCellLabel()
 {
-    int currentLabel = getMaxCellLabel();//FIXME: A loop on cells for each new label?? is it serious??
+//     int currentLabel = getMaxCellLabel();//FIXME: A loop on cells for each new label?? is it serious??
     updateReservoirLabel();
+    int currentLabel = clusters.size();
     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) {
-	    clusters.push_back(shared_ptr<PhaseCluster> (new PhaseCluster()));
-	    updateSingleCellLabelRecursion(cell,currentLabel+1,clusters.back().get());
+	    shared_ptr<PhaseCluster> clst (new PhaseCluster());
+	    clst->label=currentLabel;
+	    clusters.push_back(clst);
+	    updateSingleCellLabelRecursion(cell,clusters.back().get());
             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, PhaseCluster* cluster)
-{
-    cell->info().label=label;
-    cluster->pores.push_back(cell);
+void TwoPhaseFlowEngine::updateSingleCellLabelRecursion(CellHandle cell, PhaseCluster* cluster)
+{
+	clusterGetPore(cluster,cell);
+//     cell->info().label=label;
+//     cluster->volume+=cell->info().
+//     cluster->pores.push_back(cell);
     for (int facet = 0; facet < 4; facet ++) {
         CellHandle nCell = cell->neighbor(facet);
         if (solver->T[solver->currentTes].Triangulation().is_infinite(nCell)) continue;
@@ -677,13 +687,16 @@
 //         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,cluster);
-	else {
-		clusterGetFacet(cluster,cell,facet);
-	}
+		updateSingleCellLabelRecursion(nCell,cluster);
+	else if (nCell->info().isNWRes) clusterGetFacet(cluster,cell,facet);
     }
 }
 
+boost::python::list TwoPhaseFlowEngine::pyClusters() {
+	boost::python::list ret;
+	for(vector<shared_ptr<PhaseCluster> >::iterator it=clusters.begin(); it!=clusters.end(); ++it) ret.append(*it);
+	return ret;}
+
 void TwoPhaseFlowEngine::updatePressure()
 {
     boundaryConditions(*solver);

=== modified file 'pkg/pfv/TwoPhaseFlowEngine.hpp'
--- pkg/pfv/TwoPhaseFlowEngine.hpp	2016-11-23 14:35:05 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.hpp	2016-11-25 17:09:35 +0000
@@ -68,6 +68,9 @@
 		vector<TwoPhaseFlowEngineT::CellHandle> pores;
 		TwoPhaseFlowEngineT::RTriangulation* tri;
 		void reset() {label=entryPore=-1;volume=entryRadius=interfacialArea=0; pores.clear();}
+		vector<int> getPores() { vector<int> res;
+			for (vector<TwoPhaseFlowEngineT::CellHandle>::iterator it =  pores.begin(); it!=pores.end(); it++) res.push_back((*it)->info().id);
+			return res;}
 		
 		YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(PhaseCluster,Serializable,"Preliminary.",
 		((int,label,-1,,"Unique label of this cluster, should be reflected in pores of this cluster."))
@@ -76,6 +79,7 @@
 		((int,entryPore,-1,,"the pore of the cluster incident to the throat with smallest entry Pc."))
 		((double,interfacialArea,0,,"interfacial area of the cluster"))
 		,,,
+		.def("getPores",&PhaseCluster::getPores,"get the list of pores by index")
 		)
 };
 REGISTER_SERIALIZABLE(PhaseCluster);
@@ -127,8 +131,11 @@
 	void checkTrap(double pressure);
 	void updateReservoirLabel();
 	void updateCellLabel();
-	void updateSingleCellLabelRecursion(CellHandle cell, int label, PhaseCluster* cluster);
-	int getMaxCellLabel();
+	void updateSingleCellLabelRecursion(CellHandle cell, PhaseCluster* cluster);
+	void clusterGetFacet(PhaseCluster* cluster, CellHandle cell, int facet);//update cluster inetrfacial area and max entry radius wrt to a facet
+	void clusterGetPore(PhaseCluster* cluster, CellHandle cell);//add pore to cluster, updating flags and cluster volume
+	boost::python::list pyClusters();
+// 	int getMaxCellLabel();
 
 	void invasion2();//without-trap
 	void updateReservoirs2();
@@ -154,10 +161,7 @@
 	}
 	bool detectBridge(RTriangulation::Finite_edges_iterator& edge);
 	
-	boost::python::list pyClusters() { boost::python::list ret;
-		for(vector<shared_ptr<PhaseCluster> >::iterator it=clusters.begin(); it!=clusters.end(); ++it) ret.append(*it);
-		return ret;}
-	void clusterGetFacet(PhaseCluster* cluster, CellHandle cell, int facet);//update cluster inetrfacial area and max entry radius wrt to a facet
+
 		
 	//post-processing
 	void savePoreNetwork();
@@ -216,7 +220,8 @@
 	((bool, isImbibitionActivated, false,, "Activates imbibition."))
 
 	
-	,/*clusters.resize(2);*//*TwoPhaseFlowEngineT()*/,
+	,/*TwoPhaseFlowEngineT()*/,
+	clusters.resize(2); clusters[0]=shared_ptr<PhaseCluster>(new PhaseCluster); clusters[1]=shared_ptr<PhaseCluster>(new PhaseCluster);
 	,
 	.def("getCellIsFictious",&TwoPhaseFlowEngine::cellIsFictious,"Check the connection between pore and boundary. If true, pore throat connects the boundary.")
 	.def("setCellIsNWRes",&TwoPhaseFlowEngine::setCellIsNWRes,"set status whether 'wetting reservoir' state")