How to make some stress plates consist of 3 circles (that have radius gradually (0.5,1,2)?

Hello,

I am sorry, but I could not understand at all, what you want to achieve..
Maybe you could provide a link to a picture or scheme?

> I am trying to simulate 3 circles of stress plates that have radius dimension gradually, which has relation with contact force and the area of each circle.

what is "circle of stress plate"?
what does "having radius dimension gradually" mean?
what is "contact force"?

> In order to apply the force that corresponds to each area of circles, I need to have 3 stress plates with radius 0.5,1, and 2 in XY plane in the center of the facet.

what is "the facet"?

A few notes to the code:

> def stress_rad1():
> rad1=O.bodies[3011].shape.radius

I got "IndexError: Body id out of range."

> O.engines=O.engines+[PyRunner(command='stress_rad1()',iterPeriod=50)]

Try to avoid this structure, put the PyRunner directly to O.engines=[...] definition

> PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),

using iterPeriod or virtPeriod is preferable over realPeriod

> global boulder
> boulder=O.bodies[0]
> O.bodies[0].state.mass=2000

"global boulder" is ugly, not necessary and misleading here.
You can make a global 'boulder' variable
boulder = O.bodies[0]
after sp.toSimulation(), without "ugly" global keyword.
Or just use O.bodies[0], is not much longer and maybe even more clear what it is and you are already using it this way (inside checkUnbalanced) anyway..

> def unloadBoulder():
> if abs(O.forces.f(boulder.id)[0]):
> O.run()
> checker.command='stress_rad1()'

seems like it is not needed, O.run should not be there (the simulation is running already), so you can set checker.command='stress_rad1()' directly in checkUnbalanced

cheers
Jan

1. The following link are the images about what I want to achieve https://cutt.ly/GzVIYuA (I am trying to get stress plate due to calculate the stress increment from boulder)

2. what is "circle of stress plate"?
what does "having radius dimension gradually" mean?
what is "contact force"?

I attach the link that is the images to describe what I need to achieve

3. what is "the facet"?
I attach the link that is the images to describe what I need to achieve

4. I got "IndexError: Body id out of range."

Yes, it is because the boulder which is O.bodies[3011] has not appeared yet, when the boulder has appeared, then IndexError becomes stopped.
I have tried to change it, and I hope you could help me to check, do you think it is normal now to run it?

5. -PyRunner
-realperiod
-global boulder
-O.run
-set checker.command='stress_rad1()' directly in checkUnbalanced

I have solved all by following your directions

THANK YOU

this is my existing code now:

##Sphere Cylinder pack
from yade import pack,plot,utils,export
import math
from pylab import rand #for sand color

O.bodies.append(geom.facetBox((0,0,0),(5,2,1.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((-5,-2,-1.5),(5,2,1.5),rMean=.14,rRelFuzz=.1,num=3000)
sp.toSimulation(color=(0.6+0.15*rand(),0.5+0.15*rand(),0.15+0.15*rand()))

##Define material of the grains
b=O.bodies[0]
b.state.blockedDOFs='xyzXYZ'
b.state.vel=(0,0,0)
O.materials.append(FrictMat(young=1e7,poisson=.3,density=2700,frictionAngle=18))

##Make a floor
O.bodies.append(wall((0,0,-1.5),axis=2))

##Engines and Constitutive Law
O.engines=[ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]), #move for sphere, facet, wall
  InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()], #interaction between them
   [Ip2_FrictMat_FrictMat_FrictPhys()], # ip2 list (just one list!)
   [Law2_ScGeom_FrictPhys_CundallStrack()]), # law2 list
   NewtonIntegrator(damping=.3,gravity=[0,0,-9.81]),
   PyRunner(command='checkUnbalanced()',iterPeriod=100,label='checker'),
   PyRunner(command='stress_rad1()')
   #VTKRecorder(fileName='tes',recorders=['all'],iterPeriod=100),
   ]

O.dt=.9*PWaveTimeStep()

def checkUnbalanced():
 if O.iter<1000:return # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
 O.materials.append(FrictMat(young=1e8,poisson=0.2,density=2700,frictionAngle=25))
 O.bodies.append(utils.sphere(center=(0,0,5),radius=2,color=[0,1,1]))
 O.bodies[0].state.mass=2000
 checker.command='unloadBoulder()'
 checker.command='stress_rad1()'

def unloadBoulder():
 abs(O.forces.f(boulder.id)[0])

def stress_rad1():
 rad1=O.bodies[3011].shape.radius
 b1=O.bodies[3011]
 area1=pi*rad1**2
 area2=pi*((rad1+.8)**2-rad1**2)
 area3=pi*((rad1+1.2)**2-(rad1+.8)**2)
 ForceP=0
 StP=0
 ForceP1=0
 StP1=0
 ForceP2=0
 StP2=0
 for i in range (10,3009):
  m=O.bodies[i]
  if(m.state.pos[1]>=-2) and (m.state.pos[1]<=2):
   if (m.state.pos[0]>=-1) and (m.state.pos[0]<=1) and (m.state.pos[2]<-.5):
    if (m.state.pos[0]**2+m.state.pos[1]**2)**.5<=rad1: #1st area
      ForceP=abs(O.forces.f(m.id)[2])
      StP=ForceP/area1
    if ((m.state.pos[0]**2+m.state.pos[1]**2)**.5<=rad1+.8) and ((m.state.pos[0]**2+m.state.pos[1]**2)**.5>rad1): #2nd area
      ForceP1=abs(O.forces.f(m.id)[2])
      StP1=ForceP1/area2
    if ((m.state.pos[0]**2+m.state.pos[1]**2)**.5<=rad1+1.2) and ((m.state.pos[0]**2+m.state.pos[1]**2)**.5>rad1+.8): #3rd area
      ForceP2=abs(O.forces.f(m.id)[2])
      StP2=ForceP2/area3

 if O.iter>1000:
  plot.addData(z=b1.state.pos[2],Stress=StP,Stress1=StP1,Stress2=StP2,i=O.iter,t=O.time)

globals().update(locals())

plot.plots={'i':('z'),'t':('Stress','Stress1','Stress2')}

plot.plot(subPlots=True)

O.saveTmp()

I hope you could understand what I really want to achieve
Thank you in advance
Cheers!

Thanks for more information. I am too tired now to dive into the problem, just briefly:
You want to measure time evolution of stress at certain points at the bottom plane?
Jan

thank you for your response Mr. Stránský
I will wait whenever you're ready

The important thing is the stress (that derive from Force and Area of the certain point) and the time evolution as well (like we can put in the graph for x axis is time, y is stress) so the graph will show each of stress.

Cheers!

Thank you for your help Mr. Stránský

I am using suggestion for your easiest and general approach (last part).
Apparently I am still having many failure note in terminal.
It might be from my side, don't know how to implement my code after I use your suggestion syntax.

Thank you anyway.

> Apparently I am still having many failure note in terminal.

please provide actual code and the errors.
I wrote the code "by heart" (hence the "something like" note), there might be some problems and the need of modifications. It was meant to show the core idea..

cheers
Jan

The actual code:

##Sphere Cylinder pack
from yade import pack,plot,utils,export
import math
from pylab import rand #for sand color

facets=O.bodies.append(geom.facetBox((0,0,0),(5,2,1.5),wallMask=31))
print('the ids of the facets:',facets)

sp=pack.SpherePack()
sp.makeCloud((-5,-2,-1.5),(5,2,1.5),rMean=.14,rRelFuzz=.1,num=3000)
sp.toSimulation(color=(0.6+0.15*rand(),0.5+0.15*rand(),0.15+0.15*rand()))

##Define material of the grains
b=O.bodies[0]
b.state.blockedDOFs='xyzXYZ'
b.state.vel=(0,0,0)
O.materials.append(FrictMat(young=1e7,poisson=.3,density=2700,frictionAngle=18))

##Make a floor
bottom=wall((0,0,-1.5),axis=2)
O.bodies.append(bottom)

##Engines and Constitutive Law
O.engines=[ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]), #move for sphere, facet, wall
  InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()], #interaction between them
   [Ip2_FrictMat_FrictMat_FrictPhys()], # ip2 list (just one list!)
   [Law2_ScGeom_FrictPhys_CundallStrack()]), # law2 list
   NewtonIntegrator(damping=.3,gravity=[0,0,-9.81]),
   PyRunner(command='checkUnbalanced()',iterPeriod=100,label='checker'),
   PyRunner(command='getbottomstress()')
   #VTKRecorder(fileName='tes',recorders=['all'],iterPeriod=100),
   ]

O.dt=.9*PWaveTimeStep()

def checkUnbalanced():
 if O.iter<1000:return # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
 O.materials.append(FrictMat(young=1e8,poisson=0.2,density=2700,frictionAngle=25))
 O.bodies.append(utils.sphere(center=(0,0,5),radius=2,color=[0,1,1]))
 O.bodies[0].state.mass=2000
 checker.command='getbottomstress()'

def getbottomstress():
 b1=O.bodies[2920]
 rad1=b1.shape.radius
 area1=pi*rad1**2
 interactions = bottom.intrs()
 interactions = [i for i in interactions if (i.geom.area1).norm()]
 cp = [i.geom.area1 for i in interactions]
 fs = [i.phys.normalForce + i.phys.shearForce for i in interactions]
 StP = (fs,cp)/area1

 if O.iter>1000:
  plot.addData(z=b1.state.pos[2],Stress=StP,i=O.iter,t=O.time)

globals().update(locals())

plot.plots={'i':('z'),'t':('Stress')}

plot.plot(subPlots=True)

O.saveTmp()
##############################

What I am trying to do based on my code, is that I am trying to build a circle stress plate, then I expect to get the stress by dividing the Force and the area of circle stress plate.
I am trying to find the stress by assigning the boulder (O.bodies[2920]) to the spheres, and expecting to get the stress, and generate the value by making graph.
and if this is working, I am willing to add more circle stress plate that is slightly bigger, like this area2=pi*((rad1+.8)**2-rad1**2) , and add more area3=pi*((rad1+1.2)**2-(rad1+.8)**2)

Turns out I got the failure note:
AttributeError Traceback (most recent call last)
/usr/bin/yade in <module>

/usr/bin/yade in getbottomstress()
     94 area1=pi*rad1**2
     95 interactions = bottom.intrs()
---> 96 interactions = [i for i in interactions if (i.geom.area1).norm()]
     97 cp = [i.geom.area1 for i in interactions]
     98 fs = [i.phys.normalForce + i.phys.shearForce for i in interactions]

/usr/bin/yade in <listcomp>(.0)
     94 area1=pi*rad1**2
     95 interactions = bottom.intrs()
---> 96 interactions = [i for i in interactions if (i.geom.area1).norm()]
     97 cp = [i.geom.area1 for i in interactions]
     98 fs = [i.phys.normalForce + i.phys.shearForce for i in interactions]

AttributeError: 'ScGeom' object has no attribute 'area1'

Thank you
Cheers!

The key point here is that you should know what and how you want to compute.

Randomly modifying the code (sorry, but the last code seems like that) is rarely the way to go.

So, let's take the "last option".
You know all the interactions (i.e. essentially all the forces and their positions) acting on the bottom plate.
How to convert this discrete informations to continuum-like stress?
In general, this is not unique transformation and you can have many possible ways.
One crude option is to define a region, select forces acting in the region, sum them and divide by the area of the region.
But it will be very sensitive to the cases "at the boundary" - a force just inside or just outside, a very small difference in action point could result in a jump change in the resulting stress.

So a better option might (or might not :-) be to do some "smoothing" - weighted average with weights being distances from the region boundary.

We can help you how to code it, maybe also give some suggestions, but you have do the decision what/how to compute.

> interactions = [i for i in interactions if (i.geom.area1).norm()]

i.geom is ScGeom [1]. It does not have 'area' nor 'area1'
What should the condition "if (i.geom.area1).norm()" do?

> cp = [i.geom.area1 for i in interactions]
> fs = [i.phys.normalForce + i.phys.shearForce for i in interactions]
> StP = (fs,cp)/area1

cp and fs are lists of vectors. cp is list of contact points, fs is list of forces.
What should "(fs,cp)" do?
Now, "(fs,cp)" would be tuple of lists of vectors, something like ([Vector3(...), Vector3(...), ...], [Vector3(...), Vector3(...), ...])
Dividing this by area1 would result in error

> circle stress plate

Do you want an average stress inside the entire circle?
Or average stress at the boundary of the circle?
Or ... ?

> I am willing to add more circle stress plate

once you have one circle, more circles is easy task

cheers
Jan

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

> What should the condition "if (i.geom.area1).norm()" do?
I thought to derive the contact point we introduce the area with that syntax, in order to get the stress

> What should "(fs,cp)" do?
Now, "(fs,cp)" would be tuple of lists of vectors, something like ([Vector3(...), Vector3(...), ...], [Vector3(...), Vector3(...), ...])
Dividing this by area1 would result in error

I thought it will have such a "coordinate" like x y, represented by fs and cp value then I will get the Force value.

>Do you want an average stress inside the entire circle?
Or average stress at the boundary of the circle?
Or ... ?

the entire circle, and then after adding another circle, I want to derive also average stress inside the next circle (which has a slight bigger radius)

Thank you

> AttributeError: 'ScGeom' object has no attribute 'area1'

the error message should be clear.
Always try to get rid of errors in the very first place. The rest of the function (getBottomStress) has no meaning if you have error at the very beginning of it.
As discussed, in DEM (D = discrete), everything is discrete, so interaction.geom.area (area being not so much discrete) is likely not to be available (which actually is the case).

> if (i.geom.area1).norm()

even i.geom.area1 exists (if assuming area, why area1 and not area?), what would be the meaning of "area.norm()"? And what "if area.norm()", which would always evaluate to True?

> with that syntax

apart from syntax, code has some context and semantics. First, try to understand the code before modifying it.

> I thought ...

you can always put some prints to see what the variables are (in case of no errors, of course).

Have a try. But as discussed, this is a "crude" approach, I **personally** would add some "smoothing".
###
def getBottomStress(center,radius):
    # center ... center of the stress circle
    # radius ... radius of the circle
    interactions = bottom.intrs()
    interactions = [i for i in interactions if (i.geom.contactPoint-center).norm() < radius]
    fns = [i.phys.normalForce for i in interactions]
    fnTot = sum(fn.norm() for fn in fns)
    area = pi * radius**2
    stress = fnTot / area
    return stress
    # TODO shear stress?
###

cheers
Jan

> def getBottomStress(center,radius):

On this basis you could generate stress for each (x,y) on a grid (possibly with square domains instead of circles, with square size = grid size).

Then the average on any subdomain (e.g. a circle) would be the average on the contained grid points. That's just one of many ways.

I would say you need to separate problems to help people helping you:
1/ suppose you have a list of forces and (x,y,z) coordinates of the corresponding contact points. What do you want to do of this (mathematics)?
2/ how to obtain such a list of forces and (x,y,z) coordinates of the corresponding contact points with yade (and syntax problems)

Bruno

Thanks Jan Stránský, that solved my question.