yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #00999
[svn] r1678 - in trunk: core core/DefaultContainerImplementations extra extra/clump gui/py lib/base pkg/common/Engine/MetaEngine pkg/dem/Engine/DeusExMachina pkg/dem/Engine/StandAloneEngine
Author: eudoxos
Date: 2009-02-22 11:30:11 +0100 (Sun, 22 Feb 2009)
New Revision: 1678
Modified:
trunk/core/DefaultContainerImplementations/PhysicalActionVectorVector.cpp
trunk/core/MetaBody.cpp
trunk/core/MetaBody.hpp
trunk/extra/Brefcom.cpp
trunk/extra/clump/Shop.cpp
trunk/gui/py/yadeControl.cpp
trunk/lib/base/yadeWm3Extra.hpp
trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp
trunk/pkg/common/Engine/MetaEngine/InteractionPhysicsMetaEngine.cpp
trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp
trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessCounter.cpp
trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp
trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.hpp
Log:
1. Adapted TriaxialTest and ElasticContactLaw to BexContainer (switchable at compile-time).
2. Adapter GlobalStiffnessTimeStepper to BexContainer (GlobalStiffnessCounter code put to a function inside this one, hence that engine is not needed anymore).
3. Adapted Shop::Bex to BexContainer (probably not needed anymore?)
4. Exception is thrown in PhysicalActionVectorVector is used && built with BEX_CONTAINER?\194?\160(would most likely crash anyway).
5. transientInteractions and persistentInteractions are only references to MetaBody::interactions now. Removed extra loops in InteractionPhysicsMetaEngine and InteractionGeometryMetaEngine.
6. Remove including qglviewer into miniWm3, as it breaks compilation if using miniWm3 separately from yade (for testing purposes).
Modified: trunk/core/DefaultContainerImplementations/PhysicalActionVectorVector.cpp
===================================================================
--- trunk/core/DefaultContainerImplementations/PhysicalActionVectorVector.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/core/DefaultContainerImplementations/PhysicalActionVectorVector.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -164,6 +164,10 @@
// should be always succesfull. if it is not - you forgot to call prepare()
shared_ptr<PhysicalAction>& PhysicalActionVectorVector::find(unsigned int id , int actionIndex )
{
+ #ifdef BEX_CONTAINER
+ cerr<<"FATAL: This build of yade uses nex BexContainer instead of PhysicalActionContainer.\nFATAL: However, your simulation still uses PhysicalActionContainer.\nFATAL: Update your code, see backtrace (if in debug build) to find where the old container is used."<<endl;
+ throw std::runtime_error("Deprecated PhysicalActionContainer is not supported in this build!");
+ #endif
if( current_size <= id ) // this is very rarely executed, only at beginning.
// somebody is accesing out of bounds, make sure he will find, what he needs - a resetted PhysicalAction of his type
{
Modified: trunk/core/MetaBody.cpp
===================================================================
--- trunk/core/MetaBody.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/core/MetaBody.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -30,7 +30,7 @@
#include<unistd.h>
MetaBody::MetaBody() :
- Body(),bodies(new BodyRedirectionVector),persistentInteractions(new InteractionVecSet),transientInteractions(new InteractionVecSet),physicalActions(new PhysicalActionVectorVector)
+ Body(),bodies(new BodyRedirectionVector), interactions(new InteractionVecSet), persistentInteractions(interactions),transientInteractions(interactions),physicalActions(new PhysicalActionVectorVector)
{
engines.clear();
initializers.clear();
Modified: trunk/core/MetaBody.hpp
===================================================================
--- trunk/core/MetaBody.hpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/core/MetaBody.hpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -23,8 +23,12 @@
shared_ptr<BodyContainer> bodies;
vector<shared_ptr<Engine> > engines;
vector<shared_ptr<Engine> > initializers; // FIXME: see MovingSupport:50
- __attribute__((__deprecated__)) shared_ptr<InteractionContainer> persistentInteractions; // disappear, reappear according to physical (or any other non-spatial) criterion
- shared_ptr<InteractionContainer> transientInteractions; // disappear, reappear according to spatial criterion
+ shared_ptr<InteractionContainer> interactions;
+
+ // only aliases for interactions
+ __attribute__((__deprecated__)) shared_ptr<InteractionContainer>& persistentInteractions; // disappear, reappear according to physical (or any other non-spatial) criterion
+ shared_ptr<InteractionContainer>& transientInteractions; // disappear, reappear according to spatial criterion
+
shared_ptr<PhysicalActionContainer> physicalActions;
#ifdef BEX_CONTAINER
BexContainer bex;
Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/extra/Brefcom.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -28,7 +28,7 @@
// commulate normal strains from contacts
// get max force on contacts
Real maxContactF=0;
- FOREACH(const shared_ptr<Interaction>& I, *rb->transientInteractions){
+ FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
if(!I->isReal) continue;
shared_ptr<BrefcomContact> BC=YADE_PTR_CAST<BrefcomContact>(I->interactionPhysics); assert(BC);
maxContactF=max(maxContactF,max(BC->Fn,BC->Fs.Length()));
@@ -39,7 +39,7 @@
}
unbalancedForce=(useMaxForce?maxF:meanF)/maxContactF;
- FOREACH(const shared_ptr<Interaction>& I, *rb->transientInteractions){
+ FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
if(!I->isReal) continue;
shared_ptr<BrefcomContact> BC=YADE_PTR_CAST<BrefcomContact>(I->interactionPhysics); assert(BC);
BrefcomPhysParams* bpp1(YADE_CAST<BrefcomPhysParams*>(Body::byId(I->getId1())->physicalParameters.get()));
@@ -183,14 +183,14 @@
if(useFunctor){ // testing the functor
if(!functor) functor=shared_ptr<ef2_Spheres_Brefcom_BrefcomLaw>(new ef2_Spheres_Brefcom_BrefcomLaw);
functor->logStrain=logStrain;
- FOREACH(const shared_ptr<Interaction>& I, *rootBody->transientInteractions){
+ FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){
if(!I->isReal) continue;
functor->go(I->interactionGeometry, I->interactionPhysics, I.get(), rootBody);
}
return;
}
- FOREACH(const shared_ptr<Interaction>& I, *rootBody->transientInteractions){
+ FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){
if(!I->isReal) continue;
BC=YADE_PTR_CAST<BrefcomContact>(I->interactionPhysics);
contGeom=YADE_PTR_CAST<SpheresContactGeometry>(I->interactionGeometry);
@@ -328,7 +328,7 @@
//vector<pair<short,
vector<BodyStats> bodyStats; bodyStats.resize(rootBody->bodies->size());
assert(bodyStats[0].nCohLinks==0); // should be initialized by dfault ctor
- FOREACH(const shared_ptr<Interaction>& I, *rootBody->transientInteractions){
+ FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){
shared_ptr<BrefcomContact> BC=dynamic_pointer_cast<BrefcomContact>(I->interactionPhysics);
if(!BC || !BC->isCohesive) continue;
const body_id_t id1=I->getId1(), id2=I->getId2();
Modified: trunk/extra/clump/Shop.cpp
===================================================================
--- trunk/extra/clump/Shop.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/extra/clump/Shop.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -96,14 +96,16 @@
Shop::Bex::globalStiffnessIdx=shared_ptr<PhysicalAction>(new GlobalStiffness())->getClassIndex();
}
}
+#ifdef BEX_CONTAINER
+ Vector3r& Shop::Bex::force(body_id_t id,MetaBody* rb){ return rb->bex.force(id);}
+ Vector3r& Shop::Bex::momentum(body_id_t id,MetaBody* rb){ return rb->bex.force(id);}
+#else
+ #define __BEX_ACCESS(retType,shopBexMember,bexClass,bexIdx,bexAttribute) retType& Shop::Bex::shopBexMember(body_id_t id,MetaBody* mb){ assert(bexIdx>=0); shared_ptr<PhysicalActionContainer> pac=(mb?mb:Omega::instance().getRootBody().get())->physicalActions; /*if((long)pac->size()<=id) throw invalid_argument("No " #shopBexMember " for #"+lexical_cast<string>(id)+" (max="+lexical_cast<string>(((long)pac->size()-1))+")");*/ return static_pointer_cast<bexClass>(pac->find(id,bexIdx))->bexAttribute; }
+ __BEX_ACCESS(Vector3r,force,Force,forceIdx,force);
+ __BEX_ACCESS(Vector3r,momentum,Momentum,momentumIdx,momentum);
+ #undef __BEX_ACCESS
+#endif
-#define __BEX_ACCESS(retType,shopBexMember,bexClass,bexIdx,bexAttribute) retType& Shop::Bex::shopBexMember(body_id_t id,MetaBody* mb){ assert(bexIdx>=0); shared_ptr<PhysicalActionContainer> pac=(mb?mb:Omega::instance().getRootBody().get())->physicalActions; /*if((long)pac->size()<=id) throw invalid_argument("No " #shopBexMember " for #"+lexical_cast<string>(id)+" (max="+lexical_cast<string>(((long)pac->size()-1))+")");*/ return static_pointer_cast<bexClass>(pac->find(id,bexIdx))->bexAttribute; }
-__BEX_ACCESS(Vector3r,force,Force,forceIdx,force);
-__BEX_ACCESS(Vector3r,momentum,Momentum,momentumIdx,momentum);
-__BEX_ACCESS(Vector3r,globalStiffness,GlobalStiffness,globalStiffnessIdx,stiffness);
-__BEX_ACCESS(Vector3r,globalRStiffness,GlobalStiffness,globalStiffnessIdx,Rstiffness);
-#undef __BEX_ACCESS
-
/* Apply force on contact point to 2 bodies; the force is oriented as it applies on the first body and is reversed on the second.
*
* Shop::Bex::initCache must have been called at some point. */
Modified: trunk/gui/py/yadeControl.cpp
===================================================================
--- trunk/gui/py/yadeControl.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/gui/py/yadeControl.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -522,7 +522,7 @@
pyBodyContainer bodies_get(void){assertRootBody(); return pyBodyContainer(OMEGA.getRootBody()->bodies); }
- pyInteractionContainer interactions_get(void){assertRootBody(); return pyInteractionContainer(OMEGA.getRootBody()->transientInteractions); }
+ pyInteractionContainer interactions_get(void){assertRootBody(); return pyInteractionContainer(OMEGA.getRootBody()->interactions); }
#ifdef BEX_CONTAINER
pyBexContainer actions_get(void){return pyBexContainer();}
@@ -541,11 +541,11 @@
void interactionContainer_set(string clss){
MetaBody* rb=OMEGA.getRootBody().get();
- if(rb->transientInteractions->size()>0) throw std::runtime_error("Interaction container not empty, will not change its class.");
+ if(rb->interactions->size()>0) throw std::runtime_error("Interaction container not empty, will not change its class.");
shared_ptr<InteractionContainer> ic=dynamic_pointer_cast<InteractionContainer>(ClassFactory::instance().createShared(clss));
- rb->transientInteractions=ic;
+ rb->interactions=ic;
}
- string interactionContainer_get(string clss){ return OMEGA.getRootBody()->transientInteractions->getClassName(); }
+ string interactionContainer_get(string clss){ return OMEGA.getRootBody()->interactions->getClassName(); }
void bodyContainer_set(string clss){
MetaBody* rb=OMEGA.getRootBody().get();
Modified: trunk/lib/base/yadeWm3Extra.hpp
===================================================================
--- trunk/lib/base/yadeWm3Extra.hpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/lib/base/yadeWm3Extra.hpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -143,9 +143,8 @@
} // namespace serialization
} // namespace boost
-
+#if 0
#include<yade/lib-QGLViewer/qglviewer.h>
-
inline qglviewer::Vec toQGLViewierVec(Vector3r v){return qglviewer::Vec(v[0],v[1],v[2]);};
inline Vector3r toVec(qglviewer::Vec v){return Vector3r(v[0],v[1],v[2]);};
-
+#endif
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -5,7 +5,7 @@
* cosurgi@xxxxxxxxxx *
* *
* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
+* GNUGeneral Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
#include "InteractionGeometryMetaEngine.hpp"
@@ -43,32 +43,12 @@
void InteractionGeometryMetaEngine::action(MetaBody* ncb)
{
shared_ptr<BodyContainer>& bodies = ncb->bodies;
- FOREACH(const shared_ptr<Interaction>& interaction, *ncb->persistentInteractions){
- shared_ptr<Body>& b1=(*bodies)[interaction->getId1()];
- shared_ptr<Body>& b2=(*bodies)[interaction->getId2()];
- interaction->isReal = true;
- operator()( b1->interactingGeometry , b2->interactingGeometry , b1->physicalParameters->se3 , b2->physicalParameters->se3 , interaction );
- }
-
- FOREACH(const shared_ptr<Interaction>& interaction, *ncb->transientInteractions){
+ FOREACH(const shared_ptr<Interaction>& interaction, *ncb->interactions){
const shared_ptr<Body>& b1=(*bodies)[interaction->getId1()];
const shared_ptr<Body>& b2=(*bodies)[interaction->getId2()];
- //bool wasReal = interaction->isReal;
interaction->isReal =
b1->interactingGeometry && b2->interactingGeometry && // some bodies do not have interactingGeometry
- ( ncb->persistentInteractions->find(interaction->getId1(),interaction->getId2()) == 0 )
- &&
- operator()( b1->interactingGeometry , b2->interactingGeometry , b1->physicalParameters->se3 , b2->physicalParameters->se3 , interaction );
-
- //if(wasReal==false && interaction->isReal)
- // interaction->isNew=true;
- //cerr<<"isReal="<<interaction->isReal<<", wasReal="<<wasReal<<", isNew="<<interaction->isNew<<endl;
-
- //tmp
- //if(!(b1->interactingGeometry&&b2->interactingGeometry)){
- //cerr<<__FILE__<<":"<<__LINE__<<": no interacting geometry "<< (b1->interactingGeometry?b1->getId():-1)<<" "<<(b2->interactingGeometry?b2->getId():-1)<<endl;
- //}
-
+ operator()(b1->interactingGeometry, b2->interactingGeometry, b1->physicalParameters->se3, b2->physicalParameters->se3, interaction);
}
}
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionPhysicsMetaEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionPhysicsMetaEngine.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionPhysicsMetaEngine.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -24,32 +24,11 @@
void InteractionPhysicsMetaEngine::action(MetaBody* ncb)
{
shared_ptr<BodyContainer>& bodies = ncb->bodies;
-
- shared_ptr<InteractionContainer>& persistentInteractions = ncb->persistentInteractions;
- InteractionContainer::iterator ii = persistentInteractions->begin();
- InteractionContainer::iterator iiEnd = persistentInteractions->end();
- for( ; ii!=iiEnd ; ++ii)
- {
- const shared_ptr<Interaction> interaction = *ii;
-
+ FOREACH(const shared_ptr<Interaction>& interaction, *ncb->interactions){
shared_ptr<Body>& b1 = (*bodies)[interaction->getId1()];
shared_ptr<Body>& b2 = (*bodies)[interaction->getId2()];
- if( b1->physicalParameters && b2->physicalParameters )
- operator()( b1->physicalParameters , b2->physicalParameters , interaction );
- }
-
- shared_ptr<InteractionContainer>& transientInteractions = ncb->transientInteractions;
- ii = transientInteractions->begin();
- iiEnd = transientInteractions->end();
- for( ; ii!=iiEnd ; ++ii)
- {
- const shared_ptr<Interaction> interaction = *ii;
-
- shared_ptr<Body>& b1 = (*bodies)[interaction->getId1()];
- shared_ptr<Body>& b2 = (*bodies)[interaction->getId2()];
-
if (interaction->isReal)
- operator()( b1->physicalParameters , b2->physicalParameters , interaction );
+ operator()(b1->physicalParameters, b2->physicalParameters, interaction);
}
}
Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -168,11 +168,10 @@
void TriaxialStressController::controlExternalStress(int wall, MetaBody* ncb, Vector3r resultantForce, PhysicalParameters* p, Real wall_max_vel)
{
- Real translation=normal[wall].Dot(static_cast<Force*>( ncb->physicalActions->find(wall_id[wall],ForceClassIndex).get() )->force - resultantForce);
+ Real translation=normal[wall].Dot( getForce(ncb,wall_id[wall]) - resultantForce);
//bool log=((wall==3) && (Omega::instance().getCurrentIteration()%200==0));
const bool log=false;
- if(log) LOG_DEBUG("wall="<<wall<<" actualForce="<<static_cast<Force*>(ncb->physicalActions->find(wall_id[wall],ForceClassIndex).get() )->force<<", resultantForce="<<resultantForce<<", translation="<<translation);
- //cerr << "wall="<<wall<<" actualForce="<<static_cast<Force*>(ncb->physicalActions->find(wall_id[wall],ForceClassIndex).get() )->force<<", resultantForce="<<resultantForce<<", deltaF="<<translation << endl;
+ if(log) LOG_DEBUG("wall="<<wall<<" actualForce="<<getForce(ncb,wall_id[wall])<<", resultantForce="<<resultantForce<<", translation="<<translation);
if (translation!=0)
{
//cerr << "stiffness = " << stiffness[wall];
@@ -318,7 +317,16 @@
Real invXSurface = 1.f/(height*depth);
Real invYSurface = 1.f/(width*depth);
Real invZSurface = 1.f/(width*height);
-
+
+ force[wall_bottom]=getForce(ncb,wall_id[wall_bottom]); stress[wall_bottom]=force[wall_bottom]*invYSurface;
+ force[wall_top]= getForce(ncb,wall_id[wall_top]); stress[wall_top]= force[wall_top] *invYSurface;
+ force[wall_left]= getForce(ncb,wall_id[wall_left]); stress[wall_left]= force[wall_left] *invXSurface;
+ force[wall_right]= getForce(ncb,wall_id[wall_right]); stress[wall_right]= force[wall_right] *invXSurface;
+ force[wall_front]= getForce(ncb,wall_id[wall_front]); stress[wall_front]= force[wall_front] *invZSurface;
+ force[wall_back]= getForce(ncb,wall_id[wall_back]); stress[wall_back]= force[wall_back] *invZSurface;
+
+// remove this, it is in the 6 lines above now
+#if 0
stress[wall_bottom] = ( static_cast<Force*>( ncb->physicalActions->find(wall_id[wall_bottom],ForceClassIndex).get() )->force ) * invYSurface;
stress[wall_top] = static_cast<Force*>( ncb->physicalActions->find(wall_id[wall_top],ForceClassIndex).get() )->force * invYSurface;
stress[wall_left] = ( static_cast<Force*>( ncb->physicalActions->find(wall_id[wall_left],ForceClassIndex).get() )->force ) * invXSurface;
@@ -326,8 +334,6 @@
stress[wall_front] = ( static_cast<Force*>( ncb->physicalActions->find(wall_id[wall_front],ForceClassIndex).get() )->force ) * invZSurface;
stress[wall_back] = ( static_cast<Force*>( ncb->physicalActions->find(wall_id[wall_back],ForceClassIndex).get() )->force ) * invZSurface;
-
-
// ==================================================== jf
force[wall_bottom] = ( static_cast<Force*>( ncb->physicalActions->find(wall_id[wall_bottom],ForceClassIndex).get() )->force );
force[wall_top] = static_cast<Force*>( ncb->physicalActions->find(wall_id[wall_top],ForceClassIndex).get() )->force ;
@@ -336,8 +342,8 @@
force[wall_front] = ( static_cast<Force*>( ncb->physicalActions->find(wall_id[wall_front],ForceClassIndex).get() )->force ) ;
force[wall_back] = ( static_cast<Force*>( ncb->physicalActions->find(wall_id[wall_back],ForceClassIndex).get() )->force ) ;
// ====================================================
+#endif
-
//cerr << "stresses : " << stress[wall_bottom] << " " <<
//stress[wall_top]<< " " << stress[wall_left]<< " " << stress[wall_right]<< " "
//<< stress[wall_front] << " " << stress[wall_back] << endl;
@@ -422,7 +428,7 @@
Real f;
for( ; bi!=biEnd ; ++bi ) {
if ((*bi)->isDynamic) {
- f= (static_cast<Force*> ( ncb->physicalActions->find( (*bi)->getId() , ForceClassIndex).get() )->force).Length();
+ f=getForce(ncb,(*bi)->getId()).Length();
MeanUnbalanced += f;
if (f!=0) ++nBodies;
}
@@ -436,7 +442,7 @@
BodyContainer::iterator biEnd = bodies->end();
for( ; bi!=biEnd ; ++bi ) {
if ((*bi)->isDynamic) {
- MaxUnbalanced = std::max((static_cast<Force*> ( ncb->physicalActions->find( (*bi)->getId() , ForceClassIndex).get() )->force).Length(), MaxUnbalanced);
+ MaxUnbalanced = std::max(getForce(ncb,(*bi)->getId()).Length(),MaxUnbalanced);
}
}
if (MeanForce != 0) MaxUnbalanced = MaxUnbalanced/MeanForce;
Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -9,7 +9,7 @@
#pragma once
#include<yade/core/DeusExMachina.hpp>
-#include <Wm3Vector3.h>
+#include<yade/core/MetaBody.hpp>
#include<yade/lib-base/yadeWm3.hpp>
#define TR {if (Omega::instance().getCurrentIteration()%100==0) TRACE; }
@@ -31,6 +31,13 @@
int ForceClassIndex;
Real previousStress, previousMultiplier; //previous mean stress and size multiplier
bool firstRun;
+ inline const Vector3r getForce(MetaBody* rb, body_id_t id){
+ #ifdef BEX_CONTAINER
+ return rb->bex.force(id);
+ #else
+ return static_cast<Force*>(rb->physicalActions->find(id,ForceClassIndex).get())->force
+ #endif
+ }
public :
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -169,14 +169,18 @@
////////// PFC3d SlipModel
Vector3r f = currentContactPhysics->normalForce + shearForce;
+ #ifdef BEX_CONTAINER
+ ncb->bex.force(id1)-=f;
+ ncb->bex.force(id2)+=f;
+ ncb->bex.torque(id1)-=c1x.Cross(f);
+ ncb->bex.torque(id2)+=c2x.Cross(f);
+ #else
+ static_cast<Force*> ( ncb->physicalActions->find( id1 , actionForceIndex).get() )->force -= f;
+ static_cast<Force*> ( ncb->physicalActions->find( id2 , actionForceIndex ).get() )->force += f;
+ static_cast<Momentum*>( ncb->physicalActions->find( id1 , actionMomentumIndex ).get() )->momentum -= c1x.Cross(f);
+ static_cast<Momentum*>( ncb->physicalActions->find( id2 , actionMomentumIndex ).get() )->momentum += c2x.Cross(f);
+ #endif
- // it will be some macro( body->physicalActions, ActionType , bodyId )
- static_cast<Force*> ( ncb->physicalActions->find( id1 , actionForceIndex).get() )->force -= f;
- static_cast<Force*> ( ncb->physicalActions->find( id2 , actionForceIndex ).get() )->force += f;
-
- static_cast<Momentum*>( ncb->physicalActions->find( id1 , actionMomentumIndex ).get() )->momentum -= c1x.Cross(f);
- static_cast<Momentum*>( ncb->physicalActions->find( id2 , actionMomentumIndex ).get() )->momentum += c2x.Cross(f);
-
currentContactPhysics->prevNormal = currentContactGeometry->normal;
}
}
Modified: trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessCounter.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessCounter.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessCounter.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -54,6 +54,9 @@
#endif
void GlobalStiffnessCounter::traverseInteractions(MetaBody* ncb, const shared_ptr<InteractionContainer>& interactions){
+ #ifdef BEX_CONTAINER
+ return;
+ #endif
FOREACH(const shared_ptr<Interaction>& contact, *interactions){
if(!contact->isReal) continue;
Modified: trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -17,6 +17,7 @@
#include<yade/pkg-dem/GlobalStiffness.hpp>
CREATE_LOGGER(GlobalStiffnessTimeStepper);
+YADE_PLUGIN("GlobalStiffnessTimeStepper");
GlobalStiffnessTimeStepper::GlobalStiffnessTimeStepper() : TimeStepper() , sdecContactModel(new MacroMicroElasticRelationships), actionParameterGlobalStiffness(new GlobalStiffness)
{
@@ -54,10 +55,14 @@
// Sphere* sphere = static_cast<Sphere*>(body->geometricalModel.get());
- Vector3r& stiffness = (static_cast<GlobalStiffness*>( ncb->physicalActions->find (body->getId(), globalStiffnessClassIndex).get()))->stiffness;
- Vector3r& Rstiffness = (static_cast<GlobalStiffness*>( ncb->physicalActions->find (body->getId(), globalStiffnessClassIndex).get()))->Rstiffness;
+ #ifdef BEX_CONTAINER
+ Vector3r& stiffness= stiffnesses[body->getId()];
+ Vector3r& Rstiffness=Rstiffnesses[body->getId()];
+ #else
+ Vector3r& stiffness = (static_cast<GlobalStiffness*>( ncb->physicalActions->find (body->getId(), globalStiffnessClassIndex).get()))->stiffness;
+ Vector3r& Rstiffness = (static_cast<GlobalStiffness*>( ncb->physicalActions->find (body->getId(), globalStiffnessClassIndex).get()))->Rstiffness;
+ #endif
-
//cerr << "Vector3r& stiffness = (static_cast<GlobalStiffness*>( ncb" << endl;
if(! ( /* sphere && */ sdec && stiffness) )
return; // not possible to compute!
@@ -133,6 +138,10 @@
void GlobalStiffnessTimeStepper::computeTimeStep(MetaBody* ncb)
{
+ #ifdef BEX_CONTAINER
+ computeStiffnesses(ncb);
+ #endif
+
shared_ptr<BodyContainer>& bodies = ncb->bodies;
// shared_ptr<InteractionContainer>& transientInteractions = ncb->transientInteractions;
@@ -180,4 +189,53 @@
string(", BUT timestep is ")+lexical_cast<string>(Omega::instance().getTimeStep()))<<".");
}
-YADE_PLUGIN();
+#ifdef BEX_CONTAINER
+ void GlobalStiffnessTimeStepper::computeStiffnesses(MetaBody* rb){
+ /* check size */
+ size_t size=stiffnesses.size();
+ if(size<rb->bodies->size()){
+ size=rb->bodies->size();
+ stiffnesses.resize(size); Rstiffnesses.resize(size);
+ }
+ /* reset stored values */
+ memset(stiffnesses[0], 0,sizeof(Vector3r)*size);
+ memset(Rstiffnesses[0],0,sizeof(Vector3r)*size);
+ /* loop copied verbatim from GlobalStiffnessCounter */
+ FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
+ if(!contact->isReal) continue;
+
+ SpheresContactGeometry* geom=YADE_CAST<SpheresContactGeometry*>(contact->interactionGeometry.get()); assert(geom);
+ NormalShearInteraction* phys=YADE_CAST<NormalShearInteraction*>(contact->interactionPhysics.get()); assert(phys);
+ // all we need for getting stiffness
+ Vector3r& normal=geom->normal; Real& kn=phys->kn; Real& ks=phys->ks; Real& radius1=geom->radius1; Real& radius2=geom->radius2;
+ // FIXME? NormalShearInteraction knows nothing about whether the contact is "active" (force!=0) or not;
+ // former code: if(force==0) continue; /* disregard this interaction, it is not active */.
+ // It seems though that in such case either the interaction is accidentally at perfect equilibrium (unlikely)
+ // or it should have been deleted already. Right?
+ //ANSWER : some interactions can exist without fn, e.g. distant capillary force, wich does not contribute to the overall stiffness via kn. The test is needed.
+ Real fn = (static_cast<NormalShearInteraction *> (contact->interactionPhysics.get()))->normalForce.SquaredLength();
+
+ if (fn!=0) {
+ //Diagonal terms of the translational stiffness matrix
+ Vector3r diag_stiffness = Vector3r(std::pow(normal.X(),2),std::pow(normal.Y(),2),std::pow(normal.Z(),2));
+ diag_stiffness *= kn-ks;
+ diag_stiffness = diag_stiffness + Vector3r(1,1,1)*ks;
+
+ //diagonal terms of the rotational stiffness matrix
+ // Vector3r branch1 = currentContactGeometry->normal*currentContactGeometry->radius1;
+ // Vector3r branch2 = currentContactGeometry->normal*currentContactGeometry->radius2;
+ Vector3r diag_Rstiffness =
+ Vector3r(std::pow(normal.Y(),2)+std::pow(normal.Z(),2),
+ std::pow(normal.X(),2)+std::pow(normal.Z(),2),
+ std::pow(normal.X(),2)+std::pow(normal.Y(),2));
+ diag_Rstiffness *= ks;
+ //cerr << "diag_Rstifness=" << diag_Rstiffness << endl;
+
+ stiffnesses [contact->getId1()]+=diag_stiffness;
+ Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2);
+ stiffnesses [contact->getId2()]+=diag_stiffness;
+ Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
+ }
+ }
+ }
+#endif
Modified: trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.hpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.hpp 2009-02-21 20:14:55 UTC (rev 1677)
+++ trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.hpp 2009-02-22 10:30:11 UTC (rev 1678)
@@ -26,6 +26,12 @@
class GlobalStiffnessTimeStepper : public TimeStepper
{
private :
+ #ifdef BEX_CONTAINER
+ vector<Vector3r> stiffnesses;
+ vector<Vector3r> Rstiffnesses;
+ void computeStiffnesses(MetaBody*); // what GlobalStiffnessCounter used to do
+ #endif
+
Real newDt, previousDt;
bool computedSomething,
computedOnce;
