yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07797
[Branch ~yade-dev/yade/trunk] Rev 2887: code cleaning only (including example usage of stringstream for Ema)
------------------------------------------------------------
revno: 2887
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Mon 2011-07-18 15:53:54 +0200
message:
code cleaning only (including example usage of stringstream for Ema)
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
pkg/dem/FlowEngine.cpp
pkg/dem/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 2011-07-07 13:21:17 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2011-07-18 13:53:54 +0000
@@ -488,13 +488,6 @@
if (!cell->vertex(facetVertices[j][y])->info().isFictious) {
cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces + facetNormal*(neighbour_cell->info().p()-cell->info().p())*crossSections[j][y];
}
- ///TEST BEGING
- //Conclusion, we can't easily get ordered vertices from mirror facet...
-// for (int y=0; y<3;y++) {
-// int id1 = cell->vertex(facetVertices[j][y])->info().id();
-// int id2 = neighbour_cell->vertex(facetVertices[Tri.mirror_index(cell,j)][y])->info().id();
-// cerr <<"id1/id2 : "<<id1<<" "<<id2<<endl;}
- ///TEST END
}
}
cell->info().isvisited=!ref;
@@ -1426,38 +1419,7 @@
}
}
-// double FlowBoundingSphere::PressureProfile(const char *filename, Real& time, int& intervals)
-// {
-// RTriangulation& Tri = T[currentTes].Triangulation();
-// vector<double> Pressures;
-//
-// /** CONSOLIDATION CURVES **/
-// Cell_handle permeameter;
-// int n=0, k=0;
-// vector<double> P_ave;
-// std::ofstream consFile(filename, std::ios::out);
-//
-// double Rx = (x_max-x_min) /intervals;
-// double Ry = (y_max-y_min) /intervals;
-// double Rz = (z_max-z_min) /intervals;
-//
-// for (double Y=y_min; Y<=y_max+Ry/10; Y=Y+Ry) {
-// P_ave.push_back(0);
-// for (double X=x_min; X<=x_max+Ry/10; X=X+Rx) {
-// for (double Z=z_min; Z<=z_max+Ry/10; Z=Z+Rz) {
-// P_ave[k]+=Tri.locate(Point(X, Y, Z))->info().p();
-// n++;
-// }
-// }
-// P_ave[k]/= (n);
-// consFile<<k<<" "<<time<<" "<<P_ave[k]<<endl;
-// if (k==intervals/2) Pressures.push_back(P_ave[k]);
-// n=0; k++;
-// }
-// return P_ave[intervals/2];
-// }
-
-double FlowBoundingSphere::PressureProfile(char *filename, Real& time, int& intervals)
+double FlowBoundingSphere::PressureProfile(const char *filename, Real& time, int& intervals)
{
RTriangulation& Tri = T[currentTes].Triangulation();
vector<double> Pressures;
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2011-07-07 13:21:17 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2011-07-18 13:53:54 +0000
@@ -93,8 +93,7 @@
double Permeameter ( double P_Inf, double P_Sup, double Section, double DeltaY, const char *file );
double Sample_Permeability( double& x_Min,double& x_Max ,double& y_Min,double& y_Max,double& z_Min,double& z_Max, string key);
double Compute_HydraulicRadius (Cell_handle cell, int j );
-// double PressureProfile (const char *filename, Real& time, int& intervals );
-double PressureProfile(char *filename, Real& time, int& intervals);
+ double PressureProfile (const char *filename, Real& time, int& intervals );
double dotProduct ( Vecteur x, Vecteur y );
double Compute_EffectiveRadius(Cell_handle cell, int j);
@@ -106,7 +105,6 @@
void ComputeEdgesSurfaces();
Vector3r ComputeViscousForce(Vector3r deltaV, int edge_id);
-// Real ComputeVFacetArea(Finite_edges_iterator ed_it);
RTriangulation& Build_Triangulation ( Real x, Real y, Real z, Real radius, unsigned const id );
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2011-07-07 13:21:17 +0000
+++ pkg/dem/FlowEngine.cpp 2011-07-18 13:53:54 +0000
@@ -100,23 +100,19 @@
Real time = scene->time;
int j = scene->iter;
+ //FIXME: please Ema, put this in a separate function with python wrapper, so we keep action() code simple and usage easier
if (consolidation) {
- char file [50];
mkdir("./Consol", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
- string consol = "./Consol/"+flow->key+"%d_Consol";
-// const char* file = ;
- const char* keyconsol = consol.c_str();
-
- sprintf(file, keyconsol, j);
- char *g = file;
+ stringstream sstr; sstr<<"./Consol/"<<flow->key<<j<<"_Consol";
+ string consol = sstr.str();
timingDeltas->checkpoint("Writing cons_files");
-// MaxPressure = flow->PressureProfile(consol.c_str(), time, intervals);
- MaxPressure = flow->PressureProfile( g, time, intervals);
-
- std::ofstream max_p("pressures.txt", std::ios::app);
+ MaxPressure = flow->PressureProfile(consol.c_str(), time, intervals);
+ static bool consolidationFilesOpened=false;
+ if(!consolidationFilesOpened){
+ std::ofstream max_p("pressures.txt", std::ios::app);
+ std::ofstream settle("settle.txt", std::ios::app);
+ consolidationFilesOpened=True;}
max_p << j << " " << time << " " << MaxPressure << endl;
-
- std::ofstream settle("settle.txt", std::ios::app);
settle << j << " " << time << " " << currentStrain << endl;
}
@@ -213,8 +209,8 @@
double flux=0;
CGT::Cell_handle cell= Tri.locate(flow->imposedP[cond].first);
for (int ngb=0;ngb<4;ngb++) {
- /*if (!cell->neighbor(ngb)->info().Pcondition) */flux+= cell->info().k_norm()[ngb]*(cell->info().p()-cell->neighbor(ngb)->info().p());}
- return flux;
+ /*if (!cell->neighbor(ngb)->info().Pcondition)*/ flux+= cell->info().k_norm()[ngb]*(cell->info().p()-cell->neighbor(ngb)->info().p());}
+ return flux+cell->info().dv();
}
@@ -431,54 +427,25 @@
//AVERAGE CELL VELOCITY
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()) {
- case ( 3 ):
- for ( int g=0;g<4;g++ )
- {
- if ( !cell->vertex ( g )->info().isFictious ) {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
- for (int i=0;i<3;i++) Vel[i] = Vel[i] + sph->state->vel[i]/4;}
- }
- break;
- case ( 2 ):
- for ( int g=0;g<4;g++ )
- {
- if ( !cell->vertex ( g )->info().isFictious ) {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
- for (int i=0;i<3;i++) Vel[i] = Vel[i] + sph->state->vel[i]/4;}
- }
- break;
- case ( 1 ):
- for ( int g=0;g<4;g++ )
- {
- if ( !cell->vertex ( g )->info().isFictious ) {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
- for (int i=0;i<3;i++) Vel[i] = Vel[i] + sph->state->vel[i]/4;}
- }
- break;
- case ( 0 ) :
- for ( int g=0;g<4;g++ )
- {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
- for (int i=0;i<3;i++) Vel[i] = Vel[i] + sph->state->vel[i]/4;}
- }
- break;
-
-
- CGT::RTriangulation& Tri = flow->T[flow->currentTes].Triangulation();
- CGT::Point pos_av_facet;
- double volume_facet_translation = 0;
- CGT::Vecteur Vel_av (Vel[0], Vel[1], Vel[2]);
- for ( int i=0; i<4; i++ ) {
- volume_facet_translation = 0;
- if (!Tri.is_infinite(cell->neighbor(i))){
- CGT::Vecteur Surfk = cell->info()-cell->neighbor(i)->info();
- Real area = sqrt ( Surfk.squared_length() );
- Surfk = Surfk/area;
- CGT::Vecteur branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
- pos_av_facet = (CGT::Point) cell->info() + ( branch*Surfk ) *Surfk;
- volume_facet_translation += Vel_av*cell->info().facetSurfaces[i];
- cell->info().av_vel() = cell->info().av_vel() - volume_facet_translation/cell->info().volume() * ( pos_av_facet-CGAL::ORIGIN );}}
+ for ( int g=0;g<4;g++ ) {
+ if ( !cell->vertex ( g )->info().isFictious ) {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ for (int i=0;i<3;i++) Vel[i] = Vel[i] + sph->state->vel[i]/4;}
+ }
+ CGT::RTriangulation& Tri = flow->T[flow->currentTes].Triangulation();
+ CGT::Point pos_av_facet;
+ double volume_facet_translation = 0;
+ CGT::Vecteur Vel_av (Vel[0], Vel[1], Vel[2]);
+ for ( int i=0; i<4; i++ ) {
+ volume_facet_translation = 0;
+ if (!Tri.is_infinite(cell->neighbor(i))){
+ CGT::Vecteur Surfk = cell->info()-cell->neighbor(i)->info();
+ Real area = sqrt ( Surfk.squared_length() );
+ Surfk = Surfk/area;
+ CGT::Vecteur branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
+ pos_av_facet = (CGT::Point) cell->info() + ( branch*Surfk ) *Surfk;
+ volume_facet_translation += Vel_av*cell->info().facetSurfaces[i];
+ cell->info().av_vel() = cell->info().av_vel() - volume_facet_translation/cell->info().volume() * ( pos_av_facet-CGAL::ORIGIN );}}
}
}
@@ -504,11 +471,8 @@
dVol=cell->info().volumeSign*(newVol - cell->info().volume());
totDVol+=dVol;
eps_vol_max = max(eps_vol_max, abs(dVol/newVol));
-
- cell->info().dv() = (!cell->info().Pcondition)?dVol*invDeltaT:0;
-// cell->info().dv() = dVol*invDeltaT;
+ cell->info().dv() = dVol*invDeltaT;
cell->info().volume() = newVol;
-// if (Debug) cerr<<"v/dv : "<<cell->info().volume()<<" "<<cell->info().dv()<<" ("<<cell->info().fictious()<<")"<<endl;
}
if (Debug) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<endl;
}
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2011-07-07 13:21:17 +0000
+++ pkg/dem/FlowEngine.hpp 2011-07-18 13:53:54 +0000
@@ -32,10 +32,8 @@
Vector3r gravity;
int current_state;
Real wall_thickness;
-// bool Update_Triangulation;
bool currentTes;
int id_offset;
- // double IS;
double wall_up_y, wall_down_y;
double Eps_Vol_Cumulative;
int ReTrg;
@@ -89,7 +87,6 @@
((double, Sinus_Average, 0,,"Pressure value (average) when sinusoidal pressure is applied"))
((bool, CachedForces, true,,"Des/Activate the cached forces calculation"))
((bool, Debug, false,,"Activate debug messages"))
-// ((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"))
@@ -110,9 +107,7 @@
((double, currentStrain, 0,, "Current value of axial strain"))
((int, intervals, 30,, "Number of layers for computation average fluid pressure profiles to build consolidation curves"))
((int, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso"))
-// ((bool, liquefaction, false,,"Measure fluid pressure along the heigth of the sample"))
((double, V_d, 0,,"darcy velocity of fluid in sample"))
-// ((double, bottom_seabed_pressure,0,,"Fluid pressure measured at the bottom of the seabed on the symmetry axe"))
((bool, Flow_imposed_TOP_Boundary, true,, "if false involve pressure imposed condition"))
((bool, Flow_imposed_BOTTOM_Boundary, true,, "if false involve pressure imposed condition"))
((bool, Flow_imposed_FRONT_Boundary, true,, "if false involve pressure imposed condition"))
@@ -146,7 +141,7 @@
,
.def("imposeFlux",&FlowEngine::imposeFlux,(python::arg("pos"),python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
.def("imposePressure",&FlowEngine::imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
- .def("setImposedPressure",&FlowEngine::setImposedPressure,(python::arg("cond"),python::arg("p")),"Set pressure value at the point of index cond.")
+ .def("setImposedPressure",&FlowEngine::setImposedPressure,(python::arg("cond"),python::arg("p")),"Set pressure value at the point indexed 'cond'.")
.def("clearImposedPressure",&FlowEngine::clearImposedPressure,"Clear the list of points with pressure imposed.")
.def("getFlux",&FlowEngine::getFlux,(python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
.def("getConstrictions",&FlowEngine::getConstrictions,"Get the list of constrictions (inscribed circle) for all finite facets.")
