Paraview crashes when pressing apply buttom of the micro-strain vtk file

Hello,

> I am not sure this question is suitable to be asked here

yes

Are you using the script from [1], or something else?
What version of OS, Yade and Paraview are you using?
The resulting "strain.vtk" is a plain text file. Please copy-paste it here.

Thanks
Jan

Hi Jan,

Thanks for your quick reply.

> Are you using the script from [1], or something else?

I am using the same script from [1].

> What version of OS, Yade and Paraview are you using?

My operating system is Ubuntu 16.04 LTS. My Yade version is 2019-11-18.git-b1a2767 and Paraview version is 5.0.1 (64-bit).

> The resulting "strain.vtk" is a plain text file. Please copy-paste it here.

Yes, it's a plain text file and the content is very long. I will email it to you privately.

Thank you for your time!

Cheers,
Chien-Cheng

Hello,

I am using Ubuntu 18.04, Paraview 5.4.1
I have also the crash with the file you have sent me.
BUT I have no problem viewing the results of the script using Yade 2018.02b and 2020-04-26.git-59c7045

Maybe you are just unlucky with Yade version..

cheers
Jan

Hi Jan,

Thanks for helping me test the file.

Would you suggest me to update the Yade version to 2020-04-26.git-59c7045 and rerun the simulation?

Cheers,
Chien-Cheng

you can try 2020-04-26.git-59c7045
or you can simply try yadedaily [1] (no need of long compilation) to test the newest Yade version
cheers
Jan

Hi Jan,

I've updated my Yade version to the latest one - 2020-06-23.git-f03a37c and reran the simulation to get a new microstrain file.
But the Paraview still crashes when I apply the file.
Is it possible that Paraview cannot read the output microstrain file (a plain text file)? Or perhaps you could share me a workable microstrain file for letting me test it?
Thank you.

Cheers,
Chien-Cheng

What version is the Paraview you use?

> Is it possible that Paraview cannot read the output microstrain file (a plain text file)?

it is possible (e.g. there can be some NaN with which the export code creator did not count)

> Or perhaps you could share me a workable microstrain file for letting me test it?

sent to personal email

cheers
Jan

Hi Jan,

Thanks for your email.

I've tested your microstrain vtk file and it works with my Paraview. So I think the crash is not related to the Paraview version (btw my Paraview version is 5.0.1 (64-bit)). Maybe it's related to my code. Could you help me have a look? I'd appreciate that.

### My simplified code (direct shear simulation of a granular layer)
from yade import pack,plot,export
import numpy as np
import math

sp1=pack.SpherePack()
sp2=pack.SpherePack()
sp3=pack.SpherePack()

O.periodic=True

# dimensions of sample
RADIUS=0.25
a=15
b=1
c=1
length=a*(2*RADIUS)
height=length/b
width=length/c
thickness=RADIUS

### Guassian distribution
psdSizes=[.456,.5,.544]
psdCumm=[0,0.5,1]

# friction angles
wallFRIC=0
boundaryFRIC=0.5
spFRIC=0.5

# boundary conditions
PI=1.e5
SN=5.e6 # normal stress
RATE_NS1=1 # velocity of top plate during compaction
RATE_NS2=1 # velocity of top plate during shear
RATE_shear=1 # shear velocity
roll_stiff=0
roll_fric=0

# simulation control
DAMPSHEAR=0.

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

O.materials.append(FrictMat(density=3000,young=1e8,poisson=0.5,frictionAngle=wallFRIC,label='boxMat'))
O.materials.append(FrictMat(density=3000,young=1e8,poisson=0.5,frictionAngle=boundaryFRIC,label='boundaryMat'))
O.materials.append(FrictMat(density=3000,young=1e8,poisson=0.5,frictionAngle=spFRIC,label='sphereMat'))

upBox = utils.box(center=(length/2,2*height+thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2,height-thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,width),fixed=1,wire=False,color=(1,0,0),material='boxMat')

O.bodies.append([upBox,lowBox])

sp1.makeCloud((0,height+3*RADIUS,width),(length,2*height-3*RADIUS,2*width), psdSizes =psdSizes, psdCumm =psdCumm, periodic=True, seed =1)
sp2.makeCloud((0,height+RADIUS,width),(length,height+RADIUS-1e-10,2*width), rMean=RADIUS, periodic=True, seed =1)
sp3.makeCloud((0,2*height-RADIUS,width),(length,2*height-RADIUS-1e-10,2*width), rMean=RADIUS, periodic=True, seed =1)

sphere_id = O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='sphereMat') for s in sp1])
bottomLayer_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,0,1),material='boundaryMat') for s in sp2])
topLayer_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,0,1),material='boundaryMat') for s in sp3])

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_FrictMat_FrictMat_MindlinPhys(krot=roll_stiff,eta=roll_fric)],
  [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True)]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,PyRunner(command='fixVelocity(RATE_shear)',iterPeriod=1,label='fixVel',dead=True)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

def triaxDone():
 global phase
 volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])
 h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
 vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
 contactStress=getStress(vol)
 vol_s=Rmean=Rmax=nbSph=0
 Rmin=1e6
 x_ref=O.bodies[0].state.pos[0]
 for o in O.bodies:
  if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
   if o.shape.radius>Rmax: Rmax=o.shape.radius
   if o.shape.radius<Rmin: Rmin=o.shape.radius
   vol_s += 4.*pi/3.*(o.shape.radius)**3
 Rmean=Rmean/nbSph
 n = 1-vol_s/vol
 print ('DONE! iter=',O.iter,'| sample generated: nb spheres',nbSph,', Rmean=',Rmean,', Rratio=',Rmax/Rmin,', porosity=',n,'wall_refPos=',x_ref)
 tt=TriaxialCompressionEngine()
 #tt.setContactProperties(shearFRIC)
 triax.dead=True
 O.pause()

#### If top particle layer contact with wall then give the particle velocity, and vice versa.
def identify_topLayer():
    top = []
    for i in O.interactions.withBody(0):
        if i.id1 != 0:
            top.append(i.id1)
        if i.id2 != 0:
            top.append(i.id2)
    return top

O.run(10000000,1)

stage=0
stiff=fnPlaten=currentSN=0.
def servo():
 global stage,stiff,fnPlaten,currentSN
 if stage==0:
  currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  unbF=unbalancedForce()
  print ('SN=',SN,'| current SN = ',currentSN,' | unbF=',unbF )
  boundaryVel=copysign(min(RATE_NS1,abs(0.5*(currentSN-SN))),currentSN-SN)
  O.bodies[0].state.vel[1]=boundaryVel
  if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.005 ):
   stage+=1
   fnPlaten=O.forces.f(0)[1]
   ## the following computes the stiffness of the plate (used for stress control of the top plate)
   for i in O.interactions.withBody(O.bodies[0].id):
    stiff+=i.phys.kn
   O.pause()
 if stage==1:
  fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  boundaryVel=copysign(min(RATE_NS2,abs(0.333*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired)
  O.bodies[0].state.vel[1]=boundaryVel

O.engines = O.engines[:5]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[5:]

O.run(10000000,1)

#### fixed bottom bottomLayer_id
for i in bottomLayer_id:
    O.bodies[i].state.blockedDOFs='xyzXYZ'
    O.bodies[i].state.dynamics = False
    O.bodies[i].state.vel = Vector3(0,0,0)

### Start shearing
newton.damping=DAMPSHEAR
fixVel.dead = False

refTime = O.time

def fixVelocity(RATE_shear):
    O.bodies[0].state.vel[0] = RATE_shear
    topLayer_id = identify_topLayer()
    for i in topLayer_id:
        O.bodies[i].state.vel[0] = RATE_shear
    slip = (O.time-refTime)*RATE_shear
    h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
    ss = slip/h
    if ss > 1:
        O.pause()

O.run(10000,1)

TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()
TW.volume(10)
TW.setState(0)
O.run(100,True)
TW.setState(1)
TW.defToVtk("strain.vtk")

###

Cheers,
Chien-Cheng

>> Are you using the script from [1], or something else?
> I am using the same script from [1].
> Maybe it's related to my code.

next time please be more consistent, I have tested everything on [1]

with your code, my paraview crashes while trying to load the oustout strain.vtk file, too..
I had a look, but did not find anything suspicious..
I have no idea what could be the reason..

cheers
Jan

Hi Jan,

Sorry for my inconsistency and thank you for your time.
I will keep looking for other possibilities.

Cheers,
Chien-Cheng

Hi Jan,

I think I find the reason why the microstrain vtk file crashed.
The vtk file can be applied in Paraview when the degree of freedom of particles is not blocked.
So I guess when the particles are fixed, the microstrain vtk file can not work functionally.

Cheers,
Chien-Cheng