Output specific bodies using TW.setState

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

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

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

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

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

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

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

Thanks Karol Brzezinski, that solved my question.