Yade

How to initialize the angular velocity of a box made of pfacets

Asked by Rohit John

Hello,

I created a box using pFacets. Now I want to initialise its angular velocity. I did so by assigning the same angular velocity to all bodies in the pfacet box(nodes, cylinders and pfacet) using O.bodies[i].state.angVel, and I also assigned the appropriate linear velocities using v = r x w, where v is the velocity, r is the radius vector and w is the angular velocity. Gravity is 0 and there is no other interaction in the simulation. I ran the code and saw the box rotating in the viewport without any erratic motion. However, it is behaving differently when I plot angVel for a body in the pfacet box. It was not constant. I even tried averaging this value for all bodies in the pfacet box. It did not work. The code for this case is given below as Code A.

To check if the angVel property has an issue I tried the following. I set two spheres, positioned close to two opposite edges and moving towards the box in opposite direction, to collide with the box. This did the trick and made the box rotate. When I plotted angVel it shows a constant (almost) value. The code for this case is given below as Code B.

My question is there a way to set the angular velocity. I do not want to use angMom as the moment of inertia for all bodies in the pfacet box are different and I guess some objects are planar so they have 0 moments of inertia.

############################################################################################################## Code A

from yade import *
from yade import plot
from yade.gridpfacet import *

pfacet_rad = 0.005

young = 1e6
density = 1e3
frictAng = radians(30)
poisson = 0.3

# function --------------------------------------------------------------------------------------
## box creator function -----------------------------------------------------------
def symmetric_cuboid_pfacet(CM, length, breadth, height, radius, int_mat, ext_mat, color = [0.5, 0.5, 0.5]):
    '''
    makes a symmetric cuboid using pfacet.
    '''
    l = length
    b = breadth
    h = height
    cm1, cm2, cm3 = CM # the centre of mass of the cube
    # coordinates of the vertices
    vertices = [
        [cm1 + l/2, cm2 + b/2, cm3 + h/2],
        [cm1 + l/2, cm2 - b/2, cm3 + h/2],
        [cm1 - l/2, cm2 - b/2, cm3 + h/2],
        [cm1 - l/2, cm2 + b/2, cm3 + h/2],
        [cm1 + l/2, cm2 + b/2, cm3 - h/2],
        [cm1 + l/2, cm2 - b/2, cm3 - h/2],
        [cm1 - l/2, cm2 - b/2, cm3 - h/2],
        [cm1 - l/2, cm2 + b/2, cm3 - h/2],
        [cm1 + l/2, cm2 , cm3 ],
        [cm1 - l/2, cm2 , cm3 ],
        [cm1 , cm2 + b/2, cm3 ],
        [cm1 , cm2 - b/2, cm3 ],
        [cm1 , cm2 , cm3 + h/2],
        [cm1 , cm2 , cm3 - h/2],
    ]

    # connecting the vertices to form faces
    face = [
        [13, 2, 1], [13, 3, 2], [13, 4, 3], [13, 1, 4],
        [14, 5, 6], [14, 6, 7], [14, 7, 8], [14, 8, 5],
        [11, 1, 5], [11, 5, 8], [11, 8, 4], [11, 4, 1],
        [12, 2, 3], [12, 3, 7], [12, 7, 6], [12, 6, 2],
        [9, 1, 2], [9, 2, 6], [9, 6, 5], [9, 5, 1],
        [10, 3, 4], [10, 4, 8], [10, 8, 7], [10, 7, 3],
    ]

    # generating the grid stuff
    nodesIds = []
    cylIds = []
    pFacet = []

    for i in face:
        pfacetCreator1(
            [
                vertices[i[0] - 1],
                vertices[i[1] - 1],
                vertices[i[2] - 1]
            ],
            radius,
            nodesIds = nodesIds,
            cylIds = cylIds,
            pfIds = pFacet,
            wire = False,
            fixed = False,
            color = color,
            materialNodes = int_mat,
            material = ext_mat,
            )
    return nodesIds ,cylIds, pFacet

## cross product function ---------------------------------------------------------
def cross_prd(a,b):
    '''
    Returns the cross product of the vectors a and b
    '''
    a1,a2,a3 = a
    b1,b2,b3 = b

    c1 = a2*b3 - a3*b2
    c2 = a3*b1 - a1*b3
    c3 = a1*b2 - a2*b1

    c = [c1 ,c2, c3]
    return c

## vector difference function -----------------------------------------------------
def subtract(a,b):
    '''
    returns the difference of the vectors a and b
    '''
    a1,a2,a3 = a
    b1,b2,b3 = b

    c = [a1-b1, a2-b2, a3-b3]
    return c

## add an angular velocity to a list of bodies function ---------------------------
def body_angular_velocity_add(ang_vel, body_list, CM):
    '''
    Adds a value to the the initial angular velocity (ang_vel) of the pfacet body (body_list) about
    the centre point CM
    '''
    for i in body_list:
        pos = O.bodies[i].state.pos
        r = subtract(pos, CM)
        vel = cross_prd(r, ang_vel)
        O.bodies[i].state.vel = O.bodies[i].state.vel + vel
        O.bodies[i].state.angVel = O.bodies[i].state.angVel + ang_vel

# material ------------------------------------------------------------------------------------------
## pfacet internal material ------------------------------------------------------
O.materials.append(
    FrictMat(
        young = young,
        poisson = poisson,
        density = density,
        frictionAngle = frictAng,
        label = 'pFacet_out',
    )
)

## pfacet external material ------------------------------------------------------
O.materials.append(
    CohFrictMat(
        young = young,
        poisson = poisson,
        density = density,
        frictionAngle = frictAng,
        momentRotationLaw = True,
        normalCohesion = 1e40,
        shearCohesion = 1e40,
        label = 'pFacet_in'
    )
)

# engines -------------------------------------------------------------------------------------------
O.engines = [
    ForceResetter(),
    InsertionSortCollider([
        Bo1_GridConnection_Aabb(),
        Bo1_PFacet_Aabb()
    ]),
    InteractionLoop(
        [
            Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
            Ig2_GridNode_GridNode_GridNodeGeom6D(),
            Ig2_GridConnection_PFacet_ScGeom(),
            Ig2_PFacet_PFacet_ScGeom()
        ],
        [
            Ip2_FrictMat_FrictMat_FrictPhys(),
            Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow = True, setCohesionOnNewContacts = False)
        ],
        [
            Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
            Law2_ScGeom_FrictPhys_CundallStrack(),
            Law2_ScGeom6D_CohFrictPhys_CohesionMoment()
        ],
    ),
    NewtonIntegrator(gravity = (0,0,0),damping = 0.0),
    PyRunner(command='addPlot()',iterPeriod=10)
]
O.dt = 1e-6

# bodies --------------------------------------------------------------------------------------------
## box ----------------------------------------------------------------------------
(
    node_ids,
    cyl_ids,
    pFacet_ids
) = symmetric_cuboid_pfacet(
        [0,0,0],
        length =0.05,
        breadth =0.05,
        height =0.05,
        radius =pfacet_rad,
        int_mat = 'pFacet_in',
        ext_mat = 'pFacet_out',
        color = [0.5, 0.5, 0.5]
    )

## defining the velocity ----------------------------------------------------------
body_angular_velocity_add(
    [10,0,0],
    node_ids + cyl_ids + pFacet_ids,
    [0,0,0]
    )

# plots ---------------------------------------------------------------------------------------------
plot.plots = {'t':'wx'}

def angvelx():
    '''Returns the average angular velocity along the x axis'''
    wx = 0
    for i in (node_ids + cyl_ids + pFacet_ids):
        wx = wx + O.bodies[i].state.angVel[0]

    n = len(node_ids + cyl_ids + pFacet_ids)
    return wx/n

def addPlot():
    plot.addData(
        t = O.time,
        wx = angvelx(),
        )

plot.plot()
############################################################################################################## End of code A

############################################################################################################## Code B
# The code for box rotated by impacting it with spheres
# Code B

# new added plots
from yade import *
from yade import plot
from yade.gridpfacet import *

# set up sims ---------------------------------------------------------------------------------------

pfacet_rad = 0.005

young = 1e9
density = 1e3
frictAng = radians(30)
poisson = 0.3

# function --------------------------------------------------------------------------------------
## box creator function -----------------------------------------------------------
def symmetric_cuboid_pfacet(CM, length, breadth, height, radius, int_mat, ext_mat, color = [0.5, 0.5, 0.5]):
    '''
    makes a symmetric cuboid using pfacet.
    '''
    l = length
    b = breadth
    h = height
    cm1, cm2, cm3 = CM # the centre of mass of the cube
    # coordinates of the vertices
    vertices = [
        [cm1 + l/2, cm2 + b/2, cm3 + h/2],
        [cm1 + l/2, cm2 - b/2, cm3 + h/2],
        [cm1 - l/2, cm2 - b/2, cm3 + h/2],
        [cm1 - l/2, cm2 + b/2, cm3 + h/2],
        [cm1 + l/2, cm2 + b/2, cm3 - h/2],
        [cm1 + l/2, cm2 - b/2, cm3 - h/2],
        [cm1 - l/2, cm2 - b/2, cm3 - h/2],
        [cm1 - l/2, cm2 + b/2, cm3 - h/2],
        [cm1 + l/2, cm2 , cm3 ],
        [cm1 - l/2, cm2 , cm3 ],
        [cm1 , cm2 + b/2, cm3 ],
        [cm1 , cm2 - b/2, cm3 ],
        [cm1 , cm2 , cm3 + h/2],
        [cm1 , cm2 , cm3 - h/2],
    ]

    # connecting the vertices to form faces
    face = [
        [13, 2, 1], [13, 3, 2], [13, 4, 3], [13, 1, 4],
        [14, 5, 6], [14, 6, 7], [14, 7, 8], [14, 8, 5],
        [11, 1, 5], [11, 5, 8], [11, 8, 4], [11, 4, 1],
        [12, 2, 3], [12, 3, 7], [12, 7, 6], [12, 6, 2],
        [9, 1, 2], [9, 2, 6], [9, 6, 5], [9, 5, 1],
        [10, 3, 4], [10, 4, 8], [10, 8, 7], [10, 7, 3],
    ]

    # generating the grid stuff
    nodesIds = []
    cylIds = []
    pFacet = []

    for i in face:
        pfacetCreator1(
            [
                vertices[i[0] - 1],
                vertices[i[1] - 1],
                vertices[i[2] - 1]
            ],
            radius,
            nodesIds = nodesIds,
            cylIds = cylIds,
            pfIds = pFacet,
            wire = False,
            fixed = False,
            color = color,
            materialNodes = int_mat,
            material = ext_mat,
            )
    return nodesIds ,cylIds, pFacet

# material ------------------------------------------------------------------------------------------
## pfacet internal material ------------------------------------------------------
O.materials.append(
    FrictMat(
        young = young,
        poisson = poisson,
        density = density,
        frictionAngle = frictAng,
        label = 'pFacet_out',
    )
)

## pfacet external material ------------------------------------------------------
O.materials.append(
    CohFrictMat(
        young = 1000*young,
        poisson = poisson,
        density = density,
        frictionAngle = frictAng,
        momentRotationLaw = True,
        normalCohesion = 1e40,
        shearCohesion = 1e40,
        label = 'pFacet_in'
    )
)

## sphere material ---------------------------------------------------------------
O.materials.append(
    FrictMat(
        young = young,
        poisson = poisson,
        density = 10000*density,
        frictionAngle = frictAng,
        label = 'sphere',
    )
)
# engines -------------------------------------------------------------------------------------------
O.engines = [
    ForceResetter(),
    InsertionSortCollider([
        Bo1_Sphere_Aabb(),
        Bo1_GridConnection_Aabb(),
        Bo1_PFacet_Aabb(),
    ]),
    InteractionLoop(
        [
            Ig2_Sphere_Sphere_ScGeom(),
            Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
            Ig2_GridNode_GridNode_GridNodeGeom6D(),
            Ig2_GridConnection_PFacet_ScGeom(),
            Ig2_PFacet_PFacet_ScGeom(),
            Ig2_Sphere_PFacet_ScGridCoGeom(),
            Ig2_Sphere_GridConnection_ScGridCoGeom(),
        ],
        [
            Ip2_FrictMat_FrictMat_FrictPhys(),
            Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow = True, setCohesionOnNewContacts = False)
        ],
        [
            Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
            Law2_ScGeom_FrictPhys_CundallStrack(),
            Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
            Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
        ],
    ),
    NewtonIntegrator(gravity = (0,0,0),damping = 0.0),
    PyRunner(command='addPlot()',iterPeriod=10)
]

# bodies --------------------------------------------------------------------------------------------
## box ----------------------------------------------------------------------------
(
    node_ids,
    cyl_ids,
    pFacet_ids
) = symmetric_cuboid_pfacet(
        [0,0,0],
        length =0.05,
        breadth =0.05,
        height =0.05,
        radius =pfacet_rad,
        int_mat = 'pFacet_in',
        ext_mat = 'pFacet_out',
        color = [0.5, 0.5, 0.5]
    )

## box ----------------------------------------------------------------------------
## spheres ------------------------------------------------------------------------
s3 = utils.sphere([0, 0.04, 0.025],0.005,material='sphere')
s4 = utils.sphere([0,-0.04,-0.025],0.005,material='sphere')

s3id = O.bodies.append(s3)
s4id = O.bodies.append(s4)

O.bodies[s3id].state.vel = [0,-1,0]
O.bodies[s4id].state.vel = [0,1,0]

O.saveTmp()

# plots ---------------------------------------------------------------------------------------------
plot.plots = {'t':'wx'}

def angvelx():
    '''Returns the average angular velocity along the x axis'''
    wx = 0
    for i in (node_ids + cyl_ids + pFacet_ids):
        wx = wx + O.bodies[i].state.angVel[0]

    n = len(node_ids + cyl_ids + pFacet_ids)
    return wx/n

def addPlot():
    plot.addData(
        t = O.time,
        wx = angvelx(),
        )

plot.plot()

O.dt = utils.PWaveTimeStep()
O.saveTmp()

Question information

Language:
English Edit question
Status:
Expired
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jérôme Duriez (jduriez) said : 		#1

Hi,

Maybe your reflections can be fed from https://gitlab.com/yade-dev/trunk/-/issues/84

Revision history for this message
Launchpad Janitor (janitor) said : 		#2

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Revision history for this message
Rohit John (rohitkjohn) said : 		#3

Hello,

Sorry for the late reply. I'll see what I can do

Revision history for this message
Rohit John (rohitkjohn) said : 		#4

Hello,

I found a work around to this problem. I used the rotationEngine to initialise the angular velocity. It is active for the first few iterations and then it will be turned off using a pyrunner command.

Kind regards,
RKJ

# -------------------------------------------------------------------------------------------------------------------------------------------------- script
O.engines += [
            RotationEngine(
                rotateAroundZero = True,
                zeroPoint = (0,0,0),
                rotationAxis = axis,
                angularVelocity = -15,
                ids = target.ids,
                label = 'rotEngine'
            ),
            PyRunner(command = "switch_of_rot()", iterPeriod = 10, label = "rot_switch"),
]

def switch_of_rot():
    rotEngine.dead = True
    print("engine off")
    rot_switch.iterPeriod = 100000000
# --------------------------------------------------------------------------------------------------------------------------------------------------