# yade-users team mailing list archive

## Re: [Question #705718]: Error in the calculation of polyhedral particles

```Question #705718 on Yade changed:

The code for this polyhedron model is as follows:
******************************************************************************
import numpy
from pprint import pprint
import random
from random import uniform
from random import randint
from math import *

matP = PolyhedraMat() #the material of wall
matP.density = 7800 #kg/m^3
matP.young = 5.5E8 #10 time biger
matP.poisson = 0.21
matP.fionAngle = 0.03 #0.65
O.materials.append(matP)

steel = PolyhedraMat() #the material of plate
steel.density = 2000 #kg/m^3  in order to ignore the gravity of plate
steel.young = 1E9    #inital steel was ,try
steel.poisson = 0.21
O.materials.append(steel)

global a
a=1
#length unit: m; using stander unit

O.bodies.append(utils.wall((0,0,0),0,sense=1,material=steel))
O.bodies.append(utils.wall((0.4,0,0),0,sense=-1,material=steel))
O.bodies.append(utils.wall((0,0,0),1,sense=1,material=steel))
O.bodies.append(utils.wall((0,0.4,0),1,sense=-1,material=steel))
O.bodies.append(utils.wall((0,0,0),2,sense=1,material=steel))

'''sp=pack.SpherePack()
sp.makeCloud(minCorner=(-0.19,-0.24,0),maxCorner=(0.19,0.24,0.385),porosity=0.3,psdSizes=[0.0025,0.003,0.00395,0.00535,0.00805,0.0105],psdCumm=[0,0.0939,0.3012,0.5693,0.8842,1])
sp.toSimulation(material=spheresand)'''

ldplt1=polyhedra_utils.polyhedra(matP,v=((0.0311388231,0.269696415,0.0589928355),(0.0290950595,0.233212497,0.0434321398),(0.00477624027,0.214684841,0.0612635859),(0.000413574885,0.250807426,0.0295583737),(0.0173971545,0.25909924,0.0409667532),(0.0122290971,0.257075205,0.0368187339),(0.0185846097,0.218305976,0.0649987468),(0.00289712704,0.218589488,0.0468812254),(-0.000370547223,0.245767932,0.0267756556),(0.0226644807,0.219353472,0.0641023017),(0.0294671189,0.259766604,0.0720264969),(0.0262708557,0.253351576,0.075764224),(0.0201327597,0.265420489,0.0688631612),(0.00371496562,0.253315814,0.0610497987),(0.00134598311,0.217912015,0.0534571812),(0.00201169851,0.229893552,0.0191997213),(0.00739037819,0.226917118,0.0279797141),(0.0250662148,0.229144799,0.0378462333),(-0.00184815953,0.24418841,0.048189802),(0.0318022057,0.250408569,0.0688592611)),fixed=False, color=(0.75,0.65,0.65))
ldplt2=polyhedra_utils.polyhedra(matP,v=((0.200802937,0.00522635087,0.0764616202),(0.190203716,0.031931028,0.0616703076),(0.2163187,0.0177309575,0.0372346987),(0.199799386,0.0213384504,0.0287916627),(0.190185059,0.00565917637,0.0537674456),(0.190391663,0.00343076801,0.0752545319),(0.195803192,0.00102240847,0.0768751083),(0.195824587,0.0340955688,0.0420706728),(0.196550561,0.0193935442,0.0292134554),(0.185818083,0.0116228699,0.05176113),(0.207639423,-0.00161284012,0.0599384523),(0.211256829,-0.00130980622,0.048370807),(0.198705631,0.00420700472,0.0345769575),(0.212568236,0.0224350064,0.0613230691),(0.205961943,0.0155357745,0.0692370933),(0.191026313,0.0328246731,0.0643383676),(0.19725945,0.0266500533,0.0320846928),(0.212193909,0.00712049872,0.0514893316)),fixed=False, color=(0.75,0.65,0.65))
ldplt3=polyhedra_utils.polyhedra(matP,v=((0.0493890096,0.124625964,0.0546219955),(0.0587884523,0.101588876,0.0695327281),(0.0484898696,0.131316964,0.0671963882),(0.0851228942,0.120919536,0.0573203433),(0.0857616542,0.121410806,0.0540722961),(0.0838209627,0.116238451,0.0584039996),(0.0635240328,0.134572434,0.0328978821),(0.0662484768,0.139334204,0.0349402573),(0.0550392553,0.110181804,0.0785459614),(0.068213206,0.145236108,0.0469616925),(0.0722976776,0.143562303,0.0543968559),(0.0709269664,0.142187157,0.0668189968),(0.0611672717,0.109788495,0.0805646999),(0.0665383623,0.128782059,0.0747798246),(0.0682659455,0.114794297,0.0785211583),(0.0782640788,0.14107141,0.0605324682),(0.0748903004,0.135299948,0.0366608433),(0.0581507097,0.122073326,0.0395935881),(0.0906345309,0.128910642,0.0510022243),(0.0526194934,0.146986529,0.0615664254)),fixed=False, color=(0.75,0.65,0.65))
ldplt4=polyhedra_utils.polyhedra(matP,v=((0.269463605,0.1166371,0.0696753702),(0.259285745,0.105563545,0.0748045457),(0.251590098,0.123906883,0.0738906132),(0.236829123,0.129765341,0.0751196651),(0.247482891,0.0983998837,0.0304707525),(0.24774074,0.114695293,0.0383968612),(0.253536036,0.0795080872,0.0660853945),(0.218754376,0.093727008,0.0491631261),(0.22631638,0.0916926354,0.0704584899),(0.256512303,0.0816757632,0.0638146597),(0.258174788,0.0834932918,0.0683190929),(0.230658721,0.0761771663,0.0600711664),(0.208143222,0.126729508,0.043062315),(0.237293675,0.134463798,0.051654522),(0.241533646,0.0985797415,0.02487255),(0.281283116,0.109748354,0.0686402353),(0.213612202,0.106546334,0.0238623595),(0.207692049,0.131440181,0.0697141042),(0.217493708,0.10010838,0.0172017708)),fixed=False, color=(0.75,0.65,0.65))
ldplt5=polyhedra_utils.polyhedra(matP,v=((0.28754689,0.0542172824,0.0533610605),(0.306517053,0.06503258,0.0660064122),(0.294050907,0.0786569084,0.0519439609),(0.291634911,0.0838958271,0.0689598517),(0.272110282,0.0964089265,0.0474494643),(0.258755888,0.075240861,0.0523888492),(0.26014563,0.0582147083,0.0611699829),(0.291533876,0.0559052608,0.0578523851),(0.294005589,0.0627079764,0.0762867104),(0.305845029,0.064352088,0.0692291156),(0.273560777,0.0859867354,0.0711012417),(0.27200557,0.0815395225,0.0711780773),(0.279051003,0.0682565797,0.0778143146),(0.26416247,0.0619951397,0.0700695129),(0.279619839,0.0645552734,0.0759031198),(0.277536097,0.0962531764,0.043524177),(0.28170433,0.066615784,0.0449160177),(0.270705374,0.0918942969,0.0628764745),(0.288970399,0.0913466316,0.0610919568),(0.284065754,0.0647928491,0.0446529477),(0.260176223,0.0589487379,0.0583509404)),fixed=False, color=(0.75,0.65,0.65))
ldplt6=polyhedra_utils.polyhedra(matP,v=((0.0745047596,0.0183132421,0.0745571534),(0.0757254677,0.0181323275,0.0656490742),(0.0719328724,0.0114807741,0.0667747496),(0.0962019599,0.0161924747,0.0763739491),(0.0956007871,0.0103881815,0.0185330401),(0.109248517,0.0151778673,0.0152455953),(0.131042408,0.00726979326,0.0226276221),(0.11142901,-0.000251454602,0.0471690447),(0.11140294,-0.000788525913,0.0152358264),(0.121745495,0.00175210175,0.0124698167),(0.0998697356,0.0198069098,0.043993984),(0.069887293,0.0114109537,0.0413501293),(0.100513493,0.0175575596,0.074982852),(0.137674037,0.0201132217,0.0407585808),(0.107912936,-0.0000845798066,0.0521076373),(0.0717400859,0.000173910171,0.0314216704),(0.104058093,0.00101617836,0.0562555395),(0.139062371,0.0257933696,0.0414929237),(0.0940636122,0.00160849093,0.01615708),(0.106413253,0.0207241771,0.0688687268),(0.121563594,0.0208726091,0.0604000982)),fixed=False, color=(0.75,0.65,0.65))
ldplt7=polyhedra_utils.polyhedra(matP,v=((0.138681617,0.0709185227,0.0477621616),(0.124873262,0.0716532602,0.053608597),(0.15984327,0.0696961737,0.039230532),(0.160042149,0.0705094535,0.073073618),(0.198519816,0.0688201346,0.0556083365),(0.197752547,0.0871154535,0.0631642401),(0.183764439,0.107659447,0.0600438292),(0.186365372,0.0904045959,0.043360441),(0.194700521,0.0931648202,0.0663408142),(0.175319464,0.0657838402,0.0738609809),(0.187900689,0.0658283632,0.0481113877),(0.178458864,0.0704029275,0.0737919421),(0.204186586,0.0745222514,0.0722512729),(0.149182131,0.0874655418,0.0496051205),(0.180233883,0.108216102,0.0616678881),(0.130970361,0.0868178474,0.0532821847),(0.166580293,0.075618052,0.0378884199),(0.187937838,0.0834325988,0.0451001152),(0.177201475,0.0728739105,0.0384079861),(0.150920279,0.0949583329,0.0721540788)),fixed=False, color=(0.75,0.65,0.65))
ldplt8=polyhedra_utils.polyhedra(matP,v=((0.233399802,0.0435860937,0.052979262),(0.264677313,0.0507434894,0.0576122567),(0.259249504,0.0360805818,0.0305850682),(0.272485932,0.00296102355,0.0586477276),(0.271815712,0.0247014829,0.0309033309),(0.280258756,0.0398891032,0.0342213258),(0.23728288,0.0267067305,0.0762665959),(0.285362461,0.0284282906,0.0668103551),(0.231804952,0.0450543014,0.0570982701),(0.284392976,0.0137217487,0.0509834975),(0.289927233,0.00886335669,0.066557277),(0.285477811,0.0447512465,0.0593598616),(0.287218601,0.0339741525,0.0639914229),(0.255850049,0.0494555198,0.0594322917),(0.245291136,0.0123194777,0.0638477339),(0.254606381,-0.000785586939,0.0550352996),(0.252348928,0.00913746133,0.0490902296),(0.243556011,0.0355925646,0.0381974012),(0.285226122,0.00537147876,0.06580441)),fixed=False, color=(0.75,0.65,0.65))
ldplt9=polyhedra_utils.polyhedra(matP,v=((0.186408192,0.381826769,0.054815705),(0.168320512,0.396438463,0.0662955201),(0.188858887,0.374484895,0.0484244662),(0.169180328,0.383001389,0.0402561882),(0.18379932,0.387179458,0.0576007192),(0.172796201,0.400416415,0.0702998412),(0.177385518,0.373502408,0.030860205),(0.170410292,0.385375072,0.0635008061),(0.169411147,0.401510044,0.0711969217),(0.171920441,0.379900989,0.0610478942),(0.168292705,0.387684426,0.0435563568),(0.167664655,0.392011561,0.046027939),(0.172111028,0.371519508,0.0447943719),(0.172712281,0.370810271,0.0551308672),(0.18424394,0.375764834,0.0297553789),(0.189314681,0.374493039,0.0367433235),(0.181934706,0.395703787,0.0454660315),(0.170833253,0.372265533,0.049719743),(0.176350729,0.39624978,0.0447478883)),fixed=False, color=(0.75,0.65,0.65))
ldplt10=polyhedra_utils.polyhedra(matP,v=((0.0491463145,0.117072334,0.0664846326),(0.0462700699,0.122547202,0.0708346284),(0.0562051398,0.125717437,0.0426587167),(0.0474298193,0.117842399,0.0690315459),(0.0437157733,0.153572433,0.0633776168),(0.0448357952,0.13765534,0.0288864818),(0.0495533487,0.122075633,0.0269308703),(0.0507438469,0.119614994,0.0276565858),(0.036355923,0.153756113,0.0637762867),(0.0345060811,0.154088614,0.0471305664),(0.0468431374,0.107149963,0.0484263818),(0.0541651783,0.130510221,0.0456532818),(0.0407849899,0.141518154,0.0303559941),(0.0392811001,0.102065239,0.0339441115),(0.0448554108,0.149563669,0.0690864681),(0.0280120418,0.132804301,0.0475486932),(0.0271169116,0.14283171,0.0469772915),(0.032357532,0.11260055,0.043667408),(0.0306548108,0.124505846,0.0526176666),(0.0389902382,0.15416806,0.0674117973)),fixed=False, color=(0.75,0.65,0.65))
ldplt11=polyhedra_utils.polyhedra(matP,v=((0.0556900364,0.0909826808,0.0335290648),(0.06526648,0.0732225817,0.0264800125),(0.0591591689,0.0763219399,0.036541556),(0.0624996208,0.0764271955,0.0350806568),(0.0524131462,0.100131633,0.0221138252),(0.0580709957,0.0846045551,0.00316503935),(0.0748827804,0.103643832,0.00616873221),(0.0744504947,0.0901098979,0.000662409026),(0.0768618882,0.0838263669,0.0119797067),(0.054031248,0.109168468,0.0036971096),(0.0508836611,0.112886338,0.00720795501),(0.0511782701,0.0895539614,0.0123859429),(0.0724653024,0.0896695782,-0.000884177804),(0.0668611012,0.0731873145,0.0241323173),(0.0633945878,0.114522711,0.0139514738),(0.0620231837,0.108407498,0.0215308867),(0.0723167845,0.102558491,0.0149998265),(0.0570278753,0.0997788173,-0.00405829766),(0.0736717278,0.0906080925,0.0186850455)),fixed=False, color=(0.75,0.65,0.65))
ldplt12=polyhedra_utils.polyhedra(matP,v=((0.222452155,0.317655366,0.0749856907),(0.25827953,0.31764944,0.030858048),(0.246812449,0.316003333,0.0366106886),(0.2184549,0.305568396,0.0755400306),(0.262275929,0.317030598,0.0324706064),(0.255792817,0.349480805,0.0465823635),(0.240159154,0.348568172,0.0664519652),(0.224344425,0.320323758,0.0763528482),(0.240169224,0.34717435,0.0698085694),(0.265199721,0.312864699,0.0410596215),(0.265418055,0.311917006,0.0443562989),(0.26129048,0.308155744,0.0544564936),(0.256743737,0.305555348,0.0625472758),(0.234434135,0.300827387,0.0792163957),(0.253512449,0.354959623,0.0438396917),(0.249323628,0.353532398,0.0429945625),(0.234310634,0.3363848,0.0540814122),(0.220838934,0.301548046,0.0763120799)),fixed=False, color=(0.75,0.65,0.65))
ldplt13=polyhedra_utils.polyhedra(matP,v=((0.0574645961,0.172726493,0.0671611913),(0.0464122945,0.165099045,0.0713237631),(0.0434476564,0.157901448,0.0760789028),(0.0414561239,0.164067278,0.0478827693),(0.072804328,0.142350171,0.0317533935),(0.0951504821,0.147789042,0.0382296085),(0.0842231943,0.142831466,0.0614289276),(0.0809340468,0.155784735,0.068892435),(0.0769811439,0.146200206,0.0656535174),(0.0709411025,0.148834867,0.0690432875),(0.0633329365,0.164931365,0.0343370142),(0.0672340505,0.178422363,0.0464103079),(0.0624502757,0.154961311,0.0298507522),(0.088604781,0.163134951,0.039345871),(0.062581808,0.139811536,0.0344520801),(0.0833753609,0.183393338,0.0581430944),(0.0916609346,0.156554103,0.0383732291),(0.087481942,0.179940447,0.0500532787)),fixed=False, color=(0.75,0.65,0.65))
ldplt14=polyhedra_utils.polyhedra(matP,v=((0.137126277,0.130191443,0.0780041999),(0.118460053,0.116506782,-0.00529935144),(0.109983529,0.105402621,0.0314843693),(0.102825224,0.0982376609,0.0176778245),(0.127831748,0.0846509587,0.0262975532),(0.114121196,0.133111933,0.0567199344),(0.0981515881,0.12751338,0.0365957518),(0.140370147,0.129607361,0.0272430126),(0.136670086,0.127673405,0.019710194),(0.0932267982,0.116734381,0.000279908956),(0.128543714,0.0839454539,0.0218201911),(0.105456474,0.11389696,-0.00127562561),(0.136503244,0.134817778,0.0739578036),(0.147591793,0.130157574,0.0472411862),(0.142701519,0.13113422,0.0747800576),(0.143476285,0.0932196341,0.0493843992),(0.145793991,0.0968693037,0.0545141114),(0.143448202,0.129699738,0.0328079612),(0.133434944,0.127347647,0.0149578266),(0.100721524,0.129280015,0.0374942704)),fixed=False, color=(0.75,0.65,0.65))
ldplt15=polyhedra_utils.polyhedra(matP,v=((0.103632481,0.401145448,0.0132739138),(0.102130572,0.395039424,0.0168004721),(0.102159955,0.400176862,0.00920747126),(0.0738353629,0.399031332,-0.00153158113),(0.0859211699,0.400459717,0.00189352932),(0.0814750245,0.383795284,0.00307432998),(0.0845130656,0.385808264,0.00278164584),(0.077517455,0.389624269,0.0138655024),(0.0786593649,0.399967363,0.000980320234),(0.0818988737,0.376080941,0.0165824208),(0.0714496021,0.374358172,0.014369618),(0.0625425751,0.376119247,0.00950190513),(0.0987360356,0.39990321,0.0148076774),(0.0675537951,0.394490387,-0.00098672711),(0.0576893306,0.383223237,0.00041517103),(0.0581060072,0.374770215,0.00638029538),(0.0746463769,0.374360587,0.0158721523),(0.0741783492,0.385403956,-0.000084683578)),fixed=False, color=(0.75,0.65,0.65))
ldplt16=polyhedra_utils.polyhedra(matP,v=((0.368841395,0.185337548,0.00739735188),(0.370965142,0.18188763,0.00733450836),(0.382304377,0.179662329,0.000869860577),(0.387379664,0.180977154,0.00237854679),(0.395130067,0.206785705,0.00290562866),(0.369873982,0.186138752,0.00938902072),(0.384738174,0.207110912,0.00260589343),(0.387469572,0.196310884,0.0112643884),(0.384758465,0.193425359,-0.000819404702),(0.392235349,0.191752168,0.000146778484),(0.387902588,0.192688275,0.0125674024),(0.386829289,0.191538877,0.0124871037),(0.394821197,0.188620423,0.0112232997),(0.377181808,0.198799368,0.0080499323),(0.396630783,0.191211401,0.0106831106),(0.377540288,0.180426406,0.00153849984),(0.380690892,0.205104977,0.00110447945),(0.39947236,0.206755698,0.00533054923),(0.398738338,0.202966484,0.00724760915),(0.39861165,0.199952439,0.00250029357),(0.398114478,0.193641409,0.00820019608)),fixed=False, color=(0.75,0.65,0.65))

O.bodies.append(ldplt1)
O.bodies.append(ldplt2)
O.bodies.append(ldplt3)
O.bodies.append(ldplt4)
O.bodies.append(ldplt5)
O.bodies.append(ldplt6)
O.bodies.append(ldplt7)
O.bodies.append(ldplt8)
O.bodies.append(ldplt9)
O.bodies.append(ldplt10)
O.bodies.append(ldplt11)
O.bodies.append(ldplt12)
O.bodies.append(ldplt13)
O.bodies.append(ldplt14)
O.bodies.append(ldplt15)
O.bodies.append(ldplt16)

O.timingEnabled=True

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Sphere_Aabb()],verletDist=.000005),
InteractionLoop(
[Ig2_Wall_Polyhedra_PolyhedraGeom(),Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom(),Ig2_Sphere_Sphere_ScGeom(),Ig2_Sphere_Polyhedra_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),Ip2_FrictMat_PolyhedraMat_FrictPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(),Law2_ScGeom_FrictPhys_CundallStrack()],   # contact law -- apply forces #Bo1_Cylinder_Aabb()
),
NewtonIntegrator(damping=0.3,gravity=(0,0,-39.8)),
PyRunner(command='cycle()',iterPeriod=100,label='step'),
PyRunner(command='Compact()',iterPeriod=1,label='compact'),
PyRunner(command='save()',iterPeriod=1000,label='compact')
]

def cycle():
for b in O.bodies:
if isinstance(b.shape,Polyhedra):
if b.state.pos[0] < 0 or b.state.pos[0] > 0.4:
O.bodies.erase(b.id)
elif b.state.pos[1] <0 or b.state.pos[1] > 0.4:
O.bodies.erase(b.id)
elif b.state.pos[2] < -0.15 or b.state.pos[2] > 0.85:
O.bodies.erase(b.id)
else:
pass

def Compact():
global a
if a == 1:
l=O.bodies[6].state.pos[2]#!
if l>0.8 : return#!
ldpltheight=max([b.state.pos[2] for b in O.bodies if isinstance(b.shape,Polyhedra)])
print ("pos=%.5f"%(ldpltheight))
height=ldpltheight+0.03
ldplt=polyhedra_utils.polyhedra(steel,v=((0.4,0.4,height),(0.4,0,height),(0,0,height),(0,0.4,height),(0.4,0.4,height+0.02),(0.4,0,height+0.02),(0,0,height+0.02),(0,0.4,height+0.02)),fixed=False, color=(0.75,0.65,0.65))
O.bodies.append(ldplt)
a=a+1
elif a == 2:
Lastnum=O.bodies[-1].id
O.bodies[Lastnum].state.blockedDOFs='xyXYZ'
O.forces.setPermF(Lastnum,(0,0,-5000))
a=a+1
elif a == 3:
Lastnum=O.bodies[-1].id
plateF=O.forces.f(Lastnum)[2]
h=O.bodies[-1].state.pos[2]
print ("h=%.8f"%(h))
print ("force")
if h<0:#!
O.pause()
else:
if utils.unbalancedForce()<0.25:
a=a+1
elif a == 4:
Lastnum=O.bodies[-1].id
O.bodies[Lastnum].state.vel=[0,0,-0.25]
plateF=O.forces.f(Lastnum)[2]
h=O.bodies[-1].state.pos[2]
print ("h=%.8f"%(h))
print ("vel")
if h>0:return#!
O.pause()

def save():
export.textPolyhedra('/home/hello/wenjian/packtestlayer4-4data.dat')
pass

O.dt=5E-7

--

```