yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3953
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Mon 2016-10-24 13:17:04 +0200
message:
  Introduce Network::surfaceSolidThroatInPore() for partial solid surfaces per half-throat
  
  - Introduce Network::surfaceSolidThroatInPore() for partial solid
  surfaces per half-throat
  - consistent triangulation of clumps vs. clump members depending on
  Flow::convertClumps
  - consistent use of mask for selective triangulation of bodies
  - conditional summation of forces to skip boundaries when introduced as
  additional bodies (grep maxBodyId)
  - make consistent use of idOffset for those cases
  - introduce wall normals
modified:
  lib/triangulation/FlowBoundingSphere.hpp
  lib/triangulation/Network.hpp
  lib/triangulation/Network.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.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp	2014-10-29 16:49:20 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2016-10-24 11:17:04 +0000
@@ -30,7 +30,7 @@
 		//painfull, but we need that for templates inheritance...
 		using _N::T; using _N::xMin; using _N::xMax; using _N::yMin; using _N::yMax; using _N::zMin; using _N::zMax; using _N::Rmoy; using _N::sectionArea; using _N::Height; using _N::vTotal; using _N::currentTes; using _N::debugOut; using _N::nOfSpheres; using _N::xMinId; using _N::xMaxId; using _N::yMinId; using _N::yMaxId; using _N::zMinId; using _N::zMaxId; using _N::boundsIds; using _N::cornerMin; using _N::cornerMax;  using _N::VSolidTot; using _N::Vtotalissimo; using _N::vPoral; using _N::sSolidTot; using _N::vPoralPorosity; using _N::vTotalPorosity; using _N::boundaries; using _N::idOffset; using _N::vtkInfiniteVertices; using _N::vtkInfiniteCells; using _N::num_particles; using _N::boundingCells; using _N::facetVertices; using _N::facetNFictious;
 		//same for functions
-		using _N::defineFictiousCells; using _N::addBoundingPlanes; using _N::boundary; using _N::tesselation;
+		using _N::defineFictiousCells; using _N::addBoundingPlanes; using _N::boundary; using _N::tesselation; using _N::surfaceSolidThroatInPore;
 
 		virtual ~FlowBoundingSphere();
  		FlowBoundingSphere();

=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp	2014-06-06 15:27:15 +0000
+++ lib/triangulation/Network.hpp	2016-10-24 11:17:04 +0000
@@ -78,7 +78,8 @@
 		double volumeDoubleFictiousPore(VertexHandle SV1, VertexHandle SV2, VertexHandle SV3, Point PV1);
 		double volumeSingleFictiousPore(VertexHandle SV1, VertexHandle SV2, VertexHandle SV3, Point PV1);
 		double volumePoreVoronoiFraction ( CellHandle& cell, int& j, bool reuseFacetData=false);
-		double surfaceSolidPore( CellHandle cell, int j, bool slipBoundary, bool reuseFacetData=false);
+		double surfaceSolidThroat( CellHandle cell, int j, bool slipBoundary, bool reuseFacetData=false);
+		double surfaceSolidThroatInPore( CellHandle cell, int j, bool slipBoundary, bool reuseFacetData=false);// returns the solid area in the throat, keeping only that part of the throat in cell
 		double sphericalTriangleArea ( Sphere STA1, Sphere STA2, Sphere STA3, Point PTA1 );
 		
 		CVector surfaceDoubleFictiousFacet(VertexHandle fSV1, VertexHandle fSV2, VertexHandle SV3);

=== modified file 'lib/triangulation/Network.ipp'
--- lib/triangulation/Network.ipp	2015-03-26 09:17:44 +0000
+++ lib/triangulation/Network.ipp	2016-10-24 11:17:04 +0000
@@ -251,7 +251,7 @@
 }
 
 template<class Tesselation>
-double Network<Tesselation>::surfaceSolidPore(CellHandle cell, int j, bool slipBoundary, bool reuseFacetData)
+double Network<Tesselation>::surfaceSolidThroat(CellHandle cell, int j, bool slipBoundary, bool reuseFacetData)
 {
   if (!reuseFacetData)  facetNFictious=detectFacetFictiousVertices(cell,j);
   Point& p1 = cell->info();
@@ -358,6 +358,70 @@
 }
 
 template<class Tesselation>
+double Network<Tesselation>::surfaceSolidThroatInPore(CellHandle cell, int j, bool slipBoundary, bool reuseFacetData)
+{
+  if (!reuseFacetData)  facetNFictious=detectFacetFictiousVertices(cell,j);
+  Point& p1 = cell->info();
+  Point& p2 = cell->neighbor(j)->info();
+  double Ssolid1= 0, Ssolid2= 0, Ssolid3= 0;
+  Sphere v [3];
+  VertexHandle W [3];
+
+  for (int kk=0; kk<3; kk++) {
+	  W[kk] = cell->vertex(facetVertices[j][kk]);
+	  v[kk] = cell->vertex(facetVertices[j][kk])->point();}
+
+  switch (facetNFictious) {
+    case (0) : {
+		VertexHandle& SV1 = W[0];
+                VertexHandle& SV2 = W[1];
+                VertexHandle& SV3 = W[2];
+
+		Ssolid1 = fastSphericalTriangleArea(SV1->point(), SV2->point(), p1, SV3->point());
+                Ssolid2 = fastSphericalTriangleArea(SV2->point(),SV1->point(),p1, SV3->point());
+                Ssolid3 = fastSphericalTriangleArea(SV3->point(),SV2->point(),p1, SV1->point());
+    }; break;
+    case (1) : {
+		VertexHandle SV1 = cell->vertex(facetVertices[j][facetF1]);
+		VertexHandle SV2 = cell->vertex(facetVertices[j][facetRe1]);
+		VertexHandle SV3 = cell->vertex(facetVertices[j][facetRe2]);
+
+		Boundary &bi1 =  boundary(SV1->info().id());
+                Ssolid1 = 0;
+		if (bi1.flowCondition && ! slipBoundary) Ssolid1 = abs(0.5*CGAL::cross_product(p1-SV2->point(), SV2->point()-SV3->point())[bi1.coordinate]);
+                Ssolid2 = fastSphericalTriangleArea(SV2->point(),SV3->point(),p1, SV2->point()+bi1.normal);
+                Ssolid3 = fastSphericalTriangleArea(SV3->point(),SV2->point(),p1, SV3->point()+bi1.normal);
+    }; break;
+    case (2) : {
+		double A [3], B[3], C[3];
+		VertexHandle SV1 = cell->vertex(facetVertices[j][facetF1]);
+		VertexHandle SV2 = cell->vertex(facetVertices[j][facetF2]);
+		VertexHandle SV3 = cell->vertex(facetVertices[j][facetRe1]);
+		Boundary &bi1 =  boundary(SV1->info().id());
+                Boundary &bi2 =  boundary(SV2->info().id());
+                for (int m=0;m<3;m++) {A[m]=B[m]=C[m]= (SV3->point())[m];}
+                A[bi1.coordinate]=bi1.p[bi1.coordinate];
+                B[bi2.coordinate]=bi2.p[bi2.coordinate];
+                C[bi1.coordinate]=bi1.p[bi1.coordinate];
+                C[bi2.coordinate]=bi2.p[bi2.coordinate];
+                Point AA(A[0],A[1],A[2]);
+                Point BB(B[0],B[1],B[2]);
+                Point CC(C[0],C[1],C[2]);
+
+                Sphere A1(AA, 0);
+                Sphere B1(BB, 0);
+                Sphere C1(CC, 0);
+                //FIXME : we are computing triangle area twice here, because its computed in volume_double_fictious already -> optimize
+                Ssolid1 = 0.5*(fastSphericalTriangleArea(SV3->point(), AA, p1, p2)+ fastSphericalTriangleArea(SV3->point(), BB, p1, p2));
+                CVector p1p2v1Surface = 0.5*CGAL::cross_product(p1-p2,SV3->point()-p2);
+                if (bi1.flowCondition && ! slipBoundary) Ssolid2 = 0.5*abs(p1p2v1Surface[bi1.coordinate]);
+                if (bi2.flowCondition && ! slipBoundary) Ssolid3 = 0.5*abs(p1p2v1Surface[bi2.coordinate]); 
+    }; break;
+    }
+    return Ssolid1+Ssolid2+Ssolid3;
+}
+
+template<class Tesselation>
 CVector Network<Tesselation>::surfaceDoubleFictiousFacet(VertexHandle fSV1, VertexHandle fSV2, VertexHandle SV3)
 {
         //This function is correct only with axis-aligned boundaries

=== modified file 'pkg/pfv/FlowEngine.hpp.in'
--- pkg/pfv/FlowEngine.hpp.in	2016-02-10 18:57:40 +0000
+++ pkg/pfv/FlowEngine.hpp.in	2016-10-24 11:17:04 +0000
@@ -75,7 +75,6 @@
 		Vector3r normal [6];
 		bool currentTes;
 		bool metisForced;
-		int idOffset;
 		double epsVolCumulative;
 		int ReTrg;
 		int ellapsedIter;
@@ -165,6 +164,7 @@
 		Real volumeCellTripleFictious (Cellhandle cell);
 		template<class Cellhandle>
 		Real volumeCell (Cellhandle cell);
+		Real surfaceSolidThroatInPore(int cellId, int throatIndex) {return solver->surfaceSolidThroatInPore(solver->T[solver->currentTes].cellHandles[cellId], throatIndex, false, false);}
 		template<class Cellhandle>
 		Vector3r cellBarycenterFromHandle (Cellhandle cell) {return makeVector3r(solver->cellBarycenter(cell));}
 		Vector3r cellBarycenterFromId (unsigned long id) {
@@ -198,6 +198,11 @@
 			solver->T[solver->currentTes].cellHandles[id]->info().Pcondition=blockPressure;
 		}
 		
+		Vector3r getBoundaryNormal(int index) {
+		  if (index<0 || index>5) LOG_ERROR("index out of range (0-5)"); return normal[max(0,min(5,index))];}
+		void setBoundaryNormal(int index, Vector3r v) {
+		  if (index<0 || index>5) LOG_ERROR("index out of range (0-5)"); normal[max(0,min(5,index))]=v;}
+		
 		double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
 		double averagePressure(){return solver->averagePressure();}
 
@@ -233,6 +238,8 @@
 		YADE_CLASS_PYCLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TemplateFlowEngine_@TEMPLATE_FLOW_NAME@,@TEMPLATE_FLOW_NAME@,PartialEngine,"A generic engine from wich more specialized engines can inherit. It is defined for the sole purpose of inserting the right data classes CellInfo and VertexInfo in the triangulation, and it should not be used directly. Instead, look for specialized engines, e.g. :yref:`FlowEngine`, :yref:`PeriodicFlowEngine`, or :yref:`DFNFlowEngine`.",
 		((bool,isActivated,true,,"Activates Flow Engine"))
 		((bool,first,true,,"Controls the initialization/update phases"))
+		((int,idOffset,0,,"If the bounding walls of the fluid mesh are not walls of the scene (i.e. are not elements of O.bodies), the offset should be set equal to the size of O.bodies. If the bounding walls are bodies of the scene but are not numbered as 0-5 then offset should be the number of bodies comming before the walls. Set offset<0 to get it set equal to O.bodies.size(), it will also update :yref:`FlowEngine::wallIds`."))
+		((bool,convertClumps,true,,"If true the clumps will be temptatively converted into equivalent spheres in the triangulation, and clump members are skipped. Else clumps are ignored and spherical clump members are triangulated as independent bodies."))
 		((bool,doInterpolate,false,,"Force the interpolation of cell's info while remeshing. By default, interpolation would be done only for compressible fluids. It can be forced with this flag."))
 		((double, fluidBulkModulus, 0.,,"Bulk modulus of fluid (inverse of compressibility) K=-dP*V/dV [Pa]. Flow is compressible if fluidBulkModulus > 0, else incompressible."))
 		((Real, dt, 0,,"timestep [s]"))
@@ -258,7 +265,7 @@
 		((double,permeabilityFactor,1.0,,"permability multiplier"))
 		((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=G-Seidel, 1=Taucs, 2-Pardiso, 3-CHOLMOD"))
+		((int, useSolver, 0,, "Solver to use. 0:Gauss-Seidel, >0: Cholesky factorization (sparse CHOLMOD)"))
 		((int, xmin,0,(Attr::readonly),"Index of the boundary $x_{min}$. This index is not equal the the id of the corresponding body in general, it may be used to access the corresponding attributes (e.g. flow.bndCondValue[flow.xmin], flow.wallId[flow.xmin],...)."))
 		((int, xmax,1,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
 		((int, ymin,2,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
@@ -320,6 +327,7 @@
 		.def("getConstrictions",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getConstrictions,(boost::python::arg("all")=true),"Get the list of constriction radii (inscribed circle) for all finite facets (if all==True) or all facets not incident to a virtual bounding sphere (if all==False).  When all facets are returned, negative radii denote facet incident to one or more fictious spheres.")
 		.def("getConstrictionsFull",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getConstrictionsFull,(boost::python::arg("all")=true),"Get the list of constrictions (inscribed circle) for all finite facets (if all==True), or all facets not incident to a fictious bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres. The constrictions are returned in the format {{cell1,cell2}{rad,nx,ny,nz}}")
 		.def("edgeSize",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::edgeSize,"Return the number of interactions.")
+		.def("surfaceSolidThroatInPore",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::surfaceSolidThroatInPore,(boost::python::arg("cellId"),boost::python::arg("throatIndex")),"returns solid area in the throat (index 0-3), keeping only that part of the throat in cell.")
 		.def("OSI",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::OSI,"Return the number of interactions only between spheres.")
 		.def("saveVtk",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
 		.def("avFlVelOnSph",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::avFlVelOnSph,(boost::python::arg("idSph")),"compute a sphere-centered average fluid velocity")
@@ -360,6 +368,7 @@
 		.def("compTessVolumes",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::compTessVolumes,"Like TesselationWrapper::computeVolumes()")
 		.def("volume",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getVolume,(boost::python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
 		.def("averageVelocity",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averageVelocity,"measure the mean velocity in the period")
+		.def("setBoundaryNormal",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::setBoundaryNormal,"define the unit outward-pointing normal of a boundary (0<=index<=5).")
 		)
 };
 

=== modified file 'pkg/pfv/FlowEngine.ipp.in'
--- pkg/pfv/FlowEngine.ipp.in	2016-02-11 18:44:23 +0000
+++ pkg/pfv/FlowEngine.ipp.in	2016-10-24 11:17:04 +0000
@@ -19,7 +19,6 @@
 TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::~TemplateFlowEngine_@TEMPLATE_FLOW_NAME@() {}
 
 // YADE_PLUGIN((TFlowEng));
-
 template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
 unsigned int TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposePressure(Vector3r pos, Real p)
 {
@@ -76,21 +75,22 @@
 	Vector3r force;
 	Vector3r torque;
         FiniteVerticesIterator verticesEnd = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
+	int maxBodyId = 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);
 		torque = Vector3r(0,0,0);
                 if (shearLubrication || viscousShear){
-			force = force + solver->shearLubricationForces[vIt->info().id()];
-			torque = torque + solver->shearLubricationTorques[vIt->info().id()];
+			force = force + solver->shearLubricationForces[vId];
+			torque = torque + solver->shearLubricationTorques[vId];
 			if (pumpTorque)
-				torque = torque + solver->pumpLubricationTorques[vIt->info().id()];
+				torque = torque + solver->pumpLubricationTorques[vId];
 		}
-		if (twistTorque)
-			torque = torque + solver->twistLubricationTorques[vIt->info().id()];
-		if (normalLubrication)
-			force = force + solver-> normalLubricationForce[vIt->info().id()];
-		scene->forces.addForce ( vIt->info().id(), force);
-		scene->forces.addTorque ( vIt->info().id(), torque);
+		if (twistTorque) torque = torque + solver->twistLubricationTorques[vId];
+		if (normalLubrication) force = force + solver-> normalLubricationForce[vId];
+		if (vIt->info().id() < maxBodyId) {
+			scene->forces.addForce ( vId, force);
+			scene->forces.addTorque ( vId, torque);}
         }
         ///End compute flow and forces
         timingDeltas->checkpoint ( "Applying Forces" );
@@ -294,8 +294,7 @@
 	shared_ptr<Sphere> sph ( new Sphere );
 	const int Sph_Index = sph->getClassIndexStatic();
 	FOREACH ( const shared_ptr<Body>& b, *scene->bodies ) {
-// 		if (!b || ignoredBody==b->getId() || b->isClumpMember()) continue;
-		if (!b || (mask>0 and !b->maskCompatible(mask)) || b->isClumpMember()) continue;
+		if (!b || (mask>0 and !b->maskCompatible(mask)) || (convertClumps && b->isClumpMember()) || (b->isClump() && !convertClumps)) continue;
 		posData& dat = buffer[b->getId()];
 		dat.id=b->getId();
 		dat.pos=b->state->pos;
@@ -366,7 +365,7 @@
                 else {
 			for ( int h=0;h<3;h++ ) center[h] = buffer[*flow.boundsIds[i]].pos[h];
 // 			cerr << "id="<<*flow.boundsIds[i] <<" center="<<center[0]<<","<<center[1]<<","<<center[2]<<endl;
-                        flow.addBoundingPlane ( center, wallThickness, Normal,*flow.boundsIds[i] );
+                        flow.addBoundingPlane (center, wallThickness, Normal,*flow.boundsIds[i] );
                 }
         }
 }