Output specific bodies using TW.setState

Asked by Chien-Cheng Hung on 2021-03-19

Hi all,

Followed by my previous question [1], I found that the output micro-strain VTK plain text file from the TW.defToVtk [2] cannot be opened in the Paraview if the DOFs (degree of freedoms) of particles are blocked. For example, I wrote something like
###
for b in bottomBoundary:
    b.state.blockedDOFs='xyzXYZ'
###
to make sure part of the particles is fixed in all directions. When these particles are fixed, the output microstrain vtk text file cannot be opened in the Paraview. But when the particles are not fixed (without blockedDOFs), the vtk file works in the Paraview.

So I was wondering if I could only store the position of certain particles that are not fixed by adjusting the TW.setState [3] instead of using TW.setState(0) or TWsetState(1) to store the position of the entire body.
Thanks!

[1] https://answers.launchpad.net/yade/+question/691407
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.TesselationWrapper.defToVtk
[3] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.TesselationWrapper.setState

Cheers,
Chien-Cheng

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Karol Brzezinski
Solved:
2021-03-22
Last query:
2021-03-22
Last reply:
2021-03-21

Hi, i have assumptions on why it would be the case but nothing sure.
Could you please provide a minimal working example of the problem?
B

Chien-Cheng Hung (chiencheng) said : #2

Hi Bruno,

Thanks for your quick reply.

I've minimized my script as much as possible. I also tested it and the running time should be within 30 seconds.

Below is the MWE.

###
from yade import pack,plot,export
import numpy as np
import math

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

O.periodic=True

RADIUS1=0.25
RADIUS2=0.125
length=3*(2*RADIUS1)
height=3*(2*RADIUS1)
width=3*(2*RADIUS1)
PI=1.e5
SN=5.e6
RATE_shear=1
wallFRIC=0
boundaryFRIC=0.5 # during compaction (controls porosity)
spFRIC=0.5
DAMPSHEAR=0.

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

O.materials.append(FrictMat(density=2500,young=5.5e10,poisson=0.25,frictionAngle=wallFRIC,label='boxMat'))
O.materials.append(FrictMat(density=2500,young=5.5e10,poisson=0.25,frictionAngle=boundaryFRIC,label='boundaryMat'))
O.materials.append(FrictMat(density=2500,young=5.5e10,poisson=0.25,frictionAngle=spFRIC,label='sphereMat'))

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

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

sp1.makeCloud((0,height+1*RADIUS1,width),(length,2*height-1*RADIUS1,2*width), rMean=RADIUS2, periodic=True, seed =1)
sp2.makeCloud((0,height+0.1*RADIUS1,width),(length,height+0.1*RADIUS1-1e-10,2*width), rMean=RADIUS2, periodic=True, seed =1)
sp3.makeCloud((0,2*height-0.1*RADIUS1,width),(length,2*height-0.1*RADIUS1-1e-10,2*width), rMean=RADIUS2, 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()],
  [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True)]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,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')
 ,PyRunner(command='fixVelocity(RATE_shear)',iterPeriod=1,label='fixVel',dead=True)
 ]

def triaxDone():
 triax.dead=True
 O.pause()

O.run(10000000000,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()
  boundaryVel=copysign(min(0.1,abs(0.5*(currentSN-SN))),currentSN-SN)
  O.bodies[0].state.vel[1]=boundaryVel
  if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
   stage+=1
   fnPlaten=O.forces.f(0)[1]
   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(100,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(10000000000,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)

newton.damping=DAMPSHEAR
fixVel.dead = False

initial_particle_pos = O.bodies[topLayer_id[0]].state.pos[0]

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

O.run(1000000000,1)

### Micro-strain analysis
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

Karol Brzezinski (kbrzezinski) said : #3

Hi,

I tried to run your code and I couldn't save the file at al (Ubuntu 20.04 Yade 2020.01a )l, due to the CGAL ERROR: assertion violation.

Unfortunately, I don't know how to fix this issue, but I have some comments:
- You set something like "O.bodies[i].state.dynamics = False", I think it should be "O.bodies[i].dynamic = False" as in the example [1]
(BTW I noticed that Body.dict() doesn't show the "dynamic" value. I thought it should.)
- I am not sure, but I think that periodic boundaries may be problematic. Hypothetically: If the particle in state 0 is on the right edge and moves right it will appear on the left side (state 1). If the displacements are calculated based on positions (I am not sure about implementation) the result will be peculiar.
- If your only problem is excluding some bodies from tesselation you could save positions from state 0 and 1 only for desired bodies in external files and then compute deformations using those files instead of states (defToVtkFromPositions() instead of defToVtk() [2])

Best wishes,
Karol

[1] https://yade-dem.org/doc/user.html#motion-constraints
[2] https://yade-dem.org/doc/yade.wrapper.html?highlight=tesselationwrapper#yade.wrapper.TesselationWrapper.defToVtkFromPositions

Chien-Cheng Hung (chiencheng) said : #4

Hi Karol,

Thank you for your reply. Those comments are very helpful.

>> You set something like "O.bodies[i].state.dynamics = False", I think it should be "O.bodies[i].dynamic = False" as in the example

Thanks for pointing this out!

>> - I am not sure, but I think that periodic boundaries may be problematic.

I am not sure if the periodic boundaries could be the reason. But I guess the example offered in the user's manual [1] also applies periodic boundaries. For my current simulation, I want to see if I can see any high-strain zone using thus micro-strain methods.

>> If your only problem is excluding some bodies from tesselation, you could save positions from state 0 and 1 only for desired bodies in external files and then compute deformations using those files instead of states (defToVtkFromPositions() instead of defToVtk() [2])

This seems a feasible way to do it. But I couldn't find any example of using defToVtkFromPositions().
Could you perhaps show me an example about using it?

Cheers,
Chien-Cheng

Best Karol Brzezinski (kbrzezinski) said : #5

Hi,

here you have the two approaches compared:

####
centers = []
for x in range(3):
 for y in range(3):
  for z in range(3):
   centers.append((x,y,z))

O.bodies.append([
 sphere(center=center,radius=1,) for center in centers
])

sp = SpherePack()
sp.fromSimulation()

TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()

sp.save('sp0.txt')
TW.setState(0)

shiftBodies(list(range(18,27,1)),(0,0,0.5)) # move one layer of spheres without engines (I am just being lazy)

sp.fromSimulation()
sp.save('sp1.txt')
TW.defToVtkFromPositions( 'sp0.txt','sp1.txt' , 'strain_from_files.vtk')

TW.setState(1)
TW.defToVtk("strain_from_memory.vtk")

############

Cheers,
Karol

Hi,
I confirm TesselationWrapper may have issues with space periodicity, but not because of this:

> " Hypothetically: If the particle in state 0 is on the right edge and moves right it will appear on the left side (state 1)"

It would not appear on the left side in terms of position. Particles leaving from right and appearing on left is just a visualization trick. In memory positions remain untouched, they are all over the place instead of staying in their initial period. Hence the problem: in order to build periodic triangulation the positions should be packed into a single period, first, then additional clones should be added all around that period to close the domain. That's what we do in PeriodicFlowEngine.

Alternatively, we could use PeriodicTriangulation from CGAL directly. They implemented it after the integration in yade (partly on our request...), and it is not used currently. For the moment I don't see an out-of-the-box solution with TesselationWrapper + periodic BCs.

Bruno

Chien-Cheng Hung (chiencheng) said : #7

Hi Karol,

Thanks for showing me an example.

In your example, why do you add sp.fromSimulation() after shiftBodies(...)?
I can see the vtk file won't work if removing the sp.fromSimulation(), but I don't understand the use of it.
I am wondering if I am not using spherePack to generate granular packing but using ymport [1] from an exported sphere data, what I should add after shiftBodies(...)?
Thank you for your time again.

Cheers,
Chien-Cheng

[1] https://yade-dem.org/doc/yade.ymport.html?highlight=ymport#yade.ymport.text

Karol Brzezinski (kbrzezinski) said : #8

Hi,

I didn't use spherePack to generate granular packing. I "copied" the current state of simulation to the sphere pack using sp.fromSimulation(), so I could easily save positions of spheres to txt file. This was my first idea, not necessarily the best. Now I think that you could also use exportSpheres() function from export module [1].

Best wishes,
Karol

[1] https://yade-dem.org/doc/yade.export.html

Chien-Cheng Hung (chiencheng) said : #9

Thanks Karol Brzezinski, that solved my question.