yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2519
committer: Sergei D. <sega@think>
branch nick: trunk
timestamp: Wed 2010-10-27 22:54:33 +0400
message:
  1. Fix ViscElBasic law to handle periodicity; 2. Add some regression tests
added:
  scripts/test/facet-sphere-ViscElBasic-peri.py
  scripts/test/facet-sphere-ViscElBasic.py
  scripts/test/sphere-sphere-ViscElBasic-peri.py
modified:
  pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp
  pkg/dem/ViscoelasticPM.cpp


--
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/Ig2_Facet_Sphere_ScGeom.cpp'
--- pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp	2010-10-25 20:31:09 +0000
+++ pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp	2010-10-27 18:54:33 +0000
@@ -107,7 +107,7 @@
 		scm->radius1 = 2*sphereRadius;
 		scm->radius2 = sphereRadius;
 		if (isNew) c->geom = scm;
-		scm->precompute(state1,state2,scene,c,normal,isNew,shift2,true);
+		scm->precompute(state1,state2,scene,c,normal,isNew,shift2,false/*avoidGranularRatcheting only for sphere-sphere*/);
 		return true;
 	}
 	return false;

=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp	2010-10-25 20:31:09 +0000
+++ pkg/dem/ViscoelasticPM.cpp	2010-10-27 18:54:33 +0000
@@ -67,20 +67,17 @@
 	axis = angle*geom.normal;
 	shearForce -= shearForce.cross(axis);
 
-	//const Vector3r c1x = (geom.contactPoint - de1.pos);
-	//const Vector3r c2x = (geom.contactPoint - de2.pos);
-	const Vector3r c1x = geom.radius1*geom.normal;
-	const Vector3r c2x = -geom.radius2*geom.normal;
-	/// The following definition of c1x and c2x is to avoid "granular ratcheting" 
-	///  (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN, 
-	///   "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
-	///   ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors)
-	//Vector3r _c1x_	=  geom->radius1*geom->normal;
-	//Vector3r _c2x_	= -geom->radius2*geom->normal;
-	//Vector3r relativeVelocity		= (de2->velocity+de2->angularVelocity.Cross(_c2x_)) - (de1->velocity+de1->angularVelocity.Cross(_c1x_));
-	const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) ;
+	// Handle periodicity.
+	const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero(); 
+	const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero(); 
+
+	const Vector3r c1x = (geom.contactPoint - de1.pos);
+	const Vector3r c2x = (geom.contactPoint - de2.pos - shift2);
+	
+	const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel;
 	const Real normalVelocity	= geom.normal.dot(relativeVelocity);
 	const Vector3r shearVelocity	= relativeVelocity-normalVelocity*geom.normal;
+
 	// As Chiara Modenese suggest, we store the elastic part 
 	// and then add the viscous part if we pass the Mohr-Coulomb criterion.
 	// See http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg01391.html

=== added file 'scripts/test/facet-sphere-ViscElBasic-peri.py'
--- scripts/test/facet-sphere-ViscElBasic-peri.py	1970-01-01 00:00:00 +0000
+++ scripts/test/facet-sphere-ViscElBasic-peri.py	2010-10-27 18:54:33 +0000
@@ -0,0 +1,53 @@
+# -*- coding: utf-8
+
+# Testing facet-sphere interaction in periodic case.
+# Pass, if the sphere is rolling from left to right through the period.
+
+from yade import utils
+
+sphereRadius=0.1
+tc=0.001# collision time 
+en=0.3  # normal restitution coefficient
+es=0.3  # tangential restitution coefficient
+density=2700
+frictionAngle=radians(35)# 
+params=utils.getViscoelasticFromSpheresInteraction(tc,en,es)
+facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,**params)) 
+sphereMat=O.materials.append(ViscElMat(density=density,frictionAngle=frictionAngle,**params))
+
+#floor
+n=5.
+s=1./n
+for i in range(0,n):
+	for j in range(0,n):
+		O.bodies.append([
+			utils.facet( [(i*s,j*s,0.1),(i*s,(j+1)*s,0.1),((i+1)*s,(j+1)*s,0.1)],material=facetMat),
+			utils.facet( [(i*s,j*s,0.1),((i+1)*s,j*s,0.1),((i+1)*s,(j+1)*s,0.1)],material=facetMat),
+		])
+
+# Spheres
+sphId=O.bodies.append([
+	utils.sphere( (0.5,0.5,0.2), 0.1, material=sphereMat),
+	])
+O.bodies[sphId[-1]].state.vel=(0.5,0,0)
+
+## Engines 
+O.engines=[
+	ForceResetter(),
+	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
+	InteractionLoop(
+		[Ig2_Facet_Sphere_ScGeom()],
+		[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
+		[Law2_ScGeom_ViscElPhys_Basic()],
+	),
+	GravityEngine(gravity=[0,0,-9.81]),
+	NewtonIntegrator(damping=0),
+]
+
+O.periodic=True
+O.cell.refSize=Vector3(1,1,1)
+
+O.dt=.01*tc
+
+O.saveTmp()
+

=== added file 'scripts/test/facet-sphere-ViscElBasic.py'
--- scripts/test/facet-sphere-ViscElBasic.py	1970-01-01 00:00:00 +0000
+++ scripts/test/facet-sphere-ViscElBasic.py	2010-10-27 18:54:33 +0000
@@ -0,0 +1,48 @@
+# -*- coding: utf-8 -*-
+
+# Testing facet-sphere interaction.
+# A facet is rotated around Z axis.  Test pass, if a sphere at (0,0) position is not moving (because in this case no transfer moment from the facet to the sphere), but a sphere at facet's edge moves with the facet (for this sphere blocked the rotation DOFs in order to remove rolling).
+
+## PhysicalParameters 
+Density=2400
+frictionAngle=radians(35)
+tc = 0.001
+en = 0.3
+es = 0.3
+
+## Import wall's geometry
+params=utils.getViscoelasticFromSpheresInteraction(tc,en,es)
+facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,**params)) 
+sphereMat=O.materials.append(ViscElMat(density=Density,frictionAngle=frictionAngle,**params))
+
+facetId=O.bodies.append(utils.facet( [ (-1,0,0), (1,1,0), (1,-1,0)], material=facetMat,color=(1,0,0)))
+
+sphIds=O.bodies.append([
+	utils.sphere( (0,0,0.1),0.1, material=sphereMat,color=(0,1,0)), 
+	utils.sphere( (0.9,0,0.1),0.1, material=sphereMat,color=(0,1,0))
+	])
+
+O.bodies[sphIds[1]].state.blockedDOFs=['rx','ry','rz']
+
+## Timestep 
+O.dt=.1*tc
+
+## Engines 
+O.engines=[
+	ForceResetter(),
+	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
+	InteractionLoop(
+		[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
+		[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
+		[Law2_ScGeom_ViscElPhys_Basic()],
+	),
+	GravityEngine(gravity=[0,0,-9.81]),
+	NewtonIntegrator(damping=0),
+	RotationEngine(subscribedBodies=[facetId],rotationAxis=[0,0,1],rotateAroundZero=True,angularVelocity=0.1)
+]
+
+from yade import qt
+qt.View()
+O.saveTmp()
+#O.run()
+

=== added file 'scripts/test/sphere-sphere-ViscElBasic-peri.py'
--- scripts/test/sphere-sphere-ViscElBasic-peri.py	1970-01-01 00:00:00 +0000
+++ scripts/test/sphere-sphere-ViscElBasic-peri.py	2010-10-27 18:54:33 +0000
@@ -0,0 +1,44 @@
+# -*- coding: utf-8
+
+# Testing sphere-sphere interaction in periodic case.
+# Pass, if the spheres moves along the X axis, interacting through the period.
+
+from yade import utils
+
+sphereRadius=0.1
+tc=0.001# collision time 
+en=1  # normal restitution coefficient
+es=1  # tangential restitution coefficient
+density=2700
+frictionAngle=radians(35)# 
+params=utils.getViscoelasticFromSpheresInteraction(tc,en,es)
+sphereMat=O.materials.append(ViscElMat(density=density,frictionAngle=frictionAngle,**params))
+
+
+# Spheres
+sphId=O.bodies.append([
+	utils.sphere( (0.4,0.5,0.5), 0.1, material=sphereMat),
+	utils.sphere( (0.6,0.5,0.5), 0.1, material=sphereMat)
+	])
+O.bodies[sphId[-1]].state.vel=(0.5,0,0)
+O.bodies[sphId[0]].state.vel=(-0.5,0,0)
+
+## Engines 
+O.engines=[
+	ForceResetter(),
+	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
+	InteractionLoop(
+		[Ig2_Sphere_Sphere_ScGeom()],
+		[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
+		[Law2_ScGeom_ViscElPhys_Basic()],
+	),
+	NewtonIntegrator(damping=0),
+]
+
+O.periodic=True
+O.cell.refSize=Vector3(1,1,1)
+
+O.dt=.01*tc
+
+O.saveTmp()
+