yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #03610
[Branch ~yade-dev/yade/trunk] Rev 2087: - Verified properly re-triangulation
------------------------------------------------------------
revno: 2087
committer: Emanuele Catalano <ecatalano@r2balme>
branch nick: trunk
timestamp: Thu 2010-03-18 16:27:50 +0100
message:
- Verified properly re-triangulation
- Made corrections on Gauss-Seidel break criterion
- Assigned more severe tolerance to permeability computation
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-03-11 14:48:26 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-03-18 15:27:50 +0000
@@ -44,7 +44,7 @@
const int permut4 [4][4] = {{0,1,2,3},{1,2,3,0},{2,3,0,1},{3,0,1,2}};
vector<double> Pressures;
-
+int vtk_infinite_vertices, vtk_infinite_cells;
/*static*/
Point Corner_min;
/*static*/
@@ -271,7 +271,7 @@
}
}
-Tesselation& FlowBoundingSphere::Interpolate(Tesselation& Tes, Tesselation& NewTes)
+void FlowBoundingSphere::Interpolate(Tesselation& Tes, Tesselation& NewTes)
{
Cell_handle old_cell;
@@ -285,7 +285,7 @@
new_cell->info().p() = old_cell->info().p();
// new_cell->info().dv() = old_cell->info().dv();
}
- return NewTes;
+// return NewTes;
}
void FlowBoundingSphere::Localize()
@@ -343,7 +343,8 @@
void FlowBoundingSphere::Compute_Permeability()
{
- RTriangulation& Tri = T[currentTes].Triangulation();
+ cout << "----Computing_Permeability------" << endl;
+ RTriangulation& Tri = T[currentTes].Triangulation();
Vsolid_tot = 0, Vtotalissimo = 0, Vporale = 0, Ssolid_tot = 0;
Finite_cells_iterator cell_end = Tri.finite_cells_end();
@@ -431,6 +432,7 @@
cout<<grains<<"grains - " <<"Vtotale = " << Vtotale << " Vgrains = " << Vgrains << " Vporale1 = " << Vtotale-Vgrains << endl;
cout << "Vtotalissimo = " << Vtotalissimo << " Vsolid_tot = " << Vsolid_tot << " Vporale2 = " << Vporale << " Ssolid_tot = " << Ssolid_tot << endl<< endl;
}
+ cout << "-----Computed_Permeability-----" << endl;
}
void FlowBoundingSphere::DisplayStatistics()
@@ -479,6 +481,9 @@
cout << "There are " << Inferior << " cells INFERIOR." << endl;
cout << "There are " << Superior << " cells SUPERIOR." << endl;
cout << "There are " << Fictious << " cells FICTIOUS." << endl;
+
+ vtk_infinite_vertices = fict;
+ vtk_infinite_cells = Fictious;
}
void FlowBoundingSphere::ComputeTetrahedralForces()
@@ -733,7 +738,6 @@
// #else
else {
Real fsurf = 0.5;//= 0.5/-0.5 for defining the direction of branch vector, "0.5" is for the average pressure
- ///FIXME :all facets computed twice : the forces are x2, it needs a "visited" check before computing anything
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
for (int kk=0; kk<4;kk++) { //Loop on facets -- Compute forces
if (cell->neighbor(kk)->info().isvisited==ref) {
@@ -1515,7 +1519,7 @@
Finite_cells_iterator cell_end = Tri.finite_cells_end();
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++)
- {if (!cell->info().fictious()) cell->info().p() = P_zero;}
+ {cell->info().p() = P_zero;cell->info().dv()=0;}
for (int bound=0; bound<6;bound++) {
int& id = *boundsIds[bound];
@@ -1587,12 +1591,12 @@
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 )*/);
- // cout << "pmax " << p_max << "; pmoy : " << p_moy << endl;
+ } 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;
//Display fluxes?
@@ -1696,31 +1700,38 @@
char filename[80];
sprintf(filename,"out_%d.vtk", number++);
- basicVTKwritter vtkfile((unsigned int) T.number_of_vertices(), (unsigned int) T.number_of_finite_cells());
+ basicVTKwritter vtkfile((unsigned int) T.number_of_vertices()-vtk_infinite_vertices, (unsigned int) T.number_of_finite_cells()-vtk_infinite_cells);
+
vtkfile.open(filename,"test");
vtkfile.begin_vertices();
double x,y,z;
for (Finite_vertices_iterator v = T.finite_vertices_begin(); v != T.finite_vertices_end(); ++v) {
- x = (double)(v->point().point()[0]);
+ if (!v->info().isFictious){
+ x = (double)(v->point().point()[0]);
y = (double)(v->point().point()[1]);
z = (double)(v->point().point()[2]);
- vtkfile.write_point(x,y,z);
+ vtkfile.write_point(x,y,z);}
}
vtkfile.end_vertices();
-
+
vtkfile.begin_cells();
for (Finite_cells_iterator cell = T.finite_cells_begin(); cell != T.finite_cells_end(); ++cell) {
- vtkfile.write_cell(cell->vertex(0)->info().id(), cell->vertex(1)->info().id(), cell->vertex(2)->info().id(), cell->vertex(3)->info().id());
+ if (cell->info().fictious()==0){vtkfile.write_cell(cell->vertex(0)->info().id()-6, cell->vertex(1)->info().id()-6, cell->vertex(2)->info().id()-6, cell->vertex(3)->info().id()-6);}
}
vtkfile.end_cells();
- vtkfile.begin_data("Pressure",CELL_DATA,SCALARS,FLOAT);
-
- for (Finite_cells_iterator cell = T.finite_cells_begin(); cell != T.finite_cells_end(); ++cell) {
- vtkfile.write_data(cell->info().p());
- }
- vtkfile.end_data();
-
+
+ vtkfile.begin_data("Force",POINT_DATA,SCALARS,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]);}
+ vtkfile.end_data();
+
+ vtkfile.begin_data("Pressure",CELL_DATA,SCALARS,FLOAT);
+
+ for (Finite_cells_iterator cell = T.finite_cells_begin(); cell != T.finite_cells_end(); ++cell) {
+ if (cell->info().fictious()==0){vtkfile.write_data(cell->info().p());}
+ }
+ vtkfile.end_data();
}
void FlowBoundingSphere::MGPost(RTriangulation& Tri)
@@ -1861,7 +1872,7 @@
double P_zero = abs((boundary(y_min_id).value-boundary(y_max_id).value)/2);
Initialize_pressures( P_zero );
- GaussSeidel();
+ GaussSeidel();
char *kk;
kk = (char*) key.c_str();
=== modified file 'lib/triangulation/FlowBoundingSphere.h'
--- lib/triangulation/FlowBoundingSphere.h 2010-03-11 14:48:26 +0000
+++ lib/triangulation/FlowBoundingSphere.h 2010-03-18 15:27:50 +0000
@@ -123,7 +123,7 @@
void SliceField ( RTriangulation& Tri );
- Tesselation& Interpolate ( Tesselation& Tes, Tesselation& NewTes );
+ void Interpolate ( Tesselation& Tes, Tesselation& NewTes );
double volume_single_fictious_pore ( Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3, Point PV1 );
//Fast version, assign surface of facet for future forces calculations (pointing from PV2 to PV1)
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-03-11 14:48:26 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-03-18 15:27:50 +0000
@@ -54,9 +54,11 @@
///Compute flow and and forces here
- flow->GaussSeidel ( );
+ if (!first) flow->GaussSeidel ( );
timingDeltas->checkpoint("Gauss-Seidel");
+
+ if (save_vtk) flow->save_vtk_file(flow->T[flow->currentTes].Triangulation());
// flow->MGPost(flow->T[flow->currentTes].Triangulation());
@@ -76,6 +78,7 @@
ncb->forces.addForce ( id, f );
//ncb->forces.addTorque(id,t);
}
+
timingDeltas->checkpoint("Applying Forces");
Real time = Omega::instance().getSimulationTime();
@@ -96,6 +99,8 @@
timingDeltas->checkpoint("Storing Max Pressure");
+ first=false;
+
// }
}
}
@@ -126,14 +131,18 @@
flow->Vtotalissimo=0; flow->Vsolid_tot=0; flow->Vporale=0; flow->Ssolid_tot=0;
flow->SLIP_ON_LATERALS=slip_boundary;
flow->key = triaxialCompressionEngine->Key;
- flow->TOLERANCE=Tolerance;
+ flow->k_factor = permeability_factor;
}
- else
+ else
{
- if (compute_K) {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 );}
+ cout << "---------UPDATE PERMEABILITY VALUE--------------" << endl;
+ if (compute_K) {flow->TOLERANCE=1e-09; 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->currentTes=!flow->currentTes;
+ cout << "--------RETRIANGULATION-----------" << endl;
}
- currentTes=flow->currentTes;
+// currentTes=flow->currentTes;
+ 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 ( ncb );
@@ -144,16 +153,16 @@
flow->Localize ();
- flow->k_factor = permeability_factor;
flow->Compute_Permeability ();
if (first)
{
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) { 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 );}
+ if (compute_K) {flow->TOLERANCE=1e-09; 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 );}
Oedometer_Boundary_Conditions();
- first=!first;
+ flow->Initialize_pressures( P_zero );
+ flow->TOLERANCE=Tolerance;
}
else
{
@@ -161,7 +170,6 @@
Update_Triangulation=!Update_Triangulation;
}
- flow->Initialize_pressures( P_zero );
Initialize_volumes ( ncb );
flow->DisplayStatistics ();
}
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.hpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-03-11 14:48:26 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-03-18 15:27:50 +0000
@@ -43,14 +43,15 @@
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, slip_boundary, false, "Controls friction condition on lateral walls"))
+ ((bool,save_vtk,true,"Enable/disable vtk files creation to visualize pressure and force fields"))
+ ((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"))
((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"))
+ ((Real,loadFactor,1.1,"Load multiplicator for oedometer test"))
((double, K, 0, "Permeability of the sample"))
((bool, ComputeFlow, 1,"If false only triangulation is done, if true flow problem is solved"))
((double, MaxPressure, 0, "Maximal value of water pressure within the sample"))
