← Back to team overview

yade-dev team mailing list archive

[Branch ~yade-dev/yade/trunk] Rev 2284: - Fixed errors in walls' identification

 

------------------------------------------------------------
revno: 2284
committer: ecatalano <ecatalano@dt-rv019>
branch nick: trunk
timestamp: Tue 2010-06-08 19:16:42 +0200
message:
  - Fixed errors in walls' identification
modified:
  lib/triangulation/FlowBoundingSphere.cpp
  lib/triangulation/FlowBoundingSphere.h
  pkg/dem/Engine/PartialEngine/FlowEngine.cpp


--
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-06-04 12:07:46 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2010-06-08 17:16:42 +0000
@@ -1601,7 +1601,7 @@
 	{cell->info().p() = P_zero;cell->info().dv()=0;}
 
         for (int bound=0; bound<6;bound++) {
-                int id = boundsIds[bound];
+                int& id = *boundsIds[bound];
                 Boundary& bi = boundary(id);
                 if (!bi.flowCondition) {
                         Tesselation::Vector_Cell tmp_cells;
@@ -1987,6 +1987,7 @@
 
         char *kk;
         kk = (char*) key.c_str();
+
         return Permeameter(Tri, boundary(y_min_id).value, boundary(y_max_id).value, Section, DeltaY, kk);
 }
 
@@ -1997,7 +1998,6 @@
         Corner_min = Point(x_min, y_min, z_min);
         Corner_max = Point(x_max, y_max, z_max);
         Real min_coord = min(Extents[0],min(Extents[1],Extents[2]));
-	
         int coord=0;
         if (min_coord==Extents[0]) {
                 coord=0;
@@ -2059,27 +2059,36 @@
 {
 	Tesselation& Tes = T[currentTes];
 	
-	y_min_id = Tes.Max_id() + 1; boundsIds[0]=y_min_id;
-        y_max_id = Tes.Max_id() + 2; boundsIds[1]=y_max_id;
-        x_min_id = Tes.Max_id() + 3; boundsIds[2]=x_min_id;
-        x_max_id = Tes.Max_id() + 4; boundsIds[3]=x_max_id;
-        z_min_id = Tes.Max_id() + 5; boundsIds[4]=z_min_id;
-        z_max_id = Tes.Max_id() + 6; boundsIds[5]=z_max_id;
+	y_min_id = Tes.Max_id() + 1;
+        boundsIds[0]=&y_min_id;
+        y_max_id = Tes.Max_id() + 2;
+        boundsIds[1]=&y_max_id;
+        x_min_id = Tes.Max_id() + 3;
+        boundsIds[2]=&x_min_id;
+        x_max_id = Tes.Max_id() + 4;
+        boundsIds[3]=&x_max_id;
+        z_min_id = Tes.Max_id() + 5;
+        boundsIds[4]=&z_min_id;
+        z_max_id = Tes.Max_id() + 6;
+        boundsIds[5]=&z_max_id;
 	 
 	id_offset = Tes.Max_id() +1;//so that boundaries[vertex->id - offset] gives the ordered boundaries (also see function Boundary& boundary(int b))
 	
-	AddBoundingPlanes(y_min_id, y_max_id, x_min_id, x_max_id, z_min_id, z_max_id);
+	AddBoundingPlanes(true);
 }
 
-void FlowBoundingSphere::AddBoundingPlanes(int bottom_id, int top_id, int left_id, int right_id, int front_id, int back_id)
+void FlowBoundingSphere::AddBoundingPlanes(bool yade)
 {
         Tesselation& Tes = T[currentTes];
         Corner_min = Point(x_min, y_min, z_min);
         Corner_max = Point(x_max, y_max, z_max);
-	
-	y_min_id = bottom_id;y_max_id = top_id;x_min_id = left_id;x_max_id = right_id;z_min_id = front_id;z_max_id = back_id;
-	
-        boundsIds[0]= y_min_id;boundsIds[1]= y_max_id;boundsIds[2]= x_min_id;boundsIds[3]= x_max_id;boundsIds[4]= z_min_id;boundsIds[5]= z_max_id;
+
+        boundsIds[0]= &y_min_id;
+        boundsIds[1]= &y_max_id;
+        boundsIds[2]= &x_min_id;
+        boundsIds[3]= &x_max_id;
+        boundsIds[4]= &z_min_id;
+        boundsIds[5]= &z_max_id;
 
         Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), Corner_min.y()-FAR*(Corner_max.x()-Corner_min.x()), 0.5*(Corner_max.z()-Corner_min.z()), FAR*(Corner_max.x()-Corner_min.x()), y_min_id, true);
         boundaries[y_min_id-id_offset].p = Corner_min;

=== modified file 'lib/triangulation/FlowBoundingSphere.h'
--- lib/triangulation/FlowBoundingSphere.h	2010-06-04 12:07:46 +0000
+++ lib/triangulation/FlowBoundingSphere.h	2010-06-08 17:16:42 +0000
@@ -34,7 +34,7 @@
  		FlowBoundingSphere();
 		
 		int x_min_id, x_max_id, y_min_id, y_max_id, z_min_id, z_max_id;
-		int boundsIds [6];
+		int* boundsIds [6];
 		bool currentTes;
 		bool SLIP_ON_LATERALS;
 		double TOLERANCE;
@@ -45,13 +45,17 @@
 		int Iterations;
 		
 		Boundary boundaries [6];
-		int walls_id[6];
 		short id_offset;
  		Boundary& boundary (int b) {return boundaries[b-id_offset];}
 		
 		void mplot (RTriangulation& Tri, char *filename);
 		void Localize ();
 		void Compute_Permeability();
+		void AddBoundingPlanes();
+		
+		void AddBoundingPlanes(Real center[3], Real Extents[3], int id);
+// 		void AddBoundingPlanes(int y_min_id, int y_max_id, int x_min_id, int x_max_id, int z_min_id, int z_max_id);
+		void AddBoundingPlanes(bool yade); 
 		
 		void DisplayStatistics();
 		void GaussSeidel ( );
@@ -70,12 +74,7 @@
 		Real minPermLength; //min branch length for Poiseuille
 		
 		double P_SUP, P_INF, P_INS;
-		
 		void AddBoundingPlanes ( Tesselation& Tes, double x_Min,double x_Max ,double y_Min,double y_Max,double z_Min,double z_Max );
-		void AddBoundingPlanes(int y_min_id, int y_max_id, int x_min_id, int x_max_id, int z_min_id, int z_max_id);
-		void AddBoundingPlanes();
-		void AddBoundingPlanes(Real center[3], Real Extents[3], int id);
-		
 		void Compute_Action ( );
 		void Compute_Action ( int argc, char *argv[ ], char *envp[ ] );
 		void DisplayStatistics ( RTriangulation& Tri );

=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp	2010-06-04 12:07:46 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp	2010-06-08 17:16:42 +0000
@@ -28,8 +28,7 @@
 
 void FlowEngine::action ( )
 {
-	if (!flow) 
-	{flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);first=true;Update_Triangulation=false;}
+	if (!flow) {flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);first=true;Update_Triangulation=false;}
 	if ( !isActivated ) return;
 	else
 	{
@@ -74,8 +73,7 @@
 // 				flow->MGPost(flow->T[flow->currentTes].Triangulation());
 
 				flow->ComputeTetrahedralForces();
-//				flow->Compute_Forces();
-				
+			
 				timingDeltas->checkpoint("Compute_Forces");
 
 			///End Compute flow and forces
@@ -119,7 +117,7 @@
 				if ( Update_Triangulation ) { Build_Triangulation( );}
 				
 				timingDeltas->checkpoint("Storing Max Pressure");
-				
+
 				first=false;
 		}
 	}
@@ -173,7 +171,7 @@
 	flow->T[flow->currentTes].Compute();
 	
 	flow->Localize ();
-	flow->DisplayStatistics ();
+	flow->DisplayStatistics ();	
 
 	flow->meanK_LIMIT = meanK_correction;
 	flow->meanK_STAT = meanK_opt;
@@ -198,6 +196,7 @@
 		Update_Triangulation=!Update_Triangulation;
 		Oedometer_Boundary_Conditions();
 	}
+
 	Initialize_volumes ( );
 }
 
@@ -226,13 +225,19 @@
 			flow->y_max = max ( flow->y_max, center[1]-wall_thickness);
 			flow->z_min = min ( flow->z_min, center[2]+wall_thickness);
 			flow->z_max = max ( flow->z_max, center[2]-wall_thickness);
-
 		}
 	}
 	
 	flow->id_offset = flow->T[flow->currentTes].max_id+1;
 
-	flow->AddBoundingPlanes( triaxialCompressionEngine->wall_bottom_id, triaxialCompressionEngine->wall_top_id, triaxialCompressionEngine->wall_left_id, triaxialCompressionEngine->wall_right_id, triaxialCompressionEngine->wall_front_id, triaxialCompressionEngine->wall_back_id );
+	flow->y_min_id=triaxialCompressionEngine->wall_bottom_id;
+	flow->y_max_id=triaxialCompressionEngine->wall_top_id;
+	flow->x_max_id=triaxialCompressionEngine->wall_right_id;
+	flow->x_min_id=triaxialCompressionEngine->wall_left_id;
+	flow->z_min_id=triaxialCompressionEngine->wall_back_id;
+	flow->z_max_id=triaxialCompressionEngine->wall_front_id;
+
+	flow->AddBoundingPlanes(true);
 }
 
 void FlowEngine::Triangulate ()
@@ -340,7 +345,7 @@
 		}
 		else
 		{
-			b = cell->vertex ( y )->info().id();
+			b = cell->vertex ( y )->info().id()-flow->id_offset;
 			const shared_ptr<Body>& wll = Body::byId ( b , scene );
 			for ( int i=0;i<3;i++ ) Wall_point[i] = flow->boundaries[b].p[i];
 	Wall_point[flow->boundaries[b].coordinate] = wll->state->pos[flow->boundaries[b].coordinate]+(flow->boundaries[b].normal[flow->boundaries[b].coordinate])*wall_thickness;
@@ -362,11 +367,13 @@
 
 Real FlowEngine::Volume_cell_double_fictious ( CGT::Cell_handle cell)
 {
-// 	Real array_quattro[3] = {0, 0, 0};
-
 	Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0};
 	Real B[3]={0, 0, 0}, BS[3]={0, 0, 0}, BT[3]={0, 0, 0};
-	Real C[3]={0, 0, 0}, CS[3]={0, 0, 0}, CT[3]={0, 0, 0};
+	Real C[3]={0, 0, 0}, CS[3]={0, 0, 0}, CT[3]={0, 0, 0};	
+
+	//Real A[3], AS[3], AT[3];
+	//Real B[3], BS[3], BT[3];
+	//Real C[3], CS[3], CT[3];
 	int b[2];
 
 	Real Wall_point[2][3];
@@ -377,7 +384,7 @@
 	{
 		if ( cell->vertex ( g )->info().isFictious )
 		{
-			b[j] = cell->vertex ( g )->info().id();
+			b[j] = cell->vertex ( g )->info().id()-flow->id_offset;
 			const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
 			for ( int i=0;i<3;i++ ) Wall_point[j][i] = flow->boundaries[b[j]].p[i];
 			Wall_point[j][flow->boundaries[b[j]].coordinate] = wll->state->pos[flow->boundaries[b[j]].coordinate] +(flow->boundaries[b[j]].normal[flow->boundaries[b[j]].coordinate])*wall_thickness;
@@ -416,7 +423,8 @@
 
 Real FlowEngine::Volume_cell_triple_fictious ( CGT::Cell_handle cell)
 {
-	Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0}, AW[3]={0, 0, 0};
+	Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0}, AW[3]={0, 0, 0};	
+//Real A[3], AS[3], AT[3], AW[3];
 // 	CGT::Boundary b[3];
 	int b[3];
 	Real Wall_point[3][3];
@@ -426,7 +434,7 @@
 	{
 		if ( cell->vertex ( g )->info().isFictious )
 		{
-			b[j] = cell->vertex ( g )->info().id();
+			b[j] = cell->vertex ( g )->info().id()-flow->id_offset;
 			const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
 			for ( int i=0;i<3;i++ ) Wall_point[j][i] = flow->boundaries[b[j]].p[i];
 			Wall_point[j][flow->boundaries[b[j]].coordinate] = wll->state->pos[flow->boundaries[b[j]].coordinate]+(flow->boundaries[b[j]].normal[flow->boundaries[b[j]].coordinate])*wall_thickness;