yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07590
[Branch ~yade-dev/yade/trunk] Rev 2852: - remove many warnings
------------------------------------------------------------
revno: 2852
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Tue 2011-05-03 22:03:21 +0200
message:
- remove many warnings
- handle viscous shear force gracefully
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/Network.hpp
lib/triangulation/Tesselation.h
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-05-03 10:33:29 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2011-05-03 20:03:21 +0000
@@ -213,7 +213,7 @@
clock.top("GaussSeidel");
/** END GAUSS SEIDEL */
- char* file ="Permeability";
+ const char* file ="Permeability";
ks = Permeameter(boundary(y_min_id).value, boundary(y_max_id).value, SectionArea, Height, file);
// K_sigma << K_opt_factor << " " << ks << " "<< Iterations << endl;
clock.top("Permeameter");
@@ -222,7 +222,7 @@
clock.top("Compute_Forces");
/** VISUALISATION FILES **/
- // save_vtk_file ( );
+ // saveVtk ( );
// PressureProfile ( );
MGPost();
@@ -392,7 +392,7 @@
}}
}
-vector<Real> FlowBoundingSphere::Average_Fluid_Velocity_On_Sphere(int Id_sph)
+vector<Real> FlowBoundingSphere::Average_Fluid_Velocity_On_Sphere(unsigned int Id_sph)
{
Average_Relative_Cell_Velocity();
RTriangulation& Tri = T[currentTes].Triangulation();
@@ -408,7 +408,7 @@
for ( Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++ )
{
if (cell->info().fictious()==0){
- for (int i=0;i<4;i++){
+ for (unsigned int i=0;i<4;i++){
if (cell->vertex(i)->info().id()==Id_sph){
VelocityVolumes = VelocityVolumes + cell->info().av_vel()*cell->info().volume();
Volumes = Volumes + cell->info().volume();}}}}
@@ -1119,7 +1119,7 @@
// if (j % 10 == 0) {
// cout << "pmax " << p_max << "; pmoy : " << p_moy << "; dpmax : " << dp_max << endl;
// cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;}
-// // save_vtk_file ( Tri );
+// // saveVtk ( Tri );
// }
// if (DEBUG_OUT) {cout << "pmax " << p_max << "; pmoy : " << p_moy << endl; cout << "iteration " << j <<"; erreur : " << dp_max/p_max << " tolerance " << tolerance<<endl;}
#ifdef GS_OPEN_MP
@@ -1156,7 +1156,7 @@
-double FlowBoundingSphere::Permeameter(double P_Inf, double P_Sup, double Section, double DeltaY, char *file)
+double FlowBoundingSphere::Permeameter(double P_Inf, double P_Sup, double Section, double DeltaY, const char *file)
{
RTriangulation& Tri = T[currentTes].Triangulation();
std::ofstream kFile(file, std::ios::out);
@@ -1171,7 +1171,7 @@
Tesselation::VCell_iterator cell_up_end = Tri.incident_cells(T[currentTes].vertexHandles[y_max_id],cells_it);
for (Tesselation::VCell_iterator it = tmp_cells.begin(); it != cell_up_end; it++)
{
- Cell_handle& cell = *it;{for (int j2=0; j2<4; j2++) {if (!cell->neighbor(j2)->info().Pcondition){
+ Cell_handle& cell = *it;{for (int j2=0; j2<4; j2++) {/*if (!cell->neighbor(j2)->info().Pcondition)*/{
if ((cell->neighbor(j2)->info().p()!=cell->neighbor(j2)->info().p()) && Tri.is_infinite(cell->neighbor(j2))) cout << "oooooooooooooooh" << endl;
Q1 = Q1 + (cell->neighbor(j2)->info().k_norm())[Tri.mirror_index(cell, j2)]* (cell->neighbor(j2)->info().p()-cell->info().p());
cellQ1+=1;
@@ -1183,7 +1183,7 @@
Tesselation::VCell_iterator cell_down_end = Tri.incident_cells(T[currentTes].vertexHandles[y_min_id],cells_it);
for (Tesselation::VCell_iterator it = tmp_cells.begin(); it != cell_down_end; it++)
{
- Cell_handle& cell = *it;{for (int j2=0; j2<4; j2++) {if (!cell->neighbor(j2)->info().Pcondition){
+ Cell_handle& cell = *it;{for (int j2=0; j2<4; j2++) {/*if (!cell->neighbor(j2)->info().Pcondition)*/{
if ((cell->neighbor(j2)->info().p()!=cell->neighbor(j2)->info().p()) && Tri.is_infinite(cell->neighbor(j2))) cout << "oooooooooooooooh2" << endl;
Q2 = Q2 + (cell->neighbor(j2)->info().k_norm())[Tri.mirror_index(cell, j2)]* (cell->info().p()-cell->neighbor(j2)->info().p());
cellQ2+=1;
@@ -1275,7 +1275,7 @@
num_particles = real;
}
-void FlowBoundingSphere::save_vtk_file()
+void FlowBoundingSphere::saveVtk()
{
RTriangulation& Tri = T[currentTes].Triangulation();
static unsigned int number = 0;
@@ -1499,7 +1499,7 @@
return false;
}
-void FlowBoundingSphere::SliceField(char *filename)
+void FlowBoundingSphere::SliceField(const char *filename)
{
/** Pressure field along one cutting plane **/
RTriangulation& Tri = T[currentTes].Triangulation();
@@ -1589,7 +1589,7 @@
void FlowBoundingSphere::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++ )
{
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2011-05-03 10:33:29 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2011-05-03 20:03:21 +0000
@@ -22,12 +22,8 @@
#include "Vue3D.h" //FIXME implicit dependencies will look for this class (out of tree) even ifndef XVIEW
#endif
-
-
using namespace std;
-
-
namespace CGT{
class FlowBoundingSphere : public Network
@@ -58,6 +54,7 @@
vector <pair<int,int> > Edge_ids;
vector <Real> Edge_HydRad;
vector <Vector3r> Edge_normal;
+ vector <Vector3r> viscousShearForces;
void mplot ( char *filename);
void Localize();
@@ -87,13 +84,13 @@
/// Define forces spliting drag and buoyancy terms
void ComputeFacetForces();
void ComputeFacetForcesWithCache();
- void save_vtk_file ( );
+ void saveVtk ( );
void MGPost ( );
#ifdef XVIEW
void Dessine_Triangulation ( Vue3D &Vue, RTriangulation &T );
void Dessine_Short_Tesselation ( Vue3D &Vue, Tesselation &Tes );
#endif
- double Permeameter ( double P_Inf, double P_Sup, double Section, double DeltaY, char *file );
+ double Permeameter ( double P_Inf, double P_Sup, double Section, double DeltaY, const char *file );
double Sample_Permeability( double& x_Min,double& x_Max ,double& y_Min,double& y_Max,double& z_Min,double& z_Max, string key);
double Compute_HydraulicRadius (Cell_handle cell, int j );
double PressureProfile ( char *filename, Real& time, int& intervals );
@@ -116,7 +113,7 @@
bool isInsideSphere ( double& x, double& y, double& z );
- void SliceField (char *filename);
+ void SliceField (const char *filename);
void ComsolField();
void Interpolate ( Tesselation& Tes, Tesselation& NewTes );
@@ -125,7 +122,7 @@
void ApplySinusoidalPressure(RTriangulation& Tri, double Amplitude, double Average_Pressure, double load_intervals);
bool isOnSolid (double X, double Y, double Z);
void Measure_Pore_Pressure(double Wall_up_y, double Wall_down_y);
- vector<Real> Average_Fluid_Velocity_On_Sphere(int Id_sph);
+ vector<Real> Average_Fluid_Velocity_On_Sphere(unsigned int Id_sph);
//Solver?
int useSolver;//(0 : GaussSeidel, 1 : TAUCS, 2 : PARDISO)
=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp 2011-02-17 17:43:37 +0000
+++ lib/triangulation/Network.hpp 2011-05-03 20:03:21 +0000
@@ -25,6 +25,7 @@
{
Point p;
Vecteur normal;
+ Vector3r velocity;
int coordinate;
bool flowCondition;//flowCondition=0, pressure is imposed // flowCondition=1, flow is imposed
Real value;
=== modified file 'lib/triangulation/Tesselation.h'
--- lib/triangulation/Tesselation.h 2011-05-02 09:17:49 +0000
+++ lib/triangulation/Tesselation.h 2011-05-03 20:03:21 +0000
@@ -78,7 +78,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]; }
+ inline const Vertex_handle& vertex (unsigned int id) const { return vertexHandles[id]; }
Finite_cells_iterator finite_cells_begin(void);// {return Tri->finite_cells_begin();}
=== modified file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h 2011-05-03 10:33:29 +0000
+++ lib/triangulation/def_types.h 2011-05-03 20:03:21 +0000
@@ -13,7 +13,6 @@
#include <CGAL/number_utils.h>
#include <boost/static_assert.hpp>
-
// #include<yade/lib/base/Math.hpp> // typedef for Real
namespace CGT{
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2011-05-03 10:33:29 +0000
+++ pkg/dem/FlowEngine.cpp 2011-05-03 20:03:21 +0000
@@ -62,15 +62,17 @@
timingDeltas->checkpoint("Triangulating");
UpdateVolumes ( );
if (!first) {
- eps_vol_max=0.f;
+// eps_vol_max=0.f;//huh? in that case Eps_Vol_Cumulative will always be zero
Eps_Vol_Cumulative += eps_vol_max;
- if ((Eps_Vol_Cumulative > EpsVolPercent_RTRG || retriangulationLastIter>10000) && retriangulationLastIter>10) {
+ if ((Eps_Vol_Cumulative > EpsVolPercent_RTRG || retriangulationLastIter>1000) && retriangulationLastIter>10) {
Update_Triangulation = true;
Eps_Vol_Cumulative=0;
retriangulationLastIter=0;
ReTrg++;
} else retriangulationLastIter++;
timingDeltas->checkpoint("Update_Volumes");
+ ///Update boundary conditions
+ BoundaryConditions();
///Compute flow and and forces here
flow->GaussSeidel();
@@ -81,13 +83,16 @@
timingDeltas->checkpoint("Compute_Forces");
///Application of vicscous forces
-
+ scene->forces.sync();
if (viscousShear) ApplyViscousForces();
CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end();
for (CGT::Finite_vertices_iterator V_it = flow->T[flow->currentTes].Triangulation().finite_vertices_begin(); V_it != vertices_end; V_it++)
{
- scene->forces.addForce(V_it->info().id(), Vector3r ((V_it->info().forces)[0],V_it->info().forces[1],V_it->info().forces[2]));
+ 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])+flow->viscousShearForces[V_it->info().id()]);
}
///End Compute flow and forces
timingDeltas->checkpoint("Applying Forces");
@@ -118,7 +123,7 @@
string slice = "./Slices/Slice_"+flow->key+"_%d";
const char* keyslice = slice.c_str();
sprintf(slifile, keyslice, j);
- char *f = "slifile";
+ const char *f = "slifile";
flow->SliceField(f);
}
// if (save_vtk) {flow->save_vtk_file();}
@@ -158,12 +163,21 @@
flow->boundary ( flow->x_min_id ).value=Pressure_LEFT_Boundary;
flow->boundary ( flow->z_max_id ).value=Pressure_FRONT_Boundary;
flow->boundary ( flow->z_min_id ).value=Pressure_BACK_Boundary;
+
+ flow->boundary ( flow->y_min_id ).velocity = Vector3r::Zero();
+ flow->boundary ( flow->y_max_id ).velocity = topBoundaryVelocity;
+ flow->boundary ( flow->x_max_id ).velocity = Vector3r::Zero();
+ flow->boundary ( flow->x_min_id ).velocity = Vector3r::Zero();
+ flow->boundary ( flow->z_max_id ).velocity = Vector3r::Zero();
+ flow->boundary ( flow->z_min_id ).velocity = Vector3r::Zero();
}
void FlowEngine::imposePressure(Vector3r pos, Real p)
{
flow->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),p) );
+ //force immediate update of boundary conditions
+ Update_Triangulation=true;
}
void FlowEngine::clearImposedPressure() { flow->imposedP.clear();}
@@ -247,6 +261,7 @@
if (WaveAction) flow->ApplySinusoidalPressure(flow->T[flow->currentTes].Triangulation(), Sinus_Amplitude, Sinus_Average, 30);
}
Initialize_volumes();
+ if (viscousShear) flow->ComputeEdgesSurfaces();
}
void FlowEngine::AddBoundary ()
@@ -305,6 +320,7 @@
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 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;
@@ -318,7 +334,6 @@
flow->Corner_min = CGT::Point(flow->x_min, flow->y_min, flow->z_min);
flow->Corner_max = CGT::Point(flow->x_max, flow->y_max, flow->z_max);
-
if (Debug) {
cout << "Section area = " << flow->SectionArea << endl;
cout << "Vtotale = " << flow->Vtotale << endl;
@@ -330,6 +345,9 @@
cout << "z_min = " << flow->z_min << endl;
cout << "z_max = " << flow->z_max << endl;
cout << endl << "Adding Boundary------" << endl;}
+
+ //assign BCs types and values
+ BoundaryConditions();
double center[3];
@@ -368,11 +386,10 @@
Real x = b->state->pos[0];
Real y = b->state->pos[1];
Real z = b->state->pos[2];
-
flow->T[flow->currentTes].insert(x, y, z, rad, id);
-
}
}
+ flow->viscousShearForces.resize(flow->T[flow->currentTes].max_id+1);
}
void FlowEngine::Initialize_volumes ()
@@ -594,24 +611,40 @@
void FlowEngine::ApplyViscousForces()
{
- flow->ComputeEdgesSurfaces();
+// 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;
- vector <Vector3r> Edge_force;
- Edge_force.resize(flow->Edge_ids.size());
+ for (unsigned int k=0; k<flow->viscousShearForces.size(); k++) flow->viscousShearForces[k]=Vector3r::Zero();
+
+ const CGT::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 );
- Vector3r deltaV = (sph2->state->vel - sph1->state->vel) - (flow->Edge_normal[i].dot(sph2->state->vel - sph1->state->vel))*flow->Edge_normal[i];
+ 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(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;
+ 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;
- Edge_force.push_back(visc_f);
- scene->forces.addForce(flow->Edge_ids[i].first,Edge_force[i]);
- scene->forces.addForce(flow->Edge_ids[i].second,-Edge_force[i]);
- scene->forces.sync();
- if (Debug) cout<<"la force totale sur "<< flow->Edge_ids[i].first <<" est " << scene->forces.getForce(flow->Edge_ids[i].first) << "et la force totale sur "<< flow->Edge_ids[i].second <<" est " << scene->forces.getForce(flow->Edge_ids[i].second);
- if (Debug) cout<<"END: Application of vicscous forces"<<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)
+ flow->viscousShearForces[flow->Edge_ids[i].first]+=visc_f;
+ flow->viscousShearForces[flow->Edge_ids[i].second]-=visc_f;
}
}
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2011-05-03 10:33:29 +0000
+++ pkg/dem/FlowEngine.hpp 2011-05-03 20:03:21 +0000
@@ -56,12 +56,14 @@
void Average_real_cell_velocity();
void ApplyViscousForces();
Real getFlux(int cond);
- void saveVtk() {flow->save_vtk_file();}
- vector<Real> AvFlVelOnSph(int id_sph) {return flow->Average_Fluid_Velocity_On_Sphere(id_sph);}
+ void saveVtk() {flow->saveVtk();}
+ vector<Real> AvFlVelOnSph(unsigned int id_sph) {return flow->Average_Fluid_Velocity_On_Sphere(id_sph);}
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]);}
+ Vector3r fluidForce(unsigned int id_sph) {const CGT::Vecteur& f=flow->T[flow->currentTes].vertex(id_sph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
+ Vector3r fluidShearForce(unsigned int id_sph) {return (flow->viscousShearForces.size()>id_sph)?flow->viscousShearForces[id_sph]:Vector3r::Zero();}
+ void setBoundaryVel(Vector3r vel) {topBoundaryVelocity=vel; Update_Triangulation=true;}
virtual ~FlowEngine();
@@ -120,6 +122,7 @@
((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"))
+ ((Vector3r, topBoundaryVelocity, Vector3r::Zero(),, "velocity on top boundary, only change it using :yref:`FlowEngine::setBoundaryVel`"))
((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.)"))
@@ -145,6 +148,8 @@
.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.")
+ .def("fluidShearForce",&FlowEngine::fluidShearForce,(python::arg("Id_sph")),"Return the viscous shear force on sphere Id_sph.")
+ .def("setBoundaryVel",&FlowEngine::setBoundaryVel,(python::arg("vel")),"Change velocity on top boundary.")
)
DECLARE_LOGGER;
};
