New question #693297 on Yade: https://answers.launchpad.net/yade/+question/693297
Hi All, I am trying to obtain the volume of individual pores in my sphere packing. I found this function in Yade documentation: volume( [ (int)id=0 ] ) → float Returns the volume of Voronoi’s cell of a sphere I have few questions here: 1- Is this the correct approach to get the pore volume for each individual cell ID? 2- Does this function gives the volume of the whole tetrahedron element (which includes part of the spheres surrounding the pore as in Fig.3a in [1]) or just the pore (i.e. fluid phase)? 3- In my code below, this function is not working and I get this error: Segmentation fault (core dumped) I have tried using this function in the terminal and I found that it sometimes give a value for some cell IDs and most of the time it will crush and this error occurs. My code is copied below. The packing txt file link is here [2] I will appreciate your help in answering my questions. Thank you, Othman [1] https://onlinelibrary.wiley.com/doi/full/10.1002/nag.2198 [2] https://drive.google.com/file/d/1JqFe3W1cVFgKEa9lqU-MsCBxY7RDgY89/view?usp=sharing --------------------------------------------- from yade import plot,pack, export, ymport import numpy as np P=4347 #Pa visc=1e-3 #Pa.sec density=1000 #kg/m3 g=9.81 #m/s2 radiuscyl=.05 ########## create walls ########## allx,ally,allz,r=np.genfromtxt('S30.txt', unpack=True) #access the packed cyl file (txt) mnx=min(allx)*0.999 mny=min(ally)*0.999 mnz=min(allz)*0.999 mxx=max(allx)*1.001 mxy=max(ally)*1.001 mxz=max(allz)*1.001 mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz) walls=aabbWalls([mn,mx],thickness=0) wallIds=O.bodies.append(walls) ########## spheres ########## spheres=O.bodies.append(ymport.text('/home/auser/DEM_codes/Clogging/S30.txt')) #please change this when you download the packing file yade.qt.View() Height=max(utils.aabbDim()) dP=P/Height #Pa/m for i in spheres: body= O.bodies[i] body.state.blockedDOFs='xyzXYZ' print ('porosity = ', utils.porosity()) Height=max(utils.aabbDim()) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop" ), FlowEngine(dead=1,label="flow"), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), NewtonIntegrator(damping=0.2) ] #B. Activate flow engine and set boundary conditions in order to get permeability flow.dead=0 flow.useSolver=3 flow.permeabilityFactor=1 flow.viscosity=visc flow.bndCondIsPressure=[1,1,1,1,1,1] flow.bndCondValue=[0,0,0,0,0,P] flow.boundaryUseMaxMin=[1,1,1,1,1,1] O.dt=1e-6 O.dynDt=False flow.emulateAction() ######################################################################################################################### ################################################# Get pores volume##################################### ########################################################################################################################## cellsHit = [] numPoints = 20 xs = np.linspace(0.95*mnx,1.05*mxx,numPoints) ys = np.linspace(0.95*mny,1.05*mxy,numPoints) zs = np.linspace(0.95*mnz,1.05*mxz,numPoints) v = np.array([0,0,0]) for x,y,z in itertools.product(xs, ys, zs): cellId = flow.getCell(x,y,z) #getCell return cell id for a point position xyz if cellId in cellsHit: continue #this prevents counting a cell multiple times cellsHit.append(cellId) VV=flow.volume(cellId) -- You received this question notification because your team yade-users is an answer contact for Yade. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
