error using saveVtk() to get tesselation geometry

Asked by rhaven on 2018-06-05

Hi There,
I try to use saveVtk() to take a look at the tesselation geometry.

I have added the following when defining the engines :
FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),

and
try
 flow.saveVtk()

However I get the error: Triangulation does not exist. Sorry.

I have tried adding the tesselationwrapper as follows
 TW=TesselationWrapper()
 TW.triangulate()
 TW.computeVolumes()

 flow.saveVtk()

However I still get the same error. Is there another step missing which computes the triangulation and makes it available for saveVtk?

many thanks
Jesse

Question information

Language:
English Edit question
Status:
Invalid
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
2018-06-29
Last query:
2018-06-29
Last reply:
2018-06-29

Hi,

FlowEngine contains a triangulation/tesselation data structure which shares a lot with TesselationWrapper data structure. However if you instantiate a FlowEngine and a TW in the same script it results in two independent data structure in memory and what you do with TW is completely unrelated to the output of saveVtk().

Yes, definitely, saveVtk will only output the result of some fluid solving. If the fluid is not solved (flow.dead=True or no timestepping) the output is empty.
You may try flow.emulateAction() before saveVtk(), which triggers what normally happens in a timestep (including building the triangulation).
I actually doubt outputing a vtk is what you need, but it might be another question.
HTH
Bruno

rhaven (rhavenj) said : #2

Hi Bruno,
Thank you for the quick reply. Ive added flow.emulateAction() as well as the following boundary conditions (otherwise i would get a bus error)
 flow.dead=0
 flow.defTolerance=0.3
 flow.meshUpdateInterval=200
 flow.useSolver=3
 flow.permeabilityFactor=1
 flow.viscosity=10
 flow.bndCondIsPressure=[0,0,1,1,0,0]
 flow.bndCondValue=[0,0,1,0,0,0]
 flow.boundaryUseMaxMin=[0,0,0,0,0,0]
 O.dt=0.1e-3
 O.dynDt=False
 flow.emulateAction()
 flow.saveVtk()

The console output however from saveVtk() is:
0 : Vh==NULL!! id=0 Point=1.92862 1.65482 4.78545 rad=0.738095
1 : Vh==NULL!! id=1 Point=3.64074 2.16205 4.14654 rad=0.714286
2 : Vh==NULL!! id=2 Point=2.03866 3.15345 0.0887077 rad=0.690476
3 : Vh==NULL!! id=3 Point=3.56674 0.0664596 3.07013 rad=0.666667
4 : Vh==NULL!! id=4 Point=0.54744 0.471346 2.35463 rad=0.642857
5 : Vh==NULL!! id=5 Point=0.469151 1.40331 4.06784 rad=0.619048
6 : Vh==NULL!! id=6 Point=3.80346 1.90056 1.84533 rad=0.595238
7 : Vh==NULL!! id=7 Point=0.882771 2.24981 0.484466 rad=0.571429
8 : Vh==NULL!! id=8 Point=4.90384 3.31689 4.64228 rad=0.547619
9 : Vh==NULL!! id=9 Point=2.64444 4.79931 2.32225 rad=0.52381
10 : Vh==NULL!! id=10 Point=2.22052 1.92293 1.89341 rad=0.5
11 : Vh==NULL!! id=11 Point=1.71494 1.92762 3.09753 rad=0.47619
12 : Vh==NULL!! id=12 Point=0.988165 4.71568 0.649083 rad=0.452381
13 : Vh==NULL!! id=13 Point=3.12047 3.11159 4.68382 rad=0.428571
14 : Vh==NULL!! id=14 Point=2.58058 2.02467 2.95446 rad=0.404762
15 : Vh==NULL!! id=15 Point=4.8716 2.04733 1.46188 rad=0.380952
16 : Vh==NULL!! id=16 Point=1.43367 4.63815 3.07543 rad=0.357143
17 : Vh==NULL!! id=17 Point=4.30153 4.29709 4.99369 rad=0.333333
18 : Vh==NULL!! id=18 Point=3.33361 1.48195 0.298009 rad=0.309524
19 : Vh==NULL!! id=19 Point=1.96329 0.208186 4.32673 rad=0.285714
CHOLMOD error: invalid xtype
CHOLMOD error: argument missing
Segmentation fault (core dumped)

What does Vh=NULL mean and is this the cause of the CHOLMOD and segmentation fault?

thank you
Jesse

Vh=NULL is when a particle is refused by thriangulation algorithm. I don't think it is the result of saveVtk since it is not triangulating, it must happen before.
Main reason for this problem is excessive overlap of some particles (center of one sphere inside another sphere), are you sure you have a consistent input?
B

rhaven (rhavenj) said : #4

Hi Bruno,
I am trying with several different packings, indeed one of the testing packings has overlapping spheres. However I have also just used makeClound to generate a packing of non overlapping spheres and I have the same issue.
Ive put a MWE at https://pastebin.com/bxCKYu99

To test with the other packing would be to uncomment
 sp = ymport.textExt('/tmp/divided.txt',format='x_y_z_r_attrs',attrs=attrs)
 n = max(int(a[0]) for a in attrs)+1
 colors = [randomColor() for _ in xrange(n)]
 agglomerates = [[] for _ in xrange(n)]
 for s,a in zip(sp,attrs):
    aa = int(a[0])
    s.agglomerate = aa
    s.shape.color = colors[aa]
    agglomerates[aa].append(s)
 for a in agglomerates:
    O.bodies.appendClumped(a)

and the packing https://pastebin.com/H5Bcfn26

hope this helps to clarify
many thanks
Jesse

Hi,
In your MWE the error does not come from saveVtk() (comment line 75-> still the same error).
And in this script also there are large overlaps between the spheres and the walls (and the walls are just very large spheres, i.e. same issue), that's because you generated a periodic cloud in a finite box.
Bruno

rhaven (rhavenj) said : #6

Hi Bruno,
thank you for the reply. I try with other packings but cant get around the Vh=NULL error
for example the code below with the following input file https://pastebin.com/8a2jsCqV

I understand that the error isnt necessarily coming from saveVTK, should I open a new question instead? Is there another way to access the triangulation data? My goal is actually to be able to read the edges and vertices of the regular triangulation.

code:
        O.reset()
 sphs = ymport.text('clump' + str(clumpNum) + '.txt')
 O.bodies.append(sphs)
 print('Number of particles: ', len(sphs))
 ##########

 raw_input('Before FlowEngine saveVtk')

 O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    #InteractionLoop(
    #[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], ##triax only
  # [Ig2_Sphere_Sphere_ScGeom()],
  # [Ip2_FrictMat_FrictMat_FrictPhys()],
  # [Law2_ScGeom_FrictPhys_CundallStrack()]
    #),
    FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    #PeriTriaxController(dynCell=True,mass=0.2,maxUnbalanced=0.01,relStressTol=0.02,goal=(-1e6,-1e6,-1e6),stressMask=3,globUpdate=5,maxStrainRate=(1e10,1e10,1e10),doneHook='triaxDone()',label='compressor'),
    #NewtonIntegrator(damping=.6), ##Triaxonly
    NewtonIntegrator(damping=.2),
 ]
 factor = .5
 O.dt = factor * PWaveTimeStep()

 #B. Activate flow engine and set boundary conditions in order to get permeability
 #TW=TesselationWrapper()
 #TW.triangulate()
 #TW.computeVolumes()
 flow.dead=0
 flow.defTolerance=0.3
 flow.meshUpdateInterval=200
 flow.useSolver=3
 flow.permeabilityFactor=1
 flow.viscosity=10
 flow.bndCondIsPressure=[0,0,1,1,0,0]
 flow.bndCondValue=[0,0,1,0,0,0]
 flow.boundaryUseMaxMin=[0,0,0,0,0,0]
 O.dt=0.1e-3
 O.dynDt=False
 flow.emulateAction()
 flow.saveVtk()
 raw_input('FlowEngine saveVtk END')

Hi, better reproduce the error with makeCloud to make a true MWE please.
Bruno

p.s. Why are you using FlowEngine instead of TesselationWrapper to get list of edges? Just curious.

rhaven (rhavenj) said : #8

Hi Bruno,
I looked at the tesselationwrapper, but didnt see a way to access the edges as a list pertaining to each voronoi cell. Is this possible?
I saw saveVtk can save out the voronoi diagram and seemed to be the easier way to access the geometry...

I should really say (taking yet another step back and looking at the bigger picture), what I really want to do is obtain a pore size distribution for my packing. (The pore size I am defining to be the radius of an inscribed sphere that fits in the space between particles.)
So to calculate something that looks like this https://www.researchgate.net/profile/Frank_Werner/publication/273284325/figure/fig7/AS:294826168864779@1447303444628/An-example-of-a-Voronoi-diagram-in-Laguerre-geometry-for-a-set-of-multi-sized-circles.png
many thanks
Jesse

TesselationWrapper also has some vtk-export functions, e.g. defToVtk(), but it doesn't seem to be the question in your case.
What you want is access to internals of the triangulation and going through vtk export just for this is a bit awkward. It would be probably more convenient to manipulate the data directly in python.
FlowEngine received more attention recently and it's python interface gives more access to internals. It is thus interesting to hear your feedback on what would me missing in TesselationWrapper.

But in the end, if what you want is the set of inscribed spheres my impression is that you can get it directly from TwoPhaseFlowEngine (I don't have the precise procedure in mind right now unfortunately).

Bruno

rhaven (rhavenj) said : #10

I've tried using the TwoPhaseFlowEngine to obtain the inscribed spheres. However I had issues where the pore sizes calculated were all very very small. I guess this lead me to try to vizualise the pore structure by exporting it to vtk. Perhaps these are all the same issue. For example I have added :

TwoPhase=TwoPhaseFlowEngine()
TwoPhase.initialization()

to the script I posted on 2018-06-19 and I get the following error:
negative volume for an ordinary pore (temp warning, should still be safe)
*** stack smashing detected ***: /usr/bin/python terminated
Aborted (core dumped)

how do you suggest getting to the bottom of the triangulation/tesselation issue?
many thanks
Jesse

I can't understand this triangulation failure unless you prove it with a MWE (using makeCloud). It seems you are using some weird input, that's my guess at the moment. It is completely unrelated to the other questions like getting inscribed radius.
Bruno

rhaven (rhavenj) said : #12

Hi Bruno,
thank you again for the reply. I started making a MWE using makeCloud and it ends up that both saveVtk and the TwoPhaseFlowEngine work fine (although when I load the vtk in paraview I dont know how to see the triangulation, is there documentation for this?).

 However, when I try loading my geometry the tesselation doesnt seem work any more. Does this mean if I would like to twophaseflowEngine or the tesselation in general the geometry needs to be regenerated with makeCloud? Could O.bodies.append(sphere(c,r)) work as well?
many thanks
Jesse

Hi, loading the vtk file in paraview should already display a mesh in principle. Not the case?

Both methods work: loading a geometry or makeCloud, they both go through O.bodies.append at some point.
The problem is not the method, it is the data. There must be some inconsistency in your specific way or reloading something resulting in unexpected geometry (did you inspect visually at least? for answer #5 it was enough to open the 3D view to understand the problem).

rhaven (rhavenj) said : #14

Hi,
loading the vtk file in paraview I see a convex mesh. I dont see how to see the triangulated mesh. http://pasteall.org/pic/5aaceafb92ef561b8d5e30c463caac59, the vtk file ive pasted here http://pasteall.org/1008402/text

I've made another MWE where I load clump0.txt reset the scene and use O.bodies.append(sphere(pos,r)) to add the spheres individually.
(I also dont see anything wrong when I look at the spheres in the 3d viewport. )

When I run the twophaseflowengine I print the cell id and pore size when TwoPhaseFlowEngine::computePoreBodyRadius() runs (I modified the source code of TwoPhaseFlowEngine.cpp):
cout << endl << "Test3 Cell: " << cell->info().id << ": Pore: " << Rin;

The output is a long list of pores all the same size ending with *** stack smashing detected ***: /usr/bin/python terminated

.........
Test3 Cell: 914: Pore: 4.92e-10
Test3 Cell: 915: Pore: 4.92e-10
Test3 Cell: 916: Pore: 4.92e-10
Test3 Cell: 917: Pore: 4.92e-10
Test3 Cell: 918: Pore: 4.92e-10
Test3 Cell: 919: Pore: 4.92e-10
Test3 Cell: 920: Pore: 4.92e-10negative volume for an ordinary pore (temp warning, should still be safe)
negative volume for an ordinary pore (temp warning, should still be safe)
negative volume for an ordinary pore (temp warning, should still be safe)
*** stack smashing detected ***: /usr/bin/python terminated
Aborted (core dumped)

When I try to run saveVtk the error is as follows.

......
148 : Vh==NULL!! id=148 Point=7.42392e-07 6.92963e-07 8.02883e-07 rad=1.23e-08
149 : Vh==NULL!! id=149 Point=8.03578e-07 7.18035e-07 8.15755e-07 rad=1.23e-08
150 : Vh==NULL!! id=150 Point=7.43404e-07 7.82494e-07 7.64931e-07 rad=1.23e-08
151 : Vh==NULL!! id=151 Point=6.77428e-07 7.66925e-07 7.91695e-07 rad=1.23e-08
CHOLMOD error: invalid xtype
CHOLMOD error: argument missing
Segmentation fault (core dumped)

clump0.txt I have pasted here: https://pastebin.com/WY703ZLs

def tryvtk():
 O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    #InteractionLoop(
    #[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], ##triax only
  # [Ig2_Sphere_Sphere_ScGeom()],
  # [Ip2_FrictMat_FrictMat_FrictPhys()],
  # [Law2_ScGeom_FrictPhys_CundallStrack()]
    #),
    FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    #PeriTriaxController(dynCell=True,mass=0.2,maxUnbalanced=0.01,relStressTol=0.02,goal=(-1e6,-1e6,-1e6),stressMask=3,globUpdate=5,maxStrainRate=(1e10,1e10,1e10),doneHook='triaxDone()',label='compressor'),
    #NewtonIntegrator(damping=.6), ##Triaxonly
    NewtonIntegrator(damping=.2),
 ]
 factor = .5
 O.dt = factor * PWaveTimeStep()
 flow.dead=0
 flow.defTolerance=0.3
 flow.meshUpdateInterval=200
 flow.useSolver=3
 flow.permeabilityFactor=1
 flow.viscosity=10
 flow.bndCondIsPressure=[0,0,1,1,0,0]
 flow.bndCondValue=[0,0,1,0,0,0]
 flow.boundaryUseMaxMin=[0,0,0,0,0,0]
 O.dt=0.1e-3
 O.dynDt=False
 flow.emulateAction()
 flow.saveVtk()
 raw_input('FlowEngine saveVtk END')

def main():

 from yade import export,ymport
 import random
 random.seed(1)

 dim = Vector3(2e-6,2e-6,2e-6)

 O.materials.append(FrictMat(density=1e6,frictionAngle=0))
 O.periodic=True
 O.cell.setBox(1e-6,1e-6,1e-6)

 attrs = []

 #sp = ymport.textExt('/tmp/clump0.txt',format='x_y_z_r_attrs',attrs=attrs)
 sp = ymport.textClumps('/tmp/clump0.txt')

 spheres = []
 for b in O.bodies:
  if b.isClump == False:
   spheres.append([b.state.pos, b.shape.radius])

 raw_input('Pause1')

 O.resetThisScene() #####!!!!!!!

 raw_input('Pause1.1')
 for pos,r in spheres:
  O.bodies.append(sphere(pos,r))
  #print(pos,r)
 raw_input('Pause1.2')
 #quit()

 TwoPhase=TwoPhaseFlowEngine()
 TwoPhase.initialization()

 raw_input('Pause2')

 tryvtk()

 raw_input('Pause3')
 #quit()

 print('Number of particles = ' , len(O.bodies))
 print('Total particles inc. periodic', len(O.bodies))

 try:
  from yade import qt
  qt.View()
 except:
  pass

if __name__=="__main__":
 main()

You Vtk data file is ok, if you want to see internal mesh you need to "clip" or "slice" (and display "surface with edges" instead of wireframe).
For the second part: I guess you have overlapping spheres in the clump, and triangulation is possible only for limited overlaps.
Bruno

rhaven (rhavenj) said : #16

Thanks for the quick reply. With surface with edges and a clip I can see the internal mesh.

I made sure this time to avoid overlapping spheres. The clump is generated by using makeCloud, compressing the volume with triaxcontroller and saving only spheres which are inside a predicate.

Not sure what happens with clumps. If you can summarize the problem and reproduce in a simple way with a self-contained script (sorry but I'm tired of downloading extra files on external servers today...) I may have a look.
Bruno

rhaven (rhavenj) said : #18

Thank you for the reply.
Ive made a new MWE which also generates clump0.txt.
https://pastebin.com/1BdUpm2E

WARNING! rmin>rmax. rmin=1e-11 ,rmax=1.66533e-12
WARNING! rmin>rmax. rmin=1e-11 ,rmax=1.66533e-12
WARNING! rmin>rmax. rmin=1e-11 ,rmax=3.93018e-12
WARNING! rmin>rmax. rmin=1e-11 ,rmax=3.93018e-12
Segmentation fault

In addition, as I was seeing random lines from your last script I found:

def createAggregate(radAgg,clumpNum):
 O.reset()
        .....

then:
 for idx, radAgg in enumerate(aggRadii):
  createAggregate(radAgg, idx)

Which basically means that you keep creating new clumps, adding them to O.bodies, and finally erasing everything (including O.engines) before creating the next clump...

Another exotic thing is at the end:
if __name__=="__main__":
 main()

No clue what it's doing. Maybe it is not wrong, but for sure it is not what we expect in a MWE.

This was supposed to be about saveVtk() but it is much too confused.
I think I'll give up on this question, please really read [1] again if it can help.

B

closing

rhaven (rhavenj) said : #21

My apologies that the script is convoluted. Ive had to merge many scripts into one to make a single MWE that doesnt rely on any external files. However your description is correct. I create a clump, add it to the simulation save it. Then delete everything again to before moving on to create a new clump.

Once I have created a clump I want to have information from the tesselation. When I just create a particle cloud with makeCloud I dont have an issue running the tesselationwrapper or twophaseflowengine. However when I try to run it on the clump I get those segmentations faults.

The if __name__=="__main__":
doesnt have an effect on the script it just lets me the python file and use the functions in other scripts. Sorry for including it

I really would appreciate your continued help on this question.

> delete everything again to before moving on to create a new clump

If it crashes on clump n°1 then there is no need for the MWE to erase and restart.
If it crashes on clump n°X then send a MWE using only clump n°X.

> when I try to run it on the clump I get those segmentations faults.

May I point out that the warnings+segfault which popped up in #19 had never been mentioned before and I have no reason to believe it is related to saveVtk() or anything else discussed previously?

> The if __name__=="__main__": doesnt have an effect on the script

Why is it in the script then?

> I really would appreciate your continued help on this question.

No sorry. Please keep this question closed. When you have things cleared out maybe try another question which will follow strictly every bit of [1] (especially point 3).
Good luck
Bruno

[1] https://yade-dem.org/wiki/Howtoask