yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11431
[Branch ~yade-pkg/yade/git-trunk] Rev 3427: Merge github.com:yade/trunk into chaoUnsat
Merge authors:
Anton Gladky (gladky-anton)
Bruno Chareyre (bruno-chareyre)
Jan Stránský (honzik)
------------------------------------------------------------
revno: 3427 [merge]
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Fri 2014-06-06 17:27:15 +0200
message:
Merge github.com:yade/trunk into chaoUnsat
modified:
CMakeLists.txt
core/Body.hpp
core/ForceContainer.hpp
core/InteractionContainer.cpp
core/InteractionContainer.hpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphereLinSolv.hpp
lib/triangulation/Network.hpp
pkg/common/Collider.cpp
pkg/common/GravityEngines.cpp
pkg/common/SPHEngine.cpp
pkg/dem/CapillaryTriaxialTest.cpp
pkg/dem/CohFrictMat.hpp
pkg/dem/CohesiveFrictionalContactLaw.cpp
pkg/dem/CohesiveTriaxialTest.cpp
pkg/dem/NewtonIntegrator.cpp
pkg/dem/SimpleShear.cpp
pkg/dem/SpheresFactory.cpp
pkg/dem/TesselationWrapper.cpp
pkg/dem/TriaxialTest.cpp
pkg/dem/VTKRecorder.cpp
pkg/dem/ViscoelasticPM.hpp
pkg/pfv/FlowEngine.hpp
pkg/pfv/FlowEngine.ipp
pkg/pfv/PeriodicFlowEngine.cpp
py/export.py
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'CMakeLists.txt'
--- CMakeLists.txt 2014-05-27 11:24:16 +0000
+++ CMakeLists.txt 2014-05-30 11:58:01 +0000
@@ -434,6 +434,7 @@
ADD_LIBRARY(boot SHARED ${CMAKE_CURRENT_SOURCE_DIR}/core/main/pyboot.cpp)
SET_TARGET_PROPERTIES(boot PROPERTIES PREFIX "" LINK_FLAGS "-Wl,--as-needed" )
TARGET_LINK_LIBRARIES(yade ${Boost_LIBRARIES} ${PYTHON_LIBRARIES} ${LINKLIBS} -lrt)
+SET_TARGET_PROPERTIES(yade PROPERTIES LINK_FLAGS "-Wl,--as-needed" )
TARGET_LINK_LIBRARIES(boot yade)
IF(ENABLE_VTK)
=== modified file 'core/Body.hpp'
--- core/Body.hpp 2014-05-23 13:05:19 +0000
+++ core/Body.hpp 2014-06-05 13:19:44 +0000
@@ -20,6 +20,7 @@
#include<yade/lib/serialization/Serializable.hpp>
#include<yade/lib/multimethods/Indexable.hpp>
+
class Scene;
class Interaction;
@@ -29,6 +30,13 @@
typedef int id_t;
// internal structure to hold some interaction of a body; used by InteractionContainer;
typedef std::map<Body::id_t, shared_ptr<Interaction> > MapId2IntrT;
+ // groupMask type
+#ifdef BODY_GROUP_MASK_ARBITRARY_PRECISION
+ typedef boost::python::long_ groupMask_t;
+#else
+ typedef int groupMask_t;
+#endif
+
// bits for Body::flags
enum { FLAG_BOUNDED=1, FLAG_ASPHERICAL=2 }; /* add powers of 2 as needed */
@@ -62,13 +70,18 @@
* (otherwise, GLViewer would depend on Clump and therefore Clump would have to go to core...) */
virtual void userForcedDisplacementRedrawHook(){return;}
- boost::python::list py_intrs();
+ boost::python::list py_intrs();
Body::id_t getId() const {return id;};
unsigned int coordNumber(); // Number of neighboring particles
- int getGroupMask() {return groupMask; };
- bool maskOk(int mask){return (mask==0 || (groupMask&mask));}
+ Body::groupMask_t getGroupMask() const {return groupMask; };
+ bool maskOk(int mask) const {return (mask==0 || (bool)(groupMask & mask));}
+ bool maskCompatible(int mask) const { return (bool)(groupMask & mask); }
+#ifdef BODY_GROUP_MASK_ARBITRARY_PRECISION
+ bool maskOk(const boost::python::long_& mask) const {return (mask==0 || (bool)(groupMask & mask));}
+ bool maskCompatible(const boost::python::long_& mask) const { return (bool)(groupMask & mask); }
+#endif
// only BodyContainer can set the id of a body
friend class BodyContainer;
@@ -76,7 +89,7 @@
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Body,Serializable,"A particle, basic element of simulation; interacts with other bodies.",
((Body::id_t,id,Body::ID_NONE,Attr::readonly,"Unique id of this body."))
- ((int,groupMask,1,,"Bitmask for determining interactions."))
+ ((Body::groupMask_t,groupMask,1,,"Bitmask for determining interactions."))
((int,flags,FLAG_BOUNDED,Attr::readonly,"Bits of various body-related flags. *Do not access directly*. In c++, use isDynamic/setDynamic, isBounded/setBounded, isAspherical/setAspherical. In python, use :yref:`Body.dynamic`, :yref:`Body.bounded`, :yref:`Body.aspherical`."))
((shared_ptr<Material>,material,,,":yref:`Material` instance associated with this body."))
@@ -122,3 +135,16 @@
);
};
REGISTER_SERIALIZABLE(Body);
+
+
+#ifdef BODY_GROUP_MASK_ARBITRARY_PRECISION
+namespace boost { namespace serialization {
+template<class Archive>
+void serialize(Archive & ar, boost::python::long_ & l, const unsigned int version){
+ namespace bp = boost::python;
+ std::string value = bp::extract<std::string>(bp::str(l));
+ ar & BOOST_SERIALIZATION_NVP(value);
+ l = bp::long_(value);
+}
+}} // namespace boost::serialization
+#endif
=== modified file 'core/ForceContainer.hpp'
--- core/ForceContainer.hpp 2014-05-21 12:37:11 +0000
+++ core/ForceContainer.hpp 2014-06-04 08:44:44 +0000
@@ -65,12 +65,13 @@
assert(nThreads>omp_get_thread_num());
const Body::id_t idMaxTmp = max(id, _maxId[threadN]);
_maxId[threadN] = 0;
- if (threadN<0) resizePerm(min((size_t)1.5*(idMaxTmp+100),(size_t)(idMaxTmp+2000)));
- else if (sizeOfThreads[threadN]<=(size_t)idMaxTmp) resize(min((size_t)1.5*(idMaxTmp+100),(size_t)(idMaxTmp+2000)),threadN);
+ if (threadN<0) {
+ resizePerm(min((size_t)1.5*(idMaxTmp+100),(size_t)(idMaxTmp+2000)));
+ } else if (sizeOfThreads[threadN]<=(size_t)idMaxTmp) {
+ resize(min((size_t)1.5*(idMaxTmp+100),(size_t)(idMaxTmp+2000)),threadN);
+ }
}
-
inline void ensureSynced(){ if(!synced) throw runtime_error("ForceContainer not thread-synchronized; call sync() first!"); }
-
// dummy function to avoid template resolution failure
friend class boost::serialization::access; template<class ArchiveT> void serialize(ArchiveT & ar, unsigned int version){}
@@ -133,10 +134,17 @@
* Locks globalMutex, since one thread modifies common data (_force&_torque).
* Must be called before get* methods are used. Exception is thrown otherwise, since data are not consistent. */
inline void sync(){
+ for(int i=0; i<nThreads; i++){
+ if (_maxId[i] > 0) { synced = false;}
+ }
if(synced) return;
boost::mutex::scoped_lock lock(globalMutex);
if(synced) return; // if synced meanwhile
-
+
+ for(int i=0; i<nThreads; i++){
+ if (_maxId[i] > 0) { ensureSize(_maxId[i],i);}
+ }
+
syncSizesOfContainers();
for(long id=0; id<(long)size; id++){
@@ -283,10 +291,16 @@
lastReset=iter;
}
- void sync(){
- if (permForceUsed) for(long id=0; id<(long)size; id++)
- {_force[id]+=_permForce[id]; _torque[id]+=_permTorque[id];}
- return;}
+ void sync() {
+ if (_maxId>0) {ensureSize(_maxId); _maxId=0;}
+ if (permForceUsed) {
+ for(long id=0; id<(long)size; id++) {
+ _force[id]+=_permForce[id];
+ _torque[id]+=_permTorque[id];
+ }
+ }
+ return;
+ }
unsigned long syncCount;
// interaction in which the container was last reset; used by NewtonIntegrator to detect whether ForceResetter was not forgotten
long lastReset;
=== modified file 'core/InteractionContainer.cpp'
--- core/InteractionContainer.cpp 2014-05-23 13:05:19 +0000
+++ core/InteractionContainer.cpp 2014-06-04 08:51:38 +0000
@@ -59,7 +59,7 @@
const shared_ptr<Body>& b1((*bodies)[id1]);
const shared_ptr<Body>& b2((*bodies)[id2]);
-
+ LOG_DEBUG("InteractionContainer erase intrs id1=" << id1 << " id2=" << id2);
int linIx=-1;
if(!b1) linIx=linPos;
else {
@@ -67,6 +67,7 @@
if(I==b1->intrs.end()) linIx=linPos;
else {
linIx=I->second->linIx;
+ LOG_DEBUG("InteractionContainer linIx=" << linIx << " linPos=" << linPos);
assert(linIx==linPos);
//erase from body, we also erase from linIntrs below
b1->intrs.erase(I);
=== modified file 'core/InteractionContainer.hpp'
--- core/InteractionContainer.hpp 2014-02-22 04:56:33 +0000
+++ core/InteractionContainer.hpp 2014-06-04 08:51:38 +0000
@@ -108,29 +108,12 @@
*/
template<class T> size_t conditionalyEraseNonReal(const T& t, Scene* rb){
// beware iterators here, since erase is invalidating them. We need to iterate carefully, and keep in mind that erasing one interaction is moving the last one to the current position.
- // For the parallel flavor we build the list to be erased in parallel, then it is erased sequentially. Still significant speedup since checking bounds is the most expensive part.
- #ifndef YADE_OPENMP
size_t initSize=currSize;
for (size_t linPos=0; linPos<currSize;){
const shared_ptr<Interaction>& i=linIntrs[linPos];
if(!i->isReal() && t.shouldBeErased(i->getId1(),i->getId2(),rb)) erase(i->getId1(),i->getId2(),linPos);
else linPos++;}
return initSize-currSize;
- #else
- unsigned nThreads= omp_get_max_threads();
- assert(nThreads>0);
- std::vector<std::vector<Vector3i > >toErase;
- toErase.resize(nThreads,std::vector<Vector3i >());
- for (unsigned kk=0; kk<nThreads; kk++) toErase[kk].reserve(1000);//A smarter value than 1000?
- size_t initSize=currSize;
- #pragma omp parallel for schedule(guided,100) num_threads(nThreads)
- for (size_t linPos=0; linPos<currSize;linPos++){
- const shared_ptr<Interaction>& i=linIntrs[linPos];
- if(!i->isReal() && t.shouldBeErased(i->getId1(),i->getId2(),rb)) toErase[omp_get_thread_num()].push_back(Vector3i(i->getId1(),i->getId2(),linPos)) ;
- }
- for (unsigned int kk=0; kk<nThreads; kk++) for (size_t ii(0), jj(toErase[kk].size()); ii<jj;ii++) erase(toErase[kk][ii][0],toErase[kk][ii][1],toErase[kk][ii][2]);
- return initSize-currSize;
- #endif
}
// we must call Scene's ctor (and from Scene::postLoad), since we depend on the existing BodyContainer at that point.
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2014-05-26 23:23:29 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2014-06-04 17:19:50 +0000
@@ -31,7 +31,7 @@
//painfull, but we need that for templates inheritance...
using _N::T; using _N::xMin; using _N::xMax; using _N::yMin; using _N::yMax; using _N::zMin; using _N::zMax; using _N::Rmoy; using _N::sectionArea; using _N::Height; using _N::vTotal; using _N::currentTes; using _N::debugOut; using _N::nOfSpheres; using _N::xMinId; using _N::xMaxId; using _N::yMinId; using _N::yMaxId; using _N::zMinId; using _N::zMaxId; using _N::boundsIds; using _N::cornerMin; using _N::cornerMax; using _N::VSolidTot; using _N::Vtotalissimo; using _N::vPoral; using _N::sSolidTot; using _N::vPoralPorosity; using _N::vTotalPorosity; using _N::boundaries; using _N::idOffset; using _N::vtkInfiniteVertices; using _N::vtkInfiniteCells; using _N::num_particles; using _N::boundingCells; using _N::facetVertices; using _N::facetNFictious;
//same for functions
- using _N::defineFictiousCells; using _N::addBoundingPlanes; using _N::boundary;
+ using _N::defineFictiousCells; using _N::addBoundingPlanes; using _N::boundary; using _N::tesselation;
virtual ~FlowBoundingSphere();
FlowBoundingSphere();
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.hpp'
--- lib/triangulation/FlowBoundingSphereLinSolv.hpp 2014-04-10 05:58:01 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.hpp 2014-06-04 17:19:50 +0000
@@ -54,6 +54,7 @@
using FlowType::pressureChanged;
using FlowType::computedOnce;
using FlowType::resetNetwork;
+ using FlowType::tesselation;
//! TAUCS DECs
vector<FiniteCellsIterator> orderedCells;
=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp 2014-05-09 14:24:08 +0000
+++ lib/triangulation/Network.hpp 2014-06-06 15:27:15 +0000
@@ -46,6 +46,8 @@
Tesselation T [2];
bool currentTes;
+ Tesselation& tesselation() {return T[currentTes];};
+
double xMin, xMax, yMin, yMax, zMin, zMax, Rmoy, sectionArea, Height, vTotal;
bool debugOut;
int nOfSpheres;
=== modified file 'pkg/common/Collider.cpp'
--- pkg/common/Collider.cpp 2014-05-23 13:05:19 +0000
+++ pkg/common/Collider.cpp 2014-06-05 13:19:44 +0000
@@ -20,9 +20,9 @@
// do not collide clumps, since they are just containers, never interact
!b1->isClump() && !b2->isClump() &&
// masks must have at least 1 bit in common
- (b1->groupMask & b2->groupMask)!=0 &&
+ b1->maskCompatible(b2->groupMask) &&
// avoid contact between particles having the same mask compatible with the avoidSelfInteractionMask.
- !( (b1->groupMask == b2->groupMask) && ( avoidSelfInteractionMask & b1->groupMask ))
+ !( (b1->groupMask == b2->groupMask) && b1->maskCompatible(avoidSelfInteractionMask) )
;
}
=== modified file 'pkg/common/GravityEngines.cpp'
--- pkg/common/GravityEngines.cpp 2014-05-23 13:05:19 +0000
+++ pkg/common/GravityEngines.cpp 2014-06-05 13:19:44 +0000
@@ -22,7 +22,7 @@
YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
// skip clumps, only apply forces on their constituents
if(b->isClump()) continue;
- if(mask!=0 && (b->groupMask & mask)==0) continue;
+ if(mask!=0 && !b->maskCompatible(mask)) continue;
scene->forces.addForce(b->getId(),gravity*b->state->mass);
// work done by gravity is "negative", since the energy appears in the system from outside
if(trackEnergy) scene->energy->add(-gravity.dot(b->state->vel)*b->state->mass*dt,"gravWork",fieldWorkIx,/*non-incremental*/false);
@@ -33,7 +33,7 @@
const Vector3r& centralPos=Body::byId(centralBody)->state->pos;
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
if(b->isClump() || b->getId()==centralBody) continue; // skip clumps and central body
- if(mask!=0 && (b->groupMask & mask)==0) continue;
+ if(mask!=0 && !b->maskCompatible(mask)) continue;
Real F=accel*b->state->mass;
Vector3r toCenter=centralPos-b->state->pos; toCenter.normalize();
scene->forces.addForce(b->getId(),F*toCenter);
@@ -44,7 +44,7 @@
void AxialGravityEngine::action(){
FOREACH(const shared_ptr<Body>&b, *scene->bodies){
if(!b || b->isClump()) continue;
- if(mask!=0 && (b->groupMask & mask)==0) continue;
+ if(mask!=0 && !b->maskCompatible(mask)) continue;
/* http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html */
const Vector3r& x0=b->state->pos;
const Vector3r& x1=axisPoint;
=== modified file 'pkg/common/SPHEngine.cpp'
--- pkg/common/SPHEngine.cpp 2014-04-14 06:43:15 +0000
+++ pkg/common/SPHEngine.cpp 2014-05-30 13:24:10 +0000
@@ -186,8 +186,9 @@
Real smoothkernelLucy(const double & rr, const double & hh) {
if (rr<=hh and rr!=0 and hh!=0) {
- const Real h = 1; const Real r = rr/hh;
+ //const Real h = 1;
//return 5/(9*M_PI*pow(h,2))*(1+3*r/h)*pow((1-r/h),3);
+ const Real r = rr/hh;
return 5/(9*M_PI)*(1+3*r)*pow((1-r),3);
} else {
return 0;
@@ -196,8 +197,9 @@
Real smoothkernelLucyGrad(const double & rr, const double & hh) {
if (rr<=hh and rr!=0 and hh!=0) {
- const Real h = 1; const Real r = rr/hh;
+ //const Real h = 1;
//return -5/(9*M_PI*pow(h,2))*(-12*r/(h*h))*pow((1-r/h),2);
+ const Real r = rr/hh;
return -5/(9*M_PI)*(-12*r)*pow((1-r),2);
} else {
return 0;
@@ -206,8 +208,9 @@
Real smoothkernelLucyLapl(const double & rr, const double & hh) {
if (rr<=hh and rr!=0 and hh!=0) {
- const Real h = 1; const Real r = rr/hh;
+ //const Real h = 1;
//return 5/(9*M_PI*pow(h,2))*(-12/(h*h))*(1-r/h)*(1-3*r/h);
+ const Real r = rr/hh;
return 5/(9*M_PI)*(-12)*(1-r)*(1-3*r);
} else {
return 0;
@@ -217,7 +220,8 @@
Real smoothkernelMonaghan(const double & rr, const double & hh) {
Real ret = 0.0;
if (hh!=0) {
- const Real h = 1; const Real r = rr/hh;
+ //const Real h = 1;
+ const Real r = rr/hh;
if (rr/hh<0.5) {
//ret = 40/(7*M_PI*h*h)*(1 - 6*pow((r/h),2) + 6*pow((r/h),3));
ret = 40/(7*M_PI)*(1 - 6*pow((r),2) + 6*pow((r),3));
@@ -247,7 +251,8 @@
Real smoothkernelMonaghanLapl(const double & rr, const double & hh) {
Real ret = 0.0;
if (hh!=0) {
- const Real h = 1; const Real r = rr/hh;
+ //const Real h = 1;
+ const Real r = rr/hh;
if (rr/hh<0.5) {
//ret = 40/(7*M_PI*h*h)*( - 12/(h*h))*(1 - 3 * r/(h*h*h*h*h));
ret = 40/(7*M_PI)*( - 12)*(1 - 3 * r);
@@ -343,7 +348,7 @@
const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel;
const Real normalVelocity = geom.normal.dot(relativeVelocity);
- const Vector3r shearVelocity = relativeVelocity-normalVelocity*geom.normal;
+ //const Vector3r shearVelocity = relativeVelocity-normalVelocity*geom.normal;
// Copy-paste
//////////////////////////////////////////////////////////////////
=== modified file 'pkg/dem/CapillaryTriaxialTest.cpp'
--- pkg/dem/CapillaryTriaxialTest.cpp 2014-05-23 13:20:43 +0000
+++ pkg/dem/CapillaryTriaxialTest.cpp 2014-06-05 13:19:44 +0000
@@ -209,7 +209,7 @@
void CapillaryTriaxialTest::createSphere(shared_ptr<Body>& body, Vector3r position, Real radius, bool big, bool dynamic )
{
- body = shared_ptr<Body>(new Body); body->groupMask=2;
+ body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
shared_ptr<FrictMat> physics(new FrictMat);
shared_ptr<Aabb> aabb(new Aabb);
@@ -254,7 +254,7 @@
void CapillaryTriaxialTest::createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents, bool wire)
{
- body = shared_ptr<Body>(new Body); body->groupMask=2;
+ body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
shared_ptr<FrictMat> physics(new FrictMat);
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Box> iBox(new Box);
=== modified file 'pkg/dem/CohFrictMat.hpp'
--- pkg/dem/CohFrictMat.hpp 2014-05-26 16:49:23 +0000
+++ pkg/dem/CohFrictMat.hpp 2014-06-02 09:07:27 +0000
@@ -25,8 +25,8 @@
((Real,alphaKtw,2.0,,"Dimensionless twist stiffness."))
((Real,etaRoll,-1.,,"Dimensionless rolling (aka 'bending') strength. If negative, rolling moment will be elastic."))
((Real,etaTwist,-1.,,"Dimensionless twisting strength. If negative, twist moment will be elastic."))
- ((Real,normalCohesion,0,,"Tensile strength, homogeneous to a pressure"))
- ((Real,shearCohesion,0,,"Shear strength, homogeneous to a pressure"))
+ ((Real,normalCohesion,-1,,"Tensile strength, homogeneous to a pressure. If negative the normal force is purely elastic."))
+ ((Real,shearCohesion,-1,,"Shear strength, homogeneous to a pressure. If negative the shear force is purely elastic."))
((bool,momentRotationLaw,false,,"Use bending/twisting moment at contact. The contact will have moments only if both bodies have this flag true. See :yref:`CohFrictPhys::cohesionDisablesFriction` for details."))
,
createIndex();
=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp 2014-05-26 16:49:23 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp 2014-06-02 09:07:27 +0000
@@ -193,8 +193,8 @@
}
/// Plasticity ///
// limit rolling moment to the plastic value, if required
- Real RollMax = phys->maxRollPl*phys->normalForce.norm();
- if (RollMax>=0.){ // do we want to apply plasticity?
+ if (phys->maxRollPl>=0.){ // do we want to apply plasticity?
+ Real RollMax = phys->maxRollPl*phys->normalForce.norm();
if (!useIncrementalForm) LOG_WARN("If :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::useIncrementalForm` is false, then plasticity will not be applied correctly (the total formulation would not reproduce irreversibility).");
Real scalarRoll = phys->moment_bending.norm();
if (scalarRoll>RollMax){ // fix maximum rolling moment
@@ -206,8 +206,8 @@
}
}
// limit twisting moment to the plastic value, if required
- Real TwistMax = phys->maxTwistPl*phys->normalForce.norm();
- if (TwistMax>=0.){ // do we want to apply plasticity?
+ if (phys->maxTwistPl>=0.){ // do we want to apply plasticity?
+ Real TwistMax = phys->maxTwistPl*phys->normalForce.norm();
if (!useIncrementalForm) LOG_WARN("If :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::useIncrementalForm` is false, then plasticity will not be applied correctly (the total formulation would not reproduce irreversibility).");
Real scalarTwist= phys->moment_twist.norm();
if (scalarTwist>TwistMax){ // fix maximum rolling moment
=== modified file 'pkg/dem/CohesiveTriaxialTest.cpp'
--- pkg/dem/CohesiveTriaxialTest.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/CohesiveTriaxialTest.cpp 2014-06-05 13:19:44 +0000
@@ -200,7 +200,7 @@
void CohesiveTriaxialTest::createSphere(shared_ptr<Body>& body, Vector3r position, Real radius, bool dynamic )
{
- body = shared_ptr<Body>(new Body); body->groupMask=2;
+ body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
shared_ptr<CohFrictMat> physics(new CohFrictMat);
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Sphere> iSphere(new Sphere);
@@ -243,7 +243,7 @@
void CohesiveTriaxialTest::createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents, bool wire)
{
- body = shared_ptr<Body>(new Body); body->groupMask=2;
+ body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
shared_ptr<CohFrictMat> physics(new CohFrictMat);
shared_ptr<Aabb> aabb(new Aabb);
=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2014-06-05 13:19:44 +0000
@@ -229,7 +229,7 @@
//Reflect mean-field (periodic cell) acceleration in the velocity
if(scene->isPeriodic && homoDeform) {Vector3r dVel=dVelGrad*state->pos; state->vel+=dVel;}
- if ( (mask<=0) or ((mask>0) and ((Body::byId(id)->groupMask & mask)!=0)) ) {
+ if ( (mask<=0) or ((mask>0) and (Body::byId(id)->maskCompatible(mask))) ) {
state->pos+=state->vel*dt;
}
}
@@ -237,13 +237,13 @@
void NewtonIntegrator::leapfrogSphericalRotate(State* state, const Body::id_t& id, const Real& dt )
{
Real angle2=state->angVel.squaredNorm();
- if (angle2!=0 and ( (mask<=0) or ((mask>0) and ((Body::byId(id)->groupMask & mask)!=0)) )) {//If we have an angular velocity, we make a rotation
+ if (angle2!=0 and ( (mask<=0) or ((mask>0) and (Body::byId(id)->maskCompatible(mask))) )) {//If we have an angular velocity, we make a rotation
Real angle=sqrt(angle2);
Quaternionr q(AngleAxisr(angle*dt,state->angVel/angle));
state->ori = q*state->ori;
}
if(scene->forces.getMoveRotUsed() && scene->forces.getRot(id)!=Vector3r::Zero()
- and ( (mask<=0) or ((mask>0) and ((Body::byId(id)->groupMask & mask)!=0)) )) {
+ and ( (mask<=0) or ((mask>0) and (Body::byId(id)->maskCompatible(mask))) )) {
Vector3r r(scene->forces.getRot(id));
Real norm=r.norm(); r/=norm;
Quaternionr q(AngleAxisr(norm,r));
=== modified file 'pkg/dem/SimpleShear.cpp'
--- pkg/dem/SimpleShear.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/SimpleShear.cpp 2014-06-05 13:19:44 +0000
@@ -116,7 +116,7 @@
void SimpleShear::createSphere(shared_ptr<Body>& body, Vector3r position, Real radius)
{
- body = shared_ptr<Body>(new Body); body->groupMask=1;
+ body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(1);
shared_ptr<NormalInelasticMat> mat(new NormalInelasticMat);
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Sphere> iSphere(new Sphere);
@@ -147,7 +147,7 @@
void SimpleShear::createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents)
{
- body = shared_ptr<Body>(new Body); body->groupMask=1;
+ body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(1);
shared_ptr<NormalInelasticMat> mat(new NormalInelasticMat);
shared_ptr<Aabb> aabb(new Aabb);
// shared_ptr<BoxModel> gBox(new BoxModel);
=== modified file 'pkg/dem/SpheresFactory.cpp'
--- pkg/dem/SpheresFactory.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/SpheresFactory.cpp 2014-06-05 13:19:44 +0000
@@ -182,7 +182,7 @@
b->shape=sphere;
b->state=state;
b->material=material;
- if (mask>0) {b->groupMask=mask;}
+ if (mask>0) {b->groupMask=Body::groupMask_t(mask);}
// insert particle in the simulation
scene->bodies->insert(b);
ids.push_back(b->getId());
=== modified file 'pkg/dem/TesselationWrapper.cpp'
--- pkg/dem/TesselationWrapper.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/TesselationWrapper.cpp 2014-06-05 13:19:44 +0000
@@ -306,7 +306,7 @@
void createSphere(shared_ptr<Body>& body, Vector3r position, Real radius, bool big, bool dynamic )
{
- body = shared_ptr<Body>(new Body); body->groupMask=2;
+ body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
shared_ptr<Sphere> iSphere(new Sphere);
body->state->blockedDOFs=State::DOF_NONE;
body->state->pos=position;
=== modified file 'pkg/dem/TriaxialTest.cpp'
--- pkg/dem/TriaxialTest.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/TriaxialTest.cpp 2014-06-05 13:19:44 +0000
@@ -187,7 +187,7 @@
void TriaxialTest::createSphere(shared_ptr<Body>& body, Vector3r position, Real radius, bool big, bool dynamic )
{
- body = shared_ptr<Body>(new Body); body->groupMask=2;
+ body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Sphere> iSphere(new Sphere);
body->state->blockedDOFs=State::DOF_NONE;
@@ -213,7 +213,7 @@
void TriaxialTest::createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents, bool wire)
{
- body = shared_ptr<Body>(new Body); body->groupMask=2;
+ body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
body->state->blockedDOFs=State::DOF_ALL;
shared_ptr<Aabb> aabb(new Aabb);
aabb->color = Vector3r(1,0,0);
=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/VTKRecorder.cpp 2014-06-05 13:19:44 +0000
@@ -37,6 +37,12 @@
YADE_PLUGIN((VTKRecorder));
CREATE_LOGGER(VTKRecorder);
+#ifdef BODY_GROUP_MASK_ARBITRARY_PRECISION
+#define GET_MASK(b) boost::python::extract<int>(b->groupMask)
+#else
+#define GET_MASK(b) b->groupMask
+#endif
+
void VTKRecorder::action(){
vector<bool> recActive(REC_SENTINEL,false);
FOREACH(string& rec, recorders){
@@ -414,7 +420,7 @@
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
if (!b) continue;
- if(mask!=0 && (b->groupMask & mask)==0) continue;
+ if(mask!=0 && !b->maskCompatible(mask)) continue;
if (recActive[REC_SPHERES]){
const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
if (sphere){
@@ -425,7 +431,7 @@
spheresCells->InsertNextCell(1,pid);
radii->InsertNextValue(sphere->radius);
if (recActive[REC_ID]) spheresId->InsertNextValue(b->getId());
- if (recActive[REC_MASK]) spheresMask->InsertNextValue(b->groupMask);
+ if (recActive[REC_MASK]) spheresMask->InsertNextValue(GET_MASK(b));
if (recActive[REC_MASS]) spheresMass->InsertNextValue(b->state->mass);
if (recActive[REC_CLUMPID]) clumpId->InsertNextValue(b->clumpId);
if (recActive[REC_COLORS]){
@@ -508,7 +514,7 @@
facetsForceLen->InsertNextValue(stress.norm());
}
if (recActive[REC_MATERIALID]) facetsMaterialId->InsertNextValue(b->material->id);
- if (recActive[REC_MASK]) facetsMask->InsertNextValue(b->groupMask);
+ if (recActive[REC_MASK]) facetsMask->InsertNextValue(GET_MASK(b));
continue;
}
}
@@ -560,7 +566,7 @@
boxesForceLen->InsertNextValue(stress.norm());
}
if (recActive[REC_MATERIALID]) boxesMaterialId->InsertNextValue(b->material->id);
- if (recActive[REC_MASK]) boxesMask->InsertNextValue(b->groupMask);
+ if (recActive[REC_MASK]) boxesMask->InsertNextValue(GET_MASK(b));
}
continue;
}
@@ -885,3 +891,4 @@
};
#endif /* YADE_VTK */
+#undef GET_MASK
=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp 2014-04-10 14:04:50 +0000
+++ pkg/dem/ViscoelasticPM.hpp 2014-06-03 08:53:31 +0000
@@ -23,7 +23,7 @@
class ViscElMat : public FrictMat {
public:
virtual ~ViscElMat();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(ViscElMat,FrictMat,"Material for simple viscoelastic model of contact.\n\n.. note::\n\t ``Shop::getViscoelasticFromSpheresInteraction`` (and :yref:`yade.utils.getViscoelasticFromSpheresInteraction` in python) compute :yref:`kn<ViscElMat::kn>`, :yref:`cn<ViscElMat::cn>`, :yref:`ks<ViscElMat::ks>`, :yref:`cs<ViscElMat::cs>` from analytical solution of a pair spheres interaction problem.",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(ViscElMat,FrictMat,"Material for simple viscoelastic model of contact from analytical solution of a pair spheres interaction problem [Pournin2001]_ .",
((Real,tc,NaN,,"Contact time"))
((Real,en,NaN,,"Restitution coefficient in normal direction"))
((Real,et,NaN,,"Restitution coefficient in tangential direction"))
=== modified file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp 2014-05-27 07:50:01 +0000
+++ pkg/pfv/FlowEngine.hpp 2014-06-04 17:19:50 +0000
@@ -113,10 +113,10 @@
void updateVolumes (Solver& flow);
void initializeVolumes (Solver& flow);
void boundaryConditions(Solver& flow);
- void updateBCs ( Solver& flow ) {
- if (flow.T[flow.currentTes].maxId>0) boundaryConditions(flow);//avoids crash at iteration 0, when the packing is not bounded yet
+ void updateBCs () {
+ if (solver->tesselation().maxId>0) boundaryConditions(*solver);//avoids crash at iteration 0, when the packing is not bounded yet
else LOG_ERROR("updateBCs not applied");
- flow.pressureChanged=true;}
+ solver->pressureChanged=true;}
void imposeFlux(Vector3r pos, Real flux);
unsigned int imposePressure(Vector3r pos, Real p);
@@ -128,6 +128,7 @@
Real getBoundaryFlux(unsigned int boundary) {return solver->boundaryFlux(boundary);}
Vector3r fluidForce(unsigned int idSph) {
const CGT::CVector& f=solver->T[solver->currentTes].vertex(idSph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
+ Vector3r averageVelocity();
Vector3r shearLubForce(unsigned int id_sph) {
return (solver->shearLubricationForces.size()>id_sph)?solver->shearLubricationForces[id_sph]:Vector3r::Zero();}
@@ -353,6 +354,7 @@
#endif
.def("compTessVolumes",&TemplateFlowEngine::compTessVolumes,"Like TesselationWrapper::computeVolumes()")
.def("volume",&TemplateFlowEngine::getVolume,(boost::python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
+ .def("averageVelocity",&TemplateFlowEngine::averageVelocity,"measure the mean velocity in the period")
)
};
// Definition of functions in a separate file for clarity
=== modified file 'pkg/pfv/FlowEngine.ipp'
--- pkg/pfv/FlowEngine.ipp 2014-05-26 23:23:29 +0000
+++ pkg/pfv/FlowEngine.ipp 2014-06-04 17:19:50 +0000
@@ -222,8 +222,8 @@
flow.meanKStat = meanKStat;
flow.permeabilityMap = permeabilityMap;
flow.fluidBulkModulus = fluidBulkModulus;
-// flow.T[flow.currentTes].Clear();
- flow.T[flow.currentTes].maxId=-1;
+// flow.tesselation().Clear();
+ flow.tesselation().maxId=-1;
flow.xMin = 1000.0, flow.xMax = -10000.0, flow.yMin = 1000.0, flow.yMax = -10000.0, flow.zMin = 1000.0, flow.zMax = -10000.0;
}
@@ -256,16 +256,16 @@
addBoundary ( flow );
triangulate ( flow );
if ( debug ) cout << endl << "Tesselating------" << endl << endl;
- flow.T[flow.currentTes].compute();
+ flow.tesselation().compute();
flow.defineFictiousCells();
// For faster loops on cells define this vector
- flow.T[flow.currentTes].cellHandles.clear();
- flow.T[flow.currentTes].cellHandles.reserve(flow.T[flow.currentTes].Triangulation().number_of_finite_cells());
- FiniteCellsIterator cell_end = flow.T[flow.currentTes].Triangulation().finite_cells_end();
+ flow.tesselation().cellHandles.clear();
+ flow.tesselation().cellHandles.reserve(flow.tesselation().Triangulation().number_of_finite_cells());
+ FiniteCellsIterator cell_end = flow.tesselation().Triangulation().finite_cells_end();
int k=0;
- for ( FiniteCellsIterator cell = flow.T[flow.currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
- flow.T[flow.currentTes].cellHandles.push_back(cell);
+ for ( FiniteCellsIterator cell = flow.tesselation().Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
+ flow.tesselation().cellHandles.push_back(cell);
cell->info().id=k++;}//define unique numbering now, corresponds to position in cellHandles
flow.displayStatistics ();
flow.computePermeability();
@@ -276,9 +276,9 @@
boundaryConditions ( flow );
flow.initializePressure ( pZero );
- if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0 || doInterpolate)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
- if ( waveAction ) flow.applySinusoidalPressure ( flow.T[flow.currentTes].Triangulation(), sineMagnitude, sineAverage, 30 );
- else if (boundaryPressure.size()!=0) flow.applyUserDefinedPressure ( flow.T[flow.currentTes].Triangulation(), boundaryXPos , boundaryPressure);
+ if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0 || doInterpolate)) flow.interpolate ( flow.T[!flow.currentTes], flow.tesselation() );
+ if ( waveAction ) flow.applySinusoidalPressure ( flow.tesselation().Triangulation(), sineMagnitude, sineAverage, 30 );
+ else if (boundaryPressure.size()!=0) flow.applyUserDefinedPressure ( flow.tesselation().Triangulation(), boundaryXPos , boundaryPressure);
if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
}
template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
@@ -320,7 +320,7 @@
}
}
//FIXME idOffset must be set correctly, not the case here (always 0), then we need walls first or it will fail
- idOffset = flow.T[flow.currentTes].maxId+1;
+ idOffset = flow.tesselation().maxId+1;
flow.idOffset = idOffset;
flow.sectionArea = ( flow.xMax - flow.xMin ) * ( flow.zMax-flow.zMin );
flow.vTotal = ( flow.xMax-flow.xMin ) * ( flow.yMax-flow.yMin ) * ( flow.zMax-flow.zMin );
@@ -365,7 +365,7 @@
///Using Tesselation wrapper (faster)
// TesselationWrapper TW;
// if (TW.Tes) delete TW.Tes;
-// TW.Tes = &(flow.T[flow.currentTes]);//point to the current Tes we have in Flowengine
+// TW.Tes = &(flow.tesselation());//point to the current Tes we have in Flowengine
// TW.insertSceneSpheres();//TW is now really inserting in TemplateFlowEngine, using the faster insert(begin,end)
// TW.Tes = NULL;//otherwise, Tes would be deleted by ~TesselationWrapper() at the end of the function.
///Using one-by-one insertion
@@ -374,27 +374,27 @@
if ( !b.exists ) continue;
if ( b.isSphere ) {
if (b.id==ignoredBody) continue;
- flow.T[flow.currentTes].insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
+ flow.tesselation().insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
}
- flow.T[flow.currentTes].redirected=true;//By inserting one-by-one, we already redirected
- flow.shearLubricationForces.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.shearLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.pumpLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.twistLubricationTorques.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.shearLubricationBodyStress.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.normalLubricationForce.resize ( flow.T[flow.currentTes].maxId+1 );
- flow.normalLubricationBodyStress.resize ( flow.T[flow.currentTes].maxId+1 );
+ flow.tesselation().redirected=true;//By inserting one-by-one, we already redirected
+ flow.shearLubricationForces.resize ( flow.tesselation().maxId+1 );
+ flow.shearLubricationTorques.resize ( flow.tesselation().maxId+1 );
+ flow.pumpLubricationTorques.resize ( flow.tesselation().maxId+1 );
+ flow.twistLubricationTorques.resize ( flow.tesselation().maxId+1 );
+ flow.shearLubricationBodyStress.resize ( flow.tesselation().maxId+1 );
+ flow.normalLubricationForce.resize ( flow.tesselation().maxId+1 );
+ flow.normalLubricationBodyStress.resize ( flow.tesselation().maxId+1 );
}
template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::initializeVolumes ( Solver& flow )
{
typedef typename Solver::FiniteVerticesIterator FiniteVerticesIterator;
- FiniteVerticesIterator vertices_end = flow.T[flow.currentTes].Triangulation().finite_vertices_end();
+ FiniteVerticesIterator vertices_end = flow.tesselation().Triangulation().finite_vertices_end();
CGT::CVector Zero(0,0,0);
- for (FiniteVerticesIterator V_it = flow.T[flow.currentTes].Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero;
+ for (FiniteVerticesIterator V_it = flow.tesselation().Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero;
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles)
+ FOREACH(CellHandle& cell, flow.tesselation().cellHandles)
{
switch ( cell->info().fictious() )
{
@@ -448,12 +448,12 @@
epsVolMax=0;
Real totVol=0; Real totDVol=0;
#ifdef YADE_OPENMP
- const long size=flow.T[flow.currentTes].cellHandles.size();
+ const long size=flow.tesselation().cellHandles.size();
#pragma omp parallel for num_threads(ompThreads>0 ? ompThreads : 1)
for(long i=0; i<size; i++){
- CellHandle& cell = flow.T[flow.currentTes].cellHandles[i];
+ CellHandle& cell = flow.tesselation().cellHandles[i];
#else
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
+ FOREACH(CellHandle& cell, flow.tesselation().cellHandles){
#endif
double newVol, dVol;
switch ( cell->info().fictious() ) {
@@ -592,7 +592,7 @@
for ( unsigned int k=0; k<flow.normalLubricationBodyStress.size(); k++) flow.normalLubricationBodyStress[k]=Matrix3r::Zero();
typedef typename Solver::Tesselation Tesselation;
- const Tesselation& Tes = flow.T[flow.currentTes];
+ const Tesselation& Tes = flow.tesselation();
flow.deltaShearVel.clear(); flow.normalV.clear(); flow.deltaNormVel.clear(); flow.surfaceDistance.clear(); flow.onlySpheresInteractions.clear(); flow.normalStressInteraction.clear(); flow.shearStressInteraction.clear();
@@ -726,6 +726,24 @@
}
}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+Vector3r TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::averageVelocity()
+{
+ solver->averageRelativeCellVelocity();
+ Vector3r meanVel ( 0,0,0 );
+ Real volume=0;
+ FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
+ for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
+ //We could also define velocity using cell's center
+// if ( !cell->info().isReal() ) continue;
+ if ( cell->info().isGhost ) continue;
+ for ( int i=0;i<3;i++ )
+ meanVel[i]=meanVel[i]+ ( ( cell->info().averageVelocity() ) [i] * abs ( cell->info().volume() ) );
+ volume+=abs ( cell->info().volume() );
+ }
+ return ( meanVel/volume );
+}
+
#endif //FLOW_ENGINE
#endif /* YADE_CGAL */
=== modified file 'pkg/pfv/PeriodicFlowEngine.cpp'
--- pkg/pfv/PeriodicFlowEngine.cpp 2014-04-17 15:10:23 +0000
+++ pkg/pfv/PeriodicFlowEngine.cpp 2014-06-04 17:19:50 +0000
@@ -85,7 +85,6 @@
Real volumeCellSingleFictious (CellHandle cell);
inline void locateCell(CellHandle baseCell, unsigned int& index, int& baseIndex, FlowSolver& flow, unsigned int count=0);
- Vector3r meanVelocity();
virtual ~PeriodicFlowEngine();
@@ -232,7 +231,7 @@
void PeriodicFlowEngine::triangulate( FlowSolver& flow )
{
- Tesselation& Tes = flow.T[flow.currentTes];
+ Tesselation& Tes = flow.tesselation();
vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
FOREACH ( const posData& b, buffer ) {
if ( !b.exists || !b.isSphere || b.id==ignoredBody) continue;
@@ -332,7 +331,7 @@
PeriFlowTesselation::CellInfo& baseInfo = baseCell->info();
//already located, return FIXME: is inline working correctly? else move this test outside the function, just before the calls
if ( baseInfo.index>0 || baseInfo.isGhost ) return;
- RTriangulation& Tri = flow.T[flow.currentTes].Triangulation();
+ RTriangulation& Tri = flow.tesselation().Triangulation();
Vector3r center ( 0,0,0 );
Vector3i period;
@@ -408,23 +407,6 @@
}
}
-Vector3r PeriodicFlowEngine::meanVelocity()
-{
- solver->averageRelativeCellVelocity();
- Vector3r meanVel ( 0,0,0 );
- Real volume=0;
- FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
- for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
- //We could also define velocity using cell's center
-// if ( !cell->info().isReal() ) continue;
- if ( cell->info().isGhost ) continue;
- for ( int i=0;i<3;i++ )
- meanVel[i]=meanVel[i]+ ( ( cell->info().averageVelocity() ) [i] * abs ( cell->info().volume() ) );
- volume+=abs ( cell->info().volume() );
- }
- return ( meanVel/volume );
-}
-
void PeriodicFlowEngine::updateVolumes (FlowSolver& flow)
{
if ( debug ) cout << "Updating volumes.............." << endl;
@@ -436,7 +418,7 @@
Real totVol0=0;
Real totVol1=0;
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
+ FOREACH(CellHandle& cell, flow.tesselation().cellHandles){
switch ( cell->info().fictious() ) {
case ( 1 ) :
newVol = volumeCellSingleFictious ( cell );
@@ -463,11 +445,11 @@
void PeriodicFlowEngine::initializeVolumes (FlowSolver& flow)
{
- FiniteVerticesIterator vertices_end = flow.T[flow.currentTes].Triangulation().finite_vertices_end();
+ FiniteVerticesIterator vertices_end = flow.tesselation().Triangulation().finite_vertices_end();
CGT::CVector Zero ( 0,0,0 );
- for ( FiniteVerticesIterator V_it = flow.T[flow.currentTes].Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++ ) V_it->info().forces=Zero;
+ for ( FiniteVerticesIterator V_it = flow.tesselation().Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++ ) V_it->info().forces=Zero;
- FOREACH(CellHandle& cell, flow.T[flow.currentTes].cellHandles){
+ FOREACH(CellHandle& cell, flow.tesselation().cellHandles){
switch ( cell->info().fictious() )
{
case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
@@ -492,7 +474,7 @@
if ( debug ) cout << endl << "Added boundaries------" << endl << endl;
triangulate (flow);
if ( debug ) cout << endl << "Tesselating------" << endl << endl;
- flow.T[flow.currentTes].compute();
+ flow.tesselation().compute();
flow.defineFictiousCells();
//FIXME: this is already done in addBoundary(?)
@@ -504,7 +486,7 @@
//This must be done after boundary conditions and initialize pressure, else the indexes are not good (not accounting imposedP): FIXME
unsigned int index=0;
int baseIndex=-1;
- FlowSolver::Tesselation& Tes = flow.T[flow.currentTes];
+ FlowSolver::Tesselation& Tes = flow.tesselation();
Tes.cellHandles.resize(Tes.Triangulation().number_of_finite_cells());
const FiniteCellsIterator cellend=Tes.Triangulation().finite_cells_end();
for ( FiniteCellsIterator cell=Tes.Triangulation().finite_cells_begin(); cell!=cellend; cell++ ){
@@ -523,7 +505,7 @@
flow.displayStatistics ();
//FIXME: check interpolate() for the periodic case, at least use the mean pressure from previous step.
if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0 || doInterpolate)) flow.interpolate ( flow.T[!flow.currentTes], Tes );
-// if ( !first && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], flow.T[flow.currentTes] );
+// if ( !first && (useSolver==0 || fluidBulkModulus>0)) flow.interpolate ( flow.T[!flow.currentTes], flow.tesselation() );
if ( waveAction ) flow.applySinusoidalPressure ( Tes.Triangulation(), sineMagnitude, sineAverage, 30 );
=== modified file 'py/export.py'
--- py/export.py 2014-05-05 14:31:38 +0000
+++ py/export.py 2014-06-06 13:40:28 +0000
@@ -468,6 +468,7 @@
"""
# get list of bodies to export
bodies = self._getBodies(ids,Facet)
+ ids = [b.id for b in bodies]
if not bodies: return
nBodies = len(bodies)
if connectivityTable is None:
