yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2565
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Fri 2010-11-19 18:18:56 +0100
message:
  - handle different solvers.
  - small optimization on dt divisions.
modified:
  lib/triangulation/FlowBoundingSphere.cpp
  lib/triangulation/FlowBoundingSphere.hpp
  pkg/dem/FlowEngine.cpp
  pkg/dem/FlowEngine.hpp


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription

=== modified file 'lib/triangulation/FlowBoundingSphere.cpp'
--- lib/triangulation/FlowBoundingSphere.cpp	2010-11-18 20:01:50 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2010-11-19 17:18:56 +0000
@@ -27,7 +27,6 @@
 // #include "Vue3D.h" //FIXME implicit dependencies will look for this class (out of tree) even ifndef XVIEW
 #endif
 
-
 #define FAST
 #define TESS_BASED_FORCES
 #define FACET_BASED_FORCES 1
@@ -1033,7 +1032,7 @@
                 #endif
 		p_moy = sum_p/cell2;
                 dp_moy = sum_dp/cell2;
-		
+
 		#ifdef GS_OPEN_MP
 		#pragma omp master
 		#endif

=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2010-11-18 20:01:50 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2010-11-19 17:18:56 +0000
@@ -168,6 +168,9 @@
 // 		double surface_external_double_fictious(Cell_handle cell, Boundary b);
 // 		double surface_external_single_fictious(Cell_handle cell, Boundary b);
 
+		//Solver?
+		int useSolver;//(0 : GaussSeidel, 1 : TAUCS, 2 : PARDISO)
+
 };
 
 } //namespace CGT

=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp	2010-11-19 12:30:08 +0000
+++ pkg/dem/FlowEngine.cpp	2010-11-19 17:18:56 +0000
@@ -33,7 +33,8 @@
 		Update_Triangulation=false;/*IS=0.f;*/
 		eps_vol_max=0.f;
 		Eps_Vol_Cumulative=0.f;
-		ReTrg=1.0;
+		ReTrg=1;
+		retriangulationLastIter=0;
 	}
 	if (!isActivated) return;
 	else
@@ -63,9 +64,12 @@
 			eps_vol_max=0.f;
 			UpdateVolumes ( );
 			Eps_Vol_Cumulative += eps_vol_max;
-			if (Eps_Vol_Cumulative > ReTrg*EpsVolPercent_RTRG) {
+			if (Eps_Vol_Cumulative > EpsVolPercent_RTRG || retriangulationLastIter>100) {
 				Update_Triangulation = true;
-				ReTrg++;}
+				Eps_Vol_Cumulative=0;
+				retriangulationLastIter=0;
+				ReTrg++;
+			} else  retriangulationLastIter++;
 			timingDeltas->checkpoint("Update_Volumes");
 
 			if (!first) {
@@ -179,6 +183,7 @@
 	flow->key = triaxialCompressionEngine->Key;
 	flow->k_factor = permeability_factor;
 	flow->DEBUG_OUT = Debug;
+	flow->useSolver = useSolver;
 
 	flow->T[flow->currentTes].Clear();
 	flow->T[flow->currentTes].max_id=-1;
@@ -348,6 +353,7 @@
 {
 	if (Debug) cout << "Updating volumes.............." << endl;
 	Real deltaT = scene->dt;
+	Real invDeltaT = 1/scene->dt;
 	CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
 	for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ )
 	{
@@ -355,25 +361,25 @@
 		{
 			case ( 3 ):
 			{
-				cell->info().dv() = ( Volume_cell_triple_fictious ( cell ) - cell->info().volume() ) /deltaT;
+				cell->info().dv() = ( Volume_cell_triple_fictious ( cell ) - cell->info().volume() ) *invDeltaT;
 				eps_vol_max = max(eps_vol_max, (abs(cell->info().dv()*deltaT))/cell->info().volume());
 				cell->info().volume() = Volume_cell_triple_fictious ( cell );
 			}break;
 			case ( 2 ) :
 			{
-				cell->info().dv() = ( Volume_cell_double_fictious ( cell )-cell->info().volume() ) /deltaT;
+				cell->info().dv() = ( Volume_cell_double_fictious ( cell )-cell->info().volume() ) *invDeltaT;
 				eps_vol_max = max(eps_vol_max, (abs(cell->info().dv()*deltaT))/cell->info().volume());
 				cell->info().volume() = Volume_cell_double_fictious ( cell );
 			}break;
 			case ( 1 ) :
 			{
-				cell->info().dv() = ( Volume_cell_single_fictious ( cell )-cell->info().volume() ) /deltaT;
+				cell->info().dv() = ( Volume_cell_single_fictious ( cell )-cell->info().volume() ) *invDeltaT;
 				eps_vol_max = max(eps_vol_max, (abs(cell->info().dv()*deltaT))/cell->info().volume());
 				cell->info().volume() = Volume_cell_single_fictious ( cell );
 			}break;
 			case ( 0 ) :
 			{
-				cell->info().dv() = ( Volume_cell ( cell )-cell->info().volume() ) /deltaT;
+				cell->info().dv() = ( Volume_cell ( cell )-cell->info().volume() ) *invDeltaT;
 				eps_vol_max = max(eps_vol_max, (abs(cell->info().dv()*deltaT))/cell->info().volume());
 				cell->info().volume() = Volume_cell ( cell );
 			}break;

=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp	2010-11-17 20:07:21 +0000
+++ pkg/dem/FlowEngine.hpp	2010-11-19 17:18:56 +0000
@@ -24,6 +24,7 @@
 	private:
 		shared_ptr<TriaxialCompressionEngine> triaxialCompressionEngine;
 		shared_ptr<FlowSolver> flow;
+		int retriangulationLastIter;
 	public :
 		Vector3r gravity;
 		int current_state;
@@ -82,6 +83,7 @@
 					((double, currentStress, 0,, "Current value of axial stress"))
 					((double, currentStrain, 0,, "Current value of axial strain"))
 					((int, intervals, 30,, "Number of layers for pressure measurements"))
+					((int, useSolver, 1,, "Solver to use"))
 					((bool, Flow_imposed_TOP_Boundary, true,, "if false involve pressure imposed condition"))
 					((bool, Flow_imposed_BOTTOM_Boundary, true,, "if false involve pressure imposed condition"))
 					((bool, Flow_imposed_FRONT_Boundary, true,, "if false involve pressure imposed condition"))