Yade

question for the Clump with facet

Asked by Xin Wu

Dear Sir,
I have now some Problems with Clump. I will try Box-Experiment with Clump and write the Code. Butthe Prigramme stop wth warn Clump::updateProperties: NaNs in eigen-decomposition of inertia matrix?!. Can anyone help me. Thank you.

from yade import utils
from numpy import linspace
from numpy import arange
import gts
import itertools
from yade import pack

rMean=.0096
rRelFuzz=.0016
maxLoad=4500

frictionAngleSt=radians(35)
frictionAngleBo=radians(23.5)

a=0.05
id_Mat=O.materials.append(FrictMat(density=2650,young=175e6,poisson=0.3,frictionAngle=frictionAngleBo))
Mat=O.materials[id_Mat]
steel=O.materials.append(FrictMat(young=210e9, poisson=.25,frictionAngle=frictionAngleSt,density=8000))

from yade import pack, plot,qt

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(0.3,0.3,0.3250),rMean=.0096,rRelFuzz=.0016,periodic=False)
O.bodies.append([sphere(a,r,material=Mat) for a,r in sp])

print len(sp),' particles generated.'

def makeClumpTemplate():
    relRadList1 = [.007,.0036]
    relPosList1 = [[0,0,0],[0.009,0,0]]
    relRadList2 = [.0036,0.007,.0036]
    relPosList2 = [[0.0060,0,0],[0,0,0],[-0.0060,0,0]]
    templates= []
    templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
    templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
    O.bodies.replaceByClumps(templates,[.2,.3])

fIDSI=O.bodies.append(utils.geom.facetBox((.15,.15,.135),(.15,.15,.045),wallMask=15,material=steel))
fIDSII=O.bodies.append(utils.geom.facetBox((.15,.15,.045),(.15,.15,.045),wallMask=31,material=steel))

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   NewtonIntegrator(damping=0.7,gravity=[0,0,-9.810]),
   qt.SnapshotEngine(iterPeriod=60000,fileBase='/home/wuxin/Desktop/Clump/78Fr07-',label='snapshooter'),
   PyRunner(command='checkUnbalanced()',realPeriod=1,label='checker'),
   TranslationEngine(translationAxis=[1,0,0],velocity=0,ids=fIDSI,label='Transl'),
]
O.dt=.25*utils.PWaveTimeStep()

def AngVel():
    for b in O.bodies:
        if isinstance(b.shape,Sphere):
           b.state.angVel[1]=0

def checkUnbalanced():
    if O.iter<240000: return
    fIDSIII=O.bodies.append(utils.wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1,material=spheremat))
    global plate
    plate=O.bodies[fIDSIII]
    plate.state.vel=(0,0,-.025)
    O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=1000)]
    checker.command='unloadPlate()'
    checker.command='makeClumpTemplate()'

def unloadPlate():
   if abs(O.forces.f(plate.id)[2])>=maxLoad:
      O.engines=O.engines+[PyRunner(command='ServoContorl()',iterPeriod=1)]
      O.engines=O.engines+[PyRunner(command='AngVel()',iterPeriod=1)]
      checker.command='stopUnloading()'

plot.plots={'i':('w','Fz',),'PX':('Fx',)}
plot.plot()

qt.Controller()
qt.View()
r=qt.Renderer()

O.run()

utils.waitIfBatch()

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Christian Jakob (jakob-ifgt) said : 		#1

Hi,

Seems to be a bug. I will have a look at it.

Regards,

Christian

Revision history for this message
Xin Wu (wuxin2de) said : 		#2

Hallo Christian,
I am sure now that with utils.geom.facetBox the replaceByClumps cann't work.
i think the reason is that the replaceByClumps take the Facet as Sphere and register the length of the Facet as radius.
But how can I solve the Problem.
And the seconde Question how can I program the original C++ Code. I mean which Tool shall I use. And Which Steps shall I perform.

thank you very much

regards

Revision history for this message
Christian Jakob (jakob-ifgt) said : 		#3

Hallo,

> I am sure now that with utils.geom.facetBox the replaceByClumps cann't work.
> i think the reason is that the replaceByClumps take the Facet as Sphere and register the length of the Facet as radius.

No, definitely not [1]. But there is still a known bug in facets, so you should avoid facets.

> But how can I solve the Problem.

You can use boxes [2] or fixed particles or periodic space [3] as boundaries.

> And the seconde Question how can I program the original C++ Code. I mean which Tool shall I use.
> And Which Steps shall I perform.

You can look at/download the source code at git-hub [4]. You can install kdevelop, if you want to work on the code.

Some additional hints/comments:
1. Your script is working nethertheless (Just press play button in GUI). I think you can ignore the error message (I will have look at it tomorrow). Also for the testing case you can stay using facets.

2. Did you notice, that some particles fly away from the model? This is because replacing spheres by clumps will lead to some overlaps. You can avoid this using calm() function (see replaceByClumps example script at the end).

3. There is no need for typing "utils". You can remove it, if you want.

4. For the clumpTemplate you can use relative positions and radii, no need for "0.00x", just use "x" ;-)

Regards,

Christian

[1] https://github.com/yade/trunk/blob/master/py/wrapper/yadeWrapper.cpp#L292
[2] https://yade-dem.org/doc/yade.utils.html?highlight=box#yade.utils.box
[3] https://yade-dem.org/doc/yade.wrapper.html?highlight=periodic#yade.wrapper.Omega.periodic
[4] https://github.com/yade/trunk

Revision history for this message
Anton Gladky (gladky-anton) said : 		#4

2013/4/8 Christian Jakob <email address hidden>:
>> I am sure now that with utils.geom.facetBox the replaceByClumps cann't work.
>> i think the reason is that the replaceByClumps take the Facet as Sphere and register the length of the Facet as radius.
>
> No, definitely not [1]. But there is still a known bug in facets, so you
> should avoid facets.

const Sphere* sphere = YADE_CAST<Sphere*> (b->shape.get());

Hmm, maybe I am missing something, but where do you check in the
source code, that you have got really a sphere?

Something like that.

https://github.com/yade/trunk/blob/master/pkg/dem/VTKRecorder.cpp#L244

if(!(dynamic_cast<Sphere*>(Body::byId(I->getId1())->shape.get()))) continue;

Anton

Revision history for this message
Christian Jakob (jakob-ifgt) said : 		#5

> Hmm, maybe I am missing something, but where do you check in the
source code, that you have got really a sphere?

https://github.com/yade/trunk/blob/master/py/wrapper/yadeWrapper.cpp#L224

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#6

>if(!(dynamic_cast<Sphere*>(Body::byId(I->getId1())->shape.get()))) continue;

Side remark:
ClassIndex is probably faster than dynamic casts. More importantly, it really check spheres (if desirable). Something inheriting from sphere, like GridNode, will dynamic_cast correctly but it will have a different index.

if ( (*bi)->shape->getClassIndex() == Sph_Index )
https://github.com/yade/trunk/blob/master/pkg/dem/TesselationWrapper.cpp#L77

Revision history for this message
Xin Wu (wuxin2de) said : 		#7

Hallo Guys,
Can you send me a example to show me how can man program the original Code?
Thank you very much

Regards
Xin

Revision history for this message
Christian Jakob (jakob-ifgt) said : 		#8

> Can you send me a example to show me how can man program the original Code?

As already said, the easiest way to work on the code is to download it and use "kdevelop" ...
There is a project file called "yade.kdev4" (can be opened with kdevelop).

christian

Revision history for this message
Xin Wu (wuxin2de) said : 		#9

ok, Thank you

regards
Xin

Revision history for this message
Christian Jakob (jakob-ifgt) said : 		#10

Hi again,

I fixed this bug in replaceByClumps() method. See

https://github.com/yade/trunk/commit/164769edd61da2056434f890d7f5de5f1568aaaf

Now it is compatible with facets.

Sorry for this, I did not know, that factes are spherical ...

Cheers,

Christian

Revision history for this message
Xin Wu (wuxin2de) said : 		#11

hallo Christian Bruno Anton ,
thank you for your help.
I am trying now to organize Yade China Portal with the Professor from Agra University Beijing.
His name ist Xu Yong. Here ist the Contact Infor <email address hidden>

regards
Xin

Revision history for this message
Anton Gladky (gladky-anton) said : 		#12

2013/4/9 Christian Jakob <email address hidden>:
> Sorry for this, I did not know, that factes are spherical ...

You will not believe, but they are not :)

Anton

Can you help with this problem?

Provide an answer of your own, or ask Xin Wu for more information if necessary.

To post a message you must log in.