~ubuntu-branches/ubuntu/wily/yade/wily

template<class _CellInfo, class _VertexInfo, class _Tesselation=CGT::_Tesselation<CGT::TriangulationTypes<_VertexInfo,_CellInfo> >, class solverT=DEFAULTSOLVER >

class TemplateFlowEngine_@TEMPLATE_FLOW_NAME@ : public PartialEngine

100

void imposeFlux(Vector3r pos, Real flux);

101

unsigned int imposePressure(Vector3r pos, Real p);

102

unsigned int imposePressureFromId(long id, Real p);

102

103

void setImposedPressure(unsigned int cond, Real p);

103

104

void clearImposedPressure();

104

105

void clearImposedFlux();

165

166

Real volumeCellTripleFictious (Cellhandle cell);

166

167

template<class Cellhandle>

167

168

Real volumeCell (Cellhandle cell);

169

template<class Cellhandle>

170

Vector3r cellBarycenterFromHandle (Cellhandle cell) {return makeVector3r(solver->cellBarycenter(cell));}

171

Vector3r cellBarycenterFromId (long id) {

172

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);}

173

else return cellBarycenterFromHandle(solver->T[solver->currentTes].cellHandles[id]);}

168

174

void Oedometer_Boundary_Conditions();

169

175

void averageRealCellVelocity();

170

176

void saveVtk(const char* folder) {solver->saveVtk(folder);}

175

181

double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}

176

182

int getCell(double posX, double posY, double posZ){return solver->getCell(posX, posY, posZ);}

177

183

unsigned int nCells(){return solver->T[solver->currentTes].cellHandles.size();}

184

CELL_VECTOR_GETTER(,cellCenter)

185

CELL_VECTOR_GETTER(.averageCellVelocity,cellVelocity)

186

CELL_SCALAR_GETTER(double,.p(),cellPressure)

187

CELL_SCALAR_SETTER(double,.p(),setCellPressure)

188

CELL_SCALAR_GETTER(bool,.Pcondition,cellPImposed)

189

CELL_SCALAR_SETTER(bool,.Pcondition,setCellPImposed)

178

190

boost::python::list getVertices(unsigned int id){

179

191

boost::python::list ids;

180

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;}

192

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;}

181

193

for (unsigned int i=0;i<4;i++) ids.append(solver->T[solver->currentTes].cellHandles[id]->vertex(i)->info().id());

182

194

return ids;

183

195

}

196

void blockCell(unsigned int id, bool blockPressure){

197

if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return;}

198

solver->T[solver->currentTes].cellHandles[id]->info().blocked=true;

199

solver->T[solver->currentTes].cellHandles[id]->info().Pcondition=blockPressure;

200

}

201

184

202

double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}

185

203

double averagePressure(){return solver->averagePressure();}

186

204

213

231

if (solver->T[solver->currentTes].Volume(id) == -1) {compTessVolumes();/* LOG_WARN("Computing all volumes now, as you did not request it explicitely.");*/}

214

232

return (solver->T[solver->currentTes].Max_id() >= id) ? solver->T[solver->currentTes].Volume(id) : -1;}

215

233

216

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.",

234

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`.",

217

235

((bool,isActivated,true,,"Activates Flow Engine"))

218

236

((bool,first,true,,"Controls the initialization/update phases"))

219

237

((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."))

225

243

((double, sineMagnitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied (p )"))

226

244

((double, sineAverage, 0,,"Pressure value (average) when sinusoidal pressure is applied"))

227

245

((bool, debug, false,,"Activate debug messages"))

228

((double, wallThickness,0.001,,"Walls thickness"))

246

((double, wallThickness,0,,"Walls thickness"))

229

247

((double,pZero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model."))

230

248

((double,tolerance,1e-06,,"Gauss-Seidel tolerance"))

231

249

((double,relax,1.9,,"Gauss-Seidel relaxation"))

274

292

#endif

275

293

((vector<Real>, boundaryPressure,vector<Real>(),,"values defining pressure along x-axis for the top surface. See also :yref:`@TEMPLATE_FLOW_NAME@::boundaryXPos`"))

276

294

((vector<Real>, boundaryXPos,vector<Real>(),,"values of the x-coordinate for which pressure is defined. See also :yref:`@TEMPLATE_FLOW_NAME@::boundaryPressure`"))

295

((string,blockHook,"",,"Python command to be run after triangulation to define blocked cells (see also :yref:`TemplateFlowEngine_@TEMPLATE_FLOW_NAME@.blockCell`)"))

277

296

278

297

/*deprec*/

279

298

((meanK_opt,clampKValues,"the name changed"))

291

310

metisForced=false;

292

311

293

312

.def("imposeFlux",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::imposeFlux,(boost::python::arg("pos"),boost::python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")

313

.def("imposePressureFromId",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::imposePressureFromId,(boost::python::arg("id"),boost::python::arg("p")),"Impose pressure in cell of index 'id' (after remeshing the same condition will apply for the same location, regardless of what the new cell index is at this location). The index of the condition itself is returned (for multiple imposed pressures at different points).")

294

314

.def("imposePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::imposePressure,(boost::python::arg("pos"),boost::python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")

295

315

.def("setImposedPressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::setImposedPressure,(boost::python::arg("cond"),boost::python::arg("p")),"Set pressure value at the point indexed 'cond'.")

296

316

.def("clearImposedPressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::clearImposedPressure,"Clear the list of points with pressure imposed.")

301

321

.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}}")

302

322

.def("edgeSize",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::edgeSize,"Return the number of interactions.")

303

323

.def("OSI",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::OSI,"Return the number of interactions only between spheres.")

304

305

324

.def("saveVtk",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")

306

325

.def("avFlVelOnSph",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::avFlVelOnSph,(boost::python::arg("idSph")),"compute a sphere-centered average fluid velocity")

307

326

.def("fluidForce",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::fluidForce,(boost::python::arg("idSph")),"Return the fluid force on sphere idSph.")

322

341

.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)")

323

342

.def("emulateAction",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")

324

343

.def("getCell",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getCell,(boost::python::arg("pos")),"get id of the cell containing (X,Y,Z).")

344

.def("getCellBarycenter",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::cellBarycenterFromId,(boost::python::arg("id")),"get barycenter of cell 'id'.")

345

.def("getCellCenter",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::cellCenter,(boost::python::arg("id")),"get voronoi center of cell 'id'.")

346

.def("getCellPressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::cellPressure,(boost::python::arg("id")),"get pressure in cell 'id'.")

347

.def("setCellPressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::setCellPressure,(boost::python::arg("id"),boost::python::arg("pressure")),"set pressure in cell 'id'.")

348

.def("getCellPImposed",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::cellPImposed,(boost::python::arg("id")),"get the status of cell 'id' wrt imposed pressure.")

349

.def("setCellPImposed",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::setCellPImposed,(boost::python::arg("id"),boost::python::arg("pImposed")),"make cell 'id' assignable with imposed pressure.")

325

350

.def("nCells",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::nCells,"get the total number of finite cells in the triangulation.")

351

.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.")

326

352

.def("getVertices",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getVertices,(boost::python::arg("id")),"get the vertices of a cell")

327

353

#ifdef LINSOLV

328

354

.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

362

.def("averageVelocity",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averageVelocity,"measure the mean velocity in the period")

337

363

)

338

364

};

339

// Definition of functions in a separate file for clarity

365

340

366

341

367

class FlowCellInfo_@TEMPLATE_FLOW_NAME@ : public CGT::SimpleCellInfo {

342

368

public:

344

370

unsigned int index;

345

371

int volumeSign;

346

372

bool Pcondition;

373

bool blocked;//exclude cell from the fluid domain

347

374

Real invVoidV;

348

375

Real t;

349

376

int fict;

377

404

rayHydr.resize(4, 0);

378

405

invSumK=index=volumeSign=s=volumeVariation=pression=invVoidV=fict=0;

379

406

isFictious=false; Pcondition = false; isGhost = false;

380

isvisited = false;

407

isvisited = false;

381

408

isGhost=false;

409

blocked=false;

382

410

}

383

411

bool isGhost;

384

412

double invSumK;

414

442

inline const CVector ghostShift (void) const {return CGAL::NULL_VECTOR;}

415

443

};

416

444

445

// Implementation of functions in this separate file for clarity

417

446

#include "FlowEngine_@TEMPLATE_FLOW_NAME@.ipp"

Older »