yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #03453
[Branch ~yade-dev/yade/trunk] Rev 2040: Few corrections in flow files.
------------------------------------------------------------
revno: 2040
committer: Emanuele Catalano <ecatalano@r2balme>
branch nick: trunk
timestamp: Fri 2010-02-19 12:12:48 +0100
message:
Few corrections in flow files.
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.h
pkg/dem/Engine/PartialEngine/FlowEngine.cpp
pkg/dem/Engine/PartialEngine/FlowEngine.hpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.
=== modified file 'lib/triangulation/FlowBoundingSphere.cpp'
--- lib/triangulation/FlowBoundingSphere.cpp 2010-02-18 22:40:58 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-02-19 11:12:48 +0000
@@ -27,6 +27,7 @@
const bool DEBUG_OUT = false;
const double RELAX = 1.9;
+
const double TOLERANCE = 1e-06;
const double ONE_THIRD = 1.0/3.0;
//! Use this factor, or minLength, to reduce max permeability values (see usage below))
@@ -36,7 +37,7 @@
//! Factors including the effect of 1/2 symmetry in hydraulic radii
const Real multSym1 = 1/pow(2,0.25);
const Real multSym2 = 1/pow(4,0.25);
-const bool SLIP_ON_LATERALS = true;//no-slip/symmetry conditions on lateral boundaries
+
const double FAR = 500;
const int facetVertices [4][3] = {{1,2,3},{0,2,3},{0,1,3},{0,1,2}};
const int permut3 [3][3] = {{0,1,2},{1,2,0},{2,0,1}};
@@ -68,6 +69,7 @@
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
}
void FlowBoundingSphere::Compute_Action()
@@ -290,7 +292,7 @@
int pass = 0;
RTriangulation& Tri = T[currentTes].Triangulation();
Finite_cells_iterator cell_end = Tri.finite_cells_end();
- cout << "Localizing..." << endl;
+ cout << "Localizing..." << endl;
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
pass = 0;
for (int i=0;i<4;i++) {
@@ -301,7 +303,8 @@
Boundary& bi = boundary(V->info().id());
// Boundary& bi = boundaries [V->info().id()];
- if (bi.flowCondition) {
+// if (bi.flowCondition) {
+ if ( V->info().id() != /*(unsigned int)*/ y_min_id && V->info().id() != y_max_id) {
// Inf/Sup has priority
if (!cell->info().isSuperior && !cell->info().isInferior) {
cell->info().isLateral = true;
@@ -352,9 +355,8 @@
Vecteur n;
-
- std::ofstream oFile("Hydraulic_Radius" ,std::ios::out);
- // std::ofstream kFile ( "Permeability" ,std::ios::out );
+ std::ofstream oFile( "Hydraulic_Radius",std::ios::out);
+ std::ofstream kFile ( "Permeability" ,std::ios::out );
Real meanK=0;
Real infiniteK=1e10;
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
@@ -987,20 +989,6 @@
return sk;
}
-// double FlowBoundingSphere::epsilon_V_triple_fictious ( Cell_handle cell )
-// {
-// Boundary b;
-//
-// for (int g=0;g<4;g++)
-// {
-// if (cell->vertex(g)->info().isFictious)
-// {
-// b = boundary ( cell->vertex ( g )->info().id() );
-// double sk = surface_external_triple_fictious (cell,b);
-// }
-// }
-// }
-
void FlowBoundingSphere::Fictious_cells()
{
Finite_cells_iterator cell_end = T[currentTes].Triangulation().finite_cells_end();
@@ -1047,7 +1035,6 @@
C[bi1.coordinate]=bi1.p[bi1.coordinate];
C[bi2.coordinate]=bi2.p[bi2.coordinate];
-
// cout << "A = " << A[0] << " " << A[1] << " " << A[2] << endl;
// cout << "B = " << B[0] << " " << B[1] << " " << B[2] << endl;
// cout << "C = " << C[0] << " " << C[1] << " " << C[2] << endl;
@@ -1197,22 +1184,20 @@
double FlowBoundingSphere::Compute_HydraulicRadius(RTriangulation& Tri, Cell_handle cell, int j)
{
-
- Cell_handle neighbour_cell = cell->neighbor(j);
- if (Tri.is_infinite(neighbour_cell)) return 0;
+ if (Tri.is_infinite(cell->neighbor(j))) return 0;
Point p1 = cell->info();
- Point p2 = neighbour_cell->info();
+ Point p2 = cell->neighbor(j)->info();
Vertex_handle SV1, SV2, SV3;
double Vtot=0, Vsolid=0, Vpore=0, Vsolid1=0, Vsolid2=0;
- int F1=0, F2=0, Re1=0, Re2=0;
- //indices relative to the facet
- int facetRe1=0, facetRe2=0, facetRe3=0, facetF1=0, facetF2=0;
- int fictious_vertex=0, real_vertex=0;
+ int F1=0, F2=0, Re1=0, Re2=0; //values between 0..3, refers to one cell's fictious(F)/real(Re) vertices
+ int facetRe1=0, facetRe2=0, facetRe3=0, facetF1=0, facetF2=0; //indices relative to the facet
+ int fictious_vertex=0, real_vertex=0; //counts total fictious/real vertices
double Ssolid1=0, Ssolid1n=0, Ssolid2=0, Ssolid2n=0, Ssolid3=0, Ssolid3n=0, Ssolid=0;
+ //solid portions within a pore, for cell () and neighbour_cell (+n). Ssolid total solid surface
//bool fictious_solid = false;
@@ -1250,7 +1235,7 @@
double A [3], B[3], C[3];
/**PORE VOLUME**/
- Vpore = volume_double_fictious_pore(cell->vertex(F1), cell->vertex(F2), cell->vertex(Re1), p1,p2, cell->info().facetSurfaces[j]);
+ Vpore = volume_double_fictious_pore(cell->vertex(F1), cell->vertex(F2), cell->vertex(Re1), p1, p2, cell->info().facetSurfaces[j]);
/** **/
/**PORE SOLID SURFACE**/
@@ -1606,7 +1591,7 @@
cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;
// save_vtk_file ( Tri );
}
- } while ((dp_max/p_max) > tolerance /*( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
+ } 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;
@@ -1855,7 +1840,10 @@
P_ave[k]/= (m+n+o);
consFile<<k<<" "<<time<<" "<<P_ave[k]<<endl;
if (k==15) Pressures.push_back(P_ave[k]);
- n=0; m=0; o=0; k++;
+ n=0;
+ m=0;
+ o=0;
+ k++;
}
}
@@ -1869,7 +1857,6 @@
boundary(y_max_id).flowCondition=0;
boundary(y_min_id).value=0;
boundary(y_max_id).value=1;
-// P_SUP = 1.0; P_INF = 0.0; P_INS = 0.0;
Initialize_pressures();
@@ -1878,7 +1865,7 @@
char *kk;
kk = (char*) key.c_str();
- Permeameter(Tri, P_INF, P_SUP, Section, DeltaY, kk);
+ Permeameter(Tri, boundary(y_min_id).value, boundary(y_max_id).value, Section, DeltaY, kk);
}
void FlowBoundingSphere::AddBoundingPlanes(Real center[3], Real Extents[3], int id)
=== modified file 'lib/triangulation/FlowBoundingSphere.h'
--- lib/triangulation/FlowBoundingSphere.h 2010-02-18 22:40:58 +0000
+++ lib/triangulation/FlowBoundingSphere.h 2010-02-19 11:12:48 +0000
@@ -36,6 +36,7 @@
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;
Boundary boundaries [6];
short id_offset;
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-02-17 10:25:54 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-02-19 11:12:48 +0000
@@ -5,7 +5,6 @@
* This program is free software; it is licensed under the terms of the *
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
-#ifdef FLOW_ENGINE
#include "FlowEngine.hpp"
#include<yade/core/Scene.hpp>
@@ -15,10 +14,15 @@
#include<yade/pkg-common/Wall.hpp>
#include<yade/pkg-common/Box.hpp>
+#ifdef FLOW_ENGINE
+
YADE_REQUIRE_FEATURE (CGAL);
CREATE_LOGGER (FlowEngine);
-std::ofstream plotFile ( "plot2",std::ios::out );
+std::ofstream cons_DAMP ("cons_DAMP", std::ios::out);
+std::ofstream cons_NONDAMP ("cons_NONDAMP", std::ios::out);
+std::ofstream settle_DAMP ("settle_DAMP", std::ios::out);
+std::ofstream settle_NONDAMP ("settle_NONDAMP", std::ios::out);
FlowEngine::~FlowEngine()
{
@@ -47,7 +51,6 @@
}
current_state = triaxialCompressionEngine->currentState;
- flow->key = triaxialCompressionEngine->Key;
if ( !first && current_state==3 )
{
@@ -92,11 +95,12 @@
sprintf (file, keyconsol, j);
char *g = file;
-// string pressures = +"%d_Consol";
flow->PermeameterCurve(flow->T[currentTes].Triangulation(), g, time);
- plotFile << j << " " << flow->Pressures[cons] << endl; cons++;
-// plotFile << "replot '" << j << "_Consol' using 2:0" << endl;
+ if (damped) {cons_DAMP << j << " " << time << " " << flow->Pressures[cons] << endl; cons++;}
+ if (!damped){cons_NONDAMP << j << " " << time << " " << flow->Pressures[cons] << endl; cons++;}
+ if (damped) {settle_DAMP << j << " " << time << " " << triaxialCompressionEngine->uniaxialEpsilonCurr << endl;}
+ if (!damped) {settle_NONDAMP << j << " " << time << " " << triaxialCompressionEngine->uniaxialEpsilonCurr << endl;}
if ( Omega::instance().getCurrentIteration() % PermuteInterval == 0 )
{
@@ -128,6 +132,8 @@
Initialize ( ncb, P_zero );
flow->Vtotalissimo=0; flow->Vsolid_tot=0; flow->Vporale=0; flow->Ssolid_tot=0;
+
+ flow->SLIP_ON_LATERALS=slip_boundary;
flow->k_factor = permeability_factor;
flow->Compute_Permeability ();
@@ -143,6 +149,8 @@
}
cout << y << " deltaV initialised -----------------" << endl;
+ flow->key = triaxialCompressionEngine->Key;
+
if (compute_K) {flow->Sample_Permeability ( flow->T[currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
Oedometer_Boundary_Conditions();
@@ -150,10 +158,10 @@
flow->GaussSeidel ( );
- plotFile << "unset key" << endl;
+// plotFile << "unset key" << endl;
+
// flow->Analytical_Consolidation();
-
first = false; cons=0;
}
}
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.hpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-02-17 10:25:54 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-02-19 11:12:48 +0000
@@ -12,6 +12,8 @@
#include<yade/pkg-dem/TriaxialCompressionEngine.hpp>
#include<yade/lib-triangulation/FlowBoundingSphere.h>
+#ifdef FLOW_ENGINE
+
class FlowEngine : public PartialEngine
{
private:
@@ -20,7 +22,7 @@
public :
Vector3r gravity;
- bool first;
+// bool first;
int current_state
,previous_state
@@ -46,6 +48,9 @@
YADE_CLASS_BASE_DOC_ATTRS(FlowEngine,PartialEngine,"An engine to solve the flow problem in saturated granular media",
((bool,isActivated,true,"Activates Flow Engine "))
+ ((bool,first,true,"Controls the initialization/update phases"))
+ ((bool, damped, true, ""))
+ ((bool, slip_boundary, false, "Controls friction condition on lateral walls"))
((bool,currentTes,false,"Identifies the current triangulation/tesselation of pore space"))
((double,P_zero,0,"Initial internal pressure for oedometer test"))
((int,PermuteInterval,100000,"Pore space re-triangulation period"))
