yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3978
committer: bchareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Wed 2016-11-23 15:35:05 +0100
message:
  merging PhaseCluster code with r3975
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	2016-11-23 10:43:14 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.cpp	2016-11-23 14:35:05 +0000
@@ -17,6 +17,9 @@
 #ifdef TWOPHASEFLOW
 YADE_PLUGIN((TwoPhaseFlowEngineT));
 YADE_PLUGIN((TwoPhaseFlowEngine));
+YADE_PLUGIN((PhaseCluster));
+
+PhaseCluster::~PhaseCluster(){}
 
 void TwoPhaseFlowEngine::initialization()
 {
@@ -618,24 +621,35 @@
 {
     RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
     FiniteCellsIterator cellEnd = tri.finite_cells_end();
+    clusters[0]->reset(); clusters[1]->reset(); 
     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;
+      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);
+	      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;
     }
 }
 
-void TwoPhaseFlowEngine:: updateCellLabel()
+void TwoPhaseFlowEngine::clusterGetFacet(PhaseCluster* cluster, CellHandle cell, int facet) {
+	cell->info().hasInterface = true;
+	cluster->interfacialArea += 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;
+	}
+}
+
+void TwoPhaseFlowEngine::updateCellLabel()
 {
-    int currentLabel = getMaxCellLabel();
+    int currentLabel = getMaxCellLabel();//FIXME: A loop on cells for each new label?? is it serious??
     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) {
-	    clusters.push_back(shared_ptr<PhaseCluster>(new PhaseCluster));
-	    updateSingleCellLabelRecursion(cell,currentLabel+1,clusters.back()->pores);
+	    clusters.push_back(shared_ptr<PhaseCluster> (new PhaseCluster()));
+	    updateSingleCellLabelRecursion(cell,currentLabel+1,clusters.back().get());
             currentLabel++;
         }
     }
@@ -652,10 +666,10 @@
     return maxLabel;
 }
 
-void TwoPhaseFlowEngine:: updateSingleCellLabelRecursion(CellHandle cell, int label, std::vector<CellHandle> *outputVector)
+void TwoPhaseFlowEngine::updateSingleCellLabelRecursion(CellHandle cell, int label, PhaseCluster* cluster)
 {
     cell->info().label=label;
-    if (outputVector) outputVector->push_back(cell);
+    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;
@@ -663,7 +677,10 @@
 //         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,outputVector);
+		updateSingleCellLabelRecursion(nCell,label,cluster);
+	else {
+		clusterGetFacet(cluster,cell,facet);
+	}
     }
 }
 

=== modified file 'pkg/pfv/TwoPhaseFlowEngine.hpp'
--- pkg/pfv/TwoPhaseFlowEngine.hpp	2016-11-23 10:43:14 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.hpp	2016-11-23 14:35:05 +0000
@@ -62,22 +62,23 @@
 class PhaseCluster : public Serializable
 {
 		double totalCellVolume;
+// 		CellHandle entryPoreHandle;
 	public :
 		virtual ~PhaseCluster();
-		vector<TwoPhaseFlowEngineT::CellHandle>* pores;
+		vector<TwoPhaseFlowEngineT::CellHandle> pores;
 		TwoPhaseFlowEngineT::RTriangulation* tri;
+		void reset() {label=entryPore=-1;volume=entryRadius=interfacialArea=0; pores.clear();}
+		
 		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."))
 		((double,volume,0,,"cumulated volume of all pores."))
-		((double,entryPc,0,,"smallest entry capillary pressure."))
-		((int,entryPore,0,,"the pore of the cluster incident to the throat with smallest entry Pc."))
+		((double,entryRadius,0,,"smallest entry capillary pressure."))
+		((int,entryPore,-1,,"the pore of the cluster incident to the throat with smallest entry Pc."))
 		((double,interfacialArea,0,,"interfacial area of the cluster"))
 		,,,
 		)
 };
 REGISTER_SERIALIZABLE(PhaseCluster);
-YADE_PLUGIN((PhaseCluster));
-PhaseCluster::~PhaseCluster(){}
 
 class TwoPhaseFlowEngine : public TwoPhaseFlowEngineT
 {
@@ -126,7 +127,7 @@
 	void checkTrap(double pressure);
 	void updateReservoirLabel();
 	void updateCellLabel();
-	void updateSingleCellLabelRecursion(CellHandle cell, int label, std::vector<CellHandle> *outputVector=NULL);
+	void updateSingleCellLabelRecursion(CellHandle cell, int label, PhaseCluster* cluster);
 	int getMaxCellLabel();
 
 	void invasion2();//without-trap
@@ -153,6 +154,11 @@
 	}
 	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();
 	void saveVtk(const char* folder) {bool initT=solver->noCache; solver->noCache=false; solver->saveVtk(folder); solver->noCache=initT;}
@@ -210,7 +216,7 @@
 	((bool, isImbibitionActivated, false,, "Activates imbibition."))
 
 	
-	,/*TwoPhaseFlowEngineT()*/,
+	,/*clusters.resize(2);*//*TwoPhaseFlowEngineT()*/,
 	,
 	.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")
@@ -239,6 +245,7 @@
 	.def("computeCapillaryForce",&TwoPhaseFlowEngine::computeCapillaryForce,"Compute capillary force. ")
 	.def("saveVtk",&TwoPhaseFlowEngine::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
 	.def("getPotentialPendularSpheresPair",&TwoPhaseFlowEngine::getPotentialPendularSpheresPair,"Get the list of sphere ID pairs of potential pendular liquid bridge.")
+	.def("getClusters",&TwoPhaseFlowEngine::pyClusters/*,(boost::python::arg("folder")="./VTK")*/,"Get the list of clusters.")
 	
 	)
 	DECLARE_LOGGER;

=== modified file 'pkg/pfv/UnsaturatedEngine.cpp'
--- pkg/pfv/UnsaturatedEngine.cpp	2016-11-23 10:43:14 +0000
+++ pkg/pfv/UnsaturatedEngine.cpp	2016-11-23 14:35:05 +0000
@@ -37,27 +37,7 @@
 		double getSpecificInterfacialArea();
 
 		void printSomething();
-
-		boost::python::list getPotentialPendularSpheresPair() {
-			RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
-			boost::python::list bridgeIds;
-			FiniteEdgesIterator ed_it = Tri.finite_edges_begin();
-			for ( ; ed_it!=Tri.finite_edges_end(); ed_it++ ) {
-			  if (detectBridge(ed_it)==true) {
-			    const VertexInfo& vi1=(ed_it->first)->vertex(ed_it->second)->info();
-			    const VertexInfo& vi2=(ed_it->first)->vertex(ed_it->third)->info();
-			    const int& id1 = vi1.id();
-			    const int& id2 = vi2.id();
-			    bridgeIds.append(boost::python::make_tuple(id1,id2));}}
-			    return bridgeIds;}
-  		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;}
-				
 		virtual ~UnsaturatedEngine();
-
 		virtual void action();
 		
 		YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(UnsaturatedEngine,TwoPhaseFlowEngine,"Preliminary version engine of a drainage model for unsaturated soils. Note:Air reservoir is on the top; water reservoir is on the bottom.(deprecated engine, use TwoPhaseFlowEngine instead)",
@@ -74,8 +54,6 @@
 		.def("getWindowsSaturation",&UnsaturatedEngine::getWindowsSaturation,(boost::python::arg("windowsID"),boost::python::arg("isSideBoundaryIncluded")), "get saturation of subdomain with windowsID. If isSideBoundaryIncluded=false (default), the pores of side boundary are excluded in saturation calculating; if isSideBoundaryIncluded=true (only in isInvadeBoundary=true drainage mode), the pores of side boundary are included in saturation calculating.")
 		.def("initializeCellWindowsID",&UnsaturatedEngine::initializeCellWindowsID,"Initialize cell windows index. A temporary function for comparison with experiments, will delete soon")
 		.def("printSomething",&UnsaturatedEngine::printSomething,"print debug.")
-		.def("getPotentialPendularSpheresPair",&UnsaturatedEngine::getPotentialPendularSpheresPair,"Get the list of sphere ID pairs of potential pendular liquid bridge.")
-		.def("getClusters",&UnsaturatedEngine::pyClusters/*,(boost::python::arg("folder")="./VTK")*/,"Get the list of clusters.")
 		)
 		DECLARE_LOGGER;
 };