------------------------------------------------------------ revno: 2800 committer: Emanuele Catalano <catal...@hmg.inpg.fr branch nick: yade timestamp: Fri 2011-04-01 11:59:17 +0200 message: - Correction on facet fluid/solid surfaces calculation modified: lib/triangulation/FlowBoundingSphere.cpp lib/triangulation/FlowBoundingSphere.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-03-23 09:38:16 +0000 +++ lib/triangulation/FlowBoundingSphere.cpp 2011-04-01 09:59:17 +0000 @@ -82,6 +82,7 @@ TOLERANCE = 1e-06; RELAX = 1.9; ks=0; + V_darcy_Donia=0; distance_correction = true; meanK_LIMIT = true; meanK_STAT = false; K_opt_factor=0; @@ -328,8 +329,8 @@ cell->info().av_vel() = cell->info().av_vel() + (facet_flow_rate - volume_facet_translation) * ( pos_av_facet-CGAL::ORIGIN ); }} if (cell->info().volume()){ tVel+=cell->info().av_vel()[1]; tVol+=cell->info().volume();} - cell->info().av_vel() = cell->info().av_vel() /cell->info().volume(); - } + cell->info().av_vel() = cell->info().av_vel() /cell->info().volume(); + } } bool FlowBoundingSphere::isOnSolid (double X, double Y, double Z) @@ -685,6 +686,9 @@ Cell_handle neighbour_cell; double k=0, distance = 0, radius = 0, viscosity = VISCOSITY; + double m3=0, m1=0, m2=0, d=0, h=0; + int surfneg=0; + Real S0=0; int NEG=0, POS=0, pass=0; bool ref = Tri.finite_cells_begin()->info().isvisited; @@ -694,7 +698,7 @@ // std::ofstream kFile ( "LocalPermeabilities" ,std::ios::app ); Real meanK=0, STDEV=0, meanRadius=0, meanDistance=0; Real infiniteK=1e10; -int count_k_neg=0; + for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) { Point& p1 = cell->info(); for (int j=0; j<4; j++) { @@ -709,6 +713,7 @@ Sphere& v0 = W[0]->point(); Sphere& v1 = W[1]->point(); Sphere& v2 = W[2]->point(); + Sphere& t0=v0, t1=v1, t2=v2; #ifdef USE_FAST_MATH //FIXME : code not compiling,, do the same as in "else" assert((W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point())>=0 && (W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point())<=1); @@ -734,16 +739,46 @@ if (distance!=0) { if (minPermLength>0 && distance_correction) distance=max(minPermLength,distance); + Real fluidArea=0; + const Vecteur& Surfk = cell->info().facetSurfaces[j]; + Real area = sqrt(Surfk.squared_length()); + const Vecteur& crossSections = cell->info().facetSphereCrossSections[j]; if (areaR2Permeability){ - const Vecteur& Surfk = cell->info().facetSurfaces[j]; - Real area = sqrt(Surfk.squared_length()); - const Vecteur& crossSections = cell->info().facetSphereCrossSections[j]; - Real fluidArea=area-crossSections[0]-crossSections[1]-crossSections[2]; - if (fluidArea<0) fluidArea = -fluidArea; - k=(fluidArea * pow(radius,2)) / (8*viscosity*distance); - - } - else k = (M_PI * pow(radius,4)) / (8*viscosity*distance); + //if (cell->info().fictious()==0 && neighbour_cell->info().fictious()==0){ + + m1=sqrt((cross_product((v0-v1),v2-v1)).squared_length()/(v2-v1).squared_length()); + m2=sqrt((cross_product((v0-v1),v2-v0)).squared_length()/(v2-v0).squared_length()); + m3=sqrt((cross_product((v0-v2),v1-v0)).squared_length()/(v1-v0).squared_length()); + + if (m1<sqrt(v0.weight())) + { + d=2*sqrt((v0.weight()-m1*m1)); + h=sqrt(v0.weight())-m1; + S0=0.25*M_PI*(d*d+4*h*h);} + + if (m2<sqrt(v1.weight())) + { + d=2*sqrt((v1.weight()-m2*m2)); + h=sqrt(v1.weight())-m2; + S0=0.25*M_PI*(d*d+4*h*h);} +// + if (m3<sqrt(v2.weight())) + { + d=2*sqrt((v2.weight()-m3*m3)); + h=sqrt(v2.weight())-m3; + S0=0.25*M_PI*(d*d+4*h*h);} + +// if (S0>0) cout<<"S0= "<<S0<<endl; + + // } + fluidArea=area-crossSections[0]-crossSections[1]-crossSections[2]+S0; + k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);} + + else k = (M_PI * pow(radius,4)) / (8*viscosity*distance); + + if (k<0) {surfneg+=1; + cout<<"__ k<0 __"<<k<<" "<<area<<" "<<crossSections[0]<<" "<<crossSections[1]<<" "<<crossSections[2] <<endl;} + } else k = infiniteK;//Will be corrected in the next loop (cell->info().k_norm())[j]= k*k_factor; @@ -760,6 +795,7 @@ } cell->info().isvisited = !ref; } + cout<<"surfneg est "<<surfneg<<endl; meanK /= pass; meanRadius /= pass; meanDistance /= pass; @@ -844,8 +880,8 @@ if (!RAVERAGE) cout << "------Hydraulic Radius is used for permeability computation------" << endl << endl; else cout << "------Average Radius is used for permeability computation------" << endl << endl; - cout << "-----Computed_Permeability-----" << endl; - cout << "Negative Permeabilities = " << count_k_neg << endl; } + cout << "-----Computed_Permeability-----" << endl;} +// cout << "Negative Permeabilities = " << count_k_neg << endl; } vector<double> FlowBoundingSphere::getConstrictions() @@ -1082,7 +1118,7 @@ #ifdef GS_OPEN_MP } while (j<1500); #else - } while ((dp_max/p_max) > tolerance && j<4000 /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/); + } while ((dp_max/p_max) > tolerance /*&& j<4000*/ /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/); #endif } @@ -1161,6 +1197,7 @@ double viscosity = VISCOSITY; double gravity = 9.80665; double Vdarcy = Q1/Section; + V_darcy_Donia=Vdarcy; // double GradP = abs(P_Inf-P_Sup) /DeltaY; double DeltaP = abs(P_Inf-P_Sup); // double GradH = GradP/ (density*gravity); === modified file 'lib/triangulation/FlowBoundingSphere.hpp' --- lib/triangulation/FlowBoundingSphere.hpp 2011-03-23 09:38:16 +0000 +++ lib/triangulation/FlowBoundingSphere.hpp 2011-04-01 09:59:17 +0000 @@ -45,7 +45,7 @@ bool OUTPUT_BOUDARIES_RADII; bool noCache;//flag for checking if cached values cell->unitForceVectors have been defined vector<pair<Point,Real> > imposedP; - + double V_darcy_Donia; void initNewTri () {noCache=true; /*isLinearSystemSet=false; areCellsOrdered=false;*/}//set flags after retriangulation bool computeAllCells;//exececute computeHydraulicRadius for all facets and all spheres (double cpu time but needed for now in order to define crossSections correctly)
_______________________________________________ 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