yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11390
[Branch ~yade-pkg/yade/git-trunk] Rev 3386: -test commit
------------------------------------------------------------
revno: 3386
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Tue 2013-11-19 09:56:36 +0100
message:
-test commit
modified:
pkg/dem/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/dem/UnsaturatedEngine.cpp'
--- pkg/dem/UnsaturatedEngine.cpp 2013-08-23 16:29:53 +0000
+++ pkg/dem/UnsaturatedEngine.cpp 2013-11-19 08:56:36 +0000
@@ -47,8 +47,8 @@
setPositionsBuffer(true);//copy sphere positions in a buffer...
Build_Triangulation(P_zero,solver);//create a triangulation and initialize pressure in the elements, everything will be contained in "solver"
initializeCellIndex(solver);//initialize cell index
- getPoreRadius(solver);//save all pore radii before invade
- }
+ initializePoreRadius(solver);//save all pore radii before invade
+ }
solver->noCache = false;
}
@@ -56,7 +56,7 @@
{
//This will be used later
/*
- Build_Triangulation();initializeCellIndex();Initialize_volumes(solver);getPoreRadius();
+ Build_Triangulation();initializeCellIndex();Initialize_volumes(solver);initializePoreRadius();
*/
}
@@ -191,6 +191,18 @@
template<class Cellhandle>
double UnsaturatedEngine::computeEffPoreRadius(Cellhandle cell, int j)
{
+ double rInscribe = abs(solver->Compute_EffectiveRadius(cell, j));
+ Cell_handle cellh = Cell_handle(cell);
+ int facetNFictious = solver->detectFacetFictiousVertices (cellh,j);
+ switch (facetNFictious) {
+ case (0) : { return computeEffPoreRadiusNormal(cell,j); }; break;
+ case (1) : { return rInscribe; }; break;
+ case (2) : { return rInscribe; }; break;
+ }
+}
+template<class Cellhandle>
+double UnsaturatedEngine::computeEffPoreRadiusNormal(Cellhandle cell, int j)
+{
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
if (tri.is_infinite(cell->neighbor(j))) return 0;
@@ -209,80 +221,36 @@
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 Area_ABC = 0.5 * ((posB-posA).cross(posC-posA)).norm();
-// double Area_SolidA = 0.5*A*pow(rA,2);
-// double Area_SolidB = 0.5*B*pow(rB,2);
-// double Area_SolidC = 0.5*C*pow(rC,2);
- double rmin = max(rAB,max(rBC,rAC));
- if (rmin==0) { rmin= 1.0e-10; }
- double rmax = abs(solver->Compute_EffectiveRadius(cell, j));
-
- if((computeDeltaForce(cell,j,rmax)<0)||(rmin>rmax)) {
- double EffPoreRadius = rmax;
-// cerr<<"0";
- return EffPoreRadius;
- }
- else if((computeDeltaForce(cell,j,rmin)<0)&&(computeDeltaForce(cell,j,rmax)>0)) {
- double effPoreRadius = bisection(cell,j,rmin,rmax);
-// cerr<<"1";
- return effPoreRadius;
- }
- else if((computeDeltaForce(cell,j,rmin)>0)&&(computeDeltaForce(cell,j,rmin)<computeDeltaForce(cell,j,rmax))) {
- double EffPoreRadius = rmin;
-// cerr<<"2";
- return EffPoreRadius;
- }
- else if ((computeDeltaForce(cell,j,rmin)>computeDeltaForce(cell,j,rmax))&&(computeDeltaForce(cell,j,rmax)>0)) {
- double EffPoreRadius = rmax;
-// cerr<<"3";
- return EffPoreRadius;
- }
-// if (rmin>rmax) {cerr<<"rmin: "<< rmin << " " << "rmax" << rmax << " rA="<<rA<<" rB="<<rB <<" rC="<<rC<<endl;}
-
-// if (((Area_ABC-Area_SolidA-Area_SolidB-Area_SolidC)<0)||(rmin>rmax)) {
-// double EffPoreRadius = rmax;//for cells close to boundary spheres, the effPoreRadius set to inscribe radius.
-// cerr<<"1";
-// return EffPoreRadius;
-// }
-// if( ( computeDeltaForce(cell,j,rmin)>0 ) && ( computeDeltaForce(cell,j,rmin)<computeDeltaForce(cell,j,rmax)) ) {
-// double EffPoreRadius = rmin;
-// // cerr<<"2 ";
-// return EffPoreRadius;
-// }
-// else if( ( computeDeltaForce(cell,j,rmin)<0 ) && ( computeDeltaForce(cell,j,rmax)>0) ) {
-// double effPoreRadius = bisection(cell,j,rmin,rmax);
-// // cerr<<"3 ";
-// return effPoreRadius;
-// }
-// else if( ( computeDeltaForce(cell,j,rmin)<computeDeltaForce(cell,j,rmax) ) && ( computeDeltaForce(cell,j,rmax)<0) ) {
-// double EffPoreRadius = rmax;
-// // cerr<<"4 ";
-// return EffPoreRadius;
-// }
-// else if( ( computeDeltaForce(cell,j,rmin)>computeDeltaForce(cell,j,rmax) ) && ( computeDeltaForce(cell,j,rmax)>0) ) {
-// double EffPoreRadius = rmax;
-// // cerr<<"5 ";
-// return EffPoreRadius;
-// }
-// else if( ( computeDeltaForce(cell,j,rmin)>0 ) && ( computeDeltaForce(cell,j,rmax)<0) ) {
-// double effPoreRadius = bisection(cell,j,rmin,rmax);
-// // cerr<<"6 ";
-// return effPoreRadius;
-// }
-// else if( ( computeDeltaForce(cell,j,rmin)> computeDeltaForce(cell,j,rmax) ) && (computeDeltaForce(cell,j,rmin)<0) ) {
-// double EffPoreRadius = rmin;
-// // cerr<<"7 ";
-// return EffPoreRadius;
-// }
-// else {
-// double EffPoreRadius = rmax;
-// // cerr<<"8";
-// return EffPoreRadius;
-// }
+ double rmin = max(rAB,max(rBC,rAC)); if (rmin==0) { rmin= 1.0e-10; }
+ double rmax = abs(solver->Compute_EffectiveRadius(cell, j));//rmin>rmax ?
+ if(rmin>rmax) { cerr<<"rmin>rmax. rmin="<<rmin<<" ,rmax="<<rmax<<endl; }
+
+ double deltaForce_rMin = computeDeltaForce(cell,j,rmin);
+ double deltaForce_rMax = computeDeltaForce(cell,j,rmax);
+ if(deltaForce_rMax<0) {
+ double EffPoreRadius = rmax;
+ cerr<<"deltaForce_rMax Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;
+ return EffPoreRadius;
+ }
+ else if(deltaForce_rMin<0) {
+ double effPoreRadius = bisection(cell,j,rmin,rmax);// cerr<<"1";
+// double deltaF=computeDeltaForce(cell,j,effPoreRadius);cerr<<"deltaForce_rMin: "<<deltaForce_rMin<<" deltaForce_rMax: "<<deltaForce_rMax<<" deltaF: "<<deltaF<<endl;//test rMin<R<rMax ?
+ return effPoreRadius;
+ }
+ else if( (deltaForce_rMin>0) && (deltaForce_rMin<deltaForce_rMax) ) {
+ double EffPoreRadius = rmin;// cerr<<"2";
+// cerr<<"posA: "<<posA<<" rA: "<<rA<<" posB: "<<posB<<" rB: "<<rB<<" posC: "<<posC<<" rC: "<<rC<<endl;
+ return EffPoreRadius;
+ }
+ else if(deltaForce_rMin>deltaForce_rMax) {
+ double EffPoreRadius = rmax;
+ cerr<<"deltaForce_rMin>deltaForce_rMax. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;
+ return EffPoreRadius;
+ }
}
-template<class Cellhandle>
+template<class Cellhandle>//waiting check
double UnsaturatedEngine::bisection(Cellhandle cell, int j, double a, double b)
{
double m = 0.5*(a+b);
@@ -302,7 +270,6 @@
template<class Cellhandle>
double UnsaturatedEngine::computeDeltaForce(Cellhandle cell,int j, double rcap)
{
- double surface_tension = surfaceTension; //Surface Tension in contact with air at 20 Degrees Celsius is:0.0728(N/m)
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
if (tri.is_infinite(cell->neighbor(j))) return 0;
@@ -314,75 +281,62 @@
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 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; }
-// double Area_ABC = 0.5 * ((posB-posA).cross(posC-posA)).norm();
-// double Area_SolidA = 0.5*A*pow(rA,2);//??Area for boundary ERROR ?? FIXME
-// double Area_SolidB = 0.5*B*pow(rB,2);
-// double Area_SolidC = 0.5*C*pow(rC,2);
-
- //for triangulation ArB,rcap is the radius of sphere r; Note: (pow((rA+rcap),2)+pow(AB,2)-pow((rB+rcap),2))/(2*(rA+rcap)*AB) maybe >1, bug here!
- double _AB = (pow((rA+rcap),2)+pow(AB,2)-pow((rB+rcap),2))/(2*(rA+rcap)*AB); if(_AB>1.0) {/*cerr<<"cellID="<<cell->info().index<<" rA="<<rA<<" rB="<<rB<<" rC="<<rC<<endl;*//*cerr<<"_AB>1.0 _AB="<<_AB<<endl;*/ _AB=1.0;} if(_AB<-1.0) {cerr<<"_AB<-1.0 _AB="<<_AB<<endl; _AB=-1.0;}
+/*
+ double Area_ABC = 0.5 * ((posB-posA).cross(posC-posA)).norm();
+ double Area_SolidA = 0.5*A*pow(rA,2);//??Area for boundary ERROR ?? FIXME
+ double Area_SolidB = 0.5*B*pow(rB,2);
+ double Area_SolidC = 0.5*C*pow(rC,2);
+*/
+ //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) {/*cerr<<"cellID="<<cell->info().index<<" rA="<<rA<<" rB="<<rB<<" rC="<<rC<<endl;*//*cerr<<"_BA>1.0 _BA="<<_BA<<endl;*/ _BA=1.0;} if(_BA<-1.0) {cerr<<"_BA<-1.0 _BA="<<_BA<<endl; _BA=-1.0;}
+ 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 D=alphaAB + alphaBA + alphaArB; cerr<<D<<" ";
+// double D1=alphaAB + alphaBA + alphaArB; cerr<<D<<" ";
double length_liquidAB = alphaArB*rcap;
double AreaArB = 0.5*(rA+rcap)*(rB+rcap)*sin(alphaArB);
double Area_liquidAB = AreaArB-0.5*alphaAB*pow(rA,2)-0.5*alphaBA*pow(rB,2)-0.5*alphaArB*pow(rcap,2);
- //for triangulation ArC, rcap is the radius of sphere r;
- double _AC = (pow((rA+rcap),2)+pow(AC,2)-pow((rC+rcap),2))/(2*(rA+rcap)*AC); if(_AC>1.0) {/*cerr<<"cellID="<<cell->info().index<<" rA="<<rA<<" rB="<<rB<<" rC="<<rC<<endl;*//*cerr<<"_AC>1.0 _AC="<<_AC<<endl;*/ _AC=1.0;} if(_AC<-1.0) {cerr<<"_AC<-1.0 _AC="<<_AC<<endl; _AC=-1.0;}
+ 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) {/*cerr<<"cellID="<<cell->info().index<<" rA="<<rA<<" rB="<<rB<<" rC="<<rC<<endl;*//*cerr<<"_CA>1.0 _CA="<<_CA<<endl;*/ _CA=1.0;} if(_CA<-1.0) {cerr<<"_CA<-1.0 _CA="<<_CA<<endl; _CA=-1.0;}
+ 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 D=alphaAC + alphaCA + alphaArC; cerr<<D<<" ";
+// double D2=alphaAC + alphaCA + alphaArC; cerr<<D<<" ";
double length_liquidAC = alphaArC*rcap;
double AreaArC = 0.5*(rA+rcap)*(rC+rcap)*sin(alphaArC);
double Area_liquidAC = AreaArC-0.5*alphaAC*pow(rA,2)-0.5*alphaCA*pow(rC,2)-0.5*alphaArC*pow(rcap,2);
- //for triangulation BrC, rcap is the radius of sphere r;
- double _BC = (pow((rB+rcap),2)+pow(BC,2)-pow((rC+rcap),2))/(2*(rB+rcap)*BC); if(_BC>1.0) {/*cerr<<"cellID="<<cell->info().index<<" rA="<<rA<<" rB="<<rB<<" rC="<<rC<<endl;*//*cerr<<"_BC>1.0 _BC="<<_BC<<endl;*/ _BC=1.0;} if(_BC<-1.0) {cerr<<"_BC<-1.0 _BC="<<_BC<<endl; _BC=-1.0;}
+ 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) {/*cerr<<"cellID="<<cell->info().index<<" rA="<<rA<<" rB="<<rB<<" rC="<<rC<<endl;*//*cerr<<"_CB>1.0 _CB="<<_CB<<endl;*/ _CB=1.0;} if(_CB<-1.0) {cerr<<"_CB<-1.0 _CB="<<_CB<<endl; _CB=-1.0;}
+ 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 D=alphaBC + alphaCB + alphaBrC; cerr<<D<<" ";
+// double D3=alphaBC + alphaCB + alphaBrC; cerr<<D<<" ";
double length_liquidBC = alphaBrC*rcap;
double AreaBrC = 0.5*(rB+rcap)*(rC+rcap)*sin(alphaBrC);
double Area_liquidBC = 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 - Area_liquidAB - Area_liquidAC - Area_liquidBC;
-// double Area_pore = Area_ABC - Area_liquidAB - Area_liquidAC - Area_liquidBC - Area_SolidA - Area_SolidB - Area_SolidC;
-// if(areaCap<0) {cerr<<"areaCap="<<areaCap<<" rA="<<rA<<" rB="<<rB<<" rC="<<rC<<endl;}
- if(areaPore<0) {cerr<<"areaPore="<<areaPore<<" rA="<<rA<<" rB="<<rB<<" rC="<<rC<<endl;}
-// double Area_pore1= Area_ABC - Area_SolidA - Area_SolidB - Area_SolidC;
-// if(Area_pore1<0) {cerr<<"Area_pore1="<<Area_pore1<<" rAB="<<rAB<<" rAC="<<rAC<<" rBC="<<rBC<<" Area_ABC="<<Area_ABC<<" Area_SolidA="<<Area_SolidA<<" Area_SolidB="<<Area_SolidB<<" Area_SolidC="<<Area_SolidC<<endl;} //here because big boundary spheres
-// if((Area_pore1>0)&&(Area_pore<0))
-// { double rAB = 0.5*(AB-rA-rB);
-// double rBC = 0.5*(BC-rB-rC);
-// double rAC = 0.5*(AC-rA-rC);
-// double rInscr=abs(solver->Compute_EffectiveRadius(cell, j));
-// if((rAB>0)&&(rBC>0)&&(rAC>0))
-// { cerr<<"rA= "<<rA<<" rB="<<rB<<" rC="<<rC<<" rInscr="<<rInscr<<" Area_pore="<<Area_pore<<endl;
-// }
-// } //here also because big boundary spheres
- double Perimeter_pore = length_liquidAB + length_liquidAC + length_liquidBC + (A - alphaAB - alphaAC)*rA + (B - alphaBA - alphaBC)*rB + (C - alphaCA - alphaCB)*rC; if(Perimeter_pore<0) {/*cerr<<"cellID="<<cell->info().index<<" rA="<<rA<<" rB="<<rB<<" rC="<<rC<<endl;*//*double t1=(A - alphaAB - alphaAC)*rA; double t2= (B - alphaBA - alphaBC)*rB; double t3=(C - alphaCA - alphaCB)*rC; cerr<<"cellID="<<cell->info().index<<" Perimeter_pore<0 Perimeter_pore="<<Perimeter_pore<<" liquidAB="<<length_liquidAB<<" liquidAC="<<length_liquidAC<<" liquidBC="<<length_liquidBC<<" t1="<<t1<<" t2="<<t2<<" t3="<<t3<<endl;*/}
- if(Perimeter_pore<0) {cerr<<"Perimeter_pore="<<Perimeter_pore<<" rA= "<<rA<<" rB="<<rB<<" rC="<<rC<<endl;}
- double deltaForce = surface_tension*Perimeter_pore - areaPore*(surface_tension/rcap);
-// if((Area_pore<0)&&(deltaForce>0)) {cerr<<"rA= "<<rA<<" rB="<<rB<<" rC="<<rC<<endl;}
+ if(areaPore<0) {cerr<<"areaPore Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;}
+ //areaPore Negative, Flat facets, do nothing here ?
+ double perimeterPore = length_liquidAB + length_liquidAC + length_liquidBC + (A - alphaAB - alphaAC)*rA + (B - alphaBA - alphaBC)*rB + (C - alphaCA - alphaCB)*rC;
+ if(perimeterPore<0) {cerr<<"perimeterPore Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;}
+// double surface_tension = surfaceTension; //Surface Tension in contact with air at 20 Degrees Celsius is:0.0728(N/m)
+// double deltaForce = surface_tension*perimeterPore - areaPore*(surface_tension/rcap);
+ double deltaForce = perimeterPore - areaPore/rcap;
return deltaForce;
}
@@ -541,13 +495,6 @@
for (int k=0;k<6;k++) flow->boundary ( *flow->boundsIds[k] ).useMaxMin = boundaryUseMaxMin[k];
-// if ( flow->y_min_id>=0 ) flow->boundary ( flow->y_min_id ).useMaxMin = boundaryUseMaxMin[ymin];
-// if ( flow->y_max_id>=0 ) flow->boundary ( flow->y_max_id ).useMaxMin = boundaryUseMaxMin[ymax];
-// if ( flow->x_max_id>=0 ) flow->boundary ( flow->x_max_id ).useMaxMin = boundaryUseMaxMin[xmax];
-// if ( flow->x_min_id>=0 ) flow->boundary ( flow->x_min_id ).useMaxMin = boundaryUseMaxMin[xmin];
-// if ( flow->z_max_id>=0 ) flow->boundary ( flow->z_max_id ).useMaxMin = boundaryUseMaxMin[zmax];
-// if ( flow->z_min_id>=0 ) flow->boundary ( flow->z_min_id ).useMaxMin = boundaryUseMaxMin[zmin];
-
flow->Corner_min = CGT::Point ( flow->x_min, flow->y_min, flow->z_min );
flow->Corner_max = CGT::Point ( flow->x_max, flow->y_max, flow->z_max );
@@ -567,81 +514,81 @@
}
}
-// template<class Solver>
-// void UnsaturatedEngine::AddBoundary ( Solver& flow )
-// {
-// vector<posData>& buffer = positionBufferCurrent;
-// solver->x_min = Mathr::MAX_REAL, solver->x_max = -Mathr::MAX_REAL, solver->y_min = Mathr::MAX_REAL, solver->y_max = -Mathr::MAX_REAL, solver->z_min = Mathr::MAX_REAL, solver->z_max = -Mathr::MAX_REAL;
-// FOREACH ( const posData& b, buffer ) {
-// if ( !b.exists ) continue;
-// if ( b.isSphere ) {
-// const Real& rad = b.radius;
-// const Real& x = b.pos[0];
-// const Real& y = b.pos[1];
-// const Real& z = b.pos[2];
-// flow->x_min = min ( flow->x_min, x-rad );
-// flow->x_max = max ( flow->x_max, x+rad );
-// flow->y_min = min ( flow->y_min, y-rad );
-// flow->y_max = max ( flow->y_max, y+rad );
-// flow->z_min = min ( flow->z_min, z-rad );
-// flow->z_max = max ( flow->z_max, z+rad );
-// }
-// }
-// //FIXME id_offset must be set correctly, not the case here (always 0), then we need walls first or it will fail
-// id_offset = flow->T[flow->currentTes].max_id+1;
-// flow->id_offset = id_offset;
-// flow->SectionArea = ( flow->x_max - flow->x_min ) * ( flow->z_max-flow->z_min );
-// flow->Vtotale = ( flow->x_max-flow->x_min ) * ( flow->y_max-flow->y_min ) * ( flow->z_max-flow->z_min );
-// flow->y_min_id=wallBottomId;
-// flow->y_max_id=wallTopId;
-// flow->x_max_id=wallRightId;
-// flow->x_min_id=wallLeftId;
-// flow->z_min_id=wallBackId;
-// flow->z_max_id=wallFrontId;
-//
-// if ( flow->y_min_id>=0 ) flow->boundary ( flow->y_min_id ).useMaxMin = BOTTOM_Boundary_MaxMin;
-// if ( flow->y_max_id>=0 ) flow->boundary ( flow->y_max_id ).useMaxMin = TOP_Boundary_MaxMin;
-// if ( flow->x_max_id>=0 ) flow->boundary ( flow->x_max_id ).useMaxMin = RIGHT_Boundary_MaxMin;
-// if ( flow->x_min_id>=0 ) flow->boundary ( flow->x_min_id ).useMaxMin = LEFT_Boundary_MaxMin;
-// if ( flow->z_max_id>=0 ) flow->boundary ( flow->z_max_id ).useMaxMin = FRONT_Boundary_MaxMin;
-// if ( flow->z_min_id>=0 ) flow->boundary ( flow->z_min_id ).useMaxMin = BACK_Boundary_MaxMin;
-//
-// //FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
-// flow->boundsIds[0]= &flow->x_min_id;
-// flow->boundsIds[1]= &flow->x_max_id;
-// flow->boundsIds[2]= &flow->y_min_id;
-// flow->boundsIds[3]= &flow->y_max_id;
-// flow->boundsIds[4]= &flow->z_min_id;
-// flow->boundsIds[5]= &flow->z_max_id;
-//
-// flow->Corner_min = CGT::Point ( flow->x_min, flow->y_min, flow->z_min );
-// flow->Corner_max = CGT::Point ( flow->x_max, flow->y_max, flow->z_max );
-//
-// if ( Debug ) {
-// cout << "Section area = " << flow->SectionArea << endl;
-// cout << "Vtotale = " << flow->Vtotale << endl;
-// cout << "x_min = " << flow->x_min << endl;
-// cout << "x_max = " << flow->x_max << endl;
-// cout << "y_max = " << flow->y_max << endl;
-// cout << "y_min = " << flow->y_min << endl;
-// cout << "z_min = " << flow->z_min << endl;
-// cout << "z_max = " << flow->z_max << endl;
-// cout << endl << "Adding Boundary------" << endl;
-// }
-// //assign BCs types and values
-// BoundaryConditions ( flow );
-//
-// double center[3];
-// for ( int i=0; i<6; i++ ) {
-// if ( *flow->boundsIds[i]<0 ) continue;
-// CGT::Vecteur Normal ( normal[i].x(), normal[i].y(), normal[i].z() );
-// if ( flow->boundary ( *flow->boundsIds[i] ).useMaxMin ) flow->AddBoundingPlane ( true, Normal, *flow->boundsIds[i],5000.0 );
-// else {
-// for ( int h=0;h<3;h++ ) center[h] = buffer[*flow->boundsIds[i]].pos[h];
-// flow->AddBoundingPlane ( center, wall_thickness, Normal,*flow->boundsIds[i],5000.0 );
-// }
-// }
-// }
+/*template<class Solver>
+void UnsaturatedEngine::AddBoundary ( Solver& flow )
+{
+ vector<posData>& buffer = positionBufferCurrent;
+ solver->x_min = Mathr::MAX_REAL, solver->x_max = -Mathr::MAX_REAL, solver->y_min = Mathr::MAX_REAL, solver->y_max = -Mathr::MAX_REAL, solver->z_min = Mathr::MAX_REAL, solver->z_max = -Mathr::MAX_REAL;
+ FOREACH ( const posData& b, buffer ) {
+ if ( !b.exists ) continue;
+ if ( b.isSphere ) {
+ const Real& rad = b.radius;
+ const Real& x = b.pos[0];
+ const Real& y = b.pos[1];
+ const Real& z = b.pos[2];
+ flow->x_min = min ( flow->x_min, x-rad );
+ flow->x_max = max ( flow->x_max, x+rad );
+ flow->y_min = min ( flow->y_min, y-rad );
+ flow->y_max = max ( flow->y_max, y+rad );
+ flow->z_min = min ( flow->z_min, z-rad );
+ flow->z_max = max ( flow->z_max, z+rad );
+ }
+ }
+ //FIXME id_offset must be set correctly, not the case here (always 0), then we need walls first or it will fail
+ id_offset = flow->T[flow->currentTes].max_id+1;
+ flow->id_offset = id_offset;
+ flow->SectionArea = ( flow->x_max - flow->x_min ) * ( flow->z_max-flow->z_min );
+ flow->Vtotale = ( flow->x_max-flow->x_min ) * ( flow->y_max-flow->y_min ) * ( flow->z_max-flow->z_min );
+ flow->y_min_id=wallBottomId;
+ flow->y_max_id=wallTopId;
+ flow->x_max_id=wallRightId;
+ flow->x_min_id=wallLeftId;
+ flow->z_min_id=wallBackId;
+ flow->z_max_id=wallFrontId;
+
+ if ( flow->y_min_id>=0 ) flow->boundary ( flow->y_min_id ).useMaxMin = BOTTOM_Boundary_MaxMin;
+ if ( flow->y_max_id>=0 ) flow->boundary ( flow->y_max_id ).useMaxMin = TOP_Boundary_MaxMin;
+ if ( flow->x_max_id>=0 ) flow->boundary ( flow->x_max_id ).useMaxMin = RIGHT_Boundary_MaxMin;
+ if ( flow->x_min_id>=0 ) flow->boundary ( flow->x_min_id ).useMaxMin = LEFT_Boundary_MaxMin;
+ if ( flow->z_max_id>=0 ) flow->boundary ( flow->z_max_id ).useMaxMin = FRONT_Boundary_MaxMin;
+ if ( flow->z_min_id>=0 ) flow->boundary ( flow->z_min_id ).useMaxMin = BACK_Boundary_MaxMin;
+
+ //FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
+ flow->boundsIds[0]= &flow->x_min_id;
+ flow->boundsIds[1]= &flow->x_max_id;
+ flow->boundsIds[2]= &flow->y_min_id;
+ flow->boundsIds[3]= &flow->y_max_id;
+ flow->boundsIds[4]= &flow->z_min_id;
+ flow->boundsIds[5]= &flow->z_max_id;
+
+ flow->Corner_min = CGT::Point ( flow->x_min, flow->y_min, flow->z_min );
+ flow->Corner_max = CGT::Point ( flow->x_max, flow->y_max, flow->z_max );
+
+ if ( Debug ) {
+ cout << "Section area = " << flow->SectionArea << endl;
+ cout << "Vtotale = " << flow->Vtotale << endl;
+ cout << "x_min = " << flow->x_min << endl;
+ cout << "x_max = " << flow->x_max << endl;
+ cout << "y_max = " << flow->y_max << endl;
+ cout << "y_min = " << flow->y_min << endl;
+ cout << "z_min = " << flow->z_min << endl;
+ cout << "z_max = " << flow->z_max << endl;
+ cout << endl << "Adding Boundary------" << endl;
+ }
+ //assign BCs types and values
+ BoundaryConditions ( flow );
+
+ double center[3];
+ for ( int i=0; i<6; i++ ) {
+ if ( *flow->boundsIds[i]<0 ) continue;
+ CGT::Vecteur Normal ( normal[i].x(), normal[i].y(), normal[i].z() );
+ if ( flow->boundary ( *flow->boundsIds[i] ).useMaxMin ) flow->AddBoundingPlane ( true, Normal, *flow->boundsIds[i],5000.0 );
+ else {
+ for ( int h=0;h<3;h++ ) center[h] = buffer[*flow->boundsIds[i]].pos[h];
+ flow->AddBoundingPlane ( center, wall_thickness, Normal,*flow->boundsIds[i],5000.0 );
+ }
+ }
+}*/
template<class Solver>
void UnsaturatedEngine::Triangulate ( Solver& flow )
@@ -702,7 +649,7 @@
}
template<class Solver>
-void UnsaturatedEngine::getPoreRadius(Solver& flow)
+void UnsaturatedEngine::initializePoreRadius(Solver& flow)
{
RTriangulation& Tri = flow->T[flow->currentTes].Triangulation();
Finite_cells_iterator cell_end = Tri.finite_cells_end();
@@ -822,9 +769,7 @@
{
Initialize_volumes(solver);
Real volume = abs( cell->info().volume() ) - solver->volumeSolidPore(cell);
- if (cell->info().volume()<0) { cerr << "cell ID: " << cell->info().index << " cell volume: " << cell->info().volume() << endl; }
- if (solver->volumeSolidPore(cell)<0) { cerr << "cell ID: " << cell->info().index << " volumeSolidPore: " << solver->volumeSolidPore(cell) << endl; }
- if (volume<0) { cerr<<"cell ID: " << cell->info().index << "cell volume: " << cell->info().volume() << " volumeSolidPore: " << solver->volumeSolidPore(cell) << endl; }
+ if (volume<0) { cerr<<"volumeCapillaryCell Negative. cell ID: " << cell->info().index << "cell volume: " << cell->info().volume() << " volumeSolidPore: " << solver->volumeSolidPore(cell) << endl; }
return volume;
}
@@ -849,7 +794,7 @@
}
template<class Solver>
-void UnsaturatedEngine::saveListOfNodes(Solver& flow)
+void UnsaturatedEngine::saveListNodes(Solver& flow)
{
ofstream file;
file.open("ListOfNodes.txt");
@@ -863,10 +808,10 @@
}
template<class Solver>
-void UnsaturatedEngine::saveListOfConnections(Solver& flow)
+void UnsaturatedEngine::saveListConnection(Solver& flow)
{
ofstream file;
- file.open("ListOfConnections.txt");
+ file.open("ListConnection.txt");
file << "#List of Connections \n";
file << "Cell_ID" << " " << "NeighborCell_ID" << " " << "EntryValue" << " " << "poreRadius" << " " << "poreArea" << " " << "porePerimeter" << endl;
double surface_tension = surfaceTension ; //Surface Tension in contact with air at 20 Degrees Celsius is:0.0728(N/m)
@@ -1003,15 +948,15 @@
{
/*Here,boundsIds according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
For drainage and imbibition situations, we only care about top(3) and bottom(2) boundaries.(0-xmin,1-xmax,2-ymin,3-ymax,4-zmin,5-zmax)
- And in python script, don't forget set: Flow_imposed_TOP_Boundary=False,Flow_imposed_BOTTOM_Boundary=False, */
+ And in python script, don't forget set: bndCondIsPressure=[0,0,1,1,0,0], */
if(flow->boundingCells[3].size()==0) {
- cerr << "please set Flow_imposed_TOP_Boundary=False" << endl;
+ cerr << "please set bndCondIsPressure=[0,0,1,1,0,0]." << endl;
}
else {
vector<Cell_handle>::iterator it = flow->boundingCells[3].begin();
ofstream file;
file.open("ListAdjacentCellsTopBoundary.txt");
- file << "#List of Cells IDs adjacent top boundary \n";
+ file << "#IDs of boundary pore bodies, which are connecting water reservoir. \n";
for ( it ; it != flow->boundingCells[3].end(); it++) {
if ((*it)->info().index == 0) continue;
// cerr << (*it)->info().index << " ";
@@ -1025,13 +970,13 @@
void UnsaturatedEngine::saveListAdjCellsBottomBound(Solver& flow)
{
if(flow->boundingCells[2].size()==0) {
- cerr << "please set Flow_imposed_BOTTOM_Boundary=False"<< endl;
+ cerr << "please set bndCondIsPressure=[0,0,1,1,0,0]."<< endl;
}
else {
vector<Cell_handle>::iterator it = flow->boundingCells[2].begin();
ofstream file;
file.open("ListAdjacentCellsBottomBoundary.txt");
- file << "#List of Cells IDs adjacent bottom boundary \n";
+ file << "#IDs of boundary pore bodies, which are connecting air reservoir. \n";
for ( it ; it != flow->boundingCells[2].end(); it++) {
if ((*it)->info().index == 0) continue;
file << (*it)->info().index << endl;
@@ -1153,176 +1098,196 @@
void UnsaturatedEngine::clearImposedPressure ( Solver& flow ) { flow->imposedP.clear(); flow->IPCells.clear();}
//----tempt function for Vahid Joekar-Niasar's data----
-template<class Cellhandle>//ERROR! check later
+template<class Cellhandle >
+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 = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point());
+ Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point());
+ Vector3r posC = makeVector3r2(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; }
+
+ double rmin = max(rAB,max(rBC,rAC)); if (rmin==0) { rmin= 1.0e-10; }
+ return rmin;
+}
+template<class Solver>
+void UnsaturatedEngine::debugTemp(Solver& flow)
+{
+ ofstream file;
+ file.open("bugTemp.txt");
+// file<<"cellID "<<"deltaForce "<<"inscribeRadius*PI*tension "<<endl;
+ Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
+ for ( Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
+ if (flow->T[flow->currentTes].Triangulation().is_infinite(cell)) continue;
+ file << cell->info().index << " " <<cell->info().poreRadius[0]<<" "<<getRadiusMin(cell,0)<<" "<<abs(flow->Compute_EffectiveRadius(cell, 0))<<endl;
+ file << cell->info().index << " " <<cell->info().poreRadius[1]<<" "<<getRadiusMin(cell,1)<<" "<<abs(flow->Compute_EffectiveRadius(cell, 1))<<endl;
+ file << cell->info().index << " " <<cell->info().poreRadius[2]<<" "<<getRadiusMin(cell,2)<<" "<<abs(flow->Compute_EffectiveRadius(cell, 2))<<endl;
+ file << cell->info().index << " " <<cell->info().poreRadius[3]<<" "<<getRadiusMin(cell,3)<<" "<<abs(flow->Compute_EffectiveRadius(cell, 3))<<endl;
+ }
+ file.close();
+}
+
+template<class Cellhandle>
Real UnsaturatedEngine::computePoreArea(Cellhandle cell, int j)
{
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- if (tri.is_infinite(cell->neighbor(j))) return 0;
-
- Real rcap = cell->info().poreRadius[j];
- Vector3r posA = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point());
- Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point());
- Vector3r posC = makeVector3r2(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; }
- double Area_ABC = 0.5 * ((posB-posA).cross(posC-posA)).norm();
- double Area_SolidA = 0.5*A*pow(rA,2);
- double Area_SolidB = 0.5*B*pow(rB,2);
- double Area_SolidC = 0.5*C*pow(rC,2);
-
- double rmin = max(rAB,max(rBC,rAC));
- if (rmin==0) { rmin= 1.0e-10; }
- double rmax = abs(solver->Compute_EffectiveRadius(cell, j));
-
- Real poreArea = 0;
- if (abs(rcap-rmax)<1.0e-6) {
- poreArea = Mathr::PI*pow(rmax,2);
- return poreArea;
- }
-// else if (cell->info().poreRadius[j] > rmax) {cerr<<"poreRadius Error: "<<cell->info().poreRadius[j];}
- else {
- //for triangulation ArB,rcap is the radius of sphere r; Note: (pow((rA+rcap),2)+pow(AB,2)-pow((rB+rcap),2))/(2*(rA+rcap)*AB) maybe >1, bug here!
- double _AB = pow((rA+rcap),2)+pow(AB,2)-pow((rB+rcap),2)/(2*(rA+rcap)*AB); //ERROR HERE!
- if (_AB>1) { _AB=1.0; }
+ double rInscribe = abs(solver->Compute_EffectiveRadius(cell, j));
+ Cell_handle cellh = Cell_handle(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 = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point());
+ Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point());
+ Vector3r posC = makeVector3r2(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;}
+
+ double rcap = cell->info().poreRadius[j];
+ //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) { _BA=1.0; }
+ 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 betaAB = 0.5*Mathr::PI - alphaAB;
- double betaBA = 0.5*Mathr::PI - alphaBA;
- double length_liquidAB = (betaAB+betaBA)*rcap;
- double AreaArB = 0.5*AB*(rA+rcap)*sin(alphaAB);
- double Area_liquidAB = AreaArB-0.5*alphaAB*pow(rA,2)-0.5*alphaBA*pow(rB,2)-0.5*(betaAB+betaBA)*pow(rcap,2);
+ 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 D1=alphaAB + alphaBA + alphaArB; cerr<<D<<" ";
+ double length_liquidAB = alphaArB*rcap;
+ double AreaArB = 0.5*(rA+rcap)*(rB+rcap)*sin(alphaArB);
+ double Area_liquidAB = AreaArB-0.5*alphaAB*pow(rA,2)-0.5*alphaBA*pow(rB,2)-0.5*alphaArB*pow(rcap,2);
- //for triangulation ArC, rcap is the radius of sphere r;
- double _AC = pow((rA+rcap),2)+pow(AC,2)-pow((rC+rcap),2)/(2*(rA+rcap)*AC);
- if (_AC>1) { _AC=1.0; }
+ 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) { _CA=1.0; }
+ 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 betaAC = 0.5*Mathr::PI - alphaAC;
- double betaCA = 0.5*Mathr::PI - alphaCA;
- double length_liquidAC = (betaAC+betaCA)*rcap;
- double AreaArC = 0.5*AC*(rA+rcap)*sin(alphaAC);
- double Area_liquidAC = AreaArC-0.5*alphaAC*pow(rA,2)-0.5*alphaCA*pow(rC,2)-0.5*(betaAC+betaCA)*pow(rcap,2);
+ 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 D2=alphaAC + alphaCA + alphaArC; cerr<<D<<" ";
+ double length_liquidAC = alphaArC*rcap;
+ double AreaArC = 0.5*(rA+rcap)*(rC+rcap)*sin(alphaArC);
+ double Area_liquidAC = AreaArC-0.5*alphaAC*pow(rA,2)-0.5*alphaCA*pow(rC,2)-0.5*alphaArC*pow(rcap,2);
- //for triangulation BrC, rcap is the radius of sphere r;
- double _BC = pow((rB+rcap),2)+pow(BC,2)-pow((rC+rcap),2)/(2*(rB+rcap)*BC);
- if (_BC>1) { _BC=1.0; }
+ 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) { _CB=1.0; }
+ 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 betaBC = 0.5*Mathr::PI - alphaBC;
- double betaCB = 0.5*Mathr::PI - alphaCB;
- double length_liquidBC = (betaBC+betaCB)*rcap;
- double AreaBrC = 0.5*BC*(rB+rcap)*sin(alphaBC);
- double Area_liquidBC = AreaBrC-0.5*alphaBC*pow(rB,2)-0.5*alphaCB*pow(rC,2)-0.5*(betaBC+betaCB)*pow(rcap,2);
+ 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 D3=alphaBC + alphaCB + alphaBrC; cerr<<D<<" ";
+ double length_liquidBC = alphaBrC*rcap;
+ double AreaBrC = 0.5*(rB+rcap)*(rC+rcap)*sin(alphaBrC);
+ double Area_liquidBC = AreaBrC-0.5*alphaBC*pow(rB,2)-0.5*alphaCB*pow(rC,2)-0.5*alphaBrC*pow(rcap,2);
- double poreArea = Area_ABC - Area_liquidAB - Area_liquidAC - Area_liquidBC - Area_SolidA - Area_SolidB - Area_SolidC;
- if (poreArea<0) { poreArea=Mathr::PI*pow(rmax,2); }
- return poreArea;
- }
+ double areaCap = sqrt(cell->info().facetSurfaces[j].squared_length()) * cell->info().facetFluidSurfacesRatio[j];
+ double areaPore = areaCap - Area_liquidAB - Area_liquidAC - Area_liquidBC;
+ if(areaPore<0) {cerr<<"areaPore Negative. areaPore="<<areaPore<<endl ;areaPore=Mathr::PI*pow(rInscribe,2);}
+ return areaPore;
+ }; break;
+ case (1) : { return Mathr::PI*pow(rInscribe,2); }; break;
+ case (2) : { return Mathr::PI*pow(rInscribe,2); }; break;
+ }
}
+
template<class Cellhandle>
Real UnsaturatedEngine::computePorePerimeter(Cellhandle cell, int j)
{
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- if (tri.is_infinite(cell->neighbor(j))) return 0;
-
- Real rcap = cell->info().poreRadius[j];
- Vector3r posA = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point());
- Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point());
- Vector3r posC = makeVector3r2(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; }
- double Area_ABC = 0.5 * ((posB-posA).cross(posC-posA)).norm();
- double Area_SolidA = 0.5*A*pow(rA,2);
- double Area_SolidB = 0.5*B*pow(rB,2);
- double Area_SolidC = 0.5*C*pow(rC,2);
-
- double rmin = max(rAB,max(rBC,rAC));
- if (rmin==0) { rmin= 1.0e-10; }
- double rmax = abs(solver->Compute_EffectiveRadius(cell, j));
-
- Real porePerimeter = 0;
- if (abs(rcap-rmax)<1.0e-6) {
-// cerr<<"rcap: "<<rcap <<" "<<"rmax: "<<rmax<<endl;
- porePerimeter = 2*Mathr::PI*rmax;
- return porePerimeter;
- }
- else {
- //for triangulation ArB,rcap is the radius of sphere r; Note: (pow((rA+rcap),2)+pow(AB,2)-pow((rB+rcap),2))/(2*(rA+rcap)*AB) maybe >1, bug here!
- double _AB = pow((rA+rcap),2)+pow(AB,2)-pow((rB+rcap),2)/(2*(rA+rcap)*AB);
- if (_AB>1) { _AB=1.0; }
+ double rInscribe = abs(solver->Compute_EffectiveRadius(cell, j));
+ Cell_handle cellh = Cell_handle(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 = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point());
+ Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point());
+ Vector3r posC = makeVector3r2(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;}
+
+ double rcap = cell->info().poreRadius[j];
+ //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) { _BA=1.0; }
+ 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 betaAB = 0.5*Mathr::PI - alphaAB;
- double betaBA = 0.5*Mathr::PI - alphaBA;
- double length_liquidAB = (betaAB+betaBA)*rcap;
- double AreaArB = 0.5*AB*(rA+rcap)*sin(alphaAB);
- double Area_liquidAB = AreaArB-0.5*alphaAB*pow(rA,2)-0.5*alphaBA*pow(rB,2)-0.5*(betaAB+betaBA)*pow(rcap,2);
+ 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 D1=alphaAB + alphaBA + alphaArB; cerr<<D<<" ";
+ double length_liquidAB = alphaArB*rcap;
+ double AreaArB = 0.5*(rA+rcap)*(rB+rcap)*sin(alphaArB);
+ double Area_liquidAB = AreaArB-0.5*alphaAB*pow(rA,2)-0.5*alphaBA*pow(rB,2)-0.5*alphaArB*pow(rcap,2);
- //for triangulation ArC, rcap is the radius of sphere r;
- double _AC = pow((rA+rcap),2)+pow(AC,2)-pow((rC+rcap),2)/(2*(rA+rcap)*AC);
- if (_AC>1) { _AC=1.0; }
+ 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) { _CA=1.0; }
+ 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 betaAC = 0.5*Mathr::PI - alphaAC;
- double betaCA = 0.5*Mathr::PI - alphaCA;
- double length_liquidAC = (betaAC+betaCA)*rcap;
- double AreaArC = 0.5*AC*(rA+rcap)*sin(alphaAC);
- double Area_liquidAC = AreaArC-0.5*alphaAC*pow(rA,2)-0.5*alphaCA*pow(rC,2)-0.5*(betaAC+betaCA)*pow(rcap,2);
+ 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 D2=alphaAC + alphaCA + alphaArC; cerr<<D<<" ";
+ double length_liquidAC = alphaArC*rcap;
+ double AreaArC = 0.5*(rA+rcap)*(rC+rcap)*sin(alphaArC);
+ double Area_liquidAC = AreaArC-0.5*alphaAC*pow(rA,2)-0.5*alphaCA*pow(rC,2)-0.5*alphaArC*pow(rcap,2);
- //for triangulation BrC, rcap is the radius of sphere r;
- double _BC = pow((rB+rcap),2)+pow(BC,2)-pow((rC+rcap),2)/(2*(rB+rcap)*BC);
- if (_BC>1) { _BC=1.0; }
+ 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) { _CB=1.0; }
+ 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 betaBC = 0.5*Mathr::PI - alphaBC;
- double betaCB = 0.5*Mathr::PI - alphaCB;
- double length_liquidBC = (betaBC+betaCB)*rcap;
- double AreaBrC = 0.5*BC*(rB+rcap)*sin(alphaBC);
- double Area_liquidBC = AreaBrC-0.5*alphaBC*pow(rB,2)-0.5*alphaCB*pow(rC,2)-0.5*(betaBC+betaCB)*pow(rcap,2);
-
- double porePerimeter = length_liquidAB + length_liquidAC + length_liquidBC + (A - alphaAB - alphaAC)*rA + (B - alphaBA - alphaBC)*rB + (C - alphaCA - alphaCB)*rC;
- if (porePerimeter<0) { porePerimeter = 2*Mathr::PI*rmax; }
- return porePerimeter;
- }
+ 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 D3=alphaBC + alphaCB + alphaBrC; cerr<<D<<" ";
+ double length_liquidBC = alphaBrC*rcap;
+ double AreaBrC = 0.5*(rB+rcap)*(rC+rcap)*sin(alphaBrC);
+ double Area_liquidBC = AreaBrC-0.5*alphaBC*pow(rB,2)-0.5*alphaCB*pow(rC,2)-0.5*alphaBrC*pow(rcap,2);
+ double perimeterPore = length_liquidAB + length_liquidAC + length_liquidBC + (A - alphaAB - alphaAC)*rA + (B - alphaBA - alphaBC)*rB + (C - alphaCA - alphaCB)*rC;
+ if(perimeterPore<0) {cerr<<"perimeterPore Negative. perimeterPore="<<perimeterPore<<endl;perimeterPore=2.0*Mathr::PI*rInscribe;}
+ return perimeterPore;
+ }; break;
+ case (1) : { return 2.0*Mathr::PI*rInscribe; }; break;
+ case (2) : { return 2.0*Mathr::PI*rInscribe; }; break;
+ }
}
+
template<class Solver>
void UnsaturatedEngine::savePoreBodyInfo(Solver& flow)
{
ofstream file;
file.open("PoreBodyInfo.txt");
- file << "#List of pore bodies position (or Voronoi centers) and size (volume) \n";
+ file << "#pore bodies position (or Voronoi centers) and size (volume) \n";
file << "Cell_ID " << " x " << " y " << " z " << " volume "<< endl;
Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
for ( Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
@@ -1331,7 +1296,26 @@
}
file.close();
}
+
template<class Solver>
+void UnsaturatedEngine::savePoreThroatInfo(Solver& flow)
+{
+ ofstream file;
+ file.open("PoreThroatInfo.txt");
+ file << "#pore throat area, inscribeRadius and perimeter. \n";
+ file << "Cell_ID " << " NeighborCell_ID "<<" area " << " inscribeRadius " << " perimeter " << endl;
+ Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
+ for ( Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
+ if (flow->T[flow->currentTes].Triangulation().is_infinite(cell)) continue;
+ file << cell->info().index <<" "<< cell->neighbor(0)->info().index <<" "<< computePoreArea(cell, 0) <<" "<< abs(flow->Compute_EffectiveRadius(cell, 0)) <<" "<< computePorePerimeter(cell,0) <<endl;
+ file << cell->info().index <<" "<< cell->neighbor(1)->info().index <<" "<< computePoreArea(cell, 1) <<" "<< abs(flow->Compute_EffectiveRadius(cell, 1)) <<" "<< computePorePerimeter(cell,1) <<endl;
+ file << cell->info().index <<" "<< cell->neighbor(2)->info().index <<" "<< computePoreArea(cell, 2) <<" "<< abs(flow->Compute_EffectiveRadius(cell, 2)) <<" "<< computePorePerimeter(cell,2) <<endl;
+ file << cell->info().index <<" "<< cell->neighbor(3)->info().index <<" "<< computePoreArea(cell, 3) <<" "<< abs(flow->Compute_EffectiveRadius(cell, 3)) <<" "<< computePorePerimeter(cell,3) <<endl;
+ }
+ file.close();
+}
+
+/*template<class Solver>//test computeDeltaForce function
void UnsaturatedEngine::debugTemp(Solver& flow)
{
double surface_tension = surfaceTension; //Surface Tension in contact with air at 20 Degrees Celsius is:0.0728(N/m)
@@ -1341,58 +1325,17 @@
Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
for ( Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
if (flow->T[flow->currentTes].Triangulation().is_infinite(cell)) continue;
- double rmax0= abs(solver->Compute_EffectiveRadius(cell, 0));double D0=surface_tension*rmax0*Mathr::PI ;
- double rmax1= abs(solver->Compute_EffectiveRadius(cell, 1));double D1=surface_tension*rmax1*Mathr::PI ;
- double rmax2= abs(solver->Compute_EffectiveRadius(cell, 2));double D2=surface_tension*rmax2*Mathr::PI ;
- double rmax3= abs(solver->Compute_EffectiveRadius(cell, 3));double D3=surface_tension*rmax3*Mathr::PI ;
- file << cell->info().index << " " <<computeDeltaForce(cell,0,rmax0)<<" "<<computeDeltaMin(cell, 0)<<endl;
- file << cell->info().index << " " <<computeDeltaForce(cell,1,rmax1)<<" "<<computeDeltaMin(cell, 1)<<endl;
- file << cell->info().index << " " <<computeDeltaForce(cell,2,rmax2)<<" "<<computeDeltaMin(cell, 2)<<endl;
- file << cell->info().index << " " <<computeDeltaForce(cell,3,rmax3)<<" "<<computeDeltaMin(cell, 3)<<endl;
+ double rmax0= abs(solver->Compute_EffectiveRadius(cell, 0));double D0=rmax0*Mathr::PI ;
+ double rmax1= abs(solver->Compute_EffectiveRadius(cell, 1));double D1=rmax1*Mathr::PI ;
+ double rmax2= abs(solver->Compute_EffectiveRadius(cell, 2));double D2=rmax2*Mathr::PI ;
+ double rmax3= abs(solver->Compute_EffectiveRadius(cell, 3));double D3=rmax3*Mathr::PI ;
+ file << cell->info().index << " " <<computeDeltaForce(cell,0,rmax0)<<" "<<D0<<endl;
+ file << cell->info().index << " " <<computeDeltaForce(cell,1,rmax1)<<" "<<D1<<endl;
+ file << cell->info().index << " " <<computeDeltaForce(cell,2,rmax2)<<" "<<D2<<endl;
+ file << cell->info().index << " " <<computeDeltaForce(cell,3,rmax3)<<" "<<D3<<endl;
}
file.close();
-}
-template<class Cellhandle>
-double UnsaturatedEngine::computeDeltaMin(Cellhandle cell, int j)
-{
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- if (tri.is_infinite(cell->neighbor(j))) return 0;
-
- Vector3r posA = makeVector3r2(cell->vertex(facetVertices[j][0])->point().point());
- Vector3r posB = makeVector3r2(cell->vertex(facetVertices[j][1])->point().point());
- Vector3r posC = makeVector3r2(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; }
-// double Area_ABC = 0.5 * ((posB-posA).cross(posC-posA)).norm();
-// double Area_SolidA = 0.5*A*pow(rA,2);
-// double Area_SolidB = 0.5*B*pow(rB,2);
-// double Area_SolidC = 0.5*C*pow(rC,2);
-
- double rmax = abs(solver->Compute_EffectiveRadius(cell, j));
- double rmin = max(rAB,max(rBC,rAC));
- if (rmin==0) {/*cerr<<"1";*/
- rmin= 1.0e-10;
- }
- double deltaMin = computeDeltaForce(cell,j,rmin);
- if(deltaMin>0) {
- if((rAB!=0)&&(rBC!=0)&&(rAC!=0)) {
- /* cerr<<"rAB= "<<rAB<<" rBC="<<rBC<<" rAC="<<rAC<<" rmin="<<rmin<<endl;
- */double deltaMax = computeDeltaForce(cell,j,rmax);
-// cerr<<"deltaMax="<<deltaMax<<" deltaMin="<<deltaMin<<" rmax="<<rmax<<" rmin="<<rmin<<" rA= "<<rA<<" rB="<<rB<<" rC="<<rC<<endl;
- }
- }
- return deltaMin;
-}
+}*/
//----------end tempt function for Vahid Joekar-Niasar's data (clear later)---------------------
YADE_PLUGIN ( ( UnsaturatedEngine ) );
