GitLab's annual major release is around the corner. Along with a lot of new and exciting features, there will be a few breaking changes. Learn more here.

Commit 5fef57e3 authored by Jan Stransky's avatar Jan Stransky
Browse files

modified Clump::updateProprties to handel non-spherical bodies, added associated example

parent 48a0089f
......@@ -172,10 +172,11 @@ void Clump::updateProperties(const shared_ptr<Body>& clumpBody, unsigned int dis
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 @@ void Clump::updateProperties(const shared_ptr<Body>& clumpBody, unsigned int dis
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 @@ void Clump::updateProperties(const shared_ptr<Body>& clumpBody, unsigned int dis
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();
......
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 @@ t1.state.pos = (0,.5,0)
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 @@ O.engines=[
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)
######################################################################
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
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment