yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3514
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
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;