# 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

young = 1e6
density = 1e3
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
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 - 1],
vertices[i - 1],
vertices[i - 1]
],
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 ---------------------------
'''
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),
]
O.dt = 1e-6

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

## defining the velocity ----------------------------------------------------------
[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

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

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

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

young = 1e9
density = 1e3
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
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 - 1],
vertices[i - 1],
vertices[i - 1]
],
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),
]

# bodies --------------------------------------------------------------------------------------------
## box ----------------------------------------------------------------------------
(
node_ids,
cyl_ids,
pFacet_ids
) = symmetric_cuboid_pfacet(
[0,0,0],
length =0.05,
height =0.05,
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

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

t = O.time,
wx = angvelx(),
)

plot.plot()

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

## Question information

Language:
English Edit question
Status:
Open
For:
Assignee:
No assignee Edit question
Last query:
2020-10-15