# yade-users team mailing list archive

## Re: [Question #698572]: Is it possible to constrain the motion of a body along a circle

 Thread Previous • Date Previous • Date Next • Thread Next

```Question #698572 on Yade changed:

Rohit John posted a new comment:
Dear Jan,

I performed a simple test on your code. I tried simulating what would
happen a sphere hits the clump "pendulum". Then the angular velocity
about the origin was plotted against time. The system consists of the
clump and the sphere hitting the clump. Since, the external forces are
supposed to act on the hinge, the angular velocity about the hinge
should be conserved. But I see it is not conserved during the collision.
But it is conserved when the clump pendulum is rotating, suggesting your
idea works in that time. It could also be that, there is an error during
the integration of angular momentum. I have pasted the code below.

Kind regards,
Rohit K. John

# -------------------------------------------------------------------------------------------------------------------------------------------------------------------- Python
from yade import  plot

# ---------------------------------------------------------------------------------- bodies
sph1 = sphere((0, -2, 2), .5)
sph2 = sphere((0, 0, 0), .5)
sph3 = sphere((0, 0, 2), .5)

sph3_id = O.bodies.append(sph3)
O.bodies[sph3_id].state.vel = [0,-1,0]

# ---------------------------------------------------------------- clump
clump_id,(sph1_id,sph2_id) = O.bodies.appendClumped((sph1,sph2))
clump = O.bodies[clump_id]
center,radius = Vector3(0,0,0), sqrt(2)

clump.state.refOri = clump.state.ori # not done automatically
clump.state.blockedDOFs = "XYZ" # to force "simple" angular integration

# ---------------------------------------------------------------------------------- engines
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
PyRunner(command = "eccentricAxisEngine()", iterPeriod = 1),
NewtonIntegrator(gravity=(0,0,0), damping=0, label="newton"),
PyRunner(command = "plotter()", iterPeriod = 100)
]
O.dt = .5e-2 * PWaveTimeStep()

# ---------------------------------------------------------------------------------- additional engines
# ---------------------------------------------------------------- constraintClosesPoint
def constraintClosesPoint(pos):
delta = pos - center
direction = delta.normalized()
pos = center + direction * radius
angle = atan2(direction[2],direction[1])
ori = Quaternion()
ori.setFromTwoVectors(Vector3.UnitX,direction)
return pos,ori

# ---------------------------------------------------------------- eccentricAxisEngine
def eccentricAxisEngine():
dt = O.dt
s = clump.state
f = O.forces.f(clump_id)
a = f / s.mass + newton.gravity
pos1 = s.pos + (s.vel+a*dt)*dt
pos0,ori0 = constraintClosesPoint(pos1)
s.vel = (pos0 - s.pos) / dt - a*dt
#
v1 = s.ori * Vector3.UnitX
v2 = ori0 * Vector3.UnitX
ang = v1.cross(v2)
ang = ang / O.dt
s.angVel = ang

# ---------------------------------------------------------------- plotter
plot.plots = {'t':('ax', 'ay', 'az')}

def plotter():
body_ids    = [clump_id, sph3_id]
totalangVel = Vector3([0,0,0])
for i in body_ids:
totalangVel = totalangVel + getAngVelAboutOrigin(i)

plot.addData(t = O.time, ax = totalangVel[0],
ay = totalangVel[1], az = totalangVel[2]
)
plot.plot()

state  = O.bodies[id].state
pos    = state.pos
vel    = state.vel
angMom = state.angMom
mass   = state.mass

total_ang_mom = angMom + mass*pos.cross(vel)

--

```

 Thread Previous • Date Previous • Date Next • Thread Next