yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07105
[Branch ~yade-dev/yade/trunk] Rev 2745: - Fluid velocity calculation takes into account of facets velocities
------------------------------------------------------------
revno: 2745
committer: Emanuele Catalano <catalano@xxxxxxxxxxx
branch nick: yade
timestamp: Thu 2011-02-17 18:43:37 +0100
message:
- Fluid velocity calculation takes into account of facets velocities
- Update flow code
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/Network.hpp
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 2011-01-27 15:59:39 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2011-02-17 17:43:37 +0000
@@ -73,7 +73,7 @@
fictious_vertex = 0;
SectionArea = 0, Height=0, Vtotale=0;
vtk_infinite_vertices=0, vtk_infinite_cells=0;
-
+ VISCOSITY = 1;
tess_based_force = true;
for (int i=0;i<6;i++) boundsIds[i] = 0;
minPermLength=-1;
@@ -304,6 +304,7 @@
Point pos_av_facet;
int num_cells = 0;
double facet_flow_rate = 0;
+ double volume_facet_translation = 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();
@@ -311,6 +312,7 @@
cell->info().av_vel() =CGAL::NULL_VECTOR;
num_cells++;
for ( int i=0; i<4; i++ ) {
+ volume_facet_translation = 0;
if (!Tri.is_infinite(cell->neighbor(i))){
Vecteur Surfk = cell->info()-cell->neighbor(i)->info();
Real area = sqrt ( Surfk.squared_length() );
@@ -320,7 +322,8 @@
pos_av_facet = (Point) cell->info() + ( branch*Surfk ) *Surfk;
// 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 );
+ for (int y=0;y<3;y++) volume_facet_translation += (cell->info().facetVelocity())[i][y]*cell->info().facetSurfaces[i][y];
+ cell->info().av_vel() = cell->info().av_vel() + (facet_flow_rate - volume_facet_translation) * ( 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();
@@ -387,24 +390,31 @@
}}
}
-double FlowBoundingSphere::Measure_bottom_Pore_Pressure()
+void FlowBoundingSphere::Measure_Pore_Pressure(double Wall_up_y, double Wall_down_y)
{
RTriangulation& Tri = T[currentTes].Triangulation();
Cell_handle permeameter;
-
- int intervals = 40;
- double Rz = (z_max-z_min) /intervals;
+ std::ofstream capture ("Pressure_profile", std::ios::app);
+ int intervals = 5;
+ int captures = 6;
+ double Rz = (z_max-z_min)/intervals;
+ double Ry = (Wall_up_y-Wall_down_y)/captures;
double X=(x_max+x_min)/2;
- double Y=y_min;
+ double Y = 0;
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)) {
+// for (double Y=min(y_min,y_max); Y<=max(y_min,y_max); Y=Y+abs(Ry)) {cell=0; pressure=0.f;
+// for (double Y=Wall_down_y;Y<=Wall_up_y;Y+=Ry) {cell=0; pressure=0.f;
+ for (int i=0; i<captures; i++){
+ for (double Z=min(z_min,z_max); Z<=max(z_min,z_max); Z+=abs(Rz)) {
permeameter = Tri.locate(Point(X, Y, Z));
pressure+=permeameter->info().p();
cell++;
}
- return pressure/cell;
+ Y += Ry;
+ capture << pressure/cell << endl;}
+
}
void FlowBoundingSphere::ComputeFacetForces()
@@ -644,7 +654,7 @@
Cell_handle neighbour_cell;
- double k=0, distance = 0, radius = 0, viscosity = 1;
+ double k=0, distance = 0, radius = 0, viscosity = VISCOSITY;
int NEG=0, POS=0, pass=0;
bool ref = Tri.finite_cells_begin()->info().isvisited;
@@ -1083,7 +1093,7 @@
cout << "celle comunicanti in alto = " << cellQ1 << endl;}
double density = 1;
- double viscosity = 0.001;
+ double viscosity = VISCOSITY;
double gravity = 9.80665;
double Vdarcy = Q1/Section;
// double GradP = abs(P_Inf-P_Sup) /DeltaY;
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2011-01-27 15:25:09 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2011-02-17 17:43:37 +0000
@@ -68,7 +68,7 @@
bool tess_based_force; //allow the force computation method to be chosen from FlowEngine
Real minPermLength; //min branch length for Poiseuille
- double P_SUP, P_INF, P_INS;
+ double P_SUP, P_INF, P_INS, VISCOSITY;
Tesselation& Compute_Action ( );
Tesselation& Compute_Action ( int argc, char *argv[ ], char *envp[ ] );
@@ -112,7 +112,7 @@
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();
+ void Measure_Pore_Pressure(double Wall_up_y, double Wall_down_y);
//Solver?
int useSolver;//(0 : GaussSeidel, 1 : TAUCS, 2 : PARDISO)
=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp 2010-11-25 14:22:43 +0000
+++ lib/triangulation/Network.hpp 2011-02-17 17:43:37 +0000
@@ -52,20 +52,15 @@
short id_offset;
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
int fictious_vertex;
-// bool facet_detected;
-
-// void DisplayStatistics();
-// void AddBoundingPlanes(bool yade);
+
void AddBoundingPlanes();
void AddBoundingPlane (bool yade, Vecteur Normal, int id_wall);
void AddBoundingPlane (Real center[3], double thickness, Vecteur Normal, int id_wall );
-// void AddBoundingPlanes ( Real center[3], Real Extents[3], int id );
+
void Define_fictious_cells( );
int Detect_facet_fictious_vertices (Cell_handle& cell, int& j);
-// double Volume_Pore (Cell_handle cell);
+
double Volume_Pore_VoronoiFraction ( Cell_handle& cell, int& j);
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(const Vertex_handle& SV1, const Vertex_handle& SV2, const Vertex_handle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface);
=== modified file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h 2010-11-17 15:19:09 +0000
+++ lib/triangulation/def_types.h 2011-02-17 17:43:37 +0000
@@ -58,6 +58,7 @@
// Surface vectors of facets, pointing from outside toward inside the cell
std::vector<Vecteur> facetSurfaces;
+ std::vector<Vecteur> facetVelocities;
// Reflects the geometrical property of the cell, so that the force by cell fluid on grain "i" is pressure*unitForceVectors[i]
std::vector<Vecteur> unitForceVectors;
// Store the area of triangle-sphere intersections for each facet (used in forces definition)
@@ -75,6 +76,7 @@
cell_force.resize(4);
facetSurfaces.resize(4);
facetSphereCrossSections.resize(4);
+ facetVelocities.resize(4);
unitForceVectors.resize(4);
for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidSurfaces[k][l]=0;
RayHydr.resize(4, 0);
@@ -114,6 +116,7 @@
inline std::vector<double>& k_norm (void) {return module_permeability;}
inline std::vector< Vecteur >& facetSurf (void) {return facetSurfaces;}
+ inline std::vector< Vecteur >& facetVelocity (void) {return facetVelocities;}
inline std::vector<Vecteur>& force (void) {return cell_force;}
inline std::vector<double>& Rh (void) {return RayHydr;}
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2011-01-27 15:59:39 +0000
+++ pkg/dem/FlowEngine.cpp 2011-02-17 17:43:37 +0000
@@ -25,6 +25,8 @@
{
}
+const int facetVertices [4][3] = {{1,2,3},{0,2,3},{0,1,3},{0,1,2}};
+
void FlowEngine::action()
{
if (!isActivated) return;
@@ -130,7 +132,10 @@
}
if (velocity_profile) /*flow->FluidVelocityProfile();*/flow->Average_Fluid_Velocity();
- if (liquefaction) bottom_seabed_pressure=flow->Measure_bottom_Pore_Pressure();
+ if (first && liquefaction){
+ wall_up_y = flow->y_max;
+ wall_down_y = flow->y_min;}
+ if (liquefaction) flow->Measure_Pore_Pressure(wall_up_y, wall_down_y);
first=false;
// }
@@ -173,6 +178,7 @@
flow->k_factor = permeability_factor;
flow->DEBUG_OUT = Debug;
flow->useSolver = useSolver;
+ flow->VISCOSITY = viscosity;
flow->T[flow->currentTes].Clear();
flow->T[flow->currentTes].max_id=-1;
@@ -390,6 +396,10 @@
int w=0;
Real Wall_coordinate=0;
+
+ Vector3r Vel[3];
+ int id = 0;
+ double Vel_x=0, Vel_y=0, Vel_z=0;
for ( int y=0;y<4;y++ )
{
@@ -398,6 +408,7 @@
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];
+ Vel[w] = sph->state->pos;
w++;
}
else
@@ -406,8 +417,23 @@
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];
+ id = y;
}
}
+
+ for ( int y=0;y<4;y++ ){
+ for ( int j=0;j<3;j++ ){
+ if (cell->vertex(j)->info().isFictious){
+ Vel_x += 0.33333333333*(Vel[facetVertices[y][j]][0]);
+ Vel_y += 0.33333333333*(Vel[facetVertices[y][j]][1]);
+ Vel_z += 0.33333333333*(Vel[facetVertices[y][j]][2]);}
+ else for (int j2=0;j2<3;j2++){
+ if (!cell->vertex(j2)->info().isFictious){
+ Vel_x += 0.5*(Vel[facetVertices[y][j]][0]);
+ Vel_y += 0.5*(Vel[facetVertices[y][j]][1]);
+ Vel_z += 0.5*(Vel[facetVertices[y][j]][2]);}}}
+ CGT::Vecteur Vel_facet ( Vel_x, Vel_y, Vel_z );
+ (cell->info().facetVelocity())[y] = Vel_facet*1;}
double v1[3], v2[3];
@@ -418,40 +444,6 @@
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++ )
- {
- 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 );
- 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/2;
- }
- }
-
- 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_point[flow->boundaries[b].coordinate] );
-
- return abs ( Volume );*/
}
Real FlowEngine::Volume_cell_double_fictious ( CGT::Cell_handle cell)
@@ -465,6 +457,10 @@
int j=0;
bool first_sph=true;
+ Vector3r Vel[2];
+ vector<int> id;
+ id.resize(2);
+
for ( int g=0;g<4;g++ )
{
if ( cell->vertex ( g )->info().isFictious )
@@ -480,15 +476,24 @@
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;}
+ Vel[0] = sph1->state->vel; id[0]=g;
}
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]; }
+ Vel[1] = sph2->state->vel; id[1]=g;
}
}
+ for (int y=0;y<4;y++){
+ if ( y == id[0] ) {CGT::Vecteur Vel_facet (Vel[0][0], Vel[0][1], Vel[0][2]); (cell->info().facetVelocity())[y] = Vel_facet*1;}
+ else if ( y == id[1] ) {CGT::Vecteur Vel_facet (Vel[1][0], Vel[1][1], Vel[1][2]); (cell->info().facetVelocity())[y] = Vel_facet*1;}
+ else { CGT::Vecteur Vel_facet (0.5*(Vel[0][0]+Vel[0][1]+Vel[0][2]),0.5*(Vel[1][0]+Vel[1][1]+Vel[1][2]),0.5*(Vel[2][0]+Vel[2][1]+Vel[2][2]));
+ (cell->info().facetVelocity())[y] = Vel_facet*1;}}
+
+
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];
@@ -502,54 +507,6 @@
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)
@@ -560,6 +517,8 @@
Real Wall_coordinate[3];
int j=0;
+ Vector3r Vel;
+
for ( int g=0;g<4;g++ )
{
if ( cell->vertex ( g )->info().isFictious )
@@ -574,9 +533,13 @@
{
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];}
+ Vel = sph->state->vel;
}
}
+ CGT::Vecteur Vel_facet (Vel[0], Vel[1], Vel[2]);
+ for (int y=0;y<4;y++) (cell->info().facetVelocity())[y] = Vel_facet*1;
+
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];
@@ -591,58 +554,27 @@
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];
- int j=0;
+ Vector3r Vel[4];
+
+ double Vel_x=0, Vel_y=0, Vel_z=0;
for ( int y=0;y<4;y++ )
{
- const shared_ptr<Body>& sph = Body::byId
- ( cell->vertex ( y )->info().id(), scene );
- for ( int i=0;i<3;i++ ) A[j]=sph->state->pos;
- j++;
+ const shared_ptr<Body>& sph = Body::byId(cell->vertex ( y )->info().id(), scene);
+ A[y]=sph->state->pos;
+ Vel[y]=sph->state->vel;
}
+
+ for ( int y=0;y<4;y++ ){
+ for ( int j=0;j<3;j++ ){
+ Vel_x += 0.33333333333*(Vel[facetVertices[y][j]][0]);
+ Vel_y += 0.33333333333*(Vel[facetVertices[y][j]][1]);
+ Vel_z += 0.33333333333*(Vel[facetVertices[y][j]][2]);}
+ CGT::Vecteur Vel_facet ( Vel_x, Vel_y, Vel_z );
+ (cell->info().facetVelocity())[y] = Vel_facet*1; }
CGT::Point p1 ( ( A[0] ) [0], ( A[0] ) [1], ( A[0] ) [2] );
CGT::Point p2 ( ( A[1] ) [0], ( A[1] ) [1], ( A[1] ) [2] );
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2011-01-27 15:25:09 +0000
+++ pkg/dem/FlowEngine.hpp 2011-02-17 17:43:37 +0000
@@ -33,6 +33,7 @@
bool currentTes;
int id_offset;
// double IS;
+ double wall_up_y, wall_down_y;
double Eps_Vol_Cumulative;
int ReTrg;
void Triangulate ();
@@ -81,6 +82,7 @@
((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"))
+ ((double,viscosity,1.0,,"viscosity of fluid"))
((Real,loadFactor,1.1,,"Load multiplicator for oedometer test"))
((double, K, 0,, "Permeability of the sample"))
((double, MaxPressure, 0,, "Maximal value of water pressure within the sample"))
@@ -89,7 +91,7 @@
((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"))
+// ((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"))
