yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06085
[Branch ~yade-dev/yade/trunk] Rev 2532: 1. Make clumps subclass of Shape (instead of Body)
------------------------------------------------------------
revno: 2532
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Sat 2010-11-06 12:33:09 +0100
message:
1. Make clumps subclass of Shape (instead of Body)
2. Favicon in the batch web interface
3. Fix for clang compilation
added:
py/tests/clump.py
modified:
core/main/yade-batch.in
gui/qt4/yade-favicon.png
lib/smoothing/WeightedAverage2d.hpp
pkg/common/GravityEngines.cpp
pkg/dem/Clump.cpp
pkg/dem/Clump.hpp
pkg/dem/NewtonIntegrator.cpp
py/SConscript
py/remote.py
py/tests/__init__.py
py/tests/omega.py
py/wrapper/yadeWrapper.cpp
--
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 'core/main/yade-batch.in'
--- core/main/yade-batch.in 2010-10-18 19:27:07 +0000
+++ core/main/yade-batch.in 2010-11-06 11:33:09 +0000
@@ -128,9 +128,16 @@
from BaseHTTPServer import BaseHTTPRequestHandler,HTTPServer
import socket,re
class HttpStatsServer(BaseHTTPRequestHandler):
+ favicon=None # binary favicon, created when first requested
def do_GET(self):
if not self.path or self.path=='/': self.sendGlobal()
else:
+ if self.path=='/favicon.ico':
+ if not self.__class__.favicon:
+ import base64
+ self.__class__.favicon=base64.b64decode(yade.remote.b64favicon)
+ self.sendHttp(self.__class__.favicon,contentType='image/vnd.microsoft.icon')
+ return
jobMatch=re.match('/jobs/([0-9]+)/(.*)',self.path)
if not jobMatch:
self.send_error(404,self.path); return
=== modified file 'gui/qt4/yade-favicon.png'
Binary files gui/qt4/yade-favicon.png 2010-07-23 22:17:33 +0000 and gui/qt4/yade-favicon.png 2010-11-06 11:33:09 +0000 differ
=== modified file 'lib/smoothing/WeightedAverage2d.hpp'
--- lib/smoothing/WeightedAverage2d.hpp 2010-11-02 12:02:13 +0000
+++ lib/smoothing/WeightedAverage2d.hpp 2010-11-06 11:33:09 +0000
@@ -2,6 +2,14 @@
#pragma once
+// temporary fix:
+// clang #includes first /usr/include/c++/4.4/fenv.h, which through some obscure #ifdefs does not eventually include any code
+// why...?
+#ifdef __clang__
+ #include</usr/include/fenv.h>
+#endif
+
+
#include<iostream>
#include<vector>
#include<cstdlib>
@@ -13,9 +21,9 @@
#include<boost/lexical_cast.hpp>
#include<boost/python.hpp>
#include<yade/extra/boost_python_len.hpp>
-
#include<boost/math/distributions/normal.hpp>
+
#ifndef FOREACH
#define FOREACH BOOST_FOREACH
#endif
=== modified file 'pkg/common/GravityEngines.cpp'
--- pkg/common/GravityEngines.cpp 2010-11-05 14:13:01 +0000
+++ pkg/common/GravityEngines.cpp 2010-11-06 11:33:09 +0000
@@ -15,12 +15,6 @@
YADE_PLUGIN((GravityEngine)(CentralGravityEngine)(AxialGravityEngine)(HdapsGravityEngine));
void GravityEngine::action(){
- /* skip bodies that are within a clump;
- * even if they are marked isDynamic==false, forces applied to them are passed to the clump, which is dynamic;
- * and since clump is a body with mass equal to the sum of masses of its components, it would have gravity applied twice.
- *
- * The choice is to skip (b->isClumpMember()) or (b->isClump()). We rather skip members,
- * since that will apply smaller number of forces. */
#ifdef YADE_OPENMP
const BodyContainer& bodies=*(scene->bodies.get());
const long size=(long)bodies.size();
@@ -30,7 +24,8 @@
#else
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
#endif
- if(!b || b->isClumpMember()) continue;
+ // skip clumps, only apply forces on their constituents
+ if(!b || b->isClump()) continue;
scene->forces.addForce(b->getId(),gravity*b->state->mass);
}
}
@@ -38,7 +33,7 @@
void CentralGravityEngine::action(){
const Vector3r& centralPos=Body::byId(centralBody)->state->pos;
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
- if(!b || b->isClumpMember() || b->getId()==centralBody) continue; // skip clump members and central body
+ if(!b || b->isClump() || b->getId()==centralBody) continue; // skip clumps and central body
Real F=accel*b->state->mass;
Vector3r toCenter=centralPos-b->state->pos; toCenter.normalize();
scene->forces.addForce(b->getId(),F*toCenter);
@@ -48,7 +43,7 @@
void AxialGravityEngine::action(){
FOREACH(const shared_ptr<Body>&b, *scene->bodies){
- if(!b || b->isClumpMember()) continue;
+ if(!b || b->isClump()) continue;
/* http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html */
const Vector3r& x0=b->state->pos;
const Vector3r& x1=axisPoint;
=== modified file 'pkg/dem/Clump.cpp'
--- pkg/dem/Clump.cpp 2010-10-18 19:27:07 +0000
+++ pkg/dem/Clump.cpp 2010-11-06 11:33:09 +0000
@@ -1,4 +1,4 @@
-// (c) 2007,2009 Vaclav Smilauer <eudoxos@xxxxxxxx>
+// (c) 2007-2010 Vaclav Smilauer <eudoxos@xxxxxxxx>
#include"Clump.hpp"
#include<algorithm>
@@ -6,122 +6,60 @@
#include<yade/core/BodyContainer.hpp>
#include<yade/core/State.hpp>
-YADE_PLUGIN((Clump)/*(ClumpMemberMover)*/);
+YADE_PLUGIN((Clump));
CREATE_LOGGER(Clump);
-#if 0
-CREATE_LOGGER(ClumpMemberMover);
-/**************************************************************************************
- ************************************* ClumpMemberMover ******************************
- **************************************************************************************/
-
-/*! We only call clump's method, since it belongs there logically. It makes encapsulation of private members nicer, too.
- * _param pp passed by the dispatcher
- * _param clump passed by the dispatcher
- */
-void ClumpMemberMover::action(){
- for(BodyContainer::iterator I=scene->bodies->begin(); I!=scene->bodies->end(); ++I){
- shared_ptr<Body> b = *I;
- if(b->isClump()){
- //LOG_TRACE("Applying movement to clump #"<<b->getId());
- static_cast<Clump*>(b.get())->moveMembers();
- }
- }
- //if(!clump->isDynamic()) return; // perhaps clump that has been desactivated?!
-}
-#endif
-
-/**************************************************************************************
- ******************************************** Clump ***********************************
- **************************************************************************************/
-
-/*! @pre Body must be dynamic.
- * @pre Body must not be part or this clump already.
- * @pre Body must have valid (non-NULL) Body::physicalParameters
- * @pre Body must have id that is smaller than the clump's id (reason: processing order in NewtonIntegrator)
- */
-void Clump::add(Body::id_t subId){
- shared_ptr<Body> subBody=Body::byId(subId);
-
- // preconditions
- //assert(subBody->isDynamic);
- assert(state);
- assert(members.count(subId)==0);
- assert(subId<getId());
-
- // begin actual setup
- subBody->clumpId=getId();
- subBody->setDynamic(false);
- // for now, push just unitialized se3; will be updated by updateProperties
- members[subId]=Se3r();
-
- clumpId=getId(); // identifies a clump
-
+void Clump::add(const shared_ptr<Body>& clumpBody, const shared_ptr<Body>& subBody){
+ Body::id_t subId=subBody->getId();
+ if(subBody->clumpId!=Body::ID_NONE) throw std::invalid_argument(("Body #"+lexical_cast<string>(subId)+" is already in clump #"+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 #"+lexical_cast<string>(subId)+" is already part of this clump #"+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());
}
-/*! @pre Body with given id must be in the clump.
- */
-void Clump::del(Body::id_t subId){
- // erase the subBody; removing body that is not part of the clump is error
- assert(members.erase(subId)==1);
- // restore body's internal parameters;
- shared_ptr<Body> subBody=Body::byId(subId);
+void Clump::del(const shared_ptr<Body>& clumpBody, const shared_ptr<Body>& subBody){
+ // erase the subBody; removing body that is not part of the clump throws
+ const shared_ptr<Clump> clump=YADE_PTR_CAST<Clump>(clumpBody->shape);
+ if(clump->members.erase(subBody->id)!=1) throw std::invalid_argument(("Body #"+lexical_cast<string>(subBody->id)+" not part of clump #"+lexical_cast<string>(clumpBody->id)+"; not removing.").c_str());
subBody->clumpId=Body::ID_NONE;
- subBody->setDynamic(true);
LOG_DEBUG("Removed body #"<<subId<<" from clump #"<<getId());
}
-/*! @brief Calculate positions and orientations of members based on Clump's Se3; resets acceleration and angularAccelration to zero.
- *
- * This method is called by the ClumpMemberMover engine after each timestep.
- * @note Velocities of members are not updated, since members have isdynamic==false. It is possible, though, that someone needs to have a moving clump that is later broken apart and that liberated particle continue to move in the same way as they did within the clump. In that case, this will have to be completed.
- */
-void Clump::moveMembers(){
- //const Se3r& mySe3(physicalParameters->se3);
- //const shared_ptr<RigidBodyParameters>& myRBP=static_pointer_cast<RigidBodyParameters>(physicalParameters);
- for(Clump::memberMap::iterator I=members.begin(); I!=members.end(); I++){
- // now, I->first is Body::id_t, I->second is Se3r of that body in the clump
- shared_ptr<Body> member=Body::byId(I->first);
- //const shared_ptr<RigidBodyParameters>& subRBP(YADE_PTR_CAST<RigidBodyParameters>(member->physicalParameters));
- State* subState=member->state.get();
+void Clump::moveMembers(const shared_ptr<Body>& clumpBody, Scene* scene){
+ const shared_ptr<Clump>& clump=YADE_PTR_CAST<Clump>(clumpBody->shape);
+ const shared_ptr<State>& clumpState=clumpBody->state;
+
+ FOREACH(MemberMap::value_type& B, clump->members){
+ // B->first is Body::id_t, B->second is local Se3r of that body in the clump
+ State* subState=Body::byId(B.first,scene)->state.get();
+ Vector3r& subPos(B.second.position); Quaternionr& subOri(B.second.orientation);
//LOG_TRACE("Old #"<<I->first<<"position: "<<subRBP->se3.position);
- subState->pos=state->pos+state->ori*I->second.position;
- subState->ori=state->ori*I->second.orientation;
+
+ // position update
+ subState->pos=clumpState->pos+clumpState->ori*subPos;
+ subState->ori=clumpState->ori*subOri;
+ // velocity update
+ subState->vel=clumpState->vel+clumpState->angVel.cross(subState->pos-clumpState->pos);
+ subState->angVel=clumpState->angVel;
+
//LOG_TRACE("New #"<<I->first<<"position: "<<subRBP->se3.position);
//LOG_TRACE("Clump #"<<getId()<<" moved #"<<I->first<<".");
-
- //! FIXME: we set velocity because of damping here; but since positions are integrated after all forces applied, these velocities will be used in the NEXT step for CundallNonViscousDamping. Does that matter?!
- //subState->vel=state->vel+state->angVel.Cross(I->second.position);
- subState->vel=state->vel+state->angVel.cross(subState->pos-state->pos);
- subState->angVel=state->angVel;
}
- /* @bug Temporarily we reset acceleration and angularAcceleration of the clump here;
- * should be a new negine that will take care of that?
- */
- // const shared_ptr<RigidBodyParameters>& clumpRBP(YADE_PTR_CAST<RigidBodyParameters>(physicalParameters));
- #if 0
- if(Omega::instance().getScene()->iter%50==0){
- Real Erot=.5*clumpRBP->inertia[0]*pow(clumpRBP->angularVelocity[0],2)+.5*clumpRBP->inertia[1]*pow(clumpRBP->angularVelocity[1],2)+.5*clumpRBP->inertia[2]*pow(clumpRBP->angularVelocity[2],2);
- Real Etrans=.5*clumpRBP->mass*pow(clumpRBP->velocity.norm(),2);
- // (0,0,1) is gravity acceleration
- Real Epot=clumpRBP->se3.position.Dot(Vector3r(0,0,1))*clumpRBP->mass;
- LOG_TRACE("##"<<clumpId<<" energy "<<Erot+Etrans+Epot<<"\tv "<<Etrans<<"\tw "<<Erot<<"\tp "<<Epot);
- }
- #endif
-
- state->accel=state->angAccel=Vector3r::Zero();
}
/*! Clump's se3 will be updated (origin at centroid and axes coincident with principal inertia axes) and subSe3 modified in such a way that members positions in world coordinates will not change.
- The clump values that are changed are:
- -# Clump::members (holds position and orientation in clump's coordinate system)
- -# Clump::physicalParameters->mass (sum of masses of all members)
- -# Clump::physicalParameters->inertia (inertia of the aggregate - in clump coordinate system)
- -# Clump::physicalParameters->se3 (position and orientation of the clump; it is such that absolute positions and orientation of members will not chage)
-
- The algorithm is as follows:
+ TODO: numerical integration of inertia based on regular space sampling with relative tolerance WRT minimum sphere radius
+
+ Note: velocities and angularVelocities of constituents are zeroed.
+
+ OLD DOCS (will be cleaned up):
+
-# Clump::members values and Clump::physicalParameters::se3 are invalid from this point
-# M=0; S=vector3r(0,0,0); I=zero tensor; (ALL calculations are in world coordinates!)
-# loop over Clump::members (position x_i, mass m_i, inertia at subBody's centroid I_i) [this loop will be replaced by numerical integration (rasterization) for the intersecting case; the rest will be the same]
@@ -134,34 +72,18 @@
-# se3->orientation=quaternion(rotation_matrix); se3->position=clumpPos
-# update subSe3s
- @todo I \em think that the order of transformation of inertia is:
- - from local to global: first rotate, then translate;
- - from global to local: first translate, then rotate,
- since rotation must be done with origin at the centroid... This needs to be verified, though.
- @todo all the rest of this routine needs to be verified!
- @todo implement the loop for intersecting bodies (may cut'n'paste from slum code, but that will work for spheres only!)
+*/
- @note User is responsible for calling this function when appropriate (after adding/removing bodies and before any subsequent simulation). This function can be rather slow by virtue of numerical integration.
- @note subBodie's velocities are not taken into account. This means that clump will be at still after being created, even if its composing particles did have some velocities. If this is concern for someone, it needs to be completed in the code below. See Clump::moveMembers for complementary issue.
- @todo Needs to be tested for physical correctness
- @param intersecting if true, evaluate mass and inertia numerically; otherwise, use analytical methods (parallel axes theorem) which disregard any intersections, but are much faster. */
-void Clump::updateProperties(bool intersecting){
+void Clump::updateProperties(const shared_ptr<Body>& clumpBody, bool intersecting){
LOG_DEBUG("Updating clump #"<<getId()<<" parameters");
assert(members.size()>0);
-
- /* 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
- Matrix3r Ig(Matrix3r::Zero()), Ic(Matrix3r::Zero()); // tensors of inertia; is upper triangular, zeros instead of symmetric elements
- //Se3r& mySe3(physicalParameters->se3);
- //const shared_ptr<RigidBodyParameters>& clumpRBP(YADE_PTR_CAST<RigidBodyParameters>(physicalParameters));
-
- if(members.size()==1){
+ 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=members.begin();
+ MemberMap::iterator I=clump->members.begin();
shared_ptr<Body> subBody=Body::byId(I->first);
//const shared_ptr<RigidBodyParameters>& subRBP(YADE_PTR_CAST<RigidBodyParameters>(subBody->physicalParameters));
State* subState=subBody->state.get();
@@ -177,33 +99,38 @@
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;}
+ intersecting=false;
+ }
// begin non-intersecting loop here
if(!intersecting){
- for(memberMap::iterator I=members.begin(); I!=members.end(); I++){
- // now, I->first is Body::id_t, I->second is Se3r of that body
- shared_ptr<Body> subBody=Body::byId(I->first);
- //const shared_ptr<RigidBodyParameters>& subRBP(YADE_PTR_CAST<RigidBodyParameters>(subBody->physicalParameters));
+ 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);
State* subState=subBody->state.get();
M+=subState->mass;
Sg+=subState->mass*subState->pos;
// transform from local to global coords
- // FIXME: verify this!
- Quaternionr subState_orientation_conjugate=subState->ori.conjugate();
- Matrix3r Imatrix=Matrix3r::Zero(); for(int i=0;i<3;i++) Imatrix(i,i)=subState->inertia[i]; // eigen: Imatrix.diagonal=subState->inertia;
+ 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_orientation_conjugate),subState->mass,-1.*subState->pos);
-
+ 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);
- if(M==0){ state->mass=0; state->inertia=Vector3r(0,0,0); setDynamic(false); return; }
+ 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
@@ -219,28 +146,17 @@
// has NaNs for identity matrix!
TRWM3MAT(R_g2c);
- // these two should give the same result!
- //TRWM3MAT(Ic);
- //TRWM3MAT(Clump::inertiaTensorRotate(Ic_orientG,R_g2c));
-
// set quaternion from rotation matrix
state->ori=Quaternionr(R_g2c); state->ori.normalize();
// now Ic is diagonal
- state->inertia=Vector3r(Ic(0,0),Ic(1,1),Ic(2,2));
+ 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))){
- LOG_FATAL("EigenDecomposition gave some NaNs, we will use imaginary values for clump inertia and orientation instead. I thought this may happen only for 1-member clumps which are now treated specially. Something is broken!");
- //FIXME: since EigenDecomposition is broken, use inertia of the first body instead;
- //!!!!! note that this is HIGHLY incorrect for all non-single clumps !!!!!
- memberMap::iterator I=members.begin();
- shared_ptr<Body> subBody=Body::byId(I->first);
- state->inertia=subBody->state->inertia*10.; // 10 is arbitrary; just to have inertia of clump bigger
- // orientation of the clump is broken as well, since is result of EigenDecomposition as well (rotation matrix)
- state->ori=Quaternionr::Identity();
+ throw std::logic_error("Clump::updateProperties: NaNs in eigen-decomposition of inertia matrix?!");
}
#endif
TRWM3VEC(state->inertia);
@@ -248,17 +164,16 @@
// TODO: these might be calculated from members... but complicated... - someone needs that?!
state->vel=state->angVel=Vector3r::Zero();
- if(state->inertia[0]!=state->inertia[1] || state->inertia[0]!=state->inertia[2]) this->setAspherical(true);
+ 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
- // Conjugate is equivalent to Inverse for normalized quaternions
- for(memberMap::iterator I=members.begin(); I!=members.end(); I++){
+ 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);
+ shared_ptr<Body> subBody=Body::byId(I.first);
//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);
+ I.second.orientation=state->ori.conjugate()*subState->ori;
+ I.second.position=state->ori.conjugate()*(subState->pos-state->pos);
}
}
@@ -273,14 +188,15 @@
Matrix3r Clump::inertiaTensorTranslate(const Matrix3r& I,const Real m, const Vector3r& off){
Real ooff=off.dot(off);
Matrix3r I2=I;
- //TRWM3VEC(off); TRVAR2(ooff,m); TRWM3MAT(I);
+ // TODO: replace by nicer eigen code
+
// translation away from centroid
/* I^c_jk=I'_jk-M*(delta_jk R.R - R_j*R_k) [http://en.wikipedia.org/wiki/Moments_of_inertia#Parallel_axes_theorem] */
Matrix3r dI; dI<</* dIxx */ ooff-off[0]*off[0], /* dIxy */ -off[0]*off[1], /* dIxz */ -off[0]*off[2],
/* sym */ 0., /* dIyy */ ooff-off[1]*off[1], /* dIyz */ -off[1]*off[2],
/* sym */ 0., /* sym */ 0., /* dIzz */ ooff-off[2]*off[2];
I2+=m*dI;
- I2(1,0)=I2(0,1); I2(2,0)=I2(0,2); I2(2,1)=I2(1,2);
+ I2(1,0)=I2(0,1); I2(2,0)=I2(0,2); I2(2,1)=I2(1,2); // symmetrize
//TRWM3MAT(I2);
return I2;
}
=== modified file 'pkg/dem/Clump.hpp'
--- pkg/dem/Clump.hpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/Clump.hpp 2010-11-06 11:33:09 +0000
@@ -4,6 +4,7 @@
#include<vector>
#include<map>
+#include<stdexcept>
#include<yade/core/Body.hpp>
#include<yade/lib-base/Logging.hpp>
#include<yade/lib-base/Math.hpp>
@@ -51,22 +52,19 @@
*/
-class Clump: public Body {
- //! mapping of body IDs to their relative positions; replaces members and subSe3s;
+class Clump: public Shape {
public:
- typedef std::map<Body::id_t,Se3r> memberMap;
- virtual ~Clump(){};
- //! \brief add Body to the Clump
- void add(Body::id_t);
- //! \brief remove Body from the Clump
- void del(Body::id_t);
+ typedef std::map<Body::id_t,Se3r> MemberMap;
+
+ static void add(const shared_ptr<Body>& clump, const shared_ptr<Body>& subBody);
+ static void del(const shared_ptr<Body>& clump, const shared_ptr<Body>& subBody);
//! Recalculate physical properties of Clump.
- void updateProperties(bool intersecting);
+ static void updateProperties(const shared_ptr<Body>& clump, bool intersecting);
//! Calculate positions and orientations of members based on my own Se3.
- void moveMembers();
+ static void moveMembers(const shared_ptr<Body>& clump, Scene* scene);
//! update member positions after clump being moved by mouse (in case simulation is paused and engines will not do that).
- void userForcedDisplacementRedrawHook(){moveMembers();}
- private: // may be made public, but once properly tested...
+ void userForcedDisplacementRedrawHook(){ throw runtime_error("Clump::userForcedDisplacementRedrawHook not yet implemented (with Clump as subclass of Shape).");}
+
//! Recalculates inertia tensor of a body after translation away from (default) or towards its centroid.
static Matrix3r inertiaTensorTranslate(const Matrix3r& I,const Real m, const Vector3r& off);
//! Recalculate body's inertia tensor in rotated coordinates.
@@ -74,24 +72,12 @@
//! Recalculate body's inertia tensor in rotated coordinates.
static Matrix3r inertiaTensorRotate(const Matrix3r& I, const Quaternionr& rot);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Clump,Body,"Rigid aggregate of bodies",
- ((memberMap,members,,,"Ids and relative positions+orientations of members of the clump (should not be accessed directly)")),
- /*ctor*/setDynamic(true); /* possible source of crash is setDynamic manipulates Body::State! */
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(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)."))
+ ,/*ctor*/ createIndex();
);
DECLARE_LOGGER;
+ REGISTER_CLASS_INDEX(Clump,Shape);
};
REGISTER_SERIALIZABLE(Clump);
-
-// NewtonIntegrator calls Clump::moveMembers directly, no need for a dedicated engine
-#if 0
-/*! Update ::Clump::members positions so that the Clump behaves as a rigid body. */
-class ClumpMemberMover: public PartialEngine {
- public:
- //! Interates over scene->bodies and calls Clump::moveSubBodies() for clumps.
- virtual void action();
- virtual ~ClumpMemberMover(){};
- YADE_CLASS_BASE_DOC(ClumpMemberMover,PartialEngine,"Update Clump::members positions and orientations so that Clump behaves as rigid body. This engine is only used internally by NewtonIntegrator, not directly.");
- DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(ClumpMemberMover);
-#endif
=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2010-10-30 17:08:51 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2010-11-06 11:33:09 +0000
@@ -138,6 +138,8 @@
leapfrogAsphericalRotate(scene,state,id,dt,m);
}
} else if (b->isClump()){
+ // reset acceleration of the clump itself; computed from accels on constituents
+ state->accel=state->angAccel=Vector3r::Zero();
// clump mass forces
const Vector3r& f=scene->forces.getForce(id);
Vector3r dLinAccel=f/state->mass;
@@ -148,7 +150,7 @@
// sum force on clump memebrs
// exactAsphericalRot enabled and clump is aspherical
if (exactAsphericalRot && b->isAspherical()){
- FOREACH(Clump::memberMap::value_type mm, static_cast<Clump*>(b.get())->members){
+ FOREACH(Clump::MemberMap::value_type mm, static_cast<Clump*>(b->shape.get())->members){
const Body::id_t& memberId=mm.first;
const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(member->isClumpMember());
State* memberState=member->state.get();
@@ -163,7 +165,7 @@
Vector3r dAngAccel=M.cwise()/state->inertia;
cundallDamp(dt,M,state->angVel,dAngAccel);
state->angAccel+=dAngAccel;
- FOREACH(Clump::memberMap::value_type mm, static_cast<Clump*>(b.get())->members){
+ FOREACH(Clump::MemberMap::value_type mm, static_cast<Clump*>(b->shape.get())->members){
const Body::id_t& memberId=mm.first;
const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(member->isClumpMember());
State* memberState=member->state.get();
@@ -175,7 +177,7 @@
leapfrogTranslate(scene,state,id,dt);
leapfrogSphericalRotate(scene,state,id,dt);
}
- static_cast<Clump*>(b.get())->moveMembers();
+ Clump::moveMembers(b,scene);
}
saveMaximaVelocity(scene,id,state);
=== modified file 'py/SConscript'
--- py/SConscript 2010-10-20 15:38:06 +0000
+++ py/SConscript 2010-11-06 11:33:09 +0000
@@ -43,11 +43,7 @@
env.SharedLibrary('_customConverters',['wrapper/customConverters.cpp'],SHLIBPREFIX='',LIBS=env['LIBS']+linkPlugins(['Dispatching','MatchMaker']+(['GLDrawFunctors','OpenGLRenderer'] if 'opengl' in env['features'] else [])))
])
env.Install('$LIBDIR/py/yade/tests',[
- env.File('__init__.py','tests'),
- env.File('wrapper.py','tests'),
- env.File('omega.py','tests'),
- env.File('pbc.py','tests'),
- env.File('cohesive-chain.py','tests'),
+ env.File(env.Glob('tests/*.py'),'tests'),
])
# 3rd party modules:
=== modified file 'py/remote.py'
--- py/remote.py 2010-10-19 10:04:27 +0000
+++ py/remote.py 2010-11-06 11:33:09 +0000
@@ -18,6 +18,9 @@
plotImgFormat,plotImgMimetype='png','image/png'
#plotImgFormat,plotImgMimetype='svg','image/svg+xml'
+# yade 16x16 favicon, encoded with base64
+b64favicon='AAABAAEAEBAAAAAAAABoBQAAFgAAACgAAAAQAAAAIAAAAAEACAAAAAAAQAEAAAAAAAAAAAAAAAAAAAAAAAAAAAAA////AAiD7wCjjY8AAGhoAAwMwwBuVO0AY1RWAJG+4QAHpKMAHiRTAD88vQDW0t0AQp7XAHBwtAAMvf8ANWx0AIKK5QABycgAqZ/AADc1iwA6Oe0AM09QABskJQA+dNoAXWt/AAhAPwAAhYMApp3nAMm+/ABYWM0Avr7IAG+v5gBdS6wAint9AAdFXwApKdoAgGvnAIGBwwCjoaMANSssAEpERQByZGYAiaTmAJiQ1AAhPDkASkjfAAa0tgAGrO4APE2OABNdXgBlZdwAuqvVAA8qPAAAt9YAcmrKAF5etwAvL8YAUUm+ADNiYQAAl5QAWlrqAMPF2wCCd60ATTzNAGdPwQATc+sABjJbABKS9QCsruEAfn7VADFAQACWleUAAnV0ADea5ACIgPIAc3PfAGZk8QAVFLgAjITZAFllxAAMjpAAkpPGAFI6uwAZQ0QAZ1fXABcuSACnotkACFRaAFNT2wAIuusACI39AFpbqwAAqqwAZmbMAFFMTAB6d8YArq3SAAsuMQAApJkA0cfeAAqxqwC6t9YAN2FwAJV6gQBdXdgAhHThAFJSxQBDNucAV1FUAEQzuQB6e+AADaGrAFpOzwBkZL0AAI6PAD8+xwBIR8IAXl7iAMa88wAMMz0AAZyeAJGO2gAFvbQAcnDDAF9ezQCenesASj/FAFFR4wAJqKoAz8/bAAJ6egAWJ0UAg36wAFdX1gC7u9AAAKKiAGdrxwCDg9UAqKjcAAanmgAdQEAAYWLYAF1SUwBPTboAzcH+AAOBgQAQde8AV1fIAGNjuAAEt7kAAHFzAAGamgAFoKMAxbv3AJCQ1wCCbuYABrG0AAJqagCnpdoA0NDZAAajpAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATpA/IRNTQR5FoAAAAAAAADl0f25AVWpQGAoGnHcAAABrfDo3cSwMK0KEkSWaAAAABQtgZldPPkpEVjRoAAAAAHJ1fWF6bEMCkzUAAwAAAAA4JopcLksjWw94AAAAAAAAlR+CTE0cWDBaYiIAAAAAAA5SRiQVFDyMZUkoAAAAAACUXoZZSBEAewChMgcAAAAAi5+OaWQIAACSAIMWAAAAAG+AM4kNcGMtjZgvl20AAACbdh0gNl0AjyqeAAkaAAAAh4VnlhIEXycAKXOZgRcAAH4ZAIgbRwAAAABUeZ1RAAA9MRA7AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA/AAAABwAAAAcAAAAPAAAALwAAAD8AAAAfAAAAHwAAAo8AAANPAAAABwAAAicAAACDAAAjwwAAD/8AAP//AAA='
+
bgThreads=[] # needed to keep background threads alive
class InfoProvider:
=== modified file 'py/tests/__init__.py'
--- py/tests/__init__.py 2010-10-20 15:38:06 +0000
+++ py/tests/__init__.py 2010-11-06 11:33:09 +0000
@@ -4,7 +4,7 @@
import unittest,inspect
# add any new test suites to the list here, so that they are picked up by testAll
-allTests=['wrapper','omega','pbc','cohesive-chain']
+allTests=['wrapper','omega','pbc','clump','cohesive-chain']
# all yade modules (ugly...)
import yade.eudoxos,yade.export,yade.linterpolation,yade.log,yade.pack,yade.plot,yade.post2d,yade.timing,yade.utils,yade.ymport
=== added file 'py/tests/clump.py'
--- py/tests/clump.py 1970-01-01 00:00:00 +0000
+++ py/tests/clump.py 2010-11-06 11:33:09 +0000
@@ -0,0 +1,67 @@
+
+'''
+Various computations affected by the periodic boundary conditions.
+'''
+
+import unittest
+import random
+from yade.wrapper import *
+from miniEigen import *
+from yade._customConverters import *
+from yade import utils
+from yade import *
+from math import *
+
+class TestSimpleClump(unittest.TestCase):
+ "Test things on a simple clump composed of 2 spheres."
+ def setUp(self):
+ O.reset()
+ r1,r2,p0,p1=1,.5,Vector3.Zero,Vector3(0,0,3)
+ self.idC,(self.id1,self.id2)=O.bodies.appendClumped([
+ utils.sphere(p0,r1),
+ utils.sphere(p1,r2)
+ ])
+ def testConsistency(self):
+ "Clump: ids and flags consistency"
+ b1,b2,bC=[O.bodies[id] for id in (self.id1,self.id2,self.idC)]
+ self.assertEqual(b1.clumpId,bC.id)
+ self.assertEqual(b2.clumpId,bC.id)
+ self.assertEqual(bC.clumpId,bC.id)
+ self.assert_(bC.isClump)
+ self.assert_(b1.isClumpMember)
+ self.assert_(b2.isClumpMember)
+ self.assert_(not bC.bounded)
+ def testStaticProperties(self):
+ "Clump: mass, centroid, intertia"
+ b1,b2,bC=[O.bodies[id] for id in (self.id1,self.id2,self.idC)]
+ # mass
+ self.assertEqual(bC.state.mass,b1.state.mass+b2.state.mass)
+ # centroid
+ S=b1.state.mass*b1.state.pos+b2.state.mass*b2.state.pos
+ c=S/bC.state.mass
+ self.assertEqual(bC.state.pos,c);
+ # inertia
+ i1,i2=(8./15)*pi*b1.material.density*b1.shape.radius**5, (8./15)*pi*b2.material.density*b2.shape.radius**5 # inertia of spheres
+ iMax=i1+i2+b1.state.mass*(b1.state.pos-c).norm()**2+b2.state.mass*(b2.state.pos-c).norm()**2 # minimum principal inertia
+ iMin=i1+i2 # perpendicular to the
+ # the order of bC.state.inertia is arbitrary (though must match the orientation)
+ iC=list(bC.state.inertia); iC.sort()
+ self.assertAlmostEqual(iC[0],iMin)
+ self.assertAlmostEqual(iC[1],iMax)
+ self.assertAlmostEqual(iC[2],iMax)
+ # check orientation...?
+ #self.assertAlmostEqual
+ def testVelocity(self):
+ "Clump: velocities of member assigned by NewtonIntegrator"
+ s1,s2,sC=[O.bodies[id].state for id in (self.id1,self.id2,self.idC)]
+ O.dt=0
+ sC.vel=(1.,.2,.4)
+ sC.angVel=(0,.4,.1)
+ O.engines=[NewtonIntegrator()]; O.step() # update velocities
+ # linear velocities
+ self.assertEqual(s1.vel,sC.vel+sC.angVel.cross(s1.pos-sC.pos))
+ self.assertEqual(s2.vel,sC.vel+sC.angVel.cross(s2.pos-sC.pos))
+ # angular velocities
+ self.assertEqual(s1.angVel,sC.angVel);
+ self.assertEqual(s2.angVel,sC.angVel);
+
=== modified file 'py/tests/omega.py'
--- py/tests/omega.py 2010-10-29 10:12:44 +0000
+++ py/tests/omega.py 2010-11-06 11:33:09 +0000
@@ -20,7 +20,7 @@
class TestIO(unittest.TestCase):
def testSaveAllClasses(self):
- 'All classes can be saved and loaded with boost::serialization'
+ 'I/O: All classes can be saved and loaded with boost::serialization'
import yade.system
failed=set()
for c in yade.system.childClasses('Serializable'):
=== modified file 'py/wrapper/yadeWrapper.cpp'
--- py/wrapper/yadeWrapper.cpp 2010-11-02 12:02:13 +0000
+++ py/wrapper/yadeWrapper.cpp 2010-11-06 11:33:09 +0000
@@ -117,19 +117,25 @@
#endif
vector<Body::id_t> ret; FOREACH(shared_ptr<Body>& b, bb){ret.push_back(append(b));} return ret;
}
- python::tuple appendClump(vector<shared_ptr<Body> > bb){
- // update clump members
- vector<Body::id_t> ids(appendList(bb));
+ Body::id_t clump(vector<Body::id_t> ids){
// create and add clump itself
+ shared_ptr<Body> clumpBody=shared_ptr<Body>(new Body());
shared_ptr<Clump> clump=shared_ptr<Clump>(new Clump());
- shared_ptr<Body> clumpAsBody=static_pointer_cast<Body>(clump);
- clump->setDynamic(true);
- proxee->insert(clumpAsBody);
+ clumpBody->shape=clump;
+ clumpBody->setDynamic(true);
+ clumpBody->setBounded(false);
+ proxee->insert(clumpBody);
// add clump members to the clump
- FOREACH(Body::id_t id, ids) clump->add(id);
- // update clump
- clump->updateProperties(false);
- return python::make_tuple(clump->getId(),ids);
+ Scene* scene(Omega::instance().getScene().get());
+ FOREACH(Body::id_t id, ids) Clump::add(clumpBody,Body::byId(id,scene));
+ Clump::updateProperties(clumpBody,/*intersecting*/ false);
+ return clumpBody->getId();
+ }
+ python::tuple appendClump(vector<shared_ptr<Body> > bb){
+ // append constituent particles
+ vector<Body::id_t> ids(appendList(bb));
+ // clump them together (the clump fcn) and return
+ return python::make_tuple(clump(ids),ids);
}
vector<Body::id_t> replace(vector<shared_ptr<Body> > bb){proxee->clear(); return appendList(bb);}
long length(){return proxee->size();}
@@ -575,6 +581,7 @@
.def("append",&pyBodyContainer::append,"Append one Body instance, return its id.")
.def("append",&pyBodyContainer::appendList,"Append list of Body instance, return list of ids")
.def("appendClumped",&pyBodyContainer::appendClump,"Append given list of bodies as a clump (rigid aggregate); return list of ids.")
+ .def("clump",&pyBodyContainer::clump,"Clump given bodies together (creating a rigid aggregate); returns clump id.")
.def("clear", &pyBodyContainer::clear,"Remove all bodies (interactions not checked)")
.def("erase", &pyBodyContainer::erase,"Erase body with the given id; all interaction will be deleted by InteractionLoop in the next step.")
.def("replace",&pyBodyContainer::replace);
