Commit cacc9589 authored by Václav Šmilauer's avatar Václav Šmilauer
Browse files

1. Fix boost::sertialization class export (changed slightly in boost 1.42,...

1. Fix boost::sertialization class export (changed slightly in boost 1.42, which broke our case); that also solves the issues at maverick
2. Add framework for tracking energies (will be documented at some point), added to Law2_ScGeom_FrictPhys_CundallStrack, see scripts/test/energy.py
3. Convert energy trackers in Hertz-Mindlin to OpenMPAccumulator
4. Make OpenMPAccumulator align storage so that each thread uses one cache line (should be faster)
5. Make OpenMPAccumulator not require zero value pointer (is handled via ZeroInitializer template in lib/base/Math.hpp
6. Make OpenMPAccumulator work transparently with python (it can be used as a regular attribute now), make it boost::serializable as well; it appears transparently as a number in python.
7. Add OpenMPArrayAccumulator, for linear array of values (used in EnergyTracker)
8. Make deprecated attributes with non-g++ compilers (clang); previously, only g++>=4.3 was supported
9. Fix a few example script (not completed)
parent ed1177a2
......@@ -3,8 +3,6 @@
#include<yade/core/BodyContainer.hpp>
#include<yade/core/Body.hpp>
BOOST_CLASS_EXPORT(BodyContainer);
unsigned int BodyContainer::findFreeId(){
unsigned int max=body.size();
for(; lowestFree<max; lowestFree++){
......
#pragma once
#include<boost/python.hpp>
#include<boost/foreach.hpp>
#include<string>
#include<yade/lib-base/openmp-accu.hpp>
#include<yade/lib-serialization/Serializable.hpp>
#ifndef FOREACH
#define FOREACH BOOST_FOREACH
#endif
namespace py=boost::python;
class EnergyTracker: public Serializable{
public:
~EnergyTracker();
void findId(const std::string& name, int& id, bool reset=false, bool newIfNotFound=true){
if(id>0) return; // the caller should have checked this already
if(names.count(name)) id=names[name];
else if(newIfNotFound) {
#ifdef YADE_OPENMP
#pragma omp critical
#endif
{ energies.resize(energies.size()+1); id=energies.size()-1; resetStep.resize(id+1); resetStep[id]=reset; names[name]=id; assert(id<(int)energies.size()); assert(id>=0); }
}
}
// set value of the accumulator; note: must NOT be called from parallel sections!
void set(const Real& val, const std::string& name, int &id){
if(id<0) findId(name,id,/* do not reset value that is set directly */ false);
energies.set(id,val);
}
// add value to the accumulator; safely called from parallel sections
void add(const Real& val, const std::string& name, int &id, bool reset=false){
if(id<0) findId(name,id,reset);
energies.add(id,val);
}
Real getItem_py(const std::string& name){
int id=-1; findId(name,id,false,false);
if (id<0) {PyErr_SetString(PyExc_KeyError,("Unknown energy name '"+name+"'.").c_str()); python::throw_error_already_set(); }
return energies.get(id);
}
void setItem_py(const std::string& name, Real val){
int id=-1; set(val,name,id);
}
void clear(){ energies.clear(); names.clear(); resetStep.clear();}
py::list keys_py(){ py::list ret; FOREACH(pairStringInt p, names) ret.append(p.first); return ret; }
void resetResettables(){ size_t sz=energies.size(); for(size_t id=0; id<sz; id++){ if(resetStep[id]) energies.reset(id); } }
typedef std::map<std::string,int> mapStringInt;
typedef std::pair<std::string,int> pairStringInt;
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(EnergyTracker,Serializable,"Storage for tracing energies. Only to be used if O.traceEnergy is True.",
((OpenMPArrayAccumulator<Real>,energies,,,"Energy values, in linear array"))
((mapStringInt,names,,Attr::hidden,"Associate textual name to an index in the energies array."))
((vector<bool>,resetStep,,Attr::hidden,"Whether the respective energy value should be reset at every step."))
,/*ctor*/
,/*py*/
.def("__getitem__",&EnergyTracker::getItem_py,"Get energy value for given name.")
.def("__setitem__",&EnergyTracker::setItem_py,"Set energy value for given name (will create a non-resettable item, if it does not exist yet).")
.def("clear",&EnergyTracker::clear,"Clear all stored values.")
.def("keys",&EnergyTracker::keys_py,"Return defined energies.")
)
};
REGISTER_SERIALIZABLE(EnergyTracker);
......@@ -7,8 +7,6 @@
#include<omp.h>
#endif
BOOST_CLASS_EXPORT(InteractionContainer);
//BOOST_SERIALIZATION_FACTORY_0(InteractionContainer);
bool InteractionContainer::insert(const shared_ptr<Interaction>& i){
boost::mutex::scoped_lock lock(drawloopmutex);
......
......@@ -10,16 +10,15 @@
#pragma once
#include"Body.hpp"
#include"Cell.hpp"
#include"BodyContainer.hpp"
#include"Engine.hpp"
#include"Material.hpp"
#include"DisplayParameters.hpp"
#include"ForceContainer.hpp"
#include"InteractionContainer.hpp"
#include"GroupRelationData.hpp"
#include<yade/core/Body.hpp>
#include<yade/core/Cell.hpp>
#include<yade/core/BodyContainer.hpp>
#include<yade/core/Engine.hpp>
#include<yade/core/Material.hpp>
#include<yade/core/DisplayParameters.hpp>
#include<yade/core/ForceContainer.hpp>
#include<yade/core/InteractionContainer.hpp>
#include<yade/core/EnergyTracker.hpp>
#ifndef HOST_NAME_MAX
#define HOST_NAME_MAX 255
......@@ -73,6 +72,7 @@ class Scene: public Serializable{
((Real,stopAtRealTime,0,,"Time at which to stop the simulation"))
#endif
((bool,isPeriodic,false,Attr::readonly,"Whether periodic boundary conditions are active."))
((bool,trackEnergy,false,Attr::readonly,"Whether energies are being traced."))
((bool,needsInitializers,true,Attr::readonly,"Whether initializers will be run before the first step."))
((Body::id_t,selectedBody,-1,,"Id of body that is selected by the user"))
......@@ -81,6 +81,7 @@ class Scene: public Serializable{
((vector<shared_ptr<Engine> >,initializers,,Attr::hidden,"Engines that will be run only once, before the first step."))
((shared_ptr<BodyContainer>,bodies,new BodyContainer,Attr::hidden,"Bodies contained in the scene."))
((shared_ptr<InteractionContainer>,interactions,new InteractionContainer,Attr::hidden,"All interactions between bodies."))
((shared_ptr<EnergyTracker>,energy,new EnergyTracker,Attr::hidden,"Energy values, if energy tracking is enabled."))
((vector<shared_ptr<Material> >,materials,,Attr::hidden,"Container of shared materials. Add elements using Scene::addMaterial, not directly. Do NOT remove elements from here unless you know what you are doing!"))
((shared_ptr<Bound>,bound,,Attr::hidden,"Bounding box of the scene (only used for rendering and initialized if needed)."))
......@@ -93,5 +94,5 @@ class Scene: public Serializable{
);
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(Scene);
#include<yade/lib-factory/ClassFactory.hpp>
// make core classes known to the class factory
#include<yade/core/Body.hpp>
#include<yade/core/BodyContainer.hpp>
#include<yade/core/Bound.hpp>
#include<yade/core/Cell.hpp>
#include<yade/core/Dispatcher.hpp>
#include<yade/core/EnergyTracker.hpp>
#include<yade/core/Engine.hpp>
#include<yade/core/FileGenerator.hpp>
#include<yade/core/Functor.hpp>
#include<yade/core/GlobalEngine.hpp>
#include<yade/core/Interaction.hpp>
#include<yade/core/InteractionContainer.hpp>
#include<yade/core/IGeom.hpp>
#include<yade/core/IPhys.hpp>
#include<yade/core/Material.hpp>
......@@ -16,4 +19,14 @@
#include<yade/core/Shape.hpp>
#include<yade/core/State.hpp>
#include<yade/core/TimeStepper.hpp>
YADE_PLUGIN((Body)(Bound)(Cell)(Dispatcher)(Engine)(FileGenerator)(Functor)(GlobalEngine)(Interaction)(IGeom)(IPhys)(Material)(PartialEngine)(Shape)(State)(TimeStepper));
// these two are not accessible from python directly (though they should be in the future, perhaps)
BOOST_CLASS_EXPORT_IMPLEMENT(BodyContainer);
BOOST_CLASS_EXPORT_IMPLEMENT(InteractionContainer);
YADE_PLUGIN((Body)(Bound)(Cell)(Dispatcher)(EnergyTracker)(Engine)(FileGenerator)(Functor)(GlobalEngine)(Interaction)(IGeom)(IPhys)(Material)(PartialEngine)(Shape)(State)(TimeStepper));
EnergyTracker::~EnergyTracker(){} // vtable
//BOOST_CLASS_EXPORT(OpenMPArrayAccumulator<Real>);
//BOOST_CLASS_EXPORT(OpenMPAccumulator<Real>);
......@@ -58,9 +58,7 @@ if opts.debug:
else:
if not hasNonDebug:
print 'WARNING: non-debug build not available, running with --debug instead (try --rebuild to get the non-debug build).'
#raise RuntimeError('Non-debug build not available (run with --debug, or try --rebuild)')
else:
libDir=nonDebugLibDir
libDir=nonDebugLibDir
## remove later
## python2.5 relative module imports workaround
......@@ -119,7 +117,7 @@ if opts.updateScripts:
sys.exit(0)
if opts.manpage:
import yade.manpage
yade.manpage.generate_manpage(par,yade.config.metadata,opts.manpage,section=1,seealso='yade%s-multi (1)'%suffix)
yade.manpage.generate_manpage(par,yade.config.metadata,opts.manpage,section=1,seealso='yade%s-batch (1)'%suffix)
print 'Manual page %s generated.'%opts.manpage
sys.exit(0)
if opts.nice:
......
......@@ -50,10 +50,8 @@ colorsph2=Vector3(1,1,0);
colorsph1.normalize();
colorsph2.normalize();
colorSph=colorsph1
#O.bodies.append(utils.sphere([0,0,0],1))
for xyz in itertools.product(arange(0,numBoxes[0]),arange(0,numBoxes[1]),arange(0,numBoxes[2])):
continue
#ids_spheres=O.bodies.appendClumped(pack.regularHexa(pack.inEllipsoid((xyz[0]*(sizeBox+gapBetweenBoxes),xyz[1]*(sizeBox+gapBetweenBoxes)+sizeBox*0.5,xyz[2]*(sizeBox+gapBetweenBoxes)-radiusKnife+sizeBox*0.6),(sizeBox/2,sizeBox/2,sizeBox/2)),radius=radiusSph,gap=0,color=colorSph))
ids_spheres=O.bodies.appendClumped(pack.regularHexa(pack.inEllipsoid((xyz[0]*(sizeBox+gapBetweenBoxes),xyz[1]*(sizeBox+gapBetweenBoxes)+sizeBox*0.5,xyz[2]*(sizeBox+gapBetweenBoxes)-radiusKnife+sizeBox*0.6),(sizeBox/2,sizeBox/2,sizeBox/2)),radius=radiusSph,gap=0,color=colorSph))
if (colorSph==colorsph1):
colorSph=colorsph2
else:
......@@ -69,15 +67,13 @@ O.engines=[
InteractionLoop(
[Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_Dem3DofGeom_FrictPhys_Basic()],
[Law2_Dem3DofGeom_FrictPhys_CundallStrack()],
),
GravityEngine(gravity=(0,0,-9.8)),
TranslationEngine(translationAxis=[1,0,0],velocity=5,ids=KnifeIDs), # Buldozer motion
NewtonIntegrator(damping=.3),
#qt.SnapshotEngine(iterPeriod=100,fileBase='/tmp/bulldozer-',label='snapshooter'),
]
O.engines=[]
print len(O.bodies)
O.saveTmp()
qt.Controller()
......
......@@ -7,7 +7,7 @@ O.engines=[
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_Basic()]
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
GravityEngine(gravity=(0,0,-1000)),
NewtonIntegrator(damping=0.2)
......
......@@ -3,17 +3,14 @@
from math import *
# "instantiate" Omega, i.e. create proxy object to Omega and rootBody
o=Omega()
o.engines=[
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
IGeomDispatcher([Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()]),
IPhysDispatcher([Ip2_FrictMat_FrictMat_FrictPhys()]),
ElasticContactLaw(),
GlobalStiffnessTimeStepper(defaultDt=1e-4,active=True,timeStepUpdateInterval=500),
AxialGravityEngine(axisPoint=(0,0,0),axisDirection=(1,0,0),acceleration=1e4),
#GlobalStiffnessTimeStepper(defaultDt=1e-5,active=True,timeStepUpdateInterval=500),
AxialGravityEngine(axisPoint=(0,0,0),axisDirection=(1,0,0),acceleration=1e3),
NewtonIntegrator(damping=.4)
]
......@@ -29,18 +26,18 @@ for n in range(nMax):
O.bodies.append(utils.sphere([x,0,0],rCenter))
#lateral walls, they have wallDist gap inbetween
wLat1=utils.box([0+.5*wallDist+.5*wallThickness,0,0],[.5*wallThickness,wallSize,wallSize]); o.bodies.append(wLat1);
wLat2=utils.box([0-.5*wallDist-.5*wallThickness,0,0],[.5*wallThickness,wallSize,wallSize]); o.bodies.append(wLat2);
wLat1=utils.box([0+.5*wallDist+.5*wallThickness,0,0],[.5*wallThickness,wallSize,wallSize]); O.bodies.append(wLat1);
wLat2=utils.box([0-.5*wallDist-.5*wallThickness,0,0],[.5*wallThickness,wallSize,wallSize]); O.bodies.append(wLat2);
#angle walls, they cross at the x-axis
wAng1=utils.box([0,0,0],[.55*wallDist,.5*wallThickness,wallSize*sqrt(2)]); wAng1.state.ori=Quaternion((1,0,0),pi/4); o.bodies.append(wAng1);
wAng2=utils.box([0,0,0],[.55*wallDist,.5*wallThickness,wallSize*sqrt(2)]); wAng2.state.ori=Quaternion((1,0,0),-pi/4); o.bodies.append(wAng2)
wAng1=utils.box([0,0,0],[.55*wallDist,.5*wallThickness,wallSize*sqrt(2)]); wAng1.state.ori=Quaternion((1,0,0),pi/4); O.bodies.append(wAng1);
wAng2=utils.box([0,0,0],[.55*wallDist,.5*wallThickness,wallSize*sqrt(2)]); wAng2.state.ori=Quaternion((1,0,0),-pi/4); O.bodies.append(wAng2)
#cap
wCap=utils.box([0,0,wallSize],[.55*wallDist,wallSize,.5*wallThickness]); o.bodies.append(wCap)
wCap=utils.box([0,0,wallSize],[.55*wallDist,wallSize,.5*wallThickness]); O.bodies.append(wCap)
# all bodies up to now are fixed and only wire is will be shown
for b in o.bodies:
for b in O.bodies:
b.shape.wire=True
b.dynamic=False
......@@ -57,7 +54,9 @@ for ix in range(int(maxima[0])):
#print x,y,z,rBall
O.bodies.append(utils.sphere([x,y,z],rBall))
o.save('/tmp/a.xml')
O.dt=.1*utils.PWaveTimeStep()
O.saveTmp()
try:
from yade import qt
......@@ -65,12 +64,3 @@ try:
qt.View()
except ImportError: pass
if 0:
import os,time
os.system(yadeExecutable+' -N QtGUI -S /tmp/a.xml')
else:
o.run(30000);
o.wait()
o.save('/tmp/a_30000.xml')
#os.system(yadeExecutable+' -N QtGUI -S /tmp/a_10000.xml')
......@@ -46,7 +46,7 @@ o.engines=[
InteractionLoop(
[Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_Dem3DofGeom_FrictPhys_Basic()],
[Law2_Dem3DofGeom_FrictPhys_CundallStrack()],
),
## Apply gravity
GravityEngine(gravity=[0,-9.81,0]),
......
......@@ -22,7 +22,7 @@ O.engines=[
[Ig2_Sphere_Sphere_Dem3DofGeom(),
Ig2_Facet_Sphere_Dem3DofGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_Dem3DofGeom_FrictPhys_Basic()]
[Law2_Dem3DofGeom_FrictPhys_CundallStrack()]
),
GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(),
......@@ -34,5 +34,5 @@ O.dt=utils.PWaveTimeStep()
from yade import qt
qt.Controller()
qt.View()
O.saveTmp()
O.run(8500)
......@@ -25,6 +25,8 @@ surf=gts.read(open('horse.coarse.gts'))
O.materials.append(FrictMat(young=30e9,density=2000))
print 'surface closed',surf.is_closed(),'volume',surf.volume()
if surf.is_closed():
pred=pack.inGtsSurface(surf)
aabb=pred.aabb()
......@@ -39,7 +41,7 @@ O.engines=[
InteractionLoop(
[Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_Dem3DofGeom_FrictPhys_Basic()],
[Law2_Dem3DofGeom_FrictPhys_CundallStrack()],
),
GravityEngine(gravity=[0,0,-1e4]),
NewtonIntegrator(damping=.1)
......
......@@ -58,7 +58,7 @@ for i in (0,1): # red and blue spheres
print "Numer of grains",len(O.bodies)-len(millIds)
O.dt=utils.PWaveTimeStep()
O.dt=.5*utils.PWaveTimeStep()
O.engines=[
ForceResetter(),
......@@ -66,18 +66,18 @@ O.engines=[
InteractionLoop(
[Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_Dem3DofGeom_FrictPhys_Basic()],
[Law2_Dem3DofGeom_FrictPhys_CundallStrack()],
),
GravityEngine(gravity=(0,0,-50)), # gravity artificially high, to make it faster going ;-)
NewtonIntegrator(damping=.3),
RotationEngine(rotateAroundZero=True,zeroPoint=(0,0,0),rotationAxis=(1,0,0),angularVelocity=-20,subscribedBodies=millIds),
RotationEngine(rotateAroundZero=True,zeroPoint=(0,0,0),rotationAxis=(1,0,0),angularVelocity=-20,ids=millIds),
#SnapshotEngine(iterPeriod=30,fileBase='/tmp/mill-',viewNo=0,label='snapshooter'),
VTKRecorder(iterPeriod=100,recorders=['all'],fileName='/tmp/millVTK-',multiblock=False)
#VTKRecorder(iterPeriod=100,recorders=['all'],fileName='/tmp/millVTK-')
]
O.saveTmp()
#from yade import qt
#v=qt.View()
#v.eyePosition=(3,.8,.96); v.upVector=(-.4,-.4,.8); v.viewDir=(-.9,-.25,-.3); v.axes=True; v.sceneRadius=1.9
O.run(20000); O.wait()
from yade import qt
v=qt.View()
v.eyePosition=(3,.8,.96); v.upVector=(-.4,-.4,.8); v.viewDir=(-.9,-.25,-.3); v.axes=True; v.sceneRadius=1.9
#O.run(20000); O.wait()
#utils.encodeVideoFromFrames(snapshooter['savedSnapshots'],out='/tmp/mill.ogg',fps=30)
......@@ -86,7 +86,7 @@ O.engines=[
InteractionLoop(
[Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_Dem3DofGeom_FrictPhys_Basic()],
[Law2_Dem3DofGeom_FrictPhys_CundallStrack()],
),
GravityEngine(gravity=(1e-2,1e-2,-1000)),
NewtonIntegrator(damping=.1,exactAsphericalRot=True),
......
......@@ -12,9 +12,6 @@ mesh = 'coarse'
#mesh = 'fine'
#mesh = 'tiny'
## Omega
o=Omega()
## Import geometry
rod = O.bodies.append(ymport.stl('rod-'+mesh+'.stl',wire=True))
......@@ -34,16 +31,16 @@ for i in xrange(nbSpheres[0]):
if (i==0 or i==nbSpheres[0]-1 or j==nbSpheres[1]-1 or k==0 or k==nbSpheres[2]-1):
dynamic = False
color=[0.21,0.22,0.1]
o.bodies.append(utils.sphere([x,y,z],r,color=color,dynamic=dynamic))
O.bodies.append(utils.sphere([x,y,z],r,color=color,dynamic=dynamic))
print "done\n"
## Estimate time step
#o.dt=utils.PWaveTimeStep()
o.dt=0.0001
#O.dt=utils.PWaveTimeStep()
O.dt=0.0001
## Engines
o.engines=[
O.engines=[
## Resets forces and momenta the act on bodies
ForceResetter(),
......@@ -53,18 +50,20 @@ o.engines=[
Bo1_Facet_Aabb(),
]),
InteractionLoop(
[Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
#[Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_Dem3DofGeom_FrictPhys_Basic()],
#[Law2_Dem3DofGeom_FrictPhys_CundallStrack()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
## Apply gravity
GravityEngine(gravity=[0,-9.81,0]),
## Motion equation
NewtonIntegrator(damping=0.3),
## Apply kinematics to rod
TranslationEngine(subscribedBodies=rod,translationAxis=[0,-1,0],velocity=0.075),
TranslationEngine(ids=rod,translationAxis=[0,-1,0],velocity=0.075),
## Save force on rod
ForceRecorder(subscribedBodies=rod,file='force-'+mesh+'.dat',iterPeriod=50),
#ForceRecorder(ids=rod,file='force-'+mesh+'.dat',iterPeriod=50),
]
......@@ -76,9 +75,9 @@ nbIter=10000
from yade import qt
qt.View()
for t in xrange(2):
start=time.time();o.run(nbIter);o.wait();finish=time.time()
speed=nbIter/(finish-start); print '%g iter/sec\n'%speed
print "FINISH"
quit()
#for t in xrange(2):
# start=time.time();O.run(nbIter);O.wait();finish=time.time()
# speed=nbIter/(finish-start); print '%g iter/sec\n'%speed
#print "FINISH"
#quit()
......@@ -47,10 +47,10 @@ O.engines=[
InteractionLoop(
[Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_Dem3DofGeom_FrictPhys_Basic()],
[Law2_Dem3DofGeom_FrictPhys_CundallStrack()],
),
GravityEngine(gravity=(0,0,-1e3)), # gravity artificially high, to make it faster going ;-)
RotationEngine(rotateAroundZero=True,zeroPoint=(0,0,0),rotationAxis=(0,1,1),angularVelocity=30*(2*pi/60),subscribedBodies=cylIds,label='rotor'),
RotationEngine(rotateAroundZero=True,zeroPoint=(0,0,0),rotationAxis=(0,1,1),angularVelocity=30*(2*pi/60),ids=cylIds,label='rotor'),
NewtonIntegrator(damping=.3),
]
O.dt=utils.PWaveTimeStep()
......
......@@ -9,7 +9,7 @@ O.engines=[
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_Basic()]
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(damping=.2),
......@@ -46,7 +46,7 @@ def myAddPlotData():
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()
plot.plot(subPlots=True)
O.run(int(2./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.
......@@ -259,8 +259,10 @@ class ControllerClass(QWidget,Ui_Controller):
e=int((stopAtIter-iter)*self.iterPerSec)
self.realTimeLabel.setText('%02d:%02d:%02d (ETA %02d:%02d:%02d)'%(rt//3600,rt//60,rt%60,e//3600,e//60,e%60))
self.iterLabel.setText('#%ld / %ld, %.1f/s %s'%(O.iter,stopAtIter,self.iterPerSec,subStepInfo))
s=int(t); ms=int(t*1000)%1000; us=int(t*1000000)%1000; ns=int(t*1000000000)%1000
self.virtTimeLabel.setText(u'%03ds%03dm%03dμ%03dn'%(s,ms,us,ns))
if t!=float('inf'):
s=int(t); ms=int(t*1000)%1000; us=int(t*1000000)%1000; ns=int(t*1000000000)%1000
self.virtTimeLabel.setText(u'%03ds%03dm%03dμ%03dn'%(s,ms,us,ns))
else: self.virtTimeLabel.setText(u'[ ∞ ] ?!')
self.show3dButton.setChecked(len(views())>0)
def Generator():
......
......@@ -7,3 +7,6 @@ template<> const Real Math<Real>::TWO_PI = 2.0*Math<Real>::PI;
template<> const Real Math<Real>::HALF_PI = 0.5*Math<Real>::PI;
template<> const Real Math<Real>::DEG_TO_RAD = Math<Real>::PI/180.0;
template<> const Real Math<Real>::RAD_TO_DEG = 180.0/Math<Real>::PI;
template<> int ZeroInitializer<int>(){ return (int)0; }
template<> Real ZeroInitializer<Real>(){ return (Real)0; }
......@@ -79,6 +79,13 @@ typedef Eigen::Quaternion<Real> Quaternionr;
typedef Eigen::AngleAxis<Real> AngleAxisr;
using Eigen::AngleAxis; using Eigen::Quaternion;
// in some cases, we want to initialize types that have no default constructor (OpenMPAccumulator, for instance)
// template specialization will help us here
template<typename EigenMatrix> EigenMatrix ZeroInitializer(){ return EigenMatrix::Zero(); };
template<> int ZeroInitializer<int>();
template<> Real ZeroInitializer<Real>();
// io
template<class Scalar> std::ostream & operator<<(std::ostream &os, const VECTOR2_TEMPLATE(Scalar)& v){ os << v.x()<<" "<<v.y(); return os; };
template<class Scalar> std::ostream & operator<<(std::ostream &os, const VECTOR3_TEMPLATE(Scalar)& v){ os << v.x()<<" "<<v.y()<<" "<<v.z(); return os; };
......
// 2010 © Václav Šmilauer <eudoxos@arcig.cz>
#pragma once
// for ZeroInitializer template
#include<yade/lib-base/Math.hpp>
#include<boost/serialization/split_free.hpp>
#include<boost/lexical_cast.hpp>
#include<string>
#ifdef YADE_OPENMP
#include"omp.h"
#include<cstdlib>
#include<unistd.h>
#include<vector>
#include<stdexcept>
#include<iostream>
// O(1) access container which stores data in contiguous chunks of memory
// each chunk belonging to one thread
template<typename T>
class OpenMPArrayAccumulator{
int CLS;
size_t nThreads;
int perCL; // number of elements fitting inside cache line
std::vector<T*> chunks; // array with pointers to the chunks of memory we have allocated; each item for one thread
size_t sz; // current number of elements
size_t nCL; // current number of allocated cache lines
int nCL_for_N(size_t n){ return n/perCL+(n%perCL==0 ? 0 : 1); } // return number of cache lines to allocate for given number of elements
public:
OpenMPArrayAccumulator() : CLS(sysconf(_SC_LEVEL1_DCACHE_LINESIZE)), nThreads(omp_get_max_threads()), perCL(CLS/sizeof(T)), chunks(nThreads,NULL), sz(0), nCL(0) { }
OpenMPArrayAccumulator(size_t n): CLS(sysconf(_SC_LEVEL1_DCACHE_LINESIZE)), nThreads(omp_get_max_threads()), perCL(CLS/sizeof(T)), chunks(nThreads,NULL), sz(0), nCL(0) { resize(n); }
// change number of elements
void resize(size_t n){
if(n==sz) return; // nothing to do
size_t nCL_new=nCL_for_N(n);
if(nCL_new>nCL){
for(size_t th=0; th<nThreads; th++){
void* oldChunk=(void*)chunks[th];
// FIXME: currently we allocate 4× the memory necessary, otherwise there is crash when accessing past its half -- why?? (http://www.abclinuxu.cz/poradna/programovani/show/318324)
int succ=posix_memalign((void**)(&chunks[th]),/*alignment*/CLS,/*size*/ nCL_new*CLS);
if(succ!=0) throw std::runtime_error("OpenMPArrayAccumulator: posix_memalign failed to allocate memory.");