yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10470
[Branch ~yade-pkg/yade/git-trunk] Rev 3822: Merge pull request #40 from burak-er/master
Merge authors:
burak er (burak-er)
------------------------------------------------------------
revno: 3822 [merge]
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Mon 2014-02-17 17:37:03 +0100
message:
Merge pull request #40 from burak-er/master
Added an adaptive Runge Kutta integrator scheme.
added:
examples/adaptiveintegrator/
examples/adaptiveintegrator/simple-scene-plot-NewtonIntegrator.py
examples/adaptiveintegrator/simple-scene-plot-RungeKuttaCashKarp54.py
pkg/dem/GeneralIntegratorInsertionSortCollider.cpp
pkg/dem/GeneralIntegratorInsertionSortCollider.hpp
pkg/dem/Integrator.cpp
pkg/dem/Integrator.hpp
pkg/dem/RungeKuttaCashKarp54Integrator.cpp
pkg/dem/RungeKuttaCashKarp54Integrator.hpp
modified:
pkg/common/InsertionSortCollider.hpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== added directory 'examples/adaptiveintegrator'
=== added file 'examples/adaptiveintegrator/simple-scene-plot-NewtonIntegrator.py'
--- examples/adaptiveintegrator/simple-scene-plot-NewtonIntegrator.py 1970-01-01 00:00:00 +0000
+++ examples/adaptiveintegrator/simple-scene-plot-NewtonIntegrator.py 2014-02-17 14:31:35 +0000
@@ -0,0 +1,49 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+import matplotlib
+matplotlib.use('TkAgg')
+
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
+ [Ip2_FrictMat_FrictMat_FrictPhys()],
+ [Law2_ScGeom_FrictPhys_CundallStrack()]
+ ),
+ NewtonIntegrator(damping=0.0,gravity=(0,0,-9.81)),
+ ###
+ ### NOTE this extra engine:
+ ###
+ ### You want snapshot to be taken every 1 sec (realTimeLim) or every 50 iterations (iterLim),
+ ### whichever comes soones. virtTimeLim attribute is unset, hence virtual time period is not taken into account.
+ PyRunner(iterPeriod=20,command='myAddPlotData()')
+]
+O.bodies.append(box(center=[0,0,0],extents=[.5,.5,.5],fixed=True,color=[1,0,0]))
+O.bodies.append(sphere([0,0,2],1,color=[0,1,0]))
+O.dt=.002*PWaveTimeStep()
+
+
+############################################
+##### now the part pertaining to plots #####
+############################################
+
+from yade import plot
+## we will have 2 plots:
+## 1. t as function of i (joke test function)
+## 2. i as function of t on left y-axis ('|||' makes the separation) and z_sph, v_sph (as green circles connected with line) and z_sph_half again as function of t
+plot.plots={'i':('t'),'t':('z_sph',None,('v_sph','go-'),'z_sph_half')}
+
+## this function is called by plotDataCollector
+## it should add data with the labels that we will plot
+## if a datum is not specified (but exists), it will be NaN and will not be plotted
+def myAddPlotData():
+ sph=O.bodies[1]
+ ## store some numbers under some labels
+ plot.addData(t=O.time,i=O.iter,z_sph=sph.state.pos[2],z_sph_half=.5*sph.state.pos[2],v_sph=sph.state.vel.norm())
+print "Now calling plot.plot() to show the figures. The timestep is artificially low so that you can watch graphs being updated live."
+plot.liveInterval=.2
+plot.plot(subPlots=False)
+O.run(int(5./O.dt));
+#plot.saveGnuplot('/tmp/a')
+## you can also access the data in plot.data['i'], plot.data['t'] etc, under the labels they were saved.
=== added file 'examples/adaptiveintegrator/simple-scene-plot-RungeKuttaCashKarp54.py'
--- examples/adaptiveintegrator/simple-scene-plot-RungeKuttaCashKarp54.py 1970-01-01 00:00:00 +0000
+++ examples/adaptiveintegrator/simple-scene-plot-RungeKuttaCashKarp54.py 2014-02-17 14:31:35 +0000
@@ -0,0 +1,65 @@
+#!/usr/bin/python
+# Burak ER
+# burak.er@xxxxxxxxxx
+# github.com/burak-er
+# Mechanical Engineering Department
+# Bursa Technical University
+#
+# -*- coding: utf-8 -*-
+import matplotlib
+matplotlib.use('TkAgg')
+
+# Use an integrator engine that is derived from the interface Integrator.
+
+#RungeKuttaCashKarp54Integrator integrator performs one step of simulation for the given tolerances. Whether the time step is given, it completes it then stops.
+
+integrator=RungeKuttaCashKarp54Integrator([
+ ForceResetter(),
+ GeneralIntegratorInsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
+ [Ip2_FrictMat_FrictMat_FrictPhys()],
+ [Law2_ScGeom_FrictPhys_CundallStrack()]
+ ),
+ GravityEngine(gravity=Vector3(0,0,-9.81)),
+ PyRunner(virtPeriod=1e-99,command='myAddPlotData()')#use virtPeriod on this integrator.
+
+]);
+
+#Tolerances can be set for the optimum accuracy
+integrator.rel_err=1e-6;
+integrator.abs_err=1e-6;
+
+O.engines=[integrator,
+]
+
+O.bodies.append(box(center=[0,0,0],extents=[.5,.5,.5],fixed=True,color=[1,0,0]))
+O.bodies.append(sphere([0,0,2],1,color=[0,1,0]))
+O.dt=1e-2# this signifies the endpoint. It is not much important for the accuracy of the integration where accuracy is defined by rel_err and abs_err of the integrator.
+
+
+
+############################################
+##### now the part pertaining to plots #####
+############################################
+
+from yade import plot
+## we will have 2 plots:
+## 1. t as function of i (joke test function)
+## 2. i as function of t on left y-axis ('|||' makes the separation) and z_sph, v_sph (as green circles connected with line) and z_sph_half again as function of t
+plot.plots={'i':('t'),'t':('z_sph',None,('v_sph','go-'),'z_sph_half')}
+
+## this function is called by plotDataCollector
+## it should add data with the labels that we will plot
+## if a datum is not specified (but exists), it will be NaN and will not be plotted
+def myAddPlotData():
+ sph=O.bodies[1]
+ ## store some numbers under some labels
+ plot.addData(t=O.time,i=O.iter,z_sph=sph.state.pos[2],z_sph_half=.5*sph.state.pos[2],v_sph=sph.state.vel.norm())
+print "Now calling plot.plot() to show the figures. The timestep is artificially low so that you can watch graphs being updated live."
+plot.liveInterval=.2
+plot.plot(subPlots=False)
+print "Number of threads ", os.environ['OMP_NUM_THREADS']
+O.run(int(5./O.dt));
+#plot.saveGnuplot('/tmp/a')
+## you can also access the data in plot.data['i'], plot.data['t'] etc, under the labels they were saved.
=== modified file 'pkg/common/InsertionSortCollider.hpp'
--- pkg/common/InsertionSortCollider.hpp 2012-09-08 01:19:45 +0000
+++ pkg/common/InsertionSortCollider.hpp 2014-02-17 14:31:35 +0000
@@ -82,7 +82,14 @@
class NewtonIntegrator;
+class Integrator;
+
+class GeneralIntegratorInsertionSortCollider;// Forward decleration of child to decleare it as friend
+
class InsertionSortCollider: public Collider{
+
+ friend GeneralIntegratorInsertionSortCollider;
+
//! struct for storing bounds of bodies
struct Bounds{
//! coordinate along the given sortAxis
=== added file 'pkg/dem/GeneralIntegratorInsertionSortCollider.cpp'
--- pkg/dem/GeneralIntegratorInsertionSortCollider.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/GeneralIntegratorInsertionSortCollider.cpp 2014-02-17 14:31:35 +0000
@@ -0,0 +1,228 @@
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+
+#include"GeneralIntegratorInsertionSortCollider.hpp"
+#include<yade/core/Scene.hpp>
+#include<yade/core/Interaction.hpp>
+#include<yade/core/InteractionContainer.hpp>
+#include<yade/pkg/common/Dispatching.hpp>
+#include<yade/pkg/dem/Integrator.hpp>
+#include<yade/pkg/common/Sphere.hpp>
+
+#include<algorithm>
+#include<vector>
+#include<boost/static_assert.hpp>
+
+using namespace std;
+
+YADE_PLUGIN((GeneralIntegratorInsertionSortCollider))
+CREATE_LOGGER(GeneralIntegratorInsertionSortCollider);
+
+// STRIDE
+bool GeneralIntegratorInsertionSortCollider::isActivated(){
+ // activated if number of bodies changes (hence need to refresh collision information)
+ // or the time of scheduled run already came, or we were never scheduled yet
+ if(!strideActive) return true;
+ if(!integrator) return true;
+ if(fastestBodyMaxDist<0){fastestBodyMaxDist=0; return true;}
+ fastestBodyMaxDist=integrator->maxVelocitySq;
+ if(fastestBodyMaxDist>=1 || fastestBodyMaxDist==0) return true;
+ if((size_t)BB[0].size!=2*scene->bodies->size()) return true;
+ if(scene->interactions->dirty) return true;
+ if(scene->doSort) { scene->doSort=false; return true; }
+ return false;
+ }
+
+void GeneralIntegratorInsertionSortCollider::action(){
+ #ifdef ISC_TIMING
+ timingDeltas->start();
+ #endif
+
+ long nBodies=(long)scene->bodies->size();
+ InteractionContainer* interactions=scene->interactions.get();
+ scene->interactions->iterColliderLastRun=-1;
+
+ // periodicity changed, force reinit
+ if(scene->isPeriodic != periodic){
+ for(int i=0; i<3; i++) BB[i].vec.clear();
+ periodic=scene->isPeriodic;
+ }
+ // pre-conditions
+ // adjust storage size
+ bool doInitSort=false;
+ if (doSort) {
+ doInitSort=true;
+ doSort=false;
+ }
+ if(BB[0].size!=2*nBodies){
+ long BBsize=BB[0].size;
+ LOG_DEBUG("Resize bounds containers from "<<BBsize<<" to "<<nBodies*2<<", will std::sort.");
+ // bodies deleted; clear the container completely, and do as if all bodies were added (rather slow…)
+ // future possibility: insertion sort with such operator that deleted bodies would all go to the end, then just trim bounds
+ if(2*nBodies<BBsize){ for(int i=0; i<3; i++) BB[i].vec.clear(); }
+ // more than 100 bodies was added, do initial sort again
+ // maybe: should rather depend on ratio of added bodies to those already present...?
+ if(2*nBodies-BBsize>200 || BBsize==0) doInitSort=true;
+ assert((BBsize%2)==0);
+ for(int i=0; i<3; i++){
+ BB[i].vec.reserve(2*nBodies);
+ // add lower and upper bounds; coord is not important, will be updated from bb shortly
+ for(long id=BBsize/2; id<nBodies; id++){ BB[i].vec.push_back(Bounds(0,id,/*isMin=*/true)); BB[i].vec.push_back(Bounds(0,id,/*isMin=*/false)); }
+ BB[i].size=BB[i].vec.size();
+ }
+ }
+ if(minima.size()!=(size_t)3*nBodies){ minima.resize(3*nBodies); maxima.resize(3*nBodies); }
+ assert((size_t)BB[0].size==2*scene->bodies->size());
+
+ // update periodicity
+ assert(BB[0].axis==0); assert(BB[1].axis==1); assert(BB[2].axis==2);
+ if(periodic) for(int i=0; i<3; i++) BB[i].updatePeriodicity(scene);
+
+ // compatibility block, can be removed later
+ findBoundDispatcherInEnginesIfNoFunctorsAndWarn();
+
+ if(verletDist<0){
+ Real minR=std::numeric_limits<Real>::infinity();
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if(!b || !b->shape) continue;
+ Sphere* s=dynamic_cast<Sphere*>(b->shape.get());
+ if(!s) continue;
+ minR=min(s->radius,minR);
+ }
+ if (isinf(minR)) LOG_ERROR("verletDist is set to 0 because no spheres were found. It will result in suboptimal performances, consider setting a positive verletDist in your script.");
+ // if no spheres, disable stride
+ verletDist=isinf(minR) ? 0 : abs(verletDist)*minR;
+ }
+
+ // update bounds via boundDispatcher
+ boundDispatcher->scene=scene;
+ boundDispatcher->sweepDist=verletDist;
+ boundDispatcher->minSweepDistFactor=minSweepDistFactor;
+ boundDispatcher->targetInterv=targetInterv;
+ boundDispatcher->updatingDispFactor=updatingDispFactor;
+ boundDispatcher->action();
+
+ // if interactions are dirty, force reinitialization
+ if(scene->interactions->dirty){
+ doInitSort=true;
+ scene->interactions->dirty=false;
+ }
+
+ // STRIDE
+ if(verletDist>0){
+ // get the Integrator, to ask for the maximum velocity value
+ if(!integrator){
+ FOREACH(shared_ptr<Engine>& e, scene->engines){ integrator=dynamic_pointer_cast<Integrator>(e); if(integrator) break; }
+ if(!integrator){ throw runtime_error("InsertionSortCollider.verletDist>0, but unable to locate any Integrator within O.engines."); }
+ }
+ }
+ ISC_CHECKPOINT("init");
+
+ // STRIDE
+ // get us ready for strides, if they were deactivated
+ if(!strideActive && verletDist>0 && integrator->maxVelocitySq>=0){ // maxVelocitySq is a really computed value
+ strideActive=true;
+ }
+ if(strideActive){
+ assert(verletDist>0);
+ assert(strideActive); assert(integrator->maxVelocitySq>=0);
+ integrator->updatingDispFactor=updatingDispFactor;
+ } else { /* !strideActive */
+ boundDispatcher->sweepDist=0;
+ }
+
+ ISC_CHECKPOINT("bound");
+
+ // copy bounds along given axis into our arrays
+ for(long i=0; i<2*nBodies; i++){
+ for(int j=0; j<3; j++){
+ VecBounds& BBj=BB[j];
+ const Body::id_t id=BBj[i].id;
+ const shared_ptr<Body>& b=Body::byId(id,scene);
+ if(b){
+ const shared_ptr<Bound>& bv=b->bound;
+ // coordinate is min/max if has bounding volume, otherwise both are the position. Add periodic shift so that we are inside the cell
+ // watch out for the parentheses around ?: within ?: (there was unwanted conversion of the Reals to bools!)
+
+ BBj[i].coord=((BBj[i].flags.hasBB=((bool)bv)) ? (BBj[i].flags.isMin ? bv->min[j] : bv->max[j]) : (b->state->pos[j])) - (periodic ? BBj.cellDim*BBj[i].period : 0.);
+
+ } else { BBj[i].flags.hasBB=false; /* for vanished body, keep the coordinate as-is, to minimize inversions. */ }
+ // if initializing periodic, shift coords & record the period into BBj[i].period
+ if(doInitSort && periodic) {
+ BBj[i].coord=cellWrap(BBj[i].coord,0,BBj.cellDim,BBj[i].period);
+ }
+ }
+ }
+ // for each body, copy its minima and maxima, for quick checks of overlaps later
+ BOOST_STATIC_ASSERT(sizeof(Vector3r)==3*sizeof(Real));
+ for(Body::id_t id=0; id<nBodies; id++){
+ const shared_ptr<Body>& b=Body::byId(id,scene);
+ if(b){
+ const shared_ptr<Bound>& bv=b->bound;
+ if(bv) { memcpy(&minima[3*id],&bv->min,3*sizeof(Real)); memcpy(&maxima[3*id],&bv->max,3*sizeof(Real)); } // ⇐ faster than 6 assignments
+ else{ const Vector3r& pos=b->state->pos; memcpy(&minima[3*id],&pos,3*sizeof(Real)); memcpy(&maxima[3*id],&pos,3*sizeof(Real)); }
+ } else { memset(&minima[3*id],0,3*sizeof(Real)); memset(&maxima[3*id],0,3*sizeof(Real)); }
+ }
+
+ ISC_CHECKPOINT("copy");
+
+ // process interactions that the constitutive law asked to be erased
+// interactions->erasePending(*this,scene);
+ interactions->conditionalyEraseNonReal(*this,scene);
+
+ ISC_CHECKPOINT("erase");
+
+ // sort
+ // the regular case
+ if(!doInitSort && !sortThenCollide){
+ /* each inversion in insertionSort calls handleBoundInversion, which in turns may add/remove interaction */
+ if(!periodic) for(int i=0; i<3; i++) insertionSort(BB[i],interactions,scene);
+ else for(int i=0; i<3; i++) insertionSortPeri(BB[i],interactions,scene);
+ }
+ // create initial interactions (much slower)
+ else {
+ if(doInitSort){
+ // the initial sort is in independent in 3 dimensions, may be run in parallel; it seems that there is no time gain running in parallel, though
+ // important to reset loInx for periodic simulation (!!)
+ for(int i=0; i<3; i++) { BB[i].loIdx=0; std::sort(BB[i].vec.begin(),BB[i].vec.end()); }
+ numReinit++;
+ } else { // sortThenCollide
+ if(!periodic) for(int i=0; i<3; i++) insertionSort(BB[i],interactions,scene,false);
+ else for(int i=0; i<3; i++) insertionSortPeri(BB[i],interactions,scene,false);
+ }
+ // traverse the container along requested axis
+ assert(sortAxis==0 || sortAxis==1 || sortAxis==2);
+ VecBounds& V=BB[sortAxis];
+ // go through potential aabb collisions, create interactions as necessary
+ if(!periodic){
+ for(long i=0; i<2*nBodies; i++){
+ // start from the lower bound (i.e. skipping upper bounds)
+ // skip bodies without bbox, because they don't collide
+ if(!(V[i].flags.isMin && V[i].flags.hasBB)) continue;
+ const Body::id_t& iid=V[i].id;
+ // go up until we meet the upper bound
+ for(long j=i+1; /* handle case 2. of swapped min/max */ j<2*nBodies && V[j].id!=iid; j++){
+ const Body::id_t& jid=V[j].id;
+ // take 2 of the same condition (only handle collision [min_i..max_i]+min_j, not [min_i..max_i]+min_i (symmetric)
+ if(!V[j].flags.isMin) continue;
+ /* abuse the same function here; since it does spatial overlap check first, it is OK to use it */
+ handleBoundInversion(iid,jid,interactions,scene);
+ assert(j<2*nBodies-1);
+ }
+ }
+ } else { // periodic case: see comments above
+ for(long i=0; i<2*nBodies; i++){
+ if(!(V[i].flags.isMin && V[i].flags.hasBB)) continue;
+ const Body::id_t& iid=V[i].id;
+ long cnt=0;
+ // we might wrap over the periodic boundary here; that's why the condition is different from the aperiodic case
+ for(long j=V.norm(i+1); V[j].id!=iid; j=V.norm(j+1)){
+ const Body::id_t& jid=V[j].id;
+ if(!V[j].flags.isMin) continue;
+ handleBoundInversionPeri(iid,jid,interactions,scene);
+ if(cnt++>2*(long)nBodies){ LOG_FATAL("Uninterrupted loop in the initial sort?"); throw std::logic_error("loop??"); }
+ }
+ }
+ }
+ }
+ ISC_CHECKPOINT("sort&collide");
+}
=== added file 'pkg/dem/GeneralIntegratorInsertionSortCollider.hpp'
--- pkg/dem/GeneralIntegratorInsertionSortCollider.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/GeneralIntegratorInsertionSortCollider.hpp 2014-02-17 14:31:35 +0000
@@ -0,0 +1,39 @@
+// 2014 Burak ER <burak.er@xxxxxxxxxx>
+
+#pragma once
+#include<yade/pkg/common/InsertionSortCollider.hpp>
+
+
+/*!
+ Adaptive Integration Sort Collider:
+
+ Changing the Integrator dependence from Newton Integrator to Arbitrary Integrators. Arbitrary integrators should use the Integrator interface.
+
+*/
+
+
+
+#ifdef ISC_TIMING
+ #define ISC_CHECKPOINT(cpt) timingDeltas->checkpoint(cpt)
+#else
+ #define ISC_CHECKPOINT(cpt)
+#endif
+
+class Integrator;
+
+class GeneralIntegratorInsertionSortCollider: public InsertionSortCollider{
+
+ // we need this to find out about current maxVelocitySq
+ shared_ptr<Integrator> integrator;
+ // if False, no type of striding is used
+
+ public:
+
+ virtual bool isActivated(); //override this function to change NewtonIntegrator dependency.
+
+ virtual void action(); //override this function to change behaviour with the NewtonIntegrator dependency.
+
+ YADE_CLASS_BASE_DOC(GeneralIntegratorInsertionSortCollider,InsertionSortCollider," This class is the adaptive version of the InsertionSortCollider and changes the NewtonIntegrator dependency of the collider algorithms to the Integrator interface which is more general.");
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(GeneralIntegratorInsertionSortCollider);
=== added file 'pkg/dem/Integrator.cpp'
--- pkg/dem/Integrator.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Integrator.cpp 2014-02-17 14:31:35 +0000
@@ -0,0 +1,351 @@
+#include<yade/core/Clump.hpp>
+#include<yade/core/Scene.hpp>
+#include<yade/pkg/dem/Integrator.hpp>
+#include<boost/python.hpp>
+using namespace boost;
+#ifdef YADE_OPENMP
+ #include<omp.h>
+#endif
+
+YADE_PLUGIN((Integrator));
+
+
+
+void observer::operator()( const stateVector& x , Real t ) const
+{
+ this->integrator->scene->time=t;
+
+ this->integrator->setCurrentStates(x);
+
+}
+
+
+//! Integrator's pseudo-ctor (factory), taking nested lists of slave engines (might be moved to real ctor perhaps)
+
+
+void Integrator::action(){
+
+
+}
+
+void Integrator::system(const stateVector& currentstates, stateVector& derivatives, Real time)
+{
+
+ #ifdef YADE_OPENMP
+ //prevent https://bugs.launchpad.net/yade/+bug/923929
+ ensureSync();
+ #endif
+
+ //Calculate orientation
+
+ maxVelocitySq=-1;
+
+ setCurrentStates(currentstates);
+
+ scene->time=time;
+
+ const int size=(int)slaves.size();
+
+ for(int i=0; i<size; i++){
+ // run every slave group sequentially
+ FOREACH(const shared_ptr<Engine>& e, slaves[i]) {
+ e->scene=scene;
+ if(!e->dead && e->isActivated()) e->action();
+ }
+ }
+ derivatives=getSceneStateDot();
+
+/*
+ std::cout<<std::endl<<"Derivatives are"<<std::endl;
+ for(long int k=0;k<derivatives.size();k++)
+ std::cout<<std::endl<<derivatives[k]<<std::endl;*/
+}
+
+stateVector& Integrator::getSceneStateDot(void){
+
+ try{
+
+ long int numberofscenebodies=scene->bodies->size();
+
+ scene->forces.sync();
+
+ accumstatedotofthescene.resize(2*scene->bodies->size()*7);
+
+ YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
+
+ const Body::id_t& id=b->getId();
+
+ Vector3r force=Vector3r::Zero();
+
+ Vector3r vel_current;
+
+ Vector3r moment=Vector3r::Zero();
+
+ Vector3r angvel_current;
+
+ Quaternionr ori_current;
+
+ Quaternionr angvelquat;
+
+ Quaternionr oridot_current;
+
+ if(!b->isClumpMember()) {
+
+ Real mass=b->state->mass;
+
+ Vector3r inertia=b->state->inertia;
+
+ vel_current=b->state->vel;
+
+ angvel_current=b->state->angVel;
+
+ ori_current=b->state->ori;
+
+ // clumps forces
+ if(b->isClump()) {
+ b->shape->cast<Clump>().addForceTorqueFromMembers(b->state.get(),scene,force,moment);
+ #ifdef YADE_OPENMP
+ //it is safe here, since only one thread will read/write
+ scene->forces.getTorqueUnsynced(id)+=moment;
+ scene->forces.getForceUnsynced(id)+=force;
+ #else
+ scene->forces.addTorque(id,moment);
+ scene->forces.addForce(id,force);
+ #endif
+ }
+
+
+ force=scene->forces.getForce(id); moment=scene->forces.getTorque(id);
+
+ /*
+ * Calculation of accelerations
+ *
+ */
+
+
+ force[0]=force[0]/mass; force[1]= force[1]/mass; force[2]= force[2]/mass; //Calculate linear acceleration
+
+ moment[0]=moment[0]/inertia[0]; moment[1]= moment[1]/inertia[1]; moment[2]= moment[2]/inertia[2]; //Calculate angular acceleration
+
+ //Check for fixation
+ /*
+ This code block needs optimization
+ */
+ string str="xyzXYZ"; // Very very very hard coding!!!!! Fixation seems not handled fine by state structure, should be improved.
+
+ for(int i=0; i<3; i++) if(b->state->blockedDOFs_vec_get().find(str[i]) != std::string::npos){ force[i]=0;vel_current[i]=0;}
+ for(int i=3; i<6; i++) if(b->state->blockedDOFs_vec_get().find(str[i]) != std::string::npos){ moment[i-3]=0;angvel_current[i-3]=0;}
+
+ angvelquat = Quaternionr(0.0,angvel_current[0],angvel_current[1],angvel_current[2]);
+
+ oridot_current=0.5*angvelquat*ori_current;
+
+
+ // if (densityScaling) accel*=state->densityScaling;
+
+ }
+ else
+ {
+
+ //is clump member
+ force=Vector3r::Zero();
+ moment=Vector3r::Zero();
+ vel_current=Vector3r::Zero();
+ angvel_current=Vector3r::Zero();
+ oridot_current=Quaternionr(0,0,0,0);// zero change in quaternion with respect to time
+ }
+
+
+ /*Orientation differantion is straight forward.*/
+ accumstatedotofthescene[id*7+0]=vel_current[0]; accumstatedotofthescene[(id+numberofscenebodies)*7+0]=force[0];
+ accumstatedotofthescene[id*7+1]=vel_current[1]; accumstatedotofthescene[(id+numberofscenebodies)*7+1]=force[1];
+ accumstatedotofthescene[id*7+2]=vel_current[2]; accumstatedotofthescene[(id+numberofscenebodies)*7+2]=force[2];
+ accumstatedotofthescene[id*7+3]=oridot_current.w(); accumstatedotofthescene[(id+numberofscenebodies)*7+3]=moment[0];
+ accumstatedotofthescene[id*7+4]=oridot_current.x(); accumstatedotofthescene[(id+numberofscenebodies)*7+4]=moment[1];
+ accumstatedotofthescene[id*7+5]=oridot_current.y(); accumstatedotofthescene[(id+numberofscenebodies)*7+5]=moment[2];
+ accumstatedotofthescene[id*7+6]=oridot_current.z(); accumstatedotofthescene[(id+numberofscenebodies)*7+6]=0;
+
+
+ } YADE_PARALLEL_FOREACH_BODY_END();
+
+
+ }
+ catch(std::exception& e){
+
+ LOG_FATAL("Unhandled exception at Integrator::getSceneStateDot the exception information : "<<typeid(e).name()<<" : "<<e.what());
+ }
+
+ return accumstatedotofthescene;
+
+
+
+}
+
+
+
+stateVector& Integrator::getCurrentStates(void)
+{
+
+
+ try{
+
+ long int numberofscenebodies=scene->bodies->size();
+
+ accumstateofthescene.resize(2*scene->bodies->size()*7);
+
+ YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
+
+ const Body::id_t& id=b->getId();
+
+ Vector3r pos_current=b->state->pos;
+
+ Vector3r vel_current=b->state->vel;
+
+ Quaternionr ori=b->state->ori;
+
+ Vector3r angvel=b->state->angVel;
+
+ accumstateofthescene[id*7+0]=pos_current[0]; accumstateofthescene[(id+numberofscenebodies)*7+0]=vel_current[0];
+ accumstateofthescene[id*7+1]=pos_current[1]; accumstateofthescene[(id+numberofscenebodies)*7+1]=vel_current[1];
+ accumstateofthescene[id*7+2]=pos_current[2]; accumstateofthescene[(id+numberofscenebodies)*7+2]=vel_current[2];
+ accumstateofthescene[id*7+3]=ori.w(); accumstateofthescene[(id+numberofscenebodies)*7+3]=angvel[0];
+ accumstateofthescene[id*7+4]=ori.x(); accumstateofthescene[(id+numberofscenebodies)*7+4]=angvel[1];
+ accumstateofthescene[id*7+5]=ori.y(); accumstateofthescene[(id+numberofscenebodies)*7+5]=angvel[2];
+ accumstateofthescene[id*7+6]=ori.z(); accumstateofthescene[(id+numberofscenebodies)*7+6]=0;
+
+
+ } YADE_PARALLEL_FOREACH_BODY_END();
+
+
+ }
+ catch(std::exception& e){
+
+ LOG_FATAL("Unhandled exception at Integrator::getCurrentStates the exception information : "<<typeid(e).name()<<" : "<<e.what());
+ }
+
+ return accumstateofthescene;
+
+}
+
+bool Integrator::setCurrentStates(stateVector yscene)
+{
+
+ try{
+
+ long int numberofscenebodies=scene->bodies->size();
+
+ //Zero max velocity for each thread
+ #ifdef YADE_OPENMP
+ FOREACH(Real& thrMaxVSq, threadMaxVelocitySq) { thrMaxVSq=0; }
+ #endif
+
+ YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
+
+ if(b->isClumpMember()) continue;
+
+ const Body::id_t& id=b->getId();
+
+ Vector3r pos_current;
+
+ pos_current<<yscene[id*7+0],yscene[id*7+1],yscene[id*7+2];
+
+ Vector3r vel_current;
+
+ vel_current<<yscene[(id+numberofscenebodies)*7+0],yscene[(id+numberofscenebodies)*7+1],yscene[(id+numberofscenebodies)*7+2];
+
+ Quaternionr ori=Quaternionr(yscene[id*7+3],yscene[id*7+4],yscene[id*7+5],yscene[id*7+6]);
+
+ Vector3r angvel;
+
+ angvel<<yscene[(id+numberofscenebodies)*7+3],yscene[(id+numberofscenebodies)*7+4],yscene[(id+numberofscenebodies)*7+5];
+
+ b->state->pos=pos_current;
+
+ b->state->vel=vel_current;
+
+ b->state->ori=ori;
+
+ //std::cout<<"Setting orientation to "<<ori<<std::endl;
+
+ b->state->ori.normalize(); //Normalize orientation
+
+ //std::cout<<"Setting angvel to "<<angvel<<std::endl;
+
+ b->state->angVel=angvel;
+
+ #ifdef YADE_OPENMP
+ Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,b->state->vel.squaredNorm());
+ #else
+ maxVelocitySq=max(maxVelocitySq,b->state->vel.squaredNorm());// Set maximum velocity of the scene
+ #endif
+
+ if(b->isClump()) Clump::moveMembers(b,scene,this);
+
+ } YADE_PARALLEL_FOREACH_BODY_END();
+
+ #ifdef YADE_OPENMP
+ FOREACH(const Real& thrMaxVSq, threadMaxVelocitySq) { maxVelocitySq=max(maxVelocitySq,thrMaxVSq); }
+ #endif
+
+
+ }
+ catch(std::exception& e){
+
+ LOG_FATAL("Unhandled exception at Integrator::setCurrentStates the exception information : "<<typeid(e).name()<<" : "<<e.what());
+
+ return false;
+ }
+
+ return true;
+}
+
+
+#ifdef YADE_OPENMP
+void Integrator::ensureSync()
+{
+ if (syncEnsured) return;
+ YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
+// if(b->isClump()) continue;
+ scene->forces.addForce(b->getId(),Vector3r(0,0,0));
+ } YADE_PARALLEL_FOREACH_BODY_END();
+ syncEnsured=true;
+}
+#endif
+
+
+void Integrator::saveMaximaDisplacement(const shared_ptr<Body>& b){
+ if (!b->bound) return;//clumps for instance, have no bounds, hence not saved
+ Vector3r disp=b->state->pos-b->bound->refPos;
+ Real maxDisp=max(abs(disp[0]),max(abs(disp[1]),abs(disp[2])));
+ if (!maxDisp || maxDisp<b->bound->sweepLength) {/*b->bound->isBounding = (updatingDispFactor>0 && (updatingDispFactor*maxDisp)<b->bound->sweepLength);*/
+ maxDisp=0.5;//not 0, else it will be seen as "not updated" by the collider, but less than 1 means no colliding
+ }
+ else {/*b->bound->isBounding = false;*/ maxDisp=2;/*2 is more than 1, enough to trigger collider*/}
+
+ maxVelocitySq=max(maxVelocitySq,maxDisp);
+}
+
+void Integrator::slaves_set(const python::list& slaves2){
+std::cout<<"Adding slaves";
+ int len=python::len(slaves2);
+ slaves.clear();
+ for(int i=0; i<len; i++){
+ python::extract<std::vector<shared_ptr<Engine> > > serialGroup(slaves2[i]);
+ if (serialGroup.check()){ slaves.push_back(serialGroup()); continue; }
+ python::extract<shared_ptr<Engine> > serialAlone(slaves2[i]);
+ if (serialAlone.check()){ vector<shared_ptr<Engine> > aloneWrap; aloneWrap.push_back(serialAlone()); slaves.push_back(aloneWrap); continue; }
+ PyErr_SetString(PyExc_TypeError,"Engines that are given to Integrator should be in two cases (a) in an ordered group, (b) alone engines");
+ python::throw_error_already_set();
+ }
+}
+
+python::list Integrator::slaves_get(){
+ python::list ret;
+ FOREACH(vector<shared_ptr<Engine > >& grp, slaves){
+ if(grp.size()==1) ret.append(python::object(grp[0]));
+ else ret.append(python::object(grp));
+ }
+ return ret;
+}
+
+
=== added file 'pkg/dem/Integrator.hpp'
--- pkg/dem/Integrator.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Integrator.hpp 2014-02-17 14:31:35 +0000
@@ -0,0 +1,113 @@
+#pragma once
+
+#include<yade/core/TimeStepper.hpp>
+#include<boost/python.hpp>
+
+class Integrator;
+
+
+typedef std::vector<Real> stateVector;// Currently, we are unable to use Eigen library within odeint
+
+/*Observer used to update the state of the scene*/
+class observer
+{
+ Integrator* integrator;
+public:
+ observer(Integrator* _in):integrator(_in){}
+ void operator()( const stateVector& /* x */ , Real /* t */ ) const;
+};
+
+//[ ode_wrapper
+template< class Obj , class Mem >
+class ode_wrapper
+{
+ Obj m_obj;
+ Mem m_mem;
+
+public:
+
+ ode_wrapper( Obj obj , Mem mem ) : m_obj( obj ) , m_mem( mem ) { }
+
+ template< class State , class Deriv , class Time >
+ void operator()( const State &x , Deriv &dxdt , Time t )
+ {
+ (m_obj.*m_mem)( x , dxdt , t );
+ }
+};
+
+template< class Obj , class Mem >
+ode_wrapper< Obj , Mem > make_ode_wrapper( Obj obj , Mem mem )
+{
+ return ode_wrapper< Obj , Mem >( obj , mem );
+}
+//]
+
+
+
+class Integrator: public TimeStepper {
+
+ public:
+
+ stateVector accumstateofthescene;//pos+vel
+
+ stateVector accumstatedotofthescene;//only the accelerations
+
+ stateVector resetstate;//last state before integration attempt
+
+ Real timeresetvalue;
+
+ inline void evaluateQuaternions(const stateVector &); //evaluate quaternions after integration
+
+ typedef vector<vector<shared_ptr<Engine> > > slaveContainer;
+
+ #ifdef YADE_OPENMP
+ vector<Real> threadMaxVelocitySq;
+ #endif
+
+ virtual void action();
+
+ virtual void system(const stateVector&, stateVector&, Real); //System function to calculate the derivatives of states
+
+ virtual bool isActivated(){return true;}
+ // py access
+ boost::python::list slaves_get();
+
+ stateVector& getSceneStateDot();
+
+ bool saveCurrentState(Scene const* ourscene);//Before any integration attempt state of the scene should be saved.
+
+ bool resetLastState(void);//Before any integration attempt state of the scene should be saved.
+
+ void slaves_set(const boost::python::list& slaves);
+
+ stateVector& getCurrentStates(void);
+
+ bool setCurrentStates(stateVector);
+
+ Real updatingDispFactor;//(experimental) Displacement factor used to trigger bound update: the bound is updated only if updatingDispFactor*disp>sweepDist when >0, else all bounds are updated.
+
+ void saveMaximaDisplacement(const shared_ptr<Body>& b);
+
+ #ifdef YADE_OPENMP
+ void ensureSync(); bool syncEnsured;
+ #endif
+
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Integrator,TimeStepper,"Integration Engine Interface.",
+ ((slaveContainer,slaves,,,"[will be overridden]"))
+ ((Real,integrationsteps,,,"all integrationsteps count as all succesfull substeps"))
+ ((Real,maxVelocitySq,NaN,,"store square of max. velocity, for informative purposes; computed again at every step. |yupdate|"))
+ ,
+ /*ctor*/
+ #ifdef YADE_OPENMP
+ threadMaxVelocitySq.resize(omp_get_max_threads()); syncEnsured=false;
+ #endif
+ ,
+ /*py*/
+
+ .add_property("slaves",&Integrator::slaves_get,&Integrator::slaves_set,"List of lists of Engines to calculate the force acting on the particles; to obtain the derivatives of the states, engines inside will be run sequentially.");
+ );
+};
+REGISTER_SERIALIZABLE(Integrator);
+
+
=== added file 'pkg/dem/RungeKuttaCashKarp54Integrator.cpp'
--- pkg/dem/RungeKuttaCashKarp54Integrator.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/RungeKuttaCashKarp54Integrator.cpp 2014-02-17 14:31:35 +0000
@@ -0,0 +1,33 @@
+#include<yade/pkg/dem/RungeKuttaCashKarp54Integrator.hpp>
+#include<yade/core/Scene.hpp>
+
+YADE_PLUGIN((RungeKuttaCashKarp54Integrator));
+
+shared_ptr<Integrator> RungeKuttaCashKarp54Integrator_ctor_list(const python::list& slaves){ shared_ptr<Integrator> instance(new RungeKuttaCashKarp54Integrator); instance->slaves_set(slaves); return instance; }
+
+void RungeKuttaCashKarp54Integrator::action()
+{
+
+ Real dt=scene->dt;
+
+ Real time=scene->time;
+
+ error_checker_type rungekuttaerrorcontroller=error_checker_type(abs_err,rel_err,a_x,a_dxdt);
+
+ controlled_stepper_type rungekuttastepper=controlled_stepper_type(rungekuttaerrorcontroller);
+
+ stateVector currentstates=getCurrentStates();
+
+ resetstate.reserve(currentstates.size());
+
+ copy(currentstates.begin(),currentstates.end(),back_inserter(resetstate));//copy current state to resetstate
+
+ this->timeresetvalue=time; //set reset time to the time just before the integration
+
+ /*Try an adaptive integration*/
+
+ integrationsteps+=integrate_adaptive(rungekuttastepper,make_ode_wrapper( *((Integrator*)this) , &Integrator::system ),currentstates,time,time+dt, stepsize, observer(this));
+
+ scene->time=scene->time-dt;//Scene move next time step function already increments the time so we have to decrement it just before it.
+}
+
=== added file 'pkg/dem/RungeKuttaCashKarp54Integrator.hpp'
--- pkg/dem/RungeKuttaCashKarp54Integrator.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/RungeKuttaCashKarp54Integrator.hpp 2014-02-17 14:31:35 +0000
@@ -0,0 +1,50 @@
+#pragma once
+#include <yade/core/Scene.hpp>
+#include<yade/pkg/dem/Integrator.hpp>
+#include<boost/numeric/odeint.hpp>
+
+
+typedef boost::numeric::odeint::runge_kutta_cash_karp54< stateVector > error_stepper_type; //Runge-Kutta 54 error stepper other steppers can also be used
+
+typedef boost::numeric::odeint::controlled_runge_kutta< error_stepper_type > controlled_stepper_type;//Controlled Runge Kutta stepper
+
+typedef boost::numeric::odeint::default_error_checker< error_stepper_type::value_type,error_stepper_type::algebra_type ,error_stepper_type::operations_type > error_checker_type; //Error checker type that is redefined for initialization using different tolerance values
+
+
+shared_ptr<Integrator> RungeKuttaCashKarp54Integrator_ctor_list(const python::list& slaves);
+class RungeKuttaCashKarp54Integrator: public Integrator {
+
+ public:
+
+ error_checker_type rungekuttaerrorcontroller;
+
+ controlled_stepper_type rungekuttastepper;
+
+ void init()
+ {
+ rungekuttaerrorcontroller=error_checker_type(abs_err,rel_err,a_x,a_dxdt);
+ rungekuttastepper=controlled_stepper_type(rungekuttaerrorcontroller);
+ }
+
+
+
+ virtual void action();
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(RungeKuttaCashKarp54Integrator,Integrator,"RungeKuttaCashKarp54Integrator engine.",
+ ((Real,abs_err,1e-6,,"Relative integration tolerance"))
+ ((Real,rel_err,1e-6,,"Absolute integration tolerance"))
+ ((Real,a_x,1.0,,""))
+ ((Real,a_dxdt,1.0,,""))
+ ((Real,stepsize,1e-6,,"It is not important for an adaptive integration but important for the observer for setting the found states after integration"))
+ ,
+ /*ctor*/
+ init();
+ ,
+ .def("__init__",python::make_constructor(RungeKuttaCashKarp54Integrator_ctor_list),"Construct from (possibly nested) list of slaves.")
+ /*py*/
+ );
+};
+REGISTER_SERIALIZABLE(RungeKuttaCashKarp54Integrator);
+
+
+
