yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02900
[Branch ~yade-dev/yade/trunk] Rev 1935: 1. Preliminary attempt at https://blueprints.launchpad.net/yade/+spec/c++-numpy-arrays (for postp...
------------------------------------------------------------
revno: 1935
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Fri 2010-01-01 16:21:30 +0100
message:
1. Preliminary attempt at https://blueprints.launchpad.net/yade/+spec/c++-numpy-arrays (for postprocessing and similar), example in eudoxos.testNumpy
2. SpherePack::makeCloud can be asked for negative number of spheres (= as much as possible)
added:
lib/pyutil/
lib/pyutil/README
lib/pyutil/gil.hpp
lib/pyutil/numpy.hpp
lib/pyutil/numpy_boost.hpp
modified:
pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp
pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp
pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp
pkg/common/Engine/GlobalEngine/PeriodicPythonRunner.hpp
pkg/dem/DataClass/SpherePack.cpp
pkg/dem/DataClass/SpherePack.hpp
py/_eudoxos.cpp
py/yadeWrapper/customConverters.cpp
scripts/test/periodic-triax-velgrad.py
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.
=== added directory 'lib/pyutil'
=== added file 'lib/pyutil/README'
--- lib/pyutil/README 1970-01-01 00:00:00 +0000
+++ lib/pyutil/README 2010-01-01 15:21:30 +0000
@@ -0,0 +1,2 @@
+numpy_boost.cpp:
+ http://code.google.com/p/numpy-boost/source/checkout (svn rev 2)
=== added file 'lib/pyutil/gil.hpp'
--- lib/pyutil/gil.hpp 1970-01-01 00:00:00 +0000
+++ lib/pyutil/gil.hpp 2010-01-01 15:21:30 +0000
@@ -0,0 +1,15 @@
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+
+#include<Python.h>
+//! class (scoped lock) managing python's Global Interpreter Lock (gil)
+class gilLock{
+ PyGILState_STATE state;
+ public:
+ gilLock(){ state=PyGILState_Ensure(); }
+ ~gilLock(){ PyGILState_Release(state); }
+};
+
+//! run string as python command; locks & unlocks GIL automatically
+void pyRunString(const std::string& cmd){
+ gilLock lock; PyRun_SimpleString(cmd.c_str());
+};
=== added file 'lib/pyutil/numpy.hpp'
--- lib/pyutil/numpy.hpp 1970-01-01 00:00:00 +0000
+++ lib/pyutil/numpy.hpp 2010-01-01 15:21:30 +0000
@@ -0,0 +1,6 @@
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx
+
+#include"numpy_boost.hpp"
+
+// helper macro do assign Vector3r values to a subarray
+#define VECTOR3R_TO_NUMPY(vec,arr) arr[0]=vec[0]; arr[1]=vec[1]; arr[2]=vec[2]
=== added file 'lib/pyutil/numpy_boost.hpp'
--- lib/pyutil/numpy_boost.hpp 1970-01-01 00:00:00 +0000
+++ lib/pyutil/numpy_boost.hpp 2010-01-01 15:21:30 +0000
@@ -0,0 +1,199 @@
+/*
+Copyright (c) 2008, Michael Droettboom
+All rights reserved.
+
+Licensed under the BSD license.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above
+ copyright notice, this list of conditions and the following
+ disclaimer in the documentation and/or other materials provided
+ with the distribution.
+
+ * The names of its contributors may not be used to endorse or
+ promote products derived from this software without specific
+ prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/*
+ Documentation to come
+*/
+
+#ifndef __NUMPY_BOOST_HPP__
+#define __NUMPY_BOOST_HPP__
+
+#include <Python.h>
+#include <numpy/arrayobject.h>
+#include <boost/multi_array.hpp>
+#include <boost/cstdint.hpp>
+#include <complex>
+#include <algorithm>
+
+namespace detail {
+ template<class T>
+ class numpy_type_map {
+ public:
+ static const int typenum;
+ };
+
+ template<>
+ const int numpy_type_map<float>::typenum = NPY_FLOAT;
+
+ template<>
+ const int numpy_type_map<std::complex<float> >::typenum = NPY_CFLOAT;
+
+ template<>
+ const int numpy_type_map<double>::typenum = NPY_DOUBLE;
+
+ template<>
+ const int numpy_type_map<std::complex<double> >::typenum = NPY_CDOUBLE;
+
+ template<>
+ const int numpy_type_map<long double>::typenum = NPY_LONGDOUBLE;
+
+ template<>
+ const int numpy_type_map<std::complex<long double> >::typenum = NPY_CLONGDOUBLE;
+
+ template<>
+ const int numpy_type_map<boost::int8_t>::typenum = NPY_INT8;
+
+ template<>
+ const int numpy_type_map<boost::uint8_t>::typenum = NPY_UINT8;
+
+ template<>
+ const int numpy_type_map<boost::int16_t>::typenum = NPY_INT16;
+
+ template<>
+ const int numpy_type_map<boost::uint16_t>::typenum = NPY_UINT16;
+
+ template<>
+ const int numpy_type_map<boost::int32_t>::typenum = NPY_INT32;
+
+ template<>
+ const int numpy_type_map<boost::uint32_t>::typenum = NPY_UINT32;
+
+ template<>
+ const int numpy_type_map<boost::int64_t>::typenum = NPY_INT64;
+
+ template<>
+ const int numpy_type_map<boost::uint64_t>::typenum = NPY_INT64;
+}
+
+class numpy_boost_exception : public std::exception {
+
+};
+
+template<class T, int NDims>
+class numpy_boost : public boost::multi_array_ref<T, NDims>
+{
+public:
+ typedef numpy_boost<T, NDims> self_type;
+ typedef boost::multi_array_ref<T, NDims> super;
+ typedef typename super::size_type size_type;
+ typedef T* TPtr;
+
+private:
+ PyArrayObject* array;
+
+ void init_from_array(PyArrayObject* a) {
+ array = a;
+ super::base_ = (TPtr)PyArray_DATA(a);
+
+ // It would seem like we would want to choose C or Fortran
+ // ordering here based on the flags in the Numpy array. However,
+ // those flags are purely informational, the actually information
+ // about storage order is recorded in the strides.
+ super::storage_ = boost::c_storage_order();
+
+ boost::detail::multi_array::copy_n(PyArray_DIMS(a), NDims, super::extent_list_.begin());
+ for (size_t i = 0; i < NDims; ++i) {
+ super::stride_list_[i] = PyArray_STRIDE(a, i) / sizeof(T);
+ }
+ std::fill_n(super::index_base_list_.begin(), NDims, 0);
+ super::origin_offset_ = 0;
+ super::directional_offset_ = 0;
+ super::num_elements_ = std::accumulate(super::extent_list_.begin(),
+ super::extent_list_.end(),
+ size_type(1),
+ std::multiplies<size_type>());
+ }
+
+public:
+ numpy_boost(PyObject* obj) :
+ super(NULL, std::vector<boost::uint32_t>(NDims, 0)),
+ array(NULL)
+ {
+ PyArrayObject* a;
+
+ a = (PyArrayObject*)PyArray_FromObject(obj, detail::numpy_type_map<T>::typenum, NDims, NDims);
+ if (a == NULL) {
+ // TODO: Extract Python exception
+ throw numpy_boost_exception();
+ }
+
+ init_from_array(a);
+ }
+
+ numpy_boost(const self_type &other) :
+ super(NULL, std::vector<boost::uint32_t>(NDims, 0)),
+ array(NULL)
+ {
+ Py_INCREF(other.array);
+ init_from_array(other.array);
+ }
+
+ template<class ExtentsList>
+ explicit numpy_boost(const ExtentsList& extents) :
+ super(NULL, std::vector<boost::uint32_t>(NDims, 0)),
+ array(NULL)
+ {
+ npy_intp shape[NDims];
+ PyArrayObject* a;
+
+ boost::detail::multi_array::copy_n(extents, NDims, shape);
+
+ a = (PyArrayObject*)PyArray_SimpleNew(NDims, shape, detail::numpy_type_map<T>::typenum);
+ if (a == NULL) {
+ // TODO: Extract Python exception
+ throw numpy_boost_exception();
+ }
+
+ init_from_array(a);
+ }
+
+ ~numpy_boost() {
+ Py_DECREF(array);
+ }
+
+ void operator=(const self_type &other) {
+ Py_DECREF(array);
+ Py_INCREF(other.array);
+ init_from_array(other.array);
+ }
+
+ PyObject*
+ py_ptr() {
+ Py_INCREF(array);
+ return (PyObject*)array;
+ }
+};
+
+#endif
=== modified file 'pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp'
--- pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2009-12-25 14:46:48 +0000
+++ pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2010-01-01 15:21:30 +0000
@@ -36,7 +36,8 @@
geomDispatcher->updateScenePtr();
physDispatcher->updateScenePtr();
lawDispatcher->updateScenePtr();
- Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cell->getSize();
+ // transformed cell size
+ Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cell->trsf*scene->cell->refSize;
bool removeUnseenIntrs=(scene->interactions->iterColliderLastRun>=0 && scene->interactions->iterColliderLastRun==scene->currentIteration);
#ifdef YADE_OPENMP
const long size=scene->interactions->size();
@@ -89,8 +90,6 @@
if(!scene->isPeriodic) geomCreated=I->functorCache.geom->go(b1->shape,b2->shape, *b1->state, *b2->state, Vector3r::ZERO, /*force*/false, I);
else{ // handle periodicity
Vector3r shift2(I->cellDist[0]*cellSize[0],I->cellDist[1]*cellSize[1],I->cellDist[2]*cellSize[2]);
- // in sheared cell, apply shear on the mutual position as well
- shift2=scene->cell->shearPt(shift2);
geomCreated=I->functorCache.geom->go(b1->shape,b2->shape,*b1->state,*b2->state,shift2,/*force*/false,I);
}
if(!geomCreated){
=== modified file 'pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp'
--- pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp 2009-12-20 22:03:17 +0000
+++ pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp 2010-01-01 15:21:30 +0000
@@ -48,7 +48,7 @@
updateScenePtr();
shared_ptr<BodyContainer>& bodies = scene->bodies;
- Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cell->getSize();
+ Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cell->trsf*scene->cell->refSize;
bool removeUnseenIntrs=(scene->interactions->iterColliderLastRun>=0 && scene->interactions->iterColliderLastRun==scene->currentIteration);
#ifdef YADE_OPENMP
const long size=scene->interactions->size();
@@ -74,7 +74,6 @@
geomCreated=operator()(b1->shape, b2->shape, *b1->state, *b2->state, Vector3r::ZERO, /*force*/ false, I);
} else{
Vector3r shift2(I->cellDist[0]*cellSize[0],I->cellDist[1]*cellSize[1],I->cellDist[2]*cellSize[2]); // add periodicity to the position of the 2nd body
- shift2=scene->cell->shearPt(shift2);
geomCreated=operator()(b1->shape, b2->shape, *b1->state, *b2->state, shift2, /*force*/ false, I);
}
// reset && erase interaction that existed but now has no geometry anymore
=== modified file 'pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp'
--- pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp 2009-12-30 20:10:59 +0000
+++ pkg/common/Engine/GlobalEngine/InsertionSortCollider.cpp 2010-01-01 15:21:30 +0000
@@ -426,9 +426,6 @@
TRVAR4(pmn1,pmx1,pmn2,pmx2);
}
#endif
- /* FIXME: condition temporary disabled
- if wMn==m1 and minima[id2]<m1, then it might give spurious error; why should it be error, though?
- */
if(((pmn1!=pmx1) || (pmn2!=pmx2))){
Real span=(pmn1!=pmx1?mx1-mn1:mx2-mn2); if(span<0) span=dim-span;
LOG_FATAL("Body #"<<(pmn1!=pmx1?id1:id2)<<" spans over half of the cell size "<<dim<<" (axis="<<axis<<", min="<<(pmn1!=pmx1?mn1:mn2)<<", max="<<(pmn1!=pmx1?mx1:mx2)<<", span="<<span<<")");
=== modified file 'pkg/common/Engine/GlobalEngine/PeriodicPythonRunner.hpp'
--- pkg/common/Engine/GlobalEngine/PeriodicPythonRunner.hpp 2009-12-09 17:11:51 +0000
+++ pkg/common/Engine/GlobalEngine/PeriodicPythonRunner.hpp 2010-01-01 15:21:30 +0000
@@ -3,6 +3,7 @@
#include<yade/core/GlobalEngine.hpp>
#include<yade/core/Scene.hpp>
#include<yade/pkg-common/PeriodicEngines.hpp>
+#include<yade/lib-pyutil/gil.hpp>
/* Execute a python command (in addPlotDataCall) periodically, with defined (and adjustable) periodicity.
*
* Period constraints are iterInterval and timeInterval. When either of them is exceeded, the addPlotDataCall is run.
@@ -15,12 +16,7 @@
public :
PeriodicPythonRunner(): command("pass"){};
/* virtual bool isActivated: not overridden, StretchPeriodicEngine handles that */
- virtual void action(Scene* b){
- PyGILState_STATE gstate;
- gstate = PyGILState_Ensure();
- PyRun_SimpleString(command.c_str()); // this is suboptimal, since it has to be parsed at every execution; critical?
- PyGILState_Release(gstate);
- }
+ virtual void action(Scene* b){ pyRunString(command); }
REGISTER_ATTRIBUTES(StretchPeriodicEngine,(command));
REGISTER_CLASS_NAME(PeriodicPythonRunner);
REGISTER_BASE_CLASS_NAME(StretchPeriodicEngine);
=== modified file 'pkg/dem/DataClass/SpherePack.cpp'
--- pkg/dem/DataClass/SpherePack.cpp 2009-12-20 22:03:17 +0000
+++ pkg/dem/DataClass/SpherePack.cpp 2010-01-01 15:21:30 +0000
@@ -78,14 +78,14 @@
if(scene->isPeriodic) { cellSize=scene->cell->getSize(); }
}
-long SpherePack::makeCloud(Vector3r mn, Vector3r mx, Real rMean, Real rRelFuzz, size_t num, bool periodic){
+long SpherePack::makeCloud(Vector3r mn, Vector3r mx, Real rMean, Real rRelFuzz, int num, bool periodic){
static boost::minstd_rand randGen(TimingInfo::getNow(/* get the number even if timing is disabled globally */ true));
static boost::variate_generator<boost::minstd_rand&, boost::uniform_real<> > rnd(randGen, boost::uniform_real<>(0,1));
- const size_t maxTry=1000;
+ const int maxTry=1000;
Vector3r size=mx-mn;
if(periodic)(cellSize=size);
- for(size_t i=0; i<num; i++) {
- size_t t;
+ for(int i=0; (i<num) || (num<0); i++) {
+ int t;
Real r=(rnd()-.5)*rRelFuzz*rMean+rMean;
for(t=0; t<maxTry; ++t){
Vector3r c;
@@ -104,7 +104,7 @@
if(!overlap) { pack.push_back(Sph(c,r)); break; }
}
if (t==maxTry) {
- LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<num);
+ if(num>0) LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<num);
return i;
}
}
=== modified file 'pkg/dem/DataClass/SpherePack.hpp'
--- pkg/dem/DataClass/SpherePack.hpp 2009-11-29 11:03:25 +0000
+++ pkg/dem/DataClass/SpherePack.hpp 2010-01-01 15:21:30 +0000
@@ -50,8 +50,8 @@
void toFile(const string file) const;
void fromSimulation();
- // random generation
- long makeCloud(Vector3r min, Vector3r max, Real rMean, Real rFuzz, size_t num, bool periodic=false);
+ // random generation; if num<0, insert as many spheres as possible
+ long makeCloud(Vector3r min, Vector3r max, Real rMean, Real rFuzz, int num, bool periodic=false);
// periodic repetition
void cellRepeat(Vector3<int> count);
=== modified file 'py/_eudoxos.cpp'
--- py/_eudoxos.cpp 2009-12-18 09:48:16 +0000
+++ py/_eudoxos.cpp 2010-01-01 15:21:30 +0000
@@ -1,7 +1,10 @@
+#include<yade/lib-pyutil/numpy.hpp>
+
#include<yade/pkg-dem/ConcretePM.hpp>
#include<boost/python.hpp>
#include<yade/extra/boost_python_len.hpp>
#include<yade/pkg-dem/Shop.hpp>
+
using namespace boost::python;
using namespace std;
#ifdef YADE_LOG4CXX
@@ -126,15 +129,36 @@
return ret;
}
#endif
+
void particleConfinement(){
CpmStateUpdater::update();
}
+python::dict testNumpy(){
+ Scene* scene=Omega::instance().getScene().get();
+ int dim1[]={scene->bodies->size()};
+ int dim2[]={scene->bodies->size(),3};
+ numpy_boost<Real,1> mass(dim1);
+ numpy_boost<Real,2> vel(dim2);
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if(!b) continue;
+ mass[b->getId()]=b->state->mass;
+ VECTOR3R_TO_NUMPY(vel[b->getId()],b->state->vel);
+ }
+ python::dict ret;
+ ret["mass"]=mass;
+ ret["vel"]=vel;
+ return ret;
+}
+
+
BOOST_PYTHON_MODULE(_eudoxos){
+ import_array();
def("velocityTowardsAxis",velocityTowardsAxis,velocityTowardsAxis_overloads(args("axisPoint","axisDirection","timeToAxis","subtractDist","perturbation")));
def("yieldSigmaTMagnitude",yieldSigmaTMagnitude);
// def("spiralSphereStresses2d",spiralSphereStresses2d,(python::arg("dH_dTheta"),python::arg("axis")=2));
def("particleConfinement",particleConfinement);
+ def("testNumpy",testNumpy);
}
=== modified file 'py/yadeWrapper/customConverters.cpp'
--- py/yadeWrapper/customConverters.cpp 2009-12-18 15:04:57 +0000
+++ py/yadeWrapper/customConverters.cpp 2010-01-01 15:21:30 +0000
@@ -9,6 +9,8 @@
#include<indexing_suite/vector.hpp>
#endif
+#include<yade/lib-pyutil/numpy.hpp>
+
#include<boost/python.hpp>
#include<boost/python/class.hpp>
#include<boost/python/module.hpp>
@@ -21,6 +23,8 @@
#include<string>
#include<stdexcept>
#include<iostream>
+#include<map>
+
#include<yade/lib-base/yadeWm3.hpp>
#include<yade/lib-base/yadeWm3Extra.hpp>
@@ -126,6 +130,13 @@
+template<typename numT, int dim>
+struct custom_numpyBoost_to_py{
+ static PyObject* convert(numpy_boost<numT, dim> nb){
+ return nb.py_ptr(); // handles incref internally
+ }
+};
+
@@ -137,6 +148,8 @@
custom_Se3r_from_seq(); to_python_converter<Se3r,custom_se3_to_tuple>();
// Vector3<int> to python (not implemented the other way around yet)
custom_vector3i_to_seq(); to_python_converter<Vector3<int>,custom_vector3i_to_seq>();
+ // StrArrayMap (typedef for std::map<std::string,numpy_boost>) → python dictionary
+ //custom_StrArrayMap_to_dict();
// register from-python converter and to-python converter
custom_vector_from_seq<int>(); to_python_converter<std::vector<int>, custom_vector_to_list<int> >();
custom_vector_from_seq<Real>(); to_python_converter<std::vector<Real>, custom_vector_to_list<Real> >();
@@ -148,6 +161,13 @@
custom_vector_from_seq<shared_ptr<Material> >(); to_python_converter<std::vector<shared_ptr<Material> >, custom_vector_to_list<shared_ptr<Material> > >();
custom_vector_from_seq<shared_ptr<Serializable> >(); to_python_converter<std::vector<shared_ptr<Serializable> >, custom_vector_to_list<shared_ptr<Serializable> > >();
+ import_array();
+ to_python_converter<numpy_boost<Real,1>, custom_numpyBoost_to_py<Real,1> >();
+ to_python_converter<numpy_boost<Real,2>, custom_numpyBoost_to_py<Real,2> >();
+ to_python_converter<numpy_boost<int,1>, custom_numpyBoost_to_py<int,1> >();
+ to_python_converter<numpy_boost<int,2>, custom_numpyBoost_to_py<int,2> >();
+
+
#define VECTOR_ENGINE_UNIT(ENGINE_UNIT) custom_vector_from_seq<shared_ptr<ENGINE_UNIT> >(); to_python_converter<std::vector<shared_ptr<ENGINE_UNIT> >,custom_vector_to_list<shared_ptr<ENGINE_UNIT> > >();
VECTOR_ENGINE_UNIT(BoundFunctor)
#ifdef YADE_GEOMETRICALMODEL
=== modified file 'scripts/test/periodic-triax-velgrad.py'
--- scripts/test/periodic-triax-velgrad.py 2009-12-30 15:45:32 +0000
+++ scripts/test/periodic-triax-velgrad.py 2010-01-01 15:21:30 +0000
@@ -7,12 +7,11 @@
#log.setLevel('Shop',log.TRACE)
O.periodic=True
O.cell.refSize=Vector3(.1,.1,.1)
-#O.cell.Hsize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
O.cell.trsf=Matrix3().IDENTITY;
sp=pack.SpherePack()
radius=5e-3
-num=sp.makeCloud(Vector3().ZERO,O.cell.refSize,radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
+num=sp.makeCloud(Vector3().ZERO,O.cell.refSize,radius,.6,-1,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
O.bodies.append([utils.sphere(s[0],s[1]) for s in sp])
@@ -25,14 +24,15 @@
[SimpleElasticRelationships()],
[Law2_Dem3Dof_Elastic_Elastic()]
),
- #PeriTriaxController(goal=[-1e5,-1e5,0],stressMask=3,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='triaxDone()',label='triax'),
- NewtonIntegrator(damping=.6, homotheticCellResize=1),
- PeriodicPythonRunner(command='utils.flipCell()',iterPeriod=1000), # broken for larger strains?
+ #PeriTriaxController(goal=[-1e3,-1e3,-1e3],stressMask=7,globUpdate=5,maxStrainRate=[5.,5.,5.],label='triax'),
+ NewtonIntegrator(damping=.6, homotheticCellResize=0),
+ PeriodicPythonRunner(command='utils.flipCell()',iterPeriod=1000),
]
O.dt=0.5*utils.PWaveTimeStep()
O.run(1)
qt.View()
-O.cell.velGrad=Matrix3(0,5,0,0,0,0,0,0,0)
+O.cell.velGrad=Matrix3(0,10,0,0,0,0,0,0,0)
+#O.cell.velGrad=Matrix3(0,5,2,-5,0,0,-2,0,0)
O.saveTmp()
O.run();
rrr=qt.Renderer(); rrr['intrAllWire'],rrr['Body_interacting_geom']=True,False
