yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10937
[Branch ~yade-pkg/yade/git-trunk] Rev 4004: fix bug in FlowEngine::updateBCs() (reported by Luc Scholtes, thanks)
------------------------------------------------------------
revno: 4004
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
timestamp: Wed 2014-06-04 19:19:50 +0200
message:
fix bug in FlowEngine::updateBCs() (reported by Luc Scholtes, thanks)
modified:
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphereLinSolv.hpp
lib/triangulation/Network.hpp
pkg/pfv/FlowEngine.hpp
pkg/pfv/FlowEngine.ipp
pkg/pfv/PeriodicFlowEngine.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/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2014-05-26 23:23:29 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2014-06-04 17:19:50 +0000
@@ -31,7 +31,7 @@
//painfull, but we need that for templates inheritance...
using _N::T; using _N::xMin; using _N::xMax; using _N::yMin; using _N::yMax; using _N::zMin; using _N::zMax; using _N::Rmoy; using _N::sectionArea; using _N::Height; using _N::vTotal; using _N::currentTes; using _N::debugOut; using _N::nOfSpheres; using _N::xMinId; using _N::xMaxId; using _N::yMinId; using _N::yMaxId; using _N::zMinId; using _N::zMaxId; using _N::boundsIds; using _N::cornerMin; using _N::cornerMax; using _N::VSolidTot; using _N::Vtotalissimo; using _N::vPoral; using _N::sSolidTot; using _N::vPoralPorosity; using _N::vTotalPorosity; using _N::boundaries; using _N::idOffset; using _N::vtkInfiniteVertices; using _N::vtkInfiniteCells; using _N::num_particles; using _N::boundingCells; using _N::facetVertices; using _N::facetNFictious;
//same for functions
- using _N::defineFictiousCells; using _N::addBoundingPlanes; using _N::boundary;
+ using _N::defineFictiousCells; using _N::addBoundingPlanes; using _N::boundary; using _N::tesselation;
virtual ~FlowBoundingSphere();
FlowBoundingSphere();
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.hpp'
--- lib/triangulation/FlowBoundingSphereLinSolv.hpp 2014-04-10 05:58:01 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.hpp 2014-06-04 17:19:50 +0000
@@ -54,6 +54,7 @@
using FlowType::pressureChanged;
using FlowType::computedOnce;
using FlowType::resetNetwork;
+ using FlowType::tesselation;
//! TAUCS DECs
vector<FiniteCellsIterator> orderedCells;
=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp 2014-03-23 02:08:49 +0000
+++ lib/triangulation/Network.hpp 2014-06-04 17:19:50 +0000
@@ -46,6 +46,8 @@
Tesselation T [2];
bool currentTes;
+ Tesselation& tesselation() {return T[currentTes];};
+
double xMin, xMax, yMin, yMax, zMin, zMax, Rmoy, sectionArea, Height, vTotal;
bool debugOut;
int nOfSpheres;
=== modified file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp 2014-05-27 07:50:01 +0000
+++ pkg/pfv/FlowEngine.hpp 2014-06-04 17:19:50 +0000
@@ -113,10 +113,10 @@
void updateVolumes (Solver& flow);
void initializeVolumes (Solver& flow);
void boundaryConditions(Solver& flow);
- void updateBCs ( Solver& flow ) {
- if (flow.T[flow.currentTes].maxId>0) boundaryConditions(flow);//avoids crash at iteration 0, when the packing is not bounded yet
+ 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");
- flow.pressureChanged=true;}
+ solver->pressureChanged=true;}
void imposeFlux(Vector3r pos, Real flux);
unsigned int imposePressure(Vector3r pos, Real p);
@@ -128,6 +128,7 @@
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();}
@@ -353,6 +354,7 @@
#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
=== modified file 'pkg/pfv/FlowEngine.ipp'
--- pkg/pfv/FlowEngine.ipp 2014-05-26 23:23:29 +0000
+++ pkg/pfv/FlowEngine.ipp 2014-06-04 17:19:50 +0000
@@ -222,8 +222,8 @@
flow.meanKStat = meanKStat;
flow.permeabilityMap = permeabilityMap;
flow.fluidBulkModulus = fluidBulkModulus;
-// flow.T[flow.currentTes].Clear();
- flow.T[flow.currentTes].maxId=-1;
+// 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;
}
@@ -256,16 +256,16 @@
addBoundary ( flow );
triangulate ( flow );
if ( debug ) cout << endl << "Tesselating------" << endl << endl;
- flow.T[flow.currentTes].compute();
+ flow.tesselation().compute();
flow.defineFictiousCells();
// For faster loops on cells define this vector
- flow.T[flow.currentTes].cellHandles.clear();
- flow.T[flow.currentTes].cellHandles.reserve(flow.T[flow.currentTes].Triangulation().number_of_finite_cells());
- FiniteCellsIterator cell_end = flow.T[flow.currentTes].Triangulation().finite_cells_end();
+ 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.T[flow.currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
- flow.T[flow.currentTes].cellHandles.push_back(cell);
+ 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();
@@ -276,9 +276,9 @@
boundaryConditions ( flow );
flow.initializePressure ( pZero );
- if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0 || doInterpolate)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
- if ( waveAction ) flow.applySinusoidalPressure ( flow.T[flow.currentTes].Triangulation(), sineMagnitude, sineAverage, 30 );
- else if (boundaryPressure.size()!=0) flow.applyUserDefinedPressure ( flow.T[flow.currentTes].Triangulation(), boundaryXPos , boundaryPressure);
+ 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 >
@@ -320,7 +320,7 @@
}
}
//FIXME idOffset must be set correctly, not the case here (always 0), then we need walls first or it will fail
- idOffset = flow.T[flow.currentTes].maxId+1;
+ 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 );
@@ -365,7 +365,7 @@
///Using Tesselation wrapper (faster)
// TesselationWrapper TW;
// if (TW.Tes) delete TW.Tes;
-// TW.Tes = &(flow.T[flow.currentTes]);//point to the current Tes we have in Flowengine
+// TW.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
@@ -374,27 +374,27 @@
if ( !b.exists ) continue;
if ( b.isSphere ) {
if (b.id==ignoredBody) continue;
- flow.T[flow.currentTes].insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
+ flow.tesselation().insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
}
- flow.T[flow.currentTes].redirected=true;//By inserting one-by-one, we already redirected
- flow.shearLubricationForces.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.shearLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.pumpLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.twistLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.shearLubricationBodyStress.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.normalLubricationForce.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.normalLubricationBodyStress.resize ( flow.T[flow.currentTes].maxId+1 );
+ 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.T[flow.currentTes].Triangulation().finite_vertices_end();
+ FiniteVerticesIterator vertices_end = flow.tesselation().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;
+ 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.T[flow.currentTes].cellHandles)
+ FOREACH(CellHandle& cell, flow.tesselation().cellHandles)
{
switch ( cell->info().fictious() )
{
@@ -448,12 +448,12 @@
epsVolMax=0;
Real totVol=0; Real totDVol=0;
#ifdef YADE_OPENMP
- const long size=flow.T[flow.currentTes].cellHandles.size();
+ 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.T[flow.currentTes].cellHandles[i];
+ CellHandle& cell = flow.tesselation().cellHandles[i];
#else
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
+ FOREACH(CellHandle& cell, flow.tesselation().cellHandles){
#endif
double newVol, dVol;
switch ( cell->info().fictious() ) {
@@ -592,7 +592,7 @@
for ( unsigned int k=0; k<flow.normalLubricationBodyStress.size(); k++) flow.normalLubricationBodyStress[k]=Matrix3r::Zero();
typedef typename Solver::Tesselation Tesselation;
- const Tesselation& Tes = flow.T[flow.currentTes];
+ 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();
@@ -726,6 +726,24 @@
}
}
+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 */
=== modified file 'pkg/pfv/PeriodicFlowEngine.cpp'
--- pkg/pfv/PeriodicFlowEngine.cpp 2014-04-17 15:10:23 +0000
+++ pkg/pfv/PeriodicFlowEngine.cpp 2014-06-04 17:19:50 +0000
@@ -85,7 +85,6 @@
Real volumeCellSingleFictious (CellHandle cell);
inline void locateCell(CellHandle baseCell, unsigned int& index, int& baseIndex, FlowSolver& flow, unsigned int count=0);
- Vector3r meanVelocity();
virtual ~PeriodicFlowEngine();
@@ -232,7 +231,7 @@
void PeriodicFlowEngine::triangulate( FlowSolver& flow )
{
- Tesselation& Tes = flow.T[flow.currentTes];
+ Tesselation& Tes = flow.tesselation();
vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
FOREACH ( const posData& b, buffer ) {
if ( !b.exists || !b.isSphere || b.id==ignoredBody) continue;
@@ -332,7 +331,7 @@
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();
+ RTriangulation& Tri = flow.tesselation().Triangulation();
Vector3r center ( 0,0,0 );
Vector3i period;
@@ -408,23 +407,6 @@
}
}
-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;
@@ -436,7 +418,7 @@
Real totVol0=0;
Real totVol1=0;
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
+ FOREACH(CellHandle& cell, flow.tesselation().cellHandles){
switch ( cell->info().fictious() ) {
case ( 1 ) :
newVol = volumeCellSingleFictious ( cell );
@@ -463,11 +445,11 @@
void PeriodicFlowEngine::initializeVolumes (FlowSolver& flow)
{
- FiniteVerticesIterator vertices_end = flow.T[flow.currentTes].Triangulation().finite_vertices_end();
+ FiniteVerticesIterator vertices_end = flow.tesselation().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;
+ 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.T[flow.currentTes].cellHandles){
+ FOREACH(CellHandle& cell, flow.tesselation().cellHandles){
switch ( cell->info().fictious() )
{
case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
@@ -492,7 +474,7 @@
if ( debug ) cout << endl << "Added boundaries------" << endl << endl;
triangulate (flow);
if ( debug ) cout << endl << "Tesselating------" << endl << endl;
- flow.T[flow.currentTes].compute();
+ flow.tesselation().compute();
flow.defineFictiousCells();
//FIXME: this is already done in addBoundary(?)
@@ -504,7 +486,7 @@
//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];
+ FlowSolver::Tesselation& Tes = flow.tesselation();
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++ ){
@@ -523,7 +505,7 @@
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 || doInterpolate)) flow.interpolate ( flow.T[!flow.currentTes], Tes );
-// if ( !first && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
+// if ( !first && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], flow.tesselation() );
if ( waveAction ) flow.applySinusoidalPressure ( Tes.Triangulation(), sineMagnitude, sineAverage, 30 );
