yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06847
[Branch ~yade-dev/yade/trunk] Rev 2685: - Update flow code.
------------------------------------------------------------
revno: 2685
committer: Emanuele Catalano <catalano@xxxxxxxxxxx
branch nick: yade
timestamp: Thu 2011-01-27 16:25:09 +0100
message:
- Update flow code.
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 2010-11-25 14:22:43 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2011-01-27 15:25:09 +0000
@@ -196,7 +196,8 @@
/** END PERMEABILITY CALCULATION**/
if(DEBUG_OUT) cerr << "TOTAL VOID VOLUME: " << Vporale <<endl;
-
+ if(DEBUG_OUT) cerr << "Porosity = " << V_porale_porosity / V_totale_porosity << endl;
+
/** STATISTICS **/
DisplayStatistics();
clock.top("DisplayStatistics");
@@ -299,10 +300,16 @@
void FlowBoundingSphere::Average_Cell_Velocity()
{
+ double ranchx = (x_max - x_min)/3;
+ double ranchy = (y_max - y_min)/3;
+ double ranchz = (z_max - z_min)/3;
+
+ bool yes = false;
+
RTriangulation& Tri = T[currentTes].Triangulation();
Point pos_av_facet;
int num_cells = 0;
- double facet_flow_rate;
+ double facet_flow_rate = 0;
std::ofstream oFile ( "Average_Cells_Velocity",std::ios::app );
Real tVel=0; Real tVol=0;
Finite_cells_iterator cell_end = Tri.finite_cells_end();
@@ -310,116 +317,101 @@
cell->info().av_vel() =CGAL::NULL_VECTOR;
num_cells++;
for ( int i=0; i<4; i++ ) {
+ if (!Tri.is_infinite(cell->neighbor(i))){
Vecteur Surfk = cell->info()-cell->neighbor(i)->info();
Real area = sqrt ( Surfk.squared_length() );
Surfk = Surfk/area;
// Vecteur facetNormal = Surfk/area;
Vecteur branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
-// pos_av_facet=CGAL::ORIGIN + ((cell->vertex(facetVertices[i][0])->point() - CGAL::ORIGIN) +
-// (cell->vertex(facetVertices[i][1])->point() - CGAL::ORIGIN) +
-// (cell->vertex(facetVertices[i][2])->point() - CGAL::ORIGIN))*0.3333333333;
pos_av_facet = (Point) cell->info() + ( branch*Surfk ) *Surfk;
- facet_flow_rate = (cell->info().k_norm())[i] * (cell->neighbor (i)->info().p() - cell->info().p());
+// pos_av_facet=CGAL::ORIGIN + ((cell->vertex(facetVertices[i][0])->point() - CGAL::ORIGIN) + (cell->vertex(facetVertices[i][1])->point() - CGAL::ORIGIN) + (cell->vertex(facetVertices[i][2])->point() - CGAL::ORIGIN))*0.3333333333;
+ facet_flow_rate = (cell->info().k_norm())[i] * (cell->info().p() - cell->neighbor (i)->info().p());
cell->info().av_vel() = cell->info().av_vel() + facet_flow_rate* ( pos_av_facet-CGAL::ORIGIN );
- }
+ }}
if (cell->info().volume()){ tVel+=cell->info().av_vel()[1]; tVol+=cell->info().volume();}
cell->info().av_vel() = cell->info().av_vel() /cell->info().volume();
-// cerr << cell->info().av_vel()<<" "<<facet_flow_rate<<" "<<cell->info().volume()<<endl;
-// oFile << cell->info().av_vel() << "Cell Pressure = " << cell->info().p() << endl;
}
-// cerr <<"TOT Vol/Vel: "<<tVol<<" "<<tVel<<endl;
-}
-
-void FlowBoundingSphere::Average_Grain_Velocity()
-{
+}
+
+bool FlowBoundingSphere::isOnSolid (double X, double Y, double Z)
+{
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ Finite_cells_iterator cell_end = Tri.finite_cells_end();
+ for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+ for (int i=0; i<4; i++){
+ double radius = sqrt(cell->vertex(i)->point().weight());
+ if (X < (cell->vertex(i)->point().x()+radius) && X > (cell->vertex(i)->point().x()-radius)){
+ if (Y < (cell->vertex(i)->point().y()+radius) && Y > (cell->vertex(i)->point().y()-radius)){
+ if (Z < (cell->vertex(i)->point().z()+radius) && Z > (cell->vertex(i)->point().z()-radius)){
+ return true;}}}}}
+ return false;
+}
+
+void FlowBoundingSphere::Average_Fluid_Velocity()
+{
+ Average_Cell_Velocity();
RTriangulation& Tri = T[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++) {
- V_it->info().vel() = CGAL::NULL_VECTOR;
- V_it->info().vol_cells() = 0;}
-
- Point pos_av_facet;
- double facet_flow_rate;
+ for (Finite_vertices_iterator V_it = Tri.finite_vertices_begin(); V_it != vertices_end; V_it++) {
+ num_vertex++;}
+
+ vector<Real> Volumes;
+ vector<CGT::Vecteur> VelocityVolumes;
+ VelocityVolumes.resize(num_vertex);
+ Volumes.resize(num_vertex);
+
+ for (Finite_vertices_iterator V_it = Tri.finite_vertices_begin(); V_it != vertices_end; V_it++) {
+ VelocityVolumes[V_it->info().id()]=CGAL::NULL_VECTOR;
+ Volumes[V_it->info().id()]=0.f;
+ correction[V_it->info().id()]=false;}
+
Finite_cells_iterator cell_end = Tri.finite_cells_end();
- for ( Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++ )
+ for ( Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++ )
{
- int pass=0;
- for ( int i=0; i<4; i++ )
- {
- if(!Tri.is_infinite(cell->neighbor(i))){
- pass+=1;
- Vecteur Surfk = cell->info()-cell->neighbor(i)->info();
- Real area = sqrt ( Surfk.squared_length() );
- Surfk = Surfk/area;
- Vecteur branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
- pos_av_facet = (Point) cell->info() + ( branch*Surfk ) *Surfk;
- facet_flow_rate = (cell->info().k_norm())[i] * (cell->info().p() - cell->neighbor (i)->info().p() );
-
- for (int g=0;g<4;g++)
- {cell->vertex ( g )->info().vel() = cell->vertex ( g )->info().vel() + facet_flow_rate*( pos_av_facet-CGAL::ORIGIN );}
- }
-// else cout << "CI SONO INFINITE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
- }
- for (int g=0;g<pass;g++)
- {cell->vertex ( g )->info().vol_cells() = cell->vertex ( g )->info().vol_cells() + cell->info().volume();}
- }
-
-// 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++) {
- V_it->info().vel() = V_it->info().vel() / V_it->info().vol_cells();}
-
-// for (Finite_vertices_iterator V_it = Tri.finite_vertices_begin(); V_it != vertices_end; V_it++) {
-// V_it->info().vel() = V_it->info().vel() / sqrt ( V_it->info().vel().squared_length() );}
+ if (cell->info().fictious()==0){
+ for (int i=0;i<4;i++){
+ VelocityVolumes[cell->vertex(i)->info().id()] = VelocityVolumes[cell->vertex(i)->info().id()] + cell->info().av_vel()*cell->info().volume();
+ Volumes[cell->vertex(i)->info().id()] = Volumes[cell->vertex(i)->info().id()] + cell->info().volume();}
+ }}
+
+ std::ofstream fluid_vel ("Velocity", std::ios::out);
+ double Rx = (x_max-x_min) /10;
+ double Ry = (y_max-y_min) /12;
+ double Rz = (z_max-z_min) /20;
+ Cell_handle cellula;
+
+ Vecteur Velocity = CGAL::NULL_VECTOR;
+ int i=0;
+ for(double X=x_min+Rx;X<x_max;X+=Rx){
+ for (double Y=y_min+Ry;Y<y_max;Y+=Ry){
+ Velocity = CGAL::NULL_VECTOR; i=0;
+ for (double Z=z_min+Rz;Z<z_max;Z+=Rz){
+ cellula = Tri.locate(Point(X,Y,Z));
+ for (int y=0;y<4;y++) {if (!cellula->vertex(y)->info().isFictious) {Velocity = Velocity + (VelocityVolumes[cellula->vertex(y)->info().id()]/Volumes[cellula->vertex(y)->info().id()]);i++;}}
+ }Velocity = Velocity/i;
+ fluid_vel << X << " " << Y << " " << Velocity << endl;
+ }}
}
-void FlowBoundingSphere::vtk_average_cell_velocity(RTriangulation &Tri, int id_sphere, int num_cells )
-{
- 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++);
-
- basicVTKwritter vtkfile((unsigned int) 4*num_cells, (unsigned int) num_cells);
-
- vtkfile.open(filename,"test");
-
- Tesselation::Vector_Cell tmp_cells;
- tmp_cells.resize(10000);
- Tesselation::VCell_iterator cells_it = tmp_cells.begin();
- Tesselation::VCell_iterator cells_end = Tri.incident_cells(T[currentTes].vertexHandles[id_sphere],cells_it);
-
- vtkfile.begin_vertices();
- for (Tesselation::VCell_iterator it = tmp_cells.begin(); it != cells_end; it++)
- {
- for (int y2=0; y2<4; y2++){
- double x,y,z;
- x = (double)((*it)->vertex(y2)->point().point()[0]);
- y = (double)((*it)->vertex(y2)->point().point()[1]);
- z = (double)((*it)->vertex(y2)->point().point()[2]);
- vtkfile.write_point(x,y,z);}
- }
-
- vtkfile.end_vertices();
-
- vtkfile.begin_cells();
- for (Tesselation::VCell_iterator it = tmp_cells.begin(); it != cells_end; it++)
- {
- vtkfile.write_cell((*it)->vertex(0)->info().id()-id_offset, (*it)->vertex(1)->info().id()-id_offset, (*it)->vertex(2)->info().id()-id_offset, (*it)->vertex(3)->info().id()-id_offset);
+double FlowBoundingSphere::Measure_bottom_Pore_Pressure()
+{
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ Cell_handle permeameter;
+
+ int intervals = 40;
+ double Rz = (z_max-z_min) /intervals;
+
+ double X=(x_max+x_min)/2;
+ double Y=y_min;
+ double pressure = 0.f;
+ int cell=0;
+ for (double Z=min(z_min,z_max); Z<=max(z_min,z_max); Z=Z+abs(Rz)) {
+ permeameter = Tri.locate(Point(X, Y, Z));
+ pressure+=permeameter->info().p();
+ cell++;
}
- vtkfile.end_cells();
-
-// vtkfile.begin_data("Force",POINT_DATA,SCALARS,FLOAT);
-// vtkfile.write_data((T[currentTes].vertexHandles[id_sphere]->info().forces)[1]);
-// vtkfile.end_data();
-
- vtkfile.begin_data("Velocity",CELL_DATA,SCALARS,FLOAT);
-
- for (Tesselation::VCell_iterator it = tmp_cells.begin(); it != cells_end; it++)
- {
- vtkfile.write_data((*it)->info().av_vel()[1]);
- }
- vtkfile.end_data();
+ return pressure/cell;
}
void FlowBoundingSphere::ComputeFacetForces()
@@ -612,9 +604,8 @@
cout << "TotalForce = "<< TotalForce << endl;
}
-void FlowBoundingSphere::ApplySinusoidalPressure(RTriangulation& Tri, double Pressure, double load_intervals)
+void FlowBoundingSphere::ApplySinusoidalPressure(RTriangulation& Tri, double Amplitude, double Average_Pressure, double load_intervals)
{
-
double step = 1/load_intervals;
Tesselation::Vector_Cell tmp_cells;
tmp_cells.resize(10000);
@@ -627,35 +618,9 @@
if(!Tri.is_infinite(*it)){
Point& p1 = (*it)->info();
Cell_handle& cell = *it;
- if (p1.x()<0) cell->info().p() = Pressure;
- else if (p1.x()>x_max) cell->info().p() = 0.f;
- else if (p1.x()>(alpha*(x_max-x_min)) && p1.x()<((alpha+step)*(x_max-x_min))) cell->info().p() = (Pressure/2)*(1+cos(alpha*M_PI));
- }
- }
- }
-}
-
-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);
- }
+ if (p1.x()<x_min) cell->info().p() = Average_Pressure+Amplitude;
+ else if (p1.x()>x_max) cell->info().p() = Average_Pressure-Amplitude;
+ else if (p1.x()>(x_min+alpha*(x_max-x_min)) && p1.x()<(x_min+(alpha+step)*(x_max-x_min))) cell->info().p() = Average_Pressure + (Amplitude)*(cos(alpha*M_PI));
}
}
}
@@ -681,7 +646,7 @@
{
if (DEBUG_OUT) cout << "----Computing_Permeability------" << endl;
RTriangulation& Tri = T[currentTes].Triangulation();
- Vsolid_tot = 0, Vtotalissimo = 0, Vporale = 0, Ssolid_tot = 0;
+ Vsolid_tot = 0, Vtotalissimo = 0, Vporale = 0, Ssolid_tot = 0, V_totale_porosity=0, V_porale_porosity=0;
Finite_cells_iterator cell_end = Tri.finite_cells_end();
Cell_handle neighbour_cell;
@@ -1125,21 +1090,25 @@
cout << "celle comunicanti in alto = " << cellQ1 << endl;}
double density = 1;
- double viscosity = 1;
- double gravity = 1;
+ double viscosity = 0.001;
+ double gravity = 9.80665;
double Vdarcy = Q1/Section;
- double GradP = abs(P_Inf-P_Sup) /DeltaY;
- double GradH = GradP/ (density*gravity);
- double Ks= (Vdarcy) /GradH;
- double k= Ks*viscosity/ (density*gravity);
+// double GradP = abs(P_Inf-P_Sup) /DeltaY;
+ double DeltaP = abs(P_Inf-P_Sup);
+// double GradH = GradP/ (density*gravity);
+ double DeltaH = DeltaP/ (density*gravity);
+// double Ks= (Vdarcy) /GradH;
+// double k = Ks*viscosity/ (density*gravity);
+ double k = viscosity*Vdarcy*DeltaY / DeltaP; /**m²**/
+ double Ks = k*(density*gravity)/viscosity; /**m/s**/
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;
+ cout << "The gradient of charge is = " << DeltaH/DeltaY << " [-]" << endl;
cout << "Darcy's velocity is = " << Vdarcy << " m/s" <<endl;
cout << "The permeability of the sample is = " << k << " m^2" <<endl;}
-
+ cout << "The Darcy permeability of the sample is = " << k/9.869233e-13 << " darcys" << endl;
kFile << "the maximum superior pressure is = " << p_in_max << " the min is = " << p_in_min << endl;
kFile << "the maximum inferior pressure is = " << p_out_max << " the min is = " << p_out_min << endl;
kFile << "superior average pressure is " << p_in_moy/cellQ2 << endl;
@@ -1148,11 +1117,11 @@
kFile << "celle comunicanti in basso = " << cellQ1 << endl;
kFile << "The incoming average flow rate is = " << Q2 << " m^3/s " << endl;
kFile << "The outgoing average flow rate is = " << Q1 << " m^3/s " << endl;
- kFile << "The gradient of charge is = " << GradH << " [-]" << endl;
+ kFile << "The gradient of charge is = " << DeltaH/DeltaY << " [-]" << 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;
- // cout << "The Darcy permeability of the sample is = " << k_darcy/0.987e-12 << " darcys" << endl;
+ kFile << "The Darcy permeability of the sample is = " << k/9.869233e-13 << " darcys" << endl;
cout << endl << "The hydraulic conductivity of the sample is = " << Ks << " m/s" << endl << endl;
return Ks;
@@ -1201,7 +1170,7 @@
void FlowBoundingSphere::save_vtk_file()
{
- RTriangulation& Tri = T[currentTes].Triangulation();
+ RTriangulation& Tri = T[currentTes].Triangulation();
static unsigned int number = 0;
char filename[80];
mkdir("./VTK", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
@@ -1239,14 +1208,21 @@
}
vtkfile.end_data();
- vtkfile.begin_data("Velocity",POINT_DATA,VECTORS,FLOAT);
- for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v)
- {if (!v->info().isFictious) vtkfile.write_data((v->info().vel())[0],(v->info().vel())[1],(v->info().vel())[2]);}
- vtkfile.end_data();
+// vtkfile.begin_data("Velocity",POINT_DATA,VECTORS,FLOAT);
+// for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v)
+// {if (!v->info().isFictious) vtkfile.write_data((v->info().vel())[0],(v->info().vel())[1],(v->info().vel())[2]);}
+// vtkfile.end_data();
// vtkfile.begin_data("Velocity",CELL_DATA,VECTORS,FLOAT);
-// for (Finite_cells_iterator cell = T.finite_cells_begin(); cell != T.finite_cells_end(); ++cell) {
-// if (!cell->info().fictious()){vtkfile.write_data((cell->info().av_vel())[0],(cell->info().av_vel())[1],(cell->info().av_vel())[2]);}
+// bool control = false;
+// double zero = 0.f;
+// for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
+// if (!cell->info().fictious()==0){
+// for(int j=0;j++;j<3){if (cell->info().av_vel()[j]>0) control = true;}
+// if (!control) vtkfile.write_data((cell->info().av_vel())[0],(cell->info().av_vel())[1],(cell->info().av_vel())[2]);
+// else vtkfile.write_data(zero,zero,zero);
+// control=false;
+// }
// }
// vtkfile.end_data();
}
@@ -1321,8 +1297,8 @@
solid=false;
for (Finite_vertices_iterator V_it = Tri.finite_vertices_begin(); V_it != Tri.finite_vertices_end(); V_it++) {
- double rayon = sqrt(V_it->point().weight());
- if ((sqrt(pow((x- (V_it->point()[0])),2) +pow((y- (V_it->point()[1])),2) +pow((z- (V_it->point()[2])),2))) <= rayon) solid=true;
+ double radius = sqrt(V_it->point().weight());
+ if ((sqrt(pow((x- (V_it->point()[0])),2) +pow((y- (V_it->point()[1])),2) +pow((z- (V_it->point()[2])),2))) <= radius) solid=true;
}
if (solid) voxelfile << 1;
else voxelfile << 0;
@@ -1403,8 +1379,8 @@
{
RTriangulation& Tri = T[currentTes].Triangulation();
for (Finite_vertices_iterator V_it = Tri.finite_vertices_begin(); V_it != Tri.finite_vertices_end(); V_it++) {
- double rayon = V_it->point().weight();
- if (pow((x- (V_it->point()[0])),2) +pow((y- (V_it->point()[1])),2) +pow((z- (V_it->point()[2])),2) <= rayon) return true;
+ double radius = V_it->point().weight();
+ if (pow((x- (V_it->point()[0])),2) +pow((y- (V_it->point()[1])),2) +pow((z- (V_it->point()[2])),2) <= radius) return true;
}
return false;
}
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2010-11-25 14:22:43 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2011-01-27 15:25:09 +0000
@@ -36,15 +36,11 @@
virtual ~FlowBoundingSphere();
FlowBoundingSphere();
-// int x_min_id, x_max_id, y_min_id, y_max_id, z_min_id, z_max_id;
-// int* boundsIds [6];
-// bool currentTes;
bool SLIP_ON_LATERALS;
double TOLERANCE;
double RELAX;
double ks; //Hydraulic Conductivity
bool meanK_LIMIT, meanK_STAT, distance_correction;
-// bool DEBUG_OUT;
bool OUTPUT_BOUDARIES_RADII;
bool noCache;//flag for checking if cached values cell->unitForceVectors have been defined
@@ -55,10 +51,7 @@
int Iterations;
bool RAVERAGE;
-// Boundary boundaries [6];
int walls_id[6];
-// short id_offset;
-// Boundary& boundary (int b) {return boundaries[b-id_offset];}
void mplot ( char *filename);
void Localize();
@@ -67,13 +60,8 @@
virtual void GaussSeidel ( );
virtual void ResetNetwork();
-// void Compute_Forces ();
void Fictious_cells ( );
-// Tesselation T [2];
-
-// double x_min, x_max, y_min, y_max, z_min, z_max, Rmoy;
-// Real Vsolid_tot, Vtotalissimo, Vporale, Ssolid_tot;
double k_factor; //permeability moltiplicator
std::string key; //to give to consolidation files a name with iteration number
std::vector<double> Pressures; //for automatic write maximum pressures during consolidation
@@ -82,19 +70,11 @@
double P_SUP, P_INF, P_INS;
-// void AddBoundingPlanes ( Tesselation& Tes, double x_Min,double x_Max ,double y_Min,double y_Max,double z_Min,double z_Max );
-// void AddBoundingPlanes(bool yade);
-// void AddBoundingPlanes();
-// void AddBoundingPlanes(Real center[3], Real Extents[3], int id);
-
Tesselation& Compute_Action ( );
Tesselation& Compute_Action ( int argc, char *argv[ ], char *envp[ ] );
Tesselation& LoadPositions(int argc, char *argv[ ], char *envp[ ]);
-// Vecteur external_force_single_fictious ( Cell_handle cell );
void SpheresFileCreator ();
-// void Analytical_Consolidation ( );
void DisplayStatistics();
-// void Boundary_Conditions ( RTriangulation& Tri );
void Initialize_pressures ( double P_zero );
/// Define forces using the same averaging volumes as for permeability
void ComputeTetrahedralForces();
@@ -116,59 +96,24 @@
double Compute_EffectiveRadius(Cell_handle cell, int j);
double Compute_EquivalentRadius(Cell_handle cell, int j);
-// double crossProduct ( double x[3], double y[3] );
-
-// double surface_solid_facet ( Sphere ST1, Sphere ST2, Sphere ST3 );
-// Vecteur surface_double_fictious_facet ( Vertex_handle fSV1, Vertex_handle fSV2, Vertex_handle SV3 );
-// Vecteur surface_single_fictious_facet ( Vertex_handle fSV1, Vertex_handle SV2, Vertex_handle SV3 );
-
-// double surface_solid_double_fictious_facet ( Vertex_handle ST1, Vertex_handle ST2, Vertex_handle ST3 );
-
-// double surface_external_triple_fictious (Cell_handle cell, Boundary b );
-// 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 );
-
void GenerateVoxelFile ( );
RTriangulation& Build_Triangulation ( Real x, Real y, Real z, Real radius, unsigned const id );
void Build_Tessalation ( RTriangulation& Tri );
-// double spherical_triangle_area ( Sphere STA1, Sphere STA2, Sphere STA3, Point PTA1 );
-
-// double fast_spherical_triangle_area ( const Sphere& STA1, const Point& STA2, const Point& STA3, const Point& PTA1 );
-// Real solid_angle ( const Point& STA1, const Point& STA2, const Point& STA3, const Point& PTA1 );
-// double spherical_triangle_volume ( const Sphere& ST1, const Point& PT1, const Point& PT2, const Point& PT3 );
-// Real fast_solid_angle ( const Point& STA1, const Point& PTA1, const Point& PTA2, const Point& PTA3 );
-
bool isInsideSphere ( double& x, double& y, double& z );
void SliceField (char *filename);
void ComsolField();
void Interpolate ( Tesselation& Tes, Tesselation& NewTes );
-
-// double volume_single_fictious_pore ( Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3, Point PV1 );
- //Fast version, assign surface of facet for future forces calculations (pointing from PV2 to PV1)
-// double volume_single_fictious_pore ( const Vertex_handle& SV1, const Vertex_handle& SV2, const Vertex_handle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface);
-// double volume_double_fictious_pore ( Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3, Point PV1 );
- //Fast version, assign surface of facet for future forces calculations (pointing from PV2 to PV1)
-// double volume_double_fictious_pore (Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3, Point& PV1, Point& PV2, Vecteur& facetSurface);
-
-// double PoreVolume (RTriangulation& Tri, Cell_handle cell);
-// void Slice_Average_Cell_Velocity();
- int Average_Cell_Velocity(int id_sphere, RTriangulation& Tri);
void Average_Cell_Velocity();
- 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);
-
+ void Average_Fluid_Velocity();
+ void ApplySinusoidalPressure(RTriangulation& Tri, double Amplitude, double Average_Pressure, double load_intervals);
+ bool isOnSolid (double X, double Y, double Z);
+ double Measure_bottom_Pore_Pressure();
+
//Solver?
int useSolver;//(0 : GaussSeidel, 1 : TAUCS, 2 : PARDISO)
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2010-11-25 14:22:43 +0000
+++ pkg/dem/FlowEngine.cpp 2011-01-27 15:25:09 +0000
@@ -51,9 +51,6 @@
currentStress = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];
currentStrain = triaxialCompressionEngine->strain[1];
-// current_state = triaxialCompressionEngine->currentState;
-// if (current_state == 3){
-
timingDeltas->start();
if (first) Build_Triangulation(P_zero);
@@ -62,6 +59,7 @@
if (!first) {
eps_vol_max=0.f;
UpdateVolumes ( );
+
Eps_Vol_Cumulative += eps_vol_max;
if ((Eps_Vol_Cumulative > EpsVolPercent_RTRG || retriangulationLastIter>1000) && retriangulationLastIter>10) {
Update_Triangulation = true;
@@ -121,7 +119,7 @@
flow->SliceField(f);
}
- if (save_vtk) flow->save_vtk_file();
+ if (save_vtk) {flow->save_vtk_file();}
}
// if ( scene->iter % PermuteInterval == 0 )
// { Update_Triangulation = true; }
@@ -131,10 +129,8 @@
Update_Triangulation = false;
}
- flow->Average_Cell_Velocity();
-// flow->Average_Grain_Velocity();
-// 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);
+ if (velocity_profile) /*flow->FluidVelocityProfile();*/flow->Average_Grain_Velocity();
+ if (liquefaction) bottom_seabed_pressure=flow->Measure_bottom_Pore_Pressure();
first=false;
// }
@@ -203,8 +199,7 @@
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);
+ if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Amplitude, Sinus_Average, 30);
flow->TOLERANCE=Tolerance;
flow->RELAX=Relax;
}
@@ -218,8 +213,7 @@
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]);
- 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);
+ if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Amplitude, Sinus_Average, 30);
flow->TOLERANCE=Tolerance;
flow->RELAX=Relax;
}
@@ -389,12 +383,46 @@
}
}
-Real FlowEngine::Volume_cell_single_fictious ( CGT::Cell_handle cell)
+Real FlowEngine::Volume_cell_single_fictious ( CGT::Cell_handle cell )
{
Real V[3][3];
int b=0;
int w=0;
+ Real Wall_coordinate=0;
+
+ for ( int y=0;y<4;y++ )
+ {
+ if ( ! ( cell->vertex ( y )->info().isFictious ) )
+ {
+ const shared_ptr<Body>& sph = Body::byId
+ ( cell->vertex ( y )->info().id(), scene );
+ for ( int g=0;g<3;g++ ) V[w][g]=sph->state->pos[g];
+ w++;
+ }
+ else
+ {
+ b = cell->vertex ( y )->info().id()-flow->id_offset;
+ const shared_ptr<Body>& wll = Body::byId ( b , scene );
+ if (!flow->boundaries[b].useMaxMin) Wall_coordinate = wll->state->pos[flow->boundaries[b].coordinate]+(flow->boundaries[b].normal[flow->boundaries[b].coordinate])*wall_thickness/2;
+ else Wall_coordinate = flow->boundaries[b].p[flow->boundaries[b].coordinate];
+ }
+ }
+
+ 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 = ( CGAL::cross_product ( CGT::Vecteur ( v1[0],v1[1],v1[2] ),
+ CGT::Vecteur ( v2[0],v2[1],v2[2] ) ) *
+ flow->boundaries[b].normal ) * ( 0.33333333333* ( V[0][flow->boundaries[b].coordinate]+ V[1][flow->boundaries[b].coordinate]+ V[2][flow->boundaries[b].coordinate] ) - Wall_coordinate );
+
+ return abs ( Volume );
+ /*
+ Real V[3][3];
+ int b=0;
+ int w=0;
+
Real Wall_point[3];
for ( int y=0;y<4;y++ )
@@ -411,7 +439,7 @@
b = cell->vertex ( y )->info().id()-flow->id_offset;
const shared_ptr<Body>& wll = Body::byId ( b , scene );
for ( int i=0;i<3;i++ ) Wall_point[i] = flow->boundaries[b].p[i];
- Wall_point[flow->boundaries[b].coordinate] = wll->state->pos[flow->boundaries[b].coordinate]+(flow->boundaries[b].normal[flow->boundaries[b].coordinate])*wall_thickness;
+ Wall_point[flow->boundaries[b].coordinate] = wll->state->pos[flow->boundaries[b].coordinate]+(flow->boundaries[b].normal[flow->boundaries[b].coordinate])*wall_thickness/2;
}
}
@@ -423,29 +451,28 @@
CGT::Vecteur ( v2[0],v2[1],v2[2] ) ) *
flow->boundaries[b].normal ) * ( 0.33333333333* ( V[0][flow->boundaries[b].coordinate]+ V[1][flow->boundaries[b].coordinate]+ V[2][flow->boundaries[b].coordinate] ) - Wall_point[flow->boundaries[b].coordinate] );
- return abs ( Volume );
+ return abs ( Volume );*/
}
Real FlowEngine::Volume_cell_double_fictious ( CGT::Cell_handle cell)
{
- Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0};
+ Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0};
Real B[3]={0, 0, 0}, BS[3]={0, 0, 0}, BT[3]={0, 0, 0};
Real C[3]={0, 0, 0}, CS[3]={0, 0, 0}, CT[3]={0, 0, 0};
int b[2];
-
- Real Wall_point[2][3];
-
+ Real Wall_coordinate[2];
int j=0;
bool first_sph=true;
+
for ( int g=0;g<4;g++ )
{
if ( cell->vertex ( g )->info().isFictious )
{
b[j] = cell->vertex ( g )->info().id()-flow->id_offset;
const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
- for ( int i=0;i<3;i++ ) Wall_point[j][i] = flow->boundaries[b[j]].p[i];
- Wall_point[j][flow->boundaries[b[j]].coordinate] = wll->state->pos[flow->boundaries[b[j]].coordinate] +(flow->boundaries[b[j]].normal[flow->boundaries[b[j]].coordinate])*wall_thickness;
+ if (!flow->boundaries[b[j]].useMaxMin) Wall_coordinate[j] = wll->state->pos[flow->boundaries[b[j]].coordinate] +(flow->boundaries[b[j]].normal[flow->boundaries[b[j]].coordinate])*wall_thickness/2;
+ else Wall_coordinate[j] = flow->boundaries[b[j]].p[flow->boundaries[b[j]].coordinate];
j++;
}
else if ( first_sph )
@@ -461,9 +488,9 @@
for ( int k=0;k<3;k++ ) { B[k]=BS[k]=BT[k]=sph2->state->pos[k]; }
}
}
-
- AS[flow->boundaries[b[0]].coordinate]=BS[flow->boundaries[b[0]].coordinate] = Wall_point[0][flow->boundaries[b[0]].coordinate];
- AT[flow->boundaries[b[1]].coordinate]=BT[flow->boundaries[b[1]].coordinate] = Wall_point[1][flow->boundaries[b[1]].coordinate];
+
+ AS[flow->boundaries[b[0]].coordinate]=BS[flow->boundaries[b[0]].coordinate] = Wall_coordinate[0];
+ AT[flow->boundaries[b[1]].coordinate]=BT[flow->boundaries[b[1]].coordinate] = Wall_coordinate[1];
for ( int h=0;h<3;h++ ) {C[h]= ( A[h]+B[h] ) /2; CS[h]= ( AS[h]+BS[h] ) /2; CT[h]= ( AT[h]+BT[h] ) /2;}
@@ -475,6 +502,54 @@
Real Volume = ( CGAL::cross_product ( v1,v2 ) *flow->boundaries[b[0]].normal ) *h;
return abs ( Volume );
+
+// Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0};
+// Real B[3]={0, 0, 0}, BS[3]={0, 0, 0}, BT[3]={0, 0, 0};
+// Real C[3]={0, 0, 0}, CS[3]={0, 0, 0}, CT[3]={0, 0, 0};
+//
+// int b[2];
+//
+// Real Wall_point[2][3];
+//
+// int j=0;
+// bool first_sph=true;
+// for ( int g=0;g<4;g++ )
+// {
+// if ( cell->vertex ( g )->info().isFictious )
+// {
+// b[j] = cell->vertex ( g )->info().id()-flow->id_offset;
+// const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
+// for ( int i=0;i<3;i++ ) Wall_point[j][i] = flow->boundaries[b[j]].p[i];
+// Wall_point[j][flow->boundaries[b[j]].coordinate] = wll->state->pos[flow->boundaries[b[j]].coordinate] +(flow->boundaries[b[j]].normal[flow->boundaries[b[j]].coordinate])*wall_thickness/2;
+// j++;
+// }
+// else if ( first_sph )
+// {
+// const shared_ptr<Body>& sph1 = Body::byId
+// ( cell->vertex ( g )->info().id(), scene );
+// for ( int k=0;k<3;k++ ) { A[k]=AS[k]=AT[k]=sph1->state->pos[k]; first_sph=false;}
+// }
+// else
+// {
+// const shared_ptr<Body>& sph2 = Body::byId
+// ( cell->vertex ( g )->info().id(), scene );
+// for ( int k=0;k<3;k++ ) { B[k]=BS[k]=BT[k]=sph2->state->pos[k]; }
+// }
+// }
+//
+// AS[flow->boundaries[b[0]].coordinate]=BS[flow->boundaries[b[0]].coordinate] = Wall_point[0][flow->boundaries[b[0]].coordinate];
+// AT[flow->boundaries[b[1]].coordinate]=BT[flow->boundaries[b[1]].coordinate] = Wall_point[1][flow->boundaries[b[1]].coordinate];
+//
+// for ( int h=0;h<3;h++ ) {C[h]= ( A[h]+B[h] ) /2; CS[h]= ( AS[h]+BS[h] ) /2; CT[h]= ( AT[h]+BT[h] ) /2;}
+//
+// CGT::Vecteur v1 ( AT[0]-BT[0],AT[1]-BT[1],AT[2]-BT[2] );
+// CGT::Vecteur v2 ( C[0]-CT[0],C[1]-CT[1],C[2]-CT[2] );
+//
+// Real h = C[flow->boundaries[b[0]].coordinate]- CS[flow->boundaries[b[0]].coordinate];
+//
+// Real Volume = ( CGAL::cross_product ( v1,v2 ) *flow->boundaries[b[0]].normal ) *h;
+//
+// return abs ( Volume );
}
Real FlowEngine::Volume_cell_triple_fictious ( CGT::Cell_handle cell)
@@ -482,41 +557,80 @@
Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0}, AW[3]={0, 0, 0};
int b[3];
- Real Wall_point[3][3];
+ Real Wall_coordinate[3];
int j=0;
-
+
for ( int g=0;g<4;g++ )
{
if ( cell->vertex ( g )->info().isFictious )
{
- b[j] = cell->vertex ( g )->info().id()-flow->id_offset;
- const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
- for ( int i=0;i<3;i++ ) Wall_point[j][i] = flow->boundaries[b[j]].p[i];
- Wall_point[j][flow->boundaries[b[j]].coordinate] = wll->state->pos[flow->boundaries[b[j]].coordinate]+(flow->boundaries[b[j]].normal[flow->boundaries[b[j]].coordinate])*wall_thickness;
- j++;
+ b[j] = cell->vertex ( g )->info().id()-flow->id_offset;
+ const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
+ if (!flow->boundaries[b[j]].useMaxMin) Wall_coordinate[j] = wll->state->pos[flow->boundaries[b[j]].coordinate] + (flow->boundaries[b[j]].normal[flow->boundaries[b[j]].coordinate])*wall_thickness/2;
+ else Wall_coordinate[j] = flow->boundaries[b[j]].p[flow->boundaries[b[j]].coordinate];
+ j++;
}
else
{
- const shared_ptr<Body>& sph = Body::byId
- ( cell->vertex ( g )->info().id(), scene );
- for ( int k=0;k<3;k++ ) { A[k]=AS[k]=AT[k]=AW[k]=sph->state->pos[k];}
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ for ( int k=0;k<3;k++ ) { A[k]=AS[k]=AT[k]=AW[k]=sph->state->pos[k];}
}
}
-
- AS[flow->boundaries[b[0]].coordinate]= AT[flow->boundaries[b[0]].coordinate]= AW[flow->boundaries[b[0]].coordinate]= Wall_point[0][flow->boundaries[b[0]].coordinate];
- AT[flow->boundaries[b[1]].coordinate]= Wall_point[1][flow->boundaries[b[1]].coordinate];
- AW[flow->boundaries[b[2]].coordinate]= Wall_point[2][flow->boundaries[b[2]].coordinate];
-
+
+ AS[flow->boundaries[b[0]].coordinate]= AT[flow->boundaries[b[0]].coordinate]= AW[flow->boundaries[b[0]].coordinate]= Wall_coordinate[0];
+ AT[flow->boundaries[b[1]].coordinate]= Wall_coordinate[1];
+ AW[flow->boundaries[b[2]].coordinate]= Wall_coordinate[2];
+
CGT::Vecteur v1 ( AS[0]-AT[0],AS[1]-AT[1],AS[2]-AT[2] );
CGT::Vecteur v2 ( AS[0]-AW[0],AS[1]-AW[1],AS[2]-AW[2] );
CGT::Vecteur h ( AT[0] - A[0], AT[1] - A[1], AT[2] - A[2] );
Real Volume = ( CGAL::cross_product ( v1,v2 ) ) * h;
-
+
return abs ( Volume );
}
+// Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0}, AW[3]={0, 0, 0};
+//
+// int b[3];
+// Real Wall_point[3][3];
+// int j=0;
+
+// for ( int g=0;g<4;g++ )
+// {
+// if ( cell->vertex ( g )->info().isFictious )
+// {
+// b[j] = cell->vertex ( g )->info().id()-flow->id_offset;
+// const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
+// for ( int i=0;i<3;i++ ) Wall_point[j][i] = flow->boundaries[b[j]].p[i];
+// Wall_point[j][flow->boundaries[b[j]].coordinate] = wll->state->pos[flow->boundaries[b[j]].coordinate]+(flow->boundaries[b[j]].normal[flow->boundaries[b[j]].coordinate])*wall_thickness/2;
+// j++;
+// // cout << "id_wall = " << cell->vertex ( g )->info().id() << " Position = " << endl;
+// }
+// else
+// {
+// const shared_ptr<Body>& sph = Body::byId
+// ( cell->vertex ( g )->info().id(), scene );
+// for ( int k=0;k<3;k++ ) { A[k]=AS[k]=AT[k]=AW[k]=sph->state->pos[k];}
+//
+// }
+// }
+//
+// AS[flow->boundaries[b[0]].coordinate]= AT[flow->boundaries[b[0]].coordinate]= AW[flow->boundaries[b[0]].coordinate]= Wall_point[0][flow->boundaries[b[0]].coordinate];
+// AT[flow->boundaries[b[1]].coordinate]= Wall_point[1][flow->boundaries[b[1]].coordinate];
+// AW[flow->boundaries[b[2]].coordinate]= Wall_point[2][flow->boundaries[b[2]].coordinate];
+//
+// CGT::Vecteur v1 ( AS[0]-AT[0],AS[1]-AT[1],AS[2]-AT[2] );
+// CGT::Vecteur v2 ( AS[0]-AW[0],AS[1]-AW[1],AS[2]-AW[2] );
+//
+// CGT::Vecteur h ( AT[0] - A[0], AT[1] - A[1], AT[2] - A[2] );
+//
+// Real Volume = ( CGAL::cross_product ( v1,v2 ) ) * h;
+//
+// return abs ( Volume );
+// }
+
Real FlowEngine::Volume_cell ( CGT::Cell_handle cell)
{
Vector3r A[4];
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2010-11-25 14:22:43 +0000
+++ pkg/dem/FlowEngine.hpp 2011-01-27 15:25:09 +0000
@@ -29,7 +29,7 @@
Vector3r gravity;
int current_state;
Real wall_thickness;
- bool Update_Triangulation;
+// bool Update_Triangulation;
bool currentTes;
int id_offset;
// double IS;
@@ -59,22 +59,25 @@
((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, velocity_profile, false, ,"Enable/Disable slice velocity 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."))
+ ((double, Sinus_Amplitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied"))
+ ((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"))
+ ((bool, Update_Triangulation, 0,,"If true the medium is retriangulated"))
((int,PermuteInterval,100000,,"Pore space re-triangulation period"))
((double, eps_vol_max, 0,,"Maximal absolute volumetric strain computed at each iteration"))
((double, EpsVolPercent_RTRG,0.01,,"Percentuage of cumulate eps_vol at which retriangulation of pore space is performed"))
((double, porosity, 0,,"Porosity computed at each retriangulation"))
- ((bool,compute_K,true,,"Activates permeability measure within a granular sample"))
+ ((bool,compute_K,false,,"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"))
((double,permeability_factor,1.0,,"a permability multiplicator"))
@@ -85,6 +88,8 @@
((double, currentStrain, 0,, "Current value of axial strain"))
((int, intervals, 30,, "Number of layers for pressure measurements"))
((int, useSolver, 0,, "Solver to use"))
+ ((bool, liquefaction, false,,"Compute bottom_seabed_pressure if true, see below"))
+ ((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"))
@@ -97,7 +102,6 @@
((double, Pressure_BACK_Boundary, 0,,"Pressure imposed on back boundary"))
((double, Pressure_LEFT_Boundary, 0,, "Pressure imposed on left boundary"))
((double, Pressure_RIGHT_Boundary, 0,, "Pressure imposed on right boundary"))
- ((double, Sinus_Pressure, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied"))
((int, id_sphere, 0,, "Average velocity will be computed for all cells incident to that sphere"))
((bool, BOTTOM_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
((bool, TOP_Boundary_MaxMin, 1,,"If true bounding sphere is added as function fo max/min sphere coord, if false as function of yade wall position"))
