yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #08362
[Branch ~yade-dev/yade/trunk] Rev 3038: - fix the post processing issue of pressure field extracted from the new (uninitialized) triangul...
------------------------------------------------------------
revno: 3038
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Tue 2012-03-06 12:49:48 +0100
message:
- fix the post processing issue of pressure field extracted from the new (uninitialized) triangulation after re-triangulation
(will be mirrored to git after some verifications on gitHub)
modified:
lib/triangulation/FlowBoundingSphere.ipp
lib/triangulation/PeriodicFlow.cpp
lib/triangulation/def_types.h
pkg/dem/FlowEngine.cpp
--
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.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2012-02-14 16:58:21 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2012-03-06 11:49:48 +0000
@@ -284,7 +284,7 @@
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::Average_Relative_Cell_Velocity()
{
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
Point pos_av_facet;
int num_cells = 0;
double facet_flow_rate = 0;
@@ -329,7 +329,7 @@
void FlowBoundingSphere<Tesselation>::Average_Fluid_Velocity()
{
Average_Relative_Cell_Velocity();
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
int num_vertex = 0;
Finite_vertices_iterator vertices_end = Tri.finite_vertices_end();
for (Finite_vertices_iterator V_it = Tri.finite_vertices_begin(); V_it != vertices_end; V_it++) {
@@ -375,7 +375,7 @@
vector<Real> FlowBoundingSphere<Tesselation>::Average_Fluid_Velocity_On_Sphere(unsigned int Id_sph)
{
Average_Relative_Cell_Velocity();
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
Real Volumes; CGT::Vecteur VelocityVolumes;
vector<Real> result;
@@ -399,14 +399,14 @@
template <class Tesselation>
double FlowBoundingSphere<Tesselation>::MeasurePorePressure (double X, double Y, double Z)
{
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
Cell_handle cell = Tri.locate(Point(X,Y,Z));
return cell->info().p();
}
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::MeasurePressureProfile(double Wall_up_y, double Wall_down_y)
{
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
Cell_handle permeameter;
std::ofstream capture ("Pressure_profile", std::ios::app);
int intervals = 5;
@@ -977,7 +977,7 @@
Finite_cells_iterator cell_end = Tri.finite_cells_end();
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++){
- cell->info().p() = P_zero;cell->info().dv()=0;}
+ cell->info().p() = P_zero; cell->info().dv()=0;}
for (int bound=0; bound<6;bound++) {
int& id = *boundsIds[bound];
@@ -1277,19 +1277,25 @@
cout << "There are " << Inside << " cells INSIDE." << endl;
cout << "There are " << Fictious << " cells FICTIOUS." << endl;}
- vtk_infinite_vertices = fict;
- vtk_infinite_cells = Fictious;
num_particles = real;
}
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::saveVtk()
{
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
static unsigned int number = 0;
char filename[80];
mkdir("./VTK", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
sprintf(filename,"./VTK/out_%d.vtk", number++);
+ //count fictious vertices and cells
+ vtk_infinite_vertices=vtk_infinite_cells=0;
+ 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()) vtk_infinite_cells+=1;
+ for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v)
+ if (v->info().isFictious) vtk_infinite_vertices+=1;
+
basicVTKwritter vtkfile((unsigned int) Tri.number_of_vertices()-vtk_infinite_vertices, (unsigned int) Tri.number_of_finite_cells()-vtk_infinite_cells);
vtkfile.open(filename,"test");
@@ -1332,7 +1338,7 @@
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::MGPost()
{
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
Point P;
ofstream file("mgp.out.001");
@@ -1416,7 +1422,7 @@
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::mplot (char *filename)
{
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
std::ofstream plot (filename, std::ios::out);
Cell_handle permeameter;
int intervals = 30;
@@ -1464,7 +1470,7 @@
void FlowBoundingSphere<Tesselation>::SliceField(const char *filename)
{
/** Pressure field along one cutting plane **/
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
Cell_handle permeameter;
std::ofstream consFile(filename,std::ios::out);
@@ -1497,7 +1503,7 @@
//Compute av. velocity first, because in the following permeabilities will be overwritten with "junk" (in fact velocities from comsol)
Average_Relative_Cell_Velocity();
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
Cell_handle c;
ifstream loadFile("vx_grid_03_07_ns.txt");
ifstream loadFileY("vy_grid_03_07_ns.txt");
=== modified file 'lib/triangulation/PeriodicFlow.cpp'
--- lib/triangulation/PeriodicFlow.cpp 2012-02-29 12:46:28 +0000
+++ lib/triangulation/PeriodicFlow.cpp 2012-03-06 11:49:48 +0000
@@ -337,7 +337,7 @@
Finite_cells_iterator cell_end = Tri.finite_cells_end();
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++){
- cell->info().setP(0); P_zero;cell->info().dv()=0;}
+ cell->info().setP(P_zero); cell->info().dv()=0;}
for (int bound=0; bound<6;bound++) {
int& id = *boundsIds[bound];
=== modified file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h 2012-02-29 12:46:28 +0000
+++ lib/triangulation/def_types.h 2012-03-06 11:49:48 +0000
@@ -119,6 +119,8 @@
index=0;
volumeSign=0;
s=0;
+ VolumeVariation=0;
+ pression=0;
}
double inv_sum_k;
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2012-02-29 12:46:28 +0000
+++ pkg/dem/FlowEngine.cpp 2012-03-06 11:49:48 +0000
@@ -40,44 +40,40 @@
if ( first ) Build_Triangulation ( P_zero,solver );
timingDeltas->checkpoint ( "Triangulating" );
UpdateVolumes ( solver );
- if ( !first ) {
-// eps_vol_max=0.f;//huh? in that case Eps_Vol_Cumulative will always be zero
- 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
- solver->GaussSeidel();
- timingDeltas->checkpoint ( "Gauss-Seidel" );
- if ( save_mgpost ) solver->MGPost();
- if ( !CachedForces ) solver->ComputeFacetForces();
- else solver->ComputeFacetForcesWithCache();
- timingDeltas->checkpoint ( "Compute_Forces" );
-
- ///Application of vicscous forces
- scene->forces.sync();
- if ( viscousShear ) ApplyViscousForces ( solver );
-
- Finite_vertices_iterator vertices_end = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
- for ( Finite_vertices_iterator V_it = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); V_it != vertices_end; V_it++ ) {
- if ( !viscousShear )
- scene->forces.addForce ( V_it->info().id(), Vector3r ( ( V_it->info().forces ) [0],V_it->info().forces[1],V_it->info().forces[2] ) );
- else
- scene->forces.addForce ( V_it->info().id(), Vector3r ( ( V_it->info().forces ) [0],V_it->info().forces[1],V_it->info().forces[2] ) +solver->viscousShearForces[V_it->info().id() ] );
- }
- ///End Compute flow and forces
- timingDeltas->checkpoint ( "Applying Forces" );
+
+ 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
+ solver->GaussSeidel();
+ timingDeltas->checkpoint ( "Gauss-Seidel" );
+ if ( save_mgpost ) solver->MGPost();
+ if ( !CachedForces ) solver->ComputeFacetForces();
+ else solver->ComputeFacetForcesWithCache();
+ timingDeltas->checkpoint ( "Compute_Forces" );
+
+ ///Application of vicscous forces
+ scene->forces.sync();
+ if ( viscousShear ) ApplyViscousForces ( solver );
+
+ Finite_vertices_iterator vertices_end = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
+ for ( Finite_vertices_iterator V_it = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); V_it != vertices_end; V_it++ ) {
+ if ( !viscousShear )
+ scene->forces.addForce ( V_it->info().id(), Vector3r ( ( V_it->info().forces ) [0],V_it->info().forces[1],V_it->info().forces[2] ) );
+ else
+ scene->forces.addForce ( V_it->info().id(), Vector3r ( ( V_it->info().forces ) [0],V_it->info().forces[1],V_it->info().forces[2] ) +solver->viscousShearForces[V_it->info().id() ] );
}
- if ( scene->iter % PermuteInterval == 0 )
- { Update_Triangulation = true; }
+ ///End Compute flow and forces
+ timingDeltas->checkpoint ( "Applying Forces" );
if ( Update_Triangulation && !first ) {
- Build_Triangulation ( solver );
+ Build_Triangulation ( P_zero, solver );
Update_Triangulation = false;
}
@@ -85,9 +81,6 @@
first=false;
timingDeltas->checkpoint ( "Ending" );
-
-// if(id_sphere>=0) flow->Average_Fluid_Velocity_On_Sphere(id_sphere);
-// }
}
template<class Solver>
@@ -133,6 +126,7 @@
//force immediate update of boundary conditions
Update_Triangulation=true;
}
+
template<class Solver>
void FlowEngine::imposeFlux ( Vector3r pos, Real flux,Solver& flow )
{
@@ -190,7 +184,12 @@
flow->useSolver = useSolver;
flow->VISCOSITY = viscosity;
flow->areaR2Permeability=areaR2Permeability;
-// flow->key = key;
+ flow->TOLERANCE=Tolerance;
+ flow->RELAX=Relax;
+ flow->meanK_LIMIT = meanK_correction;
+ flow->meanK_STAT = meanK_opt;
+ flow->permeability_map = permeability_map;
+ //flow->key = key;
flow->T[flow->currentTes].Clear();
flow->T[flow->currentTes].max_id=-1;
@@ -203,21 +202,11 @@
flow->Define_fictious_cells();
flow->DisplayStatistics ();
-
- flow->meanK_LIMIT = meanK_correction;
- flow->meanK_STAT = meanK_opt;
- flow->permeability_map = permeability_map;
flow->Compute_Permeability ( );
porosity = flow->V_porale_porosity/flow->V_totale_porosity;
- flow->TOLERANCE=Tolerance;flow->RELAX=Relax;
- typedef typename Solver::element_type Flow;
- typedef typename Flow::Finite_cells_iterator Finite_cells_iterator;
- Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
- for ( 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 ) {BoundaryConditions ( flow ); K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max );}
+ if ( compute_K ) {BoundaryConditions ( flow ); K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max );}
BoundaryConditions ( flow );
flow->Initialize_pressures ( P_zero );
@@ -685,15 +674,13 @@
Real PeriodicFlowEngine::Volume_cell ( Cell_handle cell )
{
Real volume = CGT::Tetraedre (
- makeCgPoint ( Body::byId ( cell->vertex ( 0 )->info().id(), scene )->state->pos ) + cell->vertex ( 0 )->info().ghostShift(),
- makeCgPoint ( Body::byId ( cell->vertex ( 1 )->info().id(), scene )->state->pos ) + cell->vertex ( 1 )->info().ghostShift(),
- makeCgPoint ( Body::byId ( cell->vertex ( 2 )->info().id(), scene )->state->pos ) + cell->vertex ( 2 )->info().ghostShift(),
- makeCgPoint ( Body::byId ( cell->vertex ( 3 )->info().id(), scene )->state->pos ) + cell->vertex ( 3 )->info().ghostShift() )
- .volume();
-
+ makeCgPoint ( Body::byId ( cell->vertex ( 0 )->info().id(), scene )->state->pos ) + cell->vertex ( 0 )->info().ghostShift(),
+ makeCgPoint ( Body::byId ( cell->vertex ( 1 )->info().id(), scene )->state->pos ) + cell->vertex ( 1 )->info().ghostShift(),
+ makeCgPoint ( Body::byId ( cell->vertex ( 2 )->info().id(), scene )->state->pos ) + cell->vertex ( 2 )->info().ghostShift(),
+ makeCgPoint ( Body::byId ( cell->vertex ( 3 )->info().id(), scene )->state->pos ) + cell->vertex ( 3 )->info().ghostShift() )
+ .volume();
if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
return volume;
-
}
Real PeriodicFlowEngine::Volume_cell_single_fictious ( Cell_handle cell )
@@ -718,7 +705,6 @@
}
double v1[3], v2[3];
for ( int g=0;g<3;g++ ) { v1[g]=V[0][g]-V[1][g]; v2[g]=V[0][g]-V[2][g];}
-
Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
return abs ( Volume );
}
@@ -946,4 +932,4 @@
#endif //FLOW_ENGINE
#endif /* YADE_CGAL */
-// kate: indent-mode cstyle; space-indent on; indent-width 8;
+
