yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3391
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Mon 2014-01-13 17:54:44 +0100
message:
  add capillaryCellVolume in cellinfo, optimize getSaturation()
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	2013-11-20 08:40:45 +0000
+++ lib/triangulation/def_types.h	2014-01-13 16:54:44 +0000
@@ -223,12 +223,13 @@
 	UnsatCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
   	bool isWaterReservoir;
 	bool isAirReservoir;
+	Real capillaryCellVolume;//for calculating saturation
 	std::vector<double> poreRadius;
 	//pore throat radius for drainage
 	UnsatCellInfo (void)
 	{
 		poreRadius.resize(4, 0);
-		isWaterReservoir = false; isAirReservoir = false;	  
+		isWaterReservoir = false; isAirReservoir = false; capillaryCellVolume = 0;	  
 	}
 };
 

=== modified file 'pkg/dem/UnsaturatedEngine.cpp'
--- pkg/dem/UnsaturatedEngine.cpp	2014-01-13 10:57:20 +0000
+++ pkg/dem/UnsaturatedEngine.cpp	2014-01-13 16:54:44 +0000
@@ -48,6 +48,7 @@
 		Build_Triangulation(P_zero,solver);//create a triangulation and initialize pressure in the elements, everything will be contained in "solver"
 		initializeCellIndex(solver);//initialize cell index
 		initializePoreRadius(solver);//save all pore radii before invade
+		updateVolumeCapillaryCell(solver);//save capillary volume of all cells, for calculating saturation
 	}
 	solver->noCache = false;
 }
@@ -864,29 +865,35 @@
         return volume;
 }
 
-template<class Cellhandle>
-Real UnsaturatedEngine::volumeCapillaryCell ( Cellhandle cell )
+template<class Solver>
+void UnsaturatedEngine::updateVolumeCapillaryCell ( Solver& flow)
 {
-    Initialize_volumes(solver);  
-    Real volume = abs( cell->info().volume() ) - solver->volumeSolidPore(cell);
-    if (volume<0) { cerr<<"volumeCapillaryCell Negative. cell ID: " << cell->info().index << "cell volume: " << cell->info().volume() << "  volumeSolidPore: " << solver->volumeSolidPore(cell) << endl; }
-    return volume;
+    Initialize_volumes(flow);
+    RTriangulation& Tri = flow->T[flow->currentTes].Triangulation();
+    Finite_cells_iterator cell_end = Tri.finite_cells_end();
+    Cell_handle neighbour_cell;
+    for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+        cell->info().capillaryCellVolume = abs( cell->info().volume() ) - solver->volumeSolidPore(cell);
+//         if (cell->info().capillaryCellVolume<0) {cerr<<"volumeCapillaryCell Negative. cell ID: " << cell->info().index << "cell volume: " << cell->info().volume() << "  volumeSolidPore: " << solver->volumeSolidPore(cell) << endl;        }
+    }
 }
 
 template<class Solver>
 Real UnsaturatedEngine::getSaturation (Solver& flow )
 {
+    updateAirReservoir(flow);
+//     updateWaterReservoir(flow);    
     RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
-    Real capillary_volume = 0.0; //total volume
+    Real capillary_volume = 0.0; //total capillary volume
     Real air_volume = 0.0; 	//air volume
     Finite_cells_iterator cell_end = tri.finite_cells_end();
     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.has_vertex() )
-        capillary_volume = capillary_volume + volumeCapillaryCell ( cell );
-        if (cell->info().p()!=0) {
-            air_volume = air_volume + volumeCapillaryCell (cell);
+        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;
@@ -1295,7 +1302,7 @@
     Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
     for ( Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
         if (flow->T[flow->currentTes].Triangulation().is_infinite(cell)) continue;
-        file << cell->info().index << " " << cell->info() << " " << volumeCapillaryCell(cell)<< endl;
+        file << cell->info().index << " " << cell->info() << " " << cell->info().capillaryCellVolume << endl;
     }
     file.close();
 }

=== modified file 'pkg/dem/UnsaturatedEngine.hpp'
--- pkg/dem/UnsaturatedEngine.hpp	2014-01-13 10:57:20 +0000
+++ pkg/dem/UnsaturatedEngine.hpp	2014-01-13 16:54:44 +0000
@@ -97,8 +97,8 @@
 		Real Volume_cell_triple_fictious (Cellhandle cell);
 		template<class Cellhandle>
 		Real Volume_cell (Cellhandle cell);
-		template<class Cellhandle>
-		Real volumeCapillaryCell (Cellhandle cell);
+		template<class Solver>
+		void updateVolumeCapillaryCell ( Solver& flow);		
 		template<class Cellhandle>
 		double computeEffPoreRadius(Cellhandle cell, int j);
 		template<class Cellhandle>