← Back to team overview

yade-dev team mailing list archive

[Branch ~yade-pkg/yade/git-trunk] Rev 3415: -add temp function for generating sample windows, calculating responding saturation...

 

------------------------------------------------------------
revno: 3415
committer: cyuan <chaoyuan2012@xxxxxxxxx>
timestamp: Fri 2014-04-04 19:12:06 +0200
message:
  -add temp function for generating sample windows, calculating responding saturation...
modified:
  lib/triangulation/def_types.h
  pkg/dem/UnsaturatedEngine.cpp
  pkg/dem/UnsaturatedEngine.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 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h	2014-02-17 18:00:12 +0000
+++ lib/triangulation/def_types.h	2014-04-04 17:12:06 +0000
@@ -230,12 +230,14 @@
 	std::vector<double> poreRadius;//pore throat radius for drainage
 	double solidLine [4][4];//the length of intersecting line between sphere and facet. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_Lines in the facet.
 	double trapCapP;//for calculating the pressure of trapped phase, cell.Pressure_trapped = Pressure_air - trapCapP.
+	int windowsID;//a temp cell info for experiment comparison
 	UnsatCellInfo (void)
 	{
 		poreRadius.resize(4, 0);
 		isWaterReservoir = false; isAirReservoir = false; capillaryCellVolume = 0;	  
 		for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidLine[k][l]=0;
 		trapCapP = 0;
+		windowsID = 0;
 	}
 };
 

=== modified file 'pkg/dem/UnsaturatedEngine.cpp'
--- pkg/dem/UnsaturatedEngine.cpp	2014-03-29 23:32:02 +0000
+++ pkg/dem/UnsaturatedEngine.cpp	2014-04-04 17:12:06 +0000
@@ -111,10 +111,10 @@
 Real UnsaturatedEngine::getSaturation(Solver& flow )
 {
     if (isPhaseTrapped) {
-        return getSaturation1(solver);
+        return getSaturation1(flow);
     }
     else {
-        return getSaturation2(solver);
+        return getSaturation2(flow);
     }
 }
 
@@ -1158,7 +1158,7 @@
 template<class Solver>
 void UnsaturatedEngine::clearImposedPressure ( Solver& flow ) { flow->imposedP.clear(); flow->IPCells.clear();}
 
-//----tempt function for Vahid Joekar-Niasar's data----
+//----temp function for Vahid Joekar-Niasar's data----
 template<class Cellhandle >
 double UnsaturatedEngine::getRadiusMin(Cellhandle cell, int j)
 {
@@ -1379,7 +1379,92 @@
     file.close();  
 }
 
-//----------end tempt function for Vahid Joekar-Niasar's data (clear later)---------------------
+//----------end temp function for Vahid Joekar-Niasar's data (clear later)---------------------
+//----------temp functions for comparison with experiment-----------------------
+template <class Solver>
+void UnsaturatedEngine::initializeCellWindowsID(Solver&flow)
+{
+    RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
+    Finite_cells_iterator cell_end = tri.finite_cells_end();
+    for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+        for (int i=1; i<(windowsNo+1); i++) {
+            if ( (cell->info()[1]>(flow->y_min+(i-1)*(flow->y_max-flow->y_min)/windowsNo) ) && (cell->info()[1] < (flow->y_min+i*(flow->y_max-flow->y_min)/windowsNo)) )
+            {cell->info().windowsID=i; break;}
+        }
+    }
+}
+template<class Solver>
+Real UnsaturatedEngine::getWindowsSaturation(Solver&flow, int windowsID)
+{
+    RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
+    Finite_cells_iterator cell_end = tri.finite_cells_end();
+    for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+      if (cell->info().windowsID==0) {cerr<<"Please initialize windowsID"<<endl;break;}
+    }
+    if (isPhaseTrapped) {
+        return getWindowsSaturation1(flow,windowsID);
+    }
+    else {
+        return getWindowsSaturation2(flow,windowsID);
+    }
+}
+template<class Solver>
+Real UnsaturatedEngine::getWindowsSaturation1(Solver&flow, int i)
+{
+    updatePressureReservoir(flow);
+    RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
+    Real capillary_volume = 0.0; //total capillary volume
+    Real air_volume = 0.0; 	//air volume
+    Finite_cells_iterator cell_end = tri.finite_cells_end();
+
+    if (invadeBoundary==true) {
+        for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+            if (tri.is_infinite(cell)) continue;
+            if (cell->info().Pcondition) continue;//NOTE:reservoirs cells should not be included in saturation
+	    if (cell->info().windowsID != i) continue;
+            capillary_volume = capillary_volume + cell->info().capillaryCellVolume;
+            if (cell->info().isAirReservoir==true) {
+                air_volume = air_volume + cell->info().capillaryCellVolume;}}}
+    else {
+        for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+            if (tri.is_infinite(cell)) continue;
+            if (cell->info().Pcondition) continue;
+            if (cell->info().isFictious) continue;
+	    if (cell->info().windowsID != i) continue;
+            capillary_volume = capillary_volume + cell->info().capillaryCellVolume;
+            if (cell->info().isAirReservoir==true) {
+                air_volume = air_volume + cell->info().capillaryCellVolume;}}}
+    Real saturation = 1 - air_volume/capillary_volume;
+    return saturation;
+}
+template<class Solver>
+Real UnsaturatedEngine::getWindowsSaturation2(Solver& flow,int i)
+{
+    RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
+    Real capillary_volume = 0.0;
+    Real water_volume = 0.0;
+    Finite_cells_iterator cell_end = tri.finite_cells_end();
+    if (invadeBoundary==true) {
+        for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+            if (tri.is_infinite(cell)) continue;
+            if (cell->info().Pcondition) continue;
+	    if (cell->info().windowsID != i) continue;
+            capillary_volume = capillary_volume + cell->info().capillaryCellVolume;
+            if (cell->info().p()==0) {
+                water_volume = water_volume + cell->info().capillaryCellVolume;}}}
+    else {
+        for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+            if (tri.is_infinite(cell)) continue;
+            if (cell->info().Pcondition) continue;
+            if (cell->info().isFictious) continue;
+	    if (cell->info().windowsID != i) continue;
+            capillary_volume = capillary_volume + cell->info().capillaryCellVolume;
+            if (cell->info().p()==0) {
+                water_volume = water_volume + cell->info().capillaryCellVolume;}}}
+    Real saturation = water_volume/capillary_volume;
+    return saturation;
+}
+//----------end temp functions for comparison with experiment-------------------
 
 template <class Solver> 
 void UnsaturatedEngine::computeFacetPoreForcesWithCache(Solver& flow, bool onlyCache)

=== modified file 'pkg/dem/UnsaturatedEngine.hpp'
--- pkg/dem/UnsaturatedEngine.hpp	2014-03-29 15:55:08 +0000
+++ pkg/dem/UnsaturatedEngine.hpp	2014-04-04 17:12:06 +0000
@@ -63,7 +63,10 @@
 		TPL void invade (Solver& flow );
 		TPL Real getMinEntryValue (Solver& flow);
 		TPL Real getSaturation(Solver& flow);
-		
+		TPL Real getWindowsSaturation(Solver&flow, int windowsID);
+		TPL Real getWindowsSaturation1(Solver&flow, int i);
+		TPL Real getWindowsSaturation2(Solver&flow, int i);
+
 		TPL void invadeSingleCell1(Cell_handle cell, double pressure, Solver& flow);
 		TPL void invade1 (Solver& flow );
 		TPL void checkTrap(Solver& flow, double pressure);//check trapped phase, define trapCapP.	
@@ -98,7 +101,7 @@
 		TPL void savePoreBodyInfo(Solver& flow);
 		TPL void savePoreThroatInfo(Solver& flow);
 		TPL void debugTemp(Solver& flow);
-		
+		TPL void initializeCellWindowsID(Solver&flow);
 		TPL void computeSolidLine(Solver& flow);
 		TPL void computeFacetPoreForcesWithCache(Solver& flow, bool onlyCache=false);
 		
@@ -147,6 +150,7 @@
 		void		_invade() {invade(solver);}
 		Real		_getMinEntryValue() {return getMinEntryValue(solver);}
 		Real		_getSaturation() {return getSaturation(solver);}
+		Real		_getWindowsSaturation(int windowsID) {return getWindowsSaturation(solver,windowsID);}
 		void		_saveListNodes() {saveListNodes(solver);}
 		void		_saveListConnection() {saveListConnection(solver);}
  		void		_saveLatticeNodeX(double x) {saveLatticeNodeX(solver,x);}
@@ -157,6 +161,7 @@
  		void		_savePoreBodyInfo(){savePoreBodyInfo(solver);}
  		void		_savePoreThroatInfo(){savePoreThroatInfo(solver);}
  		void		_debugTemp(){debugTemp(solver);}
+ 		void		_initializeCellWindowsID(){initializeCellWindowsID(solver);}
  		void		_computeFacetPoreForcesWithCache(){computeFacetPoreForcesWithCache(solver);}
 		Vector3r 	_fluidForce(unsigned int id_sph) {return fluidForce(id_sph,solver);}
 		
@@ -192,6 +197,7 @@
 					((bool, pressureForce, true,,"Compute the pressure field and associated fluid forces. WARNING: turning off means fluid flow is not computed at all."))
 					((bool, computeForceActivated, true,,"Activate capillary force computation. WARNING: turning off means capillary force is not computed at all, but the drainage can still work."))
 					((bool, invadeBoundary, false,,"Invade from boundaries."))
+					((int, windowsNo, 10,, "Number of genrated windows/(zoomed samples)."))
 					,
 					/*deprec*/
 					,,
@@ -211,7 +217,8 @@
 					.def("getCell",&UnsaturatedEngine::_getCell,(python::arg("pos")),"get id of the cell containing (X,Y,Z).")
 					.def("testFunction",&UnsaturatedEngine::testFunction,"The playground for Chao's experiments.")
 					.def("buildTriangulation",&UnsaturatedEngine::_buildTriangulation,"Triangulate spheres of the current scene.")
-					.def("getSaturation",&UnsaturatedEngine::_getSaturation,"get saturation.")					
+					.def("getSaturation",&UnsaturatedEngine::_getSaturation,"get saturation.")
+					.def("getWindowsSaturation",&UnsaturatedEngine::_getWindowsSaturation,(python::arg("windowsID")), "get saturation of windowsID")
 					.def("getMinEntryValue",&UnsaturatedEngine::_getMinEntryValue,"get the minimum air entry pressure for the next invade step.")
 					.def("saveListNodes",&UnsaturatedEngine::_saveListNodes,"Save the list of nodes.")
 					.def("saveListConnection",&UnsaturatedEngine::_saveListConnection,"Save the connections between cells.")
@@ -223,6 +230,7 @@
 					.def("savePoreBodyInfo",&UnsaturatedEngine::_savePoreBodyInfo,"Save pore bodies positions/Voronoi centers and size/volume.")
 					.def("savePoreThroatInfo",&UnsaturatedEngine::_savePoreThroatInfo,"Save pore throat area, inscribed radius and perimeter.")
 					.def("debugTemp",&UnsaturatedEngine::_debugTemp,"debug temp file.")
+					.def("initializeCellWindowsID",&UnsaturatedEngine::_initializeCellWindowsID,"Initialize cell windows index. A temp function for comparison with experiments, will delete soon")
 					.def("invade",&UnsaturatedEngine::_invade,"Run the drainage invasion.")
 					.def("computeForce",&UnsaturatedEngine::_computeFacetPoreForcesWithCache,"Compute capillary force. ")
 					.def("fluidForce",&UnsaturatedEngine::_fluidForce,(python::arg("Id_sph")),"Return the fluid force on sphere Id_sph.")