yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01282
[svn] r1786 - in trunk: . core examples examples/collider-perf examples/concrete examples/concrete/pack extra extra/tetra gui/py lib/import pkg/common/Engine/DeusExMachina pkg/common/Engine/EngineUnit pkg/common/Engine/MetaEngine pkg/common/Engine/StandAloneEngine pkg/dem pkg/dem/DataClass/InteractionGeometry pkg/dem/Engine/DeusExMachina pkg/dem/Engine/EngineUnit pkg/dem/Engine/StandAloneEngine pkg/dem/PreProcessor pkg/dem/RenderingEngine/GLDrawCohesiveFrictionalContactInteraction pkg/dem/RenderingEngine/GLDrawElasticContactInteraction pkg/dem/RenderingEngine/GLDrawSDECLinkGeometry pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry pkg/lattice/Engine/StandAloneEngine pkg/lattice/PreProcessor pkg/mass-spring/Engine/StandAloneEngine pkg/mass-spring/PreProcessor pkg/realtime-rigidbody/Engine/StandAloneEngine pkg/snow/Engine pkg/snow/RenderingEngine scripts/test
Author: eudoxos
Date: 2009-05-29 16:35:22 +0200 (Fri, 29 May 2009)
New Revision: 1786
Added:
trunk/examples/concrete/
trunk/examples/concrete/pack/
trunk/examples/concrete/pack/mk-pack.py
Modified:
trunk/SConstruct
trunk/core/Collider.cpp
trunk/core/Interaction.cpp
trunk/core/Interaction.hpp
trunk/core/InteractionContainer.cpp
trunk/core/InteractionContainer.hpp
trunk/core/MetaEngine2D.hpp
trunk/core/yade.cpp
trunk/examples/collider-perf/perf.table
trunk/extra/Shop.cpp
trunk/extra/tetra/Tetra.cpp
trunk/gui/py/_eudoxos.cpp
trunk/gui/py/_utils.cpp
trunk/gui/py/utils.py
trunk/gui/py/yadeControl.cpp
trunk/lib/import/STLReader.hpp
trunk/pkg/common/Engine/DeusExMachina/CinemCNCEngine.cpp
trunk/pkg/common/Engine/DeusExMachina/CinemKNCEngine.cpp
trunk/pkg/common/Engine/DeusExMachina/TranslationEngine.cpp
trunk/pkg/common/Engine/EngineUnit/ElasticBodySimpleRelationship.cpp
trunk/pkg/common/Engine/MetaEngine/ConstitutiveLawDispatcher.cpp
trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp
trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.hpp
trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp
trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.hpp
trunk/pkg/common/Engine/MetaEngine/InteractionPhysicsMetaEngine.cpp
trunk/pkg/common/Engine/StandAloneEngine/DistantPersistentSAPCollider.cpp
trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp
trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp
trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.hpp
trunk/pkg/common/Engine/StandAloneEngine/PersistentTriangulationCollider.cpp
trunk/pkg/common/Engine/StandAloneEngine/SpatialQuickSortCollider.cpp
trunk/pkg/common/Engine/StandAloneEngine/SpheresFactory.cpp
trunk/pkg/dem/ConcretePM.cpp
trunk/pkg/dem/ConcretePM.hpp
trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp
trunk/pkg/dem/Engine/DeusExMachina/CapillaryRecorder.cpp
trunk/pkg/dem/Engine/DeusExMachina/CapillaryStressRecorder.cpp
trunk/pkg/dem/Engine/DeusExMachina/ContactStressRecorder.cpp
trunk/pkg/dem/Engine/DeusExMachina/ThreeDTriaxialEngine.cpp
trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.cpp
trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.hpp
trunk/pkg/dem/Engine/DeusExMachina/TriaxialStateRecorder.cpp
trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp
trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp
trunk/pkg/dem/Engine/EngineUnit/CL1Relationships.cpp
trunk/pkg/dem/Engine/EngineUnit/CohesiveFrictionalRelationships.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometryWater.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingBox4InteractionOfMyTetrahedron.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingMyTetrahedron4InteractionOfMyTetrahedron.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2BssSweptSphereLineSegment4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/EngineUnit/MacroMicroElasticRelationships.cpp
trunk/pkg/dem/Engine/EngineUnit/MacroMicroElasticRelationshipsWater.cpp
trunk/pkg/dem/Engine/EngineUnit/SimpleElasticRelationships.cpp
trunk/pkg/dem/Engine/EngineUnit/SimpleElasticRelationshipsWater.cpp
trunk/pkg/dem/Engine/EngineUnit/SimpleViscoelasticRelationships.cpp
trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/CohesiveFrictionalContactLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/ContactLaw1.cpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticCohesiveLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/ForceSnapshot.cpp
trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp
trunk/pkg/dem/Engine/StandAloneEngine/MicroMacroAnalyser.cpp
trunk/pkg/dem/Engine/StandAloneEngine/MyTetrahedronLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticSpheresInteractionRecorder.cpp
trunk/pkg/dem/Engine/StandAloneEngine/ViscousForceDamping.cpp
trunk/pkg/dem/Engine/StandAloneEngine/VolumicContactLaw.cpp
trunk/pkg/dem/PreProcessor/HydraulicTest.cpp
trunk/pkg/dem/PreProcessor/SDECLinkedSpheres.cpp
trunk/pkg/dem/PreProcessor/ThreePointBending.cpp
trunk/pkg/dem/PreProcessor/TriaxialTest.cpp
trunk/pkg/dem/PreProcessor/TriaxialTest.hpp
trunk/pkg/dem/RenderingEngine/GLDrawCohesiveFrictionalContactInteraction/GLDrawCohesiveFrictionalContactInteraction.cpp
trunk/pkg/dem/RenderingEngine/GLDrawElasticContactInteraction/GLDrawElasticContactInteraction.cpp
trunk/pkg/dem/RenderingEngine/GLDrawSDECLinkGeometry/GLDrawSDECLinkGeometry.cpp
trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp
trunk/pkg/dem/SConscript
trunk/pkg/lattice/Engine/StandAloneEngine/NonLocalInitializer.cpp
trunk/pkg/lattice/PreProcessor/LatticeExample.cpp
trunk/pkg/lattice/PreProcessor/LatticeExampleCTData.cpp
trunk/pkg/lattice/PreProcessor/LatticeExampleSimple.cpp
trunk/pkg/mass-spring/Engine/StandAloneEngine/MassSpringLaw.cpp
trunk/pkg/mass-spring/PreProcessor/HangingCloth.cpp
trunk/pkg/realtime-rigidbody/Engine/StandAloneEngine/FrictionLessElasticContactLaw.cpp
trunk/pkg/snow/Engine/Ef2_BssSnowGrain_BssSnowGrain_makeSpheresContactGeometry.cpp
trunk/pkg/snow/Engine/ElawSnowLayersDeformation.cpp
trunk/pkg/snow/RenderingEngine/Ef1_IstSnowLayersContact_glDraw.cpp
trunk/scripts/test/regular-sphere-pack.py
Log:
1. Remove Interaction::isNew, Interaction::isReal. Add Interaction::isReal() which test for presence of _both_ interactionGeometry and interactionPhysics and Interaction::isFresh() which checks the interaction creation timestamp against current iteration. Updated all code for that. Please verify that your stuff works, it is possible I made some mistakes!
2. All code saying isReal=false replaced by interaction->requestErase(id1,id2)
3. make requestErase thread-safe
4. Add interace for traversing interactions pending erase (template)
5. Add TriaxialTest::noFiles and TriaxialCompressionEngine::noFiles to not generate any files (default: off)
6. Add predicate utils.ptInAABB(p,minPt,maxPt)
7. utils.spheresToFile return number of spheres written (instead of None)
8. Add examples/concrete/pack to generate packing; other stuff in concrete will follow, with simplified version of the cpm model in ConcretePM.
Modified: trunk/SConstruct
===================================================================
--- trunk/SConstruct 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/SConstruct 2009-05-29 14:35:22 UTC (rev 1786)
@@ -129,7 +129,7 @@
ListVariable('exclude','Yade components that will not be built','none',names=['qt3','gui','extra','common','dem','fem','lattice','mass-spring','realtime-rigidbody','snow']),
EnumVariable('arcs','Whether to generate or use branch probabilities','',['','gen','use'],{'no':'','0':'','false':''},1),
# OK, dummy prevents bug in scons: if one selects all, it says all in scons.config, but without quotes, which generates error.
- ListVariable('features','Optional features that are turned on','python,log4cxx',names=['python','log4cxx','binfmt','dummy']),
+ ListVariable('features','Optional features that are turned on','python,log4cxx,gl',names=['gl','python','log4cxx','binfmt','CGAL','dummy']),
('jobs','Number of jobs to run at the same time (same as -j, but saved)',4,None,int),
('extraModules', 'Extra directories with their own SConscript files (must be in-tree) (whitespace separated)',None,None,Split),
('buildPrefix','Where to create build-[version][variant] directory for intermediary files','..'),
@@ -334,6 +334,9 @@
and conf.CheckCXXHeader(['Python.h','numpy/ndarrayobject.h'],'<>'))
if not ok: featureNotOK('python')
env.Append(CPPDEFINES=['EMBED_PYTHON'])
+ if 'CGAL' in env['features']:
+ ok=cong.CheckLibWithHeader('CGAL','CGAL/Exact_predicates_inexact_constructions_kernel.h','c++','CGAL::Exact_predicates_inexact_constructions_kernel::Point_3();')
+ if not ok: featureNotOK('CGAL')
if env['useMiniWm3']: env.Append(LIBS='miniWm3',CPPDEFINES=['MINIWM3'])
env=conf.Finish()
Modified: trunk/core/Collider.cpp
===================================================================
--- trunk/core/Collider.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/core/Collider.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -11,8 +11,9 @@
Collider::Collider(){}
Collider::~Collider(){}
-
bool Collider::handleExistingInteraction(Interaction* I){
+ throw runtime_error("handleExistingInteraction should not be called");
+#if 0
/* logically, we have 4 possibilities
* 1. real new → ¬new, keep
* 2. real ¬new → keep (same as 1.)
@@ -34,8 +35,8 @@
assert(false); // unreachable
return false;
+#endif
}
-
bool Collider::mayCollide(const Body* b1, const Body* b2){
return
// not yet implemented: only collide if at least one of the bodies is not shadow
Modified: trunk/core/Interaction.cpp
===================================================================
--- trunk/core/Interaction.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/core/Interaction.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -8,23 +8,27 @@
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
-#include "Interaction.hpp"
+#include"Interaction.hpp"
+#include<yade/core/MetaBody.hpp>
+
Interaction::Interaction(): id1(0), id2(0){ reset(); }
Interaction::Interaction(body_id_t newId1,body_id_t newId2): id1(newId1), id2(newId2){ reset(); }
+bool Interaction::isFresh(MetaBody* rb){ return iterMadeReal==rb->currentIteration;}
+
+
void Interaction::reset(){
- isNew=true;
- isReal=false;
isNeighbor = true;//NOTE : TriangulationCollider needs that
+ iterMadeReal=-1;
functorCache.geomExists=true;
//functorCache.geom=shared_ptr<InteractionGeometryEngineUnit>(); functorCache.phys=shared_ptr<InteractionPhysicsEngineUnit>(); functorCache.constLaw=shared_ptr<ConstitutiveLaw>();
}
void Interaction::swapOrder(){
- if(interactionGeometry || interactionPhysics || !isNew){
- throw std::logic_error("Bodies in interaction cannot be swapped if !isNew, have interactionGeometry or have interactionPhysics.");
+ if(interactionGeometry || interactionPhysics){
+ throw std::logic_error("Bodies in interaction cannot be swapped if they have interactionGeometry or interactionPhysics.");
}
std::swap(id1,id2);
}
Modified: trunk/core/Interaction.hpp
===================================================================
--- trunk/core/Interaction.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/core/Interaction.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -17,20 +17,26 @@
class InteractionGeometryEngineUnit;
class InteractionPhysicsEngineUnit;
class ConstitutiveLaw;
+class MetaBody;
class Interaction : public Serializable
{
private :
body_id_t id1,id2;
+ //! Step number at which the interaction was fully created (interactionGeometry and interactionPhysics).
+ //! Should be touched only by InteractionPhysicsMetaEngine and InteractionDispatchers, making them friends therefore
+ long iterMadeReal;
+ friend class InteractionPhysicsMetaEngine;
+ friend class InteractionDispatchers;
public :
- // FIXME : test if InteractionPhysics==0 and remove this flag; we can also remove this flag, if we make another container for PotetntialInteraction with only ids
- bool isNew;
- // maybe we can remove this, and check if InteractingGeometry, and InteractionPhysics are empty?
- bool isReal;
+ bool isReal() const {return (bool)interactionGeometry && (bool)interactionPhysics;}
+ //! If this interaction was just created in this step (for the constitutive law, to know that it is the first time there)
+ bool isFresh(MetaBody* rb);
+
//! phase flag to mark (for example, SpatialQuickSortCollider mark by it the stale interactions)
bool cycle;
//! NOTE : TriangulationCollider needs this (nothing else)
- bool isNeighbor;
+ bool isNeighbor;
shared_ptr<InteractionGeometry> interactionGeometry;
shared_ptr<InteractionPhysics> interactionPhysics;
@@ -61,8 +67,7 @@
REGISTER_ATTRIBUTES(/*no base*/,
(id1)
(id2)
- (isNew)
- (isReal)
+ (iterMadeReal)
(interactionGeometry)
(interactionPhysics)
);
Modified: trunk/core/InteractionContainer.cpp
===================================================================
--- trunk/core/InteractionContainer.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/core/InteractionContainer.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -9,11 +9,18 @@
*************************************************************************/
#include "InteractionContainer.hpp"
-#include "Interaction.hpp"
+void InteractionContainer::requestErase(body_id_t id1, body_id_t id2){
+ find(id1,id2)->reset(); bodyIdPair v(0,2); v.push_back(id1); v.push_back(id2);
+ boost::mutex::scoped_lock lock(pendingEraseMutex);
+ pendingErase.push_back(v);
+}
+void InteractionContainer::unconditionalErasePending(){
+ FOREACH(const bodyIdPair& p, pendingErase){ erase(p[0],p[1]); }
+ pendingErase.clear();
+}
-void InteractionContainer::requestErase(body_id_t id1, body_id_t id2){ find(id1,id2)->reset(); bodyIdPair v(0,2); v.push_back(id1); v.push_back(id2); pendingErase.push_back(v); }
void InteractionContainer::preProcessAttributes(bool deserializing)
{
Modified: trunk/core/InteractionContainer.hpp
===================================================================
--- trunk/core/InteractionContainer.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/core/InteractionContainer.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -15,8 +15,11 @@
#include<iostream>
#include<boost/range.hpp>
+// BOOST_FOREACH compatibility
+#ifndef FOREACH
+# define FOREACH BOOST_FOREACH
+#endif
-
#include<yade/core/Interaction.hpp>
using namespace boost;
@@ -81,7 +84,7 @@
class InteractionContainer : public Serializable
{
public :
- boost::mutex drawloopmutex; // FIXME a hack, containers have to be rewritten lock-free.
+ boost::mutex drawloopmutex;
InteractionContainer() { };
virtual ~InteractionContainer() {};
@@ -103,11 +106,42 @@
// std::pair is not handle by yade::serialization, use vector<body_id_t> instead
typedef vector<body_id_t> bodyIdPair;
- // Ask for erasing the interaction given (from the constitutive law); this resets the interaction (to the initial=potential state)
- // and collider should traverse pendingErase to decide whether to delete the interaction completely or keep it potential
+ //! Ask for erasing the interaction given (from the constitutive law); this resets the interaction (to the initial=potential state)
+ //! and collider should traverse pendingErase to decide whether to delete the interaction completely or keep it potential
void requestErase(body_id_t id1, body_id_t id2);
+ /*! List of pairs of interactions that will be (maybe) erased by the collider;
+
+ If accessed from within a parallel section, pendingEraseMutex must be locked (this is done inside requestErase for you)
+ If there is, at one point, a multi-threaded collider, pendingEraseMutex should be moved to the public part and used from there as well.
+ */
list<bodyIdPair> pendingErase;
+ /*! Erase all pending interactions unconditionally.
+
+ This should be called only in rare cases that collider is not used but still interactions should be erased.
+ Otherwise collider should decide on a case-by-case basis, which interaction to erase for good and which to keep in the potential state
+ (without interactionGeometry and interactionPhysics).
+
+ This function doesn't lock pendingEraseMutex, as it is (supposedly) called from no-parallel sections only once per iteration
+ */
+ void unconditionalErasePending();
+ /*! Traverse all pending interactions and erase them if the (T*)->shouldBeErased(id1,id2) return true
+ and keep it if it return false; finally, pendingErase will be clear()'ed.
+
+ Class using this interface (which is presumably a collider) must define the
+
+ bool shouldBeErased(body_id_t, body_id_t) const
+
+ method which will be called for every interaction.
+ */
+ template<class T> void erasePending(const T& t){
+ FOREACH(const vector<body_id_t>& p, pendingErase){ if(t.shouldBeErased(p[0],p[1])) erase(p[0],p[1]); }
+ pendingErase.clear();
+ }
private :
+ #ifdef YADE_OPENMP
+ // This is used only from within requestErase() for now, therefore it can be private
+ boost::mutex pendingEraseMutex;
+ #endif
// used only during serialization/deserialization
vector<shared_ptr<Interaction> > interaction;
protected :
@@ -119,10 +153,6 @@
REGISTER_SERIALIZABLE(InteractionContainer);
-// BOOST_FOREACH compatibility
-#ifndef FOREACH
-# define FOREACH BOOST_FOREACH
-#endif
namespace boost{
template<> struct range_iterator<InteractionContainer>{ typedef InteractionContainer::iterator type; };
Modified: trunk/core/MetaEngine2D.hpp
===================================================================
--- trunk/core/MetaEngine2D.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/core/MetaEngine2D.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -40,11 +40,7 @@
}
/* add functor by its literal name */
virtual void add(string euType){
- shared_ptr<EngineUnitType> eu=dynamic_pointer_cast<EngineUnitType>(ClassFactory::instance().createShared(euType));
- if(!eu){
- cerr<<__FILE__<<":"<<__LINE__<<" WARNING! dynamic cast of engine unit "<<euType<<" failed, will use static_cast. Go figure why."<<endl;
- eu=static_pointer_cast<EngineUnitType>(ClassFactory::instance().createShared(euType));
- }
+ shared_ptr<EngineUnitType> eu=static_pointer_cast<EngineUnitType>(ClassFactory::instance().createShared(euType));
add(eu);
}
Modified: trunk/core/yade.cpp
===================================================================
--- trunk/core/yade.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/core/yade.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -147,10 +147,10 @@
" NDEBUG (heavy optimizations, no assertions and debugging features)\n"
#endif
#ifdef YADE_OPENMP
- " YADE_OPENMP (supports openMP; set OMP_NUM_THREADS env. var to control parallelism.\n"
+ " YADE_OPENMP (supports openMP; set OMP_NUM_THREADS env. var to control parallelism.)\n"
#endif
#ifdef LOG4CXX
- " LOG4CXX configurable logging framework enabled (~/.yade-suffix/logging.conf)\n"
+ " LOG4CXX (configurable logging framework enabled; ~/.yade-" SUFFIX "/logging.conf)\n"
#endif
;
if(!isnan(std::numeric_limits<double>::quiet_NaN())) cerr<<
Modified: trunk/examples/collider-perf/perf.table
===================================================================
--- trunk/examples/collider-perf/perf.table 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/examples/collider-perf/perf.table 2009-05-29 14:35:22 UTC (rev 1786)
@@ -1,76 +1,77 @@
-description nSpheres collider
-128k 128000 'PersistentSAPCollider'
-96k 96000 'PersistentSAPCollider'
-64k 64000 'PersistentSAPCollider'
-56k 56000 'PersistentSAPCollider'
-48k 48000 'PersistentSAPCollider'
-40k 40000 'PersistentSAPCollider'
-36k 36000 'PersistentSAPCollider'
-32k 32000 'PersistentSAPCollider'
-28k 28000 'PersistentSAPCollider'
-24k 24000 'PersistentSAPCollider'
-20k 20000 'PersistentSAPCollider'
-18k 18000 'PersistentSAPCollider'
-16k 16000 'PersistentSAPCollider'
-14k 14000 'PersistentSAPCollider'
-12k 12000 'PersistentSAPCollider'
-10k 10000 'PersistentSAPCollider'
-9k 9000 'PersistentSAPCollider'
-8k 8000 'PersistentSAPCollider'
-7k 7000 'PersistentSAPCollider'
-6k 6000 'PersistentSAPCollider'
-5k 5000 'PersistentSAPCollider'
-4k 4000 'PersistentSAPCollider'
-3k 3000 'PersistentSAPCollider'
-2k 2000 'PersistentSAPCollider'
-1k 1000 'PersistentSAPCollider'
-128k.q 128000 'SpatialQuickSortCollider'
-96k.q 96000 'SpatialQuickSortCollider'
-64k.q 64000 'SpatialQuickSortCollider'
-56k.q 56000 'SpatialQuickSortCollider'
-48k.q 48000 'SpatialQuickSortCollider'
-40k.q 40000 'SpatialQuickSortCollider'
-36k.q 36000 'SpatialQuickSortCollider'
-32k.q 32000 'SpatialQuickSortCollider'
-28k.q 28000 'SpatialQuickSortCollider'
-24k.q 24000 'SpatialQuickSortCollider'
-20k.q 20000 'SpatialQuickSortCollider'
-18k.q 18000 'SpatialQuickSortCollider'
-16k.q 16000 'SpatialQuickSortCollider'
-14k.q 14000 'SpatialQuickSortCollider'
-12k.q 12000 'SpatialQuickSortCollider'
-10k.q 10000 'SpatialQuickSortCollider'
-9k.q 9000 'SpatialQuickSortCollider'
-8k.q 8000 'SpatialQuickSortCollider'
-7k.q 7000 'SpatialQuickSortCollider'
-6k.q 6000 'SpatialQuickSortCollider'
-5k.q 5000 'SpatialQuickSortCollider'
-4k.q 4000 'SpatialQuickSortCollider'
-3k.q 3000 'SpatialQuickSortCollider'
-2k.q 2000 'SpatialQuickSortCollider'
-1k.q 1000 'SpatialQuickSortCollider'
-128k.i 128000 'InsertionSortCollider'
-96k.i 96000 'InsertionSortCollider'
-64k.i 64000 'InsertionSortCollider'
-56k.i 56000 'InsertionSortCollider'
-48k.i 48000 'InsertionSortCollider'
-40k.i 40000 'InsertionSortCollider'
-36k.i 36000 'InsertionSortCollider'
-32k.i 32000 'InsertionSortCollider'
-28k.i 28000 'InsertionSortCollider'
-24k.i 24000 'InsertionSortCollider'
-20k.i 20000 'InsertionSortCollider'
-18k.i 18000 'InsertionSortCollider'
-16k.i 16000 'InsertionSortCollider'
-14k.i 14000 'InsertionSortCollider'
-12k.i 12000 'InsertionSortCollider'
-10k.i 10000 'InsertionSortCollider'
-9k.i 9000 'InsertionSortCollider'
-8k.i 8000 'InsertionSortCollider'
-7k.i 7000 'InsertionSortCollider'
-6k.i 6000 'InsertionSortCollider'
-5k.i 5000 'InsertionSortCollider'
-4k.i 4000 'InsertionSortCollider'
-3k.i 3000 'InsertionSortCollider'
-2k.i 2000 'InsertionSortCollider'
-1k.i 1000 'InsertionSortCollider'
+!OMP_NUM_THREADS description nSpheres collider
+3 128k 128000 'PersistentSAPCollider'
+3 96k 96000 'PersistentSAPCollider'
+3 64k 64000 'PersistentSAPCollider'
+3 56k 56000 'PersistentSAPCollider'
+3 48k 48000 'PersistentSAPCollider'
+3 40k 40000 'PersistentSAPCollider'
+3 36k 36000 'PersistentSAPCollider'
+3 32k 32000 'PersistentSAPCollider'
+3 28k 28000 'PersistentSAPCollider'
+3 24k 24000 'PersistentSAPCollider'
+3 20k 20000 'PersistentSAPCollider'
+3 18k 18000 'PersistentSAPCollider'
+3 16k 16000 'PersistentSAPCollider'
+3 14k 14000 'PersistentSAPCollider'
+3 12k 12000 'PersistentSAPCollider'
+3 10k 10000 'PersistentSAPCollider'
+3 9k 9000 'PersistentSAPCollider'
+3 8k 8000 'PersistentSAPCollider'
+3 7k 7000 'PersistentSAPCollider'
+3 6k 6000 'PersistentSAPCollider'
+3 5k 5000 'PersistentSAPCollider'
+3 4k 4000 'PersistentSAPCollider'
+3 3k 3000 'PersistentSAPCollider'
+3 2k 2000 'PersistentSAPCollider'
+3 1k 1000 'PersistentSAPCollider'
+3 128k.q 128000 'SpatialQuickSortCollider'
+3 96k.q 96000 'SpatialQuickSortCollider'
+3 64k.q 64000 'SpatialQuickSortCollider'
+3 56k.q 56000 'SpatialQuickSortCollider'
+3 48k.q 48000 'SpatialQuickSortCollider'
+3 40k.q 40000 'SpatialQuickSortCollider'
+3 36k.q 36000 'SpatialQuickSortCollider'
+3 32k.q 32000 'SpatialQuickSortCollider'
+3 28k.q 28000 'SpatialQuickSortCollider'
+3 24k.q 24000 'SpatialQuickSortCollider'
+3 20k.q 20000 'SpatialQuickSortCollider'
+3 18k.q 18000 'SpatialQuickSortCollider'
+3 16k.q 16000 'SpatialQuickSortCollider'
+3 14k.q 14000 'SpatialQuickSortCollider'
+3 12k.q 12000 'SpatialQuickSortCollider'
+3 10k.q 10000 'SpatialQuickSortCollider'
+3 9k.q 9000 'SpatialQuickSortCollider'
+3 8k.q 8000 'SpatialQuickSortCollider'
+3 7k.q 7000 'SpatialQuickSortCollider'
+3 6k.q 6000 'SpatialQuickSortCollider'
+3 5k.q 5000 'SpatialQuickSortCollider'
+3 4k.q 4000 'SpatialQuickSortCollider'
+3 3k.q 3000 'SpatialQuickSortCollider'
+3 2k.q 2000 'SpatialQuickSortCollider'
+3 1k.q 1000 'SpatialQuickSortCollider'
+1 128k.i1 128000 'InsertionSortCollider'
+3 128k.i 128000 'InsertionSortCollider'
+3 96k.i 96000 'InsertionSortCollider'
+3 64k.i 64000 'InsertionSortCollider'
+3 56k.i 56000 'InsertionSortCollider'
+3 48k.i 48000 'InsertionSortCollider'
+3 40k.i 40000 'InsertionSortCollider'
+3 36k.i 36000 'InsertionSortCollider'
+3 32k.i 32000 'InsertionSortCollider'
+3 28k.i 28000 'InsertionSortCollider'
+3 24k.i 24000 'InsertionSortCollider'
+3 20k.i 20000 'InsertionSortCollider'
+3 18k.i 18000 'InsertionSortCollider'
+3 16k.i 16000 'InsertionSortCollider'
+3 14k.i 14000 'InsertionSortCollider'
+3 12k.i 12000 'InsertionSortCollider'
+3 10k.i 10000 'InsertionSortCollider'
+3 9k.i 9000 'InsertionSortCollider'
+3 8k.i 8000 'InsertionSortCollider'
+3 7k.i 7000 'InsertionSortCollider'
+3 6k.i 6000 'InsertionSortCollider'
+3 5k.i 5000 'InsertionSortCollider'
+3 4k.i 4000 'InsertionSortCollider'
+3 3k.i 3000 'InsertionSortCollider'
+3 2k.i 2000 'InsertionSortCollider'
+3 1k.i 1000 'InsertionSortCollider'
Added: trunk/examples/concrete/pack/mk-pack.py
===================================================================
--- trunk/examples/concrete/pack/mk-pack.py 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/examples/concrete/pack/mk-pack.py 2009-05-29 14:35:22 UTC (rev 1786)
@@ -0,0 +1,73 @@
+from yade import log
+import shutil, tempfile
+#
+# Generate box or cylindrical packing that is rather compact using TriaxialTest.
+#
+
+# if non-negative, the resulting sample will be cylindrical; if < 0, generate box
+cylAxis=1
+# how many layers of average spheres to cut away from flat parts (relative to mean sphere radius); sphere center count
+cutoffFlat=3.5
+# how many layers of average spheres to cut away from round surfaces (relative to mean sphere radius); sphere centers count
+cutoffRound=2
+
+tt=TriaxialTest(
+ numberOfGrains=500,
+ radiusMean=3e-4,
+ # this is just size ratio if radiusMean is specified
+ # if you comment out the line above, it will be the corner (before compaction) and radiusMean will be set accordingly
+ upperCorner=[2,4,2],
+ radiusStdDev=0, # all spheres exactly the same radius
+
+ ##
+ ## no need to touch any the following, till the end of file
+ ##
+ noFiles=True,
+ lowerCorner=[0,0,0],
+ sigmaIsoCompaction=1e7,
+ sigmaLateralConfinement=1e3,
+ StabilityCriterion=.05,
+ strainRate=.2,
+ fast=True,
+ maxWallVelocity=.1,
+ wallOversizeFactor=2,
+ autoUnload=True, # unload after isotropic compaction
+ autoCompressionActivation=False # stop once unloaded
+)
+
+tt.load()
+log.setLevel('TriaxialCompressionEngine',log.WARN)
+O.run() ## triax stops by itself once unloaded
+O.wait()
+
+# resulting specimen geometry
+ext=utils.aabbExtrema()
+rSphere=tt['radiusMean']
+
+outFile=tempfile.NamedTemporaryFile(delete=False)
+
+if cylAxis<0: # box-shaped packing
+ aabbMin,aabbMax=tuple([ext[0][i]+rSphere*cutoffFlat for i in 0,1,2]),tuple([ext[1][i]-rSphere*cutoffFlat for i in 0,1,2])
+ def isInBox(id):
+ pos=O.bodies[id].phys.pos
+ return utils.ptInAABB(pos,aabbMin,aabbMax)
+ nSpheres=utils.spheresToFile(outFile,consider=isInBox)
+else: # cylinger packing
+ mid =[.5*(ext[1][i]+ext[0][i]) for i in [0,1,2]]
+ extent=[.5*(ext[1][i]-ext[0][i]) for i in [0,1,2]]
+ ax1,ax2=(cylAxis+1)%3,(cylAxis+2)%3
+ cylRadius=min(extent[ax1],extent[ax2])-cutoffRound*rSphere
+ cylHalfHt=extent[cylAxis]-cutoffFlat*rSphere
+ def isInCyl(id):
+ pos=O.bodies[id].phys.pos
+ axisDist=sqrt((pos[ax1]-mid[ax1])**2+(pos[ax2]-mid[ax2])**2)
+ if axisDist>cylRadius: return False
+ axialDist=abs(pos[cylAxis]-mid[cylAxis])
+ if axialDist>cylHalfHt: return False
+ return True
+ nSpheres=utils.spheresToFile(outFile,consider=isInCyl)
+
+outFile2='pack-%d-%s.sphere'%(nSpheres,'box' if cylAxis<0 else 'cyl')
+shutil.move(outFile,outFile2)
+print nSpheres,'spheres written to',outFile2
+quit()
Modified: trunk/extra/Shop.cpp
===================================================================
--- trunk/extra/Shop.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/extra/Shop.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -96,7 +96,7 @@
// get max force on contacts
Real maxContactF=0;
FOREACH(const shared_ptr<Interaction>& I, *rb->transientInteractions){
- if(!I->isReal) continue;
+ if(!I->isReal()) continue;
shared_ptr<NormalShearInteraction> nsi=YADE_PTR_CAST<NormalShearInteraction>(I->interactionPhysics); assert(nsi);
maxContactF=max(maxContactF,(nsi->normalForce+nsi->shearForce).Length());
}
Modified: trunk/extra/tetra/Tetra.cpp
===================================================================
--- trunk/extra/tetra/Tetra.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/extra/tetra/Tetra.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -49,7 +49,7 @@
shared_ptr<TetraBang> bang;
// depending whether it's a new interaction: create new one, or use the existing one.
- if (interaction->isNew) bang=shared_ptr<TetraBang>(new TetraBang());
+ if (!interaction->interactionGeometry) bang=shared_ptr<TetraBang>(new TetraBang());
else bang=YADE_PTR_CAST<TetraBang>(interaction->interactionGeometry);
interaction->interactionGeometry=bang;
@@ -375,7 +375,7 @@
{
for(InteractionContainer::iterator contactI=rootBody->transientInteractions->begin(); contactI!=rootBody->transientInteractions->end(); ++contactI){
- if (!(*contactI)->isReal) continue; // Tetra2TetraBang::go returned false for this interaction, skip it
+ if (!(*contactI)->isReal()) continue; // Tetra2TetraBang::go returned false for this interaction, skip it
const shared_ptr<TetraBang>& contactGeom(dynamic_pointer_cast<TetraBang>((*contactI)->interactionGeometry));
if(!contactGeom) continue;
Modified: trunk/gui/py/_eudoxos.cpp
===================================================================
--- trunk/gui/py/_eudoxos.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/gui/py/_eudoxos.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -39,10 +39,10 @@
MetaBody* rootBody=Omega::instance().getRootBody().get();
shared_ptr<Interaction> I;
FOREACH(I, *rootBody->transientInteractions){
- if(I->isReal) break;
+ if(I->isReal()) break;
}
Real nan=std::numeric_limits<Real>::quiet_NaN();
- if(!I->isReal) {LOG_ERROR("No real interaction found, returning NaN!"); return nan; }
+ if(!I->isReal()) {LOG_ERROR("No real interaction found, returning NaN!"); return nan; }
CpmPhys* BC=dynamic_cast<CpmPhys*>(I->interactionPhysics.get());
if(!BC) {LOG_ERROR("Interaction physics is not CpmPhys instance, returning NaN!"); return nan;}
const Real &omega(BC->omega); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle);
Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/gui/py/_utils.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -24,6 +24,8 @@
Vector3r tuple2vec(const python::tuple& t){return Vector3r(extract<double>(t[0])(),extract<double>(t[1])(),extract<double>(t[2])());}
bool isInBB(Vector3r p, Vector3r bbMin, Vector3r bbMax){return p[0]>bbMin[0] && p[0]<bbMax[0] && p[1]>bbMin[1] && p[1]<bbMax[1] && p[2]>bbMin[2] && p[2]<bbMax[2];}
+bool ptInAABB(python::tuple p, python::tuple bbMin, python::tuple bbMax){return isInBB(tuple2vec(p),tuple2vec(bbMin),tuple2vec(bbMax));}
+
/* \todo implement groupMask */
python::tuple aabbExtrema(Real cutoff=0.0, bool centers=false){
if(cutoff<0. || cutoff>1.) throw invalid_argument("Cutoff must be >=0 and <=1.");
@@ -120,7 +122,7 @@
vector<Real> cummProj(bins,0.);
shared_ptr<MetaBody> rb=Omega::instance().getRootBody();
FOREACH(const shared_ptr<Interaction>& i, *rb->transientInteractions){
- if(!i->isReal) continue;
+ if(!i->isReal()) continue;
const shared_ptr<Body>& b1=Body::byId(i->getId1(),rb), b2=Body::byId(i->getId2(),rb);
if(!b1->maskOk(mask) || !b2->maskOk(mask)) continue;
if(useBB && !isInBB(b1->physicalParameters->se3.position,bbMin,bbMax) && !isInBB(b2->physicalParameters->se3.position,bbMin,bbMax)) continue;
@@ -143,7 +145,7 @@
vector<int> bodyNumInta; bodyNumInta.resize(rb->bodies->size(),-1 /* uninitialized */);
int maxInta=0;
FOREACH(const shared_ptr<Interaction>& i, *rb->transientInteractions){
- if(!i->isReal) continue;
+ if(!i->isReal()) continue;
const body_id_t id1=i->getId1(), id2=i->getId2(); const shared_ptr<Body>& b1=Body::byId(id1,rb), b2=Body::byId(id2,rb);
if(useBB && isInBB(b1->physicalParameters->se3.position,bbMin,bbMax)) bodyNumInta[id1]=bodyNumInta[id1]>0?bodyNumInta[id1]+1:1;
if(useBB && isInBB(b2->physicalParameters->se3.position,bbMin,bbMax)) bodyNumInta[id2]=bodyNumInta[id2]>0?bodyNumInta[id2]+1:1;
@@ -334,6 +336,7 @@
def("PWaveTimeStep",PWaveTimeStep,"Get timestep accoring to the velocity of P-Wave propagation; computed from sphere radii, rigidities and masses.");
def("aabbExtrema",aabbExtrema,aabbExtrema_overloads(args("cutoff","centers"),"Return coordinates of box enclosing all bodies\n centers: do not take sphere radii in account, only their centroids (default=False)\n cutoff: 0-1 number by which the box will be scaled around its center (default=0)"));
+ def("ptInAABB",ptInAABB,"Return True/False whether the point (3-tuple) p is within box given by its min (3-tuple) and max (3-tuple) corners");
def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(args("axis","distFactor"),"Return list of ids for spheres (only) that are on extremal ends of the specimen along given axis; distFactor multiplies their radius so that sphere that do not touch the boundary coordinate can also be returned."));
def("coordsAndDisplacements",coordsAndDisplacements,coordsAndDisplacements_overloads(args("AABB"),"Return tuple of 2 same-length lists for coordinates and displacements (coordinate minus reference coordinate) along given axis (1st arg); if the AABB=((x_min,y_min,z_min),(x_max,y_max,z_max)) box is given, only bodies within this box will be considered."));
def("setRefSe3",setRefSe3,"Set reference positions and orientation of all bodies equal to their current ones.");
Modified: trunk/gui/py/utils.py
===================================================================
--- trunk/gui/py/utils.py 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/gui/py/utils.py 2009-05-29 14:35:22 UTC (rev 1786)
@@ -224,13 +224,18 @@
def spheresToFile(filename,consider=lambda id: True):
"""Save sphere coordinates into ASCII file; the format of the line is: x y z r.
- Non-spherical bodies are silently skipped."""
+ Non-spherical bodies are silently skipped.
+
+ Returns number of spheres that were written."""
o=Omega()
out=open(filename,'w')
+ count=0
for b in o.bodies:
if not b.shape or not b.shape.name=='Sphere' or not consider(b.id): continue
out.write('%g\t%g\t%g\t%g\n'%(b.phys['se3'][0],b.phys['se3'][1],b.phys['se3'][2],b.shape['radius']))
+ count+=1
out.close()
+ return count
def avgNumInteractions(cutoff=0.):
nums,counts=bodyNumInteractionsHistogram(aabbExtrema(cutoff))
Modified: trunk/gui/py/yadeControl.cpp
===================================================================
--- trunk/gui/py/yadeControl.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/gui/py/yadeControl.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -408,12 +408,12 @@
}
/* return nth _real_ iteration from the container (0-based index); this is to facilitate picking random interaction */
pyInteraction pyNth(long n){
- long i=0; FOREACH(const shared_ptr<Interaction>& I, *proxee){ if(!I->isReal) continue; if(i++==n) return pyInteraction(I); }
+ long i=0; FOREACH(const shared_ptr<Interaction>& I, *proxee){ if(!I->isReal()) continue; if(i++==n) return pyInteraction(I); }
throw invalid_argument(string("Interaction number out of range (")+lexical_cast<string>(n)+">="+lexical_cast<string>(i)+").");
}
long len(){return proxee->size();}
void clear(){proxee->clear();}
- python::list withBody(long id){ python::list ret; FOREACH(const shared_ptr<Interaction>& I, *proxee){ if(I->isReal && (I->getId1()==id || I->getId2()==id)) ret.append(pyInteraction(I));} return ret;}
+ python::list withBody(long id){ python::list ret; FOREACH(const shared_ptr<Interaction>& I, *proxee){ if(I->isReal() && (I->getId1()==id || I->getId2()==id)) ret.append(pyInteraction(I));} return ret;}
python::list withBodyAll(long id){ python::list ret; FOREACH(const shared_ptr<Interaction>& I, *proxee){ if(I->getId1()==id || I->getId2()==id) ret.append(pyInteraction(I));} return ret; }
};
Modified: trunk/lib/import/STLReader.hpp
===================================================================
--- trunk/lib/import/STLReader.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/lib/import/STLReader.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -106,19 +106,20 @@
vector<Vrtx> vcs;
set<pair<int,int> > egs;
+ int ret;
/* Read a single facet from an ASCII .STL file */
while(!feof(fp))
{
float n[3];
Vrtx v[3];
- fscanf(fp, "%*s %*s %f %f %f\n", &n[0], &n[1], &n[2]);
- fscanf(fp, "%*s %*s");
- fscanf(fp, "%*s %f %f %f\n", &v[0][0], &v[0][1], &v[0][2]);
- fscanf(fp, "%*s %f %f %f\n", &v[1][0], &v[1][1], &v[1][2]);
- fscanf(fp, "%*s %f %f %f\n", &v[2][0], &v[2][1], &v[2][2]);
- fscanf(fp, "%*s"); // end loop
- fscanf(fp, "%*s"); // end facet
+ ret=fscanf(fp, "%*s %*s %f %f %f\n", &n[0], &n[1], &n[2]);
+ ret=fscanf(fp, "%*s %*s");
+ ret=fscanf(fp, "%*s %f %f %f\n", &v[0][0], &v[0][1], &v[0][2]);
+ ret=fscanf(fp, "%*s %f %f %f\n", &v[1][0], &v[1][1], &v[1][2]);
+ ret=fscanf(fp, "%*s %f %f %f\n", &v[2][0], &v[2][1], &v[2][2]);
+ ret=fscanf(fp, "%*s"); // end loop
+ ret=fscanf(fp, "%*s"); // end facet
if(feof(fp)) break;
int vid[3];
Modified: trunk/pkg/common/Engine/DeusExMachina/CinemCNCEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/CinemCNCEngine.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/DeusExMachina/CinemCNCEngine.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -281,7 +281,7 @@
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
const shared_ptr<Interaction>& contact = *ii;
Modified: trunk/pkg/common/Engine/DeusExMachina/CinemKNCEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/CinemKNCEngine.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/DeusExMachina/CinemKNCEngine.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -274,7 +274,7 @@
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
const shared_ptr<Interaction>& contact = *ii;
Modified: trunk/pkg/common/Engine/DeusExMachina/TranslationEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/TranslationEngine.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/DeusExMachina/TranslationEngine.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -16,7 +16,7 @@
Real dt=Omega::instance().getTimeStep();
const int sign = 1; // ?
FOREACH(body_id_t id,subscribedBodies){
- assert(id<bodies->size());
+ assert(id<(body_id_t)bodies->size());
if(ParticleParameters* p = dynamic_cast<ParticleParameters*>((*bodies)[id]->physicalParameters.get())){
p->se3.position+=sign*dt*velocity*translationAxis;
p->velocity=sign*velocity*translationAxis;
Modified: trunk/pkg/common/Engine/EngineUnit/ElasticBodySimpleRelationship.cpp
===================================================================
--- trunk/pkg/common/Engine/EngineUnit/ElasticBodySimpleRelationship.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/EngineUnit/ElasticBodySimpleRelationship.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -20,7 +20,7 @@
ElasticBodyParameters* s2 = static_cast<ElasticBodyParameters*>(b2.get());
// the need to calculate this is only when the interaction is new
- if( interaction->isNew)
+ if(!interaction->interactionPhysics)
{
boost::shared_ptr<NormalInteraction> sei(new NormalInteraction); // ElasticContactInteraction
// BUG?! kn is stiffness [N], young is modulus [N/m²] !!
Modified: trunk/pkg/common/Engine/MetaEngine/ConstitutiveLawDispatcher.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/ConstitutiveLawDispatcher.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/MetaEngine/ConstitutiveLawDispatcher.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -11,9 +11,8 @@
#else
FOREACH(shared_ptr<Interaction> I, *rootBody->transientInteractions){
#endif
- if(I->isReal){
- assert(I->interactionGeometry);
- assert(I->interactionPhysics);
+ if(I->isReal()){
+ assert(I->interactionGeometry); assert(I->interactionPhysics);
operator()(I->interactionGeometry,I->interactionPhysics,I.get(),rootBody);
}
}
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -1,16 +1,25 @@
#include"InteractionDispatchers.hpp"
YADE_PLUGIN("InteractionDispatchers");
+CREATE_LOGGER(InteractionDispatchers);
InteractionDispatchers::InteractionDispatchers(){
geomDispatcher=shared_ptr<InteractionGeometryMetaEngine>(new InteractionGeometryMetaEngine);
physDispatcher=shared_ptr<InteractionPhysicsMetaEngine>(new InteractionPhysicsMetaEngine);
constLawDispatcher=shared_ptr<ConstitutiveLawDispatcher>(new ConstitutiveLawDispatcher);
+ alreadyWarnedNoCollider=false;
}
#define DISPATCH_CACHE
void InteractionDispatchers::action(MetaBody* rootBody){
+ if(rootBody->interactions->pendingErase.size()>0){
+ if(!alreadyWarnedNoCollider){
+ LOG_WARN("Interactions pending erase found, no collider being used?");
+ alreadyWarnedNoCollider=true;
+ }
+ rootBody->interactions->unconditionalErasePending();
+ }
#ifdef YADE_OPENMP
const long size=rootBody->interactions->size();
#pragma omp parallel for
@@ -24,9 +33,9 @@
const shared_ptr<Body>& b2_=Body::byId(I->getId2(),rootBody);
// we know there is no geometry functor already, take the short path
- if(!I->functorCache.geomExists) { I->isReal=false; continue; }
+ if(!I->functorCache.geomExists) { assert(!I->isReal()); continue; }
// no interaction geometry for either of bodies; no interaction possible
- if(!b1_->interactingGeometry || !b2_->interactingGeometry) { I->isReal=false; continue; }
+ if(!b1_->interactingGeometry || !b2_->interactingGeometry) { assert(!I->isReal()); continue; }
bool swap=false;
// InteractionGeometryMetaEngine
@@ -44,8 +53,8 @@
const shared_ptr<Body>& b2=Body::byId(I->getId2(),rootBody);
assert(I->functorCache.geom);
- I->isReal=I->functorCache.geom->go(b1->interactingGeometry,b2->interactingGeometry,b1->physicalParameters->se3, b2->physicalParameters->se3,I);
- if(!I->isReal) continue;
+ bool geomCreated=I->functorCache.geom->go(b1->interactingGeometry,b2->interactingGeometry,b1->physicalParameters->se3, b2->physicalParameters->se3,I);
+ if(!geomCreated) continue;
// InteractionPhysicsMetaEngine
if(!I->functorCache.phys){
@@ -54,7 +63,10 @@
}
assert(I->functorCache.phys);
I->functorCache.phys->go(b1->physicalParameters,b2->physicalParameters,I);
+ assert(I->interactionPhysics);
+ I->iterMadeReal=rootBody->currentIteration; // mark the interaction as created right now
+
// ConstitutiveLawDispatcher
// populating constLaw cache must be done after geom and physics dispatchers have been called, since otherwise the interaction
// would not have interactionGeometry and interactionPhysics yet.
@@ -64,15 +76,14 @@
}
assert(I->functorCache.constLaw);
I->functorCache.constLaw->go(I->interactionGeometry,I->interactionPhysics,I.get(),rootBody);
-
#else
const shared_ptr<Body>& b1=Body::byId(I->getId1(),rootBody);
const shared_ptr<Body>& b2=Body::byId(I->getId2(),rootBody);
// InteractionGeometryMetaEngine
- I->isReal =
+ bool geomCreated =
b1->interactingGeometry && b2->interactingGeometry && // some bodies do not have interactingGeometry
geomDispatcher->operator()(b1->interactingGeometry, b2->interactingGeometry, b1->physicalParameters->se3, b2->physicalParameters->se3,I);
- if(!I->isReal) continue;
+ if(!geomCreated) continue;
// InteractionPhysicsMetaEngine
// geom may have swapped bodies, get bodies again
physDispatcher->operator()(Body::byId(I->getId1(),rootBody)->physicalParameters, Body::byId(I->getId2(),rootBody)->physicalParameters,I);
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.hpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -6,6 +6,7 @@
#include<yade/pkg-common/ConstitutiveLawDispatcher.hpp>
class InteractionDispatchers: public StandAloneEngine {
+ bool alreadyWarnedNoCollider;
public:
InteractionDispatchers();
virtual void action(MetaBody*);
@@ -18,5 +19,6 @@
(physDispatcher)
(constLawDispatcher)
);
+ DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(InteractionDispatchers);
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -13,6 +13,8 @@
#include<yade/core/MetaBody.hpp>
+CREATE_LOGGER(InteractionGeometryMetaEngine);
+
/* Single dispatch for given pair of bodies, returning new interaction.
*
* The caller is responsible for inserting the interaction into some interaction container.
@@ -27,21 +29,30 @@
//bool op=operator()(b1->interactingGeometry,b2->interactingGeometry,b1->physicalParameters->se3,b2->physicalParameters->se3,i);
//if(!op) throw runtime_error("InteractionGeometryMetaEngine::explicitAction could not dispatch for given types ("+b1->interactingGeometry->getClassName()+","+b2->interactingGeometry->getClassName()+") or the dispatchee returned false.");
//return i;
+
+ // FIXME: not clear why it is not a good idea. If I really wnat interaction that I call this function, I also want it to say loudly that it failed.
// Seems asserts and throws in code above is not good idea.
// Below code do same (i.e. create interaction for specified bodies), but
- // without artifical exceptions. If creating interaction is fail (for
+ // without artifical exceptions. If creating interaction fails (for
// example if bodies don't have an interactionGeometry), returned
// interaction is non real, i.e. interaction->isReal==false. Sega.
shared_ptr<Interaction> interaction(new Interaction(b1->getId(),b2->getId()));
- interaction->isReal =
- b1->interactingGeometry && b2->interactingGeometry &&
- operator()( b1->interactingGeometry , b2->interactingGeometry , b1->physicalParameters->se3 , b2->physicalParameters->se3 , interaction );
+ b1->interactingGeometry && b2->interactingGeometry && operator()( b1->interactingGeometry , b2->interactingGeometry , b1->physicalParameters->se3 , b2->physicalParameters->se3 , interaction );
return interaction;
}
void InteractionGeometryMetaEngine::action(MetaBody* ncb)
{
+ // Erase interaction that were requested for erase, but not processed by the collider, if any (and warn once about that, as it is suspicious)
+ if(ncb->interactions->pendingErase.size()>0){
+ if(!alreadyWarnedNoCollider){
+ LOG_WARN("Interactions pending erase found, no collider being used?");
+ alreadyWarnedNoCollider=true;
+ }
+ ncb->interactions->unconditionalErasePending();
+ }
+
shared_ptr<BodyContainer>& bodies = ncb->bodies;
#ifdef YADE_OPENMP
const long size=ncb->transientInteractions->size();
@@ -53,9 +64,8 @@
#endif
const shared_ptr<Body>& b1=(*bodies)[interaction->getId1()];
const shared_ptr<Body>& b2=(*bodies)[interaction->getId2()];
- interaction->isReal =
- b1->interactingGeometry && b2->interactingGeometry && // some bodies do not have interactingGeometry
- operator()(b1->interactingGeometry, b2->interactingGeometry, b1->physicalParameters->se3, b2->physicalParameters->se3, interaction);
+ b1->interactingGeometry && b2->interactingGeometry && // some bodies do not have interactingGeometry
+ operator()(b1->interactingGeometry, b2->interactingGeometry, b1->physicalParameters->se3, b2->physicalParameters->se3, interaction);
}
}
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.hpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -34,13 +34,15 @@
, false // disable auto symmetry handling
>
{
+ bool alreadyWarnedNoCollider;
public :
virtual void action(MetaBody*);
shared_ptr<Interaction> explicitAction(const shared_ptr<Body>& b1, const shared_ptr<Body>& b2);
+ InteractionGeometryMetaEngine(): alreadyWarnedNoCollider(false){}
- REGISTER_CLASS_NAME(InteractionGeometryMetaEngine);
- REGISTER_BASE_CLASS_NAME(MetaEngine2D);
+ REGISTER_CLASS_AND_BASE(InteractionGeometryMetaEngine,MetaEngine2D);
REGISTER_ATTRIBUTES(MetaEngine,/* no attributes here*/ );
+ DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(InteractionGeometryMetaEngine);
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionPhysicsMetaEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionPhysicsMetaEngine.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionPhysicsMetaEngine.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -16,7 +16,7 @@
*/
void InteractionPhysicsMetaEngine::explicitAction(shared_ptr<PhysicalParameters>& pp1, shared_ptr<PhysicalParameters>& pp2, shared_ptr<Interaction>& i){
// should we throw instead of asserting?
- assert(i->isReal);
+ assert(i->isReal());
operator()(pp1,pp2,i);
}
@@ -32,10 +32,13 @@
#else
FOREACH(const shared_ptr<Interaction>& interaction, *ncb->interactions){
#endif
- shared_ptr<Body>& b1 = (*bodies)[interaction->getId1()];
- shared_ptr<Body>& b2 = (*bodies)[interaction->getId2()];
- if (interaction->isReal)
+ if(interaction->interactionGeometry){
+ shared_ptr<Body>& b1 = (*bodies)[interaction->getId1()];
+ shared_ptr<Body>& b2 = (*bodies)[interaction->getId2()];
operator()(b1->physicalParameters, b2->physicalParameters, interaction);
+ assert(interaction->interactionPhysics);
+ interaction->iterMadeReal=ncb->currentIteration;
+ }
}
}
Modified: trunk/pkg/common/Engine/StandAloneEngine/DistantPersistentSAPCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/DistantPersistentSAPCollider.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/StandAloneEngine/DistantPersistentSAPCollider.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -85,9 +85,8 @@
for( ; ii!=iiEnd ; ++ii)
{
shared_ptr<Interaction> interaction = *ii;
- // FIXME : remove this isNew flag and test if interactionPhysic ?
- if (interaction->isReal) // if a interaction was only potential then no geometry was created for it and so this time it is still a new one
- interaction->isNew = false;
+ //if (interaction->isReal()) // if a interaction was only potential then no geometry was created for it and so this time it is still a new one
+ // interaction->isNew = false;
//interaction->isReal = false;
}
@@ -200,7 +199,7 @@
if (overlap && !found)
transientInteractions->insert(body_id_t(id1),body_id_t(id2));
// removes the pair p=(id1,id2) if the two AABB do not overlapp any more and if p already exists in the overlappingBB
- else if (!overlap && found && !interaction->isReal)
+ else if (!overlap && found && !interaction->isReal())
transientInteractions->erase(body_id_t(id1),body_id_t(id2));//Bruno's hack
//else if (!overlap && found)
// transientInteractions->erase(body_id_t(id1),body_id_t(id2));
Modified: trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -16,11 +16,11 @@
CREATE_LOGGER(InsertionSortCollider);
// return true if bodies bb overlap in all 3 dimensions
-bool InsertionSortCollider::spatialOverlap(body_id_t id1, body_id_t id2){
+bool InsertionSortCollider::spatialOverlap(body_id_t id1, body_id_t id2) const {
return
- (minima[3*id1+0]<maxima[3*id2+0]) && (maxima[3*id1+0]>minima[3*id2+0]) &&
- (minima[3*id1+1]<maxima[3*id2+1]) && (maxima[3*id1+1]>minima[3*id2+1]) &&
- (minima[3*id1+2]<maxima[3*id2+2]) && (maxima[3*id1+2]>minima[3*id2+2]);
+ (minima[3*id1+0]<=maxima[3*id2+0]) && (maxima[3*id1+0]>=minima[3*id2+0]) &&
+ (minima[3*id1+1]<=maxima[3*id2+1]) && (maxima[3*id1+1]>=minima[3*id2+1]) &&
+ (minima[3*id1+2]<=maxima[3*id2+2]) && (maxima[3*id1+2]>=minima[3*id2+2]);
}
// called by the insertion sort if 2 bodies swapped their bounds
@@ -32,16 +32,16 @@
bool hasInter=(bool)I;
// interaction doesn't exist and shouldn't, or it exists and should
if(!overlap && !hasInter) return;
- if(overlap && hasInter){ /* FIXME: should check I->isNew and I->isReal; etc */ return; }
+ if(overlap && hasInter){ return; }
// create interaction if not yet existing
if(overlap && !hasInter){ // second condition only for readability
if(!Collider::mayCollide(Body::byId(id1,rb).get(),Body::byId(id2,rb).get())) return;
- LOG_TRACE("Creating new interaction #"<<id1<<"+#"<<id2);
+ // LOG_TRACE("Creating new interaction #"<<id1<<"+#"<<id2);
shared_ptr<Interaction> newI=shared_ptr<Interaction>(new Interaction(id1,id2));
interactions->insert(newI);
return;
}
- if(!overlap && hasInter){ if(!I->isReal && I->isNew) interactions->erase(id1,id2); return; }
+ if(!overlap && hasInter){ if(!I->isReal()) interactions->erase(id1,id2); return; }
assert(false); // unreachable
}
@@ -110,6 +110,8 @@
//timingDeltas->checkpoint("copy");
// process interactions that the constitutive law asked to be erased
+ interactions->erasePending(*this);
+ #if 0
FOREACH(const InteractionContainer::bodyIdPair& p, interactions->pendingErase){
// remove those that do not overlap spatially anymore
if(!spatialOverlap(p[0],p[1])){ interactions->erase(p[0],p[1]); LOG_TRACE("Deleted interaction #"<<p[0]<<"+#"<<p[1]); }
@@ -121,6 +123,7 @@
}
}
interactions->pendingErase.clear();
+ #endif
// sort
@@ -130,7 +133,16 @@
}
else {
if(doInitSort){
- std::sort(XX.begin(),XX.end()); std::sort(YY.begin(),YY.end()); std::sort(ZZ.begin(),ZZ.end());
+ // the initial sort is in independent in 3 dimensions, may be run in parallel
+ #pragma omp parallel sections
+ {
+ #pragma omp section
+ std::sort(XX.begin(),XX.end());
+ #pragma omp section
+ std::sort(YY.begin(),YY.end());
+ #pragma omp section
+ std::sort(ZZ.begin(),ZZ.end());
+ }
} else {
insertionSort(XX,interactions,rb,false); insertionSort(YY,interactions,rb,false); insertionSort(ZZ,interactions,rb,false);
}
@@ -156,18 +168,4 @@
}
}
//timingDeltas->checkpoint("sort&collide");
-
-#if 0
- // garbage collection once in a while: for interactions that were still real when the bounding boxes separated
- // the collider would never get to see them again otherwise
- if(iter%1000==0){
- typedef pair<body_id_t,body_id_t> bodyIdPair;
- list<bodyIdPair> toBeDeleted;
- FOREACH(const shared_ptr<Interaction>& I,*interactions){
- if(!I->isReal && (!I->isNew || !spatialOverlap(I->getId1(),I->getId2()))) toBeDeleted.push_back(bodyIdPair(I->getId1(),I->getId2()));
- }
- FOREACH(const bodyIdPair& p, toBeDeleted){ interactions->erase(p.first,p.second); LOG_TRACE("Deleted interaction #"<<p.first<<"+#"<<p.second); }
- }
- //timingDeltas->checkpoint("stale");
-#endif
}
Modified: trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -9,7 +9,7 @@
Insertion sort is used for sorting the bound list that is already pre-sorted from last iteration, where each inversion
calls checkOverlap which then handles either overlap (by creating interaction if necessary) or its absence (by deleting
- interaction if it exists and is only potential (!isReal && isNew).
+ interaction if it is only potential).
Bodies without bounding volume are ahndle gracefully and never collide.
@@ -41,13 +41,16 @@
*/
void insertionSort(std::vector<Bound>& v,InteractionContainer*,MetaBody*,bool doCollide=true);
void handleBoundInversion(body_id_t,body_id_t,InteractionContainer*,MetaBody*);
- bool spatialOverlap(body_id_t,body_id_t);
+ bool spatialOverlap(body_id_t,body_id_t) const;
public:
//! axis for the initial sort
int sortAxis;
//! if true, separate sorting and colliding phase; MUCH slower, but processes all interactions at every step
+ // This makes the collider non-persistent, not remembering last state
bool sortThenCollide;
+ //! Predicate called from loop within InteractionContainer::erasePending
+ bool shouldBeErased(body_id_t id1, body_id_t id2) const { return !spatialOverlap(id1,id2); }
InsertionSortCollider(): sortAxis(0), sortThenCollide(false){ /* timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);*/ }
virtual void action(MetaBody*);
Modified: trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -36,6 +36,8 @@
}
+// template void InteractionContainer::erasePending<PersistentSAPCollider>(const PersistentSAPCollider&);
+
void PersistentSAPCollider::action(MetaBody* ncb)
{
rootBody=ncb;
@@ -80,25 +82,7 @@
// timingDeltas->checkpoint("minMaxUpdate");
- typedef pair<body_id_t,body_id_t> bodyIdPair;
- list<bodyIdPair> toBeDeleted;
- FOREACH(const shared_ptr<Interaction>& I,*ncb->transientInteractions){
- // TODO: this logic will be unitedly in Collider::handleExistingInteraction
- //
- // remove interactions deleted by the constitutive law: thay are not new, but nor real either
- // to make sure, do that only with haveDistantTransient
- if(haveDistantTransient && !I->isNew && !I->isReal) { toBeDeleted.push_back(bodyIdPair(I->getId1(),I->getId2())); continue; }
- // Once the interaction has been fully created, it is not "new" anymore
- if (I->isReal) I->isNew=false;
- // OTOH if is is now real anymore, it falls back to the potential state
- if(!haveDistantTransient && !I->isReal) I->isNew=true;
- // for non-distant interactions, isReal depends on whether there is geometrical overlap; that is calculated later
- // for distant: if isReal&&!isNew means:
- // the interaction was marked (by the constitutive law) as not real anymore should be deleted
- if(!haveDistantTransient) I->isReal=false;
- //if(!I->isReal){LOG_DEBUG("Interaction #"<<I->getId1()<<"=#"<<I->getId2()<<" is not real.");}
- }
- FOREACH(const bodyIdPair& p, toBeDeleted){ transientInteractions->erase(p.first,p.second); }
+ ncb->interactions->erasePending(*this);
// timingDeltas->checkpoint("deleteInvalid");
@@ -236,6 +220,16 @@
}
}
+
+bool PersistentSAPCollider::shouldBeErased(body_id_t id1, body_id_t id2) const {
+ // if there is no bbox overlap
+ int offset1=3*id1, offset2=3*id2;
+ return (
+ maxima[offset1 ]<minima[offset2 ] || maxima[offset2 ]<minima[offset1 ] ||
+ maxima[offset1+1]<minima[offset2+1] || maxima[offset2+1]<minima[offset1+1] ||
+ maxima[offset1+2]<minima[offset2+2] || maxima[offset2+2]<minima[offset1+2] );
+}
+
/* Note that this function is called only for bodies that actually overlap along some axis */
void PersistentSAPCollider::updateOverlapingBBSet(int id1,int id2){
// look if the pair (id1,id2) already exists in the overlappingBB collection
@@ -260,7 +254,7 @@
transientInteractions->insert(body_id_t(id1),body_id_t(id2));
}
// removes the pair p=(id1,id2) if the two AABB do not overlapp any more and if p already exists in the overlappingBB
- else if(!overlap && found && (haveDistantTransient ? !interaction->isReal : true) ){
+ else if(!overlap && found && (haveDistantTransient ? !interaction->isReal() : true) ){
//LOG_DEBUG("Erasing interaction #"<<id1<<"=#"<<id2<<" (isReal="<<interaction->isReal<<")");
transientInteractions->erase(body_id_t(id1),body_id_t(id2));
}
Modified: trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.hpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -99,8 +99,7 @@
//! Don't break transient interaction once bodies don't overlap anymore; material law will be responsible for breaking it.
bool haveDistantTransient;
- //! minimum number of bodies to run updateIds in parallel secions; if 0 (default for now), never run in parallel
- //long ompBodiesMin;
+ bool shouldBeErased(body_id_t, body_id_t) const;
void registerAttributes(){
Collider::registerAttributes();
Modified: trunk/pkg/common/Engine/StandAloneEngine/PersistentTriangulationCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/PersistentTriangulationCollider.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/StandAloneEngine/PersistentTriangulationCollider.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -97,9 +97,9 @@
InteractionContainer::iterator I_end = transientInteractions->end();
for ( InteractionContainer::iterator I=transientInteractions->begin(); I!=I_end; ++I )
{
-
- if ( ( *I )->isReal ) ( *I )->isNew=false;
- if ( !haveDistantTransient ) ( *I )->isReal=false; // reset this flag, is used later... (??)
+ // FIXME: eudoxos commented out as isReal and isNew is removed...
+ // if ( ( *I )->isReal ) ( *I )->isNew=false;
+ // if ( !haveDistantTransient ) ( *I )->isReal=false; // reset this flag, is used later... (??)
( *I )->isNeighbor = false;// will be set again just below
}
@@ -132,7 +132,7 @@
I_end = transientInteractions->end();
for ( InteractionContainer::iterator I=transientInteractions->begin(); I!=I_end; ++I )
{
- if ( ( ! ( *I )->isNeighbor ) && ( haveDistantTransient ? ! ( *I )->isReal : true ) )
+ if ( ( ! ( *I )->isNeighbor ) && ( haveDistantTransient ? ! ( *I )->isReal() : true ) )
{
toErase.push_back ( pair<unsigned int,unsigned int> ( ( *I )->getId1() , ( *I )->getId2() ) );
//cerr << "to delete " << ( *I )->getId1() << "-" << ( *I )->getId2() << "(isNeighbor=" << ( *I )->isNeighbor<< endl;
Modified: trunk/pkg/common/Engine/StandAloneEngine/SpatialQuickSortCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/SpatialQuickSortCollider.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/StandAloneEngine/SpatialQuickSortCollider.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -92,7 +92,8 @@
}
interaction->cycle=true;
// if interaction !isReal it falls back to the potential state
- if (!haveDistantTransient && !interaction->isReal) interaction->isNew=true;
+ // FIXME: eudoxos: do nothing here... is that correct?
+ // if (!haveDistantTransient && !interaction->isReal()) interaction->isNew=true;
}
}
}
@@ -104,18 +105,13 @@
interaction = *ii;
if(
// if haveDistantTransient remove interactions deleted by the constitutive law
- (haveDistantTransient && !interaction->isNew && !interaction->isReal)
+ (haveDistantTransient && !interaction->isReal())
// if !haveDistantTransient remove interactions without AABB overlapping
|| (!haveDistantTransient && !interaction->cycle)
) {
transientInteractions->erase( interaction->getId1(), interaction->getId2() );
continue;
}
- // Once the interaction has been fully created, it is not "new" anymore
- if(interaction->isReal) interaction->isNew=false;
- // For non-distant interactions reset isReal (it would be calculated later);
- // for distant interactions isReal is set/unset by constitutive law
- if(!haveDistantTransient) interaction->isReal=false;
}
}
Modified: trunk/pkg/common/Engine/StandAloneEngine/SpheresFactory.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/SpheresFactory.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/common/Engine/StandAloneEngine/SpheresFactory.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -92,7 +92,7 @@
bool is_overlap=false;
for( unsigned int i=0, e=bI->probedBodies.size(); i<e; ++i)
{
- if (iGME->explicitAction(sphere,Body::byId(bI->probedBodies[i]))->isReal)
+ if (iGME->explicitAction(sphere,Body::byId(bI->probedBodies[i]))->interactionGeometry)
{
is_overlap=true;
break;
Modified: trunk/pkg/dem/ConcretePM.cpp
===================================================================
--- trunk/pkg/dem/ConcretePM.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/ConcretePM.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -16,51 +16,46 @@
CREATE_LOGGER(Ip2_CpmMat_CpmMat_CpmPhys);
void Ip2_CpmMat_CpmMat_CpmPhys::go(const shared_ptr<PhysicalParameters>& pp1, const shared_ptr<PhysicalParameters>& pp2, const shared_ptr<Interaction>& interaction){
+ if(interaction->interactionPhysics) return;
+
Dem3DofGeom* contGeom=YADE_CAST<Dem3DofGeom*>(interaction->interactionGeometry.get());
+ assert(contGeom);
- assert(contGeom); // for now, don't handle anything other than SpheresContactGeometry and Dem3DofGeom
+ const shared_ptr<BodyMacroParameters>& elast1=static_pointer_cast<BodyMacroParameters>(pp1);
+ const shared_ptr<BodyMacroParameters>& elast2=static_pointer_cast<BodyMacroParameters>(pp2);
- if(!interaction->isNew && interaction->interactionPhysics){ /* relax */ }
- else {
- interaction->isNew=false; // just in case
+ Real E12=2*elast1->young*elast2->young/(elast1->young+elast2->young); // harmonic Young's modulus average
+ //Real nu12=2*elast1->poisson*elast2->poisson/(elast1->poisson+elast2->poisson); // dtto for Poisson ratio
+ Real minRad=(contGeom->refR1<=0?contGeom->refR2:(contGeom->refR2<=0?contGeom->refR1:min(contGeom->refR1,contGeom->refR2)));
+ Real S12=Mathr::PI*pow(minRad,2); // "surface" of interaction
+ //Real E=(E12 /* was here for Kn: *S12/d0 */)*((1+alpha)/(beta*(1+nu12)+gamma*(1-alpha*nu12)));
+ //Real E=E12; // apply alpha, beta, gamma: garbage values of E !?
- const shared_ptr<BodyMacroParameters>& elast1=static_pointer_cast<BodyMacroParameters>(pp1);
- const shared_ptr<BodyMacroParameters>& elast2=static_pointer_cast<BodyMacroParameters>(pp2);
+ if(!neverDamage) { assert(!isnan(sigmaT)); }
- Real E12=2*elast1->young*elast2->young/(elast1->young+elast2->young); // harmonic Young's modulus average
- //Real nu12=2*elast1->poisson*elast2->poisson/(elast1->poisson+elast2->poisson); // dtto for Poisson ratio
- Real minRad=(contGeom->refR1<=0?contGeom->refR2:(contGeom->refR2<=0?contGeom->refR1:min(contGeom->refR1,contGeom->refR2)));
- Real S12=Mathr::PI*pow(minRad,2); // "surface" of interaction
- //Real E=(E12 /* was here for Kn: *S12/d0 */)*((1+alpha)/(beta*(1+nu12)+gamma*(1-alpha*nu12)));
- //Real E=E12; // apply alpha, beta, gamma: garbage values of E !?
+ shared_ptr<CpmPhys> contPhys(new CpmPhys());
- if(!neverDamage) { assert(!isnan(sigmaT)); }
+ contPhys->E=E12;
+ contPhys->G=E12*G_over_E;
+ contPhys->tanFrictionAngle=tan(.5*(elast1->frictionAngle+elast2->frictionAngle));
+ contPhys->undamagedCohesion=sigmaT;
+ contPhys->crossSection=S12;
+ contPhys->epsCrackOnset=epsCrackOnset;
+ contPhys->epsFracture=relDuctility*epsCrackOnset;
+ // inherited from NormalShearInteracion, used in the timestepper
+ contPhys->kn=contPhys->E*contPhys->crossSection;
+ contPhys->ks=contPhys->G*contPhys->crossSection;
- shared_ptr<CpmPhys> contPhys(new CpmPhys());
+ if(neverDamage) contPhys->neverDamage=true;
+ if(cohesiveThresholdIter<0 || Omega::instance().getCurrentIteration()<cohesiveThresholdIter) contPhys->isCohesive=true;
+ else contPhys->isCohesive=false;
+ contPhys->dmgTau=dmgTau;
+ contPhys->dmgRateExp=dmgRateExp;
+ contPhys->plTau=plTau;
+ contPhys->plRateExp=plRateExp;
+ contPhys->isoPrestress=isoPrestress;
- contPhys->E=E12;
- contPhys->G=E12*G_over_E;
- contPhys->tanFrictionAngle=tan(.5*(elast1->frictionAngle+elast2->frictionAngle));
- contPhys->undamagedCohesion=sigmaT;
- contPhys->crossSection=S12;
- contPhys->epsCrackOnset=epsCrackOnset;
- contPhys->epsFracture=relDuctility*epsCrackOnset;
- contPhys->omegaThreshold=omegaThreshold;
- // inherited from NormalShearInteracion, used in the timestepper
- contPhys->kn=contPhys->E*contPhys->crossSection;
- contPhys->ks=contPhys->G*contPhys->crossSection;
-
- if(neverDamage) contPhys->neverDamage=true;
- if(cohesiveThresholdIter<0 || Omega::instance().getCurrentIteration()<cohesiveThresholdIter) contPhys->isCohesive=true;
- else contPhys->isCohesive=false;
- contPhys->dmgTau=dmgTau;
- contPhys->dmgRateExp=dmgRateExp;
- contPhys->plTau=plTau;
- contPhys->plRateExp=plRateExp;
- contPhys->isoPrestress=isoPrestress;
-
- interaction->interactionPhysics=contPhys;
- }
+ interaction->interactionPhysics=contPhys;
}
@@ -127,21 +122,15 @@
Real Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2=1.; /* deactivated if > 0 */
Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed=1.;
Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift=0.;
+Real Law2_Dem3DofGeom_CpmPhys_Cpm::omegaThreshold=0.;
void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
//timingDeltas->start();
Dem3DofGeom* contGeom=static_cast<Dem3DofGeom*>(_geom.get());
CpmPhys* BC=static_cast<CpmPhys*>(_phys.get());
- /* kept fully damaged contacts; note that normally the contact is deleted _after_ the CPM_MATERIAL_MODEL,
- * i.e. if it is 1.0 here, omegaThreshold is >= 1.0 for sure.
- * &&'ing that just to make sure anyway ...
- */
- // if(BC->omega>=1.0 && BC->omegaThreshold>=1.0) return;
-
// shorthands
- Real& epsN(BC->epsN); Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& omegaThreshold(BC->omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& dt=Omega::instance().getTimeStep(); const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage); const Real& dmgTau(BC->dmgTau); const Real& plTau(BC->plTau); const bool& isCohesive(BC->isCohesive);
- /* const Real& transStrainCoeff(BC->transStrainCoeff); const Real& epsTrans(BC->epsTrans); const Real& xiShear(BC->xiShear); */
+ Real& epsN(BC->epsN); Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& omegaThreshold(Law2_Dem3DofGeom_CpmPhys_Cpm::omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& dt=Omega::instance().getTimeStep(); const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage); const Real& dmgTau(BC->dmgTau); const Real& plTau(BC->plTau); const bool& isCohesive(BC->isCohesive);
Real& omega(BC->omega); Real& sigmaN(BC->sigmaN); Vector3r& sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
const Real& yieldLogSpeed(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed); const int& yieldSurfType(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldSurfType);
const Real& yieldEllipseShift(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift);
@@ -171,12 +160,20 @@
//timingDeltas->checkpoint("geom");
epsN+=BC->isoPrestress/E;
- //TRVAR1(epsN);
#ifdef CPM_MATERIAL_MODEL
CPM_MATERIAL_MODEL
#else
- sigmaN=E*epsN;
- sigmaT=G*epsT;
+ // very simplified version of the constitutive law
+ kappaD=max(max(0,epsN),kappaD); // internal variable, max positive strain (non-decreasing)
+ omega=isCohesive?funcG(kappaD,epsCrackOnset,epsFracture,neverDamage):1.; // damage variable (non-decreasing, as funcG is also non-decreasing)
+ sigmaN=(1-(epsN>0?omega:0))*E*epsN; // damage taken in account in tension only
+ sigmaT=G*epsT; // trial stress
+ Real yieldSigmaT=max((Real)0.,undamagedCohesion*(1-omega)-sigmaN*tanFrictionAngle); // Mohr-Coulomb law with damage
+ if(sigmaT.SquaredLength()>yieldSigmaT*yieldSigmaT){
+ sigmaT*=yieldSigmaT/sigmaT.Length(); // stress return
+ epsPlSum+=rT*contGeom->slipToStrainTMax(rT/G); // adjust strain
+ }
+ relResidualStrength=isCohesive?(kappaD<epsCrackOnset?1.:(1-omega)*(kappaD)/epsCrackOnset):0;
#endif
sigmaN-=BC->isoPrestress;
if(contGeom->refR1<0 && Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2<=0 && epsN<Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2){
@@ -303,7 +300,7 @@
// get max force on contacts
Real maxContactF=0;
FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
- if(!I->isReal) continue;
+ if(!I->isReal()) continue;
shared_ptr<CpmPhys> BC=YADE_PTR_CAST<CpmPhys>(I->interactionPhysics); assert(BC);
maxContactF=max(maxContactF,max(BC->Fn,BC->Fs.Length()));
CpmMat* bpp1(YADE_CAST<CpmMat*>(Body::byId(I->getId1())->physicalParameters.get()));
@@ -314,7 +311,7 @@
unbalancedForce=(useMaxForce?maxF:meanF)/maxContactF;
FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
- if(!I->isReal) continue;
+ if(!I->isReal()) continue;
shared_ptr<CpmPhys> BC=YADE_PTR_CAST<CpmPhys>(I->interactionPhysics); assert(BC);
CpmMat* bpp1(YADE_CAST<CpmMat*>(Body::byId(I->getId1())->physicalParameters.get()));
CpmMat* bpp2(YADE_CAST<CpmMat*>(Body::byId(I->getId2())->physicalParameters.get()));
Modified: trunk/pkg/dem/ConcretePM.hpp
===================================================================
--- trunk/pkg/dem/ConcretePM.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/ConcretePM.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -104,10 +104,6 @@
/// since the softening law is exponential, this doesn't mean that the contact is fully damaged at this point,
/// that happens only asymptotically
epsFracture,
- //! damage after which the contact disappears (<1), since omega reaches 1 only for strain →+∞
- omegaThreshold,
- //! weight coefficient for shear strain when computing the strain semi-norm kappaD
- xiShear,
//! characteristic time for damage (if non-positive, the law without rate-dependence is used)
dmgTau,
//! exponent in the rate-dependent damage evolution
@@ -146,7 +142,7 @@
- CpmPhys(): NormalShearInteraction(),E(0), G(0), tanFrictionAngle(0), undamagedCohesion(0), crossSection(0), xiShear(0), dmgTau(-1), dmgRateExp(0), dmgStrain(0), plTau(-1), plRateExp(0), isoPrestress(0.), kappaD(0.), epsTrans(0.), epsPlSum(0.) { createIndex(); epsT=Vector3r::ZERO; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; epsPlSum=0; dmgOverstress=0; }
+ CpmPhys(): NormalShearInteraction(),E(0), G(0), tanFrictionAngle(0), undamagedCohesion(0), crossSection(0), dmgTau(-1), dmgRateExp(0), dmgStrain(0), plTau(-1), plRateExp(0), isoPrestress(0.), kappaD(0.), epsTrans(0.), epsPlSum(0.) { createIndex(); epsT=Vector3r::ZERO; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; epsPlSum=0; dmgOverstress=0; }
virtual ~CpmPhys();
REGISTER_ATTRIBUTES(NormalShearInteraction,
@@ -157,8 +153,6 @@
(crossSection)
(epsCrackOnset)
(epsFracture)
- (omegaThreshold)
- (xiShear)
(dmgTau)
(dmgRateExp)
(dmgStrain)
@@ -206,7 +200,7 @@
/* uniaxial tension resistance, bending parameter of the damage evolution law, whear weighting constant for epsT in the strain seminorm (kappa) calculation. Default to NaN so that user gets loudly notified it was not set.
*/
- Real sigmaT, epsCrackOnset, relDuctility, G_over_E, tau, expDmgRate, omegaThreshold, dmgTau, dmgRateExp, plTau, plRateExp, isoPrestress;
+ Real sigmaT, epsCrackOnset, relDuctility, G_over_E, tau, expDmgRate, dmgTau, dmgRateExp, plTau, plRateExp, isoPrestress;
//! Should new contacts be cohesive? They will before this iter#, they will not be afterwards. If 0, they will never be. If negative, they will always be created as cohesive.
long cohesiveThresholdIter;
//! Create contacts that don't receive any damage (CpmPhys::neverDamage=true); defaults to false
@@ -218,7 +212,6 @@
neverDamage=false;
cohesiveThresholdIter=-1;
dmgTau=-1; dmgRateExp=0; plTau=-1; plRateExp=-1;
- omegaThreshold=0.999;
isoPrestress=0;
}
@@ -234,7 +227,6 @@
(dmgRateExp)
(plTau)
(plRateExp)
- (omegaThreshold)
(isoPrestress)
);
@@ -258,14 +250,17 @@
int yieldSurfType;
//! scaling in the logarithmic yield surface (should be <1 for realistic results; >=0 for meaningful results)
static Real yieldLogSpeed;
+ //! horizontal scaling of the ellipse (shifts on the +x axis as interactions with +y are given)
static Real yieldEllipseShift;
+ //! damage after which the contact disappears (<1), since omega reaches 1 only for strain →+∞
+ static Real omegaThreshold;
//! HACK: limit strain on some contacts by moving body #2 in the contact; only if refR1<0 (facet); deactivated if > 0
static Real minStrain_moveBody2;
Law2_Dem3DofGeom_CpmPhys_Cpm(): logStrain(false), yieldSurfType(0) { /*timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);*/ }
void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody);
FUNCTOR2D(Dem3DofGeom,CpmPhys);
REGISTER_CLASS_AND_BASE(Law2_Dem3DofGeom_CpmPhys_Cpm,ConstitutiveLaw);
- REGISTER_ATTRIBUTES(ConstitutiveLaw,(logStrain)(yieldSurfType)(yieldLogSpeed)(yieldEllipseShift)(minStrain_moveBody2));
+ REGISTER_ATTRIBUTES(ConstitutiveLaw,(logStrain)(yieldSurfType)(yieldLogSpeed)(yieldEllipseShift)(minStrain_moveBody2)(omegaThreshold));
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_CpmPhys_Cpm);
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -48,7 +48,7 @@
Vector3r normal=facet->nf;
Real planeDist=normal.Dot(cogLine);
if(planeDist<0){normal*=-1; planeDist*=-1; }
- if(planeDist>sphereRadius && !c->isReal) { /* LOG_TRACE("Sphere too far ("<<planeDist<<") from plane"); */ return false; }
+ if(planeDist>sphereRadius && !c->isReal()) { /* LOG_TRACE("Sphere too far ("<<planeDist<<") from plane"); */ return false; }
Vector3r planarPt=cogLine-planeDist*normal; // project sphere center to the facet plane
Real normDotPt[3];
Vector3r contactPt(Vector3r::ZERO);
@@ -70,7 +70,7 @@
}
normal=cogLine-contactPt; // called normal, but it is no longer the facet's normal (for compat)
//TRVAR3(normal,contactPt,sphereRadius);
- if(!c->isReal && normal.SquaredLength()>sphereRadius*sphereRadius) { /* LOG_TRACE("Sphere too far from closest point"); */ return false; } // fast test before sqrt
+ if(!c->isReal() && normal.SquaredLength()>sphereRadius*sphereRadius) { /* LOG_TRACE("Sphere too far from closest point"); */ return false; } // fast test before sqrt
Real penetrationDepth=sphereRadius-normal.Normalize();
#else
/* This code was mostly copied from InteractingFacet2InteractinSphere4SpheresContactGeometry */
@@ -79,7 +79,7 @@
Vector3r normal=facet->nf;
Real L=normal.Dot(contactLine); // height/depth of sphere's center from facet's plane
if(L<0){normal*=-1; L*=-1;}
- if(L>sphereRadius && !c->isReal) return false; // sphere too far away from the plane
+ if(L>sphereRadius && !c->isReal()) return false; // sphere too far away from the plane
Vector3r contactPt=contactLine-L*normal; // projection of sphere's center to facet's plane (preliminary contact point)
const Vector3r* edgeNormals=facet->ne; // array[3] of edge normals (in facet plane)
@@ -130,7 +130,7 @@
// end facet-local coordinates
#endif
- if(penetrationDepth<0 && !c->isReal) return false;
+ if(penetrationDepth<0 && !c->isReal()) return false;
shared_ptr<Dem3DofGeom_FacetSphere> fs;
Vector3r normalGlob=se31.orientation*normal;
@@ -149,16 +149,6 @@
fs->se31=se31; fs->se32=se32;
fs->normal=normalGlob;
fs->contactPoint=se32.position+(-normalGlob)*(sphereRadius-penetrationDepth);
- if(c->isNew){
- //TRVAR2(planeDist,planarPt);
- //TRVAR3(normDotPt[0],normDotPt[1],normDotPt[2]);
- //TRVAR2(w,contactPt);
- TRVAR1(penetrationDepth);
- TRVAR3(fs->refLength,fs->cp1pt,fs->localFacetNormal);
- TRVAR3(fs->effR2,fs->cp2rel,fs->normal);
- TRVAR2(fs->se31.orientation,fs->se32.orientation);
- TRVAR2(fs->contPtInTgPlane1(),fs->contPtInTgPlane2());
- }
return true;
}
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -152,7 +152,7 @@
InteractingSphere *s1=static_cast<InteractingSphere*>(cm1.get()), *s2=static_cast<InteractingSphere*>(cm2.get());
Vector3r normal=se32.position-se31.position;
Real penetrationDepthSq=pow(distanceFactor*(s1->radius+s2->radius),2)-normal.SquaredLength();
- if (penetrationDepthSq<0 && !c->isReal) return false;
+ if (penetrationDepthSq<0 && !c->isReal()) return false;
Real dist=normal.Normalize(); /* Normalize() works in-place and returns length before normalization; from here, normal is unit vector */
shared_ptr<Dem3DofGeom_SphereSphere> ss;
Modified: trunk/pkg/dem/Engine/DeusExMachina/CapillaryRecorder.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/CapillaryRecorder.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/DeusExMachina/CapillaryRecorder.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -84,7 +84,7 @@
for( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal )
+ if ((*ii)->isReal() )
{
const shared_ptr<Interaction>& interaction = *ii;
Modified: trunk/pkg/dem/Engine/DeusExMachina/CapillaryStressRecorder.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/CapillaryStressRecorder.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/DeusExMachina/CapillaryStressRecorder.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -113,7 +113,7 @@
for( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
const shared_ptr<Interaction>& interaction = *ii;
Modified: trunk/pkg/dem/Engine/DeusExMachina/ContactStressRecorder.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/ContactStressRecorder.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/DeusExMachina/ContactStressRecorder.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -119,7 +119,7 @@
for( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
const shared_ptr<Interaction>& interaction = *ii;
Modified: trunk/pkg/dem/Engine/DeusExMachina/ThreeDTriaxialEngine.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/ThreeDTriaxialEngine.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/DeusExMachina/ThreeDTriaxialEngine.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -238,7 +238,7 @@
InteractionContainer::iterator ii = ncb->transientInteractions->begin();
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for( ; ii!=iiEnd ; ++ii ) {
- if (!(*ii)->isReal) continue;
+ if (!(*ii)->isReal()) continue;
const shared_ptr<BodyMacroParameters>& sdec1 = YADE_PTR_CAST<BodyMacroParameters>((*bodies)[(body_id_t) ((*ii)->getId1())]->physicalParameters);
const shared_ptr<BodyMacroParameters>& sdec2 = YADE_PTR_CAST<BodyMacroParameters>((*bodies)[(body_id_t) ((*ii)->getId2())]->physicalParameters);
//FIXME - why dynamic_cast fails here?
Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -36,6 +36,7 @@
previousState=currentState;
UnbalancedForce = 1;
Key = "";
+ noFiles=false;
Phase1End = "Compacted";
FinalIterationPhase1 = 0;
Iteration = 0;
@@ -89,6 +90,7 @@
REGISTER_ATTRIBUTE(previousSigmaIso);
REGISTER_ATTRIBUTE(sigmaLateralConfinement);
REGISTER_ATTRIBUTE(Key);
+ REGISTER_ATTRIBUTE(noFiles);
REGISTER_ATTRIBUTE(frictionAngleDegree);
REGISTER_ATTRIBUTE(epsilonMax);
REGISTER_ATTRIBUTE(uniaxialEpsilonCurr);
@@ -112,8 +114,8 @@
height0 = height; depth0 = depth; width0 = width;
//compressionActivated = true;
wall_bottom_activated=false; wall_top_activated=false;
- if(currentState==STATE_ISO_UNLOADING){ LOG_INFO("Speres -> /tmp/unloaded.spheres"); Shop::saveSpheresToFile("/tmp/unloaded.spheres"); }
- if(!firstRun) saveSimulation=true; // saving snapshot .xml will actually be done in ::applyCondition
+ if(currentState==STATE_ISO_UNLOADING && !noFiles){ LOG_INFO("Speres -> /tmp/unloaded.spheres"); Shop::saveSpheresToFile("/tmp/unloaded.spheres"); }
+ if(!firstRun && !noFiles) saveSimulation=true; // saving snapshot .xml will actually be done in ::applyCondition
Phase1End = "Unloaded";
}
else if(currentState==STATE_ISO_COMPACTION && nextState==STATE_ISO_UNLOADING){
@@ -122,17 +124,17 @@
previousSigmaIso=sigma_iso;
internalCompaction=false; // unloading will not change grain sizes
if (frictionAngleDegree>0) setContactProperties(body, frictionAngleDegree);
- if(!firstRun) saveSimulation=true;
+ if(!firstRun && !noFiles) saveSimulation=true;
Phase1End = "Compacted";
}
else if ((currentState==STATE_ISO_COMPACTION || currentState==STATE_ISO_UNLOADING) && nextState==STATE_LIMBO) {
//urrentState==STATE_DIE_COMPACTION
internalCompaction = false;
height0 = height; depth0 = depth; width0 = width;
- saveSimulation=true; // saving snapshot .xml will actually be done in ::applyCondition
+ if(!noFiles) saveSimulation=true; // saving snapshot .xml will actually be done in ::applyCondition
// stop simulation here, since nothing will happen from now on
Phase1End = (currentState==STATE_ISO_COMPACTION ? "compacted" : "unloaded");
- Shop::saveSpheresToFile("/tmp/limbo.spheres");
+ if(!noFiles) Shop::saveSpheresToFile("/tmp/limbo.spheres");
}
else if( nextState==STATE_FIXED_POROSITY_COMPACTION){
internalCompaction = false;
@@ -239,16 +241,19 @@
if ( saveSimulation )
{
- string fileName = "./"+ Key + "_" + Phase1End + "_" +
- lexical_cast<string> ( Omega::instance().getCurrentIteration() ) + "_" +
- lexical_cast<string> ( currentState ) + ".xml";
- LOG_INFO ( "saving snapshot: "<<fileName );
- Omega::instance().saveSimulation ( fileName );
- fileName="./"+ Key + "_"+Phase1End+"_"+lexical_cast<string> ( Omega::instance().getCurrentIteration() ) + "_" +
- lexical_cast<string> ( currentState ) +".spheres";
- LOG_INFO ( "saving spheres: "<<fileName );
- Shop::saveSpheresToFile ( fileName );
+ if(!noFiles){
+ string fileName = "./"+ Key + "_" + Phase1End + "_" +
+ lexical_cast<string> ( Omega::instance().getCurrentIteration() ) + "_" +
+ lexical_cast<string> ( currentState ) + ".xml";
+ LOG_INFO ( "saving snapshot: "<<fileName );
+ Omega::instance().saveSimulation ( fileName );
+ fileName="./"+ Key + "_"+Phase1End+"_"+lexical_cast<string> ( Omega::instance().getCurrentIteration() ) + "_" +
+ lexical_cast<string> ( currentState ) +".spheres";
+ LOG_INFO ( "saving spheres: "<<fileName );
+ Shop::saveSpheresToFile ( fileName );
+ }
saveSimulation = false;
+
}
if ( Omega::instance().getCurrentIteration() % testEquilibriumInterval == 0 )
@@ -333,7 +338,7 @@
InteractionContainer::iterator ii = ncb->transientInteractions->begin();
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for( ; ii!=iiEnd ; ++ii ) {
- if (!(*ii)->isReal) continue;
+ if (!(*ii)->isReal()) continue;
const shared_ptr<BodyMacroParameters>& sdec1 = YADE_PTR_CAST<BodyMacroParameters>((*bodies)[(body_id_t) ((*ii)->getId1())]->physicalParameters);
const shared_ptr<BodyMacroParameters>& sdec2 = YADE_PTR_CAST<BodyMacroParameters>((*bodies)[(body_id_t) ((*ii)->getId2())]->physicalParameters);
//FIXME - why dynamic_cast fails here?
Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -102,6 +102,8 @@
int FinalIterationPhase1, Iteration, testEquilibriumInterval;
std::string Key;//A code that is appended to file names to help distinguish between different simulations
+ //! If true, no files will be generated (.xml, .spheres)
+ bool noFiles;
// //! Is uniaxial compression currently activated?
// bool compressionActivated;
//! Auto-switch from isotropic compaction or unloading state (if sigmaLateralConfinement<sigmaIsoCompaction)
Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialStateRecorder.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialStateRecorder.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialStateRecorder.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -202,7 +202,7 @@
for ( ; ii!=iiEnd ; ++ii )
{
- if ( ( *ii )->isReal )
+ if ( ( *ii )->isReal() )
{
const shared_ptr<Interaction>& interaction = *ii;
Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -154,7 +154,7 @@
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
const shared_ptr<Interaction>& contact = *ii;
@@ -382,7 +382,7 @@
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for ( ; ii!=iiEnd ; ++ii )
{
- if ( ( *ii )->isReal )
+ if ( ( *ii )->isReal() )
{
SpheresContactGeometry* contact = static_cast<SpheresContactGeometry*> ( ( *ii )->interactionGeometry.get() );
// if ((*(ncb->bodies))[(*ii)->getId1()]->isDynamic)
@@ -412,7 +412,7 @@
InteractionContainer::iterator ii = ncb->transientInteractions->begin();
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for( ; ii!=iiEnd ; ++ii ) {
- if ((*ii)->isReal) {
+ if ((*ii)->isReal()) {
const shared_ptr<Interaction>& contact = *ii;
Real f = (static_cast<ElasticContactInteraction*> ((contact->interactionPhysics.get()))->normalForce+static_cast<ElasticContactInteraction*>(contact->interactionPhysics.get())->shearForce).SquaredLength();
if (f!=0)
Modified: trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -28,7 +28,7 @@
, const shared_ptr<PhysicalParameters>& b2 // SimpleViscoelasticBodyParameters
, const shared_ptr<Interaction>& interaction)
{
- if( !interaction->isNew ) return;
+ if(interaction->interactionPhysics) return;
SimpleViscoelasticBodyParameters* sdec1 = static_cast<SimpleViscoelasticBodyParameters*>(b1.get());
SimpleViscoelasticBodyParameters* sdec2 = static_cast<SimpleViscoelasticBodyParameters*>(b2.get());
Modified: trunk/pkg/dem/Engine/EngineUnit/CL1Relationships.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/CL1Relationships.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/CL1Relationships.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -62,7 +62,7 @@
if(interactionGeometry) // so it is SpheresContactGeometry - NON PERMANENT LINK
{
- if(interaction->isNew)
+ if(!interaction->interactionPhysics)
{
//std::cerr << " isNew, id1: " << interaction->getId1() << " id2: " << interaction->getId2() << "\n";
interaction->interactionPhysics = shared_ptr<ContactLaw1Interaction>(new ContactLaw1Interaction());
Modified: trunk/pkg/dem/Engine/EngineUnit/CohesiveFrictionalRelationships.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/CohesiveFrictionalRelationships.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/CohesiveFrictionalRelationships.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -64,7 +64,7 @@
if(interactionGeometry) // so it is SpheresContactGeometry - NON PERMANENT LINK
{
- if(interaction->isNew)
+ if(!interaction->interactionPhysics)
{
//std::cerr << " isNew, id1: " << interaction->getId1() << " id2: " << interaction->getId2() << "\n";
interaction->interactionPhysics = shared_ptr<CohesiveFrictionalContactInteraction>(new CohesiveFrictionalContactInteraction());
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometry.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -85,11 +85,12 @@
pt2 = se32.position-normal*s->radius;
shared_ptr<SpheresContactGeometry> scm;
- if (c->isNew) scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
+ bool isNew=!c->interactionGeometry;
+ if (isNew) scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
else scm = YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
#ifdef SCG_SHEAR
- if(c->isNew) { /* same as below */ scm->prevNormal=pt1-pt2; scm->prevNormal.Normalize(); }
+ if(isNew) { /* same as below */ scm->prevNormal=pt1-pt2; scm->prevNormal.Normalize(); }
else {scm->prevNormal=scm->normal;}
#endif
@@ -132,10 +133,11 @@
pt2=se32.position+cOnBox_sphere*s->radius;
shared_ptr<SpheresContactGeometry> scm;
- if (c->isNew) scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
+ bool isNew=!c->interactionGeometry;
+ if (isNew) scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
else scm = YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
#ifdef SCG_SHEAR
- if(c->isNew) { /* same as below */ scm->prevNormal=-cOnBox_sphere; }
+ if(isNew) { /* same as below */ scm->prevNormal=-cOnBox_sphere; }
else {scm->prevNormal=scm->normal;}
#endif
scm->contactPoint = 0.5*(pt1+pt2);
@@ -159,7 +161,6 @@
const Se3r& se32,
const shared_ptr<Interaction>& c)
{
- assert(c->isNew);
c->swapOrder();
return go(cm2,cm1,se32,se31,c);
}
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometryWater.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometryWater.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingBox2InteractingSphere4SpheresContactGeometryWater.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -86,7 +86,7 @@
// FIXME : remove those uncommented lines
shared_ptr<SpheresContactGeometry> scm;
- if (c->isNew)
+ if (!c->interactionGeometry)
scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
else
scm = dynamic_pointer_cast<SpheresContactGeometry>(c->interactionGeometry);
@@ -136,7 +136,7 @@
shared_ptr<SpheresContactGeometry> scm;
//cerr << "scm = " << scm << " | c = " << c << endl;
- if (c->isNew)
+ if (!c->interactionGeometry)
{ //cerr << "c->isNew";
scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
}
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -49,25 +49,8 @@
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
- contactFace=YADE_CAST<SpheresContactGeometry*>(c->interactionGeometry.get())->facetContactFace;
- // determinate contact on negative side: reverse quantities
- if(contactFace<0){normal=-normal; L=-L;}
- // indeterminate contact on negative side: set to -1; if edge, will be reset to 0 below
- else if (contactFace==0 && L<0) { normal=-normal; L=-L; contactFace=-1; }
- // indeterminate and is on the positive side: set to 1; if edge, will be reset to 0 below
- else if(contactFace==0 && L>0) {contactFace=1;}
- } else {
- 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
+ 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
Vector3r cp = cl - L*normal;
const Vector3r* ne = facet->ne;
@@ -116,7 +99,7 @@
// END everything in facet-local coordinates
//
- if (penetrationDepth>0 || c->isReal)
+ if (penetrationDepth>0 || c->isReal())
{
shared_ptr<SpheresContactGeometry> scm;
if (c->interactionGeometry)
@@ -146,7 +129,6 @@
const Se3r& se32,
const shared_ptr<Interaction>& c)
{
- assert(c->isNew);
c->swapOrder();
//LOG_WARN("Swapped interaction order for "<<c->getId2()<<"&"<<c->getId1());
return go(cm2,cm1,se32,se31,c);
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingBox4InteractionOfMyTetrahedron.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingBox4InteractionOfMyTetrahedron.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingBox4InteractionOfMyTetrahedron.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -74,7 +74,7 @@
shared_ptr<InteractionOfMyTetrahedron> imt;
// depending whether it's a new interaction: create new one, or use the existing one.
- if (c->isNew)
+ if (c->interactionGeometry)
imt = shared_ptr<InteractionOfMyTetrahedron>(new InteractionOfMyTetrahedron());
else
imt = YADE_PTR_CAST<InteractionOfMyTetrahedron>(c->interactionGeometry);
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingMyTetrahedron4InteractionOfMyTetrahedron.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingMyTetrahedron4InteractionOfMyTetrahedron.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingMyTetrahedron2InteractingMyTetrahedron4InteractionOfMyTetrahedron.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -36,7 +36,7 @@
shared_ptr<InteractionOfMyTetrahedron> imt;
// depending whether it's a new interaction: create new one, or use the existing one.
- if (c->isNew)
+ if (c->interactionGeometry)
imt = shared_ptr<InteractionOfMyTetrahedron>(new InteractionOfMyTetrahedron());
else
imt = YADE_PTR_CAST<InteractionOfMyTetrahedron>(c->interactionGeometry);
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2BssSweptSphereLineSegment4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2BssSweptSphereLineSegment4SpheresContactGeometry.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2BssSweptSphereLineSegment4SpheresContactGeometry.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -39,7 +39,7 @@
if (overlap <= 0.0)
{
shared_ptr<SpheresContactGeometry> scm;
- if (c->isNew) scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
+ if (c->interactionGeometry) scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
else scm = YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
scm->contactPoint = se32.position + proj + (ssls->radius+0.5*overlap)*ccn;
Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -25,13 +25,14 @@
InteractingSphere *s1=static_cast<InteractingSphere*>(cm1.get()), *s2=static_cast<InteractingSphere*>(cm2.get());
Vector3r normal=se32.position-se31.position;
Real penetrationDepthSq=pow(interactionDetectionFactor*(s1->radius+s2->radius),2) - normal.SquaredLength();
- if (penetrationDepthSq>0 || c->isReal){
+ if (penetrationDepthSq>0 || c->isReal()){
shared_ptr<SpheresContactGeometry> scm;
+ bool isNew=c->interactionGeometry;
if(c->interactionGeometry) scm=YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
else { scm=shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry()); c->interactionGeometry=scm; }
#ifdef SCG_SHEAR
- if(c->isNew) scm->prevNormal=normal;
+ if(isNew) scm->prevNormal=normal;
else scm->prevNormal=scm->normal;
#endif
Modified: trunk/pkg/dem/Engine/EngineUnit/MacroMicroElasticRelationships.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/MacroMicroElasticRelationships.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/MacroMicroElasticRelationships.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -67,7 +67,7 @@
contactPhysics->ks = contactPhysics->initialKs;
contactPhysics->equilibriumDistance = contactPhysics->initialEquilibriumDistance;
*/
- if( interaction->isNew)
+ if(!interaction->interactionPhysics)
{
interaction->interactionPhysics = shared_ptr<ElasticContactInteraction>(new ElasticContactInteraction());
ElasticContactInteraction* contactPhysics = YADE_CAST<ElasticContactInteraction*>(interaction->interactionPhysics.get());
Modified: trunk/pkg/dem/Engine/EngineUnit/MacroMicroElasticRelationshipsWater.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/MacroMicroElasticRelationshipsWater.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/MacroMicroElasticRelationshipsWater.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -68,7 +68,7 @@
contactPhysics->ks = contactPhysics->initialKs;
contactPhysics->equilibriumDistance = contactPhysics->initialEquilibriumDistance;
*/
- if( interaction->isNew)
+ if(!interaction->interactionPhysics)
{
interaction->interactionPhysics = shared_ptr<CapillaryParameters>(new CapillaryParameters());
CapillaryParameters* contactPhysics = static_cast<CapillaryParameters*>(interaction->interactionPhysics.get());
Modified: trunk/pkg/dem/Engine/EngineUnit/SimpleElasticRelationships.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/SimpleElasticRelationships.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/SimpleElasticRelationships.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -25,7 +25,7 @@
//if(interactionGeometry)
{
- if( interaction->isNew)
+ if(!interaction->interactionPhysics)
{
const shared_ptr<BodyMacroParameters>& sdec1 = YADE_PTR_CAST<BodyMacroParameters>(b1);
const shared_ptr<BodyMacroParameters>& sdec2 = YADE_PTR_CAST<BodyMacroParameters>(b2);
Modified: trunk/pkg/dem/Engine/EngineUnit/SimpleElasticRelationshipsWater.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/SimpleElasticRelationshipsWater.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/SimpleElasticRelationshipsWater.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -40,7 +40,7 @@
if(interactionGeometry) // so it is SpheresContactGeometry - NON PERMANENT LINK
{
//cerr << "interactionGeometry" << endl;
- if( interaction->isNew)
+ if(!interaction->interactionPhysics)
{
//cerr << "interaction->isNew" << endl;
const shared_ptr<BodyMacroParameters>& sdec1 = YADE_PTR_CAST<BodyMacroParameters>(b1);
Modified: trunk/pkg/dem/Engine/EngineUnit/SimpleViscoelasticRelationships.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/SimpleViscoelasticRelationships.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/SimpleViscoelasticRelationships.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -28,7 +28,7 @@
, const shared_ptr<PhysicalParameters>& b2 // SimpleViscoelasticBodyParameters
, const shared_ptr<Interaction>& interaction)
{
- if( !interaction->isNew ) return;
+ if(interaction->interactionPhysics) return;
SimpleViscoelasticBodyParameters* sdec1 = static_cast<SimpleViscoelasticBodyParameters*>(b1.get());
SimpleViscoelasticBodyParameters* sdec2 = static_cast<SimpleViscoelasticBodyParameters*>(b2.get());
Modified: trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -26,8 +26,9 @@
RigidBodyParameters* de2 = YADE_CAST<RigidBodyParameters*>((*bodies)[id2]->physicalParameters.get());
Vector3r& shearForce = phys->shearForce;
- if (I->isNew) shearForce=Vector3r(0,0,0);
+ if (I->isFresh(rootBody)) shearForce=Vector3r(0,0,0);
+
Real dt = Omega::instance().getTimeStep();
Vector3r axis = phys->prevNormal.Cross(geom->normal);
Modified: trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -118,7 +118,7 @@
for( ; ii!=iiEnd ; ++ii ) {
- if ((*ii)->isReal) {//FIXME : test to be removed when using DistantPersistentSAPCollider?
+ if ((*ii)->isReal()) {//FIXME : test to be removed when using DistantPersistentSAPCollider?
const shared_ptr<Interaction>& interaction = *ii;
unsigned int id1 = interaction->getId1();
@@ -268,7 +268,7 @@
for(ii= ncb->transientInteractions->begin(); ii!=iiEnd ; ++ii )
{ //cerr << "interaction " << ii << endl;
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
CapillaryParameters* currentContactPhysics = static_cast<CapillaryParameters*>((*ii)->interactionPhysics.get());
if (currentContactPhysics->meniscus)
@@ -294,7 +294,7 @@
//cerr << "id1/id2 " << (*ii)->getId1() << "/" << (*ii)->getId2() << " Fcap= " << currentContactPhysics->Fcap << endl;
}
- }
+ }
}
@@ -317,7 +317,7 @@
//Reset fusion numbers
InteractionContainer::iterator ii = ncb->transientInteractions->begin();
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
- for( ; ii!=iiEnd ; ++ii ) if ((*ii)->isReal) static_cast<CapillaryParameters*>((*ii)->interactionPhysics.get())->fusionNumber=0;
+ for( ; ii!=iiEnd ; ++ii ) if ((*ii)->isReal()) static_cast<CapillaryParameters*>((*ii)->interactionPhysics.get())->fusionNumber=0;
list< shared_ptr<Interaction> >::iterator firstMeniscus, lastMeniscus, currentMeniscus;
@@ -655,7 +655,7 @@
InteractionContainer::iterator ii = ncb->transientInteractions->begin();
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for( ; ii!=iiEnd ; ++ii ) {
- if ((*ii)->isReal) {
+ if ((*ii)->isReal()) {
if (static_cast<CapillaryParameters*>((*ii)->interactionPhysics.get())->meniscus) insert(*ii);
}
}
Modified: trunk/pkg/dem/Engine/StandAloneEngine/CohesiveFrictionalContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/CohesiveFrictionalContactLaw.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/CohesiveFrictionalContactLaw.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -90,7 +90,7 @@
for ( ; ii!=iiEnd ; ++ii )
{
//if ((*ii)->interactionGeometry && (*ii)->interactionPhysics)
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
if (detectBrokenBodies
/* FIXME - this had no effect. InteractingBox has its isBroken=false too.
@@ -119,8 +119,7 @@
Vector3r& shearForce = currentContactPhysics->shearForce;
- if (contact->isNew)
- shearForce = Vector3r::ZERO;
+ if (contact->isFresh(ncb)) shearForce = Vector3r::ZERO;
Real un = currentContactGeometry->penetrationDepth;
Real Fn = currentContactPhysics->kn*un;
@@ -136,8 +135,9 @@
//if (currentContactPhysics->cohesionBroken) {
//cerr << "broken" << endl;
- contact->isReal= false;
- currentContactPhysics->cohesionBroken = true;
+ ncb->interactions->requestErase(contact->getId1(),contact->getId2());
+ // contact->interactionPhysics was reset now; currentContactPhysics still hold the object, but is not associated with the interaction anymore
+ currentContactPhysics->cohesionBroken = true;
currentContactPhysics->normalForce = Vector3r::ZERO;
currentContactPhysics->shearForce = Vector3r::ZERO;
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ContactLaw1.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ContactLaw1.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ContactLaw1.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -51,13 +51,13 @@
for ( ; ii!=iiEnd ; ++ii )
{
//if ((*ii)->interactionGeometry && (*ii)->interactionPhysics)
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
nbreInteracTot++;
const shared_ptr<Interaction>& contact = *ii;
int id1 = contact->getId1();
int id2 = contact->getId2();
-// cout << "contact entre " << id1 << " et " << id2 << " reel ? " << contact->isReal << endl;
+// cout << "contact entre " << id1 << " et " << id2 << " reel ? " << contact->isReal() << endl;
if ( !( (*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask) )
continue; // skip other groups,
@@ -68,7 +68,7 @@
Vector3r& shearForce = currentContactPhysics->shearForce;
- if (contact->isNew)
+ if (contact->isFresh(ncb))
{
shearForce = Vector3r::ZERO;
currentContactPhysics->previousun=0.0;
@@ -114,12 +114,13 @@
//currentContactPhysics->SetBreakingState();
- contact->isReal= false;
+ ncb->interactions->requestErase(contact->getId1(),contact->getId2());
+ // probably not useful anymore
currentContactPhysics->normalForce = Vector3r::ZERO;
currentContactPhysics->shearForce = Vector3r::ZERO;
//return;
- // } else
+ // else
// currentContactPhysics->normalForce = -currentContactPhysics->normalAdhesion*currentContactGeometry->normal;
}
else
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticCohesiveLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticCohesiveLaw.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticCohesiveLaw.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -49,7 +49,7 @@
InteractionContainer::iterator iiEnd = ncb->persistentInteractions->end();
for( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
const shared_ptr<Interaction> contact2 = *ii;
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -50,7 +50,7 @@
functor->useShear=useShear;
#endif
FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){
- if(!I->isReal) continue;
+ if(!I->isReal()) continue;
#ifdef YADE_DEBUG
// these checks would be redundant in the functor (ConstitutiveLawDispatcher does that already)
if(!dynamic_cast<SpheresContactGeometry*>(I->interactionGeometry.get()) || !dynamic_cast<ElasticContactInteraction*>(I->interactionPhysics.get())) continue;
@@ -77,7 +77,7 @@
Vector3r& shearForce = currentContactPhysics->shearForce;
- if (contact->isNew) shearForce=Vector3r(0,0,0);
+ if (contact->isFresh(ncb)) shearForce=Vector3r(0,0,0);
Real un=currentContactGeometry->penetrationDepth;
TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position);
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ForceSnapshot.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ForceSnapshot.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ForceSnapshot.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -85,7 +85,7 @@
for( ;ii!=iiEnd;++ii)
{
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
const shared_ptr<Interaction>& contact = *ii;
body_id_t id1 = contact->getId1();
Modified: trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -189,7 +189,7 @@
memset(stiffnesses[0], 0,sizeof(Vector3r)*size);
memset(Rstiffnesses[0],0,sizeof(Vector3r)*size);
FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
- if(!contact->isReal) continue;
+ if(!contact->isReal()) continue;
SpheresContactGeometry* geom=YADE_CAST<SpheresContactGeometry*>(contact->interactionGeometry.get()); assert(geom);
NormalShearInteraction* phys=YADE_CAST<NormalShearInteraction*>(contact->interactionPhysics.get()); assert(phys);
Modified: trunk/pkg/dem/Engine/StandAloneEngine/MicroMacroAnalyser.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/MicroMacroAnalyser.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/MicroMacroAnalyser.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -156,7 +156,7 @@
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for ( ; ii!=iiEnd ; ++ii )
{
- if ( ( *ii )->isReal )
+ if ( ( *ii )->isReal() )
{
TriaxialState::Contact *c = new TriaxialState::Contact;
TS.contacts.push_back ( c );
Modified: trunk/pkg/dem/Engine/StandAloneEngine/MyTetrahedronLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/MyTetrahedronLaw.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/MyTetrahedronLaw.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -35,7 +35,7 @@
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal) // isReal means that InteractingMyTetrahedron2InteractingMyTetrahedron4InteractionOfMyTetrahedron returned true
+ if ((*ii)->isReal()) // isReal means that InteractingMyTetrahedron2InteractingMyTetrahedron4InteractionOfMyTetrahedron returned true
// or InteractingMyTetrahedron2InteractingBox4InteractionOfMyTetrahedron returned true
{
const shared_ptr<Interaction>& contact = *ii;
@@ -76,7 +76,6 @@
}
}
-
}
}
}
Modified: trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticSpheresInteractionRecorder.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticSpheresInteractionRecorder.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticSpheresInteractionRecorder.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -53,7 +53,7 @@
for( ; ii!=iiEnd ; ++ii )
{
const shared_ptr<Interaction>& i = *ii;
- if ( !i->isReal ) continue;
+ if ( !i->isReal() ) continue;
if ( i->interactionGeometry->getClassIndex() != interactionSphere->getClassIndex() ) continue;
if ( i->interactionPhysics->getClassIndex() != viscoelasticInteraction->getClassIndex() ) continue;
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ViscousForceDamping.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ViscousForceDamping.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ViscousForceDamping.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -51,7 +51,7 @@
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
const shared_ptr<Interaction>& contact = *ii;
int id1 = contact->getId1();
@@ -67,7 +67,7 @@
Vector3r& shearForce = currentContactPhysics->shearForce;
- if ( contact->isNew)
+ if ( contact->isFresh(ncb))
shearForce = Vector3r(0,0,0);
Real un = currentContactGeometry->penetrationDepth;
Modified: trunk/pkg/dem/Engine/StandAloneEngine/VolumicContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/VolumicContactLaw.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/Engine/StandAloneEngine/VolumicContactLaw.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -209,7 +209,7 @@
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
const shared_ptr<Interaction>& contact = *ii;
int id1 = contact->getId1();
@@ -230,7 +230,7 @@
- if (contact->isNew) // Si on a une nouvelle interpenetration detectée, on initialise
+ if (contact->isFresh(ncb)) // Si on a une nouvelle interpenetration detectée, on initialise
{
shearForce = Vector3r(0,0,0);
currentContactPhysics->previousun=0.0;
@@ -478,6 +478,7 @@
ncb->bex.addTorque(id2, c2x.Cross(f));
currentContactPhysics->prevNormal = currentContactGeometry->normal;
+
}
}
Modified: trunk/pkg/dem/PreProcessor/HydraulicTest.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/HydraulicTest.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/PreProcessor/HydraulicTest.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -23,7 +23,6 @@
#include<yade/pkg-common/AABB.hpp>
#include<yade/pkg-common/Sphere.hpp>
#include<yade/core/MetaBody.hpp>
-#include<yade/pkg-common/DistantPersistentSAPCollider.hpp>
#include<yade/lib-serialization/IOFormatManager.hpp>
#include<yade/core/Interaction.hpp>
#include<yade/pkg-common/BoundingVolumeMetaEngine.hpp>
Modified: trunk/pkg/dem/PreProcessor/SDECLinkedSpheres.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/SDECLinkedSpheres.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/PreProcessor/SDECLinkedSpheres.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -192,9 +192,6 @@
link->interactionGeometry = geometry;
link->interactionPhysics = physics;
- link->isReal = true;
- link->isNew = false;
-
rootBody->persistentInteractions->insert(link);
}
}
Modified: trunk/pkg/dem/PreProcessor/ThreePointBending.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/ThreePointBending.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/PreProcessor/ThreePointBending.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -200,9 +200,6 @@
link->interactionGeometry = geometry;
link->interactionPhysics = physics;
- link->isReal = true;
- link->isNew = false;
-
rootBody->persistentInteractions->insert(link);
}
}
Modified: trunk/pkg/dem/PreProcessor/TriaxialTest.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/TriaxialTest.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/PreProcessor/TriaxialTest.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -34,6 +34,7 @@
#include<yade/pkg-common/Sphere.hpp>
#include<yade/core/MetaBody.hpp>
#include<yade/pkg-common/PersistentSAPCollider.hpp>
+#include<yade/pkg-common/InsertionSortCollider.hpp>
#include<yade/lib-serialization/IOFormatManager.hpp>
#include<yade/core/Interaction.hpp>
#include<yade/pkg-common/BoundingVolumeMetaEngine.hpp>
@@ -159,7 +160,8 @@
isotropicCompaction=false;
fixedPorosity = 1;
- parallel=false;
+ fast=false;
+ noFiles=false;
@@ -230,7 +232,8 @@
REGISTER_ATTRIBUTE(isotropicCompaction);
REGISTER_ATTRIBUTE(fixedPorosity);
REGISTER_ATTRIBUTE(fixedBoxDims);
- REGISTER_ATTRIBUTE(parallel);
+ REGISTER_ATTRIBUTE(fast);
+ REGISTER_ATTRIBUTE(noFiles);
}
@@ -267,8 +270,7 @@
if(importFilename==""){
Vector3r dimensions=upperCorner-lowerCorner; Real volume=dimensions.X()*dimensions.Y()*dimensions.Z();
- Real really_radiusMean;
- if(radiusMean<=0) really_radiusMean=pow(volume*(1-porosity)/(Mathr::PI*(4/3.)*numberOfGrains),1/3.);
+ if(radiusMean<=0) radiusMean=pow(volume*(1-porosity)/(Mathr::PI*(4/3.)*numberOfGrains),1/3.);
else {
bool fixedDims[3];
fixedDims[0]=fixedBoxDims.find('x')!=string::npos; fixedDims[1]=fixedBoxDims.find('y')!=string::npos; fixedDims[2]=fixedBoxDims.find('z')!=string::npos;
@@ -278,9 +280,8 @@
LOG_INFO("Mean radius value of "<<radiusMean<<" requested, scaling "<<nScaled<<" dimensions by "<<boxScaleFactor);
dimensions[0]*=fixedDims[0]?1.:boxScaleFactor; dimensions[1]*=fixedDims[1]?1.:boxScaleFactor; dimensions[2]*=fixedDims[2]?1.:boxScaleFactor;
upperCorner=lowerCorner+dimensions;
- really_radiusMean=radiusMean;
}
- message+=GenerateCloud(sphere_list, lowerCorner, upperCorner, numberOfGrains, radiusStdDev, really_radiusMean, porosity);
+ message+=GenerateCloud(sphere_list, lowerCorner, upperCorner, numberOfGrains, radiusStdDev, radiusMean, porosity);
}
else {
if(radiusMean>0) LOG_WARN("radiusMean ignored, since importFilename specified.");
@@ -573,13 +574,14 @@
triaxialcompressionEngine->maxMultiplier = maxMultiplier;
triaxialcompressionEngine->finalMaxMultiplier = finalMaxMultiplier;
triaxialcompressionEngine->Key = Key;
+ triaxialcompressionEngine->noFiles=noFiles;
triaxialcompressionEngine->frictionAngleDegree = sphereFrictionDeg;
triaxialcompressionEngine->fixedPorosity = fixedPorosity;
triaxialcompressionEngine->isotropicCompaction = isotropicCompaction;
// recording global stress
- if(recordIntervalIter>0){
+ if(recordIntervalIter>0 && !noFiles){
triaxialStateRecorder = shared_ptr<TriaxialStateRecorder>(new TriaxialStateRecorder);
triaxialStateRecorder-> outputFile = WallStressRecordFile + Key;
triaxialStateRecorder-> interval = recordIntervalIter;
@@ -604,8 +606,9 @@
rootBody->engines.clear();
rootBody->engines.push_back(shared_ptr<Engine>(new PhysicalActionContainerReseter));
rootBody->engines.push_back(boundingVolumeDispatcher);
- rootBody->engines.push_back(shared_ptr<Engine>(new PersistentSAPCollider));
- if(parallel){
+ if(!fast) rootBody->engines.push_back(shared_ptr<Engine>(new PersistentSAPCollider));
+ else rootBody->engines.push_back(shared_ptr<Engine>(new InsertionSortCollider));
+ if(fast){
shared_ptr<InteractionDispatchers> ids(new InteractionDispatchers);
ids->geomDispatcher=interactionGeometryDispatcher;
ids->physDispatcher=interactionPhysicsDispatcher;
@@ -626,7 +629,7 @@
//rootBody->engines.push_back(stiffnessMatrixTimeStepper);
rootBody->engines.push_back(globalStiffnessTimeStepper);
rootBody->engines.push_back(triaxialcompressionEngine);
- if(recordIntervalIter>0) rootBody->engines.push_back(triaxialStateRecorder);
+ if(recordIntervalIter>0 && !noFiles) rootBody->engines.push_back(triaxialStateRecorder);
//rootBody->engines.push_back(gravityCondition);
shared_ptr<NewtonsDampedLaw> newton(new NewtonsDampedLaw);
@@ -714,7 +717,7 @@
" tries while generating sphere number " +
lexical_cast<string>(i+1) + "/" + lexical_cast<string>(number) + ".";
}
- return "Generated a sample with " + lexical_cast<string>(number) + "spheres inside box of dimensions: ("
+ return "Generated a sample with " + lexical_cast<string>(number) + " spheres inside box of dimensions: ("
+ lexical_cast<string>(dimensions[0]) + ","
+ lexical_cast<string>(dimensions[1]) + ","
+ lexical_cast<string>(dimensions[2]) + ").";
Modified: trunk/pkg/dem/PreProcessor/TriaxialTest.hpp
===================================================================
--- trunk/pkg/dem/PreProcessor/TriaxialTest.hpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/PreProcessor/TriaxialTest.hpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -104,10 +104,13 @@
,saveAnimationSnapshots
,biaxial2dTest
//!flag to choose an isotropic compaction until a fixed porosity choosing a same translation speed for the six walls
- ,isotropicCompaction;
+ ,isotropicCompaction
+ //! do not create any files during run (.xml, .spheres, wall stress records)
+ ,noFiles
+ ;
- //! Generate parallel simulation
- bool parallel;
+ //! Generate faster simulation: use InsertionSortCollider and InteractionDispatchers
+ bool fast;
int recordIntervalIter
,timeStepUpdateInterval
Modified: trunk/pkg/dem/RenderingEngine/GLDrawCohesiveFrictionalContactInteraction/GLDrawCohesiveFrictionalContactInteraction.cpp
===================================================================
--- trunk/pkg/dem/RenderingEngine/GLDrawCohesiveFrictionalContactInteraction/GLDrawCohesiveFrictionalContactInteraction.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/RenderingEngine/GLDrawCohesiveFrictionalContactInteraction/GLDrawCohesiveFrictionalContactInteraction.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -96,8 +96,8 @@
const shared_ptr<Body>& b2,
bool wireFrame)
{
- if(!i->isReal) return;
- isReal=i->isReal;
+ if(!i->isReal()) return;
+ isReal=i->isReal();
CohesiveFrictionalContactInteraction* ph = static_cast<CohesiveFrictionalContactInteraction*>(ih.get());
SpheresContactGeometry* sc = static_cast<SpheresContactGeometry*>(i->interactionGeometry.get());
Modified: trunk/pkg/dem/RenderingEngine/GLDrawElasticContactInteraction/GLDrawElasticContactInteraction.cpp
===================================================================
--- trunk/pkg/dem/RenderingEngine/GLDrawElasticContactInteraction/GLDrawElasticContactInteraction.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/RenderingEngine/GLDrawElasticContactInteraction/GLDrawElasticContactInteraction.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -44,7 +44,7 @@
const shared_ptr<Body>& b2,
bool wireFrame)
{
-// if(!i->isReal) return;
+// if(!i->isReal()) return;
ElasticContactInteraction* ph = static_cast<ElasticContactInteraction*>(ih.get());
SpheresContactGeometry* sc = static_cast<SpheresContactGeometry*>(i->interactionGeometry.get());
@@ -76,13 +76,13 @@
glEnd();
// draw normal
- drawArrow(cp, cp+normal*size*0.9 ,Vector3r(0,i->isReal?1:0.4,0));
+ drawArrow(cp, cp+normal*size*0.9 ,Vector3r(0,i->isReal()?1:0.4,0));
// draw prevNormal
-// drawArrow(middle, middle+ph->prevNormal*size*0.9 ,Vector3r(i->isReal?1:0.4,0,0));
+// drawArrow(middle, middle+ph->prevNormal*size*0.9 ,Vector3r(i->isReal()?1:0.4,0,0));
// draw shearForce
maxLength = std::max(maxLength,ph->shearForce.Length());
if(wireFrame) maxLength = 0.0000001;
- drawArrow(cp, cp+ph->shearForce*size*10.0/maxLength ,Vector3r(0,0,i->isReal?1:0.4));
+ drawArrow(cp, cp+ph->shearForce*size*10.0/maxLength ,Vector3r(0,0,i->isReal()?1:0.4));
// write A,B
// drawFlatText(pos1d,std::string(" A ") + boost::lexical_cast<std::string>(b1->getId()));
Modified: trunk/pkg/dem/RenderingEngine/GLDrawSDECLinkGeometry/GLDrawSDECLinkGeometry.cpp
===================================================================
--- trunk/pkg/dem/RenderingEngine/GLDrawSDECLinkGeometry/GLDrawSDECLinkGeometry.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/RenderingEngine/GLDrawSDECLinkGeometry/GLDrawSDECLinkGeometry.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -59,7 +59,7 @@
glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
// FIXME - we need a way to give parameters from outside, again.... so curerntly this scale is hardcoded here
- if( (!ip->isNew) && /*ip->isReal &&*/ ip->interactionPhysics)
+ if(ip->interactionPhysics)
{
Real force = el->normalForce.Length()/600;
forceMax = std::max(force,forceMax);
Modified: trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -62,7 +62,7 @@
// BUG: would crash with anything else than ElasticContactInteraction; use GLDrawInteractionPhysics for such things
/// ElasticContactInteraction* el = static_cast<ElasticContactInteraction*>(ip->interactionPhysics.get());
// FIXME - we need a way to give parameters from outside, again.... so curerntly this scale is hardcoded here
- if( (!ip->isNew) && ip->isReal && ip->interactionPhysics){
+ if( (!ip->isNew) && ip->isReal() && ip->interactionPhysics){
Real force = el->normalForce.Length()/600;
forceMax = std::max(force,forceMax);
Real scale = midMax*(force/forceMax)*0.3;
Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/dem/SConscript 2009-05-29 14:35:22 UTC (rev 1786)
@@ -702,6 +702,7 @@
'Sphere',
'AABB',
'PersistentSAPCollider',
+ 'InsertionSortCollider',
'MetaInteractingGeometry2AABB',
'TriaxialStressController',
'TriaxialCompressionEngine',
Modified: trunk/pkg/lattice/Engine/StandAloneEngine/NonLocalInitializer.cpp
===================================================================
--- trunk/pkg/lattice/Engine/StandAloneEngine/NonLocalInitializer.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/lattice/Engine/StandAloneEngine/NonLocalInitializer.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -106,8 +106,6 @@
nonLocal->gaussValue = std::exp( - std::pow( dist / len , 2) ) / ( len * sqPi );
- interaction->isReal = true;
- interaction->isNew = false;
interaction->interactionPhysics = nonLocal;
nonl->insert(interaction);
*/
Modified: trunk/pkg/lattice/PreProcessor/LatticeExample.cpp
===================================================================
--- trunk/pkg/lattice/PreProcessor/LatticeExample.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/lattice/PreProcessor/LatticeExample.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -1108,8 +1108,6 @@
if( dir2.Dot(beam2->otherDirection) < 0.999999 )
angularSpring->initialOffPlaneAngle2 *= -1.0;//, angularSpring->offPlaneSwap2 = true;
- interaction->isReal = true;
- interaction->isNew = false;
interaction->interactionPhysics = angularSpring;
ints->insert(interaction);
}
Modified: trunk/pkg/lattice/PreProcessor/LatticeExampleCTData.cpp
===================================================================
--- trunk/pkg/lattice/PreProcessor/LatticeExampleCTData.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/lattice/PreProcessor/LatticeExampleCTData.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -462,8 +462,6 @@
if( dir2.Dot(beam2->otherDirection) < 0.999999 )
angularSpring->initialOffPlaneAngle2 *= -1.0;//, angularSpring->offPlaneSwap2 = true;
- interaction->isReal = true;
- interaction->isNew = false;
interaction->interactionPhysics = angularSpring;
ints->insert(interaction);
}
Modified: trunk/pkg/lattice/PreProcessor/LatticeExampleSimple.cpp
===================================================================
--- trunk/pkg/lattice/PreProcessor/LatticeExampleSimple.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/lattice/PreProcessor/LatticeExampleSimple.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -414,8 +414,6 @@
if( dir2.Dot(beam2->otherDirection) < 0.999999 )
angularSpring->initialOffPlaneAngle2 *= -1.0;//, angularSpring->offPlaneSwap2 = true;
- interaction->isReal = true;
- interaction->isNew = false;
interaction->interactionPhysics = angularSpring;
ints->insert(interaction);
}
Modified: trunk/pkg/mass-spring/Engine/StandAloneEngine/MassSpringLaw.cpp
===================================================================
--- trunk/pkg/mass-spring/Engine/StandAloneEngine/MassSpringLaw.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/mass-spring/Engine/StandAloneEngine/MassSpringLaw.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -39,7 +39,7 @@
for( ; ii!=iiEnd ; ++ii )
{
shared_ptr<Interaction> spring = *ii;
- if (spring->isReal)
+ if (spring->isReal())
{
int id1 = spring->getId1();
int id2 = spring->getId2();
Modified: trunk/pkg/mass-spring/PreProcessor/HangingCloth.cpp
===================================================================
--- trunk/pkg/mass-spring/PreProcessor/HangingCloth.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/mass-spring/PreProcessor/HangingCloth.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -415,9 +415,6 @@
link->interactionGeometry = geometry;
link->interactionPhysics = physics;
- link->isReal = true;
- link->isNew = false;
-
rootBody->persistentInteractions->insert(link);
++linksNum;
}
@@ -453,9 +450,6 @@
// spring->interactionGeometry = geometry;
spring->interactionPhysics = physics;
- spring->isReal = true;
- spring->isNew = false;
-
return spring;
}
Modified: trunk/pkg/realtime-rigidbody/Engine/StandAloneEngine/FrictionLessElasticContactLaw.cpp
===================================================================
--- trunk/pkg/realtime-rigidbody/Engine/StandAloneEngine/FrictionLessElasticContactLaw.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/realtime-rigidbody/Engine/StandAloneEngine/FrictionLessElasticContactLaw.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -44,7 +44,7 @@
for( ; ii!=iiEnd; ++ii )
{
shared_ptr<Interaction> contact = *ii;
- if (contact->isReal)
+ if (contact->isReal())
{
int id1 = contact->getId1();
int id2 = contact->getId2();
Modified: trunk/pkg/snow/Engine/Ef2_BssSnowGrain_BssSnowGrain_makeSpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/snow/Engine/Ef2_BssSnowGrain_BssSnowGrain_makeSpheresContactGeometry.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/snow/Engine/Ef2_BssSnowGrain_BssSnowGrain_makeSpheresContactGeometry.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -29,7 +29,7 @@
BssSnowGrain* s1=dynamic_cast<BssSnowGrain*>(cm1.get()), *s2=dynamic_cast<BssSnowGrain*>(cm2.get());
Vector3r normal=se32.position-se31.position;
Real penetrationDepthSq=pow((s1->radius+s2->radius),2) - normal.SquaredLength();
- if (penetrationDepthSq>0 || c->isReal || assist)
+ if (penetrationDepthSq>0 || c->isReal() || assist)
{
shared_ptr<SpheresContactGeometry> scm;
if(c->interactionGeometry) scm=dynamic_pointer_cast<SpheresContactGeometry>(c->interactionGeometry);
Modified: trunk/pkg/snow/Engine/ElawSnowLayersDeformation.cpp
===================================================================
--- trunk/pkg/snow/Engine/ElawSnowLayersDeformation.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/snow/Engine/ElawSnowLayersDeformation.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -41,7 +41,7 @@
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
for ( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal)
+ if ((*ii)->isReal())
{
const shared_ptr<Interaction>& contact = *ii;
int id1 = contact->getId1();
@@ -175,7 +175,6 @@
// }
// }
}
-
}
}
}
Modified: trunk/pkg/snow/RenderingEngine/Ef1_IstSnowLayersContact_glDraw.cpp
===================================================================
--- trunk/pkg/snow/RenderingEngine/Ef1_IstSnowLayersContact_glDraw.cpp 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/pkg/snow/RenderingEngine/Ef1_IstSnowLayersContact_glDraw.cpp 2009-05-29 14:35:22 UTC (rev 1786)
@@ -99,7 +99,7 @@
const shared_ptr<Body>& b2,
bool wireFrame)
{
- if(!ip->isReal)
+ if(!ip->isReal())
return;
IstSnowLayersContact* sc = static_cast<IstSnowLayersContact*>(ig.get());
Modified: trunk/scripts/test/regular-sphere-pack.py
===================================================================
--- trunk/scripts/test/regular-sphere-pack.py 2009-05-29 12:14:18 UTC (rev 1785)
+++ trunk/scripts/test/regular-sphere-pack.py 2009-05-29 14:35:22 UTC (rev 1786)
@@ -1,3 +1,3 @@
O.bodies.append(utils.regularSphereOrthoPack([0,0,0],extents=[2,2,2],radius=.1,gap=.1,color=(1,0,1)))
-O.bodies.append(utils.regularSphereOrthoPack([0,0,4],extents=2,radius=.1,gap=.1,color=(0,1,0)))
+O.bodies.append(utils.regularSphereOrthoPack([0,0,4],extents=2,radius=.1,gap=.1,color=(0,1,0),velocity=[0,10,0]))
