------------------------------------------------------------
revno: 3514
committer: Bruno Chareyre <bruno.chare...@hmg.inpg.fr>
timestamp: Wed 2014-10-29 17:49:20 +0100
message:
FlowEngine: enhanced setter macro, detect changes in imposed pressure and
update the linear system accordingly
modified:
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphereLinSolv.hpp
lib/triangulation/FlowBoundingSphereLinSolv.ipp
lib/triangulation/PeriodicFlow.hpp
lib/triangulation/PeriodicFlowLinSolv.hpp
pkg/pfv/FlowEngine.hpp
pkg/pfv/FlowEngine.hpp.in
pkg/pfv/FlowEngine.ipp.in
pkg/pfv/PeriodicFlowEngine.cpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to
https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2014-10-15 10:12:53 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2014-10-29 16:49:20 +0000
@@ -95,8 +95,8 @@
void computePermeability();
virtual void gaussSeidel (Real dt=0);
virtual void resetNetwork();
- virtual void resetLinearSystem();
-
+ virtual void resetLinearSystem();//reset both A and B in the linear system A*P=B, done typically after updating the mesh
+ virtual void resetRHS() {};////reset only B in the linear system A*P=B, done typically after changing values of imposed pressures
double kFactor; //permeability moltiplicator
std::string key; //to give to consolidation files a name with iteration number
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.hpp'
--- lib/triangulation/FlowBoundingSphereLinSolv.hpp 2014-10-10 11:45:21 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.hpp 2014-10-29 16:49:20 +0000
@@ -53,13 +53,15 @@
using FlowType::computedOnce;
using FlowType::resetNetwork;
using FlowType::tesselation;
+ using FlowType::resetRHS;
//! TAUCS DECs
vector<FiniteCellsIterator> orderedCells;
bool isLinearSystemSet;
bool isFullLinearSystemGSSet;
bool areCellsOrdered;//true when orderedCells is filled, turn it false after retriangulation
-
+ bool updatedRHS;
+
#ifdef EIGENSPARSE_LIB
//Eigen's sparse matrix and solver
Eigen::SparseMatrix<double> A;
@@ -180,6 +182,7 @@
computedOnce=true;
}
virtual void resetLinearSystem();
+ virtual void resetRHS() {updatedRHS=false;};
};
} //namespace CGT
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.ipp'
--- lib/triangulation/FlowBoundingSphereLinSolv.ipp 2014-10-28 15:29:05 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.ipp 2014-10-29 16:49:20 +0000
@@ -71,6 +71,7 @@
isLinearSystemSet=0;
isFullLinearSystemGSSet=0;
areCellsOrdered=0;//true when orderedCells is filled, turn it false after retriangulation
+ updatedRHS=false;
ZERO=0;
#ifdef TAUCS_LIB
T_A = &SystemMatrix;
@@ -156,7 +157,7 @@
orderedCells.clear();
const FiniteCellsIterator cellEnd = Tri.finite_cells_end();
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
- orderedCells.push_back(cell);
+ orderedCells.push_back(cell); cell->info().index=0;
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>());
@@ -227,6 +228,7 @@
}
}
}
+ updatedRHS = true;
if (!isLinearSystemSet) {
if (useSolver==1 || useSolver==2){
#ifdef TAUCS_LIB
@@ -514,7 +516,7 @@
int FlowBoundingSphereLinSolv<_Tesselation,FlowType>::eigenSolve(Real dt)
{
#ifdef EIGENSPARSE_LIB
- if (!isLinearSystemSet || (isLinearSystemSet && reApplyBoundaryConditions())) ncols = setLinearSystem(dt);
+ if (!isLinearSystemSet || (isLinearSystemSet && reApplyBoundaryConditions()) || !updatedRHS) ncols = setLinearSystem(dt);
copyCellsToLin(dt);
//FIXME: we introduce new Eigen vectors, then we have to copy from/to c-arrays, can be optimized later
Eigen::VectorXd eb(ncols); Eigen::VectorXd ex(ncols);
=== modified file 'lib/triangulation/PeriodicFlow.hpp'
--- lib/triangulation/PeriodicFlow.hpp 2014-10-15 06:44:01 +0000
+++ lib/triangulation/PeriodicFlow.hpp 2014-10-29 16:49:20 +0000
@@ -28,7 +28,7 @@
//painfull, but we need that for templates inheritance...
using _N::T; using _N::xMin; using _N::xMax; using _N::yMin; using _N::yMax; using _N::zMin; using _N::zMax; using _N::Rmoy; using _N::sectionArea; using _N::Height; using _N::vTotal; using _N::currentTes; using _N::debugOut; using _N::nOfSpheres; using _N::xMinId; using _N::xMaxId; using _N::yMinId; using _N::yMaxId; using _N::zMinId; using _N::zMaxId; using _N::boundsIds; using _N::cornerMin; using _N::cornerMax; using _N::VSolidTot; using _N::Vtotalissimo; using _N::vPoral; using _N::sSolidTot; using _N::vPoralPorosity; using _N::vTotalPorosity; using _N::boundaries; using _N::idOffset; using _N::vtkInfiniteVertices; using _N::vtkInfiniteCells; using _N::num_particles; using _N::boundingCells; using _N::facetVertices; using _N::facetNFictious;
- 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;
+ 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; using BaseFlowSolver::resetRHS;
//same for functions
using _N::defineFictiousCells; using _N::addBoundingPlanes; using _N::boundary;
=== modified file 'lib/triangulation/PeriodicFlowLinSolv.hpp'
--- lib/triangulation/PeriodicFlowLinSolv.hpp 2014-07-02 16:18:24 +0000
+++ lib/triangulation/PeriodicFlowLinSolv.hpp 2014-10-29 16:49:20 +0000
@@ -28,9 +28,9 @@
//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;
+ 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; using BaseFlowSolver::resetRHS;
/// 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;
+ 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
=== modified file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp 2014-10-10 11:45:21 +0000
+++ pkg/pfv/FlowEngine.hpp 2014-10-29 16:49:20 +0000
@@ -17,10 +17,13 @@
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_SCALAR_SETTER(type, param, setterName) \
+#define CELL_SCALAR_SETTER_EXTRA(type, param, setterName, extra) \
void setterName(unsigned int id, type value){\
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()param=value;}
+ solver->T[solver->currentTes].cellHandles[id]->info()param=value;\
+ extra;}
+
+#define CELL_SCALAR_SETTER(type, param, setterName) CELL_SCALAR_SETTER_EXTRA(type, param, setterName,return)
#define CELL_VECTOR_GETTER(param, getterName) \
Vector3r getterName(unsigned int id){\
=== modified file 'pkg/pfv/FlowEngine.hpp.in'
--- pkg/pfv/FlowEngine.hpp.in 2014-10-28 15:29:05 +0000
+++ pkg/pfv/FlowEngine.hpp.in 2014-10-29 16:49:20 +0000
@@ -183,9 +183,9 @@
CELL_VECTOR_GETTER(,cellCenter)
CELL_VECTOR_GETTER(.averageCellVelocity,cellVelocity)
CELL_SCALAR_GETTER(double,.p(),cellPressure)
- CELL_SCALAR_SETTER(double,.p(),setCellPressure)
+ CELL_SCALAR_SETTER_EXTRA(double,.p(),setCellPressure,solver->resetRHS())
CELL_SCALAR_GETTER(bool,.Pcondition,cellPImposed)
- CELL_SCALAR_SETTER(bool,.Pcondition,setCellPImposed)
+ CELL_SCALAR_SETTER_EXTRA(bool,.Pcondition,setCellPImposed,solver->resetLinearSystem())
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;}
=== modified file 'pkg/pfv/FlowEngine.ipp.in'
--- pkg/pfv/FlowEngine.ipp.in 2014-10-28 15:29:05 +0000
+++ pkg/pfv/FlowEngine.ipp.in 2014-10-29 16:49:20 +0000
@@ -622,6 +622,7 @@
O1O2 = O1O2Vector.norm();
normal= (O1O2Vector/O1O2);
surfaceDist = O1O2 - r2 - r1;
+ //FIXME: what is that?
O1CVector = (O1O2/2. + (pow(r1,2) - pow(r2,2)) / (2.*O1O2))*normal;
O2CVector = -(O1O2Vector - O1CVector);
meanRad = (r2 + r1)/2.;
=== modified file 'pkg/pfv/PeriodicFlowEngine.cpp'
--- pkg/pfv/PeriodicFlowEngine.cpp 2014-10-21 21:28:01 +0000
+++ pkg/pfv/PeriodicFlowEngine.cpp 2014-10-29 16:49:20 +0000
@@ -493,7 +493,7 @@
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) {LOG_ERROR("problem here, 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;
cell->info().id=cell->info().baseIndex;
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : yade-dev@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
