yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02873
[Branch ~yade-dev/yade/trunk] Rev 1928: 1. fix size computation in cell
------------------------------------------------------------
revno: 1928
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Wed 2009-12-30 16:45:32 +0100
message:
1. fix size computation in cell
2. Save history at exit, incremental history search with arrows up/down
3. Don't return false in Ig2_Box_Sphere_ScGeom for existing intr
modified:
core/Cell.hpp
core/main/main.py.in
lib/SConscript
lib/opengl/GLUtils.cpp
lib/opengl/GLUtils.hpp
pkg/common/DataClass/Bound/Aabb.hpp
pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp
pkg/common/RenderingEngine/Gl1_Aabb.cpp
pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp
pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp
py/yadeWrapper/yadeWrapper.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.
=== modified file 'core/Cell.hpp'
--- core/Cell.hpp 2009-12-26 21:57:37 +0000
+++ core/Cell.hpp 2009-12-30 15:45:32 +0000
@@ -26,8 +26,8 @@
Vector3r getSize_copy() const { return _size; }
//! stretch ratio Λ(n) (http://en.wikipedia.org/wiki/Finite_strain_theory#Stretch_ratio)
Vector3r getStretchRatio() const { return Vector3r(trsf[0][0],trsf[1][1],trsf[2][2]); }
- //! return matrix containing cosines of shear angles
- const Matrix3r& getCosMatrix() const {return _cosMatrix;}
+ //! return vector of consines of skew angle in yz, xz, xy planes between respective transformed base vectors
+ const Vector3r& getCos() const {return _cos;}
//! transformation matrix applying puse shear (normal strain removed: ones on the diagonal)
const Matrix3r& getShearTrsf() const { return _shearTrsf; }
//! inverse of getShearTrsfMatrix().
@@ -54,14 +54,14 @@
// caches; private
private:
Matrix3r _trsfInc;
- Vector3r _size;
+ Vector3r _size, _cos;
bool _hasShear;
- Matrix3r _shearTrsf, _unshearTrsf, _cosMatrix;
+ Matrix3r _shearTrsf, _unshearTrsf;
double _glShearTrsfMatrix[16];
void fillGlShearTrsfMatrix(double m[16]){
- m[0]=1; m[4]=trsf[0][1]; m[8]=trsf[0][2]; m[12]=0;
- m[1]=trsf[1][0]; m[5]=1; m[9]=trsf[1][2]; m[13]=0;
- m[2]=trsf[2][0]; m[6]=trsf[2][1]; m[10]=1; m[14]=0;
+ m[0]=trsf[0][0]; m[4]=trsf[0][1]; m[8]=trsf[0][2]; m[12]=0;
+ m[1]=trsf[1][0]; m[5]=trsf[1][1]; m[9]=trsf[1][2]; m[13]=0;
+ m[2]=trsf[2][0]; m[6]=trsf[2][1]; m[10]=trsf[2][2];m[14]=0;
m[3]=0; m[7]=0; m[11]=0; m[15]=1;
}
@@ -69,34 +69,35 @@
//! "integrate" velGrad, update cached values used by public getter
void integrateAndUpdate(Real dt){
- #if 0
- //initialize Hsize for "lazy" default scripts, after that Hsize is free to change
- if (refSize[0]!=1 && Hsize[0][0]==0) {Hsize[0][0]=refSize[0]; Hsize[1][1]=refSize[1]; Hsize[2][2]=refSize[2];}
- //update Hsize (redundant with total transformation perhaps)
- Hsize=Hsize+_shearInc*Hsize;
- #endif
-
//incremental displacement gradient
_trsfInc=dt*velGrad;
// total transformation; M = (Id+G).M = F.M
trsf+=_trsfInc*trsf;
- #if 0
- cerr<<"velGrad "; for(int i=0; i<3; i++) for(int j=0; j<3; j++) cerr<<velGrad[i][j]<<endl;
- cerr<<"_strainInc "; for(int i=0; i<3; i++) for(int j=0; j<3; j++) cerr<<_strainInc[i][j]<<endl;
- cerr<<"strain "; for(int i=0; i<3; i++) for(int j=0; j<3; j++) cerr<<strain[i][j]<<endl;
- #endif
+ // Hsize will contain colums with transformed base vectors
+ Matrix3r Hsize(refSize[0],refSize[1],refSize[2]);
+ Hsize=trsf*Hsize;
+ // lengths of transformed cell vectors, skew cosines
+ Vector3r Hnorm[3]; // normalized transformed base vectors
+ for(int i=0; i<3; i++){
+ _size[i]=Vector3r(Hsize[i][0],Hsize[i][1],Hsize[i][2]).Length();
+ Hnorm[i]=Vector3r(Hsize[i][0],Hsize[i][1],Hsize[i][2]); Hnorm[i].Normalize();
+ };
+ //cerr<<"Hsize="<<endl<<Hsize[0][0]<<" "<<Hsize[0][1]<<" "<<Hsize[0][2]<<endl<<Hsize[1][0]<<" "<<Hsize[1][1]<<" "<<Hsize[1][2]<<endl<<Hsize[2][0]<<" "<<Hsize[2][1]<<" "<<Hsize[2][2]<<endl;
+ //cerr<<"Hnorm="<<Hnorm[0]<<" | "<<Hnorm[1]<<" | "<<Hnorm[2]<<endl;
+ // skew cosines
+ for(int i=0; i<3; i++){
+ int i1=(i+1)%3, i2=(i+2)%3;
+ // sin between axes is cos of skew
+ _cos[i]=(Hnorm[i1].Cross(Hnorm[i2])).SquaredLength();
+ }
// pure shear trsf: ones on diagonal
_shearTrsf=trsf; _shearTrsf[0][0]=_shearTrsf[1][1]=_shearTrsf[2][2]=1.;
// pure unshear transformation
_unshearTrsf=_shearTrsf.Inverse();
// some parts can branch depending on presence/absence of shear
_hasShear=(trsf[0][1]!=0 || trsf[0][2]!=0 || trsf[1][0]!=0 || trsf[1][2]!=0 || trsf[2][0]!=0 || trsf[2][1]!=0);
- // current cell size (in units on physical space axes)
- _size=diagMult(getStretchRatio(),refSize);
// OpenGL shear matrix (used frequently)
fillGlShearTrsfMatrix(_glShearTrsfMatrix);
- // precompute cosines of angles, used for enlarge factor when computing Aabb's
- for(int i=0; i<3; i++) for(int j=0; j<3; j++) _cosMatrix[i][j]=(i==j?0:cos(atan(trsf[i][j])));
}
/*! Return point inside periodic cell, even if shear is applied */
=== modified file 'core/main/main.py.in'
--- core/main/main.py.in 2009-12-20 22:03:17 +0000
+++ core/main/main.py.in 2009-12-30 15:45:32 +0000
@@ -104,9 +104,14 @@
'"\e[21~": "\C-Uyade.qt.Controller(), yade.qt.View();\C-M"', # F10
'"\e[20~": "\C-Uyade.qt.Generator();\C-M"', #F9
'"\e[19~": "\C-Uyade.plot.plot();\C-M"', #F9
+ '"\e[A": history-search-backward', '"\e[B": history-search-forward', # incremental history forward/backward
] if qtEnabled else [])
})
ipshell()
+ # save history -- a workaround for atexit handlers not being run (why?)
+ # http://lists.ipython.scipy.org/pipermail/ipython-user/2008-September/005839.html
+ import IPython.ipapi
+ IPython.ipapi.get().IP.atexit_operations()
else:
import bpython
bpython.embed()
=== modified file 'lib/SConscript'
--- lib/SConscript 2009-12-22 18:19:41 +0000
+++ lib/SConscript 2009-12-30 15:45:32 +0000
@@ -62,7 +62,7 @@
'serialization/FormatChecker.cpp','serialization/Serializable.cpp','serialization/SerializationExceptions.cpp',
]
# compile TesselationWrapper only if cgal is enabled
- +(Split('triangulation/KinematicLocalisationAnalyser.cpp triangulation/Operations.cpp triangulation/RegularTriangulation.cpp triangulation/Timer.cpp triangulation/basicVTKwritter.cpp triangulation/FlowBoundingSphere.cpp triangulation/Deformation.cpp triangulation/Empilement.cpp triangulation/stdafx.cpp triangulation/Tenseur3.cpp triangulation/Tesselation.cpp triangulation/TesselationWrapper.cpp triangulation/test.cpp triangulation/TriaxialState.cpp') if 'cgal' in env['features'] else []))
+ +(Split('triangulation/KinematicLocalisationAnalyser.cpp triangulation/Operations.cpp triangulation/RegularTriangulation.cpp triangulation/Timer.cpp triangulation/basicVTKwritter.cpp triangulation/FlowBoundingSphere.cpp triangulation/Deformation.cpp triangulation/Empilement.cpp triangulation/stdafx.cpp triangulation/Tenseur3.cpp triangulation/Tesselation.cpp triangulation/TesselationWrapper.cpp triangulation/test.cpp triangulation/TriaxialState.cpp') if 'cgal' in env['features'] else [])),
],LIBS=['dl'],CXXFLAGS=env['CXXFLAGS']+['-fPIC']
)
#)
=== modified file 'lib/opengl/GLUtils.cpp'
--- lib/opengl/GLUtils.cpp 2009-03-19 17:55:37 +0000
+++ lib/opengl/GLUtils.cpp 2009-12-30 15:45:32 +0000
@@ -1,1 +1,13 @@
#include"GLUtils.hpp"
+
+void GLUtils::Parallelepiped(const Vector3r& a, const Vector3r& b, const Vector3r& c){
+ glBegin(GL_LINE_STRIP);
+ glVertex3v(b); glVertex3v(Vector3r::ZERO); glVertex3v(a); glVertex3v(a+b); glVertex3v(a+b+c); glVertex3v(b+c); glVertex3v(b); glVertex3v(a+b);
+ glEnd();
+ glBegin(GL_LINE_STRIP);
+ glVertex3v(b+c); glVertex3v(c); glVertex3v(a+c); glVertex3v(a);
+ glEnd();
+ glBegin(GL_LINES);
+ glVertex3v(a+c); glVertex3v(a+b+c);
+ glEnd();
+}
=== modified file 'lib/opengl/GLUtils.hpp'
--- lib/opengl/GLUtils.hpp 2009-05-22 21:06:02 +0000
+++ lib/opengl/GLUtils.hpp 2009-12-30 15:45:32 +0000
@@ -12,6 +12,8 @@
#include<string>
struct GLUtils{
+ // render wire of parallelepiped with sides given by vectors a,b,c; zero corner is at origin
+ static void Parallelepiped(const Vector3r& a, const Vector3r& b, const Vector3r& c);
static void GLDrawArrow(const Vector3r& from, const Vector3r& to, const Vector3r& color=Vector3r(1,1,1)){
glEnable(GL_LIGHTING); glColor3v(color); qglviewer::Vec a(from[0],from[1],from[2]),b(to[0],to[1],to[2]); QGLViewer::drawArrow(a,b);
}
=== modified file 'pkg/common/DataClass/Bound/Aabb.hpp'
--- pkg/common/DataClass/Bound/Aabb.hpp 2009-12-18 15:04:57 +0000
+++ pkg/common/DataClass/Bound/Aabb.hpp 2009-12-30 15:45:32 +0000
@@ -22,7 +22,7 @@
virtual ~Aabb();
REGISTER_CLASS_AND_BASE(Aabb,Bound);
- REGISTER_ATTRIBUTES(Bound, /* (min)(max) */ ); // not necessary to store min and max, but it is handy for debugging/python inspection
+ REGISTER_ATTRIBUTES(Bound, (min)(max) ); // not necessary to store min and max, but it is handy for debugging/python inspection
REGISTER_CLASS_INDEX(Aabb,Bound);
};
REGISTER_SERIALIZABLE(Aabb);
=== modified file 'pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp'
--- pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp 2009-12-20 22:03:17 +0000
+++ pkg/common/Engine/Functor/Bo1_Sphere_Aabb.cpp 2009-12-30 15:45:32 +0000
@@ -15,10 +15,17 @@
Aabb* aabb = static_cast<Aabb*>(bv.get());
Vector3r halfSize = (aabbEnlargeFactor>0?aabbEnlargeFactor:1.)*Vector3r(sphere->radius,sphere->radius,sphere->radius);
// adjust box size along axes so that sphere doesn't stick out of the box even if sheared (i.e. parallelepiped)
- if(scene->isPeriodic) {
- const Matrix3r& cos=scene->cell->getCosMatrix();
- for(int i=0; i<3; i++) halfSize[i]*=(1/cos[i][(i+1)%3])*(1/cos[i][(i+2)%3]);
+ if(scene->isPeriodic && scene->cell->hasShear()) {
+ Vector3r refHalfSize(halfSize);
+ const Vector3r& cos=scene->cell->getCos();
+ for(int i=0; i<3; i++){
+ //cerr<<"cos["<<i<<"]"<<cos[i]<<" ";
+ int i1=(i+1)%3,i2=(i+2)%3;
+ halfSize[i1]+=refHalfSize[i1]*(1/cos[i]-1);
+ halfSize[i2]+=refHalfSize[i2]*(1/cos[i]-1);
+ }
}
+ //cerr<<" || "<<halfSize<<endl;
aabb->min = scene->cell->unshearPt(se3.position)-halfSize;
aabb->max = scene->cell->unshearPt(se3.position)+halfSize;
}
=== modified file 'pkg/common/RenderingEngine/Gl1_Aabb.cpp'
--- pkg/common/RenderingEngine/Gl1_Aabb.cpp 2009-12-20 22:03:17 +0000
+++ pkg/common/RenderingEngine/Gl1_Aabb.cpp 2009-12-30 15:45:32 +0000
@@ -14,10 +14,15 @@
void Gl1_Aabb::go(const shared_ptr<Bound>& bv, Scene* scene){
Aabb* aabb = static_cast<Aabb*>(bv.get());
glColor3v(bv->diffuseColor);
- glTranslatev(scene->cell->wrapShearedPt(scene->cell->shearPt(.5*(aabb->min+aabb->max))));
- // glDisable(GL_LIGHTING); // ??
- if(scene->isPeriodic) glMultMatrixd(scene->cell->getGlShearTrsfMatrix());
- glScalev(aabb->max-aabb->min);
+ // glDisable(GL_LIGHTING);
+ if(!scene->isPeriodic){
+ glTranslatev(.5*(aabb->min+aabb->max));
+ glScalev(aabb->max-aabb->min);
+ } else {
+ glTranslatev(scene->cell->shearPt(scene->cell->wrapPt(.5*(aabb->min+aabb->max))));
+ glMultMatrixd(scene->cell->getGlShearTrsfMatrix());
+ glScalev(aabb->max-aabb->min);
+ }
glutWireCube(1);
}
=== modified file 'pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp'
--- pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2009-12-25 21:58:25 +0000
+++ pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2009-12-30 15:45:32 +0000
@@ -138,12 +138,12 @@
if(!scene->isPeriodic) return;
glColor3v(Vector3r(1,1,0));
glPushMatrix();
- // order matters
- Vector3r size=scene->cell->getSize();
- if(scaleDisplacements) size=diagMult(displacementScale,size);
- glTranslatev(scene->cell->shearPt(.5*size)); // shear center (moves when sheared)
+ Vector3r refSize=scene->cell->refSize;
+ if(scaleDisplacements) refSize=diagMult(displacementScale,refSize);
+ //glTranslatev(scene->cell->shearPt(.5*size)); // shear center (moves when sheared)
+ glTranslatev(scene->cell->trsf*(.5*refSize));
glMultMatrixd(scene->cell->getGlShearTrsfMatrix());
- glScalev(size);
+ glScalev(refSize);
glutWireCube(1);
glPopMatrix();
}
=== modified file 'pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2009-12-13 20:30:13 +0000
+++ pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2009-12-30 15:45:32 +0000
@@ -110,7 +110,7 @@
Vector3r cOnBox_box = boxAxisT*cOnBox_boxLocal; // projection of sphere's center on closest box surface - relative to box's origin, but GLOBAL orientation!
Vector3r cOnBox_sphere = cOnBox_box-relPos21; // same, but origin in sphere's center
depth=s->radius-cOnBox_sphere.Length();
- if (depth<0) return false;
+ if (depth<0 && !c->isReal() && !force) return false;
/*
* +-----------------------------------+
=== modified file 'py/yadeWrapper/yadeWrapper.cpp'
--- py/yadeWrapper/yadeWrapper.cpp 2009-12-26 21:57:37 +0000
+++ py/yadeWrapper/yadeWrapper.cpp 2009-12-30 15:45:32 +0000
@@ -767,8 +767,8 @@
.def_readwrite("refSize",&Cell::refSize)
.def_readwrite("trsf",&Cell::trsf)
.def_readwrite("velGrad",&Cell::velGrad)
+ .def_readonly("size",&Cell::getSize_copy)
//.def_readwrite("Hsize",&Cell::Hsize)
- .add_property("stretch",&Cell::getStretchRatio)
//.add_property("size",&Cell::getSize,python::return_value_policy<python::return_internal_referece>()
;
=== modified file 'scripts/test/periodic-triax-velgrad.py'
--- scripts/test/periodic-triax-velgrad.py 2009-12-26 21:57:37 +0000
+++ scripts/test/periodic-triax-velgrad.py 2009-12-30 15:45:32 +0000
@@ -25,14 +25,15 @@
[SimpleElasticRelationships()],
[Law2_Dem3Dof_Elastic_Elastic()]
),
- PeriTriaxController(goal=[-1e5,-1e5,0],stressMask=3,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='triaxDone()',label='triax'),
+ #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?
]
O.dt=0.5*utils.PWaveTimeStep()
O.run(1)
qt.View()
-O.cell.velGrad=Matrix3(0,5,0,0,0,0, 0,0,-5)
+O.cell.velGrad=Matrix3(0,5,0,0,0,0,0,0,0)
+O.saveTmp()
O.run();
rrr=qt.Renderer(); rrr['intrAllWire'],rrr['Body_interacting_geom']=True,False
@@ -49,6 +50,6 @@
print 'Here we are: stress',triax['stress'],'strain',triax['strain'],'stiffness',triax['stiff']
print 'Done, pausing now.'
O.pause()
-
+
