yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05568
[Branch ~yade-dev/yade/trunk] Rev 2413: - Solved retriangulation problems in seabed simulations
------------------------------------------------------------
revno: 2413
committer: ecatalano <ecatalano@dt-rv019>
branch nick: trunk
timestamp: Wed 2010-08-25 17:38:59 +0200
message:
- Solved retriangulation problems in seabed simulations
- Code maintenance
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-07-08 09:11:42 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-08-25 15:38:59 +0000
@@ -323,7 +323,7 @@
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+cos(alpha*M_PI)*(Pressure);}
+ 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));}
}
}
}
@@ -355,7 +355,7 @@
for (int bound=0; bound<6;bound++) {
int& id = *boundsIds[bound];
Tesselation::Vector_Cell tmp_cells;
- tmp_cells.resize(1000);
+ tmp_cells.resize(10000);
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++)
@@ -1789,8 +1789,8 @@
dp_moy = sum_dp/cell2;
j++;
if (j % 1000 == 0) {
-// cout << "pmax " << p_max << "; pmoy : " << p_moy << "; dpmax : " << dp_max << endl;
-// cout << "iteration " << j <<"; erreur : " << dp_max/p_max << 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 )*/);
@@ -1949,9 +1949,9 @@
}
vtkfile.end_cells();
- vtkfile.begin_data("Force",POINT_DATA,SCALARS,FLOAT);
+ vtkfile.begin_data("Force",POINT_DATA,VECTORS,FLOAT);
for (Finite_vertices_iterator v = T.finite_vertices_begin(); v != T.finite_vertices_end(); ++v)
- {if (!v->info().isFictious) vtkfile.write_data((v->info().forces)[1]);}
+ {if (!v->info().isFictious) vtkfile.write_data((v->info().forces)[0],(v->info().forces)[1],(v->info().forces)[2]);}
vtkfile.end_data();
vtkfile.begin_data("Pressure",CELL_DATA,SCALARS,FLOAT);
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-08-16 21:31:08 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-08-25 15:38:59 +0000
@@ -60,7 +60,7 @@
if (!first) flow->GaussSeidel ( );
timingDeltas->checkpoint("Gauss-Seidel");
- if (save_mplot){int j = scene->iter;
+ if (save_mplot){int j = Omega::instance().getCurrentIteration();
char plotfile [50];
mkdir("./mplot", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
string visu_consol = "./mplot/"+flow->key+"%d_Visu_Consol";
@@ -83,25 +83,16 @@
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();
-
-// scene->forces.sync();
-// const Vector3r& fbef = scene->forces.getForce(id);
-// if (id==id_sphere) cout << "force before = " << fbef << endl;
-// if (id==id_sphere) cout << "fluid force = " << V_it->info().forces << endl;
for ( int y=0;y<3;y++ ) f[y] = ( V_it->info().forces ) [y];
scene->forces.addForce ( id, f );
-
-// scene->forces.sync();
-// const Vector3r& faft = scene->forces.getForce(id);
-// if (id==id_sphere) cout << "force after = " << faft << endl;
//scene->forces.addTorque(id,t);
}
timingDeltas->checkpoint("Applying Forces");
- Real time = scene->time;
+ Real time = Omega::instance().getSimulationTime();
- int j = scene->iter;
+ int j = Omega::instance().getCurrentIteration();
char file [50];
mkdir("./Consol", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
string consol = "./Consol/"+flow->key+"%d_Consol";
@@ -118,12 +109,13 @@
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( );}
+
first=false;
-
- if ( scene->iter>1 && scene->iter % PermuteInterval == 0 )
- { Update_Triangulation = true; }
-
- if ( Update_Triangulation ) { Build_Triangulation( );}
+ Update_Triangulation = false;
timingDeltas->checkpoint("Storing Max Pressure");
@@ -175,12 +167,10 @@
flow->T[flow->currentTes].Clear();
flow->T[flow->currentTes].max_id=-1;
flow->x_min = 1000.0, flow->x_max = -10000.0, flow->y_min = 1000.0, flow->y_max = -10000.0, flow->z_min = 1000.0, flow->z_max = -10000.0;
-AddBoundary ( );
+
+ AddBoundary ( );
Triangulate ( );
-
-
-
cout << endl << "Tesselating------" << endl << endl;
flow->T[flow->currentTes].Compute();
@@ -198,21 +188,21 @@
if (compute_K) {flow->TOLERANCE=1e-06; K = flow->Sample_Permeability ( flow->T[flow->currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
BoundaryConditions();
flow->Initialize_pressures( P_zero );
-// flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Pressure, 5);
- BoundaryConditions();
+ if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Pressure, 5);
flow->TOLERANCE=Tolerance;
flow->RELAX=Relax;
}
else
{
cout << "---------UPDATE PERMEABILITY VALUE--------------" << endl;
+ CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
+ for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){cell->info().dv() = 0; cell->info().p() = 0;}
if (compute_K) {flow->TOLERANCE=1e-06; K = flow->Sample_Permeability ( flow->T[flow->currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
- flow->TOLERANCE=Tolerance;
+ BoundaryConditions();
+ flow->Initialize_pressures( P_zero );
flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
- Update_Triangulation=!Update_Triangulation;
- BoundaryConditions();
+ if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Pressure, 5);
}
-
Initialize_volumes ( );
}
@@ -230,7 +220,7 @@
if ( b->shape->getClassIndex() == Sph_Index )
{
Sphere* s=YADE_CAST<Sphere*> ( b->shape.get() );
- const Body::id_t& id = b->getId();
+ const body_id_t& id = b->getId();
Real rad = s->radius;
Real x = b->state->pos[0];
Real y = b->state->pos[1];
@@ -258,17 +248,10 @@
if ( b->shape->getClassIndex() == Bx_Index )
{
Box* w = YADE_CAST<Box*> ( b->shape.get() );
-// const Body::id_t& id = b->getId();
+// 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]);
-
-// 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);
}
}
@@ -299,7 +282,7 @@
if ( b->shape->getClassIndex() == Sph_Index )
{
Sphere* s=YADE_CAST<Sphere*> ( b->shape.get() );
- const Body::id_t& id = b->getId();
+ const body_id_t& id = b->getId();
Real rad = s->radius;
Real x = b->state->pos[0];
Real y = b->state->pos[1];
@@ -307,13 +290,6 @@
flow->T[flow->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;
}
}
@@ -413,8 +389,6 @@
CGT::Vecteur ( v2[0],v2[1],v2[2] ) ) *
flow->boundaries[b].normal ) * ( 0.33333333333* ( V[0][flow->boundaries[b].coordinate]+ V[1][flow->boundaries[b].coordinate]+ V[2][flow->boundaries[b].coordinate] ) - Wall_point[flow->boundaries[b].coordinate] );
-// cout << "volume single cell =" << Volume << endl;
-
return abs ( Volume );
}
@@ -469,8 +443,6 @@
Real Volume = ( CGAL::cross_product ( v1,v2 ) *flow->boundaries[b[0]].normal ) *h;
-// cout << "volume double cell =" << Volume << endl;
-
return abs ( Volume );
}
@@ -512,8 +484,6 @@
Real Volume = ( CGAL::cross_product ( v1,v2 ) ) * h;
-// cout << "volume triple cell =" << Volume << endl;
-
return abs ( Volume );
}
@@ -536,8 +506,6 @@
CGT::Point p4 ( ( A[3] ) [0], ( A[3] ) [1], ( A[3] ) [2] );
Real Volume = ( std::abs ( ( CGT::Tetraedre ( p1,p2,p3,p4 ).volume() ) ) );
-//
-// cout << "volume normal cell =" << Volume << endl;
return abs ( Volume );
}
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.hpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-08-24 12:54:14 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-08-25 15:38:59 +0000
@@ -48,21 +48,22 @@
virtual void action();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(FlowEngine,PartialEngine,"An engine to solve flow problem in saturated granular media",
- ((bool,isActivated,true,,"Activates Flow Engine"))
- ((bool,first,true,,"Controls the initialization/update phases"))
- ((bool,save_vtk,false,,"Enable/disable vtk files creation for visualization"))
- ((bool,save_mplot,false,,"Enable/disable mplot files creation"))
+ ((bool,isActivated,true,"Activates Flow Engine"))
+ ((bool,first,true,"Controls the initialization/update phases"))
+ ((bool,save_vtk,false,"Enable/disable vtk files creation for visualization"))
+ ((bool,save_mplot,false,"Enable/disable mplot files creation"))
((bool, slip_boundary, true, "Controls friction condition on lateral walls"))
-// ((bool,currentTes,false,,"Identifies the current triangulation/tesselation of pore space"))
- ((double,P_zero,0,,"Initial internal pressure for oedometer test"))
- ((double,Tolerance,1e-06,,"Gauss-Seidel Tolerance"))
- ((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,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"))
+ ((bool,WaveAction, false, "Allow sinusoidal pressure condition to simulate ocean waves"))
+// ((bool,currentTes,false,"Identifies the current triangulation/tesselation of pore space"))
+ ((double,P_zero,0,"Initial internal pressure for oedometer test"))
+ ((double,Tolerance,1e-06,"Gauss-Seidel Tolerance"))
+ ((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,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"))
((double, K, 0, "Permeability of the sample"))
((double, MaxPressure, 0, "Maximal value of water pressure within the sample"))
((double, currentStress, 0, "Current value of axial stress"))
