yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11395
[Branch ~yade-pkg/yade/git-trunk] Rev 3390: -clean code
------------------------------------------------------------
revno: 3390
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Mon 2014-01-13 11:57:20 +0100
message:
-clean code
modified:
pkg/dem/UnsaturatedEngine.cpp
pkg/dem/UnsaturatedEngine.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/dem/UnsaturatedEngine.cpp'
--- pkg/dem/UnsaturatedEngine.cpp 2013-11-19 08:56:36 +0000
+++ pkg/dem/UnsaturatedEngine.cpp 2014-01-13 10:57:20 +0000
@@ -37,7 +37,7 @@
const int facetVertices [4][3] = {{1,2,3},{0,2,3},{0,1,3},{0,1,2}};
-Real UnsaturatedEngine::testFunction()
+void UnsaturatedEngine::testFunction()
{
cout<<"This is UnsaturatedEngine test program"<<endl;
RTriangulation& triangulation = solver->T[solver->currentTes].Triangulation();
@@ -188,6 +188,106 @@
}
}
+//initialize the boundingCells info().isWaterReservoir=true, on condition that those cells meet (flow->boundingCells[bound].size()!=0 && cell->info().p()==0)
+template<class Solver>
+void UnsaturatedEngine::initWaterReservoirBound(Solver& flow/*, int boundN*/)
+{
+ RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
+ Finite_cells_iterator cell_end = tri.finite_cells_end();
+ for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+ cell->info().isWaterReservoir = false;
+ }
+ for (int bound=0; bound<6; bound++) {
+ if (flow->boundingCells[bound].size()==0) continue;
+ vector<Cell_handle>::iterator it = flow->boundingCells[bound].begin();
+ for ( it ; it != flow->boundingCells[bound].end(); it++) {
+ if ((*it)->info().index == 0) continue;
+ if((*it)->info().p()!=0) continue;//FIXME: the default pressure for water is 0
+ (*it)->info().isWaterReservoir = true;
+ }
+ }
+}
+
+template<class Solver>
+void UnsaturatedEngine::updateWaterReservoir(Solver& flow)
+{
+ initWaterReservoirBound(flow);
+ for (int bound=0; bound<6; bound++) {
+ if (flow->boundingCells[bound].size()==0) continue;
+ vector<Cell_handle>::iterator it = flow->boundingCells[bound].begin();
+ for ( it ; it != flow->boundingCells[bound].end(); it++) {
+ if ((*it)->info().index == 0) continue;
+ if((*it)->info().isWaterReservoir == false) continue;
+ waterReservoirRecursion((*it),flow);
+ }
+ }
+}
+
+template<class Solver>
+void UnsaturatedEngine::waterReservoirRecursion(Cell_handle cell, Solver& flow)
+{
+ if(cell->info().isWaterReservoir == true) {
+ for (int facet = 0; facet < 4; facet ++) {
+ if (flow->T[flow->currentTes].Triangulation().is_infinite(cell->neighbor(facet))) continue;
+ if (cell->neighbor(facet)->info().p()!=0) continue;
+ if (cell->neighbor(facet)->info().isWaterReservoir==true) continue;
+ Cell_handle n_cell = cell->neighbor(facet);
+ n_cell->info().isWaterReservoir = true;
+ waterReservoirRecursion(n_cell,flow);
+ }
+ }
+}
+
+//initialize the boundingCells info().isAirReservoir=true, on condition that those cells meet (flow->boundingCells[bound].size()!=0 && cell->info().p()!=0)
+template<class Solver>
+void UnsaturatedEngine::initAirReservoirBound(Solver& flow)
+{
+ RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
+ Finite_cells_iterator cell_end = tri.finite_cells_end();
+ for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+ cell->info().isAirReservoir = false;
+ }
+ for (int bound=0; bound<6; bound++) {
+ if (flow->boundingCells[bound].size()==0) continue;
+ vector<Cell_handle>::iterator it = flow->boundingCells[bound].begin();
+ for ( it ; it != flow->boundingCells[bound].end(); it++) {
+ if ((*it)->info().index == 0) continue;
+ if((*it)->info().p() == 0) continue;//FIXME: the default pressure for water is 0
+ (*it)->info().isAirReservoir = true;
+ }
+ }
+}
+
+template<class Solver>
+void UnsaturatedEngine::updateAirReservoir(Solver& flow)
+{
+ initAirReservoirBound(flow);
+ for (int bound=0; bound<6; bound++) {
+ if (flow->boundingCells[bound].size()==0) continue;
+ vector<Cell_handle>::iterator it = flow->boundingCells[bound].begin();
+ for ( it ; it != flow->boundingCells[bound].end(); it++) {
+ if ((*it)->info().index == 0) continue;
+ if((*it)->info().isAirReservoir == false) continue;
+ airReservoirRecursion((*it),flow);
+ }
+ }
+}
+
+template<class Solver>
+void UnsaturatedEngine::airReservoirRecursion(Cell_handle cell, Solver& flow)
+{
+ if(cell->info().isAirReservoir == true) {
+ for (int facet = 0; facet < 4; facet ++) {
+ if (flow->T[flow->currentTes].Triangulation().is_infinite(cell->neighbor(facet))) continue;
+ if (cell->neighbor(facet)->info().p() == 0) continue;
+ if (cell->neighbor(facet)->info().isAirReservoir == true) continue;
+ Cell_handle n_cell = cell->neighbor(facet);
+ n_cell->info().isAirReservoir = true;
+ airReservoirRecursion(n_cell,flow);
+ }
+ }
+}
+
template<class Cellhandle>
double UnsaturatedEngine::computeEffPoreRadius(Cellhandle cell, int j)
{
@@ -234,7 +334,7 @@
return EffPoreRadius;
}
else if(deltaForce_rMin<0) {
- double effPoreRadius = bisection(cell,j,rmin,rmax);// cerr<<"1";
+ double effPoreRadius = bisection(cell,j,rmin,rmax);// cerr<<"1";//we suppose most cases should be this.
// double deltaF=computeDeltaForce(cell,j,effPoreRadius);cerr<<"deltaForce_rMin: "<<deltaForce_rMin<<" deltaForce_rMax: "<<deltaForce_rMax<<" deltaF: "<<deltaF<<endl;//test rMin<R<rMax ?
return effPoreRadius;
}
@@ -245,7 +345,7 @@
}
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;
+ cerr<<"deltaForce_rMin>deltaForce_rMax. cellID: "<<cell->info().index<<"; deltaForce_rMin="<<deltaForce_rMin<<"; deltaForce_rMax="<<deltaForce_rMax<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;
return EffPoreRadius;
}
}
@@ -956,11 +1056,12 @@
vector<Cell_handle>::iterator it = flow->boundingCells[3].begin();
ofstream file;
file.open("ListAdjacentCellsTopBoundary.txt");
- file << "#IDs of boundary pore bodies, which are connecting water reservoir. \n";
+ file << "#Check the boundingCells[3] statement of last step. boundingCells[3] connecting water reservoir(ymax), shoule be isWaterReservoir=0 in the last step. \n";
+ file << "Cell_ID"<<" Cell_Pressure"<<" isWaterReservoir"<<endl;
for ( it ; it != flow->boundingCells[3].end(); it++) {
if ((*it)->info().index == 0) continue;
// cerr << (*it)->info().index << " ";
- file << (*it)->info().index << endl;
+ file << (*it)->info().index <<" "<<(*it)->info().p()<<" "<<(*it)->info().isWaterReservoir<< endl;
}
file.close();
}
@@ -976,115 +1077,16 @@
vector<Cell_handle>::iterator it = flow->boundingCells[2].begin();
ofstream file;
file.open("ListAdjacentCellsBottomBoundary.txt");
- file << "#IDs of boundary pore bodies, which are connecting air reservoir. \n";
+ file << "#Checking the boundingCells[2] statement of last step. boundingCells[2] connecting air reservoir(ymin), should be isAirReservoir=1 in the last step.\n";
+ file << "Cell_ID"<<" Cell_Pressure"<<" isAirReservoir"<<endl;
for ( it ; it != flow->boundingCells[2].end(); it++) {
if ((*it)->info().index == 0) continue;
- file << (*it)->info().index << endl;
+ file << (*it)->info().index <<" "<<(*it)->info().p()<<" "<<(*it)->info().isAirReservoir<<endl;
}
file.close();
}
}
-//initialize the boundingCells info().isWaterReservoir=true, on condition that those cells meet (flow->boundingCells[bound].size()!=0 && cell->info().p()==0)
-template<class Solver>
-void UnsaturatedEngine::initWaterReservoirBound(Solver& flow/*, int boundN*/)
-{
- RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
- Finite_cells_iterator cell_end = tri.finite_cells_end();
- for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
- cell->info().isWaterReservoir = false;
- }
- for (int bound=0; bound<6; bound++) {
- if (flow->boundingCells[bound].size()==0) continue;
- vector<Cell_handle>::iterator it = flow->boundingCells[bound].begin();
- for ( it ; it != flow->boundingCells[bound].end(); it++) {
- if ((*it)->info().index == 0) continue;
- if((*it)->info().p()!=0) continue;//FIXME: the default pressure for water is 0
- (*it)->info().isWaterReservoir = true;
- }
- }
-}
-
-template<class Solver>
-void UnsaturatedEngine::updateWaterReservoir(Solver& flow)
-{
- initWaterReservoirBound(flow);
- for (int bound=0; bound<6; bound++) {
- if (flow->boundingCells[bound].size()==0) continue;
- vector<Cell_handle>::iterator it = flow->boundingCells[bound].begin();
- for ( it ; it != flow->boundingCells[bound].end(); it++) {
- if ((*it)->info().index == 0) continue;
- if((*it)->info().isWaterReservoir == false) continue;
- waterReservoirRecursion((*it),flow);
- }
- }
-}
-
-template<class Solver>
-void UnsaturatedEngine::waterReservoirRecursion(Cell_handle cell, Solver& flow)
-{
- if(cell->info().isWaterReservoir == true) {
- for (int facet = 0; facet < 4; facet ++) {
- if (flow->T[flow->currentTes].Triangulation().is_infinite(cell->neighbor(facet))) continue;
- if (cell->neighbor(facet)->info().p()!=0) continue;
- if (cell->neighbor(facet)->info().isWaterReservoir==true) continue;
- Cell_handle n_cell = cell->neighbor(facet);
- n_cell->info().isWaterReservoir = true;
- waterReservoirRecursion(n_cell,flow);
- }
- }
-}
-
-//initialize the boundingCells info().isAirReservoir=true, on condition that those cells meet (flow->boundingCells[bound].size()!=0 && cell->info().p()!=0)
-template<class Solver>
-void UnsaturatedEngine::initAirReservoirBound(Solver& flow)
-{
- RTriangulation& tri = flow->T[flow->currentTes].Triangulation();
- Finite_cells_iterator cell_end = tri.finite_cells_end();
- for ( Finite_cells_iterator cell = tri.finite_cells_begin(); cell != cell_end; cell++ ) {
- cell->info().isAirReservoir = false;
- }
- for (int bound=0; bound<6; bound++) {
- if (flow->boundingCells[bound].size()==0) continue;
- vector<Cell_handle>::iterator it = flow->boundingCells[bound].begin();
- for ( it ; it != flow->boundingCells[bound].end(); it++) {
- if ((*it)->info().index == 0) continue;
- if((*it)->info().p() == 0) continue;//FIXME: the default pressure for water is 0
- (*it)->info().isAirReservoir = true;
- }
- }
-}
-
-template<class Solver>
-void UnsaturatedEngine::updateAirReservoir(Solver& flow)
-{
- initAirReservoirBound(flow);
- for (int bound=0; bound<6; bound++) {
- if (flow->boundingCells[bound].size()==0) continue;
- vector<Cell_handle>::iterator it = flow->boundingCells[bound].begin();
- for ( it ; it != flow->boundingCells[bound].end(); it++) {
- if ((*it)->info().index == 0) continue;
- if((*it)->info().isAirReservoir == false) continue;
- airReservoirRecursion((*it),flow);
- }
- }
-}
-
-template<class Solver>
-void UnsaturatedEngine::airReservoirRecursion(Cell_handle cell, Solver& flow)
-{
- if(cell->info().isAirReservoir == true) {
- for (int facet = 0; facet < 4; facet ++) {
- if (flow->T[flow->currentTes].Triangulation().is_infinite(cell->neighbor(facet))) continue;
- if (cell->neighbor(facet)->info().p() == 0) continue;
- if (cell->neighbor(facet)->info().isAirReservoir == true) continue;
- Cell_handle n_cell = cell->neighbor(facet);
- n_cell->info().isAirReservoir = true;
- airReservoirRecursion(n_cell,flow);
- }
- }
-}
-
template<class Solver>
void UnsaturatedEngine::setImposedPressure ( unsigned int cond, Real p,Solver& flow )
{
@@ -1204,7 +1206,8 @@
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);}
+ if(areaPore<0) {cerr<<"areaPore Negative. cellID: "<<cell->info().index<<". posA="<<posA<<"; posB="<<posB<<"; posC="<< posC<<"; rA="<< rA<<"; rB="<< rB<<"; rC="<<rC<<endl;
+ areaPore=Mathr::PI*pow(rInscribe,2);}
return areaPore;
}; break;
case (1) : { return Mathr::PI*pow(rInscribe,2); }; break;
=== modified file 'pkg/dem/UnsaturatedEngine.hpp'
--- pkg/dem/UnsaturatedEngine.hpp 2013-11-19 08:58:13 +0000
+++ pkg/dem/UnsaturatedEngine.hpp 2014-01-13 10:57:20 +0000
@@ -41,7 +41,7 @@
vector<posData> positionBufferParallel;//keep the positions from a given step for multithread factorization
// //copy positions in a buffer for faster and/or parallel access
void setPositionsBuffer(bool current);
- Real testFunction();
+ void testFunction();
public :
// enum {wall_left=0, wall_right, wall_bottom, wall_top, wall_back, wall_front};
