Yade

flow.nCells() != "cell number" recorded by flow.saveVtk()

Asked by Lingran Zhang on 2019-11-22

Hello,

One question for FlowEngine and DFNFlowEngine: why the flow.nCells() != "cell number" recorded by flow.saveVtk()?

I would like to record the Ids of initial fractured cells and then impose different pressures within them, but I find that the ids (and vertices) are not consistent given by these two ways.

Below is an example.
Thank you in advance for your help.
Best regards,
Lingran

from yade import pack, ymport

### material
def sphereMat(): return JCFpmMat(type=1,density=3000,young=9e9,poisson=0.2,frictionAngle=radians(18))
def wallMat(): return JCFpmMat(type=0,density=3000,young=9e9,poisson=0.2,frictionAngle=radians(0))

### bodies
mn,mx=Vector3(0,0,0),Vector3(1,1,1)
O.bodies.append(aabbWalls([mn,mx],oversizeFactor=1.5,thickness=0.1,material=wallMat))
sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((0.,0.,0.),(1,1.,1.)),spheresInCell=100,radius=1/10.,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation(material=sphereMat)

dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

### engines
flow=FlowEngine(
        isActivated=1,
        ### choose solver to use (0: Gauss Seidel, 1: Taucs, 2: Pardiso, 3: CHOLMOD)
        useSolver=3, # 3 should be used by default but does not work everytime...
         boundaryUseMaxMin = [0,0,0,0,0,0],
        bndCondIsPressure = [1,1,0,0,0,0],
        bndCondValue = [1,0,0,0,0,0],
        permeabilityFactor=1,
        viscosity=0.001,
        #debug=1
)

intR=1.1
O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')])
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
 )
 ,flow
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8)
 ,NewtonIntegrator()
]

### simulation starts here

## to block the particles (better for permeability test)
for o in O.bodies:
 o.state.blockedDOFs+='xyzXYZ'

O.run(1,True)
## getBoundaryFlux get the total discharge [m3/s]
Qin = flow.getBoundaryFlux(0)
Qout = flow.getBoundaryFlux(1)

permeability = abs(Qout)*flow.viscosity*X/(Y*Z) # !!! if Pout=1, Pin=0
permeability2 = flow.averageVelocity()*flow.viscosity*X # !!! if Pout=1, Pin=0
conductivity = permeability*1000*9.82/flow.viscosity # K=rho*g*k/nu
print "Qin=",Qin," Qout=",Qout," ARE THEY EQUAL? IF NOT-> NO FLOW!"
print "Permeability [m2]=",permeability," || Hydraulic conductivity [m/s]=",conductivity
flow.saveVtk() # to see the result in Paraview

for i in range(flow.nCells()):
   print i, flow.getVertices(i)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Lingran Zhang
Solved:
2019-11-24
Last query:
2019-11-24
Last reply:
2019-11-23
Robert Caulk (rcaulk) said : #1

Hello Lingran,

I am not sure I understand the question:

flow.nCells() returns the total number of finite cells in the triangulation [1].

flow.saveVTK() prints the individual cell IDs for each cell (the ones you can index with getVertices(i)).

flow.getVertices(i) gives the vertex IDs for the particles comprising cell i.

Are you trying to say the maximum ID printed by saveVTK is not equal to the single value returned by nCells()?
What makes you say the vertices are inconsistent? What are the vertices inconsistent with? The particle body IDs?

Cheers,

Robert

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.nCells

Lingran Zhang (lingran) said : #2

Hello Robert,

Thank you.

In my simulation, I get flow.nCells()=383; flow.getVertices(0)= [5, 68, 69, 71].
It means that the particle ids comprising the the cell id=0 are 5, 68, 69, 71, right?

While when I open the VTK file printed by flow.saveVTK(), I read the following lines:
CELLS 242 1210
4 38 41 37 39
4 31 28 36 11
…
From my understanding, the VTK file tells me that I have 242 cells, each cell include 4 vertices, the particle ids comprising the cell id=0 are 38, 41, 37, 39. Am I wrong?

So my questions are: how many finite cells are produced by the triangulation, 383 or 242? What are the particle ids comprising each cell, shall I find the answer from flow.getVertices(i) or flow.saveVTK()?

Cheers,
Lingran

Robert Caulk (rcaulk) said : #3

Thanks for the clarification.

> how many finite cells are produced by the triangulation, 383 or 242?

242 withoutBoundaries (default). Try using flow.saveVTK(withBoundaries=True), the numbers should match.

>the particle ids comprising the cell id=0 are 38, 41, 37, 39. Am I wrong?

Yes, you are wrong :-) Those are not the particle IDs. Those are redirected IDs for paraview since it needs zero based vertex IDs [1]

[1]https://gitlab.com/yade-dev/trunk/blob/master/lib/triangulation/FlowBoundingSphere.ipp#L1690

Lingran Zhang (lingran) said : #4

Hi Robert,

Great, you have answered exactly my questions.
Now I should update my source code to include ‘withBoundaries’ in flow.saveVTK().
Thank you very much.

Cheers,
Lingran