yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06447
[Branch ~yade-dev/yade/trunk] Rev 2592: 1. Fixes in L3Geom, L6Geom & friends
------------------------------------------------------------
revno: 2592
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Mon 2010-12-06 23:58:31 +0100
message:
1. Fixes in L3Geom, L6Geom & friends
2. Add scripts/test/beam-l6geom.py to demonstrate L6Geom
3. Fix LawTester so that ti works with all 6 dofs reliably now
added:
scripts/test/beam-l6geom.py
modified:
pkg/dem/DomainLimiter.cpp
pkg/dem/L3Geom.cpp
scripts/test/law-test.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 'pkg/dem/DomainLimiter.cpp'
--- pkg/dem/DomainLimiter.cpp 2010-12-06 20:05:12 +0000
+++ pkg/dem/DomainLimiter.cpp 2010-12-06 22:58:31 +0000
@@ -118,12 +118,16 @@
}
// update the transformation
// the matrix is orthonormal, since axX, axY are normalized and and axZ is their corss-product
- trsf.col(0)=axX; trsf.col(1)=axY; trsf.col(2)=axZ;
+ trsf.row(0)=axX; trsf.row(1)=axY; trsf.row(2)=axZ;
} else {
trsf=l3Geom->trsf;
axX=trsf.row(0); axY=trsf.row(1); axZ=trsf.row(2);
ptGeom=l3Geom->u;
- if(l6Geom) rotGeom=l6Geom->phi;
+ if(l6Geom){
+ rotGeom=l6Geom->phi;
+ // perform allshearing by translation, as it does not induce bending
+ if(rotWeight!=0){ LOG_INFO("LawTester.rotWeight set to 0 (was"<<rotWeight<<"), since rotational DoFs are in use."); rotWeight=0; }
+ }
}
trsfQ=Quaternionr(trsf);
contPt=gsc->contactPoint;
@@ -187,8 +191,8 @@
// and the compensation is always in the -εx sense (-sign â +1 for #0, -1 for #1)
vel[i]+=-sign*axX*radius*((1-cos(arcAngleY))+(1-cos(arcAngleZ)))/scene->dt;
- // rotation
- angVel[i]+=trsf*ddPhi;
+ // rotation, convert from local to global
+ angVel[i]+=trsf.transpose()*ddPhi;
LOG_DEBUG("vel="<<vel[i]<<", angVel="<<angVel[i]<<", rotY,rotZ="<<rotY<<","<<rotZ<<", arcAngle="<<arcAngleY<<","<<arcAngleZ<<", sign="<<sign<<", weight="<<weight);
=== modified file 'pkg/dem/L3Geom.cpp'
--- pkg/dem/L3Geom.cpp 2010-12-05 17:10:06 +0000
+++ pkg/dem/L3Geom.cpp 2010-12-06 22:58:31 +0000
@@ -9,7 +9,7 @@
#include<GL/glu.h>
#endif
-YADE_PLUGIN((L3Geom)(Ig2_Sphere_Sphere_L3Geom_Inc)(Ig2_Wall_Sphere_L3Geom_Inc)(Ig2_Sphere_Sphere_L6Geom_Inc)(Law2_L3Geom_FrictPhys_ElPerfPl)(Law2_L6Geom_FrictPhys_Linear)
+YADE_PLUGIN((L3Geom)(L6Geom)(Ig2_Sphere_Sphere_L3Geom_Inc)(Ig2_Wall_Sphere_L3Geom_Inc)(Ig2_Sphere_Sphere_L6Geom_Inc)(Law2_L3Geom_FrictPhys_ElPerfPl)(Law2_L6Geom_FrictPhys_Linear)
#ifdef YADE_OPENGL
(Gl1_L3Geom)
#endif
@@ -74,11 +74,12 @@
g.refR1=r1; g.refR2=r1;
g.normal=normal; const Vector3r& locX(g.normal);
// initial local y-axis orientation, in the xz or xy plane, depending on which component is larger to avoid singularities
- Vector3r locY=g.normal.cross(g.normal[1]>g.normal[2]?Vector3r::UnitY():Vector3r::UnitZ());
- Vector3r locZ=g.normal.cross(locY);
+ Vector3r locY=normal.cross(normal[1]<normal[2]?Vector3r::UnitY():Vector3r::UnitZ());
+ Vector3r locZ=normal.cross(locY);
g.trsf.row(0)=locX; g.trsf.row(1)=locY; g.trsf.row(2)=locZ;
g.u=Vector3r(uN,0,0); // zero shear displacement
// L6Geom::phi is initialized to Vector3r::Zero() automatically
+ //cerr<<"Init trsf=\n"<<g.trsf<<endl<<"locX="<<locX<<", locY="<<locY<<", locZ="<<locZ<<endl;
return true;
}
@@ -141,7 +142,7 @@
// if requested via approxMask, just use prevTrsf
Quaternionr midTrsf=(approxMask|APPROX_NO_MID_TRSF) ? Quaternionr(prevTrsf) : Quaternionr(prevTrsf).slerp(.5,Quaternionr(currTrsf));
- // cerr<<"prevTrsf=\n"<<prevTrsf<<", currTrsf=\n"<<currTrsf<<", midTrsf=\n"<<Matrix3r(midTrsf)<<endl;
+ //cerr<<"prevTrsf=\n"<<prevTrsf<<", currTrsf=\n"<<currTrsf<<", midTrsf=\n"<<Matrix3r(midTrsf)<<endl;
// updates of geom here
=== added file 'scripts/test/beam-l6geom.py'
--- scripts/test/beam-l6geom.py 1970-01-01 00:00:00 +0000
+++ scripts/test/beam-l6geom.py 2010-12-06 22:58:31 +0000
@@ -0,0 +1,29 @@
+#
+# Demonstrate L6Geom class with suspended hook-shaped beam fixed at one side, subject to gravity
+#
+import numpy
+# radius, number and distance of spheres
+rad,num=1,6; dist=1.9999*rad
+# one arm
+O.bodies.append([utils.sphere((0,y,0),rad,wire=True) for y in numpy.arange(0,2*num-1,dist)])
+# the lateral arm
+O.bodies.append([utils.sphere((x,(num-1)*2*rad,0),rad,wire=True) for x in numpy.arange(dist,1+num/2,dist)])
+# support sphere
+O.bodies[0].state.blockedDOFs=['x','y','z','rx','ry','rz']
+# small dt to see in realtime how it swings; real critical is higher, but much less than p-wave
+O.dt=.001*utils.PWaveTimeStep()
+
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb()]),
+ InteractionLoop([Ig2_Sphere_Sphere_L6Geom_Inc()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L6Geom_FrictPhys_Linear(charLen=1)]),
+ GravityEngine(gravity=(0,0,-9.81)),
+ NewtonIntegrator(damping=0.03),
+]
+O.saveTmp()
+
+try:
+ from yade import qt
+ v=qt.View(); v.axes=True
+except: pass
+O.run()
=== modified file 'scripts/test/law-test.py'
--- scripts/test/law-test.py 2010-12-05 17:10:06 +0000
+++ scripts/test/law-test.py 2010-12-06 22:58:31 +0000
@@ -22,28 +22,28 @@
pt1=Vector3(0,0,0)
# random orientation of the interaction
normal=Vector3(random.random()-.5,random.random()-.5,random.random()-.5)
-#normal=Vector3(1,0,0)
+normal=Vector3.UnitX
O.bodies.append([
utils.sphere(pt1,r1,wire=True,color=(.7,.7,.7)),
utils.sphere(pt1+.999999*(r1+r2)*normal.normalized(),r2,wire=True,color=(0,0,0))
])
O.engines=[
+ PyRunner(iterPeriod=1,command='addPlotData()'),
ForceResetter(),
#PyRunner(iterPeriod=1,command='import time; time.sleep(.005)'),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
#[Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()] # ScGeom
- [Ig2_Sphere_Sphere_L3Geom_Inc(approxMask=63)],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L3Geom_FrictPhys_ElPerfPl(noBreak=True,noSlip=False)] # L3Geom
- #[Ig2_Sphere_Sphere_L6Geom_Inc(approxMask=63)],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L6Geom_FrictPhys_Linear(charLen=1)] # L6Geom
+ #[Ig2_Sphere_Sphere_L3Geom_Inc(approxMask=63)],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L3Geom_FrictPhys_ElPerfPl(noBreak=True,noSlip=False)] # L3Geom
+ [Ig2_Sphere_Sphere_L6Geom_Inc(approxMask=63)],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L6Geom_FrictPhys_Linear(charLen=1)] # L6Geom
),
- #LawTester(ids=[0,1],path=[(0,0,0)]*7+[(-1e-5,0,0),(-1e-5,.1,.1)],rotPath=[(0,.2,0),(0,0,0),(0,0,.2),(0,0,0),(.2,0,0),(-.2,0,0),(0,0,0)],pathSteps=[10],doneHook='tester.dead=True; O.pause();',label='tester',rotWeight=0),
- LawTester(ids=[0,1],path=[
- (-1e-5,0,0),(-.1,0,0),(-.1,.1,0),(-1e-5,.1,0), # towards, shear, back to intial normal distance
- (-.02,.1,.1),(-.02,-.1,.1),(-.02,-.1,-.1),(-.02,.1,-.1),(-.02,.1,.1), # go in square in the shear plane without changing normal deformation
- (-1e-4,0,0) # back to the origin, but keep some overlap to not delete the interaction
- ],pathSteps=[100],doneHook='tester.dead=True; O.pause()',label='tester',rotWeight=.2,idWeight=.2),
- PyRunner(iterPeriod=1,command='addPlotData()'),
+ LawTester(ids=[0,1],path=[(0,0,0)]*7+[(-1e-5,0,0),(-1e-5,.1,.1)],rotPath=[(0,.2,0),(0,0,0),(0,0,.2),(0,0,0),(.2,0,0),(-.2,0,0),(0,0,0)],pathSteps=[10],doneHook='tester.dead=True; O.pause();',label='tester',rotWeight=0),
+ #LawTester(ids=[0,1],path=[
+ # (-1e-5,0,0),(-.1,0,0),(-.1,.1,0),(-1e-5,.1,0), # towards, shear, back to intial normal distance
+ # (-.02,.1,.1),(-.02,-.1,.1),(-.02,-.1,-.1),(-.02,.1,-.1),(-.02,.1,.1), # go in square in the shear plane without changing normal deformation
+ # (-1e-4,0,0) # back to the origin, but keep some overlap to not delete the interaction
+ # ],pathSteps=[100],doneHook='tester.dead=True; O.pause()',label='tester',rotWeight=.2,idWeight=.2),
NewtonIntegrator(),
]
Follow ups