yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2096
committer: Emanuele Catalano <ecatalano@r2balme>
branch nick: trunk
timestamp: Tue 2010-03-23 12:30:28 +0100
message:
  - Fixed some error due to pointer scene* deletion
  - Added functions to write Mplotlib files
modified:
  lib/triangulation/FlowBoundingSphere.cpp
  lib/triangulation/FlowBoundingSphere.h
  pkg/dem/Engine/PartialEngine/FlowEngine.cpp
  pkg/dem/Engine/PartialEngine/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-03-18 15:27:50 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2010-03-23 11:30:28 +0000
@@ -8,6 +8,10 @@
 #include <utility>
 #include "vector"
 
+#include <sys/stat.h>
+#include <sys/types.h>
+
+
 //#define XVIEW
 #include "FlowBoundingSphere.h"//include after #define XVIEW
 
@@ -26,8 +30,6 @@
 
 const bool DEBUG_OUT = false;
 
-const double RELAX = 1.9;
-
 
 const double ONE_THIRD = 1.0/3.0;
 //! Use this factor, or minLength, to reduce max permeability values (see usage below))
@@ -71,6 +73,8 @@
         minPermLength=-1;
 	SLIP_ON_LATERALS = true;//no-slip/symmetry conditions on lateral boundaries
 	TOLERANCE = 1e-06;
+	RELAX = 1.9;
+	ks=0;
 }
 
 void FlowBoundingSphere::Compute_Action()
@@ -169,7 +173,7 @@
 
         /** END GAUSS SEIDEL */
         char* file ("Permeability");
-        double ks = Permeameter(Tri, boundary(y_min_id).value, boundary(y_max_id).value, SectionArea, H, file);
+        ks = Permeameter(Tri, boundary(y_min_id).value, boundary(y_max_id).value, SectionArea, H, file);
         clock.top("Permeameter");
 
         Compute_Forces();
@@ -1536,7 +1540,8 @@
 
 void FlowBoundingSphere::GaussSeidel()
 {
-        RTriangulation& Tri = T[currentTes].Triangulation();
+	std::ofstream iter ("Gauss_Iterations", std::ios::app);
+	RTriangulation& Tri = T[currentTes].Triangulation();
         int j = 0;
         double m, n, dp_max, p_max, sum_p, p_moy, dp_moy, dp, sum_dp;
         double tolerance = TOLERANCE;
@@ -1598,6 +1603,7 @@
         } while (((dp_max/p_max) > tolerance) /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
         cout << "pmax " << p_max << "; pmoy : " << p_moy << endl;
         cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;
+	iter << j << " " << dp_max/p_max << endl;
 
         //Display fluxes?
 //  bool ref =  Tri.finite_cells_begin()->info().isvisited;
@@ -1698,7 +1704,8 @@
 {
         static unsigned int number = 0;
         char filename[80];
-        sprintf(filename,"out_%d.vtk", number++);
+	mkdir("./VTK", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+        sprintf(filename,"./VTK/out_%d.vtk", number++);
 
         basicVTKwritter vtkfile((unsigned int) T.number_of_vertices()-vtk_infinite_vertices, (unsigned int) T.number_of_finite_cells()-vtk_infinite_cells);
 
@@ -1815,8 +1822,8 @@
 }
 
 double FlowBoundingSphere::PermeameterCurve(RTriangulation& Tri, char *filename, Real time)
-{
-        /** CONSOLIDATION CURVES **/
+{	
+	/** CONSOLIDATION CURVES **/
         Cell_handle permeameter;
         int n=0;
         int k=0;
@@ -1858,6 +1865,24 @@
 	return P_ave[intervals/2];
 }
 
+void FlowBoundingSphere::mplot (RTriangulation& Tri, char *filename)
+{
+	std::ofstream plot (filename, std::ios::out);
+	Cell_handle permeameter;
+	int intervals = 30;
+	double Rx = (x_max-x_min) /intervals;
+	double Ry = (y_max-y_min) /intervals;
+	double Z = (z_max+z_min)/2.0;
+	double Press = 0;
+	for (double Y=y_min; Y<y_max; Y=Y+Ry) {
+		for (double X=x_min; X<x_max; X=X+Rx) {
+				permeameter = Tri.locate(Point(X, Y, Z));
+				permeameter->info().p()<0? Press=0 : Press=permeameter->info().p();
+				plot << Y << " " << X << " " << Press << endl;
+		}
+	}
+}
+
 double FlowBoundingSphere::Sample_Permeability(RTriangulation& Tri, double x_Min,double x_Max ,double y_Min,double y_Max,double z_Min,double z_Max, string key)
 {
         double Section = (x_Max-x_Min) * (z_Max-z_Min);

=== modified file 'lib/triangulation/FlowBoundingSphere.h'
--- lib/triangulation/FlowBoundingSphere.h	2010-03-18 15:27:50 +0000
+++ lib/triangulation/FlowBoundingSphere.h	2010-03-23 11:30:28 +0000
@@ -38,12 +38,14 @@
 		bool currentTes;
 		bool SLIP_ON_LATERALS;
 		double TOLERANCE;
+		double RELAX;
+		double ks; //Hydraulic Conductivity
 		
 		Boundary boundaries [6];
 		short id_offset;
  		Boundary& boundary (int b) {return boundaries[b-id_offset];}
 		
-// 		void insert ( Real x, Real y, Real z, Real radius, int id );
+		void mplot (RTriangulation& Tri, char *filename);
 		void Localize ();
 		void Compute_Permeability();
 		void AddBoundingPlanes();

=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp	2010-03-20 12:40:44 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp	2010-03-23 11:30:28 +0000
@@ -14,6 +14,9 @@
 #include<yade/pkg-common/Wall.hpp>
 #include<yade/pkg-common/Box.hpp>
 
+#include <sys/stat.h>
+#include <sys/types.h>
+
 #ifdef FLOW_ENGINE
 
 YADE_REQUIRE_FEATURE (CGAL);
@@ -23,7 +26,7 @@
 {
 }
 
-void FlowEngine::action (Scene*)
+void FlowEngine::action ( )
 {
 	if (!flow) {flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);first=true;Update_Triangulation=false;}
 	if ( !isActivated ) return;
@@ -43,24 +46,31 @@
 
 		if ( current_state==3 )
 		{
-			if ( first || Update_Triangulation ) { Build_Triangulation( scene, P_zero );}
+			if ( first || Update_Triangulation ) { Build_Triangulation( P_zero );}
 // 			else
 // 			{
 				timingDeltas->checkpoint("Triangulating");
 				
-				UpdateVolumes ( scene );
+				UpdateVolumes ( );
 			
 				timingDeltas->checkpoint("Update_Volumes");
 			
 			///Compute flow and and forces here
 			
 				if (!first) flow->GaussSeidel ( );
-			
 				timingDeltas->checkpoint("Gauss-Seidel");
 				
+				if (save_mplot){int j = Omega::instance().getCurrentIteration();
+				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(flow->T[flow->currentTes].Triangulation(), gg);}
+				
 				if (save_vtk) flow->save_vtk_file(flow->T[flow->currentTes].Triangulation());
 			
-// 			flow->MGPost(flow->T[flow->currentTes].Triangulation());
+// 				flow->MGPost(flow->T[flow->currentTes].Triangulation());
 
 				flow->tess_based_force=tess_based_force;
 				flow->Compute_Forces ( );
@@ -86,7 +96,8 @@
 				int j = Omega::instance().getCurrentIteration();
 // 				int j = Omega::instance().getSimulationTime();
 				char file [50];
-				string consol = flow->key+"%d_Consol";
+				mkdir("./Consol", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+				string consol = "./Consol/"+flow->key+"%d_Consol";
 				const char* keyconsol = consol.c_str();
 				sprintf (file, keyconsol, j);
 				char *g = file;
@@ -145,8 +156,8 @@
 	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;
 
-	AddBoundary ( scene );
-	Triangulate ( scene );
+	AddBoundary ( );
+	Triangulate ( );
 	
 	cout << endl << "Tesselating------" << endl << endl;
 	flow->T[flow->currentTes].Compute();
@@ -163,6 +174,7 @@
 		Oedometer_Boundary_Conditions();
 		flow->Initialize_pressures( P_zero );
 		flow->TOLERANCE=Tolerance;
+		flow->RELAX=Relax;
 	}
 	else 
 	{
@@ -170,7 +182,7 @@
 		Update_Triangulation=!Update_Triangulation;
 	}
 
-	Initialize_volumes ( scene );
+	Initialize_volumes ( );
 	flow->DisplayStatistics ();
 }
 
@@ -247,10 +259,10 @@
 	{
 		switch ( cell->info().fictious() )
 		{
-			case ( 0 ) : cell->info().volume() = Volume_cell ( cell, scene ); break;
-			case ( 1 ) : cell->info().volume() = Volume_cell_single_fictious ( cell, scene ); break;
-			case ( 2 ) : cell->info().volume() = Volume_cell_double_fictious ( cell, scene ); break;
-			case ( 3 ) : cell->info().volume() = Volume_cell_triple_fictious ( cell, scene ); break;
+			case ( 0 ) : cell->info().volume() = Volume_cell ( cell ); break;
+			case ( 1 ) : cell->info().volume() = Volume_cell_single_fictious ( cell ); break;
+			case ( 2 ) : cell->info().volume() = Volume_cell_double_fictious ( cell ); break;
+			case ( 3 ) : cell->info().volume() = Volume_cell_triple_fictious ( cell ); break;
 		}
 	}
 	cout << "Volumes initialised." << endl;
@@ -268,23 +280,23 @@
 		{
 			case ( 3 ):
 			{
-				cell->info().dv() = ( Volume_cell_triple_fictious ( cell,scene ) - cell->info().volume() ) /deltaT;
-				cell->info().volume() = Volume_cell_triple_fictious ( cell,scene );
+				cell->info().dv() = ( Volume_cell_triple_fictious ( cell ) - cell->info().volume() ) /deltaT;
+				cell->info().volume() = Volume_cell_triple_fictious ( cell );
 			}break;
 			case ( 2 ) :
 			{
-				cell->info().dv() = ( Volume_cell_double_fictious ( cell,scene )-cell->info().volume() ) /deltaT;
-				cell->info().volume() = Volume_cell_double_fictious ( cell,scene );
+				cell->info().dv() = ( Volume_cell_double_fictious ( cell )-cell->info().volume() ) /deltaT;
+				cell->info().volume() = Volume_cell_double_fictious ( cell );
 			}break;
 			case ( 1 ) :
 			{
-				cell->info().dv() = ( Volume_cell_single_fictious ( cell,scene )-cell->info().volume() ) /deltaT;
-				cell->info().volume() = Volume_cell_single_fictious ( cell,scene );
+				cell->info().dv() = ( Volume_cell_single_fictious ( cell )-cell->info().volume() ) /deltaT;
+				cell->info().volume() = Volume_cell_single_fictious ( cell );
 			}break;
 			case ( 0 ) :
 			{
-				cell->info().dv() = ( Volume_cell ( cell,scene )-cell->info().volume() ) /deltaT;
-				cell->info().volume() = Volume_cell ( cell,scene );
+				cell->info().dv() = ( Volume_cell ( cell )-cell->info().volume() ) /deltaT;
+				cell->info().volume() = Volume_cell ( cell );
 			}break;
 		}
 	}

=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.hpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.hpp	2010-03-20 12:40:44 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.hpp	2010-03-23 11:30:28 +0000
@@ -19,7 +19,7 @@
 	private:
 		shared_ptr<TriaxialCompressionEngine> triaxialCompressionEngine;
 		shared_ptr<CGT::FlowBoundingSphere> flow;
-	public :		
+	public :	
 		Vector3r gravity;
 		int current_state;
 		Real wall_thickness;
@@ -43,11 +43,13 @@
 		YADE_CLASS_BASE_DOC_ATTRS_CTOR(FlowEngine,PartialEngine,"An engine to solve flow problem in saturated granular media",
 					((bool,isActivated,true,"Activates Flow Engine"))
 					((bool,first,true,"Controls the initialization/update phases"))
-					((bool,save_vtk,true,"Enable/disable vtk files creation to visualize pressure and force fields"))
+					((bool,save_vtk,false,"Enable/disable vtk files creation for visualization"))
+					((bool,save_mplot,false,"Enable/disable mplot files creation"))
 					((bool, slip_boundary, true, "Controls friction condition on lateral walls"))
 					((bool,currentTes,false,"Identifies the current triangulation/tesselation of pore space"))
 					((double,P_zero,0,"Initial internal pressure for oedometer test"))
 					((double,Tolerance,1e-06,"Gauss-Seidel Tolerance"))
+					((double,Relax,1.9,"Gauss-Seidel relaxation"))
 					((int,PermuteInterval,100000,"Pore space re-triangulation period"))
 					((bool,compute_K,true,"Activates permeability measure within a granular sample"))
 					((double,permeability_factor,1.0,"a permability multiplicator"))