yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10670
[Branch ~yade-pkg/yade/git-trunk] Rev 3889: Bugfixes in PeriodicFlowEngine, as discussed with Donia
------------------------------------------------------------
revno: 3889
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Mon 2014-04-07 11:33:19 +0200
message:
Bugfixes in PeriodicFlowEngine, as discussed with Donia
modified:
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphere.ipp
pkg/pfv/FlowEngine.hpp
pkg/pfv/FlowEngine.ipp
pkg/pfv/PeriodicFlowEngine.cpp
pkg/pfv/PeriodicFlowEngine.hpp
--
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 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2014-03-21 18:47:45 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2014-04-07 09:33:19 +0000
@@ -71,7 +71,7 @@
vector< vector<const Real*> > perVertexPressure;
#endif
vector <double> edgeSurfaces;
- vector <pair<int,int> > edgeIds;
+ vector <pair<const VertexInfo*,const VertexInfo*> > edgeIds;
vector <Real> edgeNormalLubF;
vector <Vector3r> shearLubricationForces;
vector <Vector3r> shearLubricationTorques;
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2014-04-03 13:07:09 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2014-04-07 09:33:19 +0000
@@ -89,7 +89,7 @@
}
template <class Tesselation>
-void FlowBoundingSphere<Tesselation>::resetNetwork() {T[0].Clear();noCache=true;}
+void FlowBoundingSphere<Tesselation>::resetNetwork() {T[currentTes].Clear();noCache=true;}
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::averageRelativeCellVelocity()
@@ -1201,41 +1201,44 @@
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::computeEdgesSurfaces()
{
- RTriangulation& Tri = T[currentTes].Triangulation();
-
- //first, copy interacting pairs and normal lub forces form prev. triangulation in a sorted structure for initializing the new lub. Forces
- vector<vector<pair<unsigned int,Real> > > lubPairs;
- lubPairs.resize(Tri.number_of_vertices()+1);
- for (unsigned int k=0; k<edgeNormalLubF.size(); k++)
- lubPairs[min(edgeIds[k].first,edgeIds[k].second)].push_back(pair<int,Real> (max(edgeIds[k].first,edgeIds[k].second),edgeNormalLubF[k]));
-
- //Now we reset the containers and initialize them
- edgeSurfaces.clear(); edgeIds.clear(); edgeNormalLubF.clear();
- FiniteEdgesIterator ed_it;
- for ( FiniteEdgesIterator ed_it = Tri.finite_edges_begin(); ed_it!=Tri.finite_edges_end();ed_it++ )
- {
- int hasFictious= (ed_it->first)->vertex(ed_it->second)->info().isFictious + (ed_it->first)->vertex(ed_it->third)->info().isFictious;
- if (hasFictious==2) continue;
- double area = T[currentTes].computeVFacetArea(ed_it);
- edgeSurfaces.push_back(area);
- unsigned int id1 = (ed_it->first)->vertex(ed_it->second)->info().id();
- unsigned int id2 = (ed_it->first)->vertex(ed_it->third)->info().id();
- edgeIds.push_back(pair<int,int>(id1,id2));
-
- //For persistant edges, we must transfer the lub. force value from the older triangulation structure
- if (id1>id2) swap(id1,id2);
- unsigned int i=0;
- //Look for the pair (id1,id2) in lubPairs
- while (i<lubPairs[id1].size()) {
- if (lubPairs[id1][i].first == id2) {
- //it's found, we copy the lub force
- edgeNormalLubF.push_back(lubPairs[id1][i].second);
- break;}
- ++i;
- }
- // not found, we initialize with zero lub force
- if (i==lubPairs[id1].size()) edgeNormalLubF.push_back(0);
- }
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ //first, copy interacting pairs and normal lub forces form prev. triangulation in a sorted structure for initializing the new lub. Forces
+ vector<vector<pair<unsigned int,Real> > > lubPairs;
+ lubPairs.resize(Tri.number_of_vertices()+1);
+ for (unsigned int k=0; k<edgeNormalLubF.size(); k++)
+ lubPairs[min(edgeIds[k].first->id(),edgeIds[k].second->id())].push_back(
+ pair<int,Real> (max(edgeIds[k].first->id(),edgeIds[k].second->id()),edgeNormalLubF[k]));
+
+ //Now we reset the containers and initialize them
+ edgeSurfaces.clear(); edgeIds.clear(); edgeNormalLubF.clear();
+ FiniteEdgesIterator ed_it;
+ for ( FiniteEdgesIterator ed_it = Tri.finite_edges_begin(); ed_it!=Tri.finite_edges_end();ed_it++ )
+ {
+ const VertexInfo& vi1=(ed_it->first)->vertex(ed_it->second)->info();
+ const VertexInfo& vi2=(ed_it->first)->vertex(ed_it->third)->info();
+
+ //We eliminate edges that would be periodic replications or involving two bounding objects, i.e. the min id must be non-ghost and non-fictious
+ if (vi1.id()<vi2.id()) {if (vi1.isFictious || vi2.isGhost) continue;}
+ else if (vi2.isFictious || vi2.isGhost) continue;
+ double area = T[currentTes].computeVFacetArea(ed_it);
+ edgeSurfaces.push_back(area);
+ unsigned int id1 = vi1.id();
+ unsigned int id2 = vi2.id();
+ edgeIds.push_back(pair<const VertexInfo*,const VertexInfo*>(&vi1,&vi2));
+ //For persistant edges, we must transfer the lub. force value from the older triangulation structure
+ if (id1>id2) swap(id1,id2);
+ unsigned int i=0;
+ //Look for the pair (id1,id2) in lubPairs
+ while (i<lubPairs[id1].size()) {
+ if (lubPairs[id1][i].first == id2) {
+ //it's found, we copy the lub force
+ edgeNormalLubF.push_back(lubPairs[id1][i].second);
+ break;}
+ ++i;
+ }
+ // not found, we initialize with zero lub force
+ if (i==lubPairs[id1].size()) edgeNormalLubF.push_back(0);
+ }
}
template <class Tesselation>
=== modified file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp 2014-04-03 16:44:33 +0000
+++ pkg/pfv/FlowEngine.hpp 2014-04-07 09:33:19 +0000
@@ -56,16 +56,18 @@
#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=CGT::FlowBoundingSphere<_Tesselation> >
+template<class _CellInfo, class _VertexInfo, class _Tesselation=CGT::_Tesselation<CGT::TriangulationTypes<_VertexInfo,_CellInfo> >, class solverT=DEFAULTSOLVER >
class TemplateFlowEngine : public PartialEngine
{
public :
- #ifdef LINSOLV
- typedef CGT::FlowBoundingSphereLinSolv<typename solverT::Tesselation,solverT> FlowSolver;
- #else
typedef solverT FlowSolver;
- #endif
typedef FlowSolver Solver;//FIXME: useless alias, search/replace then remove
typedef typename FlowSolver::Tesselation Tesselation;
typedef typename FlowSolver::RTriangulation RTriangulation;
@@ -227,8 +229,8 @@
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.");}
+ 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.",
@@ -426,7 +428,7 @@
inline CVector& force (void) {return forces;}
inline CVector& vel (void) {return grainVelocity;}
inline Real& volCells (void) {return volumeIncidentCells;}
- inline const CVector ghostShift (void) {return CGAL::NULL_VECTOR;}
+ inline const CVector ghostShift (void) const {return CGAL::NULL_VECTOR;}
};
=== modified file 'pkg/pfv/FlowEngine.ipp'
--- pkg/pfv/FlowEngine.ipp 2014-04-03 13:37:55 +0000
+++ pkg/pfv/FlowEngine.ipp 2014-04-07 09:33:19 +0000
@@ -219,7 +219,7 @@
flow.meanKStat = meanKStat;
flow.permeabilityMap = permeabilityMap;
flow.fluidBulkModulus = fluidBulkModulus;
- flow.T[flow.currentTes].Clear();
+// flow.T[flow.currentTes].Clear();
flow.T[flow.currentTes].maxId=-1;
flow.xMin = 1000.0, flow.xMax = -10000.0, flow.yMin = 1000.0, flow.yMax = -10000.0, flow.zMin = 1000.0, flow.zMax = -10000.0;
}
@@ -245,13 +245,12 @@
template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( double pZero, Solver& flow )
{
- flow.resetNetwork();
- if (first) flow.currentTes=0;
+ if (first) flow.currentTes=0;
else {
flow.currentTes=!flow.currentTes;
if (debug) cout << "--------RETRIANGULATION-----------" << endl;
}
-
+ flow.resetNetwork();
initSolver(flow);
addBoundary ( flow );
@@ -597,8 +596,10 @@
for ( int i=0; i< ( int ) flow.edgeIds.size(); i++ ) {
- const int& id1 = flow.edgeIds[i].first;
- const int& id2 = flow.edgeIds[i].second;
+ 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;
@@ -613,7 +614,7 @@
Vector3r shearLubF; Vector3r normaLubF; Vector3r pumpT; Vector3r deltaAngNormVel; Vector3r twistT; Vector3r angVel1; Vector3r angVel2;
//FIXME: if periodic and velGrad!=0, then deltaV should account for velGrad, not the case currently
if ( !hasFictious ){
- O1O2Vector = sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()) - sph1->state->pos - makeVector3r(Tes.vertex(id1)->info().ghostShift());
+ O1O2Vector = sph2->state->pos + makeVector3r(vi2.ghostShift()) - sph1->state->pos - makeVector3r(vi1.ghostShift());
O1O2 = O1O2Vector.norm();
normal= (O1O2Vector/O1O2);
surfaceDist = O1O2 - r2 - r1;
=== modified file 'pkg/pfv/PeriodicFlowEngine.cpp'
--- pkg/pfv/PeriodicFlowEngine.cpp 2014-04-03 13:37:55 +0000
+++ pkg/pfv/PeriodicFlowEngine.cpp 2014-04-07 09:33:19 +0000
@@ -9,9 +9,109 @@
#ifdef YADE_CGAL
#ifdef FLOW_ENGINE
+/// The periodic variant of FlowEngine is defined here. It should become a template class for more flexibility.
+/// It is a bit more complicated as for FlowEngine, though, because we need template inheriting from template, which breaks YADE_CLASS_XXX logic_error
+/// See below the commented exemple, for a possible solution
+
#define TEMPLATE_FLOW_NAME FlowEngine_PeriodicInfo
-#include "PeriodicFlowEngine.hpp"
-#undef TEMPLATE_FLOW_NAME
+#include <yade/pkg/pfv/FlowEngine.hpp>
+
+class PeriodicCellInfo : public FlowCellInfo
+{
+ public:
+ static CVector gradP;
+ //for real cell, baseIndex is the rank of the cell in cellHandles. For ghost cells, it is the baseIndex of the corresponding real cell.
+ //Unlike ordinary index, baseIndex is also indexing cells with imposed pressures
+ int baseIndex;
+ int period[3];
+ static CVector hSize[3];
+ static CVector deltaP;
+ int ghost;
+ Real* _pression;
+ PeriodicCellInfo (void){
+ _pression=&pression;
+ period[0]=period[1]=period[2]=0;
+ baseIndex=-1;
+ volumeSign=0;}
+ ~PeriodicCellInfo (void) {}
+
+ inline const double shiftedP (void) const {return isGhost? (*_pression)+pShift() :(*_pression) ;}
+ inline const double pShift (void) const {return deltaP[0]*period[0] + deltaP[1]*period[1] +deltaP[2]*period[2];}
+// inline const double p (void) {return shiftedP();}
+ inline void setP (const Real& p) {pression=p;}
+ bool isReal (void) {return !(isFictious || isGhost);}
+};
+
+
+class PeriodicVertexInfo : public FlowVertexInfo {
+ public:
+ PeriodicVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
+ PeriodicVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
+ PeriodicVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
+ int period[3];
+ //FIXME: the name is misleading, even non-ghost can be out of the period and therefore they need to be shifted as well
+ inline const CVector ghostShift (void) const {
+ return period[0]*PeriodicCellInfo::hSize[0]+period[1]*PeriodicCellInfo::hSize[1]+period[2]*PeriodicCellInfo::hSize[2];}
+ PeriodicVertexInfo (void) {isFictious=false; s=0; i=0; period[0]=period[1]=period[2]=0; isGhost=false;}
+ bool isReal (void) {return !(isFictious || isGhost);}
+};
+
+typedef CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> PeriFlowTriangulationTypes;
+typedef CGT::PeriodicTesselation<CGT::_Tesselation<PeriFlowTriangulationTypes> > PeriFlowTesselation;
+#ifdef LINSOLV
+#define _PeriFlowSolver CGT::PeriodicFlowLinSolv<PeriFlowTesselation>
+#else
+#define _PeriFlowSolver CGT::PeriodicFlow<PeriFlowTesselation>
+#endif
+
+typedef TemplateFlowEngine< PeriodicCellInfo,
+ PeriodicVertexInfo,
+ CGT::PeriodicTesselation<CGT::_Tesselation<CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> > >,
+ _PeriFlowSolver>
+ FlowEngine_PeriodicInfo;
+
+REGISTER_SERIALIZABLE(FlowEngine_PeriodicInfo);
+YADE_PLUGIN((FlowEngine_PeriodicInfo));
+
+
+class PeriodicFlowEngine : public FlowEngine_PeriodicInfo
+{
+ public :
+ void triangulate (FlowSolver& flow);
+ void buildTriangulation (Real pzero, FlowSolver& flow);
+ void initializeVolumes (FlowSolver& flow);
+ void updateVolumes (FlowSolver& flow);
+ Real volumeCell (CellHandle cell);
+
+ Real volumeCellSingleFictious (CellHandle cell);
+ inline void locateCell(CellHandle baseCell, unsigned int& index, int& baseIndex, FlowSolver& flow, unsigned int count=0);
+ Vector3r meanVelocity();
+
+ virtual ~PeriodicFlowEngine();
+
+ virtual void action();
+ //Cache precomputed values for pressure shifts, based on current hSize and pGrad
+ void preparePShifts();
+
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(PeriodicFlowEngine,FlowEngine_PeriodicInfo,"A variant of :yref:`FlowEngine` implementing periodic boundary conditions. The API is very similar.",
+ ((Real,duplicateThreshold, 0.06,,"distance from cell borders that will triger periodic duplication in the triangulation |yupdate|"))
+ ((Vector3r, gradP, Vector3r::Zero(),,"Macroscopic pressure gradient"))
+ ,,
+ wallIds=vector<int>(6,-1);
+ solver = shared_ptr<FlowSolver> (new FlowSolver);
+ epsVolMax=epsVolCumulative=retriangulationLastIter=0;
+ ReTrg=1;
+ first=true;
+ ,
+ //nothing special to define, we re-use FlowEngine methods
+ //.def("meanVelocity",&PeriodicFlowEngine::meanVelocity,"measure the mean velocity in the period")
+ )
+ DECLARE_LOGGER;
+
+
+};
+REGISTER_SERIALIZABLE(PeriodicFlowEngine);
+
CVector PeriodicCellInfo::hSize[]={CVector(),CVector(),CVector()};
CVector PeriodicCellInfo::deltaP=CVector();
@@ -382,11 +482,11 @@
void PeriodicFlowEngine::buildTriangulation ( double pZero, FlowSolver& flow)
{
- flow.resetNetwork();
if (first) flow.currentTes=0;
else {
flow.currentTes=!flow.currentTes;
if ( debug ) cout << "--------RETRIANGULATION-----------" << endl;}
+ flow.resetNetwork();
initSolver(flow);
addBoundary ( flow );
if ( debug ) cout << endl << "Added boundaries------" << endl << endl;
@@ -445,8 +545,6 @@
YADE_PLUGIN((PeriodicFlowEngine));
-
-
#endif //FLOW_ENGINE
#endif /* YADE_CGAL */
\ No newline at end of file
=== modified file 'pkg/pfv/PeriodicFlowEngine.hpp'
--- pkg/pfv/PeriodicFlowEngine.hpp 2014-04-03 16:44:33 +0000
+++ pkg/pfv/PeriodicFlowEngine.hpp 2014-04-07 09:33:19 +0000
@@ -6,104 +6,3 @@
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
-/// The periodic variant of FlowEngine is defined here. It should become a template class for more flexibility.
-/// It is a bit more complicated as for FlowEngine, though, because we need template inheriting from template, which breaks YADE_CLASS_XXX logic_error
-/// See below the commented exemple, for a possible solution
-
-#include <yade/pkg/pfv/FlowEngine.hpp>
-
-class PeriodicCellInfo : public FlowCellInfo
-{
- public:
- static CVector gradP;
- //for real cell, baseIndex is the rank of the cell in cellHandles. For ghost cells, it is the baseIndex of the corresponding real cell.
- //Unlike ordinary index, baseIndex is also indexing cells with imposed pressures
- int baseIndex;
- int period[3];
- static CVector hSize[3];
- static CVector deltaP;
- int ghost;
- Real* _pression;
- PeriodicCellInfo (void){
- _pression=&pression;
- period[0]=period[1]=period[2]=0;
- baseIndex=-1;
- volumeSign=0;}
- ~PeriodicCellInfo (void) {}
-
- inline const double shiftedP (void) const {return isGhost? (*_pression)+pShift() :(*_pression) ;}
- inline const double pShift (void) const {return deltaP[0]*period[0] + deltaP[1]*period[1] +deltaP[2]*period[2];}
-// inline const double p (void) {return shiftedP();}
- inline void setP (const Real& p) {pression=p;}
- bool isReal (void) {return !(isFictious || isGhost);}
-};
-
-
-class PeriodicVertexInfo : public FlowVertexInfo {
- public:
- PeriodicVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
- PeriodicVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
- PeriodicVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
- int period[3];
- //FIXME: the name is misleading, even non-ghost can be out of the period and therefore they need to be shifted as well
- inline const CVector ghostShift (void) {
- return period[0]*PeriodicCellInfo::hSize[0]+period[1]*PeriodicCellInfo::hSize[1]+period[2]*PeriodicCellInfo::hSize[2];}
- PeriodicVertexInfo (void) {isFictious=false; s=0; i=0; period[0]=period[1]=period[2]=0; isGhost=false;}
- bool isReal (void) {return !(isFictious || isGhost);}
-};
-
-typedef CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> PeriFlowTriangulationTypes;
-typedef CGT::PeriodicTesselation<CGT::_Tesselation<PeriFlowTriangulationTypes> > PeriFlowTesselation;
-#ifdef LINSOLV
-#define _PeriFlowSolver CGT::PeriodicFlowLinSolv<PeriFlowTesselation>
-#else
-#define _PeriFlowSolver CGT::PeriodicFlow<PeriFlowTesselation>
-#endif
-
-typedef TemplateFlowEngine< PeriodicCellInfo,
- PeriodicVertexInfo,
- CGT::PeriodicTesselation<CGT::_Tesselation<CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> > >,
- _PeriFlowSolver>
- FlowEngine_PeriodicInfo;
-
-REGISTER_SERIALIZABLE(FlowEngine_PeriodicInfo);
-YADE_PLUGIN((FlowEngine_PeriodicInfo));
-
-
-class PeriodicFlowEngine : public FlowEngine_PeriodicInfo
-{
- public :
- void triangulate (FlowSolver& flow);
- void buildTriangulation (Real pzero, FlowSolver& flow);
- void initializeVolumes (FlowSolver& flow);
- void updateVolumes (FlowSolver& flow);
- Real volumeCell (CellHandle cell);
-
- Real volumeCellSingleFictious (CellHandle cell);
- inline void locateCell(CellHandle baseCell, unsigned int& index, int& baseIndex, FlowSolver& flow, unsigned int count=0);
- Vector3r meanVelocity();
-
- virtual ~PeriodicFlowEngine();
-
- virtual void action();
- //Cache precomputed values for pressure shifts, based on current hSize and pGrad
- void preparePShifts();
-
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(PeriodicFlowEngine,FlowEngine_PeriodicInfo,"A variant of :yref:`FlowEngine` implementing periodic boundary conditions. The API is very similar.",
- ((Real,duplicateThreshold, 0.06,,"distance from cell borders that will triger periodic duplication in the triangulation |yupdate|"))
- ((Vector3r, gradP, Vector3r::Zero(),,"Macroscopic pressure gradient"))
- ,,
- wallIds=vector<int>(6,-1);
- solver = shared_ptr<FlowSolver> (new FlowSolver);
- epsVolMax=epsVolCumulative=retriangulationLastIter=0;
- ReTrg=1;
- first=true;
- ,
- //nothing special to define, we re-use FlowEngine methods
- //.def("meanVelocity",&PeriodicFlowEngine::meanVelocity,"measure the mean velocity in the period")
- )
- DECLARE_LOGGER;
-
-
-};
-REGISTER_SERIALIZABLE(PeriodicFlowEngine);
\ No newline at end of file
