yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2904
committer: Emanuele Catalano <catalano@xxxxxxxxxxx
branch nick: yade
timestamp: Tue 2011-08-23 12:24:56 +0200
message:
  - Consolidation files writers move to python scripts (Flow functions MeasurePorePressure and MeasurePressureProfile do what is needed)
  - Removed superfluous variables
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	2011-08-22 15:37:27 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2011-08-23 10:24:56 +0000
@@ -416,7 +416,14 @@
 	return result;
 }
 
-void FlowBoundingSphere::Measure_Pore_Pressure(double Wall_up_y, double Wall_down_y)
+double FlowBoundingSphere::MeasurePorePressure (double X, double Y, double Z)
+{
+  RTriangulation& Tri = T[currentTes].Triangulation();
+  Cell_handle cell = Tri.locate(Point(X,Y,Z));
+  return cell->info().p();
+}
+
+void FlowBoundingSphere::MeasurePressureProfile(double Wall_up_y, double Wall_down_y)
 {  
 	RTriangulation& Tri = T[currentTes].Triangulation();
         Cell_handle permeameter;
@@ -1416,37 +1423,6 @@
         }
 }
 
-double FlowBoundingSphere::PressureProfile(const char *filename, Real& time, int& intervals)
-{
-	RTriangulation& Tri = T[currentTes].Triangulation();
-	vector<double> Pressures;
-	
-	/** CONSOLIDATION CURVES **/
-        Cell_handle permeameter;
-        int n=0, k=0;
-        vector<double> P_ave;
-        std::ofstream consFile(filename, std::ios::out);
-
-        double Rx = (x_max-x_min) /intervals;
-        double Ry = (y_max-y_min) /intervals;
-        double Rz = (z_max-z_min) /intervals;
-
-	for (double Y=y_min; Y<=y_max+Ry/10; Y=Y+Ry) {
-                P_ave.push_back(0);
-		for (double X=x_min; X<=x_max+Ry/10; X=X+Rx) {
-			for (double Z=z_min; Z<=z_max+Ry/10; Z=Z+Rz) {
-                                P_ave[k]+=Tri.locate(Point(X, Y, Z))->info().p();
-				n++;
-                        }
-                }
-                P_ave[k]/= (n);
-                consFile<<k<<" "<<time<<" "<<P_ave[k]<<endl;
-                if (k==intervals/2) Pressures.push_back(P_ave[k]);
-                n=0; k++;
-	}
-	return P_ave[intervals/2];
-}
-
 void FlowBoundingSphere::mplot (char *filename)
 {
 	RTriangulation& Tri = T[currentTes].Triangulation();

=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2011-07-18 13:53:54 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2011-08-23 10:24:56 +0000
@@ -63,6 +63,7 @@
 		virtual void GaussSeidel ( );
 		virtual void ResetNetwork();
 
+		double MeasurePorePressure (double X, double Y, double Z);
 		void Fictious_cells ( );
 
 		double k_factor; //permeability moltiplicator
@@ -93,7 +94,7 @@
 		double Permeameter ( double P_Inf, double P_Sup, double Section, double DeltaY, const char *file );
 		double Sample_Permeability( double& x_Min,double& x_Max ,double& y_Min,double& y_Max,double& z_Min,double& z_Max, string key);
 		double Compute_HydraulicRadius (Cell_handle cell, int j );
-		double PressureProfile (const char *filename, Real& time, int& intervals );
+// 		double PressureProfile (const char *filename, Real& time, int& intervals );
 
 		double dotProduct ( Vecteur x, Vecteur y );
 		double Compute_EffectiveRadius(Cell_handle cell, int j);
@@ -120,7 +121,7 @@
 		void Average_Fluid_Velocity();
 		void ApplySinusoidalPressure(RTriangulation& Tri, double Amplitude, double Average_Pressure, double load_intervals);
 		bool isOnSolid  (double X, double Y, double Z);
-		void Measure_Pore_Pressure(double Wall_up_y, double Wall_down_y);
+		void MeasurePressureProfile(double Wall_up_y, double Wall_down_y);
 		vector<Real> Average_Fluid_Velocity_On_Sphere(unsigned int Id_sph);
 		//Solver?
 		int useSolver;//(0 : GaussSeidel, 1 : TAUCS, 2 : PARDISO)

=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp	2011-08-22 15:37:27 +0000
+++ pkg/dem/FlowEngine.cpp	2011-08-23 10:24:56 +0000
@@ -82,37 +82,6 @@
 		}
 		///End Compute flow and forces
 		timingDeltas->checkpoint("Applying Forces");
-
-		Real time = scene->time;
-		int j = scene->iter;
-
-		//FIXME: please Ema, put this in a separate function with python wrapper, so we keep action() code simple and usage easier
-		if (consolidation) {
-			mkdir("./Consol", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
-			stringstream sstr; sstr<<"./Consol/"<<flow->key<<j<<"_Consol";
-			string consol = sstr.str();
-			timingDeltas->checkpoint("Writing cons_files");
-			MaxPressure = flow->PressureProfile(consol.c_str(), time, intervals);
-			static bool consolidationFilesOpened=false;
-			if(!consolidationFilesOpened){
-				std::ofstream max_p("pressures.txt", std::ios::app);
-				std::ofstream settle("settle.txt", std::ios::app);
-				consolidationFilesOpened=True;}
-			max_p << j << " " << time << " " << MaxPressure << endl;
-			settle << j << " " << time << " " << currentStrain << endl;
-		}
-
-		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);
-			const char *f = "slifile";
-			flow->SliceField(f);
-		}
-// 		if (save_vtk) {flow->save_vtk_file();}
-		timingDeltas->checkpoint("Writing files");
 	}
 // 	if ( scene->iter % PermuteInterval == 0 )
 // 	{ Update_Triangulation = true; }
@@ -220,6 +189,7 @@
 	flow->useSolver = useSolver;
 	flow->VISCOSITY = viscosity;
 	flow->areaR2Permeability=areaR2Permeability;
+// 	flow->key = key;
 
 	flow->T[flow->currentTes].Clear();
 	flow->T[flow->currentTes].max_id=-1;

=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp	2011-08-22 15:37:27 +0000
+++ pkg/dem/FlowEngine.hpp	2011-08-23 10:24:56 +0000
@@ -61,7 +61,8 @@
 		Vector3r fluidForce(unsigned int id_sph) {const CGT::Vecteur& f=flow->T[flow->currentTes].vertex(id_sph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
 		Vector3r fluidShearForce(unsigned int id_sph) {return (flow->viscousShearForces.size()>id_sph)?flow->viscousShearForces[id_sph]:Vector3r::Zero();}
 		void setBoundaryVel(Vector3r vel) {topBoundaryVelocity=vel; Update_Triangulation=true;}
-		void PressureProfile(double wallUpY, double wallDownY) {return flow->Measure_Pore_Pressure(wallUpY,wallDownY);}
+		void PressureProfile(double wallUpY, double wallDownY) {return flow->MeasurePressureProfile(wallUpY,wallDownY);}
+		double MeasurePressure(double posX, double posY, double posZ){return flow->MeasurePorePressure(posX, posY, posZ);}
 
 		virtual ~FlowEngine();
 
@@ -72,13 +73,11 @@
 					((bool,first,true,,"Controls the initialization/update phases"))
 // 					((bool,save_vtk,false,,"Enable/disable vtk files creation for visualization"))
 					((bool,permeability_map,false,,"Enable/disable stocking of average permeability scalar in cell infos."))
-					((bool,save_mplot,false,,"Enable/disable mplot files creation"))
 					((bool, save_mgpost, false,,"Enable/disable mgpost file creation"))
 					((bool, slice_pressures, false, ,"Enable/Disable slice pressure measurement"))
 					((bool, velocity_profile, false, ,"Enable/Disable slice velocity measurement"))
 					((bool, consolidation,false,,"Enable/Disable storing consolidation files"))
 					((bool, slip_boundary, true,, "Controls friction condition on lateral walls"))
-					((bool, blocked_grains, false,, "Grains will/won't be moved by forces"))
 					((bool,WaveAction, false,, "Allow sinusoidal pressure condition to simulate ocean waves"))
 					((double, Sinus_Amplitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied"))
 					((double, Sinus_Average, 0,,"Pressure value (average) when sinusoidal pressure is applied"))
@@ -100,11 +99,8 @@
 					((double,viscosity,1.0,,"viscosity of fluid"))
 					((Real,loadFactor,1.1,,"Load multiplicator for oedometer test"))
 					((double, K, 0,, "Permeability of the sample"))
-					((double, MaxPressure, 0,, "Maximal value of water pressure within the sample - Case of consolidation test"))
-					((double, currentStress, 0,, "Current value of axial stress"))
-					((double, currentStrain, 0,, "Current value of axial strain"))
-					((int, intervals, 30,, "Number of layers for computation average fluid pressure profiles to build consolidation curves"))
 					((int, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso"))
+// 					((std::string,key,"",,"A string appended at the output files, use it to name simulations."))
 					((double, V_d, 0,,"darcy velocity of fluid in sample"))
 					((bool, Flow_imposed_TOP_Boundary, true,, "if false involve pressure imposed condition"))
 					((bool, Flow_imposed_BOTTOM_Boundary, true,, "if false involve pressure imposed condition"))
@@ -154,7 +150,8 @@
 					.def("fluidForce",&FlowEngine::fluidForce,(python::arg("Id_sph")),"Return the fluid force on sphere Id_sph.")
 					.def("fluidShearForce",&FlowEngine::fluidShearForce,(python::arg("Id_sph")),"Return the viscous shear force on sphere Id_sph.")
 					.def("setBoundaryVel",&FlowEngine::setBoundaryVel,(python::arg("vel")),"Change velocity on top boundary.")
-					.def("PressureProfile",&FlowEngine::PressureProfile,"Measure pore pressure in 6 equally-spaced points along the height of the sample")
+					.def("PressureProfile",&FlowEngine::PressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
+					.def("MeasurePressure",&FlowEngine::MeasurePressure,(python::arg("posX"),python::arg("posY"),python::arg("posZ")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
 					)
 		DECLARE_LOGGER;
 };