yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
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(),
 ]