------------------------------------------------------------ revno: 3973 committer: bchareyre <bruno.chare...@grenoble-inp.fr> 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; }
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp