~ubuntu-branches/ubuntu/trusty/yade/trusty

                enum {wall_xmin, wall_xmax, wall_ymin, wall_ymax, wall_zmin, wall_zmax};

                #ifdef LINSOLV

                TPL void setForceMetis (Solver& flow, bool force);

                TPL bool getForceMetis (Solver& flow);

                #endif

                TPL Vector3r shearVelocity(unsigned int interaction, Solver& flow) {

                        return (flow->deltaShearVel[interaction]);}

                TPL Vector3r normalVelocity(unsigned int interaction, Solver& flow) {

                        return (flow->deltaNormVel[interaction]);}

                TPL Matrix3r normalStressInteraction(unsigned int interaction, Solver& flow) {

                        return (flow->normalStressInteraction[interaction]);}

                TPL Matrix3r shearStressInteraction(unsigned int interaction, Solver& flow) {

                        return (flow->shearStressInteraction[interaction]);}

                TPL Vector3r normalVect(unsigned int interaction, Solver& flow) {

                        return (flow->normalV[interaction]);}

                TPL Real surfaceDistanceParticle(unsigned int interaction, Solver& flow) {

                        return (flow->surfaceDistance[interaction]);}

                TPL Real edgeSize(Solver& flow) {

                        return (flow->Edge_ids.size());}

                TPL Real OSI(Solver& flow) {

                        return (flow->onlySpheresInteractions.size());}

                TPL int onlySpheresInteractions(unsigned int interaction, Solver& flow) {

                        return (flow->onlySpheresInteractions[interaction]);}

                vector<Real> avFlVelOnSph(unsigned int id_sph) {return solver->Average_Fluid_Velocity_On_Sphere(id_sph);}

//              void setBoundaryVel(Vector3r vel) {topBoundaryVelocity=vel; updateTriangulation=true;}

                void pressureProfile(double wallUpY, double wallDownY) {return solver->measurePressureProfile(wallUpY,wallDownY);}

                double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}

                double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}

                double averagePressure(){return solver->averagePressure();}

                Vector3r        _shearVelocity(unsigned int interaction) {return shearVelocity(interaction,solver);}

                Vector3r        _normalVelocity(unsigned int interaction) {return normalVelocity(interaction,solver);}

                Matrix3r        _normalStressInteraction(unsigned int interaction) {return normalStressInteraction(interaction,solver);}

                Matrix3r        _shearStressInteraction(unsigned int interaction) {return shearStressInteraction(interaction,solver);}

                Vector3r        _normalVect(unsigned int interaction) {return normalVect(interaction,solver);}

                Real            _surfaceDistanceParticle(unsigned int interaction) {return surfaceDistanceParticle(interaction,solver);}

                int             _onlySpheresInteractions(unsigned int interaction) {return onlySpheresInteractions(interaction,solver);}

                void            _setForceMetis (bool force) {setForceMetis(solver,force);}

                bool            _getForceMetis () {return getForceMetis (solver);}

                void            cholmodStats() {

                                        cerr << cholmod_print_common("PFV Cholmod factorization",&(solver->eSolver.cholmod()))<<endl;

                                        cerr << "cholmod method:" << solver->eSolver.cholmod().selected<<endl;

                                        cerr << "METIS called:"<<solver->eSolver.cholmod().called_nd<<endl;}

                bool            metisUsed() {return bool(solver->eSolver.cholmod().called_nd);}

                YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(FlowEngine,PartialEngine,"An engine to solve flow problem in saturated granular media. Model description can be found in [Chareyre2012a]_ and [Catalano2013a]_. See the example script FluidCouplingPFV/oedometer.py. More documentation to come.\n\n.. note::Multi-threading seems to work fine for Cholesky decomposition, but it fails for the solve phase in which -j1 is the fastest, here we specify thread numbers independently using :yref:`FlowEngine::numFactorizeThreads` and :yref:`FlowEngine::numSolveThreads`. These multhreading settings are only impacting the behaviour of openblas library and are relatively independant of :yref:`FlowEngine::multithread`. However, the settings have to be globally consistent. For instance, :yref:`multithread<FlowEngine::multithread>` =True with  yref:`numFactorizeThreads<FlowEngine::numFactorizeThreads>` = yref:`numSolveThreads<FlowEngine::numSolveThreads>` = 4 implies that openblas will mobilize 8 processors at some point. If the system does not have so many procs. it will hurt performance.",

                                        ((double, sineMagnitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied (p )"))

                                        ((double, sineAverage, 0,,"Pressure value (average) when sinusoidal pressure is applied"))

                                        ((double,P_zero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model."))

                                        ((bool, updateTriangulation, 0,,"If true the medium is retriangulated. Can be switched on to force retriangulation after some events (else it will be true periodicaly based on :yref:`FlowEngine::defTolerance` and :yref:`FlowEngine::meshUpdateInterval`. Of course, it costs CPU time."))

                                        ((int,meshUpdateInterval,1000,,"Maximum number of timesteps between re-triangulation events. See also :yref:`FlowEngine::defTolerance`."))

                                        ((double, eps_vol_max, 0,(Attr::readonly),"Maximal absolute volumetric strain computed at each iteration. |yupdate|"))

                                        ((double, defTolerance,0.05,,"Cumulated deformation threshold for which retriangulation of pore space is performed. If negative, the triangulation update will occure with a fixed frequency on the basis of :yref:`FlowEngine::meshUpdateInterval`"))

                                        ((double, porosity, 0,(Attr::readonly),"Porosity computed at each retriangulation |yupdate|"))

                                        ((bool,meanKStat,false,,"report the local permeabilities' correction"))

                                        ((bool,clampKValues,true,,"If true, clamp local permeabilities in [minKdivKmean,maxKdivKmean]*globalK. This clamping can avoid singular values in the permeability matrix and may reduce numerical errors in the solve phase. It will also hide junk values if they exist, or bias all values in very heterogeneous problems. So, use this with care."))

                                        ((double,permeabilityFactor,1.0,,"permability multiplier"))

                                        ((double,viscosity,1.0,,"viscosity of the fluid"))

                                        ((double,stiffness, 10000,,"equivalent contact stiffness used in the lubrication model"))

                                        ((int, xmin,0,(Attr::readonly),"Index of the boundary $x_{min}$. This index is not equal the the id of the corresponding body in general, it may be used to access the corresponding attributes (e.g. flow.bndCondValue[flow.xmin], flow.wallId[flow.xmin],...)."))

                                        ((int, xmax,1,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))

                                        ((int, ymin,2,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))

                                        ((int, ymax,3,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))

                                        ((int, zmin,4,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))

                                        ((int, zmax,5,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))

 

                                        ((vector<bool>, bndCondIsPressure, vector<bool>(6,false),,"defines the type of boundary condition for each side. True if pressure is imposed, False for no-flux. Indexes can be retrieved with :yref:`FlowEngine::xmin` and friends."))

                                        ((vector<double>, bndCondValue, vector<double>(6,0),,"Imposed value of a boundary condition. Only applies if the boundary condition is imposed pressure, else the imposed flux is always zero presently (may be generalized to non-zero imposed fluxes in the future)."))

                                        //FIXME: to be implemented:

                                        ((vector<Vector3r>, boundaryVelocity, vector<Vector3r>(6,Vector3r::Zero()),, "velocity on top boundary, only change it using :yref:`FlowEngine::setBoundaryVel`"))

                                        ((vector<int>, wallIds,vector<int>(6),,"body ids of the boundaries (default values are ok only if aabbWalls are appended before spheres, i.e. numbered 0,...,5)"))

                                        ((vector<bool>, boundaryUseMaxMin, vector<bool>(6,true),,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))

                                        ((bool, display_force, false,,"display the lubrication force applied on particles"))

                                        ((bool, create_file, false,,"create file of velocities"))

                                        ((bool, viscousShear, false,,"Compute viscous shear terms as developped by Donia Marzougui (FIXME: ref.)"))

                                        ((bool, shearLubrication, false,,"Compute shear lubrication force as developped by Brule (FIXME: ref.) "))

                                        ((double, eps, 0.00001,,"roughness defined as a fraction of particles size, giving the minimum distance between particles in the lubrication model."))

//                                      ((bool, forceMetis, 0,,"If true, METIS is used for matrix preconditioning, else Cholmod is free to choose the best method (which may be METIS to, depending on the matrix). See ``nmethods`` in Cholmod documentation"))

                                        for (int i=0; i<6; ++i){normal[i]=Vector3r::Zero(); wallIds[i]=i;}

                                        normal[wall_ymin].y()=normal[wall_xmin].x()=normal[wall_zmin].z()=1;

                                        normal[wall_ymax].y()=normal[wall_xmax].x()=normal[wall_zmax].z()=-1;

                                        updateTriangulation=false;

                                        .def("getBoundaryFlux",&FlowEngine::_getBoundaryFlux,(python::arg("boundary")),"Get total flux through boundary defined by its body id.\n\n.. note:: The flux may be not zero even for no-flow condition. This artifact comes from cells which are incident to two or more boundaries (along the edges of the sample, typically). Such flux evaluation on impermeable boundary is just irrelevant, it does not imply that the boundary condition is not applied properly.")

                                        .def("avFlVelOnSph",&FlowEngine::avFlVelOnSph,(python::arg("Id_sph")),"Compute a sphere-centered average fluid velocity")

                                        .def("shearVelocity",&FlowEngine::_shearVelocity,(python::arg("id_sph")),"Return the shear velocity of the interaction.")

                                        .def("normalVelocity",&FlowEngine::_normalVelocity,(python::arg("id_sph")),"Return the normal velocity of the interaction.")

                                        .def("normalVect",&FlowEngine::_normalVect,(python::arg("id_sph")),"Return the normal vector between particles.")

                                        .def("surfaceDistanceParticle",&FlowEngine::_surfaceDistanceParticle,(python::arg("interaction")),"Return the distance between particles.")

                                        .def("onlySpheresInteractions",&FlowEngine::_onlySpheresInteractions,(python::arg("interaction")),"Return the id of the interaction only between spheres.")

                                        

                                        

                                        .def("pressureProfile",&FlowEngine::pressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")

                                        .def("getPorePressure",&FlowEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")

                                        .def("averageSlicePressure",&FlowEngine::averageSlicePressure,(python::arg("posY")),"Measure slice-averaged pore pressure at height posY")

                                        .def("averagePressure",&FlowEngine::averagePressure,"Measure averaged pore pressure in the entire volume")

                                        .def("emulateAction",&FlowEngine::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")

                                        .def("cholmodStats",&FlowEngine::cholmodStats,"get statistics of cholmod solver activity")

                                        .def("metisUsed",&FlowEngine::metisUsed,"check wether metis lib is effectively used")

                                        .add_property("forceMetis",&FlowEngine::_getForceMetis,&FlowEngine::_setForceMetis,"If true, METIS is used for matrix preconditioning, else Cholmod is free to choose the best method (which may be METIS to, depending on the matrix). See ``nmethods`` in Cholmod documentation")

        updateTriangulation=true;

                Matrix3r        _normalStressInteraction(unsigned int interaction) {return normalStressInteraction(interaction,solver);}

                Matrix3r        _shearStressInteraction(unsigned int interaction) {return shearStressInteraction(interaction,solver);}

                Vector3r        _shearVelocity(unsigned int id_sph) {return shearVelocity(id_sph,solver);}

                Vector3r        _normalVelocity(unsigned int id_sph) {return normalVelocity(id_sph,solver);}

                Vector3r        _normalVect(unsigned int id_sph) {return normalVect(id_sph,solver);}

                Real            _surfaceDistanceParticle(unsigned int interaction) {return surfaceDistanceParticle(interaction,solver);}

                int             _onlySpheresInteractions(unsigned int interaction) {return onlySpheresInteractions(interaction,solver);}

                Real            _edgeSize() {return edgeSize(solver);}

                Real            _OSI() {return OSI(solver);}

                double          getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}

                double          averagePressure(){return solver->averagePressure();}

                void            pressureProfile(double wallUpY, double wallDownY) {return solver->measurePressureProfile(wallUpY,wallDownY);}

                YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(PeriodicFlowEngine,FlowEngine,"A variant of :yref:`FlowEngine` implementing periodic boundary conditions. The API is very similar.",

                        wallIds=vector<int>(6,-1);

//                      wallTopId=wallBottomId=wallFrontId=wallBackId=wallLeftId=wallRightId=-1;

                        updateTriangulation=false;

                        .def("shearVelocity",&PeriodicFlowEngine::_shearVelocity,(python::arg("id_sph")),"Return the shear velocity of the interaction.")

                        .def("normalStressInteraction",&PeriodicFlowEngine::_normalStressInteraction,(python::arg("id_sph")),"Return the shear velocity of the interaction.")

                        .def("shearStressInteraction",&PeriodicFlowEngine::_shearStressInteraction,(python::arg("id_sph")),"Return the shear velocity of the interaction.")

                        .def("normalVelocity",&PeriodicFlowEngine::_normalVelocity,(python::arg("id_sph")),"Return the normal velocity of the interaction.")

                        .def("normalVect",&PeriodicFlowEngine::_normalVect,(python::arg("id_sph")),"Return the normal vector between particles.")

                        .def("surfaceDistanceParticle",&PeriodicFlowEngine::_surfaceDistanceParticle,(python::arg("interaction")),"Return the distance between particles.")

                        .def("onlySpheresInteractions",&PeriodicFlowEngine::_onlySpheresInteractions,(python::arg("interaction")),"Return the id of the interaction only between spheres.")

                        .def("edgeSize",&PeriodicFlowEngine::_edgeSize,"Return the number of interactions.")

                        .def("OSI",&PeriodicFlowEngine::_OSI,"Return the number of interactions only between spheres.")

                        .def("getPorePressure",&PeriodicFlowEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")

                        .def("averagePressure",&PeriodicFlowEngine::averagePressure,"Measure averaged pore pressure in the entire volume")

                        .def("pressureProfile",&PeriodicFlowEngine::pressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")