yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02252
[Branch ~yade-dev/yade/trunk] Rev 1796: 1. virtual facets don't have a BoundingVolume (flag noBoundingVolume in utils.facet)
------------------------------------------------------------
revno: 1796
committer: Sergei D. <sega@think>
branch nick: trunk
timestamp: Thu 2009-11-19 09:44:50 +0300
message:
1. virtual facets don't have a BoundingVolume (flag noBoundingVolume in utils.facet)
2. new engine ResetRandomPosition (another name?). This engine reset the body position to a random position in specifed (by virtual facets) region.
3. example of using of ResetRandomPosition engine: scripts/test/ResetRandomPosition.py
added:
pkg/common/Engine/StandAloneEngine/ResetRandomPosition.cpp
pkg/common/Engine/StandAloneEngine/ResetRandomPosition.hpp
scripts/test/ResetRandomPosition.py
modified:
pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
py/utils.py
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.
=== modified file 'pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp'
--- pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp 2009-11-14 12:32:32 +0000
+++ pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp 2009-11-19 06:44:50 +0000
@@ -68,7 +68,7 @@
it=XX.begin(),et=XX.end(); it < et; ++it)
{
if (it->coord > bv.max[0]) break;
- if (!it->flags.isMin) continue;
+ if (!it->flags.isMin || !it->flags.hasBB) continue;
int offset = 3*it->id;
if (!(maxima[offset] < bv.min[0] ||
minima[offset+1] > bv.max[1] ||
=== added file 'pkg/common/Engine/StandAloneEngine/ResetRandomPosition.cpp'
--- pkg/common/Engine/StandAloneEngine/ResetRandomPosition.cpp 1970-01-01 00:00:00 +0000
+++ pkg/common/Engine/StandAloneEngine/ResetRandomPosition.cpp 2009-11-19 06:44:50 +0000
@@ -0,0 +1,176 @@
+/*************************************************************************
+* Copyright (C) 2009 by Sergei Dorofeenko *
+* sega@xxxxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#include<boost/random.hpp>
+#include<yade/core/Body.hpp>
+#include<yade/pkg-common/AABB.hpp>
+#include<yade/pkg-common/InteractingSphere.hpp>
+#include<yade/pkg-common/InteractingFacet.hpp>
+#include<yade/pkg-common/InteractionDispatchers.hpp>
+#include<yade/pkg-dem/BodyMacroParameters.hpp>
+#include"ResetRandomPosition.hpp"
+#include<sstream>
+
+YADE_PLUGIN((ResetRandomPosition));
+YADE_REQUIRE_FEATURE(PYTHON)
+CREATE_LOGGER(ResetRandomPosition);
+
+namespace {
+boost::variate_generator<boost::mt19937,boost::uniform_real<> >
+ randomUnit(boost::mt19937(),boost::uniform_real<>(0,1));
+boost::variate_generator<boost::mt19937,boost::uniform_real<> >
+ randomSymmetricUnit(boost::mt19937(),boost::uniform_real<>(-1,1));
+}
+
+ResetRandomPosition::ResetRandomPosition()
+{
+ factoryFacets.clear();
+ point=Vector3r(0,0,0);
+ normal=Vector3r(0,1,0);
+ volumeSection=false;
+ maxAttempts=20;
+ velocity=Vector3r(0,0,0);
+ velocityRange=Vector3r(0,0,0);
+ angularVelocity=Vector3r(0,0,0);
+ angularVelocityRange=Vector3r(0,0,0);
+ first_run = true;
+}
+
+ResetRandomPosition::~ResetRandomPosition()
+{
+
+}
+
+void ResetRandomPosition::action(MetaBody* ncb)
+{
+ if (first_run)
+ {
+ FOREACH(shared_ptr<Engine> eng, ncb->engines)
+ {
+ bI=dynamic_cast<Collider*>(eng.get());
+ if (bI) break;
+ }
+ if (!bI)
+ {
+ LOG_FATAL("Can't find Collider." );
+ return;
+ }
+ iGME=dynamic_cast<InteractionGeometryMetaEngine*>(ncb->engineByName("InteractionGeometryMetaEngine").get());
+ if (!iGME)
+ {
+ InteractionDispatchers* iDsp=dynamic_cast<InteractionDispatchers*>(ncb->engineByName("InteractionDispatchers").get());
+ if (!iDsp)
+ {
+ LOG_FATAL("Can't find nor InteractionGeometryMetaEngine nor InteractionDispatchers." );
+ return;
+ }
+ iGME=dynamic_cast<InteractionGeometryMetaEngine*>(iDsp->geomDispatcher.get());
+ if (!iGME)
+ {
+ LOG_FATAL("Can't find InteractionGeometryMetaEngine." );
+ return;
+ }
+ }
+ first_run=false;
+ randomFacet= shared_ptr<RandomInt>(new RandomInt(boost::minstd_rand(),boost::uniform_int<>(0,factoryFacets.size()-1)));
+ }
+
+ shiftedBodies.clear();
+
+ FOREACH(int id, subscribedBodies)
+ {
+ shared_ptr<Body> b = Body::byId(id);
+ RigidBodyParameters* rb = YADE_CAST<RigidBodyParameters*>(b->physicalParameters.get());
+ Vector3r& position = rb->se3.position;
+
+ if ( (position-point).Dot(normal) < 0 )
+ {
+ Vector3r backup_pos = position;
+
+ bool is_overlap;
+
+ for (int attempt=0; attempt<maxAttempts; ++attempt)
+ {
+ position = (volumeSection) ? generatePositionInVolume():generatePositionOnSurface();
+
+ const Real r = YADE_CAST<InteractingSphere*>(b->interactingGeometry.get())->radius;
+ BoundingVolume bv;
+ bv.min = Vector3r(position[0]-r, position[1]-r, position[2]-r);
+ bv.max = Vector3r(position[0]+r, position[1]+r, position[2]+r);
+
+ is_overlap=false;
+
+ // Test overlap with already shifted bodies
+ FOREACH(shared_ptr<Body> sb, shiftedBodies)
+ {
+ if (iGME->explicitAction(b,sb)->interactionGeometry)
+ {
+ is_overlap=true;
+ break;
+ }
+ }
+ if (is_overlap) continue; // new attempt
+
+ // Test overlap with other bodies
+ if (bI->probeBoundingVolume(bv))
+ {
+ for( unsigned int i=0, e=bI->probedBodies.size(); i<e; ++i)
+ {
+ if (iGME->explicitAction(b,Body::byId(bI->probedBodies[i]))->interactionGeometry)
+ {
+ is_overlap=true;
+ break;
+ }
+ }
+ if (is_overlap) continue; // new attempt
+ }
+ break;
+ }
+ if (is_overlap)
+ {
+ LOG_WARN("Can't placing sphere during " << maxAttempts << " attemps.");
+ position=backup_pos;
+ return;
+ }
+
+ rb->velocity = Vector3r(//
+ velocity[0]+velocityRange[0]*randomSymmetricUnit(),
+ velocity[1]+velocityRange[1]*randomSymmetricUnit(),
+ velocity[2]+velocityRange[2]*randomSymmetricUnit());
+ rb->angularVelocity= Vector3r(//
+ angularVelocity[0]+angularVelocityRange[0]*randomSymmetricUnit(),
+ angularVelocity[1]+angularVelocityRange[1]*randomSymmetricUnit(),
+ angularVelocity[2]+angularVelocityRange[2]*randomSymmetricUnit());
+
+ shiftedBodies.push_back(b);
+ }
+ } // next sphere
+}
+
+Vector3r ResetRandomPosition::generatePositionOnSurface()
+{
+ body_id_t facetId = factoryFacets[(*randomFacet)()];
+ Real t1 = randomUnit();
+ Real t2 = randomUnit()*(1-t1);
+
+ shared_ptr<Body> facet = Body::byId(facetId);
+ InteractingFacet* ifacet = static_cast<InteractingFacet*>(facet->interactingGeometry.get());
+
+ return t1*(ifacet->vertices[1]-ifacet->vertices[0])+t2*(ifacet->vertices[2]-ifacet->vertices[0])+ifacet->vertices[0]+facet->physicalParameters->se3.position;
+
+}
+
+Vector3r ResetRandomPosition::generatePositionInVolume()
+{
+ Vector3r p1 = generatePositionOnSurface();
+ Vector3r p2 = generatePositionOnSurface();
+ Real t = randomUnit();
+ return p1+t*(p2-p1);
+
+}
+
=== added file 'pkg/common/Engine/StandAloneEngine/ResetRandomPosition.hpp'
--- pkg/common/Engine/StandAloneEngine/ResetRandomPosition.hpp 1970-01-01 00:00:00 +0000
+++ pkg/common/Engine/StandAloneEngine/ResetRandomPosition.hpp 2009-11-19 06:44:50 +0000
@@ -0,0 +1,96 @@
+/*************************************************************************
+* Copyright (C) 2009 by Sergei Dorofeenko *
+* sega@xxxxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+#pragma once
+
+#include <yade/pkg-common/PeriodicEngines.hpp>
+#include <yade/pkg-common/InteractionGeometryMetaEngine.hpp>
+#include <yade/core/Collider.hpp>
+#include <yade/core/MetaBody.hpp>
+#include <vector>
+#include <string>
+
+using namespace std;
+
+/// @brief Produces spheres over the course of a simulation.
+class ResetRandomPosition : public PeriodicEngine {
+public:
+
+ ResetRandomPosition();
+ virtual ~ResetRandomPosition();
+
+ /// @brief Create one sphere per call.
+ virtual void action(MetaBody*);
+
+ /// @brief The geometry of the section on which spheres will be placed.
+ vector<body_id_t> factoryFacets;
+
+ /// @brief Factory section may be a surface or volume (convex).
+ /// By default it is a surface. To make its a volume set volumeSection=true
+ bool volumeSection;
+
+ /// @brief Max attemps to place sphere.
+ /// If placing the sphere in certain random position would cause an overlap with any other physical body in the model, SpheresFactory will try to find another position. Default 20 attempts allow.
+ int maxAttempts;
+
+ /// @brief Mean velocity of spheres.
+ Vector3r velocity;
+
+ /// @brief Half size of a velocities distribution interval.
+ /// New sphere will have random velocity within the range velocity±velocityRange.
+ Vector3r velocityRange;
+
+ /// @brief Mean angularVelocity of spheres.
+ Vector3r angularVelocity;
+
+ /// @brief Half size of a angularVelocity distribution interval.
+ /// New sphere will have random angularVelocity within the range angularVelocity±angularVelocityRange.
+ Vector3r angularVelocityRange;
+
+ Vector3r point;
+ Vector3r normal;
+
+ std::vector< int > subscribedBodies;
+
+private:
+ /// @brief Pointer to Collider.
+ /// It is necessary in order to probe the bounding volume for new sphere.
+ Collider* bI;
+
+ /// @brief Pointer to InteractionGeometryMetaEngine.
+ /// It is necessary in order to detect a real overlap with other bodies.
+ InteractionGeometryMetaEngine* iGME;
+
+ std::vector<shared_ptr<Body> > shiftedBodies;
+
+ bool first_run;
+
+ //bool generateNewPosition(const shared_ptr<Body>& b, Vector3r& new_position);
+ Vector3r generatePositionOnSurface();
+ Vector3r generatePositionInVolume();
+
+ typedef boost::variate_generator<boost::minstd_rand,boost::uniform_int<> > RandomInt;
+ shared_ptr<RandomInt> randomFacet;
+
+ DECLARE_LOGGER;
+
+ REGISTER_ATTRIBUTES(PeriodicEngine,
+ (factoryFacets)
+ (subscribedBodies)
+ (point)
+ (normal)
+ (volumeSection)
+ (maxAttempts)
+ (velocity)
+ (velocityRange)
+ (angularVelocity)
+ (angularVelocityRange))
+ REGISTER_CLASS_AND_BASE(ResetRandomPosition, StandAloneEngine);
+};
+REGISTER_SERIALIZABLE(ResetRandomPosition);
+
+
=== modified file 'py/utils.py'
--- py/utils.py 2009-11-17 12:26:35 +0000
+++ py/utils.py 2009-11-19 06:44:50 +0000
@@ -124,7 +124,7 @@
b.dynamic=False
return b
-def facet(vertices,dynamic=False,wire=True,color=None,density=1,highlight=False,noInteractingGeometry=False,physParamsClass='BodyMacroParameters',**physParamsAttr):
+def facet(vertices,dynamic=False,wire=True,color=None,density=1,highlight=False,noBoundingVolume=False,physParamsClass='BodyMacroParameters',**physParamsAttr):
"""Create default facet with given parameters. Vertices are given as sequence of 3 3-tuple and they, all in global coordinates."""
b=Body()
if not color: color=randomColor()
@@ -134,12 +134,11 @@
pp.update({'se3':[center[0],center[1],center[2],1,0,0,0],'refSe3':[center[0],center[1],center[2],1,0,0,0],'inertia':[0,0,0]})
b.phys=PhysicalParameters(physParamsClass)
b.phys.updateExistingAttrs(pp)
- b.bound=BoundingVolume('AABB',diffuseColor=[0,1,0])
+ if not noBoundingVolume: b.bound=BoundingVolume('AABB',diffuseColor=[0,1,0])
b['isDynamic']=dynamic
- if not noInteractingGeometry:
- b.mold=InteractingGeometry('InteractingFacet',diffuseColor=color,wire=wire,highlight=highlight)
- b.mold['vertices']=vertices
- b.mold.postProcessAttributes(True)
+ b.mold=InteractingGeometry('InteractingFacet',diffuseColor=color,wire=wire,highlight=highlight)
+ b.mold['vertices']=vertices
+ b.mold.postProcessAttributes(True)
return b
def facetBox(center,extents,orientation=[1,0,0,0],wallMask=63,**kw):
=== added file 'scripts/test/ResetRandomPosition.py'
--- scripts/test/ResetRandomPosition.py 1970-01-01 00:00:00 +0000
+++ scripts/test/ResetRandomPosition.py 2009-11-19 06:44:50 +0000
@@ -0,0 +1,62 @@
+# -*- coding: utf-8
+
+from yade import utils,pack,export,qt
+import gts,os
+
+def Plane(v1,v2,v3,v4):
+ pts = [ [Vector3(v1),Vector3(v2),Vector3(v3),Vector3(v4)] ]
+ return pack.sweptPolylines2gtsSurface(pts,capStart=True,capEnd=True)
+
+# Parameters
+tc=0.001# collision time
+en=0.3 # normal restitution coefficient
+es=0.3 # tangential restitution coefficient
+frictionAngle=radians(35)#
+density=2700
+kw=utils.getViscoelasticFromSpheresInteraction(10e3,tc,en,es)
+
+O.dt=.2*tc # time step
+
+Rs=0.02 # mean particle radius
+Rf=0.01 # dispersion (Rs±Rf*Rs)
+nSpheres=1000# number of particles
+
+# Create geometry
+pln=Plane( (-.5, -.5, 0), (.5, -.5, -.05), (.5, .5, 0), (-.5, .5, -.05) );
+plnIds=O.bodies.append(pack.gtsSurface2Facets(pln.faces(),frictionAngle=frictionAngle,physParamsClass='SimpleViscoelasticBodyParameters',color=(0,1,0),**kw))
+
+fct=Plane( (-.25, -.25, .5), (.25, -.25, .5), (.25, .25, .5), (-.25, .25, .5) );
+fctIds=O.bodies.append(pack.gtsSurface2Facets(fct.faces(),frictionAngle=frictionAngle,physParamsClass='SimpleViscoelasticBodyParameters',color=(1,0,0),noBoundingVolume=True))
+
+# Create spheres
+sp=pack.SpherePack();
+sp.makeCloud(Vector3(-.5, -.5, 0),Vector3(.5, .5, .2), Rs, Rf, int(nSpheres), False)
+spheres=O.bodies.append([utils.sphere(s[0],s[1],color=(0.929,0.412,0.412),density=density,frictionAngle=frictionAngle,physParamsClass="SimpleViscoelasticBodyParameters") for s in sp])
+for id in spheres:
+ s=O.bodies[id]
+ p=utils.getViscoelasticFromSpheresInteraction(s.phys['mass'],tc,en,es)
+ s.phys['kn'],s.phys['cn'],s.phys['ks'],s.phys['cs']=p['kn'],p['cn'],p['ks'],p['cs']
+
+# Create engines
+O.initializers=[ BoundingVolumeMetaEngine([InteractingSphere2AABB(),InteractingFacet2AABB(),MetaInteractingGeometry2AABB()]) ]
+O.engines=[
+ BexResetter(),
+ BoundingVolumeMetaEngine([InteractingSphere2AABB(),InteractingFacet2AABB(), MetaInteractingGeometry2AABB()
+ ]),
+ InsertionSortCollider(nBins=5,sweepLength=.1*Rs),
+ InteractionDispatchers(
+ [InteractingSphere2InteractingSphere4SpheresContactGeometry(), InteractingFacet2InteractingSphere4SpheresContactGeometry()],
+ [SimpleViscoelasticRelationships()],
+ [ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw()],
+ ),
+ GravityEngine(gravity=[0,0,-9.81]),
+ NewtonsDampedLaw(damping=0),
+ ResetRandomPosition(factoryFacets=fctIds,velocity=(0,0,-2),virtPeriod=0.01,subscribedBodies=spheres,point=(0,0,-.5),normal=(0,0,1),maxAttempts=100),
+]
+
+renderer = qt.Renderer()
+qt.View()
+O.saveTmp()
+O.run()
+
+
