Vous avez reçu un message "Your GitLab account has been locked ..." ? Pas d'inquiétude : lisez cet article https://docs.gricad-pages.univ-grenoble-alpes.fr/help/unlock/

Commit 819f8a80 authored by booncw's avatar booncw
Browse files

Pushing Potential Blocks based on CWBoonEtAl2012 Computers&Geotechnics/CWBoon...

Pushing Potential Blocks based on CWBoonEtAl2012 Computers&Geotechnics/CWBoon 2013 PhD thesis, but using CLP as linear programming Solver.  CLP can be downloaded from Synaptic Package Manager: coinor-clp, coinor-libclp-dev, coinor-libclp1, coinor-libosi1v5.  To enable this, set ENABLE_POTENTIAL_BLOCKS to ON in CMakeLists.txt"
parent ca0e29ac
......@@ -44,6 +44,7 @@ INCLUDE(FindLoki)
INCLUDE(FindPythonModule)
INCLUDE(GNUInstallDirs)
INCLUDE_DIRECTORIES (${CMAKE_SOURCE_DIR})
#===========================================================
......@@ -132,6 +133,7 @@ OPTION(ENABLE_MASK_ARBITRARY "Enable arbitrary precision of bitmask variables (o
OPTION(ENABLE_PROFILING "Enable profiling, e.g. shows some more metrics, which can define bottlenecks of the code (OFF by default)" OFF)
OPTION(ENABLE_POTENTIAL_PARTICLES "Enable potential particles" OFF)
OPTION(USE_QT5 "USE Qt5 for GUI" ON)
OPTION(ENABLE_POTENTIAL_BLOCKS "Enable PotentialBlocks" OFF)
#===========================================================
# Use Eigen3 by default
......@@ -169,7 +171,7 @@ ENDIF((Boost_MAJOR_VERSION EQUAL 1) OR (Boost_MAJOR_VERSION GREATER 1) AND
((Boost_MINOR_VERSION EQUAL 53) OR (Boost_MINOR_VERSION GREATER 53)))
#===========================================================
IF(ENABLE_VTK)
FIND_PACKAGE(VTK COMPONENTS vtkCommonCore vtkIOImage vtkIOXML)
FIND_PACKAGE(VTK COMPONENTS vtkCommonCore vtkIOImage vtkIOXML vtkFiltersCore vtkImagingCore vtkRenderingCore vtkImagingGeneral vtkImagingHybrid)
IF(VTK_FOUND)
INCLUDE_DIRECTORIES(${VTK_INCLUDE_DIRS})
LINK_DIRECTORIES( ${VTK_LIBRARY_DIRS} )
......@@ -414,6 +416,30 @@ ELSE(ENABLE_POTENTIAL_PARTICLES)
ENDIF(ENABLE_POTENTIAL_PARTICLES)
#===========================================================
#===========================================================
IF(ENABLE_POTENTIAL_BLOCKS)
INCLUDE(FindCLP)
FIND_PACKAGE(OpenBlas REQUIRED)
FIND_PACKAGE(LAPACK REQUIRED)
IF(CLP_FOUND AND OPENBLAS_FOUND AND LAPACK_FOUND)
ADD_DEFINITIONS("-DYADE_POTENTIAL_BLOCKS")
INCLUDE_DIRECTORIES(${CLP_INCLUDE_DIR} ${CLP2_INCLUDE_DIR} )
SET(LINKLIBS "${LINKLIBS};${CLP_LIBRARY};${CLP2_LIBRARY};${CLP3_LIBRARY};${OPENBLAS_LIBRARY};${LAPACK_LIBRARY}")
MESSAGE(STATUS "Found CLP")
SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} PotentialBlocks")
ELSE(CLP_FOUND AND OPENBLAS_FOUND AND LAPACK_FOUND)
MESSAGE(STATUS "CLP NOT found")
SET(DISABLED_FEATS "${DISABLED_FEATS} PotentialBlocks")
SET(ENABLE_POTENTIAL_BLOCKS OFF)
ENDIF(CLP_FOUND AND OPENBLAS_FOUND AND LAPACK_FOUND)
ELSE(ENABLE_POTENTIAL_BLOCKS)
SET(DISABLED_FEATS "${DISABLED_FEATS} PotentialBlocks")
ENDIF(ENABLE_POTENTIAL_BLOCKS)
#===========================================================
IF (INSTALL_PREFIX)
SET(CMAKE_INSTALL_PREFIX ${INSTALL_PREFIX})
MESSAGE(WARNING "Use CMAKE_INSTALL_PREFIX option instead of INSTALL_PREFIX! It will be removed soon.")
......@@ -625,3 +651,5 @@ ADD_CUSTOM_COMMAND(
DEPENDS ${YADE_EXEC_PATH}/yade${SUFFIX}
)
#===========================================================
IF(CLP_INCLUDE_DIR AND CLP2_INCLUDE_DIR AND CLP_LIBRARY AND CLP2_LIBRARY AND CLP3_LIBRARY)
SET(CLP_FOUND TRUE)
ELSE(CLP_INCLUDE_DIR AND CLP2_INCLUDE_DIR AND CLP_LIBRARY AND CLP2_LIBRARY AND CLP3_LIBRARY)
FIND_LIBRARY(CLP_LIBRARY NAMES Clp
PATHS
/usr/lib/x86_64-linux-gnu
)
FIND_LIBRARY(CLP2_LIBRARY NAMES OsiClp
PATHS
/usr/lib/x86_64-linux-gnu
)
FIND_LIBRARY(CLP3_LIBRARY NAMES CoinUtils
PATHS
/usr/lib/x86_64-linux-gnu
)
FIND_PATH(
CLP_INCLUDE_DIR ClpSimplexDual.hpp
PATHS
/usr/include/coin
)
FIND_PATH(
CLP2_INCLUDE_DIR coinutils.pc
PATHS
/usr/lib/x86_64-linux-gnu/pkgconfig
)
IF(CLP_INCLUDE_DIR AND CLP2_INCLUDE_DIR AND CLP_LIBRARY AND CLP2_LIBRARY AND CLP3_LIBRARY)
SET(CLP_FOUND TRUE)
MESSAGE(STATUS "Found CLP: ${CLP_INCLUDE_DIR}, ${CLP_LIBRARY}")
ELSE(CLP_INCLUDE_DIR AND CLP2_INCLUDE_DIR AND CLP_LIBRARY AND CLP2_LIBRARY AND CLP3_LIBRARY)
SET(CLP_FOUND FALSE)
MESSAGE(STATUS "CLP not found.")
ENDIF(CLP_INCLUDE_DIR AND CLP2_INCLUDE_DIR AND CLP_LIBRARY AND CLP2_LIBRARY AND CLP3_LIBRARY)
MARK_AS_ADVANCED(CLP_INCLUDE_DIR CLP2_INCLUDE_DIR CLP_LIBRARY CLP2_LIBRARY CLP3_LIBRARY)
ENDIF(CLP_INCLUDE_DIR AND CLP2_INCLUDE_DIR AND CLP_LIBRARY AND CLP2_LIBRARY AND CLP3_LIBRARY)
......@@ -237,6 +237,125 @@ void Clump::updateProperties(const shared_ptr<Body>& clumpBody, unsigned int dis
}
}
void Clump::updatePropertiesNonSpherical(const shared_ptr<Body>& clumpBody, bool intersecting,shared_ptr<Scene> rb){ //FIXME
//LOG_DEBUG("Updating clump #"<<getId()<<" parameters");
//LOG_DEBUG("Updating clump #"<<getId()<<" parameters");
//assert(members.size()>0);
const shared_ptr<State> state(clumpBody->state);
const shared_ptr<Clump> clump(YADE_PTR_CAST<Clump>(clumpBody->shape));
// trivial case
if(clump->members.size()==1){
LOG_DEBUG("Clump of size one will be treated specially.")
MemberMap::iterator I=clump->members.begin();
shared_ptr<Body> subBody=Body::byId(I->first,rb);
//const shared_ptr<RigidBodyParameters>& subRBP(YADE_PTR_CAST<RigidBodyParameters>(subBody->physicalParameters));
State* subState=subBody->state.get();
// se3 of the clump as whole is the same as the member's se3
state->pos=subState->pos;
state->ori=subState->ori;
// relative member's se3 is identity
I->second.position=Vector3r::Zero(); I->second.orientation=Quaternionr::Identity();
state->inertia=subState->inertia;
state->mass=subState->mass;
state->vel=Vector3r::Zero();
state->angVel=Vector3r::Zero();
return;
}
/* quantities suffixed by
g: global (world) coordinates
s: local subBody's coordinates
c: local clump coordinates
*/
double M=0; // mass
Vector3r Sg(0,0,0); // static moment, for getting clump's centroid
Matrix3r Ig(Matrix3r::Zero()), Ic(Matrix3r::Zero()); // tensors of inertia; is upper triangular, zeros instead of symmetric elements
if(intersecting){
LOG_WARN("Self-intersecting clumps not yet implemented, intersections will be ignored.");
intersecting=false;
}
// begin non-intersecting loop here
if(!intersecting){
FOREACH(MemberMap::value_type& I, clump->members){
// I.first is Body::id_t, I.second is Se3r of that body
shared_ptr<Body> subBody=Body::byId(I.first,rb);
State* subState=subBody->state.get();
M+=subState->mass;
Sg+=subState->mass*subState->pos;
// transform from local to global coords
Quaternionr subState_ori_conjugate=subState->ori.conjugate();
Matrix3r Imatrix=Matrix3r::Zero(); Imatrix.diagonal()=subState->inertia;
// TRWM3MAT(Imatrix); TRWM3QUAT(subRBP_orientation_conjugate);
Ig+=Clump::inertiaTensorTranslate(Clump::inertiaTensorRotate(Imatrix,subState_ori_conjugate),subState->mass,-1.*subState->pos);
//TRWM3MAT(Clump::inertiaTensorRotate(Matrix3r(subRBP->inertia),subRBP_orientation_conjugate));
}
}
//TRVAR1(M); TRWM3MAT(Ig); TRWM3VEC(Sg);
assert(M>0);
state->pos=Sg/M; // clump's centroid
// this will calculate translation only, since rotation is zero
Matrix3r Ic_orientG=Clump::inertiaTensorTranslate(Ig, -M /* negative mass means towards centroid */, state->pos); // inertia at clump's centroid but with world orientation
//TRWM3MAT(Ic_orientG);
Matrix3r R_g2c(Matrix3r::Zero()); //rotation matrix
Ic_orientG(1,0)=Ic_orientG(0,1); Ic_orientG(2,0)=Ic_orientG(0,2); Ic_orientG(2,1)=Ic_orientG(1,2); // symmetrize
//TRWM3MAT(Ic_orientG);
matrixEigenDecomposition(Ic_orientG,R_g2c,Ic);
/*! @bug eigendecomposition might be wrong. see http://article.gmane.org/gmane.science.physics.yade.devel/99 for message. It is worked around below, however.
*/
// has NaNs for identity matrix!
//TRWM3MAT(R_g2c);
// set quaternion from rotation matrix
state->ori=Quaternionr(R_g2c); state->ori.normalize();
// now Ic is diagonal
state->inertia=Ic.diagonal();
state->mass=M;
// this block will be removed once EigenDecomposition works for diagonal matrices
//#if 1
// if(isnan(R_g2c(0,0))||isnan(R_g2c(0,1))||isnan(R_g2c(0,2))||isnan(R_g2c(1,0))||isnan(R_g2c(1,1))||isnan(R_g2c(1,2))||isnan(R_g2c(2,0))||isnan(R_g2c(2,1))||isnan(R_g2c(2,2))){
// throw std::logic_error("Clump::updateProperties: NaNs in eigen-decomposition of inertia matrix?!");
// }
//#endif
//TRWM3VEC(state->inertia);
// TODO: these might be calculated from members... but complicated... - someone needs that?!
state->vel=state->angVel=Vector3r::Zero();
clumpBody->setAspherical(state->inertia[0]!=state->inertia[1] || state->inertia[0]!=state->inertia[2]);
// update subBodySe3s; subtract clump orientation (=apply its inverse first) to subBody's orientation
FOREACH(MemberMap::value_type& I, clump->members){
// now, I->first is Body::id_t, I->second is Se3r of that body
shared_ptr<Body> subBody=Body::byId(I.first,rb);
//const shared_ptr<RigidBodyParameters>& subRBP(YADE_PTR_CAST<RigidBodyParameters>(subBody->physicalParameters));
State* subState=subBody->state.get();
I.second.orientation=state->ori.conjugate()*subState->ori;
I.second.position=state->ori.conjugate()*(subState->pos-state->pos);
}
}
void Clump::addNonSpherical(const shared_ptr<Body>& clumpBody, const shared_ptr<Body>& subBody){ //FIXME
Body::id_t subId=subBody->getId();
if(subBody->clumpId!=Body::ID_NONE) throw std::invalid_argument(("Body #"+boost::lexical_cast<string>(subId)+" is already in clump #"+boost::lexical_cast<string>(subBody->clumpId)).c_str());
const shared_ptr<Clump> clump=YADE_PTR_CAST<Clump>(clumpBody->shape);
if(clump->members.count(subId)!=0) throw std::invalid_argument(("Body #"+boost::lexical_cast<string>(subId)+" is already part of this clump #"+boost::lexical_cast<string>(clumpBody->id)).c_str());
clump->members[subId]=Se3r(); // meaningful values will be put in by Clump::updateProperties
subBody->clumpId=clumpBody->id;
clumpBody->clumpId=clumpBody->id; // just to make sure
clumpBody->setBounded(false); // disallow collisions with the clump itself
//LOG_DEBUG("Added body #"<<subId<<" to clump #"<<getId());
}
/*! @brief Recalculates inertia tensor of a body after translation away from (default) or towards its centroid.
*
* @param I inertia tensor in the original coordinates; it is assumed to be upper-triangular (elements below the diagonal are ignored).
......
......@@ -59,6 +59,7 @@ class Clump: public Shape {
static void del(const shared_ptr<Body>& clump, const shared_ptr<Body>& subBody);
//! Recalculate physical properties of Clump.
static void updateProperties(const shared_ptr<Body>& clump, unsigned int discretization);
static void updatePropertiesNonSpherical(const shared_ptr<Body>& clump, bool intersecting,shared_ptr<Scene> rb);//FIXME
//! Calculate positions and orientations of members based on relative Se3; newton pointer (if non-NULL) calls NewtonIntegrator::saveMaximaVelocity
// done as template to avoid cross-dependency between clump and newton (not necessary if all plugins are linked together)
template<class IntegratorT>
......@@ -79,7 +80,7 @@ class Clump: public Shape {
}
}
static void addNonSpherical(const shared_ptr<Body>& clump, const shared_ptr<Body>& subBody); //FIXME
//! update member positions after clump being moved by mouse (in case simulation is paused and engines will not do that).
void userForcedDisplacementRedrawHook(){ throw runtime_error("Clump::userForcedDisplacementRedrawHook not yet implemented (with Clump as subclass of Shape).");}
......@@ -98,7 +99,7 @@ class Clump: public Shape {
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Clump,Shape,"Rigid aggregate of bodies",
((MemberMap,members,,Attr::hidden,"Ids and relative positions+orientations of members of the clump (should not be accessed directly)"))
// ((vector<int>,ids,,Attr::readonly,"Ids of constituent particles (only informative; direct modifications will have no effect)."))
((vector<int>,ids,,Attr::readonly,"Ids of constituent particles (only informative; direct modifications will have no effect).")) //FIXME
,/*ctor*/ createIndex();
,/*py*/ .add_property("members",&Clump::members_get,"Return clump members as {'id1':(relPos,relOri),...}")
);
......
This diff is collapsed.
This diff is collapsed.
#!/usr/local/bin/yade-trunk -x
# -*- encoding=utf-8 -*-
# CWBoon 2015
# Some of the parameters are defunct..
# I use this for jointed rock mass, but you can use it to generate granular particles too
from yade import pack
import math
O.engines=[
ForceResetter(),
InsertionSortCollider([PotentialBlock2AABB()],verletDist=0),
InteractionLoop(
[Ig2_PB_PB_ScGeom()],
[Ip2_FrictMat_FrictMat_KnKsPBPhys(Knormal = 1e8, Kshear = 1e8,useFaceProperties=False,calJointLength=False,twoDimension=True,unitWidth2D=1.0,viscousDamping=0.7)],
[Law2_SCG_KnKsPBPhys_KnKsPBLaw(label='law',neverErase=False)]
#[Ip2_FrictMat_FrictMat_FrictPhys()],
#[Law2_ScGeom_FrictPhys_CundallStrack()]
),
#GravityEngine(gravity=[0,-10,0]),
#GlobalStiffnessTimeStepper(),
NewtonIntegrator(damping=0.0,exactAsphericalRot=False,gravity=[0,-10,0]),
#PotentialBlockVTKRecorder(fileName='/home/chiab/yadeNew/mosek/8Nov/BranchA/scripts2/boon/ComputersGeotechnics/vtk/1000PP',iterPeriod=100,sampleX=50,sampleY=50,sampleZ=50)
]
powderDensity = 10000
distanceToCentre= 0.5
meanSize = 1.0
wallThickness = 0.5*meanSize
O.materials.append(FrictMat(young=150e6,poisson=.4,frictionAngle=radians(0.0),density=powderDensity,label='frictionless'))
lengthOfBase = 9.0*meanSize
heightOfBase = 14.0*meanSize
sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),7.0*heightOfBase,0.5*(lengthOfBase-wallThickness))
sphereRad = sqrt(3.0)*0.5*meanSize
sp.makeCloud(mn,mx,sphereRad,0,100,False)
count= 0
rPP=0.05*meanSize
for s in sp:
b=Body()
radius=2.2
dynamic=True
wire=False
color=[0,0,255.0]
highlight=False
b.shape=PotentialBlock(k=0.2, r=0.05*meanSize, R=1.02*sphereRad, a=[1.0,-1.0,0.0,0.0,0.0,0.0], b=[0.0,0.0,1.0,-1.0,0.0,0.0], c=[0.0,0.0,0.0,0.0,1.0,-1.0], d=[distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP],isBoundary=False,color=color,wire=wire,highlight=highlight,minAabb=Vector3(sphereRad,sphereRad,sphereRad),maxAabb=Vector3(sphereRad,sphereRad,sphereRad),maxAabbRotated=Vector3(sphereRad,sphereRad,sphereRad),minAabbRotated=Vector3(sphereRad,sphereRad,sphereRad),AabbMinMax=False)
length=meanSize
V= 1.0
geomInert=(2./5.)*powderDensity*V*sphereRad**2
utils._commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=s[0], noBound=False, resetState=True, fixed=False)
b.state.pos = s[0] #s[0] stores center
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
O.bodies.append(b)
b.dynamic = True
count =count+1
#v1 = (0, 0, 0.2) (0, 0, 1) (0,0,2.498)
#v2 = (-0.0943, 0.1633, -0.0667) (-0.4714, 0.8165, -0.3334)
#v3 = (0.1886 0 -0.0667) (0.9428, 0, -0.3333)
#edge = 0.3266 1.6333 4.08
#volume = 0.0041 0.5132 8
r=0.1*wallThickness
bbb=Body()
wire=False
color=[0,255,0]
highlight=False
bbb.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.2*lengthOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[lengthOfBase/6.0-r,lengthOfBase/6.0-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/6.0-r,lengthOfBase/6.0-r], id=count,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),maxAabb=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),maxAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),minAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(bbb,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True, fixed=True)
bbb.dynamic=False
bbb.state.pos = [0.0,0,0]
lidID = O.bodies.append(bbb)
b1=Body()
wire=False
color=[0,255,0]
highlight=False
b1.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.2*lengthOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[lengthOfBase/6.0-r,lengthOfBase/6.0-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/6.0-r,lengthOfBase/6.0-r], id=count,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),minAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(b1,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
b1.dynamic=False
b1.state.pos = [lengthOfBase/3.0,0,lengthOfBase/3.0]
O.bodies.append(b1)
b2=Body()
wire=False
color=[0,255,0]
highlight=False
b2.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.2*lengthOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[lengthOfBase/6.0-r,lengthOfBase/6.0-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/6.0-r,lengthOfBase/6.0-r], id=count,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),minAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(b2,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
b2.dynamic=False
b2.state.pos = [-lengthOfBase/3.0,0,lengthOfBase/3.0]
O.bodies.append(b2)
b3=Body()
wire=False
color=[0,255,0]
highlight=False
b3.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.2*lengthOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[lengthOfBase/6.0-r,lengthOfBase/6.0-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/6.0-r,lengthOfBase/6.0-r], id=count,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),minAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(b3,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
b3.dynamic=False
b3.state.pos = [0.0,0,lengthOfBase/3.0]
O.bodies.append(b3)
b4=Body()
wire=False
color=[0,255,0]
highlight=False
b4.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.2*lengthOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[lengthOfBase/6.0-r,lengthOfBase/6.0-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/6.0-r,lengthOfBase/6.0-r], id=count,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),minAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(b4,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
b4.dynamic=False
b4.state.pos = [lengthOfBase/3.0,0,-lengthOfBase/3.0]
O.bodies.append(b4)
b5=Body()
wire=False
color=[0,255,0]
highlight=False
b5.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.2*lengthOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[lengthOfBase/6.0-r,lengthOfBase/6.0-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/6.0-r,lengthOfBase/6.0-r], id=count,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),minAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(b5,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
b5.dynamic=False
b5.state.pos = [0.0,0,-lengthOfBase/3.0]
O.bodies.append(b5)
b6=Body()
wire=False
color=[0,255,0]
highlight=False
b6.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.2*lengthOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[lengthOfBase/6.0-r,lengthOfBase/6.0-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/6.0-r,lengthOfBase/6.0-r], id=count,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),minAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(b6,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
b6.dynamic=False
b6.state.pos = [-lengthOfBase/3.0,0.0,-lengthOfBase/3.0]
O.bodies.append(b6)
b7=Body()
wire=False
color=[0,255,0]
highlight=False
b7.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.2*lengthOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[lengthOfBase/6.0-r,lengthOfBase/6.0-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/6.0-r,lengthOfBase/6.0-r], id=count,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),minAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(b7,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
b7.dynamic=False
b7.state.pos = [-lengthOfBase/3.0,0.0,0.0]
O.bodies.append(b7)
b8=Body()
wire=False
color=[0,255,0]
highlight=False
b8.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.2*lengthOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[lengthOfBase/6.0-r,lengthOfBase/6.0-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/6.0-r,lengthOfBase/6.0-r], id=count,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabb=1.02*Vector3(lengthOfBase/6.0,0.4*wallThickness,lengthOfBase/6.0),maxAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0),minAabbRotated=1.02*Vector3(lengthOfBase/6.0,0.5*wallThickness,lengthOfBase/6.0))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(b8,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
b8.dynamic=False
b8.state.pos = [lengthOfBase/3.0,0.0,0.0]
O.bodies.append(b8)
bA=Body()
wire=False
color=[0,255,0]
highlight=False
bA.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.5*heightOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], id=count+1,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(0.3*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),maxAabb=1.02*Vector3(0.3*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),maxAabbRotated=1.02*Vector3(0.5*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),minAabbRotated=1.02*Vector3(0.5*wallThickness,0.5*heightOfBase,0.5*lengthOfBase))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(bA,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
bA.dynamic=False
bA.state.pos = [0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bA)
bB=Body()
wire=False
color=[0,255,0]
highlight=False
bB.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.5*heightOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], id=count+2,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(0.3*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),maxAabb=1.02*Vector3(0.3*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),maxAabbRotated=1.02*Vector3(0.5*wallThickness,0.5*heightOfBase,0.5*lengthOfBase),minAabbRotated=1.02*Vector3(0.5*wallThickness,0.5*heightOfBase,0.5*lengthOfBase))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(bB,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
bB.dynamic=False
bB.state.pos = [-0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bB)
bC=Body()
wire=False
color=[0,255,0]
highlight=False
bC.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.5*heightOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*lengthOfBase-r,0.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], id=count+3,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.3*wallThickness),maxAabb=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.3*wallThickness),maxAabbRotated=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.5*wallThickness),minAabbRotated=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.5*wallThickness))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(bC,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
bC.dynamic=False
bC.state.pos = [0,0.5*heightOfBase,0.5*lengthOfBase]
O.bodies.append(bC)
bD=Body()
wire=False
color=[0,255,0]
highlight=False
bD.shape=PotentialBlock(k=0.1, r=0.1*wallThickness, R=0.5*heightOfBase,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*lengthOfBase-r,0.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], id=count+4,isBoundary=True,isBoundaryPlane=[True,True,True,True,True,True],color=color ,wire=wire,highlight=highlight,AabbMinMax=True, minAabb=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.3*wallThickness),maxAabb=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.3*wallThickness),maxAabbRotated=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.5*wallThickness),minAabbRotated=1.02*Vector3(0.5*lengthOfBase,0.5*heightOfBase,0.5*wallThickness))
length=lengthOfBase
V=lengthOfBase*lengthOfBase*wallThickness
geomInert=(1./6.)*V*length*wallThickness
utils._commonBodySetup(bD,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], noBound=False, resetState=True,fixed=True)
bD.dynamic=False
bD.state.pos = [0.0,0.5*heightOfBase,-0.5*lengthOfBase]
O.bodies.append(bD)
escapeNo=0
def myAddPlotData():
global escapeNo
global wallThickness
global meanSize
uf=utils.unbalancedForce()
if isnan(uf):
uf = 1.0
KE = utils.kineticEnergy()
for b in O.bodies:
if b.state.pos[1] < -5.0*meanSize and b.dynamic==True:
escapeNo = escapeNo+1
O.bodies.erase(b.id)
if O.iter>12000:
removeLid()
plot.addData(timeStep1=O.iter,timeStep2=O.iter,timeStep3=O.iter,timeStep4=O.iter,time=O.time,unbalancedForce=uf,kineticEn=KE,outsideNo=escapeNo)
from yade import plot
plot.plots={'timeStep1':('unbalancedForce'),'timeStep2':('kineticEn'),'time':('outsideNo')}
O.engines=O.engines+[PyRunner(iterPeriod=10,command='myAddPlotData()')]
def removeLid():
global lidID
if (O.bodies[lidID]):
O.bodies.erase(lidID)
O.engines=O.engines+[PotentialBlockVTKRecorder(fileName='/home/boon/yadeRev/trunk/examples/PotentialParticles/vtk/cubeScaled',iterPeriod=1000,sampleX=50,sampleY=50,sampleZ=50)]