yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3395
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Wed 2014-01-29 18:06:48 +0100
message:
  a test version of computing fluid force.
modified:
  lib/triangulation/FlowBoundingSphere.hpp
  lib/triangulation/FlowBoundingSphere.ipp
  lib/triangulation/Network.hpp
  lib/triangulation/Network.ipp
  pkg/dem/FlowEngine.cpp
  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/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2013-10-25 07:27:55 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2014-01-29 17:06:48 +0000
@@ -65,7 +65,7 @@
 
 		bool RAVERAGE;
 		int walls_id[6];
-		#define parallel_forces
+// 		#define parallel_forces
 		#ifdef parallel_forces
 		int ompThreads;
 		vector< vector<const Vecteur*> > perVertexUnitForce;

=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp	2013-07-26 18:16:04 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp	2014-01-29 17:06:48 +0000
@@ -94,7 +94,9 @@
 	computedOnce=false;
 	minKdivKmean=0.0001;
 	maxKdivKmean=100.;
+	#ifdef parallel_forces
 	ompThreads=1;
+	#endif 
 	errorCode=0;
 }
 

=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp	2014-01-24 18:37:37 +0000
+++ lib/triangulation/Network.hpp	2014-01-29 17:06:48 +0000
@@ -84,6 +84,7 @@
 		double surface_solid_double_fictious_facet(Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3);
 		double surface_solid_facet(Sphere ST1, Sphere ST2, Sphere ST3);
 		
+		void Line_Solid_Pore( Cell_handle cell, int j);
 		void Line_Solid_Pore( Cell_handle cell, int j, bool SLIP_ON_LATERALS, bool reuseFacetData=false);
 		double Line_solid_facet(Sphere ST1, Sphere ST2, Sphere ST3);
 

=== modified file 'lib/triangulation/Network.ipp'
--- lib/triangulation/Network.ipp	2014-01-24 18:37:37 +0000
+++ lib/triangulation/Network.ipp	2014-01-29 17:06:48 +0000
@@ -571,6 +571,11 @@
 // }
 
 template<class Tesselation>
+void Network<Tesselation>::Line_Solid_Pore(Cell_handle cell, int j)
+{
+  Line_Solid_Pore(cell, j, false, true);
+}
+template<class Tesselation>
 void Network<Tesselation>::Line_Solid_Pore(Cell_handle cell, int j, bool SLIP_ON_LATERALS, bool reuseFacetData)
 {
   if (!reuseFacetData)  facetNFictious=detectFacetFictiousVertices(cell,j);
@@ -603,7 +608,7 @@
 		Boundary &bi =  boundary(SV3->info().id());
 
 		Vecteur AB= SV1->point() - SV2->point();
-		AB[bi.coordinate] = 0;
+//		AB[bi.coordinate] = 0;
 		
 		if (bi.flowCondition && ! SLIP_ON_LATERALS) {
                         cell->info().solidLine[j][facetF1]=sqrt(AB.squared_length());

=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp	2013-10-25 07:27:55 +0000
+++ pkg/dem/FlowEngine.cpp	2014-01-29 17:06:48 +0000
@@ -55,8 +55,10 @@
 	  backgroundSolver=solver;
 	  backgroundCompleted=true;
 	}
+	#ifdef parallel_forces	
 	solver->ompThreads = ompThreads>0? ompThreads : omp_get_max_threads();
-
+	#endif
+	
         timingDeltas->checkpoint ( "Triangulating" );
 	UpdateVolumes ( solver );
         timingDeltas->checkpoint ( "Update_Volumes" );

=== modified file 'pkg/dem/UnsaturatedEngine.cpp'
--- pkg/dem/UnsaturatedEngine.cpp	2014-01-27 14:58:47 +0000
+++ pkg/dem/UnsaturatedEngine.cpp	2014-01-29 17:06:48 +0000
@@ -49,8 +49,9 @@
 		initializeCellIndex(solver);//initialize cell index
 		initializePoreRadius(solver);//save all pore radii before invade
 		updateVolumeCapillaryCell(solver);//save capillary volume of all cells, for calculating saturation
+		computeSolidLine(solver);//save cell->info().solidLine[j][y]
 	}
-	solver->noCache = false;
+	solver->noCache = true;
 }
 
 void UnsaturatedEngine::action()
@@ -1214,31 +1215,30 @@
 template <class Solver> 
 void UnsaturatedEngine::computeFacetPoreForcesWithCache(Solver& flow, bool onlyCache)
 {
-	RTriangulation& Tri = flow->T[currentTes].Triangulation();
+	RTriangulation& Tri = flow->T[solver->currentTes].Triangulation();
 	Finite_cells_iterator cell_end = Tri.finite_cells_end();
 	CGT::Vecteur nullVect(0,0,0);
 	//reset forces
 	if (!onlyCache) for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) v->info().forces=nullVect;
 
 	#ifdef parallel_forces
-	if (noCache) {
-		perVertexUnitForce.clear(); perVertexPressure.clear();
+	if (solver->noCache) {
+		solver->perVertexUnitForce.clear(); solver->perVertexPressure.clear();
 // 		vector<const Vecteur*> exf; exf.reserve(20);
 // 		vector<const Real*> exp; exp.reserve(20);
-		perVertexUnitForce.resize(Tri.number_of_vertices());
-		perVertexPressure.resize(Tri.number_of_vertices());}
+		solver->perVertexUnitForce.resize(Tri.number_of_vertices());
+		solver->perVertexPressure.resize(Tri.number_of_vertices());}
 	#endif
 
 	Cell_handle neighbour_cell;
 	Vertex_handle mirror_vertex;
 	CGT::Vecteur tempVect;
 	//FIXME : Ema, be carefull with this (noCache), it needs to be turned true after retriangulation
-	if (noCache) {for (VCell_iterator cell_it=flow->T[currentTes].cellHandles.begin(); cell_it!=flow->T[currentTes].cellHandles.end(); cell_it++){
+	if (solver->noCache) {for (VCell_iterator cell_it=flow->T[currentTes].cellHandles.begin(); cell_it!=flow->T[currentTes].cellHandles.end(); cell_it++){
 // 	if (noCache) for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
 			Cell_handle& cell = *cell_it;
 			//reset cache
 			for (int k=0;k<4;k++) cell->info().unitForceVectors[k]=nullVect;
-
 			for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j))) {
 					neighbour_cell = cell->neighbor(j);
 					const CGT::Vecteur& Surfk = cell->info().facetSurfaces[j];
@@ -1259,8 +1259,7 @@
 					/// Apply weighted forces f_k=sqRad_k/sumSqRad*f
 					CGT::Vecteur Facet_Unit_Force = -fluidSurfk*cell->info().solidLine[j][3];
 					CGT::Vecteur Facet_Force = cell->info().p()*Facet_Unit_Force;
-					
-					
+										
 					for (int y=0; y<3;y++) {
 						cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces + Facet_Force*cell->info().solidLine[j][y];
 						//add to cached value
@@ -1273,12 +1272,12 @@
 						}
 					}
 					#ifdef parallel_forces
-					perVertexUnitForce[cell->vertex(j)->info().id()].push_back(&(cell->info().unitForceVectors[j]));
-					perVertexPressure[cell->vertex(j)->info().id()].push_back(&(cell->info().p()));
+					solver->perVertexUnitForce[cell->vertex(j)->info().id()].push_back(&(cell->info().unitForceVectors[j]));
+					solver->perVertexPressure[cell->vertex(j)->info().id()].push_back(&(cell->info().p()));
 					#endif
 			}
 		}
-		noCache=false;//cache should always be defined after execution of this function
+		solver->noCache=false;//cache should always be defined after execution of this function
 		if (onlyCache) return;
 	} else {//use cached values
 		#ifndef parallel_forces
@@ -1293,8 +1292,8 @@
 			const int& id =  v->info().id();
 			CGT::Vecteur tf (0,0,0);
 			int k=0;
-			for (vector<const Real*>::iterator c = perVertexPressure[id].begin(); c != perVertexPressure[id].end(); c++)
-				tf = tf + (*(perVertexUnitForce[id][k++]))*(**c);
+			for (vector<const Real*>::iterator c = solver->perVertexPressure[id].begin(); c != solver->perVertexPressure[id].end(); c++)
+				tf = tf + (*(solver->perVertexUnitForce[id][k++]))*(**c);
 			v->info().forces = tf;
 		}
 		#endif
@@ -1307,6 +1306,47 @@
 		cout << "TotalForce = "<< TotalForce << endl;}
 }
 
+template<class Solver>
+void UnsaturatedEngine::computeSolidLine(Solver& flow)
+{
+    RTriangulation& Tri = flow->T[solver->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 j=0; j<4; j++) {
+            solver->Line_Solid_Pore(cell, j);
+        }
+    }
+}
+
+template<class Solver>//tempt test, clean later
+void UnsaturatedEngine::vertxID(Solver& flow)
+{
+    ofstream file;
+    file.open("vertexID.txt");
+    file << "vertexID	Pos	Force \n";
+    RTriangulation& Tri = flow->T[solver->currentTes].Triangulation();
+    for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v)	{
+        if (!v->info().isFictious) file<<v->info().id()<<"		"<<v->point().point()<<"		"<<v->info().forces<<endl;
+    }
+    file.close();
+}
+
+template<class Solver>//tempt test. clean later
+void UnsaturatedEngine::testSolidLine(Solver& flow)
+{
+    ofstream file;
+    file.open("solidLine.txt");
+    file << "cellID facet solidLine[j][0] solidLine[j][1] solidLine[j][2] solidLine[j][3] \n";
+    RTriangulation& Tri = flow->T[solver->currentTes].Triangulation();   
+    for (VCell_iterator cell_it=flow->T[currentTes].cellHandles.begin(); cell_it!=flow->T[currentTes].cellHandles.end(); cell_it++){
+	Cell_handle& cell = *cell_it;      
+      for(int j=0; j<4;j++) {
+	file<<cell->info().index<<" "<<j<<" "<<cell->info().solidLine[j][0]<<" "<<cell->info().solidLine[j][1]<<" "<<cell->info().solidLine[j][2]<<" "<<cell->info().solidLine[j][3]<<endl;
+      }
+    }
+    file.close();
+}
+
 YADE_PLUGIN ( ( UnsaturatedEngine ) );
 
 #endif //FLOW_ENGINE

=== modified file 'pkg/dem/UnsaturatedEngine.hpp'
--- pkg/dem/UnsaturatedEngine.hpp	2014-01-27 14:58:47 +0000
+++ pkg/dem/UnsaturatedEngine.hpp	2014-01-29 17:06:48 +0000
@@ -87,6 +87,14 @@
 		TPL void savePoreBodyInfo(Solver& flow);
 		TPL void savePoreThroatInfo(Solver& flow);
 		TPL void debugTemp(Solver& flow);
+		
+		TPL void vertxID(Solver& flow);
+		TPL void testSolidLine(Solver& flow);
+		TPL void computeSolidLine(Solver& flow);
+		TPL void computeFacetPoreForcesWithCache(Solver& flow, bool onlyCache=false);
+		
+		TPL Vector3r fluidForce(unsigned int id_sph, Solver& flow) {
+			const CGT::Vecteur& f=flow->T[flow->currentTes].vertex(id_sph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
 
 		template<class Cellhandle >
 		double getRadiusMin(Cellhandle cell, int j);
@@ -119,11 +127,6 @@
 		double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
 		TPL int getCell(double posX, double posY, double posZ, Solver& flow){return flow->getCell(posX, posY, posZ);}
 
-		bool noCache;//flag for checking if cached values cell->unitForceVectors have been defined		
-		TPL void computeFacetPoreForcesWithCache(Solver& flow, bool onlyCache=false);
-		vector< vector<const CGT::Vecteur*> > perVertexUnitForce;
-		vector< vector<const Real*> > perVertexPressure;
-		
 		void emulateAction(){
 			scene = Omega::instance().getScene().get();
 			action();}
@@ -147,6 +150,11 @@
  		void		_savePoreBodyInfo(){savePoreBodyInfo(solver);}
  		void		_savePoreThroatInfo(){savePoreThroatInfo(solver);}
  		void		_debugTemp(){debugTemp(solver);}
+ 		void		_computeFacetPoreForcesWithCache(){computeFacetPoreForcesWithCache(solver);}
+ 		void		_vertxID(){vertxID(solver);}
+ 		void		_testSolidLine(){testSolidLine(solver);}
+		Vector3r 	_fluidForce(unsigned int id_sph) {return fluidForce(id_sph,solver);}
+ 		
 		virtual ~UnsaturatedEngine();
 
 		virtual void action();
@@ -215,6 +223,10 @@
 					.def("debugTemp",&UnsaturatedEngine::_debugTemp,"debug temp file.")
 					.def("invade",&UnsaturatedEngine::_invade,"Run the drainage invasion from all cells with air pressure. ")
 					.def("invade2",&UnsaturatedEngine::_invade2,"Run the drainage invasion from all cells with air pressure.(version2,water can be trapped in cells) ")
+					.def("computeForce",&UnsaturatedEngine::_computeFacetPoreForcesWithCache,"Test computeFacetPoreForcesWithCache(). ")
+					.def("vertxID",&UnsaturatedEngine::_vertxID,"cout vertxID. ")
+					.def("testSolidLine",&UnsaturatedEngine::_testSolidLine,"For checking solidLine.")
+					.def("fluidForce",&UnsaturatedEngine::_fluidForce,(python::arg("Id_sph")),"Return the fluid force on sphere Id_sph.")					
 					)
 		DECLARE_LOGGER;
 };