yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07360
[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;
};