yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2560
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Thu 2010-11-18 21:01:50 +0100
message:
  - remove a warning in Tesselation.cpp
  - fix FlowBonding sphere for first iteration and a few other things
  - make GaussSeidel virtual in FBS
modified:
  lib/triangulation/FlowBoundingSphere.cpp
  lib/triangulation/FlowBoundingSphere.hpp
  lib/triangulation/Tesselation.cpp
  pkg/dem/FlowEngine.cpp


--
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-17 20:07:21 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2010-11-18 20:01:50 +0000
@@ -32,7 +32,7 @@
 #define TESS_BASED_FORCES
 #define FACET_BASED_FORCES 1
 
-#ifdef YADE_OPEN_MP
+#ifdef YADE_OPENMP
 //   #define GS_OPEN_MP //It should never be defined if Yade is not using openmp
 #endif
 
@@ -92,6 +92,9 @@
 	OUTPUT_BOUDARIES_RADII = false;
 	RAVERAGE = false; /** if true use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
 }
+
+void FlowBoundingSphere::ResetNetwork() {noCache=true;}
+
 Tesselation& FlowBoundingSphere::Compute_Action()
 {
         return Compute_Action(0,NULL,NULL);
@@ -1036,14 +1039,14 @@
 		#endif
 		j++;
 //                  		if (j % 100 == 0) {
-//                         // cout << "pmax " << p_max << "; pmoy : " << p_moy << "; dpmax : " << dp_max << endl;
+//                         cout << "pmax " << p_max << "; pmoy : " << p_moy << "; dpmax : " << dp_max << endl;
 //                         cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;
 //                         //     save_vtk_file ( Tri );
 //                  }
 	#ifdef GS_OPEN_MP
 	} while (j<1500);
 	#else
-	} while ((dp_max/p_max) > tolerance /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
+	} while ((dp_max/p_max) > tolerance && j<4000 /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
 	#endif
 	}
 
@@ -1171,7 +1174,6 @@
                         Fictious+=1;
                 }
         }
-
         int fict=0, real=0;
         for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
                 if (v->info().isFictious) fict+=1;

=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2010-11-17 20:07:21 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2010-11-18 20:01:50 +0000
@@ -64,7 +64,8 @@
 		void Localize();
 
 		void Compute_Permeability();
-		void GaussSeidel ( );
+		virtual void GaussSeidel ( );
+		virtual void ResetNetwork();
 
 // 		void Compute_Forces ();
 		void Fictious_cells ( );
@@ -167,17 +168,15 @@
 // 		double surface_external_double_fictious(Cell_handle cell, Boundary b);
 // 		double surface_external_single_fictious(Cell_handle cell, Boundary b);
 
-		///TAUCS
-		int SetLinearSystem();
-		int SetLinearSystemFullGS();
-		void VectorizedGaussSeidel ();
-		int TaucsSolveTest();
-		int PardisoSolveTest();
-		///END_TAUCS
-
 };
 
 } //namespace CGT
+
+#ifdef LINSOLV
+#include "FlowBoundingSphereLinSolv.hpp"
+#endif
+
+
 #endif //FLOW_ENGINE
 
 #endif

=== modified file 'lib/triangulation/Tesselation.cpp'
--- lib/triangulation/Tesselation.cpp	2010-11-02 16:23:44 +0000
+++ lib/triangulation/Tesselation.cpp	2010-11-18 20:01:50 +0000
@@ -179,7 +179,7 @@
 				pass = true;
 				if ( !Tri->is_infinite ( ( *facet ).first ) && !Tri->is_infinite ( ( *facet ).first->neighbor ( ( *facet ).second ) ) )
 				{
-					cout << "p.x()     = " << p.x() << "p.y()     = " << p.y() << "p.z()     = " << p.z() << endl;
+// 					cout << "p.x()     = " << p.x() << "p.y()     = " << p.y() << "p.z()     = " << p.z() << endl;
 					p = ( *facet ).first->info();
 					( *Coordonnes ) [k++] = p.x(); ( *Coordonnes ) [k++] = p.y(); ( *Coordonnes ) [k++] = p.z();
 					p = ( *facet ).first->neighbor ( ( *facet ).second )->info();

=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp	2010-11-17 20:07:21 +0000
+++ pkg/dem/FlowEngine.cpp	2010-11-18 20:01:50 +0000
@@ -24,81 +24,78 @@
 {
 }
 
-void FlowEngine::action ( )
+void FlowEngine::action()
 {
 	if (!flow) {
-	flow = shared_ptr<FlowSolver> (new FlowSolver);
-	first=true;Update_Triangulation=false;/*IS=0.f;*/eps_vol_max=0.f;Eps_Vol_Cumulative=0.f;ReTrg=1.0;}
-	if ( !isActivated ) return;
+		flow = shared_ptr<FlowSolver> (new FlowSolver);
+		first=true;
+		Update_Triangulation=false;/*IS=0.f;*/
+		eps_vol_max=0.f;
+		Eps_Vol_Cumulative=0.f;
+		ReTrg=1.0;
+	}
+	if (!isActivated) return;
 	else
 	{
 		timingDeltas->start();
 
-		if ( !triaxialCompressionEngine )
-		{vector<shared_ptr<Engine> >::iterator itFirst = scene->engines.begin();vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();for ( ;itFirst!=itLast; ++itFirst ){
-		if ( ( *itFirst )->getClassName() == "TriaxialCompressionEngine" ){
-		cout << "stress controller engine found - FlowEngine" << endl;
-		triaxialCompressionEngine =  YADE_PTR_CAST<TriaxialCompressionEngine> ( *itFirst );}}
-		if ( !triaxialCompressionEngine ) cout << "stress controller engine NOT found" << endl;}
+		if (!triaxialCompressionEngine)
+		{
+			vector<shared_ptr<Engine> >::iterator itFirst = scene->engines.begin();
+			vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();
+			for (;itFirst!=itLast; ++itFirst) {
+				if ((*itFirst)->getClassName() == "TriaxialCompressionEngine") {
+					cout << "stress controller engine found - FlowEngine" << endl;
+					triaxialCompressionEngine =  YADE_PTR_CAST<TriaxialCompressionEngine> (*itFirst);}}
+			if (!triaxialCompressionEngine) cout << "stress controller engine NOT found" << endl;
+		}
 
 		currentStress = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];
 		currentStrain = triaxialCompressionEngine->strain[1];
 		current_state = triaxialCompressionEngine->currentState;
 
-		if ( current_state==3 )
+// 		if ( current_state==3 )
 		{
-			if ( first ) Build_Triangulation( P_zero );
-
+			if (first) Build_Triangulation(P_zero);
 			timingDeltas->checkpoint("Triangulating");
 
-			eps_vol_max=0.f;	
+			eps_vol_max=0.f;
 			UpdateVolumes ( );
 			Eps_Vol_Cumulative += eps_vol_max;
-			if (Eps_Vol_Cumulative > ReTrg*EpsVolPercent_RTRG) {Update_Triangulation = true; ReTrg++;}
-
+			if (Eps_Vol_Cumulative > ReTrg*EpsVolPercent_RTRG) {
+				Update_Triangulation = true;
+				ReTrg++;}
 			timingDeltas->checkpoint("Update_Volumes");
 
-			///Compute flow and and forces here
-			if (!first) flow->GaussSeidel ( );
-
-			timingDeltas->checkpoint("Gauss-Seidel");
-
-			if (save_mplot) {
-				int j = scene->iter;
-				char plotfile [50];
-				mkdir("./mplot", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
-				string visu_consol = "./mplot/"+flow->key+"%d_Visu_Consol";
-				const char* key_visu_consol = visu_consol.c_str();
-				sprintf(plotfile, key_visu_consol, j);
-				char *gg = plotfile;
-				flow->mplot(gg);}
-
-			if (save_mgpost) flow->MGPost();
-
-			if (!CachedForces) flow->ComputeFacetForces();
-			else flow->ComputeFacetForcesWithCache();
-
-			timingDeltas->checkpoint("Compute_Forces");
-
-			///End Compute flow and forces
-			//int number_of_particles = flow->num_particles;
-			CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end ();
-			Vector3r f; int id;
-			//IS = 0.f;
-			//std::ofstream Istab ("Stability.txt", std::ios::app);
-			for ( CGT::Finite_vertices_iterator V_it = flow->T[flow->currentTes].Triangulation().finite_vertices_begin (); V_it !=  vertices_end; V_it++ )
-			{
-				id = V_it->info().id();
-				for ( int y=0;y<3;y++ ) f[y] = ( V_it->info().forces ) [y];
-				scene->forces.addForce ( id, f );
-				//scene->forces.sync();
-				//IS += (scene->forces.getForce(id).norm())/number_of_particles;
+			if (!first) {
+				///Compute flow and and forces here
+				flow->GaussSeidel();
+				timingDeltas->checkpoint("Gauss-Seidel");
+				if (save_mgpost) flow->MGPost();
+				if (!CachedForces) flow->ComputeFacetForces();
+				else flow->ComputeFacetForcesWithCache();
+				timingDeltas->checkpoint("Compute_Forces");
+
+				///End Compute flow and forces
+				//int number_of_particles = flow->num_particles;
+				CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
+				Vector3r f;
+				int id;
+				//IS = 0.f;
+				//std::ofstream Istab ("Stability.txt", std::ios::app);
+				for (CGT::Finite_vertices_iterator V_it = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); V_it !=  vertices_end; V_it++)
+				{
+					id = V_it->info().id();
+					for (int y=0;y<3;y++) f[y] = (V_it->info().forces)[y];
+					scene->forces.addForce(id, f);
+					//scene->forces.sync();
+					//IS += (scene->forces.getForce(id).norm())/number_of_particles;
+				}
+				timingDeltas->checkpoint("Applying Forces");
 			}
 			//cout << "STABILITY INDEX - IS = " << IS << endl;
 			//Istab << scene->time << " " << IS << endl;
 
-			timingDeltas->checkpoint("Applying Forces");
-
 			Real time = scene->time;
 			int j = scene->iter;
 
@@ -116,30 +113,28 @@
 				max_p << j << " " << time << " " << MaxPressure << endl;
 
 				std::ofstream settle("settle.txt", std::ios::app);
-				settle << j << " " << time << " " << currentStrain << endl;}
-
+				settle << j << " " << time << " " << currentStrain << endl;
+			}
 // 			if ( scene->iter % PermuteInterval == 0 )
 // 			{ Update_Triangulation = true; }
 
-			if ( Update_Triangulation && !first) { Build_Triangulation( );}
-
+			if (Update_Triangulation && !first) {
+				Build_Triangulation();
+				Update_Triangulation = false;
+			}
 			first=false;
-			Update_Triangulation = false;
-
-			timingDeltas->checkpoint("Storing Max Pressure");
-
-// 				flow->Average_Grain_Velocity();
+//			flow->Average_Grain_Velocity();
 			if (save_vtk) flow->save_vtk_file();
 
 			if (slice_pressures)
 			{
-			  char slifile [30];
-			  mkdir("./Slices", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
-			  string slice = "./Slices/Slice_"+flow->key+"_%d";
-			  const char* keyslice = slice.c_str();
-			  sprintf (slifile, keyslice, j);
-			  char *f = "slifile";
-			  flow->SliceField(f);
+				char slifile [30];
+				mkdir("./Slices", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+				string slice = "./Slices/Slice_"+flow->key+"_%d";
+				const char* keyslice = slice.c_str();
+				sprintf(slifile, keyslice, j);
+				char *f = "slifile";
+				flow->SliceField(f);
 			}
 // 				int numero = flow->Average_Cell_Velocity(id_sphere, flow->T[flow->currentTes].Triangulation());
 // 				flow->vtk_average_cell_velocity(flow->T[flow->currentTes].Triangulation(), id_sphere, numero);
@@ -172,20 +167,18 @@
 
 void FlowEngine::Build_Triangulation (double P_zero)
 {
-	if (first)
-	{
-		flow->currentTes=0;
-		flow->Vtotalissimo=0; flow->Vsolid_tot=0; flow->Vporale=0; flow->Ssolid_tot=0;
-		flow->SLIP_ON_LATERALS=slip_boundary;
-		flow->key = triaxialCompressionEngine->Key;
-		flow->k_factor = permeability_factor;
-		flow->DEBUG_OUT = Debug;
-	}
-	else
-	{
+	flow->ResetNetwork();
+	if (first) flow->currentTes=0;
+	else {
 		flow->currentTes=!flow->currentTes;
-		if (Debug) cout << "--------RETRIANGULATION-----------" << endl;
-	}
+		if (Debug) cout << "--------RETRIANGULATION-----------" << endl;}
+
+	flow->Vtotalissimo=0; flow->Vsolid_tot=0; flow->Vporale=0; flow->Ssolid_tot=0;
+	flow->SLIP_ON_LATERALS=slip_boundary;
+	flow->key = triaxialCompressionEngine->Key;
+	flow->k_factor = permeability_factor;
+	flow->DEBUG_OUT = Debug;
+
 	flow->T[flow->currentTes].Clear();
 	flow->T[flow->currentTes].max_id=-1;
 	flow->x_min = 1000.0, flow->x_max = -10000.0, flow->y_min = 1000.0, flow->y_max = -10000.0, flow->z_min = 1000.0, flow->z_max = -10000.0;
@@ -206,7 +199,7 @@
 	{
 		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++ ){cell->info().dv() = 0; cell->info().p() = 0;}
-		/*if (compute_K) {*/flow->TOLERANCE=1e-06; K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );/*}*/
+		if (compute_K) {flow->TOLERANCE=1e-06; K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
 		BoundaryConditions();
 		flow->Initialize_pressures( P_zero );
   		if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Pressure, 15);
@@ -221,14 +214,15 @@
 		for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){cell->info().dv() = 0; cell->info().p() = 0;}
 		if (compute_K) {flow->TOLERANCE=1e-06; K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
 		BoundaryConditions();
-		flow->Initialize_pressures( P_zero );// FIXME : why, if we are going to interpolate after that?
+		flow->Initialize_pressures(P_zero);// FIXME : why, if we are going to interpolate after that?
 		flow->TOLERANCE=Tolerance;//So it can be changed at run time
-		flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
-
+		flow->Interpolate (flow->T[!flow->currentTes], flow->T[flow->currentTes]);
  		if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Pressure, 15);
+		else if (TimeBC) flow->ApplySinusoidalPressure_Space_Time(flow->T[flow->currentTes].Triangulation(), Sinus_Pressure, 15, scene->time, scene->dt);
+		flow->TOLERANCE=Tolerance;
+		flow->RELAX=Relax;
 	}
-	Initialize_volumes ( );
-	flow->noCache=true;
+	Initialize_volumes();
 }
 
 void FlowEngine::AddBoundary ()