yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07514
[Branch ~yade-dev/yade/trunk] Rev 2831: - Set default values of walls ids
------------------------------------------------------------
revno: 2831
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Fri 2011-04-22 11:09:28 +0200
message:
- Set default values of walls ids
- Fix cells volume functions
- Add a vertex accessor Tesselation.vertex(body_id)
- Various improvements
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/Network.cpp
lib/triangulation/Tesselation.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-04-06 17:22:14 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2011-04-22 09:09:28 +0000
@@ -44,8 +44,9 @@
typedef vector<double> VectorR;
//! Use this factor, or minLength, to reduce max permeability values (see usage below))
-const double MAXK_DIV_KMEAN = 1;
-const double minLength = 0.20;//percentage of mean rad
+const double MAXK_DIV_KMEAN = 2;
+const double MINK_DIV_KMEAN = 0.05;
+const double minLength = 0.02;//percentage of mean rad
//! Factors including the effect of 1/2 symmetry in hydraulic radii
const Real multSym1 = 1/pow(2,0.25);
@@ -518,8 +519,13 @@
RTriangulation& Tri = T[currentTes].Triangulation();
Finite_cells_iterator cell_end = Tri.finite_cells_end();
Vecteur nullVect(0,0,0);
+ static vector<Vecteur> oldForces;
+ if (oldForces.size()<=Tri.number_of_vertices()) oldForces.resize(Tri.number_of_vertices()+1);
//reset forces
- for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) v->info().forces=nullVect;
+ for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+ if (noCache) {oldForces[v->info().id()]=nullVect; v->info().forces=nullVect;}
+ else {oldForces[v->info().id()]=v->info().forces; v->info().forces=nullVect;}
+ }
Cell_handle neighbour_cell;
Vertex_handle mirror_vertex;
@@ -567,6 +573,7 @@
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++)
for (int yy=0;yy<4;yy++) cell->vertex(yy)->info().forces = cell->vertex(yy)->info().forces + cell->info().unitForceVectors[yy]*cell->info().p();
noCache=false;//cache should always be defined after execution of this function
+ for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) v->info().forces = 0*oldForces[v->info().id()]+1*v->info().forces;
if (DEBUG_OUT) {
cout << "Facet cached scheme" <<endl;
Vecteur TotalForce = nullVect;
@@ -768,7 +775,8 @@
S0=checkSphereFacetOverlap(v0,v1,v2);
if (S0==0) S0=checkSphereFacetOverlap(v1,v2,v0);
if (S0==0) S0=checkSphereFacetOverlap(v2,v0,v1);
- fluidArea=area-crossSections[0]-crossSections[1]-crossSections[2]+S0;
+ //take absolute value, since in rare cases the surface can be negative (overlaping spheres)
+ fluidArea=abs(area-crossSections[0]-crossSections[1]-crossSections[2]+S0);
cell->info().facetFluidSurfacesRatio[j]=fluidArea/area;
k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);}
@@ -790,7 +798,7 @@
}
cell->info().isvisited = !ref;
}
- cout<<"surfneg est "<<surfneg<<endl;
+// cout<<"surfneg est "<<surfneg<<endl;
meanK /= pass;
meanRadius /= pass;
meanDistance /= pass;
@@ -809,8 +817,9 @@
neighbour_cell = cell->neighbor(j);
if (!Tri.is_infinite(neighbour_cell) && neighbour_cell->info().isvisited==ref) {
pass++;
- (cell->info().k_norm())[j] = min((cell->info().k_norm())[j], maxKdivKmean*meanK);
+ (cell->info().k_norm())[j] = max(MINK_DIV_KMEAN*meanK ,min((cell->info().k_norm())[j], maxKdivKmean*meanK));
(neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]=(cell->info().k_norm())[j];
+// cout<<(cell->info().k_norm())[j]<<endl;
// kFile << (cell->info().k_norm())[j] << endl;
}
}cell->info().isvisited = !ref;
=== modified file 'lib/triangulation/Network.cpp'
--- lib/triangulation/Network.cpp 2011-04-06 12:17:47 +0000
+++ lib/triangulation/Network.cpp 2011-04-22 09:09:28 +0000
@@ -419,13 +419,13 @@
Tesselation& Tes = T[currentTes];
//FIXME: Id's order in boundsIds is done according to the enumerotation of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
- y_min_id = Tes.Max_id() + 1;
+ y_min_id = Tes.Max_id() + 2;
boundsIds[0]=&y_min_id;
- y_max_id = Tes.Max_id() + 2;
+ y_max_id = Tes.Max_id() + 3;
boundsIds[1]=&y_max_id;
- x_min_id = Tes.Max_id() + 3;
+ x_min_id = Tes.Max_id() + 0;
boundsIds[2]=&x_min_id;
- x_max_id = Tes.Max_id() + 4;
+ x_max_id = Tes.Max_id() + 1;
boundsIds[3]=&x_max_id;
z_min_id = Tes.Max_id() + 5;
boundsIds[4]=&z_max_id;
=== modified file 'lib/triangulation/Tesselation.h'
--- lib/triangulation/Tesselation.h 2010-11-02 16:23:44 +0000
+++ lib/triangulation/Tesselation.h 2011-04-22 09:09:28 +0000
@@ -75,6 +75,7 @@
void ComputeVolumes (void);//Compute volume each voronoi cell
void ComputePorosity (void);//Compute volume and porosity of each voronoi cell
inline Real& Volume (unsigned int id) { return vertexHandles[id]->info().v(); }
+ inline Vertex_handle& vertex (unsigned int id) { return vertexHandles[id]; }
Finite_cells_iterator finite_cells_begin(void);// {return Tri->finite_cells_begin();}
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2011-04-18 18:10:20 +0000
+++ pkg/dem/FlowEngine.cpp 2011-04-22 09:09:28 +0000
@@ -21,7 +21,6 @@
CREATE_LOGGER (FlowEngine);
-
CGT::Vecteur makeCgVect (const Vector3r& yv) {return CGT::Vecteur(yv[0],yv[1],yv[2]);}
CGT::Point makeCgPoint (const Vector3r& yv) {return CGT::Point(yv[0],yv[1],yv[2]);}
@@ -50,7 +49,7 @@
vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();
for (;itFirst!=itLast; ++itFirst) {
if ((*itFirst)->getClassName() == "TriaxialCompressionEngine") {
- cout << "stress controller engine found - FlowEngine" << endl;
+// cout << "stress controller engine found - FlowEngine" << endl;
triaxialCompressionEngine = YADE_PTR_CAST<TriaxialCompressionEngine> (*itFirst);}}
if (!triaxialCompressionEngine) cout << "stress controller engine NOT found" << endl;
}
@@ -61,13 +60,11 @@
if (first) Build_Triangulation(P_zero);
timingDeltas->checkpoint("Triangulating");
-
if (!first) {
eps_vol_max=0.f;
UpdateVolumes ( );
-
Eps_Vol_Cumulative += eps_vol_max;
- if ((Eps_Vol_Cumulative > EpsVolPercent_RTRG || retriangulationLastIter>1000) && retriangulationLastIter>100) {
+ if ((Eps_Vol_Cumulative > EpsVolPercent_RTRG || retriangulationLastIter>10000) && retriangulationLastIter>10) {
Update_Triangulation = true;
Eps_Vol_Cumulative=0;
retriangulationLastIter=0;
@@ -133,6 +130,7 @@
flow->SliceField(f);
}
// if (save_vtk) {flow->save_vtk_file();}
+ timingDeltas->checkpoint("Writing files");
}
// if ( scene->iter % PermuteInterval == 0 )
// { Update_Triangulation = true; }
@@ -147,6 +145,7 @@
wall_down_y = flow->y_min;}
if (liquefaction) flow->Measure_Pore_Pressure(wall_up_y, wall_down_y);
first=false;
+ timingDeltas->checkpoint("Ending");
// if(id_sphere>=0) flow->Average_Fluid_Velocity_On_Sphere(id_sphere);
// }
@@ -291,13 +290,22 @@
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);
- flow->y_min_id=triaxialCompressionEngine->wall_bottom_id;
- flow->y_max_id=triaxialCompressionEngine->wall_top_id;
- flow->x_max_id=triaxialCompressionEngine->wall_right_id;
- flow->x_min_id=triaxialCompressionEngine->wall_left_id;
- flow->z_min_id=triaxialCompressionEngine->wall_back_id;
- flow->z_max_id=triaxialCompressionEngine->wall_front_id;
-
+ if (triaxialCompressionEngine) {
+ flow->y_min_id=triaxialCompressionEngine->wall_bottom_id;
+ flow->y_max_id=triaxialCompressionEngine->wall_top_id;
+ flow->x_max_id=triaxialCompressionEngine->wall_right_id;
+ flow->x_min_id=triaxialCompressionEngine->wall_left_id;
+ flow->z_min_id=triaxialCompressionEngine->wall_back_id;
+ flow->z_max_id=triaxialCompressionEngine->wall_front_id;
+ } else {
+ flow->y_min_id=wallBottomId;
+ flow->y_max_id=wallTopId;
+ flow->x_max_id=wallRightId;
+ flow->x_min_id=wallLeftId;
+ flow->z_min_id=wallBackId;
+ flow->z_max_id=wallFrontId;
+ }
+
flow->boundary ( flow->y_min_id ).useMaxMin = BOTTOM_Boundary_MaxMin;
flow->boundary ( flow->y_max_id ).useMaxMin = TOP_Boundary_MaxMin;
flow->boundary ( flow->x_max_id ).useMaxMin = RIGHT_Boundary_MaxMin;
@@ -305,13 +313,13 @@
flow->boundary ( flow->z_max_id ).useMaxMin = FRONT_Boundary_MaxMin;
flow->boundary ( flow->z_min_id ).useMaxMin = BACK_Boundary_MaxMin;
- //FIXME: Id's order in boundsIds is done according to the enumerotation of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
- flow->boundsIds[0]= &flow->y_min_id;
- flow->boundsIds[1]= &flow->y_max_id;
- flow->boundsIds[2]= &flow->x_min_id;
- flow->boundsIds[3]= &flow->x_max_id;
- flow->boundsIds[4]= &flow->z_max_id;
- flow->boundsIds[5]= &flow->z_min_id;
+ //FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
+ flow->boundsIds[0]= &flow->x_min_id;
+ flow->boundsIds[1]= &flow->x_max_id;
+ flow->boundsIds[2]= &flow->y_min_id;
+ flow->boundsIds[3]= &flow->y_max_id;
+ flow->boundsIds[4]= &flow->z_min_id;
+ flow->boundsIds[5]= &flow->z_max_id;
wall_thickness = triaxialCompressionEngine->thickness;
@@ -377,6 +385,10 @@
void FlowEngine::Initialize_volumes ()
{
+ CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
+ CGT::Vecteur Zero(0,0,0);
+ for (CGT::Finite_vertices_iterator V_it = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero;
+
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++ )
{
@@ -386,6 +398,7 @@
case ( 1 ) : cell->info().volume() = Volume_cell_single_fictious ( cell ); break;
case ( 2 ) : cell->info().volume() = Volume_cell_double_fictious ( cell ); break;
case ( 3 ) : cell->info().volume() = Volume_cell_triple_fictious ( cell ); break;
+ default: break;
}
}
if (Debug) cout << "Volumes initialised." << endl;
@@ -454,125 +467,102 @@
if (Debug) cout << "Updating volumes.............." << endl;
Real invDeltaT = 1/scene->dt;
CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
- double newVol=0; double dVol;
+ double newVol, dVol;
eps_vol_max=0;
+ Real totVol=0; Real totDVol=0; Real totVol0=0; Real totVol1=0; Real totVol2=0; Real totVol3=0;
for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
switch ( cell->info().fictious()) {
- case ( 3 ):
- newVol= Volume_cell_triple_fictious ( cell );
- break;
- case ( 2 ) :
- newVol = Volume_cell_double_fictious ( cell );
- break;
- case ( 1 ) :
- newVol = Volume_cell_single_fictious ( cell );
- break;
- case ( 0 ) :
- newVol = Volume_cell ( cell );
- break;
+ case ( 3 ): newVol= Volume_cell_triple_fictious ( cell ); totVol3+=newVol; break;
+ case ( 2 ) : newVol = Volume_cell_double_fictious ( cell ); totVol2+=newVol; break;
+ case ( 1 ) : newVol = Volume_cell_single_fictious ( cell ); totVol1+=newVol; break;
+ case ( 0 ) : newVol = Volume_cell ( cell ); totVol0+=newVol; break;
+ default: newVol = 0; break;
}
+ totVol+=newVol;
dVol=cell->info().volumeSign*(newVol - cell->info().volume());
+ totDVol+=dVol;
eps_vol_max = max(eps_vol_max, abs(dVol/newVol));
cell->info().dv() = (!cell->info().Pcondition)?dVol*invDeltaT:0;
cell->info().volume() = newVol;
// if (Debug) cerr<<"v/dv : "<<cell->info().volume()<<" "<<cell->info().dv()<<" ("<<cell->info().fictious()<<")"<<endl;
}
+ if (Debug) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<endl;
}
Real FlowEngine::Volume_cell_single_fictious ( CGT::Cell_handle cell )
{
- Real V[3][3];
+ Vector3r V[3];
int b=0;
int w=0;
cell->info().volumeSign=1;
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;
+ 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);
+ V[w]=sph->state->pos;
+ w++;}
+ else {
+ b = cell->vertex ( y )->info().id();
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];
- }
+ if (!flow->boundary(b).useMaxMin) Wall_coordinate = wll->state->pos[flow->boundary(b).coordinate]+(flow->boundary(b).normal[flow->boundary(b).coordinate])*wall_thickness/2;
+ else Wall_coordinate = flow->boundary(b).p[flow->boundary(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 );
+
+ Real Volume = 0.5*((V[0]-V[1]).cross(V[0]-V[2]))[flow->boundary(b).coordinate] * ( 0.33333333333* ( V[0][flow->boundary(b).coordinate]+ V[1][flow->boundary(b).coordinate]+ V[2][flow->boundary(b).coordinate] ) - Wall_coordinate );
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 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};
+ Vector3r A=Vector3r::Zero(), AS=Vector3r::Zero(),B=Vector3r::Zero(), BS=Vector3r::Zero();
cell->info().volumeSign=1;
int b[2];
+ int coord[2];
Real Wall_coordinate[2];
int j=0;
bool first_sph=true;
for ( int g=0;g<4;g++ )
{
- if ( cell->vertex ( g )->info().isFictious )
+ if ( cell->vertex(g)->info().isFictious )
{
- b[j] = cell->vertex ( g )->info().id()-flow->id_offset;
+ b[j] = cell->vertex (g)->info().id();
+ coord[j]=flow->boundary(b[j]).coordinate;
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];
+ if (!flow->boundary(b[j]).useMaxMin) Wall_coordinate[j] = wll->state->pos[coord[j]] +(flow->boundary(b[j]).normal[coord[j]])*wall_thickness/2;
+ else Wall_coordinate[j] = flow->boundary(b[j]).p[coord[j]];
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]; }
- }
+ else if ( first_sph ){
+ const shared_ptr<Body>& sph1 = Body::byId (cell->vertex ( g )->info().id(), scene );
+ A=AS=/*AT=*/(sph1->state->pos);
+ first_sph=false;}
+ else { const shared_ptr<Body>& sph2 = Body::byId(cell->vertex ( g )->info().id(), scene );
+ B=BS=/*BT=*/(sph2->state->pos);}
}
-
- 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;}
-
- 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 );
+ AS[coord[0]]=BS[coord[0]] = Wall_coordinate[0];
+
+ //first pyramid with triangular base (A,B,BS)
+ Real Vol1=0.5*((A-BS).cross(B-BS))[coord[1]]*(0.333333333*(2*B[coord[1]]+A[coord[1]])-Wall_coordinate[1]);
+ //second pyramid with triangular base (A,AS,BS)
+ Real Vol2=0.5*((AS-BS).cross(A-BS))[coord[1]]*(0.333333333*(B[coord[1]]+2*A[coord[1]])-Wall_coordinate[1]);
+ return abs ( Vol1+Vol2 );
}
Real FlowEngine::Volume_cell_triple_fictious ( CGT::Cell_handle cell)
{
- Real A[3]={0, 0, 0}, AS[3]={0, 0, 0}, AT[3]={0, 0, 0}, AW[3]={0, 0, 0};
+ Vector3r A;
int b[3];
+ int coord[3];
Real Wall_coordinate[3];
int j=0;
cell->info().volumeSign=1;
@@ -581,30 +571,20 @@
{
if ( cell->vertex ( g )->info().isFictious )
{
- b[j] = cell->vertex ( g )->info().id()-flow->id_offset;
+ b[j] = cell->vertex ( g )->info().id();
+ coord[j]=flow->boundary(b[j]).coordinate;
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];
+ if (!flow->boundary(b[j]).useMaxMin) Wall_coordinate[j] = wll->state->pos[coord[j]] + (flow->boundary(b[j]).normal[coord[j]])*wall_thickness/2;
+ else Wall_coordinate[j] = flow->boundary(b[j]).p[coord[j]];
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];}
+ A=(sph->state->pos);
}
}
-
- 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;
-
+ Real Volume = (A[coord[0]]-Wall_coordinate[0])*(A[coord[1]]-Wall_coordinate[1])*(A[coord[2]]-Wall_coordinate[2]);
return abs ( Volume );
}
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2011-04-06 16:37:32 +0000
+++ pkg/dem/FlowEngine.hpp 2011-04-22 09:09:28 +0000
@@ -20,6 +20,8 @@
#define FlowSolver CGT::FlowBoundingSphere
#endif
+
+
class FlowEngine : public PartialEngine
{
private:
@@ -58,6 +60,7 @@
python::list getConstrictions() {
vector<Real> csd=flow->getConstrictions(); python::list pycsd;
for (unsigned int k=0;k<csd.size();k++) pycsd.append(csd[k]); return pycsd;}
+ Vector3r fluidForce(int id_sph) {const CGT::Vecteur& f=flow->T[flow->currentTes].vertex(id_sph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
virtual ~FlowEngine();
@@ -91,7 +94,7 @@
((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"))
+ ((double,permeability_factor,0,,"permability multiplier"))
((double,viscosity,1.0,,"viscosity of fluid"))
((Real,loadFactor,1.1,,"Load multiplicator for oedometer test"))
((double, K, 0,, "Permeability of the sample"))
@@ -115,6 +118,12 @@
((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"))
+ ((int, wallTopId,3,,"Id of top boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
+ ((int, wallBottomId,2,,"Id of bottom boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
+ ((int, wallFrontId,5,,"Id of front boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
+ ((int, wallBackId,4,,"Id of back boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
+ ((int, wallLeftId,0,,"Id of left boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
+ ((int, wallRightId,1,,"Id of right boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
((Vector3r, id_force, 0,, "Fluid force acting on sphere with id=flow.id_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"))
@@ -132,6 +141,7 @@
.def("getConstrictions",&FlowEngine::getConstrictions,"Get the list of constrictions (inscribed circle) for all finite facets.")
.def("saveVtk",&FlowEngine::saveVtk,"Save pressure field in vtk format.")
.def("AvFlVelOnSph",&FlowEngine::AvFlVelOnSph,(python::arg("Id_sph")),"Compute a sphere-centered average fluid velocity")
+ .def("fluidForce",&FlowEngine::fluidForce,(python::arg("Id_sph")),"Return the fluid force on sphere Id_sph.")
)
DECLARE_LOGGER;
};
