yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3466
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Tue 2014-10-14 11:07:54 +0200
message:
  -add savePhaseVtk. rename function.
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	2014-10-13 12:21:20 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.cpp	2014-10-14 09:07:54 +0000
@@ -19,11 +19,8 @@
 #include "TwoPhaseFlowEngine.hpp"
 
 YADE_PLUGIN((TwoPhaseFlowEngineT));
-
 YADE_PLUGIN((TwoPhaseFlowEngine));
 
-void TwoPhaseFlowEngine::fancyFunction(Real what) {std::cerr<<"yes, I'm a new function"<<std::endl;}
-
 void TwoPhaseFlowEngine::initializeCellIndex()
 {
     int k=0;
@@ -89,7 +86,7 @@
   cell->info().p() = Pg - Pc; 
 }
 
-void TwoPhaseFlowEngine::computePoreThroatRadius()
+void TwoPhaseFlowEngine::computePoreThroatCircleRadius()
 {
   //Calculate the porethroat radii of the inscribed sphere in each pore-body. 
   RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
@@ -209,8 +206,7 @@
 	  M(0,4) = pow((r0+Rin),2);
 	  M(1,4) = pow((r1+Rin),2);
 	  M(2,4) = pow((r2+Rin),2);
-	  M(3,4) = pow((r3+Rin),2);
-	
+	  M(3,4) = pow((r3+Rin),2);	
 	
 	  if (first){
 	    first = false;
@@ -246,5 +242,56 @@
    }
 }
 
+void TwoPhaseFlowEngine::savePhaseVtk(const char* folder)
+{
+// 	RTriangulation& Tri = T[solver->noCache?(!currentTes):currentTes].Triangulation();
+	RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+        static unsigned int number = 0;
+        char filename[80];
+	mkdir(folder, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+        sprintf(filename,"%s/out_%d.vtk",folder,number++);
+	int firstReal=-1;
+
+	//count fictious vertices and cells
+	solver->vtkInfiniteVertices=solver->vtkInfiniteCells=0;
+ 	FiniteCellsIterator cellEnd = Tri.finite_cells_end();
+        for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; 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) solver->vtkInfiniteCells+=1;
+	}
+	for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+                if (!v->info().isReal()) solver->vtkInfiniteVertices+=1;
+                else if (firstReal==-1) firstReal=solver->vtkInfiniteVertices;}
+
+        basicVTKwritter vtkfile((unsigned int) Tri.number_of_vertices()-solver->vtkInfiniteVertices, (unsigned int) Tri.number_of_finite_cells()-solver->vtkInfiniteCells);
+
+        vtkfile.open(filename,"test");
+
+        vtkfile.begin_vertices();
+        double x,y,z;
+        for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+		if (v->info().isReal()){
+		x = (double)(v->point().point()[0]);
+                y = (double)(v->point().point()[1]);
+                z = (double)(v->point().point()[2]);
+                vtkfile.write_point(x,y,z);}
+        }
+        vtkfile.end_vertices();
+
+        vtkfile.begin_cells();
+        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_cell(cell->vertex(0)->info().id()-firstReal, cell->vertex(1)->info().id()-firstReal, cell->vertex(2)->info().id()-firstReal, cell->vertex(3)->info().id()-firstReal);}
+        }
+        vtkfile.end_cells();
+
+	vtkfile.begin_data("Saturation",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().saturation);}
+	}
+	vtkfile.end_data();
+}
+
 #endif //TwoPhaseFLOW
  

=== modified file 'pkg/pfv/TwoPhaseFlowEngine.hpp'
--- pkg/pfv/TwoPhaseFlowEngine.hpp	2014-10-13 12:21:20 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.hpp	2014-10-14 09:07:54 +0000
@@ -26,7 +26,6 @@
 	bool isTrapW;
 	bool isTrapNW;
 	double saturation;//the saturation of single pore (will be used in quasi-static imbibition and dynamic flow)
-	bool isImbibition;//Flag for marking pore unit which contains NW-W interface (used by Thomas)
 	double trapCapP;//for calculating the pressure of trapped pore, cell->info().p() = pressureNW- trapCapP. OR cell->info().p() = pressureW + trapCapP
 	bool hasInterface; //Indicated whether a NW-W interface is present within the pore body
 	std::vector<double> poreThroatRadius;
@@ -68,22 +67,21 @@
 	
 	//If a new function is specific to the derived engine, put it here, else go to the base TemplateFlowEngine
 	//if it is useful for everyone
-	void fancyFunction(Real what);
 	void initializeCellIndex();
 	void computePoreBodyVolume();	
 	void computePoreBodyRadius();
-	void computePoreThroatRadius();
+	void computePoreThroatCircleRadius();
 	void computePoreSatAtInterface(CellHandle cell);
 	void computePoreCapillaryPressure(CellHandle cell);
-
+	void savePhaseVtk(const char* folder);
 	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)"))
 	((bool,initialWetting,true,,"Initial wetting saturated (=true) or non-wetting saturated (=false)"))
 	((bool, isPhaseTrapped,true,,"If True, both phases can be entrapped by the other, which would correspond to snap-off. If false, both phases are always connected to their reservoirs, thus no snap-off."))
-	
+
 	,/*TwoPhaseFlowEngineT()*/,
 	,
-	.def("fancyFunction",&TwoPhaseFlowEngine::fancyFunction,(boost::python::arg("what")=0),"test function")
+	.def("savePhaseVtk",&TwoPhaseFlowEngine::savePhaseVtk,(boost::python::arg("folder")="./phaseVtk"),"Save the saturation of local pores in vtk format. Sw(NW-pore)=0, Sw(W-pore)=1. Specify a folder name for output.")
 	)
 	DECLARE_LOGGER;
 };

=== modified file 'pkg/pfv/UnsaturatedEngine.cpp'
--- pkg/pfv/UnsaturatedEngine.cpp	2014-10-13 12:12:57 +0000
+++ pkg/pfv/UnsaturatedEngine.cpp	2014-10-14 09:07:54 +0000
@@ -14,7 +14,6 @@
 #ifdef UNSATURATED_FLOW
 
 #define TWOPHASEFLOW
-#include "FlowEngine_UnsaturatedEngineT.hpp"
 #include "TwoPhaseFlowEngine.hpp"
 // #include "FlowEngine_TwoPhaseFlowEngineT.hpp"