yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3973
committer: bchareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Thu 2016-11-17 16:32:26 +0100
message:
  fix negative volumes issue in triangulation cells, consistentlyy remove the includeBoundary option (a workaround) in averageVelocity
modified:
  lib/triangulation/FlowBoundingSphere.hpp
  lib/triangulation/FlowBoundingSphere.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	2016-11-16 17:38:08 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp	2016-11-17 15:32:26 +0000
@@ -156,7 +156,7 @@
 		double getPorePressure (double X, double Y, double Z);
 		void measurePressureProfile(double WallUpy, double WallDowny);
 		double averageSlicePressure(double Y);
-		double averagePressure(bool includeBoundaries);
+		double averagePressure();
 		int getCell (double X,double Y,double Z);
 		double boundaryFlux(unsigned int boundaryId);
 		void setBlocked(CellHandle& cell);

=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp	2016-11-16 17:38:08 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp	2016-11-17 15:32:26 +0000
@@ -164,8 +164,8 @@
 	{
 	  if (cell->info().fictious()==0){
 	    for (int i=0;i<4;i++){
-	      velocityVolumes[cell->vertex(i)->info().id()] =  velocityVolumes[cell->vertex(i)->info().id()] + cell->info().averageVelocity()*cell->info().volume();
-	      volumes[cell->vertex(i)->info().id()] = volumes[cell->vertex(i)->info().id()] + cell->info().volume();}
+	      velocityVolumes[cell->vertex(i)->info().id()] =  velocityVolumes[cell->vertex(i)->info().id()] + cell->info().averageVelocity()*std::abs(cell->info().volume());
+	      volumes[cell->vertex(i)->info().id()] = volumes[cell->vertex(i)->info().id()] + std::abs(cell->info().volume());}
 	  }}	    
 	
 	std::ofstream fluid_vel ("Velocity", std::ios::out);
@@ -204,8 +204,8 @@
 	  if (cell->info().fictious()==0){
 	    for (unsigned int i=0;i<4;i++){
 	      if (cell->vertex(i)->info().id()==Id_sph){
-		velocityVolumes = velocityVolumes + cell->info().averageVelocity()*cell->info().volume();
-		volumes = volumes + cell->info().volume();}}}}
+		velocityVolumes = velocityVolumes + cell->info().averageVelocity()*std::abs(cell->info().volume());
+		volumes = volumes + std::abs(cell->info().volume());}}}}
 		
 	for (int i=0;i<3;i++) result[i] += velocityVolumes[i]/volumes;
 	return result;
@@ -272,17 +272,16 @@
   return P_ave;
 }
 template <class Tesselation> 
-double FlowBoundingSphere<Tesselation>::averagePressure(bool includeBoundaries)
+double FlowBoundingSphere<Tesselation>::averagePressure()
 {
   RTriangulation& Tri = T[currentTes].Triangulation();
   double P = 0.f, Ppond=0.f, Vpond=0.f;
   int n = 0;
   for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); cell++) {
-	if (includeBoundaries || cell->info().isReal()){
-		P+=cell->info().p();
-		n++;
-		Ppond+=cell->info().p()*cell->info().volume();
-		Vpond+=cell->info().volume();}}
+	P+=cell->info().p();
+	n++;
+	Ppond+=cell->info().p()*std::abs(cell->info().volume());
+	Vpond+=std::abs(cell->info().volume());}
   P/=n;
   Ppond/=Vpond;
   return Ppond;

=== modified file 'pkg/pfv/FlowEngine.hpp.in'
--- pkg/pfv/FlowEngine.hpp.in	2016-11-16 17:38:08 +0000
+++ pkg/pfv/FlowEngine.hpp.in	2016-11-17 15:32:26 +0000
@@ -204,7 +204,7 @@
 		  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(bool includeBoundaries){return solver->averagePressure(includeBoundaries);}
+		double averagePressure(){return solver->averagePressure();}
 
 		void emulateAction(){
 			scene = Omega::instance().getScene().get();
@@ -342,7 +342,7 @@
 		.def("pressureProfile",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::pressureProfile,(boost::python::arg("wallUpY"),boost::python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
 		.def("getPorePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getPorePressure,(boost::python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
 		.def("averageSlicePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averageSlicePressure,(boost::python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
-		.def("averagePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averagePressure,(boost::python::arg("includeBoundaries")=true),"Measure averaged pore pressure in the entire volume, the cells adjacent to the boundaries are ignored if includeBoundaries=False")
+		.def("averagePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averagePressure,"Measure averaged pore pressure in the entire volume, the cells adjacent to the boundaries are ignored if includeBoundaries=False")
 		.def("updateBCs",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::updateBCs,"tells the engine to update it's boundary conditions before running (especially useful when changing boundary pressure - should not be needed for point-wise imposed pressure)")
 		.def("emulateAction",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
 		.def("getCell",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getCell,(boost::python::arg("pos")),"get id of the cell containing (X,Y,Z).")

=== modified file 'pkg/pfv/FlowEngine.ipp.in'
--- pkg/pfv/FlowEngine.ipp.in	2016-11-16 09:53:52 +0000
+++ pkg/pfv/FlowEngine.ipp.in	2016-11-17 15:32:26 +0000
@@ -446,7 +446,7 @@
                                 CGT::CVector branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
                                 pos_av_facet = ( CGT::Point ) cell->info() + ( branch*Surfk ) *Surfk;
                                 volume_facet_translation += Vel_av*cell->info().facetSurfaces[i];
-                                cell->info().averageVelocity() = cell->info().averageVelocity() - volume_facet_translation/cell->info().volume() * ( pos_av_facet-CGAL::ORIGIN );
+                                cell->info().averageVelocity() = cell->info().averageVelocity() - volume_facet_translation/std::abs(cell->info().volume()) * ( pos_av_facet-CGAL::ORIGIN );
                         }
                 }
         }
@@ -585,7 +585,8 @@
 	const Vector3r& p1 = positionBufferCurrent[cell->vertex ( 1 )->info().id()].pos;
 	const Vector3r& p2 = positionBufferCurrent[cell->vertex ( 2 )->info().id()].pos;
 	const Vector3r& p3 = positionBufferCurrent[cell->vertex ( 3 )->info().id()].pos;
-	Real volume = inv6 * ((p0-p1).cross(p0-p2)).dot(p0-p3);
+	Real volume = -inv6 * ((p0-p1).cross(p0-p2)).dot(p0-p3);
+	if (/*debug && */volume<0) cerr<<"negative volume for an ordinary pore (temp warning, should still be safe)"<<endl; 
         if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
         return volume;
 }