cylinder larger than half of the periodic cell

Asked by Giulia Macaro on 2011-12-09

I'm trying to model a cylinder lying on a sand sample. The sample has periodic boundaries only on the planes perpendicular to the cylinder axis.

I've created a bigger periodic cell (with a very small thickness), and then, to contain the sand, 3 walls (bottom, right, and left) using the facet command. Each wall contains several facets, to overcome the "Body larger than half of the size" issues with the periodic boundaries.

My problem is now with the cylinder (generated as chainedCylinder). It's long as the thickness of the periodic cell, and it has the diameter larger than half of the size of the periodic cell.
I can divide the length of one cylinder in shorter cylinders, but I wouldn't want to divide the circumference in smaller arcs or segments to overcome the "Body larger than half of the size" issues.

Can I just comment those lines in InsertionSortCollider.cpp (see below) who give that error? The cylinder will only move horizontally and vertically, without any rotations (I'll prevent them).
So I guess if I'm using a whole single cylinder and commenting those lines, the cylinder will be reproduced in the opposite side, but exactly where the other side of the cylinder is.
These are the lines that I've commented. I'm wondering if commenting them will cause any problems.

// if(unlikely((pmn1!=pmx1) || (pmn2!=pmx2))){
// Real span=(pmn1!=pmx1?mx1-mn1:mx2-mn2); if(span<0) span=dim-span;
// LOG_FATAL("Body #"<<(pmn1!=pmx1?id1:id2)<<" spans over half of the cell size "<<dim<<" (axis="<<axis<<", min="<<(pmn1!=pmx1?mn1:mn2)<<", max="<<(pmn1!=pmx1?mx1:mx2)<<", span="<<span<<")");
// throw runtime_error(__FILE__ ": Body larger than half of the cell size encountered.");
// }

Thanks for any suggestions
Giulia

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Giulia Macaro
Solved:
2012-03-09
Last query:
2012-03-09
Last reply:
2012-03-09

This question was reopened

Try it and tell us if it works. :)
If you are very lucky, it could.

I'm not sure I understand your geometry. Looking along the cylinder's axis direction, is the cylinder's diameter larger than the size of the period?
If diameter is larger than depth but smaller than widht and height, it should not be a problem.

We will have to find a general way fix this problem one day, definitely.

Giulia Macaro (giulia-macaro) said : #2

The diameter is larger than the width of the periodic cell.
The geometry is
http://www.2shared.com/photo/LllAodFZ/frontView.html
http://www.2shared.com/photo/afOWDQeN/view.html
This is the quick version. When it works, the problem that I want to study will have a packing longer and higher than that, but with same dimensions of the particles, the cylinder and the width of the periodic cell. So basically there will be only more particles than that.

Up to now I've generated the cylinder as a unique body and it seems ok.
Should I instead generate the cylinder divided in more cylinders with same diameters and smaller length?

Thanks
Giulia

I'm trying a clean fix of this problem. The problem currentluy is that by commenting the lines, you will not crash but you may miss some interactions. Could you send me your script to test it?

The diameter is not the problem btw, since it is smaller that in-plane sizes. The problem is the total length with rounded ends.

Giulia Macaro (giulia-macaro) said : #5

Yes, you are right, the problem is that even if I divide the cylinder is shorter ones, the last section of the cylinder+the rounded end is bigger that the width of the periodic cell.

I have only created the cylinder without applying any displacement. So I have no interactions yet between it and the particles.

I'm sending the script.

It works correctly now (bzr2984). You can enable large bodies with a new collider's flag (else it returns an error as before).
I'm sending you a script where the facets are replaced by one unique large box.

There is another problem in your script: the cylinder has no bounding box. This is because a chain needs a begin and a end, while you have only a begin, so the bound is undefined. See example scripts for ending the chain.

Giulia Macaro (giulia-macaro) said : #7

Thank you very much.

I am wondering if there is any difference in using the bottom wall as facets or as a box.
I initially created it as facets to avoid that problem with the periodic cell, but now does it make any difference?

Giulia Macaro (giulia-macaro) said : #8

Thanks Chareyre, that solved my question.

It is more convenient: you don't need to write in your script boring loops creating a list of facets when all you need is a plane.
It runs faster: less bodies in total in the simulation.
It has no bug due to the double contacts on edges between facets.

Giulia Macaro (giulia-macaro) said : #10

I still have some problem with the cylinder. The bounding box for a chained cylinder is not implemented for periodic simulations.
There is only the if(!scene->isPeriodic) option:
 if(!scene->isPeriodic){
  const Vector3r& O = se3.position;
  Vector3r O2 = se3.position+cylinder->segment;
  aabb->min=aabb->max=O;
  for (int k=0;k<3;k++){
   aabb->min[k]=min(aabb->min[k],min(O[k],O2[k])-cylinder->radius);
   aabb->max[k]=max(aabb->max[k],max(O[k],O2[k])+cylinder->radius);
  }
  return;
 }
I guess I need to specify the position in relation to the periodic cell.. something like
 const Vector3r& O = scene->cell->unshearPt(se3.position);
 Vector3r O2 = scene->cell->unshearPt(se3.position)+cylinder->segment;
Is that correct? It seems working now, but I don't know if I'm forgetting something else.

I've got another question: for my simulations I need a rigid cylinder with a circular lateral surface.
I have chosen to use the chainedCylinder instead of the facetCylinder because I don't want to use the non-circular surface given when then cylinder is created with facets. But at the same time I need a non-deformable cylinder..
So I'm trying to get rid of the deformable chained state, like blocking the rotations of each part of the chain and the translations in the direction parallel to the axis of the cylinder with O.bodies[].state.blockedDOFs='xXYZ'. Would it be enough, or should I include a condition imposing the same lateral and vertical displacement for each body of the chain?

Thanks for your suggestions,
Giulia

Hi Giulia,

1/ I added the "if(!scene->isPeriodic)" at a time when there was no
option to have bodies larger than period.
In this special case, the large body is not periodic (else it would
touch itself). So, you could just comment the "if" line and it should
work (just, you must define the center of cylinder so that it is in
period (0,0,0)), no need to wrap/shear/unshear.
Could you try and confirm?

2/ For the second question, make a chain with only one cylinder, it will
be perfectly rigid.
Blocking DOFs works also but then the cylinder can't be dynamic.

Note that 1/ applies when there is only one cylinder. It is a bit less
likely to work in all cases if the chain is composed of many
smaller-than-period cylinders. Note sure, you will have to experiment.

Bruno

Giulia Macaro (giulia-macaro) said : #12

Hi Bruno,

1/ I've commented and it seems working: I have the bounding box and the interactions between cylinder and spheres. It works even if I have a chain with many cylinders.
What do you mean when saying that the centre of the cylinder have to be in (0,0,0)?
(0,0,0) isn't the origin of the periodic cell?

2/ A chain with one cylinder means without the loop for the chain, right? Then I have to close the chain anyway, is that correct?
But if I try doing it, in the 3d view I've got a wire cylinder in my domain, and a full cylinder in the periodic space.. and I can't really explain why.

Thanks,
Giulia

> >What do you mean when saying that the centre of the cylinder have to be in (0,0,0)?
Sorry, I meant an integer vector defining a period. We generate and plot
bodies in period (0,0,0), but they can move during the simulation to
other periods (say (0,+1,0) if it goes up and reach the next period).
It terms of cartesian coordinates, the cylinder must be in the box
defined by the origin (0,0,0) and hSize. It should not be at the origin,
I agree.

> >2/ A chain with one cylinder means without the loop for the chain, right? Then I have to close the chain anyway, is that correct?
> But if I try doing it, in the 3d view I've got a wire cylinder in my domain, and a full cylinder in the periodic space.. and I can't really explain why.
Oh yes, you are right, one cylinder means two nodes.
The wire cylinder is probably just because of the display trying to show
the same body in the adjacent period (you have this for spheres to).
Try turning off "ghosts" in the Display tab, and it should disappear.
Note that the length of your two-nodes cylinder can be arbitrary high,
it will not hurt. In fact you should make it very long so that spheres
moving in different periods will never touch the ends of it.

Bruno

Giulia Macaro (giulia-macaro) said : #14

I didn't explain it correctly: the wire (ghost) cylinder is also in the real domain, and the full real cylinder sometime is somehow shifted in the periodic space (where I'd expect to see only the ghost cylinder instead).
It does not happen when I run the same script without commenting the if(isPeriodic) in the Aabb for the cylinder. So I'm thinking that, even if I have interactions spheres/cylinder, commenting that line might not be the best choice. And also, now the aabb is created only on one side of the period (the one with positive coordinates, starting fomr (0,0,0)).

I'm a bit confused: the cylinder can be bigger than the cell (allowBiggerThanPeriodicCell=True) but it is still periodic, isn't it? This means that if the cylinder is longer than the width of the periodic cell, it is replicate on the other side. What happens then to the interaction cylinder vs. ghost cylinder?

This what should happen normally:
- for normal bodies, the collider is trying to find if there is a combination of shifts so that two bodies can interact
- if you create a cylinder longer than the period, then it will seen as not-periodic
- it will not be shifted (except visually if "display ghost" = True) and another body will always interact with it according to its absolute position in space
- it will not interact with itself

For instance, I made tests with sphere on inclined plane 10 times larger than period. The sphere was rolling on the plane, passing in the period cyclically. Then after it rolled on 10 period, it fells down because it reached the end of the plane.
It is a bit strange to imagine a non-periodic body in periodic space I agree, but it is in fact just an easy way to simulate an infinite body. That is why you want a very long cylinder, so that particles will never reach the ends of it.
We may find a better way so that any body larger than period is really seen as infinite, but it is not done yet (and a bit more difficult).

Maybe send a script to show unexpected behavior if it doesn't work as expected.

I forgot to ask: what version are you using? Because I made changes very recently (bzr3035).

Giulia Macaro (giulia-macaro) said : #17

I'm using 3036, but I've commented the 'if(scene->Periodic)' condition in the Bo1_ChainedCylinder_Aabb (line 328).
I'm sending you my script.

1. Thanks for the detailed explanation. The example of the sphere on the inclined plane is quite useful.
But how many cells are reproduced during the simulation? Is it an infinite number? I assume that, if I have a sphere falling under gravity inside a periodic cell without any plane, it keeps on going out and in the primary cell, in loop for ever, is that right?

2. I've run some more tests, and now it works fine with a large number of chains (Ne=30 in generatePipe() of my script) and blocking the rotations and lat translations.
If then I make the cylinder with only one chain (Ne=1), and still block xXYZ, the cylinder doesn't move (as expected). Plus it seems to be periodic only on one side.
If I remove the constrain to the DOFs, and make the cylinder with only one chain, I've Segmentation Fault.
Could I just use the long chain with blockDOFs='xXYZ'? The cylinder seems to behave still dynamically..

3. Is there any difference between creating a periodic cell (O.periodic=True and O.cell.setBox()) and creating a periodic packing (particleSD(periodic=True))?
My script works fine with O.periodic=True and particleSD(periodic=False), but gives weird behaviour when also particleSD(periodic=True).

4. In the periodicSandPile example, it's not clear to me why the bottom floor is not periodic. In my simulation I always have a periodic bottom floor, and so the cylinder. How do I make them non-periodic?

Thank you

1. Yes, right.
2. I recommend unique cylinder (+ending node). Why do you say it seems periodic on only one side?
3. Yes, periodic BCs and and periodic cloud generation are independant things. I don't know particleSD. With makeCloud it works correctly (see the periodicSandPile example)
4. It is not periodic in the x,z directions because it's bigger than period, you have no choice here. It is still periodic in y (if one particle go through the floor, it will hit it again in the (0,-1,0) period).

Please send scripts showing the problems. You should have a attach button in the question interface.

Giulia Macaro (giulia-macaro) said : #19

I'm sorry, but I don't see any attach button.. I think it's only for the bugs pages. I've sent the script to your email yesterday. I'm sending it again.

2. I mean that the ghost cylinder (unique cylinder and blockDOFs='xXYZ') is only on one side, at least this is what I see from the 3D view. Instead the unique cylinder without blockDOFs gives segmentation fault or sometime more than one cylinder generated outside the cell.
3. What is then their difference? I can't find an exhaustive answer in the documentation. I have different results even if I use makeCloud and periodic=True or False. In periodicSandPile, instead, I don't see any difference if I change periodic=True/False in makeCloud.
4. I still don't get it. In my script, even if I have allowBiggerThanPeriod=True, I have ghost parts of the bottom floor and of the cylinder. Is it reasonable?

I found your script, and I understand a bit more.
With only one cylinder and fixed xXYZ, it seems to work correctly.
The 3D display is misleading, because it will always display the
position of the long cylinder inside the base period. And the "position"
of a cylinder is at a node, not a mid-range between nodes. If you impose
a velocity on x to the cylinder, you will see this effect clearly.
This display problem aside, I think the simulation is correct. I'm
sending you a script that - I think - gives the correct behaviour.
I have only one doubt on the ending node, which is small and may
therefore be seen as a periodic object. I'm not sure it is a problem but
I have not much time to investigate.
For sure, you should not trust what you "see" in the 3D view, instead
you should inspect the results of the simulation (contact forces, etc.)
and see if it is consistent.

Bruno

On 08/03/12 17:40, Giulia Macaro wrote:
> Question #181411 on Yade changed:
> https://answers.launchpad.net/yade/+question/181411
>
> Status: Answered => Open
>
> Giulia Macaro is still having a problem:
> I'm sorry, but I don't see any attach button.. I think it's only for the
> bugs pages. I've sent the script to your email yesterday. I'm sending it
> again.
>
> 2. I mean that the ghost cylinder (unique cylinder and blockDOFs='xXYZ') is only on one side, at least this is what I see from the 3D view. Instead the unique cylinder without blockDOFs gives segmentation fault or sometime more than one cylinder generated outside the cell.
> 3. What is then their difference? I can't find an exhaustive answer in the documentation. I have different results even if I use makeCloud and periodic=True or False. In periodicSandPile, instead, I don't see any difference if I change periodic=True/False in makeCloud.
> 4. I still don't get it. In my script, even if I have allowBiggerThanPeriod=True, I have ghost parts of the bottom floor and of the cylinder. Is it reasonable?
>

--
_______________
Bruno Chareyre
Associate Professor
ENSE³ - Grenoble INP
11, rue des Mathématiques
BP 46
38402 St Martin d'Hères, France
Tél : +33 4 56 52 86 21
Fax : +33 4 76 82 70 43
________________

Giulia Macaro (giulia-macaro) said : #21

Thank you very much. I'll run more tests starting from that script and will investigate the results obtained.

Giulia