Pore size

Hi,

I am working on a similar problem as explained in the question [1]. I am trying implement the solution suggested by Robert in #7 by using print all the vertices [2].
I have compacted spherical particles into a cylinder shape, and want estimates/measure pore size in the compact. I was running the suggestion by print all vertices..
When I am using this I am getting following error message:

"ArgumentError: Python argument types in
    FlowEngineT.printVertices()
did not match C++ signature:
    printVertices(yade::TemplateFlowEngine_FlowEngineT<yade::FlowCellInfo_FlowEngineT, yade::FlowVertexInfo_FlowEngineT, yade::CGT::_Tesselation<yade::CGT::TriangulationTypes<yade::FlowVertexInfo_FlowEngineT, yade::FlowCellInfo_FlowEngineT> >, yade::CGT::FlowBoundingSphereLinSolv<yade::CGT::_Tesselation<yade::CGT::TriangulationTypes<yade::FlowVertexInfo_FlowEngineT, yade::FlowCellInfo_FlowEngineT> >, yade::CGT::FlowBoundingSphere<yade::CGT::_Tesselation<yade::CGT::TriangulationTypes<yade::FlowVertexInfo_FlowEngineT, yade::FlowCellInfo_FlowEngineT> > > > > {lvalue})
"

My question how should implement this code to extract the vertices and type of the void/pore space? What should I use as an input?

Best regards,
Mithu

[1] https://answers.launchpad.net/yade/+question/700963
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.printVertices

MWE:

# -*- coding: utf-8 -*-
"""
Created on Thu Jul 15 16:29:55 2021

@author: njb18198
"""

#!/usr/bin/env python
#encoding: ascii

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test
# The reference paper [Haustein2017]
from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export
from scipy.interpolate import interp1d
from csv import writer
import os
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import csv
from matplotlib.pyplot import figure
from pylab import *
from scipy.optimize import curve_fit
readParamsFromTable(comp_press=0.54e8,h_tab=2.05,m_tab=0.1974, r_tab=5.02,tab_porosity=20,tab_height=1,save=0)
from yade.params.table import *

o = Omega()
save=save
# Physical parameters
fr = 0.41
rho = 1561
Diameter = 7.9e-5
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 10000
kp = 140000
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

if save==0:
    g=-9.81
elif save==1:
    g=0

o.dt = 1.0e-8
particleMass = (4.0/3.0)*math.pi*r1*r1*r1*rho
Tab_rad=0.001
r_tab=r_tab*1e-3 #real size
h_tab=h_tab*1e-3
m_tab=m_tab*1e-3
v_tab=math.pi*(r_tab**2)*h_tab
v_1mm=math.pi*(Tab_rad**2)*(0.4e-3)
#m_tot=4e-06
m_tot=(v_1mm/v_tab)*m_tab
no_p=m_tot/particleMass
#no_p=8000
Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

Cyl_height=0.005
cross_area=math.pi*(Tab_rad**2)

Comp_press_up= comp_press
Comp_force_up=Comp_press_up*cross_area

Comp_press_lp= comp_press
Comp_force_lp=Comp_press_lp*cross_area

compression_data_save=[]
sc_por_15=2
#sc_por_2=2
#sc_por_1=1
data_to_save=[comp_press/1e6,round(no_p),rho]
compression_data_save.append(data_to_save)

# Add material
mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression and walls
sp=pack.SpherePack()
sp.makeCloud((-7.5*Diameter,-7.5*Diameter,-30*Diameter),(7.5*Diameter,7.5*Diameter,5.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=round(no_p))
sp.toSimulation(material=mat1)
walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))

vtkRecorder = VTKRecorder(fileName='vtkRecorder_'+str(tab_porosity),recorders=['all'])
tab_porosity=tab_porosity
tab_height=tab_height

##Single particle swelling model
def model(r,t,Q_max,rho_t,rho_w,r_0,Diff):
    Q=((rho_w*(r**3))/(rho_t*(r_0**3)))-(rho_w/rho_t)+1;
    drdt =((Diff*rho_t)/(r*rho_w))*((Q_max-Q)/Q);
    return drdt
P=[1.45,rho,1000,396.39e-12]

# Add engines
o.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
                         Bo1_Wall_Aabb(),
                         Bo1_Facet_Aabb()
                         ]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
     Ig2_Facet_Sphere_ScGeom(),
     Ig2_Wall_Sphere_ScGeom()],
    [Ip2_LudingMat_LudingMat_LudingPhys()],
    [Law2_ScGeom_LudingPhys_Basic()]
  ),
  NewtonIntegrator(damping=0.1, gravity=[0, 0, g]),
  PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
  #DeformControl(label="DefControl")
]

def checkForce():
    if O.iter < 3000000:
        return
    if unbalancedForce() > 1:
        return
    global upper_punch
    upper_punch=O.bodies.append(geom.facetCylinder((0,0,((-Cyl_height/2)+0.0001)+utils.aabbDim()[2]),Tab_rad-.00001,0,segmentsNumber=50,wallMask=1))
    for i in upper_punch:
        body= O.bodies[i]
        body.state.vel = (0,0,-0.02)
    global lower_punch
    lower_punch= O.bodies.append(geom.facetCylinder((0,0,(-Cyl_height/2)-0.0001),Tab_rad-.00001,0,segmentsNumber=50,wallMask=1))
    for n in lower_punch:
        body= O.bodies[n]
        body.state.vel = (0,0,0.02)
    O.engines = O.engines + [PyRunner(command='storeData()', iterPeriod=25000)]+ [PyRunner(command='saveData()', iterPeriod=100000)]
    fCheck.command = 'unloadPlate()'

def unloadPlate():
    force_up=0
    for i in upper_punch:
        body= O.bodies[i]
        force_up=force_up+abs(O.forces.f(body.id)[2])
    force_lp=0
    for n in lower_punch:
        body = O.bodies[n]
        force_lp = force_lp + abs(O.forces.f(body.id)[2])
    if ((force_up > Comp_force_up) and (force_lp > Comp_force_lp)):
        for i in upper_punch:
            body= O.bodies[i]
            body.state.vel = (0,0,0.04)
        for n in lower_punch:
            body= O.bodies[n]
            body.state.vel = (0,0,-0.04)
        fCheck.command = 'stopUnloading()'

def stopUnloading():
    #force_up=0
    #for i in upper_punch:
        #body= O.bodies[i]
        #force_up=force_up+(O.forces.f(body.id)[2])
    force_lp=0
    for n in lower_punch:
        body = O.bodies[n]
        force_lp = force_lp + abs(O.forces.f(body.id)[2])
    if force_lp==0:
        for i in lower_punch:
            body= O.bodies[i]
            body.state.vel = (0,0,0)
    #if ((force_up==0) and (force_lp==0)):
    for i in upper_punch:
        body=O.bodies[i]
        pos_up=body.state.pos
    for i in lower_punch:
        body=O.bodies[i]
        pos_lp=body.state.pos
    if pos_up[2]> pos_lp[2]+utils.aabbDim()[2]+0.0002:
        for j in upper_punch: O.bodies.erase(j)
        fCheck.command = 'Savecheck()'

def Savecheck():
    if save==1:
        utils.saveVars('lower_punch',lower_punch=lower_punch)
        save_filename='PH101_'+str(tab_porosity)+'_'+str(tab_height)+'mm.xml'
        o.save(save_filename)
        o.pause()
    if save==0:
        o.engines = o.engines+[PyRunner(command='ParticleSwelling()', iterPeriod=100000)]
        fCheck.dead = True # (!!!)
        storeData.dead=True
        saveData.dead=True

if save==1:
    O.run()
    # when running with yade-batch, the script must not finish until the simulation is done fully
    # this command will wait for that (has no influence in the non-batch mode)
    waitIfBatch()
    g=-9.81

if save==0:
    read_filename='PH101_'+str(tab_porosity)+'_'+str(tab_height)+'mm.xml'
    o.load(read_filename)
    utils.loadVars('lower_punch')
    from yade.params.lower_punch import *
    for b in O.bodies:
        if isinstance(b.shape, Sphere):
            #b.state.blockedDOFs = 'xyXY'

            r=b.shape.radius
            oldm=b.state.mass
            oldI=b.state.inertia

            m=oldm*3./4./r
            b.state.mass=m

            b.state.inertia[0] = 15./16./r*oldI[0] #inertia with respect to x and y axes are not used and the computation here is wrong
            b.state.inertia[1] = 15./16./r*oldI[1] #inertia with respect to x and y axes are not used and the computation here is wrong
            b.state.inertia[2] = 15./16./r*oldI[2] #only inertia with respect to z axis is usefull
    o.dt=1e-6))
    for j in lower_punch: O.bodies.erase(j)
    for b in O.bodies:
        if isinstance(b.shape, Sphere):
            continue
        else:
            O.bodies.erase(b.id)
    g=0
    walls=O.bodies.append(yade.geom.facetCylinder((0,0,-0.00095),radius=Tab_rad+0.0008,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))
    FlowEngine(dead=1,label="flow")
    FlowEngine.printVertices()