yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06293
[Branch ~yade-dev/yade/trunk] Rev 2572: - new function for adding bounding planes
------------------------------------------------------------
revno: 2572
committer: Emanuele Catalano <catalano@xxxxxxxxxxx
branch nick: yade
timestamp: Thu 2010-11-25 15:22:43 +0100
message:
- new function for adding bounding planes
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/Network.cpp
lib/triangulation/Network.hpp
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.hpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'lib/triangulation/FlowBoundingSphere.cpp'
--- lib/triangulation/FlowBoundingSphere.cpp 2010-11-19 17:18:56 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-11-25 14:22:43 +0000
@@ -180,13 +180,11 @@
//GenerateVoxelFile(); ///*GENERATION OF A VOXEL FILE *///
/** INITIALIZATION OF VOLUMES AND PRESSURES **/
- Real totV=0;
Finite_cells_iterator cell_end = Tri.finite_cells_end();
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
- cell->info().volume() = Volume_Pore(cell); totV+=cell->info().volume();
+ cell->info().volume() = ( std::abs ( ( CGT::Tetraedre ( cell->vertex(0)->point(),cell->vertex(1)->point(),cell->vertex(2)->point(),cell->vertex(3)->point()).volume() ) ) );
cell->info().dv() = 0;
}
- if(DEBUG_OUT) cerr << "TOTAL VOID VOLUME: "<<totV<<endl;
clock.top("initializing delta_volumes");
@@ -196,7 +194,9 @@
Compute_Permeability();
clock.top("Compute_Permeability");
/** END PERMEABILITY CALCULATION**/
-
+
+ if(DEBUG_OUT) cerr << "TOTAL VOID VOLUME: " << Vporale <<endl;
+
/** STATISTICS **/
DisplayStatistics();
clock.top("DisplayStatistics");
@@ -627,7 +627,9 @@
if(!Tri.is_infinite(*it)){
Point& p1 = (*it)->info();
Cell_handle& cell = *it;
- if (p1.x()>(alpha*(x_max-x_min)) && p1.x()<((alpha+step)*(x_max-x_min))){cell->info().p() = (Pressure/2)*(1+cos(alpha*M_PI));}
+ if (p1.x()<0) cell->info().p() = Pressure;
+ else if (p1.x()>x_max) cell->info().p() = 0.f;
+ else if (p1.x()>(alpha*(x_max-x_min)) && p1.x()<((alpha+step)*(x_max-x_min))) cell->info().p() = (Pressure/2)*(1+cos(alpha*M_PI));
}
}
}
@@ -688,12 +690,11 @@
int NEG=0, POS=0, pass=0;
bool ref = Tri.finite_cells_begin()->info().isvisited;
-
Vecteur n;
// std::ofstream oFile( "Radii",std::ios::out);
// std::ofstream fFile( "Radii_Fictious",std::ios::out);
// std::ofstream kFile ( "LocalPermeabilities" ,std::ios::app );
- Real meanK=0, STDEV=0;
+ Real meanK=0, STDEV=0, meanRadius=0, meanDistance=0;
Real infiniteK=1e10;
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
@@ -741,7 +742,8 @@
(cell->info().k_norm())[j]= k*k_factor;
// (cell->info().facetSurf())[j]= k*n;
(neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]= k*k_factor;
-
+ meanDistance += distance;
+ meanRadius += radius;
meanK += (cell->info().k_norm())[j];
// if (!meanK_LIMIT) kFile << ( cell->info().k_norm() )[j] << endl;
@@ -752,11 +754,15 @@
cell->info().isvisited = !ref;
}
meanK /= pass;
+ meanRadius /= pass;
+ meanDistance /= pass;
double maxKdivKmean=MAXK_DIV_KMEAN;
if (DEBUG_OUT) {
cout << "PassCompK = " << pass << endl;
cout << "meanK = " << meanK << endl;
cout << "maxKdivKmean = " << maxKdivKmean << endl;
+ cout << "meanTubesRadius = " << meanRadius << endl;
+ cout << "meanDistance = " << meanDistance << endl;
}
ref = Tri.finite_cells_begin()->info().isvisited;
pass=0;
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2010-11-19 17:18:56 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2010-11-25 14:22:43 +0000
@@ -158,6 +158,7 @@
// double volume_double_fictious_pore (Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3, Point& PV1, Point& PV2, Vecteur& facetSurface);
// double PoreVolume (RTriangulation& Tri, Cell_handle cell);
+// void Slice_Average_Cell_Velocity();
int Average_Cell_Velocity(int id_sphere, RTriangulation& Tri);
void Average_Cell_Velocity();
void Average_Grain_Velocity();
=== modified file 'lib/triangulation/Network.cpp'
--- lib/triangulation/Network.cpp 2010-11-17 20:07:21 +0000
+++ lib/triangulation/Network.cpp 2010-11-25 14:22:43 +0000
@@ -31,30 +31,13 @@
const double FAR = 500;
using namespace std;
-using namespace boost;
+// using namespace boost;
namespace CGT {
Network::~Network(){}
Network::Network(){}
-double Network::Volume_Pore ( Cell_handle cell)
-{
- RTriangulation& Tri = T[currentTes].Triangulation();
- Real v;
- switch (cell->info().fictious()) {
- case(0) :
- v=std::abs(Tri.tetrahedron(cell).volume());
- for (int k=0;k<4;k++) v-= spherical_triangle_volume(cell->vertex(permut4[k][0])->point(),
- cell->vertex(permut4[k][1])->point(),
- cell->vertex(permut4[k][2])->point(),
- cell->vertex(permut4[k][3])->point());
- return v;
- case(1) : return 0;
- case(2) : return 0;
- case(3) : return 0;}
-}
-
int Network::Detect_facet_fictious_vertices (Cell_handle& cell, int& j)
{
fictious_vertex = 0;
@@ -104,6 +87,8 @@
Vsolid_tot += Vsolid1 + Vsolid2;
Vporale += Vtot - (Vsolid1 + Vsolid2);
+ V_porale_porosity += Vtot - (Vsolid1 + Vsolid2);
+ V_totale_porosity += Vtot;
/**Vpore**/ return Vtot - (Vsolid1 + Vsolid2);
}; break;
@@ -116,20 +101,20 @@
double A [3], B[3];
Boundary &bi1 = boundaries [SV1->info().id()];
-
+
for (int m=0;m<3;m++) {A[m]= (SV2->point())[m];}
for (int m=0;m<3;m++) {B[m]= (SV3->point())[m];}
-
+
A[bi1.coordinate]=bi1.p[bi1.coordinate];
B[bi1.coordinate]=bi1.p[bi1.coordinate];
-
+
Point AA(A[0],A[1],A[2]);
Point BB(B[0],B[1],B[2]);
facetSurface = surface_single_fictious_facet(SV1,SV2,SV3);
if (facetSurface*(PV2-PV1) > 0) facetSurface = -1.0*facetSurface;
Real Vtot=ONE_THIRD*abs(facetSurface*(PV1-PV2));
Vtotalissimo += Vtot;
-
+
Sphere A1(AA, 0);
Sphere B1(BB, 0);
Sphere& SW2 = SV2->point();
@@ -137,7 +122,7 @@
Real Vsolid1 = spherical_triangle_volume(SW2, AA, PV1, PV2)+spherical_triangle_volume(SW2, SW3, PV1, PV2);
Real Vsolid2 = spherical_triangle_volume(SW3, BB, PV1, PV2)+spherical_triangle_volume(SW3, SW2, PV1, PV2);
-
+
Vsolid_tot += Vsolid1 + Vsolid2;
Vporale += Vtot - (Vsolid1 + Vsolid2);
@@ -431,7 +416,8 @@
void Network::AddBoundingPlanes()
{
Tesselation& Tes = T[currentTes];
-
+
+ //FIXME: Id's order in boundsIds is done according to the enumerotation of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
y_min_id = Tes.Max_id() + 1;
boundsIds[0]=&y_min_id;
y_max_id = Tes.Max_id() + 2;
@@ -441,74 +427,71 @@
x_max_id = Tes.Max_id() + 4;
boundsIds[3]=&x_max_id;
z_min_id = Tes.Max_id() + 5;
- boundsIds[4]=&z_min_id;
+ boundsIds[4]=&z_max_id;
z_max_id = Tes.Max_id() + 6;
- boundsIds[5]=&z_max_id;
-
+ boundsIds[5]=&z_min_id;
+
+ Corner_min = Point(x_min, y_min, z_min);
+ Corner_max = Point(x_max, y_max, z_max);
+
id_offset = Tes.Max_id() +1;//so that boundaries[vertex->id - offset] gives the ordered boundaries (also see function Boundary& boundary(int b))
-
- AddBoundingPlanes(true);
-}
-
-void Network::AddBoundingPlanes(bool yade)
-{
- Tesselation& Tes = T[currentTes];
- Corner_min = Point(x_min, y_min, z_min);
- Corner_max = Point(x_max, y_max, z_max);
-
- boundsIds[0]= &y_min_id;
- boundsIds[1]= &y_max_id;
- boundsIds[2]= &x_min_id;
- boundsIds[3]= &x_max_id;
- boundsIds[4]= &z_min_id;
- boundsIds[5]= &z_max_id;
-
- Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), Corner_min.y()-FAR*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.z()+Corner_min.z()), FAR*(Corner_max.y()-Corner_min.y()), y_min_id, true);
- boundaries[y_min_id-id_offset].p = Corner_min;
- boundaries[y_min_id-id_offset].normal = Vecteur(0,1,0);
- boundaries[y_min_id-id_offset].coordinate = 1;
-
- if(DEBUG_OUT) cout << "Bottom boundary has been created. ID = " << y_min_id << endl;
-
- Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), Corner_max.y() +FAR*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.z()+Corner_min.z()), FAR*(Corner_max.y()-Corner_min.y()), y_max_id, true);
- boundaries[y_max_id-id_offset].p = Corner_max;
- boundaries[y_max_id-id_offset].normal = Vecteur(0,-1,0);
- boundaries[y_max_id-id_offset].coordinate = 1;
-
- if(DEBUG_OUT) cout << "Top boundary has been created. ID = " << y_max_id << endl;
-
- Tes.insert(Corner_min.x()-FAR*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.y()+Corner_min.y()), 0.5*(Corner_max.z()+Corner_min.z()), FAR*(Corner_max.y()-Corner_min.y()), x_min_id, true);
- boundaries[x_min_id-id_offset].p = Corner_min;
- boundaries[x_min_id-id_offset].normal = Vecteur(1,0,0);
- boundaries[x_min_id-id_offset].coordinate = 0;
-
- if(DEBUG_OUT) cout << "Left boundary has been created. ID = " << x_min_id << endl;
-
- Tes.insert(Corner_max.x() +FAR*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.y()+Corner_min.y()), 0.5*(Corner_max.z()+Corner_min.z()), FAR*(Corner_max.y()-Corner_min.y()), x_max_id, true);
- boundaries[x_max_id-id_offset].p = Corner_max;
- boundaries[x_max_id-id_offset].normal = Vecteur(-1,0,0);
- boundaries[x_max_id-id_offset].coordinate = 0;
-
- if(DEBUG_OUT) cout << "Right boundary has been created. ID = " << x_max_id << endl;
-
- Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), 0.5*(Corner_max.y()+Corner_min.y()), Corner_min.z()-FAR*(Corner_max.y()-Corner_min.y()), FAR*(Corner_max.y()-Corner_min.y()), z_min_id, true);
- boundaries[z_min_id-id_offset].p = Corner_min;
- boundaries[z_min_id-id_offset].normal = Vecteur(0,0,1);
- boundaries[z_min_id-id_offset].coordinate = 2;
-
- if(DEBUG_OUT) cout << "Front boundary has been created. ID = " << z_min_id << endl;
-
- Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), 0.5*(Corner_max.y()+Corner_min.y()), Corner_max.z() +FAR*(Corner_max.y()-Corner_min.y()), FAR*(Corner_max.y()-Corner_min.y()), z_max_id, true);
- boundaries[z_max_id-id_offset].p = Corner_max;
- boundaries[z_max_id-id_offset].normal = Vecteur(0,0,-1);
- boundaries[z_max_id-id_offset].coordinate = 2;
-
- if(DEBUG_OUT) cout << "Back boundary has been created. ID = " << z_max_id << endl;
-
- for (int k=0;k<6;k++) {
- boundaries[k].flowCondition = 1;
- boundaries[k].value = 0;
- }
+
+ AddBoundingPlane (true, Vecteur(0,1,0) , y_min_id);
+ AddBoundingPlane (true, Vecteur(0,-1,0) , y_max_id);
+ AddBoundingPlane (true, Vecteur(-1,0,0) , x_max_id);
+ AddBoundingPlane (true, Vecteur(1,0,0) , x_min_id);
+ AddBoundingPlane (true, Vecteur(0,0,1) , z_min_id);
+ AddBoundingPlane (true, Vecteur(0,0,-1) , z_max_id);
+
+// AddBoundingPlanes(true);
+}
+
+void Network::AddBoundingPlane (bool yade, Vecteur Normal, int id_wall)
+{
+ Tesselation& Tes = T[currentTes];
+
+ int Coordinate = abs(Normal[0])*0 + abs(Normal[1])*1 + abs(Normal[2])*2;
+
+ double pivot = Normal[Coordinate]<0 ?
+ Corner_max.x()*abs(Normal[0])+Corner_max.y()*abs(Normal[1])+Corner_max.z()*abs(Normal[2]) : Corner_min.x()*abs(Normal[0])+Corner_min.y()*abs(Normal[1])+Corner_min.z()*abs(Normal[2]);
+
+ Tes.insert(0.5*(Corner_min.x() +Corner_max.x())*(1-abs(Normal[0]))+(pivot-Normal[0]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]),
+ 0.5*(Corner_max.y() +Corner_min.y())*(1-abs(Normal[1]))+(pivot-Normal[1]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]),
+ 0.5*(Corner_max.z() +Corner_min.z())*(1-abs(Normal[2]))+(pivot-Normal[2]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]),
+ FAR*(Corner_max.y()-Corner_min.y()), id_wall, true);
+
+ if (Normal[Coordinate]<0) boundaries[id_wall-id_offset].p = Corner_max;
+ else boundaries[id_wall-id_offset].p = Corner_min;
+
+ boundaries[id_wall-id_offset].normal = Normal;
+ boundaries[id_wall-id_offset].coordinate = Coordinate;
+
+ boundaries[id_wall-id_offset].flowCondition = 1;
+ boundaries[id_wall-id_offset].value = 0;
+
+ if(DEBUG_OUT) cout << "A boundary -max/min-has been created. ID = " << id_wall << " position = " << 0.5*(Corner_min.x() +Corner_max.x())*(1-abs(Normal[0]))+(pivot-Normal[0]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]) << " , " << 0.5*(Corner_max.y() +Corner_min.y())*(1-abs(Normal[1]))+(pivot-Normal[1]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]) << " , " << 0.5*(Corner_max.z() +Corner_min.z())*(1-abs(Normal[2]))+(pivot-Normal[2]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]) << ". Radius = " << FAR*(Corner_max.y()-Corner_min.y()) << endl;
+}
+
+void Network::AddBoundingPlane (Real center[3], double thickness, Vecteur Normal, int id_wall )
+{
+ Tesselation& Tes = T[currentTes];
+
+ int Coordinate = abs(Normal[0])*0 + abs(Normal[1])*1 + abs(Normal[2])*2;
+
+ Tes.insert((center[0]+Normal[0]*thickness/2)*(1-abs(Normal[0])) + (center[0]+Normal[0]*thickness/2-Normal[0]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]),
+ (center[1]+Normal[1]*thickness/2)*(1-abs(Normal[1])) + (center[1]+Normal[1]*thickness/2-Normal[1]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]),
+ (center[2]+Normal[2]*thickness/2)*(1-abs(Normal[2])) + (center[2]+Normal[2]*thickness/2-Normal[2]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]),
+ FAR*(Corner_max.y()-Corner_min.y()), id_wall, true);
+
+ Point P (center[0],center[1],center[2]);
+ boundaries[id_wall-id_offset].p = P;
+ boundaries[id_wall-id_offset].normal = Normal;
+ boundaries[id_wall-id_offset].coordinate = Coordinate;
+
+ boundaries[id_wall-id_offset].flowCondition = 1;
+ boundaries[id_wall-id_offset].value = 0;
+
+ if(DEBUG_OUT) cout << "A boundary -center/thick- has been created. ID = " << id_wall << " position = " << center[0]+Normal[0]*thickness/2 + (center[0]+Normal[0]*thickness/2-Normal[0]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]) << " , " << center[1]+Normal[0]*thickness/2 + (center[1]+Normal[1]*thickness/2-Normal[1]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]) << " , " << center[2]+Normal[0]*thickness/2 + (center[2]+Normal[2]*thickness/2-Normal[2]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]) << ". Radius = " << FAR*(Corner_max.y()-Corner_min.y()) << endl;
}
void Network::Define_fictious_cells()
@@ -531,6 +514,8 @@
(cell->info().fictious())+=1;
}
}
+
+ if(DEBUG_OUT) cout << "Fictious cell defined" << endl;
}
// double Network::spherical_triangle_area ( Sphere STA1, Sphere STA2, Sphere STA3, Point PTA1 )
=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp 2010-11-17 20:07:21 +0000
+++ lib/triangulation/Network.hpp 2010-11-25 14:22:43 +0000
@@ -28,6 +28,7 @@
int coordinate;
bool flowCondition;//flowCondition=0, pressure is imposed // flowCondition=1, flow is imposed
Real value;
+ bool useMaxMin;
};
class Network
@@ -45,7 +46,7 @@
int* boundsIds [6];
Point Corner_min;
Point Corner_max;
- Real Vsolid_tot, Vtotalissimo, Vporale, Ssolid_tot;
+ Real Vsolid_tot, Vtotalissimo, Vporale, Ssolid_tot, V_porale_porosity, V_totale_porosity;
Boundary boundaries [6];
Boundary& boundary (int b) {return boundaries[b-id_offset];}
short id_offset;
@@ -57,11 +58,14 @@
// bool facet_detected;
// void DisplayStatistics();
- void AddBoundingPlanes(bool yade);
+// void AddBoundingPlanes(bool yade);
void AddBoundingPlanes();
+ void AddBoundingPlane (bool yade, Vecteur Normal, int id_wall);
+ void AddBoundingPlane (Real center[3], double thickness, Vecteur Normal, int id_wall );
+// void AddBoundingPlanes ( Real center[3], Real Extents[3], int id );
void Define_fictious_cells( );
int Detect_facet_fictious_vertices (Cell_handle& cell, int& j);
- double Volume_Pore (Cell_handle cell);
+// double Volume_Pore (Cell_handle cell);
double Volume_Pore_VoronoiFraction ( Cell_handle& cell, int& j);
double volume_single_fictious_pore(const Vertex_handle& SV1, const Vertex_handle& SV2, const Vertex_handle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface);
double volume_double_fictious_pore(const Vertex_handle& SV1, const Vertex_handle& SV2, const Vertex_handle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface);
@@ -86,4 +90,4 @@
#endif
-#endif //FLOW_ENGINE
\ No newline at end of file
+#endif //FLOW_ENGINE
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2010-11-25 12:16:32 +0000
+++ pkg/dem/FlowEngine.cpp 2010-11-25 14:22:43 +0000
@@ -25,128 +25,119 @@
{
}
-
void FlowEngine::action()
{
+ if (!isActivated) return;
if (!flow) {
flow = shared_ptr<FlowSolver> (new FlowSolver);
first=true;
cerr <<"first = true"<<endl;
- Update_Triangulation=false;/*IS=0.f;*/
+ Update_Triangulation=false;
eps_vol_max=0.f;
Eps_Vol_Cumulative=0.f;
ReTrg=1;
retriangulationLastIter=0;
}
- if (!isActivated) return;
- else
+ if (!triaxialCompressionEngine)
{
- timingDeltas->start();
-
- if (!triaxialCompressionEngine)
- {
- vector<shared_ptr<Engine> >::iterator itFirst = scene->engines.begin();
- vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();
- for (;itFirst!=itLast; ++itFirst) {
- if ((*itFirst)->getClassName() == "TriaxialCompressionEngine") {
- cout << "stress controller engine found - FlowEngine" << endl;
- triaxialCompressionEngine = YADE_PTR_CAST<TriaxialCompressionEngine> (*itFirst);}}
- if (!triaxialCompressionEngine) cout << "stress controller engine NOT found" << endl;
- }
-
- currentStress = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];
- currentStrain = triaxialCompressionEngine->strain[1];
- current_state = triaxialCompressionEngine->currentState;
-
-// if ( current_state==3 )
- {
- if (first) Build_Triangulation(P_zero);
- timingDeltas->checkpoint("Triangulating");
-
- eps_vol_max=0.f;
- UpdateVolumes ( );
- Eps_Vol_Cumulative += eps_vol_max;
- if ((Eps_Vol_Cumulative > EpsVolPercent_RTRG || retriangulationLastIter>1000) && retriangulationLastIter>10) {
- Update_Triangulation = true;
- Eps_Vol_Cumulative=0;
- retriangulationLastIter=0;
- ReTrg++;
- } else retriangulationLastIter++;
- timingDeltas->checkpoint("Update_Volumes");
-
- if (!first) {
- ///Compute flow and and forces here
- flow->GaussSeidel();
- timingDeltas->checkpoint("Gauss-Seidel");
- if (save_mgpost) flow->MGPost();
- if (!CachedForces) flow->ComputeFacetForces();
- else flow->ComputeFacetForcesWithCache();
- timingDeltas->checkpoint("Compute_Forces");
-
- ///End Compute flow and forces
- //int number_of_particles = flow->num_particles;
- CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
- Vector3r f;
- int id;
- //IS = 0.f;
- //std::ofstream Istab ("Stability.txt", std::ios::app);
- for (CGT::Finite_vertices_iterator V_it = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); V_it != vertices_end; V_it++)
- {
- id = V_it->info().id();
- for (int y=0;y<3;y++) f[y] = (V_it->info().forces)[y];
- scene->forces.addForce(id, f);
- //scene->forces.sync();
- //IS += (scene->forces.getForce(id).norm())/number_of_particles;
- }
- timingDeltas->checkpoint("Applying Forces");
- }
- //cout << "STABILITY INDEX - IS = " << IS << endl;
- //Istab << scene->time << " " << IS << endl;
-
- Real time = scene->time;
- int j = scene->iter;
-
- if (consolidation) {
- char file [50];
- mkdir("./Consol", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
- string consol = "./Consol/"+flow->key+"%d_Consol";
- const char* keyconsol = consol.c_str();
- sprintf(file, keyconsol, j);
- char *g = file;
- timingDeltas->checkpoint("Writing cons_files");
- MaxPressure = flow->PressureProfile(g, time, intervals);
-
- std::ofstream max_p("pressures.txt", std::ios::app);
- max_p << j << " " << time << " " << MaxPressure << endl;
-
- std::ofstream settle("settle.txt", std::ios::app);
- settle << j << " " << time << " " << currentStrain << endl;
- }
-// if ( scene->iter % PermuteInterval == 0 )
-// { Update_Triangulation = true; }
-
- if (Update_Triangulation && !first) {
- Build_Triangulation();
- Update_Triangulation = false;
- }
- first=false;
-// flow->Average_Grain_Velocity();
- if (save_vtk) flow->save_vtk_file();
-
- if (slice_pressures)
- {
- char slifile [30];
- mkdir("./Slices", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
- string slice = "./Slices/Slice_"+flow->key+"_%d";
- const char* keyslice = slice.c_str();
- sprintf(slifile, keyslice, j);
- char *f = "slifile";
- flow->SliceField(f);
- }
-// int numero = flow->Average_Cell_Velocity(id_sphere, flow->T[flow->currentTes].Triangulation());
-// flow->vtk_average_cell_velocity(flow->T[flow->currentTes].Triangulation(), id_sphere, numero);
- }
- }
+ vector<shared_ptr<Engine> >::iterator itFirst = scene->engines.begin();
+ vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();
+ for (;itFirst!=itLast; ++itFirst) {
+ if ((*itFirst)->getClassName() == "TriaxialCompressionEngine") {
+ cout << "stress controller engine found - FlowEngine" << endl;
+ triaxialCompressionEngine = YADE_PTR_CAST<TriaxialCompressionEngine> (*itFirst);}}
+ if (!triaxialCompressionEngine) cout << "stress controller engine NOT found" << endl;
+ }
+ currentStress = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];
+ currentStrain = triaxialCompressionEngine->strain[1];
+
+// current_state = triaxialCompressionEngine->currentState;
+// if (current_state == 3){
+
+ timingDeltas->start();
+
+ if (first) Build_Triangulation(P_zero);
+ timingDeltas->checkpoint("Triangulating");
+
+ if (!first) {
+ eps_vol_max=0.f;
+ UpdateVolumes ( );
+ Eps_Vol_Cumulative += eps_vol_max;
+ if ((Eps_Vol_Cumulative > EpsVolPercent_RTRG || retriangulationLastIter>1000) && retriangulationLastIter>10) {
+ Update_Triangulation = true;
+ Eps_Vol_Cumulative=0;
+ retriangulationLastIter=0;
+ ReTrg++;
+ } else retriangulationLastIter++;
+ timingDeltas->checkpoint("Update_Volumes");
+
+ ///Compute flow and and forces here
+ flow->GaussSeidel();
+ timingDeltas->checkpoint("Gauss-Seidel");
+ if (save_mgpost) flow->MGPost();
+ if (!CachedForces) flow->ComputeFacetForces();
+ else flow->ComputeFacetForcesWithCache();
+ timingDeltas->checkpoint("Compute_Forces");
+
+ ///End Compute flow and forces
+ CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
+ Vector3r f;
+ int id;
+ for (CGT::Finite_vertices_iterator V_it = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); V_it != vertices_end; V_it++)
+ {
+ id = V_it->info().id();
+ for (int y=0;y<3;y++) f[y] = (V_it->info().forces)[y];
+ scene->forces.addForce(id, f);
+ }
+ timingDeltas->checkpoint("Applying Forces");
+
+ Real time = scene->time;
+ int j = scene->iter;
+
+ if (consolidation) {
+ char file [50];
+ mkdir("./Consol", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ string consol = "./Consol/"+flow->key+"%d_Consol";
+ const char* keyconsol = consol.c_str();
+ sprintf(file, keyconsol, j);
+ char *g = file;
+ timingDeltas->checkpoint("Writing cons_files");
+ MaxPressure = flow->PressureProfile(g, time, intervals);
+
+ std::ofstream max_p("pressures.txt", std::ios::app);
+ max_p << j << " " << time << " " << MaxPressure << endl;
+
+ std::ofstream settle("settle.txt", std::ios::app);
+ settle << j << " " << time << " " << currentStrain << endl;
+ }
+
+ if (slice_pressures){
+ char slifile [30];
+ mkdir("./Slices", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ string slice = "./Slices/Slice_"+flow->key+"_%d";
+ const char* keyslice = slice.c_str();
+ sprintf(slifile, keyslice, j);
+ char *f = "slifile";
+ flow->SliceField(f);
+ }
+
+ if (save_vtk) flow->save_vtk_file();
+ }
+// if ( scene->iter % PermuteInterval == 0 )
+// { Update_Triangulation = true; }
+
+ if (Update_Triangulation && !first) {
+ Build_Triangulation();
+ Update_Triangulation = false;
+ }
+
+ flow->Average_Cell_Velocity();
+// flow->Average_Grain_Velocity();
+// int numero = flow->Average_Cell_Velocity(id_sphere, flow->T[flow->currentTes].Triangulation());
+// flow->vtk_average_cell_velocity(flow->T[flow->currentTes].Triangulation(), id_sphere, numero);
+
+ first=false;
+// }
}
void FlowEngine::BoundaryConditions()
@@ -195,13 +186,15 @@
Triangulate ( );
if (Debug) cout << endl << "Tesselating------" << endl << endl;
flow->T[flow->currentTes].Compute();
-
+
flow->Define_fictious_cells();
flow->DisplayStatistics ();
flow->meanK_LIMIT = meanK_correction;
flow->meanK_STAT = meanK_opt;
flow->Compute_Permeability ( );
+
+ porosity = flow->V_porale_porosity/flow->V_totale_porosity;
if (first)
{
@@ -236,21 +229,19 @@
void FlowEngine::AddBoundary ()
{
shared_ptr<Sphere> sph ( new Sphere );
-
int Sph_Index = sph->getClassIndexStatic();
- int contator = 0;
FOREACH ( const shared_ptr<Body>& b, *scene->bodies )
{
if ( !b ) continue;
if ( b->shape->getClassIndex() == Sph_Index )
{
- Sphere* s=YADE_CAST<Sphere*> ( b->shape.get() );
+ Sphere* s=YADE_CAST<Sphere*> ( b->shape.get() );
//const Body::id_t& id = b->getId();
- Real rad = s->radius;
- Real x = b->state->pos[0];
- Real y = b->state->pos[1];
- Real z = b->state->pos[2];
+ Real& rad = s->radius;
+ Real& x = b->state->pos[0];
+ Real& y = b->state->pos[1];
+ Real& z = b->state->pos[2];
flow->x_min = min ( flow->x_min, x-rad);
flow->x_max = max ( flow->x_max, x+rad);
@@ -258,14 +249,46 @@
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);
-
- contator+=1;
}
}
+
+ 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=triaxialCompressionEngine->wall_bottom_id;
+ flow->y_max_id=triaxialCompressionEngine->wall_top_id;
+ flow->x_max_id=triaxialCompressionEngine->wall_right_id;
+ flow->x_min_id=triaxialCompressionEngine->wall_left_id;
+ flow->z_min_id=triaxialCompressionEngine->wall_back_id;
+ flow->z_max_id=triaxialCompressionEngine->wall_front_id;
+
+ flow->boundary ( flow->y_min_id ).useMaxMin = BOTTOM_Boundary_MaxMin;
+ flow->boundary ( flow->y_max_id ).useMaxMin = TOP_Boundary_MaxMin;
+ flow->boundary ( flow->x_max_id ).useMaxMin = RIGHT_Boundary_MaxMin;
+ flow->boundary ( flow->x_min_id ).useMaxMin = LEFT_Boundary_MaxMin;
+ flow->boundary ( flow->z_max_id ).useMaxMin = FRONT_Boundary_MaxMin;
+ flow->boundary ( flow->z_min_id ).useMaxMin = BACK_Boundary_MaxMin;
+
+ //FIXME: Id's order in boundsIds is done according to the enumerotation of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
+ flow->boundsIds[0]= &flow->y_min_id;
+ flow->boundsIds[1]= &flow->y_max_id;
+ flow->boundsIds[2]= &flow->x_min_id;
+ flow->boundsIds[3]= &flow->x_max_id;
+ flow->boundsIds[4]= &flow->z_max_id;
+ flow->boundsIds[5]= &flow->z_min_id;
+
+ wall_thickness = triaxialCompressionEngine->thickness;
+
+ 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 (flow->DEBUG_OUT) {cout << "Section area = " << flow->SectionArea << endl;
+ if (Debug) {
+ cout << "Section area = " << flow->SectionArea << endl;
cout << "Vtotale = " << flow->Vtotale << endl;
// cout << "Rmoy " << Rmoy << endl;
cout << "x_min = " << flow->x_min << endl;
@@ -273,37 +296,22 @@
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;}
-
- if (Debug) cout << "Adding Boundary------" << endl;
-
-// shared_ptr<Box> bx ( new Box );
-// int Bx_Index = bx->getClassIndexStatic();
-//
-// FOREACH ( const shared_ptr<Body>& b, *scene->bodies )
-// {
-// if ( !b ) continue;
-// if ( b->shape->getClassIndex() == Bx_Index )
-// {
-// Box* w = YADE_CAST<Box*> ( b->shape.get() );
-// // const Body::id_t& id = b->getId();
-// Real center [3], Extent[3];
-// for ( int h=0;h<3;h++ ) {center[h] = b->state->pos[h]; Extent[h] = w->extents[h];}
-// wall_thickness = min(min(Extent[0],Extent[1]),Extent[2]);
-// }
-// }
- id_offset = flow->T[flow->currentTes].max_id+1;
-
- flow->id_offset = id_offset;
-
- flow->y_min_id=triaxialCompressionEngine->wall_bottom_id;
- flow->y_max_id=triaxialCompressionEngine->wall_top_id;
- flow->x_max_id=triaxialCompressionEngine->wall_right_id;
- flow->x_min_id=triaxialCompressionEngine->wall_left_id;
- flow->z_min_id=triaxialCompressionEngine->wall_back_id;
- flow->z_max_id=triaxialCompressionEngine->wall_front_id;
-
- flow->AddBoundingPlanes(true);
+ cout << "z_max = " << flow->z_max << endl;
+ cout << endl << "Adding Boundary------" << endl;}
+
+ double center[3];
+
+ for (int i=0; i<6; i++)
+ {
+ CGT::Vecteur Normal (triaxialCompressionEngine->normal[i].x(), triaxialCompressionEngine->normal[i].y(), triaxialCompressionEngine->normal[i].z());
+ if (flow->boundary(*flow->boundsIds[i]).useMaxMin) flow->AddBoundingPlane (true, Normal, *flow->boundsIds[i]);
+ else {
+ const shared_ptr<Body>& wll = Body::byId ( *flow->boundsIds[i] , scene );
+ for ( int h=0;h<3;h++ ){center[h] = wll->state->pos[h];}
+ flow->AddBoundingPlane (center, wall_thickness, Normal,*flow->boundsIds[i]);}}
+
+// flow->AddBoundingPlanes(true);
+
}
void FlowEngine::Triangulate ()
@@ -354,7 +362,6 @@
void FlowEngine::UpdateVolumes ()
{
if (Debug) cout << "Updating volumes.............." << endl;
- Real deltaT = scene->dt;
Real invDeltaT = 1/scene->dt;
CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
double newVol=0; double dVol;
@@ -424,10 +431,7 @@
Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0};
Real B[3]={0, 0, 0}, BS[3]={0, 0, 0}, BT[3]={0, 0, 0};
Real C[3]={0, 0, 0}, CS[3]={0, 0, 0}, CT[3]={0, 0, 0};
-
- //Real A[3], AS[3], AT[3];
- //Real B[3], BS[3], BT[3];
- //Real C[3], CS[3], CT[3];
+
int b[2];
Real Wall_point[2][3];
@@ -476,8 +480,7 @@
Real FlowEngine::Volume_cell_triple_fictious ( CGT::Cell_handle cell)
{
Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0}, AW[3]={0, 0, 0};
-//Real A[3], AS[3], AT[3], AW[3];
-// CGT::Boundary b[3];
+
int b[3];
Real Wall_point[3][3];
int j=0;
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2010-11-19 17:18:56 +0000
+++ pkg/dem/FlowEngine.hpp 2010-11-25 14:22:43 +0000
@@ -73,6 +73,7 @@
((int,PermuteInterval,100000,,"Pore space re-triangulation period"))
((double, eps_vol_max, 0,,"Maximal absolute volumetric strain computed at each iteration"))
((double, EpsVolPercent_RTRG,0.01,,"Percentuage of cumulate eps_vol at which retriangulation of pore space is performed"))
+ ((double, porosity, 0,,"Porosity computed at each retriangulation"))
((bool,compute_K,true,,"Activates permeability measure within a granular sample"))
((bool,meanK_correction,true,,"Local permeabilities' correction through meanK threshold"))
((bool,meanK_opt,false,,"Local permeabilities' correction through an optimized threshold"))
@@ -83,7 +84,7 @@
((double, currentStress, 0,, "Current value of axial stress"))
((double, currentStrain, 0,, "Current value of axial strain"))
((int, intervals, 30,, "Number of layers for pressure measurements"))
- ((int, useSolver, 1,, "Solver to use"))
+ ((int, useSolver, 0,, "Solver to use"))
((bool, Flow_imposed_TOP_Boundary, true,, "if false involve pressure imposed condition"))
((bool, Flow_imposed_BOTTOM_Boundary, true,, "if false involve pressure imposed condition"))
((bool, Flow_imposed_FRONT_Boundary, true,, "if false involve pressure imposed condition"))
@@ -98,6 +99,12 @@
((double, Pressure_RIGHT_Boundary, 0,, "Pressure imposed on right boundary"))
((double, Sinus_Pressure, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied"))
((int, id_sphere, 0,, "Average velocity will be computed for all cells incident to that sphere"))
+ ((bool, BOTTOM_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
+ ((bool, TOP_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
+ ((bool, RIGHT_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
+ ((bool, LEFT_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
+ ((bool, FRONT_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
+ ((bool, BACK_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
,timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas));
DECLARE_LOGGER;
};
