yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10664
[Branch ~yade-pkg/yade/git-trunk] Rev 3886: - introduce SoluteFlowEngine (from Thomas Sweijen)
------------------------------------------------------------
revno: 3886
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
timestamp: Thu 2014-04-03 15:37:55 +0200
message:
- introduce SoluteFlowEngine (from Thomas Sweijen)
- move Flow classes to a separate folder
removed:
pkg/dem/DFNFlow.cpp
pkg/dem/DummyFlowEngine.cpp
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.hpp
pkg/dem/FlowEngine.ipp
pkg/dem/PeriodicFlowEngine.cpp
pkg/dem/PeriodicFlowEngine.hpp
added:
pkg/pfv/
pkg/pfv/DFNFlow.cpp
pkg/pfv/DummyFlowEngine.cpp
pkg/pfv/FlowEngine.cpp
pkg/pfv/FlowEngine.hpp
pkg/pfv/FlowEngine.ipp
pkg/pfv/PeriodicFlowEngine.cpp
pkg/pfv/PeriodicFlowEngine.hpp
pkg/pfv/SoluteFlowEngine.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
=== removed file 'pkg/dem/DFNFlow.cpp'
--- pkg/dem/DFNFlow.cpp 2014-04-01 13:18:38 +0000
+++ pkg/dem/DFNFlow.cpp 1970-01-01 00:00:00 +0000
@@ -1,85 +0,0 @@
-
-/*************************************************************************
-* Copyright (C) 2014 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-#ifdef FLOW_ENGINE
-
-//keep this #ifdef for commited versions unless you really have stable version that should be compiled by default
-//it will save compilation time for everyone else
-//when you want it compiled, you can pass -DDFNFLOW to cmake, or just uncomment the following line
-// #define DFNFLOW
-#ifdef DFNFLOW
-#define TEMPLATE_FLOW_NAME DFNFlowEngineT
-#include <yade/pkg/dem/FlowEngine.hpp>
-
-class DFNCellInfo : public FlowCellInfo
-{
- public:
- Real anotherVariable;
- void anotherFunction() {};
-};
-
-
-class DFNVertexInfo : public FlowVertexInfo {
- public:
- //same here if needed
-};
-
-typedef TemplateFlowEngine<DFNCellInfo,DFNVertexInfo> DFNFlowEngineT;
-REGISTER_SERIALIZABLE(DFNFlowEngineT);
-YADE_PLUGIN((DFNFlowEngineT));
-
-class DFNFlowEngine : public DFNFlowEngineT
-{
- public :
- void trickPermeability();
- void trickPermeability (RTriangulation::Facet_circulator& facet,Real somethingBig);
- void trickPermeability (RTriangulation::Finite_edges_iterator& edge,Real somethingBig);
-
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DFNFlowEngine,DFNFlowEngineT,"documentation here",
- ((Real, myNewAttribute, 0,,"useless example"))
- ,/*DFNFlowEngineT()*/,
- ,
- )
- DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(DFNFlowEngine);
-YADE_PLUGIN((DFNFlowEngine));
-
-void DFNFlowEngine::trickPermeability (RTriangulation::Facet_circulator& facet, Real somethingBig)
-{
- const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
- const CellHandle& cell1 = facet->first;
- const CellHandle& cell2 = facet->first->neighbor(facet->second);
- if ( Tri.is_infinite(cell1) || Tri.is_infinite(cell2)) cerr<<"Infinite cell found in trickPermeability, should be handled somehow, maybe"<<endl;
- cell1->info().kNorm()[facet->second] = somethingBig;
- cell2->info().kNorm()[Tri.mirror_index(cell1, facet->second)] = somethingBig;
-}
-
-void DFNFlowEngine::trickPermeability(RTriangulation::Finite_edges_iterator& edge, Real somethingBig)
-{
- const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
- RTriangulation::Facet_circulator facet1 = Tri.incident_facets(*edge);
- RTriangulation::Facet_circulator facet0=facet1++;
- trickPermeability(facet0,somethingBig);
- while ( facet1!=facet0 ) {trickPermeability(facet1,somethingBig); facet1++;}
-}
-
-
-void DFNFlowEngine::trickPermeability()
-{
- Real somethingBig=10;//is that big??
- const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
- //We want to change permeability perpendicular to the 10th edge, let's say.
- //in the end this function should have a loop on all edges I guess
- FiniteEdgesIterator edge = Tri.finite_edges_begin();
- for(int k=10; k>0;--k) edge++;
- trickPermeability(edge,somethingBig);
-}
-
-
-#endif //DFNFLOW
-#endif //FLOWENGINE
\ No newline at end of file
=== removed file 'pkg/dem/DummyFlowEngine.cpp'
--- pkg/dem/DummyFlowEngine.cpp 2014-04-01 13:18:38 +0000
+++ pkg/dem/DummyFlowEngine.cpp 1970-01-01 00:00:00 +0000
@@ -1,60 +0,0 @@
-
-/*************************************************************************
-* Copyright (C) 2014 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-// This is an example of how to derive a new FlowEngine with additional data and possibly completely new behaviour.
-// Every functions of the base engine can be overloaded, and new functions can be added
-
-//keep this #ifdef as long as you don't really want to realize a final version publicly, it will save compilation time for everyone else
-//when you want it compiled, you can pass -DDFNFLOW to cmake, or just uncomment the following line
-// #define DUMMYFLOW
-#ifdef DUMMYFLOW
-#define TEMPLATE_FLOW_NAME DummyFlowEngineT
-#include <yade/pkg/dem/FlowEngine.hpp>
-
-/// We can add data to the Info types by inheritance
-class DummyCellInfo : public FlowCellInfo
-{
- public:
- Real anotherVariable;
- void anotherFunction() {};
-};
-
-class DummyVertexInfo : public FlowVertexInfo {
- public:
- //same here if needed
-};
-
-typedef TemplateFlowEngine<DummyCellInfo,DummyVertexInfo> TEMPLATE_FLOW_NAME;
-REGISTER_SERIALIZABLE(TEMPLATE_FLOW_NAME);
-YADE_PLUGIN((TEMPLATE_FLOW_NAME));
-
-class DummyFlowEngine : public TEMPLATE_FLOW_NAME
-{
- public :
- //We can overload every functions of the base engine to make it behave differently
- //if we overload action() like this, this engine is doing nothing in a standard timestep, it can still have useful functions
- virtual void action() {};
-
- //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); {cerr<<"yes, I'm a new function"<<end;}
-
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DummyFlowEngine,DummyFlowEngineT,"documentation here",
- ((Real, myNewAttribute, 0,,"useless example"))
- ,/*DummyFlowEngineT()*/,
- ,
- .def("fancyFunction",&DummyFlowEngine::fancyFunction,(python::arg("what")=0),"test function")
- )
- DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(DummyFlowEngine);
-YADE_PLUGIN((DummyFlowEngine));
-
-void DummyFlowEngine::fancyFunction(Real what) {cerr<<"yes, I'm a new function"<<end;}
-
-#endif //DummyFLOW
\ No newline at end of file
=== removed file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2014-04-01 13:18:38 +0000
+++ pkg/dem/FlowEngine.cpp 1970-01-01 00:00:00 +0000
@@ -1,40 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
-* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-#ifdef YADE_CGAL
-#ifdef FLOW_ENGINE
-
-#define TEMPLATE_FLOW_NAME FlowEngineT
-#include "FlowEngine.hpp"
-
-// To register properly, we need to first instantiate an intermediate class, then inherit from it with correct class names in YADE_CLASS macro
-// The intermediate one would be seen with the name "TemplateFlowEngine" by python, thus it would not work when more than one class are derived, they would all
-// be named "TemplateFlowEngine" ...
-typedef TemplateFlowEngine<FlowCellInfo,FlowVertexInfo> FlowEngineT;
-REGISTER_SERIALIZABLE(FlowEngineT);
-
-class FlowEngine : public FlowEngineT
-{
- public :
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(FlowEngine,FlowEngineT,"An engine to solve flow problem in saturated granular media. Model description can be found in [Chareyre2012a]_ and [Catalano2014a]_. 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.",
- ,,
- ,
- //nothing special to define here, we simply re-use FlowEngine methods
- //.def("meanVelocity",&PeriodicFlowEngine::meanVelocity,"measure the mean velocity in the period")
- )
- DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(FlowEngine);
-
-YADE_PLUGIN((FlowEngineT));
-CREATE_LOGGER(FlowEngine );
-YADE_PLUGIN((FlowEngine));
-
-#endif //FLOW_ENGINE
-
-#endif /* YADE_CGAL */
-
=== removed file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2014-04-03 13:07:09 +0000
+++ pkg/dem/FlowEngine.hpp 1970-01-01 00:00:00 +0000
@@ -1,433 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
-* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-/*!
-
-FlowEngine is an interface between Yade and the fluid solvers defined in lib/triangulation and using the PFV scheme.
-There are also a few non-trivial functions defined here, such has building triangulation and computating elements volumes.
-The code strongly relies on CGAL library for building triangulations on the top of sphere packings.
-CGAL's trangulations introduce the data associated to each element of the triangulation through template parameters (Cell_info and Vertex_info).
-FlowEngine is thus defined via the template class TemplateFlowEngine(Cellinfo,VertexInfo), for easier addition of variables in various couplings.
-PeriodicFlowEngine is a variant for periodic boundary conditions.
-
-Which solver will be actually used internally to obtain pore pressure will depend partly on compile time flags (libcholmod available and #define LINSOLV),
-and on runtime settings (useSolver=0: Gauss-Seidel (iterative), useSolver=1: CHOLMOD if available (direct sparse cholesky)).
-
-The files defining lower level classes are in yade/lib. The code uses CGAL::Triangulation3 for managing the mesh and storing data. Eigen3::Sparse, suitesparse::cholmod, and metis are used for solving the linear systems with a direct method (Cholesky). An iterative method (Gauss-Seidel) is implemented directly in Yade and can be used as a fallback (see FlowEngine::useSolver).
-
-Most classes in lib/triangulation are templates, and are therefore completely defined in header files.
-A pseudo hpp/cpp split is reproduced for clarity, with hpp/ipp extensions (but again, in the end they are all include files).
-The files are
-- RegularTriangulation.h : declaration of info types embeded in the mesh (template parameters of CGAL::Triangulation3), and instanciation of RTriangulation classes
-- Tesselation.h/ipp : encapsulate RegularTriangulation and adds functions to manipulate the dual (Voronoi) graph of the triangulation
-- Network.hpp/ipp : specialized for PFV model (the two former are used independently by TesselationWrapper), a set of functions to determine volumes and surfaces of intersections between spheres and various subdomains. Contains two triangulations for smooth transitions while remeshing - e.g. interpolating values in the new mesh using the previous one.
-- FlowBoundingSphere.hpp/ipp and PeriodicFlow.hpp/ipp + LinSolv variants: implement the solver in itself (mesh, boundary conditions, solving, defining fluid-particles interactions)
-- FlowEngine.hpp/ipp/cpp (this file)
-
-Variants for periodic boundary conditions are also present.
-
-*/
-
-#pragma once
-#include<yade/core/PartialEngine.hpp>
-#include<yade/pkg/dem/TriaxialCompressionEngine.hpp>
-#include<yade/pkg/dem/TesselationWrapper.hpp>
-#include<yade/lib/triangulation/FlowBoundingSphere.hpp>
-#include<yade/lib/triangulation/PeriodicFlow.hpp>
-
-/// Frequently used:
-typedef CGT::CVector CVector;
-typedef CGT::Point Point;
-
-/// Converters for Eigen and CGAL vectors
-inline CVector makeCgVect ( const Vector3r& yv ) {return CVector ( yv[0],yv[1],yv[2] );}
-inline Point makeCgPoint ( const Vector3r& yv ) {return Point ( yv[0],yv[1],yv[2] );}
-inline Vector3r makeVector3r ( const Point& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
-inline Vector3r makeVector3r ( const CVector& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
-
-/// C++ templates and YADE_CLASS_... macro family are not very compatible, this #define is a hack to make it work
-/// TEMPLATE_FLOW_NAME will be the name of a serializable TemplateFlowEngine<...> instance, which can in turn be
-/// inherited from. The instance itself will be useless for actual simulations.
-#ifndef TEMPLATE_FLOW_NAME
-#error You must define TEMPLATE_FLOW_NAME in your source file before including FlowEngine.hpp
-#endif
-
-template<class _CellInfo, class _VertexInfo, class _Tesselation=CGT::_Tesselation<CGT::TriangulationTypes<_VertexInfo,_CellInfo> >, class solverT=CGT::FlowBoundingSphere<_Tesselation> >
-class TemplateFlowEngine : public PartialEngine
-{
- public :
- #ifdef LINSOLV
- typedef CGT::FlowBoundingSphereLinSolv<typename solverT::Tesselation,solverT> FlowSolver;
- #else
- typedef solverT FlowSolver;
- #endif
- typedef FlowSolver Solver;//FIXME: useless alias, search/replace then remove
- typedef typename FlowSolver::Tesselation Tesselation;
- typedef typename FlowSolver::RTriangulation RTriangulation;
- typedef typename FlowSolver::FiniteVerticesIterator FiniteVerticesIterator;
- typedef typename FlowSolver::FiniteCellsIterator FiniteCellsIterator;
- typedef typename FlowSolver::CellHandle CellHandle;
- typedef typename FlowSolver::FiniteEdgesIterator FiniteEdgesIterator;
- typedef typename FlowSolver::VertexHandle VertexHandle;
- typedef typename RTriangulation::Triangulation_data_structure::Cell::Info CellInfo;
- typedef typename RTriangulation::Triangulation_data_structure::Vertex::Info VertexInfo;
-
- protected:
- shared_ptr<FlowSolver> solver;
- shared_ptr<FlowSolver> backgroundSolver;
- volatile bool backgroundCompleted;
- Cell cachedCell;
- struct posData {Body::id_t id; Vector3r pos; Real radius; bool isSphere; bool exists; posData(){exists=0;}};
- vector<posData> positionBufferCurrent;//reflect last known positions before we start computations
- vector<posData> positionBufferParallel;//keep the positions from a given step for multithread factorization
- //copy positions in a buffer for faster and/or parallel access
- void setPositionsBuffer(bool current);
- virtual void trickPermeability() {};
-
- public :
- int retriangulationLastIter;
- enum {wall_xmin, wall_xmax, wall_ymin, wall_ymax, wall_zmin, wall_zmax};
- Vector3r normal [6];
- bool currentTes;
- bool metisForced;
- int idOffset;
- double epsVolCumulative;
- int ReTrg;
- int ellapsedIter;
- void initSolver (FlowSolver& flow);
- #ifdef LINSOLV
- void setForceMetis (bool force);
- bool getForceMetis ();
- #endif
- void triangulate (Solver& flow);
- void addBoundary (Solver& flow);
- void buildTriangulation (double pZero, Solver& flow);
- void buildTriangulation (Solver& flow);
- void updateVolumes (Solver& flow);
- void initializeVolumes (Solver& flow);
- void boundaryConditions(Solver& flow);
- void updateBCs ( Solver& flow ) {
- if (flow.T[flow.currentTes].maxId>0) boundaryConditions(flow);//avoids crash at iteration 0, when the packing is not bounded yet
- else LOG_ERROR("updateBCs not applied");
- flow.pressureChanged=true;}
-
- void imposeFlux(Vector3r pos, Real flux);
- unsigned int imposePressure(Vector3r pos, Real p);
- void setImposedPressure(unsigned int cond, Real p);
- void clearImposedPressure();
- void clearImposedFlux();
- void computeViscousForces(Solver& flow);
- Real getCellFlux(unsigned int cond);
- Real getBoundaryFlux(unsigned int boundary) {return solver->boundaryFlux(boundary);}
- Vector3r fluidForce(unsigned int idSph) {
- const CGT::CVector& f=solver->T[solver->currentTes].vertex(idSph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
-
- Vector3r shearLubForce(unsigned int id_sph) {
- return (solver->shearLubricationForces.size()>id_sph)?solver->shearLubricationForces[id_sph]:Vector3r::Zero();}
- Vector3r shearLubTorque(unsigned int id_sph) {
- return (solver->shearLubricationTorques.size()>id_sph)?solver->shearLubricationTorques[id_sph]:Vector3r::Zero();}
- Vector3r pumpLubTorque(unsigned int id_sph) {
- return (solver->pumpLubricationTorques.size()>id_sph)?solver->pumpLubricationTorques[id_sph]:Vector3r::Zero();}
- Vector3r twistLubTorque(unsigned int id_sph) {
- return (solver->twistLubricationTorques.size()>id_sph)?solver->twistLubricationTorques[id_sph]:Vector3r::Zero();}
- Vector3r normalLubForce(unsigned int id_sph) {
- return (solver->normalLubricationForce.size()>id_sph)?solver->normalLubricationForce[id_sph]:Vector3r::Zero();}
- Matrix3r bodyShearLubStress(unsigned int id_sph) {
- return (solver->shearLubricationBodyStress.size()>id_sph)?solver->shearLubricationBodyStress[id_sph]:Matrix3r::Zero();}
- Matrix3r bodyNormalLubStress(unsigned int id_sph) {
- return (solver->normalLubricationBodyStress.size()>id_sph)?solver->normalLubricationBodyStress[id_sph]:Matrix3r::Zero();}
- Vector3r shearVelocity(unsigned int interaction) {
- return (solver->deltaShearVel[interaction]);}
- Vector3r normalVelocity(unsigned int interaction) {
- return (solver->deltaNormVel[interaction]);}
- Matrix3r normalStressInteraction(unsigned int interaction) {
- return (solver->normalStressInteraction[interaction]);}
- Matrix3r shearStressInteraction(unsigned int interaction) {
- return (solver->shearStressInteraction[interaction]);}
- Vector3r normalVect(unsigned int interaction) {
- return (solver->normalV[interaction]);}
- Real surfaceDistanceParticle(unsigned int interaction) {
- return (solver->surfaceDistance[interaction]);}
- Real edgeSize() {
- return (solver->edgeIds.size());}
- Real OSI() {
- return (solver->onlySpheresInteractions.size());}
- int onlySpheresInteractions(unsigned int interaction) {
- return (solver->onlySpheresInteractions[interaction]);}
- python::list getConstrictions(bool all) {
- vector<Real> csd=solver->getConstrictions(); python::list pycsd;
- for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k]!=0) || csd[k]>0) pycsd.append(csd[k]); return pycsd;}
- python::list getConstrictionsFull(bool all) {
- vector<Constriction> csd=solver->getConstrictionsFull(); python::list pycsd;
- for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k].second[0]!=0) || csd[k].second[0]>0) {
- python::list cons;
- cons.append(csd[k].first.first);
- cons.append(csd[k].first.second);
- cons.append(csd[k].second[0]);
- cons.append(csd[k].second[1]);
- cons.append(csd[k].second[2]);
- cons.append(csd[k].second[3]);
- pycsd.append(cons);}
- return pycsd;}
-
- template<class Cellhandle>
- Real volumeCellSingleFictious (Cellhandle cell);
- template<class Cellhandle>
- Real volumeCellDoubleFictious (Cellhandle cell);
- template<class Cellhandle>
- Real volumeCellTripleFictious (Cellhandle cell);
- template<class Cellhandle>
- Real volumeCell (Cellhandle cell);
- void Oedometer_Boundary_Conditions();
- void averageRealCellVelocity();
- void saveVtk(const char* folder) {solver->saveVtk(folder);}
- vector<Real> avFlVelOnSph(unsigned int idSph) {return solver->averageFluidVelocityOnSphere(idSph);}
-
-// 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]);}
- int getCell(double posX, double posY, double posZ){return solver->getCell(posX, posY, posZ);}
- unsigned int nCells(){return solver->T[solver->currentTes].cellHandles.size();}
- python::list getVertices(unsigned int id){
- python::list ids;
- if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return ids;}
- for (unsigned int i=0;i<4;i++) ids.append(solver->T[solver->currentTes].cellHandles[id]->vertex(i)->info().id());
- return ids;
- }
- double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
- double averagePressure(){return solver->averagePressure();}
-
- void emulateAction(){
- scene = Omega::instance().getScene().get();
- action();}
- #ifdef LINSOLV
- void exportMatrix(string filename) {if (useSolver==3) solver->exportMatrix(filename.c_str());
- else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
- void exportTriplets(string filename) {if (useSolver==3) solver->exportTriplets(filename.c_str());
- else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
- void cholmodStats() {
- cerr << cholmod_print_common((char*)string("PFV Cholmod factorization").c_str(),&(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);}
- #endif
-
- virtual ~TemplateFlowEngine();
- virtual void action();
- virtual void backgroundAction();
-
- //commodities
- void compTessVolumes() {
- solver->T[solver->currentTes].compute();
- solver->T[solver->currentTes].computeVolumes();
- }
- Real getVolume (Body::id_t id) {
- if (solver->T[solver->currentTes].Max_id() <= 0) {emulateAction(); LOG_WARN("Not triangulated yet, emulating action");}
- if (solver->T[solver->currentTes].Volume(id) == -1) {compTessVolumes(); LOG_WARN("Computing all volumes now, as you did not request it explicitely.");}
- return (solver->T[solver->currentTes].Max_id() >= id) ? solver->T[solver->currentTes].Volume(id) : -1;}
-
- YADE_CLASS_PYCLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TemplateFlowEngine,TEMPLATE_FLOW_NAME,PartialEngine,"An engine to solve flow problem in saturated granular media. Model description can be found in [Chareyre2012a]_ and [Catalano2014a]_. 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.",
- ((bool,isActivated,true,,"Activates Flow Engine"))
- ((bool,first,true,,"Controls the initialization/update phases"))
- ((double, fluidBulkModulus, 0.,,"Bulk modulus of fluid (inverse of compressibility) K=-dP*V/dV [Pa]. Flow is compressible if fluidBulkModulus > 0, else incompressible."))
- ((Real, dt, 0,,"timestep [s]"))
- ((bool,permeabilityMap,false,,"Enable/disable stocking of average permeability scalar in cell infos."))
- ((bool, slipBoundary, true,, "Controls friction condition on lateral walls"))
- ((bool,waveAction, false,, "Allow sinusoidal pressure condition to simulate ocean waves"))
- ((double, sineMagnitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied (p )"))
- ((double, sineAverage, 0,,"Pressure value (average) when sinusoidal pressure is applied"))
- ((bool, debug, false,,"Activate debug messages"))
- ((double, wallThickness,0.001,,"Walls thickness"))
- ((double,pZero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model."))
- ((double,tolerance,1e-06,,"Gauss-Seidel tolerance"))
- ((double,relax,1.9,,"Gauss-Seidel relaxation"))
- ((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, epsVolMax, 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."))
- ((Real,minKdivKmean,0.0001,,"define the min K value (see :yref:`FlowEngine::clampKValues`)"))
- ((Real,maxKdivKmean,100,,"define the max K value (see :yref:`FlowEngine::clampKValues`)"))
- ((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, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso, 3-CHOLMOD"))
- ((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`"))
- ((int, ignoredBody,-1,,"Id of a sphere to exclude from the triangulation.)"))
- ((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.) "))
- ((bool, pumpTorque, false,,"Compute pump torque applied on particles "))
- ((bool, twistTorque, false,,"Compute twist torque applied on particles "))
- ((double, eps, 0.00001,,"roughness defined as a fraction of particles size, giving the minimum distance between particles in the lubrication model."))
- ((bool, pressureForce, true,,"compute the pressure field and associated fluid forces. WARNING: turning off means fluid flow is not computed at all."))
- ((bool, normalLubrication, false,,"compute normal lubrication force as developped by Brule"))
- ((bool, viscousNormalBodyStress, false,,"compute normal viscous stress applied on each body"))
- ((bool, viscousShearBodyStress, false,,"compute shear viscous stress applied on each body"))
- ((bool, multithread, false,,"Build triangulation and factorize in the background (multi-thread mode)"))
- #ifdef EIGENSPARSE_LIB
- ((int, numSolveThreads, 1,,"number of openblas threads in the solve phase."))
- ((int, numFactorizeThreads, 1,,"number of openblas threads in the factorization phase"))
- #endif
- ,
- /*deprec*/
- ((meanK_opt,clampKValues,"the name changed"))
- ,,
- timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
- 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;
- solver = shared_ptr<FlowSolver> (new FlowSolver);
- first=true;
- epsVolMax=epsVolCumulative=retriangulationLastIter=0;
- ReTrg=1;
- backgroundCompleted=true;
- ellapsedIter=0;
- metisForced=false;
- ,
- .def("imposeFlux",&TemplateFlowEngine::imposeFlux,(python::arg("pos"),python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
- .def("imposePressure",&TemplateFlowEngine::imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
- .def("setImposedPressure",&TemplateFlowEngine::setImposedPressure,(python::arg("cond"),python::arg("p")),"Set pressure value at the point indexed 'cond'.")
- .def("clearImposedPressure",&TemplateFlowEngine::clearImposedPressure,"Clear the list of points with pressure imposed.")
- .def("clearImposedFlux",&TemplateFlowEngine::clearImposedFlux,"Clear the list of points with flux imposed.")
- .def("getCellFlux",&TemplateFlowEngine::getCellFlux,(python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
- .def("getBoundaryFlux",&TemplateFlowEngine::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("getConstrictions",&TemplateFlowEngine::getConstrictions,(python::arg("all")=true),"Get the list of constriction radii (inscribed circle) for all finite facets (if all==True) or all facets not incident to a virtual bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres.")
- .def("getConstrictionsFull",&TemplateFlowEngine::getConstrictionsFull,(python::arg("all")=true),"Get the list of constrictions (inscribed circle) for all finite facets (if all==True), or all facets not incident to a fictious bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres. The constrictions are returned in the format {{cell1,cell2}{rad,nx,ny,nz}}")
- .def("edgeSize",&TemplateFlowEngine::edgeSize,"Return the number of interactions.")
- .def("OSI",&TemplateFlowEngine::OSI,"Return the number of interactions only between spheres.")
-
- .def("saveVtk",&TemplateFlowEngine::saveVtk,(python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
- .def("avFlVelOnSph",&TemplateFlowEngine::avFlVelOnSph,(python::arg("idSph")),"compute a sphere-centered average fluid velocity")
- .def("fluidForce",&TemplateFlowEngine::fluidForce,(python::arg("idSph")),"Return the fluid force on sphere idSph.")
- .def("shearLubForce",&TemplateFlowEngine::shearLubForce,(python::arg("idSph")),"Return the shear lubrication force on sphere idSph.")
- .def("shearLubTorque",&TemplateFlowEngine::shearLubTorque,(python::arg("idSph")),"Return the shear lubrication torque on sphere idSph.")
- .def("normalLubForce",&TemplateFlowEngine::normalLubForce,(python::arg("idSph")),"Return the normal lubrication force on sphere idSph.")
- .def("bodyShearLubStress",&TemplateFlowEngine::bodyShearLubStress,(python::arg("idSph")),"Return the shear lubrication stress on sphere idSph.")
- .def("bodyNormalLubStress",&TemplateFlowEngine::bodyNormalLubStress,(python::arg("idSph")),"Return the normal lubrication stress on sphere idSph.")
- .def("shearVelocity",&TemplateFlowEngine::shearVelocity,(python::arg("idSph")),"Return the shear velocity of the interaction.")
- .def("normalVelocity",&TemplateFlowEngine::normalVelocity,(python::arg("idSph")),"Return the normal velocity of the interaction.")
- .def("normalVect",&TemplateFlowEngine::normalVect,(python::arg("idSph")),"Return the normal vector between particles.")
- .def("surfaceDistanceParticle",&TemplateFlowEngine::surfaceDistanceParticle,(python::arg("interaction")),"Return the distance between particles.")
- .def("onlySpheresInteractions",&TemplateFlowEngine::onlySpheresInteractions,(python::arg("interaction")),"Return the id of the interaction only between spheres.")
- .def("pressureProfile",&TemplateFlowEngine::pressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
- .def("getPorePressure",&TemplateFlowEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
- .def("averageSlicePressure",&TemplateFlowEngine::averageSlicePressure,(python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
- .def("averagePressure",&TemplateFlowEngine::averagePressure,"Measure averaged pore pressure in the entire volume")
- .def("updateBCs",&TemplateFlowEngine::updateBCs,"tells the engine to update it's boundary conditions before running (especially useful when changing boundary pressure - should not be needed for point-wise imposed pressure)")
- .def("emulateAction",&TemplateFlowEngine::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
- .def("getCell",&TemplateFlowEngine::getCell,(python::arg("pos")),"get id of the cell containing (X,Y,Z).")
- .def("nCells",&TemplateFlowEngine::nCells,"get the total number of finite cells in the triangulation.")
- .def("getVertices",&TemplateFlowEngine::getVertices,(python::arg("id")),"get the vertices of a cell")
- #ifdef LINSOLV
- .def("exportMatrix",&TemplateFlowEngine::exportMatrix,(python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
- .def("exportTriplets",&TemplateFlowEngine::exportTriplets,(python::arg("filename")="triplets"),"Export system matrix to a file with only non-zero entries.")
- .def("cholmodStats",&TemplateFlowEngine::cholmodStats,"get statistics of cholmod solver activity")
- .def("metisUsed",&TemplateFlowEngine::metisUsed,"check wether metis lib is effectively used")
- .add_property("forceMetis",&TemplateFlowEngine::getForceMetis,&TemplateFlowEngine::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")
- #endif
- .def("compTessVolumes",&TemplateFlowEngine::compTessVolumes,"Like TesselationWrapper::computeVolumes()")
- .def("volume",&TemplateFlowEngine::getVolume,(python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
- )
-};
-// Definition of functions in a separate file for clarity
-#include<yade/pkg/dem/FlowEngine.ipp>
-
-class FlowCellInfo : public CGT::SimpleCellInfo {
- public:
- //For vector storage of all cells
- unsigned int index;
- int volumeSign;
- bool Pcondition;
- Real invVoidV;
- Real t;
- int fict;
- Real volumeVariation;
- double pression;
- //average relative (fluid - facet translation) velocity defined for a single cell as 1/Volume * SUM_ON_FACETS(x_average_facet*average_facet_flow_rate)
- CVector averageCellVelocity;
- // Surface vectors of facets, pointing from outside toward inside the cell
- std::vector<CVector> facetSurfaces;
- //Ratio between fluid surface and facet surface
- std::vector<Real> facetFluidSurfacesRatio;
- // Reflects the geometrical property of the cell, so that the force by cell fluid on grain "i" is pressure*unitForceVectors[i]
- std::vector<CVector> unitForceVectors;
- // Store the area of triangle-sphere intersections for each facet (used in forces definition)
- std::vector<CVector> facetSphereCrossSections;
- std::vector<CVector> cellForce;
- std::vector<double> rayHydr;
- std::vector<double> modulePermeability;
- // Partial surfaces of spheres in the double-tetrahedron linking two voronoi centers. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_surfaces in the facet.
- double solidSurfaces [4][4];
-
- FlowCellInfo (void)
- {
- modulePermeability.resize(4, 0);
- cellForce.resize(4);
- facetSurfaces.resize(4);
- facetFluidSurfacesRatio.resize(4);
- facetSphereCrossSections.resize(4);
- unitForceVectors.resize(4);
- for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidSurfaces[k][l]=0;
- rayHydr.resize(4, 0);
- invSumK=index=volumeSign=s=volumeVariation=pression=invVoidV=fict=0;
- isFictious=false; Pcondition = false; isGhost = false;
- isvisited = false;
- isGhost=false;
- }
- bool isGhost;
- double invSumK;
- bool isvisited;
-
- inline Real& volume (void) {return t;}
- inline const Real& invVoidVolume (void) const {return invVoidV;}
- inline Real& invVoidVolume (void) {return invVoidV;}
- inline Real& dv (void) {return volumeVariation;}
- inline int& fictious (void) {return fict;}
- inline double& p (void) {return pression;}
- inline const double shiftedP (void) const {return pression;} //For compatibility with the periodic case
- inline const std::vector<double>& kNorm (void) const {return modulePermeability;}
- inline std::vector<double>& kNorm (void) {return modulePermeability;}
- inline std::vector< CVector >& facetSurf (void) {return facetSurfaces;}
- inline std::vector<CVector>& force (void) {return cellForce;}
- inline std::vector<double>& Rh (void) {return rayHydr;}
- inline CVector& averageVelocity (void) {return averageCellVelocity;}
- //used for transfering values between two triangulations, overload with more variables in derived classes (see e.g. SoluteFlow)
- inline void getInfo(const FlowCellInfo& otherCellInfo) {p()=otherCellInfo.shiftedP();}
-};
-
-class FlowVertexInfo : public CGT::SimpleVertexInfo {
- CVector grainVelocity;
- Real volumeIncidentCells;
-public:
- CVector forces;
- bool isGhost;
- FlowVertexInfo (void) {isGhost=false;}
- inline CVector& force (void) {return forces;}
- inline CVector& vel (void) {return grainVelocity;}
- inline Real& volCells (void) {return volumeIncidentCells;}
- inline const CVector ghostShift (void) {return CGAL::NULL_VECTOR;}
-};
-
-
-
=== removed file 'pkg/dem/FlowEngine.ipp'
--- pkg/dem/FlowEngine.ipp 2014-04-01 13:18:38 +0000
+++ pkg/dem/FlowEngine.ipp 1970-01-01 00:00:00 +0000
@@ -1,730 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
-* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-#ifdef YADE_CGAL
-
-#ifdef FLOW_ENGINE
-#include<yade/core/Scene.hpp>
-#include<yade/lib/base/Math.hpp>
-#include<yade/pkg/dem/TesselationWrapper.hpp>
-#include<yade/pkg/common/Sphere.hpp>
-#include<yade/pkg/common/Wall.hpp>
-#include<yade/pkg/common/Box.hpp>
-#include <sys/stat.h>
-#include <sys/types.h>
-#include <boost/thread.hpp>
-#include <boost/date_time.hpp>
-#include <boost/date_time/posix_time/posix_time.hpp>
-
-#ifdef LINSOLV
-#include <cholmod.h>
-#endif
-
-#include "FlowEngine.hpp"
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::~TemplateFlowEngine() {}
-
-// YADE_PLUGIN((TFlowEng));
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-unsigned int TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposePressure(Vector3r pos, Real p)
-{
-// if (!flow) LOG_ERROR("no flow defined yet, run at least one iter");
- solver->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),p) );
- //force immediate update of boundary conditions
- updateTriangulation=true;
- return solver->imposedP.size()-1;
-}
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::action()
-{
- if ( !isActivated ) return;
- timingDeltas->start();
- setPositionsBuffer(true);
- timingDeltas->checkpoint ( "Position buffer" );
- if (first) {
- if (multithread) setPositionsBuffer(false);
- buildTriangulation(pZero,*solver);
- initializeVolumes(*solver);
- backgroundSolver=solver;
- backgroundCompleted=true;
- }
- solver->ompThreads = ompThreads>0? ompThreads : omp_get_max_threads();
-
- timingDeltas->checkpoint ( "Triangulating" );
- updateVolumes ( *solver );
- timingDeltas->checkpoint ( "Update_Volumes" );
-
- epsVolCumulative += epsVolMax;
- retriangulationLastIter++;
- if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
- (defTolerance>0 && epsVolCumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval;
-
- ///compute flow and and forces here
- if (pressureForce){
- solver->gaussSeidel(scene->dt);
- timingDeltas->checkpoint ( "Gauss-Seidel (includes matrix construct and factorization in single-thread mode)" );
- solver->computeFacetForcesWithCache();}
- timingDeltas->checkpoint ( "compute_Forces" );
- ///Application of vicscous forces
- scene->forces.sync();
- timingDeltas->checkpoint ( "forces.sync()" );
- computeViscousForces ( *solver );
- timingDeltas->checkpoint ( "viscous forces" );
- Vector3r force;
- Vector3r torque;
- FiniteVerticesIterator verticesEnd = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
- for ( FiniteVerticesIterator vIt = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); vIt != verticesEnd; vIt++ ) {
- force = pressureForce ? Vector3r ( vIt->info().forces[0],vIt->info().forces[1],vIt->info().forces[2] ): Vector3r(0,0,0);
- torque = Vector3r(0,0,0);
- if (shearLubrication || viscousShear){
- force = force + solver->shearLubricationForces[vIt->info().id()];
- torque = torque + solver->shearLubricationTorques[vIt->info().id()];
- if (pumpTorque)
- torque = torque + solver->pumpLubricationTorques[vIt->info().id()];
- }
- if (twistTorque)
- torque = torque + solver->twistLubricationTorques[vIt->info().id()];
- if (normalLubrication)
- force = force + solver-> normalLubricationForce[vIt->info().id()];
- scene->forces.addForce ( vIt->info().id(), force);
- scene->forces.addTorque ( vIt->info().id(), torque);
- }
- ///End compute flow and forces
- timingDeltas->checkpoint ( "Applying Forces" );
- int sleeping = 0;
- if (multithread && !first) {
- while (updateTriangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/
- sleeping++;
- boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
- if (debug && sleeping) cerr<<"sleeping..."<<sleeping<<endl;
- if (updateTriangulation || (ellapsedIter>(0.5*meshUpdateInterval) && backgroundCompleted)) {
- if (debug) cerr<<"switch flow solver"<<endl;
- if (useSolver==0) LOG_ERROR("background calculations not available for Gauss-Seidel");
- if (fluidBulkModulus>0) solver->interpolate (solver->T[solver->currentTes], backgroundSolver->T[backgroundSolver->currentTes]);
- solver=backgroundSolver;
- backgroundSolver = shared_ptr<FlowSolver> (new FlowSolver);
- if (metisForced) {backgroundSolver->eSolver.cholmod().nmethods=1; backgroundSolver->eSolver.cholmod().method[0].ordering=CHOLMOD_METIS;}
- //Copy imposed pressures/flow from the old solver
- backgroundSolver->imposedP = vector<pair<CGT::Point,Real> >(solver->imposedP);
- backgroundSolver->imposedF = vector<pair<CGT::Point,Real> >(solver->imposedF);
- if (debug) cerr<<"switched"<<endl;
- setPositionsBuffer(false);//set "parallel" buffer for background calculation
- backgroundCompleted=false;
- retriangulationLastIter=ellapsedIter;
- updateTriangulation=false;
- epsVolCumulative=0;
- ellapsedIter=0;
- boost::thread workerThread(&TemplateFlowEngine::backgroundAction,this);
- workerThread.detach();
- if (debug) cerr<<"backgrounded"<<endl;
- initializeVolumes(*solver);
- computeViscousForces(*solver);
- if (debug) cerr<<"volumes initialized"<<endl;
- }
- else {
- if (debug && !backgroundCompleted) cerr<<"still computing solver in the background, ellapsedIter="<<ellapsedIter<<endl;
- ellapsedIter++;
- }
- } else {
- if (updateTriangulation && !first) {
- buildTriangulation (pZero, *solver);
- initializeVolumes(*solver);
- computeViscousForces(*solver);
- updateTriangulation = false;
- epsVolCumulative=0;
- retriangulationLastIter=0;
- ReTrg++;}
- }
- first=false;
- timingDeltas->checkpoint ( "triangulate + init volumes" );
-}
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::backgroundAction()
-{
- if (useSolver<1) {LOG_ERROR("background calculations not available for Gauss-Seidel"); return;}
- buildTriangulation ( pZero,*backgroundSolver );
- //FIXME: GS is computing too much, we need only matrix factorization in fact
- backgroundSolver->gaussSeidel(scene->dt);
- //FIXME(2): and here we need only cached variables, not forces
- backgroundSolver->computeFacetForcesWithCache(/*onlyCache?*/ true);
-// boost::this_thread::sleep(boost::posix_time::seconds(5));
- backgroundCompleted = true;
-}
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::boundaryConditions ( Solver& flow )
-{
- for (int k=0;k<6;k++) {
- flow.boundary (wallIds[k]).flowCondition=!bndCondIsPressure[k];
- flow.boundary (wallIds[k]).value=bndCondValue[k];
- flow.boundary (wallIds[k]).velocity = boundaryVelocity[k];//FIXME: needs correct implementation, maybe update the cached pos/vel?
- }
-}
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setImposedPressure ( unsigned int cond, Real p)
-{
- if ( cond>=solver->imposedP.size() ) LOG_ERROR ( "Setting p with cond higher than imposedP size." );
- solver->imposedP[cond].second=p;
- //force immediate update of boundary conditions
- solver->pressureChanged=true;
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposeFlux ( Vector3r pos, Real flux){
- solver->imposedF.push_back ( pair<CGT::Point,Real> ( CGT::Point ( pos[0],pos[1],pos[2] ), flux ) );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedPressure () { solver->imposedP.clear(); solver->IPCells.clear();}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedFlux () { solver->imposedF.clear(); solver->IFCells.clear();}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::getCellFlux ( unsigned int cond)
-{
- if ( cond>=solver->imposedP.size() ) {LOG_ERROR ( "Getting flux with cond higher than imposedP size." ); return 0;}
- double flux=0;
- typename Solver::CellHandle& cell= solver->IPCells[cond];
- for ( int ngb=0;ngb<4;ngb++ ) {
- /*if (!cell->neighbor(ngb)->info().Pcondition)*/
- flux+= cell->info().kNorm() [ngb]* ( cell->info().p()-cell->neighbor ( ngb )->info().p() );
- }
- return flux+cell->info().dv();
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::initSolver ( FlowSolver& flow )
-{
- flow.Vtotalissimo=0; flow.VSolidTot=0; flow.vPoral=0; flow.sSolidTot=0;
- flow.slipBoundary=slipBoundary;
- flow.kFactor = permeabilityFactor;
- flow.debugOut = debug;
- flow.useSolver = useSolver;
- #ifdef EIGENSPARSE_LIB
- flow.numSolveThreads = numSolveThreads;
- flow.numFactorizeThreads = numFactorizeThreads;
- #endif
- flow.meanKStat = meanKStat;
- flow.viscosity = viscosity;
- flow.tolerance=tolerance;
- flow.relax=relax;
- flow.clampKValues = clampKValues;
- flow.maxKdivKmean = maxKdivKmean;
- flow.minKdivKmean = minKdivKmean;
- flow.meanKStat = meanKStat;
- flow.permeabilityMap = permeabilityMap;
- flow.fluidBulkModulus = fluidBulkModulus;
- flow.T[flow.currentTes].Clear();
- flow.T[flow.currentTes].maxId=-1;
- flow.xMin = 1000.0, flow.xMax = -10000.0, flow.yMin = 1000.0, flow.yMax = -10000.0, flow.zMin = 1000.0, flow.zMax = -10000.0;
-}
-
-#ifdef LINSOLV
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setForceMetis ( bool force )
-{
- if (force) {
- metisForced=true;
- solver->eSolver.cholmod().nmethods=1;
- solver->eSolver.cholmod().method[0].ordering=CHOLMOD_METIS;
- } else {cholmod_defaults(&(solver->eSolver.cholmod())); metisForced=false;}
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-bool TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::getForceMetis () {return (solver->eSolver.cholmod().nmethods==1);}
-#endif
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( Solver& flow )
-{
- buildTriangulation ( 0.f,flow );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( double pZero, Solver& flow )
-{
- flow.resetNetwork();
- if (first) flow.currentTes=0;
- else {
- flow.currentTes=!flow.currentTes;
- if (debug) cout << "--------RETRIANGULATION-----------" << endl;
- }
-
- initSolver(flow);
-
- addBoundary ( flow );
- triangulate ( flow );
- if ( debug ) cout << endl << "Tesselating------" << endl << endl;
- flow.T[flow.currentTes].compute();
-
- flow.defineFictiousCells();
- // For faster loops on cells define this vector
- flow.T[flow.currentTes].cellHandles.clear();
- flow.T[flow.currentTes].cellHandles.reserve(flow.T[flow.currentTes].Triangulation().number_of_finite_cells());
- FiniteCellsIterator cell_end = flow.T[flow.currentTes].Triangulation().finite_cells_end();
- int k=0;
- for ( FiniteCellsIterator cell = flow.T[flow.currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
- flow.T[flow.currentTes].cellHandles.push_back(cell);
- cell->info().id=k++;}//define unique numbering now, corresponds to position in cellHandles
- flow.displayStatistics ();
- flow.computePermeability();
- //This virtual function does nothing yet, derived class may overload it to make permeability different (see DFN engine)
- trickPermeability();
- porosity = flow.vPoralPorosity/flow.vTotalPorosity;
-
- boundaryConditions ( flow );
- flow.initializePressure ( pZero );
-
- if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
- if ( waveAction ) flow.applySinusoidalPressure ( flow.T[flow.currentTes].Triangulation(), sineMagnitude, sineAverage, 30 );
- if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setPositionsBuffer(bool current)
-{
- vector<posData>& buffer = current? positionBufferCurrent : positionBufferParallel;
- buffer.clear();
- buffer.resize(scene->bodies->size());
- shared_ptr<Sphere> sph ( new Sphere );
- const int Sph_Index = sph->getClassIndexStatic();
- FOREACH ( const shared_ptr<Body>& b, *scene->bodies ) {
- if (!b || ignoredBody==b->getId()) continue;
- posData& dat = buffer[b->getId()];
- dat.id=b->getId();
- dat.pos=b->state->pos;
- dat.isSphere= (b->shape->getClassIndex() == Sph_Index);
- if (dat.isSphere) dat.radius = YADE_CAST<Sphere*>(b->shape.get())->radius;
- dat.exists=true;
- }
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::addBoundary ( Solver& flow )
-{
- vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
- solver->xMin = Mathr::MAX_REAL, solver->xMax = -Mathr::MAX_REAL, solver->yMin = Mathr::MAX_REAL, solver->yMax = -Mathr::MAX_REAL, solver->zMin = Mathr::MAX_REAL, solver->zMax = -Mathr::MAX_REAL;
- FOREACH ( const posData& b, buffer ) {
- if ( !b.exists ) continue;
- if ( b.isSphere ) {
- const Real& rad = b.radius;
- const Real& x = b.pos[0];
- const Real& y = b.pos[1];
- const Real& z = b.pos[2];
- flow.xMin = min ( flow.xMin, x-rad );
- flow.xMax = max ( flow.xMax, x+rad );
- flow.yMin = min ( flow.yMin, y-rad );
- flow.yMax = max ( flow.yMax, y+rad );
- flow.zMin = min ( flow.zMin, z-rad );
- flow.zMax = max ( flow.zMax, z+rad );
- }
- }
- //FIXME idOffset must be set correctly, not the case here (always 0), then we need walls first or it will fail
- idOffset = flow.T[flow.currentTes].maxId+1;
- flow.idOffset = idOffset;
- flow.sectionArea = ( flow.xMax - flow.xMin ) * ( flow.zMax-flow.zMin );
- flow.vTotal = ( flow.xMax-flow.xMin ) * ( flow.yMax-flow.yMin ) * ( flow.zMax-flow.zMin );
- flow.yMinId=wallIds[ymin];
- flow.yMaxId=wallIds[ymax];
- flow.xMaxId=wallIds[xmax];
- flow.xMinId=wallIds[xmin];
- flow.zMinId=wallIds[zmin];
- flow.zMaxId=wallIds[zmax];
-
- //FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
- flow.boundsIds[0]= &flow.xMinId;
- flow.boundsIds[1]= &flow.xMaxId;
- flow.boundsIds[2]= &flow.yMinId;
- flow.boundsIds[3]= &flow.yMaxId;
- flow.boundsIds[4]= &flow.zMinId;
- flow.boundsIds[5]= &flow.zMaxId;
-
- for (int k=0;k<6;k++) flow.boundary ( *flow.boundsIds[k] ).useMaxMin = boundaryUseMaxMin[k];
-
- flow.cornerMin = CGT::Point ( flow.xMin, flow.yMin, flow.zMin );
- flow.cornerMax = CGT::Point ( flow.xMax, flow.yMax, flow.zMax );
-
- //assign BCs types and values
- boundaryConditions ( flow );
-
- double center[3];
- for ( int i=0; i<6; i++ ) {
- if ( *flow.boundsIds[i]<0 ) continue;
- CGT::CVector Normal ( normal[i].x(), normal[i].y(), normal[i].z() );
- if ( flow.boundary ( *flow.boundsIds[i] ).useMaxMin ) flow.addBoundingPlane(Normal, *flow.boundsIds[i] );
- else {
- for ( int h=0;h<3;h++ ) center[h] = buffer[*flow.boundsIds[i]].pos[h];
-// cerr << "id="<<*flow.boundsIds[i] <<" center="<<center[0]<<","<<center[1]<<","<<center[2]<<endl;
- flow.addBoundingPlane ( center, wallThickness, Normal,*flow.boundsIds[i] );
- }
- }
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::triangulate ( Solver& flow )
-{
-///Using Tesselation wrapper (faster)
-// TesselationWrapper TW;
-// if (TW.Tes) delete TW.Tes;
-// TW.Tes = &(flow.T[flow.currentTes]);//point to the current Tes we have in Flowengine
-// TW.insertSceneSpheres();//TW is now really inserting in TemplateFlowEngine, using the faster insert(begin,end)
-// TW.Tes = NULL;//otherwise, Tes would be deleted by ~TesselationWrapper() at the end of the function.
-///Using one-by-one insertion
- vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
- FOREACH ( const posData& b, buffer ) {
- if ( !b.exists ) continue;
- if ( b.isSphere ) {
- if (b.id==ignoredBody) continue;
- flow.T[flow.currentTes].insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
- }
- flow.T[flow.currentTes].redirected=true;//By inserting one-by-one, we already redirected
- flow.shearLubricationForces.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.shearLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.pumpLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.twistLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.shearLubricationBodyStress.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.normalLubricationForce.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.normalLubricationBodyStress.resize ( flow.T[flow.currentTes].maxId+1 );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::initializeVolumes ( Solver& flow )
-{
- typedef typename Solver::FiniteVerticesIterator FiniteVerticesIterator;
-
- FiniteVerticesIterator vertices_end = flow.T[flow.currentTes].Triangulation().finite_vertices_end();
- CGT::CVector Zero(0,0,0);
- for (FiniteVerticesIterator V_it = flow.T[flow.currentTes].Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero;
-
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles)
- {
- switch ( cell->info().fictious() )
- {
- case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
- case ( 1 ) : cell->info().volume() = volumeCellSingleFictious ( cell ); break;
- case ( 2 ) : cell->info().volume() = volumeCellDoubleFictious ( cell ); break;
- case ( 3 ) : cell->info().volume() = volumeCellTripleFictious ( cell ); break;
- default: break;
- }
- if (flow.fluidBulkModulus>0) { cell->info().invVoidVolume() = 1 / ( abs(cell->info().volume()) - flow.volumeSolidPore(cell) ); }
- }
- if (debug) cout << "Volumes initialised." << endl;
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::averageRealCellVelocity()
-{
- solver->averageRelativeCellVelocity();
- Vector3r Vel ( 0,0,0 );
- //AVERAGE CELL VELOCITY
- FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
- for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
- for ( int g=0;g<4;g++ ) {
- if ( !cell->vertex ( g )->info().isFictious ) {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
- for ( int i=0;i<3;i++ ) Vel[i] = Vel[i] + sph->state->vel[i]/4;
- }
- }
- RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
- CGT::Point pos_av_facet;
- double volume_facet_translation = 0;
- CGT::CVector Vel_av ( Vel[0], Vel[1], Vel[2] );
- for ( int i=0; i<4; i++ ) {
- volume_facet_translation = 0;
- if ( !Tri.is_infinite ( cell->neighbor ( i ) ) ) {
- CGT::CVector Surfk = cell->info()-cell->neighbor ( i )->info();
- Real area = sqrt ( Surfk.squared_length() );
- Surfk = Surfk/area;
- CGT::CVector branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
- pos_av_facet = ( CGT::Point ) cell->info() + ( branch*Surfk ) *Surfk;
- volume_facet_translation += Vel_av*cell->info().facetSurfaces[i];
- cell->info().averageVelocity() = cell->info().averageVelocity() - volume_facet_translation/cell->info().volume() * ( pos_av_facet-CGAL::ORIGIN );
- }
- }
- }
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::updateVolumes ( Solver& flow )
-{
- if ( debug ) cout << "Updating volumes.............." << endl;
- Real invDeltaT = 1/scene->dt;
- epsVolMax=0;
- Real totVol=0; Real totDVol=0;
- #ifdef YADE_OPENMP
- const long size=flow.T[flow.currentTes].cellHandles.size();
- #pragma omp parallel for num_threads(ompThreads>0 ? ompThreads : 1)
- for(long i=0; i<size; i++){
- CellHandle& cell = flow.T[flow.currentTes].cellHandles[i];
- #else
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
- #endif
- double newVol, dVol;
- switch ( cell->info().fictious() ) {
- case ( 3 ) : newVol = volumeCellTripleFictious ( cell ); break;
- case ( 2 ) : newVol = volumeCellDoubleFictious ( cell ); break;
- case ( 1 ) : newVol = volumeCellSingleFictious ( cell ); break;
- case ( 0 ) : newVol = volumeCell (cell ); break;
- default: newVol = 0; break;}
- dVol=cell->info().volumeSign* ( newVol - cell->info().volume() );
- cell->info().dv() = dVol*invDeltaT;
- cell->info().volume() = newVol;
- if (defTolerance>0) { //if the criterion is not used, then we skip these updates a save a LOT of time when Nthreads > 1
- #pragma omp atomic
- totVol+=newVol;
- #pragma omp atomic
- totDVol+=abs(dVol);}
- }
- if (defTolerance>0) epsVolMax = totDVol/totVol;
- for (unsigned int n=0; n<flow.imposedF.size();n++) {
- flow.IFCells[n]->info().dv()+=flow.imposedF[n].second;
- flow.IFCells[n]->info().Pcondition=false;}
- if ( debug ) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<endl;
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-template<class Cellhandle>
-Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellSingleFictious ( Cellhandle cell )
-{
- Vector3r V[3];
- int b=0;
- int w=0;
- cell->info().volumeSign=1;
- Real Wall_coordinate=0;
-
- for ( int y=0;y<4;y++ ) {
- if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
- V[w]=positionBufferCurrent[cell->vertex ( y )->info().id()].pos;
- w++;
- } else {
- b = cell->vertex ( y )->info().id();
- const shared_ptr<Body>& wll = Body::byId ( b , scene );
- if ( !solver->boundary ( b ).useMaxMin ) Wall_coordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wallThickness/2.;
- else Wall_coordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate];
- }
- }
- Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
- return abs ( Volume );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-template<class Cellhandle>
-Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellDoubleFictious ( Cellhandle cell )
-{
- Vector3r A=Vector3r::Zero(), AS=Vector3r::Zero(),B=Vector3r::Zero(), BS=Vector3r::Zero();
-
- cell->info().volumeSign=1;
- int b[2];
- int coord[2];
- Real Wall_coordinate[2];
- int j=0;
- bool first_sph=true;
-
- for ( int g=0;g<4;g++ ) {
- if ( cell->vertex ( g )->info().isFictious ) {
- b[j] = cell->vertex ( g )->info().id();
- coord[j]=solver->boundary ( b[j] ).coordinate;
- if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = positionBufferCurrent[b[j]].pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2.;
- else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
- j++;
- } else if ( first_sph ) {
- A=AS=/*AT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;
- first_sph=false;
- } else {
- B=BS=/*BT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;;
- }
- }
- AS[coord[0]]=BS[coord[0]] = Wall_coordinate[0];
-
- //first pyramid with triangular base (A,B,BS)
- Real Vol1=0.5* ( ( A-BS ).cross ( B-BS ) ) [coord[1]]* ( 0.333333333* ( 2*B[coord[1]]+A[coord[1]] )-Wall_coordinate[1] );
- //second pyramid with triangular base (A,AS,BS)
- Real Vol2=0.5* ( ( AS-BS ).cross ( A-BS ) ) [coord[1]]* ( 0.333333333* ( B[coord[1]]+2*A[coord[1]] )-Wall_coordinate[1] );
- return abs ( Vol1+Vol2 );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-template<class Cellhandle>
-Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellTripleFictious ( Cellhandle cell )
-{
- Vector3r A;
-
- int b[3];
- int coord[3];
- Real Wall_coordinate[3];
- int j=0;
- cell->info().volumeSign=1;
-
- for ( int g=0;g<4;g++ ) {
- if ( cell->vertex ( g )->info().isFictious ) {
- b[j] = cell->vertex ( g )->info().id();
- coord[j]=solver->boundary ( b[j] ).coordinate;
- const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
- if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = wll->state->pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2.;
- else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
- j++;
- } else {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
- A= ( sph->state->pos );
- }
- }
- Real Volume = ( A[coord[0]]-Wall_coordinate[0] ) * ( A[coord[1]]-Wall_coordinate[1] ) * ( A[coord[2]]-Wall_coordinate[2] );
- return abs ( Volume );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-template<class Cellhandle>
-Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCell ( Cellhandle cell )
-{
- static const Real inv6 = 1/6.;
- const Vector3r& p0 = positionBufferCurrent[cell->vertex ( 0 )->info().id()].pos;
- const Vector3r& p1 = positionBufferCurrent[cell->vertex ( 1 )->info().id()].pos;
- const Vector3r& p2 = positionBufferCurrent[cell->vertex ( 2 )->info().id()].pos;
- const Vector3r& p3 = positionBufferCurrent[cell->vertex ( 3 )->info().id()].pos;
- Real volume = inv6 * ((p0-p1).cross(p0-p2)).dot(p0-p3);
- if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
- return volume;
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::computeViscousForces ( Solver& flow )
-{
- if (normalLubrication || shearLubrication || viscousShear){
- if ( debug ) cout << "Application of viscous forces" << endl;
- if ( debug ) cout << "Number of edges = " << flow.edgeIds.size() << endl;
- for ( unsigned int k=0; k<flow.shearLubricationForces.size(); k++ ) flow.shearLubricationForces[k]=Vector3r::Zero();
- for ( unsigned int k=0; k<flow.shearLubricationTorques.size(); k++ ) flow.shearLubricationTorques[k]=Vector3r::Zero();
- for ( unsigned int k=0; k<flow.pumpLubricationTorques.size(); k++ ) flow.pumpLubricationTorques[k]=Vector3r::Zero();
- for ( unsigned int k=0; k<flow.twistLubricationTorques.size(); k++ ) flow.twistLubricationTorques[k]=Vector3r::Zero();
- for ( unsigned int k=0; k<flow.shearLubricationBodyStress.size(); k++) flow.shearLubricationBodyStress[k]=Matrix3r::Zero();
- for ( unsigned int k=0; k<flow.normalLubricationForce.size(); k++ ) flow.normalLubricationForce[k]=Vector3r::Zero();
- for ( unsigned int k=0; k<flow.normalLubricationBodyStress.size(); k++) flow.normalLubricationBodyStress[k]=Matrix3r::Zero();
-
- typedef typename Solver::Tesselation Tesselation;
- const Tesselation& Tes = flow.T[flow.currentTes];
- flow.deltaShearVel.clear(); flow.normalV.clear(); flow.deltaNormVel.clear(); flow.surfaceDistance.clear(); flow.onlySpheresInteractions.clear(); flow.normalStressInteraction.clear(); flow.shearStressInteraction.clear();
-
-
- for ( int i=0; i< ( int ) flow.edgeIds.size(); i++ ) {
- const int& id1 = flow.edgeIds[i].first;
- const int& id2 = flow.edgeIds[i].second;
-
- int hasFictious= Tes.vertex ( id1 )->info().isFictious + Tes.vertex ( id2 )->info().isFictious;
- if (hasFictious>0 or id1==id2) continue;
- const shared_ptr<Body>& sph1 = Body::byId ( id1, scene );
- const shared_ptr<Body>& sph2 = Body::byId ( id2, scene );
- Sphere* s1=YADE_CAST<Sphere*> ( sph1->shape.get() );
- Sphere* s2=YADE_CAST<Sphere*> ( sph2->shape.get() );
- const Real& r1 = s1->radius;
- const Real& r2 = s2->radius;
- Vector3r deltaV; Real deltaNormV; Vector3r deltaShearV;
- Vector3r O1O2Vector; Real O1O2; Vector3r normal; Real surfaceDist; Vector3r O1CVector; Vector3r O2CVector;Real meanRad ;Real Rh; Vector3r deltaAngVel; Vector3r deltaShearAngVel;
- Vector3r shearLubF; Vector3r normaLubF; Vector3r pumpT; Vector3r deltaAngNormVel; Vector3r twistT; Vector3r angVel1; Vector3r angVel2;
- //FIXME: if periodic and velGrad!=0, then deltaV should account for velGrad, not the case currently
- if ( !hasFictious ){
- O1O2Vector = sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()) - sph1->state->pos - makeVector3r(Tes.vertex(id1)->info().ghostShift());
- O1O2 = O1O2Vector.norm();
- normal= (O1O2Vector/O1O2);
- surfaceDist = O1O2 - r2 - r1;
- O1CVector = (O1O2/2. + (pow(r1,2) - pow(r2,2)) / (2.*O1O2))*normal;
- O2CVector = -(O1O2Vector - O1CVector);
- meanRad = (r2 + r1)/2.;
- Rh = (r1 < r2)? surfaceDist + 0.45 * r1 : surfaceDist + 0.45 * r2;
- deltaV = (sph2->state->vel + sph2->state->angVel.cross(-r2 * normal)) - (sph1->state->vel+ sph1->state->angVel.cross(r1 * normal));
- angVel1 = sph1->state->angVel;
- angVel2 = sph2->state->angVel;
- deltaAngVel = sph2->state->angVel - sph1->state->angVel;
-
- } else {
- if ( hasFictious==1 ) {//for the fictious sphere, use velocity of the boundary, not of the body
- bool v1fictious = Tes.vertex ( id1 )->info().isFictious;
- int bnd = v1fictious? id1 : id2;
- int coord = flow.boundary(bnd).coordinate;
- O1O2 = v1fictious ? abs((sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]) : abs((sph1->state->pos + makeVector3r(Tes.vertex(id1)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]);
- if(v1fictious)
- normal = makeVector3r(flow.boundary(id1).normal);
- else
- normal = -makeVector3r(flow.boundary(id2).normal);
- O1O2Vector = O1O2 * normal;
- meanRad = v1fictious ? r2:r1;
- surfaceDist = O1O2 - meanRad;
- if (v1fictious){
- O1CVector = Vector3r::Zero();
- O2CVector = - O1O2Vector;}
- else{
- O1CVector = O1O2Vector;
- O2CVector = Vector3r::Zero();}
-
- Rh = surfaceDist + 0.45 * meanRad;
- Vector3r v1 = ( Tes.vertex ( id1 )->info().isFictious ) ? flow.boundary ( id1 ).velocity:sph1->state->vel + sph1->state->angVel.cross(r1 * normal);
- Vector3r v2 = ( Tes.vertex ( id2 )->info().isFictious ) ? flow.boundary ( id2 ).velocity:sph2->state->vel + sph2->state->angVel.cross(-r2 * (normal));
- deltaV = v2-v1;
- angVel1 = ( Tes.vertex ( id1 )->info().isFictious ) ? Vector3r::Zero() : sph1->state->angVel;
- angVel2 = ( Tes.vertex ( id2 )->info().isFictious ) ? Vector3r::Zero() : sph2->state->angVel;
- deltaAngVel = angVel2 - angVel1;
- }
- }
- deltaShearV = deltaV - ( normal.dot ( deltaV ) ) *normal;
- deltaShearAngVel = deltaAngVel - ( normal.dot ( deltaAngVel ) ) *normal;
- flow.deltaShearVel.push_back(deltaShearV);
- flow.normalV.push_back(normal);
- flow.surfaceDistance.push_back(max(surfaceDist, 0.) + eps*meanRad);
-
- /// Compute the shear Lubrication force and torque on each particle
-
- if (shearLubrication)
- shearLubF = flow.computeShearLubricationForce(deltaShearV,surfaceDist,i,eps,O1O2,meanRad);
- else if (viscousShear)
- shearLubF = flow.computeViscousShearForce ( deltaShearV, i , Rh);
-
- if (viscousShear || shearLubrication){
-
- flow.shearLubricationForces[id1]+=shearLubF;
- flow.shearLubricationForces[id2]+=(-shearLubF);
- flow.shearLubricationTorques[id1]+=O1CVector.cross(shearLubF);
- flow.shearLubricationTorques[id2]+=O2CVector.cross(-shearLubF);
-
- /// Compute the pump Lubrication torque on each particle
-
- if (pumpTorque){
- pumpT = flow.computePumpTorque(deltaShearAngVel, surfaceDist, i, eps, meanRad );
- flow.pumpLubricationTorques[id1]+=(-pumpT);
- flow.pumpLubricationTorques[id2]+=pumpT;}
-
- /// Compute the twist Lubrication torque on each particle
-
- if (twistTorque){
- deltaAngNormVel = (normal.dot(deltaAngVel))*normal ;
- twistT = flow.computeTwistTorque(deltaAngNormVel, surfaceDist, i, eps, meanRad );
- flow.twistLubricationTorques[id1]+=(-twistT);
- flow.twistLubricationTorques[id2]+=twistT;
- }
- }
- /// Compute the viscous shear stress on each particle
-
- if (viscousShearBodyStress){
- flow.shearLubricationBodyStress[id1] += shearLubF * O1CVector.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
- flow.shearLubricationBodyStress[id2] += (-shearLubF) * O2CVector.transpose()/ (4.0/3.0 *3.14* pow(r2,3));
- flow.shearStressInteraction.push_back(shearLubF * O1O2Vector.transpose()/(4.0/3.0 *3.14* pow(r1,3)));
- }
-
- /// Compute the normal lubrication force applied on each particle
-
- if (normalLubrication){
- deltaNormV = normal.dot(deltaV);
- flow.deltaNormVel.push_back(deltaNormV * normal);
- normaLubF = flow.computeNormalLubricationForce (deltaNormV, surfaceDist, i,eps,stiffness,scene->dt,meanRad)*normal;
- flow.normalLubricationForce[id1]+=normaLubF;
- flow.normalLubricationForce[id2]+=(-normaLubF);
-
- /// Compute the normal lubrication stress on each particle
-
- if (viscousNormalBodyStress){
- flow.normalLubricationBodyStress[id1] += normaLubF * O1CVector.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
- flow.normalLubricationBodyStress[id2] += (-normaLubF) *O2CVector.transpose() / (4.0/3.0 *3.14* pow(r2,3));
- flow.normalStressInteraction.push_back(normaLubF * O1O2Vector.transpose()/(4.0/3.0 *3.14* pow(r1,3)));
- }
- }
-
- if (!hasFictious)
- flow.onlySpheresInteractions.push_back(i);
-
- }
- }
-}
-
-#endif //FLOW_ENGINE
-
-#endif /* YADE_CGAL */
-
=== removed file 'pkg/dem/PeriodicFlowEngine.cpp'
--- pkg/dem/PeriodicFlowEngine.cpp 2014-04-01 13:18:38 +0000
+++ pkg/dem/PeriodicFlowEngine.cpp 1970-01-01 00:00:00 +0000
@@ -1,452 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* Copyright (C) 2012 by Donia Marzougui <donia.marzougui@xxxxxxxxxxxxxxx>*
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#ifdef YADE_CGAL
-#ifdef FLOW_ENGINE
-
-#define TEMPLATE_FLOW_NAME FlowEngine_PeriodicInfo
-#include "PeriodicFlowEngine.hpp"
-#undef TEMPLATE_FLOW_NAME
-
-CVector PeriodicCellInfo::hSize[]={CVector(),CVector(),CVector()};
-CVector PeriodicCellInfo::deltaP=CVector();
-CVector PeriodicCellInfo::gradP=CVector();
-
-CREATE_LOGGER ( PeriodicFlowEngine );
-
-PeriodicFlowEngine::~PeriodicFlowEngine(){}
-
-void PeriodicFlowEngine:: action()
-{
- if ( !isActivated ) return;
- timingDeltas->start();
- preparePShifts();
- setPositionsBuffer(true);
- if (first) {
- if (multithread) setPositionsBuffer(false);
- cachedCell= Cell(*(scene->cell));
- buildTriangulation(pZero,*solver);
- if (solver->errorCode>0) {LOG_INFO("triangulation error, pausing"); Omega::instance().pause(); return;}
- initializeVolumes(*solver); backgroundSolver=solver; backgroundCompleted=true;}
-// if ( first ) {buildTriangulation ( pZero ); updateTriangulation = false; initializeVolumes();}
-
- timingDeltas->checkpoint("Triangulating");
- updateVolumes (*solver);
- epsVolCumulative += epsVolMax;
- retriangulationLastIter++;
- if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
- (defTolerance>0 && epsVolCumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval;
-
- timingDeltas->checkpoint("Update_Volumes");
-
- ///compute flow and and forces here
- if (pressureForce){
- solver->gaussSeidel(scene->dt);
- timingDeltas->checkpoint("Gauss-Seidel");
- solver->computeFacetForcesWithCache();}
- timingDeltas->checkpoint("compute_Pressure_Forces");
-
- ///compute vicscous forces
- scene->forces.sync();
- computeViscousForces(*solver);
- timingDeltas->checkpoint("compute_Viscous_Forces");
- Vector3r force;
- Vector3r torque;
- const Tesselation& Tes = solver->T[solver->currentTes];
- for (int id=0; id<=Tes.maxId; id++){
- assert (Tes.vertexHandles[id] != NULL);
- const Tesselation::VertexInfo& v_info = Tes.vertexHandles[id]->info();
- force =(pressureForce) ? Vector3r ( ( v_info.forces ) [0],v_info.forces[1],v_info.forces[2] ) : Vector3r(0,0,0);
- torque = Vector3r(0,0,0);
- if (shearLubrication || viscousShear){
- force = force +solver->shearLubricationForces[v_info.id()];
- torque = torque +solver->shearLubricationTorques[v_info.id()];
- if (pumpTorque)
- torque = torque +solver->pumpLubricationTorques[v_info.id()];
- if (twistTorque)
- torque = torque +solver->twistLubricationTorques[v_info.id()];
- }
-
- if (normalLubrication)
- force = force + solver->normalLubricationForce[v_info.id()];
- scene->forces.addForce ( v_info.id(), force);
- scene->forces.addTorque ( v_info.id(), torque);
- }
- ///End Compute flow and forces
- timingDeltas->checkpoint("Applying Forces");
- if (multithread && !first) {
- while (updateTriangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/ boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
- if (updateTriangulation || ellapsedIter>(0.5*meshUpdateInterval)) {
- if (useSolver==0) LOG_ERROR("background calculations not available for Gauss-Seidel");
- if (fluidBulkModulus>0) solver->interpolate (solver->T[solver->currentTes], backgroundSolver->T[backgroundSolver->currentTes]);
- solver=backgroundSolver;
- backgroundSolver = shared_ptr<FlowSolver> (new FlowSolver);
- //Copy imposed pressures/flow from the old solver
- backgroundSolver->imposedP = vector<pair<CGT::Point,Real> >(solver->imposedP);
- backgroundSolver->imposedF = vector<pair<CGT::Point,Real> >(solver->imposedF);
- setPositionsBuffer(false);
- cachedCell= Cell(*(scene->cell));
- backgroundCompleted=false;
- retriangulationLastIter=ellapsedIter;
- ellapsedIter=0;
- epsVolCumulative=0;
- boost::thread workerThread(&PeriodicFlowEngine::backgroundAction,this);
- workerThread.detach();
- initializeVolumes(*solver);
- computeViscousForces(*solver);
- }
- else if (debug && !first) {
- if (debug && !backgroundCompleted) cerr<<"still computing solver in the background"<<endl;
- ellapsedIter++;}
- } else {
- if (updateTriangulation && !first) {
- cachedCell= Cell(*(scene->cell));
- buildTriangulation (pZero, *solver);
- initializeVolumes(*solver);
- computeViscousForces(*solver);
- updateTriangulation = false;
- epsVolCumulative=0;
- retriangulationLastIter=0;
- ReTrg++;}
- }
- first=false;
- timingDeltas->checkpoint("Ending");
-}
-
-
-// void PeriodicFlowEngine::backgroundAction()
-// {
-// if (useSolver<1) {LOG_ERROR("background calculations not available for Gauss-Seidel"); return;}
-// buildTriangulation (pZero,*backgroundSolver);
-// //FIXME: GS is computing too much, we need only matrix factorization in fact
-// backgroundSolver->gaussSeidel(scene->dt);
-// backgroundSolver->computeFacetForcesWithCache(/*onlyCache?*/ true);
-// // boost::this_thread::sleep(boost::posix_time::seconds(10));
-// backgroundCompleted = true;
-// }
-
-void PeriodicFlowEngine::triangulate( FlowSolver& flow )
-{
- Tesselation& Tes = flow.T[flow.currentTes];
- vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
- FOREACH ( const posData& b, buffer ) {
- if ( !b.exists || !b.isSphere || b.id==ignoredBody) continue;
- Vector3i period; Vector3r wpos;
- // FIXME: use "sheared" variant if the cell is sheared
- wpos=cachedCell.wrapPt ( b.pos,period );
- const Body::id_t& id = b.id;
- const Real& rad = b.radius;
- const Real& x = wpos[0];
- const Real& y = wpos[1];
- const Real& z = wpos[2];
- VertexHandle vh0=Tes.insert ( x, y, z, rad, id );
-// VertexHandle vh0=Tes.insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );
- if (vh0==NULL) {
- flow.errorCode = 2;
- LOG_ERROR("Vh NULL in PeriodicFlowEngine::triangulate(), check input data"); continue;}
- for ( int k=0;k<3;k++ ) vh0->info().period[k]=-period[k];
- const Vector3r cellSize ( cachedCell.getSize() );
- //FIXME: if hasShear, comment in
-// wpos=scene->cell->unshearPt ( wpos );
- // traverse all periodic cells around the body, to see if any of them touches
- Vector3r halfSize= ( rad+duplicateThreshold ) *Vector3r ( 1,1,1 );
- Vector3r pmin,pmax;
- Vector3i i;
- for ( i[0]=-1; i[0]<=1; i[0]++ )
- for ( i[1]=-1;i[1]<=1; i[1]++ )
- for ( i[2]=-1; i[2]<=1; i[2]++ ) {
- if ( i[0]!=0 || i[1]!=0 || i[2]!=0 ) { // middle; already rendered above
- Vector3r pos2=wpos+Vector3r ( cellSize[0]*i[0],cellSize[1]*i[1],cellSize[2]*i[2] ); // shift, but without shear!
- pmin=pos2-halfSize;
- pmax=pos2+halfSize;
- if ( (pmin[0]<=cellSize[0]) && (pmax[0]>=0) && (pmin[1]<=cellSize[1]) && (pmax[1]>=0) && (pmin[2]<=cellSize[2]) && (pmax[2]>=0) ) {
- //with shear:
- //Vector3r pt=scene->cell->shearPt ( pos2 );
- //without shear:
- const Vector3r& pt= pos2;
- VertexHandle vh=Tes.insert ( pt[0],pt[1],pt[2],rad,id,false,id );
- for ( int k=0;k<3;k++ ) vh->info().period[k]=i[k]-period[k];}}
- }
- //re-assign the original vertex pointer since duplicates may have overwrite it
- Tes.vertexHandles[id]=vh0;
- }
- Tes.redirected=true;//By inserting one-by-one, we already redirected
- flow.shearLubricationForces.resize ( Tes.maxId+1 );
- flow.shearLubricationTorques.resize ( Tes.maxId+1 );
- flow.pumpLubricationTorques.resize ( Tes.maxId+1 );
- flow.twistLubricationTorques.resize ( Tes.maxId+1 );
- flow.shearLubricationBodyStress.resize ( Tes.maxId+1 );
- flow.normalLubricationForce.resize ( Tes.maxId+1 );
- flow.normalLubricationBodyStress.resize ( Tes.maxId+1 );
-}
-
-
-Real PeriodicFlowEngine::volumeCell ( CellHandle cell )
-{
- static const Real inv6 = 1/6.;
- const Vector3r p0 = positionBufferCurrent[cell->vertex(0)->info().id()].pos + makeVector3r(cell->vertex(0)->info().ghostShift());
- const Vector3r p1 = positionBufferCurrent[cell->vertex(1)->info().id()].pos + makeVector3r(cell->vertex(1)->info().ghostShift());
- const Vector3r p2 = positionBufferCurrent[cell->vertex(2)->info().id()].pos + makeVector3r(cell->vertex(2)->info().ghostShift());
- const Vector3r p3 = positionBufferCurrent[cell->vertex(3)->info().id()].pos + makeVector3r(cell->vertex(3)->info().ghostShift());
- Real volume = inv6*((p0-p1).cross(p0-p2)).dot(p0-p3);
- if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
- return volume;
-}
-
-Real PeriodicFlowEngine::volumeCellSingleFictious ( CellHandle cell )
-{
- Vector3r V[3];
- int b=0;
- int w=0;
- cell->info().volumeSign=1;
- Real Wall_coordinate=0;
-
- for ( int y=0;y<4;y++ ) {
- if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( y )->info().id(), scene );
- V[w]=sph->state->pos+ makeVector3r ( cell->vertex ( y )->info().ghostShift() );
- w++;
- } else {
- b = cell->vertex ( y )->info().id();
- const shared_ptr<Body>& wll = Body::byId ( b,scene );
- if ( !solver->boundary ( b ).useMaxMin ) Wall_coordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wallThickness/2.;
- else Wall_coordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate];
- }
- }
- Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
- return abs ( Volume );
-}
-
-
-void PeriodicFlowEngine::locateCell ( CellHandle baseCell, unsigned int& index, int& baseIndex, FlowSolver& flow, unsigned int count)
-{
- if (count>10) {
- LOG_ERROR("More than 10 attempts to locate a cell, duplicateThreshold may be too small, resulting in periodicity inconsistencies.");
- flow.errorCode=1; return;
- }
- PeriFlowTesselation::CellInfo& baseInfo = baseCell->info();
- //already located, return FIXME: is inline working correctly? else move this test outside the function, just before the calls
- if ( baseInfo.index>0 || baseInfo.isGhost ) return;
- RTriangulation& Tri = flow.T[flow.currentTes].Triangulation();
- Vector3r center ( 0,0,0 );
- Vector3i period;
-
- if (baseCell->info().fictious()==0)
- for ( int k=0;k<4;k++ ) center+= 0.25*makeVector3r (baseCell->vertex(k)->point());
- else {
-
- Real boundPos=0; int coord=0;
- for ( int k=0;k<4;k++ ) {
- if ( !baseCell->vertex ( k )->info().isFictious ) center+= 0.3333333333*makeVector3r ( baseCell->vertex ( k )->point() );
- else {
- coord=flow.boundary ( baseCell->vertex ( k )->info().id() ).coordinate;
- boundPos=flow.boundary ( baseCell->vertex ( k )->info().id() ).p[coord];}
- }
- center[coord]=boundPos;
- }
- Vector3r wdCenter= cachedCell.wrapPt ( center,period );
- if ( period[0]!=0 || period[1]!=0 || period[2]!=0 ) {
- if ( baseCell->info().index>0 ) {
- cout<<"indexed cell is found ghost!"<<baseInfo.index <<endl;
- baseInfo.isGhost=false;
- return;
- }
- CellHandle ch= Tri.locate ( CGT::Point ( wdCenter[0],wdCenter[1],wdCenter[2] )
-// ,/*hint*/ v0
- );
- baseInfo.period[0]=period[0];
- baseInfo.period[1]=period[1];
- baseInfo.period[2]=period[2];
- //call recursively, since the returned cell could be also a ghost (especially if baseCell is a non-periodic type from the external contour
- locateCell ( ch,index,baseIndex,flow,++count );
- if ( ch==baseCell ) cerr<<"WTF!!"<<endl;
- //check consistency
- bool checkC=false;
- for (int kk=0; kk<4;kk++) if ((!baseCell->vertex(kk)->info().isGhost) && ((!baseCell->vertex(kk)->info().isFictious))) checkC = true;
- if (checkC) {
- bool checkV=true;
- for (int kk=0; kk<4;kk++) {
- checkV=false;
- for (int jj=0; jj<4;jj++)
- if (baseCell->vertex(kk)->info().id() == ch->vertex(jj)->info().id()) checkV = true;
- if (!checkV) {cerr <<"periodicity is broken"<<endl;
- for (int jj=0; jj<4;jj++) cerr<<baseCell->vertex(jj)->info().id()<<" ";
- cerr<<" vs. ";
- for (int jj=0; jj<4;jj++) cerr<<ch->vertex(jj)->info().id()<<" ";
- cerr<<endl;}
- }
- } else {
-// bool checkV=true;
-// for (int kk=0; kk<4;kk++) {
-// checkV=false;
-// for (int jj=0; jj<4;jj++)
-// if (baseCell->vertex(kk)->info().id() == ch->vertex(jj)->info().id()) checkV = true;
-// if (!checkV) {cerr <<"periodicity is broken (that's ok probably)"<<endl;
-// for (int jj=0; jj<4;jj++) cerr<<baseCell->vertex(jj)->info().id()<<" ";
-// cerr<<" vs. ";
-// for (int jj=0; jj<4;jj++) cerr<<ch->vertex(jj)->info().id()<<" ";
-// cerr<<endl;}
-// }
- }
-
- baseInfo.isGhost=true;
- baseInfo._pression=& ( ch->info().p() );
- baseInfo.index=ch->info().index;
- baseInfo.baseIndex=ch->info().baseIndex;
- baseInfo.Pcondition=ch->info().Pcondition;
- } else {
- baseInfo.isGhost=false;
- //index is 1-based, if it is zero it is not initialized, we define it here
- if ( baseInfo.baseIndex<0 ){
- baseInfo.baseIndex=++baseIndex;
- if (!baseInfo.Pcondition) baseInfo.index=++index;}
- }
-}
-
-Vector3r PeriodicFlowEngine::meanVelocity()
-{
- solver->averageRelativeCellVelocity();
- Vector3r meanVel ( 0,0,0 );
- Real volume=0;
- FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
- for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
- //We could also define velocity using cell's center
-// if ( !cell->info().isReal() ) continue;
- if ( cell->info().isGhost ) continue;
- for ( int i=0;i<3;i++ )
- meanVel[i]=meanVel[i]+ ( ( cell->info().averageVelocity() ) [i] * abs ( cell->info().volume() ) );
- volume+=abs ( cell->info().volume() );
- }
- return ( meanVel/volume );
-}
-
-void PeriodicFlowEngine::updateVolumes (FlowSolver& flow)
-{
- if ( debug ) cout << "Updating volumes.............." << endl;
- Real invDeltaT = 1/scene->dt;
- double newVol, dVol;
- epsVolMax=0;
- Real totVol=0;
- Real totDVol=0;
- Real totVol0=0;
- Real totVol1=0;
-
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
- switch ( cell->info().fictious() ) {
- case ( 1 ) :
- newVol = volumeCellSingleFictious ( cell );
- totVol1+=newVol;
- break;
- case ( 0 ) :
- newVol = volumeCell ( cell );
- totVol0+=newVol;
- break;
- default:
- newVol = 0;
- break;
- }
- totVol+=newVol;
- dVol=cell->info().volumeSign * ( newVol - cell->info().volume() );
- totDVol+=dVol;
- epsVolMax = max ( epsVolMax, abs ( dVol/newVol ) );
- cell->info().dv() = dVol * invDeltaT;
- cell->info().volume() = newVol;
- }
- if ( debug ) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<" "<< totVol0<<" "<< totVol1<<endl;
-}
-
-
-void PeriodicFlowEngine::initializeVolumes (FlowSolver& flow)
-{
- FiniteVerticesIterator vertices_end = flow.T[flow.currentTes].Triangulation().finite_vertices_end();
- CGT::CVector Zero ( 0,0,0 );
- for ( FiniteVerticesIterator V_it = flow.T[flow.currentTes].Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++ ) V_it->info().forces=Zero;
-
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
- switch ( cell->info().fictious() )
- {
- case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
- case ( 1 ) : cell->info().volume() = volumeCellSingleFictious ( cell ); break;
- default: cell->info().volume() = 0; break;
- }
- //FIXME: the void volume is negative sometimes, hence crashing...
- if (flow.fluidBulkModulus>0) { cell->info().invVoidVolume() = 1. / (max(0.1*cell->info().volume(),abs(cell->info().volume()) - flow.volumeSolidPore(cell)) ); }
- }
- if ( debug ) cout << "Volumes initialised." << endl;
-}
-
-void PeriodicFlowEngine::buildTriangulation ( double pZero, FlowSolver& flow)
-{
- flow.resetNetwork();
- if (first) flow.currentTes=0;
- else {
- flow.currentTes=!flow.currentTes;
- if ( debug ) cout << "--------RETRIANGULATION-----------" << endl;}
- initSolver(flow);
- addBoundary ( flow );
- if ( debug ) cout << endl << "Added boundaries------" << endl << endl;
- triangulate (flow);
- if ( debug ) cout << endl << "Tesselating------" << endl << endl;
- flow.T[flow.currentTes].compute();
- flow.defineFictiousCells();
-
- //FIXME: this is already done in addBoundary(?)
- boundaryConditions ( flow );
- if ( debug ) cout << endl << "boundaryConditions------" << endl << endl;
- flow.initializePressure ( pZero );
- if ( debug ) cout << endl << "initializePressure------" << endl << endl;
- // Define the ghost cells and add indexes to the cells inside the period (the ones that will contain the pressure unknowns)
- //This must be done after boundary conditions and initialize pressure, else the indexes are not good (not accounting imposedP): FIXME
- unsigned int index=0;
- int baseIndex=-1;
- FlowSolver::Tesselation& Tes = flow.T[flow.currentTes];
- Tes.cellHandles.resize(Tes.Triangulation().number_of_finite_cells());
- const FiniteCellsIterator cellend=Tes.Triangulation().finite_cells_end();
- for ( FiniteCellsIterator cell=Tes.Triangulation().finite_cells_begin(); cell!=cellend; cell++ ){
- locateCell ( cell,index,baseIndex,flow );
- if (flow.errorCode>0) return;
- //Fill this vector than can be later used to speedup loops
- if (!cell->info().isGhost) Tes.cellHandles[cell->info().baseIndex]=cell;
- }
- Tes.cellHandles.resize(baseIndex+1);
-
- if ( debug ) cout << endl << "locateCell------" << endl << endl;
- flow.computePermeability ( );
- porosity = flow.vPoralPorosity/flow.vTotalPorosity;
- flow.tolerance=tolerance;flow.relax=relax;
-
- flow.displayStatistics ();
- //FIXME: check interpolate() for the periodic case, at least use the mean pressure from previous step.
- if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], Tes );
-// if ( !first && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
-
- if ( waveAction ) flow.applySinusoidalPressure ( Tes.Triangulation(), sineMagnitude, sineAverage, 30 );
-
- if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
- if ( debug ) cout << endl << "end buildTri------" << endl << endl;
-}
-
-void PeriodicFlowEngine::preparePShifts()
-{
- CellInfo::gradP = makeCgVect ( gradP );
- CellInfo::hSize[0] = makeCgVect ( scene->cell->hSize.col ( 0 ) );
- CellInfo::hSize[1] = makeCgVect ( scene->cell->hSize.col ( 1 ) );
- CellInfo::hSize[2] = makeCgVect ( scene->cell->hSize.col ( 2 ) );
- CellInfo::deltaP=CGT::CVector (
- CellInfo::hSize[0]*CellInfo::gradP,
- CellInfo::hSize[1]*CellInfo::gradP,
- CellInfo::hSize[2]*CellInfo::gradP );
-}
-
-
-YADE_PLUGIN((PeriodicFlowEngine));
-
-
-#endif //FLOW_ENGINE
-
-#endif /* YADE_CGAL */
\ No newline at end of file
=== removed file 'pkg/dem/PeriodicFlowEngine.hpp'
--- pkg/dem/PeriodicFlowEngine.hpp 2014-04-03 13:07:09 +0000
+++ pkg/dem/PeriodicFlowEngine.hpp 1970-01-01 00:00:00 +0000
@@ -1,109 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* Copyright (C) 2012 by Donia Marzougui <donia.marzougui@xxxxxxxxxxxxxxx>*
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-/// The periodic variant of FlowEngine is defined here. It should become a template class for more flexibility.
-/// It is a bit more complicated as for FlowEngine, though, because we need template inheriting from template, which breaks YADE_CLASS_XXX logic_error
-/// See below the commented exemple, for a possible solution
-
-#include <yade/pkg/dem/FlowEngine.hpp>
-
-class PeriodicCellInfo : public FlowCellInfo
-{
- public:
- static CVector gradP;
- //for real cell, baseIndex is the rank of the cell in cellHandles. For ghost cells, it is the baseIndex of the corresponding real cell.
- //Unlike ordinary index, baseIndex is also indexing cells with imposed pressures
- int baseIndex;
- int period[3];
- static CVector hSize[3];
- static CVector deltaP;
- int ghost;
- Real* _pression;
- PeriodicCellInfo (void){
- _pression=&pression;
- period[0]=period[1]=period[2]=0;
- baseIndex=-1;
- volumeSign=0;}
- ~PeriodicCellInfo (void) {}
-
- inline const double shiftedP (void) const {return isGhost? (*_pression)+pShift() :(*_pression) ;}
- inline const double pShift (void) const {return deltaP[0]*period[0] + deltaP[1]*period[1] +deltaP[2]*period[2];}
-// inline const double p (void) {return shiftedP();}
- inline void setP (const Real& p) {pression=p;}
- bool isReal (void) {return !(isFictious || isGhost);}
-};
-
-
-class PeriodicVertexInfo : public FlowVertexInfo {
- public:
- PeriodicVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
- PeriodicVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
- PeriodicVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
- int period[3];
- //FIXME: the name is misleading, even non-ghost can be out of the period and therefore they need to be shifted as well
- inline const CVector ghostShift (void) {
- return period[0]*PeriodicCellInfo::hSize[0]+period[1]*PeriodicCellInfo::hSize[1]+period[2]*PeriodicCellInfo::hSize[2];}
- PeriodicVertexInfo (void) {isFictious=false; s=0; i=0; period[0]=period[1]=period[2]=0; isGhost=false;}
- bool isReal (void) {return !(isFictious || isGhost);}
-};
-
-typedef CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> PeriFlowTriangulationTypes;
-typedef CGT::PeriodicTesselation<CGT::_Tesselation<PeriFlowTriangulationTypes> > PeriFlowTesselation;
-#ifdef LINSOLV
-#define _PeriFlowSolver CGT::PeriodicFlowLinSolv<PeriFlowTesselation>
-#else
-#define _PeriFlowSolver CGT::PeriodicFlow<PeriFlowTesselation>
-#endif
-
-typedef TemplateFlowEngine< PeriodicCellInfo,
- PeriodicVertexInfo,
- CGT::PeriodicTesselation<CGT::_Tesselation<CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> > >,
- _PeriFlowSolver>
- FlowEngine_PeriodicInfo;
-
-REGISTER_SERIALIZABLE(FlowEngine_PeriodicInfo);
-YADE_PLUGIN((FlowEngine_PeriodicInfo));
-
-
-class PeriodicFlowEngine : public FlowEngine_PeriodicInfo
-{
- public :
- void triangulate (FlowSolver& flow);
- void buildTriangulation (Real pzero, FlowSolver& flow);
- void initializeVolumes (FlowSolver& flow);
- void updateVolumes (FlowSolver& flow);
- Real volumeCell (CellHandle cell);
-
- Real volumeCellSingleFictious (CellHandle cell);
- inline void locateCell(CellHandle baseCell, unsigned int& index, int& baseIndex, FlowSolver& flow, unsigned int count=0);
- Vector3r meanVelocity();
-
- virtual ~PeriodicFlowEngine();
-
- virtual void action();
- //Cache precomputed values for pressure shifts, based on current hSize and pGrad
- void preparePShifts();
-
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(PeriodicFlowEngine,FlowEngine_PeriodicInfo,"A variant of :yref:`FlowEngine` implementing periodic boundary conditions. The API is very similar.",
- ((Real,duplicateThreshold, 0.06,,"distance from cell borders that will triger periodic duplication in the triangulation |yupdate|"))
- ((Vector3r, gradP, Vector3r::Zero(),,"Macroscopic pressure gradient"))
- ,,
- wallIds=vector<int>(6,-1);
- solver = shared_ptr<FlowSolver> (new FlowSolver);
- epsVolMax=epsVolCumulative=retriangulationLastIter=0;
- ReTrg=1;
- first=true;
- ,
- //nothing special to define, we re-use FlowEngine methods
- //.def("meanVelocity",&PeriodicFlowEngine::meanVelocity,"measure the mean velocity in the period")
- )
- DECLARE_LOGGER;
-
-
-};
-REGISTER_SERIALIZABLE(PeriodicFlowEngine);
\ No newline at end of file
=== added directory 'pkg/pfv'
=== added file 'pkg/pfv/DFNFlow.cpp'
--- pkg/pfv/DFNFlow.cpp 1970-01-01 00:00:00 +0000
+++ pkg/pfv/DFNFlow.cpp 2014-04-03 13:37:55 +0000
@@ -0,0 +1,85 @@
+
+/*************************************************************************
+* Copyright (C) 2014 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+#ifdef FLOW_ENGINE
+
+//keep this #ifdef for commited versions unless you really have stable version that should be compiled by default
+//it will save compilation time for everyone else
+//when you want it compiled, you can pass -DDFNFLOW to cmake, or just uncomment the following line
+// #define DFNFLOW
+#ifdef DFNFLOW
+#define TEMPLATE_FLOW_NAME DFNFlowEngineT
+#include <yade/pkg/dem/FlowEngine.hpp>
+
+class DFNCellInfo : public FlowCellInfo
+{
+ public:
+ Real anotherVariable;
+ void anotherFunction() {};
+};
+
+
+class DFNVertexInfo : public FlowVertexInfo {
+ public:
+ //same here if needed
+};
+
+typedef TemplateFlowEngine<DFNCellInfo,DFNVertexInfo> DFNFlowEngineT;
+REGISTER_SERIALIZABLE(DFNFlowEngineT);
+YADE_PLUGIN((DFNFlowEngineT));
+
+class DFNFlowEngine : public DFNFlowEngineT
+{
+ public :
+ void trickPermeability();
+ void trickPermeability (RTriangulation::Facet_circulator& facet,Real somethingBig);
+ void trickPermeability (RTriangulation::Finite_edges_iterator& edge,Real somethingBig);
+
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DFNFlowEngine,DFNFlowEngineT,"documentation here",
+ ((Real, myNewAttribute, 0,,"useless example"))
+ ,/*DFNFlowEngineT()*/,
+ ,
+ )
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(DFNFlowEngine);
+YADE_PLUGIN((DFNFlowEngine));
+
+void DFNFlowEngine::trickPermeability (RTriangulation::Facet_circulator& facet, Real somethingBig)
+{
+ const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ const CellHandle& cell1 = facet->first;
+ const CellHandle& cell2 = facet->first->neighbor(facet->second);
+ if ( Tri.is_infinite(cell1) || Tri.is_infinite(cell2)) cerr<<"Infinite cell found in trickPermeability, should be handled somehow, maybe"<<endl;
+ cell1->info().kNorm()[facet->second] = somethingBig;
+ cell2->info().kNorm()[Tri.mirror_index(cell1, facet->second)] = somethingBig;
+}
+
+void DFNFlowEngine::trickPermeability(RTriangulation::Finite_edges_iterator& edge, Real somethingBig)
+{
+ const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ RTriangulation::Facet_circulator facet1 = Tri.incident_facets(*edge);
+ RTriangulation::Facet_circulator facet0=facet1++;
+ trickPermeability(facet0,somethingBig);
+ while ( facet1!=facet0 ) {trickPermeability(facet1,somethingBig); facet1++;}
+}
+
+
+void DFNFlowEngine::trickPermeability()
+{
+ Real somethingBig=10;//is that big??
+ const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ //We want to change permeability perpendicular to the 10th edge, let's say.
+ //in the end this function should have a loop on all edges I guess
+ FiniteEdgesIterator edge = Tri.finite_edges_begin();
+ for(int k=10; k>0;--k) edge++;
+ trickPermeability(edge,somethingBig);
+}
+
+
+#endif //DFNFLOW
+#endif //FLOWENGINE
\ No newline at end of file
=== added file 'pkg/pfv/DummyFlowEngine.cpp'
--- pkg/pfv/DummyFlowEngine.cpp 1970-01-01 00:00:00 +0000
+++ pkg/pfv/DummyFlowEngine.cpp 2014-04-03 13:37:55 +0000
@@ -0,0 +1,60 @@
+
+/*************************************************************************
+* Copyright (C) 2014 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+// This is an example of how to derive a new FlowEngine with additional data and possibly completely new behaviour.
+// Every functions of the base engine can be overloaded, and new functions can be added
+
+//keep this #ifdef as long as you don't really want to realize a final version publicly, it will save compilation time for everyone else
+//when you want it compiled, you can pass -DDFNFLOW to cmake, or just uncomment the following line
+// #define DUMMYFLOW
+#ifdef DUMMYFLOW
+#define TEMPLATE_FLOW_NAME DummyFlowEngineT
+#include <yade/pkg/dem/FlowEngine.hpp>
+
+/// We can add data to the Info types by inheritance
+class DummyCellInfo : public FlowCellInfo
+{
+ public:
+ Real anotherVariable;
+ void anotherFunction() {};
+};
+
+class DummyVertexInfo : public FlowVertexInfo {
+ public:
+ //same here if needed
+};
+
+typedef TemplateFlowEngine<DummyCellInfo,DummyVertexInfo> TEMPLATE_FLOW_NAME;
+REGISTER_SERIALIZABLE(TEMPLATE_FLOW_NAME);
+YADE_PLUGIN((TEMPLATE_FLOW_NAME));
+
+class DummyFlowEngine : public TEMPLATE_FLOW_NAME
+{
+ public :
+ //We can overload every functions of the base engine to make it behave differently
+ //if we overload action() like this, this engine is doing nothing in a standard timestep, it can still have useful functions
+ virtual void action() {};
+
+ //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); {cerr<<"yes, I'm a new function"<<end;}
+
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DummyFlowEngine,TEMPLATE_FLOW_NAME,"documentation here",
+ ((Real, myNewAttribute, 0,,"useless example"))
+ ,/*DummyFlowEngineT()*/,
+ ,
+ .def("fancyFunction",&DummyFlowEngine::fancyFunction,(python::arg("what")=0),"test function")
+ )
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(DummyFlowEngine);
+YADE_PLUGIN((DummyFlowEngine));
+
+void DummyFlowEngine::fancyFunction(Real what) {cerr<<"yes, I'm a new function"<<end;}
+#undef TEMPLATE_FLOW_NAME DummyFlowEngineT //To be sure it will not conflict, maybe not needed
+#endif //DummyFLOW
\ No newline at end of file
=== added file 'pkg/pfv/FlowEngine.cpp'
--- pkg/pfv/FlowEngine.cpp 1970-01-01 00:00:00 +0000
+++ pkg/pfv/FlowEngine.cpp 2014-04-03 13:37:55 +0000
@@ -0,0 +1,40 @@
+/*************************************************************************
+* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
+* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+#ifdef YADE_CGAL
+#ifdef FLOW_ENGINE
+
+#define TEMPLATE_FLOW_NAME FlowEngineT
+#include "FlowEngine.hpp"
+
+// To register properly, we need to first instantiate an intermediate class, then inherit from it with correct class names in YADE_CLASS macro
+// The intermediate one would be seen with the name "TemplateFlowEngine" by python, thus it would not work when more than one class are derived, they would all
+// be named "TemplateFlowEngine" ...
+typedef TemplateFlowEngine<FlowCellInfo,FlowVertexInfo> FlowEngineT;
+REGISTER_SERIALIZABLE(FlowEngineT);
+
+class FlowEngine : public FlowEngineT
+{
+ public :
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(FlowEngine,FlowEngineT,"An engine to solve flow problem in saturated granular media. Model description can be found in [Chareyre2012a]_ and [Catalano2014a]_. 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.",
+ ,,
+ ,
+ //nothing special to define here, we simply re-use FlowEngine methods
+ //.def("meanVelocity",&PeriodicFlowEngine::meanVelocity,"measure the mean velocity in the period")
+ )
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(FlowEngine);
+
+YADE_PLUGIN((FlowEngineT));
+CREATE_LOGGER(FlowEngine );
+YADE_PLUGIN((FlowEngine));
+
+#endif //FLOW_ENGINE
+
+#endif /* YADE_CGAL */
+
=== added file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp 1970-01-01 00:00:00 +0000
+++ pkg/pfv/FlowEngine.hpp 2014-04-03 13:37:55 +0000
@@ -0,0 +1,433 @@
+/*************************************************************************
+* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
+* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+/*!
+
+FlowEngine is an interface between Yade and the fluid solvers defined in lib/triangulation and using the PFV scheme.
+There are also a few non-trivial functions defined here, such has building triangulation and computating elements volumes.
+The code strongly relies on CGAL library for building triangulations on the top of sphere packings.
+CGAL's trangulations introduce the data associated to each element of the triangulation through template parameters (Cell_info and Vertex_info).
+FlowEngine is thus defined via the template class TemplateFlowEngine(Cellinfo,VertexInfo), for easier addition of variables in various couplings.
+PeriodicFlowEngine is a variant for periodic boundary conditions.
+
+Which solver will be actually used internally to obtain pore pressure will depend partly on compile time flags (libcholmod available and #define LINSOLV),
+and on runtime settings (useSolver=0: Gauss-Seidel (iterative), useSolver=1: CHOLMOD if available (direct sparse cholesky)).
+
+The files defining lower level classes are in yade/lib. The code uses CGAL::Triangulation3 for managing the mesh and storing data. Eigen3::Sparse, suitesparse::cholmod, and metis are used for solving the linear systems with a direct method (Cholesky). An iterative method (Gauss-Seidel) is implemented directly in Yade and can be used as a fallback (see FlowEngine::useSolver).
+
+Most classes in lib/triangulation are templates, and are therefore completely defined in header files.
+A pseudo hpp/cpp split is reproduced for clarity, with hpp/ipp extensions (but again, in the end they are all include files).
+The files are
+- RegularTriangulation.h : declaration of info types embeded in the mesh (template parameters of CGAL::Triangulation3), and instanciation of RTriangulation classes
+- Tesselation.h/ipp : encapsulate RegularTriangulation and adds functions to manipulate the dual (Voronoi) graph of the triangulation
+- Network.hpp/ipp : specialized for PFV model (the two former are used independently by TesselationWrapper), a set of functions to determine volumes and surfaces of intersections between spheres and various subdomains. Contains two triangulations for smooth transitions while remeshing - e.g. interpolating values in the new mesh using the previous one.
+- FlowBoundingSphere.hpp/ipp and PeriodicFlow.hpp/ipp + LinSolv variants: implement the solver in itself (mesh, boundary conditions, solving, defining fluid-particles interactions)
+- FlowEngine.hpp/ipp/cpp (this file)
+
+Variants for periodic boundary conditions are also present.
+
+*/
+
+#pragma once
+#include<yade/core/PartialEngine.hpp>
+#include<yade/pkg/dem/TriaxialCompressionEngine.hpp>
+#include<yade/pkg/dem/TesselationWrapper.hpp>
+#include<yade/lib/triangulation/FlowBoundingSphere.hpp>
+#include<yade/lib/triangulation/PeriodicFlow.hpp>
+
+/// Frequently used:
+typedef CGT::CVector CVector;
+typedef CGT::Point Point;
+
+/// Converters for Eigen and CGAL vectors
+inline CVector makeCgVect ( const Vector3r& yv ) {return CVector ( yv[0],yv[1],yv[2] );}
+inline Point makeCgPoint ( const Vector3r& yv ) {return Point ( yv[0],yv[1],yv[2] );}
+inline Vector3r makeVector3r ( const Point& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
+inline Vector3r makeVector3r ( const CVector& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
+
+/// C++ templates and YADE_CLASS_... macro family are not very compatible, this #define is a hack to make it work
+/// TEMPLATE_FLOW_NAME will be the name of a serializable TemplateFlowEngine<...> instance, which can in turn be
+/// inherited from. The instance itself will be useless for actual simulations.
+#ifndef TEMPLATE_FLOW_NAME
+#error You must define TEMPLATE_FLOW_NAME in your source file before including FlowEngine.hpp
+#endif
+
+template<class _CellInfo, class _VertexInfo, class _Tesselation=CGT::_Tesselation<CGT::TriangulationTypes<_VertexInfo,_CellInfo> >, class solverT=CGT::FlowBoundingSphere<_Tesselation> >
+class TemplateFlowEngine : public PartialEngine
+{
+ public :
+ #ifdef LINSOLV
+ typedef CGT::FlowBoundingSphereLinSolv<typename solverT::Tesselation,solverT> FlowSolver;
+ #else
+ typedef solverT FlowSolver;
+ #endif
+ typedef FlowSolver Solver;//FIXME: useless alias, search/replace then remove
+ typedef typename FlowSolver::Tesselation Tesselation;
+ typedef typename FlowSolver::RTriangulation RTriangulation;
+ typedef typename FlowSolver::FiniteVerticesIterator FiniteVerticesIterator;
+ typedef typename FlowSolver::FiniteCellsIterator FiniteCellsIterator;
+ typedef typename FlowSolver::CellHandle CellHandle;
+ typedef typename FlowSolver::FiniteEdgesIterator FiniteEdgesIterator;
+ typedef typename FlowSolver::VertexHandle VertexHandle;
+ typedef typename RTriangulation::Triangulation_data_structure::Cell::Info CellInfo;
+ typedef typename RTriangulation::Triangulation_data_structure::Vertex::Info VertexInfo;
+
+ protected:
+ shared_ptr<FlowSolver> solver;
+ shared_ptr<FlowSolver> backgroundSolver;
+ volatile bool backgroundCompleted;
+ Cell cachedCell;
+ struct posData {Body::id_t id; Vector3r pos; Real radius; bool isSphere; bool exists; posData(){exists=0;}};
+ vector<posData> positionBufferCurrent;//reflect last known positions before we start computations
+ vector<posData> positionBufferParallel;//keep the positions from a given step for multithread factorization
+ //copy positions in a buffer for faster and/or parallel access
+ void setPositionsBuffer(bool current);
+ virtual void trickPermeability() {};
+
+ public :
+ int retriangulationLastIter;
+ enum {wall_xmin, wall_xmax, wall_ymin, wall_ymax, wall_zmin, wall_zmax};
+ Vector3r normal [6];
+ bool currentTes;
+ bool metisForced;
+ int idOffset;
+ double epsVolCumulative;
+ int ReTrg;
+ int ellapsedIter;
+ void initSolver (FlowSolver& flow);
+ #ifdef LINSOLV
+ void setForceMetis (bool force);
+ bool getForceMetis ();
+ #endif
+ void triangulate (Solver& flow);
+ void addBoundary (Solver& flow);
+ void buildTriangulation (double pZero, Solver& flow);
+ void buildTriangulation (Solver& flow);
+ void updateVolumes (Solver& flow);
+ void initializeVolumes (Solver& flow);
+ void boundaryConditions(Solver& flow);
+ void updateBCs ( Solver& flow ) {
+ if (flow.T[flow.currentTes].maxId>0) boundaryConditions(flow);//avoids crash at iteration 0, when the packing is not bounded yet
+ else LOG_ERROR("updateBCs not applied");
+ flow.pressureChanged=true;}
+
+ void imposeFlux(Vector3r pos, Real flux);
+ unsigned int imposePressure(Vector3r pos, Real p);
+ void setImposedPressure(unsigned int cond, Real p);
+ void clearImposedPressure();
+ void clearImposedFlux();
+ void computeViscousForces(Solver& flow);
+ Real getCellFlux(unsigned int cond);
+ Real getBoundaryFlux(unsigned int boundary) {return solver->boundaryFlux(boundary);}
+ Vector3r fluidForce(unsigned int idSph) {
+ const CGT::CVector& f=solver->T[solver->currentTes].vertex(idSph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
+
+ Vector3r shearLubForce(unsigned int id_sph) {
+ return (solver->shearLubricationForces.size()>id_sph)?solver->shearLubricationForces[id_sph]:Vector3r::Zero();}
+ Vector3r shearLubTorque(unsigned int id_sph) {
+ return (solver->shearLubricationTorques.size()>id_sph)?solver->shearLubricationTorques[id_sph]:Vector3r::Zero();}
+ Vector3r pumpLubTorque(unsigned int id_sph) {
+ return (solver->pumpLubricationTorques.size()>id_sph)?solver->pumpLubricationTorques[id_sph]:Vector3r::Zero();}
+ Vector3r twistLubTorque(unsigned int id_sph) {
+ return (solver->twistLubricationTorques.size()>id_sph)?solver->twistLubricationTorques[id_sph]:Vector3r::Zero();}
+ Vector3r normalLubForce(unsigned int id_sph) {
+ return (solver->normalLubricationForce.size()>id_sph)?solver->normalLubricationForce[id_sph]:Vector3r::Zero();}
+ Matrix3r bodyShearLubStress(unsigned int id_sph) {
+ return (solver->shearLubricationBodyStress.size()>id_sph)?solver->shearLubricationBodyStress[id_sph]:Matrix3r::Zero();}
+ Matrix3r bodyNormalLubStress(unsigned int id_sph) {
+ return (solver->normalLubricationBodyStress.size()>id_sph)?solver->normalLubricationBodyStress[id_sph]:Matrix3r::Zero();}
+ Vector3r shearVelocity(unsigned int interaction) {
+ return (solver->deltaShearVel[interaction]);}
+ Vector3r normalVelocity(unsigned int interaction) {
+ return (solver->deltaNormVel[interaction]);}
+ Matrix3r normalStressInteraction(unsigned int interaction) {
+ return (solver->normalStressInteraction[interaction]);}
+ Matrix3r shearStressInteraction(unsigned int interaction) {
+ return (solver->shearStressInteraction[interaction]);}
+ Vector3r normalVect(unsigned int interaction) {
+ return (solver->normalV[interaction]);}
+ Real surfaceDistanceParticle(unsigned int interaction) {
+ return (solver->surfaceDistance[interaction]);}
+ Real edgeSize() {
+ return (solver->edgeIds.size());}
+ Real OSI() {
+ return (solver->onlySpheresInteractions.size());}
+ int onlySpheresInteractions(unsigned int interaction) {
+ return (solver->onlySpheresInteractions[interaction]);}
+ python::list getConstrictions(bool all) {
+ vector<Real> csd=solver->getConstrictions(); python::list pycsd;
+ for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k]!=0) || csd[k]>0) pycsd.append(csd[k]); return pycsd;}
+ python::list getConstrictionsFull(bool all) {
+ vector<Constriction> csd=solver->getConstrictionsFull(); python::list pycsd;
+ for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k].second[0]!=0) || csd[k].second[0]>0) {
+ python::list cons;
+ cons.append(csd[k].first.first);
+ cons.append(csd[k].first.second);
+ cons.append(csd[k].second[0]);
+ cons.append(csd[k].second[1]);
+ cons.append(csd[k].second[2]);
+ cons.append(csd[k].second[3]);
+ pycsd.append(cons);}
+ return pycsd;}
+
+ template<class Cellhandle>
+ Real volumeCellSingleFictious (Cellhandle cell);
+ template<class Cellhandle>
+ Real volumeCellDoubleFictious (Cellhandle cell);
+ template<class Cellhandle>
+ Real volumeCellTripleFictious (Cellhandle cell);
+ template<class Cellhandle>
+ Real volumeCell (Cellhandle cell);
+ void Oedometer_Boundary_Conditions();
+ void averageRealCellVelocity();
+ void saveVtk(const char* folder) {solver->saveVtk(folder);}
+ vector<Real> avFlVelOnSph(unsigned int idSph) {return solver->averageFluidVelocityOnSphere(idSph);}
+
+// 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]);}
+ int getCell(double posX, double posY, double posZ){return solver->getCell(posX, posY, posZ);}
+ unsigned int nCells(){return solver->T[solver->currentTes].cellHandles.size();}
+ python::list getVertices(unsigned int id){
+ python::list ids;
+ if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return ids;}
+ for (unsigned int i=0;i<4;i++) ids.append(solver->T[solver->currentTes].cellHandles[id]->vertex(i)->info().id());
+ return ids;
+ }
+ double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
+ double averagePressure(){return solver->averagePressure();}
+
+ void emulateAction(){
+ scene = Omega::instance().getScene().get();
+ action();}
+ #ifdef LINSOLV
+ void exportMatrix(string filename) {if (useSolver==3) solver->exportMatrix(filename.c_str());
+ else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
+ void exportTriplets(string filename) {if (useSolver==3) solver->exportTriplets(filename.c_str());
+ else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
+ void cholmodStats() {
+ cerr << cholmod_print_common((char*)string("PFV Cholmod factorization").c_str(),&(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);}
+ #endif
+
+ virtual ~TemplateFlowEngine();
+ virtual void action();
+ virtual void backgroundAction();
+
+ //commodities
+ void compTessVolumes() {
+ solver->T[solver->currentTes].compute();
+ solver->T[solver->currentTes].computeVolumes();
+ }
+ Real getVolume (Body::id_t id) {
+ if (solver->T[solver->currentTes].Max_id() <= 0) {emulateAction(); LOG_WARN("Not triangulated yet, emulating action");}
+ if (solver->T[solver->currentTes].Volume(id) == -1) {compTessVolumes(); LOG_WARN("Computing all volumes now, as you did not request it explicitely.");}
+ return (solver->T[solver->currentTes].Max_id() >= id) ? solver->T[solver->currentTes].Volume(id) : -1;}
+
+ YADE_CLASS_PYCLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TemplateFlowEngine,TEMPLATE_FLOW_NAME,PartialEngine,"An engine to solve flow problem in saturated granular media. Model description can be found in [Chareyre2012a]_ and [Catalano2014a]_. 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.",
+ ((bool,isActivated,true,,"Activates Flow Engine"))
+ ((bool,first,true,,"Controls the initialization/update phases"))
+ ((double, fluidBulkModulus, 0.,,"Bulk modulus of fluid (inverse of compressibility) K=-dP*V/dV [Pa]. Flow is compressible if fluidBulkModulus > 0, else incompressible."))
+ ((Real, dt, 0,,"timestep [s]"))
+ ((bool,permeabilityMap,false,,"Enable/disable stocking of average permeability scalar in cell infos."))
+ ((bool, slipBoundary, true,, "Controls friction condition on lateral walls"))
+ ((bool,waveAction, false,, "Allow sinusoidal pressure condition to simulate ocean waves"))
+ ((double, sineMagnitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied (p )"))
+ ((double, sineAverage, 0,,"Pressure value (average) when sinusoidal pressure is applied"))
+ ((bool, debug, false,,"Activate debug messages"))
+ ((double, wallThickness,0.001,,"Walls thickness"))
+ ((double,pZero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model."))
+ ((double,tolerance,1e-06,,"Gauss-Seidel tolerance"))
+ ((double,relax,1.9,,"Gauss-Seidel relaxation"))
+ ((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, epsVolMax, 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."))
+ ((Real,minKdivKmean,0.0001,,"define the min K value (see :yref:`FlowEngine::clampKValues`)"))
+ ((Real,maxKdivKmean,100,,"define the max K value (see :yref:`FlowEngine::clampKValues`)"))
+ ((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, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso, 3-CHOLMOD"))
+ ((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`"))
+ ((int, ignoredBody,-1,,"Id of a sphere to exclude from the triangulation.)"))
+ ((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.) "))
+ ((bool, pumpTorque, false,,"Compute pump torque applied on particles "))
+ ((bool, twistTorque, false,,"Compute twist torque applied on particles "))
+ ((double, eps, 0.00001,,"roughness defined as a fraction of particles size, giving the minimum distance between particles in the lubrication model."))
+ ((bool, pressureForce, true,,"compute the pressure field and associated fluid forces. WARNING: turning off means fluid flow is not computed at all."))
+ ((bool, normalLubrication, false,,"compute normal lubrication force as developped by Brule"))
+ ((bool, viscousNormalBodyStress, false,,"compute normal viscous stress applied on each body"))
+ ((bool, viscousShearBodyStress, false,,"compute shear viscous stress applied on each body"))
+ ((bool, multithread, false,,"Build triangulation and factorize in the background (multi-thread mode)"))
+ #ifdef EIGENSPARSE_LIB
+ ((int, numSolveThreads, 1,,"number of openblas threads in the solve phase."))
+ ((int, numFactorizeThreads, 1,,"number of openblas threads in the factorization phase"))
+ #endif
+ ,
+ /*deprec*/
+ ((meanK_opt,clampKValues,"the name changed"))
+ ,,
+ timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
+ 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;
+ solver = shared_ptr<FlowSolver> (new FlowSolver);
+ first=true;
+ epsVolMax=epsVolCumulative=retriangulationLastIter=0;
+ ReTrg=1;
+ backgroundCompleted=true;
+ ellapsedIter=0;
+ metisForced=false;
+ ,
+ .def("imposeFlux",&TemplateFlowEngine::imposeFlux,(python::arg("pos"),python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
+ .def("imposePressure",&TemplateFlowEngine::imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
+ .def("setImposedPressure",&TemplateFlowEngine::setImposedPressure,(python::arg("cond"),python::arg("p")),"Set pressure value at the point indexed 'cond'.")
+ .def("clearImposedPressure",&TemplateFlowEngine::clearImposedPressure,"Clear the list of points with pressure imposed.")
+ .def("clearImposedFlux",&TemplateFlowEngine::clearImposedFlux,"Clear the list of points with flux imposed.")
+ .def("getCellFlux",&TemplateFlowEngine::getCellFlux,(python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
+ .def("getBoundaryFlux",&TemplateFlowEngine::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("getConstrictions",&TemplateFlowEngine::getConstrictions,(python::arg("all")=true),"Get the list of constriction radii (inscribed circle) for all finite facets (if all==True) or all facets not incident to a virtual bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres.")
+ .def("getConstrictionsFull",&TemplateFlowEngine::getConstrictionsFull,(python::arg("all")=true),"Get the list of constrictions (inscribed circle) for all finite facets (if all==True), or all facets not incident to a fictious bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres. The constrictions are returned in the format {{cell1,cell2}{rad,nx,ny,nz}}")
+ .def("edgeSize",&TemplateFlowEngine::edgeSize,"Return the number of interactions.")
+ .def("OSI",&TemplateFlowEngine::OSI,"Return the number of interactions only between spheres.")
+
+ .def("saveVtk",&TemplateFlowEngine::saveVtk,(python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
+ .def("avFlVelOnSph",&TemplateFlowEngine::avFlVelOnSph,(python::arg("idSph")),"compute a sphere-centered average fluid velocity")
+ .def("fluidForce",&TemplateFlowEngine::fluidForce,(python::arg("idSph")),"Return the fluid force on sphere idSph.")
+ .def("shearLubForce",&TemplateFlowEngine::shearLubForce,(python::arg("idSph")),"Return the shear lubrication force on sphere idSph.")
+ .def("shearLubTorque",&TemplateFlowEngine::shearLubTorque,(python::arg("idSph")),"Return the shear lubrication torque on sphere idSph.")
+ .def("normalLubForce",&TemplateFlowEngine::normalLubForce,(python::arg("idSph")),"Return the normal lubrication force on sphere idSph.")
+ .def("bodyShearLubStress",&TemplateFlowEngine::bodyShearLubStress,(python::arg("idSph")),"Return the shear lubrication stress on sphere idSph.")
+ .def("bodyNormalLubStress",&TemplateFlowEngine::bodyNormalLubStress,(python::arg("idSph")),"Return the normal lubrication stress on sphere idSph.")
+ .def("shearVelocity",&TemplateFlowEngine::shearVelocity,(python::arg("idSph")),"Return the shear velocity of the interaction.")
+ .def("normalVelocity",&TemplateFlowEngine::normalVelocity,(python::arg("idSph")),"Return the normal velocity of the interaction.")
+ .def("normalVect",&TemplateFlowEngine::normalVect,(python::arg("idSph")),"Return the normal vector between particles.")
+ .def("surfaceDistanceParticle",&TemplateFlowEngine::surfaceDistanceParticle,(python::arg("interaction")),"Return the distance between particles.")
+ .def("onlySpheresInteractions",&TemplateFlowEngine::onlySpheresInteractions,(python::arg("interaction")),"Return the id of the interaction only between spheres.")
+ .def("pressureProfile",&TemplateFlowEngine::pressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
+ .def("getPorePressure",&TemplateFlowEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
+ .def("averageSlicePressure",&TemplateFlowEngine::averageSlicePressure,(python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
+ .def("averagePressure",&TemplateFlowEngine::averagePressure,"Measure averaged pore pressure in the entire volume")
+ .def("updateBCs",&TemplateFlowEngine::updateBCs,"tells the engine to update it's boundary conditions before running (especially useful when changing boundary pressure - should not be needed for point-wise imposed pressure)")
+ .def("emulateAction",&TemplateFlowEngine::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
+ .def("getCell",&TemplateFlowEngine::getCell,(python::arg("pos")),"get id of the cell containing (X,Y,Z).")
+ .def("nCells",&TemplateFlowEngine::nCells,"get the total number of finite cells in the triangulation.")
+ .def("getVertices",&TemplateFlowEngine::getVertices,(python::arg("id")),"get the vertices of a cell")
+ #ifdef LINSOLV
+ .def("exportMatrix",&TemplateFlowEngine::exportMatrix,(python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
+ .def("exportTriplets",&TemplateFlowEngine::exportTriplets,(python::arg("filename")="triplets"),"Export system matrix to a file with only non-zero entries.")
+ .def("cholmodStats",&TemplateFlowEngine::cholmodStats,"get statistics of cholmod solver activity")
+ .def("metisUsed",&TemplateFlowEngine::metisUsed,"check wether metis lib is effectively used")
+ .add_property("forceMetis",&TemplateFlowEngine::getForceMetis,&TemplateFlowEngine::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")
+ #endif
+ .def("compTessVolumes",&TemplateFlowEngine::compTessVolumes,"Like TesselationWrapper::computeVolumes()")
+ .def("volume",&TemplateFlowEngine::getVolume,(python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
+ )
+};
+// Definition of functions in a separate file for clarity
+#include<yade/pkg/dem/FlowEngine.ipp>
+
+class FlowCellInfo : public CGT::SimpleCellInfo {
+ public:
+ //For vector storage of all cells
+ unsigned int index;
+ int volumeSign;
+ bool Pcondition;
+ Real invVoidV;
+ Real t;
+ int fict;
+ Real volumeVariation;
+ double pression;
+ //average relative (fluid - facet translation) velocity defined for a single cell as 1/Volume * SUM_ON_FACETS(x_average_facet*average_facet_flow_rate)
+ CVector averageCellVelocity;
+ // Surface vectors of facets, pointing from outside toward inside the cell
+ std::vector<CVector> facetSurfaces;
+ //Ratio between fluid surface and facet surface
+ std::vector<Real> facetFluidSurfacesRatio;
+ // Reflects the geometrical property of the cell, so that the force by cell fluid on grain "i" is pressure*unitForceVectors[i]
+ std::vector<CVector> unitForceVectors;
+ // Store the area of triangle-sphere intersections for each facet (used in forces definition)
+ std::vector<CVector> facetSphereCrossSections;
+ std::vector<CVector> cellForce;
+ std::vector<double> rayHydr;
+ std::vector<double> modulePermeability;
+ // Partial surfaces of spheres in the double-tetrahedron linking two voronoi centers. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_surfaces in the facet.
+ double solidSurfaces [4][4];
+
+ FlowCellInfo (void)
+ {
+ modulePermeability.resize(4, 0);
+ cellForce.resize(4);
+ facetSurfaces.resize(4);
+ facetFluidSurfacesRatio.resize(4);
+ facetSphereCrossSections.resize(4);
+ unitForceVectors.resize(4);
+ for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidSurfaces[k][l]=0;
+ rayHydr.resize(4, 0);
+ invSumK=index=volumeSign=s=volumeVariation=pression=invVoidV=fict=0;
+ isFictious=false; Pcondition = false; isGhost = false;
+ isvisited = false;
+ isGhost=false;
+ }
+ bool isGhost;
+ double invSumK;
+ bool isvisited;
+
+ inline Real& volume (void) {return t;}
+ inline const Real& invVoidVolume (void) const {return invVoidV;}
+ inline Real& invVoidVolume (void) {return invVoidV;}
+ inline Real& dv (void) {return volumeVariation;}
+ inline int& fictious (void) {return fict;}
+ inline double& p (void) {return pression;}
+ inline const double shiftedP (void) const {return pression;} //For compatibility with the periodic case
+ inline const std::vector<double>& kNorm (void) const {return modulePermeability;}
+ inline std::vector<double>& kNorm (void) {return modulePermeability;}
+ inline std::vector< CVector >& facetSurf (void) {return facetSurfaces;}
+ inline std::vector<CVector>& force (void) {return cellForce;}
+ inline std::vector<double>& Rh (void) {return rayHydr;}
+ inline CVector& averageVelocity (void) {return averageCellVelocity;}
+ //used for transfering values between two triangulations, overload with more variables in derived classes (see e.g. SoluteFlow)
+ inline void getInfo(const FlowCellInfo& otherCellInfo) {p()=otherCellInfo.shiftedP();}
+};
+
+class FlowVertexInfo : public CGT::SimpleVertexInfo {
+ CVector grainVelocity;
+ Real volumeIncidentCells;
+public:
+ CVector forces;
+ bool isGhost;
+ FlowVertexInfo (void) {isGhost=false;}
+ inline CVector& force (void) {return forces;}
+ inline CVector& vel (void) {return grainVelocity;}
+ inline Real& volCells (void) {return volumeIncidentCells;}
+ inline const CVector ghostShift (void) {return CGAL::NULL_VECTOR;}
+};
+
+
+
=== added file 'pkg/pfv/FlowEngine.ipp'
--- pkg/pfv/FlowEngine.ipp 1970-01-01 00:00:00 +0000
+++ pkg/pfv/FlowEngine.ipp 2014-04-03 13:37:55 +0000
@@ -0,0 +1,730 @@
+/*************************************************************************
+* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
+* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+#ifdef YADE_CGAL
+
+#ifdef FLOW_ENGINE
+#include<yade/core/Scene.hpp>
+#include<yade/lib/base/Math.hpp>
+#include<yade/pkg/dem/TesselationWrapper.hpp>
+#include<yade/pkg/common/Sphere.hpp>
+#include<yade/pkg/common/Wall.hpp>
+#include<yade/pkg/common/Box.hpp>
+#include <sys/stat.h>
+#include <sys/types.h>
+#include <boost/thread.hpp>
+#include <boost/date_time.hpp>
+#include <boost/date_time/posix_time/posix_time.hpp>
+
+#ifdef LINSOLV
+#include <cholmod.h>
+#endif
+
+#include "FlowEngine.hpp"
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::~TemplateFlowEngine() {}
+
+// YADE_PLUGIN((TFlowEng));
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+unsigned int TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposePressure(Vector3r pos, Real p)
+{
+// if (!flow) LOG_ERROR("no flow defined yet, run at least one iter");
+ solver->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),p) );
+ //force immediate update of boundary conditions
+ updateTriangulation=true;
+ return solver->imposedP.size()-1;
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::action()
+{
+ if ( !isActivated ) return;
+ timingDeltas->start();
+ setPositionsBuffer(true);
+ timingDeltas->checkpoint ( "Position buffer" );
+ if (first) {
+ if (multithread) setPositionsBuffer(false);
+ buildTriangulation(pZero,*solver);
+ initializeVolumes(*solver);
+ backgroundSolver=solver;
+ backgroundCompleted=true;
+ }
+ solver->ompThreads = ompThreads>0? ompThreads : omp_get_max_threads();
+
+ timingDeltas->checkpoint ( "Triangulating" );
+ updateVolumes ( *solver );
+ timingDeltas->checkpoint ( "Update_Volumes" );
+
+ epsVolCumulative += epsVolMax;
+ retriangulationLastIter++;
+ if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
+ (defTolerance>0 && epsVolCumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval;
+
+ ///compute flow and and forces here
+ if (pressureForce){
+ solver->gaussSeidel(scene->dt);
+ timingDeltas->checkpoint ( "Gauss-Seidel (includes matrix construct and factorization in single-thread mode)" );
+ solver->computeFacetForcesWithCache();}
+ timingDeltas->checkpoint ( "compute_Forces" );
+ ///Application of vicscous forces
+ scene->forces.sync();
+ timingDeltas->checkpoint ( "forces.sync()" );
+ computeViscousForces ( *solver );
+ timingDeltas->checkpoint ( "viscous forces" );
+ Vector3r force;
+ Vector3r torque;
+ FiniteVerticesIterator verticesEnd = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
+ for ( FiniteVerticesIterator vIt = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); vIt != verticesEnd; vIt++ ) {
+ force = pressureForce ? Vector3r ( vIt->info().forces[0],vIt->info().forces[1],vIt->info().forces[2] ): Vector3r(0,0,0);
+ torque = Vector3r(0,0,0);
+ if (shearLubrication || viscousShear){
+ force = force + solver->shearLubricationForces[vIt->info().id()];
+ torque = torque + solver->shearLubricationTorques[vIt->info().id()];
+ if (pumpTorque)
+ torque = torque + solver->pumpLubricationTorques[vIt->info().id()];
+ }
+ if (twistTorque)
+ torque = torque + solver->twistLubricationTorques[vIt->info().id()];
+ if (normalLubrication)
+ force = force + solver-> normalLubricationForce[vIt->info().id()];
+ scene->forces.addForce ( vIt->info().id(), force);
+ scene->forces.addTorque ( vIt->info().id(), torque);
+ }
+ ///End compute flow and forces
+ timingDeltas->checkpoint ( "Applying Forces" );
+ int sleeping = 0;
+ if (multithread && !first) {
+ while (updateTriangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/
+ sleeping++;
+ boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
+ if (debug && sleeping) cerr<<"sleeping..."<<sleeping<<endl;
+ if (updateTriangulation || (ellapsedIter>(0.5*meshUpdateInterval) && backgroundCompleted)) {
+ if (debug) cerr<<"switch flow solver"<<endl;
+ if (useSolver==0) LOG_ERROR("background calculations not available for Gauss-Seidel");
+ if (fluidBulkModulus>0) solver->interpolate (solver->T[solver->currentTes], backgroundSolver->T[backgroundSolver->currentTes]);
+ solver=backgroundSolver;
+ backgroundSolver = shared_ptr<FlowSolver> (new FlowSolver);
+ if (metisForced) {backgroundSolver->eSolver.cholmod().nmethods=1; backgroundSolver->eSolver.cholmod().method[0].ordering=CHOLMOD_METIS;}
+ //Copy imposed pressures/flow from the old solver
+ backgroundSolver->imposedP = vector<pair<CGT::Point,Real> >(solver->imposedP);
+ backgroundSolver->imposedF = vector<pair<CGT::Point,Real> >(solver->imposedF);
+ if (debug) cerr<<"switched"<<endl;
+ setPositionsBuffer(false);//set "parallel" buffer for background calculation
+ backgroundCompleted=false;
+ retriangulationLastIter=ellapsedIter;
+ updateTriangulation=false;
+ epsVolCumulative=0;
+ ellapsedIter=0;
+ boost::thread workerThread(&TemplateFlowEngine::backgroundAction,this);
+ workerThread.detach();
+ if (debug) cerr<<"backgrounded"<<endl;
+ initializeVolumes(*solver);
+ computeViscousForces(*solver);
+ if (debug) cerr<<"volumes initialized"<<endl;
+ }
+ else {
+ if (debug && !backgroundCompleted) cerr<<"still computing solver in the background, ellapsedIter="<<ellapsedIter<<endl;
+ ellapsedIter++;
+ }
+ } else {
+ if (updateTriangulation && !first) {
+ buildTriangulation (pZero, *solver);
+ initializeVolumes(*solver);
+ computeViscousForces(*solver);
+ updateTriangulation = false;
+ epsVolCumulative=0;
+ retriangulationLastIter=0;
+ ReTrg++;}
+ }
+ first=false;
+ timingDeltas->checkpoint ( "triangulate + init volumes" );
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::backgroundAction()
+{
+ if (useSolver<1) {LOG_ERROR("background calculations not available for Gauss-Seidel"); return;}
+ buildTriangulation ( pZero,*backgroundSolver );
+ //FIXME: GS is computing too much, we need only matrix factorization in fact
+ backgroundSolver->gaussSeidel(scene->dt);
+ //FIXME(2): and here we need only cached variables, not forces
+ backgroundSolver->computeFacetForcesWithCache(/*onlyCache?*/ true);
+// boost::this_thread::sleep(boost::posix_time::seconds(5));
+ backgroundCompleted = true;
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::boundaryConditions ( Solver& flow )
+{
+ for (int k=0;k<6;k++) {
+ flow.boundary (wallIds[k]).flowCondition=!bndCondIsPressure[k];
+ flow.boundary (wallIds[k]).value=bndCondValue[k];
+ flow.boundary (wallIds[k]).velocity = boundaryVelocity[k];//FIXME: needs correct implementation, maybe update the cached pos/vel?
+ }
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setImposedPressure ( unsigned int cond, Real p)
+{
+ if ( cond>=solver->imposedP.size() ) LOG_ERROR ( "Setting p with cond higher than imposedP size." );
+ solver->imposedP[cond].second=p;
+ //force immediate update of boundary conditions
+ solver->pressureChanged=true;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposeFlux ( Vector3r pos, Real flux){
+ solver->imposedF.push_back ( pair<CGT::Point,Real> ( CGT::Point ( pos[0],pos[1],pos[2] ), flux ) );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedPressure () { solver->imposedP.clear(); solver->IPCells.clear();}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedFlux () { solver->imposedF.clear(); solver->IFCells.clear();}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::getCellFlux ( unsigned int cond)
+{
+ if ( cond>=solver->imposedP.size() ) {LOG_ERROR ( "Getting flux with cond higher than imposedP size." ); return 0;}
+ double flux=0;
+ typename Solver::CellHandle& cell= solver->IPCells[cond];
+ for ( int ngb=0;ngb<4;ngb++ ) {
+ /*if (!cell->neighbor(ngb)->info().Pcondition)*/
+ flux+= cell->info().kNorm() [ngb]* ( cell->info().p()-cell->neighbor ( ngb )->info().p() );
+ }
+ return flux+cell->info().dv();
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::initSolver ( FlowSolver& flow )
+{
+ flow.Vtotalissimo=0; flow.VSolidTot=0; flow.vPoral=0; flow.sSolidTot=0;
+ flow.slipBoundary=slipBoundary;
+ flow.kFactor = permeabilityFactor;
+ flow.debugOut = debug;
+ flow.useSolver = useSolver;
+ #ifdef EIGENSPARSE_LIB
+ flow.numSolveThreads = numSolveThreads;
+ flow.numFactorizeThreads = numFactorizeThreads;
+ #endif
+ flow.meanKStat = meanKStat;
+ flow.viscosity = viscosity;
+ flow.tolerance=tolerance;
+ flow.relax=relax;
+ flow.clampKValues = clampKValues;
+ flow.maxKdivKmean = maxKdivKmean;
+ flow.minKdivKmean = minKdivKmean;
+ flow.meanKStat = meanKStat;
+ flow.permeabilityMap = permeabilityMap;
+ flow.fluidBulkModulus = fluidBulkModulus;
+ flow.T[flow.currentTes].Clear();
+ flow.T[flow.currentTes].maxId=-1;
+ flow.xMin = 1000.0, flow.xMax = -10000.0, flow.yMin = 1000.0, flow.yMax = -10000.0, flow.zMin = 1000.0, flow.zMax = -10000.0;
+}
+
+#ifdef LINSOLV
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setForceMetis ( bool force )
+{
+ if (force) {
+ metisForced=true;
+ solver->eSolver.cholmod().nmethods=1;
+ solver->eSolver.cholmod().method[0].ordering=CHOLMOD_METIS;
+ } else {cholmod_defaults(&(solver->eSolver.cholmod())); metisForced=false;}
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+bool TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::getForceMetis () {return (solver->eSolver.cholmod().nmethods==1);}
+#endif
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( Solver& flow )
+{
+ buildTriangulation ( 0.f,flow );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( double pZero, Solver& flow )
+{
+ flow.resetNetwork();
+ if (first) flow.currentTes=0;
+ else {
+ flow.currentTes=!flow.currentTes;
+ if (debug) cout << "--------RETRIANGULATION-----------" << endl;
+ }
+
+ initSolver(flow);
+
+ addBoundary ( flow );
+ triangulate ( flow );
+ if ( debug ) cout << endl << "Tesselating------" << endl << endl;
+ flow.T[flow.currentTes].compute();
+
+ flow.defineFictiousCells();
+ // For faster loops on cells define this vector
+ flow.T[flow.currentTes].cellHandles.clear();
+ flow.T[flow.currentTes].cellHandles.reserve(flow.T[flow.currentTes].Triangulation().number_of_finite_cells());
+ FiniteCellsIterator cell_end = flow.T[flow.currentTes].Triangulation().finite_cells_end();
+ int k=0;
+ for ( FiniteCellsIterator cell = flow.T[flow.currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
+ flow.T[flow.currentTes].cellHandles.push_back(cell);
+ cell->info().id=k++;}//define unique numbering now, corresponds to position in cellHandles
+ flow.displayStatistics ();
+ flow.computePermeability();
+ //This virtual function does nothing yet, derived class may overload it to make permeability different (see DFN engine)
+ trickPermeability();
+ porosity = flow.vPoralPorosity/flow.vTotalPorosity;
+
+ boundaryConditions ( flow );
+ flow.initializePressure ( pZero );
+
+ if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
+ if ( waveAction ) flow.applySinusoidalPressure ( flow.T[flow.currentTes].Triangulation(), sineMagnitude, sineAverage, 30 );
+ if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setPositionsBuffer(bool current)
+{
+ vector<posData>& buffer = current? positionBufferCurrent : positionBufferParallel;
+ buffer.clear();
+ buffer.resize(scene->bodies->size());
+ shared_ptr<Sphere> sph ( new Sphere );
+ const int Sph_Index = sph->getClassIndexStatic();
+ FOREACH ( const shared_ptr<Body>& b, *scene->bodies ) {
+ if (!b || ignoredBody==b->getId()) continue;
+ posData& dat = buffer[b->getId()];
+ dat.id=b->getId();
+ dat.pos=b->state->pos;
+ dat.isSphere= (b->shape->getClassIndex() == Sph_Index);
+ if (dat.isSphere) dat.radius = YADE_CAST<Sphere*>(b->shape.get())->radius;
+ dat.exists=true;
+ }
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::addBoundary ( Solver& flow )
+{
+ vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
+ solver->xMin = Mathr::MAX_REAL, solver->xMax = -Mathr::MAX_REAL, solver->yMin = Mathr::MAX_REAL, solver->yMax = -Mathr::MAX_REAL, solver->zMin = Mathr::MAX_REAL, solver->zMax = -Mathr::MAX_REAL;
+ FOREACH ( const posData& b, buffer ) {
+ if ( !b.exists ) continue;
+ if ( b.isSphere ) {
+ const Real& rad = b.radius;
+ const Real& x = b.pos[0];
+ const Real& y = b.pos[1];
+ const Real& z = b.pos[2];
+ flow.xMin = min ( flow.xMin, x-rad );
+ flow.xMax = max ( flow.xMax, x+rad );
+ flow.yMin = min ( flow.yMin, y-rad );
+ flow.yMax = max ( flow.yMax, y+rad );
+ flow.zMin = min ( flow.zMin, z-rad );
+ flow.zMax = max ( flow.zMax, z+rad );
+ }
+ }
+ //FIXME idOffset must be set correctly, not the case here (always 0), then we need walls first or it will fail
+ idOffset = flow.T[flow.currentTes].maxId+1;
+ flow.idOffset = idOffset;
+ flow.sectionArea = ( flow.xMax - flow.xMin ) * ( flow.zMax-flow.zMin );
+ flow.vTotal = ( flow.xMax-flow.xMin ) * ( flow.yMax-flow.yMin ) * ( flow.zMax-flow.zMin );
+ flow.yMinId=wallIds[ymin];
+ flow.yMaxId=wallIds[ymax];
+ flow.xMaxId=wallIds[xmax];
+ flow.xMinId=wallIds[xmin];
+ flow.zMinId=wallIds[zmin];
+ flow.zMaxId=wallIds[zmax];
+
+ //FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
+ flow.boundsIds[0]= &flow.xMinId;
+ flow.boundsIds[1]= &flow.xMaxId;
+ flow.boundsIds[2]= &flow.yMinId;
+ flow.boundsIds[3]= &flow.yMaxId;
+ flow.boundsIds[4]= &flow.zMinId;
+ flow.boundsIds[5]= &flow.zMaxId;
+
+ for (int k=0;k<6;k++) flow.boundary ( *flow.boundsIds[k] ).useMaxMin = boundaryUseMaxMin[k];
+
+ flow.cornerMin = CGT::Point ( flow.xMin, flow.yMin, flow.zMin );
+ flow.cornerMax = CGT::Point ( flow.xMax, flow.yMax, flow.zMax );
+
+ //assign BCs types and values
+ boundaryConditions ( flow );
+
+ double center[3];
+ for ( int i=0; i<6; i++ ) {
+ if ( *flow.boundsIds[i]<0 ) continue;
+ CGT::CVector Normal ( normal[i].x(), normal[i].y(), normal[i].z() );
+ if ( flow.boundary ( *flow.boundsIds[i] ).useMaxMin ) flow.addBoundingPlane(Normal, *flow.boundsIds[i] );
+ else {
+ for ( int h=0;h<3;h++ ) center[h] = buffer[*flow.boundsIds[i]].pos[h];
+// cerr << "id="<<*flow.boundsIds[i] <<" center="<<center[0]<<","<<center[1]<<","<<center[2]<<endl;
+ flow.addBoundingPlane ( center, wallThickness, Normal,*flow.boundsIds[i] );
+ }
+ }
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::triangulate ( Solver& flow )
+{
+///Using Tesselation wrapper (faster)
+// TesselationWrapper TW;
+// if (TW.Tes) delete TW.Tes;
+// TW.Tes = &(flow.T[flow.currentTes]);//point to the current Tes we have in Flowengine
+// TW.insertSceneSpheres();//TW is now really inserting in TemplateFlowEngine, using the faster insert(begin,end)
+// TW.Tes = NULL;//otherwise, Tes would be deleted by ~TesselationWrapper() at the end of the function.
+///Using one-by-one insertion
+ vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
+ FOREACH ( const posData& b, buffer ) {
+ if ( !b.exists ) continue;
+ if ( b.isSphere ) {
+ if (b.id==ignoredBody) continue;
+ flow.T[flow.currentTes].insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
+ }
+ flow.T[flow.currentTes].redirected=true;//By inserting one-by-one, we already redirected
+ flow.shearLubricationForces.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.shearLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.pumpLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.twistLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.shearLubricationBodyStress.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.normalLubricationForce.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.normalLubricationBodyStress.resize ( flow.T[flow.currentTes].maxId+1 );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::initializeVolumes ( Solver& flow )
+{
+ typedef typename Solver::FiniteVerticesIterator FiniteVerticesIterator;
+
+ FiniteVerticesIterator vertices_end = flow.T[flow.currentTes].Triangulation().finite_vertices_end();
+ CGT::CVector Zero(0,0,0);
+ for (FiniteVerticesIterator V_it = flow.T[flow.currentTes].Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero;
+
+ FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles)
+ {
+ switch ( cell->info().fictious() )
+ {
+ case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
+ case ( 1 ) : cell->info().volume() = volumeCellSingleFictious ( cell ); break;
+ case ( 2 ) : cell->info().volume() = volumeCellDoubleFictious ( cell ); break;
+ case ( 3 ) : cell->info().volume() = volumeCellTripleFictious ( cell ); break;
+ default: break;
+ }
+ if (flow.fluidBulkModulus>0) { cell->info().invVoidVolume() = 1 / ( abs(cell->info().volume()) - flow.volumeSolidPore(cell) ); }
+ }
+ if (debug) cout << "Volumes initialised." << endl;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::averageRealCellVelocity()
+{
+ solver->averageRelativeCellVelocity();
+ Vector3r Vel ( 0,0,0 );
+ //AVERAGE CELL VELOCITY
+ FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
+ for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
+ for ( int g=0;g<4;g++ ) {
+ if ( !cell->vertex ( g )->info().isFictious ) {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ for ( int i=0;i<3;i++ ) Vel[i] = Vel[i] + sph->state->vel[i]/4;
+ }
+ }
+ RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ CGT::Point pos_av_facet;
+ double volume_facet_translation = 0;
+ CGT::CVector Vel_av ( Vel[0], Vel[1], Vel[2] );
+ for ( int i=0; i<4; i++ ) {
+ volume_facet_translation = 0;
+ if ( !Tri.is_infinite ( cell->neighbor ( i ) ) ) {
+ CGT::CVector Surfk = cell->info()-cell->neighbor ( i )->info();
+ Real area = sqrt ( Surfk.squared_length() );
+ Surfk = Surfk/area;
+ CGT::CVector branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
+ pos_av_facet = ( CGT::Point ) cell->info() + ( branch*Surfk ) *Surfk;
+ volume_facet_translation += Vel_av*cell->info().facetSurfaces[i];
+ cell->info().averageVelocity() = cell->info().averageVelocity() - volume_facet_translation/cell->info().volume() * ( pos_av_facet-CGAL::ORIGIN );
+ }
+ }
+ }
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::updateVolumes ( Solver& flow )
+{
+ if ( debug ) cout << "Updating volumes.............." << endl;
+ Real invDeltaT = 1/scene->dt;
+ epsVolMax=0;
+ Real totVol=0; Real totDVol=0;
+ #ifdef YADE_OPENMP
+ const long size=flow.T[flow.currentTes].cellHandles.size();
+ #pragma omp parallel for num_threads(ompThreads>0 ? ompThreads : 1)
+ for(long i=0; i<size; i++){
+ CellHandle& cell = flow.T[flow.currentTes].cellHandles[i];
+ #else
+ FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
+ #endif
+ double newVol, dVol;
+ switch ( cell->info().fictious() ) {
+ case ( 3 ) : newVol = volumeCellTripleFictious ( cell ); break;
+ case ( 2 ) : newVol = volumeCellDoubleFictious ( cell ); break;
+ case ( 1 ) : newVol = volumeCellSingleFictious ( cell ); break;
+ case ( 0 ) : newVol = volumeCell (cell ); break;
+ default: newVol = 0; break;}
+ dVol=cell->info().volumeSign* ( newVol - cell->info().volume() );
+ cell->info().dv() = dVol*invDeltaT;
+ cell->info().volume() = newVol;
+ if (defTolerance>0) { //if the criterion is not used, then we skip these updates a save a LOT of time when Nthreads > 1
+ #pragma omp atomic
+ totVol+=newVol;
+ #pragma omp atomic
+ totDVol+=abs(dVol);}
+ }
+ if (defTolerance>0) epsVolMax = totDVol/totVol;
+ for (unsigned int n=0; n<flow.imposedF.size();n++) {
+ flow.IFCells[n]->info().dv()+=flow.imposedF[n].second;
+ flow.IFCells[n]->info().Pcondition=false;}
+ if ( debug ) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<endl;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellSingleFictious ( Cellhandle cell )
+{
+ Vector3r V[3];
+ int b=0;
+ int w=0;
+ cell->info().volumeSign=1;
+ Real Wall_coordinate=0;
+
+ for ( int y=0;y<4;y++ ) {
+ if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
+ V[w]=positionBufferCurrent[cell->vertex ( y )->info().id()].pos;
+ w++;
+ } else {
+ b = cell->vertex ( y )->info().id();
+ const shared_ptr<Body>& wll = Body::byId ( b , scene );
+ if ( !solver->boundary ( b ).useMaxMin ) Wall_coordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wallThickness/2.;
+ else Wall_coordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate];
+ }
+ }
+ Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
+ return abs ( Volume );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellDoubleFictious ( Cellhandle cell )
+{
+ Vector3r A=Vector3r::Zero(), AS=Vector3r::Zero(),B=Vector3r::Zero(), BS=Vector3r::Zero();
+
+ cell->info().volumeSign=1;
+ int b[2];
+ int coord[2];
+ Real Wall_coordinate[2];
+ int j=0;
+ bool first_sph=true;
+
+ for ( int g=0;g<4;g++ ) {
+ if ( cell->vertex ( g )->info().isFictious ) {
+ b[j] = cell->vertex ( g )->info().id();
+ coord[j]=solver->boundary ( b[j] ).coordinate;
+ if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = positionBufferCurrent[b[j]].pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2.;
+ else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
+ j++;
+ } else if ( first_sph ) {
+ A=AS=/*AT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;
+ first_sph=false;
+ } else {
+ B=BS=/*BT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;;
+ }
+ }
+ AS[coord[0]]=BS[coord[0]] = Wall_coordinate[0];
+
+ //first pyramid with triangular base (A,B,BS)
+ Real Vol1=0.5* ( ( A-BS ).cross ( B-BS ) ) [coord[1]]* ( 0.333333333* ( 2*B[coord[1]]+A[coord[1]] )-Wall_coordinate[1] );
+ //second pyramid with triangular base (A,AS,BS)
+ Real Vol2=0.5* ( ( AS-BS ).cross ( A-BS ) ) [coord[1]]* ( 0.333333333* ( B[coord[1]]+2*A[coord[1]] )-Wall_coordinate[1] );
+ return abs ( Vol1+Vol2 );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellTripleFictious ( Cellhandle cell )
+{
+ Vector3r A;
+
+ int b[3];
+ int coord[3];
+ Real Wall_coordinate[3];
+ int j=0;
+ cell->info().volumeSign=1;
+
+ for ( int g=0;g<4;g++ ) {
+ if ( cell->vertex ( g )->info().isFictious ) {
+ b[j] = cell->vertex ( g )->info().id();
+ coord[j]=solver->boundary ( b[j] ).coordinate;
+ const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
+ if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = wll->state->pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2.;
+ else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
+ j++;
+ } else {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ A= ( sph->state->pos );
+ }
+ }
+ Real Volume = ( A[coord[0]]-Wall_coordinate[0] ) * ( A[coord[1]]-Wall_coordinate[1] ) * ( A[coord[2]]-Wall_coordinate[2] );
+ return abs ( Volume );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCell ( Cellhandle cell )
+{
+ static const Real inv6 = 1/6.;
+ const Vector3r& p0 = positionBufferCurrent[cell->vertex ( 0 )->info().id()].pos;
+ const Vector3r& p1 = positionBufferCurrent[cell->vertex ( 1 )->info().id()].pos;
+ const Vector3r& p2 = positionBufferCurrent[cell->vertex ( 2 )->info().id()].pos;
+ const Vector3r& p3 = positionBufferCurrent[cell->vertex ( 3 )->info().id()].pos;
+ Real volume = inv6 * ((p0-p1).cross(p0-p2)).dot(p0-p3);
+ if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
+ return volume;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::computeViscousForces ( Solver& flow )
+{
+ if (normalLubrication || shearLubrication || viscousShear){
+ if ( debug ) cout << "Application of viscous forces" << endl;
+ if ( debug ) cout << "Number of edges = " << flow.edgeIds.size() << endl;
+ for ( unsigned int k=0; k<flow.shearLubricationForces.size(); k++ ) flow.shearLubricationForces[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.shearLubricationTorques.size(); k++ ) flow.shearLubricationTorques[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.pumpLubricationTorques.size(); k++ ) flow.pumpLubricationTorques[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.twistLubricationTorques.size(); k++ ) flow.twistLubricationTorques[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.shearLubricationBodyStress.size(); k++) flow.shearLubricationBodyStress[k]=Matrix3r::Zero();
+ for ( unsigned int k=0; k<flow.normalLubricationForce.size(); k++ ) flow.normalLubricationForce[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.normalLubricationBodyStress.size(); k++) flow.normalLubricationBodyStress[k]=Matrix3r::Zero();
+
+ typedef typename Solver::Tesselation Tesselation;
+ const Tesselation& Tes = flow.T[flow.currentTes];
+ flow.deltaShearVel.clear(); flow.normalV.clear(); flow.deltaNormVel.clear(); flow.surfaceDistance.clear(); flow.onlySpheresInteractions.clear(); flow.normalStressInteraction.clear(); flow.shearStressInteraction.clear();
+
+
+ for ( int i=0; i< ( int ) flow.edgeIds.size(); i++ ) {
+ const int& id1 = flow.edgeIds[i].first;
+ const int& id2 = flow.edgeIds[i].second;
+
+ int hasFictious= Tes.vertex ( id1 )->info().isFictious + Tes.vertex ( id2 )->info().isFictious;
+ if (hasFictious>0 or id1==id2) continue;
+ const shared_ptr<Body>& sph1 = Body::byId ( id1, scene );
+ const shared_ptr<Body>& sph2 = Body::byId ( id2, scene );
+ Sphere* s1=YADE_CAST<Sphere*> ( sph1->shape.get() );
+ Sphere* s2=YADE_CAST<Sphere*> ( sph2->shape.get() );
+ const Real& r1 = s1->radius;
+ const Real& r2 = s2->radius;
+ Vector3r deltaV; Real deltaNormV; Vector3r deltaShearV;
+ Vector3r O1O2Vector; Real O1O2; Vector3r normal; Real surfaceDist; Vector3r O1CVector; Vector3r O2CVector;Real meanRad ;Real Rh; Vector3r deltaAngVel; Vector3r deltaShearAngVel;
+ Vector3r shearLubF; Vector3r normaLubF; Vector3r pumpT; Vector3r deltaAngNormVel; Vector3r twistT; Vector3r angVel1; Vector3r angVel2;
+ //FIXME: if periodic and velGrad!=0, then deltaV should account for velGrad, not the case currently
+ if ( !hasFictious ){
+ O1O2Vector = sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()) - sph1->state->pos - makeVector3r(Tes.vertex(id1)->info().ghostShift());
+ O1O2 = O1O2Vector.norm();
+ normal= (O1O2Vector/O1O2);
+ surfaceDist = O1O2 - r2 - r1;
+ O1CVector = (O1O2/2. + (pow(r1,2) - pow(r2,2)) / (2.*O1O2))*normal;
+ O2CVector = -(O1O2Vector - O1CVector);
+ meanRad = (r2 + r1)/2.;
+ Rh = (r1 < r2)? surfaceDist + 0.45 * r1 : surfaceDist + 0.45 * r2;
+ deltaV = (sph2->state->vel + sph2->state->angVel.cross(-r2 * normal)) - (sph1->state->vel+ sph1->state->angVel.cross(r1 * normal));
+ angVel1 = sph1->state->angVel;
+ angVel2 = sph2->state->angVel;
+ deltaAngVel = sph2->state->angVel - sph1->state->angVel;
+
+ } else {
+ if ( hasFictious==1 ) {//for the fictious sphere, use velocity of the boundary, not of the body
+ bool v1fictious = Tes.vertex ( id1 )->info().isFictious;
+ int bnd = v1fictious? id1 : id2;
+ int coord = flow.boundary(bnd).coordinate;
+ O1O2 = v1fictious ? abs((sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]) : abs((sph1->state->pos + makeVector3r(Tes.vertex(id1)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]);
+ if(v1fictious)
+ normal = makeVector3r(flow.boundary(id1).normal);
+ else
+ normal = -makeVector3r(flow.boundary(id2).normal);
+ O1O2Vector = O1O2 * normal;
+ meanRad = v1fictious ? r2:r1;
+ surfaceDist = O1O2 - meanRad;
+ if (v1fictious){
+ O1CVector = Vector3r::Zero();
+ O2CVector = - O1O2Vector;}
+ else{
+ O1CVector = O1O2Vector;
+ O2CVector = Vector3r::Zero();}
+
+ Rh = surfaceDist + 0.45 * meanRad;
+ Vector3r v1 = ( Tes.vertex ( id1 )->info().isFictious ) ? flow.boundary ( id1 ).velocity:sph1->state->vel + sph1->state->angVel.cross(r1 * normal);
+ Vector3r v2 = ( Tes.vertex ( id2 )->info().isFictious ) ? flow.boundary ( id2 ).velocity:sph2->state->vel + sph2->state->angVel.cross(-r2 * (normal));
+ deltaV = v2-v1;
+ angVel1 = ( Tes.vertex ( id1 )->info().isFictious ) ? Vector3r::Zero() : sph1->state->angVel;
+ angVel2 = ( Tes.vertex ( id2 )->info().isFictious ) ? Vector3r::Zero() : sph2->state->angVel;
+ deltaAngVel = angVel2 - angVel1;
+ }
+ }
+ deltaShearV = deltaV - ( normal.dot ( deltaV ) ) *normal;
+ deltaShearAngVel = deltaAngVel - ( normal.dot ( deltaAngVel ) ) *normal;
+ flow.deltaShearVel.push_back(deltaShearV);
+ flow.normalV.push_back(normal);
+ flow.surfaceDistance.push_back(max(surfaceDist, 0.) + eps*meanRad);
+
+ /// Compute the shear Lubrication force and torque on each particle
+
+ if (shearLubrication)
+ shearLubF = flow.computeShearLubricationForce(deltaShearV,surfaceDist,i,eps,O1O2,meanRad);
+ else if (viscousShear)
+ shearLubF = flow.computeViscousShearForce ( deltaShearV, i , Rh);
+
+ if (viscousShear || shearLubrication){
+
+ flow.shearLubricationForces[id1]+=shearLubF;
+ flow.shearLubricationForces[id2]+=(-shearLubF);
+ flow.shearLubricationTorques[id1]+=O1CVector.cross(shearLubF);
+ flow.shearLubricationTorques[id2]+=O2CVector.cross(-shearLubF);
+
+ /// Compute the pump Lubrication torque on each particle
+
+ if (pumpTorque){
+ pumpT = flow.computePumpTorque(deltaShearAngVel, surfaceDist, i, eps, meanRad );
+ flow.pumpLubricationTorques[id1]+=(-pumpT);
+ flow.pumpLubricationTorques[id2]+=pumpT;}
+
+ /// Compute the twist Lubrication torque on each particle
+
+ if (twistTorque){
+ deltaAngNormVel = (normal.dot(deltaAngVel))*normal ;
+ twistT = flow.computeTwistTorque(deltaAngNormVel, surfaceDist, i, eps, meanRad );
+ flow.twistLubricationTorques[id1]+=(-twistT);
+ flow.twistLubricationTorques[id2]+=twistT;
+ }
+ }
+ /// Compute the viscous shear stress on each particle
+
+ if (viscousShearBodyStress){
+ flow.shearLubricationBodyStress[id1] += shearLubF * O1CVector.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
+ flow.shearLubricationBodyStress[id2] += (-shearLubF) * O2CVector.transpose()/ (4.0/3.0 *3.14* pow(r2,3));
+ flow.shearStressInteraction.push_back(shearLubF * O1O2Vector.transpose()/(4.0/3.0 *3.14* pow(r1,3)));
+ }
+
+ /// Compute the normal lubrication force applied on each particle
+
+ if (normalLubrication){
+ deltaNormV = normal.dot(deltaV);
+ flow.deltaNormVel.push_back(deltaNormV * normal);
+ normaLubF = flow.computeNormalLubricationForce (deltaNormV, surfaceDist, i,eps,stiffness,scene->dt,meanRad)*normal;
+ flow.normalLubricationForce[id1]+=normaLubF;
+ flow.normalLubricationForce[id2]+=(-normaLubF);
+
+ /// Compute the normal lubrication stress on each particle
+
+ if (viscousNormalBodyStress){
+ flow.normalLubricationBodyStress[id1] += normaLubF * O1CVector.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
+ flow.normalLubricationBodyStress[id2] += (-normaLubF) *O2CVector.transpose() / (4.0/3.0 *3.14* pow(r2,3));
+ flow.normalStressInteraction.push_back(normaLubF * O1O2Vector.transpose()/(4.0/3.0 *3.14* pow(r1,3)));
+ }
+ }
+
+ if (!hasFictious)
+ flow.onlySpheresInteractions.push_back(i);
+
+ }
+ }
+}
+
+#endif //FLOW_ENGINE
+
+#endif /* YADE_CGAL */
+
=== added file 'pkg/pfv/PeriodicFlowEngine.cpp'
--- pkg/pfv/PeriodicFlowEngine.cpp 1970-01-01 00:00:00 +0000
+++ pkg/pfv/PeriodicFlowEngine.cpp 2014-04-03 13:37:55 +0000
@@ -0,0 +1,452 @@
+/*************************************************************************
+* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
+* Copyright (C) 2012 by Donia Marzougui <donia.marzougui@xxxxxxxxxxxxxxx>*
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#ifdef YADE_CGAL
+#ifdef FLOW_ENGINE
+
+#define TEMPLATE_FLOW_NAME FlowEngine_PeriodicInfo
+#include "PeriodicFlowEngine.hpp"
+#undef TEMPLATE_FLOW_NAME
+
+CVector PeriodicCellInfo::hSize[]={CVector(),CVector(),CVector()};
+CVector PeriodicCellInfo::deltaP=CVector();
+CVector PeriodicCellInfo::gradP=CVector();
+
+CREATE_LOGGER ( PeriodicFlowEngine );
+
+PeriodicFlowEngine::~PeriodicFlowEngine(){}
+
+void PeriodicFlowEngine:: action()
+{
+ if ( !isActivated ) return;
+ timingDeltas->start();
+ preparePShifts();
+ setPositionsBuffer(true);
+ if (first) {
+ if (multithread) setPositionsBuffer(false);
+ cachedCell= Cell(*(scene->cell));
+ buildTriangulation(pZero,*solver);
+ if (solver->errorCode>0) {LOG_INFO("triangulation error, pausing"); Omega::instance().pause(); return;}
+ initializeVolumes(*solver); backgroundSolver=solver; backgroundCompleted=true;}
+// if ( first ) {buildTriangulation ( pZero ); updateTriangulation = false; initializeVolumes();}
+
+ timingDeltas->checkpoint("Triangulating");
+ updateVolumes (*solver);
+ epsVolCumulative += epsVolMax;
+ retriangulationLastIter++;
+ if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
+ (defTolerance>0 && epsVolCumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval;
+
+ timingDeltas->checkpoint("Update_Volumes");
+
+ ///compute flow and and forces here
+ if (pressureForce){
+ solver->gaussSeidel(scene->dt);
+ timingDeltas->checkpoint("Gauss-Seidel");
+ solver->computeFacetForcesWithCache();}
+ timingDeltas->checkpoint("compute_Pressure_Forces");
+
+ ///compute vicscous forces
+ scene->forces.sync();
+ computeViscousForces(*solver);
+ timingDeltas->checkpoint("compute_Viscous_Forces");
+ Vector3r force;
+ Vector3r torque;
+ const Tesselation& Tes = solver->T[solver->currentTes];
+ for (int id=0; id<=Tes.maxId; id++){
+ assert (Tes.vertexHandles[id] != NULL);
+ const Tesselation::VertexInfo& v_info = Tes.vertexHandles[id]->info();
+ force =(pressureForce) ? Vector3r ( ( v_info.forces ) [0],v_info.forces[1],v_info.forces[2] ) : Vector3r(0,0,0);
+ torque = Vector3r(0,0,0);
+ if (shearLubrication || viscousShear){
+ force = force +solver->shearLubricationForces[v_info.id()];
+ torque = torque +solver->shearLubricationTorques[v_info.id()];
+ if (pumpTorque)
+ torque = torque +solver->pumpLubricationTorques[v_info.id()];
+ if (twistTorque)
+ torque = torque +solver->twistLubricationTorques[v_info.id()];
+ }
+
+ if (normalLubrication)
+ force = force + solver->normalLubricationForce[v_info.id()];
+ scene->forces.addForce ( v_info.id(), force);
+ scene->forces.addTorque ( v_info.id(), torque);
+ }
+ ///End Compute flow and forces
+ timingDeltas->checkpoint("Applying Forces");
+ if (multithread && !first) {
+ while (updateTriangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/ boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
+ if (updateTriangulation || ellapsedIter>(0.5*meshUpdateInterval)) {
+ if (useSolver==0) LOG_ERROR("background calculations not available for Gauss-Seidel");
+ if (fluidBulkModulus>0) solver->interpolate (solver->T[solver->currentTes], backgroundSolver->T[backgroundSolver->currentTes]);
+ solver=backgroundSolver;
+ backgroundSolver = shared_ptr<FlowSolver> (new FlowSolver);
+ //Copy imposed pressures/flow from the old solver
+ backgroundSolver->imposedP = vector<pair<CGT::Point,Real> >(solver->imposedP);
+ backgroundSolver->imposedF = vector<pair<CGT::Point,Real> >(solver->imposedF);
+ setPositionsBuffer(false);
+ cachedCell= Cell(*(scene->cell));
+ backgroundCompleted=false;
+ retriangulationLastIter=ellapsedIter;
+ ellapsedIter=0;
+ epsVolCumulative=0;
+ boost::thread workerThread(&PeriodicFlowEngine::backgroundAction,this);
+ workerThread.detach();
+ initializeVolumes(*solver);
+ computeViscousForces(*solver);
+ }
+ else if (debug && !first) {
+ if (debug && !backgroundCompleted) cerr<<"still computing solver in the background"<<endl;
+ ellapsedIter++;}
+ } else {
+ if (updateTriangulation && !first) {
+ cachedCell= Cell(*(scene->cell));
+ buildTriangulation (pZero, *solver);
+ initializeVolumes(*solver);
+ computeViscousForces(*solver);
+ updateTriangulation = false;
+ epsVolCumulative=0;
+ retriangulationLastIter=0;
+ ReTrg++;}
+ }
+ first=false;
+ timingDeltas->checkpoint("Ending");
+}
+
+
+// void PeriodicFlowEngine::backgroundAction()
+// {
+// if (useSolver<1) {LOG_ERROR("background calculations not available for Gauss-Seidel"); return;}
+// buildTriangulation (pZero,*backgroundSolver);
+// //FIXME: GS is computing too much, we need only matrix factorization in fact
+// backgroundSolver->gaussSeidel(scene->dt);
+// backgroundSolver->computeFacetForcesWithCache(/*onlyCache?*/ true);
+// // boost::this_thread::sleep(boost::posix_time::seconds(10));
+// backgroundCompleted = true;
+// }
+
+void PeriodicFlowEngine::triangulate( FlowSolver& flow )
+{
+ Tesselation& Tes = flow.T[flow.currentTes];
+ vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
+ FOREACH ( const posData& b, buffer ) {
+ if ( !b.exists || !b.isSphere || b.id==ignoredBody) continue;
+ Vector3i period; Vector3r wpos;
+ // FIXME: use "sheared" variant if the cell is sheared
+ wpos=cachedCell.wrapPt ( b.pos,period );
+ const Body::id_t& id = b.id;
+ const Real& rad = b.radius;
+ const Real& x = wpos[0];
+ const Real& y = wpos[1];
+ const Real& z = wpos[2];
+ VertexHandle vh0=Tes.insert ( x, y, z, rad, id );
+// VertexHandle vh0=Tes.insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );
+ if (vh0==NULL) {
+ flow.errorCode = 2;
+ LOG_ERROR("Vh NULL in PeriodicFlowEngine::triangulate(), check input data"); continue;}
+ for ( int k=0;k<3;k++ ) vh0->info().period[k]=-period[k];
+ const Vector3r cellSize ( cachedCell.getSize() );
+ //FIXME: if hasShear, comment in
+// wpos=scene->cell->unshearPt ( wpos );
+ // traverse all periodic cells around the body, to see if any of them touches
+ Vector3r halfSize= ( rad+duplicateThreshold ) *Vector3r ( 1,1,1 );
+ Vector3r pmin,pmax;
+ Vector3i i;
+ for ( i[0]=-1; i[0]<=1; i[0]++ )
+ for ( i[1]=-1;i[1]<=1; i[1]++ )
+ for ( i[2]=-1; i[2]<=1; i[2]++ ) {
+ if ( i[0]!=0 || i[1]!=0 || i[2]!=0 ) { // middle; already rendered above
+ Vector3r pos2=wpos+Vector3r ( cellSize[0]*i[0],cellSize[1]*i[1],cellSize[2]*i[2] ); // shift, but without shear!
+ pmin=pos2-halfSize;
+ pmax=pos2+halfSize;
+ if ( (pmin[0]<=cellSize[0]) && (pmax[0]>=0) && (pmin[1]<=cellSize[1]) && (pmax[1]>=0) && (pmin[2]<=cellSize[2]) && (pmax[2]>=0) ) {
+ //with shear:
+ //Vector3r pt=scene->cell->shearPt ( pos2 );
+ //without shear:
+ const Vector3r& pt= pos2;
+ VertexHandle vh=Tes.insert ( pt[0],pt[1],pt[2],rad,id,false,id );
+ for ( int k=0;k<3;k++ ) vh->info().period[k]=i[k]-period[k];}}
+ }
+ //re-assign the original vertex pointer since duplicates may have overwrite it
+ Tes.vertexHandles[id]=vh0;
+ }
+ Tes.redirected=true;//By inserting one-by-one, we already redirected
+ flow.shearLubricationForces.resize ( Tes.maxId+1 );
+ flow.shearLubricationTorques.resize ( Tes.maxId+1 );
+ flow.pumpLubricationTorques.resize ( Tes.maxId+1 );
+ flow.twistLubricationTorques.resize ( Tes.maxId+1 );
+ flow.shearLubricationBodyStress.resize ( Tes.maxId+1 );
+ flow.normalLubricationForce.resize ( Tes.maxId+1 );
+ flow.normalLubricationBodyStress.resize ( Tes.maxId+1 );
+}
+
+
+Real PeriodicFlowEngine::volumeCell ( CellHandle cell )
+{
+ static const Real inv6 = 1/6.;
+ const Vector3r p0 = positionBufferCurrent[cell->vertex(0)->info().id()].pos + makeVector3r(cell->vertex(0)->info().ghostShift());
+ const Vector3r p1 = positionBufferCurrent[cell->vertex(1)->info().id()].pos + makeVector3r(cell->vertex(1)->info().ghostShift());
+ const Vector3r p2 = positionBufferCurrent[cell->vertex(2)->info().id()].pos + makeVector3r(cell->vertex(2)->info().ghostShift());
+ const Vector3r p3 = positionBufferCurrent[cell->vertex(3)->info().id()].pos + makeVector3r(cell->vertex(3)->info().ghostShift());
+ Real volume = inv6*((p0-p1).cross(p0-p2)).dot(p0-p3);
+ if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
+ return volume;
+}
+
+Real PeriodicFlowEngine::volumeCellSingleFictious ( CellHandle cell )
+{
+ Vector3r V[3];
+ int b=0;
+ int w=0;
+ cell->info().volumeSign=1;
+ Real Wall_coordinate=0;
+
+ for ( int y=0;y<4;y++ ) {
+ if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( y )->info().id(), scene );
+ V[w]=sph->state->pos+ makeVector3r ( cell->vertex ( y )->info().ghostShift() );
+ w++;
+ } else {
+ b = cell->vertex ( y )->info().id();
+ const shared_ptr<Body>& wll = Body::byId ( b,scene );
+ if ( !solver->boundary ( b ).useMaxMin ) Wall_coordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wallThickness/2.;
+ else Wall_coordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate];
+ }
+ }
+ Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
+ return abs ( Volume );
+}
+
+
+void PeriodicFlowEngine::locateCell ( CellHandle baseCell, unsigned int& index, int& baseIndex, FlowSolver& flow, unsigned int count)
+{
+ if (count>10) {
+ LOG_ERROR("More than 10 attempts to locate a cell, duplicateThreshold may be too small, resulting in periodicity inconsistencies.");
+ flow.errorCode=1; return;
+ }
+ PeriFlowTesselation::CellInfo& baseInfo = baseCell->info();
+ //already located, return FIXME: is inline working correctly? else move this test outside the function, just before the calls
+ if ( baseInfo.index>0 || baseInfo.isGhost ) return;
+ RTriangulation& Tri = flow.T[flow.currentTes].Triangulation();
+ Vector3r center ( 0,0,0 );
+ Vector3i period;
+
+ if (baseCell->info().fictious()==0)
+ for ( int k=0;k<4;k++ ) center+= 0.25*makeVector3r (baseCell->vertex(k)->point());
+ else {
+
+ Real boundPos=0; int coord=0;
+ for ( int k=0;k<4;k++ ) {
+ if ( !baseCell->vertex ( k )->info().isFictious ) center+= 0.3333333333*makeVector3r ( baseCell->vertex ( k )->point() );
+ else {
+ coord=flow.boundary ( baseCell->vertex ( k )->info().id() ).coordinate;
+ boundPos=flow.boundary ( baseCell->vertex ( k )->info().id() ).p[coord];}
+ }
+ center[coord]=boundPos;
+ }
+ Vector3r wdCenter= cachedCell.wrapPt ( center,period );
+ if ( period[0]!=0 || period[1]!=0 || period[2]!=0 ) {
+ if ( baseCell->info().index>0 ) {
+ cout<<"indexed cell is found ghost!"<<baseInfo.index <<endl;
+ baseInfo.isGhost=false;
+ return;
+ }
+ CellHandle ch= Tri.locate ( CGT::Point ( wdCenter[0],wdCenter[1],wdCenter[2] )
+// ,/*hint*/ v0
+ );
+ baseInfo.period[0]=period[0];
+ baseInfo.period[1]=period[1];
+ baseInfo.period[2]=period[2];
+ //call recursively, since the returned cell could be also a ghost (especially if baseCell is a non-periodic type from the external contour
+ locateCell ( ch,index,baseIndex,flow,++count );
+ if ( ch==baseCell ) cerr<<"WTF!!"<<endl;
+ //check consistency
+ bool checkC=false;
+ for (int kk=0; kk<4;kk++) if ((!baseCell->vertex(kk)->info().isGhost) && ((!baseCell->vertex(kk)->info().isFictious))) checkC = true;
+ if (checkC) {
+ bool checkV=true;
+ for (int kk=0; kk<4;kk++) {
+ checkV=false;
+ for (int jj=0; jj<4;jj++)
+ if (baseCell->vertex(kk)->info().id() == ch->vertex(jj)->info().id()) checkV = true;
+ if (!checkV) {cerr <<"periodicity is broken"<<endl;
+ for (int jj=0; jj<4;jj++) cerr<<baseCell->vertex(jj)->info().id()<<" ";
+ cerr<<" vs. ";
+ for (int jj=0; jj<4;jj++) cerr<<ch->vertex(jj)->info().id()<<" ";
+ cerr<<endl;}
+ }
+ } else {
+// bool checkV=true;
+// for (int kk=0; kk<4;kk++) {
+// checkV=false;
+// for (int jj=0; jj<4;jj++)
+// if (baseCell->vertex(kk)->info().id() == ch->vertex(jj)->info().id()) checkV = true;
+// if (!checkV) {cerr <<"periodicity is broken (that's ok probably)"<<endl;
+// for (int jj=0; jj<4;jj++) cerr<<baseCell->vertex(jj)->info().id()<<" ";
+// cerr<<" vs. ";
+// for (int jj=0; jj<4;jj++) cerr<<ch->vertex(jj)->info().id()<<" ";
+// cerr<<endl;}
+// }
+ }
+
+ baseInfo.isGhost=true;
+ baseInfo._pression=& ( ch->info().p() );
+ baseInfo.index=ch->info().index;
+ baseInfo.baseIndex=ch->info().baseIndex;
+ baseInfo.Pcondition=ch->info().Pcondition;
+ } else {
+ baseInfo.isGhost=false;
+ //index is 1-based, if it is zero it is not initialized, we define it here
+ if ( baseInfo.baseIndex<0 ){
+ baseInfo.baseIndex=++baseIndex;
+ if (!baseInfo.Pcondition) baseInfo.index=++index;}
+ }
+}
+
+Vector3r PeriodicFlowEngine::meanVelocity()
+{
+ solver->averageRelativeCellVelocity();
+ Vector3r meanVel ( 0,0,0 );
+ Real volume=0;
+ FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
+ for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
+ //We could also define velocity using cell's center
+// if ( !cell->info().isReal() ) continue;
+ if ( cell->info().isGhost ) continue;
+ for ( int i=0;i<3;i++ )
+ meanVel[i]=meanVel[i]+ ( ( cell->info().averageVelocity() ) [i] * abs ( cell->info().volume() ) );
+ volume+=abs ( cell->info().volume() );
+ }
+ return ( meanVel/volume );
+}
+
+void PeriodicFlowEngine::updateVolumes (FlowSolver& flow)
+{
+ if ( debug ) cout << "Updating volumes.............." << endl;
+ Real invDeltaT = 1/scene->dt;
+ double newVol, dVol;
+ epsVolMax=0;
+ Real totVol=0;
+ Real totDVol=0;
+ Real totVol0=0;
+ Real totVol1=0;
+
+ FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
+ switch ( cell->info().fictious() ) {
+ case ( 1 ) :
+ newVol = volumeCellSingleFictious ( cell );
+ totVol1+=newVol;
+ break;
+ case ( 0 ) :
+ newVol = volumeCell ( cell );
+ totVol0+=newVol;
+ break;
+ default:
+ newVol = 0;
+ break;
+ }
+ totVol+=newVol;
+ dVol=cell->info().volumeSign * ( newVol - cell->info().volume() );
+ totDVol+=dVol;
+ epsVolMax = max ( epsVolMax, abs ( dVol/newVol ) );
+ cell->info().dv() = dVol * invDeltaT;
+ cell->info().volume() = newVol;
+ }
+ if ( debug ) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<" "<< totVol0<<" "<< totVol1<<endl;
+}
+
+
+void PeriodicFlowEngine::initializeVolumes (FlowSolver& flow)
+{
+ FiniteVerticesIterator vertices_end = flow.T[flow.currentTes].Triangulation().finite_vertices_end();
+ CGT::CVector Zero ( 0,0,0 );
+ for ( FiniteVerticesIterator V_it = flow.T[flow.currentTes].Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++ ) V_it->info().forces=Zero;
+
+ FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
+ switch ( cell->info().fictious() )
+ {
+ case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
+ case ( 1 ) : cell->info().volume() = volumeCellSingleFictious ( cell ); break;
+ default: cell->info().volume() = 0; break;
+ }
+ //FIXME: the void volume is negative sometimes, hence crashing...
+ if (flow.fluidBulkModulus>0) { cell->info().invVoidVolume() = 1. / (max(0.1*cell->info().volume(),abs(cell->info().volume()) - flow.volumeSolidPore(cell)) ); }
+ }
+ if ( debug ) cout << "Volumes initialised." << endl;
+}
+
+void PeriodicFlowEngine::buildTriangulation ( double pZero, FlowSolver& flow)
+{
+ flow.resetNetwork();
+ if (first) flow.currentTes=0;
+ else {
+ flow.currentTes=!flow.currentTes;
+ if ( debug ) cout << "--------RETRIANGULATION-----------" << endl;}
+ initSolver(flow);
+ addBoundary ( flow );
+ if ( debug ) cout << endl << "Added boundaries------" << endl << endl;
+ triangulate (flow);
+ if ( debug ) cout << endl << "Tesselating------" << endl << endl;
+ flow.T[flow.currentTes].compute();
+ flow.defineFictiousCells();
+
+ //FIXME: this is already done in addBoundary(?)
+ boundaryConditions ( flow );
+ if ( debug ) cout << endl << "boundaryConditions------" << endl << endl;
+ flow.initializePressure ( pZero );
+ if ( debug ) cout << endl << "initializePressure------" << endl << endl;
+ // Define the ghost cells and add indexes to the cells inside the period (the ones that will contain the pressure unknowns)
+ //This must be done after boundary conditions and initialize pressure, else the indexes are not good (not accounting imposedP): FIXME
+ unsigned int index=0;
+ int baseIndex=-1;
+ FlowSolver::Tesselation& Tes = flow.T[flow.currentTes];
+ Tes.cellHandles.resize(Tes.Triangulation().number_of_finite_cells());
+ const FiniteCellsIterator cellend=Tes.Triangulation().finite_cells_end();
+ for ( FiniteCellsIterator cell=Tes.Triangulation().finite_cells_begin(); cell!=cellend; cell++ ){
+ locateCell ( cell,index,baseIndex,flow );
+ if (flow.errorCode>0) return;
+ //Fill this vector than can be later used to speedup loops
+ if (!cell->info().isGhost) Tes.cellHandles[cell->info().baseIndex]=cell;
+ }
+ Tes.cellHandles.resize(baseIndex+1);
+
+ if ( debug ) cout << endl << "locateCell------" << endl << endl;
+ flow.computePermeability ( );
+ porosity = flow.vPoralPorosity/flow.vTotalPorosity;
+ flow.tolerance=tolerance;flow.relax=relax;
+
+ flow.displayStatistics ();
+ //FIXME: check interpolate() for the periodic case, at least use the mean pressure from previous step.
+ if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], Tes );
+// if ( !first && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
+
+ if ( waveAction ) flow.applySinusoidalPressure ( Tes.Triangulation(), sineMagnitude, sineAverage, 30 );
+
+ if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
+ if ( debug ) cout << endl << "end buildTri------" << endl << endl;
+}
+
+void PeriodicFlowEngine::preparePShifts()
+{
+ CellInfo::gradP = makeCgVect ( gradP );
+ CellInfo::hSize[0] = makeCgVect ( scene->cell->hSize.col ( 0 ) );
+ CellInfo::hSize[1] = makeCgVect ( scene->cell->hSize.col ( 1 ) );
+ CellInfo::hSize[2] = makeCgVect ( scene->cell->hSize.col ( 2 ) );
+ CellInfo::deltaP=CGT::CVector (
+ CellInfo::hSize[0]*CellInfo::gradP,
+ CellInfo::hSize[1]*CellInfo::gradP,
+ CellInfo::hSize[2]*CellInfo::gradP );
+}
+
+
+YADE_PLUGIN((PeriodicFlowEngine));
+
+
+#endif //FLOW_ENGINE
+
+#endif /* YADE_CGAL */
\ No newline at end of file
=== added file 'pkg/pfv/PeriodicFlowEngine.hpp'
--- pkg/pfv/PeriodicFlowEngine.hpp 1970-01-01 00:00:00 +0000
+++ pkg/pfv/PeriodicFlowEngine.hpp 2014-04-03 13:37:55 +0000
@@ -0,0 +1,109 @@
+/*************************************************************************
+* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
+* Copyright (C) 2012 by Donia Marzougui <donia.marzougui@xxxxxxxxxxxxxxx>*
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+/// The periodic variant of FlowEngine is defined here. It should become a template class for more flexibility.
+/// It is a bit more complicated as for FlowEngine, though, because we need template inheriting from template, which breaks YADE_CLASS_XXX logic_error
+/// See below the commented exemple, for a possible solution
+
+#include <yade/pkg/dem/FlowEngine.hpp>
+
+class PeriodicCellInfo : public FlowCellInfo
+{
+ public:
+ static CVector gradP;
+ //for real cell, baseIndex is the rank of the cell in cellHandles. For ghost cells, it is the baseIndex of the corresponding real cell.
+ //Unlike ordinary index, baseIndex is also indexing cells with imposed pressures
+ int baseIndex;
+ int period[3];
+ static CVector hSize[3];
+ static CVector deltaP;
+ int ghost;
+ Real* _pression;
+ PeriodicCellInfo (void){
+ _pression=&pression;
+ period[0]=period[1]=period[2]=0;
+ baseIndex=-1;
+ volumeSign=0;}
+ ~PeriodicCellInfo (void) {}
+
+ inline const double shiftedP (void) const {return isGhost? (*_pression)+pShift() :(*_pression) ;}
+ inline const double pShift (void) const {return deltaP[0]*period[0] + deltaP[1]*period[1] +deltaP[2]*period[2];}
+// inline const double p (void) {return shiftedP();}
+ inline void setP (const Real& p) {pression=p;}
+ bool isReal (void) {return !(isFictious || isGhost);}
+};
+
+
+class PeriodicVertexInfo : public FlowVertexInfo {
+ public:
+ PeriodicVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
+ PeriodicVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
+ PeriodicVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
+ int period[3];
+ //FIXME: the name is misleading, even non-ghost can be out of the period and therefore they need to be shifted as well
+ inline const CVector ghostShift (void) {
+ return period[0]*PeriodicCellInfo::hSize[0]+period[1]*PeriodicCellInfo::hSize[1]+period[2]*PeriodicCellInfo::hSize[2];}
+ PeriodicVertexInfo (void) {isFictious=false; s=0; i=0; period[0]=period[1]=period[2]=0; isGhost=false;}
+ bool isReal (void) {return !(isFictious || isGhost);}
+};
+
+typedef CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> PeriFlowTriangulationTypes;
+typedef CGT::PeriodicTesselation<CGT::_Tesselation<PeriFlowTriangulationTypes> > PeriFlowTesselation;
+#ifdef LINSOLV
+#define _PeriFlowSolver CGT::PeriodicFlowLinSolv<PeriFlowTesselation>
+#else
+#define _PeriFlowSolver CGT::PeriodicFlow<PeriFlowTesselation>
+#endif
+
+typedef TemplateFlowEngine< PeriodicCellInfo,
+ PeriodicVertexInfo,
+ CGT::PeriodicTesselation<CGT::_Tesselation<CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> > >,
+ _PeriFlowSolver>
+ FlowEngine_PeriodicInfo;
+
+REGISTER_SERIALIZABLE(FlowEngine_PeriodicInfo);
+YADE_PLUGIN((FlowEngine_PeriodicInfo));
+
+
+class PeriodicFlowEngine : public FlowEngine_PeriodicInfo
+{
+ public :
+ void triangulate (FlowSolver& flow);
+ void buildTriangulation (Real pzero, FlowSolver& flow);
+ void initializeVolumes (FlowSolver& flow);
+ void updateVolumes (FlowSolver& flow);
+ Real volumeCell (CellHandle cell);
+
+ Real volumeCellSingleFictious (CellHandle cell);
+ inline void locateCell(CellHandle baseCell, unsigned int& index, int& baseIndex, FlowSolver& flow, unsigned int count=0);
+ Vector3r meanVelocity();
+
+ virtual ~PeriodicFlowEngine();
+
+ virtual void action();
+ //Cache precomputed values for pressure shifts, based on current hSize and pGrad
+ void preparePShifts();
+
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(PeriodicFlowEngine,FlowEngine_PeriodicInfo,"A variant of :yref:`FlowEngine` implementing periodic boundary conditions. The API is very similar.",
+ ((Real,duplicateThreshold, 0.06,,"distance from cell borders that will triger periodic duplication in the triangulation |yupdate|"))
+ ((Vector3r, gradP, Vector3r::Zero(),,"Macroscopic pressure gradient"))
+ ,,
+ wallIds=vector<int>(6,-1);
+ solver = shared_ptr<FlowSolver> (new FlowSolver);
+ epsVolMax=epsVolCumulative=retriangulationLastIter=0;
+ ReTrg=1;
+ first=true;
+ ,
+ //nothing special to define, we re-use FlowEngine methods
+ //.def("meanVelocity",&PeriodicFlowEngine::meanVelocity,"measure the mean velocity in the period")
+ )
+ DECLARE_LOGGER;
+
+
+};
+REGISTER_SERIALIZABLE(PeriodicFlowEngine);
\ No newline at end of file
=== added file 'pkg/pfv/SoluteFlowEngine.cpp'
--- pkg/pfv/SoluteFlowEngine.cpp 1970-01-01 00:00:00 +0000
+++ pkg/pfv/SoluteFlowEngine.cpp 2014-04-03 13:37:55 +0000
@@ -0,0 +1,197 @@
+/*************************************************************************
+* Copyright (C) 2013 by T. Sweijen (T.sweijen@xxxxx) *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#ifdef YADE_CGAL
+#ifdef FLOW_ENGINE
+
+#define SOLUTE_FLOW
+#ifdef SOLUTE_FLOW
+
+#define TEMPLATE_FLOW_NAME SoluteFlowEngineT
+#include <yade/pkg/dem/FlowEngine.hpp>
+#undef TEMPLATE_FLOW_NAME
+
+class SoluteCellInfo : public FlowCellInfo
+{
+ public:
+ Real solute_concentration;
+ SoluteCellInfo (void) : FlowCellInfo() {solute_concentration=0;}
+ inline Real& solute (void) {return solute_concentration;}
+ inline const Real& solute (void) const {return solute_concentration;}
+ inline void getInfo (const SoluteCellInfo& otherCellInfo) {FlowCellInfo::getInfo(otherCellInfo); solute()=otherCellInfo.solute();}
+};
+
+typedef TemplateFlowEngine<SoluteCellInfo,FlowVertexInfo> SoluteFlowEngineT;
+REGISTER_SERIALIZABLE(SoluteFlowEngineT);
+YADE_PLUGIN((SoluteFlowEngineT));
+
+class SoluteFlowEngine : public SoluteFlowEngineT
+{
+ public :
+ void initializeSoluteTransport();
+ void soluteTransport (double deltatime, double D);
+ double getConcentration(unsigned int id){return solver->T[solver->currentTes].cellHandles[id]->info().solute(); }
+ double insertConcentration(unsigned int id,double conc){
+ solver->T[solver->currentTes].cellHandles[id]->info().solute() = conc;
+ return conc;}
+ void soluteBC(unsigned int bc_id1, unsigned int bc_id2, double bc_concentration1, double bc_concentration2,unsigned int s);
+ double getConcentrationPlane (double Yobs,double Yr, int xyz);
+ ///Elaborate the description as you wish
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(SoluteFlowEngine,SoluteFlowEngineT,"A variant of :yref:`FlowEngine` with solute transport).",
+ ///No additional variable yet, else input here
+// ((Vector3r, gradP, Vector3r::Zero(),,"Macroscopic pressure gradient"))
+ ,,,
+ .def("soluteTransport",&SoluteFlowEngine::soluteTransport,(python::arg("deltatime"),python::arg("D")),"Solute transport (advection and diffusion) engine for diffusion use a diffusion coefficient (D) other than 0.")
+ .def("getConcentration",&SoluteFlowEngine::getConcentration,(python::arg("id")),"get concentration of pore with ID")
+ .def("insertConcentration",&SoluteFlowEngine::insertConcentration,(python::arg("id"),python::arg("conc")),"Insert Concentration (ID, Concentration)")
+ .def("solute_BC",&SoluteFlowEngine::soluteBC,(python::arg("bc_id1"),python::arg("bc_id2"),python::arg("bc_concentration1"),python::arg("bc_concentration2"),python::arg("s")),"Enter X,Y,Z for concentration observation'.")
+ //I guess there is getConcentrationPlane missing here, but it is not on github
+ )
+};
+REGISTER_SERIALIZABLE(SoluteFlowEngine);
+
+
+// PeriodicFlowEngine::~PeriodicFlowEngine(){}
+
+void SoluteFlowEngine::initializeSoluteTransport ()
+{
+ //Prepare a list with concentrations in range of 0,nCells. Or resets the list of concentrations to 0
+// typedef typename Solver::element_type Flow;
+// typedef typename Flow::Finite_vertices_iterator Finite_vertices_iterator;
+// typedef typename Solver::element_type Flow;
+
+ FOREACH(CellHandle& cell, solver->T[solver->currentTes].cellHandles)
+ {
+ cell->info().solute() = 0.0;
+ }
+}
+
+void SoluteFlowEngine::soluteTransport (double deltatime, double D)
+{
+ //soluteTransport is a function to solve transport of solutes for advection (+diffusion).
+ //Call this function with a pyRunner in the python script, implement a diffusion coefficient and a dt.
+ // Optimalization has to be done such that the coefficient matrix is not solved for each time step,only after triangulation.
+ // Extensive testing has to be done to check its ability to simulate a deforming porous media.
+ double coeff = 0.00;
+ double coeff1 = 0.00; //Coefficient for off-diagonal element
+ double coeff2 = 0.00; //Coefficient for diagonal element
+ double qin = 0.00; //Flux into the pore per pore throat
+ double Qout=0.0; //Total Flux out of the pore
+ int ncells = 0; //Number of cells
+ int ID = 0;
+ //double D = 0.0;
+ double invdistance = 0.0; //Fluid facet area divided by pore throat length for each pore throat
+ double invdistancelocal = 0.0; //Sum of invdistance
+ ncells=solver->T[solver->currentTes].cellHandles.size();
+
+ Eigen::SparseMatrix<double, Eigen::ColMajor> Aconc(ncells,ncells);
+ typedef Eigen::Triplet<double> ETriplet2;
+ std::vector<ETriplet2> tripletList2;
+ Eigen::SparseLU<Eigen::SparseMatrix<double,Eigen::ColMajor>,Eigen::COLAMDOrdering<int> > eSolver2;
+
+ // Prepare (copy) concentration vector
+ Eigen::VectorXd eb2(ncells); Eigen::VectorXd ex2(ncells);
+ FOREACH(CellHandle& cell, solver->T[solver->currentTes].cellHandles){
+ eb2[cell->info().id]=cell->info().solute();
+ }
+
+ // Fill coefficient matrix
+ FOREACH(CellHandle& cell, solver->T[solver->currentTes].cellHandles){
+ cell->info().invVoidVolume() = 1 / ( abs(cell->info().volume()) - abs(solver->volumeSolidPore(cell) ) );
+ invdistance=0.0;
+
+
+ unsigned int i=cell->info().id;
+ for (unsigned int ngb=0;ngb<4;ngb++){
+ CGT::Point& p2 = cell ->neighbor(ngb)->info();
+ CGT::Point& p1 = cell->info();
+ CGT::CVector l = p1-p2;
+ CGT::Real fluidSurf = sqrt(cell->info().facetSurfaces[ngb].squared_length())*cell->info().facetFluidSurfacesRatio[ngb];
+ invdistancelocal = (fluidSurf/sqrt(l.squared_length()));
+ invdistance+=(fluidSurf/sqrt(l.squared_length()));
+ coeff = deltatime*cell->info().invVoidVolume();
+ ID = cell->neighbor(ngb)->info().id;
+ qin=abs(cell->info().kNorm() [ngb])* ( cell->neighbor ( ngb )->info().p()-cell->info().p());
+ Qout=Qout+max(qin,0.0);
+ coeff1=-1*coeff*(abs(max(qin,0.0))-(D*invdistancelocal));
+ if (coeff1 != 0.0){
+ tripletList2.push_back(ETriplet2(i,ID,coeff1));
+ }
+
+ }
+ coeff2=1.0+(coeff*abs(Qout))+(coeff*D*invdistance);
+ tripletList2.push_back(ETriplet2(i,i,coeff2));
+ Qout=0.0;
+ }
+ //Solve Matrix
+ Aconc.setFromTriplets(tripletList2.begin(), tripletList2.end());
+ //if (eSolver2.signDeterminant() < 1){cerr << "determinant is negative!!!!!!! " << eSolver2.signDeterminant()<<endl;}
+ //eSolver2.setPivotThreshold(10e-8);
+ eSolver2.analyzePattern(Aconc);
+ eSolver2.factorize(Aconc);
+ eSolver2.compute(Aconc);
+ ex2 = eSolver2.solve(eb2);
+
+
+ double averageConc = 0.0;
+ FOREACH(CellHandle& cell, solver->T[solver->currentTes].cellHandles){
+ averageConc+=cell->info().solute();}
+
+ //Copy data to concentration array
+ FOREACH(CellHandle& cell, solver->T[solver->currentTes].cellHandles){
+ cell->info().solute()= ex2[cell->info().id];
+
+ }
+
+ tripletList2.clear();
+
+ }
+
+void SoluteFlowEngine::soluteBC(unsigned int bcid1, unsigned int bcid2, double bcconcentration1, double bcconcentration2, unsigned int s)
+{
+ //Boundary conditions according to soluteTransport.
+ //It simply assigns boundary concentrations to cells with a common vertices (e.g. infinite large sphere which makes up the boundary condition in flowEngine)
+ //s is a switch, if 0 only bc_id1 is used (advection only). If >0 than both bc_id1 and bc_id2 are used.
+ FOREACH(CellHandle& cell, solver->T[solver->currentTes].cellHandles)
+ {
+ for (unsigned int ngb=0;ngb<4;ngb++){
+ if (cell->vertex(ngb)->info().id() == bcid1){cell->info().solute() = bcconcentration1;}
+ if (s > 0){if (cell->vertex(ngb)->info().id() == bcid2){cell->info().solute() = bcconcentration2;}}
+ }
+ }
+}
+
+double SoluteFlowEngine::getConcentrationPlane (double Yobs,double Yr, int xyz)
+{
+ //Get the concentration within a certain plane (Y_obs), whilst the cells are located in a small volume around this plane
+ //The concentration in cells are weighed for their distance to the observation point
+ //for a point on the x-axis (xyz=0), point on the y-axis (xyz=1), point on the z-axis (xyz=2)
+ double sumConcentration = 0.0;
+ double sumFraction=0.0;
+ double concentration=0.0;
+ //Find cells within designated volume
+
+ FOREACH(CellHandle& cell, solver->T[solver->currentTes].cellHandles)
+ {
+ CGT::Point& p1 = cell->info();
+ if (abs(p1[xyz]) < abs(abs(Yobs) + abs(Yr))){
+ if(abs(p1[xyz]) > abs(abs(Yobs) - abs(Yr))){
+ sumConcentration += cell->info().solute()*(1-(abs(p1[xyz])-abs(Yobs))/abs(Yr));
+ sumFraction += (1-(abs(p1[xyz])-abs(Yobs))/abs(Yr));
+ }
+ }
+ }
+ concentration = sumConcentration / sumFraction;
+ return concentration;
+}
+
+YADE_PLUGIN ( ( SoluteFlowEngine ) );
+
+#endif //SOLUTE_FLOW
+#endif //FLOW_ENGINE
+
+#endif /* YADE_CGAL */
\ No newline at end of file
