yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #03398
[Branch ~yade-dev/yade/trunk] Rev 2028: - Made corrections on flow code
------------------------------------------------------------
revno: 2028
committer: Emanuele Catalano <ecatalano@r2balme>
branch nick: trunk
timestamp: Mon 2010-02-15 18:58:41 +0100
message:
- Made corrections on flow code
- Class registration via new macros introduced
modified:
lib/triangulation/FlowBoundingSphere.cpp
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-02-14 22:32:05 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-02-15 17:58:41 +0000
@@ -226,8 +226,6 @@
// Boundary_Conditions ( Tri );
- P_SUP = 1, P_INF = 0, P_INS=0;
-
Initialize_pressures();
GaussSeidel();
@@ -237,7 +235,7 @@
/** END GAUSS SEIDEL */
char *file = "Permeability";
- Permeameter ( Tri, P_INF, P_SUP, SectionArea, H, file );
+ Permeameter ( Tri, boundary ( y_min_id ).value, boundary ( y_max_id ).value, SectionArea, H, file );
clock.top ( "Permeameter" );
Compute_Forces();
@@ -382,30 +380,26 @@
for ( Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++ )
{
-// cell->info().fictious() = 0;
+ cell->info().fictious() = 0;
pass = 0;
for ( int i=0;i<4;i++ )
{
V = cell->vertex ( i );
if ( V->info().isFictious )
{
-// cell->info().fictious() ++;
+ cell->info().fictious() ++;
pass = 1;
- //FIXME : remove the isFictious flag and use cell->info().fictious() instead
-// cell->info().isFictious = true;
- Boundary& bi = boundary ( V->info().id() );
- // Boundary& bi = boundaries [V->info().id()];
-
- if ( bi.flowCondition )
+ // Boundary& bi = boundary ( V->info().id() );
+ // if ( bi.flowCondition )
+ if ( V->info().id() != y_min_id && V->info().id() != y_max_id)
{
- // Inf/Sup has priority
if ( !cell->info().isSuperior && !cell->info().isInferior )
{
- cell->info().isLateral = true;
+ cell->info().isLateral = true;cell->info().isInside=false; cell->info().isInferior=false; cell->info().isSuperior=false;
#ifdef XVIEW
Vue1.SetSpheresColor ( 0.1,0.95,0.1,1 );
- Vue1.Dessine_Sphere ( cell->info().x(), cell->info().y(), cell->info().z(), 0.02 , 4 );
+ Vue1.Dessine_Sphere (cell->info().x(), cell->info().y(), cell->info().z(), 0.02 , 4 );
#endif
}
}
@@ -413,9 +407,12 @@
{
if ( V->info().id() == y_min_id )
{
- cell->info().isInferior=true; cell->info().isLateral=false; cell->info().isSuperior=false;
- }
- else {cell->info().isSuperior=true;cell->info().isLateral=false;cell->info().isInferior=false;}
+ cell->info().isInferior=true; cell->info().isLateral=false; cell->info().isSuperior=false;cell->info().isInside=false;
+ }
+ else{
+ cell->info().isSuperior=true;cell->info().isLateral=false;
+ cell->info().isInferior=false;cell->info().isInside=false;
+ }
#ifdef XVIEW
Vue1.SetSpheresColor ( 1,0.1,0.1,1 );
Vue1.Dessine_Sphere ( cell->info().x(), cell->info().y(), cell->info().z(), 0.02 , 4 );
@@ -423,7 +420,7 @@
}
}
}
- //if (!pass) cell->info().isInside=true;
+ if (!pass) cell->info().isInside=true;
}
cout << "Localised -------------" << endl;
}
@@ -491,8 +488,6 @@
// else if ( Tri.is_infinite ( neighbour_cell )) k = 0;//connection to an infinite cells
}
cell->info().isvisited = !ref;
- // for (int y=0;y<4;y++) cout << "Permeability " << y << " = " << (cell->info().k_norm())[y] << endl;
- // for ( int y=0;y<4;y++ ) kFile << y << " = " << ( cell->info().k_norm() ) [y] << endl;
}
cout << "POS = " << POS << " NEG = " << NEG << " pass = " << pass << endl;
@@ -825,7 +820,7 @@
{
if ( !v->info().isFictious )
{
- cout << "id = " << v->info().id() << " force = " << v->info().forces << endl;
+// cout << "id = " << v->info().id() << " force = " << v->info().forces << endl;
TotalForce = TotalForce + v->info().forces;
}
else
@@ -1480,7 +1475,7 @@
// else cell->info().p() =0;//FIXME : assign better values for faster convergence?
// }
- for ( int bound=0; bound<6;bound++ )
+ for ( int bound=0;bound<6;bound++ )
{
int& id = *boundsIds[bound];
Boundary& bi = boundary ( id );
@@ -1490,7 +1485,8 @@
tmp_cells.resize ( 1000 );
Tesselation::VCell_iterator cells_it = tmp_cells.begin();
Tesselation::VCell_iterator cells_end = Tri.incident_cells ( T[currentTes].vertexHandles[id],cells_it );
- for ( Tesselation::VCell_iterator it = tmp_cells.begin(); it != cells_end; it++ ) ( *it )->info().p() = bi.value;
+ for ( Tesselation::VCell_iterator it = tmp_cells.begin(); it != cells_end; it++ )
+ ( *it )->info().p() = bi.value;
}
}
}
@@ -1554,12 +1550,12 @@
if ( j % 1000 == 0 )
{
- cout << "pmax " << p_max << "; pmoy : " << p_moy << endl;
+ cout << "pmax " << p_max << "; pmoy : " << p_moy << "; dpmax : " << dp_max << endl;
cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;
// save_vtk_file ( Tri );
}
}
- while ( ( dp_max/p_max ) > tolerance /*( dp_max > tolerance )*//* &&*/ /*( j<50 )*/ );
+ while ( ( dp_max/*/p_max*/ ) > tolerance /*( dp_max > tolerance )*//* &&*/ /*( j<50 )*/ );
// cout << "pmax " << p_max << "; pmoy : " << p_moy << endl;
// cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;
}
@@ -1587,7 +1583,7 @@
{
//PRESSION_EXT=P_Inf;
//Qout = Qout + ( cell->info().k_vector() ) [j2]* ( cell->info().p()-PRESSION_EXT );
- cout << "outFlow : "<< ( cell->info().k_norm() ) [j2]* ( neighbour_cell->info().p()-cell->info().p() ) <<endl;
+// cout << "outFlow : "<< ( cell->info().k_norm() ) [j2]* ( neighbour_cell->info().p()-cell->info().p() ) <<endl;
Qout = Qout + ( cell->info().k_norm() ) [j2]* ( neighbour_cell->info().p()-cell->info().p() );
cellQout+=1;/*( cell->info().k_vector() ) [j2]* ( cell->info().p()-PRESSION_EXT )*/;
p_out_max = std::max ( cell->info().p(), p_out_max );
@@ -1838,7 +1834,11 @@
double DeltaY = y_Max-y_Min;
- P_SUP = 1.0; P_INF = 0.0; P_INS = 0.0;
+ boundary ( y_min_id ).flowCondition=0;
+ boundary ( y_max_id ).flowCondition=0;
+ boundary ( y_min_id ).value=0;
+ boundary ( y_max_id ).value=1;
+// P_SUP = 1.0; P_INF = 0.0; P_INS = 0.0;
Initialize_pressures();
@@ -1847,7 +1847,7 @@
char *kk;
kk = ( char* ) key.c_str();
- Permeameter ( Tri, P_INF, P_SUP, Section, DeltaY, kk );
+ Permeameter ( Tri, boundary ( y_min_id ).value, boundary ( y_max_id ).value, Section, DeltaY, kk );
}
void FlowBoundingSphere::AddBoundingPlanes ( Real center[3], Real Extents[3], int id )
@@ -1951,31 +1951,37 @@
boundaries[0].p = Corner_min;
boundaries[0].normal = Vecteur ( 0,1,0 );
boundaries[0].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;
+ 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;
+ 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;
+ 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;
+ 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;
+ cout << "Back boundary has been created. ID = " << z_max_id << endl;
for ( int k=0;k<6;k++ )
{
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-02-14 22:32:05 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-02-15 17:58:41 +0000
@@ -20,27 +20,13 @@
std::ofstream plotFile ( "plot2",std::ios::out );
-FlowEngine::FlowEngine() : gravity ( Vector3r::ZERO ), isActivated ( true )
-{
- first = true;
- currentTes = 0;
- P_zero=0;
- PermuteInterval = 100000;
- permeability_factor=1.0;
- loadFactor=1.0;
- compute_K=true;
- unload=false;
- tess_based_force=true;
- flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);
-}
-
-
FlowEngine::~FlowEngine()
{
}
void FlowEngine::applyCondition ( Scene* ncb )
{
+ if (!flow) {flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);first=true;}
if ( !isActivated ) return;
else
{
@@ -141,57 +127,34 @@
{
Initialize ( ncb, P_zero );
- cout << "ok2" << endl;
-
- Oedometer_Boundary_Conditions();
-
- cout << "ok3" << endl;
-
flow->Vtotalissimo=0; flow->Vsolid_tot=0; flow->Vporale=0; flow->Ssolid_tot=0;
-
+
flow->k_factor = permeability_factor;
-
flow->Compute_Permeability ();
- cout << "Vtotalissimo = " << flow->Vtotalissimo << " Vsolid_tot = " << flow->Vsolid_tot << " Vporale2 = " << flow->Vporale << " Ssolid_tot = " << flow->Ssolid_tot << endl << endl;
-
flow->DisplayStatistics ();
CGT::Finite_cells_iterator cell_end = flow->T[currentTes].Triangulation().finite_cells_end();
-
int y=0;
for ( CGT::Finite_cells_iterator cell = flow->T[currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ )
{
-// cell->info().p() = P_zero;
- cell->info().dv() = 0;
+ cell->info().dv() = 0; cell->info().p() = 0;
y++;
}
cout << y << " deltaV initialised -----------------" << endl;
if (compute_K) {flow->Sample_Permeability ( flow->T[currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
- double P_ext=0, P_int=0.0;
-
- flow->P_SUP=P_ext; flow->P_INF=P_ext; flow->P_INS=P_int;
-
- flow->Initialize_pressures( );
+ Oedometer_Boundary_Conditions();
+ flow->Initialize_pressures();
flow->GaussSeidel ( );
plotFile << "unset key" << endl;
- triaxialCompressionEngine->wall_left_activated=0;
- triaxialCompressionEngine->wall_right_activated=0;
- triaxialCompressionEngine->wall_front_activated=0;
- triaxialCompressionEngine->wall_back_activated=0;
- triaxialCompressionEngine->wall_top_activated=1;
- triaxialCompressionEngine->wall_bottom_activated=1;
-
- triaxialCompressionEngine->sigma_iso=(triaxialCompressionEngine->sigma_iso)*loadFactor;
-
// flow->Analytical_Consolidation();
- first = false;cons=0;
+ first = false; cons=0;
}
}
}
@@ -201,7 +164,16 @@
flow->boundary ( flow->y_min_id ).flowCondition=0;
flow->boundary ( flow->y_max_id ).flowCondition=0;
flow->boundary ( flow->y_min_id ).value=0;
- flow->boundary ( flow->y_max_id ).value=1;
+ flow->boundary ( flow->y_max_id ).value=0;
+
+ triaxialCompressionEngine->wall_left_activated=0;
+ triaxialCompressionEngine->wall_right_activated=0;
+ triaxialCompressionEngine->wall_front_activated=0;
+ triaxialCompressionEngine->wall_back_activated=0;
+ triaxialCompressionEngine->wall_top_activated=1;
+ triaxialCompressionEngine->wall_bottom_activated=1;
+
+ triaxialCompressionEngine->sigma_iso=(triaxialCompressionEngine->sigma_iso)*loadFactor;
}
void FlowEngine::Initialize ( Scene* ncb, double P_zero )
@@ -217,13 +189,11 @@
cout << endl << "Tesselating------" << endl << endl;
flow->T[currentTes].Compute();
- flow->Fictious_cells();
+// flow->Fictious_cells();
flow->Localize ();
Initialize_volumes ( ncb );
-
- cout << "ok1" << endl;
}
void FlowEngine::AddBoundary ( Scene* ncb )
@@ -233,38 +203,27 @@
shared_ptr<Box> bx ( new Box );
int Bx_Index = bx->getClassIndexStatic();
- int contator=0;
-
FOREACH ( const shared_ptr<Body>& b, *ncb->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();
- cout << "is it a wall?? id == " << id << endl;
-
Real center [3], Extent[3];
-
for ( int h=0;h<3;h++ ) {center[h] = b->state->pos[h]; Extent[h] = w->extents[h];}
// flow->AddBoundingPlanes ( center, Extent, id );
-
- flow->x_min = min ( flow->x_min, center[0]);
- flow->x_max = max ( flow->x_max, center[0]);
- flow->y_min = min ( flow->y_min, center[1]);
- flow->y_max = max ( flow->y_max, center[1]);
- flow->z_min = min ( flow->z_min, center[2]);
- flow->z_max = max ( flow->z_max, center[2]);
-
- contator+=1;
+ flow->x_min = min ( flow->x_min, center[0]-wall_thickness);
+ flow->x_max = max ( flow->x_max, center[0]+wall_thickness);
+ flow->y_min = min ( flow->y_min, center[1]-wall_thickness);
+ 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();
-
- cout << contator << " walls inserted -------- ADDED BOUNDING PLANES" << endl;
}
void FlowEngine::Triangulate ( Scene* ncb )
@@ -275,8 +234,6 @@
int Sph_Index = sph->getClassIndexStatic();
int contator = 0;
-
-// std::list<CGT::Point> input;
FOREACH ( const shared_ptr<Body>& b, *ncb->bodies )
{
@@ -292,21 +249,14 @@
flow->T[currentTes].insert(x, y, z, rad, id);
-// 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 );
-
contator+=1;
}
}
cout << contator << "spheres inserted " << endl;
- double SectionArea = ( flow->x_max-flow->x_min ) * ( flow->z_max-flow->z_min );
+ double SectionArea = ( flow->x_max - flow->x_min ) * ( flow->z_max-flow->z_min );
- cout << "section area = " << SectionArea << endl;
+ cout << "Section area = " << SectionArea << endl;
// cout << "Rmoy " << Rmoy << endl;
cout << "x_min = " << flow->x_min << endl;
cout << "x_max = " << flow->x_max << endl;
@@ -314,7 +264,6 @@
cout << "y_min = " << flow->y_min << endl;
cout << "z_min = " << flow->z_min << endl;
cout << "z_max = " << flow->z_max << endl;
-// cout << "SectionArea = " << SectionArea << endl;
}
void FlowEngine::NewTriangulation ( Scene* ncb )
@@ -355,7 +304,6 @@
void FlowEngine::Initialize_volumes ( Scene* ncb )
{
CGT::Finite_cells_iterator cell_end = flow->T[currentTes].Triangulation().finite_cells_end();
-
for ( CGT::Finite_cells_iterator cell = flow->T[currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ )
{
switch ( cell->info().fictious() )
@@ -374,39 +322,31 @@
cout << "Updating volumes.............." << endl;
Real deltaT = ncb->dt;
-
- cout << "deltaT = " << deltaT << endl;
-
CGT::Finite_cells_iterator cell_end = flow->T[currentTes].Triangulation().finite_cells_end();
-
for ( CGT::Finite_cells_iterator cell = flow->T[currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ )
{
switch ( cell->info().fictious() )
{
- case ( 3 ) :
+ case ( 3 ):
{
cell->info().dv() = ( Volume_cell_triple_fictious ( cell,ncb ) - cell->info().volume() ) /deltaT;
cell->info().volume() = Volume_cell_triple_fictious ( cell,ncb );
- }
- break;
+ }break;
case ( 2 ) :
{
cell->info().dv() = ( Volume_cell_double_fictious ( cell,ncb )-cell->info().volume() ) /deltaT;
cell->info().volume() = Volume_cell_double_fictious ( cell,ncb );
- }
- break;
+ }break;
case ( 1 ) :
{
cell->info().dv() = ( Volume_cell_single_fictious ( cell,ncb )-cell->info().volume() ) /deltaT;
cell->info().volume() = Volume_cell_single_fictious ( cell,ncb );
- }
- break;
+ }break;
case ( 0 ) :
{
cell->info().dv() = ( Volume_cell ( cell,ncb )-cell->info().volume() ) /deltaT;
cell->info().volume() = Volume_cell ( cell,ncb );
- }
- break;
+ }break;
}
}
}
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.hpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-02-14 22:32:05 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-02-15 17:58:41 +0000
@@ -6,6 +6,10 @@
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
+// #define FLOW_ENGINE
+
+#ifdef FLOW_ENGINE
+
#pragma once
#include<yade/core/PartialEngine.hpp>
@@ -17,29 +21,16 @@
private:
shared_ptr<TriaxialCompressionEngine> triaxialCompressionEngine;
shared_ptr<CGT::FlowBoundingSphere> flow;
-// FlowBoundingSphere* flow;
- //Tesselation* Tes;
public :
Vector3r gravity;
- bool isActivated;
bool first;
- int PermuteInterval,
- current_state;
- double permeability_factor;
- int previous_state;
- bool currentTes
- ,compute_K
- ,unload
- ,tess_based_force;
- Real loadFactor;
- int cons;
-
+ int current_state
+ ,previous_state
+ ,cons;
Real wall_thickness;
-
- double P_zero;
void Triangulate ( Scene* ncb );
void AddBoundary ( Scene* ncb );
@@ -53,21 +44,23 @@
void NewTriangulation ( Scene* ncb );
void Oedometer_Boundary_Conditions();
- FlowEngine();
virtual ~FlowEngine();
virtual void applyCondition(Scene*);
- REGISTER_ATTRIBUTES(PartialEngine,(isActivated)(first)(currentTes)(P_zero)(PermuteInterval)(compute_K)(permeability_factor)(loadFactor)(unload)(tess_based_force));
-
- protected :
- REGISTER_CLASS_NAME(FlowEngine);
- REGISTER_BASE_CLASS_NAME(PartialEngine);
-
+ YADE_CLASS_BASE_DOC_ATTRS(FlowEngine,PartialEngine,"An engine to solve the flow problem in saturated granular media",
+ ((bool,isActivated,true,"Activates Flow Engine "))
+ ((bool,currentTes,false,"Identifies the current triangulation/tesselation of pore space"))
+ ((double,P_zero,0,"Initial internal pressure for oedometer test"))
+ ((int,PermuteInterval,100000,"Pore space re-triangulation period"))
+ ((bool,compute_K,true,"Activates permeability measure within a granular sample"))
+ ((double,permeability_factor,1.0,"a permability multiplicator"))
+ ((Real,loadFactor,1.5,"Load multiplicator for oedometer test"))
+ ((bool,unload,false,"Remove the load in oedometer test"))
+ ((bool,tess_based_force,true,"true=force computation based on tessalation, false=force computation based on triangulation")));
DECLARE_LOGGER;
};
-
REGISTER_SERIALIZABLE(FlowEngine);
-
+#endif
\ No newline at end of file
