yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11445
[Branch ~yade-pkg/yade/git-trunk] Rev 3441: -rename variables for computePoreRadius(); clean code.
------------------------------------------------------------
revno: 3441
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Fri 2014-07-04 16:02:36 +0200
message:
-rename variables for computePoreRadius(); 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-03 13:35:28 +0000
+++ pkg/pfv/UnsaturatedEngine.cpp 2014-07-04 14:02:36 +0000
@@ -57,12 +57,13 @@
void updatePoreRadius();
void updateTotalCellVolume();
void updateVolumeCapillaryCell();
- double computeCellInterfacialArea(CellHandle cell, int j, double rCap);
+ double computeCellInterfacialArea(CellHandle cell, int j, double rC);
double computeEffPoreRadius(CellHandle cell, int j);
double computeEffPoreRadiusFine(CellHandle cell, int j);
double bisection(CellHandle cell, int j, double a, double b);
- double computeDeltaForce(CellHandle cell,int j, double rCap);
-
+ double computeTriRadian(double a, double b, double c);
+ double computeDeltaForce(CellHandle cell,int j, double rC);
+
void computeSolidLine();
void computeFacetPoreForcesWithCache(bool onlyCache=false);
void computeCapillaryForce() {computeFacetPoreForcesWithCache(false);}
@@ -87,21 +88,18 @@
void invade2();
void updateReservoirs2();
+ void saveVtk(const char* folder) {bool initT=solver->noCache; solver->noCache=false; solver->saveVtk(folder); solver->noCache=initT;}
+ void savePhaseVtk(const char* folder);
+
//record and test functions
void checkCellsConnection();
void checkEntryCapillaryPressure();
void checkLatticeNodeY(double y);
void checkReservoirInfo(int boundN);
void checkBoundingCellsInfo();
- void saveVtk(const char* folder) {bool initT=solver->noCache; solver->noCache=false; solver->saveVtk(folder); solver->noCache=initT;}
- void savePhaseVtk(const char* folder);
//temp functions
void initializeCellWindowsID();
double getWindowsSaturation(int i);
-// double getWindowsSaturation1(int i);
-// double getWindowsSaturation2(int i);
- double getRadiusMin(CellHandle cell, int j);
- void debugTemp();
bool checknoCache() {return solver->noCache;}
double getRMin(CellHandle cell, int j);
@@ -138,7 +136,6 @@
.def("testFunction",&UnsaturatedEngine::testFunction,"The playground for Chao's experiments.")
.def("checknoCache",&UnsaturatedEngine::checknoCache,"check noCache. (temp func.)")
.def("getWindowsSaturation",&UnsaturatedEngine::getWindowsSaturation,(boost::python::arg("windowsID")), "get saturation of windowsID")
- .def("debugTemp",&UnsaturatedEngine::debugTemp,"debug temp file.(temporary)")
.def("initializeCellWindowsID",&UnsaturatedEngine::initializeCellWindowsID,"Initialize cell windows index. A temp function for comparison with experiments, will delete soon")
.def("checkRCompare",&UnsaturatedEngine::checkRCompare,"debug R RMin RMax.")
)
@@ -514,77 +511,70 @@
return interfacialArea/totalCellVolume;
}
-double UnsaturatedEngine::computeCellInterfacialArea(CellHandle cell, int j, double rCap)
+double UnsaturatedEngine::computeCellInterfacialArea(CellHandle cell, int j, double rC)
{
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;
- }
+
+ if(facetNFictious==0) {
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ if (tri.is_infinite(cell->neighbor(j))) return 0;
+
+ Vector3r pos[3]; //solid pos
+ double r[3]; //solid radius
+ double rRc[3]; //r[i] + rC (rC: capillary radius)
+ double e[3]; //edges of triangulation
+ double rad[4][3]; //angle in radian
+
+ for (int i=0; i<3; i++) {
+ pos[i] = makeVector3r(cell->vertex(facetVertices[j][i])->point().point());
+ r[i] = sqrt(cell->vertex(facetVertices[j][i])->point().weight());
+ rRc[i] = r[i]+rC;
+ }
+
+ e[0] = (pos[1]-pos[2]).norm();
+ e[1] = (pos[2]-pos[0]).norm();
+ e[2] = (pos[1]-pos[0]).norm();
+
+ rad[3][0]=acos(((pos[1]-pos[0]).dot(pos[2]-pos[0]))/(e[2]*e[1]));
+ rad[3][1]=acos(((pos[2]-pos[1]).dot(pos[0]-pos[1]))/(e[0]*e[2]));
+ rad[3][2]=acos(((pos[0]-pos[2]).dot(pos[1]-pos[2]))/(e[1]*e[0]));
+
+ rad[0][0]=computeTriRadian(e[0],rRc[1],rRc[2]);
+ rad[0][1]=computeTriRadian(rRc[2],e[0],rRc[1]);
+ rad[0][2]=computeTriRadian(rRc[1],rRc[2],e[0]);
+
+ rad[1][0]=computeTriRadian(rRc[2],e[1],rRc[0]);
+ rad[1][1]=computeTriRadian(e[1],rRc[0],rRc[2]);
+ rad[1][2]=computeTriRadian(rRc[0],rRc[2],e[1]);
+
+ rad[2][0]=computeTriRadian(rRc[1],e[2],rRc[0]);
+ rad[2][1]=computeTriRadian(rRc[0],rRc[1],e[2]);
+ rad[2][2]=computeTriRadian(e[2],rRc[0],rRc[1]);
+
+ double sW0=0.5*rRc[1]*rRc[2]*sin(rad[0][0])-0.5*rad[0][0]*pow(rC,2)-0.5*rad[0][1]*pow(r[1],2)-0.5*rad[0][2]*pow(r[2],2) ;
+ double sW1=0.5*rRc[2]*rRc[0]*sin(rad[1][1])-0.5*rad[1][1]*pow(rC,2)-0.5*rad[1][2]*pow(r[2],2)-0.5*rad[1][0]*pow(r[0],2) ;
+ double sW2=0.5*rRc[0]*rRc[1]*sin(rad[2][2])-0.5*rad[2][2]*pow(rC,2)-0.5*rad[2][0]*pow(r[0],2)-0.5*rad[2][1]*pow(r[1],2) ;
+ double sW=sW0+sW1+sW2;
+ double sVoid=sqrt(cell->info().facetSurfaces[j].squared_length()) * cell->info().facetFluidSurfacesRatio[j];
+ double sInterface=sVoid-sW;
+ return sInterface;
+ }
+ else {
+ return Mathr::PI*pow(rInscribe,2);
+ }
}
double UnsaturatedEngine::computeEffPoreRadius(CellHandle cell, int j)
{
- double rInscribe = abs(solver->computeEffectiveRadius(cell, j));
+ double rInscribe = abs(solver->computeEffectiveRadius(cell, j));
CellHandle cellh = CellHandle(cell);
int facetNFictious = solver->detectFacetFictiousVertices (cellh,j);
- switch (facetNFictious) {
- case (0) : { return computeEffPoreRadiusFine(cell,j); }; break;
- case (1) : { return rInscribe; }; break;
- case (2) : { return rInscribe; }; break;
- }
+ double r;
+ if(facetNFictious==0) {r=computeEffPoreRadiusFine(cell,j);}
+ else r=rInscribe;
+ return r;
}
double UnsaturatedEngine::computeEffPoreRadiusFine(CellHandle cell, int j)
@@ -592,41 +582,49 @@
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 ?
+ Vector3r pos[3]; //solid pos
+ double r[3]; //solid radius
+ double e[3]; //edges of triangulation
+ double g[3]; //gap radius between solid
+
+ for (int i=0; i<3; i++) {
+ pos[i] = makeVector3r(cell->vertex(facetVertices[j][i])->point().point());
+ r[i] = sqrt(cell->vertex(facetVertices[j][i])->point().weight());
+ }
+
+ e[0] = (pos[1]-pos[2]).norm();
+ e[1] = (pos[2]-pos[0]).norm();
+ e[2] = (pos[1]-pos[0]).norm();
+ g[0] = ((e[0]-r[1]-r[2])>0) ? 0.5*(e[0]-r[1]-r[2]):0 ;
+ g[1] = ((e[1]-r[2]-r[0])>0) ? 0.5*(e[1]-r[2]-r[0]):0 ;
+ g[2] = ((e[2]-r[0]-r[1])>0) ? 0.5*(e[2]-r[0]-r[1]):0 ;
+
+ double rmin= (max(g[0],max(g[1],g[2]))==0) ? 1.0e-10:max(g[0],max(g[1],g[2])) ;
+ double rmax = abs(solver->computeEffectiveRadius(cell, j));
// 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 effPoreRadius;
+
+ if(deltaForceRMin>deltaForceRMax) {
+ effPoreRadius=rmax;
+ //cerr<<"ERROR! deltaForceRMin>deltaForceRMax. cellID: "<<cell->info().index<<"; deltaForceRMin="<<deltaForceRMin<<"; deltaForceRMax="<<deltaForceRMax<<". posA="<<pos[0]<<"; posB="<<pos[1]<<"; posC="<< pos[2]<<"; rA="<< r[0]<<"; rB="<< r[1]<<"; rC="<<r[2]<<endl;
+ }
+ else if(deltaForceRMax<0) {
+ effPoreRadius=rmax;
+ //cerr<<"deltaForceRMax<0. cellID: "<<cell->info().index<<"; deltaForceRMin="<<deltaForceRMin<<"; deltaForceRMax="<<deltaForceRMax<<". posA="<<pos[0]<<"; posB="<<pos[1]<<"; posC="<< pos[2]<<"; rA="<< r[0]<<"; rB="<< r[1]<<"; rC="<<r[2]<<endl;
+ }
+ else if(deltaForceRMin>0) {
+ effPoreRadius=rmin;
+ //cerr<<"0";
+ }
+ else {
+ effPoreRadius=bisection(cell,j,rmin,rmax);
+ //cerr<<"1";
+ }
+ return effPoreRadius;
}
-
double UnsaturatedEngine::bisection(CellHandle cell, int j, double a, double b)
{
double m = 0.5*(a+b);
@@ -640,6 +638,158 @@
else return m;
}
+//calculate radian with law of cosines. (solve $\alpha$)
+double UnsaturatedEngine::computeTriRadian(double a, double b, double c)
+{
+ double cosAlpha = (pow(b,2) + pow(c,2) - pow(a,2))/(2*b*c);
+ if (cosAlpha>1.0) {cosAlpha=1.0;} if (cosAlpha<-1.0) {cosAlpha=-1.0;}
+ double alpha = acos(cosAlpha);
+ return alpha;
+}
+
+double UnsaturatedEngine::computeDeltaForce(CellHandle cell,int j, double rC)
+{
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ if (tri.is_infinite(cell->neighbor(j))) return 0;
+
+ Vector3r pos[3]; //solid pos
+ double r[3]; //solid radius
+ double rRc[3]; //r[i] + rC (rC: capillary radius)
+ double e[3]; //edges of triangulation
+ double rad[4][3]; //angle in radian
+
+ for (int i=0; i<3; i++) {
+ pos[i] = makeVector3r(cell->vertex(facetVertices[j][i])->point().point());
+ r[i] = sqrt(cell->vertex(facetVertices[j][i])->point().weight());
+ rRc[i] = r[i]+rC;
+ }
+
+ e[0] = (pos[1]-pos[2]).norm();
+ e[1] = (pos[2]-pos[0]).norm();
+ e[2] = (pos[1]-pos[0]).norm();
+
+ rad[3][0]=acos(((pos[1]-pos[0]).dot(pos[2]-pos[0]))/(e[2]*e[1]));
+ rad[3][1]=acos(((pos[2]-pos[1]).dot(pos[0]-pos[1]))/(e[0]*e[2]));
+ rad[3][2]=acos(((pos[0]-pos[2]).dot(pos[1]-pos[2]))/(e[1]*e[0]));
+
+ rad[0][0]=computeTriRadian(e[0],rRc[1],rRc[2]);
+ rad[0][1]=computeTriRadian(rRc[2],e[0],rRc[1]);
+ rad[0][2]=computeTriRadian(rRc[1],rRc[2],e[0]);
+
+ rad[1][0]=computeTriRadian(rRc[2],e[1],rRc[0]);
+ rad[1][1]=computeTriRadian(e[1],rRc[0],rRc[2]);
+ rad[1][2]=computeTriRadian(rRc[0],rRc[2],e[1]);
+
+ rad[2][0]=computeTriRadian(rRc[1],e[2],rRc[0]);
+ rad[2][1]=computeTriRadian(rRc[0],rRc[1],e[2]);
+ rad[2][2]=computeTriRadian(e[2],rRc[0],rRc[1]);
+
+ double lNW = (rad[0][0]+rad[1][1]+rad[2][2])*rC;
+ double lNS = (rad[3][0]-rad[1][0]-rad[2][0])*r[0] + (rad[3][1]-rad[2][1]-rad[0][1])*r[1] + (rad[3][2]-rad[1][2]-rad[0][2])*r[2] ;
+ double lInterface=lNW+lNS;
+
+ double sW0=0.5*rRc[1]*rRc[2]*sin(rad[0][0])-0.5*rad[0][0]*pow(rC,2)-0.5*rad[0][1]*pow(r[1],2)-0.5*rad[0][2]*pow(r[2],2) ;
+ double sW1=0.5*rRc[2]*rRc[0]*sin(rad[1][1])-0.5*rad[1][1]*pow(rC,2)-0.5*rad[1][2]*pow(r[2],2)-0.5*rad[1][0]*pow(r[0],2) ;
+ double sW2=0.5*rRc[0]*rRc[1]*sin(rad[2][2])-0.5*rad[2][2]*pow(rC,2)-0.5*rad[2][0]*pow(r[0],2)-0.5*rad[2][1]*pow(r[1],2) ;
+ double sW=sW0+sW1+sW2;
+ double sVoid=sqrt(cell->info().facetSurfaces[j].squared_length()) * cell->info().facetFluidSurfacesRatio[j];
+ double sInterface=sVoid-sW;
+
+ double deltaF = lInterface - sInterface/rC;//deltaF=surfaceTension*(perimeterPore - areaPore/rCap)
+ return deltaF;
+}
+
+void UnsaturatedEngine::savePhaseVtk(const char* folder)
+{
+// RTriangulation& Tri = T[solver->noCache?(!currentTes):currentTes].Triangulation();
+ RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ static unsigned int number = 0;
+ char filename[80];
+ mkdir(folder, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ sprintf(filename,"%s/out_%d.vtk",folder,number++);
+ int firstReal=-1;
+
+ //count fictious vertices and cells
+ solver->vtkInfiniteVertices=solver->vtkInfiniteCells=0;
+ FiniteCellsIterator cellEnd = Tri.finite_cells_end();
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
+ bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
+ if (!isDrawable) solver->vtkInfiniteCells+=1;
+ }
+ for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+ if (!v->info().isReal()) solver->vtkInfiniteVertices+=1;
+ else if (firstReal==-1) firstReal=solver->vtkInfiniteVertices;}
+
+ basicVTKwritter vtkfile((unsigned int) Tri.number_of_vertices()-solver->vtkInfiniteVertices, (unsigned int) Tri.number_of_finite_cells()-solver->vtkInfiniteCells);
+
+ vtkfile.open(filename,"test");
+
+ vtkfile.begin_vertices();
+ double x,y,z;
+ for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+ if (v->info().isReal()){
+ x = (double)(v->point().point()[0]);
+ y = (double)(v->point().point()[1]);
+ z = (double)(v->point().point()[2]);
+ vtkfile.write_point(x,y,z);}
+ }
+ vtkfile.end_vertices();
+
+ vtkfile.begin_cells();
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
+ bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
+ if (isDrawable){vtkfile.write_cell(cell->vertex(0)->info().id()-firstReal, cell->vertex(1)->info().id()-firstReal, cell->vertex(2)->info().id()-firstReal, cell->vertex(3)->info().id()-firstReal);}
+ }
+ vtkfile.end_cells();
+
+ vtkfile.begin_data("Phase",CELL_DATA,SCALARS,FLOAT);
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
+ bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
+ if (isDrawable){vtkfile.write_data(!cell->info().isAirReservoir);}
+ }
+ 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();
@@ -706,27 +856,90 @@
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)
{
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- if (tri.is_infinite(cell->neighbor(j))) {return 0;cerr<<"tri.is_infinite(cell->neighbor(j)"<<endl;}
-
- 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; }
+ if (tri.is_infinite(cell->neighbor(j))) return 0;
+ Vector3r pos[3]; //solid pos
+ double r[3]; //solid radius
+ double e[3]; //edges of triangulation
+ double g[3]; //gap radius between solid
+
+ for (int i=0; i<3; i++) {
+ pos[i] = makeVector3r(cell->vertex(facetVertices[j][i])->point().point());
+ r[i] = sqrt(cell->vertex(facetVertices[j][i])->point().weight());
+ }
+
+ e[0] = (pos[1]-pos[2]).norm();
+ e[1] = (pos[2]-pos[0]).norm();
+ e[2] = (pos[1]-pos[0]).norm();
+ g[0] = ((e[0]-r[1]-r[2])>0) ? 0.5*(e[0]-r[1]-r[2]):0 ;
+ g[1] = ((e[1]-r[2]-r[0])>0) ? 0.5*(e[1]-r[2]-r[0]):0 ;
+ g[2] = ((e[2]-r[0]-r[1])>0) ? 0.5*(e[2]-r[0]-r[1]):0 ;
+
+ double rmin= (max(g[0],max(g[1],g[2]))==0) ? 1.0e-10:max(g[0],max(g[1],g[2])) ;
return rmin;
}
double UnsaturatedEngine::getRMax(CellHandle cell, int j)
@@ -744,10 +957,9 @@
file<<"cellID "<<"j "<<"rmin "<<"r "<<"rmax "<<endl;
FiniteCellsIterator cellEnd = tri.finite_cells_end();
for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
- file << cell->info().index << " " <<"0"<< " " <<getRMin(cell,0) << " " << computeEffPoreRadius(cell,0)<< " " <<getRMax(cell,0) << endl;
- file << cell->info().index << " " <<"1"<< " " <<getRMin(cell,1) << " " << computeEffPoreRadius(cell,1)<< " " <<getRMax(cell,1) << endl;
- file << cell->info().index << " " <<"2"<< " " <<getRMin(cell,2) << " " << computeEffPoreRadius(cell,2)<< " " <<getRMax(cell,2) << endl;
- file << cell->info().index << " " <<"3"<< " " <<getRMin(cell,3) << " " << computeEffPoreRadius(cell,3)<< " " <<getRMax(cell,3) << endl;
+ for (int facet=0; facet<4; facet++) {
+ file << cell->info().index << " " <<facet<< " " <<getRMin(cell,facet) << " " << computeEffPoreRadius(cell,facet)<< " " <<getRMax(cell,facet) << endl;
+ }
}
file.close();
@@ -755,17 +967,9 @@
file<<"cellID "<<"j "<<"rmin "<<"r "<<"rmax "<<endl;
FiniteCellsIterator cellEnd1 = tri.finite_cells_end();
for ( FiniteCellsIterator cell1 = tri.finite_cells_begin(); cell1 != cellEnd1; cell1++ ) {
- if(getRMin(cell1,0)==computeEffPoreRadius(cell1,0)) {
- file << cell1->info().index << " " <<"0"<< " " <<getRMin(cell1,0) << " " << computeEffPoreRadius(cell1,0)<< " " <<getRMax(cell1,0) << endl;
- }
- if(getRMin(cell1,1)==computeEffPoreRadius(cell1,1)) {
- file << cell1->info().index << " " <<"1"<< " " <<getRMin(cell1,1) << " " << computeEffPoreRadius(cell1,1)<< " " <<getRMax(cell1,1) << endl;
- }
- if(getRMin(cell1,2)==computeEffPoreRadius(cell1,2)) {
- file << cell1->info().index << " " <<"2"<< " " <<getRMin(cell1,2) << " " << computeEffPoreRadius(cell1,2)<< " " <<getRMax(cell1,2) << endl;
- }
- if(getRMin(cell1,3)==computeEffPoreRadius(cell1,3)) {
- file << cell1->info().index << " " <<"3"<< " " <<getRMin(cell1,3) << " " << computeEffPoreRadius(cell1,3)<< " " <<getRMax(cell1,3) << endl;
+ for (int facet=0; facet<4; facet++) {
+ if(getRMin(cell1,facet)==computeEffPoreRadius(cell1,facet))
+ file << cell1->info().index << " " <<facet<< " " <<getRMin(cell1,facet) << " " << computeEffPoreRadius(cell1,facet)<< " " <<getRMax(cell1,facet) << endl;
}
}
file.close();
@@ -774,17 +978,9 @@
file<<"cellID "<<"j "<<"rmin "<<"r "<<"rmax "<<endl;
FiniteCellsIterator cellEnd2 = tri.finite_cells_end();
for ( FiniteCellsIterator cell2 = tri.finite_cells_begin(); cell2 != cellEnd2; cell2++ ) {
- if(getRMax(cell2,0)==computeEffPoreRadius(cell2,0)) {
- file << cell2->info().index << " " <<"0"<< " " <<getRMin(cell2,0) << " " << computeEffPoreRadius(cell2,0)<< " " <<getRMax(cell2,0) << endl;
- }
- if(getRMax(cell2,1)==computeEffPoreRadius(cell2,1)) {
- file << cell2->info().index << " " <<"1"<< " " <<getRMin(cell2,1) << " " << computeEffPoreRadius(cell2,1)<< " " <<getRMax(cell2,1) << endl;
- }
- if(getRMax(cell2,2)==computeEffPoreRadius(cell2,2)) {
- file << cell2->info().index << " " <<"2"<< " " <<getRMin(cell2,2) << " " << computeEffPoreRadius(cell2,2)<< " " <<getRMax(cell2,2) << endl;
- }
- if(getRMax(cell2,3)==computeEffPoreRadius(cell2,3)) {
- file << cell2->info().index << " " <<"3"<< " " <<getRMin(cell2,3) << " " << computeEffPoreRadius(cell2,3)<< " " <<getRMax(cell2,3) << endl;
+ for (int facet=0; facet<4; facet++) {
+ if(getRMax(cell2,facet)==computeEffPoreRadius(cell2,facet))
+ file << cell2->info().index << " " <<facet<< " " <<getRMin(cell2,facet) << " " << computeEffPoreRadius(cell2,facet)<< " " <<getRMax(cell2,facet) << endl;
}
}
file.close();
@@ -793,18 +989,11 @@
file<<"cellID "<<"j "<<"rmin "<<"r "<<"rmax "<<endl;
FiniteCellsIterator cellEnd3 = tri.finite_cells_end();
for ( FiniteCellsIterator cell3 = tri.finite_cells_begin(); cell3 != cellEnd3; cell3++ ) {
- if( (getRMax(cell3,0)>computeEffPoreRadius(cell3,0)) && (getRMin(cell3,0)<computeEffPoreRadius(cell3,0)) ) {
- file << cell3->info().index << " " <<"0"<< " " <<getRMin(cell3,0) << " " << computeEffPoreRadius(cell3,0)<< " " <<getRMax(cell3,0) << endl;
- }
- if( (getRMax(cell3,1)>computeEffPoreRadius(cell3,1)) && (getRMin(cell3,1)<computeEffPoreRadius(cell3,1)) ) {
- file << cell3->info().index << " " <<"1"<< " " <<getRMin(cell3,1) << " " << computeEffPoreRadius(cell3,1)<< " " <<getRMax(cell3,1) << endl;
- }
- if( (getRMax(cell3,2)>computeEffPoreRadius(cell3,2)) && (getRMin(cell3,2)<computeEffPoreRadius(cell3,2)) ) {
- file << cell3->info().index << " " <<"2"<< " " <<getRMin(cell3,2) << " " << computeEffPoreRadius(cell3,2)<< " " <<getRMax(cell3,2) << endl;
- }
- if( (getRMax(cell3,3)>computeEffPoreRadius(cell3,3)) && (getRMin(cell3,3)<computeEffPoreRadius(cell3,3)) ) {
- file << cell3->info().index << " " <<"3"<< " " <<getRMin(cell3,3) << " " << computeEffPoreRadius(cell3,3)<< " " <<getRMax(cell3,3) << endl;
- }
+ for (int facet=0; facet<4; facet++) {
+ if( (getRMax(cell3,facet)>computeEffPoreRadius(cell3,facet)) && (getRMin(cell3,facet)<computeEffPoreRadius(cell3,facet)) ) {
+ file << cell3->info().index << " " <<facet<< " " <<getRMin(cell3,facet) << " " << computeEffPoreRadius(cell3,facet)<< " " <<getRMax(cell3,facet) << endl;
+ }
+ }
}
file.close();
@@ -812,24 +1001,14 @@
file<<"cellID "<<"j "<<"rmin "<<"r "<<"rmax "<<endl;
FiniteCellsIterator cellEnd4 = tri.finite_cells_end();
for ( FiniteCellsIterator cell4 = tri.finite_cells_begin(); cell4 != cellEnd4; cell4++ ) {
- if(getRMax(cell4,0)<getRMin(cell4,0)) {
- file << cell4->info().index << " " <<"0"<< " " <<getRMin(cell4,0) << " " << computeEffPoreRadius(cell4,0)<< " " <<getRMax(cell4,0) << endl;
- }
- if(getRMax(cell4,1)<getRMin(cell4,1)) {
- file << cell4->info().index << " " <<"1"<< " " <<getRMin(cell4,1) << " " << computeEffPoreRadius(cell4,1)<< " " <<getRMax(cell4,1) << endl;
- }
- if(getRMax(cell4,2)<getRMin(cell4,2)) {
- file << cell4->info().index << " " <<"2"<< " " <<getRMin(cell4,2) << " " << computeEffPoreRadius(cell4,2)<< " " <<getRMax(cell4,2) << endl;
- }
- if(getRMax(cell4,3)<getRMin(cell4,3)) {
- file << cell4->info().index << " " <<"3"<< " " <<getRMin(cell4,3) << " " << computeEffPoreRadius(cell4,3)<< " " <<getRMax(cell4,3) << endl;
- }
+ for(int facet=0; facet<4; facet++) {
+ if(getRMax(cell4,facet)<getRMin(cell4,facet)) {
+ file << cell4->info().index << " " <<facet<< " " <<getRMin(cell4,facet) << " " << computeEffPoreRadius(cell4,facet)<< " " <<getRMax(cell4,facet) << endl;
+ }
+ }
}
file.close();
-
}
-// --------------------------------------
-
void UnsaturatedEngine::checkCellsConnection()
{
@@ -854,10 +1033,9 @@
FiniteCellsIterator cellEnd = tri.finite_cells_end();
for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
if (tri.is_infinite(cell)) continue;
- file << cell->info().index << " " <<cell->neighbor(0)->info().index << " " << surfaceTension/cell->info().poreRadius[0] << endl;
- file << cell->info().index << " " <<cell->neighbor(1)->info().index << " " << surfaceTension/cell->info().poreRadius[1] << endl;
- file << cell->info().index << " " <<cell->neighbor(2)->info().index << " " << surfaceTension/cell->info().poreRadius[2] << endl;
- file << cell->info().index << " " <<cell->neighbor(3)->info().index << " " << surfaceTension/cell->info().poreRadius[3] << endl;
+ for (int facet=0; facet<4; facet++) {
+ file << cell->info().index << " " <<cell->neighbor(facet)->info().index << " " << surfaceTension/cell->info().poreRadius[facet] << endl;
+ }
}
file.close();
}
@@ -945,95 +1123,6 @@
}
}
-void UnsaturatedEngine::savePhaseVtk(const char* folder)
-{
-// RTriangulation& Tri = T[solver->noCache?(!currentTes):currentTes].Triangulation();
- RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
- static unsigned int number = 0;
- char filename[80];
- mkdir(folder, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
- sprintf(filename,"%s/out_%d.vtk",folder,number++);
- int firstReal=-1;
-
- //count fictious vertices and cells
- solver->vtkInfiniteVertices=solver->vtkInfiniteCells=0;
- FiniteCellsIterator cellEnd = Tri.finite_cells_end();
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
- bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
- if (!isDrawable) solver->vtkInfiniteCells+=1;
- }
- for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
- if (!v->info().isReal()) solver->vtkInfiniteVertices+=1;
- else if (firstReal==-1) firstReal=solver->vtkInfiniteVertices;}
-
- basicVTKwritter vtkfile((unsigned int) Tri.number_of_vertices()-solver->vtkInfiniteVertices, (unsigned int) Tri.number_of_finite_cells()-solver->vtkInfiniteCells);
-
- vtkfile.open(filename,"test");
-
- vtkfile.begin_vertices();
- double x,y,z;
- for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
- if (v->info().isReal()){
- x = (double)(v->point().point()[0]);
- y = (double)(v->point().point()[1]);
- z = (double)(v->point().point()[2]);
- vtkfile.write_point(x,y,z);}
- }
- vtkfile.end_vertices();
-
- vtkfile.begin_cells();
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
- bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
- if (isDrawable){vtkfile.write_cell(cell->vertex(0)->info().id()-firstReal, cell->vertex(1)->info().id()-firstReal, cell->vertex(2)->info().id()-firstReal, cell->vertex(3)->info().id()-firstReal);}
- }
- vtkfile.end_cells();
-
- vtkfile.begin_data("Phase",CELL_DATA,SCALARS,FLOAT);
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
- bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
- if (isDrawable){vtkfile.write_data(!cell->info().isAirReservoir);}
- }
- vtkfile.end_data();
-}
-
-//----temp function for Vahid Joekar-Niasar's data----
-double UnsaturatedEngine::getRadiusMin(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; }
- return rmin;
-}
-
-void UnsaturatedEngine::debugTemp()
-{
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- ofstream file;
- file.open("bugTemp.txt");
- FiniteCellsIterator cellEnd = tri.finite_cells_end();
- for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
- if (tri.is_infinite(cell)) continue;
- for(int i=0; i<4; i++) {
- file << cell->info().index << " " <<cell->info().solidLine[i][0]<<" " <<cell->info().solidLine[i][1]<<" " <<cell->info().solidLine[i][2]<<" " <<cell->info().solidLine[i][3]<<endl;
- }
- }
- file.close();
-}//----------end temp function for Vahid Joekar-Niasar's data (clear later)---------------------
-
//----------temp functions for comparison with experiment-----------------------
void UnsaturatedEngine::initializeCellWindowsID()
{
