yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11346
[Branch ~yade-pkg/yade/git-trunk] Rev 4181: selective blocking of cells of the mesh in FlowEngines (preliminary steps)
------------------------------------------------------------
revno: 4181
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
timestamp: Tue 2014-10-07 20:41:54 +0200
message:
selective blocking of cells of the mesh in FlowEngines (preliminary steps)
modified:
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphere.ipp
lib/triangulation/FlowBoundingSphereLinSolv.ipp
lib/triangulation/PeriodicFlow.hpp
pkg/pfv/FlowEngine.hpp
pkg/pfv/FlowEngine.hpp.in
pkg/pfv/FlowEngine.ipp.in
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2014-07-02 16:18:24 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2014-10-07 18:41:54 +0000
@@ -51,6 +51,8 @@
vector<CellHandle> IPCells;
vector<pair<Point,Real> > imposedF;
vector<CellHandle> IFCells;
+ //Blocked cells, where pressure may be computed in undrained condition
+ vector<CellHandle> blockedCells;
//Pointers to vectors used for user defined boundary pressure
vector<Real> *pxpos, *ppval;
@@ -155,6 +157,7 @@
double averagePressure();
int getCell (double X,double Y,double Z);
double boundaryFlux(unsigned int boundaryId);
+ void setBlocked(CellHandle& cell);
vector<Real> averageFluidVelocityOnSphere(unsigned int Id_sph);
//Solver?
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2014-07-16 16:49:41 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2014-10-07 18:41:54 +0000
@@ -467,6 +467,19 @@
}
template <class Tesselation>
+void FlowBoundingSphere<Tesselation>::setBlocked(CellHandle& cell)
+{
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ cerr<<"skipping a blocked cell"<<endl;
+ if (cell->info().Pcondition=true) cell->info().p() = 0;
+ else blockedCells.push_back(cell);
+ for (int j=0; j<4; j++) {
+ (cell->info().kNorm())[j]= 0;
+ (cell->neighbor(j)->info().kNorm())[Tri.mirror_index(cell, j)]= 0;}
+}
+
+
+template <class Tesselation>
void FlowBoundingSphere<Tesselation>::computePermeability()
{
if (debugOut) cout << "----Computing_Permeability------" << endl;
@@ -486,6 +499,9 @@
for (VCellIterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
CellHandle& cell = *cellIt;
+ if (cell->info().blocked) {
+ setBlocked(cell);
+ cell->info().isvisited = !ref;}
Point& p1 = cell->info();
for (int j=0; j<4; j++) {
neighbourCell = cell->neighbor(j);
@@ -504,6 +520,8 @@
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())));
+ //FIXME: it should be possible to skip completely blocked cells, currently the problem is it segfault for undefined areas
+// if (cell->info().blocked) continue;//We don't need permeability for blocked cells, it will be set to zero anyway
pass+=1;
CVector l = p1 - p2;
@@ -846,7 +864,7 @@
int bb=-1;
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
bb++;
- if ( !cell->info().Pcondition ) {
+ if ( !cell->info().Pcondition && !cell->info().blocked) {
cell2++;
#endif
if (compressible && j==0) { previousP[bb]=cell->info().p(); }
@@ -913,7 +931,6 @@
} while ((dp_max/p_max) > tolerance /*&& j<4000*/ /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
#endif
}
-
if (debugOut) {cout << "pmax " << p_max << "; pmoy : " << p_moy << endl;
cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;}
computedOnce=true;
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.ipp'
--- lib/triangulation/FlowBoundingSphereLinSolv.ipp 2014-07-02 16:18:24 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.ipp 2014-10-07 18:41:54 +0000
@@ -156,7 +156,7 @@
const FiniteCellsIterator cellEnd = Tri.finite_cells_end();
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
orderedCells.push_back(cell);
- if (!cell->info().Pcondition) ++ncols;}
+ if (!cell->info().Pcondition && !cell->info().blocked) ++ncols;}
// //Segfault on 14.10, and useless overall since SuiteSparse has preconditionners (including metis)
// spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort<RTriangulation>());
@@ -189,7 +189,7 @@
FiniteCellsIterator& cell = orderedCells[i];
///Non-ordered cells
// for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
- if (!cell->info().Pcondition) {
+ if (!cell->info().Pcondition && !cell->info().blocked) {
index=cell->info().index;
if (index==0) {
T_cells[++T_index]=cell;
@@ -325,7 +325,7 @@
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
orderedCells.push_back(cell);
- if (!cell->info().Pcondition) ++ncols;
+ if (!cell->info().Pcondition && !cell->info().blocked) ++ncols;
}
//FIXME: does it really help? test by commenting this "sorting" line
spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort<RTriangulation>());
@@ -361,7 +361,7 @@
FiniteCellsIterator& cell = orderedCells[i];
///Non-ordered cells
// for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
- if (!cell->info().Pcondition) {
+ if (!cell->info().Pcondition && !cell->info().blocked) {
if (cell->info().index==0) {
T_cells[++T_index]=cell;
cell->info().index=T_index;
@@ -403,7 +403,7 @@
FiniteCellsIterator& cell = orderedCells[i];
///Non-ordered cells
// for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
- if (!cell->info().Pcondition) for (int j=0; j<4; j++) {
+ if (!cell->info().Pcondition && !cell->info().blocked) for (int j=0; j<4; j++) {
CellHandle neighbourCell = cell->neighbor(j);
if (!Tri.is_infinite(neighbourCell) && neighbourCell->info().Pcondition)
gsB[cell->info().index]+=cell->info().kNorm()[j]*neighbourCell->info().p();
=== modified file 'lib/triangulation/PeriodicFlow.hpp'
--- lib/triangulation/PeriodicFlow.hpp 2014-09-15 10:27:11 +0000
+++ lib/triangulation/PeriodicFlow.hpp 2014-10-07 18:41:54 +0000
@@ -179,6 +179,9 @@
for (VCellIterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
CellHandle& cell = *cellIt;
Point& p1 = cell->info();
+ if (cell->info().blocked) {
+ setBlocked(cell);
+ continue;}
if (cell->info().isGhost) {cerr<<"skipping a ghost"<<endl; continue;}
for (int j=0; j<4; j++){
neighbourCell = cell->neighbor(j);
=== modified file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp 2014-06-20 15:28:21 +0000
+++ pkg/pfv/FlowEngine.hpp 2014-10-07 18:41:54 +0000
@@ -1,4 +1,5 @@
#pragma once
+
/// Frequently used:
typedef CGT::CVector CVector;
typedef CGT::Point Point;
@@ -9,6 +10,18 @@
inline Vector3r makeVector3r ( const Point& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
inline Vector3r makeVector3r ( const CVector& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
+/// The following macros can be used to expose CellInfo members.
+/// The syntax is CELL_SCALAR_GETTER(double,.p(),pressure), note the "." before member name, data members would be without the "()"
+#define CELL_SCALAR_GETTER(type, param, getterName) \
+type getterName(unsigned int 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 0;}\
+ return solver->T[solver->currentTes].cellHandles[id]->info()param;}
+
+#define CELL_VECTOR_GETTER(param, getterName) \
+Vector3r getterName(unsigned int 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 Vector3r(0,0,0);}\
+ return makeVector3r(solver->T[solver->currentTes].cellHandles[id]->info()param);}
+
#ifdef LINSOLV
#define DEFAULTSOLVER CGT::FlowBoundingSphereLinSolv<_Tesselation>
#else
=== modified file 'pkg/pfv/FlowEngine.hpp.in'
--- pkg/pfv/FlowEngine.hpp.in 2014-09-03 15:59:54 +0000
+++ pkg/pfv/FlowEngine.hpp.in 2014-10-07 18:41:54 +0000
@@ -175,12 +175,21 @@
double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
int getCell(double posX, double posY, double posZ){return solver->getCell(posX, posY, posZ);}
unsigned int nCells(){return solver->T[solver->currentTes].cellHandles.size();}
+ CELL_VECTOR_GETTER(,cellCenter)
+ CELL_VECTOR_GETTER(.averageCellVelocity,cellVelocity)
+ CELL_SCALAR_GETTER(double,.p(),pressure)
boost::python::list getVertices(unsigned int id){
boost::python::list ids;
- 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;}
+ 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;
}
+ void blockCell(unsigned int id, bool blockPressure){
+ if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return;}
+ solver->T[solver->currentTes].cellHandles[id]->info().blocked=true;
+ solver->T[solver->currentTes].cellHandles[id]->info().Pcondition=blockPressure;
+ }
+
double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
double averagePressure(){return solver->averagePressure();}
@@ -213,7 +222,7 @@
if (solver->T[solver->currentTes].Volume(id) == -1) {compTessVolumes();/* LOG_WARN("Computing all volumes now, as you did not request it explicitely.");*/}
return (solver->T[solver->currentTes].Max_id() >= id) ? solver->T[solver->currentTes].Volume(id) : -1;}
- YADE_CLASS_PYCLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TemplateFlowEngine_@TEMPLATE_FLOW_NAME@,@TEMPLATE_FLOW_NAME@,PartialEngine,"An engine to solve flow problem in saturated granular media. Model description can be found in [Chareyre2012a]_ and [Catalano2014a]_. See the example script FluidCouplingPFV/oedometer.py. More documentation to come.\n\n.. note::Multi-threading seems to work fine for Cholesky decomposition, but it fails for the solve phase in which -j1 is the fastest, here we specify thread numbers independently using :yref:`FlowEngine::numFactorizeThreads` and :yref:`FlowEngine::numSolveThreads`. These multhreading settings are only impacting the behaviour of openblas library and are relatively independant of :yref:`FlowEngine::multithread`. However, the settings have to be globally consistent. For instance, :yref:`multithread<FlowEngine::multithread>` =True with yref:`numFactorizeThreads<FlowEngine::numFactorizeThreads>` = yref:`numSolveThreads<FlowEngine::numSolveThreads>` = 4 implies that openblas will mobilize 8 processors at some point. If the system does not have so many procs. it will hurt performance.",
+ YADE_CLASS_PYCLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TemplateFlowEngine_@TEMPLATE_FLOW_NAME@,@TEMPLATE_FLOW_NAME@,PartialEngine,"A generic engine from wich more specialized engines can inherit. It is defined for the sole purpose of inserting the right data classes CellInfo and VertexInfo in the triangulation, and it should not be used directly. Instead, look for specialized engines, e.g. :yref:`FlowEngine`, :yref:`PeriodicFlowEngine`, or :yref:`DFNFlowEngine`.",
((bool,isActivated,true,,"Activates Flow Engine"))
((bool,first,true,,"Controls the initialization/update phases"))
((bool,doInterpolate,false,,"Force the interpolation of cell's info while remeshing. By default, interpolation would be done only for compressible fluids. It can be forced with this flag."))
@@ -274,6 +283,7 @@
#endif
((vector<Real>, boundaryPressure,vector<Real>(),,"values defining pressure along x-axis for the top surface. See also :yref:`@TEMPLATE_FLOW_NAME@::boundaryXPos`"))
((vector<Real>, boundaryXPos,vector<Real>(),,"values of the x-coordinate for which pressure is defined. See also :yref:`@TEMPLATE_FLOW_NAME@::boundaryPressure`"))
+ ((string,blockHook,"",,"Python command to be run after triangulation to define blocked cells (see also :yref:`TemplateFlowEngine_@TEMPLATE_FLOW_NAME@.blockCell`)"))
,
/*deprec*/
((meanK_opt,clampKValues,"the name changed"))
@@ -301,7 +311,6 @@
.def("getConstrictionsFull",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getConstrictionsFull,(boost::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_@TEMPLATE_FLOW_NAME@::edgeSize,"Return the number of interactions.")
.def("OSI",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::OSI,"Return the number of interactions only between spheres.")
-
.def("saveVtk",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
.def("avFlVelOnSph",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::avFlVelOnSph,(boost::python::arg("idSph")),"compute a sphere-centered average fluid velocity")
.def("fluidForce",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::fluidForce,(boost::python::arg("idSph")),"Return the fluid force on sphere idSph.")
@@ -322,7 +331,9 @@
.def("updateBCs",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::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_@TEMPLATE_FLOW_NAME@::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
.def("getCell",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getCell,(boost::python::arg("pos")),"get id of the cell containing (X,Y,Z).")
+ .def("getCellCenter",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::cellCenter,(boost::python::arg("id")),"get voronoi center of cell 'id'.")
.def("nCells",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::nCells,"get the total number of finite cells in the triangulation.")
+ .def("blockCell",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::blockCell,(boost::python::arg("id"),boost::python::arg("blockPressure")),"block cell 'id'. The cell will be excluded from the fluid flow problem and the conductivity of all incident facets will be null. If blockPressure=False, deformation is reflected in the pressure, else it is constantly 0.")
.def("getVertices",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getVertices,(boost::python::arg("id")),"get the vertices of a cell")
#ifdef LINSOLV
.def("exportMatrix",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::exportMatrix,(boost::python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
@@ -336,7 +347,7 @@
.def("averageVelocity",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averageVelocity,"measure the mean velocity in the period")
)
};
-// Definition of functions in a separate file for clarity
+
class FlowCellInfo_@TEMPLATE_FLOW_NAME@ : public CGT::SimpleCellInfo {
public:
@@ -344,6 +355,7 @@
unsigned int index;
int volumeSign;
bool Pcondition;
+ bool blocked;//exclude cell from the fluid domain
Real invVoidV;
Real t;
int fict;
@@ -377,8 +389,9 @@
rayHydr.resize(4, 0);
invSumK=index=volumeSign=s=volumeVariation=pression=invVoidV=fict=0;
isFictious=false; Pcondition = false; isGhost = false;
- isvisited = false;
+ isvisited = false;
isGhost=false;
+ blocked=false;
}
bool isGhost;
double invSumK;
@@ -414,4 +427,5 @@
inline const CVector ghostShift (void) const {return CGAL::NULL_VECTOR;}
};
+// Implementation of functions in this separate file for clarity
#include "FlowEngine_@TEMPLATE_FLOW_NAME@.ipp"
=== modified file 'pkg/pfv/FlowEngine.ipp.in'
--- pkg/pfv/FlowEngine.ipp.in 2014-08-01 11:16:21 +0000
+++ pkg/pfv/FlowEngine.ipp.in 2014-10-07 18:41:54 +0000
@@ -12,6 +12,8 @@
#ifdef LINSOLV
#include <cholmod.h>
#endif
+//For pyRunString
+#include<yade/lib/pyutil/gil.hpp>
template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::~TemplateFlowEngine_@TEMPLATE_FLOW_NAME@() {}
@@ -211,6 +213,7 @@
flow.fluidBulkModulus = fluidBulkModulus;
// flow.tesselation().Clear();
flow.tesselation().maxId=-1;
+ flow.blockedCells.clear();
flow.xMin = 1000.0, flow.xMax = -10000.0, flow.yMin = 1000.0, flow.yMax = -10000.0, flow.zMin = 1000.0, flow.zMax = -10000.0;
}
@@ -255,6 +258,7 @@
flow.tesselation().cellHandles.push_back(cell);
cell->info().id=k++;}//define unique numbering now, corresponds to position in cellHandles
flow.displayStatistics ();
+ if(!blockHook.empty()){ LOG_INFO("Running blockHook: "<<blockHook); pyRunString(blockHook); }
flow.computePermeability();
//This virtual function does nothing yet, derived class may overload it to make permeability different (see DFN engine)
trickPermeability();
