I have 4 different polyhedra shapes with vertices given by t1,t2,t3, and t4. I want to perform a gravity deposition of these polyhedra in a box by deposition of 9 random polyhedra shapes. This is the code that i have:
# -*- coding: utf-8 -*-
"""
Created on Thu May 18 15:08:40 2023
@author: Carine
"""
from __future__ import print_function
from yade import plot, polyhedra_utils, qt
import numpy as np
wire = False
######################################################################
# Need to define polyhedral material for Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys() to work
m = PolyhedraMat() #defines polyhedral properties
m.density = 2600 #kg/m^3
m.Ks = 130E3
m.Kn = 160E3 #N/m
m.frictionAngle = 0.54 #tan(31)=0.6, so 31 is 0.54 rad and the input is in rad
O.materials.append(m) #adds material to yade
######################################################################
# Defining polyhedra shapes in terms of their vertices so technically defining vertices, dimaters are 1 mm
t1_coords = (
(0.00033091155911506355, -0.0003242136881048748, -0.00033985697995026614),
(0.00033091155911506355, 0.00038088860029783984, -0.00033985697995026614),
(-0.0003662726180924643, -0.0003242136881048748, -0.00033985697995026614),
(-0.00008498849969921465, -0.0003242136881048748, 0.0003440300800075135),
(0.00033091155911506355, 0.0000956749508999913, 0.0003055457315613278),
(0.00033091155911506355, 0.00038088860029783984, 0.00025920106521738856),
(-0.0003662726180924643, 0.00038088860029783984, 0.00020493952283110147),
(-0.0003662726180924643, -0.0003242136881048748, 0.00032033663438746235)
)
t2_coords = (
(-0.00014309681589806823, -0.0002272260744496193, 0.00038401015387691325),
(0.000358010967088001, -0.0002272260744496193, 0.00038401015387691325),
(0.0006643612513710932, -0.0002272260744496193, -0.00011312310165154156),
(0.00009399535164170355, -0.0002272260744496193, -0.0004453308199780648),
(0.00028983044257111084, 0.00006624457321674196, -0.0002657890485909965),
(-0.0003263752923764808, 0.00009806342306798776, -0.00021177708099387714),
(-0.00014309681589806823, 0.0005961347452425924, 0.00005063802822300939),
(-0.00027303541351821494, 0.0002667732486219123, 0.00018065628098092602),
(0.00037721690590733186, -0.0000964985942102914, 0.0003057083518969001),
(-0.00047649283929694435, -0.0002272260744496193, -0.00011312310165154156)
)
t3_coords = (
(0.0002600238357760947, 0.0002851848128988758, -0.00024701557197894438),
(0.00016012830358156934, 0.0003484895606572876, -0.00018513953373961826),
(-0.00016008801057746095, 0.0003484895606572876, -0.00018513953373961826),
(-0.00026000224348968266, 0.0002851848128988758, -0.00024701557197894438),
(-0.00038471948905778405, 0.0002025288577517229, -0.0000915820081211398),
(-0.00038471948905778405, -0.00020256899286437123, -0.0000915820081211398),
(-0.00023030565260400093, -0.00040365230950484085, 0.0002613542406089658),
(-0.00016008801057746095, -0.000500774842940389, 0.00018311089846328836),
(0.00016012830358156934, -0.000500774842940389, 0.00018311089846328836),
(0.00023034594560810932, -0.00040365230950484085, 0.0002613542406089658),
(0.00038475978206189245, -0.00020256899286437123, -0.0000915820081211398),
(0.00038475978206189245, 0.0002025288577517229, -0.0000915820081211398),
(0.00023034594560810932, 0.00040361217439219253, 0.0002613542406089658),
(0.00016012830358156934, 0.000500734707827741, 0.00018311089846328836),
(-0.00016008801057746095, 0.000500734707827741, 0.00018311089846328836),
(-0.00023030565260400093, 0.00040361217439219253, 0.0002613542406089658),
(-0.00016008801057746095, 0.0002025288577517229, 0.0004286664300572351),
(-0.00016008801057746095, -0.00020256899286437123, 0.0004286664300572351),
(0.00016012830358156934, -0.00020256899286437123, 0.0004286664300572351),
(0.00016012830358156934, 0.0002025288577517229, 0.0004286664300572351),
(0.00016012830358156934, 0.0002025288577517229, -0.0004030812475869826),
(-0.00016008801057746095, 0.0002025288577517229, -0.0004030812475869826),
(-0.00016008801057746095, -0.00020256899286437123, -0.0004030812475869826),
(0.00016012830358156934, -0.00020256899286437123, -0.0004030812475869826),
(0.0002600238357760947, -0.0002852249480115241, -0.00024701557197894438),
(0.00016012830358156934, -0.000348529695769936, -0.00018513953373961826),
(-0.00016008801057746095, -0.000348529695769936, -0.00018513953373961826),
(-0.00026000224348968266, -0.0002852249480115241, -0.00024701557197894438)
)
t4_coords = (
(0.000425451497911189, 0.000410663686410634, 0.000414950888597802),
(0.000568476283629213, 0.000328088287631944, 0.00033237549181911),
(0.000568476283629213, -0.00032837342517883, 0.00033237549181911),
(0.000425451497911189, -0.000410948821957522, 0.000414950888597802),
(-3.62362424167e-08, -0.000656604280584219, -1.05368442558058e-05),
(-0.000406908599377877, -0.00042169641187843, -0.000417409208691264),
(-0.000568548755514047, -0.000328373423178833, -0.000324086220991665),
(-0.000568548755514047, 0.000328088287631944, -0.000324086220991665),
(-0.000406908599377877, 0.000421411275331543, -0.000417409208691264),
(-3.62362424167e-08, 0.000656319144037332, -1.05368442558058e-05),
(0.000568476283629213, 0.000328088287631944, -0.000324086220991665),
(0.000568476283629213, -0.000328373423178833, -0.000324086220991665),
(0.000406836127493042, -0.00042169641187843, -0.000417409208691264),
(0.000126867164394245, -0.000328373423178833, -0.000579049363827436),
(-0.000126939637279078, -0.000328373423178833, -0.000579049363827436),
(-3.62362424167e-08, 3.71246810048347e-05, -0.000652317076397053),
(-0.000126939637279078, 0.000328088287631944, -0.000579049363827436),
(0.000126867164394245, 0.000328088287631944, -0.000579049363827436),
(0.000406836127493042, 0.000421411275331543, -0.000417409208691264),
(0.000177725303575118, 0.000328088287631944, 0.000557975674315826),
(-3.62362424167e-08, 3.71246810048347e-05, 0.000660606347224498),
(0.000177725303575118, -0.000328373423178833, 0.000557975674315826),
(-0.000177797776459951, -0.000328373423178833, 0.000557975674315826),
(-0.000425523970796022, -0.000410948821957522, 0.000414950888597802),
(-0.000568548755514047, -0.000328373423178833, 0.00033237549181911),
(-0.000568548755514047, 0.000328088287631944, 0.00033237549181911),
(-0.000425523970796022, 0.000410663686410634, 0.000414950888597802),
(-0.000177797776459951, 0.000328088287631944, 0.000557975674315826)
)
shapes = [t1_coords, t2_coords, t3_coords, t4_coords]
bottom = polyhedron(((-5E-3,-0.2E-3,25E-3),(50E-3,-0.2E-3,25E-3),(50E-3,0E-3,25E-3),(-5E-3,0,25E-3),(-5E-3,-0.2E-3,-15E-3),(50E-3,-0.2E-3,-15E-3),(50E-3,0,-15E-3),(-5E-3,0,-15E-3)),fixed=True,color=(0.1,0.2,0.5),wire=True,material=m)
leftside = polyhedron(((-0.2E-3,0,10E-3),(0,0,10E-3),(0,8E-3,10E-3),(-0.E-3,8E-3,10E-3),(-0.2E-3,0,0),(0,0,0),(0,8E-3,0),(-0.2E-3,8E-3,0)),fixed=True,color=(0.1,0.2,0.5),wire=True,material=m)
rightside = polyhedron(((10E-3,0,10E-3),(10.2E-3,0,10E-3),(10.2E-3,8E-3,10E-3),(10E-3,8E-3,10E-3),(10E-3,0,0),(10.2E-3,0,0),(10.2E-3,8E-3,0),(10E-3,8E-3,0)),fixed=True,color=(0.1,0.2,0.5),wire=True,material=m) #side that will be removed
back = polyhedron(((0,0,10.2E-3),(10E-3,0,10.2E-3),(10E-3,8E-3,10.2E-3),(0,8E-3,10.2E-3),(0,0,10E-3),(10E-3,0,10E-3),(10E-3,8E-3,10E-3),(0,8E-3,10E-3)),fixed=True,color=(0.1,0.2,0.5),wire=True,material=m)
front = polyhedron(((0,0,0),(10E-3,0,0),(10E-3,8E-3,0),(0,8E-3,0),(0,0,-0.2E-3),(10E-3,0,-0.2E-3),(10E-3,8E-3,-0.2E-3),(0,8E-3,-0.2E-3)),fixed=True,color=(0.1,0.2,0.5),wire=True,material=m)
######################################################################
# Setting polyhedra positions in terms of their centers of mass
part_list = [] # Creates an empty list
i = 0 # counter
for x in np.linspace(1.5E-3, 8.5E-3, num=3): # (min coordinate, max coordinate, number of points) box width is 10 mm
for y in np.linspace(10E-3,14E-3, num=1): # box height is 8 mm and we are dropping into the box so drop from a heigher height
for z in np.linspace(1.5E-3, 8.5E-3, num=1): # width is 10 mm
shape_id = np.random.randint(0,4) # chooses a number between zero and 3
shape = shapes[shape_id] # uses the random number to pick one of the 4 shapes
part_list += [polyhedron(shape,color=(1,1,1),wire=wire,material=m)] # adds a new particle to the particle list
part_list[i].state.pos = (x,y,z)
i = i + 1
# centroid location
bottom.state.pos = (22.5E-3,-0.1E-3,5E-3)
leftside.state.pos = (-0.1E-3,4.1E-3,5E-3)
rightside.state.pos = (10.1E-3,4.1E-3,5E-3)
back.state.pos = (5E-3,4.1E-3,10.2E-3)
front.state.pos = (5E-3,4.1E-3,-0.2E-3)
part_list = part_list + [bottom] + [leftside] + [back] + [front]
######################################################################
# Adding the bodies to the simulation, to see how to assign material properties to different polyhedra check O.bodies.append used in free-fall.py
O.bodies.append(part_list)
ThisID = O.bodies.append(rightside)
######################################################################
# Simulation loop
O.engines=[
# Reset forces
ForceResetter(),
# Appropriate collision detection, create interactions
InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb()]),
# Handle interaction
InteractionLoop(
[Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()] # contact law -- apply forces
),
# GravityEngine(gravity=(0,0,-9.81)), and Newton Integrator which updates position using Newton's equation
NewtonIntegrator(damping=0.10,gravity=(0,-9.810,0)),
# save data from Yade's own 3d view
qt.SnapshotEngine(fileBase='pic', iterPeriod=10000, label='snapshot'),
#VTKRecorder(fileName='./vtk/small',label='vtkRecorder',iterPeriod=10000),
# this engine will be called after 20000 steps, only once
PyRunner(command='finish()', iterPeriod=500000)
]
######################################################################
# Time step
O.dt = 4.03E-7
from yade import qt
qt.Controller()
v=qt.View()
O.saveTmp()
Is there a more efficient way of writing it?
