yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2846
author: Donia Marzougui <marzougui.donia@xxxxxxxxxx>
committer: Emanuele Catalano <catalano@xxxxxxxxxxx
branch nick: yade
timestamp: Mon 2011-05-02 11:17:49 +0200
message:
  - Add viscous shear force definition (disabled by default) (Donia+Bruno)
modified:
  lib/triangulation/FlowBoundingSphere.cpp
  lib/triangulation/FlowBoundingSphere.hpp
  lib/triangulation/Tesselation.cpp
  lib/triangulation/Tesselation.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-04-27 15:15:04 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp	2011-05-02 09:17:49 +0000
@@ -695,8 +695,6 @@
 	Real m=(cross_product(v0-v1,v2-v1)).squared_length()/v1v2;
 	if (m<v0.weight()) {
 		Real d=2*sqrt((v0.weight()-m));
-// 		h=sqrt(v0.weight())-sqrt(m);
-// 		Real S0=0.25*M_PI*(d*d+4*h*h);
 		Real teta=2*acos(sqrt(m/v0.weight()));
 		return 0.5*(teta*v0.weight()-d*sqrt(m));//this is S0, we use crossSection to avoid computing an "asin"
 // 		return crossSection-m*d;
@@ -1268,6 +1266,8 @@
 
 void FlowBoundingSphere::save_vtk_file()
 {
+	ComputeEdgesSurfaces();
+	
 	RTriangulation& Tri = T[currentTes].Triangulation();
         static unsigned int number = 0;
         char filename[80];
@@ -1570,6 +1570,34 @@
 	cerr << "meanCmsVel "<<meanCmsVel/totCmsPoints<<" mean diff "<<diff/kk<<endl;
 }
 
+void FlowBoundingSphere::ComputeEdgesSurfaces()
+{
+  RTriangulation& Tri = T[currentTes].Triangulation();
+ 
+  Finite_edges_iterator ed_it;
+  for ( Finite_edges_iterator ed_it = Tri.finite_edges_begin(); ed_it!=Tri.finite_edges_end();ed_it++ )
+  {
+    Real Rh;
+    if (((ed_it->first)->vertex(ed_it->second)->info().isFictious) && ((ed_it->first)->vertex(ed_it->third)->info().isFictious)) continue;
+    int id1 = (ed_it->first)->vertex(ed_it->second)->info().id();
+    int id2 = (ed_it->first)->vertex(ed_it->third)->info().id();
+    double area = T[currentTes].ComputeVFacetArea(ed_it);
+    Edge_Surfaces.push_back(area);
+    Edge_ids.push_back(pair<int,int>(id1,id2));
+    double radius1 = sqrt((ed_it->first)->vertex(ed_it->second)->point().weight());
+    double radius2 = sqrt((ed_it->first)->vertex(ed_it->third)->point().weight());
+    Vecteur x = (ed_it->first)->vertex(ed_it->third)->point().point()- (ed_it->first)->vertex(ed_it->second)->point().point();
+    Vecteur n = x / sqrt(x.squared_length());
+    Edge_normal.push_back(Vector3r(n[0],n[1],n[2]));
+    double d = x*n - radius1 - radius2;
+    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;
+    
+  }
+}
+
 } //namespace CGT
 
 #endif //FLOW_ENGINE

=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2011-04-06 16:37:32 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2011-05-02 09:17:49 +0000
@@ -53,6 +53,11 @@
 
 		bool RAVERAGE;
 		int walls_id[6];
+		vector <double> Edge_Surfaces;
+		vector <pair<int,int> > Edge_ids;
+		vector <Real> Edge_HydRad;
+		vector <Vector3r> Edge_normal;
+		
 		void mplot ( char *filename);
 		void Localize();
 
@@ -99,6 +104,9 @@
 		vector<double> getConstrictions();
 
 		void GenerateVoxelFile ( );
+		
+		void ComputeEdgesSurfaces();
+// 		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/Tesselation.cpp'
--- lib/triangulation/Tesselation.cpp	2010-11-18 20:01:50 +0000
+++ lib/triangulation/Tesselation.cpp	2011-05-02 09:17:49 +0000
@@ -12,9 +12,10 @@
 //  std::cout << "Tesselation(void)" << std::endl;
 	Tri = new RTriangulation;
 	Tes = Tri;
-	computed = false;
+	computed=false;
 	max_id = -1;
 	TotalFiniteVoronoiVolume=0;
+	area=0;
 	TotalInternalVoronoiPorosity=0;
 	TotalInternalVoronoiVolume=0;
 	redirected = false;
@@ -351,6 +352,31 @@
 	return Segment ( f_it->first->info(), ( f_it->first->neighbor ( f_it->second ) )->info() );
 }
 
+double Tesselation::ComputeVFacetArea ( Finite_edges_iterator ed_it )
+{
+	Cell_circulator cell0 = Tri->incident_cells ( *ed_it );
+	Cell_circulator cell2 = cell0;
+
+	if ( Tri->is_infinite ( cell2 ) )
+	{
+		++cell2;
+		while ( Tri->is_infinite ( cell2 ) && cell2!=cell0 ) ++cell2;
+		if ( cell2==cell0 ) return 0;
+	}
+	cell0=cell2++;
+	Cell_circulator cell1=cell2++;
+	Real area = 0;
+	
+	while ( cell2!=cell0 )
+	{
+	  	area+= sqrt(std::abs (( Triangle ( cell0->info(), cell1->info(), cell2->info() ) ).squared_area())) ;
+		++cell1; 
+		++cell2;
+	}
+	
+	return area;
+}
+
 void Tesselation::AssignPartialVolume ( Finite_edges_iterator& ed_it )
 {
 	//Edge_iterator ed_it

=== modified file 'lib/triangulation/Tesselation.h'
--- lib/triangulation/Tesselation.h	2011-04-22 09:09:28 +0000
+++ lib/triangulation/Tesselation.h	2011-05-02 09:17:49 +0000
@@ -36,6 +36,7 @@
 	bool computed;
 public:
 	Real TotalFiniteVoronoiVolume;
+	Real area; 
 	Real TotalInternalVoronoiVolume;
 	Real TotalInternalVoronoiPorosity;
 	Vector_Vertex vertexHandles;//This is a redirection vector to get vertex pointers by spheres id
@@ -71,6 +72,8 @@
 	static Segment  Dual	(Finite_facets_iterator &facet);	//G�n�re le segment dual d'une facette finie
 	static Real	Volume	(Finite_cells_iterator cell);
 	inline void 	AssignPartialVolume	(Finite_edges_iterator& ed_it);
+// 	inline void 	ComputeVFacetArea	(Finite_edges_iterator& ed_it);
+	double ComputeVFacetArea (Finite_edges_iterator ed_it);
 	void		ResetVCellVolumes	(void);
 	void		ComputeVolumes		(void);//Compute volume each voronoi cell
 	void		ComputePorosity		(void);//Compute volume and porosity of each voronoi cell
@@ -95,7 +98,6 @@
 															//du graph de Voronoi
 	long New_liste_short_edges2	( Real** Coordonnes );
 	long New_liste_adjacent_edges ( Vertex_handle vertex0, Real** Coordonnes );
-	
 };
 
 

=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp	2011-04-27 15:15:04 +0000
+++ pkg/dem/FlowEngine.cpp	2011-05-02 09:17:49 +0000
@@ -79,6 +79,43 @@
 		if (!CachedForces) flow->ComputeFacetForces();
 		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;
+		  
+		}
+		}
 
 		///End Compute flow and forces
 		CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end();

=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp	2011-04-27 15:15:04 +0000
+++ pkg/dem/FlowEngine.hpp	2011-05-02 09:17:49 +0000
@@ -94,7 +94,7 @@
 					((bool,compute_K,false,,"Activates permeability measure within a granular sample"))
 					((bool,meanK_correction,true,,"Local permeabilities' correction through meanK threshold"))
 					((bool,meanK_opt,false,,"Local permeabilities' correction through an optimized threshold"))
-					((double,permeability_factor,0,,"permability multiplier"))
+					((double,permeability_factor,0.0,,"permability multiplier"))
 					((double,viscosity,1.0,,"viscosity of fluid"))
 					((Real,loadFactor,1.1,,"Load multiplicator for oedometer test"))
 					((double, K, 0,, "Permeability of the sample"))
@@ -132,6 +132,7 @@
 					((bool, FRONT_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
 					((bool, BACK_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
 					((bool, areaR2Permeability, 1,,"Use corrected formula for permeabilities calculation in flowboundingsphere (areaR2permeability variable)"))
+					((bool, viscousShear, false,,"Compute viscous shear terms as developped by Donia Marzougui"))
 					,,
 					timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
 					,