yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10269
[Branch ~yade-pkg/yade/git-trunk] Rev 3752: Merge branch 'master' of github.com:yade/trunk
Merge authors:
Anton Gladky (gladky-anton)
Anton Gladky (gladky-anton)
Bruno Chareyre (bruno-chareyre)
Jan Stránský (honzik)
------------------------------------------------------------
revno: 3752 [merge]
committer: Jerome Duriez <jerome.duriez@xxxxxxxxxxxxxxx>
timestamp: Thu 2013-11-28 10:45:15 +0100
message:
Merge branch 'master' of github.com:yade/trunk
added:
examples/baraban/BicyclePedalEngine.py
pkg/common/Grid_GUI.cpp
scripts/ppa/
scripts/ppa/createtar.py
modified:
CMakeLists.txt
README.rst
doc/sphinx/installation.rst
lib/triangulation/FlowBoundingSphere.ipp
lib/triangulation/FlowBoundingSphereLinSolv.ipp
lib/triangulation/PeriodicFlow.cpp
lib/triangulation/def_types.h
pkg/common/Grid.cpp
pkg/common/KinematicEngines.cpp
pkg/common/KinematicEngines.hpp
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.hpp
py/_utils.cpp
scripts/checks-and-tests/checks/DEM-PFV-check.py
--
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
=== modified file 'CMakeLists.txt'
--- CMakeLists.txt 2013-10-30 07:54:30 +0000
+++ CMakeLists.txt 2013-11-21 22:57:34 +0000
@@ -40,15 +40,15 @@
INCLUDE(GNUInstallDirs)
#===========================================================
# HACK!!! If the version of gcc is 4.8 or greater, we add -ftrack-macro-expansion=0
-# to compiler to reduce the memory consumption during compilation step.
+# and -save-temps into compiler to reduce the memory consumption during compilation.
# See http://bugs.debian.org/726009 for more information
# Can be removed later, if gcc fixes its regression
# Taken from http://stackoverflow.com/questions/4058565/check-gcc-minor-in-cmake
EXECUTE_PROCESS(COMMAND ${CMAKE_C_COMPILER} -dumpversion OUTPUT_VARIABLE GCC_VERSION)
IF (GCC_VERSION VERSION_GREATER 4.8 OR GCC_VERSION VERSION_EQUAL 4.8)
- MESSAGE(STATUS "GCC Version >= 4.8. Adding -ftrack-macro-expansion=0")
- SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ftrack-macro-expansion=0 ")
+ MESSAGE(STATUS "GCC Version >= 4.8. Adding -ftrack-macro-expansion=0 and -save-temps")
+ SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ftrack-macro-expansion=0 -save-temps")
ENDIF()
#===========================================================
@@ -259,7 +259,7 @@
MESSAGE(STATUS "Found Metis")
SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} LinSolv")
ELSE(CHOLMOD_FOUND AND OPENBLAS_FOUND AND METIS_FOUND)
- MESSAGE(STATUS "CHOLMOD NOT found, LINSOLV disabled")
+ MESSAGE(STATUS "Missing dependency for LINSOLV, disabled")
SET(DISABLED_FEATS "${DISABLED_FEATS} LinSolv")
SET(ENABLE_LINSOLV OFF)
ENDIF(CHOLMOD_FOUND AND OPENBLAS_FOUND AND METIS_FOUND)
=== modified file 'README.rst'
--- README.rst 2013-04-26 20:00:56 +0000
+++ README.rst 2013-11-10 21:23:19 +0000
@@ -18,7 +18,6 @@
- Packages for Ubuntu:
- - Stable-builds: https://launchpad.net/~yade-pkg/+archive/stable
- Daily-builds: https://launchpad.net/~yade-pkg/+archive/snapshots
- External packages: https://launchpad.net/~yade-users/+archive/external
=== modified file 'doc/sphinx/installation.rst'
--- doc/sphinx/installation.rst 2013-11-01 06:25:52 +0000
+++ doc/sphinx/installation.rst 2013-11-26 14:29:06 +0000
@@ -82,7 +82,8 @@
Prerequisites
^^^^^^^^^^^^^
-Yade relies on a number of external software to run; its installation is checked before the compilation starts.
+Yade relies on a number of external software to run; they are checked before the compilation starts.
+Some of them are only optional. The last ones are only relevant for using the fluid coupling module (:yref:`FlowEngine`).
* `cmake <http://www.cmake.org/>`_ build system
* `gcc <http://www.gcc.gnu.org>`_ compiler (g++); other compilers will not work; you need g++>=4.2 for openMP support
@@ -98,13 +99,20 @@
* `Loki <http://loki-lib.sf.net>`_ library
* `VTK <http://www.vtk.org/>`_ library (optional but recommended)
* `CGAL <http://www.cgal.org/>`_ library (optional)
+* `SuiteSparse <http://www.cise.ufl.edu/research/sparse/SuiteSparse/>`_ sparse algebra library (fluid coupling, optional, requires eigen>=3.1)
+* `OpenBLAS <http://www.openblas.net/>`_ optimized and parallelized alternative to the standard blas+lapack (fluid coupling, optional)
+* `Metis <http://glaros.dtc.umn.edu/gkhome/metis/metis/overview/>`_ matrix preconditioning (fluid coupling, optional)
Most of the list above is very likely already packaged for your distribution.
The following commands have to be executed in command line of corresponding
distributions. Just copy&paste to the terminal. To perform commands you
should have root privileges
- * **Ubuntu**, **Debian** and their derivatives::
+.. warning:: if you have Ubuntu 12.10 or older, you need to install libqglviewer-qt4-dev
+ package instead of libqglviewer-dev.
+
+
+* **Ubuntu**, **Debian** and their derivatives::
sudo apt-get install cmake git freeglut3-dev libloki-dev \
libboost-all-dev fakeroot dpkg-dev build-essential g++ \
@@ -112,17 +120,36 @@
libgts-dev python-pygraphviz libvtk5-dev python-scientific libeigen3-dev \
python-xlib python-qt4 pyqt4-dev-tools gtk2-engines-pixbuf python-argparse \
libqglviewer-dev python-imaging libjs-jquery python-sphinx python-git python-bibtex \
- libxmu-dev libxi-dev libcgal-dev help2man libsuitesparse-dev \
- libopenblas-dev libmetis-dev python-gts libbz2-dev zlib1g-dev
-
-If you are using other distribuition, than Debian or its derivatives, you should
-install the software, which is listed above. Their names can differ from the
-names of Debian-packages.
+ libxmu-dev libxi-dev libcgal-dev help2man libbz2-dev zlib1g-dev \
+
Some of packages (for example, cmake, eigen3) are mandatory, some of them
are optional. Watch for notes and warnings/errors, which are shown
by cmake during configuration step. If the missing package is optional,
some of Yade features will be disabled (see the messages at the end of configuration).
+
+Additional packages, which can become mandatory later::
+
+ sudo apt-get install python-gts python-minieigen \
+
+For effective usage of direct solvers in the fluid coupling, the following libraries are recommended, together with eigen>=3.1: blas, lapack, suitesparse, and metis.
+All four of them are available in many different versions. Different combinations are possible and not all of them will work. The following was found to be effective on ubuntu 12.04 and debian wheezy.
+(openblas provides its own version of lapack, and suitesparse-metis will trigger the installation of parmetis)::
+
+ sudo apt-get install libopenblas-dev libsuitesparse-metis-dev \
+
+Some packages listed here are relatively new and they can be absent
+in your distribution (for example, libmetis-dev or python-gts). They can be
+installed from our `external PPA <https://launchpad.net/~yade-users/+archive/external/>`_
+or just ignored. In this case some features can be disabled.
+
+If you are using other distribuition, than Debian or its derivatives, you should
+install the softwares listed above. Their names can differ from the
+names of Debian-packages.
+
+
+
+
Compilation
=== added file 'examples/baraban/BicyclePedalEngine.py'
--- examples/baraban/BicyclePedalEngine.py 1970-01-01 00:00:00 +0000
+++ examples/baraban/BicyclePedalEngine.py 2013-11-11 16:22:15 +0000
@@ -0,0 +1,62 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+import time
+
+## PhysicalParameters
+Density=2400
+frictionAngle=radians(35)
+tc = 0.001
+en = 0.3
+es = 0.3
+
+## Import wall's geometry
+params=getViscoelasticFromSpheresInteraction(tc,en,es)
+facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,**params)) # **params sets kn, cn, ks, cs
+sphereMat=O.materials.append(ViscElMat(density=Density,frictionAngle=frictionAngle,**params))
+from yade import ymport
+fctIds=O.bodies.append(ymport.stl('baraban.stl',color=(1,0,0),material=facetMat))
+## Spheres
+sphereRadius = 0.2
+nbSpheres = (10,10,10)
+#nbSpheres = (50,50,50)
+for i in xrange(nbSpheres[0]):
+ for j in xrange(nbSpheres[1]):
+ for k in xrange(nbSpheres[2]):
+ x = (i*2 - nbSpheres[0])*sphereRadius*1.1
+ y = (j*2 - nbSpheres[1])*sphereRadius*1.1
+ z = (k*2 - nbSpheres[2])*sphereRadius*1.1
+ s=sphere([x,y,z],sphereRadius,material=sphereMat)
+ O.bodies.append(s)
+
+## Timestep
+O.dt=.2*tc
+
+## Engines
+O.engines=[
+ ## Resets forces and momenta the act on bodies
+ ForceResetter(),
+ ## Using bounding boxes find possible body collisions.
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
+ ## Interactions
+ InteractionLoop(
+ ## Create geometry information about each potential collision.
+ [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
+ ## Create physical information about the interaction.
+ [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
+ ## Constitutive law
+ [Law2_ScGeom_ViscElPhys_Basic()],
+ ),
+ ## Apply gravity
+ ## Cundall damping must been disabled!
+ NewtonIntegrator(damping=0,gravity=[0,-9.81,0]),
+## Saving results
+ #VTKRecorder(virtPeriod=0.04,fileName='/tmp/stlimp-',recorders=['spheres','facets']),
+ ## Apply kinematics to walls
+ BicyclePedalEngine(ids=fctIds,rotationAxis=[0,0,1],radius = 1.2,angularVelocity=4.0)
+]
+
+from yade import qt
+qt.View()
+#O.saveTmp()
+#O.run()
+
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2013-07-26 18:16:04 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2013-11-21 01:28:17 +0000
@@ -416,7 +416,7 @@
if (noCache) {cerr<<"Triangulation does not exist. Waht did you do?!"<<endl; return -1;}
RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
Cell_handle cell = Tri.locate(Point(X,Y,Z));
- return cell->info().index;
+ return cell->info().id;
}
template <class Tesselation>
@@ -485,7 +485,6 @@
void FlowBoundingSphere<Tesselation>::ComputeFacetForcesWithCache(bool onlyCache)
{
RTriangulation& Tri = T[currentTes].Triangulation();
- Finite_cells_iterator cell_end = Tri.finite_cells_end();
Vecteur nullVect(0,0,0);
//reset forces
if (!onlyCache) for (Finite_vertices_iterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) v->info().forces=nullVect;
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.ipp'
--- lib/triangulation/FlowBoundingSphereLinSolv.ipp 2013-06-27 11:29:17 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.ipp 2013-11-26 15:37:26 +0000
@@ -161,7 +161,8 @@
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
orderedCells.push_back(cell);
if (!cell->info().Pcondition) ++ncols;}
- spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort<RTriangulation>());
+// //Segfault on 14.10, and useless overall since SuiteSparse has preconditionners (including metis)
+// spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort<RTriangulation>());
T_cells.clear();
T_index=0;
=== modified file 'lib/triangulation/PeriodicFlow.cpp'
--- lib/triangulation/PeriodicFlow.cpp 2013-05-03 16:45:52 +0000
+++ lib/triangulation/PeriodicFlow.cpp 2013-11-21 01:28:17 +0000
@@ -47,7 +47,6 @@
void PeriodicFlow::ComputeFacetForcesWithCache(bool onlyCache)
{
RTriangulation& Tri = T[currentTes].Triangulation();
- Finite_cells_iterator cell_end = Tri.finite_cells_end();
Vecteur nullVect(0,0,0);
static vector<Vecteur> oldForces;
if (oldForces.size()<=Tri.number_of_vertices()) oldForces.resize(Tri.number_of_vertices()+1);
=== modified file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h 2013-08-22 14:32:01 +0000
+++ lib/triangulation/def_types.h 2013-11-21 01:28:17 +0000
@@ -44,6 +44,9 @@
class SimpleCellInfo : public Point {
public:
+ //"id": unique identifier of each cell, independant of other numberings used in the fluid types.
+ // Care to initialize it, there is no magic numbering to rely on
+ unsigned int id;
Real s;
bool isFictious;
SimpleCellInfo (void) {isFictious=false; s=0;}
=== modified file 'pkg/common/Grid.cpp'
--- pkg/common/Grid.cpp 2013-08-23 15:21:20 +0000
+++ pkg/common/Grid.cpp 2013-11-12 08:29:26 +0000
@@ -6,9 +6,6 @@
*************************************************************************/
#include "Grid.hpp"
-#ifdef YADE_OPENGL
- #include<yade/lib/opengl/OpenGLWrapper.hpp>
-#endif
//!################## SHAPES #####################
@@ -580,61 +577,3 @@
}
YADE_PLUGIN((Bo1_GridConnection_Aabb));
-
-#ifdef YADE_OPENGL
-//!################## Rendering #####################
-
-bool Gl1_GridConnection::wire;
-bool Gl1_GridConnection::glutNormalize;
-int Gl1_GridConnection::glutSlices;
-int Gl1_GridConnection::glutStacks;
-
-void Gl1_GridConnection::out( Quaternionr q )
-{
- AngleAxisr aa(q);
- std::cout << " axis: " << aa.axis()[0] << " " << aa.axis()[1] << " " << aa.axis()[2] << ", angle: " << aa.angle() << " | ";
-}
-
-void Gl1_GridConnection::go(const shared_ptr<Shape>& cm, const shared_ptr<State>& st ,bool wire2, const GLViewInfo&)
-{
- GridConnection *GC=static_cast<GridConnection*>(cm.get());
- Real r=GC->radius;
- Real length=GC->getLength();
- const shared_ptr<Interaction> intr = scene->interactions->find((int)GC->node1->getId(),(int)GC->node2->getId());
- Vector3r segt = GC->node2->state->pos - GC->node1->state->pos;
- if (scene->isPeriodic && intr) segt+=scene->cell->intrShiftPos(intr->cellDist);
- //glMaterialv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, Vector3f(cm->color[0],cm->color[1],cm->color[2]));
-
- glColor3v(cm->color);
- if(glutNormalize) glPushAttrib(GL_NORMALIZE);
-// glPushMatrix();
- Quaternionr shift;
- shift.setFromTwoVectors(Vector3r::UnitZ(),segt);
- if(intr){drawCylinder(wire || wire2, r,length,shift);}
-// if (intr && scene->isPeriodic) { glTranslatef(-segt[0],-segt[1],-segt[2]); drawCylinder(wire || wire2, r,length,-shift);}
- if(glutNormalize) glPopAttrib();
-// glPopMatrix();
- return;
-}
-
-void Gl1_GridConnection::drawCylinder(bool wire, Real radius, Real length, const Quaternionr& shift)
-{
- glPushMatrix();
- GLUquadricObj *quadObj = gluNewQuadric();
- gluQuadricDrawStyle(quadObj, (GLenum) (wire ? GLU_SILHOUETTE : GLU_FILL));
- gluQuadricNormals(quadObj, (GLenum) GLU_SMOOTH);
- gluQuadricOrientation(quadObj, (GLenum) GLU_OUTSIDE);
- AngleAxisr aa(shift);
- glRotatef(aa.angle()*180.0/Mathr::PI,aa.axis()[0],aa.axis()[1],aa.axis()[2]);
- gluCylinder(quadObj, radius, radius, length, glutSlices,glutStacks);
- gluQuadricOrientation(quadObj, (GLenum) GLU_INSIDE);
- //glutSolidSphere(radius,glutSlices,glutStacks);
- glTranslatef(0.0,0.0,length);
-
- //glutSolidSphere(radius,glutSlices,glutStacks);
-// gluDisk(quadObj,0.0,radius,glutSlices,_loops);
- gluDeleteQuadric(quadObj);
- glPopMatrix();
-}
-YADE_PLUGIN((Gl1_GridConnection));
-#endif
=== added file 'pkg/common/Grid_GUI.cpp'
--- pkg/common/Grid_GUI.cpp 1970-01-01 00:00:00 +0000
+++ pkg/common/Grid_GUI.cpp 2013-11-12 08:29:26 +0000
@@ -0,0 +1,67 @@
+/*************************************************************************
+* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
+* Copyright (C) 2012 by Bruno Chareyre bruno.chareyre@xxxxxxxxxxx *
+* 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 "Grid.hpp"
+#ifdef YADE_OPENGL
+ #include<yade/lib/opengl/OpenGLWrapper.hpp>
+
+//!################## Rendering #####################
+
+bool Gl1_GridConnection::wire;
+bool Gl1_GridConnection::glutNormalize;
+int Gl1_GridConnection::glutSlices;
+int Gl1_GridConnection::glutStacks;
+
+void Gl1_GridConnection::out( Quaternionr q )
+{
+ AngleAxisr aa(q);
+ std::cout << " axis: " << aa.axis()[0] << " " << aa.axis()[1] << " " << aa.axis()[2] << ", angle: " << aa.angle() << " | ";
+}
+
+void Gl1_GridConnection::go(const shared_ptr<Shape>& cm, const shared_ptr<State>& st ,bool wire2, const GLViewInfo&)
+{
+ GridConnection *GC=static_cast<GridConnection*>(cm.get());
+ Real r=GC->radius;
+ Real length=GC->getLength();
+ const shared_ptr<Interaction> intr = scene->interactions->find((int)GC->node1->getId(),(int)GC->node2->getId());
+ Vector3r segt = GC->node2->state->pos - GC->node1->state->pos;
+ if (scene->isPeriodic && intr) segt+=scene->cell->intrShiftPos(intr->cellDist);
+ //glMaterialv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, Vector3f(cm->color[0],cm->color[1],cm->color[2]));
+
+ glColor3v(cm->color);
+ if(glutNormalize) glPushAttrib(GL_NORMALIZE);
+// glPushMatrix();
+ Quaternionr shift;
+ shift.setFromTwoVectors(Vector3r::UnitZ(),segt);
+ if(intr){drawCylinder(wire || wire2, r,length,shift);}
+// if (intr && scene->isPeriodic) { glTranslatef(-segt[0],-segt[1],-segt[2]); drawCylinder(wire || wire2, r,length,-shift);}
+ if(glutNormalize) glPopAttrib();
+// glPopMatrix();
+ return;
+}
+
+void Gl1_GridConnection::drawCylinder(bool wire, Real radius, Real length, const Quaternionr& shift)
+{
+ glPushMatrix();
+ GLUquadricObj *quadObj = gluNewQuadric();
+ gluQuadricDrawStyle(quadObj, (GLenum) (wire ? GLU_SILHOUETTE : GLU_FILL));
+ gluQuadricNormals(quadObj, (GLenum) GLU_SMOOTH);
+ gluQuadricOrientation(quadObj, (GLenum) GLU_OUTSIDE);
+ AngleAxisr aa(shift);
+ glRotatef(aa.angle()*180.0/Mathr::PI,aa.axis()[0],aa.axis()[1],aa.axis()[2]);
+ gluCylinder(quadObj, radius, radius, length, glutSlices,glutStacks);
+ gluQuadricOrientation(quadObj, (GLenum) GLU_INSIDE);
+ //glutSolidSphere(radius,glutSlices,glutStacks);
+ glTranslatef(0.0,0.0,length);
+
+ //glutSolidSphere(radius,glutSlices,glutStacks);
+// gluDisk(quadObj,0.0,radius,glutSlices,_loops);
+ gluDeleteQuadric(quadObj);
+ glPopMatrix();
+}
+YADE_PLUGIN((Gl1_GridConnection));
+#endif
=== modified file 'pkg/common/KinematicEngines.cpp'
--- pkg/common/KinematicEngines.cpp 2013-10-16 13:48:29 +0000
+++ pkg/common/KinematicEngines.cpp 2013-11-11 16:22:15 +0000
@@ -4,7 +4,7 @@
#include<yade/pkg/dem/Shop.hpp>
#include<yade/lib/smoothing/LinearInterpolate.hpp>
-YADE_PLUGIN((KinematicEngine)(CombinedKinematicEngine)(TranslationEngine)(HarmonicMotionEngine)(RotationEngine)(HelixEngine)(InterpolatingHelixEngine)(HarmonicRotationEngine)(ServoPIDController));
+YADE_PLUGIN((KinematicEngine)(CombinedKinematicEngine)(TranslationEngine)(HarmonicMotionEngine)(RotationEngine)(HelixEngine)(InterpolatingHelixEngine)(HarmonicRotationEngine)(ServoPIDController)(BicyclePedalEngine));
CREATE_LOGGER(KinematicEngine);
@@ -109,7 +109,6 @@
}
}
-
void RotationEngine::apply(const vector<Body::id_t>& ids){
if (ids.size()>0) {
#ifdef YADE_OPENMP
@@ -183,3 +182,33 @@
TranslationEngine::apply(ids);
}
+
+void BicyclePedalEngine::apply(const vector<Body::id_t>& ids){
+ if (ids.size()>0) {
+ Quaternionr qRotateZVec(Quaternionr().setFromTwoVectors(Vector3r(0,0,1), rotationAxis));
+
+ Vector3r newPos = Vector3r(cos(fi + angularVelocity*scene->dt)*radius, sin(fi + angularVelocity*scene->dt)*radius, 0.0);
+ Vector3r oldPos = Vector3r(cos(fi)*radius, sin(fi)*radius, 0.0);
+
+ Vector3r newVel = (oldPos - newPos)/scene->dt;
+
+ fi += angularVelocity*scene->dt;
+ newVel = qRotateZVec*newVel;
+
+ #ifdef YADE_OPENMP
+ const long size=ids.size();
+ #pragma omp parallel for schedule(static)
+ for(long i=0; i<size; i++){
+ const Body::id_t& id=ids[i];
+ #else
+ FOREACH(Body::id_t id,ids){
+ #endif
+ assert(id<(Body::id_t)scene->bodies->size());
+ Body* b=Body::byId(id,scene).get();
+ if(!b) continue;
+ b->state->vel+=newVel;
+ }
+ } else {
+ LOG_WARN("The list of ids is empty! Can't move any body.");
+ }
+}
=== modified file 'pkg/common/KinematicEngines.hpp'
--- pkg/common/KinematicEngines.hpp 2013-07-04 17:36:34 +0000
+++ pkg/common/KinematicEngines.hpp 2013-11-11 16:22:15 +0000
@@ -115,3 +115,15 @@
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(ServoPIDController);
+
+struct BicyclePedalEngine: public KinematicEngine{
+ virtual void apply(const vector<Body::id_t>& ids);
+ void postLoad(BicyclePedalEngine&){ rotationAxis.normalize(); }
+ YADE_CLASS_BASE_DOC_ATTRS(BicyclePedalEngine,KinematicEngine,"Engine applying the linear motion of ``bicycle pedal`` e.g. moving points around the axis without rotation",
+ ((Real,angularVelocity,0,,"Angular velocity. [rad/s]"))
+ ((Vector3r,rotationAxis,Vector3r::UnitX(),Attr::triggerPostLoad,"Axis of rotation (direction); will be normalized automatically."))
+ ((Real,radius,-1.0,,"Rotation radius. [m]"))
+ ((Real,fi,Mathr::PI/2.0,,"Initial phase [radians]"))
+ );
+};
+REGISTER_SERIALIZABLE(BicyclePedalEngine);
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2013-10-25 07:27:55 +0000
+++ pkg/dem/FlowEngine.cpp 2013-11-25 02:04:22 +0000
@@ -62,15 +62,9 @@
timingDeltas->checkpoint ( "Update_Volumes" );
Eps_Vol_Cumulative += eps_vol_max;
- if ( (defTolerance>0 && Eps_Vol_Cumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval) {
- updateTriangulation = true;
- Eps_Vol_Cumulative=0;
- retriangulationLastIter=0;
- ReTrg++;
- } else {
- updateTriangulation = false;
- retriangulationLastIter++;}
-
+ retriangulationLastIter++;
+ if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
+ (defTolerance>0 && Eps_Vol_Cumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval;
///Compute flow and and forces here
if (pressureForce){
@@ -117,6 +111,7 @@
backgroundCompleted=false;
retriangulationLastIter=ellapsedIter;
updateTriangulation=false;
+ Eps_Vol_Cumulative=0;
ellapsedIter=0;
boost::thread workerThread(&FlowEngine::backgroundAction,this);
workerThread.detach();
@@ -134,7 +129,10 @@
Build_Triangulation (P_zero, solver);
Initialize_volumes(solver);
ComputeViscousForces(*solver);
- updateTriangulation = false;}
+ updateTriangulation = false;
+ Eps_Vol_Cumulative=0;
+ retriangulationLastIter=0;
+ ReTrg++;}
}
first=false;
timingDeltas->checkpoint ( "Triangulate + init volumes" );
@@ -265,8 +263,10 @@
flow->T[flow->currentTes].cellHandles.clear();
flow->T[flow->currentTes].cellHandles.reserve(flow->T[flow->currentTes].Triangulation().number_of_finite_cells());
Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
- for ( Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ )
+ int k=0;
+ for ( Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
flow->T[flow->currentTes].cellHandles.push_back(cell);
+ cell->info().id=k++;}//define unique numbering now, corresponds to position in cellHandles
flow->DisplayStatistics ();
flow->Compute_Permeability();
porosity = flow->V_porale_porosity/flow->V_totale_porosity;
@@ -730,14 +730,10 @@
timingDeltas->checkpoint("Triangulating");
UpdateVolumes (solver);
Eps_Vol_Cumulative += eps_vol_max;
- if ( (defTolerance>0 && Eps_Vol_Cumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval ) {
- updateTriangulation = true;
- Eps_Vol_Cumulative=0;
- retriangulationLastIter=0;
- ReTrg++;
- } else {
- updateTriangulation = false;
- retriangulationLastIter++;}
+ retriangulationLastIter++;
+ if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
+ (defTolerance>0 && Eps_Vol_Cumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval;
+
timingDeltas->checkpoint("Update_Volumes");
///Compute flow and and forces here
@@ -782,6 +778,7 @@
backgroundCompleted=false;
retriangulationLastIter=ellapsedIter;
ellapsedIter=0;
+ Eps_Vol_Cumulative=0;
boost::thread workerThread(&PeriodicFlowEngine::backgroundAction,this);
workerThread.detach();
Initialize_volumes(solver);
@@ -796,7 +793,10 @@
Build_Triangulation (P_zero, solver);
Initialize_volumes(solver);
ComputeViscousForces(*solver);
- updateTriangulation = false;}
+ updateTriangulation = false;
+ Eps_Vol_Cumulative=0;
+ retriangulationLastIter=0;
+ ReTrg++;}
}
first=false;
timingDeltas->checkpoint("Ending");
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2013-10-25 07:27:55 +0000
+++ pkg/dem/FlowEngine.hpp 2013-11-25 00:52:12 +0000
@@ -146,6 +146,13 @@
void pressureProfile(double wallUpY, double wallDownY) {return solver->measurePressureProfile(wallUpY,wallDownY);}
double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
TPL int getCell(double posX, double posY, double posZ, Solver& flow){return flow->getCell(posX, posY, posZ);}
+ TPL unsigned int nCells(Solver& flow){return flow->T[flow->currentTes].cellHandles.size();}
+ TPL python::list getVertices(unsigned int id, Solver& flow){
+ python::list ids;
+ if (id>=flow->T[flow->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<flow->T[flow->currentTes].cellHandles.size()); return ids;}
+ for (unsigned int i=0;i<4;i++) ids.append(flow->T[flow->currentTes].cellHandles[id]->vertex(i)->info().id());
+ return ids;
+ }
double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
double averagePressure(){return solver->averagePressure();}
#ifdef LINSOLV
@@ -181,6 +188,8 @@
Real _getCellFlux(unsigned int cond) {return getCellFlux(cond,solver);}
Real _getBoundaryFlux(unsigned int boundary) {return getBoundaryFlux(boundary,solver);}
int _getCell(Vector3r pos) {return getCell(pos[0],pos[1],pos[2],solver);}
+ unsigned int _nCells() {return nCells(solver);}
+ python::list _getVertices(unsigned int id) {return getVertices(id,solver);}
#ifdef LINSOLV
void _exportMatrix(string filename) {exportMatrix(filename,solver);}
void _exportTriplets(string filename) {exportTriplets(filename,solver);}
@@ -267,7 +276,6 @@
normal[wall_ymax].y()=normal[wall_xmax].x()=normal[wall_zmax].z()=-1;
solver = shared_ptr<FlowSolver> (new FlowSolver);
first=true;
- updateTriangulation=false;
eps_vol_max=Eps_Vol_Cumulative=retriangulationLastIter=0;
ReTrg=1;
backgroundCompleted=true;
@@ -304,6 +312,8 @@
.def("updateBCs",&FlowEngine::_updateBCs,"tells the engine to update it's boundary conditions before running (especially useful when changing boundary pressure - should not be needed for point-wise imposed pressure)")
.def("emulateAction",&FlowEngine::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
.def("getCell",&FlowEngine::_getCell,(python::arg("pos")),"get id of the cell containing (X,Y,Z).")
+ .def("nCells",&FlowEngine::_nCells,"get the total number of finite cells in the triangulation.")
+ .def("getVertices",&FlowEngine::_getVertices,(python::arg("id")),"get the vertices of a cell")
#ifdef LINSOLV
.def("exportMatrix",&FlowEngine::_exportMatrix,(python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
.def("exportTriplets",&FlowEngine::_exportTriplets,(python::arg("filename")="triplets"),"Export system matrix to a file with only non-zero entries.")
@@ -426,7 +436,6 @@
wallIds=vector<int>(6,-1);
// wallTopId=wallBottomId=wallFrontId=wallBackId=wallLeftId=wallRightId=-1;
solver = shared_ptr<FlowSolver> (new FlowSolver);
- updateTriangulation=false;
eps_vol_max=Eps_Vol_Cumulative=retriangulationLastIter=0;
ReTrg=1;
first=true;
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2013-11-28 09:45:10 +0000
+++ py/_utils.cpp 2013-11-28 09:45:15 +0000
@@ -355,35 +355,37 @@
upper part" (lower and upper parts with respect to the plane's normal).
(This could be easily extended to return sum of only normal forces or only of shear forces.)
- // FIXME: adapt to ScGeom or remove
*/
-// Vector3r forcesOnPlane(const Vector3r& planePt, const Vector3r& normal){
-// Vector3r ret(Vector3r::Zero());
-// Scene* scene=Omega::instance().getScene().get();
-// FOREACH(const shared_ptr<Interaction>&I, *scene->interactions){
-// if(!I->isReal()) continue;
-// NormShearPhys* nsi=dynamic_cast<NormShearPhys*>(I->phys.get());
-// if(!nsi) continue;
-// Vector3r pos1,pos2;
-// Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(I->geom.get()); // Dem3DofGeom has copy of se3 in itself, otherwise we have to look up the bodies
-// if(d3dg){ pos1=d3dg->se31.position; pos2=d3dg->se32.position; }
-// else{ pos1=Body::byId(I->getId1(),scene)->state->pos; pos2=Body::byId(I->getId2(),scene)->state->pos; }
-// Real dot1=(pos1-planePt).dot(normal), dot2=(pos2-planePt).dot(normal);
-// if(dot1*dot2>0) continue; // both (centers of) bodies are on the same side of the plane=> this interaction has to be disregarded
-// // if pt1 is on the negative plane side, d3dg->normal.Dot(normal)>0, the force is well oriented;
-// // otherwise, reverse its contribution. So that we return finally
-// // Sum [ ( normal(plane) dot normal(interaction= from 1 to 2) ) "nsi->force" ]
-// ret+=(dot1<0.?1.:-1.)*(nsi->normalForce+nsi->shearForce);
-// }
-// return ret;
-// }
+Vector3r forcesOnPlane(const Vector3r& planePt, const Vector3r& normal){
+ Vector3r ret(Vector3r::Zero());
+ Scene* scene=Omega::instance().getScene().get();
+ FOREACH(const shared_ptr<Interaction>&I, *scene->interactions){
+ if(!I->isReal()) continue;
+ NormShearPhys* nsi=dynamic_cast<NormShearPhys*>(I->phys.get());
+ if(!nsi) continue;
+ Vector3r pos1,pos2;
+ /*
+ Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(I->geom.get()); // Dem3DofGeom has copy of se3 in itself, otherwise we have to look up the bodies
+ if(d3dg){ pos1=d3dg->se31.position; pos2=d3dg->se32.position; }
+ else{ pos1=Body::byId(I->getId1(),scene)->state->pos; pos2=Body::byId(I->getId2(),scene)->state->pos; }
+ */
+ pos1=Body::byId(I->getId1(),scene)->state->pos; pos2=Body::byId(I->getId2(),scene)->state->pos;
+ Real dot1=(pos1-planePt).dot(normal), dot2=(pos2-planePt).dot(normal);
+ if(dot1*dot2>0) continue; // both (centers of) bodies are on the same side of the plane=> this interaction has to be disregarded
+ // if pt1 is on the negative plane side, d3dg->normal.Dot(normal)>0, the force is well oriented;
+ // otherwise, reverse its contribution. So that we return finally
+ // Sum [ ( normal(plane) dot normal(interaction= from 1 to 2) ) "nsi->force" ]
+ ret+=(dot1<0.?1.:-1.)*(nsi->normalForce+nsi->shearForce);
+ }
+ return ret;
+}
/* Less general than forcesOnPlane, computes force on plane perpendicular to axis, passing through coordinate coord. */
-// /*Vector3r forcesOnCoordPlane(Real coord, int axis){
-// Vector3r planePt(Vector3r::Zero()); planePt[axis]=coord;
-// Vector3r normal(Vector3r::Zero()); normal[axis]=1;
-// return forcesOnPlane(planePt,normal);
-// }*/
+Vector3r forcesOnCoordPlane(Real coord, int axis){
+ Vector3r planePt(Vector3r::Zero()); planePt[axis]=coord;
+ Vector3r normal(Vector3r::Zero()); normal[axis]=1;
+ return forcesOnPlane(planePt,normal);
+}
py::tuple spiralProject(const Vector3r& pt, Real dH_dTheta, int axis=2, Real periodStart=std::numeric_limits<Real>::quiet_NaN(), Real theta0=0){
@@ -517,8 +519,8 @@
py::def("sumForces",sumForces,(py::arg("ids"),py::arg("direction")),"Return summary force on bodies with given *ids*, projected on the *direction* vector.");
py::def("sumTorques",sumTorques,(py::arg("ids"),py::arg("axis"),py::arg("axisPt")),"Sum forces and torques on bodies given in *ids* with respect to axis specified by a point *axisPt* and its direction *axis*.");
py::def("sumFacetNormalForces",sumFacetNormalForces,(py::arg("ids"),py::arg("axis")=-1),"Sum force magnitudes on given bodies (must have :yref:`shape<Body.shape>` of the :yref:`Facet` type), considering only part of forces perpendicular to each :yref:`facet's<Facet>` face; if *axis* has positive value, then the specified axis (0=x, 1=y, 2=z) will be used instead of facet's normals.");
-// py::def("forcesOnPlane",forcesOnPlane,(py::arg("planePt"),py::arg("normal")),"Find all interactions deriving from :yref:`NormShearPhys` that cross given plane and sum forces (both normal and shear) on them.\n\n:param Vector3 planePt: a point on the plane\n:param Vector3 normal: plane normal (will be normalized).\n");
-// py::def("forcesOnCoordPlane",forcesOnCoordPlane);
+ py::def("forcesOnPlane",forcesOnPlane,(py::arg("planePt"),py::arg("normal")),"Find all interactions deriving from :yref:`NormShearPhys` that cross given plane and sum forces (both normal and shear) on them.\n\n:param Vector3 planePt: a point on the plane\n:param Vector3 normal: plane normal (will be normalized).\n");
+ py::def("forcesOnCoordPlane",forcesOnCoordPlane);
py::def("totalForceInVolume",Shop__totalForceInVolume,"Return summed forces on all interactions and average isotropic stiffness, as tuple (Vector3,float)");
py::def("createInteraction",Shop__createExplicitInteraction,(py::arg("id1"),py::arg("id2")),"Create interaction between given bodies by hand.\n\nCurrent engines are searched for :yref:`IGeomDispatcher` and :yref:`IPhysDispatcher` (might be both hidden in :yref:`InteractionLoop`). Geometry is created using ``force`` parameter of the :yref:`geometry dispatcher<IGeomDispatcher>`, wherefore the interaction will exist even if bodies do not spatially overlap and the functor would return ``false`` under normal circumstances. \n\n.. warning:: This function will very likely behave incorrectly for periodic simulations (though it could be extended it to handle it farily easily).");
py::def("spiralProject",spiralProject,(py::arg("pt"),py::arg("dH_dTheta"),py::arg("axis")=2,py::arg("periodStart")=std::numeric_limits<Real>::quiet_NaN(),py::arg("theta0")=0));
=== modified file 'scripts/checks-and-tests/checks/DEM-PFV-check.py'
--- scripts/checks-and-tests/checks/DEM-PFV-check.py 2013-08-30 10:33:17 +0000
+++ scripts/checks-and-tests/checks/DEM-PFV-check.py 2013-11-26 16:22:14 +0000
@@ -79,7 +79,7 @@
e22=e22-triax.strain[1]
modulus = 1000./abs(e22)
- target=249064.586653
+ target=252759.905803
if abs((modulus-target)/target)>tolerance :
print "DEM-PFV: difference in bulk modulus:", modulus, "vs. target ",target
errors+=1
@@ -105,7 +105,7 @@
print "DEM-PFV: unbalanced Qin vs. Qout"
errors+=1
- target=0.0408678245942
+ target=0.040399916554
if abs((permeability-target)/target)>tolerance :
print "DEM-PFV: difference in permeability:",permeability," vs. target ",target
errors+=1
@@ -123,18 +123,17 @@
from yade import timing
O.run(3000,1)
- target=637.268936033
+ target=628.314160434
if abs((flow.getPorePressure((0.5,0.1,0.5))-target)/target)>tolerance :
print "DEM-PFV: difference in final pressure:",flow.getPorePressure((0.5,0.1,0.5))," vs. target ",target
errors+=1
- target=0.00260892345196
+ target=0.00258113045083
if abs((triax.strain[1]-zeroe22-target)/target)>tolerance :
print "DEM-PFV: difference in final deformation",triax.strain[1]-zeroe22," vs. target ",target
errors+=1
if (float(flow.execTime)/float(sum([e.execTime for e in O.engines])))>0.6 :
- print "DEM-PFV: More than 60\% of cpu time in FlowEngine (",100.*(float(flow.execTime)/float(sum([e.execTime for e in O.engines]))) ,"%). Should not happen with efficient libraries (check blas/lapack/cholmod implementations)"
- errors+=1
+ print "(INFO) DEM-PFV: More than 60\% of cpu time in FlowEngine (",100.*(float(flow.execTime)/float(sum([e.execTime for e in O.engines]))) ,"%). Should not happen with efficient libraries (check blas/lapack/cholmod implementations)"
flow.forceMetis=True
O.run(201,1)
=== added directory 'scripts/ppa'
=== added file 'scripts/ppa/createtar.py'
--- scripts/ppa/createtar.py 1970-01-01 00:00:00 +0000
+++ scripts/ppa/createtar.py 2013-11-20 19:50:17 +0000
@@ -0,0 +1,98 @@
+#!/usr/bin/env python
+import argparse, os, git, shutil
+
+
+jobsNumber = 6
+
+parser = argparse.ArgumentParser(description='Process some integers.')
+parser.add_argument("g", help="git debian-directory")
+parser.add_argument("p", help="direcrory, where will be created pbuilder tarballs")
+parser.add_argument("u", help="git upstream-direcrory")
+parser.add_argument("-u","--update", help="update tarballs, if they are exist", action='store_true')
+args = parser.parse_args()
+
+gitdebdir = args.g
+gitupsdir = args.u
+pbdir = args.p
+
+if not(os.path.isdir(gitdebdir)):
+ raise RuntimeError('Git-debian-directory does not exists')
+
+if not(os.path.isdir(gitupsdir)):
+ raise RuntimeError('Git-upstream-directory does not exists')
+
+if not(os.path.isdir(pbdir)):
+ print ('Pbuilder-directory does not exists. Creating')
+ os.mkdir(pbdir)
+
+repodeb = git.Repo(gitdebdir)
+repoups = git.Repo(gitupsdir)
+assert repodeb.bare == False
+assert repoups.bare == False
+
+if (repodeb.is_dirty()):
+ raise RuntimeError('Git-debian-repo has an uncommitted changes. Exiting.')
+
+if (repoups.is_dirty()):
+ raise RuntimeError('Git-upstream-repo has an uncommitted changes. Exiting.')
+
+for branch in repodeb.branches:
+ branchstr = str(branch)
+ if (branchstr<>'master'):
+ print "Switching to branch %s"%(branch)
+ repodeb.git.checkout(branch)
+ infile = open(gitdebdir+"/pbuilder")
+ lines = infile.readlines()
+ mirror = lines[0].strip()
+ components = lines[1].strip()
+ archs = lines[2].split()
+ keyringuse = lines[3].strip()
+ othermirror = lines[4].strip()
+ infile.close()
+ for a in archs:
+ tarball = "%s/%s_%s.tgz"%(pbdir, branchstr, a.strip())
+ addAllowuntrusted = ""
+ if (othermirror!="#"):
+ addAllowuntrusted = " --allow-untrusted "
+ if not(os.path.isfile(tarball)):
+ createPbTar = ('sudo pbuilder --create --distribution %s --mirror %s --components "%s" --architecture %s --debootstrapopts "--keyring=%s" --basetgz %s'%
+ (branchstr, mirror, components, a, keyringuse, tarball))
+ if (othermirror!="#"):
+ createPbTar += ' --othermirror "' + othermirror + '"'
+ addAllowuntrusted = " --allow-untrusted "
+ print createPbTar
+
+ print "Creating tarball %s"%(tarball)
+ os.system(createPbTar)
+ else:
+ print "Tarball %s exists"%(tarball)
+ if (args.update):
+ print "Updating %s as requested" %(tarball)
+ updatePbTar = ('sudo pbuilder --update --basetgz %s'%(tarball))
+ os.system(updatePbTar)
+
+ # Creating dir for building
+ builddirup="%s_%s/"%(branchstr, a)
+ builddirdeb="%s/build/"%(builddirup)
+ builddirres="%s/result/"%(builddirup)
+ shutil.rmtree(builddirdeb,ignore_errors=True)
+ shutil.copytree(gitupsdir, builddirdeb )
+ shutil.rmtree(builddirdeb+".git")
+ # Get package version
+ versiondebian = repoups.git.describe() + "~" + branchstr
+
+ # Get package name
+ infilepname = open(gitdebdir+"/changelog"); sourcePackName = infilepname.readlines()[0].split()[0]
+ print sourcePackName
+
+ os.system('cd %s; apack %s_%s.orig.tar.xz build'%(builddirup,sourcePackName,versiondebian))
+ shutil.copytree(gitdebdir, builddirdeb+"/debian")
+
+ os.system('sed -i.bak -e s/VERSION/%s/ -e s/DISTRIBUTION/%s/ %s/debian/changelog'%(versiondebian,branch,builddirdeb))
+ os.system('cd %s; dpkg-source -b -I build'%(builddirup))
+ os.mkdir(builddirres)
+ print "Building package %s_%s"%(sourcePackName, versiondebian)
+ buildPackage = ('sudo pbuilder --build --architecture %s --basetgz %s %s --logfile %s/pbuilder.log --debbuildopts "-j%d" --buildresult %s %s/*.dsc'%
+ (a, tarball, addAllowuntrusted, builddirup, jobsNumber, builddirres, builddirup))
+ print buildPackage
+ os.system(buildPackage)
Follow ups
