yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10593
[Branch ~yade-pkg/yade/git-trunk] Rev 3847: make FlowEngine a template yade class
------------------------------------------------------------
revno: 3847
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Fri 2014-03-21 19:47:45 +0100
message:
make FlowEngine a template yade class
added:
pkg/dem/FlowEngine.ipp
--
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
=== added file 'pkg/dem/FlowEngine.ipp'
--- pkg/dem/FlowEngine.ipp 1970-01-01 00:00:00 +0000
+++ pkg/dem/FlowEngine.ipp 2014-03-21 18:47:45 +0000
@@ -0,0 +1,729 @@
+/*************************************************************************
+* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
+* Copyright (C) 2009 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+#ifdef YADE_CGAL
+
+#ifdef FLOW_ENGINE
+#include<yade/core/Scene.hpp>
+#include<yade/lib/base/Math.hpp>
+#include<yade/pkg/dem/TesselationWrapper.hpp>
+#include<yade/pkg/common/Sphere.hpp>
+#include<yade/pkg/common/Wall.hpp>
+#include<yade/pkg/common/Box.hpp>
+#include <sys/stat.h>
+#include <sys/types.h>
+#include <boost/thread.hpp>
+#include <boost/date_time.hpp>
+#include <boost/date_time/posix_time/posix_time.hpp>
+
+#ifdef LINSOLV
+#include <cholmod.h>
+#endif
+
+#include "FlowEngine.hpp"
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::~TemplateFlowEngine() {}
+
+// YADE_PLUGIN((TFlowEng));
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+unsigned int TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposePressure(Vector3r pos, Real p)
+{
+// if (!flow) LOG_ERROR("no flow defined yet, run at least one iter");
+ solver->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),p) );
+ //force immediate update of boundary conditions
+ updateTriangulation=true;
+ return solver->imposedP.size()-1;
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::action()
+{
+ if ( !isActivated ) return;
+ timingDeltas->start();
+ setPositionsBuffer(true);
+ timingDeltas->checkpoint ( "Position buffer" );
+ if (first) {
+ if (multithread) setPositionsBuffer(false);
+ buildTriangulation(pZero,*solver);
+ initializeVolumes(*solver);
+ backgroundSolver=solver;
+ backgroundCompleted=true;
+ }
+ solver->ompThreads = ompThreads>0? ompThreads : omp_get_max_threads();
+
+ timingDeltas->checkpoint ( "Triangulating" );
+ updateVolumes ( *solver );
+ timingDeltas->checkpoint ( "Update_Volumes" );
+
+ epsVolCumulative += epsVolMax;
+ retriangulationLastIter++;
+ if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
+ (defTolerance>0 && epsVolCumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval;
+
+ ///compute flow and and forces here
+ if (pressureForce){
+ solver->gaussSeidel(scene->dt);
+ timingDeltas->checkpoint ( "Gauss-Seidel (includes matrix construct and factorization in single-thread mode)" );
+ solver->computeFacetForcesWithCache();}
+ timingDeltas->checkpoint ( "compute_Forces" );
+ ///Application of vicscous forces
+ scene->forces.sync();
+ timingDeltas->checkpoint ( "forces.sync()" );
+ computeViscousForces ( *solver );
+ timingDeltas->checkpoint ( "viscous forces" );
+ Vector3r force;
+ Vector3r torque;
+ FiniteVerticesIterator verticesEnd = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
+ for ( FiniteVerticesIterator vIt = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); vIt != verticesEnd; vIt++ ) {
+ force = pressureForce ? Vector3r ( vIt->info().forces[0],vIt->info().forces[1],vIt->info().forces[2] ): Vector3r(0,0,0);
+ torque = Vector3r(0,0,0);
+ if (shearLubrication || viscousShear){
+ force = force + solver->shearLubricationForces[vIt->info().id()];
+ torque = torque + solver->shearLubricationTorques[vIt->info().id()];
+ if (pumpTorque)
+ torque = torque + solver->pumpLubricationTorques[vIt->info().id()];
+ }
+ if (twistTorque)
+ torque = torque + solver->twistLubricationTorques[vIt->info().id()];
+ if (normalLubrication)
+ force = force + solver-> normalLubricationForce[vIt->info().id()];
+ scene->forces.addForce ( vIt->info().id(), force);
+ scene->forces.addTorque ( vIt->info().id(), torque);
+ }
+ ///End compute flow and forces
+ timingDeltas->checkpoint ( "Applying Forces" );
+ int sleeping = 0;
+ if (multithread && !first) {
+ while (updateTriangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/
+ sleeping++;
+ boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
+ if (debug && sleeping) cerr<<"sleeping..."<<sleeping<<endl;
+ if (updateTriangulation || (ellapsedIter>(0.5*meshUpdateInterval) && backgroundCompleted)) {
+ if (debug) cerr<<"switch flow solver"<<endl;
+ if (useSolver==0) LOG_ERROR("background calculations not available for Gauss-Seidel");
+ if (fluidBulkModulus>0) solver->interpolate (solver->T[solver->currentTes], backgroundSolver->T[backgroundSolver->currentTes]);
+ solver=backgroundSolver;
+ backgroundSolver = shared_ptr<FlowSolver> (new FlowSolver);
+ //Copy imposed pressures/flow from the old solver
+ backgroundSolver->imposedP = vector<pair<CGT::Point,Real> >(solver->imposedP);
+ backgroundSolver->imposedF = vector<pair<CGT::Point,Real> >(solver->imposedF);
+ if (debug) cerr<<"switched"<<endl;
+ setPositionsBuffer(false);//set "parallel" buffer for background calculation
+ backgroundCompleted=false;
+ retriangulationLastIter=ellapsedIter;
+ updateTriangulation=false;
+ epsVolCumulative=0;
+ ellapsedIter=0;
+ boost::thread workerThread(&TemplateFlowEngine::backgroundAction,this);
+ workerThread.detach();
+ if (debug) cerr<<"backgrounded"<<endl;
+ initializeVolumes(*solver);
+ computeViscousForces(*solver);
+ if (debug) cerr<<"volumes initialized"<<endl;
+ }
+ else {
+ if (debug && !backgroundCompleted) cerr<<"still computing solver in the background, ellapsedIter="<<ellapsedIter<<endl;
+ ellapsedIter++;
+ }
+ } else {
+ if (updateTriangulation && !first) {
+ buildTriangulation (pZero, *solver);
+ initializeVolumes(*solver);
+ computeViscousForces(*solver);
+ updateTriangulation = false;
+ epsVolCumulative=0;
+ retriangulationLastIter=0;
+ ReTrg++;}
+ }
+ first=false;
+ timingDeltas->checkpoint ( "triangulate + init volumes" );
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::backgroundAction()
+{
+ if (useSolver<1) {LOG_ERROR("background calculations not available for Gauss-Seidel"); return;}
+ buildTriangulation ( pZero,*backgroundSolver );
+ //FIXME: GS is computing too much, we need only matrix factorization in fact
+ backgroundSolver->gaussSeidel(scene->dt);
+ //FIXME(2): and here we need only cached variables, not forces
+ backgroundSolver->computeFacetForcesWithCache(/*onlyCache?*/ true);
+// boost::this_thread::sleep(boost::posix_time::seconds(5));
+ backgroundCompleted = true;
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::boundaryConditions ( Solver& flow )
+{
+
+ for (int k=0;k<6;k++) {
+ flow.boundary (wallIds[k]).flowCondition=!bndCondIsPressure[k];
+ flow.boundary (wallIds[k]).value=bndCondValue[k];
+ flow.boundary (wallIds[k]).velocity = boundaryVelocity[k];//FIXME: needs correct implementation, maybe update the cached pos/vel?
+ }
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setImposedPressure ( unsigned int cond, Real p)
+{
+ if ( cond>=solver->imposedP.size() ) LOG_ERROR ( "Setting p with cond higher than imposedP size." );
+ solver->imposedP[cond].second=p;
+ //force immediate update of boundary conditions
+ solver->pressureChanged=true;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposeFlux ( Vector3r pos, Real flux){
+ solver->imposedF.push_back ( pair<CGT::Point,Real> ( CGT::Point ( pos[0],pos[1],pos[2] ), flux ) );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedPressure () { solver->imposedP.clear(); solver->IPCells.clear();}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedFlux () { solver->imposedF.clear(); solver->IFCells.clear();}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::getCellFlux ( unsigned int cond)
+{
+ if ( cond>=solver->imposedP.size() ) {LOG_ERROR ( "Getting flux with cond higher than imposedP size." ); return 0;}
+ double flux=0;
+ typename Solver::CellHandle& cell= solver->IPCells[cond];
+ for ( int ngb=0;ngb<4;ngb++ ) {
+ /*if (!cell->neighbor(ngb)->info().Pcondition)*/
+ flux+= cell->info().kNorm() [ngb]* ( cell->info().p()-cell->neighbor ( ngb )->info().p() );
+ }
+ return flux+cell->info().dv();
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::initSolver ( FlowSolver& flow )
+{
+ flow.Vtotalissimo=0; flow.VSolidTot=0; flow.vPoral=0; flow.sSolidTot=0;
+ flow.slipBoundary=slipBoundary;
+ flow.kFactor = permeabilityFactor;
+ flow.debugOut = debug;
+ flow.useSolver = useSolver;
+ #ifdef EIGENSPARSE_LIB
+ flow.numSolveThreads = numSolveThreads;
+ flow.numFactorizeThreads = numFactorizeThreads;
+ #endif
+ flow.meanKStat = meanKStat;
+ flow.viscosity = viscosity;
+ flow.tolerance=tolerance;
+ flow.relax=relax;
+ flow.clampKValues = clampKValues;
+ flow.maxKdivKmean = maxKdivKmean;
+ flow.minKdivKmean = minKdivKmean;
+ flow.meanKStat = meanKStat;
+ flow.permeabilityMap = permeabilityMap;
+ flow.fluidBulkModulus = fluidBulkModulus;
+ flow.T[flow.currentTes].Clear();
+ flow.T[flow.currentTes].maxId=-1;
+ flow.xMin = 1000.0, flow.xMax = -10000.0, flow.yMin = 1000.0, flow.yMax = -10000.0, flow.zMin = 1000.0, flow.zMax = -10000.0;
+}
+
+#ifdef LINSOLV
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setForceMetis ( Solver& flow, bool force )
+{
+ if (force) {
+ flow.eSolver.cholmod().nmethods=1;
+ flow.eSolver.cholmod().method[0].ordering=CHOLMOD_METIS;
+ } else cholmod_defaults(&(flow.eSolver.cholmod()));
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+bool TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::getForceMetis ( Solver& flow ) {return (flow.eSolver.cholmod().nmethods==1);}
+#endif
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( Solver& flow )
+{
+ buildTriangulation ( 0.f,flow );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( double pZero, Solver& flow )
+{
+ flow.resetNetwork();
+ if (first) flow.currentTes=0;
+ else {
+ flow.currentTes=!flow.currentTes;
+ if (debug) cout << "--------RETRIANGULATION-----------" << endl;
+ }
+
+ initSolver(flow);
+
+ addBoundary ( flow );
+ triangulate ( flow );
+ if ( debug ) cout << endl << "Tesselating------" << endl << endl;
+ flow.T[flow.currentTes].compute();
+
+ flow.defineFictiousCells();
+ // For faster loops on cells define this vector
+ flow.T[flow.currentTes].cellHandles.clear();
+ flow.T[flow.currentTes].cellHandles.reserve(flow.T[flow.currentTes].Triangulation().number_of_finite_cells());
+ FiniteCellsIterator cell_end = flow.T[flow.currentTes].Triangulation().finite_cells_end();
+ int k=0;
+ for ( FiniteCellsIterator cell = flow.T[flow.currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
+ flow.T[flow.currentTes].cellHandles.push_back(cell);
+ cell->info().id=k++;}//define unique numbering now, corresponds to position in cellHandles
+ flow.displayStatistics ();
+ flow.computePermeability();
+ //This virtual function does nothing yet, derived class may overload it to make permeability different (see DFN engine)
+ trickPermeability();
+ porosity = flow.vPoralPorosity/flow.vTotalPorosity;
+
+ boundaryConditions ( flow );
+ flow.initializePressure ( pZero );
+
+ if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
+ if ( waveAction ) flow.applySinusoidalPressure ( flow.T[flow.currentTes].Triangulation(), sineMagnitude, sineAverage, 30 );
+ if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setPositionsBuffer(bool current)
+{
+ vector<posData>& buffer = current? positionBufferCurrent : positionBufferParallel;
+ buffer.clear();
+ buffer.resize(scene->bodies->size());
+ shared_ptr<Sphere> sph ( new Sphere );
+ const int Sph_Index = sph->getClassIndexStatic();
+ FOREACH ( const shared_ptr<Body>& b, *scene->bodies ) {
+ if (!b || ignoredBody==b->getId()) continue;
+ posData& dat = buffer[b->getId()];
+ dat.id=b->getId();
+ dat.pos=b->state->pos;
+ dat.isSphere= (b->shape->getClassIndex() == Sph_Index);
+ if (dat.isSphere) dat.radius = YADE_CAST<Sphere*>(b->shape.get())->radius;
+ dat.exists=true;
+ }
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::addBoundary ( Solver& flow )
+{
+ vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
+ solver->xMin = Mathr::MAX_REAL, solver->xMax = -Mathr::MAX_REAL, solver->yMin = Mathr::MAX_REAL, solver->yMax = -Mathr::MAX_REAL, solver->zMin = Mathr::MAX_REAL, solver->zMax = -Mathr::MAX_REAL;
+ FOREACH ( const posData& b, buffer ) {
+ if ( !b.exists ) continue;
+ if ( b.isSphere ) {
+ const Real& rad = b.radius;
+ const Real& x = b.pos[0];
+ const Real& y = b.pos[1];
+ const Real& z = b.pos[2];
+ flow.xMin = min ( flow.xMin, x-rad );
+ flow.xMax = max ( flow.xMax, x+rad );
+ flow.yMin = min ( flow.yMin, y-rad );
+ flow.yMax = max ( flow.yMax, y+rad );
+ flow.zMin = min ( flow.zMin, z-rad );
+ flow.zMax = max ( flow.zMax, z+rad );
+ }
+ }
+ //FIXME idOffset must be set correctly, not the case here (always 0), then we need walls first or it will fail
+ idOffset = flow.T[flow.currentTes].maxId+1;
+ flow.idOffset = idOffset;
+ flow.sectionArea = ( flow.xMax - flow.xMin ) * ( flow.zMax-flow.zMin );
+ flow.vTotal = ( flow.xMax-flow.xMin ) * ( flow.yMax-flow.yMin ) * ( flow.zMax-flow.zMin );
+ flow.yMinId=wallIds[ymin];
+ flow.yMaxId=wallIds[ymax];
+ flow.xMaxId=wallIds[xmax];
+ flow.xMinId=wallIds[xmin];
+ flow.zMinId=wallIds[zmin];
+ flow.zMaxId=wallIds[zmax];
+
+ //FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
+ flow.boundsIds[0]= &flow.xMinId;
+ flow.boundsIds[1]= &flow.xMaxId;
+ flow.boundsIds[2]= &flow.yMinId;
+ flow.boundsIds[3]= &flow.yMaxId;
+ flow.boundsIds[4]= &flow.zMinId;
+ flow.boundsIds[5]= &flow.zMaxId;
+
+ for (int k=0;k<6;k++) flow.boundary ( *flow.boundsIds[k] ).useMaxMin = boundaryUseMaxMin[k];
+
+ flow.cornerMin = CGT::Point ( flow.xMin, flow.yMin, flow.zMin );
+ flow.cornerMax = CGT::Point ( flow.xMax, flow.yMax, flow.zMax );
+
+ //assign BCs types and values
+ boundaryConditions ( flow );
+
+ double center[3];
+ for ( int i=0; i<6; i++ ) {
+ if ( *flow.boundsIds[i]<0 ) continue;
+ CGT::CVector Normal ( normal[i].x(), normal[i].y(), normal[i].z() );
+ if ( flow.boundary ( *flow.boundsIds[i] ).useMaxMin ) flow.addBoundingPlane(Normal, *flow.boundsIds[i] );
+ else {
+ for ( int h=0;h<3;h++ ) center[h] = buffer[*flow.boundsIds[i]].pos[h];
+// cerr << "id="<<*flow.boundsIds[i] <<" center="<<center[0]<<","<<center[1]<<","<<center[2]<<endl;
+ flow.addBoundingPlane ( center, wallThickness, Normal,*flow.boundsIds[i] );
+ }
+ }
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::triangulate ( Solver& flow )
+{
+///Using Tesselation wrapper (faster)
+// TesselationWrapper TW;
+// if (TW.Tes) delete TW.Tes;
+// TW.Tes = &(flow.T[flow.currentTes]);//point to the current Tes we have in Flowengine
+// TW.insertSceneSpheres();//TW is now really inserting in TemplateFlowEngine, using the faster insert(begin,end)
+// TW.Tes = NULL;//otherwise, Tes would be deleted by ~TesselationWrapper() at the end of the function.
+///Using one-by-one insertion
+ vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
+ FOREACH ( const posData& b, buffer ) {
+ if ( !b.exists ) continue;
+ if ( b.isSphere ) {
+ if (b.id==ignoredBody) continue;
+ flow.T[flow.currentTes].insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
+ }
+ flow.T[flow.currentTes].redirected=true;//By inserting one-by-one, we already redirected
+ flow.shearLubricationForces.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.shearLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.pumpLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.twistLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.shearLubricationBodyStress.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.normalLubricationForce.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.normalLubricationBodyStress.resize ( flow.T[flow.currentTes].maxId+1 );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::initializeVolumes ( Solver& flow )
+{
+ typedef typename Solver::FiniteVerticesIterator FiniteVerticesIterator;
+
+ FiniteVerticesIterator vertices_end = flow.T[flow.currentTes].Triangulation().finite_vertices_end();
+ CGT::CVector Zero(0,0,0);
+ for (FiniteVerticesIterator V_it = flow.T[flow.currentTes].Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero;
+
+ FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles)
+ {
+ switch ( cell->info().fictious() )
+ {
+ case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
+ case ( 1 ) : cell->info().volume() = volumeCellSingleFictious ( cell ); break;
+ case ( 2 ) : cell->info().volume() = volumeCellDoubleFictious ( cell ); break;
+ case ( 3 ) : cell->info().volume() = volumeCellTripleFictious ( cell ); break;
+ default: break;
+ }
+ if (flow.fluidBulkModulus>0) { cell->info().invVoidVolume() = 1 / ( abs(cell->info().volume()) - flow.volumeSolidPore(cell) ); }
+ }
+ if (debug) cout << "Volumes initialised." << endl;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::averageRealCellVelocity()
+{
+ solver->averageRelativeCellVelocity();
+ Vector3r Vel ( 0,0,0 );
+ //AVERAGE CELL VELOCITY
+ FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
+ for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
+ for ( int g=0;g<4;g++ ) {
+ if ( !cell->vertex ( g )->info().isFictious ) {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ for ( int i=0;i<3;i++ ) Vel[i] = Vel[i] + sph->state->vel[i]/4;
+ }
+ }
+ RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ CGT::Point pos_av_facet;
+ double volume_facet_translation = 0;
+ CGT::CVector Vel_av ( Vel[0], Vel[1], Vel[2] );
+ for ( int i=0; i<4; i++ ) {
+ volume_facet_translation = 0;
+ if ( !Tri.is_infinite ( cell->neighbor ( i ) ) ) {
+ CGT::CVector Surfk = cell->info()-cell->neighbor ( i )->info();
+ Real area = sqrt ( Surfk.squared_length() );
+ Surfk = Surfk/area;
+ CGT::CVector branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
+ pos_av_facet = ( CGT::Point ) cell->info() + ( branch*Surfk ) *Surfk;
+ volume_facet_translation += Vel_av*cell->info().facetSurfaces[i];
+ cell->info().averageVelocity() = cell->info().averageVelocity() - volume_facet_translation/cell->info().volume() * ( pos_av_facet-CGAL::ORIGIN );
+ }
+ }
+ }
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::updateVolumes ( Solver& flow )
+{
+ if ( debug ) cout << "Updating volumes.............." << endl;
+ Real invDeltaT = 1/scene->dt;
+ epsVolMax=0;
+ Real totVol=0; Real totDVol=0;
+ #ifdef YADE_OPENMP
+ const long size=flow.T[flow.currentTes].cellHandles.size();
+ #pragma omp parallel for num_threads(ompThreads>0 ? ompThreads : 1)
+ for(long i=0; i<size; i++){
+ CellHandle& cell = flow.T[flow.currentTes].cellHandles[i];
+ #else
+ FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
+ #endif
+ double newVol, dVol;
+ switch ( cell->info().fictious() ) {
+ case ( 3 ) : newVol = volumeCellTripleFictious ( cell ); break;
+ case ( 2 ) : newVol = volumeCellDoubleFictious ( cell ); break;
+ case ( 1 ) : newVol = volumeCellSingleFictious ( cell ); break;
+ case ( 0 ) : newVol = volumeCell (cell ); break;
+ default: newVol = 0; break;}
+ dVol=cell->info().volumeSign* ( newVol - cell->info().volume() );
+ cell->info().dv() = dVol*invDeltaT;
+ cell->info().volume() = newVol;
+ if (defTolerance>0) { //if the criterion is not used, then we skip these updates a save a LOT of time when Nthreads > 1
+ #pragma omp atomic
+ totVol+=newVol;
+ #pragma omp atomic
+ totDVol+=abs(dVol);}
+ }
+ if (defTolerance>0) epsVolMax = totDVol/totVol;
+ for (unsigned int n=0; n<flow.imposedF.size();n++) {
+ flow.IFCells[n]->info().dv()+=flow.imposedF[n].second;
+ flow.IFCells[n]->info().Pcondition=false;}
+ if ( debug ) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<endl;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellSingleFictious ( Cellhandle cell )
+{
+ Vector3r V[3];
+ int b=0;
+ int w=0;
+ cell->info().volumeSign=1;
+ Real Wall_coordinate=0;
+
+ for ( int y=0;y<4;y++ ) {
+ if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
+ V[w]=positionBufferCurrent[cell->vertex ( y )->info().id()].pos;
+ w++;
+ } else {
+ b = cell->vertex ( y )->info().id();
+ const shared_ptr<Body>& wll = Body::byId ( b , scene );
+ if ( !solver->boundary ( b ).useMaxMin ) Wall_coordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wallThickness/2.;
+ else Wall_coordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate];
+ }
+ }
+ Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
+ return abs ( Volume );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellDoubleFictious ( Cellhandle cell )
+{
+ Vector3r A=Vector3r::Zero(), AS=Vector3r::Zero(),B=Vector3r::Zero(), BS=Vector3r::Zero();
+
+ cell->info().volumeSign=1;
+ int b[2];
+ int coord[2];
+ Real Wall_coordinate[2];
+ int j=0;
+ bool first_sph=true;
+
+ for ( int g=0;g<4;g++ ) {
+ if ( cell->vertex ( g )->info().isFictious ) {
+ b[j] = cell->vertex ( g )->info().id();
+ coord[j]=solver->boundary ( b[j] ).coordinate;
+ if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = positionBufferCurrent[b[j]].pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2.;
+ else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
+ j++;
+ } else if ( first_sph ) {
+ A=AS=/*AT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;
+ first_sph=false;
+ } else {
+ B=BS=/*BT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;;
+ }
+ }
+ AS[coord[0]]=BS[coord[0]] = Wall_coordinate[0];
+
+ //first pyramid with triangular base (A,B,BS)
+ Real Vol1=0.5* ( ( A-BS ).cross ( B-BS ) ) [coord[1]]* ( 0.333333333* ( 2*B[coord[1]]+A[coord[1]] )-Wall_coordinate[1] );
+ //second pyramid with triangular base (A,AS,BS)
+ Real Vol2=0.5* ( ( AS-BS ).cross ( A-BS ) ) [coord[1]]* ( 0.333333333* ( B[coord[1]]+2*A[coord[1]] )-Wall_coordinate[1] );
+ return abs ( Vol1+Vol2 );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellTripleFictious ( Cellhandle cell )
+{
+ Vector3r A;
+
+ int b[3];
+ int coord[3];
+ Real Wall_coordinate[3];
+ int j=0;
+ cell->info().volumeSign=1;
+
+ for ( int g=0;g<4;g++ ) {
+ if ( cell->vertex ( g )->info().isFictious ) {
+ b[j] = cell->vertex ( g )->info().id();
+ coord[j]=solver->boundary ( b[j] ).coordinate;
+ const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
+ if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = wll->state->pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2.;
+ else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
+ j++;
+ } else {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ A= ( sph->state->pos );
+ }
+ }
+ Real Volume = ( A[coord[0]]-Wall_coordinate[0] ) * ( A[coord[1]]-Wall_coordinate[1] ) * ( A[coord[2]]-Wall_coordinate[2] );
+ return abs ( Volume );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCell ( Cellhandle cell )
+{
+ static const Real inv6 = 1/6.;
+ const Vector3r& p0 = positionBufferCurrent[cell->vertex ( 0 )->info().id()].pos;
+ const Vector3r& p1 = positionBufferCurrent[cell->vertex ( 1 )->info().id()].pos;
+ const Vector3r& p2 = positionBufferCurrent[cell->vertex ( 2 )->info().id()].pos;
+ const Vector3r& p3 = positionBufferCurrent[cell->vertex ( 3 )->info().id()].pos;
+ Real volume = inv6 * ((p0-p1).cross(p0-p2)).dot(p0-p3);
+ if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
+ return volume;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::computeViscousForces ( Solver& flow )
+{
+ if (normalLubrication || shearLubrication || viscousShear){
+ if ( debug ) cout << "Application of viscous forces" << endl;
+ if ( debug ) cout << "Number of edges = " << flow.edgeIds.size() << endl;
+ for ( unsigned int k=0; k<flow.shearLubricationForces.size(); k++ ) flow.shearLubricationForces[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.shearLubricationTorques.size(); k++ ) flow.shearLubricationTorques[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.pumpLubricationTorques.size(); k++ ) flow.pumpLubricationTorques[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.twistLubricationTorques.size(); k++ ) flow.twistLubricationTorques[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.shearLubricationBodyStress.size(); k++) flow.shearLubricationBodyStress[k]=Matrix3r::Zero();
+ for ( unsigned int k=0; k<flow.normalLubricationForce.size(); k++ ) flow.normalLubricationForce[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.normalLubricationBodyStress.size(); k++) flow.normalLubricationBodyStress[k]=Matrix3r::Zero();
+
+ typedef typename Solver::Tesselation Tesselation;
+ const Tesselation& Tes = flow.T[flow.currentTes];
+ flow.deltaShearVel.clear(); flow.normalV.clear(); flow.deltaNormVel.clear(); flow.surfaceDistance.clear(); flow.onlySpheresInteractions.clear(); flow.normalStressInteraction.clear(); flow.shearStressInteraction.clear();
+
+
+ for ( int i=0; i< ( int ) flow.edgeIds.size(); i++ ) {
+ const int& id1 = flow.edgeIds[i].first;
+ const int& id2 = flow.edgeIds[i].second;
+
+ int hasFictious= Tes.vertex ( id1 )->info().isFictious + Tes.vertex ( id2 )->info().isFictious;
+ if (hasFictious>0 or id1==id2) continue;
+ const shared_ptr<Body>& sph1 = Body::byId ( id1, scene );
+ const shared_ptr<Body>& sph2 = Body::byId ( id2, scene );
+ Sphere* s1=YADE_CAST<Sphere*> ( sph1->shape.get() );
+ Sphere* s2=YADE_CAST<Sphere*> ( sph2->shape.get() );
+ const Real& r1 = s1->radius;
+ const Real& r2 = s2->radius;
+ Vector3r deltaV; Real deltaNormV; Vector3r deltaShearV;
+ Vector3r O1O2Vector; Real O1O2; Vector3r normal; Real surfaceDist; Vector3r O1CVector; Vector3r O2CVector;Real meanRad ;Real Rh; Vector3r deltaAngVel; Vector3r deltaShearAngVel;
+ Vector3r shearLubF; Vector3r normaLubF; Vector3r pumpT; Vector3r deltaAngNormVel; Vector3r twistT; Vector3r angVel1; Vector3r angVel2;
+ //FIXME: if periodic and velGrad!=0, then deltaV should account for velGrad, not the case currently
+ if ( !hasFictious ){
+ O1O2Vector = sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()) - sph1->state->pos - makeVector3r(Tes.vertex(id1)->info().ghostShift());
+ O1O2 = O1O2Vector.norm();
+ normal= (O1O2Vector/O1O2);
+ surfaceDist = O1O2 - r2 - r1;
+ O1CVector = (O1O2/2. + (pow(r1,2) - pow(r2,2)) / (2.*O1O2))*normal;
+ O2CVector = -(O1O2Vector - O1CVector);
+ meanRad = (r2 + r1)/2.;
+ Rh = (r1 < r2)? surfaceDist + 0.45 * r1 : surfaceDist + 0.45 * r2;
+ deltaV = (sph2->state->vel + sph2->state->angVel.cross(-r2 * normal)) - (sph1->state->vel+ sph1->state->angVel.cross(r1 * normal));
+ angVel1 = sph1->state->angVel;
+ angVel2 = sph2->state->angVel;
+ deltaAngVel = sph2->state->angVel - sph1->state->angVel;
+
+ } else {
+ if ( hasFictious==1 ) {//for the fictious sphere, use velocity of the boundary, not of the body
+ bool v1fictious = Tes.vertex ( id1 )->info().isFictious;
+ int bnd = v1fictious? id1 : id2;
+ int coord = flow.boundary(bnd).coordinate;
+ O1O2 = v1fictious ? abs((sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]) : abs((sph1->state->pos + makeVector3r(Tes.vertex(id1)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]);
+ if(v1fictious)
+ normal = makeVector3r(flow.boundary(id1).normal);
+ else
+ normal = -makeVector3r(flow.boundary(id2).normal);
+ O1O2Vector = O1O2 * normal;
+ meanRad = v1fictious ? r2:r1;
+ surfaceDist = O1O2 - meanRad;
+ if (v1fictious){
+ O1CVector = Vector3r::Zero();
+ O2CVector = - O1O2Vector;}
+ else{
+ O1CVector = O1O2Vector;
+ O2CVector = Vector3r::Zero();}
+
+ Rh = surfaceDist + 0.45 * meanRad;
+ Vector3r v1 = ( Tes.vertex ( id1 )->info().isFictious ) ? flow.boundary ( id1 ).velocity:sph1->state->vel + sph1->state->angVel.cross(r1 * normal);
+ Vector3r v2 = ( Tes.vertex ( id2 )->info().isFictious ) ? flow.boundary ( id2 ).velocity:sph2->state->vel + sph2->state->angVel.cross(-r2 * (normal));
+ deltaV = v2-v1;
+ angVel1 = ( Tes.vertex ( id1 )->info().isFictious ) ? Vector3r::Zero() : sph1->state->angVel;
+ angVel2 = ( Tes.vertex ( id2 )->info().isFictious ) ? Vector3r::Zero() : sph2->state->angVel;
+ deltaAngVel = angVel2 - angVel1;
+ }
+ }
+ deltaShearV = deltaV - ( normal.dot ( deltaV ) ) *normal;
+ deltaShearAngVel = deltaAngVel - ( normal.dot ( deltaAngVel ) ) *normal;
+ flow.deltaShearVel.push_back(deltaShearV);
+ flow.normalV.push_back(normal);
+ flow.surfaceDistance.push_back(max(surfaceDist, 0.) + eps*meanRad);
+
+ /// Compute the shear Lubrication force and torque on each particle
+
+ if (shearLubrication)
+ shearLubF = flow.computeShearLubricationForce(deltaShearV,surfaceDist,i,eps,O1O2,meanRad);
+ else if (viscousShear)
+ shearLubF = flow.computeViscousShearForce ( deltaShearV, i , Rh);
+
+ if (viscousShear || shearLubrication){
+
+ flow.shearLubricationForces[id1]+=shearLubF;
+ flow.shearLubricationForces[id2]+=(-shearLubF);
+ flow.shearLubricationTorques[id1]+=O1CVector.cross(shearLubF);
+ flow.shearLubricationTorques[id2]+=O2CVector.cross(-shearLubF);
+
+ /// Compute the pump Lubrication torque on each particle
+
+ if (pumpTorque){
+ pumpT = flow.computePumpTorque(deltaShearAngVel, surfaceDist, i, eps, meanRad );
+ flow.pumpLubricationTorques[id1]+=(-pumpT);
+ flow.pumpLubricationTorques[id2]+=pumpT;}
+
+ /// Compute the twist Lubrication torque on each particle
+
+ if (twistTorque){
+ deltaAngNormVel = (normal.dot(deltaAngVel))*normal ;
+ twistT = flow.computeTwistTorque(deltaAngNormVel, surfaceDist, i, eps, meanRad );
+ flow.twistLubricationTorques[id1]+=(-twistT);
+ flow.twistLubricationTorques[id2]+=twistT;
+ }
+ }
+ /// Compute the viscous shear stress on each particle
+
+ if (viscousShearBodyStress){
+ flow.shearLubricationBodyStress[id1] += shearLubF * O1CVector.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
+ flow.shearLubricationBodyStress[id2] += (-shearLubF) * O2CVector.transpose()/ (4.0/3.0 *3.14* pow(r2,3));
+ flow.shearStressInteraction.push_back(shearLubF * O1O2Vector.transpose()/(4.0/3.0 *3.14* pow(r1,3)));
+ }
+
+ /// Compute the normal lubrication force applied on each particle
+
+ if (normalLubrication){
+ deltaNormV = normal.dot(deltaV);
+ flow.deltaNormVel.push_back(deltaNormV * normal);
+ normaLubF = flow.computeNormalLubricationForce (deltaNormV, surfaceDist, i,eps,stiffness,scene->dt,meanRad)*normal;
+ flow.normalLubricationForce[id1]+=normaLubF;
+ flow.normalLubricationForce[id2]+=(-normaLubF);
+
+ /// Compute the normal lubrication stress on each particle
+
+ if (viscousNormalBodyStress){
+ flow.normalLubricationBodyStress[id1] += normaLubF * O1CVector.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
+ flow.normalLubricationBodyStress[id2] += (-normaLubF) *O2CVector.transpose() / (4.0/3.0 *3.14* pow(r2,3));
+ flow.normalStressInteraction.push_back(normaLubF * O1O2Vector.transpose()/(4.0/3.0 *3.14* pow(r1,3)));
+ }
+ }
+
+ if (!hasFictious)
+ flow.onlySpheresInteractions.push_back(i);
+
+ }
+ }
+}
+
+#endif //FLOW_ENGINE
+
+#endif /* YADE_CGAL */
+
