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

Asked by Rohit John on 2020-10-15

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:
Open
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-10-15
Last reply:
Jérôme Duriez (jduriez) said : #1

Hi,

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

Can you help with this problem?

Provide an answer of your own, or ask Rohit John for more information if necessary.

To post a message you must log in.