yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3576
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Sat 2015-02-07 18:49:29 +0100
message:
  modified Clump::updateProprties to handel non-spherical bodies, added associated example
added:
  examples/test/clump-facet.py
modified:
  core/Clump.cpp
  examples/polyhedra/clump.py


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'core/Clump.cpp'
--- core/Clump.cpp	2014-10-15 06:44:01 +0000
+++ core/Clump.cpp	2015-02-07 17:49:29 +0000
@@ -172,10 +172,11 @@
 							dens = subBody->material->density;
 							const Sphere* sphere = YADE_CAST<Sphere*> (subBody->shape.get());
 							if((x-subBody->state->pos).squaredNorm() < pow(sphere->radius,2)){
-								M += dv;
-								Sg += dv*x;
+								Real m = dens*dv;
+								M += m;
+								Sg += m*x;
 								//inertia I = sum_i( mass_i*dist^2 + I_s) )	//steiners theorem
-								Ig += dv*( x.dot(x)*Matrix3r::Identity()-x*x.transpose())/*dist^2*/+Matrix3r(Vector3r::Constant(dv*pow(dx,2)/6.).asDiagonal())/*I_s/m = d^2: along princial axes of dv; perhaps negligible?*/;
+								Ig += m*( x.dot(x)*Matrix3r::Identity()-x*x.transpose())/*dist^2*/+Matrix3r(Vector3r::Constant(dv*pow(dx,2)/6.).asDiagonal())/*I_s/m = d^2: along princial axes of dv; perhaps negligible?*/;
 								break;
 							}
 						}
@@ -193,9 +194,19 @@
 				State* subState=subBody->state.get();
 				const Sphere* sphere = YADE_CAST<Sphere*> (subBody->shape.get());
 				Real vol = (4./3.)*Mathr::PI*pow(sphere->radius,3.);
-				M+=vol;
-				Sg+=vol*subState->pos;
-				Ig+=Clump::inertiaTensorTranslate(Vector3r::Constant((2/5.)*vol*pow(sphere->radius,2)).asDiagonal(),vol,-1.*subState->pos);
+				Real m = dens*vol;
+				M+=m;
+				Sg+=m*subState->pos;
+				Ig+=Clump::inertiaTensorTranslate(Vector3r::Constant((2/5.)*m*pow(sphere->radius,2)).asDiagonal(),m,-1.*subState->pos);
+			} else { // non-spherical bodies
+				State* subState = subBody->state.get();
+				const Real& m = subState->mass;
+				const Vector3r& inertia = subState->inertia;
+				const Vector3r& pos = subState->pos;
+				const Quaternionr& ori = subState->ori;
+				M += m;
+				Sg += m*pos;
+				Ig += inertiaTensorTranslate(inertiaTensorRotate(inertia.asDiagonal(),ori),m,-pos);
 			}
 		}
 	}
@@ -212,8 +223,8 @@
 	LOG_TRACE("R_g2c=\n"<<R_g2c);
 	// set quaternion from rotation matrix
 	state->ori=Quaternionr(R_g2c); state->ori.normalize();
-	state->inertia=dens*decomposed.eigenvalues();
-	state->mass=dens*M;
+	state->inertia=decomposed.eigenvalues();
+	state->mass=M;
 	
 	// TODO: these might be calculated from members... but complicated... - someone needs that?!
 	state->vel=state->angVel=Vector3r::Zero();

=== modified file 'examples/polyhedra/clump.py'
--- examples/polyhedra/clump.py	2014-10-22 10:13:18 +0000
+++ examples/polyhedra/clump.py	2015-02-07 17:49:29 +0000
@@ -1,20 +1,12 @@
-
-def clumpPolyhedra(polyhedra):
-	clumpId,junk = O.bodies.appendClumped(polyhedra)
-	return
-	clump = O.bodies[clumpId]
-	clump.state.pos = sum((t.state.pos*t.state.mass for t in polyhedra),Vector3.Zero)/sum(t.state.mass for t in polyhedra)
-	print clump.state.pos
-	clump.state.mass = sum(t.state.mass for t in polyhedra)
-	clump.state.inertia = sum((t.state.inertia for t in polyhedra),Vector3.Zero)	
-
 m = PolyhedraMat()
 m.density = 2600 #kg/m^3 
 m.Ks = 20000
-m.Kn = 1E6 #Pa
+m.Kn = 1E9 #Pa
 m.frictionAngle = 0.6 #rad
 O.materials.append(m)
 
+######################################################################
+# basic test
 t1 = polyhedron(((0,0,0),(0,0,1),(0,1,0),(0,1,1),(1,0,0),(1,0,1),(1,1,0),(1,1,1)))
 t2 = polyhedron(((0,0,1),(0,0,2),(0,1,1),(0,1,2),(2,0,1),(2,0,2),(2,1,1),(2,1,2)))
 
@@ -25,8 +17,29 @@
 t2.state.ori = ((1,2,3),1)
 t2.state.pos = (sqrt(2),.5,-sqrt(2))
 
-clumpPolyhedra((t1,t2))
-
+clumpId,polys = O.bodies.appendClumped((t1,t2))
+clump = O.bodies[clumpId]
+s = clump.state
+s2 = t2.state
+# if t1.mass and inertia are set to almost zero, clump properties should equal t2 properties
+print s.pos, s2.pos
+print s.mass, s2.mass
+print s.inertia, s2.inertia
+######################################################################
+
+
+######################################################################
+# something more real
+wire = False
+O.bodies.clear()
+t1 = polyhedron(((0,0,0),(0,0,1),(0,1,0),(0,1,1),(1,0,0),(1,0,1),(1,1,0),(1,1,1)),color=(0,1,0),wire=wire)
+t2 = polyhedron(((0,0,1),(0,0,2),(0,1,1),(0,1,2),(2,0,1),(2,0,2),(2,1,1),(2,1,2)),color=(0,1,0),wire=wire)
+bottom = polyhedron(((-5,-5,0),(5,-5,0),(5,5,0),(-5,5,0),(-5,-5,-1),(5,-5,-1),(5,5,-1),(-5,5,-1)),fixed=True,color=(0,1,1),wire=wire)
+O.bodies.append(bottom)
+clumpId,polys = O.bodies.appendClumped((t1,t2))
+s = O.bodies[clumpId].state
+s.pos = (0,0,5)
+s.ori = Quaternion((1,1,.3),1.5)
 
 O.engines=[
    ForceResetter(),
@@ -34,16 +47,18 @@
    InteractionLoop(
       [Ig2_Polyhedra_Polyhedra_PolyhedraGeom()],
       [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
-      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()]
+      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()],
    ),
-   NewtonIntegrator()
+   NewtonIntegrator(label='newton'),
 ]
 
-O.dt=0.00025
-#O.step()
+O.dt=0.00001
 
 try:
 	from yade import qt
 	qt.View()
 except:
 	pass
+O.step()
+newton.gravity = (0,0,-9.81)
+######################################################################

=== added file 'examples/test/clump-facet.py'
--- examples/test/clump-facet.py	1970-01-01 00:00:00 +0000
+++ examples/test/clump-facet.py	2015-02-07 17:49:29 +0000
@@ -0,0 +1,36 @@
+O.engines=[
+	ForceResetter(),
+	InsertionSortCollider([
+		Bo1_Sphere_Aabb(),
+		Bo1_Facet_Aabb(),
+	]),
+	InteractionLoop(
+		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
+		[Ip2_FrictMat_FrictMat_FrictPhys()],
+		[Law2_ScGeom_FrictPhys_CundallStrack()]
+	),
+	NewtonIntegrator(),
+]
+
+s1 = sphere((2,2,0),1)
+s2 = sphere((-2,2,0),1,fixed=True)
+f1 = facet(((0,0,-2),(0,0,2),(+5,0,0)),wire=False)
+f2 = facet(((0,0,-2),(0,0,2),(-5,0,0)),wire=False)
+# needs to assign (any) nonzero mass and inertia to facets. Total mass and inertia of clump itself can be assigned independently
+f1.state.mass = f2.state.mass = 1
+f2.state.inertia = f2.state.inertia = (1,1,1)
+O.bodies.append((s1,s2))
+clumpId,facetsId = O.bodies.appendClumped((f1,f2))
+s = O.bodies[clumpId].state
+# now we can assign arbitrary inertia of the overall clump
+s.mass = 5000
+s.inertia = 10000*Vector3(1,1,1)
+s.blockedDOFs = 'xyzXY'
+O.forces.addT(clumpId,(0,0,1e1),True)
+
+O.dt = 5e-5
+try:
+	from yade import qt
+	qt.View()
+except:
+	pass