yade-users team mailing list archive

[Rohit John <question698572@xxxxxxxxxxxxxxxxxxxxx>]

Question #698572 on Yade changed:
https://answers.launchpad.net/yade/+question/698572

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 *
from yade import *
from yade import plot


# ---------------------------------------------------------------------------------- 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)
]

# ---------------------------------------------------------------------------------- additional engines
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()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.