yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2848
committer: Emanuele Catalano <catalano@xxxxxxxxxxx
branch nick: yade
timestamp: Tue 2011-05-03 12:33:29 +0200
message:
  - Reordering new function for viscous force calculation
  - Added computation of average cell permeability for visualization
modified:
  lib/triangulation/FlowBoundingSphere.cpp
  lib/triangulation/FlowBoundingSphere.hpp
  lib/triangulation/def_types.h
  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-02 09:17:49 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2011-05-03 10:33:29 +0000
@@ -92,6 +92,7 @@
 	OUTPUT_BOUDARIES_RADII = false;
 	RAVERAGE = false; /** if true use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
 	areaR2Permeability=true;
+	permeability_map = false;
 }
 
 void FlowBoundingSphere::ResetNetwork() {noCache=true;}
@@ -723,6 +724,8 @@
 	Real meanK=0, STDEV=0, meanRadius=0, meanDistance=0;
 	Real infiniteK=1e10;
 
+	double volume_sub_pore = 0.f;
+
 	for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
 		Point& p1 = cell->info();
 		for (int j=0; j<4; j++) {
@@ -760,10 +763,10 @@
 					cout << "INS-INS PROBLEM!!!!!!!" << endl;
 				}
 // 				Real h,d;
+				Real fluidArea=0;
 				int test=0;
 				if (distance!=0) {
 					if (minPermLength>0 && distance_correction) distance=max(minPermLength,distance);
-					Real fluidArea=0;
 					const Vecteur& Surfk = cell->info().facetSurfaces[j];
 					Real area = sqrt(Surfk.squared_length());
 					const Vecteur& crossSections = cell->info().facetSphereCrossSections[j];
@@ -792,9 +795,17 @@
 				meanDistance += distance;
 				meanRadius += radius;
 				meanK += (cell->info().k_norm())[j];
+				
+				if(permeability_map){
+				  Cell_handle c = cell;
+				  cell->info().s = cell->info().s + k*distance/fluidArea*Volume_Pore_VoronoiFraction (c,j);
+				  volume_sub_pore += Volume_Pore_VoronoiFraction (c,j);}
+				
 			}
 		}
 		cell->info().isvisited = !ref;
+		cell->info().s = cell->info().s/volume_sub_pore;
+		volume_sub_pore = 0.f;
 	}
 	if (DEBUG_OUT) cout<<"surfneg est "<<surfneg<<endl;
 	meanK /= pass;
@@ -1266,8 +1277,6 @@
 
 void FlowBoundingSphere::save_vtk_file()
 {
-	ComputeEdgesSurfaces();
-	
 	RTriangulation& Tri = T[currentTes].Triangulation();
         static unsigned int number = 0;
         char filename[80];
@@ -1300,11 +1309,18 @@
 // 	{if (!v->info().isFictious) vtkfile.write_data((v->info().forces)[0],(v->info().forces)[1],(v->info().forces)[2]);}
 // 	vtkfile.end_data();
 
+	if (permeability_map){
+	vtkfile.begin_data("Permeability",CELL_DATA,SCALARS,FLOAT);
+	for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
+		if (!cell->info().fictious()){vtkfile.write_data(cell->info().s);}
+	}
+	vtkfile.end_data();}
+	else{
 	vtkfile.begin_data("Pressure",CELL_DATA,SCALARS,FLOAT);
 	for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
 		if (!cell->info().fictious()){vtkfile.write_data(cell->info().p());}
 	}
-	vtkfile.end_data();
+	vtkfile.end_data();}
 
 // 	vtkfile.begin_data("Velocity",POINT_DATA,VECTORS,FLOAT);
 // 	for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v)
@@ -1593,11 +1609,19 @@
     if (radius1<radius2)  Rh = d + 0.45 * radius1;
     else  Rh = d + 0.45 * radius2;
     Edge_HydRad. push_back(Rh);
-    cout<<"id1= "<<id1<<", id2= "<<id2<<", area= "<<area<<", R1= "<<radius1<<", R2= "<<radius2<<" x= "<<x<<", n= "<<n<<", Rh= "<<Rh<<endl;
+    if (DEBUG_OUT) cout<<"id1= "<<id1<<", id2= "<<id2<<", area= "<<area<<", R1= "<<radius1<<", R2= "<<radius2<<" x= "<<x<<", n= "<<n<<", Rh= "<<Rh<<endl;
     
   }
 }
 
+Vector3r FlowBoundingSphere::ComputeViscousForce(Vector3r deltaV, int edge_id)
+{
+    Vector3r tau = deltaV*VISCOSITY/Edge_HydRad[edge_id];
+    return tau * Edge_Surfaces[edge_id];
+}
+
+
+
 } //namespace CGT
 
 #endif //FLOW_ENGINE

=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2011-05-02 09:17:49 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2011-05-03 10:33:29 +0000
@@ -46,6 +46,7 @@
 		bool noCache;//flag for checking if cached values cell->unitForceVectors have been defined
 		vector<pair<Point,Real> > imposedP;
 		void initNewTri () {noCache=true; /*isLinearSystemSet=false; areCellsOrdered=false;*/}//set flags after retriangulation
+		bool permeability_map;
 
 		bool computeAllCells;//exececute computeHydraulicRadius for all facets and all spheres (double cpu time but needed for now in order to define crossSections correctly)
 		double K_opt_factor;
@@ -106,6 +107,7 @@
 		void GenerateVoxelFile ( );
 		
 		void ComputeEdgesSurfaces();
+		Vector3r ComputeViscousForce(Vector3r deltaV, int edge_id);
 // 		Real ComputeVFacetArea(Finite_edges_iterator ed_it);
 
 		RTriangulation& Build_Triangulation ( Real x, Real y, Real z, Real radius, unsigned const id );

=== modified file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h	2011-04-08 11:15:34 +0000
+++ lib/triangulation/def_types.h	2011-05-03 10:33:29 +0000
@@ -87,6 +87,7 @@
 		isFictious=false; Pcondition = false; isInferior = false; isSuperior = false; isLateral = false; isvisited = false; isExternal=false;
 		index=0;
 		volumeSign=0;
+		s=0;
 	}	
 
 	double inv_sum_k;

=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp	2011-05-02 09:17:49 +0000
+++ pkg/dem/FlowEngine.cpp	2011-05-03 10:33:29 +0000
@@ -80,61 +80,16 @@
 		else flow->ComputeFacetForcesWithCache();
 		timingDeltas->checkpoint("Compute_Forces");
 		
-		if (viscousShear){
 		///Application of vicscous forces
-		int length_ids = flow -> Edge_ids.size();
-		cout <<"length_ids est "<<length_ids<< endl;
-		vector <Vector3r> Edge_force;
-		for (int i=0;i<length_ids;i++)
-		{
-// 		  cout <<"length_ids est "<<length_ids<< endl;
-		  cout<<"Application of vicscous forces"<<endl;
-		  
-		  const shared_ptr<Body>& sph1 = Body::byId( flow->Edge_ids[i].first, scene );
-		  const shared_ptr<Body>& sph2 = Body::byId( flow->Edge_ids[i].second, scene );
-		  Vector3r deltaV = (sph2->state->vel - sph1->state->vel) - (flow->Edge_normal[i].dot(sph2->state->vel - sph1->state->vel))*flow->Edge_normal[i];
-		  Vector3r tau = deltaV*viscosity/flow->Edge_HydRad[i];
-		  Vector3r visc_f = tau * flow->Edge_Surfaces[i];
-		  cout<<"la force visqueuse entre "<<flow->Edge_ids[i].first<<" et "<<flow->Edge_ids[i].second<<"est "<<visc_f<< endl;
-		  Edge_force.push_back(visc_f);
-	
-		///Compute total force
-      
-		 
-		  int id1 = flow->Edge_ids[i].first;
-		  int id2 = flow->Edge_ids[i].second;
-// 		  
-		  scene->forces.addForce(id1,Edge_force[i]);
-		  scene->forces.addForce(id2,-Edge_force[i]);
-		  
-		  scene->forces.sync();
-		  Vector3r F1=scene->forces.getForce(id1);
-		  Vector3r F2=scene->forces.getForce(id2);
-		  cout<<"la force totale sur "<<id1<<" est "<<F1<< "et la force totale sur "<<id2<<" est "<<F2;
 		
-		  cout<<"END: Application of vicscous forces"<<endl;
-		  
-		}
-		}
+		if (viscousShear) ApplyViscousForces();
 
-		///End Compute flow and forces
 		CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
-		Vector3r f;
-		int id;
 		for (CGT::Finite_vertices_iterator V_it = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); V_it !=  vertices_end; V_it++)
 		{
-			id = V_it->info().id();
-			for (int y=0;y<3;y++) f[y] = (V_it->info().forces)[y];
-// 			if (id==id_sphere)
-// 			  id_force = f;
-// 			  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++ ){
-// 			    for (int j=0;j<4;j++){
-// 			      if (cell->vertex(j)->info().id() == id_sphere){
-// 				cout << endl << "Cell 
-			    
-			scene->forces.addForce(id, f);
+			scene->forces.addForce(V_it->info().id(), Vector3r ((V_it->info().forces)[0],V_it->info().forces[1],V_it->info().forces[2]));
 		}
+		///End Compute flow and forces
 		timingDeltas->checkpoint("Applying Forces");
 
 		Real time = scene->time;
@@ -260,6 +215,7 @@
 
 	flow->meanK_LIMIT = meanK_correction;
 	flow->meanK_STAT = meanK_opt;
+	flow->permeability_map = permeability_map;
 	flow->Compute_Permeability ( );
 
 	porosity = flow->V_porale_porosity/flow->V_totale_porosity;
@@ -636,6 +592,29 @@
 	return volume;
 }
 
+void FlowEngine::ApplyViscousForces()
+{
+  flow->ComputeEdgesSurfaces();
+  if (Debug) cout << "Application of viscous forces" << endl;
+  if (Debug) cout << "Number of edges = " << flow->Edge_ids.size() << endl;
+  vector <Vector3r> Edge_force;
+  Edge_force.resize(flow->Edge_ids.size());
+  for (int i=0; i<(int)flow->Edge_ids.size(); i++)
+  {
+    const shared_ptr<Body>& sph1 = Body::byId( flow->Edge_ids[i].first, scene );
+    const shared_ptr<Body>& sph2 = Body::byId( flow->Edge_ids[i].second, scene );
+    Vector3r deltaV = (sph2->state->vel - sph1->state->vel) - (flow->Edge_normal[i].dot(sph2->state->vel - sph1->state->vel))*flow->Edge_normal[i];
+    Vector3r visc_f = flow->ComputeViscousForce(deltaV, i);
+    if (Debug) cout << "la force visqueuse entre " << flow->Edge_ids[i].first << " et " << flow->Edge_ids[i].second << "est " << visc_f << endl;
+    Edge_force.push_back(visc_f);
+    scene->forces.addForce(flow->Edge_ids[i].first,Edge_force[i]);
+    scene->forces.addForce(flow->Edge_ids[i].second,-Edge_force[i]);
+    scene->forces.sync();
+    if (Debug) cout<<"la force totale sur "<< flow->Edge_ids[i].first <<" est " << scene->forces.getForce(flow->Edge_ids[i].first) << "et la force totale sur "<< flow->Edge_ids[i].second <<" est " << scene->forces.getForce(flow->Edge_ids[i].second);
+    if (Debug) cout<<"END: Application of vicscous forces"<<endl;
+  }
+}
+
 YADE_PLUGIN ((FlowEngine));
 #endif //FLOW_ENGINE
 

=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp	2011-05-02 09:17:49 +0000
+++ pkg/dem/FlowEngine.hpp	2011-05-03 10:33:29 +0000
@@ -54,6 +54,7 @@
 		void imposePressure(Vector3r pos, Real p);
 		void clearImposedPressure();
 		void Average_real_cell_velocity();
+		void ApplyViscousForces();
 		Real getFlux(int cond);
 		void saveVtk() {flow->save_vtk_file();}
 		vector<Real> AvFlVelOnSph(int id_sph) {return flow->Average_Fluid_Velocity_On_Sphere(id_sph);}
@@ -70,6 +71,7 @@
 					((bool,isActivated,true,,"Activates Flow Engine"))
 					((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"))