← Back to team overview

yade-dev team mailing list archive

[Branch ~yade-dev/yade/trunk] Rev 2794: - Update flow code

 

------------------------------------------------------------
revno: 2794
committer: Emanuele Catalano <catalano@xxxxxxxxxxx
branch nick: yade
timestamp: Wed 2011-03-23 10:38:16 +0100
message:
  - Update flow code
modified:
  doc/yade-conferences.bib
  lib/triangulation/FlowBoundingSphere.cpp
  lib/triangulation/FlowBoundingSphere.hpp
  lib/triangulation/Network.cpp
  pkg/dem/FlowEngine.cpp
  pkg/dem/FlowEngine.hpp


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'doc/yade-conferences.bib'
--- doc/yade-conferences.bib	2010-09-09 14:02:49 +0000
+++ doc/yade-conferences.bib	2011-03-23 09:38:16 +0000
@@ -230,3 +230,49 @@
 	month={June}
 }
 
+@InProceedings{ Catalano1,
+	title = {Couplage solide - fluide dans les mod{\`e}les discrets. Approche m{\'e}soscopique.},
+	author = {E. Catalano, B. Chareyre, E. Barth{\'e}l{\'e}my},
+	month = {Jun},
+	year = {2009},
+	booktitle = {GdR MeGe}
+	address = {La Rochelle, France}
+}
+
+@InProceedings{ Catalano2,
+	title = {Couplage solide-fluide par volumes finis {\`a} l'{\'e}chelle porale},
+	author = {B. Chareyre, E. Catalano, E. Barth{\'e}l{\'e}my},
+	month = {Nov},
+	year = {2009},
+	booktitle = {GdR MeGe},
+	address = {Montpellier, France}
+}
+
+@InProceedings{ Catalano3,
+	title = {Fluid-solid coupling in discrete models},
+	author = {E. Catalano, B. Chareyre, E. Barth{\'e}l{\'e}my},
+	month = {Oct},
+	year = {2009},
+	booktitle = {Alert Geomaterials Workshop 2009},
+	address = {Aussois, France}
+}
+
+@InProceedings{ Catalano4,
+	title = {A coupled model for fluid-solid interaction analysis in geomaterials},
+	author = {E. Catalano, B. Chareyre, E. Barth{\'e}l{\'e}my},
+	month = {Oct},
+	year = {2010},
+	booktitle = {Alert Geomaterials Workshop 2010},
+	address = {Aussois, France}
+}
+
+@InProceedings{ Catalano5,
+	title = {Pore scale modelling of Stokes flow},
+	author = {E. Catalano, B. Chareyre, E. Barth{\'e}l{\'e}my},
+	editor = {GdR MeGe},
+	month = {Dec},
+	year = {2010},
+	booktitle = {GdR MeGe},
+	address = {Nantes, France}
+}
+

=== modified file 'lib/triangulation/FlowBoundingSphere.cpp'
--- lib/triangulation/FlowBoundingSphere.cpp	2011-03-22 18:32:04 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2011-03-23 09:38:16 +0000
@@ -67,6 +67,7 @@
 
 FlowBoundingSphere::FlowBoundingSphere()
 {
+	id_Sphere=0;
 	x_min = 1000.0, x_max = -10000.0, y_min = 1000.0, y_max = -10000.0, z_min = 1000.0, z_max = -10000.0;
 	currentTes = 0;
 	nOfSpheres = 0;
@@ -90,6 +91,7 @@
 	RAVERAGE = false; /** use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
 	OUTPUT_BOUDARIES_RADII = false;
 	RAVERAGE = false; /** if true use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
+	areaR2Permeability=true;
 }
 
 void FlowBoundingSphere::ResetNetwork() {noCache=true;}
@@ -390,6 +392,31 @@
 	  }}
 }
 
+vector<Real> FlowBoundingSphere::Average_Fluid_Velocity_On_Sphere(int Id_sph)
+{
+	Average_Cell_Velocity();
+	RTriangulation& Tri = T[currentTes].Triangulation();
+	
+	Real Volumes; CGT::Vecteur VelocityVolumes;
+	vector<Real> result;
+	result.resize(3);
+	
+	VelocityVolumes=CGAL::NULL_VECTOR;
+	Volumes=0.f;
+	
+	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().fictious()==0){
+	    for (int i=0;i<4;i++){
+	      if (cell->vertex(i)->info().id()==Id_sph){
+		VelocityVolumes = VelocityVolumes + cell->info().av_vel()*cell->info().volume();
+		Volumes = Volumes + cell->info().volume();}}}}
+		
+	for (int i=0;i<3;i++) result[i] += VelocityVolumes[i]/Volumes;
+	return result;
+}
+
 void FlowBoundingSphere::Measure_Pore_Pressure(double Wall_up_y, double Wall_down_y)
 {  
 	RTriangulation& Tri = T[currentTes].Triangulation();
@@ -439,6 +466,7 @@
 				Vecteur facetNormal = Surfk/area;
 				const std::vector<Vecteur>& crossSections = cell->info().facetSphereCrossSections;
 				Real fluidSurfRatio = (area-crossSections[j][0]-crossSections[j][1]-crossSections[j][2])/area;
+				if (fluidSurfRatio<0) fluidSurfRatio=-fluidSurfRatio;
 				Vecteur fluidSurfk = cell->info().facetSurfaces[j]*fluidSurfRatio;
 				/// handle fictious vertex since we can get the projected surface easily here
 				if (cell->vertex(j)->info().isFictious) {
@@ -522,6 +550,8 @@
 					/// Apply weighted forces f_k=sqRad_k/sumSqRad*f
 					Vecteur Facet_Unit_Force = -fluidSurfk*cell->info().solidSurfaces[j][3];
 					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().solidSurfaces[j][y];
 						//add to cached value
@@ -533,7 +563,6 @@
 							cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]-facetNormal*crossSections[j][y];
 						}
 					}
-
 				}
 		}
 	else //use cached values
@@ -665,7 +694,7 @@
 //         std::ofstream kFile ( "LocalPermeabilities" ,std::ios::app );
 	Real meanK=0, STDEV=0, meanRadius=0, meanDistance=0;
 	Real infiniteK=1e10;
-
+int count_k_neg=0;
 	for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
 		Point& p1 = cell->info();
 		for (int j=0; j<4; j++) {
@@ -705,7 +734,16 @@
 
 				if (distance!=0) {
 					if (minPermLength>0 && distance_correction) distance=max(minPermLength,distance);
-					k = (M_PI * pow(radius,4)) / (8*viscosity*distance);
+					if (areaR2Permeability){
+						const Vecteur& Surfk = cell->info().facetSurfaces[j];
+						Real area = sqrt(Surfk.squared_length());
+						const Vecteur& crossSections = cell->info().facetSphereCrossSections[j];
+						Real fluidArea=area-crossSections[0]-crossSections[1]-crossSections[2];
+						if (fluidArea<0) fluidArea = -fluidArea;
+						k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);
+						
+					}
+					else k = (M_PI * pow(radius,4)) / (8*viscosity*distance);
 				} else  k = infiniteK;//Will be corrected in the next loop
 
 				(cell->info().k_norm())[j]= k*k_factor;
@@ -801,12 +839,13 @@
 		}
 	}
 	if (DEBUG_OUT) {
-		cout<<grains<<"grains - " <<"Vtotale = " << 2*Vtotale << " Vgrains = " << 2*Vgrains << " Vporale1 = " << 2*(Vtotale-Vgrains) << endl;
-		cout << "Vtotalissimo = " << Vtotalissimo << " Vsolid_tot = " << Vsolid_tot << " Vporale2 = " << Vporale  << " Ssolid_tot = " << Ssolid_tot << endl<< endl;
+		cout<<grains<<"grains - " <<"Vtotale = " << Vtotale << " Vgrains = " << Vgrains << " Vporale1 = " << (Vtotale-Vgrains) << endl;
+		cout << "Vtotalissimo = " << Vtotalissimo/2 << " Vsolid_tot = " << Vsolid_tot/2 << " Vporale2 = " << Vporale/2  << " Ssolid_tot = " << Ssolid_tot << endl<< endl;
 
 		if (!RAVERAGE) cout << "------Hydraulic Radius is used for permeability computation------" << endl << endl;
 		else cout << "------Average Radius is used for permeability computation------" << endl << endl;
-		cout << "-----Computed_Permeability-----" << endl;}
+		cout << "-----Computed_Permeability-----" << endl;
+	cout << "Negative Permeabilities = " << count_k_neg << endl; }
 }
 
 vector<double> FlowBoundingSphere::getConstrictions()
@@ -982,6 +1021,12 @@
                                 for (int j2=0; j2<4; j2++) {
 					if (!Tri.is_infinite(cell->neighbor(j2))) {
                                                 m += (cell->info().k_norm())[j2] * cell->neighbor(j2)->info().p();
+						if ( isinf(m) && j<10 ) cout << "(cell->info().k_norm())[j2] = " << (cell->info().k_norm())[j2] << " cell->neighbor(j2)->info().p() = " << cell->neighbor(j2)->info().p() << endl;
+// 						if (m!=m) 
+// 						{cout << "cell->neighbor(j2)->info().p() =" << cell->neighbor(j2)->info().p() << endl;
+// 						for (int i =0;i<4;i++) cout << "vertex " << i << " = " << cell->neighbor(j2)->vertex(i)->point().point()[0] << " " << cell->neighbor(j2)->vertex(i)->point().point()[1] << " " << cell->neighbor(j2)->vertex(i)->point().point()[2] << endl;}
+// 						else if ( isinf(cell->neighbor(j2)->info().p()) || cell->neighbor(j2)->info().p()> 1e150) cout << "(cell->info().k_norm())[j2] = " << (cell->info().k_norm())[j2] << " cell->neighbor(j2)->info().p() = " << cell->neighbor(j2)->info().p() << endl;
+						
                                                 if (j==0) n += (cell->info().k_norm())[j2];}
                                 }
                                 dp = cell->info().p();
@@ -989,7 +1034,11 @@
                                         //     cell->info().p() =   - ( cell->info().dv() - m ) / ( n );
 					if (j==0) cell->info().inv_sum_k=1/n;
 // 					cell->info().p() = (- (cell->info().dv() - m) / (n) - cell->info().p()) * relax + cell->info().p();
+// if (isinf(cell->info().p())) cout << " PRIMA-- cell->info().dv() = " << cell->info().dv()<< " m =" << m<<" cell->info().inv_sum_k =" << cell->info().inv_sum_k<< " n = " << n << endl;
 					cell->info().p() = (- (cell->info().dv() - m) * cell->info().inv_sum_k - cell->info().p()) * relax + cell->info().p();
+					
+// 					if (isinf(cell->info().p())) cout << " DOPO-- cell->info().dv() = " << cell->info().dv()<< " m =" << m<<" cell->info().inv_sum_k =" << cell->info().inv_sum_k<< " n = " << n << endl;
+					
 					#ifdef GS_OPEN_MP
 // 					double r = sqrt(sqrt(sqrt(cell->info().p())/(1+sqrt(cell->info().p()))));
 // 					if (j % 100 == 0) cout<<"cell->info().p() "<<cell->info().p()<<" vs. "<< (- (cell->info().dv() - m) / (n) - cell->info().p())* relax<<endl;
@@ -1025,9 +1074,9 @@
 		#pragma omp master
 		#endif
 		j++;
-//                  		if (j % 100 == 0) {
+//                  		if (j % 10 == 0) {
 //                         cout << "pmax " << p_max << "; pmoy : " << p_moy << "; dpmax : " << dp_max << endl;
-//                         cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;
+//                         cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;}
 //                         //     save_vtk_file ( Tri );
 //                  }
 	#ifdef GS_OPEN_MP
@@ -1080,6 +1129,7 @@
   for (Tesselation::VCell_iterator it = tmp_cells.begin(); it != cell_up_end; it++)
   {
     Cell_handle& cell = *it;{for (int j2=0; j2<4; j2++) {if (!cell->neighbor(j2)->info().Pcondition){
+      if ((cell->neighbor(j2)->info().p()!=cell->neighbor(j2)->info().p()) && Tri.is_infinite(cell->neighbor(j2))) cout << "oooooooooooooooh" << endl;
     Q1 = Q1 + (cell->neighbor(j2)->info().k_norm())[Tri.mirror_index(cell, j2)]* (cell->neighbor(j2)->info().p()-cell->info().p());
     cellQ1+=1;
     p_out_max = std::max(cell->neighbor(j2)->info().p(), p_out_max);
@@ -1091,6 +1141,7 @@
   for (Tesselation::VCell_iterator it = tmp_cells.begin(); it != cell_down_end; it++)
   {
     Cell_handle& cell = *it;{for (int j2=0; j2<4; j2++) {if (!cell->neighbor(j2)->info().Pcondition){
+      if ((cell->neighbor(j2)->info().p()!=cell->neighbor(j2)->info().p()) && Tri.is_infinite(cell->neighbor(j2))) cout << "oooooooooooooooh2" << endl;
     Q2 = Q2 + (cell->neighbor(j2)->info().k_norm())[Tri.mirror_index(cell, j2)]* (cell->info().p()-cell->neighbor(j2)->info().p());
     cellQ2+=1;
     p_in_max = std::max(cell->neighbor(j2)->info().p(), p_in_max);
@@ -1099,8 +1150,8 @@
   }}}
 
 	if (DEBUG_OUT){
-	cout << "the maximum superior pressure is = " << p_out_max << " the min is = " << p_in_min << endl;
-	cout << "the maximum inferior pressure is = " << p_in_max << " the min is = " << p_out_min << endl;
+	cout << "the maximum superior pressure is = " << p_out_max << " the min is = " << p_out_min << endl;
+	cout << "the maximum inferior pressure is = " << p_in_max << " the min is = " << p_in_min << endl;
 	cout << "superior average pressure is " << p_out_moy/cellQ1 << endl;
         cout << "inferior average pressure is " << p_in_moy/cellQ2 << endl;
         cout << "celle comunicanti in basso = " << cellQ2 << endl;

=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2011-03-22 18:32:04 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2011-03-23 09:38:16 +0000
@@ -37,6 +37,7 @@
  		FlowBoundingSphere();
 
 		bool SLIP_ON_LATERALS;
+		bool areaR2Permeability;
 		double TOLERANCE;
 		double RELAX;
 		double ks; //Hydraulic Conductivity
@@ -53,6 +54,8 @@
 
 		bool RAVERAGE;
 		int walls_id[6];
+		
+		int id_Sphere;
 
 		void mplot ( char *filename);
 		void Localize();
@@ -116,7 +119,7 @@
 		void ApplySinusoidalPressure(RTriangulation& Tri, double Amplitude, double Average_Pressure, double load_intervals);
 		bool isOnSolid  (double X, double Y, double Z);
 		void Measure_Pore_Pressure(double Wall_up_y, double Wall_down_y);
-		
+		vector<Real> Average_Fluid_Velocity_On_Sphere(int Id_sph);
 		//Solver?
 		int useSolver;//(0 : GaussSeidel, 1 : TAUCS, 2 : PARDISO)
 

=== modified file 'lib/triangulation/Network.cpp'
--- lib/triangulation/Network.cpp	2011-03-22 18:32:04 +0000
+++ lib/triangulation/Network.cpp	2011-03-23 09:38:16 +0000
@@ -481,7 +481,7 @@
 	  
 	  Tes.insert((center[0]+Normal[0]*thickness/2)*(1-abs(Normal[0])) + (center[0]+Normal[0]*thickness/2-Normal[0]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]),
 		     (center[1]+Normal[1]*thickness/2)*(1-abs(Normal[1])) + (center[1]+Normal[1]*thickness/2-Normal[1]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]),
-		     (center[2]+Normal[2]*thickness/2)*(1-abs(Normal[2])) + (center[2]+Normal[2]*thickness/2-Normal[2]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]), 
+		     (center[2]+Normal[2]*thickness/2)*(1-abs(Normal[2])) + (center[2]+Normal[2]*thickness/2-Normal[2]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]),
 		     FAR*(Corner_max.y()-Corner_min.y()), id_wall, true);
 	  
  	  Point P (center[0],center[1],center[2]);
@@ -492,7 +492,7 @@
           boundaries[id_wall-id_offset].flowCondition = 1;
           boundaries[id_wall-id_offset].value = 0;
 	  
-	  if(DEBUG_OUT) cout << "A boundary -center/thick- has been created. ID = " << id_wall << " position = " << center[0]+Normal[0]*thickness/2 + (center[0]+Normal[0]*thickness/2-Normal[0]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]) << " , " << center[1]+Normal[0]*thickness/2 + (center[1]+Normal[1]*thickness/2-Normal[1]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]) << " , " <<  center[2]+Normal[0]*thickness/2 + (center[2]+Normal[2]*thickness/2-Normal[2]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]) << ". Radius = " << FAR*(Corner_max.y()-Corner_min.y()) << endl;
+	  if(DEBUG_OUT) cout << "A boundary -center/thick- has been created. ID = " << id_wall << " position = " << (center[0]+Normal[0]*thickness/2)*(1-abs(Normal[0])) + (center[0]+Normal[0]*thickness/2-Normal[0]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]) << " , " << (center[1]+Normal[1]*thickness/2)*(1-abs(Normal[1])) + (center[1]+Normal[1]*thickness/2-Normal[1]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]) << " , " <<  (center[2]+Normal[2]*thickness/2)*(1-abs(Normal[2])) + (center[2]+Normal[2]*thickness/2-Normal[2]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]) << ". Radius = " << FAR*(Corner_max.y()-Corner_min.y()) << endl;
 }
 
 void Network::Define_fictious_cells()

=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp	2011-03-22 18:32:04 +0000
+++ pkg/dem/FlowEngine.cpp	2011-03-23 09:38:16 +0000
@@ -21,6 +21,12 @@
 
 CREATE_LOGGER (FlowEngine);
 
+BOOST_PYTHON_MODULE(_vtkFile){
+        python::scope().attr("__doc__")="SaveVTKFile";
+        YADE_SET_DOCSTRING_OPTS;
+        python::class_<FlowEngine>("FlowEngine","Create file for parallel visualization" )
+                .def("saveVtk",&FlowEngine::saveVtk,"Save VTK File, each dd iterations");}
+
 FlowEngine::~FlowEngine()
 {
 }
@@ -87,6 +93,14 @@
 		{
 			id = V_it->info().id();
 			for (int y=0;y<3;y++) f[y] = (V_it->info().forces)[y];
+			if (id==id_sphere)
+			  id_force = f;
+// 			  CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
+// 			  for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
+// 			    for (int j=0;j<4;j++){
+// 			      if (cell->vertex(j)->info().id() == id_sphere){
+// 				cout << endl << "Cell 
+			    
 			scene->forces.addForce(id, f);
 		}
 		timingDeltas->checkpoint("Applying Forces");
@@ -135,6 +149,8 @@
 	  wall_down_y = flow->y_min;}
 	if (liquefaction) flow->Measure_Pore_Pressure(wall_up_y, wall_down_y);
 	first=false;
+	
+// 	if(id_sphere>=0) flow->Average_Fluid_Velocity_On_Sphere(id_sphere);
 // }
 }
 
@@ -194,6 +210,8 @@
 	flow->DEBUG_OUT = Debug;
 	flow->useSolver = useSolver;
 	flow->VISCOSITY = viscosity;
+	
+	flow->id_Sphere=id_sphere;
 
 	flow->T[flow->currentTes].Clear();
 	flow->T[flow->currentTes].max_id=-1;

=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp	2011-03-22 18:32:04 +0000
+++ pkg/dem/FlowEngine.hpp	2011-03-23 09:38:16 +0000
@@ -12,6 +12,7 @@
 #include<yade/lib/triangulation/FlowBoundingSphere.hpp>
 #include<yade/pkg/dem/TesselationWrapper.hpp>
 
+class Flow;
 class TesselationWrapper;
 #ifdef LINSOLV
 #define FlowSolver CGT::FlowBoundingSphereLinSolv
@@ -52,6 +53,7 @@
 		void clearImposedPressure();
 		Real getFlux(int cond);
 		void saveVtk() {flow->save_vtk_file();}
+		vector<Real> AvFlVelOnSph(int id_sph) {flow->Average_Fluid_Velocity_On_Sphere(id_sph);}
 		python::list getConstrictions() {
 			vector<Real> csd=flow->getConstrictions(); python::list pycsd;
 			for (unsigned int k=0;k<csd.size();k++) pycsd.append(csd[k]); return pycsd;}
@@ -92,12 +94,12 @@
 					((double,viscosity,1.0,,"viscosity of fluid"))
 					((Real,loadFactor,1.1,,"Load multiplicator for oedometer test"))
 					((double, K, 0,, "Permeability of the sample"))
-					((double, MaxPressure, 0,, "Maximal value of water pressure within the sample"))
+					((double, MaxPressure, 0,, "Maximal value of water pressure within the sample - Case of consolidation test"))
 					((double, currentStress, 0,, "Current value of axial stress"))
 					((double, currentStrain, 0,, "Current value of axial strain"))
-					((int, intervals, 30,, "Number of layers for pressure measurements"))
-					((int, useSolver, 0,, "Solver to use"))
-					((bool, liquefaction, false,,"Compute bottom_seabed_pressure if true, see below"))
+					((int, intervals, 30,, "Number of layers for computation average fluid pressure profiles to build consolidation curves"))
+					((int, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso"))
+					((bool, liquefaction, false,,"Measure fluid pressure along the heigth of the sample"))
 // 					((double, bottom_seabed_pressure,0,,"Fluid pressure measured at the bottom of the seabed on the symmetry axe"))
 					((bool, Flow_imposed_TOP_Boundary, true,, "if false involve pressure imposed condition"))
 					((bool, Flow_imposed_BOTTOM_Boundary, true,, "if false involve pressure imposed condition"))
@@ -112,6 +114,7 @@
 					((double, Pressure_LEFT_Boundary,  0,, "Pressure imposed on left boundary"))
 					((double, Pressure_RIGHT_Boundary,  0,, "Pressure imposed on right boundary"))
 					((int, id_sphere, 0,, "Average velocity will be computed for all cells incident to that sphere"))
+					((Vector3r, id_force, 0,, "Fluid force acting on sphere with id=flow.id_sphere"))
 					((bool, BOTTOM_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
 					((bool, TOP_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
 					((bool, RIGHT_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
@@ -126,6 +129,7 @@
 					.def("getFlux",&FlowEngine::getFlux,(python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
 					.def("getConstrictions",&FlowEngine::getConstrictions,"Get the list of constrictions (inscribed circle) for all finite facets.")
 					.def("saveVtk",&FlowEngine::saveVtk,"Save pressure field in vtk format.")
+					.def("AvFlVelOnSph",&FlowEngine::AvFlVelOnSph,(python::arg("Id_sph")),"Compute a sphere-centered average fluid velocity")
 					)
 		DECLARE_LOGGER;
 };