# yade-users team mailing list archive

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

```Question #698572 on Yade changed:

Rohit John posted a new comment:
Hello all,

I  managed to make the code I had posted earlier faster (#6) by getting
an analytical solution for the non linear equations. The code is given
below. It affects the linear motion. I believe the angular motion will
require more complicated equations

Kind regards,
Rohit K. John

from math import *

# ---------------------------------------------------------------------------------- bodies
sph1 = sphere((2, 0, 2), .5)
sph2 = sphere(center=(0, 0, 0), radius=.5, fixed=True)

sph1_id = O.bodies.append(sph1)
sph2_id = O.bodies.append(sph2)

g      =  Vector3([0,0,-10])
mass_1 =  O.bodies[sph1_id].state.mass
length =  O.bodies[sph1_id].state.pos.norm()
# ---------------------------------------------------------------------------------- engines
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],  # collision geometry
[Ip2_FrictMat_FrictMat_FrictPhys()],  # collision "physics"
[Law2_ScGeom_FrictPhys_CundallStrack()]  # contact law -- apply forces
),
PyRunner(command = "constrainingForceEngine2()", iterPeriod = 1),
NewtonIntegrator(gravity=g, damping=0.0)
]

plot.plots = {"t": ("fx", "fy", "fz")}

def constrainingForceEngine2():
m  = O.bodies[sph1_id].state.mass
dt = O.dt
l  = length
x0,  y0,  z0   = O.bodies[sph1_id].state.pos
vx0, vy0, vz0  = O.bodies[sph1_id].state.vel
F0x, F0y, F0z  = O.forces.f(sph1_id) +  m*g

a = dt**2*F0x + dt*m*vx0 + m*x0
b = dt**2*F0z + dt*m*vz0 + m*z0

z = b*l/(sqrt(b**2+a**2))
x = a*l/(sqrt(b**2+a**2))

vx = (x  - x0)/dt
vz = (z  - z0)/dt

if x !=0:
lamda = (m*(vx - vx0)/dt - F0x)/x
else:
lamda = (m*(vz - vz0)/dt - F0z)/z

# Assiging the force
Fx = lamda * x
Fz = lamda * z
force = Vector3([Fx, 0, Fz])
O.forces.addF(id = sph1_id, f = force)

if O.iter % 1000 == 0:
plot.addData(t = O.time, fx = force[0], fy = force[1], fz = force[2])

plot.plot()
# ---------------------------------------------------------------------------------- sim setup
O.dt = 1e-5
O.saveTmp()

--