yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #03517
[Branch ~yade-dev/yade/trunk] Rev 2057: - Modified the type of some function in FlowBoundingSphere
------------------------------------------------------------
revno: 2057
committer: Emanuele Catalano <ecatalano@r2balme>
branch nick: trunk
timestamp: Mon 2010-03-01 18:55:59 +0100
message:
- Modified the type of some function in FlowBoundingSphere
- Adapted FlowEngine to that modifications
- Adjusted the use of "currentTes" in FlowEngine
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-02-19 11:12:48 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-03-01 17:55:59 +0000
@@ -167,7 +167,7 @@
/** END GAUSS SEIDEL */
char* file ("Permeability");
- Permeameter(Tri, boundary(y_min_id).value, boundary(y_max_id).value, SectionArea, H, file);
+ double ks = Permeameter(Tri, boundary(y_min_id).value, boundary(y_max_id).value, SectionArea, H, file);
clock.top("Permeameter");
Compute_Forces();
@@ -355,8 +355,8 @@
Vecteur n;
- std::ofstream oFile( "Hydraulic_Radius",std::ios::out);
- std::ofstream kFile ( "Permeability" ,std::ios::out );
+// std::ofstream oFile( "Hydraulic_Radius",std::ios::out);
+// std::ofstream kFile ( "Permeability" ,std::ios::out );
Real meanK=0;
Real infiniteK=1e10;
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
@@ -371,9 +371,9 @@
else POS++;
(cell->info().Rh())[j]=Rhv;
(neighbour_cell->info().Rh())[Tri.mirror_index(cell, j)]= Rhv;
- if (DEBUG_OUT)
- oFile << pass << " " << Rhv << " "<<cell->vertex(facetVertices[j][0])->point()
- << " "<<cell->vertex(facetVertices[j][1])->point() << " "<<cell->vertex(facetVertices[j][2])->point()<<endl;
+// if (DEBUG_OUT)
+// oFile << pass << " " << Rhv << " "<<cell->vertex(facetVertices[j][0])->point()
+// << " "<<cell->vertex(facetVertices[j][1])->point() << " "<<cell->vertex(facetVertices[j][2])->point()<<endl;
Vecteur l = p1 - p2;
distance = sqrt(l.squared_length());
n = l/distance;
@@ -1608,7 +1608,7 @@
}
-void FlowBoundingSphere::Permeameter(RTriangulation &Tri, double P_Inf, double P_Sup, double Section, double DeltaY, char *file)
+double FlowBoundingSphere::Permeameter(RTriangulation &Tri, double P_Inf, double P_Sup, double Section, double DeltaY, char *file)
{
std::ofstream kFile(file, std::ios::out);
@@ -1686,6 +1686,8 @@
cout << "The permeability of the sample is = " << k << " m^2" <<endl;
kFile << "The permeability of the sample is = " << k << " m^2" <<endl;
// cout << "The Darcy permeability of the sample is = " << k_darcy/0.987e-12 << " darcys" << endl;
+
+ return Ks;
}
void FlowBoundingSphere::save_vtk_file(RTriangulation &T)
@@ -1801,7 +1803,7 @@
}
}
-void FlowBoundingSphere::PermeameterCurve(RTriangulation& Tri, char *filename, Real time)
+double FlowBoundingSphere::PermeameterCurve(RTriangulation& Tri, char *filename, Real time)
{
/** CONSOLIDATION CURVES **/
Cell_handle permeameter;
@@ -1839,15 +1841,13 @@
// cout << "P_ave[" << k << "] = " << P_ave[k]/ ( m+n+o ) << " n = " << n << " m = " << m << " o = " << o << endl;
P_ave[k]/= (m+n+o);
consFile<<k<<" "<<time<<" "<<P_ave[k]<<endl;
- if (k==15) Pressures.push_back(P_ave[k]);
- n=0;
- m=0;
- o=0;
- k++;
- }
+ if (k==intervals/2) Pressures.push_back(P_ave[k]);
+ n=0, m=0, o=0; k++;
+ }
+ return P_ave[intervals/2];
}
-void FlowBoundingSphere::Sample_Permeability(RTriangulation& Tri, double x_Min,double x_Max ,double y_Min,double y_Max,double z_Min,double z_Max, string key)
+double FlowBoundingSphere::Sample_Permeability(RTriangulation& Tri, double x_Min,double x_Max ,double y_Min,double y_Max,double z_Min,double z_Max, string key)
{
double Section = (x_Max-x_Min) * (z_Max-z_Min);
@@ -1865,7 +1865,7 @@
char *kk;
kk = (char*) key.c_str();
- Permeameter(Tri, boundary(y_min_id).value, boundary(y_max_id).value, Section, DeltaY, kk);
+ return 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)
=== modified file 'lib/triangulation/FlowBoundingSphere.h'
--- lib/triangulation/FlowBoundingSphere.h 2010-02-19 11:12:48 +0000
+++ lib/triangulation/FlowBoundingSphere.h 2010-03-01 17:55:59 +0000
@@ -84,10 +84,10 @@
void Dessine_Triangulation ( Vue3D &Vue, RTriangulation &T );
void Dessine_Short_Tesselation ( Vue3D &Vue, Tesselation &Tes );
#endif
- void Permeameter ( RTriangulation& Tri, double P_Inf, double P_Sup, double Section, double DeltaY, char *file );
- void Sample_Permeability ( RTriangulation& Tri, double x_Min,double x_Max ,double y_Min,double y_Max,double z_Min,double z_Max, std::string key );
+ double Permeameter ( RTriangulation& Tri, double P_Inf, double P_Sup, double Section, double DeltaY, char *file );
+ double Sample_Permeability ( RTriangulation& Tri, double x_Min,double x_Max ,double y_Min,double y_Max,double z_Min,double z_Max, std::string key );
double Compute_HydraulicRadius ( RTriangulation& Tri, Cell_handle cell, int j );
- void PermeameterCurve ( RTriangulation& Tri, char *filename, Real time );
+ double PermeameterCurve ( RTriangulation& Tri, char *filename, Real time );
double dotProduct ( Vecteur x, Vecteur y );
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-02-19 13:31:48 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-03-01 17:55:59 +0000
@@ -22,6 +22,12 @@
FlowEngine::~FlowEngine()
{
}
+
+std::ofstream plot ("pressurees.txt", std::ios::out);
+// std::ofstream cons_NONDAMP ("cons_NONDAMP_SLIP", std::ios::out);
+// std::ofstream settle_DAMP ("settle_DAMP_SLIP", std::ios::out);
+// std::ofstream settle_NONDAMP ("settle_NONDAMP_SLIP", std::ios::out);
+
void FlowEngine::applyCondition ( Scene* ncb )
{
if (!flow) {flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);first=true;}
@@ -44,6 +50,7 @@
if ( !triaxialCompressionEngine ) cout << "stress controller engine NOT found" << endl;
}
+ currentStress = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];
current_state = triaxialCompressionEngine->currentState;
if ( !first && current_state==3 )
@@ -58,19 +65,19 @@
flow->GaussSeidel ( );
- flow->MGPost(flow->T[currentTes].Triangulation());
+// flow->MGPost(flow->T[flow->currentTes].Triangulation());
flow->tess_based_force=tess_based_force;
flow->Compute_Forces ( );
///End Compute flow and forces
- CGT::Finite_vertices_iterator vertices_end = flow->T[currentTes].Triangulation().finite_vertices_end ();
+ CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end ();
Vector3r f;
int id;
- for ( CGT::Finite_vertices_iterator V_it = flow->T[currentTes].Triangulation().finite_vertices_begin (); V_it != vertices_end; V_it++ )
+ 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];
@@ -89,11 +96,17 @@
sprintf (file, keyconsol, j);
char *g = file;
- flow->PermeameterCurve(flow->T[currentTes].Triangulation(), g, time);
+ MaxPressure = flow->PermeameterCurve(flow->T[flow->currentTes].Triangulation(), g, time);
+ plot << j << " " << MaxPressure << endl;
+
+// if (damped) {cons_DAMP << j << " " << time << " " << flow->Pressures[cons] << endl; cons++;}
+// if (!damped){cons_NONDAMP << j << " " << time << " " << flow->Pressures[cons] << endl; cons++;}
+// if (damped) {settle_DAMP << j << " " << time << " " << triaxialCompressionEngine->uniaxialEpsilonCurr << endl;}
+// if (!damped) {settle_NONDAMP << j << " " << time << " " << triaxialCompressionEngine->uniaxialEpsilonCurr << endl;}
if ( Omega::instance().getCurrentIteration() % PermuteInterval == 0 )
{
-// 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 );
+// 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 );
cout << endl << "---NEW TRIANGULATION---" << endl << endl;
@@ -105,10 +118,10 @@
std::ofstream PFile ( "NewTriangulation_Pressures",std::ios::out );
- CGT::Finite_cells_iterator cell_end = flow->T[currentTes].Triangulation().finite_cells_end();
+ CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
int j=0;
- for ( CGT::Finite_cells_iterator cell = flow->T[currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ )
+ for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ )
{
PFile << ++j << " " << cell->info().p() << endl;
}
@@ -129,9 +142,9 @@
flow->DisplayStatistics ();
- CGT::Finite_cells_iterator cell_end = flow->T[currentTes].Triangulation().finite_cells_end();
+ CGT::Finite_cells_iterator cell_end = flow->T[flow->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++ )
+ 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;
y++;
@@ -140,7 +153,7 @@
flow->key = triaxialCompressionEngine->Key;
- 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 );}
+ 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 );}
Oedometer_Boundary_Conditions();
flow->Initialize_pressures();
@@ -175,8 +188,8 @@
void FlowEngine::Initialize ( Scene* ncb, double P_zero )
{
- currentTes=0;
- flow->currentTes=currentTes;
+ flow->currentTes=0;
+ currentTes=flow->currentTes;
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;
@@ -184,9 +197,7 @@
Triangulate ( ncb );
cout << endl << "Tesselating------" << endl << endl;
- flow->T[currentTes].Compute();
-
-// flow->Fictious_cells();
+ flow->T[flow->currentTes].Compute();
flow->Localize ();
@@ -244,7 +255,7 @@
Real y = b->state->pos[1];
Real z = b->state->pos[2];
- flow->T[currentTes].insert(x, y, z, rad, id);
+ flow->T[flow->currentTes].insert(x, y, z, rad, id);
contator+=1;
}
@@ -267,18 +278,17 @@
{
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;
-// flow->currentTes=!currentTes;
- currentTes=!currentTes;
- flow->currentTes=currentTes;
+ flow->currentTes=!flow->currentTes;
+ currentTes = flow->currentTes;
- flow->T[currentTes].Clear();
+ flow->T[flow->currentTes].Clear();
Triangulate ( ncb );
AddBoundary ( ncb );
// flow->Tesselate();
- flow->T[currentTes].Compute();
+ flow->T[flow->currentTes].Compute();
// flow->Fictious_cells();
@@ -290,18 +300,18 @@
// flow->Sample_Permeability ( t2.Triangulation(), x_min, x_max, z_min, z_max, y_max, y_min );
- flow->Interpolate ( flow->T[!currentTes], flow->T[currentTes] );
-
-// currentTes = !currentTes;
-
-// flow->T[currentTes] = flow->T[currentTes];
-// t2 = flow->T[!currentTes];
+ flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
+
+// flow->currentTes = !flow->currentTes;
+
+// flow->T[flow->currentTes] = flow->T[flow->currentTes];
+// t2 = flow->T[!flow->currentTes];
}
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++ )
+ 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++ )
{
switch ( cell->info().fictious() )
{
@@ -319,8 +329,8 @@
cout << "Updating volumes.............." << endl;
Real deltaT = ncb->dt;
- 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++ )
+ 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++ )
{
switch ( cell->info().fictious() )
{
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.hpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-02-19 13:31:48 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-03-01 17:55:59 +0000
@@ -58,6 +58,10 @@
((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"))
+ ((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"))
+ ((double, currentStress, 0, "Current value of normal stress"))
((bool,tess_based_force,true,"true=force computation based on tessalation, false=force computation based on triangulation")));
DECLARE_LOGGER;
};
