yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10989
[Branch ~yade-pkg/yade/git-trunk] Rev 4033: PFV code refactoring.
------------------------------------------------------------
revno: 4033
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Fri 2014-06-20 17:28:21 +0200
message:
PFV code refactoring.
Implement generaton of hpp-files for each cpp to escape
function redifinitions.
http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg10232.html
added:
pkg/pfv/FlowEngine.hpp.in
pkg/pfv/FlowEngine.ipp.in
modified:
CMakeLists.txt
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/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
=== modified file 'CMakeLists.txt'
--- CMakeLists.txt 2014-06-17 10:20:23 +0000
+++ CMakeLists.txt 2014-06-20 15:28:21 +0000
@@ -499,6 +499,17 @@
CONFIGURE_FILE(py/__init__.py.in "${CMAKE_BINARY_DIR}/__init__.py")
#===========================================================
+# Create header files for PFV from FlowEngine.hpp.in-template.
+# All @TEMPLATE_FLOW_NAME@ are replacing by a given names
+
+SET (TEMPLATE_FLOW_NAMES DFNFlowEngineT DummyFlowEngineT FlowEngineT FlowEngine_PeriodicInfo SoluteFlowEngineT)
+FOREACH(TF ${TEMPLATE_FLOW_NAMES})
+ SET (TEMPLATE_FLOW_NAME ${TF})
+ CONFIGURE_FILE(pkg/pfv/FlowEngine.hpp.in "${CMAKE_BINARY_DIR}/pkg/pfv/FlowEngine_${TF}.hpp" @ONLY)
+ CONFIGURE_FILE(pkg/pfv/FlowEngine.ipp.in "${CMAKE_BINARY_DIR}/pkg/pfv/FlowEngine_${TF}.ipp" @ONLY)
+ENDFOREACH(TF)
+INCLUDE_DIRECTORIES("${CMAKE_BINARY_DIR}/pkg/pfv/")
+#===========================================================
INSTALL(PROGRAMS "${CMAKE_BINARY_DIR}/bins/yade${SUFFIX}-batch" DESTINATION ${YADE_EXEC_PATH}/)
INSTALL(PROGRAMS "${CMAKE_BINARY_DIR}/bins/yade${SUFFIX}" DESTINATION ${YADE_EXEC_PATH}/)
=== modified file 'pkg/pfv/DFNFlow.cpp'
--- pkg/pfv/DFNFlow.cpp 2014-05-27 10:33:52 +0000
+++ pkg/pfv/DFNFlow.cpp 2014-06-20 15:28:21 +0000
@@ -11,11 +11,11 @@
//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/pfv/FlowEngine.hpp>
+#include "FlowEngine_DFNFlowEngineT.hpp"
-class DFNCellInfo : public FlowCellInfo
+class DFNCellInfo : public FlowCellInfo_DFNFlowEngineT
{
public:
Real anotherVariable;
@@ -23,12 +23,12 @@
};
-class DFNVertexInfo : public FlowVertexInfo {
+class DFNVertexInfo : public FlowVertexInfo_DummyFlowEngineT {
public:
//same here if needed
};
-typedef TemplateFlowEngine<DFNCellInfo,DFNVertexInfo> DFNFlowEngineT;
+typedef TemplateFlowEngine_DFNFlowEngineT<DFNCellInfo,DFNVertexInfo> DFNFlowEngineT;
REGISTER_SERIALIZABLE(DFNFlowEngineT);
YADE_PLUGIN((DFNFlowEngineT));
@@ -104,4 +104,4 @@
#endif //DFNFLOW
-#endif //FLOWENGINE
\ No newline at end of file
+#endif //FLOWENGINE
=== modified file 'pkg/pfv/DummyFlowEngine.cpp'
--- pkg/pfv/DummyFlowEngine.cpp 2014-05-23 13:20:43 +0000
+++ pkg/pfv/DummyFlowEngine.cpp 2014-06-20 15:28:21 +0000
@@ -13,27 +13,27 @@
//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/pfv/FlowEngine.hpp>
+
+#include "FlowEngine_DummyFlowEngineT.hpp"
/// We can add data to the Info types by inheritance
-class DummyCellInfo : public FlowCellInfo
+class DummyCellInfo : public FlowCellInfo_DummyFlowEngineT
{
public:
Real anotherVariable;
void anotherFunction() {};
};
-class DummyVertexInfo : public FlowVertexInfo {
+class DummyVertexInfo : public FlowVertexInfo_DummyFlowEngineT {
public:
//same here if needed
};
-typedef TemplateFlowEngine<DummyCellInfo,DummyVertexInfo> TEMPLATE_FLOW_NAME;
-REGISTER_SERIALIZABLE(TEMPLATE_FLOW_NAME);
-YADE_PLUGIN((TEMPLATE_FLOW_NAME));
+typedef TemplateFlowEngine_DummyFlowEngineT<DummyCellInfo,DummyVertexInfo> DummyFlowEngineT;
+REGISTER_SERIALIZABLE(DummyFlowEngineT);
+YADE_PLUGIN((DummyFlowEngineT));
-class DummyFlowEngine : public TEMPLATE_FLOW_NAME
+class DummyFlowEngine : public DummyFlowEngineT
{
public :
//We can overload every functions of the base engine to make it behave differently
@@ -42,9 +42,9 @@
//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;}
+ void fancyFunction(Real what);
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DummyFlowEngine,TEMPLATE_FLOW_NAME,"documentation here",
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DummyFlowEngine,DummyFlowEngineT,"documentation here",
((Real, myNewAttribute, 0,,"useless example"))
,/*DummyFlowEngineT()*/,
,
@@ -55,6 +55,5 @@
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
+void DummyFlowEngine::fancyFunction(Real what) {std::cerr<<"yes, I'm a new function"<<std::endl;}
#endif //DummyFLOW
=== modified file 'pkg/pfv/FlowEngine.cpp'
--- pkg/pfv/FlowEngine.cpp 2014-04-03 13:37:55 +0000
+++ pkg/pfv/FlowEngine.cpp 2014-06-20 15:28:21 +0000
@@ -8,13 +8,12 @@
#ifdef YADE_CGAL
#ifdef FLOW_ENGINE
-#define TEMPLATE_FLOW_NAME FlowEngineT
-#include "FlowEngine.hpp"
+#include "FlowEngine_FlowEngineT.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;
+typedef TemplateFlowEngine_FlowEngineT<FlowCellInfo_FlowEngineT,FlowVertexInfo_DummyFlowEngineT> FlowEngineT;
REGISTER_SERIALIZABLE(FlowEngineT);
class FlowEngine : public FlowEngineT
=== modified file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp 2014-06-04 17:19:50 +0000
+++ pkg/pfv/FlowEngine.hpp 2014-06-20 15:28:21 +0000
@@ -1,45 +1,4 @@
-/*************************************************************************
-* 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;
@@ -50,391 +9,8 @@
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
-
#ifdef LINSOLV
#define DEFAULTSOLVER CGT::FlowBoundingSphereLinSolv<_Tesselation>
#else
#define DEFAULTSOLVER CGT::FlowBoundingSphere<_Tesselation>
#endif
-
-template<class _CellInfo, class _VertexInfo, class _Tesselation=CGT::_Tesselation<CGT::TriangulationTypes<_VertexInfo,_CellInfo> >, class solverT=DEFAULTSOLVER >
-class TemplateFlowEngine : public PartialEngine
-{
- public :
- typedef solverT FlowSolver;
- 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
- virtual 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 () {
- if (solver->tesselation().maxId>0) boundaryConditions(*solver);//avoids crash at iteration 0, when the packing is not bounded yet
- else LOG_ERROR("updateBCs not applied");
- solver->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 averageVelocity();
-
- 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]);}
- boost::python::list getConstrictions(bool all) {
- vector<Real> csd=solver->getConstrictions(); boost::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;}
- boost::python::list getConstrictionsFull(bool all) {
- vector<Constriction> csd=solver->getConstrictionsFull(); boost::python::list pycsd;
- for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k].second[0]!=0) || csd[k].second[0]>0) {
- boost::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();}
- boost::python::list getVertices(unsigned int id){
- boost::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"))
- ((bool,doInterpolate,false,,"Force the interpolation of cell's info while remeshing. By default, interpolation would be done only for compressible fluids. It can be forced with this flag."))
- ((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
- ((vector<Real>, boundaryPressure,vector<Real>(),,"values defining pressure along x-axis for the top surface. See also :yref:`TEMPLATE_FLOW_NAME::boundaryXPos`"))
- ((vector<Real>, boundaryXPos,vector<Real>(),,"values of the x-coordinate for which pressure is defined. See also :yref:`TEMPLATE_FLOW_NAME::boundaryPressure`"))
- ,
- /*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,(boost::python::arg("pos"),boost::python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
- .def("imposePressure",&TemplateFlowEngine::imposePressure,(boost::python::arg("pos"),boost::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,(boost::python::arg("cond"),boost::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,(boost::python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
- .def("getBoundaryFlux",&TemplateFlowEngine::getBoundaryFlux,(boost::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,(boost::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,(boost::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,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
- .def("avFlVelOnSph",&TemplateFlowEngine::avFlVelOnSph,(boost::python::arg("idSph")),"compute a sphere-centered average fluid velocity")
- .def("fluidForce",&TemplateFlowEngine::fluidForce,(boost::python::arg("idSph")),"Return the fluid force on sphere idSph.")
- .def("shearLubForce",&TemplateFlowEngine::shearLubForce,(boost::python::arg("idSph")),"Return the shear lubrication force on sphere idSph.")
- .def("shearLubTorque",&TemplateFlowEngine::shearLubTorque,(boost::python::arg("idSph")),"Return the shear lubrication torque on sphere idSph.")
- .def("normalLubForce",&TemplateFlowEngine::normalLubForce,(boost::python::arg("idSph")),"Return the normal lubrication force on sphere idSph.")
- .def("bodyShearLubStress",&TemplateFlowEngine::bodyShearLubStress,(boost::python::arg("idSph")),"Return the shear lubrication stress on sphere idSph.")
- .def("bodyNormalLubStress",&TemplateFlowEngine::bodyNormalLubStress,(boost::python::arg("idSph")),"Return the normal lubrication stress on sphere idSph.")
- .def("shearVelocity",&TemplateFlowEngine::shearVelocity,(boost::python::arg("idSph")),"Return the shear velocity of the interaction.")
- .def("normalVelocity",&TemplateFlowEngine::normalVelocity,(boost::python::arg("idSph")),"Return the normal velocity of the interaction.")
- .def("normalVect",&TemplateFlowEngine::normalVect,(boost::python::arg("idSph")),"Return the normal vector between particles.")
- .def("surfaceDistanceParticle",&TemplateFlowEngine::surfaceDistanceParticle,(boost::python::arg("interaction")),"Return the distance between particles.")
- .def("onlySpheresInteractions",&TemplateFlowEngine::onlySpheresInteractions,(boost::python::arg("interaction")),"Return the id of the interaction only between spheres.")
- .def("pressureProfile",&TemplateFlowEngine::pressureProfile,(boost::python::arg("wallUpY"),boost::python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
- .def("getPorePressure",&TemplateFlowEngine::getPorePressure,(boost::python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
- .def("averageSlicePressure",&TemplateFlowEngine::averageSlicePressure,(boost::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,(boost::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,(boost::python::arg("id")),"get the vertices of a cell")
- #ifdef LINSOLV
- .def("exportMatrix",&TemplateFlowEngine::exportMatrix,(boost::python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
- .def("exportTriplets",&TemplateFlowEngine::exportTriplets,(boost::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,(boost::python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
- .def("averageVelocity",&TemplateFlowEngine::averageVelocity,"measure the mean velocity in the period")
- )
-};
-// Definition of functions in a separate file for clarity
-#include<yade/pkg/pfv/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.doInterpolate
- double solidSurfaces [4][4];
-
- FlowCellInfo (void)
- {
- modulePermeability.resize(4, 0);
- cellForce.resize(4,CGAL::NULL_VECTOR);
- facetSurfaces.resize(4,CGAL::NULL_VECTOR);
- facetFluidSurfacesRatio.resize(4,0);
- facetSphereCrossSections.resize(4,CGAL::NULL_VECTOR);
- unitForceVectors.resize(4,CGAL::NULL_VECTOR);
- 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) const {return CGAL::NULL_VECTOR;}
-};
-
-
-
=== added file 'pkg/pfv/FlowEngine.hpp.in'
--- pkg/pfv/FlowEngine.hpp.in 1970-01-01 00:00:00 +0000
+++ pkg/pfv/FlowEngine.hpp.in 2014-06-20 15:28:21 +0000
@@ -0,0 +1,416 @@
+/*************************************************************************
+* 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>
+#include<yade/pkg/pfv/FlowEngine.hpp>
+
+template<class _CellInfo, class _VertexInfo, class _Tesselation=CGT::_Tesselation<CGT::TriangulationTypes<_VertexInfo,_CellInfo> >, class solverT=DEFAULTSOLVER >
+class TemplateFlowEngine_@TEMPLATE_FLOW_NAME@ : public PartialEngine
+{
+ public :
+ typedef solverT FlowSolver;
+ 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
+ virtual 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 () {
+ if (solver->tesselation().maxId>0) boundaryConditions(*solver);//avoids crash at iteration 0, when the packing is not bounded yet
+ else LOG_ERROR("updateBCs not applied");
+ solver->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 averageVelocity();
+
+ 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]);}
+ boost::python::list getConstrictions(bool all) {
+ vector<Real> csd=solver->getConstrictions(); boost::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;}
+ boost::python::list getConstrictionsFull(bool all) {
+ vector<Constriction> csd=solver->getConstrictionsFull(); boost::python::list pycsd;
+ for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k].second[0]!=0) || csd[k].second[0]>0) {
+ boost::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();}
+ boost::python::list getVertices(unsigned int id){
+ boost::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_@TEMPLATE_FLOW_NAME@();
+ 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@,@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"))
+ ((bool,doInterpolate,false,,"Force the interpolation of cell's info while remeshing. By default, interpolation would be done only for compressible fluids. It can be forced with this flag."))
+ ((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
+ ((vector<Real>, boundaryPressure,vector<Real>(),,"values defining pressure along x-axis for the top surface. See also :yref:`@TEMPLATE_FLOW_NAME@::boundaryXPos`"))
+ ((vector<Real>, boundaryXPos,vector<Real>(),,"values of the x-coordinate for which pressure is defined. See also :yref:`@TEMPLATE_FLOW_NAME@::boundaryPressure`"))
+ ,
+ /*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_@TEMPLATE_FLOW_NAME@::imposeFlux,(boost::python::arg("pos"),boost::python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
+ .def("imposePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::imposePressure,(boost::python::arg("pos"),boost::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_@TEMPLATE_FLOW_NAME@::setImposedPressure,(boost::python::arg("cond"),boost::python::arg("p")),"Set pressure value at the point indexed 'cond'.")
+ .def("clearImposedPressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::clearImposedPressure,"Clear the list of points with pressure imposed.")
+ .def("clearImposedFlux",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::clearImposedFlux,"Clear the list of points with flux imposed.")
+ .def("getCellFlux",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getCellFlux,(boost::python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
+ .def("getBoundaryFlux",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getBoundaryFlux,(boost::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_@TEMPLATE_FLOW_NAME@::getConstrictions,(boost::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_@TEMPLATE_FLOW_NAME@::getConstrictionsFull,(boost::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_@TEMPLATE_FLOW_NAME@::edgeSize,"Return the number of interactions.")
+ .def("OSI",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::OSI,"Return the number of interactions only between spheres.")
+
+ .def("saveVtk",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
+ .def("avFlVelOnSph",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::avFlVelOnSph,(boost::python::arg("idSph")),"compute a sphere-centered average fluid velocity")
+ .def("fluidForce",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::fluidForce,(boost::python::arg("idSph")),"Return the fluid force on sphere idSph.")
+ .def("shearLubForce",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::shearLubForce,(boost::python::arg("idSph")),"Return the shear lubrication force on sphere idSph.")
+ .def("shearLubTorque",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::shearLubTorque,(boost::python::arg("idSph")),"Return the shear lubrication torque on sphere idSph.")
+ .def("normalLubForce",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::normalLubForce,(boost::python::arg("idSph")),"Return the normal lubrication force on sphere idSph.")
+ .def("bodyShearLubStress",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::bodyShearLubStress,(boost::python::arg("idSph")),"Return the shear lubrication stress on sphere idSph.")
+ .def("bodyNormalLubStress",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::bodyNormalLubStress,(boost::python::arg("idSph")),"Return the normal lubrication stress on sphere idSph.")
+ .def("shearVelocity",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::shearVelocity,(boost::python::arg("idSph")),"Return the shear velocity of the interaction.")
+ .def("normalVelocity",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::normalVelocity,(boost::python::arg("idSph")),"Return the normal velocity of the interaction.")
+ .def("normalVect",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::normalVect,(boost::python::arg("idSph")),"Return the normal vector between particles.")
+ .def("surfaceDistanceParticle",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::surfaceDistanceParticle,(boost::python::arg("interaction")),"Return the distance between particles.")
+ .def("onlySpheresInteractions",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::onlySpheresInteractions,(boost::python::arg("interaction")),"Return the id of the interaction only between spheres.")
+ .def("pressureProfile",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::pressureProfile,(boost::python::arg("wallUpY"),boost::python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
+ .def("getPorePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getPorePressure,(boost::python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
+ .def("averageSlicePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averageSlicePressure,(boost::python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
+ .def("averagePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averagePressure,"Measure averaged pore pressure in the entire volume")
+ .def("updateBCs",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::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_@TEMPLATE_FLOW_NAME@::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
+ .def("getCell",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getCell,(boost::python::arg("pos")),"get id of the cell containing (X,Y,Z).")
+ .def("nCells",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::nCells,"get the total number of finite cells in the triangulation.")
+ .def("getVertices",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getVertices,(boost::python::arg("id")),"get the vertices of a cell")
+ #ifdef LINSOLV
+ .def("exportMatrix",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::exportMatrix,(boost::python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
+ .def("exportTriplets",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::exportTriplets,(boost::python::arg("filename")="triplets"),"Export system matrix to a file with only non-zero entries.")
+ .def("cholmodStats",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::cholmodStats,"get statistics of cholmod solver activity")
+ .def("metisUsed",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::metisUsed,"check wether metis lib is effectively used")
+ .add_property("forceMetis",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getForceMetis,&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::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_@TEMPLATE_FLOW_NAME@::compTessVolumes,"Like TesselationWrapper::computeVolumes()")
+ .def("volume",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getVolume,(boost::python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
+ .def("averageVelocity",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averageVelocity,"measure the mean velocity in the period")
+ )
+};
+// Definition of functions in a separate file for clarity
+
+class FlowCellInfo_@TEMPLATE_FLOW_NAME@ : 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.doInterpolate
+ double solidSurfaces [4][4];
+
+ FlowCellInfo_@TEMPLATE_FLOW_NAME@ (void)
+ {
+ modulePermeability.resize(4, 0);
+ cellForce.resize(4,CGAL::NULL_VECTOR);
+ facetSurfaces.resize(4,CGAL::NULL_VECTOR);
+ facetFluidSurfacesRatio.resize(4,0);
+ facetSphereCrossSections.resize(4,CGAL::NULL_VECTOR);
+ unitForceVectors.resize(4,CGAL::NULL_VECTOR);
+ 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_@TEMPLATE_FLOW_NAME@& otherCellInfo) {p()=otherCellInfo.shiftedP();}
+};
+
+class FlowVertexInfo_@TEMPLATE_FLOW_NAME@ : public CGT::SimpleVertexInfo {
+ CVector grainVelocity;
+ Real volumeIncidentCells;
+public:
+ CVector forces;
+ bool isGhost;
+ FlowVertexInfo_@TEMPLATE_FLOW_NAME@ (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) const {return CGAL::NULL_VECTOR;}
+};
+
+#include "FlowEngine_@TEMPLATE_FLOW_NAME@.ipp"
=== modified file 'pkg/pfv/FlowEngine.ipp'
--- pkg/pfv/FlowEngine.ipp 2014-06-04 17:19:50 +0000
+++ pkg/pfv/FlowEngine.ipp 2014-06-20 15:28:21 +0000
@@ -1,750 +1,1 @@
-/*************************************************************************
-* 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" );
- #ifdef LINSOLV
- 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 || doInterpolate) 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
- #endif
- {
- 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.tesselation().Clear();
- flow.tesselation().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 )
-{
- if (first) flow.currentTes=0;
- else { flow.currentTes=!flow.currentTes; if (debug) cout << "--------RETRIANGULATION-----------" << endl;}
- flow.resetNetwork();
- initSolver(flow);
-
- addBoundary ( flow );
- triangulate ( flow );
- if ( debug ) cout << endl << "Tesselating------" << endl << endl;
- flow.tesselation().compute();
-
- flow.defineFictiousCells();
- // For faster loops on cells define this vector
- flow.tesselation().cellHandles.clear();
- flow.tesselation().cellHandles.reserve(flow.tesselation().Triangulation().number_of_finite_cells());
- FiniteCellsIterator cell_end = flow.tesselation().Triangulation().finite_cells_end();
- int k=0;
- for ( FiniteCellsIterator cell = flow.tesselation().Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
- flow.tesselation().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 || doInterpolate)) flow.interpolate ( flow.T[!flow.currentTes], flow.tesselation() );
- if ( waveAction ) flow.applySinusoidalPressure ( flow.tesselation().Triangulation(), sineMagnitude, sineAverage, 30 );
- else if (boundaryPressure.size()!=0) flow.applyUserDefinedPressure ( flow.tesselation().Triangulation(), boundaryXPos , boundaryPressure);
- 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.tesselation().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.tesselation());//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.tesselation().insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
- }
- flow.tesselation().redirected=true;//By inserting one-by-one, we already redirected
- flow.shearLubricationForces.resize ( flow.tesselation().maxId+1 );
- flow.shearLubricationTorques.resize ( flow.tesselation().maxId+1 );
- flow.pumpLubricationTorques.resize ( flow.tesselation().maxId+1 );
- flow.twistLubricationTorques.resize ( flow.tesselation().maxId+1 );
- flow.shearLubricationBodyStress.resize ( flow.tesselation().maxId+1 );
- flow.normalLubricationForce.resize ( flow.tesselation().maxId+1 );
- flow.normalLubricationBodyStress.resize ( flow.tesselation().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.tesselation().Triangulation().finite_vertices_end();
- CGT::CVector Zero(0,0,0);
- for (FiniteVerticesIterator V_it = flow.tesselation().Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero;
-
- FOREACH(CellHandle& cell, flow.tesselation().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.tesselation().cellHandles.size();
- #pragma omp parallel for num_threads(ompThreads>0 ? ompThreads : 1)
- for(long i=0; i<size; i++){
- CellHandle& cell = flow.tesselation().cellHandles[i];
- #else
- FOREACH(CellHandle& cell, flow.tesselation().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+=cell->info().volumeSign*newVol;
- #pragma omp atomic
- totDVol+=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.tesselation();
- 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 VertexInfo& vi1 = *flow.edgeIds[i].first;
- const VertexInfo& vi2 = *flow.edgeIds[i].second;
- const int& id1 = vi1.id();
- const int& id2 = vi2.id();
-
- 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(vi2.ghostShift()) - sph1->state->pos - makeVector3r(vi1.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);
-
- }
- }
-}
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-Vector3r TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::averageVelocity()
-{
- 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 );
-}
-
-#endif //FLOW_ENGINE
-
-#endif /* YADE_CGAL */
=== added file 'pkg/pfv/FlowEngine.ipp.in'
--- pkg/pfv/FlowEngine.ipp.in 1970-01-01 00:00:00 +0000
+++ pkg/pfv/FlowEngine.ipp.in 2014-06-20 15:28:21 +0000
@@ -0,0 +1,737 @@
+/*************************************************************************
+* 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
+
+#ifdef LINSOLV
+#include <cholmod.h>
+#endif
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::~TemplateFlowEngine_@TEMPLATE_FLOW_NAME@() {}
+
+// YADE_PLUGIN((TFlowEng));
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+unsigned int TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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" );
+ #ifdef LINSOLV
+ 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 || doInterpolate) 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_@TEMPLATE_FLOW_NAME@::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
+ #endif
+ {
+ 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_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedPressure () { solver->imposedP.clear(); solver->IPCells.clear();}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedFlux () { solver->imposedF.clear(); solver->IFCells.clear();}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+Real TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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.tesselation().Clear();
+ flow.tesselation().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_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::getForceMetis () {return (solver->eSolver.cholmod().nmethods==1);}
+#endif
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( Solver& flow )
+{
+ buildTriangulation ( 0.f,flow );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( double pZero, Solver& flow )
+{
+ if (first) flow.currentTes=0;
+ else { flow.currentTes=!flow.currentTes; if (debug) cout << "--------RETRIANGULATION-----------" << endl;}
+ flow.resetNetwork();
+ initSolver(flow);
+
+ addBoundary ( flow );
+ triangulate ( flow );
+ if ( debug ) cout << endl << "Tesselating------" << endl << endl;
+ flow.tesselation().compute();
+
+ flow.defineFictiousCells();
+ // For faster loops on cells define this vector
+ flow.tesselation().cellHandles.clear();
+ flow.tesselation().cellHandles.reserve(flow.tesselation().Triangulation().number_of_finite_cells());
+ FiniteCellsIterator cell_end = flow.tesselation().Triangulation().finite_cells_end();
+ int k=0;
+ for ( FiniteCellsIterator cell = flow.tesselation().Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
+ flow.tesselation().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 || doInterpolate)) flow.interpolate ( flow.T[!flow.currentTes], flow.tesselation() );
+ if ( waveAction ) flow.applySinusoidalPressure ( flow.tesselation().Triangulation(), sineMagnitude, sineAverage, 30 );
+ else if (boundaryPressure.size()!=0) flow.applyUserDefinedPressure ( flow.tesselation().Triangulation(), boundaryXPos , boundaryPressure);
+ if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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.tesselation().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_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::triangulate ( Solver& flow )
+{
+///Using Tesselation wrapper (faster)
+// TesselationWrapper TW;
+// if (TW.Tes) delete TW.Tes;
+// TW.Tes = &(flow.tesselation());//point to the current Tes we have in Flowengine
+// TW.insertSceneSpheres();//TW is now really inserting in TemplateFlowEngine_@TEMPLATE_FLOW_NAME@, 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.tesselation().insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
+ }
+ flow.tesselation().redirected=true;//By inserting one-by-one, we already redirected
+ flow.shearLubricationForces.resize ( flow.tesselation().maxId+1 );
+ flow.shearLubricationTorques.resize ( flow.tesselation().maxId+1 );
+ flow.pumpLubricationTorques.resize ( flow.tesselation().maxId+1 );
+ flow.twistLubricationTorques.resize ( flow.tesselation().maxId+1 );
+ flow.shearLubricationBodyStress.resize ( flow.tesselation().maxId+1 );
+ flow.normalLubricationForce.resize ( flow.tesselation().maxId+1 );
+ flow.normalLubricationBodyStress.resize ( flow.tesselation().maxId+1 );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::initializeVolumes ( Solver& flow )
+{
+ typedef typename Solver::FiniteVerticesIterator FiniteVerticesIterator;
+
+ FiniteVerticesIterator vertices_end = flow.tesselation().Triangulation().finite_vertices_end();
+ CGT::CVector Zero(0,0,0);
+ for (FiniteVerticesIterator V_it = flow.tesselation().Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero;
+
+ FOREACH(CellHandle& cell, flow.tesselation().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_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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.tesselation().cellHandles.size();
+ #pragma omp parallel for num_threads(ompThreads>0 ? ompThreads : 1)
+ for(long i=0; i<size; i++){
+ CellHandle& cell = flow.tesselation().cellHandles[i];
+ #else
+ FOREACH(CellHandle& cell, flow.tesselation().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+=cell->info().volumeSign*newVol;
+ #pragma omp atomic
+ totDVol+=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_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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_@TEMPLATE_FLOW_NAME@<_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.tesselation();
+ 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 VertexInfo& vi1 = *flow.edgeIds[i].first;
+ const VertexInfo& vi2 = *flow.edgeIds[i].second;
+ const int& id1 = vi1.id();
+ const int& id2 = vi2.id();
+
+ 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(vi2.ghostShift()) - sph1->state->pos - makeVector3r(vi1.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);
+
+ }
+ }
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+Vector3r TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::averageVelocity()
+{
+ 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 );
+}
+
+#endif //FLOW_ENGINE
+
+#endif /* YADE_CGAL */
+
=== modified file 'pkg/pfv/PeriodicFlowEngine.cpp'
--- pkg/pfv/PeriodicFlowEngine.cpp 2014-06-04 17:19:50 +0000
+++ pkg/pfv/PeriodicFlowEngine.cpp 2014-06-20 15:28:21 +0000
@@ -13,10 +13,9 @@
/// 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
-#define TEMPLATE_FLOW_NAME FlowEngine_PeriodicInfo
-#include <yade/pkg/pfv/FlowEngine.hpp>
+#include "FlowEngine_FlowEngine_PeriodicInfo.hpp"
-class PeriodicCellInfo : public FlowCellInfo
+class PeriodicCellInfo : public FlowCellInfo_FlowEngine_PeriodicInfo
{
public:
static CVector gradP;
@@ -43,7 +42,7 @@
};
-class PeriodicVertexInfo : public FlowVertexInfo {
+class PeriodicVertexInfo : public FlowVertexInfo_DummyFlowEngineT {
public:
PeriodicVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
PeriodicVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
@@ -64,7 +63,7 @@
#define _PeriFlowSolver CGT::PeriodicFlow<PeriFlowTesselation>
#endif
-typedef TemplateFlowEngine< PeriodicCellInfo,
+typedef TemplateFlowEngine_FlowEngine_PeriodicInfo< PeriodicCellInfo,
PeriodicVertexInfo,
CGT::PeriodicTesselation<CGT::_Tesselation<CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> > >,
_PeriFlowSolver>
@@ -529,4 +528,4 @@
YADE_PLUGIN((PeriodicFlowEngine));
#endif //FLOW_ENGINE
-#endif /* YADE_CGAL */
\ No newline at end of file
+#endif /* YADE_CGAL */
=== modified file 'pkg/pfv/SoluteFlowEngine.cpp'
--- pkg/pfv/SoluteFlowEngine.cpp 2014-05-23 13:20:43 +0000
+++ pkg/pfv/SoluteFlowEngine.cpp 2014-06-20 15:28:21 +0000
@@ -11,22 +11,21 @@
// #define SOLUTE_FLOW
#ifdef SOLUTE_FLOW
-#define TEMPLATE_FLOW_NAME SoluteFlowEngineT
-#include <yade/pkg/pfv/FlowEngine.hpp>
+#include "FlowEngine_SoluteFlowEngineT.hpp"
#include <Eigen/Sparse>
-class SoluteCellInfo : public FlowCellInfo
+class SoluteCellInfo : public FlowCellInfo_SoluteFlowEngineT
{
public:
Real solute_concentration;
- SoluteCellInfo (void) : FlowCellInfo() {solute_concentration=0;}
+ SoluteCellInfo (void) : FlowCellInfo_SoluteFlowEngineT() {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();}
+ inline void getInfo (const SoluteCellInfo& otherCellInfo) {FlowCellInfo_SoluteFlowEngineT::getInfo(otherCellInfo); solute()=otherCellInfo.solute();}
};
-typedef TemplateFlowEngine<SoluteCellInfo,FlowVertexInfo> SoluteFlowEngineT;
+typedef TemplateFlowEngine_SoluteFlowEngineT<SoluteCellInfo,FlowVertexInfo_DummyFlowEngineT> SoluteFlowEngineT;
REGISTER_SERIALIZABLE(SoluteFlowEngineT);
YADE_PLUGIN((SoluteFlowEngineT));
