yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2918
committer: Emanuele Catalano <catalano@xxxxxxxxxxx
branch nick: yade
timestamp: Wed 2011-09-07 11:46:39 +0200
message:
  - Add function to measure slice-averaged pressure values
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-23 10:24:56 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2011-09-07 09:46:39 +0000
@@ -450,6 +450,24 @@
 	
 }
 
+double FlowBoundingSphere::MeasureAveragedPressure(double Y)
+{
+  RTriangulation& Tri = T[currentTes].Triangulation();
+  double P_ave = 0.f;
+  int n = 0;
+  double Ry = (y_max-y_min)/30;
+  double Rx = (x_max-x_min)/30;
+  double Rz = (z_max-z_min)/30;
+  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+=Tri.locate(Point(X, Y, Z))->info().p();
+	  n++;
+	}
+  }
+  P_ave/=n;
+  return P_ave;
+}
+
 void FlowBoundingSphere::ComputeFacetForces()
 {
 	RTriangulation& Tri = T[currentTes].Triangulation();
@@ -1442,7 +1460,7 @@
 	}
 }
 
-double FlowBoundingSphere::Sample_Permeability(double& x_Min,double& x_Max ,double& y_Min,double& y_Max,double& z_Min,double& z_Max, string key)
+double FlowBoundingSphere::Sample_Permeability(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);
         double DeltaY = y_Max-y_Min;
@@ -1454,8 +1472,7 @@
         Initialize_pressures( P_zero );
 	GaussSeidel();
 
-        char *kk;
-        kk = (char*) key.c_str();
+	const char *kk = "Permeability";
         return Permeameter(boundary(y_min_id).value, boundary(y_max_id).value, Section, DeltaY, kk);
 }
 

=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2011-08-23 13:22:46 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2011-09-07 09:46:39 +0000
@@ -92,7 +92,7 @@
 		void Dessine_Short_Tesselation ( Vue3D &Vue, Tesselation &Tes );
 #endif
 		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 Sample_Permeability( double& x_Min,double& x_Max ,double& y_Min,double& y_Max,double& z_Min,double& z_Max);
 		double Compute_HydraulicRadius (Cell_handle cell, int j );
 
 		double dotProduct ( Vecteur x, Vecteur y );
@@ -121,6 +121,7 @@
 		void ApplySinusoidalPressure(RTriangulation& Tri, double Amplitude, double Average_Pressure, double load_intervals);
 		bool isOnSolid  (double X, double Y, double Z);
 		void MeasurePressureProfile(double Wall_up_y, double Wall_down_y);
+		double MeasureAveragedPressure(double 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-23 10:24:56 +0000
+++ pkg/dem/FlowEngine.cpp	2011-09-07 09:46:39 +0000
@@ -214,7 +214,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) {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) {K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max );}
 	
 	BoundaryConditions();
 	flow->Initialize_pressures(P_zero);

=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp	2011-08-23 13:22:46 +0000
+++ pkg/dem/FlowEngine.hpp	2011-09-07 09:46:39 +0000
@@ -63,6 +63,7 @@
 		void setBoundaryVel(Vector3r vel) {topBoundaryVelocity=vel; Update_Triangulation=true;}
 		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);}
+		double MeasureAveragedPressure(double posY){return flow->MeasureAveragedPressure(posY);}
 
 		virtual ~FlowEngine();
 
@@ -153,6 +154,7 @@
 					.def("setBoundaryVel",&FlowEngine::setBoundaryVel,(python::arg("vel")),"Change velocity on top boundary.")
 					.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]")
+					.def("MeasureAveragedPressure",&FlowEngine::MeasureAveragedPressure,(python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
 					)
 		DECLARE_LOGGER;
 };