yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #08326
[Branch ~yade-dev/yade/trunk] Rev 3028: - TWrapper: python binding for local deformations
------------------------------------------------------------
revno: 3028
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Tue 2012-02-14 17:58:21 +0100
message:
- TWrapper: python binding for local deformations
- Flow: fixes in the periodic BCs
- Flow: improvements in BCs handling and imposed fluxes, toward per-step update of the RHS + code cleaning
- Periodic flow: support for optimized solvers (not in trunk) by introducing cells indexes
- Periodic flow: optimization of periodicity by precomputed Hsize[k]*gradP products
- SpherePack: more code documentation in makeCloud
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphere.ipp
lib/triangulation/Network.hpp
lib/triangulation/PeriodicFlow.cpp
lib/triangulation/PeriodicFlow.hpp
lib/triangulation/Tesselation.ipp
lib/triangulation/TriaxialState.h
lib/triangulation/def_types.h
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.hpp
pkg/dem/SpherePack.cpp
pkg/dem/TesselationWrapper.cpp
pkg/dem/TesselationWrapper.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 2012-01-24 17:14:09 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2012-02-14 16:58:21 +0000
@@ -11,13 +11,15 @@
namespace CGT {
Vecteur PeriodicCellInfo::gradP;
Vecteur PeriodicCellInfo::hSize[3];
+Vecteur PeriodicCellInfo::deltaP;
}
+
//Forcing instanciation of the template to avoid linkage problems
typedef CGT::FlowBoundingSphere<FlowTesselation> FlowBoundingSphere;
FlowBoundingSphere ex;
#ifdef LINSOLV
-typedef CGT::FlowBoundingSphereLinSolv<FlowTesselation> FlowBoundingSphereLinSolv;
+typedef CGT::FlowBoundingSphereLinSolv<FlowBoundingSphere> FlowBoundingSphereLinSolv;
FlowBoundingSphereLinSolv exls;
#endif
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2012-02-01 14:42:03 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2012-02-14 16:58:21 +0000
@@ -29,7 +29,7 @@
DECLARE_TESSELATION_TYPES(Network<Tesselation>)
//painfull, but we need that for templates inheritance...
- using _N::T; using _N::x_min; using _N::x_max; using _N::y_min; using _N::y_max; using _N::z_min; using _N::z_max; using _N::Rmoy; using _N::SectionArea; using _N::Height; using _N::Vtotale; using _N::currentTes; using _N::DEBUG_OUT; using _N::nOfSpheres; using _N::x_min_id; using _N::x_max_id; using _N::y_min_id; using _N::y_max_id; using _N::z_min_id; using _N::z_max_id; using _N::boundsIds; using _N::Corner_min; using _N::Corner_max; using _N::Vsolid_tot; using _N::Vtotalissimo; using _N::Vporale; using _N::Ssolid_tot; using _N::V_porale_porosity; using _N::V_totale_porosity; using _N::boundaries; using _N::id_offset; using _N::vtk_infinite_vertices; using _N::vtk_infinite_cells; using _N::num_particles; using _N::fictious_vertex;
+ using _N::T; using _N::x_min; using _N::x_max; using _N::y_min; using _N::y_max; using _N::z_min; using _N::z_max; using _N::Rmoy; using _N::SectionArea; using _N::Height; using _N::Vtotale; using _N::currentTes; using _N::DEBUG_OUT; using _N::nOfSpheres; using _N::x_min_id; using _N::x_max_id; using _N::y_min_id; using _N::y_max_id; using _N::z_min_id; using _N::z_max_id; using _N::boundsIds; using _N::Corner_min; using _N::Corner_max; using _N::Vsolid_tot; using _N::Vtotalissimo; using _N::Vporale; using _N::Ssolid_tot; using _N::V_porale_porosity; using _N::V_totale_porosity; using _N::boundaries; using _N::id_offset; using _N::vtk_infinite_vertices; using _N::vtk_infinite_cells; using _N::num_particles; using _N::fictious_vertex; using _N::boundingCells;
//same for functions
using _N::Define_fictious_cells; using _N::AddBoundingPlanes; using _N::boundary;
@@ -87,6 +87,7 @@
void SpheresFileCreator ();
void DisplayStatistics();
void Initialize_pressures ( double P_zero );
+ void reApplyBoundaryConditions ();
/// Define forces using the same averaging volumes as for permeability
void ComputeTetrahedralForces();
/// Define forces spliting drag and buoyancy terms
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2012-01-25 10:02:04 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2012-02-14 16:58:21 +0000
@@ -981,16 +981,21 @@
for (int bound=0; bound<6;bound++) {
int& id = *boundsIds[bound];
+ boundingCells[bound].clear();
+ if (id<0) continue;
Boundary& bi = boundary(id);
if (!bi.flowCondition) {
Vector_Cell tmp_cells;
tmp_cells.resize(10000);
VCell_iterator cells_it = tmp_cells.begin();
VCell_iterator cells_end = Tri.incident_cells(T[currentTes].vertexHandles[id],cells_it);
- for (VCell_iterator it = tmp_cells.begin(); it != cells_end; it++)
- {(*it)->info().p() = bi.value;(*it)->info().Pcondition=true;}
+ for (VCell_iterator it = tmp_cells.begin(); it != cells_end; it++){
+ (*it)->info().p() = bi.value;(*it)->info().Pcondition=true;
+ boundingCells[bound].push_back(*it);
+ }
}
}
+ IPCells.clear();
for (unsigned int n=0; n<imposedP.size();n++) {
Cell_handle cell=Tri.locate(imposedP[n].first);
IPCells.push_back(cell);
@@ -1003,6 +1008,31 @@
cell->info().p()=imposedP[n].second;
cell->info().Pcondition=true;}
}
+
+template <class Tesselation>
+void FlowBoundingSphere<Tesselation>::reApplyBoundaryConditions()
+{
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ Finite_cells_iterator cell_end = Tri.finite_cells_end();
+
+ for (int bound=0; bound<6;bound++) {
+ int& id = *boundsIds[bound];
+ if (id<0) continue;
+ Boundary& bi = boundary(id);
+ if (!bi.flowCondition) {
+ for (VCell_iterator it = boundingCells[bound].begin(); it != boundingCells[bound].end(); it++){
+ (*it)->info().p() = bi.value;(*it)->info().Pcondition=true;
+ }
+ }
+ }
+ for (unsigned int n=0; n<imposedP.size();n++) {
+ IPCells[n]->info().p()=imposedP[n].second;
+ IPCells[n]->info().Pcondition=true;}
+}
+
+
+
+
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::GaussSeidel()
{
=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp 2012-01-20 17:31:56 +0000
+++ lib/triangulation/Network.hpp 2012-02-14 16:58:21 +0000
@@ -44,6 +44,7 @@
int nOfSpheres;
int x_min_id, x_max_id, y_min_id, y_max_id, z_min_id, z_max_id;
int* boundsIds [6];
+ vector<Cell_handle> boundingCells [6];
Point Corner_min;
Point Corner_max;
Real Vsolid_tot, Vtotalissimo, Vporale, Ssolid_tot, V_porale_porosity, V_totale_porosity;
=== modified file 'lib/triangulation/PeriodicFlow.cpp'
--- lib/triangulation/PeriodicFlow.cpp 2012-02-01 14:42:03 +0000
+++ lib/triangulation/PeriodicFlow.cpp 2012-02-14 16:58:21 +0000
@@ -83,6 +83,10 @@
}
}
//use cached values
+ //First define products that will be used below for all cells:
+ Real pDeltas [3];
+ for (unsigned int k=0; k<3;k++) pDeltas[k]=PeriodicCellInfo::hSize[k]*PeriodicCellInfo::gradP;
+ //Then compute the forces
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++){
if (cell->info().isGhost) continue;
for (int yy=0;yy<4;yy++) {
@@ -91,7 +95,7 @@
// cerr <<"p="<<cell->info().p()<<" shifted="<<cell->info().shiftedP()<<endl;
//the pressure translated to a ghost cell adjacent to the non-ghost vertex
//FIXME: the hsize[k]*gradP products could be computed only once for all cells and facets
- if (vhi.isGhost) unshiftedP -= (PeriodicCellInfo::hSize[0]*PeriodicCellInfo::gradP)*vhi.period[0] + (PeriodicCellInfo::hSize[1]*PeriodicCellInfo::gradP)*vhi.period[1] +(PeriodicCellInfo::hSize[2]*PeriodicCellInfo::gradP)*vhi.period[2];
+ if (vhi.isGhost) unshiftedP -= pDeltas[0]*vhi.period[0] + pDeltas[1]*vhi.period[1] +pDeltas[2]*vhi.period[2];
vhi.forces = vhi.forces + cell->info().unitForceVectors[yy]*unshiftedP;
// cerr <<"unshifted="<<unshiftedP<<endl;
}
@@ -348,10 +352,17 @@
{(*it)->info().setP(bi.value);(*it)->info().Pcondition=true;}
}
}
+ IPCells.clear();
for (unsigned int n=0; n<imposedP.size();n++) {
Cell_handle cell=Tri.locate(imposedP[n].first);
+
+ //check redundancy
+ for (unsigned int kk=0;kk<IPCells.size();kk++){
+ if (cell==IPCells[kk]) cerr<<"Two imposed pressures fall in the same cell."<<endl;
+ else if (cell->info().Pcondition) cerr<<"Imposed pressure fall in a boundary condition."<<endl;}
// cerr<<"cell found : "<<cell->vertex(0)->point()<<" "<<cell->vertex(1)->point()<<" "<<cell->vertex(2)->point()<<" "<<cell->vertex(3)->point()<<endl;
// assert(cell);
+ IPCells.push_back(cell);
cell->info().setP(imposedP[n].second);
cell->info().Pcondition=true;}
}
@@ -406,7 +417,7 @@
dp_moy = sum_dp/cell2;
j++;
if (j%100==0) cerr <<"j="<<j<<" p_moy="<<p_moy<<" dp="<< dp_moy<<" p_max="<<p_max<<" dp_max="<<dp_max<<endl;
- if (j>=40000) cerr<<"GS not converged after 40k iterations, break"<<endl;
+ if (j>=40000) cerr<<"\r GS not converged after 40k iterations, break";
} while ((dp_max/p_max) > tolerance && j<40000 /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
@@ -476,9 +487,9 @@
int num_cells = 0;
double facet_flow_rate = 0;
double volume_facet_translation = 0;
- Real tVel=0; Real tVol=0;
Finite_cells_iterator cell_end = Tri.finite_cells_end();
for ( Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++ ) {
+ if (cell->info().isGhost) continue;
cell->info().av_vel() =CGAL::NULL_VECTOR;
num_cells++;
for ( int i=0; i<4; i++ ) {
@@ -490,12 +501,10 @@
// Vecteur facetNormal = Surfk/area;
Vecteur branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
pos_av_facet = (Point) cell->info() + ( branch*Surfk ) *Surfk;
-// pos_av_facet=CGAL::ORIGIN + ((cell->vertex(facetVertices[i][0])->point() - CGAL::ORIGIN) + (cell->vertex(facetVertices[i][1])->point() - CGAL::ORIGIN) + (cell->vertex(facetVertices[i][2])->point() - CGAL::ORIGIN))*0.3333333333;
facet_flow_rate = (cell->info().k_norm())[i] * (cell->info().shiftedP() - cell->neighbor (i)->info().shiftedP());
- cell->info().av_vel() = cell->info().av_vel() + (facet_flow_rate) * ( pos_av_facet-CGAL::ORIGIN );
+ cell->info().av_vel() = cell->info().av_vel() + (facet_flow_rate) * ( pos_av_facet-CGAL::ORIGIN );
}}
- if (cell->info().volume()){ tVel+=cell->info().av_vel()[1]; tVol+=cell->info().volume();}
- cell->info().av_vel() = cell->info().av_vel() /cell->info().volume();
+ cell->info().av_vel() = cell->info().av_vel() /abs(cell->info().volume());
}
}
@@ -504,7 +513,7 @@
#endif //FLOW_ENGINE
-/*
+
#ifdef LINSOLV
-#include "FlowBoundingSphereLinSolv.cpp"
-#endif*/
+#include "PeriodicFlowLinSolv.ipp"
+#endif
=== modified file 'lib/triangulation/PeriodicFlow.hpp'
--- lib/triangulation/PeriodicFlow.hpp 2012-01-24 17:14:09 +0000
+++ lib/triangulation/PeriodicFlow.hpp 2012-02-14 16:58:21 +0000
@@ -14,7 +14,6 @@
// #include "Vue3D.h" //FIXME implicit dependencies will look for this class (out of tree) even ifndef XVIEW
// #endif
-// using namespace std;
namespace CGT{
@@ -41,9 +40,9 @@
};
}
-// #ifdef LINSOLV
-// #include "FlowBoundingSphereLinSolv.hpp"
-// #endif
+#ifdef LINSOLV
+#include "PeriodicFlowLinSolv.hpp"
+#endif
#endif //FLOW_ENGINE
=== modified file 'lib/triangulation/Tesselation.ipp'
--- lib/triangulation/Tesselation.ipp 2012-01-24 17:14:09 +0000
+++ lib/triangulation/Tesselation.ipp 2012-02-14 16:58:21 +0000
@@ -402,49 +402,50 @@
}
cell0=cell2++;
Cell_circulator cell1=cell2++;
- //std::cout << "edge : " << ed_it->first->vertex ( ed_it->second )->info().id() << "-" << ed_it->first->vertex ( ed_it->third )->info().id() << std::endl;
+// std::cout << "edge : " << ed_it->first->vertex ( ed_it->second )->info().id() << "-" << ed_it->first->vertex ( ed_it->third )->info().id() << std::endl;
bool isFictious1 = ( ed_it->first )->vertex ( ed_it->second )->info().isFictious;
bool isFictious2 = ( ed_it->first )->vertex ( ed_it->third )->info().isFictious;
+// cout<<"fictious "<<isFictious1<<" "<<isFictious2<<endl;
Real r;
- //cout << "cell0 : " << cell0->vertex(0)->info().id() << " "
- // << cell0->vertex(1)->info().id() << " "
- // << cell0->vertex(2)->info().id() << " "
- // << cell0->vertex(3)->info().id() << "(center : " << (Point) cell0->info() << ")" << endl;
+// cout << "cell0 : " << cell0->vertex(0)->info().id() << " "
+// << cell0->vertex(1)->info().id() << " "
+// << cell0->vertex(2)->info().id() << " "
+// << cell0->vertex(3)->info().id() << "(center : " << (Point) cell0->info() << ")" << endl;
while ( cell2!=cell0 )
{
if ( !Tri->is_infinite ( cell1 ) && !Tri->is_infinite ( cell2 ) )
{
- // cout << "cell1 : " << cell1->vertex(0)->info().id() << " "
- // << cell1->vertex(1)->info().id() << " "
- // << cell1->vertex(2)->info().id() << " "
- // << cell1->vertex(3)->info().id() << "(center : " << (Point) cell1->info() << ")" << endl;
- // cout << "cell2 : " << cell2->vertex(0)->info().id() << " "
- // << cell2->vertex(1)->info().id() << " "
- // << cell2->vertex(2)->info().id() << " "
- // << cell2->vertex(3)->info().id() << "(center : " << (Point) cell2->info() << ")" << endl;
+// cout << "cell1 : " << cell1->vertex(0)->info().id() << " "
+// << cell1->vertex(1)->info().id() << " "
+// << cell1->vertex(2)->info().id() << " "
+// << cell1->vertex(3)->info().id() << "(center : " << (Point) cell1->info() << ")" << endl;
+// cout << "cell2 : " << cell2->vertex(0)->info().id() << " "
+// << cell2->vertex(1)->info().id() << " "
+// << cell2->vertex(2)->info().id() << " "
+// << cell2->vertex(3)->info().id() << "(center : " << (Point) cell2->info() << ")" << endl;
- //std::cout << "assign tetra : (" << ed_it->first->vertex ( ed_it->second )->point() << "),(" << cell0->info() << "),(" << cell1->info() << "),(" <<cell2->info() << ")" << std::endl;
+// std::cout << "assign tetra : (" << ed_it->first->vertex ( ed_it->second )->point() << "),(" << cell0->info() << "),(" << cell1->info() << "),(" <<cell2->info() << ")" << std::endl;
if ( !isFictious1 )
{
r = std::abs ( ( Tetraedre ( ed_it->first->vertex ( ed_it->second )->point(), cell0->info(), cell1->info(), cell2->info() ) ).volume() );
- //std::cout << "assigned1=" << r << " on " << ed_it->first->vertex ( ed_it->second )->info().id() << std::endl;
+// std::cout << "assigned1=" << r << " on " << ed_it->first->vertex ( ed_it->second )->info().id() << std::endl;
( ed_it->first )->vertex ( ed_it->second )->info().v() += r;
TotalFiniteVoronoiVolume+=r;
}
- //std::cout << "assign tetra : (" << ed_it->first->vertex ( ed_it->third )->point() << "),(" << cell0->info() << "),(" << cell1->info() << "),(" <<cell2->info() << ")" << std::endl;
+// std::cout << "assign tetra : (" << ed_it->first->vertex ( ed_it->third )->point() << "),(" << cell0->info() << "),(" << cell1->info() << "),(" <<cell2->info() << ")" << std::endl;
if ( !isFictious2 )
{
r = std::abs ( ( Tetraedre ( ed_it->first->vertex ( ed_it->third )->point(), cell0->info(), cell1->info(), cell2->info() ) ).volume() );
- //std::cout << "assigned2=" << r << " on " << ed_it->first->vertex ( ed_it->third )->info().id() << std::endl;
+// std::cout << "assigned2=" << r << " on " << ed_it->first->vertex ( ed_it->third )->info().id() << std::endl;
ed_it->first->vertex ( ed_it->third )->info().v() +=r;
TotalFiniteVoronoiVolume+=r;
}
}
++cell1; ++cell2;
}
- //std::cout << "fin AssignPartialVolume" << std::endl;
+// std::cout << "fin AssignPartialVolume,total " << TotalFiniteVoronoiVolume<< std::endl;
}
template<class TT>
=== modified file 'lib/triangulation/TriaxialState.h'
--- lib/triangulation/TriaxialState.h 2012-01-20 17:31:56 +0000
+++ lib/triangulation/TriaxialState.h 2012-02-14 16:58:21 +0000
@@ -25,8 +25,6 @@
namespace CGT {
using namespace std;
-
-
class TriaxialState
{
public:
=== modified file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h 2012-01-24 17:14:09 +0000
+++ lib/triangulation/def_types.h 2012-02-14 16:58:21 +0000
@@ -71,6 +71,7 @@
inline Real& f (void) {return s;}
inline Real& v (void) {return vol;}
inline const unsigned int& id (void) const {return i;}
+ SimpleVertexInfo (void) {isFictious=false; s=0; i=0;}
};
@@ -163,20 +164,24 @@
static Vecteur gradP;
int period[3];
static Vecteur hSize[3];
+ static Vecteur deltaP;
int ghost;
Real* _pression;
bool isGhost;
PeriodicCellInfo (void)
{
_pression=&pression;
+ period[0]=period[1]=period[2]=0;
+ isGhost=false;
}
~PeriodicCellInfo (void) {}
PeriodicCellInfo& operator= (const Point &p) { Point::operator= (p); return *this; }
PeriodicCellInfo& operator= (const float &scalar) { s=scalar; return *this; }
- inline const double shiftedP (void) {return isGhost? (*_pression)+(hSize[0]*gradP)*period[0] + (hSize[1]*gradP)*period[1] +(hSize[2]*gradP)*period[2] :(*_pression) ;}
+ inline const double shiftedP (void) const {return isGhost? (*_pression)+pShift() :(*_pression) ;}
+ inline const double pShift (void) const {return deltaP[0]*period[0] + deltaP[1]*period[1] +deltaP[2]*period[2];}
// inline const double p (void) {return shiftedP();}
- inline double setP (Real p) {pression=p;}
+ inline void setP (const Real& p) {pression=p;}
};
class PeriodicVertexInfo : public FlowVertexInfo {
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2012-01-26 08:46:25 +0000
+++ pkg/dem/FlowEngine.cpp 2012-02-14 16:58:21 +0000
@@ -15,7 +15,6 @@
#include<yade/pkg/common/Sphere.hpp>
#include<yade/pkg/common/Wall.hpp>
#include<yade/pkg/common/Box.hpp>
-
#include <sys/stat.h>
#include <sys/types.h>
@@ -138,16 +137,16 @@
template<class Solver>
void FlowEngine::imposeFlux(Vector3r pos, Real flux,Solver& flow)
{
- RTriangulation& Tri = flow->T[flow->currentTes].Triangulation();
+ typename Solver::RTriangulation& Tri = flow.T[flow.currentTes].Triangulation();
double flux_base=0.f;
double perm_base=0.f;
- Cell_handle cell = Tri.locate(CGT::Point(pos[0],pos[1],pos[2]));
+ typename Solver::Cell_handle cell = Tri.locate(CGT::Point(pos[0],pos[1],pos[2]));
for (int ngb=0;ngb<4;ngb++) {
if (!cell->neighbor(ngb)->info().Pcondition) {
flux_base += cell->info().k_norm()[ngb]*(cell->neighbor(ngb)->info().p());
perm_base += cell->info().k_norm()[ngb];}}
- flow->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),(flux_base-flux)/perm_base));
+ flow.imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),(flux_base-flux)/perm_base));
//force immediate update of boundary conditions
Update_Triangulation=true;
}
@@ -587,17 +586,25 @@
void PeriodicFlowEngine:: action()
{
+ if (!isActivated) return;
CGT::PeriodicCellInfo::gradP = makeCgVect(gradP);
CGT::PeriodicCellInfo::hSize[0] = makeCgVect(scene->cell->hSize.col(0));
CGT::PeriodicCellInfo::hSize[1] = makeCgVect(scene->cell->hSize.col(1));
CGT::PeriodicCellInfo::hSize[2] = makeCgVect(scene->cell->hSize.col(2));
- if (!isActivated) return;
+ CGT::PeriodicCellInfo::deltaP=CGT::Vecteur(
+ CGT::PeriodicCellInfo::hSize[0]*CGT::PeriodicCellInfo::gradP,
+ CGT::PeriodicCellInfo::hSize[1]*CGT::PeriodicCellInfo::gradP,
+ CGT::PeriodicCellInfo::hSize[2]*CGT::PeriodicCellInfo::gradP);
+/*
+ CGT::PeriodicCellInfo::hSize[2]
+ (hSize[0]*gradP)*period[0] + (hSize[1]*gradP)*period[1] +(hSize[2]*gradP)*period[2]*/
+
// timingDeltas->start();
- if (first) Build_Triangulation(P_zero);
+ if (first) {Build_Triangulation(P_zero); Update_Triangulation = false;}
// timingDeltas->checkpoint("Triangulating");
UpdateVolumes ( );
- if (!first) {
+// if (!first) {
// 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>1000) && retriangulationLastIter>10) {
@@ -628,7 +635,7 @@
}
///End Compute flow and forces
// timingDeltas->checkpoint("Applying Forces");
- }
+// }
// if ( scene->iter % PermuteInterval == 0 )
// { Update_Triangulation = true; }
@@ -687,10 +694,6 @@
// cerr <<"number of vtx"<< solver->T[solver -> currentTes].Triangulation().number_of_vertices()<<endl;
}
}
- // Define the ghost cells
- Finite_cells_iterator cellend=solver->T[solver->currentTes].Triangulation().finite_cells_end();
- for (Finite_cells_iterator cell=solver -> T[solver -> currentTes].Triangulation().finite_cells_begin(); cell!=cellend; cell++)
- locateCell(cell);
solver -> viscousShearForces.resize(solver -> T[solver -> currentTes].max_id+1);
}
@@ -750,21 +753,46 @@
}
}
-void PeriodicFlowEngine::locateCell(Cell_handle baseCell)
+void PeriodicFlowEngine::locateCell(Cell_handle baseCell, unsigned int& index)
{
+ PeriFlowTesselation::Cell_Info& base_info = baseCell->info();
+ 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;
for (int k=0;k<4;k++) center+= 0.25*makeVector3r(baseCell->vertex(k)->point());
Vector3r wdCenter= scene->cell->wrapPt(center,period);
if (period[0]!=0 || period[1]!=0 || period[2]!=0) {
- baseCell->info().isGhost=true;
- Cell_handle ch= Tri.locate(CGT::Point(wdCenter[0],wdCenter[1],wdCenter[2]));
- baseCell->info().period[0]=period[0];
- baseCell->info().period[1]=period[1];
- baseCell->info().period[2]=period[2];
- baseCell->info()._pression=&(ch->info().p());
- } else baseCell->info().isGhost=false;
+ if (baseCell->info().index>0) {
+ cout<<"indexed cell is found ghost!"<<base_info.index <<endl;
+ base_info.isGhost=false;
+ return;
+ }
+// Cell_handle ch= Tri.locate(CGT::Point(wdCenter[0],wdCenter[1],wdCenter[2]),Tes.vertexHandles[baseCell->vertex(0)->info().id()]);//T[currentTes].vertexHandles[id]
+ Cell_handle ch= Tri.locate(CGT::Point(wdCenter[0],wdCenter[1],wdCenter[2]));//T[currentTes].vertexHandles[id]
+ base_info.period[0]=period[0];
+ base_info.period[1]=period[1];
+ base_info.period[2]=period[2];
+ //call recursively, since the returned cell could be also a ghost (especially if baseCell is a non-periodic type from the external contour
+ locateCell(ch,index);
+ if (ch==baseCell) cerr<<"WTF!!"<<endl;
+
+ base_info.isGhost=true;
+ //index is 1-based, if it is zero it is not initialized, we define it here
+// if (!ch->info().Pcondition && ch->info().index==0) {
+// ch->info().index=++index;
+// ch->info().isGhost=false;
+// }
+ //we make the ghost point to the real cell
+ base_info._pression=&(ch->info().p());
+ base_info.index=ch->info().index;
+ base_info.Pcondition=ch->info().Pcondition;
+ } else {
+ base_info.isGhost=false;
+ //index is 1-based, if it is zero it is not initialized, we define it here
+ if (!base_info.Pcondition && base_info.index==0) base_info.index=++index;
+ }
}
Vector3r PeriodicFlowEngine::MeanVelocity()
@@ -776,8 +804,8 @@
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] * cell->info().volume());
- volume+=cell->info().volume();
+ MeanVel[i]=MeanVel[i]+((cell->info().av_vel())[i] * abs(cell->info().volume()));
+ volume+=abs(cell->info().volume());
}
return (MeanVel/volume);
}
@@ -814,7 +842,7 @@
// totVol1+=newVol;
// break;
// case ( 0 ) :
- newVol = Volume_cell ( cell );
+// newVol = Volume_cell ( cell );
// totVol0+=newVol;
// break;
// default:
@@ -824,7 +852,7 @@
newVol = Volume_cell ( cell );
totVol+=newVol;
dVol=cell->info().volumeSign * (newVol - cell->info().volume());
- totDVol+=dVol;
+ totDVol+=dVol;
eps_vol_max = max(eps_vol_max, abs(dVol/newVol));
cell->info().dv() = dVol * invDeltaT;
cell->info().volume() = newVol;
@@ -885,7 +913,6 @@
solver->T[solver->currentTes].Compute();
// solver->Define_fictious_cells();
- solver->DisplayStatistics ();
solver->meanK_LIMIT = meanK_correction;
solver->meanK_STAT = meanK_opt;
@@ -899,9 +926,22 @@
for ( Finite_cells_iterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){cell->info().dv() = 0; cell->info().p() = 0;}
BoundaryConditions(solver);
+
+
solver->Initialize_pressures(P_zero);
- if (!first && useSolver==0) solver->Interpolate (solver->T[!solver->currentTes], solver->T[solver->currentTes]);
+ // Define the ghost cells and add indexes to the cells inside the period (the ones that will contain the pressure unknowns)
+ //This must be done after boundary conditions and before initialize pressure, else the indexes are not good (not accounting imposedP):
+ Finite_cells_iterator cellend=solver->T[solver->currentTes].Triangulation().finite_cells_end();
+ unsigned int index=0;
+ for (Finite_cells_iterator cell=solver -> T[solver -> currentTes].Triangulation().finite_cells_begin(); cell!=cellend; cell++)
+ locateCell(cell,index);
+
+ solver->DisplayStatistics ();
+
+ //FIXME: check interpolate() for the periodic case, at least use the mean pressure from previous step.
+// if (!first && useSolver==0) solver->Interpolate (solver->T[!solver->currentTes], solver->T[solver->currentTes]);
+
if (WaveAction) solver->ApplySinusoidalPressure(solver->T[solver->currentTes].Triangulation(), Sinus_Amplitude, Sinus_Average, 30);
Initialize_volumes();
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2012-01-26 08:46:25 +0000
+++ pkg/dem/FlowEngine.hpp 2012-02-14 16:58:21 +0000
@@ -16,7 +16,7 @@
class Flow;
class TesselationWrapper;
#ifdef LINSOLV
-#define _FlowSolver CGT::FlowBoundingSphereLinSolv<FlowTesselation>
+#define _FlowSolver CGT::FlowBoundingSphereLinSolv<CGT::FlowBoundingSphere<FlowTesselation> >
#else
#define _FlowSolver CGT::FlowBoundingSphere<FlowTesselation>
#endif
@@ -95,7 +95,7 @@
//Instanciation of templates for python binding
Vector3r _fluidShearForce(unsigned int id_sph) {return fluidShearForce(id_sph,solver);}
- void _imposeFlux(Vector3r pos, Real flux) {return imposeFlux(pos,flux,solver);}
+ void _imposeFlux(Vector3r pos, Real flux) {return imposeFlux(pos,flux,*solver);}
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);}
@@ -211,21 +211,17 @@
REGISTER_SERIALIZABLE(FlowEngine);
-// #ifdef LINSOLV
-// #define _PeriFlowSolver CGT::FlowBoundingSphereLinSolv<FlowTesselation>
-// #else
-// #define _PeriFlowSolver CGT::FlowBoundingSphere<PeriodicFlowTesselation>
-// #endif
+#ifdef LINSOLV
+#define _PeriFlowSolver CGT::PeriodicFlowLinSolv
+#else
+#define _PeriFlowSolver CGT::PeriodicFlow
+#endif
class PeriodicFlowEngine : public FlowEngine
{
- private:
-// shared_ptr<FlowSolver> flow;
-// int retriangulationLastIter;
-
public :
public :
- typedef CGT::PeriodicFlow FlowSolver;
+ typedef _PeriFlowSolver FlowSolver;
typedef PeriFlowTesselation Tesselation;
typedef FlowSolver::RTriangulation RTriangulation;
typedef FlowSolver::Finite_vertices_iterator Finite_vertices_iterator;
@@ -243,7 +239,7 @@
void UpdateVolumes ();
Real Volume_cell (Cell_handle cell);
void ApplyViscousForces();
- void locateCell(Cell_handle baseCell);
+ void locateCell(Cell_handle baseCell, unsigned int& index);
Vector3r MeanVelocity();
virtual ~PeriodicFlowEngine();
@@ -252,7 +248,7 @@
//Instanciation of templates for python binding
Vector3r _fluidShearForce(unsigned int id_sph) {return fluidShearForce(id_sph,solver);}
-// void _imposeFlux(Vector3r pos, Real flux) {return imposeFlux(pos,flux,solver);}
+ void _imposeFlux(Vector3r pos, Real flux) {return imposeFlux(pos,flux,*solver);}
unsigned int _imposePressure(Vector3r pos, Real p) {return imposePressure(pos,p,this->solver);}
// void _setImposedPressure(unsigned int cond, Real p) {setImposedPressure(cond,p,solver);}
// void _clearImposedPressure() {clearImposedPressure(solver);}
=== modified file 'pkg/dem/SpherePack.cpp'
--- pkg/dem/SpherePack.cpp 2012-02-01 18:06:57 +0000
+++ pkg/dem/SpherePack.cpp 2012-02-14 16:58:21 +0000
@@ -115,6 +115,7 @@
if (!area) throw invalid_argument("The box defined has null volume AND null surface. Define at least maxCorner of the box, or hSize if periodic.");
rMean=pow(area*(1-porosity)/(Mathr::PI*(1+rRelFuzz*rRelFuzz)*num),0.5);}
}
+ // transform sizes and cummulated fractions values in something convenient for the generation process
if(psdSizes.size()>0){
err=(mode>=0); mode=RDIST_PSD;
if(psdSizes.size()!=psdCumm.size()) throw invalid_argument(("SpherePack.makeCloud: psdSizes and psdCumm must have same dimensions ("+lexical_cast<string>(psdSizes.size())+"!="+lexical_cast<string>(psdCumm.size())).c_str());
@@ -137,8 +138,6 @@
appliedPsdScaling=1;
if(distributeMass) {
if (psdCumm2[psdSizes.size()-1]<num) appliedPsdScaling=pow(psdCumm2[psdSizes.size()-1]/num,1./3.);
- //Normalize psdCumm2 so it's between 0 and 1
-// for(size_t i=1; i<psdSizes.size(); i++) psdCumm2[i]/=psdCumm2[psdSizes.size()-1];
} else {
double totVol=0;
for(size_t i=1; i<psdSizes.size(); i++) totVol+= 4/3*Mathr::PI*(psdCumm[i]-psdCumm[i-1])*num*
@@ -158,13 +157,13 @@
Real r=0;
for(int i=0; (i<num) || (num<0); i++) {
Real norm, rand;
- //Determine radius of the next sphere we will attempt to place in space. If (num>0), generate radii the deterministic way, in decreasing order, else radii are stochastic since we don't know what the final number will be
+ //Determine radius of the next sphere that will be placed in space. If (num>0), generate radii the deterministic way, in decreasing order, else radii are stochastic since we don't know what the final number will be
if (num>0) rand = ((Real)num-(Real)i+0.5)/((Real)num+1.);
else rand = rnd();
int t;
switch(mode){
case RDIST_RMEAN:
- //FIXME : r is never defined, it will be zero at first iteration, but it will have values in the next ones. Some magic?
+ //FIXME : r is never defined, it will be zero at first iteration, but it will have values in the next ones. Some magic?
case RDIST_NUM:
if(distributeMass) r=pow3Interp(rand,rMean*(1-rRelFuzz),rMean*(1+rRelFuzz));
else r=rMean*(2*(rand-.5)*rRelFuzz+1); // uniform distribution in rMean*(1±rRelFuzz)
@@ -205,6 +204,7 @@
if (t==maxTry) {
if(num>0) {
if (mode!=RDIST_RMEAN) {
+ //if rMean is not imposed, then we call makeCloud recursively, scaling the PSD down until the target num is obtained
Real nextPoro = porosity+(1-porosity)/10.;
LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<num<<". Trying again with porosity "<<nextPoro<<". The size distribution is being scaled down");
pack.clear();
=== modified file 'pkg/dem/TesselationWrapper.cpp'
--- pkg/dem/TesselationWrapper.cpp 2012-01-23 14:43:54 +0000
+++ pkg/dem/TesselationWrapper.cpp 2012-02-14 16:58:21 +0000
@@ -158,7 +158,7 @@
// clock.top("Triangulation");
}
-double TesselationWrapper::Volume(unsigned int id) {return ((unsigned int) Tes->Max_id() > id) ? Tes->Volume(id) : -1;}
+double TesselationWrapper::Volume(unsigned int id) {return ((unsigned int) Tes->Max_id() >= id) ? Tes->Volume(id) : -1;}
bool TesselationWrapper::insert(double x, double y, double z, double rad, unsigned int id)
{
@@ -205,6 +205,7 @@
void TesselationWrapper::ComputeVolumes(void)
{
+ if (!bounded) AddBoundingPlanes();
ComputeTesselation();
Tes->ComputeVolumes();
}
=== modified file 'pkg/dem/TesselationWrapper.hpp'
--- pkg/dem/TesselationWrapper.hpp 2012-01-20 17:31:56 +0000
+++ pkg/dem/TesselationWrapper.hpp 2012-02-14 16:58:21 +0000
@@ -77,8 +77,13 @@
///Compute Voronoi vertices + volumes of all cells
///use ComputeTesselation to force update, e.g. after spheres positions have been updated
void ComputeVolumes (void);
+ void computeDeformations (void) {mma.analyser->ComputeParticlesDeformation();}
///Get volume of the sphere inserted with indentifier "id""
double Volume (unsigned int id);
+ double deformation (unsigned int id,unsigned int i,unsigned int j) {
+ if (!mma.analyser->ParticleDeformation.size()) {LOG_ERROR("Compute deformations first"); return 0;}
+ if (mma.analyser->ParticleDeformation.size()<id) {LOG_ERROR("id out of bounds"); return 0;}
+ return mma.analyser->ParticleDeformation[id](i,j);}
/// number of facets in the tesselation (finite branches of the triangulation)
unsigned int NumberOfFacets(bool initIters=false);
@@ -121,13 +126,15 @@
.def("triangulate",&TesselationWrapper::insertSceneSpheres,(python::arg("reset")=true),"triangulate spheres of the packing")
.def("setState",&TesselationWrapper::setState,(python::arg("state")=0),"Make the current state of the simulation the initial (0) or final (1) configuration for the definition of displacement increments, use only state=0 if you just want to get volmumes and porosity.")
.def("loadState",&TesselationWrapper::loadState,(python::arg("stateNumber")=0),"Load a file with positions to define state 0 or 1.")
- .def("saveState",&TesselationWrapper::saveState,(python::arg("outputFile")="state",python::arg("bz2")=true),"Save a file with positions, can be later reloaded in order to define state 0 or 1.")
+ .def("saveState",&TesselationWrapper::saveState,(python::arg("outputFile")="state",python::arg("state")=0,python::arg("bz2")=true),"Save a file with positions, can be later reloaded in order to define state 0 or 1.")
.def("volume",&TesselationWrapper::Volume,(python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
.def("defToVtk",&TesselationWrapper::defToVtk,(python::arg("outputFile")="def.vtk"),"Write local deformations in vtk format from states 0 and 1.")
- .def("defToVtkFromStates",&TesselationWrapper::defToVtkFromStates,(python::arg("input1")="state1",python::arg("input2")="state2",python::arg("outputFile")="def.vtk",python::arg("bz2")=false),"Write local deformations in vtk format from state files (since the file format is very special, consider using defToVtkFromPositions instead).")
- .def("defToVtkFromPositions",&TesselationWrapper::defToVtkFromPositions,(python::arg("input1")="pos1",python::arg("input2")="pos2",python::arg("outputFile")="def.vtk",python::arg("bz2")=false),"Write local deformations in vtk format from state files (since the file format is very special, consider using defToVtkFromPositions instead).")
+ .def("defToVtkFromStates",&TesselationWrapper::defToVtkFromStates,(python::arg("input1")="state1",python::arg("input2")="state2",python::arg("outputFile")="def.vtk",python::arg("bz2")=false),"Write local deformations in vtk format from state files (since the file format is very special, consider using defToVtkFromPositions if the input files were not generated by TesselationWrapper).")
+ .def("defToVtkFromPositions",&TesselationWrapper::defToVtkFromPositions,(python::arg("input1")="pos1",python::arg("input2")="pos2",python::arg("outputFile")="def.vtk",python::arg("bz2")=false),"Write local deformations in vtk format from positions files (one sphere per line, with x,y,z,rad separated by spaces).")
.def("computeVolumes",&TesselationWrapper::ComputeVolumes,"Compute volumes of all Voronoi's cells.")
.def("getVolPoroDef",&TesselationWrapper::getVolPoroDef,(python::arg("deformation")=false),"Return a table with per-sphere computed quantities. Include deformations on the increment defined by states 0 and 1 if deformation=True (make sure to define states 0 and 1 consistently).")
+ .def("ComputeDeformations",&TesselationWrapper::computeDeformations,"Compute per-particle deformation. Get it with :yref:`TesselationWrapper.deformation`(id,i,j).")
+ .def("deformation",&TesselationWrapper::deformation,(python::arg("id"),python::arg("i"),python::arg("j")),"Get particle deformation")
);
DECLARE_LOGGER;
};
