yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #03875
[Branch ~yade-dev/yade/trunk] Rev 2131: - Adjustments to the code
------------------------------------------------------------
revno: 2131
committer: Emanuele Catalano <ecatalano@r2balme>
branch nick: trunk
timestamp: Wed 2010-04-07 17:03:09 +0200
message:
- Adjustments to the code
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-03-23 11:30:28 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-04-07 15:03:09 +0000
@@ -363,8 +363,8 @@
Vecteur n;
// std::ofstream oFile( "Hydraulic_Radius",std::ios::out);
-// std::ofstream kFile ( "Permeability" ,std::ios::out );
- Real meanK=0;
+ std::ofstream kFile ( "LocalPermeabilities" ,std::ios::app );
+ Real meanK=0, k_moy;
Real infiniteK=1e10;
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
p1 = cell->info();
@@ -393,29 +393,36 @@
if (distance!=0) {
if (minPermLength>0) distance=max(minPermLength,distance);
k = (M_PI * pow(radius,4)) / (8*viscosity*distance);
- meanK += k;
+// meanK += k;
} else k = infiniteK;//Will be corrected in the next loop
// cout << "k="<<k<<endl;
(cell->info().k_norm())[j]= k*k_factor;
// (cell->info().facetSurf())[j]= k*n;
(neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]= k*k_factor;
+
+ meanK += 1/k;
+
// (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
}
cell->info().isvisited = !ref;
- // for (int y=0;y<4;y++) cout << "Permeability " << y << " = " << (cell->info().k_norm())[y] << endl;
- // for ( int y=0;y<4;y++ ) kFile << y << " = " << ( cell->info().k_norm() ) [y] << endl;
+// for (int y=0;y<4;y++) cout << "Permeability " << y << " = " << (cell->info().k_norm())[y] << endl;
+ for ( int y=0;y<4;y++ ) kFile << ( cell->info().k_norm() ) [y] << endl;
}
// A loop to reduce K maxima, needs a better equation : the mean value is influenced by the very big K
+ kFile << "----------reduction reduction reduction---------" << endl;
meanK /= pass;
+ k_moy = 1/meanK;
Real maxKdivKmean = MAXK_DIV_KMEAN;
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
for (int j=0; j<4; j++) {
// if (cell->info().k_norm()[j]>maxKdivKmean*meanK) cerr <<"Adjusting permeability : " <<cell->info().k_norm()[j]<<" > "<<maxKdivKmean<<" * "<<meanK<<endl;
- (cell->info().k_norm())[j] = min(cell->info().k_norm()[j], maxKdivKmean*meanK);
+// (cell->info().k_norm())[j] = min(cell->info().k_norm()[j], maxKdivKmean*meanK);
+ (cell->info().k_norm())[j] = min(cell->info().k_norm()[j], k_moy);
+ kFile << (cell->info().k_norm())[j] << endl;
}
}
if (DEBUG_OUT) cout << "POS = " << POS << " NEG = " << NEG << " pass = " << pass << endl;
@@ -436,6 +443,7 @@
cout<<grains<<"grains - " <<"Vtotale = " << Vtotale << " Vgrains = " << Vgrains << " Vporale1 = " << Vtotale-Vgrains << endl;
cout << "Vtotalissimo = " << Vtotalissimo << " Vsolid_tot = " << Vsolid_tot << " Vporale2 = " << Vporale << " Ssolid_tot = " << Ssolid_tot << endl<< endl;
}
+
cout << "-----Computed_Permeability-----" << endl;
}
@@ -1538,7 +1546,7 @@
}
}
-void FlowBoundingSphere::GaussSeidel()
+void FlowBoundingSphere::GaussSeidel ()
{
std::ofstream iter ("Gauss_Iterations", std::ios::app);
RTriangulation& Tri = T[currentTes].Triangulation();
@@ -1566,7 +1574,7 @@
m=0, n=0;
// if ((!cell->info().isSuperior && !cell->info().isInferior)) {
- if (!cell->info().fictious() || cell->info().isLateral) {
+ if ( !cell->info().fictious() || cell->info().isLateral ) {
cell2++;
for (int j2=0; j2<4; j2++) {
if (!Tri.is_infinite(cell->neighbor(j2))) {
@@ -1821,8 +1829,8 @@
}
}
-double FlowBoundingSphere::PermeameterCurve(RTriangulation& Tri, char *filename, Real time)
-{
+double FlowBoundingSphere::PermeameterCurve(RTriangulation& Tri, char *filename, Real time, int intervals)
+{
/** CONSOLIDATION CURVES **/
Cell_handle permeameter;
int n=0;
@@ -1831,16 +1839,16 @@
int o=0;
vector<double> P_ave;
std::ofstream consFile(filename, std::ios::out);
- consFile << "j " << " time " << "p_ave " << endl;
- int intervals = 30;
+// consFile << "j " << " time " << "p_ave " << endl;
+// int intervals = 30;
double Rx = (x_max-x_min) /intervals;
double Ry = (y_max-y_min) /intervals;
double Rz = (z_max-z_min) /intervals;
- for (double Y=y_min; Y<y_max; Y=Y+Ry) {
+ for (double Y=y_min; Y<=y_max+Ry/10; Y=Y+Ry) {
P_ave.push_back(0);
- for (double X=x_min; X<x_max; X=X+Rx) {
- for (double Z=z_min; Z<z_max; Z=Z+Rz) {
+ for (double X=x_min; X<=x_max+Ry/10; X=X+Rx) {
+ for (double Z=z_min; Z<=z_max+Ry/10; Z=Z+Rz) {
permeameter = Tri.locate(Point(X, Y, Z));
if (permeameter->info().isInside) {
n++;
@@ -1856,7 +1864,7 @@
}
}
}
- // cout << "P_ave[" << k << "] = " << P_ave[k]/ ( m+n+o ) << " n = " << n << " m = " << m << " o = " << o << endl;
+// cout << "P_ave[" << k << "] = " << P_ave[k]/ ( m+n+o ) << " n = " << n << " m = " << m << " o = " << o << endl;
P_ave[k]/= (m+n+o);
consFile<<k<<" "<<time<<" "<<P_ave[k]<<endl;
if (k==intervals/2) Pressures.push_back(P_ave[k]);
=== modified file 'lib/triangulation/FlowBoundingSphere.h'
--- lib/triangulation/FlowBoundingSphere.h 2010-03-23 11:30:28 +0000
+++ lib/triangulation/FlowBoundingSphere.h 2010-04-07 15:03:09 +0000
@@ -90,7 +90,7 @@
double Permeameter ( RTriangulation& Tri, double P_Inf, double P_Sup, double Section, double DeltaY, char *file );
double Sample_Permeability ( RTriangulation& Tri, double x_Min,double x_Max ,double y_Min,double y_Max,double z_Min,double z_Max, std::string key );
double Compute_HydraulicRadius ( RTriangulation& Tri, Cell_handle cell, int j );
- double PermeameterCurve ( RTriangulation& Tri, char *filename, Real time );
+ double PermeameterCurve ( RTriangulation& Tri, char *filename, Real time, int intervals );
double dotProduct ( Vecteur x, Vecteur y );
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-03-23 11:30:28 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-04-07 15:03:09 +0000
@@ -42,6 +42,7 @@
if ( !triaxialCompressionEngine ) cout << "stress controller engine NOT found" << endl;}
currentStress = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];
+ currentStrain = triaxialCompressionEngine->strain[1];
current_state = triaxialCompressionEngine->currentState;
if ( current_state==3 )
@@ -103,7 +104,7 @@
char *g = file;
timingDeltas->checkpoint("Writing cons_files");
- MaxPressure = flow->PermeameterCurve(flow->T[flow->currentTes].Triangulation(), g, time);
+ MaxPressure = flow->PermeameterCurve(flow->T[flow->currentTes].Triangulation(), g, time, intervals);
if ( Omega::instance().getCurrentIteration() % PermuteInterval == 0 )
{ Update_Triangulation = true; }
@@ -149,6 +150,7 @@
cout << "---------UPDATE PERMEABILITY VALUE--------------" << endl;
if (compute_K) {flow->TOLERANCE=1e-09; K = flow->Sample_Permeability ( flow->T[flow->currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
flow->currentTes=!flow->currentTes;
+ flow->TOLERANCE=Tolerance;
cout << "--------RETRIANGULATION-----------" << endl;
}
// currentTes=flow->currentTes;
@@ -178,6 +180,7 @@
}
else
{
+ Oedometer_Boundary_Conditions();
flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
Update_Triangulation=!Update_Triangulation;
}
@@ -202,14 +205,15 @@
// 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]);
// flow->AddBoundingPlanes ( center, Extent, id );
- flow->x_min = min ( flow->x_min, center[0]-wall_thickness);
- flow->x_max = max ( flow->x_max, center[0]+wall_thickness);
- flow->y_min = min ( flow->y_min, center[1]-wall_thickness);
- flow->y_max = max ( flow->y_max, center[1]+wall_thickness);
- flow->z_min = min ( flow->z_min, center[2]-wall_thickness);
- flow->z_max = max ( flow->z_max, center[2]+wall_thickness);
+ flow->x_min = min ( flow->x_min, center[0]+wall_thickness);
+ flow->x_max = max ( flow->x_max, center[0]-wall_thickness);
+ flow->y_min = min ( flow->y_min, center[1]+wall_thickness);
+ flow->y_max = max ( flow->y_max, center[1]-wall_thickness);
+ flow->z_min = min ( flow->z_min, center[2]+wall_thickness);
+ flow->z_max = max ( flow->z_max, center[2]-wall_thickness);
}
}
@@ -324,7 +328,7 @@
b = cell->vertex ( y )->info().id();
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]-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;
}
}
@@ -359,7 +363,7 @@
b[j] = cell->vertex ( g )->info().id();
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] - wall_thickness;
+ 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++;
}
else if ( first_sph )
@@ -408,7 +412,7 @@
b[j] = cell->vertex ( g )->info().id();
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]-wall_thickness;
+ 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++;
}
else
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.hpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-03-23 11:30:28 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-04-07 15:03:09 +0000
@@ -1,6 +1,6 @@
/*************************************************************************
* Copyright (C) 2009 by Emanuele Catalano *
-* emanuele.catalanog@xxxxxxxxxxx *
+* emanuele.catalano@xxxxxxxxxxx *
* *
* 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. *
@@ -55,9 +55,10 @@
((double,permeability_factor,1.0,"a permability multiplicator"))
((Real,loadFactor,1.1,"Load multiplicator for oedometer test"))
((double, K, 0, "Permeability of the sample"))
- ((bool, ComputeFlow, 1,"If false only triangulation is done, if true flow problem is solved"))
((double, MaxPressure, 0, "Maximal value of water pressure within the sample"))
- ((double, currentStress, 0, "Current value of normal stress"))
+ ((double, currentStress, 0, "Current value of axial stress"))
+ ((double, currentStrain, 0, "Current value of axial strain"))
+ ((int, intervals, 30, "Number of layers for pressure measurements"))
((bool,tess_based_force,true,"true=force computation based on tessalation, false=force computation based on triangulation")),
timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas));
DECLARE_LOGGER;
