yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11446
[Branch ~yade-pkg/yade/git-trunk] Rev 3442: -clean code
------------------------------------------------------------
revno: 3442
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Fri 2014-07-04 16:08:45 +0200
message:
-clean code
modified:
pkg/pfv/UnsaturatedEngine.cpp
--
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 'pkg/pfv/UnsaturatedEngine.cpp'
--- pkg/pfv/UnsaturatedEngine.cpp 2014-07-04 14:02:36 +0000
+++ pkg/pfv/UnsaturatedEngine.cpp 2014-07-04 14:08:45 +0000
@@ -750,173 +750,6 @@
vtkfile.end_data();
}
-/*
-double UnsaturatedEngine::computeEffPoreRadiusFine(CellHandle cell, int j)
-{
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- if (tri.is_infinite(cell->neighbor(j))) return 0;
-
- Vector3r posA = makeVector3r(cell->vertex(facetVertices[j][0])->point().point());
- Vector3r posB = makeVector3r(cell->vertex(facetVertices[j][1])->point().point());
- Vector3r posC = makeVector3r(cell->vertex(facetVertices[j][2])->point().point());
- double rA = sqrt(cell->vertex(facetVertices[j][0])->point().weight());
- double rB = sqrt(cell->vertex(facetVertices[j][1])->point().weight());
- double rC = sqrt(cell->vertex(facetVertices[j][2])->point().weight());
- double AB = (posA-posB).norm();
- double AC = (posA-posC).norm();
- double BC = (posB-posC).norm();
- double rAB = 0.5*(AB-rA-rB); if (rAB<0) { rAB=0; }
- double rBC = 0.5*(BC-rB-rC); if (rBC<0) { rBC=0; }
- double rAC = 0.5*(AC-rA-rC); if (rAC<0) { rAC=0; }
-
- double rmin = max(rAB,max(rBC,rAC)); if (rmin==0) { rmin= 1.0e-10; }
- double rmax = abs(solver->computeEffectiveRadius(cell, j));//rmin>rmax ?
-// if(rmin>rmax) { cerr<<"WARNING! rmin>rmax. rmin="<<rmin<<" ,rmax="<<rmax<<endl; }
-
- double deltaForceRMin = computeDeltaForce(cell,j,rmin);
- double deltaForceRMax = computeDeltaForce(cell,j,rmax);
- if(deltaForceRMax<0) {
- double effPoreRadius = rmax;
-// cerr<<"deltaForceRMax Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;
- return effPoreRadius;}
- else if(deltaForceRMin<0) {
- double effPoreRadius = bisection(cell,j,rmin,rmax);// cerr<<"1";//we suppose most cases should be this.
- return effPoreRadius;}
- else if( (deltaForceRMin>0) && (deltaForceRMin<deltaForceRMax) ) {
- double effPoreRadius = rmin;// cerr<<"2";
- return effPoreRadius;}
- else if(deltaForceRMin>deltaForceRMax) {
- double effPoreRadius = rmax;
-// cerr<<"WARNING! deltaForceRMin>deltaForceRMax. cellID: "<<cell->info().index<<"; deltaForceRMin="<<deltaForceRMin<<"; deltaForceRMax="<<deltaForceRMax<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;
- return effPoreRadius;}
-}
-double UnsaturatedEngine::computeDeltaForce(CellHandle cell,int j, double rCap)
-{
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- if (tri.is_infinite(cell->neighbor(j))) return 0;
-
- Vector3r posA = makeVector3r(cell->vertex(facetVertices[j][0])->point().point());
- Vector3r posB = makeVector3r(cell->vertex(facetVertices[j][1])->point().point());
- Vector3r posC = makeVector3r(cell->vertex(facetVertices[j][2])->point().point());
- double rA = sqrt(cell->vertex(facetVertices[j][0])->point().weight());
- double rB = sqrt(cell->vertex(facetVertices[j][1])->point().weight());
- double rC = sqrt(cell->vertex(facetVertices[j][2])->point().weight());
- double AB = (posA-posB).norm();
- double AC = (posA-posC).norm();
- double BC = (posB-posC).norm();
- double A = acos(((posB-posA).dot(posC-posA))/((posB-posA).norm()*(posC-posA).norm()));
- double B = acos(((posA-posB).dot(posC-posB))/((posA-posB).norm()*(posC-posB).norm()));
- double C = acos(((posA-posC).dot(posB-posC))/((posA-posC).norm()*(posB-posC).norm()));
- double rAB = 0.5*(AB-rA-rB); if (rAB<0) { rAB=0; }
- double rBC = 0.5*(BC-rB-rC); if (rBC<0) { rBC=0; }
- double rAC = 0.5*(AC-rA-rC); if (rAC<0) { rAC=0; }
-
- //In triangulation ArB,rCap is the radius of sphere r;
- double _AB = (pow((rA+rCap),2)+pow(AB,2)-pow((rB+rCap),2))/(2*(rA+rCap)*AB); if(_AB>1.0) {_AB=1.0;} if(_AB<-1.0) {_AB=-1.0;}
- double alphaAB = acos(_AB);
- double _BA = (pow((rB+rCap),2)+pow(AB,2)-pow((rA+rCap),2))/(2*(rB+rCap)*AB); if(_BA>1.0) {_BA=1.0;} if(_BA<-1.0) {_BA=-1.0;}
- double alphaBA = acos(_BA);
- double _ArB = (pow((rA+rCap),2)+pow((rB+rCap),2)-pow(AB,2))/(2*(rA+rCap)*(rB+rCap)); if(_ArB>1.0) {_ArB=1.0;} if(_ArB<-1.0) {_ArB=-1.0;}
- double alphaArB = acos(_ArB);
-
- double lengthLiquidAB = alphaArB*rCap;
- double AreaArB = 0.5*(rA+rCap)*(rB+rCap)*sin(alphaArB);
- double areaLiquidAB = AreaArB-0.5*alphaAB*pow(rA,2)-0.5*alphaBA*pow(rB,2)-0.5*alphaArB*pow(rCap,2);
-
- double _AC = (pow((rA+rCap),2)+pow(AC,2)-pow((rC+rCap),2))/(2*(rA+rCap)*AC); if(_AC>1.0) {_AC=1.0;} if(_AC<-1.0) {_AC=-1.0;}
- double alphaAC = acos(_AC);
- double _CA = (pow((rC+rCap),2)+pow(AC,2)-pow((rA+rCap),2))/(2*(rC+rCap)*AC); if(_CA>1.0) {_CA=1.0;} if(_CA<-1.0) {_CA=-1.0;}
- double alphaCA = acos(_CA);
- double _ArC = (pow((rA+rCap),2)+pow((rC+rCap),2)-pow(AC,2))/(2*(rA+rCap)*(rC+rCap)); if(_ArC>1.0) {_ArC=1.0;} if(_ArC<-1.0) {_ArC=-1.0;}
- double alphaArC = acos(_ArC);
-
- double lengthLiquidAC = alphaArC*rCap;
- double AreaArC = 0.5*(rA+rCap)*(rC+rCap)*sin(alphaArC);
- double areaLiquidAC = AreaArC-0.5*alphaAC*pow(rA,2)-0.5*alphaCA*pow(rC,2)-0.5*alphaArC*pow(rCap,2);
-
- double _BC = (pow((rB+rCap),2)+pow(BC,2)-pow((rC+rCap),2))/(2*(rB+rCap)*BC); if(_BC>1.0) {_BC=1.0;} if(_BC<-1.0) {_BC=-1.0;}
- double alphaBC = acos(_BC);
- double _CB = (pow((rC+rCap),2)+pow(BC,2)-pow((rB+rCap),2))/(2*(rC+rCap)*BC); if(_CB>1.0) {_CB=1.0;} if(_CB<-1.0) {_CB=-1.0;}
- double alphaCB = acos(_CB);
- double _BrC = (pow((rB+rCap),2)+pow((rC+rCap),2)-pow(BC,2))/(2*(rB+rCap)*(rC+rCap)); if(_BrC>1.0) {_BrC=1.0;} if(_BrC<-1.0) {_BrC=-1.0;}
- double alphaBrC = acos(_BrC);
-
- double lengthLiquidBC = alphaBrC*rCap;
- double AreaBrC = 0.5*(rB+rCap)*(rC+rCap)*sin(alphaBrC);
- double areaLiquidBC = AreaBrC-0.5*alphaBC*pow(rB,2)-0.5*alphaCB*pow(rC,2)-0.5*alphaBrC*pow(rCap,2);
-
- double areaCap = sqrt(cell->info().facetSurfaces[j].squared_length()) * cell->info().facetFluidSurfacesRatio[j];
- double areaPore = areaCap - areaLiquidAB - areaLiquidAC - areaLiquidBC;
-
- //FIXME:rethink here,areaPore Negative, Flat facets, do nothing ?
-// if(areaPore<0) {cerr<<"ERROR! areaPore Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;}
- double perimeterPore = lengthLiquidAB + lengthLiquidAC + lengthLiquidBC + (A - alphaAB - alphaAC)*rA + (B - alphaBA - alphaBC)*rB + (C - alphaCA - alphaCB)*rC;
-// if(perimeterPore<0) {cerr<<"ERROR! perimeterPore Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;}
-
- double deltaForce = perimeterPore - areaPore/rCap;//deltaForce=surfaceTension*(perimeterPore - areaPore/rCap)
- return deltaForce;
-}
-double UnsaturatedEngine::computeCellInterfacialArea(CellHandle cell, int j, double rCap)
-{
- double rInscribe = abs(solver->computeEffectiveRadius(cell, j));
- CellHandle cellh = CellHandle(cell);
- int facetNFictious = solver->detectFacetFictiousVertices (cellh,j);
- switch (facetNFictious) {
- case (0) : {
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- if (tri.is_infinite(cell->neighbor(j))) return 0;
-
- Vector3r posA = makeVector3r(cell->vertex(facetVertices[j][0])->point().point());
- Vector3r posB = makeVector3r(cell->vertex(facetVertices[j][1])->point().point());
- Vector3r posC = makeVector3r(cell->vertex(facetVertices[j][2])->point().point());
- double rA = sqrt(cell->vertex(facetVertices[j][0])->point().weight());
- double rB = sqrt(cell->vertex(facetVertices[j][1])->point().weight());
- double rC = sqrt(cell->vertex(facetVertices[j][2])->point().weight());
- double AB = (posA-posB).norm();
- double AC = (posA-posC).norm();
- double BC = (posB-posC).norm();
- double rAB = 0.5*(AB-rA-rB); if (rAB<0) { rAB=0; }
- double rBC = 0.5*(BC-rB-rC); if (rBC<0) { rBC=0; }
- double rAC = 0.5*(AC-rA-rC); if (rAC<0) { rAC=0; }
-
- //In triangulation ArB,rCap is the radius of sphere r;
- double _AB = (pow((rA+rCap),2)+pow(AB,2)-pow((rB+rCap),2))/(2*(rA+rCap)*AB); if(_AB>1.0) {_AB=1.0;} if(_AB<-1.0) {_AB=-1.0;}
- double alphaAB = acos(_AB);
- double _BA = (pow((rB+rCap),2)+pow(AB,2)-pow((rA+rCap),2))/(2*(rB+rCap)*AB); if(_BA>1.0) {_BA=1.0;} if(_BA<-1.0) {_BA=-1.0;}
- double alphaBA = acos(_BA);
- double _ArB = (pow((rA+rCap),2)+pow((rB+rCap),2)-pow(AB,2))/(2*(rA+rCap)*(rB+rCap)); if(_ArB>1.0) {_ArB=1.0;} if(_ArB<-1.0) {_ArB=-1.0;}
- double alphaArB = acos(_ArB);
- double AreaArB = 0.5*(rA+rCap)*(rB+rCap)*sin(alphaArB);
- double areaLiquidAB = AreaArB-0.5*alphaAB*pow(rA,2)-0.5*alphaBA*pow(rB,2)-0.5*alphaArB*pow(rCap,2);
-
- double _AC = (pow((rA+rCap),2)+pow(AC,2)-pow((rC+rCap),2))/(2*(rA+rCap)*AC); if(_AC>1.0) {_AC=1.0;} if(_AC<-1.0) {_AC=-1.0;}
- double alphaAC = acos(_AC);
- double _CA = (pow((rC+rCap),2)+pow(AC,2)-pow((rA+rCap),2))/(2*(rC+rCap)*AC); if(_CA>1.0) {_CA=1.0;} if(_CA<-1.0) {_CA=-1.0;}
- double alphaCA = acos(_CA);
- double _ArC = (pow((rA+rCap),2)+pow((rC+rCap),2)-pow(AC,2))/(2*(rA+rCap)*(rC+rCap)); if(_ArC>1.0) {_ArC=1.0;} if(_ArC<-1.0) {_ArC=-1.0;}
- double alphaArC = acos(_ArC);
- double AreaArC = 0.5*(rA+rCap)*(rC+rCap)*sin(alphaArC);
- double areaLiquidAC = AreaArC-0.5*alphaAC*pow(rA,2)-0.5*alphaCA*pow(rC,2)-0.5*alphaArC*pow(rCap,2);
-
- double _BC = (pow((rB+rCap),2)+pow(BC,2)-pow((rC+rCap),2))/(2*(rB+rCap)*BC); if(_BC>1.0) {_BC=1.0;} if(_BC<-1.0) {_BC=-1.0;}
- double alphaBC = acos(_BC);
- double _CB = (pow((rC+rCap),2)+pow(BC,2)-pow((rB+rCap),2))/(2*(rC+rCap)*BC); if(_CB>1.0) {_CB=1.0;} if(_CB<-1.0) {_CB=-1.0;}
- double alphaCB = acos(_CB);
- double _BrC = (pow((rB+rCap),2)+pow((rC+rCap),2)-pow(BC,2))/(2*(rB+rCap)*(rC+rCap)); if(_BrC>1.0) {_BrC=1.0;} if(_BrC<-1.0) {_BrC=-1.0;}
- double alphaBrC = acos(_BrC);
- double AreaBrC = 0.5*(rB+rCap)*(rC+rCap)*sin(alphaBrC);
- double areaLiquidBC = AreaBrC-0.5*alphaBC*pow(rB,2)-0.5*alphaCB*pow(rC,2)-0.5*alphaBrC*pow(rCap,2);
-
- double areaCap = sqrt(cell->info().facetSurfaces[j].squared_length()) * cell->info().facetFluidSurfacesRatio[j];
- double areaPore = areaCap - areaLiquidAB - areaLiquidAC - areaLiquidBC;//FIXME:areaPore may be out of range
- if (areaPore<0) cerr<<"areaPore:"<<areaPore<<" posA:"<<posA<<" rA:"<<rA<<" posB"<<posB<<" rB"<<rB<<" posC:"<<posC<<" rC"<<rC<<endl;
- if (areaPore>10) cerr<<"areaPore:"<<areaPore<<" posA:"<<posA<<" rA:"<<rA<<" posB"<<posB<<" rB"<<rB<<" posC:"<<posC<<" rC"<<rC<<endl;
- return areaPore;}; break;
- case (1) : { return Mathr::PI*pow(rInscribe,2); }; break;
- case (2) : { return Mathr::PI*pow(rInscribe,2); }; break;
- }
-}
-*/
// ------------------for checking----------
double UnsaturatedEngine::getRMin(CellHandle cell, int j)
{
