yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10595
[Branch ~yade-pkg/yade/git-trunk] Rev 3849: PFV: simplify (c++ wise) value assignement to Info types. More refactoring. Final step of code cl...
------------------------------------------------------------
revno: 3849
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Fri 2014-03-21 19:47:45 +0100
message:
PFV: simplify (c++ wise) value assignement to Info types. More refactoring. Final step of code cleaning.
modified:
lib/triangulation/RegularTriangulation.h
lib/triangulation/Tesselation.ipp
pkg/dem/DFNFlow.cpp
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.hpp
pkg/dem/TesselationWrapper.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 'lib/triangulation/RegularTriangulation.h'
--- lib/triangulation/RegularTriangulation.h 2014-03-21 18:47:45 +0000
+++ lib/triangulation/RegularTriangulation.h 2014-03-21 18:47:45 +0000
@@ -52,14 +52,14 @@
Real s;
bool isFictious;
SimpleCellInfo (void) {isFictious=false; s=0;}
- SimpleCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
- SimpleCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
+ SimpleCellInfo& setPoint(const Point &p) { Point::operator= (p); return *this; }
+ SimpleCellInfo& setScalar(const Real &scalar) { s=scalar; return *this; }
inline Real x (void) {return Point::x();}
inline Real y (void) {return Point::y();}
inline Real z (void) {return Point::z();}
inline Real& f (void) {return s;}
//virtual function that will be defined for all classes, allowing shared function (e.g. for display of periodic and non-periodic with the unique function saveVTK)
- virtual bool isReal (void) {return !isFictious;}
+ bool isReal (void) {return !isFictious;}
};
class SimpleVertexInfo : public CVector {
@@ -69,9 +69,9 @@
Real vol;
public:
bool isFictious;
- SimpleVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
- SimpleVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
- SimpleVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
+ SimpleVertexInfo& setVector(const CVector &u) { CVector::operator= (u); return *this; }
+ SimpleVertexInfo& setFloat(const float &scalar) { s=scalar; return *this; }
+ SimpleVertexInfo& setId(const unsigned int &id) { i= id; return *this; }
inline Real ux (void) {return CVector::x();}
inline Real uy (void) {return CVector::y();}
inline Real uz (void) {return CVector::z();}
@@ -80,145 +80,10 @@
inline const unsigned int& id (void) const {return i;}
SimpleVertexInfo (void) {isFictious=false; s=0; i=0; vol=-1;}
//virtual function that will be defined for all classes, allowing shared function (e.g. for display)
- virtual bool isReal (void) {return !isFictious;}
+ bool isReal (void) {return !isFictious;}
};
-// class FlowCellInfo : public SimpleCellInfo {
-//
-// public:
-// //For vector storage of all cells
-// unsigned int index;
-// int volumeSign;
-// bool Pcondition;
-// Real invVoidV;
-// Real t;
-// int fict;
-// Real volumeVariation;
-// double pression;
-// //average relative (fluid - facet translation) velocity defined for a single cell as 1/Volume * SUM_ON_FACETS(x_average_facet*average_facet_flow_rate)
-// CVector averageCellVelocity;
-// // Surface vectors of facets, pointing from outside toward inside the cell
-// std::vector<CVector> facetSurfaces;
-// //Ratio between fluid surface and facet surface
-// std::vector<Real> facetFluidSurfacesRatio;
-// // Reflects the geometrical property of the cell, so that the force by cell fluid on grain "i" is pressure*unitForceVectors[i]
-// std::vector<CVector> unitForceVectors;
-// // Store the area of triangle-sphere intersections for each facet (used in forces definition)
-// std::vector<CVector> facetSphereCrossSections;
-// std::vector<CVector> cellForce;
-// std::vector<double> rayHydr;
-// std::vector<double> modulePermeability;
-// // Partial surfaces of spheres in the double-tetrahedron linking two voronoi centers. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_surfaces in the facet.
-// double solidSurfaces [4][4];
-//
-// FlowCellInfo (void)
-// {
-// modulePermeability.resize(4, 0);
-// cellForce.resize(4);
-// facetSurfaces.resize(4);
-// facetFluidSurfacesRatio.resize(4);
-// facetSphereCrossSections.resize(4);
-// unitForceVectors.resize(4);
-// for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidSurfaces[k][l]=0;
-// rayHydr.resize(4, 0);
-// // isInside = false;
-// invSumK=0;
-// isFictious=false; Pcondition = false; isGhost = false;
-// // isInferior = false; isSuperior = false; isLateral = false; isExternal=false;
-// isvisited = false;
-// index=0;
-// volumeSign=0;
-// s=0;
-// volumeVariation=0;
-// pression=0;
-// invVoidV=0;
-// fict=0;
-// isGhost=false;
-// }
-// bool isGhost;
-// double invSumK;
-// bool isvisited;
-//
-// FlowCellInfo& operator= (const std::vector<double> &v) { for (int i=0; i<4;i++) modulePermeability[i]= v[i]; return *this; }
-// FlowCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
-// FlowCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
-//
-// 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;}
-// //For compatibility with the periodic case
-// inline const double shiftedP (void) const {return pression;}
-// 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;}
-// };
-//
-// class FlowVertexInfo : public SimpleVertexInfo {
-// CVector grainVelocity;
-// Real volumeIncidentCells;
-// public:
-// FlowVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
-// FlowVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
-// FlowVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
-// CVector forces;
-// bool isGhost;
-// FlowVertexInfo (void) {isGhost=false;}
-// inline CVector& force (void) {return forces;}
-// inline CVector& vel (void) {return grainVelocity;}
-// inline Real& volCells (void) {return volumeIncidentCells;}
-// inline const CVector ghostShift (void) {return CGAL::NULL_VECTOR;}
-// };
-
-// 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) {}
-// PeriodicCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
-// PeriodicCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
-//
-// 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;}
-// virtual 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;}
-// virtual bool isReal (void) {return !(isFictious || isGhost);}
-// };
-
-
-///
-
+
template<class vertex_info, class cell_info>
class TriangulationTypes {
=== modified file 'lib/triangulation/Tesselation.ipp'
--- lib/triangulation/Tesselation.ipp 2014-03-21 18:47:45 +0000
+++ lib/triangulation/Tesselation.ipp 2014-03-21 18:47:45 +0000
@@ -56,7 +56,7 @@
Vh = Tri->insert(Sphere(Point(x,y,z),pow(rad,2)));
if ( Vh!=NULL )
{
- Vh->info() = id;
+ Vh->info().setId(id);
Vh->info().isFictious = isFictious;
assert (vertexHandles.size()>id);
vertexHandles[id] = Vh;
@@ -74,7 +74,7 @@
if ( Vh!=NULL )
{
vertexHandles[id] = Vh;
- Vh->info() = id;
+ Vh->info().setId(id);
Vh->info().isFictious = fictious;
}
else cerr << "Vh==NULL" << " id=" << id << " Point=" << Point ( x,y,z ) << " rad=" << rad << endl;
@@ -158,7 +158,7 @@
S2.point().x(), S2.point().y(), S2.point().z(), S2.weight(),
S3.point().x(), S3.point().y(), S3.point().z(), S3.weight(),
x, y, z );
- cell->info() =Point ( x,y,z );
+ cell->info().setPoint(Point(x,y,z));
}
computed = true;
}
=== modified file 'pkg/dem/DFNFlow.cpp'
--- pkg/dem/DFNFlow.cpp 2014-03-21 18:47:45 +0000
+++ pkg/dem/DFNFlow.cpp 2014-03-21 18:47:45 +0000
@@ -10,7 +10,7 @@
//keep this #ifdef for commited versions unless you really have stable version that should be compiled by default
//it will save compilation time for everyone else
//when you want it compiled, you can pass -DDFNFLOW to cmake, or just uncomment the following line
-// #define DFNFLOW
+#define DFNFLOW
#ifdef DFNFLOW
#include <yade/pkg/dem/FlowEngine.hpp>
@@ -35,7 +35,7 @@
class DFNFlowEngine : public DFNFlowEngineT
{
public :
- void pyTrickPermeability(Real somethingBig);
+ void trickPermeability();
void trickPermeability (RTriangulation::Facet_circulator& facet,Real somethingBig);
void trickPermeability (RTriangulation::Finite_edges_iterator& edge,Real somethingBig);
@@ -43,7 +43,6 @@
((Real, myNewAttribute, 0,,"useless example"))
,/*DFNFlowEngineT()*/,
,
- .def("trickPermeability",&DFNFlowEngine::pyTrickPermeability,(python::arg("somethingBig")=10.),"change permeability inside a plane")
)
DECLARE_LOGGER;
};
@@ -70,8 +69,9 @@
}
-void DFNFlowEngine::pyTrickPermeability(Real somethingBig)
+void DFNFlowEngine::trickPermeability()
{
+ Real somethingBig=10;//is that big??
const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
//We want to change permeability perpendicular to the 10th edge, let's say.
//in the end this function should have a loop on all edges I guess
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2014-03-21 18:47:45 +0000
+++ pkg/dem/FlowEngine.cpp 2014-03-21 18:47:45 +0000
@@ -1,7 +1,6 @@
/*************************************************************************
* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* Copyright (C) 2012 by Donia Marzougui <donia.marzougui@xxxxxxxxxxxxxxx>*
* *
* This program is free software; it is licensed under the terms of the *
* GNU General Public License v2 or later. See file LICENSE for details. *
@@ -9,483 +8,14 @@
#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"
-
-
-
-CREATE_LOGGER ( FlowEngine );
-CREATE_LOGGER ( PeriodicFlowEngine );
+YADE_PLUGIN((FlowEngineT));
+CREATE_LOGGER(FlowEngine );
YADE_PLUGIN((FlowEngine));
-//______________________________________________________________
-
-//___________________ PERIODIC VERSION _________________________
-//______________________________________________________________
-
-PeriodicFlowEngine::~PeriodicFlowEngine(){}
-
-void PeriodicFlowEngine:: action()
-{
- if ( !isActivated ) return;
- timingDeltas->start();
- preparePShifts();
- setPositionsBuffer(true);
- if (first) {
- if (multithread) setPositionsBuffer(false);
- cachedCell= Cell(*(scene->cell));
- buildTriangulation(pZero,*solver);
- if (solver->errorCode>0) {LOG_INFO("triangulation error, pausing"); Omega::instance().pause(); return;}
- initializeVolumes(*solver); backgroundSolver=solver; backgroundCompleted=true;}
-// if ( first ) {buildTriangulation ( pZero ); updateTriangulation = false; initializeVolumes();}
-
- timingDeltas->checkpoint("Triangulating");
- updateVolumes (*solver);
- epsVolCumulative += epsVolMax;
- retriangulationLastIter++;
- if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
- (defTolerance>0 && epsVolCumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval;
-
- timingDeltas->checkpoint("Update_Volumes");
-
- ///compute flow and and forces here
- if (pressureForce){
- solver->gaussSeidel(scene->dt);
- timingDeltas->checkpoint("Gauss-Seidel");
- solver->computeFacetForcesWithCache();}
- timingDeltas->checkpoint("compute_Pressure_Forces");
-
- ///compute vicscous forces
- scene->forces.sync();
- computeViscousForces(*solver);
- timingDeltas->checkpoint("compute_Viscous_Forces");
- Vector3r force;
- Vector3r torque;
- const Tesselation& Tes = solver->T[solver->currentTes];
- for (int id=0; id<=Tes.maxId; id++){
- assert (Tes.vertexHandles[id] != NULL);
- const Tesselation::VertexInfo& v_info = Tes.vertexHandles[id]->info();
- force =(pressureForce) ? Vector3r ( ( v_info.forces ) [0],v_info.forces[1],v_info.forces[2] ) : Vector3r(0,0,0);
- torque = Vector3r(0,0,0);
- if (shearLubrication || viscousShear){
- force = force +solver->shearLubricationForces[v_info.id()];
- torque = torque +solver->shearLubricationTorques[v_info.id()];
- if (pumpTorque)
- torque = torque +solver->pumpLubricationTorques[v_info.id()];
- if (twistTorque)
- torque = torque +solver->twistLubricationTorques[v_info.id()];
- }
-
- if (normalLubrication)
- force = force + solver->normalLubricationForce[v_info.id()];
- scene->forces.addForce ( v_info.id(), force);
- scene->forces.addTorque ( v_info.id(), torque);
- }
- ///End Compute flow and forces
- timingDeltas->checkpoint("Applying Forces");
- if (multithread && !first) {
- while (updateTriangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/ boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
- if (updateTriangulation || ellapsedIter>(0.5*meshUpdateInterval)) {
- if (useSolver==0) LOG_ERROR("background calculations not available for Gauss-Seidel");
- if (fluidBulkModulus>0) solver->interpolate (solver->T[solver->currentTes], backgroundSolver->T[backgroundSolver->currentTes]);
- solver=backgroundSolver;
- backgroundSolver = shared_ptr<FlowSolver> (new FlowSolver);
- //Copy imposed pressures/flow from the old solver
- backgroundSolver->imposedP = vector<pair<CGT::Point,Real> >(solver->imposedP);
- backgroundSolver->imposedF = vector<pair<CGT::Point,Real> >(solver->imposedF);
- setPositionsBuffer(false);
- cachedCell= Cell(*(scene->cell));
- backgroundCompleted=false;
- retriangulationLastIter=ellapsedIter;
- ellapsedIter=0;
- epsVolCumulative=0;
- boost::thread workerThread(&PeriodicFlowEngine::backgroundAction,this);
- workerThread.detach();
- initializeVolumes(*solver);
- computeViscousForces(*solver);
- }
- else if (debug && !first) {
- if (debug && !backgroundCompleted) cerr<<"still computing solver in the background"<<endl;
- ellapsedIter++;}
- } else {
- if (updateTriangulation && !first) {
- cachedCell= Cell(*(scene->cell));
- buildTriangulation (pZero, *solver);
- initializeVolumes(*solver);
- computeViscousForces(*solver);
- updateTriangulation = false;
- epsVolCumulative=0;
- retriangulationLastIter=0;
- ReTrg++;}
- }
- first=false;
- timingDeltas->checkpoint("Ending");
-}
-
-
-// void PeriodicFlowEngine::backgroundAction()
-// {
-// if (useSolver<1) {LOG_ERROR("background calculations not available for Gauss-Seidel"); return;}
-// buildTriangulation (pZero,*backgroundSolver);
-// //FIXME: GS is computing too much, we need only matrix factorization in fact
-// backgroundSolver->gaussSeidel(scene->dt);
-// backgroundSolver->computeFacetForcesWithCache(/*onlyCache?*/ true);
-// // boost::this_thread::sleep(boost::posix_time::seconds(10));
-// backgroundCompleted = true;
-// }
-
-void PeriodicFlowEngine::triangulate( FlowSolver& flow )
-{
- Tesselation& Tes = flow.T[flow.currentTes];
- vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
- FOREACH ( const posData& b, buffer ) {
- if ( !b.exists || !b.isSphere || b.id==ignoredBody) continue;
- Vector3i period; Vector3r wpos;
- // FIXME: use "sheared" variant if the cell is sheared
- wpos=cachedCell.wrapPt ( b.pos,period );
- const Body::id_t& id = b.id;
- const Real& rad = b.radius;
- const Real& x = wpos[0];
- const Real& y = wpos[1];
- const Real& z = wpos[2];
- VertexHandle vh0=Tes.insert ( x, y, z, rad, id );
-// VertexHandle vh0=Tes.insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );
- if (vh0==NULL) {
- flow.errorCode = 2;
- LOG_ERROR("Vh NULL in PeriodicFlowEngine::triangulate(), check input data"); continue;}
- for ( int k=0;k<3;k++ ) vh0->info().period[k]=-period[k];
- const Vector3r cellSize ( cachedCell.getSize() );
- //FIXME: if hasShear, comment in
-// wpos=scene->cell->unshearPt ( wpos );
- // traverse all periodic cells around the body, to see if any of them touches
- Vector3r halfSize= ( rad+duplicateThreshold ) *Vector3r ( 1,1,1 );
- Vector3r pmin,pmax;
- Vector3i i;
- for ( i[0]=-1; i[0]<=1; i[0]++ )
- for ( i[1]=-1;i[1]<=1; i[1]++ )
- for ( i[2]=-1; i[2]<=1; i[2]++ ) {
- if ( i[0]!=0 || i[1]!=0 || i[2]!=0 ) { // middle; already rendered above
- Vector3r pos2=wpos+Vector3r ( cellSize[0]*i[0],cellSize[1]*i[1],cellSize[2]*i[2] ); // shift, but without shear!
- pmin=pos2-halfSize;
- pmax=pos2+halfSize;
- if ( (pmin[0]<=cellSize[0]) && (pmax[0]>=0) && (pmin[1]<=cellSize[1]) && (pmax[1]>=0) && (pmin[2]<=cellSize[2]) && (pmax[2]>=0) ) {
- //with shear:
- //Vector3r pt=scene->cell->shearPt ( pos2 );
- //without shear:
- const Vector3r& pt= pos2;
- VertexHandle vh=Tes.insert ( pt[0],pt[1],pt[2],rad,id,false,id );
- for ( int k=0;k<3;k++ ) vh->info().period[k]=i[k]-period[k];}}
- }
- //re-assign the original vertex pointer since duplicates may have overwrite it
- Tes.vertexHandles[id]=vh0;
- }
- Tes.redirected=true;//By inserting one-by-one, we already redirected
- flow.shearLubricationForces.resize ( Tes.maxId+1 );
- flow.shearLubricationTorques.resize ( Tes.maxId+1 );
- flow.pumpLubricationTorques.resize ( Tes.maxId+1 );
- flow.twistLubricationTorques.resize ( Tes.maxId+1 );
- flow.shearLubricationBodyStress.resize ( Tes.maxId+1 );
- flow.normalLubricationForce.resize ( Tes.maxId+1 );
- flow.normalLubricationBodyStress.resize ( Tes.maxId+1 );
-}
-
-
-Real PeriodicFlowEngine::volumeCell ( CellHandle cell )
-{
- static const Real inv6 = 1/6.;
- const Vector3r p0 = positionBufferCurrent[cell->vertex(0)->info().id()].pos + makeVector3r(cell->vertex(0)->info().ghostShift());
- const Vector3r p1 = positionBufferCurrent[cell->vertex(1)->info().id()].pos + makeVector3r(cell->vertex(1)->info().ghostShift());
- const Vector3r p2 = positionBufferCurrent[cell->vertex(2)->info().id()].pos + makeVector3r(cell->vertex(2)->info().ghostShift());
- const Vector3r p3 = positionBufferCurrent[cell->vertex(3)->info().id()].pos + makeVector3r(cell->vertex(3)->info().ghostShift());
- Real volume = inv6*((p0-p1).cross(p0-p2)).dot(p0-p3);
- if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
- return volume;
-}
-
-Real PeriodicFlowEngine::volumeCellSingleFictious ( CellHandle cell )
-{
- Vector3r V[3];
- int b=0;
- int w=0;
- cell->info().volumeSign=1;
- Real Wall_coordinate=0;
-
- for ( int y=0;y<4;y++ ) {
- if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( y )->info().id(), scene );
- V[w]=sph->state->pos+ makeVector3r ( cell->vertex ( y )->info().ghostShift() );
- w++;
- } else {
- b = cell->vertex ( y )->info().id();
- const shared_ptr<Body>& wll = Body::byId ( b,scene );
- if ( !solver->boundary ( b ).useMaxMin ) Wall_coordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wallThickness/2.;
- else Wall_coordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate];
- }
- }
- Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
- return abs ( Volume );
-}
-
-
-void PeriodicFlowEngine::locateCell ( CellHandle baseCell, unsigned int& index, int& baseIndex, FlowSolver& flow, unsigned int count)
-{
- if (count>10) {
- LOG_ERROR("More than 10 attempts to locate a cell, duplicateThreshold may be too small, resulting in periodicity inconsistencies.");
- flow.errorCode=1; return;
- }
- PeriFlowTesselation::CellInfo& baseInfo = baseCell->info();
- //already located, return FIXME: is inline working correctly? else move this test outside the function, just before the calls
- if ( baseInfo.index>0 || baseInfo.isGhost ) return;
- RTriangulation& Tri = flow.T[flow.currentTes].Triangulation();
- Vector3r center ( 0,0,0 );
- Vector3i period;
-
- if (baseCell->info().fictious()==0)
- for ( int k=0;k<4;k++ ) center+= 0.25*makeVector3r (baseCell->vertex(k)->point());
- else {
-
- Real boundPos=0; int coord=0;
- for ( int k=0;k<4;k++ ) {
- if ( !baseCell->vertex ( k )->info().isFictious ) center+= 0.3333333333*makeVector3r ( baseCell->vertex ( k )->point() );
- else {
- coord=flow.boundary ( baseCell->vertex ( k )->info().id() ).coordinate;
- boundPos=flow.boundary ( baseCell->vertex ( k )->info().id() ).p[coord];}
- }
- center[coord]=boundPos;
- }
- Vector3r wdCenter= cachedCell.wrapPt ( center,period );
- if ( period[0]!=0 || period[1]!=0 || period[2]!=0 ) {
- if ( baseCell->info().index>0 ) {
- cout<<"indexed cell is found ghost!"<<baseInfo.index <<endl;
- baseInfo.isGhost=false;
- return;
- }
- CellHandle ch= Tri.locate ( CGT::Point ( wdCenter[0],wdCenter[1],wdCenter[2] )
-// ,/*hint*/ v0
- );
- baseInfo.period[0]=period[0];
- baseInfo.period[1]=period[1];
- baseInfo.period[2]=period[2];
- //call recursively, since the returned cell could be also a ghost (especially if baseCell is a non-periodic type from the external contour
- locateCell ( ch,index,baseIndex,flow,++count );
- if ( ch==baseCell ) cerr<<"WTF!!"<<endl;
- //check consistency
- bool checkC=false;
- for (int kk=0; kk<4;kk++) if ((!baseCell->vertex(kk)->info().isGhost) && ((!baseCell->vertex(kk)->info().isFictious))) checkC = true;
- if (checkC) {
- bool checkV=true;
- for (int kk=0; kk<4;kk++) {
- checkV=false;
- for (int jj=0; jj<4;jj++)
- if (baseCell->vertex(kk)->info().id() == ch->vertex(jj)->info().id()) checkV = true;
- if (!checkV) {cerr <<"periodicity is broken"<<endl;
- for (int jj=0; jj<4;jj++) cerr<<baseCell->vertex(jj)->info().id()<<" ";
- cerr<<" vs. ";
- for (int jj=0; jj<4;jj++) cerr<<ch->vertex(jj)->info().id()<<" ";
- cerr<<endl;}
- }
- } else {
-// bool checkV=true;
-// for (int kk=0; kk<4;kk++) {
-// checkV=false;
-// for (int jj=0; jj<4;jj++)
-// if (baseCell->vertex(kk)->info().id() == ch->vertex(jj)->info().id()) checkV = true;
-// if (!checkV) {cerr <<"periodicity is broken (that's ok probably)"<<endl;
-// for (int jj=0; jj<4;jj++) cerr<<baseCell->vertex(jj)->info().id()<<" ";
-// cerr<<" vs. ";
-// for (int jj=0; jj<4;jj++) cerr<<ch->vertex(jj)->info().id()<<" ";
-// cerr<<endl;}
-// }
- }
-
- baseInfo.isGhost=true;
- baseInfo._pression=& ( ch->info().p() );
- baseInfo.index=ch->info().index;
- baseInfo.baseIndex=ch->info().baseIndex;
- baseInfo.Pcondition=ch->info().Pcondition;
- } else {
- baseInfo.isGhost=false;
- //index is 1-based, if it is zero it is not initialized, we define it here
- if ( baseInfo.baseIndex<0 ){
- baseInfo.baseIndex=++baseIndex;
- if (!baseInfo.Pcondition) baseInfo.index=++index;}
- }
-}
-
-Vector3r PeriodicFlowEngine::meanVelocity()
-{
- solver->averageRelativeCellVelocity();
- Vector3r meanVel ( 0,0,0 );
- Real volume=0;
- FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
- for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
- //We could also define velocity using cell's center
-// if ( !cell->info().isReal() ) continue;
- if ( cell->info().isGhost ) continue;
- for ( int i=0;i<3;i++ )
- meanVel[i]=meanVel[i]+ ( ( cell->info().averageVelocity() ) [i] * abs ( cell->info().volume() ) );
- volume+=abs ( cell->info().volume() );
- }
- return ( meanVel/volume );
-}
-
-void PeriodicFlowEngine::updateVolumes (FlowSolver& flow)
-{
- if ( debug ) cout << "Updating volumes.............." << endl;
- Real invDeltaT = 1/scene->dt;
- double newVol, dVol;
- epsVolMax=0;
- Real totVol=0;
- Real totDVol=0;
- Real totVol0=0;
- Real totVol1=0;
-
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
- switch ( cell->info().fictious() ) {
- case ( 1 ) :
- newVol = volumeCellSingleFictious ( cell );
- totVol1+=newVol;
- break;
- case ( 0 ) :
- newVol = volumeCell ( cell );
- totVol0+=newVol;
- break;
- default:
- newVol = 0;
- break;
- }
- totVol+=newVol;
- dVol=cell->info().volumeSign * ( newVol - cell->info().volume() );
- totDVol+=dVol;
- epsVolMax = max ( epsVolMax, abs ( dVol/newVol ) );
- cell->info().dv() = dVol * invDeltaT;
- cell->info().volume() = newVol;
- }
- if ( debug ) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<" "<< totVol0<<" "<< totVol1<<endl;
-}
-
-
-void PeriodicFlowEngine::initializeVolumes (FlowSolver& flow)
-{
- FiniteVerticesIterator vertices_end = flow.T[flow.currentTes].Triangulation().finite_vertices_end();
- CGT::CVector Zero ( 0,0,0 );
- for ( FiniteVerticesIterator V_it = flow.T[flow.currentTes].Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++ ) V_it->info().forces=Zero;
-
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
- switch ( cell->info().fictious() )
- {
- case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
- case ( 1 ) : cell->info().volume() = volumeCellSingleFictious ( cell ); break;
- default: cell->info().volume() = 0; break;
- }
- //FIXME: the void volume is negative sometimes, hence crashing...
- if (flow.fluidBulkModulus>0) { cell->info().invVoidVolume() = 1. / (max(0.1*cell->info().volume(),abs(cell->info().volume()) - flow.volumeSolidPore(cell)) ); }
- }
- if ( debug ) cout << "Volumes initialised." << endl;
-}
-
-void PeriodicFlowEngine::buildTriangulation ( double pZero, FlowSolver& flow)
-{
- flow.resetNetwork();
- if (first) flow.currentTes=0;
- else {
- flow.currentTes=!flow.currentTes;
- if ( debug ) cout << "--------RETRIANGULATION-----------" << endl;}
- initSolver(flow);
- addBoundary ( flow );
- if ( debug ) cout << endl << "Added boundaries------" << endl << endl;
- triangulate (flow);
- if ( debug ) cout << endl << "Tesselating------" << endl << endl;
- flow.T[flow.currentTes].compute();
- flow.defineFictiousCells();
-
- //FIXME: this is already done in addBoundary(?)
- boundaryConditions ( flow );
- if ( debug ) cout << endl << "boundaryConditions------" << endl << endl;
- flow.initializePressure ( pZero );
- if ( debug ) cout << endl << "initializePressure------" << endl << endl;
- // Define the ghost cells and add indexes to the cells inside the period (the ones that will contain the pressure unknowns)
- //This must be done after boundary conditions and initialize pressure, else the indexes are not good (not accounting imposedP): FIXME
- unsigned int index=0;
- int baseIndex=-1;
- FlowSolver::Tesselation& Tes = flow.T[flow.currentTes];
- Tes.cellHandles.resize(Tes.Triangulation().number_of_finite_cells());
- const FiniteCellsIterator cellend=Tes.Triangulation().finite_cells_end();
- for ( FiniteCellsIterator cell=Tes.Triangulation().finite_cells_begin(); cell!=cellend; cell++ ){
- locateCell ( cell,index,baseIndex,flow );
- if (flow.errorCode>0) return;
- //Fill this vector than can be later used to speedup loops
- if (!cell->info().isGhost) Tes.cellHandles[cell->info().baseIndex]=cell;
- }
- Tes.cellHandles.resize(baseIndex+1);
-
- if ( debug ) cout << endl << "locateCell------" << endl << endl;
- flow.computePermeability ( );
- porosity = flow.vPoralPorosity/flow.vTotalPorosity;
- flow.tolerance=tolerance;flow.relax=relax;
-
- flow.displayStatistics ();
- //FIXME: check interpolate() for the periodic case, at least use the mean pressure from previous step.
- if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], Tes );
-// if ( !first && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
-
- if ( waveAction ) flow.applySinusoidalPressure ( Tes.Triangulation(), sineMagnitude, sineAverage, 30 );
-
- if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
- if ( debug ) cout << endl << "end buildTri------" << endl << endl;
-}
-
-void PeriodicFlowEngine::preparePShifts()
-{
- CellInfo::gradP = makeCgVect ( gradP );
- CellInfo::hSize[0] = makeCgVect ( scene->cell->hSize.col ( 0 ) );
- CellInfo::hSize[1] = makeCgVect ( scene->cell->hSize.col ( 1 ) );
- CellInfo::hSize[2] = makeCgVect ( scene->cell->hSize.col ( 2 ) );
- CellInfo::deltaP=CGT::CVector (
- CellInfo::hSize[0]*CellInfo::gradP,
- CellInfo::hSize[1]*CellInfo::gradP,
- CellInfo::hSize[2]*CellInfo::gradP );
-}
-
-
-YADE_PLUGIN((PeriodicFlowEngine));
-
-
-///========================================================================================================================
-
-/// Definition of base types for the basic version of FlowEngine
-
-// #ifdef LINSOLV
-// #define _FlowSolver CGT::FlowBoundingSphereLinSolv<CGT::FlowBoundingSphere<FlowTesselation> >
-// #else
-// #define _FlowSolver CGT::FlowBoundingSphere<FlowTesselation>
-// #endif
-//
-// typedef CGT::FlowBoundingSphereLinSolv<CGT::FlowBoundingSphere<FlowTesselation> > tparam;
-// typedef TemplateFlowEngine<_FlowSolver> TFlowEng;
-// REGISTER_SERIALIZABLE(TFlowEng);
-//
-// // template<>
-// // TFlowEng::~TFlowEng()
-// // {
-// // }
-// YADE_PLUGIN((TFlowEng));
+#include <yade/pkg/dem/PeriodicFlowEngine.cpp>
+#include <yade/pkg/dem/DFNFlow.cpp>
#endif //FLOW_ENGINE
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2014-03-21 18:47:45 +0000
+++ pkg/dem/FlowEngine.hpp 2014-03-21 18:47:45 +0000
@@ -1,7 +1,7 @@
/*************************************************************************
* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* Copyright (C) 2012 by Donia Marzougui <donia.marzougui@xxxxxxxxxxxxxxx>*
+
* *
* This program is free software; it is licensed under the terms of the *
* GNU General Public License v2 or later. See file LICENSE for details. *
@@ -43,28 +43,16 @@
CGT::Point makeCgPoint ( const Vector3r& yv ) {return CGT::Point ( yv[0],yv[1],yv[2] );}
Vector3r makeVector3r ( const CGT::Point& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
Vector3r makeVector3r ( const CGT::CVector& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
-/*
-
-typedef TriangulationTypes<FlowVertexInfo,FlowCellInfo> FlowTriangulationTypes;
-typedef CGT::_Tesselation<FlowTriangulationTypes> FlowTesselation;
-
-
-#ifdef LINSOLV
-#define _FlowSolver CGT::FlowBoundingSphereLinSolv<CGT::FlowBoundingSphere<FlowTesselation> >
-#else
-#define _FlowSolver CGT::FlowBoundingSphere<FlowTesselation>
-#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=CGT::FlowBoundingSphere<_Tesselation> >
class TemplateFlowEngine : public PartialEngine
{
public :
#ifdef LINSOLV
- typedef CGT::FlowBoundingSphereLinSolv<typename solverT::Tesselation,solverT> FlowSolver;
+ typedef CGT::FlowBoundingSphereLinSolv<typename solverT::Tesselation,solverT> FlowSolver;
#else
- typedef solverT FlowSolver;
+ typedef solverT FlowSolver;
#endif
-// typedef solverT FlowSolver;
typedef FlowSolver Solver;//FIXME: useless alias, search/replace then remove
typedef typename FlowSolver::Tesselation Tesselation;
typedef typename FlowSolver::RTriangulation RTriangulation;
@@ -86,6 +74,7 @@
vector<posData> positionBufferParallel;//keep the positions from a given step for multithread factorization
//copy positions in a buffer for faster and/or parallel access
void setPositionsBuffer(bool current);
+ virtual void trickPermeability() {};
public :
int retriangulationLastIter;
@@ -403,10 +392,6 @@
double invSumK;
bool isvisited;
- FlowCellInfo& operator= (const std::vector<double> &v) { for (int i=0; i<4;i++) modulePermeability[i]= v[i]; return *this; }
- FlowCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
- FlowCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
-
inline Real& volume (void) {return t;}
inline const Real& invVoidVolume (void) const {return invVoidV;}
inline Real& invVoidVolume (void) {return invVoidV;}
@@ -427,9 +412,6 @@
CVector grainVelocity;
Real volumeIncidentCells;
public:
- FlowVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
- FlowVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
- FlowVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
CVector forces;
bool isGhost;
FlowVertexInfo (void) {isGhost=false;}
@@ -439,127 +421,22 @@
inline const CVector ghostShift (void) {return CGAL::NULL_VECTOR;}
};
-typedef CGT::TriangulationTypes<FlowVertexInfo,FlowCellInfo> FlowTriangulationTypes;
-typedef CGT::_Tesselation<FlowTriangulationTypes> FlowTesselation;
-
-
-#ifdef LINSOLV
-#define _FlowSolver CGT::FlowBoundingSphereLinSolv<FlowTesselation,CGT::FlowBoundingSphere<FlowTesselation> >
-#else
-#define _FlowSolver CGT::FlowBoundingSphere<FlowTesselation>
-#endif
-
-// typedef CGT::FlowBoundingSphereLinSolv<CGT::FlowBoundingSphere<FlowTesselation> > tparam;
-typedef TemplateFlowEngine<FlowCellInfo,FlowVertexInfo> FlowEngine;
-REGISTER_SERIALIZABLE(FlowEngine);
-
-// YADE_PLUGIN((FlowEngine));
-
-///==========================================
-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) {}
- PeriodicCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
- PeriodicCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
-
- 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;}
- virtual bool isReal (void) {return !(isFictious || isGhost);}
-};
-CVector PeriodicCellInfo::hSize[]={CVector(),CVector(),CVector()};
-CVector PeriodicCellInfo::deltaP=CVector();
-CVector PeriodicCellInfo::gradP=CVector();
-
-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;}
- virtual 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
-//CGT::PeriodicFlowLinSolv<CGT::PeriodicTesselation<CGT::_Tesselation<TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> > > >
-
-typedef TemplateFlowEngine< PeriodicCellInfo,
- PeriodicVertexInfo,
- CGT::PeriodicTesselation<CGT::_Tesselation<CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> > >,
- _PeriFlowSolver>
- FlowEngine_PeriodicInfo;
-REGISTER_SERIALIZABLE(FlowEngine_PeriodicInfo);
-
-// template<>
-// TFlowEng::~TFlowEng()
-// {
-// }
-YADE_PLUGIN((FlowEngine_PeriodicInfo));
-
-
-class PeriodicFlowEngine : public FlowEngine_PeriodicInfo
+// To register properly, we need to first instantiate an intermediate class, then inherit from it with correct class names in YADE_CLASS macro
+// The intermediate one would be seen with the name "TemplateFlowEngine" by python, thus it would not work when more than one class are derived, they would all
+// be named "TemplateFlowEngine" ...
+typedef TemplateFlowEngine<FlowCellInfo,FlowVertexInfo> FlowEngineT;
+REGISTER_SERIALIZABLE(FlowEngineT);
+
+class FlowEngine : public FlowEngineT
{
public :
- 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"))
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(FlowEngine,FlowEngineT,"An engine to solve flow problem in saturated granular media. Model description can be found in [Chareyre2012a]_ and [Catalano2013a]_. See the example script FluidCouplingPFV/oedometer.py. More documentation to come.\n\n.. note::Multi-threading seems to work fine for Cholesky decomposition, but it fails for the solve phase in which -j1 is the fastest, here we specify thread numbers independently using :yref:`FlowEngine::numFactorizeThreads` and :yref:`FlowEngine::numSolveThreads`. These multhreading settings are only impacting the behaviour of openblas library and are relatively independant of :yref:`FlowEngine::multithread`. However, the settings have to be globally consistent. For instance, :yref:`multithread<FlowEngine::multithread>` =True with yref:`numFactorizeThreads<FlowEngine::numFactorizeThreads>` = yref:`numSolveThreads<FlowEngine::numSolveThreads>` = 4 implies that openblas will mobilize 8 processors at some point. If the system does not have so many procs. it will hurt performance.",
,,
- 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
+ //nothing special to define here, we simply re-use FlowEngine methods
//.def("meanVelocity",&PeriodicFlowEngine::meanVelocity,"measure the mean velocity in the period")
)
DECLARE_LOGGER;
-
-
};
-
-REGISTER_SERIALIZABLE(PeriodicFlowEngine);
-
-// #endif
+REGISTER_SERIALIZABLE(FlowEngine);
+
=== modified file 'pkg/dem/TesselationWrapper.cpp'
--- pkg/dem/TesselationWrapper.cpp 2014-03-21 18:47:45 +0000
+++ pkg/dem/TesselationWrapper.cpp 2014-03-21 18:47:45 +0000
@@ -109,7 +109,7 @@
if (v==RTriangulation::Vertex_handle())
hint=c;
else {
- v->info() = (const unsigned int) p->second;
+ v->info().setId((const unsigned int) p->second);
//Vh->info().isFictious = false;//false is the default
Tes.maxId = std::max(Tes.maxId,(int) p->second);
Tes.vertexHandles[p->second]=v;
