yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06192
[Branch ~yade-dev/yade/trunk] Rev 2555: - Retriangulation controlled by a volumetric deformation threshold
------------------------------------------------------------
revno: 2555
committer: ecatalano <ecatalano@dt-rv019>
branch nick: yade
timestamp: Tue 2010-11-16 14:10:35 +0100
message:
- Retriangulation controlled by a volumetric deformation threshold
- Substituted pkg-dem in pkg/dem (and similar) to included files
- Removed #define flow_engine line from def_type, definition is moved to the scons profile
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.h
lib/triangulation/Network.cpp
lib/triangulation/Network.h
lib/triangulation/def_types.h
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 2010-10-20 11:28:25 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-11-16 13:10:35 +0000
@@ -21,6 +21,7 @@
#include "Network.h"
+
// #define XVIEW
#include "FlowBoundingSphere.h"//include after #define XVIEW
@@ -63,17 +64,17 @@
FlowBoundingSphere::FlowBoundingSphere()
{
- x_min = 1000.0, x_max = -10000.0, y_min = 1000.0, y_max = -10000.0, z_min = 1000.0, z_max = -10000.0;
- currentTes = 0;
- nOfSpheres = 0;
- fictious_vertex = 0;
- SectionArea = 0, Height=0, Vtotale=0;
- vtk_infinite_vertices=0, vtk_infinite_cells=0;
-
- tess_based_force = true;
- for (int i=0;i<6;i++) boundsIds[i] = 0;
- minPermLength=-1;
- SLIP_ON_LATERALS = false;//no-slip/symmetry conditions on lateral boundaries
+ x_min = 1000.0, x_max = -10000.0, y_min = 1000.0, y_max = -10000.0, z_min = 1000.0, z_max = -10000.0;
+ currentTes = 0;
+ nOfSpheres = 0;
+ fictious_vertex = 0;
+ SectionArea = 0, Height=0, Vtotale=0;
+ vtk_infinite_vertices=0, vtk_infinite_cells=0;
+
+ tess_based_force = true;
+ for (int i=0;i<6;i++) boundsIds[i] = 0;
+ minPermLength=-1;
+ SLIP_ON_LATERALS = true;//no-slip/symmetry conditions on lateral boundaries
TOLERANCE = 1e-06;
RELAX = 1.9;
ks=0;
@@ -82,11 +83,11 @@
meanK_STAT = false; K_opt_factor=0;
noCache=true;
computeAllCells=true;//might be turned false IF the code is reorganized (we can make a separate function to compute unitForceVectors outside Compute_Permeability) AND it really matters for CPU time
- DEBUG_OUT = true;
+ DEBUG_OUT = false;
RAVERAGE = false; /** if true use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
}
-Tesselation& FlowBoundingSphere::Compute_Action()
+void FlowBoundingSphere::Compute_Action()
{
Compute_Action(0,NULL,NULL);
}
@@ -554,6 +555,32 @@
}
}
+void FlowBoundingSphere::ApplySinusoidalPressure_Space_Time(RTriangulation& Tri, double Pressure, double load_intervals, double time, double dt)
+{
+ //FIXME : rivedere!!
+
+ double step = 1/load_intervals;
+ Tesselation::Vector_Cell tmp_cells;
+ tmp_cells.resize(1000);
+ Tesselation::VCell_iterator cells_it = tmp_cells.begin();
+ for (double alpha=0; alpha<1.001; alpha+=step)
+ {
+ Tesselation::VCell_iterator cells_end = Tri.incident_cells(T[currentTes].vertexHandles[y_max_id],cells_it);
+ for (Tesselation::VCell_iterator it = tmp_cells.begin(); it != cells_end; it++)
+ {
+ if(!Tri.is_infinite(*it)){
+ Point& p1 = (*it)->info();
+ Cell_handle& cell = *it;
+ if (p1.x()>(alpha*(x_max-x_min)) && p1.x()<((alpha+step)*(x_max-x_min)))
+ {
+ if (alpha<0.5) cell->info().p() = (Pressure/2)*(1+cos(alpha*M_PI)-(1-cos(time/(20*dt)))*M_PI);
+ if (alpha>0.5) cell->info().p() = (Pressure/2)*(1+cos(alpha*M_PI)+(1-cos(time/(20*dt)))*M_PI);
+ }
+ }
+ }
+ }
+}
+
void FlowBoundingSphere::Interpolate(Tesselation& Tes, Tesselation& NewTes)
{
Cell_handle old_cell;
@@ -572,7 +599,7 @@
void FlowBoundingSphere::Compute_Permeability()
{
- if (DEBUG_OUT) cout << "----Computing_Permeability------" << endl;
+ if (DEBUG_OUT) 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();
@@ -586,9 +613,9 @@
Vecteur n;
- std::ofstream oFile( "Radii",std::ios::out);
- std::ofstream fFile( "Radii_Fictious",std::ios::out);
- std::ofstream kFile ( "LocalPermeabilities" ,std::ios::app );
+// std::ofstream oFile( "Radii",std::ios::out);
+// std::ofstream fFile( "Radii_Fictious",std::ios::out);
+// std::ofstream kFile ( "LocalPermeabilities" ,std::ios::app );
Real meanK=0, STDEV=0;
Real infiniteK=1e10;
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
@@ -613,6 +640,7 @@
W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0-v2)*(v1-v2)/sqrt((v1-v2).squared_length()*(v2-v0).squared_length())));
#endif
pass+=1;
+// cerr << "test1" << endl;
Rhv = Compute_HydraulicRadius(cell, j);
if (Rhv<0) NEG++;
else POS++;
@@ -623,10 +651,10 @@
// if ((reff+requiv)/2 > 2*Rhv) {cout << "Here it is ---------- " << endl << "p1 = <" << cell->info() << ">; p2 = <" << neighbour_cell->info() << ">;" << endl << "Rav = " << (reff+requiv)/2 << " Rhydr = " << 2*Rhv << endl << endl << "Vpore = " << Vpore << " Ssolid = " << Ssolid << " cell_fictious = " << cell->info().fictious() << " neighbour_cell_fictious = " << neighbour_cell->info().fictious() << endl << endl << "facet vertices = <" << cell->vertex(facetVertices[j][0])->point() << ">;,<" << cell->vertex(facetVertices[j][1])->point() << ">;,<" << cell->vertex(facetVertices[j][2])->point() << ">;" << endl << endl;}
// else cout << "NORMAL CASE" << endl << "p1 = <" << cell->info() << ">; p2 = <" << neighbour_cell->info() << ">;" << endl << "Rav = " << (reff+requiv)/2 << " Rhydr = " << 2*Rhv << endl << endl << "Vpore = " << Vpore << " Ssolid = " << Ssolid << " cell_fictious = " << cell->info().fictious() << " neighbour_cell_fictious = " << neighbour_cell->info().fictious() << endl << endl << "facet vertices = <" << cell->vertex(facetVertices[j][0])->point() << ">;,<" << cell->vertex(facetVertices[j][1])->point() << ">;,<" << cell->vertex(facetVertices[j][2])->point() << ">;" << endl << endl;
- if (cell->info().fictious()==0){
- oFile << 2*Rhv << " ";
- oFile << (reff+requiv)/2 << endl;}
- else {fFile << 2*Rhv << " ";fFile << (reff+requiv)/2 << endl;}
+// if (cell->info().fictious()==0){
+// oFile << 2*Rhv << " ";
+// oFile << (reff+requiv)/2 << endl;}
+// else {fFile << 2*Rhv << " ";fFile << (reff+requiv)/2 << endl;}
(cell->info().Rh())[j]=Rhv;
(neighbour_cell->info().Rh())[Tri.mirror_index(cell, j)]= Rhv;
@@ -655,7 +683,7 @@
meanK += (cell->info().k_norm())[j];
- if (!meanK_LIMIT) kFile << ( cell->info().k_norm() )[j] << endl;
+// if (!meanK_LIMIT) kFile << ( cell->info().k_norm() )[j] << endl;
// (neighbour_cell->info().facetSurf())[Tri.mirror_index(cell, j)]= (-k) *n;
}
// else if ( Tri.is_infinite ( neighbour_cell )) k = 0;//connection to an infinite cells
@@ -667,9 +695,9 @@
if (DEBUG_OUT) cout << "PassCompK = " << pass << endl;
if (meanK_LIMIT) {
- cout << "meanK = " << meanK << endl;
+ if (DEBUG_OUT) cout << "meanK = " << meanK << endl;
Real maxKdivKmean = MAXK_DIV_KMEAN;
- cout << "maxKdivKmean = " << maxKdivKmean << endl;
+ if (DEBUG_OUT) cout << "maxKdivKmean = " << maxKdivKmean << endl;
ref = Tri.finite_cells_begin()->info().isvisited; pass=0;
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
for (int j=0; j<4; j++) {
@@ -678,7 +706,7 @@
pass++;
(cell->info().k_norm())[j] = min((cell->info().k_norm())[j], maxKdivKmean*meanK);
(neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]=(cell->info().k_norm())[j];
- kFile << (cell->info().k_norm())[j] << endl;
+// kFile << (cell->info().k_norm())[j] << endl;
}
}cell->info().isvisited = !ref;
}
@@ -726,12 +754,12 @@
if (DEBUG_OUT) {
cout<<grains<<"grains - " <<"Vtotale = " << 2*Vtotale << " Vgrains = " << 2*Vgrains << " Vporale1 = " << 2*(Vtotale-Vgrains) << endl;
cout << "Vtotalissimo = " << Vtotalissimo << " Vsolid_tot = " << Vsolid_tot << " Vporale2 = " << Vporale << " Ssolid_tot = " << Ssolid_tot << endl<< endl;
- }
+
if (!RAVERAGE) cout << "------Hydraulic Radius is used for permeability computation------" << endl << endl;
else cout << "------Average Radius is used for permeability computation------" << endl << endl;
- cout << "-----Computed_Permeability-----" << endl;
+ cout << "-----Computed_Permeability-----" << endl;}
}
double FlowBoundingSphere::Compute_EffectiveRadius(Cell_handle cell, int j)
@@ -786,10 +814,12 @@
double FlowBoundingSphere::Compute_HydraulicRadius(Cell_handle cell, int j)
{
+// cerr << "test11" << endl;
RTriangulation& Tri = T[currentTes].Triangulation();
if (Tri.is_infinite(cell->neighbor(j))) return 0;
-
+// cerr << "test12" << endl;
/*double */Vpore = Volume_Pore_VoronoiFraction(cell, j);
+// cerr << "test13" << endl;
/*double */Ssolid = Surface_Solid_Pore(cell, j, SLIP_ON_LATERALS);
#ifdef FACET_BASED_FORCES
@@ -797,7 +827,7 @@
// if (FACET_BASED_FORCES) cell->info().facetSurfaces[j]=cell->info().facetSurfaces[j]*((facetSurfaces[j].length()-facetSphereCrossSections[j][0]-facetSphereCrossSections[j][1]-facetSphereCrossSections[j][2])/facetSurfaces[j].length());
#endif
// if (DEBUG_OUT) std::cerr << "total facet surface "<< cell->info().facetSurfaces[j] << " with solid sectors : " << cell->info().facetSphereCrossSections[j][0] << " " << cell->info().facetSphereCrossSections[j][1] << " " << cell->info().facetSphereCrossSections[j][2] << " difference "<<sqrt(cell->info().facetSurfaces[j].squared_length())-cell->info().facetSphereCrossSections[j][0]-cell->info().facetSphereCrossSections[j][2]-cell->info().facetSphereCrossSections[j][1]<<endl;
-
+// cerr << "test14" << endl;
//handle symmetry (tested ok)
if (/*SLIP_ON_LATERALS &&*/ fictious_vertex>0)
{
@@ -805,6 +835,7 @@
Real mult= fictious_vertex==1 ? multSym1 : multSym2;
return Vpore/Ssolid*mult;
}
+// cerr << "test15" << endl;
return Vpore/Ssolid;
}
@@ -832,7 +863,7 @@
void FlowBoundingSphere::GaussSeidel ()
{
- std::ofstream iter ("Gauss_Iterations", std::ios::app);
+// std::ofstream iter ("Gauss_Iterations", std::ios::app);
std::ofstream p_av ("P_moyenne", std::ios::app);
RTriangulation& Tri = T[currentTes].Triangulation();
int j = 0;
@@ -840,8 +871,8 @@
double tolerance = TOLERANCE;
double relax = RELAX;
- cout << "tolerance = " << tolerance << endl;
- cout << "relax = " << relax << endl;
+ if(DEBUG_OUT){ cout << "tolerance = " << tolerance << endl;
+ cout << "relax = " << relax << endl;}
Finite_cells_iterator cell_end = Tri.finite_cells_end();
@@ -881,7 +912,7 @@
} 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;
- iter << j << " " << dp_max/p_max << endl;
+// iter << j << " " << dp_max/p_max << endl;
int cel=0;
double Pav=0;
@@ -939,12 +970,13 @@
p_in_moy += cell->neighbor(j2)->info().p();}
}}}
+ if (DEBUG_OUT){
cout << "the maximum superior pressure is = " << p_out_max << " the min is = " << p_in_min << endl;
cout << "the maximum inferior pressure is = " << p_in_max << " the min is = " << p_out_min << endl;
cout << "superior average pressure is " << p_out_moy/cellQ1 << endl;
cout << "inferior average pressure is " << p_in_moy/cellQ2 << endl;
cout << "celle comunicanti in basso = " << cellQ2 << endl;
- cout << "celle comunicanti in alto = " << cellQ1 << endl;
+ cout << "celle comunicanti in alto = " << cellQ1 << endl;}
double density = 1;
double viscosity = 1;
@@ -955,6 +987,7 @@
double Ks= (Vdarcy) /GradH;
double k= Ks*viscosity/ (density*gravity);
+ if (DEBUG_OUT){
cout << "The incoming average flow rate is = " << Q2 << " m^3/s " << endl;
cout << "The outgoing average flow rate is = " << Q1 << " m^3/s " << endl;
cout << "The gradient of charge is = " << GradH << " [-]" << endl;
@@ -973,7 +1006,7 @@
kFile << "The gradient of charge is = " << GradH << " [-]" << endl;
kFile << "Darcy's velocity is = " << Vdarcy << " m/s" <<endl;
kFile << "The hydraulic conductivity of the sample is = " << Ks << " m/s" <<endl;
- kFile << "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;
@@ -1164,7 +1197,7 @@
double FlowBoundingSphere::Sample_Permeability(double& x_Min,double& x_Max ,double& y_Min,double& y_Max,double& z_Min,double& z_Max, string key)
{
- RTriangulation& Tri = T[currentTes].Triangulation();
+// RTriangulation& Tri = T[currentTes].Triangulation();
double Section = (x_Max-x_Min) * (z_Max-z_Min);
double DeltaY = y_Max-y_Min;
boundary(y_min_id).flowCondition=0;
@@ -1190,32 +1223,33 @@
return false;
}
-void FlowBoundingSphere::SliceField()
+void FlowBoundingSphere::SliceField(char *filename)
{
/** Pressure field along one cutting plane **/
RTriangulation& Tri = T[currentTes].Triangulation();
Cell_handle permeameter;
- std::ofstream consFile("slices",std::ios::out);
+ std::ofstream consFile(filename,std::ios::out);
int intervals = 400;
double Rx = 20* (x_max-x_min) /intervals;
double Ry = (y_max-y_min) /intervals;
double Rz = (z_max-z_min) /intervals;
- cout<<Rx<<" "<<Ry<<" "<<Rz<<" "<<z_max<<" "<<z_min<<" "<<y_max<<" "<<y_min<<" "<<x_max<<" "<<x_min<<endl;
+// cout<<Rx<<" "<<Ry<<" "<<Rz<<" "<<z_max<<" "<<z_min<<" "<<y_max<<" "<<y_min<<" "<<x_max<<" "<<x_min<<endl;
// for (double X=min(x_min,x_max); X<=max(x_min,x_max); X=X+abs(Rx)) {
double X=0.5;
for (double Y=min(y_max,y_min); Y<=max(y_max,y_min); Y=Y+abs(Ry)) {
for (double Z=min(z_min,z_max); Z<=max(z_min,z_max); Z=Z+abs(Rz)) {
- if (!isInsideSphere(X,Y,Z)) {
- permeameter = Tri.locate(Point(X, Y, Z));
- consFile << permeameter->info().p() <<" ";
- //cout <<"valeur trouvée";
- } else consFile << "Nan ";
+ permeameter = Tri.locate(Point(X, Y, Z));
+ consFile << permeameter->info().p() <<" ";
+// if (!isInsideSphere(X,Y,Z)) {
+// permeameter = Tri.locate(Point(X, Y, Z));
+// consFile << permeameter->info().p() <<" ";
+// //cout <<"valeur trouvée";
+// } else consFile << "Nan ";
}
- consFile << endl;
- }
+ consFile << endl;}
consFile << endl;
// }
consFile.close();
=== modified file 'lib/triangulation/FlowBoundingSphere.h'
--- lib/triangulation/FlowBoundingSphere.h 2010-10-20 11:28:25 +0000
+++ lib/triangulation/FlowBoundingSphere.h 2010-11-16 13:10:35 +0000
@@ -42,7 +42,7 @@
double RELAX;
double ks; //Hydraulic Conductivity
bool meanK_LIMIT, meanK_STAT, distance_correction;
- bool DEBUG_OUT;
+// bool DEBUG_OUT;
bool noCache;//flag for checking if cached values cell->unitForceVectors have been defined
bool computeAllCells;//exececute computeHydraulicRadius for all facets and all spheres (double cpu time but needed for now in order to define crossSections correctly)
double K_opt_factor;
@@ -55,7 +55,7 @@
// Boundary& boundary (int b) {return boundaries[b-id_offset];}
void mplot ( char *filename);
- void Localize ();
+ void Localize();
void Compute_Permeability();
@@ -82,7 +82,7 @@
// void AddBoundingPlanes();
// void AddBoundingPlanes(Real center[3], Real Extents[3], int id);
- Tesselation& Compute_Action ( );
+ void Compute_Action ( );
Tesselation& Compute_Action ( int argc, char *argv[ ], char *envp[ ] );
// Vecteur external_force_single_fictious ( Cell_handle cell );
void SpheresFileCreator ();
@@ -139,7 +139,7 @@
bool isInsideSphere ( double& x, double& y, double& z );
- void SliceField ();
+ void SliceField (char *filename);
void ComsolField();
void Interpolate ( Tesselation& Tes, Tesselation& NewTes );
@@ -157,7 +157,7 @@
void Average_Grain_Velocity();
void vtk_average_cell_velocity(RTriangulation &T, int id_sphere, int num_cells);
void ApplySinusoidalPressure(RTriangulation& Tri, double Pressure, double load_intervals);
-
+ void ApplySinusoidalPressure_Space_Time(RTriangulation& Tri, double Pressure, double load_intervals, double time, double dt);
// double surface_external_triple_fictious(Real position[3], Cell_handle cell, Boundary b);
// double surface_external_double_fictious(Cell_handle cell, Boundary b);
// double surface_external_single_fictious(Cell_handle cell, Boundary b);
=== modified file 'lib/triangulation/Network.cpp'
--- lib/triangulation/Network.cpp 2010-10-20 12:10:06 +0000
+++ lib/triangulation/Network.cpp 2010-11-16 13:10:35 +0000
@@ -1,3 +1,4 @@
+
#ifdef FLOW_ENGINE
#include "def_types.h"
@@ -15,7 +16,6 @@
#include "Network.h"
-
#define FAST
const double ONE_THIRD = 1.0/3.0;
@@ -82,13 +82,10 @@
{
Point& p1 = cell->info();
Point& p2 = cell->neighbor(j)->info();
-
fictious_vertex = Detect_facet_fictious_vertices (cell,j);
-
Sphere v [3];
Vertex_handle W [3];
for (int kk=0; kk<3; kk++) {W[kk] = cell->vertex(facetVertices[j][kk]);v[kk] = cell->vertex(facetVertices[j][kk])->point();}
-
switch (fictious_vertex) {
case (0) : {
Vertex_handle& SV1 = W[0];
@@ -111,8 +108,7 @@
/**Vpore**/ return Vtot - (Vsolid1 + Vsolid2);
}; break;
case (1) : {return volume_single_fictious_pore(cell->vertex(facetVertices[j][facetF1]), cell->vertex(facetVertices[j][facetRe1]), cell->vertex(facetVertices[j][facetRe2]), p1,p2, cell->info().facetSurfaces[j]);}; break;
- case (2) : {return volume_double_fictious_pore(cell->vertex(facetVertices[j][facetF1]), cell->vertex(facetVertices[j][facetF2]), cell->vertex(facetVertices[j][facetRe1]), p1,p2, cell->info().facetSurfaces[j]);}; break;
- }
+ case (2) : {return volume_double_fictious_pore(cell->vertex(facetVertices[j][facetF1]), cell->vertex(facetVertices[j][facetF2]), cell->vertex(facetVertices[j][facetRe1]), p1,p2, cell->info().facetSurfaces[j]);}; break;}
}
double Network::volume_single_fictious_pore(const Vertex_handle& SV1, const Vertex_handle& SV2, const Vertex_handle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface)
@@ -265,7 +261,7 @@
{
if (!facet_detected) fictious_vertex=Detect_facet_fictious_vertices(cell,j);
- RTriangulation& Tri = T[currentTes].Triangulation();
+// RTriangulation& Tri = T[currentTes].Triangulation();
Point& p1 = cell->info();
Point& p2 = cell->neighbor(j)->info();
@@ -469,41 +465,41 @@
boundsIds[4]= &z_min_id;
boundsIds[5]= &z_max_id;
- Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), Corner_min.y()-FAR*(Corner_max.x()-Corner_min.x()), 0.5*(Corner_max.z()-Corner_min.z()), FAR*(Corner_max.x()-Corner_min.x()), y_min_id, true);
+ Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), Corner_min.y()-FAR*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.z()+Corner_min.z()), FAR*(Corner_max.y()-Corner_min.y()), y_min_id, true);
boundaries[y_min_id-id_offset].p = Corner_min;
boundaries[y_min_id-id_offset].normal = Vecteur(0,1,0);
boundaries[y_min_id-id_offset].coordinate = 1;
- cout << "Bottom boundary has been created. ID = " << y_min_id << endl;
+ if(DEBUG_OUT) cout << "Bottom boundary has been created. ID = " << y_min_id << endl;
- Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), Corner_max.y() +FAR*(Corner_max.x()-Corner_min.x()), 0.5*(Corner_max.z()-Corner_min.z()), FAR*(Corner_max.x()-Corner_min.x()), y_max_id, true);
+ Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), Corner_max.y() +FAR*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.z()+Corner_min.z()), FAR*(Corner_max.y()-Corner_min.y()), y_max_id, true);
boundaries[y_max_id-id_offset].p = Corner_max;
boundaries[y_max_id-id_offset].normal = Vecteur(0,-1,0);
boundaries[y_max_id-id_offset].coordinate = 1;
- cout << "Top boundary has been created. ID = " << y_max_id << endl;
+ if(DEBUG_OUT) cout << "Top boundary has been created. ID = " << y_max_id << endl;
- Tes.insert(Corner_min.x()-FAR*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.z()-Corner_min.z()), FAR*(Corner_max.y()-Corner_min.y()), x_min_id, true);
+ Tes.insert(Corner_min.x()-FAR*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.y()+Corner_min.y()), 0.5*(Corner_max.z()+Corner_min.z()), FAR*(Corner_max.y()-Corner_min.y()), x_min_id, true);
boundaries[x_min_id-id_offset].p = Corner_min;
boundaries[x_min_id-id_offset].normal = Vecteur(1,0,0);
boundaries[x_min_id-id_offset].coordinate = 0;
- cout << "Left boundary has been created. ID = " << x_min_id << endl;
+ if(DEBUG_OUT) cout << "Left boundary has been created. ID = " << x_min_id << endl;
- Tes.insert(Corner_max.x() +FAR*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.z()-Corner_min.z()), FAR*(Corner_max.y()-Corner_min.y()), x_max_id, true);
+ Tes.insert(Corner_max.x() +FAR*(Corner_max.y()-Corner_min.y()), 0.5*(Corner_max.y()+Corner_min.y()), 0.5*(Corner_max.z()+Corner_min.z()), FAR*(Corner_max.y()-Corner_min.y()), x_max_id, true);
boundaries[x_max_id-id_offset].p = Corner_max;
boundaries[x_max_id-id_offset].normal = Vecteur(-1,0,0);
boundaries[x_max_id-id_offset].coordinate = 0;
- cout << "Right boundary has been created. ID = " << x_max_id << endl;
+ if(DEBUG_OUT) cout << "Right boundary has been created. ID = " << x_max_id << endl;
- Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), 0.5*(Corner_max.y()-Corner_min.y()), Corner_min.z()-FAR*(Corner_max.y()-Corner_min.y()), FAR*(Corner_max.y()-Corner_min.y()), z_min_id, true);
+ Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), 0.5*(Corner_max.y()+Corner_min.y()), Corner_min.z()-FAR*(Corner_max.y()-Corner_min.y()), FAR*(Corner_max.y()-Corner_min.y()), z_min_id, true);
boundaries[z_min_id-id_offset].p = Corner_min;
boundaries[z_min_id-id_offset].normal = Vecteur(0,0,1);
boundaries[z_min_id-id_offset].coordinate = 2;
- cout << "Front boundary has been created. ID = " << z_min_id << endl;
+ if(DEBUG_OUT) cout << "Front boundary has been created. ID = " << z_min_id << endl;
Tes.insert(0.5*(Corner_min.x() +Corner_max.x()), 0.5*(Corner_max.y()-Corner_min.y()), Corner_max.z() +FAR*(Corner_max.y()-Corner_min.y()), FAR*(Corner_max.y()-Corner_min.y()), z_max_id, true);
boundaries[z_max_id-id_offset].p = Corner_max;
boundaries[z_max_id-id_offset].normal = Vecteur(0,0,-1);
boundaries[z_max_id-id_offset].coordinate = 2;
- cout << "Back boundary has been created. ID = " << z_max_id << endl;
+ if(DEBUG_OUT) cout << "Back boundary has been created. ID = " << z_max_id << endl;
for (int k=0;k<6;k++) {
boundaries[k].flowCondition = 1;
@@ -571,6 +567,7 @@
vtk_infinite_vertices = fict;
vtk_infinite_cells = Fictious;
+ num_particles = real;
}
// double Network::spherical_triangle_area ( Sphere STA1, Sphere STA2, Sphere STA3, Point PTA1 )
=== modified file 'lib/triangulation/Network.h'
--- lib/triangulation/Network.h 2010-10-20 11:28:25 +0000
+++ lib/triangulation/Network.h 2010-11-16 13:10:35 +0000
@@ -39,6 +39,7 @@
Tesselation T [2];
bool currentTes;
double x_min, x_max, y_min, y_max, z_min, z_max, Rmoy, SectionArea, Height, Vtotale;
+ bool DEBUG_OUT;
int nOfSpheres;
int x_min_id, x_max_id, y_min_id, y_max_id, z_min_id, z_max_id;
int* boundsIds [6];
@@ -48,7 +49,7 @@
Boundary boundaries [6];
Boundary& boundary (int b) {return boundaries[b-id_offset];}
short id_offset;
- int vtk_infinite_vertices, vtk_infinite_cells;
+ int vtk_infinite_vertices, vtk_infinite_cells, num_particles;
// int F1, F2, Re1, Re2; //values between 0..3, refers to one cell's fictious(F)/real(Re) vertices
// int facetRe1, facetRe2, facetRe3, facetF1, facetF2; //indices relative to the facet
=== modified file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h 2010-11-12 19:22:10 +0000
+++ lib/triangulation/def_types.h 2010-11-16 13:10:35 +0000
@@ -16,8 +16,6 @@
#include<yade/lib/base/Math.hpp> // typedef for Real
-//#define FLOW_ENGINE
-
namespace CGT{
//Robust kernel
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2010-11-07 11:46:20 +0000
+++ pkg/dem/FlowEngine.cpp 2010-11-16 13:10:35 +0000
@@ -18,7 +18,6 @@
#include <sys/stat.h>
#include <sys/types.h>
-
YADE_REQUIRE_FEATURE (CGAL);
CREATE_LOGGER (FlowEngine);
@@ -29,8 +28,8 @@
void FlowEngine::action ( )
{
if (!flow) {
- flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);
- first=true;Update_Triangulation=false;}
+ flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);
+ 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
{
@@ -52,79 +51,100 @@
if ( first ) Build_Triangulation( P_zero );
timingDeltas->checkpoint("Triangulating");
-
+
+ eps_vol_max=0.f;
+ std::ofstream eps_vol ("Epsilon_volume.txt", std::ios::app);
UpdateVolumes ( );
+ eps_vol << eps_vol_max << endl;
+ Eps_Vol_Cumulative += eps_vol_max;
+ 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);}
+ timingDeltas->checkpoint("Gauss-Seidel");
- flow->MGPost();
+ 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");
+ if (!CachedForces) flow->ComputeFacetForces();
+ else flow->ComputeFacetForcesWithCache();
+
+ timingDeltas->checkpoint("Compute_Forces");
///End Compute flow and forces
-
- 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[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.addTorque(id,t);
- }
-
- timingDeltas->checkpoint("Applying Forces");
-
- Real time = scene->time;
-
- int j = scene->iter;
- char file [50];
- mkdir("./Consol", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
- string consol = "./Consol/"+flow->key+"%d_Consol";
- const char* keyconsol = consol.c_str();
- sprintf (file, keyconsol, j);
- char *g = file;
- timingDeltas->checkpoint("Writing cons_files");
-
- MaxPressure = flow->PressureProfile( g, time, intervals);
-
- std::ofstream max_p ("pressures.txt", std::ios::app);
- max_p << j << " " << time << " " << MaxPressure << endl;
-
- std::ofstream settle ("settle.txt", std::ios::app);
- settle << j << " " << time << " " << currentStrain << endl;
-
- if ( scene->iter % PermuteInterval == 0 )
- { Update_Triangulation = true; }
-
- if ( Update_Triangulation && !first) { Build_Triangulation( );}
-
- first=false;
- Update_Triangulation = false;
-
- timingDeltas->checkpoint("Storing Max Pressure");
-
- flow->Average_Grain_Velocity();
- if (save_vtk) flow->save_vtk_file();
-
+
+ //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;
+ }
+ //cout << "STABILITY INDEX - IS = " << IS << endl;
+ //Istab << scene->time << " " << IS << endl;
+
+ timingDeltas->checkpoint("Applying Forces");
+
+ Real time = scene->time;
+ int j = scene->iter;
+
+ if (consolidation) {
+ char file [50];
+ mkdir("./Consol", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ string consol = "./Consol/"+flow->key+"%d_Consol";
+ const char* keyconsol = consol.c_str();
+ sprintf (file, keyconsol, j);
+ char *g = file;
+ timingDeltas->checkpoint("Writing cons_files");
+ MaxPressure = flow->PressureProfile( g, time, intervals);
+
+ std::ofstream max_p ("pressures.txt", std::ios::app);
+ max_p << j << " " << time << " " << MaxPressure << endl;
+
+ std::ofstream settle ("settle.txt", std::ios::app);
+ settle << j << " " << time << " " << currentStrain << endl;}
+
+// if ( scene->iter % PermuteInterval == 0 )
+// { Update_Triangulation = true; }
+
+ if ( Update_Triangulation && !first) { Build_Triangulation( );}
+
+ first=false;
+ Update_Triangulation = false;
+
+ timingDeltas->checkpoint("Storing Max Pressure");
+
+// 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);
+ }
// 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);
}
}
@@ -162,11 +182,12 @@
flow->SLIP_ON_LATERALS=slip_boundary;
flow->key = triaxialCompressionEngine->Key;
flow->k_factor = permeability_factor;
+ flow->DEBUG_OUT = Debug;
}
else
{
flow->currentTes=!flow->currentTes;
- if (flow->DEBUG_OUT) cout << "--------RETRIANGULATION-----------" << endl;
+ if (Debug) cout << "--------RETRIANGULATION-----------" << endl;
}
flow->T[flow->currentTes].Clear();
flow->T[flow->currentTes].max_id=-1;
@@ -175,7 +196,7 @@
AddBoundary ( );
Triangulate ( );
- if (flow->DEBUG_OUT) cout << endl << "Tesselating------" << endl << endl;
+ if (Debug) cout << endl << "Tesselating------" << endl << endl;
flow->T[flow->currentTes].Compute();
flow->Define_fictious_cells();
@@ -189,16 +210,17 @@
{
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);
+ 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;
}
else
{
- if (flow->DEBUG_OUT) cout << "---------UPDATE PERMEABILITY VALUE--------------" << endl;
+ if (Debug && compute_K) cout << "---------UPDATE PERMEABILITY VALUE--------------" << endl;
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 );}
@@ -214,7 +236,6 @@
void FlowEngine::AddBoundary ()
{
-
shared_ptr<Sphere> sph ( new Sphere );
int Sph_Index = sph->getClassIndexStatic();
@@ -242,7 +263,7 @@
contator+=1;
}
}
- flow->SectionArea = ( flow->x_max - flow->x_min ) * ( flow->z_max-flow->z_min );
+ flow->SectionArea = ( flow->x_max - flow->x_min ) * ( flow->z_max-flow->z_min );
flow->Vtotale = (flow->x_max-flow->x_min) * (flow->y_max-flow->y_min) * (flow->z_max-flow->z_min);
if (flow->DEBUG_OUT) {cout << "Section area = " << flow->SectionArea << endl;
@@ -255,23 +276,23 @@
cout << "z_min = " << flow->z_min << endl;
cout << "z_max = " << flow->z_max << endl;}
- if (flow->DEBUG_OUT) cout << "Adding Boundary------" << endl;
-
- shared_ptr<Box> bx ( new Box );
- int Bx_Index = bx->getClassIndexStatic();
-
- FOREACH ( const shared_ptr<Body>& b, *scene->bodies )
- {
- if ( !b ) continue;
- if ( b->shape->getClassIndex() == Bx_Index )
- {
- Box* w = YADE_CAST<Box*> ( b->shape.get() );
-// const Body::id_t& id = b->getId();
- Real center [3], Extent[3];
- for ( int h=0;h<3;h++ ) {center[h] = b->state->pos[h]; Extent[h] = w->extents[h];}
- wall_thickness = min(min(Extent[0],Extent[1]),Extent[2]);
- }
- }
+ if (Debug) cout << "Adding Boundary------" << endl;
+
+// shared_ptr<Box> bx ( new Box );
+// int Bx_Index = bx->getClassIndexStatic();
+//
+// FOREACH ( const shared_ptr<Body>& b, *scene->bodies )
+// {
+// if ( !b ) continue;
+// if ( b->shape->getClassIndex() == Bx_Index )
+// {
+// Box* w = YADE_CAST<Box*> ( b->shape.get() );
+// // const Body::id_t& id = b->getId();
+// Real center [3], Extent[3];
+// for ( int h=0;h<3;h++ ) {center[h] = b->state->pos[h]; Extent[h] = w->extents[h];}
+// wall_thickness = min(min(Extent[0],Extent[1]),Extent[2]);
+// }
+// }
id_offset = flow->T[flow->currentTes].max_id+1;
@@ -290,31 +311,30 @@
void FlowEngine::Triangulate ()
{
///Using Tesselation wrapper (faster)
- TesselationWrapper TW;
- if (TW.Tes) delete TW.Tes;
- TW.Tes = &(flow->T[flow->currentTes]);//point to the current Tes we have in Flowengine
- TW.insertSceneSpheres();//TW is now really inserting in FlowEngine, using the faster insert(begin,end)
- TW.Tes = NULL;//otherwise, Tes would be deleted by ~TesselationWrapper() at the end of the function.
+// TesselationWrapper TW;
+// if (TW.Tes) delete TW.Tes;
+// TW.Tes = &(flow->T[flow->currentTes]);//point to the current Tes we have in Flowengine
+// TW.insertSceneSpheres();//TW is now really inserting in FlowEngine, using the faster insert(begin,end)
+// TW.Tes = NULL;//otherwise, Tes would be deleted by ~TesselationWrapper() at the end of the function.
///Using one-by-one insertion
-// shared_ptr<Sphere> sph ( new Sphere );
-// int Sph_Index = sph->getClassIndexStatic();
-// FOREACH ( const shared_ptr<Body>& b, *scene->bodies )
-// {
-// if ( !b ) continue;
-// if ( b->shape->getClassIndex() == Sph_Index )
-// {
-// Sphere* s=YADE_CAST<Sphere*> ( b->shape.get() );
-// const Body::id_t& id = b->getId();
-// Real rad = s->radius;
-// Real x = b->state->pos[0];
-// Real y = b->state->pos[1];
-// Real z = b->state->pos[2];
-//
-// flow->T[flow->currentTes].insert(x, y, z, rad, id);
-//
-// contator+=1;
-// }
-// }
+ shared_ptr<Sphere> sph ( new Sphere );
+ int Sph_Index = sph->getClassIndexStatic();
+ FOREACH ( const shared_ptr<Body>& b, *scene->bodies )
+ {
+ if ( !b ) continue;
+ if ( b->shape->getClassIndex() == Sph_Index )
+ {
+ Sphere* s=YADE_CAST<Sphere*> ( b->shape.get() );
+ const Body::id_t& id = b->getId();
+ Real rad = s->radius;
+ Real x = b->state->pos[0];
+ Real y = b->state->pos[1];
+ Real z = b->state->pos[2];
+
+ flow->T[flow->currentTes].insert(x, y, z, rad, id);
+
+ }
+ }
}
void FlowEngine::Initialize_volumes ()
@@ -330,13 +350,12 @@
case ( 3 ) : cell->info().volume() = Volume_cell_triple_fictious ( cell ); break;
}
}
- if (flow->DEBUG_OUT) cout << "Volumes initialised." << endl;
+ if (Debug) cout << "Volumes initialised." << endl;
}
void FlowEngine::UpdateVolumes ()
{
- if (flow->DEBUG_OUT) cout << "Updating volumes.............." << endl;
-
+ if (Debug) cout << "Updating volumes.............." << endl;
Real deltaT = scene->dt;
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++ )
@@ -346,21 +365,25 @@
case ( 3 ):
{
cell->info().dv() = ( Volume_cell_triple_fictious ( cell ) - cell->info().volume() ) /deltaT;
+ eps_vol_max = max(eps_vol_max, (abs(cell->info().dv()*deltaT))/cell->info().volume());
cell->info().volume() = Volume_cell_triple_fictious ( cell );
}break;
case ( 2 ) :
{
cell->info().dv() = ( Volume_cell_double_fictious ( cell )-cell->info().volume() ) /deltaT;
+ eps_vol_max = max(eps_vol_max, (abs(cell->info().dv()*deltaT))/cell->info().volume());
cell->info().volume() = Volume_cell_double_fictious ( cell );
}break;
case ( 1 ) :
{
cell->info().dv() = ( Volume_cell_single_fictious ( cell )-cell->info().volume() ) /deltaT;
+ eps_vol_max = max(eps_vol_max, (abs(cell->info().dv()*deltaT))/cell->info().volume());
cell->info().volume() = Volume_cell_single_fictious ( cell );
}break;
case ( 0 ) :
{
cell->info().dv() = ( Volume_cell ( cell )-cell->info().volume() ) /deltaT;
+ eps_vol_max = max(eps_vol_max, (abs(cell->info().dv()*deltaT))/cell->info().volume());
cell->info().volume() = Volume_cell ( cell );
}break;
}
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2010-11-07 11:46:20 +0000
+++ pkg/dem/FlowEngine.hpp 2010-11-16 13:10:35 +0000
@@ -11,8 +11,9 @@
#include<yade/core/PartialEngine.hpp>
#include<yade/pkg/dem/TriaxialCompressionEngine.hpp>
#include<yade/lib/triangulation/FlowBoundingSphere.h>
-
-
+#include<yade/pkg/dem/TesselationWrapper.hpp>
+
+class TesselationWrapper;
class FlowEngine : public PartialEngine
{
@@ -26,7 +27,11 @@
bool Update_Triangulation;
bool currentTes;
int id_offset;
-
+ // double IS;
+ double eps_vol_max;
+ double Eps_Vol_Cumulative;
+ int ReTrg;
+
void Triangulate ();
void AddBoundary ();
void Build_Triangulation (double P_zero );
@@ -49,15 +54,21 @@
((bool,first,true,,"Controls the initialization/update phases"))
((bool,save_vtk,false,,"Enable/disable vtk files creation for visualization"))
((bool,save_mplot,false,,"Enable/disable mplot files creation"))
+ ((bool, save_mgpost, false,,"Enable/disable mgpost file creation"))
+ ((bool, slice_pressures, false, ,"Enable/Disable slice pressure measurement"))
+ ((bool, consolidation,false,,"Enable/Disable storing consolidation files"))
((bool, slip_boundary, true,, "Controls friction condition on lateral walls"))
((bool, blocked_grains, false,, "Grains will/won't be moved by forces"))
((bool,WaveAction, false,, "Allow sinusoidal pressure condition to simulate ocean waves"))
+ ((bool, TimeBC, false,,"Activate evolution in time of pressure B.C."))
((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"))
((int,PermuteInterval,100000,,"Pore space re-triangulation period"))
+ ((double,EpsVolPercent_RTRG,0.01,,"Percentuage of cumulate eps_vol at which retriangulation of pore space is performed"))
((bool,compute_K,true,,"Activates permeability measure within a granular sample"))
((bool,meanK_correction,true,,"Local permeabilities' correction through meanK threshold"))
((bool,meanK_opt,false,,"Local permeabilities' correction through an optimized threshold"))
Follow ups
