yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 4008
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Tue 2017-03-07 11:48:59 +0100
message:
  add functions to get pore throat radius (by cells or spheres).
modified:
  pkg/pfv/TwoPhaseFlowEngine.cpp
  pkg/pfv/TwoPhaseFlowEngine.hpp


--
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/TwoPhaseFlowEngine.cpp'
--- pkg/pfv/TwoPhaseFlowEngine.cpp	2017-01-25 16:17:13 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.cpp	2017-03-07 10:48:59 +0000
@@ -1122,6 +1122,163 @@
 	}
 }
 
+bool TwoPhaseFlowEngine::isCellNeighbor(unsigned int cell1, unsigned int cell2)
+{
+  bool neighbor=false;
+  for (unsigned int i=0;i<4;i++) {
+    if (solver->T[solver->currentTes].cellHandles[cell1]->neighbor(i)->info().id==cell2)
+    {neighbor=true;break;}
+  }
+  return neighbor;
+}
+
+void TwoPhaseFlowEngine::setPoreThroatRadius(unsigned int cell1, unsigned int cell2, double radius)
+{
+    if (isCellNeighbor(cell1,cell2)==false) {
+        cout<<"cell1 and cell2 are not neighbors."<<endl;}
+    else {
+        for (unsigned int i=0; i<4; i++) {
+            if (solver->T[solver->currentTes].cellHandles[cell1]->neighbor(i)->info().id==cell2)
+                solver->T[solver->currentTes].cellHandles[cell1]->info().poreThroatRadius[i]=radius;}}
+}
+double TwoPhaseFlowEngine::getPoreThroatRadius(unsigned int cell1, unsigned int cell2)
+{
+    double r =-1.;
+    if (isCellNeighbor(cell1,cell2)==false) {
+        cout<<"cell1 and cell2 are not neighbors."<<endl;}
+    else {
+        for (unsigned int i=0; i<4; i++) {
+            if (solver->T[solver->currentTes].cellHandles[cell1]->neighbor(i)->info().id==cell2)
+                r = solver->T[solver->currentTes].cellHandles[cell1]->info().poreThroatRadius[i];}}
+    return r;
+}
+
+////compute entry pc independently
+double TwoPhaseFlowEngine::computeEffRcByPosRadius(const Vector3r& posA, const double& rA, const Vector3r& posB, const double& rB, const Vector3r& posC, const double& rC)
+{
+	CVector B = makeCgVect(posB - posA);
+	CVector x = B/sqrt(B.squared_length());
+	CVector C = makeCgVect(posC - posA);
+	CVector z = CGAL::cross_product(x,C);
+	CVector y = CGAL::cross_product(x,z);
+	y = y/sqrt(y.squared_length());
+
+	double b1[2]; b1[0] = B*x; b1[1] = B*y;
+	double c1[2]; c1[0] = C*x; c1[1] = C*y;
+
+	double A = ((pow(rA,2))*(1-c1[0]/b1[0])+((pow(rB,2)*c1[0])/b1[0])-pow(rC,2)+pow(c1[0],2)+pow(c1[1],2)-((pow(b1[0],2)+pow(b1[1],2))*c1[0]/b1[0]))/(2*c1[1]-2*b1[1]*c1[0]/b1[0]);
+	double BB = (rA-rC-((rA-rB)*c1[0]/b1[0]))/(c1[1]-b1[1]*c1[0]/b1[0]);
+	double CC = (pow(rA,2)-pow(rB,2)+pow(b1[0],2)+pow(b1[1],2))/(2*b1[0]);
+	double D = (rA-rB)/b1[0];
+	double E = b1[1]/b1[0];
+	double F = pow(CC,2)+pow(E,2)*pow(A,2)-2*CC*E*A;
+
+	double c = -F-pow(A,2)+pow(rA,2);
+	double b = 2*rA-2*(D-BB*E)*(CC-E*A)-2*A*BB;
+	double a = 1-pow((D-BB*E),2)-pow(BB,2);
+
+	if ((pow(b,2)-4*a*c)<0){cout << "NEGATIVE DETERMINANT" << endl; }
+	double reff = (-b+sqrt(pow(b,2)-4*a*c))/(2*a);
+	if (reff<0) return 0;//happens very rarely, with bounding spheres most probably
+	//if the facet involves one ore more bounding sphere, we return R with a minus sign
+// 	if (cell->vertex(facetVertices[j][2])->info().isFictious || cell->vertex(facetVertices[j][1])->info().isFictious || cell->vertex(facetVertices[j][2])->info().isFictious) return -reff;
+	else return reff;
+}
+double TwoPhaseFlowEngine::computeMSPRcByPosRadius(const Vector3r& posA, const double& rA, const Vector3r& posB, const double& rB, const Vector3r& posC, const double& rC)
+{
+//     Vector3r pos[3]; //solid pos
+//     double r[3]; //solid radius
+    double e[3]; //edges of triangulation
+    double g[3]; //gap radius between solid
+        
+    e[0] = (posB-posC).norm();
+    e[1] = (posC-posA).norm();
+    e[2] = (posB-posA).norm();
+    g[0] = ((e[0]-rB-rC)>0) ? 0.5*(e[0]-rB-rC):0 ;
+    g[1] = ((e[1]-rC-rA)>0) ? 0.5*(e[1]-rC-rA):0 ;
+    g[2] = ((e[2]-rA-rB)>0) ? 0.5*(e[2]-rA-rB):0 ;
+    
+    double rmin= (std::max(g[0],std::max(g[1],g[2]))==0) ? 1.0e-10:std::max(g[0],std::max(g[1],g[2])) ;
+    double rmax = std::abs(computeEffRcByPosRadius(posA, rA, posB, rB, posC, rC));
+//     if(rmin>rmax) { cerr<<"WARNING! rmin>rmax. rmin="<<rmin<<" ,rmax="<<rmax<<endl; }
+    
+    double deltaForceRMin = computeDeltaForceTmp(posA, rA, posB, rB, posC, rC, rmin);
+    double deltaForceRMax = computeDeltaForceTmp(posA, rA, posB, rB, posC, rC, rmax);
+    double effPoreRadius;
+    
+    if(deltaForceRMin>deltaForceRMax) { effPoreRadius=rmax; }
+    else if(deltaForceRMax<0) { effPoreRadius=rmax; }
+    else if(deltaForceRMin>0) { effPoreRadius=rmin; }
+    else { effPoreRadius=bisectionTmp(posA, rA, posB, rB, posC, rC, rmin,rmax); }
+    return effPoreRadius;
+}
+double TwoPhaseFlowEngine::bisectionTmp(const Vector3r& posA, const double& rA, const Vector3r& posB, const double& rB, const Vector3r& posC, const double& rC, double a, double b)
+{
+    double m = 0.5*(a+b);
+    if (std::abs(b-a)>std::abs((computeEffRcByPosRadius(posA, rA, posB, rB, posC, rC)*1.0e-6))) {
+        if ( computeDeltaForceTmp(posA, rA, posB, rB, posC, rC,m) * computeDeltaForceTmp(posA, rA, posB, rB, posC, rC,a) < 0 ) {
+            b = m;
+//             if(solve->debugOut) {cerr<<"deltaF(m)="<<computeDeltaForceTmp(posA, rA, posB, rB, posC, rC,m)<<",m="<<m<<endl;}
+            return bisectionTmp(posA, rA, posB, rB, posC, rC, a,b);}
+        else {
+            a = m;
+//             if(solve->debugOut) {cerr<<"deltaF(m)="<<computeDeltaForceTmp(posA, rA, posB, rB, posC, rC,m)<<",m="<<m<<endl;}
+            return bisectionTmp(posA, rA, posB, rB, posC, rC, a,b);}
+    }
+    else {
+//         if(solve->debugOut) {cerr<<"deltaF(m)="<<computeDeltaForceTmp(posA, rA, posB, rB, posC, rC,m)<<",m="<<m<<endl;}
+        return m;}
+}
+double TwoPhaseFlowEngine::computeDeltaForceTmp(const Vector3r& posA, const double& rA, const Vector3r& posB, const double& rB, const Vector3r& posC, const double& rC, double r)
+{
+//     Vector3r pos[3]; //solid pos
+//     double r[3]; //solid radius
+    double rRc[3]; //r[i] + r (r: capillary radius)
+    double e[3]; //edges of triangulation
+    double rad[4][3]; //angle in radian
+    
+    rRc[0] = rA+r;
+    rRc[1] = rB+r;
+    rRc[2] = rC+r;
+        
+    e[0] = (posB-posC).norm();
+    e[1] = (posC-posA).norm();
+    e[2] = (posB-posA).norm();
+    
+    rad[3][0]=acos(((posB-posA).dot(posC-posA))/(e[2]*e[1]));
+    rad[3][1]=acos(((posC-posB).dot(posA-posB))/(e[0]*e[2]));
+    rad[3][2]=acos(((posA-posC).dot(posB-posC))/(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])*r;
+    double lNS = (rad[3][0]-rad[1][0]-rad[2][0])*rA + (rad[3][1]-rad[2][1]-rad[0][1])*rB + (rad[3][2]-rad[1][2]-rad[0][2])*rC ;
+    double lInterface=lNW+lNS;
+    
+    double sW0=0.5*rRc[1]*rRc[2]*sin(rad[0][0])-0.5*rad[0][0]*pow(r,2)-0.5*rad[0][1]*pow(rB,2)-0.5*rad[0][2]*pow(rC,2) ;
+    double sW1=0.5*rRc[2]*rRc[0]*sin(rad[1][1])-0.5*rad[1][1]*pow(r,2)-0.5*rad[1][2]*pow(rC,2)-0.5*rad[1][0]*pow(rA,2) ;
+    double sW2=0.5*rRc[0]*rRc[1]*sin(rad[2][2])-0.5*rad[2][2]*pow(r,2)-0.5*rad[2][0]*pow(rA,2)-0.5*rad[2][1]*pow(rB,2) ;
+    double sW=sW0+sW1+sW2;
+    
+    
+    CVector facetSurface = 0.5*CGAL::cross_product(makeCgVect(posA-posC),makeCgVect(posB-posC));
+    double sVoid = sqrt(facetSurface.squared_length()) - (0.5*rad[3][0]*pow(rA,2) + 0.5*rad[3][1]*pow(rB,2) + 0.5*rad[3][2]*pow(rC,2));
+    double sInterface=sVoid-sW;
+
+    double deltaF = lInterface - sInterface/r;//deltaF=surfaceTension*(perimeterPore - areaPore/rCap)
+    return deltaF;
+}
+
 
 #endif //TwoPhaseFLOW
  

=== modified file 'pkg/pfv/TwoPhaseFlowEngine.hpp'
--- pkg/pfv/TwoPhaseFlowEngine.hpp	2017-01-25 18:08:06 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.hpp	2017-03-07 10:48:59 +0000
@@ -201,6 +201,17 @@
 	//FIXME, I removed indexing cells from inside UnsatEngine (SoluteEngine shouldl be ok (?)) in order to get pressure computed, problem is they are not indexed at all if flow is not calculated
 	void computeOnePhaseFlow() {scene = Omega::instance().getScene().get(); if (!solver) cerr<<"no solver!"<<endl; solver->gaussSeidel(scene->dt);initSolver(*solver);}
 	
+	///manipulate/get/set on pore geometry
+	bool isCellNeighbor(unsigned int cell1, unsigned int cell2);
+	void setPoreThroatRadius(unsigned int cell1, unsigned int cell2, double radius);
+	double getPoreThroatRadius(unsigned int cell1, unsigned int cell2);
+
+	double computeEffRcByPosRadius(const Vector3r& posA, const double& rA, const Vector3r& posB, const double& rB, const Vector3r& posC, const double& rC);
+	double computeMSPRcByPosRadius(const Vector3r& posA, const double& rA, const Vector3r& posB, const double& rB, const Vector3r& posC, const double& rC);
+	double bisectionTmp(const Vector3r& posA, const double& rA, const Vector3r& posB, const double& rB, const Vector3r& posC, const double& rC, double a, double b);
+	double computeDeltaForceTmp(const Vector3r& posA, const double& rA, const Vector3r& posB, const double& rB, const Vector3r& posC, const double& rC, double r);
+
+	
 	CELL_SCALAR_GETTER(bool,.isWRes,cellIsWRes)
 	CELL_SCALAR_GETTER(bool,.isNWRes,cellIsNWRes)
 	CELL_SCALAR_GETTER(bool,.isTrapW,cellIsTrapW)
@@ -215,6 +226,7 @@
 	CELL_SCALAR_GETTER(Real,.porosity,cellPorosity)
 	CELL_SCALAR_SETTER(bool,.hasInterface,setCellHasInterface) //Temporary function to allow for simulations in Python
 	CELL_SCALAR_GETTER(int,.label,cellLabel)
+	CELL_SCALAR_GETTER(Real,.volume(),cellVolume) //Temporary function to allow for simulations in Python	
 
 	YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(TwoPhaseFlowEngine,TwoPhaseFlowEngineT,"documentation here",
 	((double,surfaceTension,0.0728,,"Water Surface Tension in contact with air at 20 Degrees Celsius is: 0.0728(N/m)"))
@@ -247,7 +259,7 @@
 	.def("computeOnePhaseFlow",&TwoPhaseFlowEngine::computeOnePhaseFlow,"compute pressure and fluxes in the W-phase")
 	.def("initialization",&TwoPhaseFlowEngine::initialization,"Initialize invasion setup. Build network, compute pore geometry info and initialize reservoir boundary conditions. ")
 	.def("computePoreSatAtInterface",&TwoPhaseFlowEngine::computePoreSatAtInterface,(boost::python::arg("ID")),"compute pressure and fluxes in the W-phase")
-	.def("getPoreThroatRadius",&TwoPhaseFlowEngine::cellporeThroatRadius,"get 4 pore throat radii")
+	.def("getPoreThroatRadiusList",&TwoPhaseFlowEngine::cellporeThroatRadius,(boost::python::arg("cell_ID")),"get 4 pore throat radii of a cell.")
 	.def("getNeighbors",&TwoPhaseFlowEngine::getNeighbors,"get 4 neigboring cells")
 	.def("getCellHasInterface",&TwoPhaseFlowEngine::cellHasInterface,"indicates whether a NW-W interface is present within the cell")
 	.def("getCellInSphereRadius",&TwoPhaseFlowEngine::cellInSphereRadius,"get the radius of the inscribed sphere in a pore unit")
@@ -264,6 +276,12 @@
 	.def("saveVtk",&TwoPhaseFlowEngine::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
 	.def("getPotentialPendularSpheresPair",&TwoPhaseFlowEngine::getPotentialPendularSpheresPair,"Get the list of sphere ID pairs of potential pendular liquid bridge.")
 	.def("getClusters",&TwoPhaseFlowEngine::pyClusters/*,(boost::python::arg("folder")="./VTK")*/,"Get the list of clusters.")
+	.def("getCellVolume",&TwoPhaseFlowEngine::cellVolume,"get the volume of each cell")
+	.def("isCellNeighbor",&TwoPhaseFlowEngine::isCellNeighbor,(boost::python::arg("cell1_ID"), boost::python::arg("cell2_ID")),"check if cell1 and cell2 are neigbors.")
+	.def("setPoreThroatRadius",&TwoPhaseFlowEngine::setPoreThroatRadius, (boost::python::arg("cell1_ID"), boost::python::arg("cell2_ID"), boost::python::arg("radius")), "set the pore throat radius between cell1 and cell2.")
+	.def("getPoreThroatRadius",&TwoPhaseFlowEngine::getPoreThroatRadius, (boost::python::arg("cell1_ID"), boost::python::arg("cell2_ID")), "get the pore throat radius between cell1 and cell2.")
+	.def("getEffRcByPosRadius",&TwoPhaseFlowEngine::computeEffRcByPosRadius, (boost::python::arg("position1"),boost::python::arg("radius1"),boost::python::arg("position2"),boost::python::arg("radius2"),boost::python::arg("position3"),boost::python::arg("radius3")), "get effective radius by three spheres position and radius.(inscribed sphere)")
+	.def("getMSPRcByPosRadius",&TwoPhaseFlowEngine::computeMSPRcByPosRadius, (boost::python::arg("position1"),boost::python::arg("radius1"),boost::python::arg("position2"),boost::python::arg("radius2"),boost::python::arg("position3"),boost::python::arg("radius3")), "get entry radius wrt MSP method by three spheres position and radius.")
 	
 	)
 	DECLARE_LOGGER;