← Back to team overview

yade-dev team mailing list archive

  1. yade-dev team
  2. Mailing list archive
  3. Message #07747

[Branch ~yade-dev/yade/trunk] Rev 2881: - Removed correction on smallest permeabilities

  Thread Previous • Date Previous • • Thread Next 
------------------------------------------------------------
revno: 2881
committer: Emanuele Catalano <catalano@xxxxxxxxxxx
branch nick: yade
timestamp: Thu 2011-07-07 15:21:17 +0200
message:
  - Removed correction on smallest permeabilities
  - Porosity computation neglects boundaries
  - New functions ImposeFlux() and GetFlux() for one-point fluid injection
modified:
  lib/triangulation/FlowBoundingSphere.cpp
  lib/triangulation/FlowBoundingSphere.hpp
  lib/triangulation/Network.cpp
  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-05-03 20:03:21 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2011-07-07 13:21:17 +0000
@@ -499,20 +499,20 @@
 			}
 		cell->info().isvisited=!ref;
 	}
-// 	if (DEBUG_OUT) {
-// 		cout << "Facet scheme" <<endl;
-// 		Vecteur TotalForce = nullVect;
-// 		for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
-// 			if (!v->info().isFictious) {
-// 				TotalForce = TotalForce + v->info().forces;
-// 				cout << "real_id = " << v->info().id() << " force = " << v->info().forces << endl;
-// 			} else {
-// 				if (boundary(v->info().id()).flowCondition==1) TotalForce = TotalForce + v->info().forces;
-// 				cout << "fictious_id = " << v->info().id() << " force = " << v->info().forces << endl;
-// 			}
-// 		}
-// 		cout << "TotalForce = "<< TotalForce << endl;
-// 	}
+	if (DEBUG_OUT) {
+		cout << "Facet scheme" <<endl;
+		Vecteur TotalForce = nullVect;
+		for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+			if (!v->info().isFictious) {
+				TotalForce = TotalForce + v->info().forces;
+				cout << "real_id = " << v->info().id() << " force = " << v->info().forces << endl;
+			} else {
+				if (boundary(v->info().id()).flowCondition==1) TotalForce = TotalForce + v->info().forces;
+				cout << "fictious_id = " << v->info().id() << " force = " << v->info().forces << endl;
+			}
+		}
+		cout << "TotalForce = "<< TotalForce << endl;
+	}
 }
 
 void FlowBoundingSphere::ComputeFacetForcesWithCache()
@@ -582,10 +582,10 @@
 		{
 			if (!v->info().isFictious) {
 				TotalForce = TotalForce + v->info().forces;
-// 				if (DEBUG_OUT) cout << "real_id = " << v->info().id() << " force = " << v->info().forces << endl;
+				if (DEBUG_OUT) cout << "real_id = " << v->info().id() << " force = " << v->info().forces << endl;
 			} else {
 				if (boundary(v->info().id()).flowCondition==1) TotalForce = TotalForce + v->info().forces;
-// 				if (DEBUG_OUT) cout << "fictious_id = " << v->info().id() << " force = " << v->info().forces << endl;
+				if (DEBUG_OUT) cout << "fictious_id = " << v->info().id() << " force = " << v->info().forces << endl;
 			}
 		}
 		cout << "TotalForce = "<< TotalForce << endl;
@@ -677,6 +677,7 @@
 		if (new_cell->info().Pcondition) continue;
                 old_cell = Tri.locate((Point) new_cell->info());
                 new_cell->info().p() = old_cell->info().p();
+// 		new_cell->info().dv() = old_cell->info().dv();
         }
 	Tes.Clear();
 }
@@ -804,8 +805,8 @@
 			}
 		}
 		cell->info().isvisited = !ref;
-		cell->info().s = cell->info().s/volume_sub_pore;
-		volume_sub_pore = 0.f;
+		if(permeability_map){cell->info().s = cell->info().s/volume_sub_pore;
+		volume_sub_pore = 0.f;}
 	}
 	if (DEBUG_OUT) cout<<"surfneg est "<<surfneg<<endl;
 	meanK /= pass;
@@ -826,7 +827,8 @@
 			neighbour_cell = cell->neighbor(j);
 			if (!Tri.is_infinite(neighbour_cell) && neighbour_cell->info().isvisited==ref) {
 				pass++;
-				(cell->info().k_norm())[j] = max(MINK_DIV_KMEAN*meanK ,min((cell->info().k_norm())[j], maxKdivKmean*meanK));
+				(cell->info().k_norm())[j] = min((cell->info().k_norm())[j], maxKdivKmean*meanK);
+// 				(cell->info().k_norm())[j] = max(MINK_DIV_KMEAN*meanK ,min((cell->info().k_norm())[j], maxKdivKmean*meanK));
 				(neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]=(cell->info().k_norm())[j];
 // 				cout<<(cell->info().k_norm())[j]<<endl;
 // 				kFile << (cell->info().k_norm())[j] << endl;
@@ -1171,26 +1173,28 @@
   Tesselation::VCell_iterator cell_up_end = Tri.incident_cells(T[currentTes].vertexHandles[y_max_id],cells_it);
   for (Tesselation::VCell_iterator it = tmp_cells.begin(); it != cell_up_end; it++)
   {
-    Cell_handle& cell = *it;{for (int j2=0; j2<4; j2++) {/*if (!cell->neighbor(j2)->info().Pcondition)*/{
-      if ((cell->neighbor(j2)->info().p()!=cell->neighbor(j2)->info().p()) && Tri.is_infinite(cell->neighbor(j2))) cout << "oooooooooooooooh" << endl;
-    Q1 = Q1 + (cell->neighbor(j2)->info().k_norm())[Tri.mirror_index(cell, j2)]* (cell->neighbor(j2)->info().p()-cell->info().p());
-    cellQ1+=1;
-    p_out_max = std::max(cell->neighbor(j2)->info().p(), p_out_max);
-    p_out_min = std::min(cell->neighbor(j2)->info().p(), p_out_min);
-    p_out_moy += cell->neighbor(j2)->info().p();}
-  }}}
+    Cell_handle& cell = *it;
+    for (int j2=0; j2<4; j2++) {
+      if (!cell->neighbor(j2)->info().Pcondition){
+	Q1 = Q1 + (cell->neighbor(j2)->info().k_norm())[Tri.mirror_index(cell, j2)]* (cell->neighbor(j2)->info().p()-cell->info().p());
+	cellQ1+=1;
+	p_out_max = std::max(cell->neighbor(j2)->info().p(), p_out_max);
+	p_out_min = std::min(cell->neighbor(j2)->info().p(), p_out_min);
+	p_out_moy += cell->neighbor(j2)->info().p();}
+  }}
 
   Tesselation::VCell_iterator cell_down_end = Tri.incident_cells(T[currentTes].vertexHandles[y_min_id],cells_it);
   for (Tesselation::VCell_iterator it = tmp_cells.begin(); it != cell_down_end; it++)
   {
-    Cell_handle& cell = *it;{for (int j2=0; j2<4; j2++) {/*if (!cell->neighbor(j2)->info().Pcondition)*/{
-      if ((cell->neighbor(j2)->info().p()!=cell->neighbor(j2)->info().p()) && Tri.is_infinite(cell->neighbor(j2))) cout << "oooooooooooooooh2" << endl;
-    Q2 = Q2 + (cell->neighbor(j2)->info().k_norm())[Tri.mirror_index(cell, j2)]* (cell->info().p()-cell->neighbor(j2)->info().p());
-    cellQ2+=1;
-    p_in_max = std::max(cell->neighbor(j2)->info().p(), p_in_max);
-    p_in_min = std::min(cell->neighbor(j2)->info().p(), p_in_min);
-    p_in_moy += cell->neighbor(j2)->info().p();}
-  }}}
+    Cell_handle& cell = *it;
+    for (int j2=0; j2<4; j2++){
+      if (!cell->neighbor(j2)->info().Pcondition){
+	Q2 = Q2 + (cell->neighbor(j2)->info().k_norm())[Tri.mirror_index(cell, j2)]* (cell->info().p()-cell->neighbor(j2)->info().p());
+	cellQ2+=1;
+	p_in_max = std::max(cell->neighbor(j2)->info().p(), p_in_max);
+	p_in_min = std::min(cell->neighbor(j2)->info().p(), p_in_min);
+	p_in_moy += cell->neighbor(j2)->info().p();}
+  }}
 
 	if (DEBUG_OUT){
 	cout << "the maximum superior pressure is = " << p_out_max << " the min is = " << p_out_min << endl;
@@ -1201,8 +1205,8 @@
         cout << "celle comunicanti in alto = " << cellQ1 << endl;}
 
         double density = 1;
-        double viscosity = VISCOSITY;
-        double gravity = 9.80665;
+        double viscosity = 1.0;
+        double gravity = 1;
         double Vdarcy = Q1/Section;
 	double DeltaP = abs(P_Inf-P_Sup);
 	double DeltaH = DeltaP/ (density*gravity);
@@ -1422,11 +1426,42 @@
         }
 }
 
+// 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];
+// }
+
 double FlowBoundingSphere::PressureProfile(char *filename, Real& time, int& intervals)
 {
 	RTriangulation& Tri = T[currentTes].Triangulation();
 	vector<double> Pressures;
-
+	
 	/** CONSOLIDATION CURVES **/
         Cell_handle permeameter;
         int n=0, k=0;
@@ -1620,8 +1655,6 @@
     return tau * Edge_Surfaces[edge_id];
 }
 
-
-
 } //namespace CGT
 
 #endif //FLOW_ENGINE

=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2011-05-03 20:03:21 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2011-07-07 13:21:17 +0000
@@ -55,7 +55,7 @@
 		vector <Real> Edge_HydRad;
 		vector <Vector3r> Edge_normal;
 		vector <Vector3r> viscousShearForces;
-		
+	
 		void mplot ( char *filename);
 		void Localize();
 
@@ -93,7 +93,8 @@
 		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 ( char *filename, Real& time, int& intervals );
+// 		double PressureProfile (const char *filename, Real& time, int& intervals );
+double PressureProfile(char *filename, Real& time, int& intervals);
 
 		double dotProduct ( Vecteur x, Vecteur y );
 		double Compute_EffectiveRadius(Cell_handle cell, int j);
@@ -125,7 +126,6 @@
 		vector<Real> Average_Fluid_Velocity_On_Sphere(unsigned int Id_sph);
 		//Solver?
 		int useSolver;//(0 : GaussSeidel, 1 : TAUCS, 2 : PARDISO)
-
 };
 
 } //namespace CGT

=== modified file 'lib/triangulation/Network.cpp'
--- lib/triangulation/Network.cpp	2011-04-22 09:09:28 +0000
+++ lib/triangulation/Network.cpp	2011-07-07 13:21:17 +0000
@@ -87,8 +87,12 @@
 
 		Vsolid_tot += Vsolid1 + Vsolid2;
 		Vporale += Vtot - (Vsolid1 + Vsolid2);
-		V_porale_porosity += Vtot - (Vsolid1 + Vsolid2);
-		V_totale_porosity += Vtot;
+		
+		bool border=false;
+		for (int i=0;i<4;i++){
+		  if (cell->neighbor(i)->info().fictious()!=0) border=true;}
+		if (!border) {V_porale_porosity += Vtot - (Vsolid1 + Vsolid2);
+		    V_totale_porosity += Vtot;}
 
 		/**Vpore**/ return Vtot - (Vsolid1 + Vsolid2);
     }; break;

=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp	2011-05-18 07:06:02 +0000
+++ pkg/dem/FlowEngine.cpp	2011-07-07 13:21:17 +0000
@@ -104,11 +104,14 @@
 			char file [50];
 			mkdir("./Consol", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
 			string consol = "./Consol/"+flow->key+"%d_Consol";
+// 			const char* file = ;
 			const char* keyconsol = consol.c_str();
+			
 			sprintf(file, keyconsol, j);
 			char *g = file;
 			timingDeltas->checkpoint("Writing cons_files");
-			MaxPressure = flow->PressureProfile(g, time, intervals);
+// 			MaxPressure = flow->PressureProfile(consol.c_str(), time, intervals);
+			MaxPressure = flow->PressureProfile( g, time, intervals);
 
 			std::ofstream max_p("pressures.txt", std::ios::app);
 			max_p << j << " " << time << " " << MaxPressure << endl;
@@ -137,10 +140,7 @@
 		Update_Triangulation = false;}
 		
 	if (velocity_profile) /*flow->FluidVelocityProfile();*/flow->Average_Fluid_Velocity();
-	if (first && liquefaction){
-	  wall_up_y = flow->y_max;
-	  wall_down_y = flow->y_min;}
-	if (liquefaction) flow->Measure_Pore_Pressure(wall_up_y, wall_down_y);
+	
 	first=false;
 	timingDeltas->checkpoint("Ending");
 	
@@ -189,6 +189,22 @@
 	Update_Triangulation=true;
 }
 
+void FlowEngine::imposeFlux(Vector3r pos, Real flux)
+{
+	CGT::RTriangulation& Tri = flow->T[flow->currentTes].Triangulation();
+	double flux_base=0.f;
+	double perm_base=0.f;
+	CGT::Cell_handle cell = Tri.locate(CGT::Point(pos[0],pos[1],pos[2]));
+	for (int ngb=0;ngb<4;ngb++) {
+		if (!cell->neighbor(ngb)->info().Pcondition) {
+		  flux_base += cell->info().k_norm()[ngb]*(cell->neighbor(ngb)->info().p());
+		  perm_base += cell->info().k_norm()[ngb];}}
+	
+	flow->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),(flux_base-flux)/perm_base));
+	//force immediate update of boundary conditions
+	Update_Triangulation=true;
+}
+
 void FlowEngine::clearImposedPressure() { flow->imposedP.clear();}
 
 Real FlowEngine::getFlux(unsigned int cond) {
@@ -197,7 +213,7 @@
 	double flux=0;
 	CGT::Cell_handle cell= Tri.locate(flow->imposedP[cond].first);
 	for (int ngb=0;ngb<4;ngb++) {
-		if (!cell->neighbor(ngb)->info().Pcondition) flux+= cell->info().k_norm()[ngb]*(cell->info().p()-cell->neighbor(ngb)->info().p());}
+		/*if (!cell->neighbor(ngb)->info().Pcondition) */flux+= cell->info().k_norm()[ngb]*(cell->info().p()-cell->neighbor(ngb)->info().p());}
 	return flux;
 }
 
@@ -243,32 +259,18 @@
 
 	porosity = flow->V_porale_porosity/flow->V_totale_porosity;
 
-	if (first)
-	{
-		flow->TOLERANCE=Tolerance;
-		flow->RELAX=Relax;
-		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 );}
-		BoundaryConditions();
-		flow->Initialize_pressures( P_zero );
-  		if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Amplitude, Sinus_Average, 30);
-		
-
-	}
-	else
-	{
-		flow->TOLERANCE=Tolerance;
-		flow->RELAX=Relax;
-		if (Debug && compute_K) cout << "---------UPDATE PERMEABILITY VALUE--------------" << endl;
-		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 );}
-		BoundaryConditions();
-		flow->Initialize_pressures(P_zero);// FIXME : why, if we are going to interpolate after that?
-		flow->Interpolate (flow->T[!flow->currentTes], flow->T[flow->currentTes]);
- 		if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Amplitude, Sinus_Average, 30);
-	}
+	flow->TOLERANCE=Tolerance;flow->RELAX=Relax;
+	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 );}
+	
+	BoundaryConditions();
+	flow->Initialize_pressures(P_zero);
+	
+	if (!first) flow->Interpolate (flow->T[!flow->currentTes], flow->T[flow->currentTes]);
+	if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Amplitude, Sinus_Average, 30);
+	
 	Initialize_volumes();
 	if (viscousShear) flow->ComputeEdgesSurfaces();
 }
@@ -485,12 +487,14 @@
 	if (Debug) cout << "Updating volumes.............." << endl;
 	Real invDeltaT = 1/scene->dt;
 	CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
+	
 	double newVol, dVol;
 	eps_vol_max=0;
 	Real totVol=0; Real totDVol=0; Real totVol0=0; Real totVol1=0; Real totVol2=0; Real totVol3=0;
+	
 	for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
 		switch ( cell->info().fictious()) {
-			case ( 3 ): newVol= Volume_cell_triple_fictious ( cell ); totVol3+=newVol; break;
+			case ( 3 ) : newVol = Volume_cell_triple_fictious ( cell ); totVol3+=newVol; break;
 			case ( 2 ) : newVol = Volume_cell_double_fictious ( cell ); totVol2+=newVol; break;
 			case ( 1 ) : newVol = Volume_cell_single_fictious ( cell ); totVol1+=newVol; break;
 			case ( 0 ) : newVol = Volume_cell ( cell ); totVol0+=newVol; break;
@@ -502,6 +506,7 @@
 		eps_vol_max = max(eps_vol_max, abs(dVol/newVol));
 
 		cell->info().dv() = (!cell->info().Pcondition)?dVol*invDeltaT:0;
+// 		cell->info().dv() = dVol*invDeltaT;
 		cell->info().volume() = newVol;
 // 		if (Debug) cerr<<"v/dv : "<<cell->info().volume()<<" "<<cell->info().dv()<<" ("<<cell->info().fictious()<<")"<<endl;
 	}

=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp	2011-05-17 17:48:28 +0000
+++ pkg/dem/FlowEngine.hpp	2011-07-07 13:21:17 +0000
@@ -51,6 +51,7 @@
 		Real Volume_cell (CGT::Cell_handle cell);
 		void Oedometer_Boundary_Conditions();
 		void BoundaryConditions();
+		void imposeFlux(Vector3r pos, Real flux);
 		unsigned int imposePressure(Vector3r pos, Real p);
 		void setImposedPressure(unsigned int cond, Real p);
 		void clearImposedPressure();
@@ -65,6 +66,7 @@
 		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);}
 
 		virtual ~FlowEngine();
 
@@ -108,7 +110,7 @@
 					((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"))
-					((bool, liquefaction, false,,"Measure fluid pressure along the heigth of the sample"))
+// 					((bool, liquefaction, false,,"Measure fluid pressure along the heigth of the sample"))
 					((double, V_d, 0,,"darcy velocity of fluid in sample"))
 // 					((double, bottom_seabed_pressure,0,,"Fluid pressure measured at the bottom of the seabed on the symmetry axe"))
 					((bool, Flow_imposed_TOP_Boundary, true,, "if false involve pressure imposed condition"))
@@ -142,6 +144,7 @@
 					,,
 					timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
 					,
+					.def("imposeFlux",&FlowEngine::imposeFlux,(python::arg("pos"),python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
 					.def("imposePressure",&FlowEngine::imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
 					.def("setImposedPressure",&FlowEngine::setImposedPressure,(python::arg("cond"),python::arg("p")),"Set pressure value at the point of index cond.")
 					.def("clearImposedPressure",&FlowEngine::clearImposedPressure,"Clear the list of points with pressure imposed.")
@@ -152,6 +155,7 @@
 					.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")
 					)
 		DECLARE_LOGGER;
 };
  Thread Previous • Date Previous • • Thread Next