yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01220
[svn] r1766 - in trunk: extra extra/clump gui gui/py pkg/dem pkg/dem/DataClass pkg/dem/DataClass/InteractionGeometry pkg/dem/Engine/DeusExMachina pkg/dem/Engine/EngineUnit pkg/dem/Engine/StandAloneEngine pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry scripts/test
Author: eudoxos
Date: 2009-05-05 07:15:26 +0200 (Tue, 05 May 2009)
New Revision: 1766
Added:
trunk/pkg/dem/DataClass/Clump.cpp
trunk/pkg/dem/DataClass/Clump.hpp
trunk/pkg/dem/Engine/EngineUnit/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.cpp
trunk/pkg/dem/Engine/EngineUnit/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.hpp
trunk/scripts/test/clump.py
Removed:
trunk/extra/clump/Clump.cpp
trunk/extra/clump/Clump.hpp
trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.cpp
trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.hpp
Modified:
trunk/extra/Brefcom.cpp
trunk/extra/BrefcomTestGen.cpp
trunk/extra/SConscript
trunk/extra/clump/Shop.cpp
trunk/gui/SConscript
trunk/gui/py/_eudoxos.cpp
trunk/gui/py/_utils.cpp
trunk/gui/py/yadeControl.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.hpp
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp
trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp
trunk/pkg/dem/SConscript
Log:
1. Add python wrapper for clumps, see scripts/test/clump.py on how to use it (O.bodies.appendClumped)
2. Fix clump support in NewtonsDampedLaw (damping is quite wrong, it seems...; Bruno, any ideas on that?)
3. Move clump to pkg-dem (even though it depends on extra/Shop)
4. Remove absolute shear computation from SpheresContactGeometry, as Dem3Dof works fine now and Brefcom will not use SCG anymore; adapt a few other classes to that; remove ElasticContactLaw2 that was using it.
5.
Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/extra/Brefcom.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -139,6 +139,7 @@
}
#endif
+
CREATE_LOGGER(ef2_Spheres_Brefcom_BrefcomLaw);
long BrefcomContact::cummBetaIter=0, BrefcomContact::cummBetaCount=0;
@@ -295,7 +296,7 @@
void GLDrawBrefcomContact::go(const shared_ptr<InteractionPhysics>& ip, const shared_ptr<Interaction>& i, const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool wireFrame){
const shared_ptr<BrefcomContact>& BC=static_pointer_cast<BrefcomContact>(ip);
- const shared_ptr<SpheresContactGeometry>& geom=YADE_PTR_CAST<SpheresContactGeometry>(i->interactionGeometry);
+ const shared_ptr<Dem3DofGeom>& geom=YADE_PTR_CAST<Dem3DofGeom>(i->interactionGeometry);
//Vector3r lineColor(BC->omega,1-BC->omega,0.0); /* damaged links red, undamaged green */
Vector3r lineColor=Shop::scalarOnColorScale(1.-BC->relResidualStrength);
@@ -329,11 +330,11 @@
glPopMatrix();
}
- const Vector3r& cp=static_pointer_cast<SpheresContactGeometry>(i->interactionGeometry)->contactPoint;
+ const Vector3r& cp=static_pointer_cast<Dem3DofGeom>(i->interactionGeometry)->contactPoint;
if(epsT){
Real maxShear=(BC->undamagedCohesion-BC->sigmaN*BC->tanFrictionAngle)/BC->G;
Real relShear=BC->epsT.Length()/maxShear;
- Real scale=.5*geom->d0;
+ Real scale=.5*geom->refLength;
Vector3r dirShear=BC->epsT; dirShear.Normalize();
if(epsTAxes){
GLUtils::GLDrawLine(cp-Vector3r(scale,0,0),cp+Vector3r(scale,0,0));
Modified: trunk/extra/BrefcomTestGen.cpp
===================================================================
--- trunk/extra/BrefcomTestGen.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/extra/BrefcomTestGen.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -17,7 +17,6 @@
#include<yade/pkg-common/MetaInteractingGeometry2AABB.hpp>
#include<yade/pkg-common/InteractionGeometryMetaEngine.hpp>
#include<yade/pkg-common/InteractionPhysicsMetaEngine.hpp>
-#include<yade/pkg-dem/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp>
#include<yade/pkg-common/PhysicalActionApplier.hpp>
#include<yade/pkg-common/PhysicalParametersMetaEngine.hpp>
#include<yade/pkg-common/NewtonsForceLaw.hpp>
@@ -28,6 +27,7 @@
#include<yade/pkg-common/ConstitutiveLawDispatcher.hpp>
#include<yade/pkg-dem/PositionOrientationRecorder.hpp>
#include<yade/pkg-dem/GlobalStiffnessTimeStepper.hpp>
+#include<yade/pkg-dem/Dem3DofGeom_SphereSphere.hpp>
#include<yade/extra/UniaxialStrainControlledTest.hpp>
@@ -51,9 +51,8 @@
rootBody->engines.push_back(collider);
shared_ptr<InteractionGeometryMetaEngine> igeomDispatcher(new InteractionGeometryMetaEngine);
- shared_ptr<InteractingSphere2InteractingSphere4SpheresContactGeometry> is2is4scg(new InteractingSphere2InteractingSphere4SpheresContactGeometry);
- is2is4scg->hasShear=true;
- igeomDispatcher->add(is2is4scg);
+ shared_ptr<ef2_Sphere_Sphere_Dem3DofGeom> ef2ssd3d(new ef2_Sphere_Sphere_Dem3DofGeom);
+ igeomDispatcher->add(ef2ssd3d);
rootBody->engines.push_back(igeomDispatcher);
shared_ptr<InteractionPhysicsMetaEngine> iphysDispatcher(new InteractionPhysicsMetaEngine);
@@ -83,11 +82,6 @@
shared_ptr<BrefcomDamageColorizer> dmg(new BrefcomDamageColorizer);
rootBody->engines.push_back(dmg);
- shared_ptr<PositionOrientationRecorder> rec(new PositionOrientationRecorder);
- rec->outputFile="/tmp/brefcom-test";
- rec->interval=500;
- rec->saveRgb=true;
- rootBody->engines.push_back(rec);
#if 0
shared_ptr<BrefcomStiffnessCounter> bsc(new BrefcomStiffnessCounter);
bsc->interval=100;
Modified: trunk/extra/SConscript
===================================================================
--- trunk/extra/SConscript 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/extra/SConscript 2009-05-05 05:15:26 UTC (rev 1766)
@@ -28,10 +28,6 @@
'yade-base',
]),
- env.SharedLibrary('Clump',
- ['clump/Clump.cpp'],
- LIBS=env['LIBS']+['Shop']),
-
env.SharedLibrary('Tetra',['tetra/Tetra.cpp'],LIBS=env['LIBS']+['Tetrahedron',
'ParticleParameters',
'RigidBodyParameters',
@@ -45,7 +41,7 @@
env.SharedLibrary('Brefcom',['Brefcom.cpp'],CXXFLAGS=env['CXXFLAGS']+brefcomInclude,LIBS=env['LIBS']+['Shop','InteractingSphere2InteractingSphere4SpheresContactGeometry','DemXDofGeom']),
- env.SharedLibrary('BrefcomTestGen',['BrefcomTestGen.cpp'],LIBS=env['LIBS']+['Shop','UniaxialStrainControlledTest','PositionOrientationRecorder','Brefcom','ConstitutiveLawDispatcher']),
+ env.SharedLibrary('BrefcomTestGen',['BrefcomTestGen.cpp'],LIBS=env['LIBS']+['Shop','UniaxialStrainControlledTest','PositionOrientationRecorder','Brefcom','ConstitutiveLawDispatcher','Dem3DofGeom_SphereSphere']),
env.SharedLibrary('SimpleScene',['SimpleScene.cpp'],LIBS=env['LIBS']+['Shop','SimpleElasticRelationships']),
Deleted: trunk/extra/clump/Clump.cpp
===================================================================
--- trunk/extra/clump/Clump.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/extra/clump/Clump.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -1,428 +0,0 @@
-// (c) 2007 Vaclav Smilauer <eudoxos@xxxxxxxx>
-
-#include"Clump.hpp"
-#include<algorithm>
-
-YADE_PLUGIN("Clump","ClumpMemberMover","ClumpTestGen");
-
-CREATE_LOGGER(Clump);
-CREATE_LOGGER(ClumpMemberMover);
-CREATE_LOGGER(ClumpTestGen);
-
-/**************************************************************************************
- ************************************* ClumpMemberMover ******************************
- **************************************************************************************/
-
-// Constructor must be in the .cpp file (?)
-ClumpMemberMover::ClumpMemberMover(){/*createIndex();*/ }
-
-/*! 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::applyCondition(MetaBody* rootBody){
- for(BodyContainer::iterator I=rootBody->bodies->begin(); I!=rootBody->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?!
-}
-
-/**************************************************************************************
- ******************************************** Clump ***********************************
- **************************************************************************************/
-
-/*! Create zero'ed RigidBodyParameters; they should not be manipulated directly, since they are all calculated in Clump::update.
- * @todo do we need to do the same for GeomtricalModel, InteractingGeometry and BoundingVolume? They will never be used. Sort that out for sure.
- * @bug setting Clump::isDynamic in constructor is not enough (as if it were modified somewhere), must be set explicitly by the user after construction (why?)
- */
-Clump::Clump(): Body(){
- isDynamic=true;
- physicalParameters=shared_ptr<RigidBodyParameters>(new RigidBodyParameters);
-
- // these will not be defined for the moment...
-#if 0
- boundingVolume=shared_ptr<AABB>(new AABB);
- boundingVolume->diffuseColor=Vector3r(1,0,0);
-
- interactingGeometry=shared_ptr<InteractingGeometry>(new InteractingGeometry);
- interactingGeometry->diffuseColor=Vector3r(0,0,0);
-
- geometricalModel=shared_ptr<GeometricalModel>(new GeometricalModel);
- geometricalModel->diffuseColor=Vector3r(0,0,0); geometricalModel->wire=false; geometricalModel->visible=false; geometricalModel->shadowCaster=false;
-#endif
-
-}
-
-/*! @pre Body must be dynamic.
- * @pre Body must not be part or this clump already.
- * @pre Body must have valid (non-NULL) Body::physicalParameters
- * @todo se3 calculation is not tested yet
- */
-void Clump::add(body_id_t subId){
- shared_ptr<Body> subBody=Body::byId(subId);
-
- // preconditions
- assert(subBody->isDynamic);
- assert(physicalParameters);
- assert(members.count(subId)==0);
-
- // begin actual setup
- subBody->clumpId=getId();
- subBody->isDynamic=false;
- // for now, push just unitialized se3; will be updated by updateProperties
- members[subId]=Se3r();
-
- clumpId=getId(); // identifies a clump
-
- 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);
- subBody->clumpId=Body::ID_NONE;
- subBody->isDynamic=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));
- //LOG_TRACE("Old #"<<I->first<<"position: "<<subRBP->se3.position);
- subRBP->se3.position=mySe3.position+mySe3.orientation*I->second.position;
- subRBP->se3.orientation=mySe3.orientation*I->second.orientation;
- //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?!
- subRBP->velocity=myRBP->velocity+myRBP->angularVelocity.Cross(I->second.position);
- subRBP->angularVelocity=myRBP->angularVelocity;
- }
- /* @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().getCurrentIteration()%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.Length(),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
-
- clumpRBP->acceleration=Vector3r(0,0,0);
- clumpRBP->angularAcceleration=Vector3r(0,0,0);
-
-}
-
-/*! 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:
- -# 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]
- - M+=m_i
- - S+=m_i*x_i (local static moments are zero (centroid)
- - get inertia tensor of subBody in world coordinates, by rotating the principal (local) tensor against subBody->se3->orientation; then translate it to world origin (parallel axes theorem), then I+=I_i_world
- -# clumpPos=S/M
- -# translate aggregate's inertia tensor; parallel axes on I (R=clumpPos): 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]
- -# eigen decomposition of I, get principal inertia and rotation matrix of the clump
- -# 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){
- 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(true /* fill with zeros */ ), Ic(true); // 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){
- LOG_DEBUG("Clump of size one will be treated specially.")
- memberMap::iterator I=members.begin();
- shared_ptr<Body> subBody=Body::byId(I->first);
- const shared_ptr<RigidBodyParameters>& subRBP(YADE_PTR_CAST<RigidBodyParameters>(subBody->physicalParameters));
- // se3 of the clump as whole is the same as the member's se3
- mySe3.position=subRBP->se3.position;
- mySe3.orientation=subRBP->se3.orientation;
- // relative member's se3 is identity
- I->second.position=Vector3r::ZERO; I->second.orientation=Quaternionr::IDENTITY;
- clumpRBP->inertia=subRBP->inertia;
- clumpRBP->mass=subRBP->mass;
- clumpRBP->velocity=Vector3r::ZERO;
- clumpRBP->angularVelocity=Vector3r::ZERO;
- return;
- }
-
- if(intersecting){
- LOG_WARN("Self-intersecting clumps not yet implemented, intersections will be ignored.");
- 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));
- M+=subRBP->mass;
- Sg+=subRBP->mass*subRBP->se3.position;
- // transform from local to global coords
- // FIXME: verify this!
- Quaternionr subRBP_orientation_conjugate=subRBP->se3.orientation.Conjugate();
- Matrix3r Imatrix(subRBP->inertia[0],subRBP->inertia[1],subRBP->inertia[2]);
- // TRWM3MAT(Imatrix); TRWM3QUAT(subRBP_orientation_conjugate);
- Ig+=Clump::inertiaTensorTranslate(Clump::inertiaTensorRotate(Imatrix,subRBP_orientation_conjugate),subRBP->mass,-1.*subRBP->se3.position);
-
- //TRWM3MAT(Clump::inertiaTensorRotate(Matrix3r(subRBP->inertia),subRBP_orientation_conjugate));
- }
- }
- TRVAR1(M);
- TRWM3MAT(Ig);
- TRWM3VEC(Sg);
-
- mySe3.position=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 */, mySe3.position); // inertia at clump's centroid but with world orientation
- TRWM3MAT(Ic_orientG);
-
- Matrix3r R_g2c(true); //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);
- Ic_orientG.EigenDecomposition(R_g2c,Ic);
- /*! @bug: eigendecomposition is wrong. see http://article.gmane.org/gmane.science.physics.yade.devel/99 for message. */
- // 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
- mySe3.orientation.FromRotationMatrix(R_g2c);
- // now Ic is diagonal
- clumpRBP->inertia=Vector3r(Ic(0,0),Ic(1,1),Ic(2,2));
- clumpRBP->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);
- const shared_ptr<RigidBodyParameters>& subRBP(YADE_PTR_CAST<RigidBodyParameters>(subBody->physicalParameters));
- clumpRBP->inertia=subRBP->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)
- mySe3.orientation.FromRotationMatrix(Matrix3r(1,0,0,0,1,0,0,0,1));
- }
- #endif
- TRWM3VEC(clumpRBP->inertia);
-
- // TODO: these might be calculated from members... but complicated... - someone needs that?!
- clumpRBP->velocity=Vector3r(0,0,0);
- clumpRBP->angularVelocity=Vector3r(0,0,0);
-
- // 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++){
- // 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));
- I->second.orientation=mySe3.orientation.Conjugate()*subRBP->se3.orientation;
- I->second.position=mySe3.orientation.Conjugate()*(subRBP->se3.position-mySe3.position);
- }
-
-}
-
-/*! @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).
- * @param m mass of the body; if positive, translation is away from the centroid; if negative, towards centroid.
- * @param off offset of the new origin from the original origin
- * @return inertia tensor in the new coordinate system; the matrix is symmetric.
- */
-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);
- // 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] */
- I2+=m*Matrix3r(/* 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(1,0)=I2(0,1); I2(2,0)=I2(0,2); I2(2,1)=I2(1,2);
- //TRWM3MAT(I2);
- return I2;
-}
-
-/*! @brief Recalculate body's inertia tensor in rotated coordinates.
- *
- * @param I inertia tensor in old coordinates
- * @param T rotation matrix from old to new coordinates
- * @return inertia tensor in new coordinates
- */
-Matrix3r Clump::inertiaTensorRotate(const Matrix3r& I,const Matrix3r& T){
- /* [http://www.kwon3d.com/theory/moi/triten.html] */
- //TRWM3MAT(I); TRWM3MAT(T);
- return T.Transpose()*I*T;
-}
-
-/*! @brief Recalculate body's inertia tensor in rotated coordinates.
- *
- * @param I inertia tensor in old coordinates
- * @param rot quaternion that describes rotation from old to new coordinates
- * @return inertia tensor in new coordinates
- */
-Matrix3r Clump::inertiaTensorRotate(const Matrix3r& I, const Quaternionr& rot){
- Matrix3r T;
- rot.ToRotationMatrix(T);
- return inertiaTensorRotate(I,T);
-}
-
-
-
-/**************************************************************************************
- ********************* ClumpTestGen ***************************************************
- **************************************************************************************/
-
-#include<yade/core/MetaBody.hpp>
-#include<yade/extra/Shop.hpp>
-
-bool ClumpTestGen::generate()
-{
- //Shop::setDefault("param_pythonRunExpr",string("if S.i%50==0 and S.i<1000 and S.i>500:\n\tprint S.i,len(S.sel),B[1].x, B[1].E"));
-
- rootBody=Shop::rootBody();
- Shop::rootBodyActors(rootBody);
- // clumps do not need to subscribe currently (that will most likely change, though)
- rootBody->engines.push_back(shared_ptr<ClumpMemberMover>(new ClumpMemberMover));
-
-
- shared_ptr<MetaBody> oldRootBody=Omega::instance().getRootBody();
- Omega::instance().setRootBody(rootBody);
-
- shared_ptr<Body> ground=Shop::box(Vector3r(0,0,-1),Vector3r(3,3,.2));
- ground->isDynamic=false;
- // revert random colors for this single case...
- ground->geometricalModel->diffuseColor=Vector3r(.6,.6,.6);
- ground->interactingGeometry->diffuseColor=Vector3r(.3,.3,.3);
- rootBody->bodies->insert(ground);
-
- vector<Vector3r> relPos; vector<Real> radii; Vector3r clumpPos;
-
- // standalone (non-clump!) sphere as well
- shared_ptr<Body> sphere=Shop::sphere(Vector3r(0,0,0),.5);
- rootBody->bodies->insert(sphere);
- // one-sphere clump
- clumpPos=Vector3r(-2,0,0);
- relPos.push_back(Vector3r(0,0,0)); radii.push_back(.5);
- createOneClump(rootBody,clumpPos,relPos,radii);
- relPos.clear(); radii.clear();
- // two-sphere clump
- clumpPos=Vector3r(2,0,0);
- relPos.push_back(Vector3r(0,-.5,.5)); radii.push_back(.5);
- relPos.push_back(Vector3r(0,.5,0)); radii.push_back(.5);
- createOneClump(rootBody,clumpPos,relPos,radii);
- relPos.clear(); radii.clear();
- // three-sphere slump
- clumpPos=Vector3r(0,2,0);
- relPos.push_back(Vector3r(0,-.5,.5)); radii.push_back(.5);
- relPos.push_back(Vector3r(0,.5,0)); radii.push_back(.5);
- relPos.push_back(Vector3r(.5,0,0)); radii.push_back(.5);
- createOneClump(rootBody,clumpPos,relPos,radii);
- relPos.clear(); radii.clear();
- // four-sphere slump
- clumpPos=Vector3r(0,-2,0);
- relPos.push_back(Vector3r(0,0,0)); radii.push_back(.5);
- relPos.push_back(Vector3r(.5,0,0)); radii.push_back(.5);
- relPos.push_back(Vector3r(0,.5,0)); radii.push_back(.5);
- relPos.push_back(Vector3r(0,0,.5)); radii.push_back(.5);
- createOneClump(rootBody,clumpPos,relPos,radii);
- relPos.clear(); radii.clear();
-
- // restore Omega
- Omega::instance().setRootBody(oldRootBody);
-
- message="OK";
- return true;
-}
-
-/*! \brief Generate clump of spheres, the result will be inserted into rootBody.
- *
- * To create a clump, first the clump itself needs to be instantiated \em and inserted into rootBody (this will assign an Body::id).
- * In order for this to work, Omega::roootBody must have been assigned; within generators, use Omega::setRootBody for this.
- *
- * The body to add to clump must have been also created and added to the rootBody (so that it has id, again).
- *
- * Finally, call Clump::updateProperties to get physical properties physically right (inertia, position, orientation, mass, ...).
- *
- * @param clumpPos Center of the clump (not necessarily centroid); serves merely as reference for sphere positions.
- * @param relPos Relative positions of individual spheres' centers.
- * @param radii Radii of composing spheres. Must have the same length as relPos.
- */
-void ClumpTestGen::createOneClump(shared_ptr<MetaBody>& rootBody, Vector3r clumpPos, vector<Vector3r> relPos, vector<Real> radii)
-{
- assert(relPos.size()==radii.size());
-
- // empty clump
- shared_ptr<Clump> clump=shared_ptr<Clump>(new Clump());
- shared_ptr<Body> clumpAsBody(static_pointer_cast<Body>(clump));
- rootBody->bodies->insert(clumpAsBody);
-
- clump->isDynamic=true;
- // if subscribedBodies work some day: clumpMover->subscribedBodies.push_back(clump->getId());
-
- for(size_t i=0; i<relPos.size(); i++){
- shared_ptr<Body> sphere=Shop::sphere(clumpPos+relPos[i],radii[i]);
- body_id_t lastId=(body_id_t)rootBody->bodies->insert(sphere);
- clump->add(lastId);
- LOG_TRACE("Generated clumped sphere #"<<lastId);
- }
- clump->updateProperties(false);
-}
-
Deleted: trunk/extra/clump/Clump.hpp
===================================================================
--- trunk/extra/clump/Clump.hpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/extra/clump/Clump.hpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -1,148 +0,0 @@
-// (c) 2007 Vaclav Smilauer <eudoxos@xxxxxxxx>
-
-/* To compile this class, you MUST
- 1. use recent version of scons
- 2. say extraModules=yade-extra/clump (this will cause SConscript in this directory to be processed)
- Any further documentation is in doxygen comments.
-*/
-
-#pragma once
-
-#include<vector>
-#include<map>
-#include<yade/core/Body.hpp>
-#include<yade/core/MetaBody.hpp>
-#include<yade/core/FileGenerator.hpp>
-#include<yade/core/DeusExMachina.hpp>
-#include<yade/lib-factory/Factorable.hpp>
-#include<yade/pkg-common/PhysicalParametersEngineUnit.hpp>
-#include<yade/pkg-common/RigidBodyParameters.hpp>
-#include<yade/pkg-common/AABB.hpp>
-#include<yade/lib-base/Logging.hpp>
-#include<yade/lib-base/yadeWm3Extra.hpp>
-
-
-/*! Body representing clump (rigid aggregate) composed by other existing bodies.
-
- Clump is one of bodies that reside in rootBody->bodies.
- When an existing body is added to ::Clump, it's ::Body::isDynamic flag is set to false
- (it is still subscribed to all its engines, to make it possible to remove it from the clump again).
- All forces acting on Clump::members are made to act on the clump itself, which will ensure that they
- influence all Clump::members as if the clump were a rigid particle.
-
- What are clump requirements so that they function?
- -# Given any body, tell
- - if it is a clump member: Body::isClumpMember()
- - if it is a clump: Body:: isClump(). (Correct result is assured at each call to Clump::add).
- (we could use RTTI instead? Would that be more reliable?)
- - if it is a standalone Body: Body::isStandalone()
- - what is it's clump id (Body::clumpId)
- -# given the root body, tell
- - what clumps it contains (enumerate all bodies and filter clumps, see above)
- -# given a clump, tell
- - what bodies it contains (keys of ::Clump::members)
- - what are se3 of these bodies (values of ::Clump::members)
- -# add/delete bodies from/to clump (::Clump::add, ::Clump::del)
- - This includes saving se3 of the subBody: it \em must be in clump's local coordinates so that it is constant. The transformation from global to local is given by clump's se3 at the moment of addition. Clump's se3 is initially (origin,identity)
- -# Update clump's physical properties (Clump::updateProperties)
- - This \em must reposition members so that they have the same se3 globally
- -# Apply forces acting on members to the clump instead (done in NewtonsForceLaw, NewtonsMomentumLaw) - uses world coordinates to calculate effect on the clump's centroid
- -# Integrate position and orientation of the clump
- - LeapFrogPositionIntegrator and LeapFrogOrientationIntegrator move clump as whole
- - clump members are skipped, since they have Body::isDynamic==false.
- - ClumpMemberMover is an engine that updates positions of the clump memebers in each timestep (calls Clump::moveSubBodies internally)
-
- Some more information can be found http://beta.arcig.cz/~eudoxos/phd/index.cgi/YaDe/HighLevelClumps
-
- For an example how to generate a clump, see ClumpTestGen::createOneClump.
-
- @todo GravityEngine should be applied to members, not to clump as such?! Still not sure. Perhaps Clumps should have mass and inertia set to zeros so that engines unaware of clumps do not act on it. It would have some private mass and insertia that would be used in NewtonsForceLaw etc for clumps specially...
-
- @note PersistentSAPCollider bypass Clumps explicitly. This no longer depends on the absence of boundingVolume.
- @note Clump relies on its id being assigned (as well as id of its components); therefore, only bodies that have already been inserted to the container may be added to Clump which has been itself already added to the container.
-
- */
-
-class Clump: public Body {
- //! mapping of body IDs to their relative positions; replaces members and subSe3s;
- public:
- typedef std::map<body_id_t,Se3r> memberMap;
- memberMap members;
-
- Clump();
- virtual ~Clump(){};
- //! \brief add Body to the Clump
- void add(body_id_t);
- //! \brief remove Body from the Clump
- void del(body_id_t);
- //! Recalculate physical properties of Clump.
- void updateProperties(bool intersecting);
- //! Calculate positions and orientations of members based on my own Se3.
- void moveMembers();
- //! 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...
- //! 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.
- static Matrix3r inertiaTensorRotate(const Matrix3r& I, const Matrix3r& T);
- //! Recalculate body's inertia tensor in rotated coordinates.
- static Matrix3r inertiaTensorRotate(const Matrix3r& I, const Quaternionr& rot);
-
- void registerAttributes(){Body::registerAttributes(); REGISTER_ATTRIBUTE(members);}
- REGISTER_CLASS_NAME(Clump);
- REGISTER_BASE_CLASS_NAME(Body);
- // REGISTER_CLASS_INDEX(Clump,Body);
- DECLARE_LOGGER;
-};
-
-REGISTER_SERIALIZABLE(Clump);
-
-/*! Update ::Clump::members positions so that the Clump behaves as a rigid body.
- *
- *
-*/
-class ClumpMemberMover: public DeusExMachina {
- public:
- //! Interates over rootBody->bodies and calls Clump::moveSubBodies() for clumps.
- virtual void applyCondition(MetaBody* rootBody);
- ClumpMemberMover();
- virtual ~ClumpMemberMover(){};
-
- REGISTER_CLASS_NAME(ClumpMemberMover);
- REGISTER_BASE_CLASS_NAME(DeusExMachina);
- // REGISTER_CLASS_INDEX(ClumpMemberMover,PhysicalParametersEngineUnit);
- DECLARE_LOGGER;
-};
-
-REGISTER_SERIALIZABLE(ClumpMemberMover);
-
-/*! \brief Test some basic clump functionality; show how to use clumps as well. */
-class ClumpTestGen : public FileGenerator {
- DECLARE_LOGGER;
- private :
- //Vector3r nbTetrahedrons,groundSize,gravity;
- //Real minSize,density,maxSize,dampingForce,disorder,dampingMomentum,youngModulus;
- //int timeStepUpdateInterval;
- //bool rotationBlocked;
- //void createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents);
- void createOneClump(shared_ptr<MetaBody>& rootBody, Vector3r clumpPos, vector<Vector3r> relPos, vector<Real> radii);
- //shared_ptr<Body> createOneSphere(Vector3r position, Real radius);
- //void createActors(shared_ptr<MetaBody>& rootBody);
- //void positionRootBody(shared_ptr<MetaBody>& rootBody);
- //void calculatePropertiesAndReposition(const shared_ptr<SlumShape>& slum, shared_ptr<ElasticBodyParameters>& rbp, Real density);
- //void makeTet(shared_ptr<Tetrahedron>& tet, Real radius);
- shared_ptr<ClumpMemberMover> clumpMover;
- public :
- ClumpTestGen (){};
- ~ClumpTestGen (){};
- bool generate();
- protected :
- virtual void postProcessAttributes(bool deserializing){};
- void registerAttributes(){};
- REGISTER_CLASS_NAME(ClumpTestGen);
- REGISTER_BASE_CLASS_NAME(FileGenerator);
-};
-
-REGISTER_SERIALIZABLE(ClumpTestGen);
-
Modified: trunk/extra/clump/Shop.cpp
===================================================================
--- trunk/extra/clump/Shop.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/extra/clump/Shop.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -66,10 +66,7 @@
#include<yade/extra/Tetra.hpp>
-//#include<yade/extra/Clump.hpp>
-//#include <yade/pkg-dem/BodyMacroParameters.hpp>
-
#define _SPEC_CAST(orig,cast) template<> void Shop::setDefault<orig>(string key, orig val){setDefault(key,cast(val));}
_SPEC_CAST(const char*,string);
_SPEC_CAST(char*,string);
Modified: trunk/gui/SConscript
===================================================================
--- trunk/gui/SConscript 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/gui/SConscript 2009-05-05 05:15:26 UTC (rev 1766)
@@ -69,7 +69,8 @@
'ConstitutiveLawDispatcher',
'InteractionDispatchers',
'STLImporter',
- 'ParallelEngine'
+ 'ParallelEngine',
+ 'Clump'
],
),
env.SharedLibrary('_utils',['py/_utils.cpp'],SHLIBPREFIX='',
Modified: trunk/gui/py/_eudoxos.cpp
===================================================================
--- trunk/gui/py/_eudoxos.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/gui/py/_eudoxos.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -8,7 +8,6 @@
#endif
-
# if 0
Real elasticEnergyDensityInAABB(python::tuple AABB){
Vector3r bbMin=tuple2vec(python::extract<python::tuple>(AABB[0])()), bbMax=tuple2vec(python::extract<python::tuple>(AABB[1])()); Vector3r box=bbMax-bbMin;
Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/gui/py/_utils.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -5,6 +5,7 @@
#include<yade/core/Omega.hpp>
#include<yade/pkg-common/Sphere.hpp>
#include<yade/pkg-dem/SpheresContactGeometry.hpp>
+#include<yade/pkg-dem/DemXDofGeom.hpp>
#include<yade/pkg-dem/SimpleViscoelasticBodyParameters.hpp>
#include<yade/pkg-common/NormalShearInteractions.hpp>
#include<cmath>
@@ -81,7 +82,7 @@
FOREACH(const shared_ptr<Interaction>&i, *rb->transientInteractions){
if(!i->interactionPhysics) continue;
shared_ptr<NormalShearInteraction> bc=dynamic_pointer_cast<NormalShearInteraction>(i->interactionPhysics); if(!bc) continue;
- shared_ptr<SpheresContactGeometry> geom=dynamic_pointer_cast<SpheresContactGeometry>(i->interactionGeometry); if(!bc){LOG_ERROR("NormalShearInteraction contact doesn't have SpheresContactGeomety associated?!"); continue;}
+ shared_ptr<Dem3DofGeom> geom=dynamic_pointer_cast<Dem3DofGeom>(i->interactionGeometry); if(!bc){LOG_ERROR("NormalShearInteraction contact doesn't have SpheresContactGeomety associated?!"); continue;}
const shared_ptr<Body>& b1=Body::byId(i->getId1(),rb), b2=Body::byId(i->getId2(),rb);
bool isIn1=isInBB(b1->physicalParameters->se3.position,bbMin,bbMax), isIn2=isInBB(b2->physicalParameters->se3.position,bbMin,bbMax);
if(!isIn1 && !isIn2) continue;
@@ -96,7 +97,7 @@
//cerr<<"Interaction crosses AABB boundary, weight is "<<weight<<endl;
//LOG_DEBUG("Interaction crosses AABB boundary, weight is "<<weight);
} else {assert(isIn1 && isIn2); /* cerr<<"Interaction inside, weight is "<<weight<<endl;*/ /*LOG_DEBUG("Interaction inside, weight is "<<weight);*/}
- E+=geom->d0*weight*(.5*bc->kn*pow(geom->epsN(),2)+.5*bc->ks*pow(geom->epsT().Length(),2));
+ E+=geom->refLength*weight*(.5*bc->kn*pow(geom->strainN(),2)+.5*bc->ks*pow(geom->strainT().Length(),2));
}
return E;
}
Modified: trunk/gui/py/yadeControl.cpp
===================================================================
--- trunk/gui/py/yadeControl.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/gui/py/yadeControl.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -62,6 +62,7 @@
#include<yade/pkg-common/ConstitutiveLaw.hpp>
#include<yade/extra/Shop.hpp>
+#include<yade/pkg-dem/Clump.hpp>
using namespace boost;
using namespace std;
@@ -338,6 +339,16 @@
}
body_id_t insert(pyBody b){return proxee->insert(b.proxee);}
python::list insertList(python::list bb){python::list ret; for(int i=0; i<len(bb); i++){ret.append(insert(python::extract<pyBody>(bb[i])()));} return ret;}
+ python::tuple insertClump(python::list bb){/*clump: first add constitutents, then add clump, then add constitutents to the clump, then update clump props*/
+ python::list ids=insertList(bb);
+ shared_ptr<Clump> clump=shared_ptr<Clump>(new Clump());
+ shared_ptr<Body> clumpAsBody=static_pointer_cast<Body>(clump);
+ clump->isDynamic=true;
+ proxee->insert(clumpAsBody);
+ for(int i=0; i<len(ids); i++){clump->add(python::extract<body_id_t>(ids[i])());}
+ clump->updateProperties(false);
+ return python::make_tuple(clump->getId(),ids);
+ }
python::list replace(python::list bb){proxee->clear(); return insertList(bb);}
long length(){return proxee->size();}
void clear(){proxee->clear();}
@@ -702,6 +713,7 @@
.def("__len__",&pyBodyContainer::length)
.def("append",&pyBodyContainer::insert)
.def("append",&pyBodyContainer::insertList)
+ .def("appendClumped",&pyBodyContainer::insertClump)
.def("clear", &pyBodyContainer::clear)
.def("replace",&pyBodyContainer::replace);
boost::python::class_<pyInteractionContainer>("InteractionContainer",python::init<pyInteractionContainer&>())
Copied: trunk/pkg/dem/DataClass/Clump.cpp (from rev 1763, trunk/extra/clump/Clump.cpp)
===================================================================
--- trunk/extra/clump/Clump.cpp 2009-05-03 07:16:04 UTC (rev 1763)
+++ trunk/pkg/dem/DataClass/Clump.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -0,0 +1,427 @@
+// (c) 2007,2009 Vaclav Smilauer <eudoxos@xxxxxxxx>
+
+#include"Clump.hpp"
+#include<algorithm>
+
+YADE_PLUGIN("Clump","ClumpMemberMover","ClumpTestGen");
+
+CREATE_LOGGER(Clump);
+CREATE_LOGGER(ClumpMemberMover);
+CREATE_LOGGER(ClumpTestGen);
+
+/**************************************************************************************
+ ************************************* ClumpMemberMover ******************************
+ **************************************************************************************/
+
+// Constructor must be in the .cpp file (?)
+ClumpMemberMover::ClumpMemberMover(){/*createIndex();*/ }
+
+/*! 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::applyCondition(MetaBody* rootBody){
+ for(BodyContainer::iterator I=rootBody->bodies->begin(); I!=rootBody->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?!
+}
+
+/**************************************************************************************
+ ******************************************** Clump ***********************************
+ **************************************************************************************/
+
+/*! Create zero'ed RigidBodyParameters; they should not be manipulated directly, since they are all calculated in Clump::update.
+ * @todo do we need to do the same for GeomtricalModel, InteractingGeometry and BoundingVolume? They will never be used. Sort that out for sure.
+ * @bug setting Clump::isDynamic in constructor is not enough (as if it were modified somewhere), must be set explicitly by the user after construction (why?)
+ */
+Clump::Clump(): Body(){
+ isDynamic=true;
+ physicalParameters=shared_ptr<RigidBodyParameters>(new RigidBodyParameters);
+
+ // these will not be defined for the moment...
+#if 0
+ boundingVolume=shared_ptr<AABB>(new AABB);
+ boundingVolume->diffuseColor=Vector3r(1,0,0);
+
+ interactingGeometry=shared_ptr<InteractingGeometry>(new InteractingGeometry);
+ interactingGeometry->diffuseColor=Vector3r(0,0,0);
+
+ geometricalModel=shared_ptr<GeometricalModel>(new GeometricalModel);
+ geometricalModel->diffuseColor=Vector3r(0,0,0); geometricalModel->wire=false; geometricalModel->visible=false; geometricalModel->shadowCaster=false;
+#endif
+
+}
+
+/*! @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 NewtonsDampedLaw)
+ */
+void Clump::add(body_id_t subId){
+ shared_ptr<Body> subBody=Body::byId(subId);
+
+ // preconditions
+ assert(subBody->isDynamic);
+ assert(physicalParameters);
+ assert(members.count(subId)==0);
+ assert(subId<getId());
+
+ // begin actual setup
+ subBody->clumpId=getId();
+ subBody->isDynamic=false;
+ // for now, push just unitialized se3; will be updated by updateProperties
+ members[subId]=Se3r();
+
+ clumpId=getId(); // identifies a clump
+
+ 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);
+ subBody->clumpId=Body::ID_NONE;
+ subBody->isDynamic=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));
+ //LOG_TRACE("Old #"<<I->first<<"position: "<<subRBP->se3.position);
+ subRBP->se3.position=mySe3.position+mySe3.orientation*I->second.position;
+ subRBP->se3.orientation=mySe3.orientation*I->second.orientation;
+ //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?!
+ subRBP->velocity=myRBP->velocity+myRBP->angularVelocity.Cross(I->second.position);
+ subRBP->angularVelocity=myRBP->angularVelocity;
+ }
+ /* @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().getCurrentIteration()%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.Length(),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
+
+ clumpRBP->acceleration=Vector3r(0,0,0);
+ clumpRBP->angularAcceleration=Vector3r(0,0,0);
+
+}
+
+/*! 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:
+ -# 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]
+ - M+=m_i
+ - S+=m_i*x_i (local static moments are zero (centroid)
+ - get inertia tensor of subBody in world coordinates, by rotating the principal (local) tensor against subBody->se3->orientation; then translate it to world origin (parallel axes theorem), then I+=I_i_world
+ -# clumpPos=S/M
+ -# translate aggregate's inertia tensor; parallel axes on I (R=clumpPos): 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]
+ -# eigen decomposition of I, get principal inertia and rotation matrix of the clump
+ -# 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){
+ 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(true /* fill with zeros */ ), Ic(true); // 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){
+ LOG_DEBUG("Clump of size one will be treated specially.")
+ memberMap::iterator I=members.begin();
+ shared_ptr<Body> subBody=Body::byId(I->first);
+ const shared_ptr<RigidBodyParameters>& subRBP(YADE_PTR_CAST<RigidBodyParameters>(subBody->physicalParameters));
+ // se3 of the clump as whole is the same as the member's se3
+ mySe3.position=subRBP->se3.position;
+ mySe3.orientation=subRBP->se3.orientation;
+ // relative member's se3 is identity
+ I->second.position=Vector3r::ZERO; I->second.orientation=Quaternionr::IDENTITY;
+ clumpRBP->inertia=subRBP->inertia;
+ clumpRBP->mass=subRBP->mass;
+ clumpRBP->velocity=Vector3r::ZERO;
+ clumpRBP->angularVelocity=Vector3r::ZERO;
+ return;
+ }
+
+ if(intersecting){
+ LOG_WARN("Self-intersecting clumps not yet implemented, intersections will be ignored.");
+ 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));
+ M+=subRBP->mass;
+ Sg+=subRBP->mass*subRBP->se3.position;
+ // transform from local to global coords
+ // FIXME: verify this!
+ Quaternionr subRBP_orientation_conjugate=subRBP->se3.orientation.Conjugate();
+ Matrix3r Imatrix(subRBP->inertia[0],subRBP->inertia[1],subRBP->inertia[2]);
+ // TRWM3MAT(Imatrix); TRWM3QUAT(subRBP_orientation_conjugate);
+ Ig+=Clump::inertiaTensorTranslate(Clump::inertiaTensorRotate(Imatrix,subRBP_orientation_conjugate),subRBP->mass,-1.*subRBP->se3.position);
+
+ //TRWM3MAT(Clump::inertiaTensorRotate(Matrix3r(subRBP->inertia),subRBP_orientation_conjugate));
+ }
+ }
+ TRVAR1(M);
+ TRWM3MAT(Ig);
+ TRWM3VEC(Sg);
+
+ mySe3.position=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 */, mySe3.position); // inertia at clump's centroid but with world orientation
+ TRWM3MAT(Ic_orientG);
+
+ Matrix3r R_g2c(true); //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);
+ Ic_orientG.EigenDecomposition(R_g2c,Ic);
+ /*! @bug: eigendecomposition is wrong. see http://article.gmane.org/gmane.science.physics.yade.devel/99 for message. */
+ // 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
+ mySe3.orientation.FromRotationMatrix(R_g2c);
+ // now Ic is diagonal
+ clumpRBP->inertia=Vector3r(Ic(0,0),Ic(1,1),Ic(2,2));
+ clumpRBP->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);
+ const shared_ptr<RigidBodyParameters>& subRBP(YADE_PTR_CAST<RigidBodyParameters>(subBody->physicalParameters));
+ clumpRBP->inertia=subRBP->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)
+ mySe3.orientation.FromRotationMatrix(Matrix3r(1,0,0,0,1,0,0,0,1));
+ }
+ #endif
+ TRWM3VEC(clumpRBP->inertia);
+
+ // TODO: these might be calculated from members... but complicated... - someone needs that?!
+ clumpRBP->velocity=Vector3r(0,0,0);
+ clumpRBP->angularVelocity=Vector3r(0,0,0);
+
+ // 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++){
+ // 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));
+ I->second.orientation=mySe3.orientation.Conjugate()*subRBP->se3.orientation;
+ I->second.position=mySe3.orientation.Conjugate()*(subRBP->se3.position-mySe3.position);
+ }
+
+}
+
+/*! @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).
+ * @param m mass of the body; if positive, translation is away from the centroid; if negative, towards centroid.
+ * @param off offset of the new origin from the original origin
+ * @return inertia tensor in the new coordinate system; the matrix is symmetric.
+ */
+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);
+ // 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] */
+ I2+=m*Matrix3r(/* 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(1,0)=I2(0,1); I2(2,0)=I2(0,2); I2(2,1)=I2(1,2);
+ //TRWM3MAT(I2);
+ return I2;
+}
+
+/*! @brief Recalculate body's inertia tensor in rotated coordinates.
+ *
+ * @param I inertia tensor in old coordinates
+ * @param T rotation matrix from old to new coordinates
+ * @return inertia tensor in new coordinates
+ */
+Matrix3r Clump::inertiaTensorRotate(const Matrix3r& I,const Matrix3r& T){
+ /* [http://www.kwon3d.com/theory/moi/triten.html] */
+ //TRWM3MAT(I); TRWM3MAT(T);
+ return T.Transpose()*I*T;
+}
+
+/*! @brief Recalculate body's inertia tensor in rotated coordinates.
+ *
+ * @param I inertia tensor in old coordinates
+ * @param rot quaternion that describes rotation from old to new coordinates
+ * @return inertia tensor in new coordinates
+ */
+Matrix3r Clump::inertiaTensorRotate(const Matrix3r& I, const Quaternionr& rot){
+ Matrix3r T;
+ rot.ToRotationMatrix(T);
+ return inertiaTensorRotate(I,T);
+}
+
+
+
+/**************************************************************************************
+ ********************* ClumpTestGen ***************************************************
+ **************************************************************************************/
+
+#include<yade/core/MetaBody.hpp>
+#include<yade/extra/Shop.hpp>
+
+bool ClumpTestGen::generate()
+{
+ //Shop::setDefault("param_pythonRunExpr",string("if S.i%50==0 and S.i<1000 and S.i>500:\n\tprint S.i,len(S.sel),B[1].x, B[1].E"));
+
+ rootBody=Shop::rootBody();
+ Shop::rootBodyActors(rootBody);
+ // clumps do not need to subscribe currently (that will most likely change, though)
+ rootBody->engines.push_back(shared_ptr<ClumpMemberMover>(new ClumpMemberMover));
+
+
+ shared_ptr<MetaBody> oldRootBody=Omega::instance().getRootBody();
+ Omega::instance().setRootBody(rootBody);
+
+ shared_ptr<Body> ground=Shop::box(Vector3r(0,0,-1),Vector3r(3,3,.2));
+ ground->isDynamic=false;
+ // revert random colors for this single case...
+ ground->geometricalModel->diffuseColor=Vector3r(.6,.6,.6);
+ ground->interactingGeometry->diffuseColor=Vector3r(.3,.3,.3);
+ rootBody->bodies->insert(ground);
+
+ vector<Vector3r> relPos; vector<Real> radii; Vector3r clumpPos;
+
+ // standalone (non-clump!) sphere as well
+ shared_ptr<Body> sphere=Shop::sphere(Vector3r(0,0,0),.5);
+ rootBody->bodies->insert(sphere);
+ // one-sphere clump
+ clumpPos=Vector3r(-2,0,0);
+ relPos.push_back(Vector3r(0,0,0)); radii.push_back(.5);
+ createOneClump(rootBody,clumpPos,relPos,radii);
+ relPos.clear(); radii.clear();
+ // two-sphere clump
+ clumpPos=Vector3r(2,0,0);
+ relPos.push_back(Vector3r(0,-.5,.5)); radii.push_back(.5);
+ relPos.push_back(Vector3r(0,.5,0)); radii.push_back(.5);
+ createOneClump(rootBody,clumpPos,relPos,radii);
+ relPos.clear(); radii.clear();
+ // three-sphere slump
+ clumpPos=Vector3r(0,2,0);
+ relPos.push_back(Vector3r(0,-.5,.5)); radii.push_back(.5);
+ relPos.push_back(Vector3r(0,.5,0)); radii.push_back(.5);
+ relPos.push_back(Vector3r(.5,0,0)); radii.push_back(.5);
+ createOneClump(rootBody,clumpPos,relPos,radii);
+ relPos.clear(); radii.clear();
+ // four-sphere slump
+ clumpPos=Vector3r(0,-2,0);
+ relPos.push_back(Vector3r(0,0,0)); radii.push_back(.5);
+ relPos.push_back(Vector3r(.5,0,0)); radii.push_back(.5);
+ relPos.push_back(Vector3r(0,.5,0)); radii.push_back(.5);
+ relPos.push_back(Vector3r(0,0,.5)); radii.push_back(.5);
+ createOneClump(rootBody,clumpPos,relPos,radii);
+ relPos.clear(); radii.clear();
+
+ // restore Omega
+ Omega::instance().setRootBody(oldRootBody);
+
+ message="OK";
+ return true;
+}
+
+/*! \brief Generate clump of spheres, the result will be inserted into rootBody.
+ *
+ * Attention here: clump's id must be greater than id of any of its constituents; therefore
+ * 1. create bodies that will be clumped, add them to bodies (at that moment they get their id) and save their ids in clumpMembers
+ * 2. create (empty) clump and add it to bodies
+ * 3. add bodies to be clumped to the clump
+ * 4. call Clump::updateProperties to get physical properties physically right (inertia, position, orientation, mass, ...).
+ *
+ * @param clumpPos Center of the clump (not necessarily centroid); serves merely as reference for sphere positions.
+ * @param relPos Relative positions of individual spheres' centers.
+ * @param radii Radii of composing spheres. Must have the same length as relPos.
+ */
+void ClumpTestGen::createOneClump(shared_ptr<MetaBody>& rootBody, Vector3r clumpPos, vector<Vector3r> relPos, vector<Real> radii)
+{
+ assert(relPos.size()==radii.size());
+ list<body_id_t> clumpMembers;
+ for(size_t i=0; i<relPos.size(); i++){
+ shared_ptr<Body> sphere=Shop::sphere(clumpPos+relPos[i],radii[i]);
+ body_id_t lastId=(body_id_t)rootBody->bodies->insert(sphere);
+ clumpMembers.push_back(lastId);
+ LOG_TRACE("Generated (not yet) clumped sphere #"<<lastId);
+ }
+ shared_ptr<Clump> clump=shared_ptr<Clump>(new Clump());
+ shared_ptr<Body> clumpAsBody=static_pointer_cast<Body>(clump);
+ clump->isDynamic=true;
+ rootBody->bodies->insert(clumpAsBody);
+ FOREACH(body_id_t id, clumpMembers){
+ clump->add(id);
+ }
+ clump->updateProperties(false);
+}
+
Property changes on: trunk/pkg/dem/DataClass/Clump.cpp
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: trunk/pkg/dem/DataClass/Clump.hpp (from rev 1763, trunk/extra/clump/Clump.hpp)
===================================================================
--- trunk/extra/clump/Clump.hpp 2009-05-03 07:16:04 UTC (rev 1763)
+++ trunk/pkg/dem/DataClass/Clump.hpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -0,0 +1,115 @@
+// (c) 2007 Vaclav Smilauer <eudoxos@xxxxxxxx>
+
+#pragma once
+
+#include<vector>
+#include<map>
+#include<yade/core/Body.hpp>
+#include<yade/core/MetaBody.hpp>
+#include<yade/core/FileGenerator.hpp>
+#include<yade/core/DeusExMachina.hpp>
+#include<yade/lib-factory/Factorable.hpp>
+#include<yade/pkg-common/PhysicalParametersEngineUnit.hpp>
+#include<yade/pkg-common/RigidBodyParameters.hpp>
+#include<yade/pkg-common/AABB.hpp>
+#include<yade/lib-base/Logging.hpp>
+#include<yade/lib-base/yadeWm3Extra.hpp>
+
+
+/*! Body representing clump (rigid aggregate) composed by other existing bodies.
+
+ Clump is one of bodies that reside in rootBody->bodies.
+ When an existing body is added to ::Clump, it's ::Body::isDynamic flag is set to false
+ (it is still subscribed to all its engines, to make it possible to remove it from the clump again).
+ All forces acting on Clump::members are made to act on the clump itself, which will ensure that they
+ influence all Clump::members as if the clump were a rigid particle.
+
+ What are clump requirements so that they function?
+ -# Given any body, tell
+ - if it is a clump member: Body::isClumpMember()
+ - if it is a clump: Body:: isClump(). (Correct result is assured at each call to Clump::add).
+ (we could use RTTI instead? Would that be more reliable?)
+ - if it is a standalone Body: Body::isStandalone()
+ - what is it's clump id (Body::clumpId)
+ -# given the root body, tell
+ - what clumps it contains (enumerate all bodies and filter clumps, see above)
+ -# given a clump, tell
+ - what bodies it contains (keys of ::Clump::members)
+ - what are se3 of these bodies (values of ::Clump::members)
+ -# add/delete bodies from/to clump (::Clump::add, ::Clump::del)
+ - This includes saving se3 of the subBody: it \em must be in clump's local coordinates so that it is constant. The transformation from global to local is given by clump's se3 at the moment of addition. Clump's se3 is initially (origin,identity)
+ -# Update clump's physical properties (Clump::updateProperties)
+ - This \em must reposition members so that they have the same se3 globally
+ -# Apply forces acting on members to the clump instead (done in NewtonsForceLaw, NewtonsMomentumLaw) - uses world coordinates to calculate effect on the clump's centroid
+ -# Integrate position and orientation of the clump
+ - LeapFrogPositionIntegrator and LeapFrogOrientationIntegrator move clump as whole
+ - clump members are skipped, since they have Body::isDynamic==false.
+ - ClumpMemberMover is an engine that updates positions of the clump memebers in each timestep (calls Clump::moveSubBodies internally)
+
+ Some more information can be found http://beta.arcig.cz/~eudoxos/phd/index.cgi/YaDe/HighLevelClumps
+
+ For an example how to generate a clump, see ClumpTestGen::createOneClump.
+
+ @todo GravityEngine should be applied to members, not to clump as such?! Still not sure. Perhaps Clumps should have mass and inertia set to zeros so that engines unaware of clumps do not act on it. It would have some private mass and insertia that would be used in NewtonsForceLaw etc for clumps specially...
+
+ @note PersistentSAPCollider bypass Clumps explicitly. This no longer depends on the absence of boundingVolume.
+ @note Clump relies on its id being assigned (as well as id of its components); therefore, only bodies that have already been inserted to the container may be added to Clump which has been itself already added to the container. We further requier that clump id is greater than ids of clumped bodies
+
+ */
+
+class Clump: public Body {
+ //! mapping of body IDs to their relative positions; replaces members and subSe3s;
+ public:
+ typedef std::map<body_id_t,Se3r> memberMap;
+ memberMap members;
+
+ Clump();
+ virtual ~Clump(){};
+ //! \brief add Body to the Clump
+ void add(body_id_t);
+ //! \brief remove Body from the Clump
+ void del(body_id_t);
+ //! Recalculate physical properties of Clump.
+ void updateProperties(bool intersecting);
+ //! Calculate positions and orientations of members based on my own Se3.
+ void moveMembers();
+ //! 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...
+ //! 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.
+ static Matrix3r inertiaTensorRotate(const Matrix3r& I, const Matrix3r& T);
+ //! Recalculate body's inertia tensor in rotated coordinates.
+ static Matrix3r inertiaTensorRotate(const Matrix3r& I, const Quaternionr& rot);
+ REGISTER_ATTRIBUTES(Body,(members));
+ REGISTER_CLASS_AND_BASE(Clump,Body);
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Clump);
+
+/*! Update ::Clump::members positions so that the Clump behaves as a rigid body. */
+class ClumpMemberMover: public DeusExMachina {
+ public:
+ //! Interates over rootBody->bodies and calls Clump::moveSubBodies() for clumps.
+ virtual void applyCondition(MetaBody* rootBody);
+ ClumpMemberMover();
+ virtual ~ClumpMemberMover(){};
+ REGISTER_CLASS_AND_BASE(ClumpMemberMover,DeusExMachina);
+ REGISTER_ATTRIBUTES(DeusExMachina,/*nothing here*/);
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(ClumpMemberMover);
+
+/*! \brief Test some basic clump functionality; show how to use clumps as well. */
+class ClumpTestGen : public FileGenerator {
+ void createOneClump(shared_ptr<MetaBody>& rootBody, Vector3r clumpPos, vector<Vector3r> relPos, vector<Real> radii);
+ shared_ptr<ClumpMemberMover> clumpMover;
+ public :
+ bool generate();
+ DECLARE_LOGGER;
+ REGISTER_CLASS_AND_BASE(ClumpTestGen,FileGenerator);
+ REGISTER_ATTRIBUTES(FileGenerator,/*nothing here*/);
+};
+REGISTER_SERIALIZABLE(ClumpTestGen);
+
Property changes on: trunk/pkg/dem/DataClass/Clump.hpp
___________________________________________________________________
Name: svn:mergeinfo
+
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -117,106 +117,20 @@
}
-
-
-Vector3r SpheresContactGeometry::relRotVector() const{
- Quaternionr relOri12=ori1.Conjugate()*ori2;
- Quaternionr oriDiff=initRelOri12.Conjugate()*relOri12;
- Vector3r axis; Real angle;
- oriDiff.ToAxisAngle(axis,angle);
- if(angle>Mathr::PI)angle-=Mathr::TWO_PI;
- return angle*axis;
-}
-
-void SpheresContactGeometry::bendingTorsionAbs(Vector3r& bend, Real& tors){
- Vector3r relRot=relRotVector();
- tors=relRot.Dot(normal);
- bend=relRot-tors*normal;
-}
-
-/* Set contact points on both spheres such that their projection is the one given
- * (should be on the plane passing through origin and oriented with normal; not checked!)
- */
-void SpheresContactGeometry::setTgPlanePts(Vector3r p1new, Vector3r p2new){
- assert(hasShear);
- cp1rel=ori1.Conjugate()*rollPlanePtToSphere(p1new,d1,normal);
- cp2rel=ori2.Conjugate()*rollPlanePtToSphere(p2new,d2,-normal);
-}
-
-
-/* Move contact point on both spheres in such way that their relative position (displacementT) is the same;
- * this should be done regularly to ensure that the angle doesn't go over π, since then quaternion would
- * flip axis and the point would project on other side of the tangent plane piece.
- */
-void SpheresContactGeometry::relocateContactPoints(){
- assert(hasShear);
- relocateContactPoints(contPtInTgPlane1(),contPtInTgPlane2());
-}
-
-/* Like SpheresContactGeometry::relocateContactPoints(), but use already computed tangent plane points.
- *
- *
- */
-void SpheresContactGeometry::relocateContactPoints(const Vector3r& p1, const Vector3r& p2){
- Vector3r midPt=(d1/(d1+d2))*(p1+p2); // proportionally to radii, so that angle would be the same
- if((p1.SquaredLength()>pow(d1,2) || p2.SquaredLength()>pow(d2,2)) && midPt.SquaredLength()>pow(min(d1,d2),2)){
- //cerr<<"RELOCATION with displacementT="<<displacementT(); // should be the same before and after relocation
- setTgPlanePts(p1-midPt,p2-midPt);
- //cerr<<" → "<<displacementT()<<endl;
+/* keep this for reference; declarations in the header */
+#if 0
+ Vector3r SpheresContactGeometry::relRotVector() const{
+ Quaternionr relOri12=ori1.Conjugate()*ori2;
+ Quaternionr oriDiff=initRelOri12.Conjugate()*relOri12;
+ Vector3r axis; Real angle;
+ oriDiff.ToAxisAngle(axis,angle);
+ if(angle>Mathr::PI)angle-=Mathr::TWO_PI;
+ return angle*axis;
}
-}
-
-/*! Perform slip of the projected contact points so that their distance becomes equal (or remains smaller) than the given one.
- * The slipped distance is returned.
- */
-Real SpheresContactGeometry::slipToDisplacementTMax(Real displacementTMax){
- assert(hasShear);
- // very close, reset shear
- if(displacementTMax<=0.){ setTgPlanePts(Vector3r(0,0,0),Vector3r(0,0,0)); return displacementTMax;}
- // otherwise
- Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
- Real currDistSq=(p2-p1).SquaredLength();
- if(currDistSq<pow(displacementTMax,2)) return 0; // close enough, no slip needed
- Vector3r diff=.5*(sqrt(currDistSq)/displacementTMax-1)*(p2-p1);
- setTgPlanePts(p1+diff,p2-diff);
- return 2*diff.Length();
-}
-
-/*! Project point from sphere surface to tangent plane,
- * such that the angle of shortest arc from (1,0,0) pt on the sphere to the point itself is the same
- * as the angle of segment of the same length on the tangent plane.
- *
- * This function is (or should be) inverse of SpheresContactGeometry::rollPlanePtToSphere.
- *
- * @param fromXtoPtOri gives orientation of the vector from sphere center to the sphere point from the global +x axis.
- * @param radius the distance from sphere center to the contact plane
- * @param planeNormal unit vector pointing away from the sphere center, determining plane orientation on which the projected point lies.
- * @returns The projected point coordinates (with origin at the contact point).
- */
-Vector3r SpheresContactGeometry::unrollSpherePtToPlane(const Quaternionr& fromXtoPtOri, const Real& radius, const Vector3r& planeNormal){
- Quaternionr normal2pt; normal2pt.Align(planeNormal,fromXtoPtOri*Vector3r::UNIT_X);
- Vector3r axis; Real angle; normal2pt.ToAxisAngle(axis,angle);
- return (angle*radius) /* length */ *(axis.Cross(planeNormal)) /* direction: both are unit vectors */;
-}
-
-/*! Project point from tangent plane to the sphere.
- *
- * This function is (or should be) inverse of SpheresContactGeometry::unrollSpherePtToPlane.
- *
- * @param planePt point on the tangent plane, with origin at the contact point (i.e. at sphere center + normal*radius)
- * @param radius sphere radius
- * @param planeNorma _unit_ vector pointing away from sphere center
- * @returns orientation that transforms +x axis to the vector between sphere center and point on the sphere that corresponds to planePt.
- *
- * @note It is not checked whether planePt relly lies on the tangent plane. If not, result will be incorrect.
- */
-Quaternionr SpheresContactGeometry::rollPlanePtToSphere(const Vector3r& planePt, const Real& radius, const Vector3r& planeNormal){
- Vector3r axis=planeNormal.Cross(planePt); axis.Normalize();
- Real angle=planePt.Length()/radius;
- Quaternionr normal2pt(axis,angle);
- Quaternionr ret; ret.Align(Vector3r::UNIT_X,normal2pt*planeNormal);
- return ret;
-}
-
-
+ void SpheresContactGeometry::bendingTorsionAbs(Vector3r& bend, Real& tors){
+ Vector3r relRot=relRotVector();
+ tors=relRot.Dot(normal);
+ bend=relRot-tors*normal;
+ }
+#endif
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -29,70 +29,17 @@
Vector3r updateShear(const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting=true);
#endif
- // all the rest here will hopefully disappear at some point...
+ /* keep this for reference, until it is implemented elsewhere, like Dem6DofGeom, if it ever exists */
+ #if 0
+ // initial relative orientation, used for bending and torsion computation
+ Quaternionr initRelOri12;
+ void bendingTorsionAbs(Vector3r& bend, Real& tors);
+ void bendingTorsionRel(Vector3r& bend, Real& tors){ bendingTorsionAbs(bend,tors); bend/=d0; tors/=d0;}
+ Vector3r relRotVector() const;
+ #endif
- // begin abusive storage
- //! Whether the original contact was on the positive (1) or negative (-1) facet side; this is to permit repulsion to the right side even if the sphere passes, under extreme pressure, geometrically to the other facet's side. This is used only in InteractingFacet2IteractingSphere4SpheresContactGeometry. If at any point the contact is with the edge or vertex instead of face, this attribute is reset so that contact with the other face is possible.
- int facetContactFace;
- // end abusive storage
- bool hasShear; // whether the exact rotation code is being used -- following fields are needed for that
- //! positions and orientations of both spheres -- must be updated at every iteration
- Vector3r pos1, pos2; Quaternionr ori1, ori2;
- /*! Orientation of the contact point relative to each sphere-local coordinates.+
- * Those fields are almost constant, except for a few cases
- * (a) plastic slip and
- * (b) spheres mutually rolling without slipping with big angle, when the contact point must be moved since quaternions
- * describe only rotations up to π (if the rotation is bigger, then the shorter path from the other side is
- * taken instead)
- */
- Quaternionr cp1rel, cp2rel;
- // interaction "radii" and total length; this is _really_ constant throughout the interaction life
- // d1 is really distance from the sphere1 center to the contact plane, it may not correspond to the sphere radius!
- // therefore, d1+d2=d0 (distance at which the contact was created)
- // d0fixup is added to d0 when computing normal strain; it should fix problems with sphere-facet interactions never getting enough compressed
- Real d1, d2, d0, d0fixup;
- // initial relative orientation, used for bending and torsion computation
- Quaternionr initRelOri12;
-
- // auxiliary functions; the quaternion magic is all in here
- static Vector3r unrollSpherePtToPlane(const Quaternionr& fromXtoPtOri, const Real& radius, const Vector3r& normal);
- static Quaternionr rollPlanePtToSphere(const Vector3r& planePt, const Real& radius, const Vector3r& normal);
-
- void setTgPlanePts(Vector3r p1new, Vector3r p2new);
-
- Vector3r contPtInTgPlane1() const {assert(hasShear); return unrollSpherePtToPlane(ori1*cp1rel,d1,normal);}
- Vector3r contPtInTgPlane2() const {assert(hasShear); return unrollSpherePtToPlane(ori2*cp2rel,d2,-normal);}
- Vector3r contPt() const {return contactPoint; /*pos1+(d1/d0)*(pos2-pos1)*/}
-
- Real displacementN() const {assert(hasShear); return (pos2-pos1).Length()-(d0+d0fixup);}
- Real epsN() const {return displacementN()*(1./(d0+d0fixup));}
- Vector3r displacementT() { assert(hasShear);
- // enabling automatic relocation decreases overall simulation speed by about 3%
- // perhaps: bool largeStrains ... ?
- #if 1
- Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
- relocateContactPoints(p1,p2);
- return p2-p1; // shear before relocation, but that is OK
- #else
- return contPtInTgPlane2()-contPtInTgPlane1();
- #endif
- }
- Vector3r epsT() {return displacementT()*(1./(d0+d0fixup));}
-
- Real slipToDisplacementTMax(Real displacementTMax);
- //! slip to epsTMax if current epsT is greater; return the relative slip magnitude
- Real slipToStrainTMax(Real epsTMax){ return (1/d0)*slipToDisplacementTMax(epsTMax*d0);}
-
- void relocateContactPoints();
- void relocateContactPoints(const Vector3r& tgPlanePt1, const Vector3r& tgPlanePt2);
-
- void bendingTorsionAbs(Vector3r& bend, Real& tors);
- void bendingTorsionRel(Vector3r& bend, Real& tors){ bendingTorsionAbs(bend,tors); bend/=d0; tors/=d0;}
-
- Vector3r relRotVector() const;
-
- SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0.),radius2(0.),facetContactFace(0.),hasShear(false),pos1(Vector3r::ZERO),pos2(Vector3r::ZERO),ori1(Quaternionr::IDENTITY),ori2(Quaternionr::IDENTITY),cp1rel(Quaternionr::IDENTITY),cp2rel(Quaternionr::IDENTITY),d1(0),d2(0),d0(0),d0fixup(0),initRelOri12(Quaternionr::IDENTITY){createIndex();
+ SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0.),radius2(0.){createIndex();
#ifdef SCG_SHEAR
shear=Vector3r::ZERO; prevNormal=Vector3r::ZERO /*initialized to proper value by geom functor*/;
#endif
@@ -110,20 +57,7 @@
(shear)
(prevNormal)
#endif
- (facetContactFace)
- // hasShear
- (hasShear)
- (pos1)
- (pos2)
- (ori1)
- (ori2)
- (cp1rel)
- (cp2rel)
- (d1)
- (d2)
- (d0)
- (d0fixup)
- (initRelOri12));
+ );
REGISTER_CLASS_AND_BASE(SpheresContactGeometry,InteractionGeometry);
REGISTER_CLASS_INDEX(SpheresContactGeometry,InteractionGeometry);
};
Modified: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -10,6 +10,7 @@
#include<yade/core/MetaBody.hpp>
#include<yade/pkg-common/RigidBodyParameters.hpp>
#include<yade/lib-base/yadeWm3Extra.hpp>
+#include<yade/pkg-dem/Clump.hpp>
YADE_PLUGIN("NewtonsDampedLaw");
@@ -29,7 +30,8 @@
{
ncb->bex.sync();
FOREACH(const shared_ptr<Body>& b, *ncb->bodies){
- if (!b->isDynamic) continue;
+ // clump members are non-dynamic; they skip the rest of loop once their forces are properly taken into account, however
+ if (!b->isDynamic && !b->isClumpMember()) continue;
RigidBodyParameters* rb = YADE_CAST<RigidBodyParameters*>(b->physicalParameters.get());
unsigned int id = b->getId();
@@ -39,45 +41,40 @@
Real dt = Omega::instance().getTimeStep();
//Newtons mometum law :
- if ( b->isStandalone() ) rb->angularAcceleration=diagDiv ( m,rb->inertia );
- else if ( b->isClump() ) rb->angularAcceleration+=diagDiv ( m,rb->inertia );
- else
- { // isClumpMember()
- const shared_ptr<Body>& clump ( Body::byId ( b->clumpId ) );
- RigidBodyParameters* clumpRBP=YADE_CAST<RigidBodyParameters*> ( clump->physicalParameters.get() );
- /* angularAcceleration is reset by ClumpMemberMover engine */
- clumpRBP->angularAcceleration+=diagDiv ( m,clumpRBP->inertia );
+ if (b->isStandalone()){
+ rb->angularAcceleration=diagDiv(m,rb->inertia);
+ rb->acceleration=f/rb->mass;
}
-
- // Newtons force law
- if ( b->isStandalone() ) rb->acceleration=f/rb->mass;
- else if ( b->isClump() )
- {
- // accel for clump reset in Clump::moveMembers, called by ClumpMemberMover engine
- rb->acceleration+=f/rb->mass;
+ else if (b->isClump()){
+ // at this point, forces from clump members are already summed up, this is just for forces applied to clump proper, if there are such (does it have some physical meaning?)
+ rb->angularAcceleration+=diagDiv(m,rb->inertia);
+ rb->acceleration+=f/rb->mass; // accel for clump will be reset in Clump::moveMembers, called from the clump itself
}
- else
- { // force applied to a clump member is applied to clump itself
+ else {
+ assert(b->isClumpMember());
+ assert(b->clumpId>b->id);
const shared_ptr<Body>& clump ( Body::byId ( b->clumpId ) );
RigidBodyParameters* clumpRBP=YADE_CAST<RigidBodyParameters*> ( clump->physicalParameters.get() );
- // accels reset by Clump::moveMembers in last iteration
+ /* angularAcceleration is reset by ClumpMemberMover engine */
+ clumpRBP->angularAcceleration+=diagDiv ( m,clumpRBP->inertia );
clumpRBP->acceleration+=f/clumpRBP->mass;
clumpRBP->angularAcceleration+=diagDiv ( ( rb->se3.position-clumpRBP->se3.position ).Cross ( f ),clumpRBP->inertia ); //acceleration from torque generated by the force WRT particle centroid on the clump centroid
+ continue;
}
- if (!b->isClump()) {
+ // damping: applied to both clumps and non-clumps
+ // FIXME: clumps with just 1 body are not dumped? See sphere jumping forever in scripts/test/clump.py
+ // other clumps seem to be damped in some strange way (again, see the same script).
+ // The damping also looks different from ClumpTestGen, which uses separate engines instead of NewtonsDampedLaw
+ for ( int i=0; i<3; ++i ) {
//Damping moments
- for ( int i=0; i<3; ++i )
- {
- rb->angularAcceleration[i] *= 1 - damping*Mathr::Sign ( m[i]*(rb->angularVelocity[i] + (Real) 0.5 *dt*rb->angularAcceleration[i]) );
- }
+ rb->angularAcceleration[i] *= 1 - damping*Mathr::Sign ( m[i]*(rb->angularVelocity[i] + (Real) 0.5 *dt*rb->angularAcceleration[i]) );
//Damping force :
- for ( int i=0; i<3; ++i )
- {
- rb->acceleration[i] *= 1 - damping*Mathr::Sign ( f[i]*(rb->velocity[i] + (Real) 0.5 *dt*rb->acceleration[i]) );
- }
+ rb->acceleration[i] *= 1 - damping*Mathr::Sign ( f[i]*(rb->velocity[i] + (Real) 0.5 *dt*rb->acceleration[i]) );
}
+
+ // blocking DOFs
if(rb->blockedDOFs==0){ /* same as: rb->blockedDOFs==PhysicalParameters::DOF_NONE */
rb->angularVelocity=rb->angularVelocity+dt*rb->angularAcceleration;
rb->velocity=rb->velocity+dt*rb->acceleration;
@@ -98,6 +95,8 @@
rb->se3.orientation.Normalize();
rb->se3.position += rb->velocity*dt;
+
+ if(b->isClump()) static_cast<Clump*>(b.get())->moveMembers();
}
}
@@ -106,8 +105,7 @@
[16:09:37] eudoxos2: for(int state=LOOP1; state!=END; state++){
[16:09:37] eudoxos2: FOREACH(const shared_ptr<Body>& b, rootBody->bodies){
[16:09:38] eudoxos2: if(b->isClumpMember() && LOOP1){ [[apply that on b->clumpId]] }
-[16:09:38] eudoxos2: if((b->isStandalone && LOOP1) || (b->isClump && LOOP2){ [[damping, newton, integrate]] }
-[16:09:38] eudoxos2: if(b->isClump() && LOOP 2){ b->moveMembers(); }
+[16:09:38] eudoxos2: if((b->isStandalone && LOOP1) || (b->isClump && LOOP2){ [[damping, newton, integrate]]; b->moveClumpMembers(); }
[16:09:40] eudoxos2: }
[16:09:42] eudoxos2: }
[16:09:44] eudoxos2: }*/
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -18,14 +18,12 @@
InteractingFacet2InteractingSphere4SpheresContactGeometry::InteractingFacet2InteractingSphere4SpheresContactGeometry()
{
shrinkFactor=0;
- hasShear=false;
}
void InteractingFacet2InteractingSphere4SpheresContactGeometry::registerAttributes()
{
InteractionGeometryEngineUnit::registerAttributes();
REGISTER_ATTRIBUTE(shrinkFactor);
- REGISTER_ATTRIBUTE(hasShear);
}
bool InteractingFacet2InteractingSphere4SpheresContactGeometry::go(const shared_ptr<InteractingGeometry>& cm1,
@@ -49,7 +47,9 @@
Vector3r normal = facet->nf;
Real L = normal.Dot(cl);
+ if (L<0) {normal=-normal; L=-L; }
+#if 0
int contactFace=0; // temp to save what will be maybe needed for new contact
//assert((c->interactionGeometry&&c->isReal)||(!c->interactionGeometry&&!c->isReal));
if(c->interactionGeometry){ // contact already exists, use old data here
@@ -64,6 +64,7 @@
if (L<0) { normal=-normal; L=-L; contactFace=-1;} // new contact on the negative face, reverse and save that information so that since now this contact is always reversed
else contactFace=1;
}
+#endif
Real sphereRadius = static_cast<InteractingSphere*>(cm2.get())->radius;
if (L>sphereRadius && !c->isReal) return false; // no contact, but only if there was no previous contact; ortherwise, the constitutive law is responsible for setting Interaction::isReal=false
@@ -100,8 +101,6 @@
}
else
{
- // edge or vertex contact: become indeterminate with respect to face
- contactFace=0;
cp = cp + ne[m]*(icr-bm);
if (cp.Dot(ne[(m-1<0)?2:m-1])>icr) // contact with vertex m
// cp = facet->vertices[m];
@@ -123,8 +122,7 @@
if (c->interactionGeometry)
scm = YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
else
- { scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
- scm->facetContactFace=contactFace; }
+ scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
normal = facetAxisT*normal; // in global orientation
scm->contactPoint = se32.position - (sphereRadius-0.5*penetrationDepth)*normal;
@@ -136,23 +134,6 @@
if (!c->interactionGeometry)
c->interactionGeometry = scm;
- if(hasShear){
- scm->hasShear=true;
- // fictive center of the sphere representing the facet in the sphere-sphere contact
- scm->pos1=scm->contactPoint-(scm->radius1-.5*penetrationDepth)*normal; scm->pos2=se32.position;
- scm->ori1=se31.orientation; scm->ori2=se32.orientation;
- if(c->isNew){
- scm->d0=scm->radius1+scm->radius2-penetrationDepth;
- scm->d1=scm->radius1-.5*penetrationDepth; scm->d2=scm->radius2-.5*penetrationDepth;
- scm->d0fixup=-scm->radius1-.5*penetrationDepth;
- // quasi-constants
- scm->cp1rel.Align(Vector3r::UNIT_X,se31.orientation.Conjugate()*normal);
- scm->cp2rel.Align(Vector3r::UNIT_X,se32.orientation.Conjugate()*(-normal));
- scm->cp1rel.Normalize(); scm->cp2rel.Normalize();
- }
- TRVAR2(scm->contPtInTgPlane1(),scm->contPtInTgPlane2());
- }
-
return true;
}
return false;
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.hpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.hpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -44,8 +44,6 @@
// multiplied by 2*shrinkFactor. By default shrinkFactor=0.
Real shrinkFactor;
- bool hasShear;
-
protected :
virtual void registerAttributes();
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -14,7 +14,6 @@
InteractingSphere2InteractingSphere4SpheresContactGeometry::InteractingSphere2InteractingSphere4SpheresContactGeometry()
{
interactionDetectionFactor = 1;
- hasShear=false;
}
bool InteractingSphere2InteractingSphere4SpheresContactGeometry::go( const shared_ptr<InteractingGeometry>& cm1,
@@ -42,28 +41,10 @@
scm->penetrationDepth=penetrationDepth;
scm->radius1=s1->radius;
scm->radius2=s2->radius;
- if(hasShear){
- scm->hasShear=true;
- scm->pos1=se31.position; scm->pos2=se32.position;
- scm->ori1=se31.orientation; scm->ori2=se32.orientation;
- if(c->isNew){
- // contact constants
- scm->d0=(se32.position-se31.position).Length();
- scm->d1=s1->radius-.5*penetrationDepth; scm->d2=s2->radius-.5*penetrationDepth;
- scm->initRelOri12=se31.orientation.Conjugate()*se32.orientation;
- // quasi-constants
- scm->cp1rel.Align(Vector3r::UNIT_X,se31.orientation.Conjugate()*normal);
- scm->cp2rel.Align(Vector3r::UNIT_X,se32.orientation.Conjugate()*(-normal));
- scm->cp1rel.Normalize(); scm->cp2rel.Normalize();
- }
- // testing only
- #if 0
- if((Omega::instance().getCurrentIteration()%10000)==0) scm->relocateContactPoints();
- if((Omega::instance().getCurrentIteration()%10000)==0) {
- Real slip=scm->slipToEpsNMax(1.); if(slip>0) cerr<<"SLIP by Δε_N "<<slip<<" → ε_N="<<scm->epsN()<<endl;
- }
- #endif
- }
+ // keep this for reference on how to compute bending and torsion from relative orientation; parts in SpheresContactGeometry header
+ #if 0
+ scm->initRelOri12=se31.orientation.Conjugate()*se32.orientation;
+ #endif
return true;
}
return false;
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -25,11 +25,9 @@
* @note This parameter is functionally coupled with InteractinSphere2AABB::aabbEnlargeFactor,
* which will create larger bounding boxes and should be of the same value. */
double interactionDetectionFactor;
- /*! Whether we create SpheresContactGeometry with data necessary for (exact) shear computation. By default false */
- bool hasShear;
REGISTER_CLASS_AND_BASE(InteractingSphere2InteractingSphere4SpheresContactGeometry,InteractionGeometryEngineUnit);
- REGISTER_ATTRIBUTES(InteractionGeometryEngineUnit,(interactionDetectionFactor)(hasShear));
+ REGISTER_ATTRIBUTES(InteractionGeometryEngineUnit,(interactionDetectionFactor));
FUNCTOR2D(InteractingSphere,InteractingSphere);
// needed for the dispatcher, even if it is symmetric
DEFINE_FUNCTOR_ORDER_2D(InteractingSphere,InteractingSphere);
Copied: trunk/pkg/dem/Engine/EngineUnit/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.cpp (from rev 1763, trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.cpp)
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.cpp 2009-05-03 07:16:04 UTC (rev 1763)
+++ trunk/pkg/dem/Engine/EngineUnit/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -0,0 +1,17 @@
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+#include"ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.hpp"
+#include<yade/pkg-dem/DemXDofGeom.hpp>
+#include<yade/pkg-common/NormalShearInteractions.hpp>
+YADE_PLUGIN("ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw");
+
+void ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
+ Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(_geom.get());
+ NormalShearInteraction* phys=static_cast<NormalShearInteraction*>(_phys.get());
+
+ // if (geom->displacementN()>0) return; // non-cohesive behavior
+ phys->normalForce=(geom->displacementN()*phys->kn)*geom->normal;
+ phys->shearForce=geom->displacementT()*phys->ks;
+
+ applyForceAtContactPoint(phys->normalForce+phys->shearForce, geom->contactPoint, I->getId1(), geom->se31.position, I->getId2(), geom->se32.position, rootBody);
+}
+
Property changes on: trunk/pkg/dem/Engine/EngineUnit/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.cpp
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: trunk/pkg/dem/Engine/EngineUnit/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.hpp (from rev 1763, trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.hpp)
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.hpp 2009-05-03 07:16:04 UTC (rev 1763)
+++ trunk/pkg/dem/Engine/EngineUnit/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.hpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -0,0 +1,10 @@
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+#include<yade/pkg-common/ConstitutiveLaw.hpp>
+/* Experimental constitutive law using the ConstitutiveLawDispatcher.
+ * Has only purely elastic normal and shear components. */
+class ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw: public ConstitutiveLaw {
+ virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*, MetaBody*);
+ FUNCTOR2D(Dem3DofGeom,NormalShearInteraction);
+ REGISTER_CLASS_AND_BASE(ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw,ConstitutiveLaw);
+};
+REGISTER_SERIALIZABLE(ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw);
Property changes on: trunk/pkg/dem/Engine/EngineUnit/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.hpp
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -1,17 +0,0 @@
-// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
-#include"ef2_Spheres_NormalShear_ElasticFrictionalLaw.hpp"
-#include<yade/pkg-dem/SpheresContactGeometry.hpp>
-#include<yade/pkg-common/NormalShearInteractions.hpp>
-YADE_PLUGIN("ef2_Spheres_NormalShear_ElasticFrictionalLaw");
-
-void ef2_Spheres_NormalShear_ElasticFrictionalLaw::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
- SpheresContactGeometry* geom=static_cast<SpheresContactGeometry*>(_geom.get());
- NormalShearInteraction* phys=static_cast<NormalShearInteraction*>(_phys.get());
-
- // if (geom->displacementN()>0) return; // non-cohesive behavior
- phys->normalForce=(geom->displacementN()*phys->kn)*geom->normal;
- phys->shearForce=geom->displacementT()*phys->ks;
-
- applyForceAtContactPoint(phys->normalForce+phys->shearForce, geom->contPt(), I->getId1(), geom->pos1, I->getId2(), geom->pos2, rootBody);
-}
-
Deleted: trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.hpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.hpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -1,10 +0,0 @@
-// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
-#include<yade/pkg-common/ConstitutiveLaw.hpp>
-/* Experimental constitutive law using the ConstitutiveLawDispatcher.
- * Has only purely elastic normal and shear components. */
-class ef2_Spheres_NormalShear_ElasticFrictionalLaw: public ConstitutiveLaw {
- virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*, MetaBody*);
- FUNCTOR2D(SpheresContactGeometry,NormalShearInteraction);
- REGISTER_CLASS_AND_BASE(ef2_Spheres_NormalShear_ElasticFrictionalLaw,ConstitutiveLaw);
-};
-REGISTER_SERIALIZABLE(ef2_Spheres_NormalShear_ElasticFrictionalLaw);
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -17,40 +17,9 @@
#include<yade/extra/Shop.hpp>
-YADE_PLUGIN("ElasticContactLaw2","ef2_Spheres_Elastic_ElasticLaw","ef2_Dem3Dof_Elastic_ElasticLaw","ElasticContactLaw");
+YADE_PLUGIN("ef2_Spheres_Elastic_ElasticLaw","ef2_Dem3Dof_Elastic_ElasticLaw","ElasticContactLaw");
-ElasticContactLaw2::ElasticContactLaw2(){
- isCohesive=true;
-}
-ElasticContactLaw2::~ElasticContactLaw2(){}
-
-void ElasticContactLaw2::action(MetaBody* rb){
- Real /* bending stiffness */ kb=1e7, /* torsion stiffness */ ktor=1e8;
- FOREACH(shared_ptr<Interaction> i, *rb->transientInteractions){
- if(!i->isReal) continue;
- shared_ptr<SpheresContactGeometry> contGeom=YADE_PTR_CAST<SpheresContactGeometry>(i->interactionGeometry);
- shared_ptr<ElasticContactInteraction> contPhys=YADE_PTR_CAST<ElasticContactInteraction>(i->interactionPhysics);
- assert(contGeom); assert(contPhys);
- if(!contGeom->hasShear) throw runtime_error("SpheresContactGeometry::hasShear must be true for ElasticContactLaw2");
- Real Fn=contPhys->kn*contGeom->displacementN(); // scalar normal force; displacementN()>=0 ≡ elongation of the contact
- if(!isCohesive && contGeom->displacementN()>0){ cerr<<"deleting"<<endl; /* delete the interaction */ i->isReal=false; continue;}
- contPhys->normalForce=Fn*contGeom->normal;
- //contGeom->relocateContactPoints();
- //contGeom->slipToDisplacementTMax(max(0.,(-Fn*contPhys->tangensOfFrictionAngle)/contPhys->ks)); // limit shear displacement -- Coulomb criterion
- contPhys->shearForce=contPhys->ks*contGeom->displacementT();
- Vector3r force=contPhys->shearForce+contPhys->normalForce;
- Shop::applyForceAtContactPoint(force,contGeom->contactPoint,i->getId1(),contGeom->pos1,i->getId2(),contGeom->pos2,rb);
-
- Vector3r bendAbs; Real torsionAbs; contGeom->bendingTorsionAbs(bendAbs,torsionAbs);
- rb->bex.addTorque(i->getId1(), contGeom->normal*torsionAbs*ktor+bendAbs*kb);
- rb->bex.addTorque(i->getId2(),-contGeom->normal*torsionAbs*ktor-bendAbs*kb);
- }
-}
-
-
-
-
ElasticContactLaw::ElasticContactLaw() : InteractionSolver()
{
sdecGroupMask=1;
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -18,22 +18,6 @@
#include <boost/tuple/tuple.hpp>
-class ElasticContactLaw2: public InteractionSolver{
- public:
- //! this should really be property of the interaction, but for simplicity keep it here now...
- bool isCohesive;
- ElasticContactLaw2();
- virtual ~ElasticContactLaw2();
- void action(MetaBody*);
- void registerAttributes(){
- InteractionSolver::registerAttributes();
- REGISTER_ATTRIBUTE(isCohesive);
- }
- REGISTER_CLASS_NAME(ElasticContactLaw2);
- REGISTER_BASE_CLASS_NAME(InteractionSolver);
-};
-REGISTER_SERIALIZABLE(ElasticContactLaw2);
-
class ef2_Spheres_Elastic_ElasticLaw: public ConstitutiveLaw{
public:
virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody);
Modified: trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp 2009-05-05 05:15:26 UTC (rev 1766)
@@ -77,43 +77,5 @@
#endif
}
- if(sc->hasShear){
- Vector3r pos1=sc->pos1, pos2=sc->pos2, contPt=sc->contPt();
- #if 0
- GLUtils::GLDrawArrow(contPt,contPt+sc->normal*.5*sc->d0); // normal of the contact
- #endif
- #if 0
- //Vector3r contPt=se31.position+(sc->d1/sc->d0)*(se32.position-se31.position); // must be recalculated to not be unscaled if scaling displacements ...
- GLUtils::GLDrawLine(pos1,pos2,Vector3r(.5,.5,.5));
- Vector3r bend; Real tors;
- sc->bendingTorsionRel(bend,tors);
- GLUtils::GLDrawLine(contPt,contPt+10*sc->radius1*(bend+sc->normal*tors),Vector3r(1,0,0));
- #endif
- #if 0
- GLUtils::GLDrawNum(bend[0],contPt-.2*sc->normal*sc->radius1,Vector3r(1,0,0));
- GLUtils::GLDrawNum(bend[1],contPt,Vector3r(0,1,0));
- GLUtils::GLDrawNum(bend[2],contPt+.2*sc->normal*sc->radius1,Vector3r(0,0,1));
- GLUtils::GLDrawNum(tors,contPt+.5*sc->normal*sc->radius2,Vector3r(1,1,0));
- #endif
- // sphere center to point on the sphere
- //GLUtils::GLDrawLine(pos1,pos1+(sc->ori1*sc->cp1rel*Vector3r::UNIT_X*sc->d1),Vector3r(0,.5,1));
- //GLUtils::GLDrawLine(pos2,pos2+(sc->ori2*sc->cp2rel*Vector3r::UNIT_X*sc->d2),Vector3r(0,1,.5));
- //TRVAR4(pos1,sc->ori1,pos2,sc->ori2);
- //TRVAR2(sc->cp2rel,pos2+(sc->ori2*sc->cp2rel*Vector3r::UNIT_X*sc->d2));
- // contact point to projected points
- #if 1
- Vector3r ptTg1=sc->contPtInTgPlane1(), ptTg2=sc->contPtInTgPlane2();
- //TRVAR3(ptTg1,ptTg2,sc->normal)
- GLUtils::GLDrawLine(contPt,contPt+ptTg1,Vector3r(0,.5,1)); GLUtils::GLDrawLine(pos1,contPt+ptTg1,Vector3r(0,.5,1));
- GLUtils::GLDrawLine(contPt,contPt+ptTg2,Vector3r(0,1,.5)); GLUtils::GLDrawLine(pos2,contPt+ptTg2,Vector3r(0,1,.5));
- // projected shear
- //GLUtils::GLDrawLine(contPt+ptTg1,contPt+ptTg2,Vector3r(1,1,1));
- //
- //GLUtils::GLDrawNum(sc->epsN(),contPt,Vector3r(1,1,1));
- #endif
- }
-
-
-
}
Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/pkg/dem/SConscript 2009-05-05 05:15:26 UTC (rev 1766)
@@ -33,6 +33,7 @@
env.SharedLibrary('FacetTopologyAnalyzer',['Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp'],LIBS=env['LIBS']+['InteractingFacet']),
+ env.SharedLibrary('Clump',['DataClass/Clump.cpp'],LIBS=env['LIBS']+['Shop']),
env.SharedLibrary('SQLiteRecorder',
['Engine/StandAloneEngine/SQLiteRecorder.cpp'],
@@ -949,7 +950,8 @@
LIBS=env['LIBS']+['yade-serialization',
'yade-base',
'yade-multimethods',
- 'RigidBodyParameters'
+ 'RigidBodyParameters',
+ 'Clump'
]),
env.SharedLibrary('SimpleViscoelasticBodyParameters'
@@ -1049,9 +1051,9 @@
,'RotationEngine'
]),
- env.SharedLibrary('ef2_Spheres_NormalShear_ElasticFrictionalLaw',
- ['Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.cpp'],
- LIBS=env['LIBS']+['ConstitutiveLawDispatcher','NormalShearInteractions','SpheresContactGeometry']),
+ env.SharedLibrary('ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw',
+ ['Engine/EngineUnit/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.cpp'],
+ LIBS=env['LIBS']+['ConstitutiveLawDispatcher','NormalShearInteractions','DemXDofGeom']),
env.SharedLibrary('ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw',
['Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp'],
Added: trunk/scripts/test/clump.py
===================================================================
--- trunk/scripts/test/clump.py 2009-05-04 14:31:07 UTC (rev 1765)
+++ trunk/scripts/test/clump.py 2009-05-05 05:15:26 UTC (rev 1766)
@@ -0,0 +1,33 @@
+O.initializers=[BoundingVolumeMetaEngine([InteractingSphere2AABB(),InteractingBox2AABB(),MetaInteractingGeometry2AABB()])]
+O.engines=[
+ PhysicalActionContainerReseter(),
+ BoundingVolumeMetaEngine([
+ InteractingSphere2AABB(),
+ InteractingBox2AABB(),
+ MetaInteractingGeometry2AABB()
+ ]),
+ PersistentSAPCollider(),
+ InteractionDispatchers(
+ [InteractingSphere2InteractingSphere4SpheresContactGeometry(),InteractingBox2InteractingSphere4SpheresContactGeometry()],
+ [SimpleElasticRelationships()],
+ [ef2_Spheres_Elastic_ElasticLaw()]
+ ),
+ GravityEngine(gravity=[0,0,-9.81]),
+ NewtonsDampedLaw(damping=.2)
+]
+# support
+O.bodies.append(utils.box([0,0,-1.5],[3,3,.2],dynamic=False))
+# stand-alone sphere
+O.bodies.append(utils.sphere([0,0,0],.5,density=1000,color=[1,1,1]))
+# clumps
+relPos=[(0,-.5,-.5),(0,.5,0),(.5,0,0),(0,0,.5)]
+coords=[(-2,0,0),(2,0,0),(0,2,0),(0,-2,0)]
+for i,cc in enumerate(coords):
+ # This shorthand command does something like this:
+ # O.bodies.appendClumped([utils.sphere(...),utils.sphere(...),utils.sphere(...)])
+ # and returns tuple of clumpId,[bodyId1,bodyId2,bodyId3]
+ clump,spheres=O.bodies.appendClumped([utils.sphere([relPos[j][0]+coords[i][0],relPos[j][1]+coords[i][1],relPos[j][2]+coords[i][2]],.5,density=1000) for j in range(0,i+1)])
+ print clump,spheres
+O.dt=utils.PWaveTimeStep()
+O.saveTmp('init')
+
