yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2145
committer: Emanuele Catalano <ecatalano@r2balme>
branch nick: trunk
timestamp: Thu 2010-04-15 12:41:33 +0200
message:
  - Flow Code Maintenance
  - Introduced local permeabilities' correction agent
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-04-08 17:49:39 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2010-04-15 10:41:33 +0000
@@ -33,7 +33,7 @@
 
 const double ONE_THIRD = 1.0/3.0;
 //! Use this factor, or minLength, to reduce max permeability values (see usage below))
-const double MAXK_DIV_KMEAN = 5000;
+const double MAXK_DIV_KMEAN = 1;
 const double minLength = 0.20;//percentage of mean rad
 
 //! Factors including the effect of 1/2 symmetry in hydraulic radii
@@ -75,6 +75,8 @@
 	TOLERANCE = 1e-06;
 	RELAX = 1.9;
 	ks=0;
+	meanK_LIMIT= false;
+	meanK_STAT = false; K_opt_factor=0;
 }
 
 void FlowBoundingSphere::Compute_Action()
@@ -157,6 +159,7 @@
         /** PERMEABILITY **/
         /** START PERMEABILITY CALCULATION**/
         k_factor = 1;
+	
         Compute_Permeability();
         clock.top("Compute_Permeability");
         /** END PERMEABILITY CALCULATION**/
@@ -365,7 +368,7 @@
 
 //         std::ofstream oFile( "Hydraulic_Radius",std::ios::out);
         std::ofstream kFile ( "LocalPermeabilities" ,std::ios::app );
-        Real meanK=0, k_moy;
+	Real meanK=0, STDEV=0;
         Real infiniteK=1e10;
         for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
                 p1 = cell->info();
@@ -401,33 +404,64 @@
 //     (cell->info().facetSurf())[j]= k*n;
                                 (neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]= k*k_factor;
 				
-				meanK += 1/k;
-				
+				meanK += (cell->info().k_norm())[j];
+				kFile << ( cell->info().k_norm() )[j] << endl;
 //     (neighbour_cell->info().facetSurf())[Tri.mirror_index(cell, j)]= (-k) *n;
                         }
                         //    else if ( Tri.is_infinite ( neighbour_cell )) k = 0;//connection to an infinite cells
                 }
                 cell->info().isvisited = !ref;
-
-//                 for (int y=0;y<4;y++) cout << "Permeability " << y << " = " << (cell->info().k_norm())[y] << endl;
-		for ( int y=0;y<4;y++ ) kFile << ( cell->info().k_norm() ) [y] << endl;
         }
+	meanK /= pass;
+	if (DEBUG_OUT) cout << "PassCompK = " << pass << endl;
 
         // A loop to reduce K maxima, needs a better equation : the mean value is influenced by the very big K
-	kFile << "----------reduction reduction reduction---------" << endl;
-        meanK /= pass;
-	k_moy = 1/meanK;
+
         Real maxKdivKmean = MAXK_DIV_KMEAN;
+	ref = Tri.finite_cells_begin()->info().isvisited; pass=0;
         for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
                 for (int j=0; j<4; j++) {
-//    if (cell->info().k_norm()[j]>maxKdivKmean*meanK) cerr <<"Adjusting permeability : " <<cell->info().k_norm()[j]<<" > "<<maxKdivKmean<<" * "<<meanK<<endl;
-//                         (cell->info().k_norm())[j] = min(cell->info().k_norm()[j], maxKdivKmean*meanK);
-			(cell->info().k_norm())[j] = min(cell->info().k_norm()[j], k_moy);
-			kFile << (cell->info().k_norm())[j] << endl;
+			neighbour_cell = cell->neighbor(j);
+			if (!Tri.is_infinite(neighbour_cell) && neighbour_cell->info().isvisited==ref){
+			pass++;
+			if (meanK_LIMIT) {
+				kFile << "--Correction--MEANKTHRESHOLD--" << endl;
+				(cell->info().k_norm())[j] = min(cell->info().k_norm()[j], maxKdivKmean*meanK);
+				kFile << (cell->info().k_norm())[j] << endl;}
+				STDEV += pow(((cell->info().k_norm())[j]-meanK),2);
+			}
                 }
         }
-        if (DEBUG_OUT) cout << "POS = " << POS << " NEG = " << NEG << " pass = " << pass << endl;
+	STDEV = sqrt(STDEV/pass);
+	
+	if (DEBUG_OUT) cout << "PassKcorrect = " << pass << endl;
+	if (DEBUG_OUT) cout << "POS = " << POS << " NEG = " << NEG << " pass = " << pass << endl;
+	
+	if (meanK_STAT)
+	{
+		cout << "STATISTIC K" << endl;
+		double k_min = 0;
+		double k_max = meanK + K_opt_factor*STDEV;
 
+		cout << "Kmoy = " << meanK << " Standard Deviation = " << STDEV << endl;
+		cout << "kmin = " << k_min << " kmax = " << k_max << endl;
+	
+		ref = Tri.finite_cells_begin()->info().isvisited;
+		pass=0;
+		for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+			for (int j=0; j<4; j++) {neighbour_cell = cell->neighbor(j);
+				if (!Tri.is_infinite(neighbour_cell) && neighbour_cell->info().isvisited==ref){
+					pass+=1;
+// 					if ((cell->info().k_norm())[j]<k_min) 
+// 					{(cell->info().k_norm())[j]=k_min;(neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]= (cell->info().k_norm())[j];}
+					if ((cell->info().k_norm())[j]>k_max)
+					{(cell->info().k_norm())[j]=k_max;(neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]= (cell->info().k_norm())[j];}
+				}
+			}cell->info().isvisited=!ref;
+		}
+		if (DEBUG_OUT) cout << "PassKopt = " << pass << endl;
+	}
+	
         Finite_vertices_iterator vertices_end = Tri.finite_vertices_end();
         Real Vgrains = 0;
         int grains=0;
@@ -532,6 +566,17 @@
                         }
                 cell->info().isvisited=!ref;
         }
+	cout << "tetrahedral 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;
+		} 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;
+		}
+	}
+	cout << "TotalForce = "<< TotalForce << endl;
 }
 
 void FlowBoundingSphere::Compute_Forces()

=== modified file 'lib/triangulation/FlowBoundingSphere.h'
--- lib/triangulation/FlowBoundingSphere.h	2010-04-07 15:03:09 +0000
+++ lib/triangulation/FlowBoundingSphere.h	2010-04-15 10:41:33 +0000
@@ -40,6 +40,8 @@
 		double TOLERANCE;
 		double RELAX;
 		double ks; //Hydraulic Conductivity
+		bool meanK_LIMIT, meanK_STAT;
+		double K_opt_factor;
 		
 		Boundary boundaries [6];
 		short id_offset;

=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp	2010-04-08 17:49:39 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp	2010-04-15 10:41:33 +0000
@@ -48,8 +48,7 @@
 		if ( current_state==3 )
 		{
 			if ( first || Update_Triangulation ) { Build_Triangulation( P_zero );}
-// 			else
-// 			{
+
 				timingDeltas->checkpoint("Triangulating");
 				
 				UpdateVolumes ( );
@@ -73,8 +72,8 @@
 			
 // 				flow->MGPost(flow->T[flow->currentTes].Triangulation());
 
-				flow->tess_based_force=tess_based_force;
-				flow->Compute_Forces ( );
+				flow->Compute_Forces();
+// 				flow->ComputeTetrahedralForces();
 			
 				timingDeltas->checkpoint("Compute_Forces");
 
@@ -107,14 +106,11 @@
 				MaxPressure = flow->PermeameterCurve(flow->T[flow->currentTes].Triangulation(), g, time, intervals);
 				
 				if ( Omega::instance().getCurrentIteration() % PermuteInterval == 0 )
-// 				if (Retriangulation)
 				{ Update_Triangulation = true; }
 				
 				timingDeltas->checkpoint("Storing Max Pressure");
 				
 				first=false;
-				
-// 			}
 		}
 	}
 }
@@ -163,13 +159,16 @@
 	
 	flow->Localize ();
 
+	flow->meanK_LIMIT = meanK_correction;
+	flow->meanK_STAT = meanK_opt;
 	flow->Compute_Permeability ();
 
 	if (first)
 	{
 		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-09; K = flow->Sample_Permeability ( flow->T[flow->currentTes].Triangulation(), 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->T[flow->currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
+		Oedometer_Boundary_Conditions();
 		flow->Initialize_pressures( P_zero );
 		flow->TOLERANCE=Tolerance;
 		flow->RELAX=Relax;
@@ -177,14 +176,13 @@
 	else 
 	{
 		cout << "---------UPDATE PERMEABILITY VALUE--------------" << endl;
-		if (compute_K) {flow->TOLERANCE=1e-09; K = flow->Sample_Permeability ( flow->T[flow->currentTes].Triangulation(), 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->T[flow->currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
 		flow->TOLERANCE=Tolerance;
 		flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
 		Update_Triangulation=!Update_Triangulation;
-		Retriangulation=false;
+		Oedometer_Boundary_Conditions();
 	}
 
-	Oedometer_Boundary_Conditions();
 	Initialize_volumes ( );
 	flow->DisplayStatistics ();
 }

=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.hpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.hpp	2010-04-08 17:49:39 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.hpp	2010-04-15 10:41:33 +0000
@@ -24,6 +24,7 @@
 		int current_state;
 		Real wall_thickness;
 		bool Update_Triangulation;
+		bool currentTes;
 		
 		void Triangulate ();
 		void AddBoundary ();
@@ -46,21 +47,21 @@
 					((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"))
+// 					((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"))
+					((bool,meanK_correction,false,"Local permeabilities' correction through meanK threshold"))
+					((bool,meanK_opt,false,"Local permeabilities' correction through an optimized threshold"))
 					((double,permeability_factor,1.0,"a permability multiplicator"))
 					((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"))
 					((double, currentStress, 0, "Current value of axial stress"))
 					((double, currentStrain, 0, "Current value of axial strain"))
-					((bool, Retriangulation, 0, "Enable/Disable sample retriangulation"))
-					((int, intervals, 30, "Number of layers for pressure measurements"))
-					((bool,tess_based_force,true,"true=force computation based on tessalation, false=force computation based on triangulation")),
+					((int, intervals, 30, "Number of layers for pressure measurements")),
 					timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas));
 		DECLARE_LOGGER;
 };