yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06211
[Branch ~yade-dev/yade/trunk] Rev 2560: - remove a warning in Tesselation.cpp
------------------------------------------------------------
revno: 2560
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Thu 2010-11-18 21:01:50 +0100
message:
- remove a warning in Tesselation.cpp
- fix FlowBonding sphere for first iteration and a few other things
- make GaussSeidel virtual in FBS
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/Tesselation.cpp
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.cpp'
--- lib/triangulation/FlowBoundingSphere.cpp 2010-11-17 20:07:21 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-11-18 20:01:50 +0000
@@ -32,7 +32,7 @@
#define TESS_BASED_FORCES
#define FACET_BASED_FORCES 1
-#ifdef YADE_OPEN_MP
+#ifdef YADE_OPENMP
// #define GS_OPEN_MP //It should never be defined if Yade is not using openmp
#endif
@@ -92,6 +92,9 @@
OUTPUT_BOUDARIES_RADII = false;
RAVERAGE = false; /** if true use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
}
+
+void FlowBoundingSphere::ResetNetwork() {noCache=true;}
+
Tesselation& FlowBoundingSphere::Compute_Action()
{
return Compute_Action(0,NULL,NULL);
@@ -1036,14 +1039,14 @@
#endif
j++;
// if (j % 100 == 0) {
-// // cout << "pmax " << p_max << "; pmoy : " << p_moy << "; dpmax : " << dp_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 );
// }
#ifdef GS_OPEN_MP
} while (j<1500);
#else
- } while ((dp_max/p_max) > tolerance /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
+ } while ((dp_max/p_max) > tolerance && j<4000 /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
#endif
}
@@ -1171,7 +1174,6 @@
Fictious+=1;
}
}
-
int fict=0, real=0;
for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
if (v->info().isFictious) fict+=1;
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2010-11-17 20:07:21 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2010-11-18 20:01:50 +0000
@@ -64,7 +64,8 @@
void Localize();
void Compute_Permeability();
- void GaussSeidel ( );
+ virtual void GaussSeidel ( );
+ virtual void ResetNetwork();
// void Compute_Forces ();
void Fictious_cells ( );
@@ -167,17 +168,15 @@
// double surface_external_double_fictious(Cell_handle cell, Boundary b);
// double surface_external_single_fictious(Cell_handle cell, Boundary b);
- ///TAUCS
- int SetLinearSystem();
- int SetLinearSystemFullGS();
- void VectorizedGaussSeidel ();
- int TaucsSolveTest();
- int PardisoSolveTest();
- ///END_TAUCS
-
};
} //namespace CGT
+
+#ifdef LINSOLV
+#include "FlowBoundingSphereLinSolv.hpp"
+#endif
+
+
#endif //FLOW_ENGINE
#endif
=== modified file 'lib/triangulation/Tesselation.cpp'
--- lib/triangulation/Tesselation.cpp 2010-11-02 16:23:44 +0000
+++ lib/triangulation/Tesselation.cpp 2010-11-18 20:01:50 +0000
@@ -179,7 +179,7 @@
pass = true;
if ( !Tri->is_infinite ( ( *facet ).first ) && !Tri->is_infinite ( ( *facet ).first->neighbor ( ( *facet ).second ) ) )
{
- cout << "p.x() = " << p.x() << "p.y() = " << p.y() << "p.z() = " << p.z() << endl;
+// cout << "p.x() = " << p.x() << "p.y() = " << p.y() << "p.z() = " << p.z() << endl;
p = ( *facet ).first->info();
( *Coordonnes ) [k++] = p.x(); ( *Coordonnes ) [k++] = p.y(); ( *Coordonnes ) [k++] = p.z();
p = ( *facet ).first->neighbor ( ( *facet ).second )->info();
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2010-11-17 20:07:21 +0000
+++ pkg/dem/FlowEngine.cpp 2010-11-18 20:01:50 +0000
@@ -24,81 +24,78 @@
{
}
-void FlowEngine::action ( )
+void FlowEngine::action()
{
if (!flow) {
- flow = shared_ptr<FlowSolver> (new FlowSolver);
- first=true;Update_Triangulation=false;/*IS=0.f;*/eps_vol_max=0.f;Eps_Vol_Cumulative=0.f;ReTrg=1.0;}
- if ( !isActivated ) return;
+ flow = shared_ptr<FlowSolver> (new FlowSolver);
+ first=true;
+ Update_Triangulation=false;/*IS=0.f;*/
+ eps_vol_max=0.f;
+ Eps_Vol_Cumulative=0.f;
+ ReTrg=1.0;
+ }
+ if (!isActivated) return;
else
{
timingDeltas->start();
- if ( !triaxialCompressionEngine )
- {vector<shared_ptr<Engine> >::iterator itFirst = scene->engines.begin();vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();for ( ;itFirst!=itLast; ++itFirst ){
- if ( ( *itFirst )->getClassName() == "TriaxialCompressionEngine" ){
- cout << "stress controller engine found - FlowEngine" << endl;
- triaxialCompressionEngine = YADE_PTR_CAST<TriaxialCompressionEngine> ( *itFirst );}}
- if ( !triaxialCompressionEngine ) cout << "stress controller engine NOT found" << endl;}
+ if (!triaxialCompressionEngine)
+ {
+ vector<shared_ptr<Engine> >::iterator itFirst = scene->engines.begin();
+ vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();
+ for (;itFirst!=itLast; ++itFirst) {
+ if ((*itFirst)->getClassName() == "TriaxialCompressionEngine") {
+ cout << "stress controller engine found - FlowEngine" << endl;
+ triaxialCompressionEngine = YADE_PTR_CAST<TriaxialCompressionEngine> (*itFirst);}}
+ if (!triaxialCompressionEngine) cout << "stress controller engine NOT found" << endl;
+ }
currentStress = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];
currentStrain = triaxialCompressionEngine->strain[1];
current_state = triaxialCompressionEngine->currentState;
- if ( current_state==3 )
+// if ( current_state==3 )
{
- if ( first ) Build_Triangulation( P_zero );
-
+ if (first) Build_Triangulation(P_zero);
timingDeltas->checkpoint("Triangulating");
- eps_vol_max=0.f;
+ eps_vol_max=0.f;
UpdateVolumes ( );
Eps_Vol_Cumulative += eps_vol_max;
- if (Eps_Vol_Cumulative > ReTrg*EpsVolPercent_RTRG) {Update_Triangulation = true; ReTrg++;}
-
+ if (Eps_Vol_Cumulative > ReTrg*EpsVolPercent_RTRG) {
+ Update_Triangulation = true;
+ ReTrg++;}
timingDeltas->checkpoint("Update_Volumes");
- ///Compute flow and and forces here
- if (!first) flow->GaussSeidel ( );
-
- timingDeltas->checkpoint("Gauss-Seidel");
-
- if (save_mplot) {
- int j = scene->iter;
- char plotfile [50];
- mkdir("./mplot", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
- string visu_consol = "./mplot/"+flow->key+"%d_Visu_Consol";
- const char* key_visu_consol = visu_consol.c_str();
- sprintf(plotfile, key_visu_consol, j);
- char *gg = plotfile;
- flow->mplot(gg);}
-
- if (save_mgpost) flow->MGPost();
-
- if (!CachedForces) flow->ComputeFacetForces();
- else flow->ComputeFacetForcesWithCache();
-
- timingDeltas->checkpoint("Compute_Forces");
-
- ///End Compute flow and forces
- //int number_of_particles = flow->num_particles;
- CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end ();
- Vector3r f; int id;
- //IS = 0.f;
- //std::ofstream Istab ("Stability.txt", std::ios::app);
- 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();
- for ( int y=0;y<3;y++ ) f[y] = ( V_it->info().forces ) [y];
- scene->forces.addForce ( id, f );
- //scene->forces.sync();
- //IS += (scene->forces.getForce(id).norm())/number_of_particles;
+ if (!first) {
+ ///Compute flow and and forces here
+ flow->GaussSeidel();
+ timingDeltas->checkpoint("Gauss-Seidel");
+ if (save_mgpost) flow->MGPost();
+ if (!CachedForces) flow->ComputeFacetForces();
+ else flow->ComputeFacetForcesWithCache();
+ timingDeltas->checkpoint("Compute_Forces");
+
+ ///End Compute flow and forces
+ //int number_of_particles = flow->num_particles;
+ CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
+ Vector3r f;
+ int id;
+ //IS = 0.f;
+ //std::ofstream Istab ("Stability.txt", std::ios::app);
+ 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();
+ for (int y=0;y<3;y++) f[y] = (V_it->info().forces)[y];
+ scene->forces.addForce(id, f);
+ //scene->forces.sync();
+ //IS += (scene->forces.getForce(id).norm())/number_of_particles;
+ }
+ timingDeltas->checkpoint("Applying Forces");
}
//cout << "STABILITY INDEX - IS = " << IS << endl;
//Istab << scene->time << " " << IS << endl;
- timingDeltas->checkpoint("Applying Forces");
-
Real time = scene->time;
int j = scene->iter;
@@ -116,30 +113,28 @@
max_p << j << " " << time << " " << MaxPressure << endl;
std::ofstream settle("settle.txt", std::ios::app);
- settle << j << " " << time << " " << currentStrain << endl;}
-
+ settle << j << " " << time << " " << currentStrain << endl;
+ }
// if ( scene->iter % PermuteInterval == 0 )
// { Update_Triangulation = true; }
- if ( Update_Triangulation && !first) { Build_Triangulation( );}
-
+ if (Update_Triangulation && !first) {
+ Build_Triangulation();
+ Update_Triangulation = false;
+ }
first=false;
- Update_Triangulation = false;
-
- timingDeltas->checkpoint("Storing Max Pressure");
-
-// flow->Average_Grain_Velocity();
+// flow->Average_Grain_Velocity();
if (save_vtk) flow->save_vtk_file();
if (slice_pressures)
{
- char slifile [30];
- mkdir("./Slices", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
- string slice = "./Slices/Slice_"+flow->key+"_%d";
- const char* keyslice = slice.c_str();
- sprintf (slifile, keyslice, j);
- char *f = "slifile";
- flow->SliceField(f);
+ char slifile [30];
+ mkdir("./Slices", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ string slice = "./Slices/Slice_"+flow->key+"_%d";
+ const char* keyslice = slice.c_str();
+ sprintf(slifile, keyslice, j);
+ char *f = "slifile";
+ flow->SliceField(f);
}
// int numero = flow->Average_Cell_Velocity(id_sphere, flow->T[flow->currentTes].Triangulation());
// flow->vtk_average_cell_velocity(flow->T[flow->currentTes].Triangulation(), id_sphere, numero);
@@ -172,20 +167,18 @@
void FlowEngine::Build_Triangulation (double P_zero)
{
- if (first)
- {
- flow->currentTes=0;
- 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->k_factor = permeability_factor;
- flow->DEBUG_OUT = Debug;
- }
- else
- {
+ flow->ResetNetwork();
+ if (first) flow->currentTes=0;
+ else {
flow->currentTes=!flow->currentTes;
- if (Debug) cout << "--------RETRIANGULATION-----------" << endl;
- }
+ if (Debug) cout << "--------RETRIANGULATION-----------" << endl;}
+
+ 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->k_factor = permeability_factor;
+ flow->DEBUG_OUT = Debug;
+
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;
@@ -206,7 +199,7 @@
{
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->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-06; K = flow->Sample_Permeability ( 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 );
if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Pressure, 15);
@@ -221,14 +214,15 @@
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->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 );// FIXME : why, if we are going to interpolate after that?
+ flow->Initialize_pressures(P_zero);// FIXME : why, if we are going to interpolate after that?
flow->TOLERANCE=Tolerance;//So it can be changed at run time
- flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
-
+ flow->Interpolate (flow->T[!flow->currentTes], flow->T[flow->currentTes]);
if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Pressure, 15);
+ else if (TimeBC) flow->ApplySinusoidalPressure_Space_Time(flow->T[flow->currentTes].Triangulation(), Sinus_Pressure, 15, scene->time, scene->dt);
+ flow->TOLERANCE=Tolerance;
+ flow->RELAX=Relax;
}
- Initialize_volumes ( );
- flow->noCache=true;
+ Initialize_volumes();
}
void FlowEngine::AddBoundary ()