yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10591
[Branch ~yade-pkg/yade/git-trunk] Rev 3845: PFV: First step in refactoring templates
------------------------------------------------------------
revno: 3845
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Fri 2014-03-21 19:47:45 +0100
message:
PFV: First step in refactoring templates
removed:
lib/triangulation/PeriodicFlow.cpp
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphere.ipp
lib/triangulation/FlowBoundingSphereLinSolv.hpp
lib/triangulation/FlowBoundingSphereLinSolv.ipp
lib/triangulation/KinematicLocalisationAnalyser.cpp
lib/triangulation/KinematicLocalisationAnalyser.hpp
lib/triangulation/Network.hpp
lib/triangulation/Network.ipp
lib/triangulation/PeriodicFlow.hpp
lib/triangulation/PeriodicFlowLinSolv.hpp
lib/triangulation/PeriodicFlowLinSolv.ipp
lib/triangulation/RegularTriangulation.h
lib/triangulation/Tesselation.h
lib/triangulation/Tesselation.ipp
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.hpp
pkg/dem/MicroMacroAnalyser.cpp
pkg/dem/MicroMacroAnalyser.hpp
pkg/dem/TesselationWrapper.cpp
pkg/dem/TesselationWrapper.hpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'lib/triangulation/FlowBoundingSphere.cpp'
--- lib/triangulation/FlowBoundingSphere.cpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2014-03-21 18:47:45 +0000
@@ -4,25 +4,25 @@
* 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 FLOW_ENGINE
-
-#include "yade/lib/triangulation/FlowBoundingSphere.hpp"
-
-namespace CGT {
-CVector PeriodicCellInfo::gradP;
-CVector PeriodicCellInfo::hSize[3];
-CVector PeriodicCellInfo::deltaP;
-}
-
-
-//Forcing instanciation of the template to avoid linkage problems
-typedef CGT::FlowBoundingSphere<FlowTesselation> FlowBoundingSphere;
-FlowBoundingSphere ex;
-#ifdef LINSOLV
-typedef CGT::FlowBoundingSphereLinSolv<FlowBoundingSphere> FlowBoundingSphereLinSolv;
-FlowBoundingSphereLinSolv exls;
-#endif
-
-
-#endif //FLOW_ENGINE
+// #ifdef FLOW_ENGINE
+//
+// #include "yade/lib/triangulation/FlowBoundingSphere.hpp"
+//
+// namespace CGT {
+// CVector PeriodicCellInfo::gradP;
+// CVector PeriodicCellInfo::hSize[3];
+// CVector PeriodicCellInfo::deltaP;
+// }
+//
+//
+// //Forcing instanciation of the template to avoid linkage problems
+// typedef CGT::FlowBoundingSphere<FlowTesselation> FlowBoundingSphere;
+// FlowBoundingSphere ex;
+// #ifdef LINSOLV
+// typedef CGT::FlowBoundingSphereLinSolv<FlowBoundingSphere> FlowBoundingSphereLinSolv;
+// FlowBoundingSphereLinSolv exls;
+// #endif
+//
+//
+// #endif //FLOW_ENGINE
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2014-03-21 18:47:45 +0000
@@ -29,29 +29,29 @@
DECLARE_TESSELATION_TYPES(Network<Tesselation>)
//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::vTotalePorosity; 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;
+ 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;
virtual ~FlowBoundingSphere();
FlowBoundingSphere();
- bool slipOnLaterals;
+ bool slipBoundary;
double tolerance;
double relax;
double ks; //Hydraulic Conductivity
- bool clampKValues, meanKStat, distance_correction;
+ bool clampKValues, meanKStat, distanceCorrection;
bool OUTPUT_BOUDARIES_RADII;
bool noCache;//flag for checking if cached values cell->unitForceVectors have been defined
bool computedOnce;//flag for checking if current triangulation has been computed at least once
bool pressureChanged;//are imposed pressures modified (on python side)? When it happens, we have to reApplyBoundaryConditions
int errorCode;
- //Handling imposed pressures/fluxes on elements in the form of {point,value} pairs, ipCells contains the cell handles corresponding to point
+ //Handling imposed pressures/fluxes on elements in the form of {point,value} pairs, IPCells contains the cell handles corresponding to point
vector<pair<Point,Real> > imposedP;
- vector<CellHandle> ipCells;
+ vector<CellHandle> IPCells;
vector<pair<Point,Real> > imposedF;
- vector<CellHandle> ifCells;
+ vector<CellHandle> IFCells;
void initNewTri () {noCache=true; /*isLinearSystemSet=false; areCellsOrdered=false;*/}//set flags after retriangulation
bool permeabilityMap;
@@ -62,7 +62,7 @@
double maxKdivKmean;
int Iterations;
- bool RAVERAGE;
+ bool rAverage;
int walls_id[6];
#define parallel_forces
#ifdef parallel_forces
@@ -104,10 +104,10 @@
double fluidBulkModulus;
void displayStatistics();
- void initializePressures ( double pZero );
+ void initializePressure ( double pZero );
bool reApplyBoundaryConditions ();
void computeFacetForcesWithCache(bool onlyCache=false);
- void saveVtk (const char* folder );
+ void saveVtk (const char* folder);
#ifdef XVIEW
void dessineTriangulation ( Vue3D &Vue, RTriangulation &T );
void dessineShortTesselation ( Vue3D &Vue, Tesselation &Tes );
@@ -160,12 +160,12 @@
};
} //namespace CGT
-
+#include <yade/lib/triangulation/FlowBoundingSphere.ipp>
#ifdef LINSOLV
#include "yade/lib/triangulation/FlowBoundingSphereLinSolv.hpp"
#endif
/// _____ Template Implementation ____
-#include "yade/lib/triangulation/FlowBoundingSphere.ipp"
+// #include "yade/lib/triangulation/FlowBoundingSphereLinSolv.ipp"
#endif //FLOW_ENGINE
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2014-03-21 18:47:45 +0000
@@ -65,20 +65,20 @@
tessBasedForce = true;
for (int i=0;i<6;i++) boundsIds[i] = 0;
minPermLength=-1;
- slipOnLaterals = false;//no-slip/symmetry conditions on lateral boundaries
+ slipBoundary = false;//no-slip/symmetry conditions on lateral boundaries
tolerance = 1e-07;
relax = 1.9;
ks=0;
- distance_correction = true;
+ distanceCorrection = true;
clampKValues = true;
meanKStat = true; KOptFactor=0;
noCache=true;
pressureChanged=false;
- computeAllCells=true;//might be turned false IF the code is reorganized (we can make a separate function to compute unitForceVectors outside Compute_Permeability) AND it really matters for CPU time
+ computeAllCells=true;//might be turned false IF the code is reorganized (we can make a separate function to compute unitForceVectors outside compute_Permeability) AND it really matters for CPU time
debugOut = true;
- RAVERAGE = false; /** use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
+ rAverage = false; /** use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
OUTPUT_BOUDARIES_RADII = false;
- RAVERAGE = false; /** if true use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
+ rAverage = false; /** if true use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
// areaR2Permeability=true;
permeabilityMap = false;
computedOnce=false;
@@ -447,7 +447,7 @@
{
if (debugOut) cout << "----Computing_Permeability------" << endl;
RTriangulation& Tri = T[currentTes].Triangulation();
- VSolidTot = 0, Vtotalissimo = 0, vPoral = 0, sSolidTot = 0, vTotalePorosity=0, vPoralPorosity=0;
+ VSolidTot = 0, Vtotalissimo = 0, vPoral = 0, sSolidTot = 0, vTotalPorosity=0, vPoralPorosity=0;
FiniteCellsIterator cellEnd = Tri.finite_cells_end();
CellHandle neighbourCell;
@@ -467,7 +467,7 @@
neighbourCell = cell->neighbor(j);
Point& p2 = neighbourCell->info();
if (!Tri.is_infinite(neighbourCell) && (neighbourCell->info().isvisited==ref || computeAllCells)) {
- //Compute and store the area of sphere-facet intersections for later use
+ //compute and store the area of sphere-facet intersections for later use
VertexHandle W [3];
for (int kk=0; kk<3; kk++) {
W[kk] = cell->vertex(facetVertices[j][kk]);
@@ -484,7 +484,7 @@
pass+=1;
CVector l = p1 - p2;
distance = sqrt(l.squared_length());
- if (!RAVERAGE) radius = 2* computeHydraulicRadius(cell, j);
+ if (!rAverage) radius = 2* computeHydraulicRadius(cell, j);
else radius = (computeEffectiveRadius(cell, j)+computeEquivalentRadius(cell,j))*0.5;
if (radius<0) NEG++;
else POS++;
@@ -493,7 +493,7 @@
}
Real fluidArea=0;
if (distance!=0) {
- if (minPermLength>0 && distance_correction) distance=max(minPermLength,distance);
+ if (minPermLength>0 && distanceCorrection) distance=max(minPermLength,distance);
const CVector& Surfk = cell->info().facetSurfaces[j];
Real area = sqrt(Surfk.squared_length());
const CVector& crossSections = cell->info().facetSphereCrossSections[j];
@@ -596,9 +596,9 @@
Vgrains += 1.33333333 * M_PI * pow(vIt->point().weight(),1.5);}}
cout<<grains<<"grains - " <<"vTotal = " << vTotal << " Vgrains = " << Vgrains << " vPoral1 = " << (vTotal-Vgrains) << endl;
cout << "Vtotalissimo = " << Vtotalissimo/2 << " VSolidTot = " << VSolidTot/2 << " vPoral2 = " << vPoral/2 << " sSolidTot = " << sSolidTot << endl<< endl;
- if (!RAVERAGE) cout << "------Hydraulic Radius is used for permeability computation------" << endl << endl;
+ if (!rAverage) cout << "------Hydraulic Radius is used for permeability computation------" << endl << endl;
else cout << "------Average Radius is used for permeability computation------" << endl << endl;
- cout << "-----Computed_Permeability-----" << endl;}
+ cout << "-----computed_Permeability-----" << endl;}
}
template <class Tesselation>
@@ -690,17 +690,17 @@
RTriangulation& Tri = T[currentTes].Triangulation();
if (Tri.is_infinite(cell->neighbor(j))) return 0;
double Vpore = this->volumePoreVoronoiFraction(cell, j);
- double Ssolid = this->surfaceSolidPore(cell, j, slipOnLaterals, /*reuse the same facet data*/ true);
+ double Ssolid = this->surfaceSolidPore(cell, j, slipBoundary, /*reuse the same facet data*/ true);
//handle symmetry (tested ok)
- if (slipOnLaterals && facetNFictious>0) {
+ if (slipBoundary && facetNFictious>0) {
//! Include a multiplier so that permeability will be K/2 or K/4 in symmetry conditions
Real mult= facetNFictious==1 ? multSym1 : multSym2;
return Vpore/Ssolid*mult;}
return Vpore/Ssolid;
}
template <class Tesselation>
-void FlowBoundingSphere<Tesselation>::initializePressures( double pZero )
+void FlowBoundingSphere<Tesselation>::initializePressure( double pZero )
{
RTriangulation& Tri = T[currentTes].Triangulation();
FiniteCellsIterator cellEnd = Tri.finite_cells_end();
@@ -724,26 +724,26 @@
}
}
}
- ipCells.clear();
+ IPCells.clear();
for (unsigned int n=0; n<imposedP.size();n++) {
CellHandle cell=Tri.locate(imposedP[n].first);
//check redundancy
- for (unsigned int kk=0;kk<ipCells.size();kk++){
- if (cell==ipCells[kk]) cerr<<"Two imposed pressures fall in the same cell."<<endl;
+ for (unsigned int kk=0;kk<IPCells.size();kk++){
+ if (cell==IPCells[kk]) cerr<<"Two imposed pressures fall in the same cell."<<endl;
else if (cell->info().Pcondition) cerr<<"Imposed pressure fall in a boundary condition."<<endl;}
- ipCells.push_back(cell);
+ IPCells.push_back(cell);
cell->info().p()=imposedP[n].second;
cell->info().Pcondition=true;}
pressureChanged=false;
- ifCells.clear();
+ IFCells.clear();
for (unsigned int n=0; n<imposedF.size();n++) {
CellHandle cell=Tri.locate(imposedF[n].first);
//check redundancy
- for (unsigned int kk=0;kk<ipCells.size();kk++){
- if (cell==ipCells[kk]) cerr<<"Both flux and pressure are imposed in the same cell."<<endl;
+ for (unsigned int kk=0;kk<IPCells.size();kk++){
+ if (cell==IPCells[kk]) cerr<<"Both flux and pressure are imposed in the same cell."<<endl;
else if (cell->info().Pcondition) cerr<<"Imposed flux fall in a pressure boundary condition."<<endl;}
- ifCells.push_back(cell);
+ IFCells.push_back(cell);
cell->info().Pcondition=false;}
}
@@ -763,8 +763,8 @@
}
}
for (unsigned int n=0; n<imposedP.size();n++) {
- ipCells[n]->info().p()=imposedP[n].second;
- ipCells[n]->info().Pcondition=true;}
+ IPCells[n]->info().p()=imposedP[n].second;
+ IPCells[n]->info().Pcondition=true;}
pressureChanged=false;
return true;
}
@@ -907,8 +907,9 @@
{
const CellHandle& cell = *it;
if (cell->info().isGhost) continue;
+ Q1 -= cell->info().dv();
for (int j2=0; j2<4; j2++)
- Q1 += (cell->neighbor(j2)->info().kNorm())[Tri.mirror_index(cell, j2)]* (cell->neighbor(j2)->info().p()-cell->info().p());
+ Q1 += (cell->info().kNorm())[j2]* (cell->neighbor(j2)->info().shiftedP()-cell->info().shiftedP());
}
return Q1;
}
@@ -1113,7 +1114,7 @@
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::dessineShortTesselation(Vue3D &Vue, Tesselation &Tes)
{
- if (!Tes.Computed()) Tes.Compute();
+ if (!Tes.computed()) Tes.compute();
double* Segments = NULL;
long N_seg = Tes.newListeShortEdges(&Segments);
Vue.Dessine_Segment(Segments, N_seg);
@@ -1158,7 +1159,7 @@
boundary(yMinId).value=0;
boundary(yMaxId).value=1;
double pZero = abs((boundary(yMinId).value-boundary(yMaxId).value)/2);
- initializePressures( pZero );
+ initializePressure( pZero );
gaussSeidel();
const char *kk = "Permeability";
return permeameter(boundary(yMinId).value, boundary(yMaxId).value, Section, DeltaY, kk);
@@ -1214,7 +1215,7 @@
{
int hasFictious= (ed_it->first)->vertex(ed_it->second)->info().isFictious + (ed_it->first)->vertex(ed_it->third)->info().isFictious;
if (hasFictious==2) continue;
- double area = T[currentTes].ComputeVFacetArea(ed_it);
+ double area = T[currentTes].computeVFacetArea(ed_it);
edgeSurfaces.push_back(area);
unsigned int id1 = (ed_it->first)->vertex(ed_it->second)->info().id();
unsigned int id2 = (ed_it->first)->vertex(ed_it->third)->info().id();
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.hpp'
--- lib/triangulation/FlowBoundingSphereLinSolv.hpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.hpp 2014-03-21 18:47:45 +0000
@@ -34,7 +34,7 @@
namespace CGT {
-template<class FlowType>
+template<class _Tesselation, class FlowType=FlowBoundingSphere<_Tesselation> >
class FlowBoundingSphereLinSolv : public FlowType
{
public:
@@ -53,6 +53,7 @@
using FlowType::reApplyBoundaryConditions;
using FlowType::pressureChanged;
using FlowType::computedOnce;
+ using FlowType::resetNetwork;
//! TAUCS DECs
vector<FiniteCellsIterator> orderedCells;
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.ipp'
--- lib/triangulation/FlowBoundingSphereLinSolv.ipp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.ipp 2014-03-21 18:47:45 +0000
@@ -60,15 +60,15 @@
#ifdef XVIEW
Vue3D Vue1;
#endif
-template<class FlowType>
-FlowBoundingSphereLinSolv<FlowType>::~FlowBoundingSphereLinSolv()
+template<class _Tesselation, class FlowType>
+FlowBoundingSphereLinSolv<_Tesselation,FlowType>::~FlowBoundingSphereLinSolv()
{
#ifdef TAUCS_LIB
if (Fccs) taucs_ccs_free(Fccs);
#endif
}
-template<class FlowType>
-FlowBoundingSphereLinSolv<FlowType>::FlowBoundingSphereLinSolv(): FlowType() {
+template<class _Tesselation, class FlowType>
+FlowBoundingSphereLinSolv<_Tesselation,FlowType>::FlowBoundingSphereLinSolv(): FlowType() {
useSolver=0;
isLinearSystemSet=0;
isFullLinearSystemGSSet=0;
@@ -90,10 +90,10 @@
}
-template<class FlowType>
-void FlowBoundingSphereLinSolv<FlowType>::swapFwd (double* v, int i) {double temp = v[i]; v[i]=v[i+1]; v[i+1]=temp;}
-template<class FlowType>
-void FlowBoundingSphereLinSolv<FlowType>::swapFwd (int* v, int i) {double temp = v[i]; v[i]=v[i+1]; v[i+1]=temp;}
+template<class _Tesselation, class FlowType>
+void FlowBoundingSphereLinSolv<_Tesselation,FlowType>::swapFwd (double* v, int i) {double temp = v[i]; v[i]=v[i+1]; v[i+1]=temp;}
+template<class _Tesselation, class FlowType>
+void FlowBoundingSphereLinSolv<_Tesselation,FlowType>::swapFwd (int* v, int i) {double temp = v[i]; v[i]=v[i+1]; v[i+1]=temp;}
//spatial sort traits to use with a pair of CGAL::sphere pointers and integer.
//template<class _Triangulation>
@@ -117,8 +117,8 @@
Less_z_3 less_z_3_object() const {return Less_z_3();}
};
-template<class FlowType>
-void FlowBoundingSphereLinSolv<FlowType>::resetNetwork() {
+template<class _Tesselation, class FlowType>
+void FlowBoundingSphereLinSolv<_Tesselation,FlowType>::resetNetwork() {
FlowType::resetNetwork();
isLinearSystemSet=false;
isFullLinearSystemGSSet=false;
@@ -140,8 +140,8 @@
}
#endif
}
-template<class FlowType>
-int FlowBoundingSphereLinSolv<FlowType>::setLinearSystem(Real dt)
+template<class _Tesselation, class FlowType>
+int FlowBoundingSphereLinSolv<_Tesselation,FlowType>::setLinearSystem(Real dt)
{
RTriangulation& Tri = T[currentTes].Triangulation();
@@ -285,11 +285,11 @@
return ncols;
}
-template<class FlowType>
-void FlowBoundingSphereLinSolv<FlowType>::copyGsToCells() {for (int ii=1; ii<=ncols; ii++) T_cells[ii]->info().p()=gsP[ii];}
+template<class _Tesselation, class FlowType>
+void FlowBoundingSphereLinSolv<_Tesselation,FlowType>::copyGsToCells() {for (int ii=1; ii<=ncols; ii++) T_cells[ii]->info().p()=gsP[ii];}
-template<class FlowType>
-void FlowBoundingSphereLinSolv<FlowType>::copyCellsToGs (Real dt)
+template<class _Tesselation, class FlowType>
+void FlowBoundingSphereLinSolv<_Tesselation,FlowType>::copyCellsToGs (Real dt)
{
for (int ii=1; ii<=ncols; ii++){
gsP[ii]=T_cells[ii]->info().p();
@@ -298,11 +298,11 @@
}
}
-template<class FlowType>
-void FlowBoundingSphereLinSolv<FlowType>::copyLinToCells() {for (int ii=1; ii<=ncols; ii++) {T_cells[ii]->info().p()=T_x[ii-1];} }
+template<class _Tesselation, class FlowType>
+void FlowBoundingSphereLinSolv<_Tesselation,FlowType>::copyLinToCells() {for (int ii=1; ii<=ncols; ii++) {T_cells[ii]->info().p()=T_x[ii-1];} }
-template<class FlowType>
-void FlowBoundingSphereLinSolv<FlowType>::copyCellsToLin (Real dt)
+template<class _Tesselation, class FlowType>
+void FlowBoundingSphereLinSolv<_Tesselation,FlowType>::copyCellsToLin (Real dt)
{
for (int ii=1; ii<=ncols; ii++) {
T_bv[ii-1]=T_b[ii-1]-T_cells[ii]->info().dv();
@@ -310,9 +310,9 @@
}
/// For Gauss Seidel, we need the full matrix
-template<class FlowType>
-// int FlowBoundingSphereLinSolv<FlowType>::SetLinearSystemFullGS()
-int FlowBoundingSphereLinSolv<FlowType>::setLinearSystemFullGS(Real dt)
+template<class _Tesselation, class FlowType>
+// int FlowBoundingSphereLinSolv<_Tesselation,FlowType>::SetLinearSystemFullGS()
+int FlowBoundingSphereLinSolv<_Tesselation,FlowType>::setLinearSystemFullGS(Real dt)
{
//WARNING : boundary conditions (Pcondition, p values) must have been set for a correct definition
RTriangulation& Tri = T[currentTes].Triangulation();
@@ -333,7 +333,7 @@
spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort<RTriangulation>());
// double pZero=0;
-// if (y_min_id>=0 and y_max_id>y_min_id) pZero = abs((boundary(y_min_id).value-boundary(y_max_id).value)/2);
+// if (yMinId>=0 and yMaxId>yMinId) pZero = abs((boundary(yMinId).value-boundary(yMaxId).value)/2);
gsP.resize(ncols+1);
// _gsP.resize(ncols+1);
gsB.resize(ncols+1);
@@ -415,8 +415,8 @@
return ncols;
}
-template<class FlowType>
-void FlowBoundingSphereLinSolv<FlowType>::vectorizedGaussSeidel(Real dt)
+template<class _Tesselation, class FlowType>
+void FlowBoundingSphereLinSolv<_Tesselation,FlowType>::vectorizedGaussSeidel(Real dt)
{
// cout<<"VectorizedGaussSeidel"<<endl;
if (!isFullLinearSystemGSSet || (isFullLinearSystemGSSet && reApplyBoundaryConditions())) setLinearSystemFullGS(dt);
@@ -498,8 +498,8 @@
if (debugOut) cerr <<"GS iterations : "<<j-1<<endl;
}
-template<class FlowType>
-void FlowBoundingSphereLinSolv<FlowType>::sortV(int k1, int k2, int* is, double* ds){
+template<class _Tesselation, class FlowType>
+void FlowBoundingSphereLinSolv<_Tesselation,FlowType>::sortV(int k1, int k2, int* is, double* ds){
for (int k=k1; k<k2; k++) {
int kk=k;
while (kk>=k1 && is[kk]>is[kk+1]) {
@@ -509,8 +509,8 @@
}
}
-template<class FlowType>
-int FlowBoundingSphereLinSolv<FlowType>::eigenSolve(Real dt)
+template<class _Tesselation, class FlowType>
+int FlowBoundingSphereLinSolv<_Tesselation,FlowType>::eigenSolve(Real dt)
{
#ifdef EIGENSPARSE_LIB
if (!isLinearSystemSet || (isLinearSystemSet && reApplyBoundaryConditions())) ncols = setLinearSystem(dt);
@@ -541,8 +541,8 @@
}
-template<class FlowType>
-int FlowBoundingSphereLinSolv<FlowType>::taucsSolve(Real dt)
+template<class _Tesselation, class FlowType>
+int FlowBoundingSphereLinSolv<_Tesselation,FlowType>::taucsSolve(Real dt)
{
#ifdef TAUCS_LIB
if (debugOut) cerr <<endl<<"TAUCS solve"<<endl;
@@ -597,8 +597,8 @@
#endif
return 0;
}
-template<class FlowType>
-int FlowBoundingSphereLinSolv<FlowType>::pardisoSolve(Real dt)
+template<class _Tesselation, class FlowType>
+int FlowBoundingSphereLinSolv<_Tesselation,FlowType>::pardisoSolve(Real dt)
{
cerr <<endl<<"PardisoSolve solve"<<endl;
#ifndef PARDISO
@@ -738,8 +738,8 @@
}
-template<class FlowType>
-int FlowBoundingSphereLinSolv<FlowType>::pardisoSolveTest()
+template<class _Tesselation, class FlowType>
+int FlowBoundingSphereLinSolv<_Tesselation,FlowType>::pardisoSolveTest()
{
#ifndef PARDISO
return 0;
@@ -886,8 +886,8 @@
return 0;
#endif
}
-template<class FlowType>
-int FlowBoundingSphereLinSolv<FlowType>::taucsSolveTest()
+template<class _Tesselation, class FlowType>
+int FlowBoundingSphereLinSolv<_Tesselation,FlowType>::taucsSolveTest()
{
#ifdef TAUCS_LIB
cout <<endl<<"TAUCS solve test"<<endl;
=== modified file 'lib/triangulation/KinematicLocalisationAnalyser.cpp'
--- lib/triangulation/KinematicLocalisationAnalyser.cpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/KinematicLocalisationAnalyser.cpp 2014-03-21 18:47:45 +0000
@@ -75,7 +75,7 @@
Delta_epsilon(2,2) = TS1->eps2 - TS0->eps2;
}
-const vector<Tenseur3>& KinematicLocalisationAnalyser::ComputeParticlesDeformation(const char* state_file1, const char* state_file0, bool usebz2)
+const vector<Tenseur3>& KinematicLocalisationAnalyser::computeParticlesDeformation(const char* state_file1, const char* state_file0, bool usebz2)
{
consecutive = false;
bz2 = usebz2;
@@ -85,7 +85,7 @@
Delta_epsilon(3,3) = TS1->eps3 - TS0->eps3;
Delta_epsilon(1,1) = TS1->eps1 - TS0->eps1;
Delta_epsilon(2,2) = TS1->eps2 - TS0->eps2;
- return ComputeParticlesDeformation();
+ return computeParticlesDeformation();
}
@@ -193,7 +193,7 @@
bool KinematicLocalisationAnalyser::DefToFile(const char* output_file_name)
{
- ComputeParticlesDeformation();
+ computeParticlesDeformation();
Tesselation& Tes = TS1->tesselation();
RTriangulation& Tri = Tes.Triangulation();
basicVTKwritter vtk(n_real_vertices, n_real_cells);
@@ -762,7 +762,7 @@
if (vol_divide) T/= Tesselation::Volume(cell);
}
-const vector<Tenseur3>& KinematicLocalisationAnalyser::ComputeParticlesDeformation(void)
+const vector<Tenseur3>& KinematicLocalisationAnalyser::computeParticlesDeformation(void)
{
Tesselation& Tes = TS1->tesselation();
RTriangulation& Tri = Tes.Triangulation();
@@ -777,8 +777,8 @@
Delta_epsilon(1,1) = TS1->eps1 - TS0->eps1;
Delta_epsilon(2,2) = TS1->eps2 - TS0->eps2;
- //Compute Voronoi tesselation (i.e. voronoi center of each cell)
- if (!Tes.Computed()) Tes.Compute();
+ //compute Voronoi tesselation (i.e. voronoi center of each cell)
+ if (!Tes.computed) Tes.compute();
if (ParticleDeformation.size() != (unsigned int)(Tes.Max_id() + 1)) {
ParticleDeformation.clear();
ParticleDeformation.resize(Tes.Max_id() + 1);
@@ -794,7 +794,7 @@
Finite_cells_iterator cell = Tri.finite_cells_begin();
Finite_cells_iterator cell0 = Tri.finite_cells_end();
- //Compute grad_u and volumes of all cells in the triangulation, and assign them to each of the vertices ( volume*grad_u is added here rather than grad_u, the weighted average is computed later )
+ //compute grad_u and volumes of all cells in the triangulation, and assign them to each of the vertices ( volume*grad_u is added here rather than grad_u, the weighted average is computed later )
//define the number of non-fictious cells, i.e. not in contact with a boundary
n_real_cells=0;
for (; cell != cell0; cell++) { // calcule la norme du d�viateur dans chaque cellule
@@ -839,7 +839,7 @@
return ParticleDeformation;
}
-Real KinematicLocalisationAnalyser::ComputeMacroPorosity(void)
+Real KinematicLocalisationAnalyser::computeMacroPorosity(void)
{
return (1-v_solid_total/(TS1->haut*TS1->larg*TS1->prof));
}
=== modified file 'lib/triangulation/KinematicLocalisationAnalyser.hpp'
--- lib/triangulation/KinematicLocalisationAnalyser.hpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/KinematicLocalisationAnalyser.hpp 2014-03-21 18:47:45 +0000
@@ -6,7 +6,7 @@
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
-/*! \brief Computes statistics of micro-variables assuming axi-symetry.
+/*! \brief computes statistics of micro-variables assuming axi-symetry.
*/
@@ -53,7 +53,7 @@
bool DistribsToFile (const char* output_file_name);
- ///Write the averaged deformation on each grain in a file (vertices and cells lists included in the file), no need to call ComputeParticlesDeformation()
+ ///Write the averaged deformation on each grain in a file (vertices and cells lists included in the file), no need to call computeParticlesDeformation()
bool DefToFile (const char* output_file_name = "deformations");
bool DefToFile (const char* state_file1, const char* state_file0, const char* output_file_name="deformation.vtk", bool usebz2=false);
///Save/Load states using bz2 compression
@@ -79,14 +79,14 @@
///Add surface*displacement to T
void Grad_u ( Finite_cells_iterator cell, int facet, CVector &V, Tenseur3& T );
- ///Compute grad_u in cell (by default, T= average grad_u in cell, if !vol_divide, T=grad_u*volume
+ ///compute grad_u in cell (by default, T= average grad_u in cell, if !vol_divide, T=grad_u*volume
void Grad_u ( Finite_cells_iterator cell, Tenseur3& T, bool vol_divide=true );
- /// Compute grad_u for all particles, by summing grad_u of all adjaent cells using current states
- const vector<Tenseur3>& ComputeParticlesDeformation (void);
+ /// compute grad_u for all particles, by summing grad_u of all adjaent cells using current states
+ const vector<Tenseur3>& computeParticlesDeformation (void);
/// Do everything in one step by giving some final (file1) and initial (file0) positions
- const vector<Tenseur3>& ComputeParticlesDeformation(const char* state_file1, const char* state_file0, bool usebz2 = true);
- ///Compute porisity from cumulated spheres volumes and positions of boxes
- Real ComputeMacroPorosity (void );
+ const vector<Tenseur3>& computeParticlesDeformation(const char* state_file1, const char* state_file0, bool usebz2 = true);
+ ///compute porisity from cumulated spheres volumes and positions of boxes
+ Real computeMacroPorosity (void );
CVector Deplacement ( Cell_handle cell ); //donne le d�placement d'un sommet de voronoi
CVector Deplacement ( Finite_cells_iterator cell, int facet ); //mean displacement on a facet
=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/Network.hpp 2014-03-21 18:47:45 +0000
@@ -54,7 +54,7 @@
vector<CellHandle> boundingCells [6];
Point cornerMin;
Point cornerMax;
- Real VSolidTot, Vtotalissimo, vPoral, sSolidTot, vPoralPorosity, vTotalePorosity;
+ Real VSolidTot, Vtotalissimo, vPoral, sSolidTot, vPoralPorosity, vTotalPorosity;
Boundary boundaries [6];
Boundary& boundary (int b) {return boundaries[b-idOffset];}
short idOffset;
@@ -76,7 +76,7 @@
double volumeDoubleFictiousPore(VertexHandle SV1, VertexHandle SV2, VertexHandle SV3, Point PV1);
double volumeSingleFictiousPore(VertexHandle SV1, VertexHandle SV2, VertexHandle SV3, Point PV1);
double volumePoreVoronoiFraction ( CellHandle& cell, int& j, bool reuseFacetData=false);
- double surfaceSolidPore( CellHandle cell, int j, bool slipOnLaterals, bool reuseFacetData=false);
+ double surfaceSolidPore( CellHandle cell, int j, bool slipBoundary, bool reuseFacetData=false);
double sphericalTriangleArea ( Sphere STA1, Sphere STA2, Sphere STA3, Point PTA1 );
CVector surfaceDoubleFictiousFacet(VertexHandle fSV1, VertexHandle fSV2, VertexHandle SV3);
=== modified file 'lib/triangulation/Network.ipp'
--- lib/triangulation/Network.ipp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/Network.ipp 2014-03-21 18:47:45 +0000
@@ -84,7 +84,7 @@
for (int i=0;i<4;i++){
if (cell->neighbor(i)->info().fictious()!=0) border=true;}
if (!border) {vPoralPorosity += Vtot - (Vsolid1 + Vsolid2);
- vTotalePorosity += Vtot;}
+ vTotalPorosity += Vtot;}
/**Vpore**/ return Vtot - (Vsolid1 + Vsolid2);
}; break;
@@ -251,7 +251,7 @@
}
template<class Tesselation>
-double Network<Tesselation>::surfaceSolidPore(CellHandle cell, int j, bool slipOnLaterals, bool reuseFacetData)
+double Network<Tesselation>::surfaceSolidPore(CellHandle cell, int j, bool slipBoundary, bool reuseFacetData)
{
if (!reuseFacetData) facetNFictious=detectFacetFictiousVertices(cell,j);
Point& p1 = cell->info();
@@ -291,7 +291,7 @@
Boundary &bi1 = boundary(SV1->info().id());
Ssolid1 = 0;
- if (bi1.flowCondition && ! slipOnLaterals) {
+ if (bi1.flowCondition && ! slipBoundary) {
Ssolid1 = abs(0.5*CGAL::cross_product(p1-p2, SV2->point()-SV3->point())[bi1.coordinate]);
cell->info().solidSurfaces[j][facetF1]=Ssolid1;
}
@@ -332,13 +332,13 @@
cell->info().solidSurfaces[j][facetRe1]=Ssolid1+Ssolid1n;
//area vector of triangle (p1,sphere,p2)
CVector p1p2v1Surface = 0.5*CGAL::cross_product(p1-p2,SV3->point()-p2);
- if (bi1.flowCondition && ! slipOnLaterals) {
+ if (bi1.flowCondition && ! slipBoundary) {
//projection on boundary 1
Ssolid2 = abs(p1p2v1Surface[bi1.coordinate]);
cell->info().solidSurfaces[j][facetF1]=Ssolid2;
} else cell->info().solidSurfaces[j][facetF1]=0;
- if (bi2.flowCondition && ! slipOnLaterals) {
+ if (bi2.flowCondition && ! slipBoundary) {
//projection on boundary 2
Ssolid3 = abs(p1p2v1Surface[bi2.coordinate]);
cell->info().solidSurfaces[j][facetF2]=Ssolid3;
=== removed file 'lib/triangulation/PeriodicFlow.cpp'
--- lib/triangulation/PeriodicFlow.cpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/PeriodicFlow.cpp 1970-01-01 00:00:00 +0000
@@ -1,518 +0,0 @@
-#ifdef FLOW_ENGINE
-
-#include"PeriodicFlow.hpp"
-
-#ifdef YADE_OPENMP
-// #define GS_OPEN_MP //It should never be defined if Yade is not using openmp
-#endif
-
-
-namespace CGT {
-
-void PeriodicFlow::interpolate(Tesselation& Tes, Tesselation& NewTes)
-{
- CellHandle oldCell;
- RTriangulation& Tri = Tes.Triangulation();
- for (VectorCell::iterator cellIt=NewTes.cellHandles.begin(); cellIt!=NewTes.cellHandles.end(); cellIt++){
- CellHandle& newCell = *cellIt;
- if (newCell->info().Pcondition || newCell->info().isGhost) continue;
- CVector center ( 0,0,0 );
- if (newCell->info().fictious()==0) for ( int k=0;k<4;k++ ) center= center + 0.25* (Tes.vertex(newCell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
- else {
- Real boundPos=0; int coord=0;
- for ( int k=0;k<4;k++ ) {
- if (!newCell->vertex (k)->info().isFictious) center= center+0.3333333333*(Tes.vertex(newCell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
- else {
- coord=boundary (newCell->vertex(k)->info().id()).coordinate;
- boundPos=boundary (newCell->vertex(k)->info().id()).p[coord];
- }
- }
- center=CVector(coord==0?boundPos:center[0],coord==1?boundPos:center[1],coord==2?boundPos:center[2]);
- }
- oldCell = Tri.locate(Point(center[0],center[1],center[2]));
- newCell->info().p() = oldCell->info().shiftedP();
- }
-// Tes.Clear();//Don't reset to avoid segfault when getting pressure in scripts just after interpolation
-}
-
-
-
-void PeriodicFlow::computeFacetForcesWithCache(bool onlyCache)
-{
- RTriangulation& Tri = T[currentTes].Triangulation();
- CVector nullVect(0,0,0);
- static vector<CVector> oldForces;
- if (oldForces.size()<=Tri.number_of_vertices()) oldForces.resize(Tri.number_of_vertices()+1);
- //reset forces
- for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
- if (noCache) {oldForces[v->info().id()]=nullVect; v->info().forces=nullVect;}
- else {oldForces[v->info().id()]=v->info().forces; v->info().forces=nullVect;}
- }
- CellHandle neighbourCell;
- VertexHandle mirrorVertex;
- CVector tempVect;
-
- //FIXME : Ema, be carefull with this (noCache), it needs to be turned true after retriangulation
- if (noCache) {
- for (VectorCell::iterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
- CellHandle& cell = *cellIt;
- for (int k=0;k<4;k++) cell->info().unitForceVectors[k]=nullVect;
- for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j))) {
-// #ifdef EIGENSPARSE_LIB
-// if (!cell->info().Pcondition) ++nPCells;
-// #endif
- neighbourCell = cell->neighbor(j);
- const CVector& Surfk = cell->info().facetSurfaces[j];
- //FIXME : later compute that fluidSurf only once in hydraulicRadius, for now keep full surface not modified in cell->info for comparison with other forces schemes
- //The ratio void surface / facet surface
- Real area = sqrt(Surfk.squared_length());
- CVector facetNormal = Surfk/area;
- const std::vector<CVector>& crossSections = cell->info().facetSphereCrossSections;
- CVector fluidSurfk = cell->info().facetSurfaces[j]*cell->info().facetFluidSurfacesRatio[j];
- /// handle fictious vertex since we can get the projected surface easily here
- if (cell->vertex(j)->info().isFictious) {
- Real projSurf=abs(Surfk[boundary(cell->vertex(j)->info().id()).coordinate]);
- tempVect=-projSurf*boundary(cell->vertex(j)->info().id()).normal;
- //define the cached value for later use with cache*p
- cell->info().unitForceVectors[j]=cell->info().unitForceVectors[j]+ tempVect;
- }
- /// Apply weighted forces f_k=sqRad_k/sumSqRad*f
- CVector facetUnitForce = -fluidSurfk*cell->info().solidSurfaces[j][3];
- for (int y=0; y<3;y++) {
- //add to cached value
- cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]+facetUnitForce*cell->info().solidSurfaces[j][y];
- //uncomment to get total force / comment to get only viscous forces (Bruno)
- if (!cell->vertex(facetVertices[j][y])->info().isFictious) {
- //add to cached value
- cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]-facetNormal*crossSections[j][y];
- }
- }
- }
- }
-
- noCache=false;//cache should always be defined after execution of this function
- if (onlyCache) return;
- }// end if(noCache)
-
- //use cached values
- //First define products that will be used below for all cells:
- Real pDeltas [3];
- for (unsigned int k=0; k<3;k++) pDeltas[k]=PeriodicCellInfo::hSize[k]*PeriodicCellInfo::gradP;
- //Then compute the forces
- for (VectorCell::iterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
- const CellHandle& cell = *cellIt;
- for (int yy=0;yy<4;yy++) {
- VertexInfo& vhi = cell->vertex(yy)->info();
- Real unshiftedP = cell->info().p();
- //the pressure translated to a ghost cell adjacent to the non-ghost vertex
- unshiftedP -= pDeltas[0]*vhi.period[0] + pDeltas[1]*vhi.period[1] +pDeltas[2]*vhi.period[2];
- T[currentTes].vertexHandles[vhi.id()]->info().forces=T[currentTes].vertexHandles[vhi.id()]->info().forces + cell->info().unitForceVectors[yy]*unshiftedP;
- }
- }
- if (debugOut) {
- CVector totalForce = nullVect;
- for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); v++)
- {
- if (!v->info().isFictious /*&& !v->info().isGhost*/ ){
- totalForce = totalForce + v->info().forces;}
- else /*if (!v->info().isGhost)*/{
- if (boundary(v->info().id()).flowCondition==1) totalForce = totalForce + v->info().forces;
- }
- }
- cout << "totalForce = "<< totalForce << endl;}
-}
-
-void PeriodicFlow::computePermeability()
-{
- if (debugOut) cout << "----Computing_Permeability (Periodic)------" << endl;
- RTriangulation& Tri = T[currentTes].Triangulation();
- VSolidTot = 0, Vtotalissimo = 0, vPoral = 0, sSolidTot = 0, vTotalePorosity=0, vPoralPorosity=0;
- CellHandle neighbourCell;
-
- double k=0, distance = 0, radius = 0;
- int surfneg=0; int NEG=0, POS=0, pass=0;
- Real meanK=0, STDEV=0, meanRadius=0, meanDistance=0;
- Real infiniteK=1e3;
- double volume_sub_pore = 0.f;
- VectorCell& cellHandles= T[currentTes].cellHandles;
- for (VectorCell::iterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
- CellHandle& cell = *cellIt;
- Point& p1 = cell->info();
- for (int j=0; j<4; j++){
- neighbourCell = cell->neighbor(j);
- Point& p2 = neighbourCell->info();
- if (!Tri.is_infinite(neighbourCell) /*&& (neighbour_cell->info().isvisited==ref || computeAllCells)*/) {
- //Compute and store the area of sphere-facet intersections for later use
- VertexHandle W [3];
- for (int kk=0; kk<3; kk++) {
- W[kk] = cell->vertex(facetVertices[j][kk]);
- }
- Sphere& v0 = W[0]->point();
- Sphere& v1 = W[1]->point();
- Sphere& v2 = W[2]->point();
-#ifdef USE_FAST_MATH
- //FIXME : code not compiling,, do the same as in "else"
- assert((W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point())>=0 && (W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point())<=1);
- for (int jj=0;jj<3;jj++)
- cell->info().facetSphereCrossSections[j][jj]=0.5*W[jj]->point().weight()*Wm3::FastInvCos1((W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point()));
-#else
- cell->info().facetSphereCrossSections[j]=CVector(
- W[0]->info().isFictious ? 0 : 0.5*v0.weight()*acos((v1-v0)*(v2-v0)/sqrt((v1-v0).squared_length()*(v2-v0).squared_length())),
- W[1]->info().isFictious ? 0 : 0.5*v1.weight()*acos((v0-v1)*(v2-v1)/sqrt((v1-v0).squared_length()*(v2-v1).squared_length())),
- W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0-v2)*(v1-v2)/sqrt((v1-v2).squared_length()*(v2-v0).squared_length())));
-#endif
- pass+=1;
- CVector l = p1 - p2;
- distance = sqrt(l.squared_length());
- if (!RAVERAGE) radius = 2* computeHydraulicRadius(cell, j);
- else radius = (computeEffectiveRadius(cell, j)+computeEquivalentRadius(cell,j))*0.5;
- if (radius<0) NEG++;
- else POS++;
- if (radius==0) {
- cout << "INS-INS PROBLEM!!!!!!!" << endl;
- }
- Real fluidArea=0;
- int test=0;
- if (distance!=0) {
- if (minPermLength>0 && distance_correction) distance=max(minPermLength,distance);
- const CVector& Surfk = cell->info().facetSurfaces[j];
- Real area = sqrt(Surfk.squared_length());
- const CVector& crossSections = cell->info().facetSphereCrossSections[j];
- Real S0=0;
- S0=checkSphereFacetOverlap(v0,v1,v2);
- if (S0==0) S0=checkSphereFacetOverlap(v1,v2,v0);
- if (S0==0) S0=checkSphereFacetOverlap(v2,v0,v1);
- //take absolute value, since in rare cases the surface can be negative (overlaping spheres)
- fluidArea=abs(area-crossSections[0]-crossSections[1]-crossSections[2]+S0);
- cell->info().facetFluidSurfacesRatio[j]=fluidArea/area;
- k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);
- meanDistance += distance;
- meanRadius += radius;
- meanK += k*kFactor;
-
- if (k<0 && debugOut) {surfneg+=1;
- cout<<"__ k<0 __"<<k<<" "<<" fluidArea "<<fluidArea<<" area "<<area<<" "<<crossSections[0]<<" "<<crossSections[1]<<" "<<crossSections[2] <<" "<<W[0]->info().id()<<" "<<W[1]->info().id()<<" "<<W[2]->info().id()<<" "<<p1<<" "<<p2<<" test "<<test<<endl;}
-
- } else cout <<"infinite K2! surfaces will be missing (FIXME)"<<endl; k = infiniteK;
- //Will be corrected in the next loop
- (cell->info().kNorm())[j]= k*kFactor;
- if (!neighbourCell->info().isGhost) (neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)]= (cell->info().kNorm())[j];
- //The following block is correct but not very usefull, since all values are clamped below with MIN and MAX, skip for now
-// else {//find the real neighbor connected to our cell through periodicity
-// CellHandle true_neighbour_cell = cellHandles[neighbour_cell->info().baseIndex];
-// for (int ii=0;ii<4;ii++)
-// if (true_neighbour_cell->neighbor(ii)->info().index == cell->info().index){
-// (true_neighbour_cell->info().kNorm())[ii]=(cell->info().kNorm())[j]; break;
-// }
-// }
- if(permeabilityMap){
- CellHandle c = cell;
- cell->info().s = cell->info().s + k*distance/fluidArea*volumePoreVoronoiFraction (c,j);
- volume_sub_pore += volumePoreVoronoiFraction (c,j);
- }
- }
- }
-// cell->info().isvisited = !ref;
- if(permeabilityMap){
- cell->info().s = cell->info().s/volume_sub_pore;
- volume_sub_pore = 0.f;
- }
-// }
- }
-
-
- if (debugOut) cout<<"surfneg est "<<surfneg<<endl;
- meanK /= pass;
- meanRadius /= pass;
- meanDistance /= pass;
- Real globalK=kFactor*meanDistance*vPoral/(sSolidTot*8.*viscosity);//An approximate value of macroscopic permeability, for clamping local values below
- if (debugOut) {
- cout << "PassCompK = " << pass << endl;
- cout << "meanK = " << meanK << endl;
- cout << "globalK = " << globalK << endl;
- cout << "maxKdivKmean*globalK = " << maxKdivKmean*globalK << endl;
- cout << "minKdivKmean*globalK = " << minKdivKmean*globalK << endl;
- cout << "meanTubesRadius = " << meanRadius << endl;
- cout << "meanDistance = " << meanDistance << endl;
- }
- pass=0;
- if (clampKValues) for (VectorCell::iterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
- CellHandle& cell = *cellIt;
- for (int j=0; j<4; j++) {
- neighbourCell = cell->neighbor(j);
- if (!Tri.is_infinite(neighbourCell) /*&& neighbour_cell->info().isvisited==ref*/) {
- pass++;
- (cell->info().kNorm())[j] = max(minKdivKmean*globalK, min((cell->info().kNorm())[j], maxKdivKmean*globalK));
- (neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)]=(cell->info().kNorm())[j];
- if (!neighbourCell->info().isGhost) (neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)]= (cell->info().kNorm())[j];
- else {//find the real neighbor connected to our cell through periodicity, as we want exactly the same permeability without rounding errors
- CellHandle& true_neighbourCell = cellHandles[neighbourCell->info().baseIndex];
- for (int ii=0;ii<4;ii++)
- if (true_neighbourCell->neighbor(ii)->info().index == cell->info().index){
- (true_neighbourCell->info().kNorm())[ii]=(cell->info().kNorm())[j]; break;
- }
- }
-
- }
- }
- }
- if (debugOut) cout << "PassKcorrect = " << pass << endl;
-
- if (debugOut) cout << "POS = " << POS << " NEG = " << NEG << " pass = " << pass << endl;
-
-// A loop to compute the standard deviation of the local K distribution, and use it to include/exclude K values higher then (meanK +/- K_opt_factor*STDEV)
- if (meanKStat)
- {
- std::ofstream k_opt_file("k_stdev.txt" ,std::ios::out);
- pass=0;
- FiniteCellsIterator cellEnd = Tri.finite_cells_end();
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
- for (int j=0; j<4; j++) {
- neighbourCell = cell->neighbor(j);
- if (!Tri.is_infinite(neighbourCell) /*&& neighbour_cell->info().isvisited==ref*/) {
- pass++;
- STDEV += pow(((cell->info().kNorm())[j]-meanK),2);
- }
- }
- }
- STDEV = sqrt(STDEV/pass);
- if (debugOut) cout << "PassSTDEV = " << pass << endl;
- cout << "STATISTIC K" << endl;
- double k_min = 0, k_max = meanK + KOptFactor*STDEV;
- cout << "Kmoy = " << meanK << " Standard Deviation = " << STDEV << endl;
- cout << "kmin = " << k_min << " kmax = " << k_max << endl;
- pass=0;
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
- for (int j=0; j<4; j++) {
- neighbourCell = cell->neighbor(j);
- if (!Tri.is_infinite(neighbourCell) /*&& neighbour_cell->info().isvisited==ref*/) {
- pass+=1;
- if ((cell->info().kNorm())[j]>k_max) {
- (cell->info().kNorm())[j]=k_max;
- (neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)]= (cell->info().kNorm())[j];
- }
- k_opt_file << KOptFactor << " " << (cell->info().kNorm())[j] << endl;
- }
- }
- }
- if (debugOut) cout << "PassKopt = " << pass << endl;
- }
- if (debugOut) {
- FiniteVerticesIterator verticesEnd = Tri.finite_vertices_end();
- Real Vgrains = 0;
- int grains=0;
-
- for (FiniteVerticesIterator vIt = Tri.finite_vertices_begin(); vIt != verticesEnd; vIt++) {
- if (!vIt->info().isFictious && !vIt->info().isGhost) {
- grains +=1;
- Vgrains += 1.33333333 * M_PI * pow(vIt->point().weight(),1.5);
- }
- }
- cout<<grains<<"grains - " <<"vTotal = " << vTotal << " Vgrains = " << Vgrains << " vPoral1 = " << (vTotal-Vgrains) << endl;
- cout << "Vtotalissimo = " << Vtotalissimo/2 << " VSolidTot = " << VSolidTot/2 << " vPoral2 = " << vPoral/2 << " sSolidTot = " << sSolidTot << endl<< endl;
-
- if (!RAVERAGE) cout << "------Hydraulic Radius is used for permeability computation------" << endl << endl;
- else cout << "------Average Radius is used for permeability computation------" << endl << endl;
- cout << "-----Computed_Permeability Periodic-----" << endl;
- }
-}
-
-
-
-
-void PeriodicFlow::gaussSeidel(Real dt)
-{
- RTriangulation& Tri = T[currentTes].Triangulation();
- int j = 0;
- double m, n, dp_max, p_max, sum_p, dp, sum_dp;
- double compFlowFactor=0;
- vector<Real> previousP;
- previousP.resize(Tri.number_of_finite_cells());
- const int num_threads=1;
- bool compressible= fluidBulkModulus>0;
-
- vector<Real> t_sum_p, t_dp_max, t_sum_dp, t_p_max;
- t_sum_dp.resize(num_threads);
- t_dp_max.resize(num_threads);
- t_p_max.resize(num_threads);
- t_sum_p.resize(num_threads);
- FiniteCellsIterator cellEnd = Tri.finite_cells_end();
- do {
- int cell2=0;
- dp_max = 0;
- p_max = 0;
- sum_p=0;
- sum_dp=0;
- int bb=-1;
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
- bb++;
- if ( !cell->info().Pcondition && !cell->info().isGhost) {
- cell2++;
- if (compressible && j==0) previousP[bb]=cell->info().shiftedP();
- m=0, n=0;
- for (int j2=0; j2<4; j2++) {
- if (!Tri.is_infinite(cell->neighbor(j2))) {
- if ( compressible ) {
- compFlowFactor = fluidBulkModulus*dt*cell->info().invVoidVolume();
- m += compFlowFactor*(cell->info().kNorm())[j2]*cell->neighbor(j2)->info().shiftedP();
- if (j==0) n += compFlowFactor*(cell->info().kNorm())[j2];
- } else {
- m += (cell->info().kNorm())[j2]*cell->neighbor(j2)->info().shiftedP();
- if ( isinf(m) && j<10 ) cout << "(cell->info().kNorm())[j2] = " << (cell->info().kNorm())[j2] << " cell->neighbor(j2)->info().shiftedP() = " << cell->neighbor(j2)->info().shiftedP() << endl;
- if (j==0) n += (cell->info().kNorm())[j2];
- }
- }
- }
- dp = cell->info().p();
- if (n!=0 || j!=0) {
- if (j==0) { if (compressible) cell->info().invSumK=1/(1+n); else cell->info().invSumK=1/n; }
- if ( compressible ) {
- cell->info().setP( ( ((previousP[bb] - ((fluidBulkModulus*dt*cell->info().invVoidVolume())*(cell->info().dv()))) + m) * cell->info().invSumK - cell->info().shiftedP()) * relax + cell->info().shiftedP());
- } else {
- cell->info().setP((-(cell->info().dv()-m)*cell->info().invSumK-cell->info().p())*relax+cell->info().shiftedP());
- }
- }
- dp -= cell->info().p();
- dp_max = max(dp_max, std::abs(dp));
- p_max = max(p_max, std::abs(cell->info().shiftedP()));
- sum_p += cell->info().shiftedP();
- sum_dp += std::abs(dp);
- }
- }
- j++;
-
- if (j>=40000) cerr<<"\r GS not converged after 40k iterations, break";
-
- } while ((dp_max/p_max) > tolerance && j<40000 /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
-
- int cel=0;
- double Pav=0;
- for (FiniteCellsIterator cell = Tri.finite_cells_begin();
- cell != cellEnd;
- cell++) {
- cel++;
- Pav+=cell->info().shiftedP();
- }
- Pav/=cel;
-}
-
-void PeriodicFlow::displayStatistics()
-{
- RTriangulation& Tri = T[currentTes].Triangulation();
- int Zero =0, Inside=0, Fictious=0, ghostC=0,realC=0, ghostV=0, realV=0;
- FiniteCellsIterator cellEnd = Tri.finite_cells_end();
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
- int zeros =0;
- for (int j=0; j!=4; j++) {
- if ((cell->info().kNorm())[j]==0) {
- zeros+=1;
- }
- }
- if (zeros==4) {
- Zero+=1;
- }
- if (!cell->info().fictious()) {
- Inside+=1;
- }
- else {
- Fictious+=1;
- }
- if (cell->info().isGhost) ghostC+=1; else realC+=1;
- }
- int fict=0, real=0;
- for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
- if (v->info().isFictious) fict+=1;
- else real+=1;
- }
- long Vertices = Tri.number_of_vertices();
- long Cells = Tri.number_of_finite_cells();
- long Facets = Tri.number_of_finite_facets();
- if(debugOut) {
- cout << "zeros = " << Zero << endl;
- cout << "There are " << Vertices << " vertices, dont " << fict << " fictious et " << real << " reeeeeel" << std::endl;
- cout << "There are " << ghostV+realV << " vertices, dont " << ghostV << " ghost et " << realV << " reeeeeel" << std::endl;
- cout << "There are " << ghostC+realC << " cells, dont " << ghostC << " ghost et " << realC << " reeeeeel" << std::endl;
- cout << "There are " << Cells << " cells " << std::endl;
- cout << "There are " << Facets << " facets " << std::endl;
- cout << "There are " << Inside << " cells INSIDE." << endl;
- cout << "There are " << Fictious << " cells FICTIOUS." << endl;
- }
-
- vtkInfiniteVertices = fict;
- vtkInfiniteCells = Fictious;
- num_particles = real;
-}
-
-double PeriodicFlow::boundaryFlux(unsigned int boundaryId)
-{
- RTriangulation& Tri = T[currentTes].Triangulation();
- double Q1=0;
-
- VectorCell tmpCells;
- tmpCells.resize(10000);
- VCellIterator cells_it = tmpCells.begin();
-
- VCellIterator cell_up_end = Tri.incident_cells(T[currentTes].vertexHandles[boundaryId],cells_it);
- for (VCellIterator it = tmpCells.begin(); it != cell_up_end; it++)
- {
- const CellHandle& cell = *it;
- if (cell->info().isGhost) continue;
- Q1 -= cell->info().dv();
- for (int j2=0; j2<4; j2++)
- Q1 += (cell->info().kNorm())[j2]* (cell->neighbor(j2)->info().shiftedP()-cell->info().shiftedP());
- }
- return Q1;
-}
-
-void PeriodicFlow::computeEdgesSurfaces()
-{
- RTriangulation& Tri = T[currentTes].Triangulation();
-//first, copy interacting pairs and normal lub forces form prev. triangulation in a sorted structure for initializing the new lub. Forces
- vector<vector<pair<unsigned int,Real> > > lubPairs;
- lubPairs.resize(Tri.number_of_vertices()+1);
- for (unsigned int k=0; k<edgeNormalLubF.size(); k++)
- lubPairs[min(edgeIds[k].first,edgeIds[k].second)].push_back(pair<int,Real> (max(edgeIds[k].first,edgeIds[k].second),edgeNormalLubF[k]));
-
- //Now we reset the containers and initialize them
-
- edgeSurfaces.clear(); edgeIds.clear(); edgeNormalLubF.clear();
- FiniteEdgesIterator ed_it;
- for ( FiniteEdgesIterator ed_it = Tri.finite_edges_begin(); ed_it!=Tri.finite_edges_end();ed_it++ )
- {
- int hasFictious= (ed_it->first)->vertex(ed_it->second)->info().isFictious + (ed_it->first)->vertex(ed_it->third)->info().isFictious;
- if (hasFictious==2) continue;
- if (((ed_it->first)->vertex(ed_it->second)->info().isGhost) && ((ed_it->first)->vertex(ed_it->third)->info().isGhost)) continue;
- if (((ed_it->first)->vertex(ed_it->second)->info().isGhost) && ((ed_it->first)->vertex(ed_it->third)->info().isFictious)) continue;
- if (((ed_it->first)->vertex(ed_it->second)->info().isFictious) && ((ed_it->first)->vertex(ed_it->third)->info().isGhost)) continue;
-
- unsigned int id1 = (ed_it->first)->vertex(ed_it->second)->info().id();
- unsigned int id2 = (ed_it->first)->vertex(ed_it->third)->info().id();
- double area = T[currentTes].ComputeVFacetArea(ed_it);
- edgeSurfaces.push_back(area);
- edgeIds.push_back(pair<int,int>(id1,id2));
-
- //For persistant edges, we must transfer the lub. force value from the older triangulation structure
- if (id1>id2) swap(id1,id2);
- unsigned int i=0;
- //Look for the pair (id1,id2) in lubPairs
- while (i<lubPairs[id1].size()) {
- if (lubPairs[id1][i].first == id2) {
- //it's found, we copy the lub force
- edgeNormalLubF.push_back(lubPairs[id1][i].second);
- break;}
- ++i;
- }
- // not found, we initialize with zero lub force
- if (i==lubPairs[id1].size()) edgeNormalLubF.push_back(0);
-
- }
-}
-
-} //namespace CGT
-
-#endif //FLOW_ENGINE
-
-
-#ifdef LINSOLV
-#include "PeriodicFlowLinSolv.ipp"
-#endif
=== modified file 'lib/triangulation/PeriodicFlow.hpp'
--- lib/triangulation/PeriodicFlow.hpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/PeriodicFlow.hpp 2014-03-21 18:47:45 +0000
@@ -16,17 +16,28 @@
namespace CGT{
- typedef CGT::FlowBoundingSphere<PeriFlowTesselation> PeriodicFlowBoundingSphere;
- class PeriodicFlow : public PeriodicFlowBoundingSphere
+// typedef CGT::FlowBoundingSphere<PeriFlowTesselation> PeriodicFlowBoundingSphere;
+ template<class _Tesselation>
+ class PeriodicFlow : public CGT::FlowBoundingSphere<_Tesselation>
{
public:
+ typedef _Tesselation Tesselation;
+ typedef Network<Tesselation> _N;
+ DECLARE_TESSELATION_TYPES(Network<Tesselation>)
+ typedef CGT::FlowBoundingSphere<_Tesselation> BaseFlowSolver;
+
+ //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;
+ using BaseFlowSolver::noCache; using BaseFlowSolver::rAverage; using BaseFlowSolver::distanceCorrection; using BaseFlowSolver::minPermLength; using BaseFlowSolver::checkSphereFacetOverlap; using BaseFlowSolver::viscosity; using BaseFlowSolver::kFactor; using BaseFlowSolver::permeabilityMap; using BaseFlowSolver::maxKdivKmean; using BaseFlowSolver::clampKValues; using BaseFlowSolver::KOptFactor; using BaseFlowSolver::meanKStat; using BaseFlowSolver::fluidBulkModulus; using BaseFlowSolver::relax; using BaseFlowSolver::tolerance; using BaseFlowSolver::minKdivKmean;
+
+ //same for functions
+ using _N::defineFictiousCells; using _N::addBoundingPlanes; using _N::boundary;
+
void interpolate(Tesselation& Tes, Tesselation& NewTes);
void computeFacetForcesWithCache(bool onlyCache=false);
void computePermeability();
void gaussSeidel(Real dt=0);
void displayStatistics();
- void computeEdgesSurfaces();
- double boundaryFlux(unsigned int boundaryId);
#ifdef EIGENSPARSE_LIB
//Eigen's sparse matrix for forces computation
// Eigen::SparseMatrix<double> FIntegral;
@@ -35,7 +46,435 @@
// Eigen::VectorXd forces;
#endif
};
-}
+
+template<class _Tesselation>
+void PeriodicFlow<_Tesselation>::interpolate(Tesselation& Tes, Tesselation& NewTes)
+{
+ CellHandle oldCell;
+ RTriangulation& Tri = Tes.Triangulation();
+ for (VCellIterator cellIt=NewTes.cellHandles.begin(); cellIt!=NewTes.cellHandles.end(); cellIt++){
+ CellHandle& newCell = *cellIt;
+ if (newCell->info().Pcondition || newCell->info().isGhost) continue;
+ CVector center ( 0,0,0 );
+ if (newCell->info().fictious()==0) for ( int k=0;k<4;k++ ) center= center + 0.25* (Tes.vertex(newCell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
+ else {
+ Real boundPos=0; int coord=0;
+ for ( int k=0;k<4;k++ ) {
+ if (!newCell->vertex (k)->info().isFictious) center= center+0.3333333333*(Tes.vertex(newCell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
+ else {
+ coord=boundary (newCell->vertex(k)->info().id()).coordinate;
+ boundPos=boundary (newCell->vertex(k)->info().id()).p[coord];
+ }
+ }
+ center=CVector(coord==0?boundPos:center[0],coord==1?boundPos:center[1],coord==2?boundPos:center[2]);
+ }
+ oldCell = Tri.locate(Point(center[0],center[1],center[2]));
+ newCell->info().p() = oldCell->info().shiftedP();
+ }
+// Tes.Clear();//Don't reset to avoid segfault when getting pressure in scripts just after interpolation
+}
+
+
+template<class _Tesselation>
+void PeriodicFlow<_Tesselation>::computeFacetForcesWithCache(bool onlyCache)
+{
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ CVector nullVect(0,0,0);
+ static vector<CVector> oldForces;
+ if (oldForces.size()<=Tri.number_of_vertices()) oldForces.resize(Tri.number_of_vertices()+1);
+ //reset forces
+ for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+ if (noCache) {oldForces[v->info().id()]=nullVect; v->info().forces=nullVect;}
+ else {oldForces[v->info().id()]=v->info().forces; v->info().forces=nullVect;}
+ }
+ CellHandle neighbourCell;
+ VertexHandle mirrorVertex;
+ CVector tempVect;
+
+ //FIXME : Ema, be carefull with this (noCache), it needs to be turned true after retriangulation
+ if (noCache) {
+ for (VCellIterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
+ CellHandle& cell = *cellIt;
+ for (int k=0;k<4;k++) cell->info().unitForceVectors[k]=nullVect;
+ for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j))) {
+// #ifdef EIGENSPARSE_LIB
+// if (!cell->info().Pcondition) ++nPCells;
+// #endif
+ neighbourCell = cell->neighbor(j);
+ const CVector& Surfk = cell->info().facetSurfaces[j];
+ //FIXME : later compute that fluidSurf only once in hydraulicRadius, for now keep full surface not modified in cell->info for comparison with other forces schemes
+ //The ratio void surface / facet surface
+ Real area = sqrt(Surfk.squared_length());
+ CVector facetNormal = Surfk/area;
+ const std::vector<CVector>& crossSections = cell->info().facetSphereCrossSections;
+ CVector fluidSurfk = cell->info().facetSurfaces[j]*cell->info().facetFluidSurfacesRatio[j];
+ /// handle fictious vertex since we can get the projected surface easily here
+ if (cell->vertex(j)->info().isFictious) {
+ Real projSurf=abs(Surfk[boundary(cell->vertex(j)->info().id()).coordinate]);
+ tempVect=-projSurf*boundary(cell->vertex(j)->info().id()).normal;
+ //define the cached value for later use with cache*p
+ cell->info().unitForceVectors[j]=cell->info().unitForceVectors[j]+ tempVect;
+ }
+ /// Apply weighted forces f_k=sqRad_k/sumSqRad*f
+ CVector facetUnitForce = -fluidSurfk*cell->info().solidSurfaces[j][3];
+ for (int y=0; y<3;y++) {
+ //add to cached value
+ cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]+facetUnitForce*cell->info().solidSurfaces[j][y];
+ //uncomment to get total force / comment to get only viscous forces (Bruno)
+ if (!cell->vertex(facetVertices[j][y])->info().isFictious) {
+ //add to cached value
+ cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]-facetNormal*crossSections[j][y];
+ }
+ }
+ }
+ }
+
+ noCache=false;//cache should always be defined after execution of this function
+ if (onlyCache) return;
+ }// end if(noCache)
+
+ //use cached values
+ //First define products that will be used below for all cells:
+ Real pDeltas [3];
+ for (unsigned int k=0; k<3;k++) pDeltas[k]=CellInfo::hSize[k]*CellInfo::gradP;
+ //Then compute the forces
+ for (VCellIterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
+ const CellHandle& cell = *cellIt;
+ for (int yy=0;yy<4;yy++) {
+ VertexInfo& vhi = cell->vertex(yy)->info();
+ Real unshiftedP = cell->info().p();
+ //the pressure translated to a ghost cell adjacent to the non-ghost vertex
+ unshiftedP -= pDeltas[0]*vhi.period[0] + pDeltas[1]*vhi.period[1] +pDeltas[2]*vhi.period[2];
+ T[currentTes].vertexHandles[vhi.id()]->info().forces=T[currentTes].vertexHandles[vhi.id()]->info().forces + cell->info().unitForceVectors[yy]*unshiftedP;
+ }
+ }
+ if (debugOut) {
+ CVector totalForce = nullVect;
+ for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); v++)
+ {
+ if (!v->info().isFictious /*&& !v->info().isGhost*/ ){
+ totalForce = totalForce + v->info().forces;}
+ else /*if (!v->info().isGhost)*/{
+ if (boundary(v->info().id()).flowCondition==1) totalForce = totalForce + v->info().forces;
+ }
+ }
+ cout << "totalForce = "<< totalForce << endl;}
+}
+
+template<class _Tesselation>
+void PeriodicFlow<_Tesselation>::computePermeability()
+{
+ if (debugOut) cout << "----Computing_Permeability (Periodic)------" << endl;
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ VSolidTot = 0, Vtotalissimo = 0, vPoral = 0, sSolidTot = 0, vTotalPorosity=0, vPoralPorosity=0;
+ CellHandle neighbourCell;
+
+ double k=0, distance = 0, radius = 0;
+ int surfneg=0; int NEG=0, POS=0, pass=0;
+ Real meanK=0, STDEV=0, meanRadius=0, meanDistance=0;
+ Real infiniteK=1e3;
+ double volume_sub_pore = 0.f;
+ VectorCell& cellHandles= T[currentTes].cellHandles;
+ for (VCellIterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
+ CellHandle& cell = *cellIt;
+ Point& p1 = cell->info();
+ for (int j=0; j<4; j++){
+ neighbourCell = cell->neighbor(j);
+ Point& p2 = neighbourCell->info();
+ if (!Tri.is_infinite(neighbourCell) /*&& (neighbour_cell->info().isvisited==ref || computeAllCells)*/) {
+ //compute and store the area of sphere-facet intersections for later use
+ VertexHandle W [3];
+ for (int kk=0; kk<3; kk++) {
+ W[kk] = cell->vertex(facetVertices[j][kk]);
+ }
+ Sphere& v0 = W[0]->point();
+ Sphere& v1 = W[1]->point();
+ Sphere& v2 = W[2]->point();
+#ifdef USE_FAST_MATH
+ //FIXME : code not compiling,, do the same as in "else"
+ assert((W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point())>=0 && (W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point())<=1);
+ for (int jj=0;jj<3;jj++)
+ cell->info().facetSphereCrossSections[j][jj]=0.5*W[jj]->point().weight()*Wm3::FastInvCos1((W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point()));
+#else
+ cell->info().facetSphereCrossSections[j]=CVector(
+ W[0]->info().isFictious ? 0 : 0.5*v0.weight()*acos((v1-v0)*(v2-v0)/sqrt((v1-v0).squared_length()*(v2-v0).squared_length())),
+ W[1]->info().isFictious ? 0 : 0.5*v1.weight()*acos((v0-v1)*(v2-v1)/sqrt((v1-v0).squared_length()*(v2-v1).squared_length())),
+ W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0-v2)*(v1-v2)/sqrt((v1-v2).squared_length()*(v2-v0).squared_length())));
+#endif
+ pass+=1;
+ CVector l = p1 - p2;
+ distance = sqrt(l.squared_length());
+ if (!rAverage) radius = 2* computeHydraulicRadius(cell, j);
+ else radius = (computeEffectiveRadius(cell, j)+computeEquivalentRadius(cell,j))*0.5;
+ if (radius<0) NEG++;
+ else POS++;
+ if (radius==0) {
+ cout << "INS-INS PROBLEM!!!!!!!" << endl;
+ }
+ Real fluidArea=0;
+ int test=0;
+ if (distance!=0) {
+ if (minPermLength>0 && distanceCorrection) distance=max(minPermLength,distance);
+ const CVector& Surfk = cell->info().facetSurfaces[j];
+ Real area = sqrt(Surfk.squared_length());
+ const CVector& crossSections = cell->info().facetSphereCrossSections[j];
+ Real S0=0;
+ S0=checkSphereFacetOverlap(v0,v1,v2);
+ if (S0==0) S0=checkSphereFacetOverlap(v1,v2,v0);
+ if (S0==0) S0=checkSphereFacetOverlap(v2,v0,v1);
+ //take absolute value, since in rare cases the surface can be negative (overlaping spheres)
+ fluidArea=abs(area-crossSections[0]-crossSections[1]-crossSections[2]+S0);
+ cell->info().facetFluidSurfacesRatio[j]=fluidArea/area;
+ k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);
+ meanDistance += distance;
+ meanRadius += radius;
+ meanK += k*kFactor;
+
+ if (k<0 && debugOut) {surfneg+=1;
+ cout<<"__ k<0 __"<<k<<" "<<" fluidArea "<<fluidArea<<" area "<<area<<" "<<crossSections[0]<<" "<<crossSections[1]<<" "<<crossSections[2] <<" "<<W[0]->info().id()<<" "<<W[1]->info().id()<<" "<<W[2]->info().id()<<" "<<p1<<" "<<p2<<" test "<<test<<endl;}
+
+ } else cout <<"infinite K2! surfaces will be missing (FIXME)"<<endl; k = infiniteK;
+ //Will be corrected in the next loop
+ (cell->info().kNorm())[j]= k*kFactor;
+ if (!neighbourCell->info().isGhost) (neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)]= (cell->info().kNorm())[j];
+ //The following block is correct but not very usefull, since all values are clamped below with MIN and MAX, skip for now
+// else {//find the real neighbor connected to our cell through periodicity
+// CellHandle true_neighbour_cell = cellHandles[neighbour_cell->info().baseIndex];
+// for (int ii=0;ii<4;ii++)
+// if (true_neighbour_cell->neighbor(ii)->info().index == cell->info().index){
+// (true_neighbour_cell->info().kNorm())[ii]=(cell->info().kNorm())[j]; break;
+// }
+// }
+ if(permeabilityMap){
+ CellHandle c = cell;
+ cell->info().s = cell->info().s + k*distance/fluidArea*volumePoreVoronoiFraction (c,j);
+ volume_sub_pore += volumePoreVoronoiFraction (c,j);
+ }
+ }
+ }
+// cell->info().isvisited = !ref;
+ if(permeabilityMap){
+ cell->info().s = cell->info().s/volume_sub_pore;
+ volume_sub_pore = 0.f;
+ }
+// }
+ }
+
+
+ if (debugOut) cout<<"surfneg est "<<surfneg<<endl;
+ meanK /= pass;
+ meanRadius /= pass;
+ meanDistance /= pass;
+ Real globalK=kFactor*meanDistance*vPoral/(sSolidTot*8.*viscosity);//An approximate value of macroscopic permeability, for clamping local values below
+ if (debugOut) {
+ cout << "PassCompK = " << pass << endl;
+ cout << "meanK = " << meanK << endl;
+ cout << "globalK = " << globalK << endl;
+ cout << "maxKdivKmean*globalK = " << maxKdivKmean*globalK << endl;
+ cout << "minKdivKmean*globalK = " << minKdivKmean*globalK << endl;
+ cout << "meanTubesRadius = " << meanRadius << endl;
+ cout << "meanDistance = " << meanDistance << endl;
+ }
+ pass=0;
+ if (clampKValues) for (VCellIterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
+ CellHandle& cell = *cellIt;
+ for (int j=0; j<4; j++) {
+ neighbourCell = cell->neighbor(j);
+ if (!Tri.is_infinite(neighbourCell) /*&& neighbour_cell->info().isvisited==ref*/) {
+ pass++;
+ (cell->info().kNorm())[j] = max(minKdivKmean*globalK, min((cell->info().kNorm())[j], maxKdivKmean*globalK));
+ (neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)]=(cell->info().kNorm())[j];
+ if (!neighbourCell->info().isGhost) (neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)]= (cell->info().kNorm())[j];
+ else {//find the real neighbor connected to our cell through periodicity, as we want exactly the same permeability without rounding errors
+ CellHandle& true_neighbourCell = cellHandles[neighbourCell->info().baseIndex];
+ for (int ii=0;ii<4;ii++)
+ if (true_neighbourCell->neighbor(ii)->info().index == cell->info().index){
+ (true_neighbourCell->info().kNorm())[ii]=(cell->info().kNorm())[j]; break;
+ }
+ }
+
+ }
+ }
+ }
+ if (debugOut) cout << "PassKcorrect = " << pass << endl;
+
+ if (debugOut) cout << "POS = " << POS << " NEG = " << NEG << " pass = " << pass << endl;
+
+// A loop to compute the standard deviation of the local K distribution, and use it to include/exclude K values higher then (meanK +/- K_opt_factor*STDEV)
+ if (meanKStat)
+ {
+ std::ofstream k_opt_file("k_stdev.txt" ,std::ios::out);
+ pass=0;
+ FiniteCellsIterator cellEnd = Tri.finite_cells_end();
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
+ for (int j=0; j<4; j++) {
+ neighbourCell = cell->neighbor(j);
+ if (!Tri.is_infinite(neighbourCell) /*&& neighbour_cell->info().isvisited==ref*/) {
+ pass++;
+ STDEV += pow(((cell->info().kNorm())[j]-meanK),2);
+ }
+ }
+ }
+ STDEV = sqrt(STDEV/pass);
+ if (debugOut) cout << "PassSTDEV = " << pass << endl;
+ cout << "STATISTIC K" << endl;
+ double k_min = 0, k_max = meanK + KOptFactor*STDEV;
+ cout << "Kmoy = " << meanK << " Standard Deviation = " << STDEV << endl;
+ cout << "kmin = " << k_min << " kmax = " << k_max << endl;
+ pass=0;
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
+ for (int j=0; j<4; j++) {
+ neighbourCell = cell->neighbor(j);
+ if (!Tri.is_infinite(neighbourCell) /*&& neighbour_cell->info().isvisited==ref*/) {
+ pass+=1;
+ if ((cell->info().kNorm())[j]>k_max) {
+ (cell->info().kNorm())[j]=k_max;
+ (neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)]= (cell->info().kNorm())[j];
+ }
+ k_opt_file << KOptFactor << " " << (cell->info().kNorm())[j] << endl;
+ }
+ }
+ }
+ if (debugOut) cout << "PassKopt = " << pass << endl;
+ }
+ if (debugOut) {
+ FiniteVerticesIterator verticesEnd = Tri.finite_vertices_end();
+ Real Vgrains = 0;
+ int grains=0;
+
+ for (FiniteVerticesIterator vIt = Tri.finite_vertices_begin(); vIt != verticesEnd; vIt++) {
+ if (!vIt->info().isFictious && !vIt->info().isGhost) {
+ grains +=1;
+ Vgrains += 1.33333333 * M_PI * pow(vIt->point().weight(),1.5);
+ }
+ }
+ cout<<grains<<"grains - " <<"vTotal = " << vTotal << " Vgrains = " << Vgrains << " vPoral1 = " << (vTotal-Vgrains) << endl;
+ cout << "Vtotalissimo = " << Vtotalissimo/2 << " VSolidTot = " << VSolidTot/2 << " vPoral2 = " << vPoral/2 << " sSolidTot = " << sSolidTot << endl<< endl;
+
+ if (!rAverage) cout << "------Hydraulic Radius is used for permeability computation------" << endl << endl;
+ else cout << "------Average Radius is used for permeability computation------" << endl << endl;
+ cout << "-----computed_Permeability Periodic-----" << endl;
+ }
+}
+
+
+
+
+template<class _Tesselation>
+void PeriodicFlow<_Tesselation>::gaussSeidel(Real dt)
+{
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ int j = 0;
+ double m, n, dp_max, p_max, sum_p, dp, sum_dp;
+ double compFlowFactor=0;
+ vector<Real> previousP;
+ previousP.resize(Tri.number_of_finite_cells());
+ const int num_threads=1;
+ bool compressible= fluidBulkModulus>0;
+
+ vector<Real> t_sum_p, t_dp_max, t_sum_dp, t_p_max;
+ t_sum_dp.resize(num_threads);
+ t_dp_max.resize(num_threads);
+ t_p_max.resize(num_threads);
+ t_sum_p.resize(num_threads);
+ FiniteCellsIterator cellEnd = Tri.finite_cells_end();
+ do {
+ int cell2=0;
+ dp_max = 0;
+ p_max = 0;
+ sum_p=0;
+ sum_dp=0;
+ int bb=-1;
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
+ bb++;
+ if ( !cell->info().Pcondition && !cell->info().isGhost) {
+ cell2++;
+ if (compressible && j==0) previousP[bb]=cell->info().shiftedP();
+ m=0, n=0;
+ for (int j2=0; j2<4; j2++) {
+ if (!Tri.is_infinite(cell->neighbor(j2))) {
+ if ( compressible ) {
+ compFlowFactor = fluidBulkModulus*dt*cell->info().invVoidVolume();
+ m += compFlowFactor*(cell->info().kNorm())[j2]*cell->neighbor(j2)->info().shiftedP();
+ if (j==0) n += compFlowFactor*(cell->info().kNorm())[j2];
+ } else {
+ m += (cell->info().kNorm())[j2]*cell->neighbor(j2)->info().shiftedP();
+ if ( isinf(m) && j<10 ) cout << "(cell->info().kNorm())[j2] = " << (cell->info().kNorm())[j2] << " cell->neighbor(j2)->info().shiftedP() = " << cell->neighbor(j2)->info().shiftedP() << endl;
+ if (j==0) n += (cell->info().kNorm())[j2];
+ }
+ }
+ }
+ dp = cell->info().p();
+ if (n!=0 || j!=0) {
+ if (j==0) { if (compressible) cell->info().invSumK=1/(1+n); else cell->info().invSumK=1/n; }
+ if ( compressible ) {
+ cell->info().setP( ( ((previousP[bb] - ((fluidBulkModulus*dt*cell->info().invVoidVolume())*(cell->info().dv()))) + m) * cell->info().invSumK - cell->info().shiftedP()) * relax + cell->info().shiftedP());
+ } else {
+ cell->info().setP((-(cell->info().dv()-m)*cell->info().invSumK-cell->info().p())*relax+cell->info().shiftedP());
+ }
+ }
+ dp -= cell->info().p();
+ dp_max = max(dp_max, std::abs(dp));
+ p_max = max(p_max, std::abs(cell->info().shiftedP()));
+ sum_p += cell->info().shiftedP();
+ sum_dp += std::abs(dp);
+ }
+ }
+ j++;
+
+ if (j>=40000) cerr<<"\r GS not converged after 40k iterations, break";
+
+ } while ((dp_max/p_max) > tolerance && j<40000 /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
+
+ int cel=0;
+ double Pav=0;
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin();
+ cell != cellEnd;
+ cell++) {
+ cel++;
+ Pav+=cell->info().shiftedP();
+ }
+ Pav/=cel;
+}
+
+template<class _Tesselation>
+void PeriodicFlow<_Tesselation>::displayStatistics()
+{
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ int Zero =0, Inside=0, Fictious=0, ghostC=0,realC=0, ghostV=0, realV=0;
+ FiniteCellsIterator cellEnd = Tri.finite_cells_end();
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
+ int zeros =0;
+ for (int j=0; j!=4; j++)
+ if ((cell->info().kNorm())[j]==0) zeros+=1;
+ if (zeros==4) Zero+=1;
+ if (!cell->info().fictious()) Inside+=1; else Fictious+=1;
+ if (cell->info().isGhost) ghostC+=1; else realC+=1;
+ }
+ int fict=0, real=0;
+ for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+ if (v->info().isFictious) fict+=1;
+ else real+=1;
+ }
+ long Vertices = Tri.number_of_vertices();
+ long Cells = Tri.number_of_finite_cells();
+ long Facets = Tri.number_of_finite_facets();
+ if(debugOut) {
+ cout << "zeros = " << Zero << endl;
+ cout << "There are " << Vertices << " vertices, dont " << fict << " fictious et " << real << " reeeeeel" << std::endl;
+ cout << "There are " << ghostV+realV << " vertices, dont " << ghostV << " ghost et " << realV << " reeeeeel" << std::endl;
+ cout << "There are " << ghostC+realC << " cells, dont " << ghostC << " ghost et " << realC << " reeeeeel" << std::endl;
+ cout << "There are " << Cells << " cells " << std::endl;
+ cout << "There are " << Facets << " facets " << std::endl;
+ cout << "There are " << Inside << " cells INSIDE." << endl;
+ cout << "There are " << Fictious << " cells FICTIOUS." << endl;
+ }
+ vtkInfiniteVertices = fict;
+ vtkInfiniteCells = Fictious;
+ num_particles = real;
+}
+
+} //END NAMESPACE
#ifdef LINSOLV
#include "PeriodicFlowLinSolv.hpp"
=== modified file 'lib/triangulation/PeriodicFlowLinSolv.hpp'
--- lib/triangulation/PeriodicFlowLinSolv.hpp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/PeriodicFlowLinSolv.hpp 2014-03-21 18:47:45 +0000
@@ -15,12 +15,25 @@
namespace CGT {
-typedef FlowBoundingSphereLinSolv<PeriodicFlow> LinSolver;
-
-class PeriodicFlowLinSolv : public LinSolver
+template<class _Tesselation>
+class PeriodicFlowLinSolv : public FlowBoundingSphereLinSolv<_Tesselation,PeriodicFlow<_Tesselation> >
{
public:
- typedef PeriodicFlow FlowType;
+ typedef _Tesselation Tesselation;
+ typedef Network<Tesselation> _N;
+ DECLARE_TESSELATION_TYPES(Network<Tesselation>)
+ typedef FlowBoundingSphereLinSolv<_Tesselation,PeriodicFlow<_Tesselation> > BaseFlowSolver;
+ typedef typename BaseFlowSolver::ETriplet ETriplet;
+
+ ///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 BaseFlowSolver::noCache; using BaseFlowSolver::rAverage; using BaseFlowSolver::distanceCorrection; using BaseFlowSolver::minPermLength; using BaseFlowSolver::checkSphereFacetOverlap; using BaseFlowSolver::viscosity; using BaseFlowSolver::kFactor; using BaseFlowSolver::permeabilityMap; using BaseFlowSolver::maxKdivKmean; using BaseFlowSolver::clampKValues; using BaseFlowSolver::KOptFactor; using BaseFlowSolver::meanKStat; using BaseFlowSolver::fluidBulkModulus; using BaseFlowSolver::relax; using BaseFlowSolver::tolerance; using BaseFlowSolver::minKdivKmean;
+ /// More members from LinSolv variant
+ using BaseFlowSolver::areCellsOrdered; using BaseFlowSolver::T_nnz; using BaseFlowSolver::ncols; using BaseFlowSolver::T_cells; using BaseFlowSolver::T_index; using BaseFlowSolver::orderedCells; using BaseFlowSolver::isLinearSystemSet; using BaseFlowSolver::T_x; using BaseFlowSolver::T_b; using BaseFlowSolver::T_bv; using BaseFlowSolver::bodv; using BaseFlowSolver::xodv; using BaseFlowSolver::errorCode; using BaseFlowSolver::useSolver; using BaseFlowSolver::tripletList; using BaseFlowSolver::A; using BaseFlowSolver::gsP; using BaseFlowSolver::gsB; using BaseFlowSolver::fullAcolumns; using BaseFlowSolver::fullAvalues; using BaseFlowSolver::isFullLinearSystemGSSet; using BaseFlowSolver::gsdV;
+
vector<int> indices;//redirection vector containing the rank of cell so that T_cells[indices[cell->info().index]]=cell
virtual ~PeriodicFlowLinSolv();
@@ -32,5 +45,5 @@
};
} //namespace CGTF
-
+#include "PeriodicFlowLinSolv.ipp"
#endif //FLOW_ENGINE
=== modified file 'lib/triangulation/PeriodicFlowLinSolv.ipp'
--- lib/triangulation/PeriodicFlowLinSolv.ipp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/PeriodicFlowLinSolv.ipp 2014-03-21 18:47:45 +0000
@@ -41,14 +41,15 @@
double *, int *, int *, int *, int *, int *,
int *, double *, double *, int *, double *);
#endif
-
-PeriodicFlowLinSolv::~PeriodicFlowLinSolv()
+template<class _Tesselation>
+PeriodicFlowLinSolv<_Tesselation>::~PeriodicFlowLinSolv()
{
}
-
-PeriodicFlowLinSolv::PeriodicFlowLinSolv(): LinSolver() {}
-
-int PeriodicFlowLinSolv::setLinearSystem(Real dt)
+template<class _Tesselation>
+PeriodicFlowLinSolv<_Tesselation>::PeriodicFlowLinSolv(): BaseFlowSolver() {}
+
+template<class _Tesselation>
+int PeriodicFlowLinSolv<_Tesselation>::setLinearSystem(Real dt)
{
//WARNING : boundary conditions (Pcondition, p values) must have been set for a correct definition of
RTriangulation& Tri = T[currentTes].Triangulation();
@@ -135,66 +136,8 @@
}
}
if (!isLinearSystemSet) {
- if (useSolver==1 || useSolver==2){
- #ifdef TAUCS_LIB
- clen.resize(ncols+1);
- T_jn.resize(ncols+1);
- T_A->colptr = &T_jn[0];
- T_ia.resize(T_nnz);
- T_A->rowind = &T_ia[0];
- T_A->flags = (TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER);
- T_an.resize(T_nnz);
- T_A->values.d = &T_an[0];
- T_A->n = ncols;
- T_A->m = ncols;
- int i,j,k;
- for (j=0; j<ncols; j++) clen[j] = 0;
- for (k=0; k<T_nnz; k++) {
- i = is[k]-1; /* make it 1-based */
- j = js[k]-1; /* make it 1-based */
- (clen[j])++;
- }
- /* now compute column pointers */
- k = 0;
- for (j=0; j<ncols; j++) {
- int tmp;
- tmp = clen[j];
- clen[j] = (T_A->colptr[j]) = k;
- k += tmp;
- }
- clen[ncols] = (T_A->colptr[ncols]) = k;
-
- /* now read matrix into data structure */
- for (k=0; k<T_nnz; k++) {
- i = is[k] - 1; /* make it 1-based */
- j = js[k] - 1; /* make it 1-based */
- assert(i < ncols);
- assert(j < ncols);
- (T_A->taucs_values)[clen[j]] = vs[k];
- (T_A->rowind)[clen[j]] = i;
- clen[j] ++;
- }
- #else
- cerr<<"yade compiled without Taucs, FlowEngine.useSolver="<< useSolver <<" not supported"<<endl;
- #endif //TAUCS_LIB
- } else if (useSolver==3){
+ if (useSolver>0){
#ifdef EIGENSPARSE_LIB
-// //here the matrix can be exported in in MatrixMarket format for experiments
-// static int mn=0;
-// ofstream file; ofstream file2; ofstream file3;
-// stringstream ss,ss2,ss3;
-// ss<<"matrix"<<mn;
-// ss3<<"matrixf2"<<mn;
-// ss2<<"matrixf"<<mn++;
-// file.open(ss.str().c_str());
-// file2.open(ss2.str().c_str());
-// file3.open(ss3.str().c_str());
-// file <<"%%MatrixMarket matrix coordinate real symmetric"<<endl;
-// file2 <<"%%MatrixMarket matrix coordinate real symmetric"<<endl;
-// file3 <<"%%MatrixMarket matrix coordinate real symmetric"<<endl;
-// file <<ncols<<" "<<ncols<<" "<<T_nnz<<" -1"<<endl;
-// file2 <<ncols<<" "<<ncols<<" "<<T_nnz<<" -1"<<endl;
-// file3 <<ncols<<" "<<ncols<<" "<<T_nnz<<" -1"<<endl;
tripletList.clear(); tripletList.resize(T_nnz);
for(int k=0;k<T_nnz;k++) {
tripletList[k]=ETriplet(is[k]-1,js[k]-1,vs[k]);
@@ -202,7 +145,7 @@
A.resize(ncols,ncols);
A.setFromTriplets(tripletList.begin(), tripletList.end());
#else
- cerr<<"yade compiled without CHOLMOD, FlowEngine.useSolver="<< useSolver <<" not supported"<<endl;
+ cerr<<"yade compiled without CHOLMOD, FlowEngine.useSolver="<< useSolver <<">0 not supported"<<endl;
#endif
}
isLinearSystemSet=true;
@@ -213,8 +156,8 @@
/// For Gauss Seidel, we need the full matrix
-
-int PeriodicFlowLinSolv::setLinearSystemFullGS(Real dt)
+template<class _Tesselation>
+int PeriodicFlowLinSolv<_Tesselation>::setLinearSystemFullGS(Real dt)
{
//WARNING : boundary conditions (Pcondition, p values) must have been set for a correct definition
RTriangulation& Tri = T[currentTes].Triangulation();
=== modified file 'lib/triangulation/RegularTriangulation.h'
--- lib/triangulation/RegularTriangulation.h 2014-03-21 18:45:24 +0000
+++ lib/triangulation/RegularTriangulation.h 2014-03-21 18:47:45 +0000
@@ -22,6 +22,8 @@
//This include from yade let us use Eigen types
#include <yade/lib/base/Math.hpp>
+const int facetVertices [4][3] = {{1,2,3},{0,2,3},{0,1,3},{0,1,2}};
+
namespace CGT {
//Robust kernel
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
@@ -81,138 +83,138 @@
virtual 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);}
-};
+// 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);}
+// };
///
@@ -249,7 +251,5 @@
};
typedef TriangulationTypes<SimpleVertexInfo,SimpleCellInfo> SimpleTriangulationTypes;
-typedef TriangulationTypes<FlowVertexInfo,FlowCellInfo> FlowTriangulationTypes;
-typedef TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> PeriFlowTriangulationTypes;
} // namespace CGT
=== modified file 'lib/triangulation/Tesselation.h'
--- lib/triangulation/Tesselation.h 2014-03-21 18:45:24 +0000
+++ lib/triangulation/Tesselation.h 2014-03-21 18:47:45 +0000
@@ -69,7 +69,7 @@
protected:
RTriangulation* Tri;
RTriangulation* Tes; //=NULL or Tri depending on the constructor used.
- bool computed;
+
public:
Real TotalFiniteVoronoiVolume;
Real area;
@@ -77,10 +77,10 @@
Real TotalInternalVoronoiPorosity;
VectorVertex vertexHandles;//This is a redirection vector to get vertex pointers by spheres id
VectorCell cellHandles;//for speedup of global loops, iterating on this vector is faster than cellIterator++
- bool redirected;//is vertexHandles filled with current vertex pointers?
+ bool redirected;//is vertexHandles filled with current vertex pointers?
+ bool computed;
-public:
- _Tesselation(void);
+ _Tesselation(void);
_Tesselation(RTriangulation &T);// : Tri(&T) { Calcule(); }
~_Tesselation(void);
@@ -94,9 +94,9 @@
bool remove (unsigned int id);
int Max_id (void) {return maxId;}
- void Compute (); //Calcule le centres de Voronoi pour chaque cellule
+ void compute (); //Calcule le centres de Voronoi pour chaque cellule
void Invalidate () {computed=false;} //Set the tesselation as "not computed" (computed=false), this will launch //tesselation internaly when using functions like computeVolumes())
- // N.B : Compute() must be executed before the functions below are used
+ // N.B : compute() must be executed before the functions below are used
void Clear(void);
static Point Dual (const CellHandle &cell);
@@ -104,10 +104,10 @@
static Segment Dual (FiniteFacetsIterator &facet); //G�n�re le segment dual d'une facette finie
static Real Volume (FiniteCellsIterator cell);
inline void AssignPartialVolume (FiniteEdgesIterator& ed_it);
- double ComputeVFacetArea (FiniteEdgesIterator ed_it);
+ double computeVFacetArea (FiniteEdgesIterator ed_it);
void ResetVCellVolumes (void);
- void ComputeVolumes (void);//Compute volume each voronoi cell
- void ComputePorosity (void);//Compute volume and porosity of each voronoi cell
+ void computeVolumes (void);//compute volume each voronoi cell
+ void computePorosity (void);//compute volume and porosity of each voronoi cell
inline Real& Volume (unsigned int id) { return vertexHandles[id]->info().v(); }
inline const VertexHandle& vertex (unsigned int id) const { return vertexHandles[id]; }
@@ -117,7 +117,7 @@
void voisins (VertexHandle v, VectorVertex& Output_vector);// {Tri->incident_vertices(v, back_inserter(Output_vector));}
RTriangulation& Triangulation (void);// {return *Tri;}
- bool Computed (void) {return computed;}
+// bool computed (void) {return computed;}
bool is_short ( FiniteFacetsIterator f_it );
inline bool is_internal ( FiniteFacetsIterator &facet );//
@@ -152,7 +152,6 @@
//Explicit instanciation
typedef CGT::_Tesselation<CGT::SimpleTriangulationTypes> SimpleTesselation;
-typedef CGT::_Tesselation<CGT::FlowTriangulationTypes> FlowTesselation;
-typedef CGT::PeriodicTesselation<CGT::_Tesselation<CGT::PeriFlowTriangulationTypes> > PeriFlowTesselation;
+
=== modified file 'lib/triangulation/Tesselation.ipp'
--- lib/triangulation/Tesselation.ipp 2014-03-21 18:45:24 +0000
+++ lib/triangulation/Tesselation.ipp 2014-03-21 18:47:45 +0000
@@ -32,7 +32,7 @@
_Tesselation<TT>::_Tesselation ( RTriangulation &T ) : Tri ( &T ), Tes ( &T ), computed ( false )
{
std::cout << "Tesselation(RTriangulation &T)" << std::endl;
- Compute();
+ compute();
}
template<class TT>
@@ -140,7 +140,7 @@
}
template<class TT>
-void _Tesselation<TT>::Compute ()
+void _Tesselation<TT>::compute ()
{
if (!redirected) redirect();
FiniteCellsIterator cellEnd = Tri->finite_cells_end();
@@ -170,7 +170,7 @@
}
template<class TT>
-double _Tesselation<TT>::ComputeVFacetArea ( FiniteEdgesIterator ed_it )
+double _Tesselation<TT>::computeVFacetArea ( FiniteEdgesIterator ed_it )
{
CellCirculator cell0 = Tri->incident_cells ( *ed_it );
CellCirculator cell2 = cell0;
@@ -239,9 +239,9 @@
}
template<class TT>
-void _Tesselation<TT>::ComputeVolumes ( void )
+void _Tesselation<TT>::computeVolumes ( void )
{
- if ( !computed ) Compute();
+ if ( !computed ) compute();
ResetVCellVolumes();
for ( FiniteEdgesIterator ed_it=Tri->finite_edges_begin(); ed_it!=Tri->finite_edges_end();ed_it++ )
{
@@ -251,9 +251,9 @@
}
template<class TT>
-void _Tesselation<TT>::ComputePorosity ( void ) //WARNING : This function will erase real volumes of cells
+void _Tesselation<TT>::computePorosity ( void ) //WARNING : This function will erase real volumes of cells
{
- ComputeVolumes();
+ computeVolumes();
FiniteVerticesIterator verticesEnd = Tri->finite_vertices_end ();
for ( FiniteVerticesIterator vIt = Tri->finite_vertices_begin (); vIt != verticesEnd; vIt++ )
{
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2014-03-21 18:45:24 +0000
+++ pkg/dem/FlowEngine.cpp 2014-03-21 18:47:45 +0000
@@ -27,704 +27,12 @@
#include "FlowEngine.hpp"
+
+
+
CREATE_LOGGER ( FlowEngine );
CREATE_LOGGER ( PeriodicFlowEngine );
-
-CGT::CVector makeCgVect ( const Vector3r& yv ) {return CGT::CVector ( 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] );}
-
-
-FlowEngine::~FlowEngine()
-{
-}
-
-const int facetVertices [4][3] = {{1,2,3},{0,2,3},{0,1,3},{0,1,2}};
-
-void FlowEngine::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()" );
- computeLubricationContributions ( *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(&FlowEngine::backgroundAction,this);
- workerThread.detach();
- if (debug) cerr<<"backgrounded"<<endl;
- initializeVolumes(solver);
- computeLubricationContributions(*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);
- computeLubricationContributions(*solver);
- updateTriangulation = false;
- epsVolCumulative=0;
- retriangulationLastIter=0;
- ReTrg++;}
- }
- first=false;
- timingDeltas->checkpoint ( "triangulate + init volumes" );
-}
-
-void FlowEngine::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 Solver>
-
-void FlowEngine::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 Solver>
-void FlowEngine::setImposedPressure ( unsigned int cond, Real p,Solver& flow )
-{
- if ( cond>=flow->imposedP.size() ) LOG_ERROR ( "Setting p with cond higher than imposedP size." );
- flow->imposedP[cond].second=p;
- //force immediate update of boundary conditions
- flow->pressureChanged=true;
-}
-
-template<class Solver>
-void FlowEngine::imposeFlux ( Vector3r pos, Real flux,Solver& flow ){
- flow.imposedF.push_back ( pair<CGT::Point,Real> ( CGT::Point ( pos[0],pos[1],pos[2] ), flux ) );
-}
-
-template<class Solver>
-void FlowEngine::clearImposedPressure ( Solver& flow ) { flow->imposedP.clear(); flow->ipCells.clear();}
-template<class Solver>
-void FlowEngine::clearImposedFlux ( Solver& flow ) { flow->imposedF.clear(); flow->ifCells.clear();}
-
-template<class Solver>
-Real FlowEngine::getCellFlux ( unsigned int cond, const shared_ptr<Solver>& flow )
-{
- if ( cond>=flow->imposedP.size() ) {LOG_ERROR ( "Getting flux with cond higher than imposedP size." ); return 0;}
- double flux=0;
- typename Solver::CellHandle& cell= flow->ipCells[cond];
- for ( int ngb=0;ngb<4;ngb++ ) {
- flux+= cell->info().kNorm() [ngb]* ( cell->info().p()-cell->neighbor ( ngb )->info().p() );
- }
- return flux+cell->info().dv();
-}
-
-template<class Solver>
-void FlowEngine::initSolver ( Solver& flow )
-{
- flow->Vtotalissimo=0; flow->VSolidTot=0; flow->vPoral=0; flow->sSolidTot=0;
- flow->slipOnLaterals = 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 Solver>
-void FlowEngine::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 Solver>
-bool FlowEngine::getForceMetis ( Solver& flow ) {return (flow->eSolver.cholmod().nmethods==1);}
-#endif
-
-template<class Solver>
-void FlowEngine::buildTriangulation ( Solver& flow )
-{
- buildTriangulation ( 0.f,flow );
-}
-
-template<class Solver>
-void FlowEngine::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 cellEnd = flow->T[flow->currentTes].Triangulation().finite_cells_end();
- int k=0;
- for ( FiniteCellsIterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cellEnd; 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();
- porosity = flow->vPoralPorosity/flow->vTotalePorosity;
-
- boundaryConditions ( flow );
- flow->initializePressures ( 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();
-}
-
-void FlowEngine::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 Solver>
-void FlowEngine::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 id_offset 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 Solver>
-void FlowEngine::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 FlowEngine, 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 Solver>
-void FlowEngine::initializeVolumes ( Solver& flow )
-{
- typedef typename Solver::element_type Flow;
- typedef typename Flow::FiniteVerticesIterator FiniteVerticesIterator;
- typedef typename Solver::element_type Flow;
-
- FiniteVerticesIterator verticesEnd = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
- CGT::CVector Zero(0,0,0);
- for (FiniteVerticesIterator vIt = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); vIt!= verticesEnd; vIt++) vIt->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;
-}
-
-void FlowEngine::averageRealCellVelocity()
-{
- solver->averageRelativeCellVelocity();
- Vector3r Vel ( 0,0,0 );
- //AVERAGE CELL VELOCITY
- FiniteCellsIterator cellEnd = solver->T[solver->currentTes].Triangulation().finite_cells_end();
- for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cellEnd; 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 posAvFacet;
- double volumeFacetTranslation = 0;
- CGT::CVector velAv ( Vel[0], Vel[1], Vel[2] );
- for ( int i=0; i<4; i++ ) {
- volumeFacetTranslation = 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();
- posAvFacet = ( CGT::Point ) cell->info() + ( branch*Surfk ) *Surfk;
- volumeFacetTranslation += velAv*cell->info().facetSurfaces[i];
- cell->info().averageVelocity() = cell->info().averageVelocity() - volumeFacetTranslation/cell->info().volume() * ( posAvFacet-CGAL::ORIGIN );
- }
- }
- }
-}
-template<class Solver>
-void FlowEngine::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 Cellhandle>
-Real FlowEngine::volumeCellSingleFictious ( Cellhandle cell )
-{
- Vector3r V[3];
- int b=0;
- int w=0;
- cell->info().volumeSign=1;
- Real wallCoordinate=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 ) wallCoordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wallThickness/2;
- else wallCoordinate = 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] ) - wallCoordinate );
- return abs ( volume );
-}
-template<class Cellhandle>
-Real FlowEngine::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 wallCoordinate[2];
- int j=0;
- bool firstSph=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 ) wallCoordinate[j] = positionBufferCurrent[b[j]].pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2;
- else wallCoordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
- j++;
- } else if ( firstSph ) {
- A=AS=/*AT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;
- firstSph=false;
- } else {
- B=BS=/*BT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;;
- }
- }
- AS[coord[0]]=BS[coord[0]] = wallCoordinate[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]] )-wallCoordinate[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]] )-wallCoordinate[1] );
- return abs ( Vol1+Vol2 );
-}
-template<class Cellhandle>
-Real FlowEngine::volumeCellTripleFictious ( Cellhandle cell )
-{
- Vector3r A;
-
- int b[3];
- int coord[3];
- Real wallCoordinate[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 ) wallCoordinate[j] = wll->state->pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2;
- else wallCoordinate[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]]-wallCoordinate[0] ) * ( A[coord[1]]-wallCoordinate[1] ) * ( A[coord[2]]-wallCoordinate[2] );
- return abs ( volume );
-}
-template<class Cellhandle>
-Real FlowEngine::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 Solver>
-void FlowEngine::computeLubricationContributions ( 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);
- }
- }
-}
-
-YADE_PLUGIN ( ( FlowEngine ) );
+YADE_PLUGIN((FlowEngine));
//______________________________________________________________
@@ -742,13 +50,13 @@
if (first) {
if (multithread) setPositionsBuffer(false);
cachedCell= Cell(*(scene->cell));
- buildTriangulation(pZero,solver);
+ buildTriangulation(pZero,*solver);
if (solver->errorCode>0) {LOG_INFO("triangulation error, pausing"); Omega::instance().pause(); return;}
- initializeVolumes(solver); backgroundSolver=solver; backgroundCompleted=true;}
-// if ( first ) {Build_Triangulation ( pZero ); updateTriangulation = false; initializeVolumes();}
+ initializeVolumes(*solver); backgroundSolver=solver; backgroundCompleted=true;}
+// if ( first ) {buildTriangulation ( pZero ); updateTriangulation = false; initializeVolumes();}
timingDeltas->checkpoint("Triangulating");
- updateVolumes (solver);
+ updateVolumes (*solver);
epsVolCumulative += epsVolMax;
retriangulationLastIter++;
if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
@@ -756,17 +64,17 @@
timingDeltas->checkpoint("Update_Volumes");
- ///Compute flow and and forces here
+ ///compute flow and and forces here
if (pressureForce){
solver->gaussSeidel(scene->dt);
timingDeltas->checkpoint("Gauss-Seidel");
solver->computeFacetForcesWithCache();}
- timingDeltas->checkpoint("Compute_Pressure_Forces");
+ timingDeltas->checkpoint("compute_Pressure_Forces");
- ///Compute vicscous forces
+ ///compute vicscous forces
scene->forces.sync();
- computeLubricationContributions(*solver);
- timingDeltas->checkpoint("Compute_Viscous_Forces");
+ computeViscousForces(*solver);
+ timingDeltas->checkpoint("compute_Viscous_Forces");
Vector3r force;
Vector3r torque;
const Tesselation& Tes = solver->T[solver->currentTes];
@@ -809,8 +117,8 @@
epsVolCumulative=0;
boost::thread workerThread(&PeriodicFlowEngine::backgroundAction,this);
workerThread.detach();
- initializeVolumes(solver);
- computeLubricationContributions(*solver);
+ initializeVolumes(*solver);
+ computeViscousForces(*solver);
}
else if (debug && !first) {
if (debug && !backgroundCompleted) cerr<<"still computing solver in the background"<<endl;
@@ -818,9 +126,9 @@
} else {
if (updateTriangulation && !first) {
cachedCell= Cell(*(scene->cell));
- buildTriangulation (pZero, solver);
- initializeVolumes(solver);
- computeLubricationContributions(*solver);
+ buildTriangulation (pZero, *solver);
+ initializeVolumes(*solver);
+ computeViscousForces(*solver);
updateTriangulation = false;
epsVolCumulative=0;
retriangulationLastIter=0;
@@ -831,20 +139,20 @@
}
-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::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( shared_ptr<FlowSolver>& flow )
+void PeriodicFlowEngine::triangulate( FlowSolver& flow )
{
- Tesselation& Tes = flow->T[flow->currentTes];
+ 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;
@@ -859,7 +167,7 @@
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;
+ 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() );
@@ -888,13 +196,13 @@
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 );
+ 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 );
}
@@ -916,7 +224,7 @@
int b=0;
int w=0;
cell->info().volumeSign=1;
- Real wallCoordinate=0;
+ Real Wall_coordinate=0;
for ( int y=0;y<4;y++ ) {
if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
@@ -926,25 +234,25 @@
} else {
b = cell->vertex ( y )->info().id();
const shared_ptr<Body>& wll = Body::byId ( b,scene );
- if ( !solver->boundary ( b ).useMaxMin ) wallCoordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wallThickness/2;
- else wallCoordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate];
+ 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] ) - wallCoordinate );
+ 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, shared_ptr<FlowSolver>& flow, unsigned int count)
+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;
+ 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();
+ RTriangulation& Tri = flow.T[flow.currentTes].Triangulation();
Vector3r center ( 0,0,0 );
Vector3i period;
@@ -956,8 +264,8 @@
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];}
+ coord=flow.boundary ( baseCell->vertex ( k )->info().id() ).coordinate;
+ boundPos=flow.boundary ( baseCell->vertex ( k )->info().id() ).p[coord];}
}
center[coord]=boundPos;
}
@@ -1025,8 +333,8 @@
solver->averageRelativeCellVelocity();
Vector3r meanVel ( 0,0,0 );
Real volume=0;
- FiniteCellsIterator cellEnd = solver->T[solver->currentTes].Triangulation().finite_cells_end();
- for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cellEnd; cell++ ) {
+ 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;
@@ -1037,7 +345,7 @@
return ( meanVel/volume );
}
-void PeriodicFlowEngine::updateVolumes (shared_ptr<FlowSolver>& flow)
+void PeriodicFlowEngine::updateVolumes (FlowSolver& flow)
{
if ( debug ) cout << "Updating volumes.............." << endl;
Real invDeltaT = 1/scene->dt;
@@ -1048,7 +356,7 @@
Real totVol0=0;
Real totVol1=0;
- FOREACH(CellHandle& cell, flow->T[flow->currentTes].cellHandles){
+ FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
switch ( cell->info().fictious() ) {
case ( 1 ) :
newVol = volumeCellSingleFictious ( cell );
@@ -1073,13 +381,13 @@
}
-void PeriodicFlowEngine::initializeVolumes (shared_ptr<FlowSolver>& flow)
+void PeriodicFlowEngine::initializeVolumes (FlowSolver& flow)
{
- FiniteVerticesIterator verticesEnd = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
+ FiniteVerticesIterator vertices_end = flow.T[flow.currentTes].Triangulation().finite_vertices_end();
CGT::CVector Zero ( 0,0,0 );
- for ( FiniteVerticesIterator vIt = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); vIt!= verticesEnd; vIt++ ) vIt->info().forces=Zero;
+ 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){
+ FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
switch ( cell->info().fictious() )
{
case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
@@ -1087,77 +395,98 @@
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 (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, shared_ptr<FlowSolver>& flow)
+void PeriodicFlowEngine::buildTriangulation ( double pZero, FlowSolver& flow)
{
- flow->resetNetwork();
- if (first) flow->currentTes=0;
+ flow.resetNetwork();
+ if (first) flow.currentTes=0;
else {
- flow->currentTes=!flow->currentTes;
+ 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();
+ flow.T[flow.currentTes].compute();
+ flow.defineFictiousCells();
//FIXME: this is already done in addBoundary(?)
boundaryConditions ( flow );
if ( debug ) cout << endl << "boundaryConditions------" << endl << endl;
- flow->initializePressures ( pZero );
- if ( debug ) cout << endl << "initializePressures------" << 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];
+ 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;
+ 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->vTotalePorosity;
- flow->tolerance=tolerance;flow->relax=relax;
+ flow.computePermeability ( );
+ porosity = flow.vPoralPorosity/flow.vTotalPorosity;
+ flow.tolerance=tolerance;flow.relax=relax;
- flow->displayStatistics ();
+ 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 ( !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 ( waveAction ) flow.applySinusoidalPressure ( Tes.Triangulation(), sineMagnitude, sineAverage, 30 );
- if (normalLubrication || shearLubrication || viscousShear) flow->computeEdgesSurfaces();
+ if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
if ( debug ) cout << endl << "end buildTri------" << endl << endl;
}
void PeriodicFlowEngine::preparePShifts()
{
- CGT::PeriodicCellInfo::gradP = makeCgVect ( gradP );
- CGT::PeriodicCellInfo::hSize[0] = makeCgVect ( scene->cell->hSize.col ( 0 ) );
- CGT::PeriodicCellInfo::hSize[1] = makeCgVect ( scene->cell->hSize.col ( 1 ) );
- CGT::PeriodicCellInfo::hSize[2] = makeCgVect ( scene->cell->hSize.col ( 2 ) );
- CGT::PeriodicCellInfo::deltaP=CGT::CVector (
- CGT::PeriodicCellInfo::hSize[0]*CGT::PeriodicCellInfo::gradP,
- CGT::PeriodicCellInfo::hSize[1]*CGT::PeriodicCellInfo::gradP,
- CGT::PeriodicCellInfo::hSize[2]*CGT::PeriodicCellInfo::gradP );
+ 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));
+
#endif //FLOW_ENGINE
#endif /* YADE_CGAL */
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2014-03-21 18:45:24 +0000
+++ pkg/dem/FlowEngine.hpp 2014-03-21 18:47:45 +0000
@@ -9,18 +9,19 @@
/*!
-FlowEngine is mainly an interface to the fluid solvers defined in lib/triangulation and using the PFV scheme.
+FlowEngine is an interface between Yade and the fluid solvers defined in lib/triangulation and using the PFV scheme.
There are also a few non-trivial functions defined here, such has building triangulation and computating of elements volumes.
+FlowEngine is defined via the template class TemplateFlowEngine, allowing flexible definition of cell info and vertex info types for handling different types of models and couplings, via inheritance and/or template specialization.
PeriodicFlowEngine is a variant for periodic boundary conditions.
-Which solver will be actually used internally to obtain pore pressure will depend on compile time flags (libcholmod available and #define LINSOLV),
+Which solver will be actually used internally to obtain pore pressure will depend partly on compile time flags (libcholmod available and #define LINSOLV),
and on runtime settings (useSolver=0: Gauss-Seidel (iterative), useSolver=1: CHOLMOD if available (direct sparse cholesky)).
The files defining lower level classes are in yade/lib. The code uses CGAL::Triangulation3 for managing the mesh and storing data; eigen3::Sparse, suitesparse::cholmod, and metis for direct resolution of the linear systems. The Gauss-Seidel iterative solver OTOH is implemented directly in Yade.
Most classes in lib/triangulation are templates, and therefore completely defined in header files.
-A pseudo hpp/cpp split is reproduced for clarity, with hpp/ipp extensions (but again, in the end they are all included).
+A pseudo hpp/cpp split is reproduced for clarity, with hpp/ipp extensions (but again, in the end they are all include files).
The files are
- RegularTriangulation.h : declaration of info types embeded in the mesh (template parameters of CGAL::Triangulation3), and instanciation of RTriangulation classes
- Tesselation.h/ipp : encapsulate RegularTriangulation and adds functions to manipulate the dual (Voronoi) graph of the triangulation
@@ -37,31 +38,44 @@
#include<yade/lib/triangulation/FlowBoundingSphere.hpp>
#include<yade/lib/triangulation/PeriodicFlow.hpp>
-class Flow;
-class TesselationWrapper;
+/// Converters for Eigen and CGAL vectors
+CGT::CVector makeCgVect ( const Vector3r& yv ) {return CGT::CVector ( yv[0],yv[1],yv[2] );}
+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
-
-#define TPL template<class Solver>
-
-
-class FlowEngine : public PartialEngine
-{
+#endif*/
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation=CGT::_Tesselation<CGT::TriangulationTypes<_VertexInfo,_CellInfo> >, class solverT=CGT::FlowBoundingSphere<_Tesselation> >
+class TemplateFlowEngine : public PartialEngine
+{
public :
- typedef _FlowSolver FlowSolver;
- typedef FlowTesselation Tesselation;
- typedef FlowSolver::RTriangulation RTriangulation;
- typedef FlowSolver::FiniteVerticesIterator FiniteVerticesIterator;
- typedef FlowSolver::FiniteCellsIterator FiniteCellsIterator;
- typedef FlowSolver::CellHandle CellHandle;
- typedef RTriangulation::Finite_edges_iterator FiniteEdgesIterator;
- typedef FlowSolver::VertexHandle VertexHandle;
-
-
-
+ #ifdef LINSOLV
+ typedef CGT::FlowBoundingSphereLinSolv<typename solverT::Tesselation,solverT> FlowSolver;
+ #else
+ 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;
+ typedef typename FlowSolver::FiniteVerticesIterator FiniteVerticesIterator;
+ typedef typename FlowSolver::FiniteCellsIterator FiniteCellsIterator;
+ typedef typename FlowSolver::CellHandle CellHandle;
+ typedef typename FlowSolver::FiniteEdgesIterator FiniteEdgesIterator;
+ typedef typename FlowSolver::VertexHandle VertexHandle;
+ typedef typename RTriangulation::Triangulation_data_structure::Cell::Info CellInfo;
+ typedef typename RTriangulation::Triangulation_data_structure::Vertex::Info VertexInfo;
+
protected:
shared_ptr<FlowSolver> solver;
shared_ptr<FlowSolver> backgroundSolver;
@@ -75,77 +89,78 @@
public :
int retriangulationLastIter;
- enum {wallxMin, wallxMax, wallyMin, wallyMax, wallzMin, wallzMax};
+ enum {wall_xmin, wall_xmax, wall_ymin, wall_ymax, wall_zmin, wall_zmax};
Vector3r normal [6];
bool currentTes;
int idOffset;
double epsVolCumulative;
int ReTrg;
int ellapsedIter;
- TPL void initSolver (Solver& flow);
+ void initSolver (FlowSolver& flow);
#ifdef LINSOLV
- TPL void setForceMetis (Solver& flow, bool force);
- TPL bool getForceMetis (Solver& flow);
+ void setForceMetis (Solver& flow, bool force);
+ bool getForceMetis (Solver& flow);
#endif
- TPL void triangulate (Solver& flow);
- TPL void addBoundary (Solver& flow);
- TPL void buildTriangulation (double pZero, Solver& flow);
- TPL void buildTriangulation (Solver& flow);
- TPL void updateVolumes (Solver& flow);
- TPL void initializeVolumes (Solver& flow);
- TPL void boundaryConditions(Solver& flow);
- TPL 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 triangulate (Solver& flow);
+ void addBoundary (Solver& flow);
+ void buildTriangulation (double pZero, Solver& flow);
+ void buildTriangulation (Solver& flow);
+ void updateVolumes (Solver& flow);
+ void initializeVolumes (Solver& flow);
+ void boundaryConditions(Solver& flow);
+ void updateBCs ( Solver& flow ) {
+ if (flow.T[flow.currentTes].maxId>0) boundaryConditions(flow);//avoids crash at iteration 0, when the packing is not bounded yet
else LOG_ERROR("updateBCs not applied");
- flow->pressureChanged=true;}
+ flow.pressureChanged=true;}
- TPL void imposeFlux(Vector3r pos, Real flux,Solver& flow);
- TPL unsigned int imposePressure(Vector3r pos, Real p,Solver& flow);
- TPL void setImposedPressure(unsigned int cond, Real p,Solver& flow);
- TPL void clearImposedPressure(Solver& flow);
- TPL void clearImposedFlux(Solver& flow);
- TPL void computeLubricationContributions(Solver& flow);
- TPL Real getCellFlux(unsigned int cond, const shared_ptr<Solver>& flow);
- TPL Real getBoundaryFlux(unsigned int boundary,Solver& flow) {return flow->boundaryFlux(boundary);}
- TPL Vector3r fluidForce(unsigned int id_sph, Solver& flow) {
- const CGT::CVector& f=flow->T[flow->currentTes].vertex(id_sph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
- TPL Vector3r shearLubForce(unsigned int id_sph, Solver& flow) {
- return (flow->shearLubricationForces.size()>id_sph)?flow->shearLubricationForces[id_sph]:Vector3r::Zero();}
- TPL Vector3r shearLubTorque(unsigned int id_sph, Solver& flow) {
- return (flow->shearLubricationTorques.size()>id_sph)?flow->shearLubricationTorques[id_sph]:Vector3r::Zero();}
- TPL Vector3r pumpLubTorque(unsigned int id_sph, Solver& flow) {
- return (flow->pumpLubricationTorques.size()>id_sph)?flow->pumpLubricationTorques[id_sph]:Vector3r::Zero();}
- TPL Vector3r twistLubTorque(unsigned int id_sph, Solver& flow) {
- return (flow->twistLubricationTorques.size()>id_sph)?flow->twistLubricationTorques[id_sph]:Vector3r::Zero();}
- TPL Vector3r normalLubForce(unsigned int id_sph, Solver& flow) {
- return (flow->normalLubricationForce.size()>id_sph)?flow->normalLubricationForce[id_sph]:Vector3r::Zero();}
- TPL Matrix3r bodyShearLubStress(unsigned int id_sph, Solver& flow) {
- return (flow->shearLubricationBodyStress.size()>id_sph)?flow->shearLubricationBodyStress[id_sph]:Matrix3r::Zero();}
- TPL Matrix3r bodyNormalLubStress(unsigned int id_sph, Solver& flow) {
- return (flow->normalLubricationBodyStress.size()>id_sph)?flow->normalLubricationBodyStress[id_sph]:Matrix3r::Zero();}
- TPL Vector3r shearVelocity(unsigned int interaction, Solver& flow) {
- return (flow->deltaShearVel[interaction]);}
- TPL Vector3r normalVelocity(unsigned int interaction, Solver& flow) {
- return (flow->deltaNormVel[interaction]);}
- TPL Matrix3r normalStressInteraction(unsigned int interaction, Solver& flow) {
- return (flow->normalStressInteraction[interaction]);}
- TPL Matrix3r shearStressInteraction(unsigned int interaction, Solver& flow) {
- return (flow->shearStressInteraction[interaction]);}
- TPL Vector3r normalVect(unsigned int interaction, Solver& flow) {
- return (flow->normalV[interaction]);}
- TPL Real surfaceDistanceParticle(unsigned int interaction, Solver& flow) {
- return (flow->surfaceDistance[interaction]);}
- TPL Real edgeSize(Solver& flow) {
- return (flow->edgeIds.size());}
- TPL Real OSI(Solver& flow) {
- return (flow->onlySpheresInteractions.size());}
- TPL int onlySpheresInteractions(unsigned int interaction, Solver& flow) {
- return (flow->onlySpheresInteractions[interaction]);}
- TPL python::list getConstrictions(bool all, Solver& flow) {
- vector<Real> csd=flow->getConstrictions(); python::list pycsd;
+ void imposeFlux(Vector3r pos, Real flux);
+ unsigned int imposePressure(Vector3r pos, Real p);
+ void setImposedPressure(unsigned int cond, Real p);
+ void clearImposedPressure();
+ void clearImposedFlux();
+ void computeViscousForces(Solver& flow);
+ Real getCellFlux(unsigned int cond);
+ Real getBoundaryFlux(unsigned int boundary) {return solver->boundaryFlux(boundary);}
+ Vector3r fluidForce(unsigned int idSph) {
+ const CGT::CVector& f=solver->T[solver->currentTes].vertex(idSph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
+
+ Vector3r shearLubForce(unsigned int id_sph) {
+ return (solver->shearLubricationForces.size()>id_sph)?solver->shearLubricationForces[id_sph]:Vector3r::Zero();}
+ Vector3r shearLubTorque(unsigned int id_sph) {
+ return (solver->shearLubricationTorques.size()>id_sph)?solver->shearLubricationTorques[id_sph]:Vector3r::Zero();}
+ Vector3r pumpLubTorque(unsigned int id_sph) {
+ return (solver->pumpLubricationTorques.size()>id_sph)?solver->pumpLubricationTorques[id_sph]:Vector3r::Zero();}
+ Vector3r twistLubTorque(unsigned int id_sph) {
+ return (solver->twistLubricationTorques.size()>id_sph)?solver->twistLubricationTorques[id_sph]:Vector3r::Zero();}
+ Vector3r normalLubForce(unsigned int id_sph) {
+ return (solver->normalLubricationForce.size()>id_sph)?solver->normalLubricationForce[id_sph]:Vector3r::Zero();}
+ Matrix3r bodyShearLubStress(unsigned int id_sph) {
+ return (solver->shearLubricationBodyStress.size()>id_sph)?solver->shearLubricationBodyStress[id_sph]:Matrix3r::Zero();}
+ Matrix3r bodyNormalLubStress(unsigned int id_sph) {
+ return (solver->normalLubricationBodyStress.size()>id_sph)?solver->normalLubricationBodyStress[id_sph]:Matrix3r::Zero();}
+ Vector3r shearVelocity(unsigned int interaction) {
+ return (solver->deltaShearVel[interaction]);}
+ Vector3r normalVelocity(unsigned int interaction) {
+ return (solver->deltaNormVel[interaction]);}
+ Matrix3r normalStressInteraction(unsigned int interaction) {
+ return (solver->normalStressInteraction[interaction]);}
+ Matrix3r shearStressInteraction(unsigned int interaction) {
+ return (solver->shearStressInteraction[interaction]);}
+ Vector3r normalVect(unsigned int interaction) {
+ return (solver->normalV[interaction]);}
+ Real surfaceDistanceParticle(unsigned int interaction) {
+ return (solver->surfaceDistance[interaction]);}
+ Real edgeSize() {
+ return (solver->edgeIds.size());}
+ Real OSI() {
+ return (solver->onlySpheresInteractions.size());}
+ int onlySpheresInteractions(unsigned int interaction) {
+ return (solver->onlySpheresInteractions[interaction]);}
+ python::list getConstrictions(bool all) {
+ vector<Real> csd=solver->getConstrictions(); python::list pycsd;
for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k]!=0) || csd[k]>0) pycsd.append(csd[k]); return pycsd;}
- TPL python::list getConstrictionsFull(bool all, Solver& flow) {
- vector<Constriction> csd=flow->getConstrictionsFull(); python::list pycsd;
+ python::list getConstrictionsFull(bool all) {
+ vector<Constriction> csd=solver->getConstrictionsFull(); python::list pycsd;
for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k].second[0]!=0) || csd[k].second[0]>0) {
python::list cons;
cons.append(csd[k].first.first);
@@ -165,344 +180,380 @@
Real volumeCellTripleFictious (Cellhandle cell);
template<class Cellhandle>
Real volumeCell (Cellhandle cell);
+ void Oedometer_Boundary_Conditions();
void averageRealCellVelocity();
- void saveVtk() {solver->saveVtk();}
- vector<Real> avFlVelOnSph(unsigned int id_sph) {return solver->averageFluidVelocityOnSphere(id_sph);}
+ void saveVtk(const char* folder) {solver->saveVtk(folder);}
+ vector<Real> avFlVelOnSph(unsigned int idSph) {return solver->averageFluidVelocityOnSphere(idSph);}
+
+// void setBoundaryVel(Vector3r vel) {topBoundaryVelocity=vel; updateTriangulation=true;}
void pressureProfile(double wallUpY, double wallDownY) {return solver->measurePressureProfile(wallUpY,wallDownY);}
double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
- TPL int getCell(double posX, double posY, double posZ, Solver& flow){return flow->getCell(posX, posY, posZ);}
- TPL unsigned int nCells(Solver& flow){return flow->T[flow->currentTes].cellHandles.size();}
- TPL python::list getVertices(unsigned int id, Solver& flow){
+ int getCell(double posX, double posY, double posZ){return solver->getCell(posX, posY, posZ);}
+ unsigned int nCells(){return solver->T[solver->currentTes].cellHandles.size();}
+ python::list getVertices(unsigned int id){
python::list ids;
- if (id>=flow->T[flow->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<flow->T[flow->currentTes].cellHandles.size()); return ids;}
- for (unsigned int i=0;i<4;i++) ids.append(flow->T[flow->currentTes].cellHandles[id]->vertex(i)->info().id());
+ if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return ids;}
+ for (unsigned int i=0;i<4;i++) ids.append(solver->T[solver->currentTes].cellHandles[id]->vertex(i)->info().id());
return ids;
}
double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
double averagePressure(){return solver->averagePressure();}
#ifdef LINSOLV
- TPL void exportMatrix(string filename,Solver& flow) {if (useSolver>0) flow->exportMatrix(filename.c_str());
- else LOG_ERROR("available for Cholmod solver (useSolver>0)");}
- TPL void exportTriplets(string filename,Solver& flow) {if (useSolver>0) flow->exportTriplets(filename.c_str());
- else LOG_ERROR("available for Cholmod solver (useSolver>0)");}
+ void exportMatrix(string filename) {if (useSolver==3) solver->exportMatrix(filename.c_str());
+ else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
+ void exportTriplets(string filename) {if (useSolver==3) solver->exportTriplets(filename.c_str());
+ else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
#endif
void emulateAction(){
scene = Omega::instance().getScene().get();
action();}
- //Instanciation of templates for python binding
- Vector3r _shearLubForce(unsigned int id_sph) {return shearLubForce(id_sph,solver);}
- Vector3r _shearLubTorque(unsigned int id_sph) {return shearLubTorque(id_sph,solver);}
- Vector3r _pumpLubTorque(unsigned int id_sph) {return pumpLubTorque(id_sph,solver);}
- Vector3r _twistLubTorque(unsigned int id_sph) {return twistLubTorque(id_sph,solver);}
- Vector3r _normalLubForce(unsigned int id_sph) {return normalLubForce(id_sph,solver);}
- Matrix3r _bodyShearLubStress(unsigned int id_sph) {return bodyShearLubStress(id_sph,solver);}
- Matrix3r _bodyNormalLubStress(unsigned int id_sph) {return bodyNormalLubStress(id_sph,solver);}
- Vector3r _fluidForce(unsigned int id_sph) {return fluidForce(id_sph,solver);}
- Vector3r _shearVelocity(unsigned int interaction) {return shearVelocity(interaction,solver);}
- Vector3r _normalVelocity(unsigned int interaction) {return normalVelocity(interaction,solver);}
- Matrix3r _normalStressInteraction(unsigned int interaction) {return normalStressInteraction(interaction,solver);}
- Matrix3r _shearStressInteraction(unsigned int interaction) {return shearStressInteraction(interaction,solver);}
- Vector3r _normalVect(unsigned int interaction) {return normalVect(interaction,solver);}
- Real _surfaceDistanceParticle(unsigned int interaction) {return surfaceDistanceParticle(interaction,solver);}
- int _onlySpheresInteractions(unsigned int interaction) {return onlySpheresInteractions(interaction,solver);}
- void _imposeFlux(Vector3r pos, Real flux) {return imposeFlux(pos,flux,*solver);}
- unsigned int _imposePressure(Vector3r pos, Real p) {return imposePressure(pos,p,solver);}
- void _setImposedPressure(unsigned int cond, Real p) {setImposedPressure(cond,p,solver);}
- void _clearImposedPressure() {clearImposedPressure(solver);}
- void _clearImposedFlux() {clearImposedFlux(solver);}
- void _updateBCs() {updateBCs(solver);}
- Real _getCellFlux(unsigned int cond) {return getCellFlux(cond,solver);}
- Real _getBoundaryFlux(unsigned int boundary) {return getBoundaryFlux(boundary,solver);}
- int _getCell(Vector3r pos) {return getCell(pos[0],pos[1],pos[2],solver);}
- unsigned int _nCells() {return nCells(solver);}
- python::list _getVertices(unsigned int id) {return getVertices(id,solver);}
#ifdef LINSOLV
- void _exportMatrix(string filename) {exportMatrix(filename,solver);}
- void _exportTriplets(string filename) {exportTriplets(filename,solver);}
- void _setForceMetis (bool force) {setForceMetis(solver,force);}
- bool _getForceMetis () {return getForceMetis (solver);}
void cholmodStats() {
cerr << cholmod_print_common("PFV Cholmod factorization",&(solver->eSolver.cholmod()))<<endl;
cerr << "cholmod method:" << solver->eSolver.cholmod().selected<<endl;
cerr << "METIS called:"<<solver->eSolver.cholmod().called_nd<<endl;}
bool metisUsed() {return bool(solver->eSolver.cholmod().called_nd);}
#endif
- python::list _getConstrictions(bool all) {return getConstrictions(all,solver);}
- python::list _getConstrictionsFull(bool all) {return getConstrictionsFull(all,solver);}
-
- virtual ~FlowEngine();
-
+
+ virtual ~TemplateFlowEngine();
virtual void action();
virtual void backgroundAction();
-
- YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(FlowEngine,PartialEngine,"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.",
- ((bool,isActivated,true,,"Activates Flow Engine"))
- ((bool,first,true,,"Controls the initialization/update phases"))
- ((double, fluidBulkModulus, 0.,,"Bulk modulus of fluid (inverse of compressibility) K=-dP*V/dV [Pa]. Flow is compressible if fluidBulkModulus > 0, else incompressible."))
- ((Real, dt, 0,,"timestep [s]"))
- ((bool,permeabilityMap,false,,"Enable/disable stocking of average permeability scalar in cell infos."))
- ((bool, slipBoundary, true,, "Controls friction condition on lateral walls"))
- ((bool,waveAction, false,, "Allow sinusoidal pressure condition to simulate ocean waves"))
- ((double, sineMagnitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied (p )"))
- ((double, sineAverage, 0,,"Pressure value (average) when sinusoidal pressure is applied"))
- ((bool, debug, false,,"Activate debug messages"))
- ((double, wallThickness,0.001,,"Walls thickness"))
- ((double,pZero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model."))
- ((double,tolerance,1e-06,,"Gauss-Seidel Tolerance"))
- ((double,relax,1.9,,"Gauss-Seidel relaxation"))
- ((bool, updateTriangulation, 0,,"If true the medium is retriangulated. Can be switched on to force retriangulation after some events (else it will be true periodicaly based on :yref:`FlowEngine::defTolerance` and :yref:`FlowEngine::meshUpdateInterval`. Of course, it costs CPU time."))
- ((int,meshUpdateInterval,1000,,"Maximum number of timesteps between re-triangulation events. See also :yref:`FlowEngine::defTolerance`."))
- ((double, epsVolMax, 0,(Attr::readonly),"Maximal absolute volumetric strain computed at each iteration. |yupdate|"))
- ((double, defTolerance,0.05,,"Cumulated deformation threshold for which retriangulation of pore space is performed. If negative, the triangulation update will occure with a fixed frequency on the basis of :yref:`FlowEngine::meshUpdateInterval`"))
- ((double, porosity, 0,(Attr::readonly),"Porosity computed at each retriangulation |yupdate|"))
- ((bool,meanKStat,false,,"report the local permeabilities' correction"))
- ((bool,clampKValues,true,,"If true, clamp local permeabilities in [minKdivKmean,maxKdivKmean]*globalK. This clamping can avoid singular values in the permeability matrix and may reduce numerical errors in the solve phase. It will also hide junk values if they exist, or bias all values in very heterogeneous problems. So, use this with care."))
- ((Real,minKdivKmean,0.0001,,"define the min K value (see :yref:`FlowEngine::clampKValues`)"))
- ((Real,maxKdivKmean,100,,"define the max K value (see :yref:`FlowEngine::clampKValues`)"))
- ((double,permeabilityFactor,1.0,,"permability multiplier"))
- ((double,viscosity,1.0,,"viscosity of the fluid"))
- ((double,stiffness, 10000,,"equivalent contact stiffness used in the lubrication model"))
- ((int, useSolver, 0,, "Solver to use: 0=Gauss-Seidel (iterative), 1=CHOLMOD (direct sparse cholesky)"))
- ((int, xmin,0,(Attr::readonly),"Index of the boundary $x_{min}$. This index is not equal the the id of the corresponding body in general, it may be used to access the corresponding attributes (e.g. flow.bndCondValue[flow.xmin], flow.wallId[flow.xmin],...)."))
- ((int, xmax,1,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
- ((int, ymin,2,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
- ((int, ymax,3,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
- ((int, zmin,4,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
- ((int, zmax,5,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
-
- ((vector<bool>, bndCondIsPressure, vector<bool>(6,false),,"defines the type of boundary condition for each side. True if pressure is imposed, False for no-flux. Indexes can be retrieved with :yref:`FlowEngine::xmin` and friends."))
- ((vector<double>, bndCondValue, vector<double>(6,0),,"Imposed value of a boundary condition. Only applies if the boundary condition is imposed pressure, else the imposed flux is always zero presently (may be generalized to non-zero imposed fluxes in the future)."))
- //FIXME: to be implemented:
- ((vector<Vector3r>, boundaryVelocity, vector<Vector3r>(6,Vector3r::Zero()),, "velocity on top boundary, only change it using :yref:`FlowEngine::setBoundaryVel`"))
- ((int, ignoredBody,-1,,"Id of a sphere to exclude from the triangulation.)"))
- ((vector<int>, wallIds,vector<int>(6),,"body ids of the boundaries (default values are ok only if aabbWalls are appended before spheres, i.e. numbered 0,...,5)"))
- ((vector<bool>, boundaryUseMaxMin, vector<bool>(6,true),,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
- ((bool, viscousShear, false,,"Compute viscous shear terms as developped by Donia Marzougui (FIXME: ref.)"))
- ((bool, shearLubrication, false,,"Compute shear lubrication force as developped by Brule (FIXME: ref.) "))
- ((bool, pumpTorque, false,,"Compute pump torque applied on particles "))
- ((bool, twistTorque, false,,"Compute twist torque applied on particles "))
- ((double, eps, 0.00001,,"roughness defined as a fraction of particles size, giving the minimum distance between particles in the lubrication model."))
- ((bool, pressureForce, true,,"Compute the pressure field and associated fluid forces. WARNING: turning off means fluid flow is not computed at all."))
- ((bool, normalLubrication, false,,"Compute normal lubrication force as developped by Brule"))
- ((bool, viscousNormalBodyStress, false,,"Compute normal viscous stress applied on each body"))
- ((bool, viscousShearBodyStress, false,,"Compute shear viscous stress applied on each body"))
- ((bool, multithread, false,,"Build triangulation and factorize in the background (multi-thread mode)"))
- #ifdef EIGENSPARSE_LIB
- ((int, numSolveThreads, 1,,"number of openblas threads in the solve phase."))
- ((int, numFactorizeThreads, 1,,"number of openblas threads in the factorization phase"))
- #endif
- ,
- /*deprec*/
- ((meanK_opt,clampKValues,"the name changed"))
- ,,
- timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
- for (int i=0; i<6; ++i){normal[i]=Vector3r::Zero(); wallIds[i]=i;}
- normal[wallyMin].y()=normal[wallxMin].x()=normal[wallzMin].z()=1;
- normal[wallyMax].y()=normal[wallxMax].x()=normal[wallzMax].z()=-1;
- solver = shared_ptr<FlowSolver> (new FlowSolver);
- first=true;
- epsVolMax=epsVolCumulative=retriangulationLastIter=0;
- ReTrg=1;
- backgroundCompleted=true;
- ellapsedIter=0;
- ,
- .def("imposeFlux",&FlowEngine::_imposeFlux,(python::arg("pos"),python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
- .def("imposePressure",&FlowEngine::_imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
- .def("setImposedPressure",&FlowEngine::_setImposedPressure,(python::arg("cond"),python::arg("p")),"Set pressure value at the point indexed 'cond'.")
- .def("clearImposedPressure",&FlowEngine::_clearImposedPressure,"Clear the list of points with pressure imposed.")
- .def("clearImposedFlux",&FlowEngine::_clearImposedFlux,"Clear the list of points with flux imposed.")
- .def("getCellFlux",&FlowEngine::_getCellFlux,(python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
- .def("getBoundaryFlux",&FlowEngine::_getBoundaryFlux,(python::arg("boundary")),"Get total flux through boundary defined by its body id.\n\n.. note:: The flux may be not zero even for no-flow condition. This artifact comes from cells which are incident to two or more boundaries (along the edges of the sample, typically). Such flux evaluation on impermeable boundary is just irrelevant, it does not imply that the boundary condition is not applied properly.")
- .def("getConstrictions",&FlowEngine::_getConstrictions,(python::arg("all")=true),"Get the list of constriction radii (inscribed circle) for all finite facets (if all==True) or all facets not incident to a virtual bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres.")
- .def("getConstrictionsFull",&FlowEngine::_getConstrictionsFull,(python::arg("all")=true),"Get the list of constrictions (inscribed circle) for all finite facets (if all==True), or all facets not incident to a fictious bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres. The constrictions are returned in the format {{cell1,cell2}{rad,nx,ny,nz}}")
- .def("saveVtk",&FlowEngine::saveVtk,(python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
- .def("avFlVelOnSph",&FlowEngine::avFlVelOnSph,(python::arg("Id_sph")),"Compute a sphere-centered average fluid velocity")
- .def("fluidForce",&FlowEngine::_fluidForce,(python::arg("Id_sph")),"Return the fluid force on sphere Id_sph.")
- .def("shearLubForce",&FlowEngine::_shearLubForce,(python::arg("Id_sph")),"Return the shear lubrication force on sphere Id_sph.")
- .def("shearLubTorque",&FlowEngine::_shearLubTorque,(python::arg("Id_sph")),"Return the shear lubrication torque on sphere Id_sph.")
- .def("pumpLubTorque",&FlowEngine::_pumpLubTorque,(python::arg("Id_sph")),"Return the pump torque on sphere Id_sph.")
- .def("twistLubTorque",&FlowEngine::_twistLubTorque,(python::arg("Id_sph")),"Return the twist torque on sphere Id_sph.")
- .def("normalLubForce",&FlowEngine::_normalLubForce,(python::arg("Id_sph")),"Return the normal lubrication force on sphere Id_sph.")
- .def("bodyShearLubStress",&FlowEngine::_bodyShearLubStress,(python::arg("Id_sph")),"Return the shear lubrication stress on sphere Id_sph.")
- .def("bodyNormalLubStress",&FlowEngine::_bodyNormalLubStress,(python::arg("Id_sph")),"Return the normal lubrication stress on sphere Id_sph.")
- .def("shearVelocity",&FlowEngine::_shearVelocity,(python::arg("id_sph")),"Return the shear velocity of the interaction.")
- .def("normalVelocity",&FlowEngine::_normalVelocity,(python::arg("id_sph")),"Return the normal velocity of the interaction.")
- .def("normalVect",&FlowEngine::_normalVect,(python::arg("id_sph")),"Return the normal vector between particles.")
- .def("surfaceDistanceParticle",&FlowEngine::_surfaceDistanceParticle,(python::arg("interaction")),"Return the distance between particles.")
- .def("onlySpheresInteractions",&FlowEngine::_onlySpheresInteractions,(python::arg("interaction")),"Return the id of the interaction only between spheres.")
-
-
- .def("pressureProfile",&FlowEngine::pressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
- .def("getPorePressure",&FlowEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
- .def("averageSlicePressure",&FlowEngine::averageSlicePressure,(python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
- .def("averagePressure",&FlowEngine::averagePressure,"Measure averaged pore pressure in the entire volume")
- .def("updateBCs",&FlowEngine::_updateBCs,"tells the engine to update it's boundary conditions before running (especially useful when changing boundary pressure - should not be needed for point-wise imposed pressure)")
- .def("emulateAction",&FlowEngine::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
- .def("getCell",&FlowEngine::_getCell,(python::arg("pos")),"get id of the cell containing (X,Y,Z).")
- .def("nCells",&FlowEngine::_nCells,"get the total number of finite cells in the triangulation.")
- .def("getVertices",&FlowEngine::_getVertices,(python::arg("id")),"get the vertices of a cell")
- #ifdef LINSOLV
- .def("exportMatrix",&FlowEngine::_exportMatrix,(python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
- .def("exportTriplets",&FlowEngine::_exportTriplets,(python::arg("filename")="triplets"),"Export system matrix to a file with only non-zero entries.")
- .def("cholmodStats",&FlowEngine::cholmodStats,"get statistics of cholmod solver activity")
- .def("metisUsed",&FlowEngine::metisUsed,"check wether metis lib is effectively used")
- .add_property("forceMetis",&FlowEngine::_getForceMetis,&FlowEngine::_setForceMetis,"If true, METIS is used for matrix preconditioning, else Cholmod is free to choose the best method (which may be METIS to, depending on the matrix). See ``nmethods`` in Cholmod documentation")
- #endif
- )
- DECLARE_LOGGER;
-};
-
-
-template<class Solver>
-unsigned int FlowEngine::imposePressure(Vector3r pos, Real p,Solver& flow)
-{
- if (!flow) LOG_ERROR("no flow defined yet, run at least one iter");
- flow->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 flow->imposedP.size()-1;
-}
-
-
-
+
+ //commodities
+ void compTessVolumes() {
+ solver->T[solver->currentTes].compute();
+ solver->T[solver->currentTes].computeVolumes();
+ }
+ Real getVolume (Body::id_t id) {
+ if (solver->T[solver->currentTes].Max_id() <= 0) {emulateAction(); LOG_WARN("Not triangulated yet, emulating action");}
+ if (solver->T[solver->currentTes].Volume(id) == -1) {compTessVolumes(); LOG_WARN("Computing all volumes now, as you did not request it explicitely.");}
+ return (solver->T[solver->currentTes].Max_id() >= id) ? solver->T[solver->currentTes].Volume(id) : -1;}
+
+ YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TemplateFlowEngine,PartialEngine,"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.",
+ ((bool,isActivated,true,,"Activates Flow Engine"))
+ ((bool,first,true,,"Controls the initialization/update phases"))
+ ((double, fluidBulkModulus, 0.,,"Bulk modulus of fluid (inverse of compressibility) K=-dP*V/dV [Pa]. Flow is compressible if fluidBulkModulus > 0, else incompressible."))
+ ((Real, dt, 0,,"timestep [s]"))
+ ((bool,permeabilityMap,false,,"Enable/disable stocking of average permeability scalar in cell infos."))
+ ((bool, slipBoundary, true,, "Controls friction condition on lateral walls"))
+ ((bool,waveAction, false,, "Allow sinusoidal pressure condition to simulate ocean waves"))
+ ((double, sineMagnitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied (p )"))
+ ((double, sineAverage, 0,,"Pressure value (average) when sinusoidal pressure is applied"))
+ ((bool, debug, false,,"Activate debug messages"))
+ ((double, wallThickness,0.001,,"Walls thickness"))
+ ((double,pZero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model."))
+ ((double,tolerance,1e-06,,"Gauss-Seidel tolerance"))
+ ((double,relax,1.9,,"Gauss-Seidel relaxation"))
+ ((bool, updateTriangulation, 0,,"If true the medium is retriangulated. Can be switched on to force retriangulation after some events (else it will be true periodicaly based on :yref:`FlowEngine::defTolerance` and :yref:`FlowEngine::meshUpdateInterval`. Of course, it costs CPU time."))
+ ((int,meshUpdateInterval,1000,,"Maximum number of timesteps between re-triangulation events. See also :yref:`FlowEngine::defTolerance`."))
+ ((double, epsVolMax, 0,(Attr::readonly),"Maximal absolute volumetric strain computed at each iteration. |yupdate|"))
+ ((double, defTolerance,0.05,,"Cumulated deformation threshold for which retriangulation of pore space is performed. If negative, the triangulation update will occure with a fixed frequency on the basis of :yref:`FlowEngine::meshUpdateInterval`"))
+ ((double, porosity, 0,(Attr::readonly),"Porosity computed at each retriangulation |yupdate|"))
+ ((bool,meanKStat,false,,"report the local permeabilities' correction"))
+ ((bool,clampKValues,true,,"If true, clamp local permeabilities in [minKdivKmean,maxKdivKmean]*globalK. This clamping can avoid singular values in the permeability matrix and may reduce numerical errors in the solve phase. It will also hide junk values if they exist, or bias all values in very heterogeneous problems. So, use this with care."))
+ ((Real,minKdivKmean,0.0001,,"define the min K value (see :yref:`FlowEngine::clampKValues`)"))
+ ((Real,maxKdivKmean,100,,"define the max K value (see :yref:`FlowEngine::clampKValues`)"))
+ ((double,permeabilityFactor,1.0,,"permability multiplier"))
+ ((double,viscosity,1.0,,"viscosity of the fluid"))
+ ((double,stiffness, 10000,,"equivalent contact stiffness used in the lubrication model"))
+ ((int, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso, 3-CHOLMOD"))
+ ((int, xmin,0,(Attr::readonly),"Index of the boundary $x_{min}$. This index is not equal the the id of the corresponding body in general, it may be used to access the corresponding attributes (e.g. flow.bndCondValue[flow.xmin], flow.wallId[flow.xmin],...)."))
+ ((int, xmax,1,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, ymin,2,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, ymax,3,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, zmin,4,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, zmax,5,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+
+ ((vector<bool>, bndCondIsPressure, vector<bool>(6,false),,"defines the type of boundary condition for each side. True if pressure is imposed, False for no-flux. Indexes can be retrieved with :yref:`FlowEngine::xmin` and friends."))
+ ((vector<double>, bndCondValue, vector<double>(6,0),,"Imposed value of a boundary condition. Only applies if the boundary condition is imposed pressure, else the imposed flux is always zero presently (may be generalized to non-zero imposed fluxes in the future)."))
+ //FIXME: to be implemented:
+ ((vector<Vector3r>, boundaryVelocity, vector<Vector3r>(6,Vector3r::Zero()),, "velocity on top boundary, only change it using :yref:`FlowEngine::setBoundaryVel`"))
+ ((int, ignoredBody,-1,,"Id of a sphere to exclude from the triangulation.)"))
+ ((vector<int>, wallIds,vector<int>(6),,"body ids of the boundaries (default values are ok only if aabbWalls are appended before spheres, i.e. numbered 0,...,5)"))
+ ((vector<bool>, boundaryUseMaxMin, vector<bool>(6,true),,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
+// ((bool, display_force, false,,"display the lubrication force applied on particles"))
+// ((bool, create_file, false,,"create file of velocities"))
+ ((bool, viscousShear, false,,"compute viscous shear terms as developped by Donia Marzougui (FIXME: ref.)"))
+ ((bool, shearLubrication, false,,"compute shear lubrication force as developped by Brule (FIXME: ref.) "))
+ ((bool, pumpTorque, false,,"Compute pump torque applied on particles "))
+ ((bool, twistTorque, false,,"Compute twist torque applied on particles "))
+ ((double, eps, 0.00001,,"roughness defined as a fraction of particles size, giving the minimum distance between particles in the lubrication model."))
+ ((bool, pressureForce, true,,"compute the pressure field and associated fluid forces. WARNING: turning off means fluid flow is not computed at all."))
+ ((bool, normalLubrication, false,,"compute normal lubrication force as developped by Brule"))
+ ((bool, viscousNormalBodyStress, false,,"compute normal viscous stress applied on each body"))
+ ((bool, viscousShearBodyStress, false,,"compute shear viscous stress applied on each body"))
+ ((bool, multithread, false,,"Build triangulation and factorize in the background (multi-thread mode)"))
+ #ifdef EIGENSPARSE_LIB
+ ((int, numSolveThreads, 1,,"number of openblas threads in the solve phase."))
+ ((int, numFactorizeThreads, 1,,"number of openblas threads in the factorization phase"))
+ #endif
+ ,
+ /*deprec*/
+ ((meanK_opt,clampKValues,"the name changed"))
+ ,,
+ timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
+ for (int i=0; i<6; ++i){normal[i]=Vector3r::Zero(); wallIds[i]=i;}
+ normal[wall_ymin].y()=normal[wall_xmin].x()=normal[wall_zmin].z()=1;
+ normal[wall_ymax].y()=normal[wall_xmax].x()=normal[wall_zmax].z()=-1;
+ solver = shared_ptr<FlowSolver> (new FlowSolver);
+ first=true;
+ epsVolMax=epsVolCumulative=retriangulationLastIter=0;
+ ReTrg=1;
+ backgroundCompleted=true;
+ ellapsedIter=0;
+ ,
+ .def("imposeFlux",&TemplateFlowEngine::imposeFlux,(python::arg("pos"),python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
+ .def("imposePressure",&TemplateFlowEngine::imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
+ .def("setImposedPressure",&TemplateFlowEngine::setImposedPressure,(python::arg("cond"),python::arg("p")),"Set pressure value at the point indexed 'cond'.")
+ .def("clearImposedPressure",&TemplateFlowEngine::clearImposedPressure,"Clear the list of points with pressure imposed.")
+ .def("clearImposedFlux",&TemplateFlowEngine::clearImposedFlux,"Clear the list of points with flux imposed.")
+ .def("getCellFlux",&TemplateFlowEngine::getCellFlux,(python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
+ .def("getBoundaryFlux",&TemplateFlowEngine::getBoundaryFlux,(python::arg("boundary")),"Get total flux through boundary defined by its body id.\n\n.. note:: The flux may be not zero even for no-flow condition. This artifact comes from cells which are incident to two or more boundaries (along the edges of the sample, typically). Such flux evaluation on impermeable boundary is just irrelevant, it does not imply that the boundary condition is not applied properly.")
+ .def("getConstrictions",&TemplateFlowEngine::getConstrictions,(python::arg("all")=true),"Get the list of constriction radii (inscribed circle) for all finite facets (if all==True) or all facets not incident to a virtual bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres.")
+ .def("getConstrictionsFull",&TemplateFlowEngine::getConstrictionsFull,(python::arg("all")=true),"Get the list of constrictions (inscribed circle) for all finite facets (if all==True), or all facets not incident to a fictious bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres. The constrictions are returned in the format {{cell1,cell2}{rad,nx,ny,nz}}")
+ .def("edgeSize",&TemplateFlowEngine::edgeSize,"Return the number of interactions.")
+ .def("OSI",&TemplateFlowEngine::OSI,"Return the number of interactions only between spheres.")
+
+ .def("saveVtk",&TemplateFlowEngine::saveVtk,(python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
+ .def("avFlVelOnSph",&TemplateFlowEngine::avFlVelOnSph,(python::arg("idSph")),"compute a sphere-centered average fluid velocity")
+ .def("fluidForce",&TemplateFlowEngine::fluidForce,(python::arg("idSph")),"Return the fluid force on sphere idSph.")
+ .def("shearLubForce",&TemplateFlowEngine::shearLubForce,(python::arg("idSph")),"Return the shear lubrication force on sphere idSph.")
+ .def("shearLubTorque",&TemplateFlowEngine::shearLubTorque,(python::arg("idSph")),"Return the shear lubrication torque on sphere idSph.")
+ .def("normalLubForce",&TemplateFlowEngine::normalLubForce,(python::arg("idSph")),"Return the normal lubrication force on sphere idSph.")
+ .def("bodyShearLubStress",&TemplateFlowEngine::bodyShearLubStress,(python::arg("idSph")),"Return the shear lubrication stress on sphere idSph.")
+ .def("bodyNormalLubStress",&TemplateFlowEngine::bodyNormalLubStress,(python::arg("idSph")),"Return the normal lubrication stress on sphere idSph.")
+ .def("shearVelocity",&TemplateFlowEngine::shearVelocity,(python::arg("idSph")),"Return the shear velocity of the interaction.")
+ .def("normalVelocity",&TemplateFlowEngine::normalVelocity,(python::arg("idSph")),"Return the normal velocity of the interaction.")
+ .def("normalVect",&TemplateFlowEngine::normalVect,(python::arg("idSph")),"Return the normal vector between particles.")
+ .def("surfaceDistanceParticle",&TemplateFlowEngine::surfaceDistanceParticle,(python::arg("interaction")),"Return the distance between particles.")
+ .def("onlySpheresInteractions",&TemplateFlowEngine::onlySpheresInteractions,(python::arg("interaction")),"Return the id of the interaction only between spheres.")
+
+
+ .def("pressureProfile",&TemplateFlowEngine::pressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
+ .def("getPorePressure",&TemplateFlowEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
+ .def("averageSlicePressure",&TemplateFlowEngine::averageSlicePressure,(python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
+ .def("averagePressure",&TemplateFlowEngine::averagePressure,"Measure averaged pore pressure in the entire volume")
+ .def("updateBCs",&TemplateFlowEngine::updateBCs,"tells the engine to update it's boundary conditions before running (especially useful when changing boundary pressure - should not be needed for point-wise imposed pressure)")
+ .def("emulateAction",&TemplateFlowEngine::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
+ .def("getCell",&TemplateFlowEngine::getCell,(python::arg("pos")),"get id of the cell containing (X,Y,Z).")
+ .def("nCells",&TemplateFlowEngine::nCells,"get the total number of finite cells in the triangulation.")
+ .def("getVertices",&TemplateFlowEngine::getVertices,(python::arg("id")),"get the vertices of a cell")
+ #ifdef LINSOLV
+ .def("exportMatrix",&TemplateFlowEngine::exportMatrix,(python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
+ .def("exportTriplets",&TemplateFlowEngine::exportTriplets,(python::arg("filename")="triplets"),"Export system matrix to a file with only non-zero entries.")
+ .def("cholmodStats",&TemplateFlowEngine::cholmodStats,"get statistics of cholmod solver activity")
+ .def("metisUsed",&TemplateFlowEngine::metisUsed,"check wether metis lib is effectively used")
+ .add_property("forceMetis",&TemplateFlowEngine::getForceMetis,&TemplateFlowEngine::setForceMetis,"If true, METIS is used for matrix preconditioning, else Cholmod is free to choose the best method (which may be METIS to, depending on the matrix). See ``nmethods`` in Cholmod documentation")
+ #endif
+ .def("compTessVolumes",&TemplateFlowEngine::compTessVolumes,"Like TesselationWrapper::computeVolumes()")
+ .def("volume",&TemplateFlowEngine::getVolume,(python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
+ )
+};
+// Definition of functions in a separate file for clarity
+#include<yade/pkg/dem/FlowEngine.ipp>
+
+// using namespace CGT;
+typedef CGT::CVector CVector;
+typedef CGT::Point Point;
+
+class FlowCellInfo : public CGT::SimpleCellInfo {
+ public:
+ //For vector storage of all cells
+ unsigned int index;
+ int volumeSign;
+ bool Pcondition;
+ Real invVoidV;
+ Real t;
+ int fict;
+ Real volumeVariation;
+ double pression;
+ //average relative (fluid - facet translation) velocity defined for a single cell as 1/Volume * SUM_ON_FACETS(x_average_facet*average_facet_flow_rate)
+ CVector averageCellVelocity;
+ // Surface vectors of facets, pointing from outside toward inside the cell
+ std::vector<CVector> facetSurfaces;
+ //Ratio between fluid surface and facet surface
+ std::vector<Real> facetFluidSurfacesRatio;
+ // Reflects the geometrical property of the cell, so that the force by cell fluid on grain "i" is pressure*unitForceVectors[i]
+ std::vector<CVector> unitForceVectors;
+ // Store the area of triangle-sphere intersections for each facet (used in forces definition)
+ std::vector<CVector> facetSphereCrossSections;
+ std::vector<CVector> cellForce;
+ std::vector<double> rayHydr;
+ std::vector<double> modulePermeability;
+ // Partial surfaces of spheres in the double-tetrahedron linking two voronoi centers. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_surfaces in the facet.
+ 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);
+ invSumK=index=volumeSign=s=volumeVariation=pression=invVoidV=fict=0;
+ isFictious=false; Pcondition = false; isGhost = false;
+ isvisited = false;
+ isGhost=false;
+ }
+ bool isGhost;
+ double invSumK;
+ bool isvisited;
+
+ 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 CGT::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;}
+};
+
+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
+#define _PeriFlowSolver CGT::PeriodicFlowLinSolv<PeriFlowTesselation>
#else
-#define _PeriFlowSolver CGT::PeriodicFlow
+#define _PeriFlowSolver CGT::PeriodicFlow<PeriFlowTesselation>
#endif
-
-class PeriodicFlowEngine : public FlowEngine
+//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
{
public :
- public :
- typedef _PeriFlowSolver FlowSolver;
- typedef PeriFlowTesselation Tesselation;
- typedef FlowSolver::RTriangulation RTriangulation;
- typedef FlowSolver::FiniteVerticesIterator FiniteVerticesIterator;
- typedef FlowSolver::FiniteCellsIterator FiniteCellsIterator;
- typedef FlowSolver::CellHandle CellHandle;
- typedef RTriangulation::Finite_edges_iterator FiniteEdgesIterator;
- typedef RTriangulation::Vertex_handle VertexHandle;
-
- shared_ptr<FlowSolver> solver;
- shared_ptr<FlowSolver> backgroundSolver;
-
- void triangulate (shared_ptr<FlowSolver>& flow);
- void buildTriangulation (Real pzero, shared_ptr<FlowSolver>& flow);
- void initializeVolumes (shared_ptr<FlowSolver>& flow);
- void updateVolumes (shared_ptr<FlowSolver>& flow);
+ 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, shared_ptr<FlowSolver>& flow, unsigned int count=0);
+ inline void locateCell(CellHandle baseCell, unsigned int& index, int& baseIndex, FlowSolver& flow, unsigned int count=0);
Vector3r meanVelocity();
virtual ~PeriodicFlowEngine();
virtual void action();
- void backgroundAction();
//Cache precomputed values for pressure shifts, based on current hSize and pGrad
void preparePShifts();
- //(re)instanciation of templates and others, for python binding
- void saveVtk(const char* folder) {solver->saveVtk(folder);}
- Vector3r _shearLubForce(unsigned int id_sph) {return shearLubForce(id_sph,solver);}
- Vector3r _shearLubTorque(unsigned int id_sph) {return shearLubTorque(id_sph,solver);}
- Vector3r _pumpLubTorque(unsigned int id_sph) {return pumpLubTorque(id_sph,solver);}
- Vector3r _twistLubTorque(unsigned int id_sph) {return twistLubTorque(id_sph,solver);}
- Vector3r _normalLubForce(unsigned int id_sph) {return normalLubForce(id_sph,solver);}
- Matrix3r _bodyShearLubStress(unsigned int id_sph) {return bodyShearLubStress(id_sph,solver);}
- Matrix3r _bodyNormalLubStress(unsigned int id_sph) {return bodyNormalLubStress(id_sph,solver);}
- Matrix3r _normalStressInteraction(unsigned int interaction) {return normalStressInteraction(interaction,solver);}
- Matrix3r _shearStressInteraction(unsigned int interaction) {return shearStressInteraction(interaction,solver);}
- Vector3r _shearVelocity(unsigned int id_sph) {return shearVelocity(id_sph,solver);}
- Vector3r _normalVelocity(unsigned int id_sph) {return normalVelocity(id_sph,solver);}
- Vector3r _normalVect(unsigned int id_sph) {return normalVect(id_sph,solver);}
- Real _surfaceDistanceParticle(unsigned int interaction) {return surfaceDistanceParticle(interaction,solver);}
- int _onlySpheresInteractions(unsigned int interaction) {return onlySpheresInteractions(interaction,solver);}
- Real _edgeSize() {return edgeSize(solver);}
- Real _OSI() {return OSI(solver);}
-
- Vector3r _fluidForce(unsigned int id_sph) {return fluidForce(id_sph, solver);}
- void _imposeFlux(Vector3r pos, Real flux) {return imposeFlux(pos,flux,*solver);}
- unsigned int _imposePressure(Vector3r pos, Real p) {return imposePressure(pos,p,this->solver);}
- Real _getBoundaryFlux(unsigned int boundary) {return getBoundaryFlux(boundary,solver);}
-
- void _updateBCs() {updateBCs(solver);}
- double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
- double averagePressure(){return solver->averagePressure();}
- void pressureProfile(double wallUpY, double wallDownY) {return solver->measurePressureProfile(wallUpY,wallDownY);}
-
- int _getCell(Vector3r pos) {return getCell(pos[0],pos[1],pos[2],solver);}
- #ifdef LINSOLV
- void _exportMatrix(string filename) {exportMatrix(filename,solver);}
- void _exportTriplets(string filename) {exportTriplets(filename,solver);}
- #endif
- Real _getCellFlux(unsigned int cond) {return getCellFlux(cond,solver);}
- python::list _getConstrictions(bool all) {return getConstrictions(all,solver);}
- python::list _getConstrictionsFull(bool all) {return getConstrictionsFull(all,solver);}
-
- //commodities
- void compTessVolumes() {
- solver->T[solver->currentTes].Compute();
- solver->T[solver->currentTes].ComputeVolumes();
- }
- Real getVolume (Body::id_t id) {
- if (solver->T[solver->currentTes].Max_id() <= 0) {emulateAction(); LOG_WARN("Not triangulated yet, emulating action");}
- if (solver->T[solver->currentTes].Volume(id) == -1) {compTessVolumes(); LOG_WARN("Computing all volumes now, as you did not request it explicitely.");}
- return (solver->T[solver->currentTes].Max_id() >= id) ? solver->T[solver->currentTes].Volume(id) : -1;}
-
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(PeriodicFlowEngine,FlowEngine,"A variant of :yref:`FlowEngine` implementing periodic boundary conditions. The API is very similar.",
- ((Real,duplicateThreshold, 0.06,,"distance from cell borders that will triger periodic duplication in the triangulation |yupdate|"))
- ((Vector3r, gradP, Vector3r::Zero(),,"Macroscopic pressure gradient"))
- ,,
- wallIds=vector<int>(6,-1);
- solver = shared_ptr<FlowSolver> (new FlowSolver);
- epsVolMax=epsVolCumulative=retriangulationLastIter=0;
- ReTrg=1;
- first=true;
- ,
- .def("meanVelocity",&PeriodicFlowEngine::meanVelocity,"measure the mean velocity in the period")
- .def("fluidForce",&PeriodicFlowEngine::_fluidForce,(python::arg("Id_sph")),"Return the fluid force on sphere Id_sph.")
- .def("shearLubForce",&PeriodicFlowEngine::_shearLubForce,(python::arg("Id_sph")),"Return the shear lubrication force on sphere Id_sph.")
- .def("shearLubTorque",&PeriodicFlowEngine::_shearLubTorque,(python::arg("Id_sph")),"Return the shear lubrication torque on sphere Id_sph.")
- .def("pumpLubTorque",&PeriodicFlowEngine::_pumpLubTorque,(python::arg("Id_sph")),"Return the pump torque on sphere Id_sph.")
- .def("twistLubTorque",&PeriodicFlowEngine::_twistLubTorque,(python::arg("Id_sph")),"Return the twist torque on sphere Id_sph.")
- .def("normalLubForce",&PeriodicFlowEngine::_normalLubForce,(python::arg("Id_sph")),"Return the normal lubrication force on sphere Id_sph.")
- .def("bodyShearLubStress",&PeriodicFlowEngine::_bodyShearLubStress,(python::arg("Id_sph")),"Return the shear lubrication stress on sphere Id_sph.")
- .def("bodyNormalLubStress",&PeriodicFlowEngine::_bodyNormalLubStress,(python::arg("Id_sph")),"Return the normal lubrication stress on sphere Id_sph.")
- .def("shearVelocity",&PeriodicFlowEngine::_shearVelocity,(python::arg("id_sph")),"Return the shear velocity of the interaction.")
- .def("normalStressInteraction",&PeriodicFlowEngine::_normalStressInteraction,(python::arg("id_sph")),"Return the shear velocity of the interaction.")
- .def("shearStressInteraction",&PeriodicFlowEngine::_shearStressInteraction,(python::arg("id_sph")),"Return the shear velocity of the interaction.")
- .def("normalVelocity",&PeriodicFlowEngine::_normalVelocity,(python::arg("id_sph")),"Return the normal velocity of the interaction.")
- .def("normalVect",&PeriodicFlowEngine::_normalVect,(python::arg("id_sph")),"Return the normal vector between particles.")
- .def("surfaceDistanceParticle",&PeriodicFlowEngine::_surfaceDistanceParticle,(python::arg("interaction")),"Return the distance between particles.")
- .def("onlySpheresInteractions",&PeriodicFlowEngine::_onlySpheresInteractions,(python::arg("interaction")),"Return the id of the interaction only between spheres.")
- .def("edgeSize",&PeriodicFlowEngine::_edgeSize,"Return the number of interactions.")
- .def("OSI",&PeriodicFlowEngine::_OSI,"Return the number of interactions only between spheres.")
-
-// .def("imposeFlux",&FlowEngine::_imposeFlux,(python::arg("pos"),python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
- .def("saveVtk",&PeriodicFlowEngine::saveVtk,(python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
- .def("imposePressure",&PeriodicFlowEngine::_imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
- .def("getBoundaryFlux",&PeriodicFlowEngine::_getBoundaryFlux,(python::arg("boundary")),"Get total flux through boundary defined by its body id.")
- .def("getPorePressure",&PeriodicFlowEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
- .def("averagePressure",&PeriodicFlowEngine::averagePressure,"Measure averaged pore pressure in the entire volume")
- .def("pressureProfile",&PeriodicFlowEngine::pressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
-
- .def("updateBCs",&PeriodicFlowEngine::_updateBCs,"tells the engine to update it's boundary conditions before running (especially useful when changing boundary pressure - should not be needed for point-wise imposed pressure)")
-
- .def("getCell",&PeriodicFlowEngine::_getCell,python::arg("pos"),"get id of the cell containing 'pos'.")
- .def("getConstrictions",&PeriodicFlowEngine::_getConstrictions,(python::arg("all")=true),"Get the list of constriction radii (inscribed circle) for all finite facets (if all==True) or all facets not incident to a virtual bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres.")
- .def("getConstrictionsFull",&PeriodicFlowEngine::_getConstrictionsFull,(python::arg("all")=true),"Get the list of constrictions (inscribed circle) for all finite facets (if all==True), or all facets not incident to a fictious bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres. The constrictions are returned in the format {{cell1,cell2}{rad,nx,ny,nz}}")
- #ifdef LINSOLV
- .def("exportMatrix",&PeriodicFlowEngine::_exportMatrix,(python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
- .def("exportTriplets",&PeriodicFlowEngine::_exportTriplets,(python::arg("filename")="triplets"),"Export system matrix to a file with only non-zero entries.")
- #endif
- .def("compTessVolumes",&PeriodicFlowEngine::compTessVolumes,"Like TesselationWrapper::computeVolumes()")
- .def("volume",&PeriodicFlowEngine::getVolume,(python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(PeriodicFlowEngine,FlowEngine_PeriodicInfo,"A variant of :yref:`FlowEngine` implementing periodic boundary conditions. The API is very similar.",
+ ((Real,duplicateThreshold, 0.06,,"distance from cell borders that will triger periodic duplication in the triangulation |yupdate|"))
+ ((Vector3r, gradP, Vector3r::Zero(),,"Macroscopic pressure gradient"))
+ ,,
+ wallIds=vector<int>(6,-1);
+ solver = shared_ptr<FlowSolver> (new FlowSolver);
+ epsVolMax=epsVolCumulative=retriangulationLastIter=0;
+ ReTrg=1;
+ first=true;
+ ,
+ //nothing special to define, we re-use FlowEngine methods
+ //.def("meanVelocity",&PeriodicFlowEngine::meanVelocity,"measure the mean velocity in the period")
)
DECLARE_LOGGER;
@@ -511,5 +562,4 @@
REGISTER_SERIALIZABLE(PeriodicFlowEngine);
-
// #endif
=== modified file 'pkg/dem/MicroMacroAnalyser.cpp'
--- pkg/dem/MicroMacroAnalyser.cpp 2014-03-21 18:45:24 +0000
+++ pkg/dem/MicroMacroAnalyser.cpp 2014-03-21 18:47:45 +0000
@@ -58,7 +58,7 @@
} else if (scene->iter % interval == 0) {
setState(2, true, compIncrt);
if (compDeformation) {
- analyser->ComputeParticlesDeformation();
+ analyser->computeParticlesDeformation();
//for (int i=0; i<analyser->ParticleDeformation.size();i++) cerr<< analyser->ParticleDeformation[i]<<endl;
ostringstream oss;
oss<<"deformation"<<incrtNumber++<<".vtk";
@@ -214,8 +214,8 @@
TS.haut = triaxialCompressionEngine->height;//find_parameter("haut=", Statefile);
TS.larg = triaxialCompressionEngine->width;//find_parameter("larg=", Statefile);
TS.prof = triaxialCompressionEngine->depth;//find_parameter("prof=", Statefile);
- TS.porom = 0/*analyser->ComputeMacroPorosity() crasher?*/;//find_parameter("porom=", Statefile);
- TS.ratio_f = triaxialCompressionEngine-> ComputeUnbalancedForce(scene); //find_parameter("ratio_f=", Statefile);
+ TS.porom = 0/*analyser->computeMacroPorosity() crasher?*/;//find_parameter("porom=", Statefile);
+ TS.ratio_f = triaxialCompressionEngine->ComputeUnbalancedForce(scene); //find_parameter("ratio_f=", Statefile);
} else TS.wszzh=TS.wsxxd=TS.wsyyfa=TS.eps3=TS.eps1=TS.eps2=TS.haut=TS.larg=TS.prof=TS.porom=TS.ratio_f=0;
if (filename!=NULL) TS.to_file(filename);
return TS;
@@ -223,7 +223,7 @@
// const vector<CGT::Tenseur3>& MicroMacroAnalyser::makeDeformationArray(const char* state_file1, const char* state_file0)
// {
-// return analyser->ComputeParticlesDeformation(state_file1, state_file0);
+// return analyser->computeParticlesDeformation(state_file1, state_file0);
// }
#endif /* YADE_CGAL */
=== modified file 'pkg/dem/MicroMacroAnalyser.hpp'
--- pkg/dem/MicroMacroAnalyser.hpp 2012-01-13 11:24:44 +0000
+++ pkg/dem/MicroMacroAnalyser.hpp 2014-03-21 18:47:45 +0000
@@ -15,7 +15,7 @@
#include <string>
#include <fstream>
-/*! \brief Compute fabric tensor, local porosity, local deformation, and other micromechanicaly defined quantities based on triangulation/tesselation of the packing.
+/*! \brief compute fabric tensor, local porosity, local deformation, and other micromechanicaly defined quantities based on triangulation/tesselation of the packing.
*/
@@ -45,7 +45,7 @@
void postLoad(MicroMacroAnalyser&);
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(MicroMacroAnalyser,GlobalEngine,"Compute fabric tensor, local porosity, local deformation, and other micromechanicaly defined quantities based on triangulation/tesselation of the packing.",
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(MicroMacroAnalyser,GlobalEngine,"compute fabric tensor, local porosity, local deformation, and other micromechanicaly defined quantities based on triangulation/tesselation of the packing.",
((unsigned int,stateNumber,0,,"A number incremented and appended at the end of output files to reflect increment number."))
((unsigned int,incrtNumber,1,,""))
((std::string,outputFile,"MicroMacroAnalysis",,"Base name for increment analysis output file."))
=== modified file 'pkg/dem/TesselationWrapper.cpp'
--- pkg/dem/TesselationWrapper.cpp 2014-03-21 18:45:24 +0000
+++ pkg/dem/TesselationWrapper.cpp 2014-03-21 18:47:45 +0000
@@ -189,26 +189,26 @@
}
}
-void TesselationWrapper::ComputeTesselation(void)
+void TesselationWrapper::computeTesselation(void)
{
if (!rad_divided) {
mean_radius /= n_spheres;
rad_divided = true;
}
- Tes->Compute();
+ Tes->compute();
}
-void TesselationWrapper::ComputeTesselation(double pminx, double pmaxx, double pminy, double pmaxy, double pminz, double pmaxz, double dt)
+void TesselationWrapper::computeTesselation(double pminx, double pmaxx, double pminy, double pmaxy, double pminz, double pmaxz, double dt)
{
addBoundingPlanes(pminx, pmaxx, pminy, pmaxy, pminz, pmaxz, dt);
- ComputeTesselation();
+ computeTesselation();
}
-void TesselationWrapper::ComputeVolumes(void)
+void TesselationWrapper::computeVolumes(void)
{
if (!bounded) addBoundingPlanes();
- ComputeTesselation();
- Tes->ComputeVolumes();
+ computeTesselation();
+ Tes->computeVolumes();
}
unsigned int TesselationWrapper::NumberOfFacets(bool initIters)
{
@@ -349,7 +349,7 @@
delete Tes;
CGT::TriaxialState* ts;
if (deformation){//use the final state to compute volumes
- /*const vector<CGT::Tenseur3>& def =*/ mma.analyser->ComputeParticlesDeformation();
+ /*const vector<CGT::Tenseur3>& def =*/ mma.analyser->computeParticlesDeformation();
Tes = &mma.analyser->TS1->tesselation();
ts = mma.analyser->TS1;
}
@@ -358,7 +358,7 @@
RTriangulation& Tri = Tes->Triangulation();
Pmin=ts->box.base; Pmax=ts->box.sommet;
//if (!scene->isPeriodic) addBoundingPlanes();
- ComputeVolumes();
+ computeVolumes();
int bodiesDim = Tes->Max_id() + 1; //=scene->bodies->size();
cerr<<"bodiesDim="<<bodiesDim<<endl;
int dim1[]={bodiesDim};
=== modified file 'pkg/dem/TesselationWrapper.hpp'
--- pkg/dem/TesselationWrapper.hpp 2014-03-21 18:45:24 +0000
+++ pkg/dem/TesselationWrapper.hpp 2014-03-21 18:47:45 +0000
@@ -71,18 +71,18 @@
/// Force boudaries at positions not equal to precomputed ones
void addBoundingPlanes(double pminx, double pmaxx, double pminy, double pmaxy, double pminz, double pmaxz, double dt);
void RemoveBoundingPlanes (void);
- ///Compute voronoi centers then stop (don't compute anything else)
- void ComputeTesselation (void);
- void ComputeTesselation( double pminx, double pmaxx, double pminy, double pmaxy, double pminz, double pmaxz, double dt);
+ ///compute voronoi centers then stop (don't compute anything else)
+ void computeTesselation (void);
+ void computeTesselation( double pminx, double pmaxx, double pminy, double pmaxy, double pminz, double pmaxz, double dt);
- ///Compute Voronoi vertices + volumes of all cells
- ///use ComputeTesselation to force update, e.g. after spheres positions have been updated
- void ComputeVolumes (void);
- void computeDeformations (void) {mma.analyser->ComputeParticlesDeformation();}
+ ///compute Voronoi vertices + volumes of all cells
+ ///use computeTesselation to force update, e.g. after spheres positions have been updated
+ void computeVolumes (void);
+ void computeDeformations (void) {mma.analyser->computeParticlesDeformation();}
///Get volume of the sphere inserted with indentifier "id""
double Volume (unsigned int id);
double deformation (unsigned int id,unsigned int i,unsigned int j) {
- if (!mma.analyser->ParticleDeformation.size()) {LOG_ERROR("Compute deformations first"); return 0;}
+ if (!mma.analyser->ParticleDeformation.size()) {LOG_ERROR("compute deformations first"); return 0;}
if (mma.analyser->ParticleDeformation.size()<id) {LOG_ERROR("id out of bounds"); return 0;}
return mma.analyser->ParticleDeformation[id](i,j);}
@@ -132,9 +132,9 @@
.def("defToVtk",&TesselationWrapper::defToVtk,(python::arg("outputFile")="def.vtk"),"Write local deformations in vtk format from states 0 and 1.")
.def("defToVtkFromStates",&TesselationWrapper::defToVtkFromStates,(python::arg("input1")="state1",python::arg("input2")="state2",python::arg("outputFile")="def.vtk",python::arg("bz2")=true),"Write local deformations in vtk format from state files (since the file format is very special, consider using defToVtkFromPositions if the input files were not generated by TesselationWrapper).")
.def("defToVtkFromPositions",&TesselationWrapper::defToVtkFromPositions,(python::arg("input1")="pos1",python::arg("input2")="pos2",python::arg("outputFile")="def.vtk",python::arg("bz2")=false),"Write local deformations in vtk format from positions files (one sphere per line, with x,y,z,rad separated by spaces).")
- .def("computeVolumes",&TesselationWrapper::ComputeVolumes,"Compute volumes of all Voronoi's cells.")
+ .def("computeVolumes",&TesselationWrapper::computeVolumes,"compute volumes of all Voronoi's cells.")
.def("getVolPoroDef",&TesselationWrapper::getVolPoroDef,(python::arg("deformation")=false),"Return a table with per-sphere computed quantities. Include deformations on the increment defined by states 0 and 1 if deformation=True (make sure to define states 0 and 1 consistently).")
- .def("ComputeDeformations",&TesselationWrapper::computeDeformations,"Compute per-particle deformation. Get it with :yref:`TesselationWrapper::deformation` (id,i,j).")
+ .def("computeDeformations",&TesselationWrapper::computeDeformations,"compute per-particle deformation. Get it with :yref:`TesselationWrapper::deformation` (id,i,j).")
.def("deformation",&TesselationWrapper::deformation,(python::arg("id"),python::arg("i"),python::arg("j")),"Get particle deformation")
);
DECLARE_LOGGER;
