yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #08368
[Branch ~yade-dev/yade/trunk] Rev 3041: - compile error
------------------------------------------------------------
revno: 3041
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Wed 2012-03-07 20:22:27 +0100
message:
- compile error
- remove warnings
- add last diffs from Donia
------------- This line and the following will be ignored --------------
modified:
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/KinematicLocalisationAnalyser.cpp
lib/triangulation/PeriodicFlow.cpp
lib/triangulation/PeriodicFlow.hpp
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.hpp
modified:
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/KinematicLocalisationAnalyser.cpp
lib/triangulation/PeriodicFlow.cpp
lib/triangulation/PeriodicFlow.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.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2012-02-14 16:58:21 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2012-03-07 19:22:27 +0000
@@ -21,10 +21,11 @@
using namespace std;
namespace CGT {
-template<class Tesselation>
-class FlowBoundingSphere : public Network<Tesselation>
+template<class _Tesselation>
+class FlowBoundingSphere : public Network<_Tesselation>
{
public:
+ typedef _Tesselation Tesselation;
typedef Network<Tesselation> _N;
DECLARE_TESSELATION_TYPES(Network<Tesselation>)
=== modified file 'lib/triangulation/KinematicLocalisationAnalyser.cpp'
--- lib/triangulation/KinematicLocalisationAnalyser.cpp 2012-02-14 17:04:44 +0000
+++ lib/triangulation/KinematicLocalisationAnalyser.cpp 2012-03-07 19:22:27 +0000
@@ -188,6 +188,7 @@
TS1->from_file(state_file1,/*use bz2?*/ bz2);
TS0->from_file(state_file0,/*use bz2?*/ bz2);
DefToFile(output_file_name);
+ return 0;
}
bool KinematicLocalisationAnalyser::DefToFile(const char* output_file_name)
=== modified file 'lib/triangulation/PeriodicFlow.cpp'
--- lib/triangulation/PeriodicFlow.cpp 2012-03-06 11:49:48 +0000
+++ lib/triangulation/PeriodicFlow.cpp 2012-03-07 19:22:27 +0000
@@ -28,6 +28,7 @@
// #define GS_OPEN_MP //It should never be defined if Yade is not using openmp
#endif
+
namespace CGT {
void PeriodicFlow::ComputeFacetForcesWithCache()
@@ -52,7 +53,6 @@
if (cell->info().isGhost) continue;
for (int k=0;k<4;k++) cell->info().unitForceVectors[k]=nullVect;
for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j))) {
-// if (cell->neighbor(j)->info().isGhost == true) continue;
neighbour_cell = cell->neighbor(j);
const Vecteur& Surfk = cell->info().facetSurfaces[j];
//FIXME : later compute that fluidSurf only once in hydraulicRadius, for now keep full surface not modified in cell->info for comparison with other forces schemes
@@ -508,7 +508,37 @@
cell->info().av_vel() = cell->info().av_vel() /abs(cell->info().volume());
}
}
-
+void PeriodicFlow::ComputeEdgesSurfaces()
+{
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ Edge_normal.clear(); Edge_Surfaces.clear(); Edge_ids.clear(); Edge_HydRad.clear();
+ Finite_edges_iterator ed_it;
+ for ( Finite_edges_iterator ed_it = Tri.finite_edges_begin(); ed_it!=Tri.finite_edges_end();ed_it++ )
+ {
+ Real Rh;
+ if (((ed_it->first)->vertex(ed_it->second)->info().isFictious) && ((ed_it->first)->vertex(ed_it->third)->info().isFictious)) continue;
+ else if (((ed_it->first)->vertex(ed_it->second)->info().isFictious) && ((ed_it->first)->vertex(ed_it->third)->info().isGhost)) continue;
+ else if (((ed_it->first)->vertex(ed_it->second)->info().isGhost) && ((ed_it->first)->vertex(ed_it->third)->info().isFictious)) continue;
+ else if (((ed_it->first)->vertex(ed_it->second)->info().isGhost) && ((ed_it->first)->vertex(ed_it->third)->info().isGhost)) continue;
+ int id1 = (ed_it->first)->vertex(ed_it->second)->info().id();
+ int id2 = (ed_it->first)->vertex(ed_it->third)->info().id();
+ double area = T[currentTes].ComputeVFacetArea(ed_it);
+ Edge_Surfaces.push_back(area);
+ Edge_ids.push_back(pair<int,int>(id1,id2));
+ double radius1 = sqrt((ed_it->first)->vertex(ed_it->second)->point().weight());
+ double radius2 = sqrt((ed_it->first)->vertex(ed_it->third)->point().weight());
+ Vecteur x = (ed_it->first)->vertex(ed_it->third)->point().point() - (ed_it->first)->vertex(ed_it->second)->point().point();
+ Vecteur n = x / sqrt(x.squared_length());
+ Edge_normal.push_back(Vector3r(n[0],n[1],n[2]));
+ double d = x*n - radius1 - radius2;
+ if (radius1<radius2) Rh = d + 0.45 * radius1;
+ else Rh = d + 0.45 * radius2;
+ Edge_HydRad.push_back(Rh);
+// if (DEBUG_OUT) cout<<"id1= "<<id1<<", id2= "<<id2<<", area= "<<area<<", R1= "<<radius1<<", R2= "<<radius2<<" x= "<<x<<", n= "<<n<<", Rh= "<<Rh<<endl;
+
+ }
+ if (DEBUG_OUT)cout << "size of Edge_ids "<< Edge_ids.size()<< ", size of area "<< Edge_Surfaces.size() << ", size of Rh "<< Edge_HydRad.size()<< ", size of normal "<<Edge_normal.size() << endl;
+}
} //namespace CGT
=== modified file 'lib/triangulation/PeriodicFlow.hpp'
--- lib/triangulation/PeriodicFlow.hpp 2012-02-14 16:58:21 +0000
+++ lib/triangulation/PeriodicFlow.hpp 2012-03-07 19:22:27 +0000
@@ -37,6 +37,7 @@
void GaussSeidel();
void DisplayStatistics();
void Average_Relative_Cell_Velocity();
+ void ComputeEdgesSurfaces();
};
}
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2012-03-06 19:22:30 +0000
+++ pkg/dem/FlowEngine.cpp 2012-03-07 19:22:27 +0000
@@ -60,7 +60,7 @@
///Application of vicscous forces
scene->forces.sync();
- if ( viscousShear ) ApplyViscousForces ( solver );
+ if ( viscousShear ) ApplyViscousForces ( *solver );
Finite_vertices_iterator vertices_end = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
for ( Finite_vertices_iterator V_it = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); V_it != vertices_end; V_it++ ) {
@@ -84,6 +84,7 @@
}
template<class Solver>
+
void FlowEngine::BoundaryConditions ( Solver& flow )
{
if ( flow->y_min_id>=0 ) {
@@ -152,7 +153,6 @@
Real FlowEngine::getFlux ( unsigned int cond,Solver& flow )
{
if ( cond>=flow->imposedP.size() ) LOG_ERROR ( "Getting flux with cond higher than imposedP size." );
- RTriangulation& Tri = flow->T[flow->currentTes].Triangulation();
double flux=0;
Cell_handle& cell= flow->IPCells[cond];
for ( int ngb=0;ngb<4;ngb++ ) {
@@ -189,7 +189,6 @@
flow->meanK_LIMIT = meanK_correction;
flow->meanK_STAT = meanK_opt;
flow->permeability_map = permeability_map;
- //flow->key = key;
flow->T[flow->currentTes].Clear();
flow->T[flow->currentTes].max_id=-1;
@@ -526,29 +525,34 @@
{
// flow->ComputeEdgesSurfaces(); //only done in buildTriangulation
if ( Debug ) cout << "Application of viscous forces" << endl;
- if ( Debug ) cout << "Number of edges = " << flow->Edge_ids.size() << endl;
- for ( unsigned int k=0; k<flow->viscousShearForces.size(); k++ ) flow->viscousShearForces[k]=Vector3r::Zero();
+ if ( Debug ) cout << "Number of edges = " << flow.Edge_ids.size() << endl;
+ for ( unsigned int k=0; k<flow.viscousShearForces.size(); k++ ) flow.viscousShearForces[k]=Vector3r::Zero();
- const Tesselation& Tes = flow->T[flow->currentTes];
- for ( int i=0; i< ( int ) flow->Edge_ids.size(); i++ ) {
- int hasFictious= Tes.vertex ( flow->Edge_ids[i].first )->info().isFictious + Tes.vertex ( flow->Edge_ids[i].second )->info().isFictious;
- const shared_ptr<Body>& sph1 = Body::byId ( flow->Edge_ids[i].first, scene );
- const shared_ptr<Body>& sph2 = Body::byId ( flow->Edge_ids[i].second, scene );
+ typedef typename Solver::Tesselation Tesselation;
+ const Tesselation& Tes = flow.T[flow.currentTes];
+ for ( int i=0; i< ( int ) flow.Edge_ids.size(); i++ ) {
+ int hasFictious= Tes.vertex ( flow.Edge_ids[i].first )->info().isFictious + Tes.vertex ( flow.Edge_ids[i].second )->info().isFictious;
+ const shared_ptr<Body>& sph1 = Body::byId ( flow.Edge_ids[i].first, scene );
+ const shared_ptr<Body>& sph2 = Body::byId ( flow.Edge_ids[i].second, scene );
+ Sphere* s1=YADE_CAST<Sphere*> ( sph1->shape.get() );
+ Sphere* s2=YADE_CAST<Sphere*> ( sph2->shape.get() );
+ Vector3r x = sph1->state->pos - sph2->state->pos;
+ Vector3r n = x / sqrt(makeCgVect(x).squared_length());
Vector3r deltaV;
if ( !hasFictious )
- deltaV = ( sph2->state->vel - sph1->state->vel );
+ deltaV = (sph2->state->vel + sph2->state->angVel.cross(s2->radius * n)) - (sph1->state->vel+ sph1->state->angVel.cross(s1->radius * n));
else {
if ( hasFictious==1 ) {//for the fictious sphere, use velocity of the boundary, not of the body
- Vector3r v1 = ( Tes.vertex ( flow->Edge_ids[i].first )->info().isFictious ) ? flow->boundary ( flow->Edge_ids[i].first ).velocity:sph1->state->vel;
- Vector3r v2 = ( Tes.vertex ( flow->Edge_ids[i].second )->info().isFictious ) ? flow->boundary ( flow->Edge_ids[i].second ).velocity:sph2->state->vel;
+ Vector3r v1 = ( Tes.vertex ( flow.Edge_ids[i].first )->info().isFictious ) ? flow.boundary ( flow.Edge_ids[i].first ).velocity:sph1->state->vel + sph1->state->angVel.cross(s1->radius * n);
+ Vector3r v2 = ( Tes.vertex ( flow.Edge_ids[i].second )->info().isFictious ) ? flow.boundary ( flow.Edge_ids[i].second ).velocity:sph2->state->vel + sph2->state->angVel.cross(s2->radius * n);
deltaV = v2-v1;
} else {//both fictious, ignore
deltaV = Vector3r::Zero();
}
}
- deltaV = deltaV - ( flow->Edge_normal[i].dot ( deltaV ) ) *flow->Edge_normal[i];
- Vector3r visc_f = flow->ComputeViscousForce ( deltaV, i );
- if ( Debug ) cout << "la force visqueuse entre " << flow->Edge_ids[i].first << " et " << flow->Edge_ids[i].second << "est " << visc_f << endl;
+ deltaV = deltaV - ( flow.Edge_normal[i].dot ( deltaV ) ) *flow.Edge_normal[i];
+ Vector3r visc_f = flow.ComputeViscousForce ( deltaV, i );
+// if ( Debug ) cout << "la force visqueuse entre " << flow->Edge_ids[i].first << " et " << flow->Edge_ids[i].second << "est " << visc_f << endl;
/// //(1) directement sur le body Yade...
// scene->forces.addForce(flow->Edge_ids[i].first,visc_f);
// scene->forces.addForce(flow->Edge_ids[i].second,-visc_f);
@@ -556,9 +560,10 @@
// Tes.vertex(flow->Edge_ids[i].first)->info().forces=Tes.vertex(flow->Edge_ids[i].first)->info().forces+makeCgVect(visc_f);
// Tes.vertex(flow->Edge_ids[i].second)->info().forces=Tes.vertex(flow->Edge_ids[i].second)->info().forces+makeCgVect(visc_f);
/// //(3) ou dans un vecteur séparé (rapide)
- flow->viscousShearForces[flow->Edge_ids[i].first]+=visc_f;
- flow->viscousShearForces[flow->Edge_ids[i].second]-=visc_f;
+ flow.viscousShearForces[flow.Edge_ids[i].first]+=visc_f;
+ flow.viscousShearForces[flow.Edge_ids[i].second]-=visc_f;
}
+ if(Debug) cout<<"number of viscousShearForce"<<flow.viscousShearForces.size()<<endl;
}
YADE_PLUGIN ( ( FlowEngine ) );
@@ -605,14 +610,16 @@
///Application of vicscous forces
scene->forces.sync();
- if ( viscousShear ) ApplyViscousForces(solver);
+ if ( viscousShear ) ApplyViscousForces(*solver);
Finite_vertices_iterator vertices_end = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
for ( Finite_vertices_iterator V_it = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); V_it != vertices_end; V_it++ ) {
if ( !viscousShear )
scene->forces.addForce ( V_it->info().id(), Vector3r ( ( V_it->info().forces ) [0],V_it->info().forces[1],V_it->info().forces[2] ) );
else
- scene->forces.addForce ( V_it->info().id(), Vector3r ( ( V_it->info().forces ) [0],V_it->info().forces[1],V_it->info().forces[2] ) +solver->viscousShearForces[V_it->info().id() ] );
+ if (V_it->info().isGhost) continue;
+ else
+ scene->forces.addForce ( V_it->info().id(), Vector3r ( ( V_it->info().forces ) [0],V_it->info().forces[1],V_it->info().forces[2] ) +solver->viscousShearForces[V_it->info().id() ] );
}
///End Compute flow and forces
// timingDeltas->checkpoint("Applying Forces");
@@ -710,46 +717,6 @@
}
-
-// void PeriodicFlowEngine::ApplyViscousForces()
-// {
-// // flow->ComputeEdgesSurfaces(); //only done in buildTriangulation
-// if ( Debug ) cout << "Application of viscous forces" << endl;
-// if ( Debug ) cout << "Number of edges = " << solver->Edge_ids.size() << endl;
-// for ( unsigned int k=0; k<solver->viscousShearForces.size(); k++ ) solver->viscousShearForces[k]=Vector3r::Zero();
-//
-// const Tesselation& Tes = solver->T[solver->currentTes];
-// for ( int i=0; i< ( int ) solver->Edge_ids.size(); i++ ) {
-// int hasFictious= Tes.vertex ( solver->Edge_ids[i].first )->info().isFictious + Tes.vertex ( solver->Edge_ids[i].second )->info().isFictious;
-// const shared_ptr<Body>& sph1 = Body::byId ( solver->Edge_ids[i].first, scene );
-// const shared_ptr<Body>& sph2 = Body::byId ( solver->Edge_ids[i].second, scene );
-// Vector3r deltaV;
-// if ( !hasFictious )
-// deltaV = ( sph2->state->vel - sph1->state->vel );
-// else {
-// if ( hasFictious==1 ) {//for the fictious sphere, use velocity of the boundary, not of the body
-// Vector3r v1 = ( Tes.vertex ( solver->Edge_ids[i].first )->info().isFictious ) ? solver->boundary ( solver->Edge_ids[i].first ).velocity:sph1->state->vel;
-// Vector3r v2 = ( Tes.vertex ( solver->Edge_ids[i].second )->info().isFictious ) ? solver->boundary ( solver->Edge_ids[i].second ).velocity:sph2->state->vel;
-// deltaV = v2-v1;
-// } else {//both fictious, ignore
-// deltaV = Vector3r::Zero();
-// }
-// }
-// deltaV = deltaV - ( solver->Edge_normal[i].dot ( deltaV ) ) *solver->Edge_normal[i];
-// Vector3r visc_f = solver->ComputeViscousForce ( deltaV, i );
-// if ( Debug ) cout << "la force visqueuse entre " << solver->Edge_ids[i].first << " et " << solver->Edge_ids[i].second << "est " << visc_f << endl;
-// /// //(1) directement sur le body Yade...
-// // scene->forces.addForce(flow->Edge_ids[i].first,visc_f);
-// // scene->forces.addForce(flow->Edge_ids[i].second,-visc_f);
-// /// //(2) ou dans CGAL? On a le choix (on pourrait même avoir info->viscousF pour faire la différence entre les deux types de forces... mais ça prend un peu plus de mémoire et de temps de calcul)
-// // Tes.vertex(flow->Edge_ids[i].first)->info().forces=Tes.vertex(flow->Edge_ids[i].first)->info().forces+makeCgVect(visc_f);
-// // Tes.vertex(flow->Edge_ids[i].second)->info().forces=Tes.vertex(flow->Edge_ids[i].second)->info().forces+makeCgVect(visc_f);
-// /// //(3) ou dans un vecteur séparé (rapide)
-// solver->viscousShearForces[solver->Edge_ids[i].first]+=visc_f;
-// solver->viscousShearForces[solver->Edge_ids[i].second]-=visc_f;
-// }
-// }
-
void PeriodicFlowEngine::locateCell ( Cell_handle baseCell, unsigned int& index, unsigned int count)
{
if (count>10) LOG_ERROR("More than 10 attempts to locate a cell, duplicateThreshold may be too small, resulting in periodicity inconsistencies.");
@@ -757,13 +724,12 @@
//already located, return FIXME: is inline working correctly? else move this test outside the function, just before the calls
if ( base_info.index>0 || base_info.isGhost ) return;
RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
- Tesselation& Tes = solver->T[solver->currentTes];
Vector3r center ( 0,0,0 );
Vector3i period;
if (baseCell->info().fictious()==0)
for ( int k=0;k<4;k++ ) center+= 0.25*makeVector3r (baseCell->vertex(k)->point());
else {
- Real boundPos; int coord;
+ Real boundPos=0; int coord=0;
for ( int k=0;k<4;k++ ) {
if ( !baseCell->vertex ( k )->info().isFictious ) center+= 0.3333333333*makeVector3r ( baseCell->vertex ( k )->point() );
else {
@@ -797,19 +763,19 @@
}
}
-Vector3r PeriodicFlowEngine::MeanVelocity()
+Vector3r PeriodicFlowEngine::meanVelocity()
{
solver->Average_Relative_Cell_Velocity();
- Vector3r MeanVel ( 0,0,0 );
+ Vector3r meanVel ( 0,0,0 );
Real volume=0;
Finite_cells_iterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
for ( Finite_cells_iterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
if ( cell->info().isGhost ) continue;
for ( int i=0;i<3;i++ )
- MeanVel[i]=MeanVel[i]+ ( ( cell->info().av_vel() ) [i] * abs ( cell->info().volume() ) );
+ meanVel[i]=meanVel[i]+ ( ( cell->info().av_vel() ) [i] * abs ( cell->info().volume() ) );
volume+=abs ( cell->info().volume() );
}
- return ( MeanVel/volume );
+ return ( meanVel/volume );
}
void PeriodicFlowEngine::UpdateVolumes ()
@@ -824,10 +790,8 @@
Real totDVol=0;
Real totVol0=0;
Real totVol1=0;
- Real totVol2=0;
- Real totVol3=0;
- for ( Finite_cells_iterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
+ for ( Finite_cells_iterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
if ( cell->info().isGhost ) continue;
switch ( cell->info().fictious() ) {
case ( 1 ) :
@@ -850,7 +814,6 @@
cell->info().volume() = newVol;
}
if ( Debug ) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<" "<< totVol0<<" "<< totVol1<<endl;
-
}
void PeriodicFlowEngine::Initialize_volumes ()
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2012-03-06 19:22:30 +0000
+++ pkg/dem/FlowEngine.hpp 2012-03-07 19:22:27 +0000
@@ -99,7 +99,7 @@
unsigned int _imposePressure(Vector3r pos, Real p) {return imposePressure(pos,p,solver);}
void _setImposedPressure(unsigned int cond, Real p) {setImposedPressure(cond,p,solver);}
void _clearImposedPressure() {clearImposedPressure(solver);}
- Real _getFlux(unsigned int cond) {getFlux(cond,solver);}
+ Real _getFlux(unsigned int cond) {return getFlux(cond,solver);}
virtual ~FlowEngine();
@@ -238,10 +238,10 @@
void Initialize_volumes ();
void UpdateVolumes ();
Real Volume_cell (Cell_handle cell);
+
Real Volume_cell_single_fictious (Cell_handle cell);
-// void ApplyViscousForces();
inline void locateCell(Cell_handle baseCell, unsigned int& index, unsigned int count=0);
- Vector3r MeanVelocity();
+ Vector3r meanVelocity();
virtual ~PeriodicFlowEngine();
@@ -266,7 +266,7 @@
eps_vol_max=Eps_Vol_Cumulative=retriangulationLastIter=0;
ReTrg=1;
,
- .def("MeanVelocity",&PeriodicFlowEngine::MeanVelocity,"measure the mean velocity in the period")
+ .def("meanVelocity",&PeriodicFlowEngine::meanVelocity,"measure the mean velocity in the period")
// .def("imposeFlux",&FlowEngine::_imposeFlux,(python::arg("pos"),python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
.def("imposePressure",&PeriodicFlowEngine::_imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
// .def("setImposedPressure",&FlowEngine::_setImposedPressure,(python::arg("cond"),python::arg("p")),"Set pressure value at the point indexed 'cond'.")
