yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04758
[Branch ~yade-dev/yade/trunk] Rev 2276: - Walls proper identification is no more ID dependent
------------------------------------------------------------
revno: 2276
committer: ecatalano <ecatalano@dt-rv019>
branch nick: trunk
timestamp: Fri 2010-06-04 14:07:46 +0200
message:
- Walls proper identification is no more ID dependent
- New output files for fluid pressure and settlement evolution
- Fixed some warning came out with new compiler version
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.h
pkg/dem/Engine/PartialEngine/FlowEngine.cpp
pkg/dem/Engine/PartialEngine/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-04-23 09:35:10 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-06-04 12:07:46 +0000
@@ -1601,7 +1601,7 @@
{cell->info().p() = P_zero;cell->info().dv()=0;}
for (int bound=0; bound<6;bound++) {
- int& id = *boundsIds[bound];
+ int id = boundsIds[bound];
Boundary& bi = boundary(id);
if (!bi.flowCondition) {
Tesselation::Vector_Cell tmp_cells;
@@ -1682,7 +1682,7 @@
iter << j << " " << dp_max/p_max << endl;
int cel=0;
- double Pav;
+ double Pav=0;
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
cel++;
Pav+=cell->info().p();
@@ -1987,7 +1987,6 @@
char *kk;
kk = (char*) key.c_str();
-
return Permeameter(Tri, boundary(y_min_id).value, boundary(y_max_id).value, Section, DeltaY, kk);
}
@@ -1998,6 +1997,7 @@
Corner_min = Point(x_min, y_min, z_min);
Corner_max = Point(x_max, y_max, z_max);
Real min_coord = min(Extents[0],min(Extents[1],Extents[2]));
+
int coord=0;
if (min_coord==Extents[0]) {
coord=0;
@@ -2057,59 +2057,64 @@
void FlowBoundingSphere::AddBoundingPlanes()
{
+ Tesselation& Tes = T[currentTes];
+
+ y_min_id = Tes.Max_id() + 1; boundsIds[0]=y_min_id;
+ y_max_id = Tes.Max_id() + 2; boundsIds[1]=y_max_id;
+ x_min_id = Tes.Max_id() + 3; boundsIds[2]=x_min_id;
+ x_max_id = Tes.Max_id() + 4; boundsIds[3]=x_max_id;
+ z_min_id = Tes.Max_id() + 5; boundsIds[4]=z_min_id;
+ z_max_id = Tes.Max_id() + 6; boundsIds[5]=z_max_id;
+
+ id_offset = Tes.Max_id() +1;//so that boundaries[vertex->id - offset] gives the ordered boundaries (also see function Boundary& boundary(int b))
+
+ AddBoundingPlanes(y_min_id, y_max_id, x_min_id, x_max_id, z_min_id, z_max_id);
+}
+
+void FlowBoundingSphere::AddBoundingPlanes(int bottom_id, int top_id, int left_id, int right_id, int front_id, int back_id)
+{
Tesselation& Tes = T[currentTes];
Corner_min = Point(x_min, y_min, z_min);
Corner_max = Point(x_max, y_max, z_max);
-
- y_min_id = Tes.Max_id() +1;
- boundsIds[0]=&y_min_id;
- y_max_id = Tes.Max_id() +2;
- boundsIds[1]=&y_max_id;
- x_min_id = Tes.Max_id() +3;
- boundsIds[2]=&x_min_id;
- x_max_id = Tes.Max_id() +4;
- boundsIds[3]=&x_max_id;
- z_min_id = Tes.Max_id() +5;
- boundsIds[4]=&z_min_id;
- z_max_id = Tes.Max_id() +6;
- boundsIds[5]=&z_max_id;
-
- id_offset = Tes.Max_id() +1;//so that boundaries[vertex->id - offset] gives the ordered boundaries (also see function Boundary& boundary(int b))
+
+ y_min_id = bottom_id;y_max_id = top_id;x_min_id = left_id;x_max_id = right_id;z_min_id = front_id;z_max_id = back_id;
+
+ 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.x()-Corner_min.x()), 0.5*(Corner_max.z()-Corner_min.z()), FAR*(Corner_max.x()-Corner_min.x()), y_min_id, true);
- boundaries[0].p = Corner_min;
- boundaries[0].normal = Vecteur(0,1,0);
- boundaries[0].coordinate = 1;
+ 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;
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.x()-Corner_min.x()), 0.5*(Corner_max.z()-Corner_min.z()), FAR*(Corner_max.x()-Corner_min.x()), y_max_id, true);
- boundaries[1].p = Corner_max;
- boundaries[1].normal = Vecteur(0,-1,0);
- boundaries[1].coordinate = 1;
+ 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;
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[2].p = Corner_min;
- boundaries[2].normal = Vecteur(1,0,0);
- boundaries[2].coordinate = 0;
+ 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;
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[3].p = Corner_max;
- boundaries[3].normal = Vecteur(-1,0,0);
- boundaries[3].coordinate = 0;
+ 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;
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[4].p = Corner_min;
- boundaries[4].normal = Vecteur(0,0,1);
- boundaries[4].coordinate = 2;
+ 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;
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[5].p = Corner_max;
- boundaries[5].normal = Vecteur(0,0,-1);
- boundaries[5].coordinate = 2;
+ 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;
cout << "Back boundary has been created. ID = " << z_max_id << endl;
for (int k=0;k<6;k++) {
=== modified file 'lib/triangulation/FlowBoundingSphere.h'
--- lib/triangulation/FlowBoundingSphere.h 2010-04-23 09:35:10 +0000
+++ lib/triangulation/FlowBoundingSphere.h 2010-06-04 12:07:46 +0000
@@ -34,7 +34,7 @@
FlowBoundingSphere();
int x_min_id, x_max_id, y_min_id, y_max_id, z_min_id, z_max_id;
- int* boundsIds [6];
+ int boundsIds [6];
bool currentTes;
bool SLIP_ON_LATERALS;
double TOLERANCE;
@@ -45,15 +45,13 @@
int Iterations;
Boundary boundaries [6];
+ int walls_id[6];
short id_offset;
Boundary& boundary (int b) {return boundaries[b-id_offset];}
void mplot (RTriangulation& Tri, char *filename);
void Localize ();
void Compute_Permeability();
- void AddBoundingPlanes();
-
- void AddBoundingPlanes(Real center[3], Real Extents[3], int id);
void DisplayStatistics();
void GaussSeidel ( );
@@ -72,7 +70,12 @@
Real minPermLength; //min branch length for Poiseuille
double P_SUP, P_INF, P_INS;
+
void AddBoundingPlanes ( Tesselation& Tes, double x_Min,double x_Max ,double y_Min,double y_Max,double z_Min,double z_Max );
+ void AddBoundingPlanes(int y_min_id, int y_max_id, int x_min_id, int x_max_id, int z_min_id, int z_max_id);
+ void AddBoundingPlanes();
+ void AddBoundingPlanes(Real center[3], Real Extents[3], int id);
+
void Compute_Action ( );
void Compute_Action ( int argc, char *argv[ ], char *envp[ ] );
void DisplayStatistics ( RTriangulation& Tri );
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-04-23 09:35:10 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-06-04 12:07:46 +0000
@@ -28,7 +28,8 @@
void FlowEngine::action ( )
{
- if (!flow) {flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);first=true;Update_Triangulation=false;}
+ if (!flow)
+ {flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);first=true;Update_Triangulation=false;}
if ( !isActivated ) return;
else
{
@@ -47,7 +48,7 @@
if ( current_state==3 )
{
- if ( first ) { Build_Triangulation( P_zero );first=false;}
+ if ( first ) { Build_Triangulation( P_zero );}
timingDeltas->checkpoint("Triangulating");
@@ -73,7 +74,8 @@
// flow->MGPost(flow->T[flow->currentTes].Triangulation());
flow->ComputeTetrahedralForces();
-
+// flow->Compute_Forces();
+
timingDeltas->checkpoint("Compute_Forces");
///End Compute flow and forces
@@ -104,12 +106,21 @@
MaxPressure = flow->PermeameterCurve(flow->T[flow->currentTes].Triangulation(), g, time, intervals);
+ std::ofstream max_p ("pressures.txt", std::ios::app);
+ MaxPressure = flow->PermeameterCurve(flow->T[flow->currentTes].Triangulation(), g, time, intervals);
+ max_p << j << " " << time << " " << MaxPressure << endl;
+
+ std::ofstream settle ("settle.txt", std::ios::app);
+ settle << j << " " << time << " " << currentStrain << endl;
+
if ( Omega::instance().getCurrentIteration() % PermuteInterval == 0 )
{ Update_Triangulation = true; }
if ( Update_Triangulation ) { Build_Triangulation( );}
timingDeltas->checkpoint("Storing Max Pressure");
+
+ first=false;
}
}
}
@@ -162,6 +173,7 @@
flow->T[flow->currentTes].Compute();
flow->Localize ();
+ flow->DisplayStatistics ();
flow->meanK_LIMIT = meanK_correction;
flow->meanK_STAT = meanK_opt;
@@ -186,9 +198,7 @@
Update_Triangulation=!Update_Triangulation;
Oedometer_Boundary_Conditions();
}
-
Initialize_volumes ( );
- flow->DisplayStatistics ();
}
void FlowEngine::AddBoundary ()
@@ -216,10 +226,13 @@
flow->y_max = max ( flow->y_max, center[1]-wall_thickness);
flow->z_min = min ( flow->z_min, center[2]+wall_thickness);
flow->z_max = max ( flow->z_max, center[2]-wall_thickness);
+
}
}
- flow->AddBoundingPlanes();
+ flow->id_offset = flow->T[flow->currentTes].max_id+1;
+
+ flow->AddBoundingPlanes( triaxialCompressionEngine->wall_bottom_id, triaxialCompressionEngine->wall_top_id, triaxialCompressionEngine->wall_left_id, triaxialCompressionEngine->wall_right_id, triaxialCompressionEngine->wall_front_id, triaxialCompressionEngine->wall_back_id );
}
void FlowEngine::Triangulate ()
@@ -349,9 +362,11 @@
Real FlowEngine::Volume_cell_double_fictious ( CGT::Cell_handle cell)
{
- Real A[3], AS[3], AT[3];
- Real B[3], BS[3], BT[3];
- Real C[3], CS[3], CT[3];
+// Real array_quattro[3] = {0, 0, 0};
+
+ 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};
int b[2];
Real Wall_point[2][3];
@@ -401,7 +416,7 @@
Real FlowEngine::Volume_cell_triple_fictious ( CGT::Cell_handle cell)
{
- Real A[3], AS[3], AT[3], AW[3];
+ Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0}, AW[3]={0, 0, 0};
// CGT::Boundary b[3];
int b[3];
Real Wall_point[3][3];
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.hpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-04-23 09:35:10 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-06-04 12:07:46 +0000
@@ -54,7 +54,7 @@
((double,Relax,1.9,"Gauss-Seidel relaxation"))
((int,PermuteInterval,100000,"Pore space re-triangulation period"))
((bool,compute_K,true,"Activates permeability measure within a granular sample"))
- ((bool,meanK_correction,false,"Local permeabilities' correction through meanK threshold"))
+ ((bool,meanK_correction,true,"Local permeabilities' correction through meanK threshold"))
((bool,meanK_opt,false,"Local permeabilities' correction through an optimized threshold"))
((double,permeability_factor,1.0,"a permability multiplicator"))
((Real,loadFactor,1.1,"Load multiplicator for oedometer test"))
