Yade

How to compute average contact radius

Asked by 张文卿

I want to get avgrage contact radius to my next simulation.I have write C++ file of phy to get contactradius for each contact.And I have achieved it via python.Here is my script:
id1=[]
id2=[]
for i in O.interactions:
                          id1.append(i.id1)
                          id2.append(i.id2)
for k in rang (0,len(id1)):
                          Acontactradius=O.interactions[id1[k],id2[k]].phys.contactradius+Acontactradius
Acontactradius=Acontactradius/(len(id1)

But this script costs a lot of cpu to compute.It make my computer breakdown.So what if I write c++ file to get average contactradius,will it cost litter cpu?
And how can I get it by C++ file?(Ps: I have modified '_utils file'.But it don't work)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
Last query:
Last reply:
Revision history for this message
Best Jérôme Duriez (jduriez) said : 		#1

Hi,

I got somewhat confused in your question between the Python and the C++ parts, but it is clear your (Python ?) commands could be optimized. See e.g. this more "Pythonic" way to compute the average of local stiffnesses:

knList = [i.phys.kn for i in O.interactions] # note "for in in O.interactions" which is obviously the best manner to loop over interactions in YADE.. [*]
sum(knList)/len(knList) # returns the average stiffness

This takes an infinitesimal time on my PC, for 85519 interactions.

[*] https://yade-dem.org/doc/yade.wrapper.html#interactioncontainer

Revision history for this message
Jan Stránský (honzik) said : 		#2

Hello,

to get a relevant answer, please read and follow [1], I think point 6 is very close to this question. Specifically:
- please give more information on your problem. E.g. I use 1,000 or 10,000,000 particles etc.
- please provide a complete MWE such that we can try it
- I guess contactRadius is added by you in C++ files? if so please also provide the source code, the problem might be there
- information "But it don't work" has no value. Please give more information how did you modify _utils.cpp and why it didn't work (some errors? still took a lot of time?)

short source files / pieces of code can be copy-paste to the message, but larger files is better to upload somewhere and send a link

> And how can I get it by C++ file?(Ps: I have modified '_utils file'.But it don't work)

the easiest way is to write the loop inside _utils.cpp, here we are waiting to your response :-)

cheers
Jan

[1] https://yade-dem.org/wiki/Howtoask

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

There is at least one thing in the above code which is bizzare.
This:

for i in O.interactions:
                          id1.append(i.id1); id2.append(i.id2)
for k in rang (0,len(id1)):
                          Acontactradius=O.interactions[id1[k],id2[k]].phys.contactradius+...

It is literally equivalent to (in pseudo code):

for i in O.interactions:
    get(id1,id2)
for all (id1,id2):
    find i=O.interaction(id1,id2)

Which is obviously pointless.
Jérôme (#1) was mainly speaking about syntax ([i.phys.kn for i in O.interactions] is exactly like "for i in O.interactions: append" in terms of operations) but what he suggested also replaces these two loops by one.

Bruno

Revision history for this message
张文卿 (sandedadi) said : 		#4

Thanks,Bruno,Jan.I have understood where I was wrong. That helps a lot: )

Revision history for this message
张文卿 (sandedadi) said : 		#5

Thanks Jérôme Duriez, that solved my question.