yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3971
committer: bchareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Wed 2016-11-16 10:53:52 +0100
message:
  Implement Darcy scale permeability in FlowEngine, introduce alphaShapes (commented), some code cleaning and smallfixes
modified:
  lib/triangulation/FlowBoundingSphere.ipp
  lib/triangulation/FlowBoundingSphereLinSolv.ipp
  lib/triangulation/PeriodicFlow.hpp
  lib/triangulation/RegularTriangulation.h
  lib/triangulation/Tesselation.h
  lib/triangulation/Tesselation.ipp
  pkg/pfv/FlowEngine.hpp.in
  pkg/pfv/FlowEngine.ipp.in


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp	2016-11-02 15:14:59 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp	2016-11-16 09:53:52 +0000
@@ -578,7 +578,9 @@
 	meanK /= pass;
 	meanRadius /= pass;
 	meanDistance /= pass;
-	Real globalK=kFactor*meanDistance*vPoral/(sSolidTot*8.*viscosity);//An approximate value of macroscopic permeability, for clamping local values below
+	Real globalK;
+	if (kFactor>0) globalK=kFactor*meanDistance*vPoral/(sSolidTot*8.*viscosity);//An approximate value of macroscopic permeability, for clamping local values below
+	else globalK=meanK;
 	if (debugOut) {
 		cout << "PassCompK = " << pass << endl;
 		cout << "meanK = " << meanK << endl;

=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.ipp'
--- lib/triangulation/FlowBoundingSphereLinSolv.ipp	2015-09-17 14:08:41 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.ipp	2016-11-16 09:53:52 +0000
@@ -160,7 +160,6 @@
 			if (!cell->info().Pcondition && !cell->info().blocked) ++ncols;}
 //		//Segfault on 14.10, and useless overall since SuiteSparse has preconditionners (including metis)
 // 		spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort<RTriangulation>());
-
 		T_cells.clear();
 		T_index=0;
 		isLinearSystemSet=false;
@@ -273,12 +272,8 @@
 		#ifdef EIGENSPARSE_LIB
 		} else if (useSolver==3){
 			tripletList.clear(); tripletList.resize(T_nnz);
-			for(int k=0;k<T_nnz;k++) {
-// 				if (is[k]<js[k]) cerr<<"not the good relation"<<endl;
-// 				else cerr<< "comp " <<is[k]-1<<" "<<js[k]-1<<" "<<vs[k]<<endl;
-				tripletList[k]=ETriplet(is[k]-1,js[k]-1,vs[k]);
-			}
-			A.resize(ncols,ncols);
+			for(int k=0;k<T_nnz;k++) tripletList[k]=ETriplet(is[k]-1,js[k]-1,vs[k]);
+ 			A.resize(ncols,ncols);
 			A.setFromTriplets(tripletList.begin(), tripletList.end());
 		#endif
 		}
@@ -526,7 +521,7 @@
 		eSolver.compute(A);
 		//Check result
 		if (eSolver.cholmod().status>0) {
-			cerr << "something went wrong in Cholesky factorization, use LDLt as fallback this time" << endl;
+			cerr << "something went wrong in Cholesky factorization, use LDLt as fallback this time" << eSolver.cholmod().status << endl;
 			eSolver.setMode(Eigen::CholmodLDLt);
 			eSolver.compute(A);
 		}

=== modified file 'lib/triangulation/PeriodicFlow.hpp'
--- lib/triangulation/PeriodicFlow.hpp	2016-11-02 15:14:59 +0000
+++ lib/triangulation/PeriodicFlow.hpp	2016-11-16 09:53:52 +0000
@@ -14,6 +14,11 @@
 // #endif
 
 
+/* TODO:
+- remove computePermeability(), mostly a duplicate of the non-periodic version
+
+*/
+
 namespace CGT{
 
 // 	typedef CGT::FlowBoundingSphere<PeriFlowTesselation> PeriodicFlowBoundingSphere;

=== modified file 'lib/triangulation/RegularTriangulation.h'
--- lib/triangulation/RegularTriangulation.h	2014-10-15 06:44:01 +0000
+++ lib/triangulation/RegularTriangulation.h	2016-11-16 09:53:52 +0000
@@ -12,6 +12,12 @@
 #include <CGAL/Cartesian.h>
 #include <CGAL/Regular_triangulation_3.h>
 #include <CGAL/Regular_triangulation_euclidean_traits_3.h>
+#define ALPHASHAPES
+#ifdef ALPHASHAPES
+	#include <CGAL/Alpha_shape_vertex_base_3.h>
+	#include <CGAL/Alpha_shape_cell_base_3.h>
+	#include <CGAL/Alpha_shape_3.h>
+#endif
 #include <CGAL/Triangulation_vertex_base_with_info_3.h>
 #include <CGAL/Triangulation_cell_base_with_info_3.h>
 #include <CGAL/Delaunay_triangulation_3.h>
@@ -92,11 +98,19 @@
 typedef cell_info								Cell_Info;
 typedef CGAL::Triangulation_vertex_base_with_info_3<Vertex_Info, Traits>	Vb_info;
 typedef CGAL::Triangulation_cell_base_with_info_3<Cell_Info, Traits>		Cb_info;
+#ifdef ALPHASHAPES
+typedef CGAL::Alpha_shape_vertex_base_3<Traits,Vb_info> Vb;
+typedef CGAL::Alpha_shape_cell_base_3<Traits,Cb_info>   Fb;
+typedef CGAL::Triangulation_data_structure_3<Vb, Fb>				Tds;
+#else
 typedef CGAL::Triangulation_data_structure_3<Vb_info, Cb_info>			Tds;
+#endif
 
 typedef CGAL::Triangulation_3<K>						Triangulation;
 typedef CGAL::Regular_triangulation_3<Traits, Tds>				RTriangulation;
-
+#ifdef ALPHASHAPES
+typedef CGAL::Alpha_shape_3<RTriangulation>  					AlphaShape;
+#endif
 typedef typename RTriangulation::Vertex_iterator                    		VertexIterator;
 typedef typename RTriangulation::Vertex_handle                      		VertexHandle;
 typedef typename RTriangulation::Finite_vertices_iterator                    	FiniteVerticesIterator;

=== modified file 'lib/triangulation/Tesselation.h'
--- lib/triangulation/Tesselation.h	2014-03-21 18:47:45 +0000
+++ lib/triangulation/Tesselation.h	2016-11-16 09:53:52 +0000
@@ -13,6 +13,7 @@
 //Since template inheritance does not automatically give access to the members of the base class, this macro can be used to declare all members at once. 
 #define DECLARE_TESSELATION_TYPES(baseType)\
 		typedef typename baseType::RTriangulation		 	RTriangulation;\
+		typedef typename baseType::AlphaShape		 		AlphaShape;\
 		typedef typename baseType::VertexInfo				VertexInfo;\
 		typedef typename baseType::CellInfo				CellInfo;\
 		typedef typename baseType::VertexIterator			VertexIterator;\
@@ -43,6 +44,7 @@
 {
 public:
 	typedef typename TT::RTriangulation							RTriangulation;
+	typedef typename TT::AlphaShape						 		AlphaShape;
 	typedef typename TT::Vertex_Info							VertexInfo;
 	typedef typename TT::Cell_Info								CellInfo;
 	typedef typename RTriangulation::Vertex_iterator		 			VertexIterator;

=== modified file 'lib/triangulation/Tesselation.ipp'
--- lib/triangulation/Tesselation.ipp	2015-05-19 19:03:44 +0000
+++ lib/triangulation/Tesselation.ipp	2016-11-16 09:53:52 +0000
@@ -166,6 +166,16 @@
 		cell->info().setPoint(Point(x,y,z));
 	}
 	computed = true;
+	
+//	AlphaShape as (*Tri);
+//	as.set_alpha(as.find_alpha_solid());
+//	cerr << "Alpha shape computed" <<endl;
+//	std::list<CellHandle> cells,cells2;
+//	as.get_alpha_shape_cells(std::back_inserter(cells),
+  //             AlphaShape::EXTERIOR);
+//	as.get_alpha_shape_cells(std::back_inserter(cells2),
+  //             AlphaShape::INTERIOR);
+//	std::cerr<< "num exterior cells "<< cells.size() <<" vs. "<<cells2.size() <<std::endl;
 }
 
 template<class TT>

=== modified file 'pkg/pfv/FlowEngine.hpp.in'
--- pkg/pfv/FlowEngine.hpp.in	2016-10-24 11:17:04 +0000
+++ pkg/pfv/FlowEngine.hpp.in	2016-11-16 09:53:52 +0000
@@ -262,7 +262,7 @@
 		((bool,clampKValues,true,,"If true, clamp local permeabilities in [minKdivKmean,maxKdivKmean]*globalK. This clamping can avoid singular values in the permeability matrix and may reduce numerical errors in the solve phase. It will also hide junk values if they exist, or bias all values in very heterogeneous problems. So, use this with care."))
 		((Real,minKdivKmean,0.0001,,"define the min K value (see :yref:`FlowEngine::clampKValues`)"))
 		((Real,maxKdivKmean,100,,"define the max K value (see :yref:`FlowEngine::clampKValues`)"))
-		((double,permeabilityFactor,1.0,,"permability multiplier"))
+		((double,permeabilityFactor,1.0,,"Permability multiplier ($m$): $m=1$ (default) attempts to predicty the actual hydraulic conductivity using a Poiseuille equation; $m>0$ multiplies the default values by $m$; $m<0$ defines the conductivity independently of particle size and viscosity as if the material was a homogeneous continuum of conductivity $-m$"))
 		((double,viscosity,1.0,,"viscosity of the fluid"))
 		((double,stiffness, 10000,,"equivalent contact stiffness used in the lubrication model"))
 		((int, useSolver, 0,, "Solver to use. 0:Gauss-Seidel, >0: Cholesky factorization (sparse CHOLMOD)"))
@@ -275,14 +275,11 @@
 
 		((vector<bool>, bndCondIsPressure, vector<bool>(6,false),,"defines the type of boundary condition for each side. True if pressure is imposed, False for no-flux. Indexes can be retrieved with :yref:`FlowEngine::xmin` and friends."))
 		((vector<double>, bndCondValue, vector<double>(6,0),,"Imposed value of a boundary condition. Only applies if the boundary condition is imposed pressure, else the imposed flux is always zero presently (may be generalized to non-zero imposed fluxes in the future)."))
-		//FIXME: to be implemented:
 		((vector<Vector3r>, boundaryVelocity, vector<Vector3r>(6,Vector3r::Zero()),, "velocity on top boundary, only change it using :yref:`FlowEngine::setBoundaryVel`"))
 		((int, ignoredBody,-1,,"DEPRECATED, USE MASK - Id of a sphere to exclude from the triangulation.)"))
 		((int,mask,0,,"If mask defined, only bodies with corresponding groupMask will be affected by this engine. If 0, all bodies will be affected."))
 		((vector<int>, wallIds,vector<int>(6),,"body ids of the boundaries (default values are ok only if aabbWalls are appended before spheres, i.e. numbered 0,...,5)"))
 		((vector<bool>, boundaryUseMaxMin, vector<bool>(6,true),,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
-// 					((bool, display_force, false,,"display the lubrication force applied on particles"))
-// 					((bool, create_file, false,,"create file of velocities"))
 		((bool, viscousShear, false,,"compute viscous shear terms as developped by Donia Marzougui (FIXME: ref.)"))
 		((bool, shearLubrication, false,,"compute shear lubrication force as developped by Brule (FIXME: ref.) "))
 		((bool, pumpTorque, false,,"Compute pump torque applied on particles "))

=== modified file 'pkg/pfv/FlowEngine.ipp.in'
--- pkg/pfv/FlowEngine.ipp.in	2016-10-24 11:17:04 +0000
+++ pkg/pfv/FlowEngine.ipp.in	2016-11-16 09:53:52 +0000
@@ -75,7 +75,7 @@
 	Vector3r force;
 	Vector3r torque;
         FiniteVerticesIterator verticesEnd = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
-	int maxBodyId = scene->bodies->size();
+	int bodiesLength = scene->bodies->size();
         for ( FiniteVerticesIterator vIt = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); vIt !=  verticesEnd; vIt++ ) {
 		int vId = vIt->info().id();
 		force = pressureForce ? Vector3r ( vIt->info().forces[0],vIt->info().forces[1],vIt->info().forces[2] ): Vector3r(0,0,0);
@@ -88,7 +88,7 @@
 		}
 		if (twistTorque) torque = torque + solver->twistLubricationTorques[vId];
 		if (normalLubrication) force = force + solver-> normalLubricationForce[vId];
-		if (vIt->info().id() < maxBodyId) {
+		if (vIt->info().id() < bodiesLength) {
 			scene->forces.addForce ( vId, force);
 			scene->forces.addTorque ( vId, torque);}
         }
@@ -329,8 +329,11 @@
                         flow.zMax = max ( flow.zMax, z+rad );
                 }
         }
-	//FIXME idOffset must be set correctly, not the case here (always 0), then we need walls first or it will fail
-        idOffset = flow.tesselation().maxId+1;
+	//set idOffset if<0, in that case overwrite wallIds and assign values out of range scene->bodies
+        if (idOffset<0) {
+	  idOffset = scene->bodies->size();
+	  for (int i=0; i<6; ++i){wallIds[i]=i+idOffset;}
+	}
         flow.idOffset = idOffset;
         flow.sectionArea = ( flow.xMax - flow.xMin ) * ( flow.zMax-flow.zMin );
         flow.vTotal = ( flow.xMax-flow.xMin ) * ( flow.yMax-flow.yMin ) * ( flow.zMax-flow.zMin );