Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Bruno Chareyre]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Status: Open => Answered

Bruno Chareyre proposed the following answer:
> How do you get the normal component of the relative velocity?

Something like this (it would need another term from the velocity gradient, you 
can get inspiration from scgeom code):
(state2->vel - state1->vel).dot(normal) 

>from a physical point of view I believe that the question of the
contact point is not that simple and is fundamental.

Of course it is not really half-way between non-symmetric objects (different 
sizes, different properties, ...) but the real question is where should it be 
if it is not half-way?
Fundamentally, in fact, the notion of contact point/force is in itself 
irrelevant, that is a reason why I don't like very much the grammar " c1x = 
(geom.contactPoint - de1.pos)". I prefer to write equations without contact 
points, even if of course it can be mathematically equivalent.

The difference between the two sets of equation is that the first one uses 
positions, the second doesn't.
The problem in the first one is this (I didn't realize immediatly): if one body 
is larger than one period it is not clear how to define I->cellDist since it is 
not unique. If positions are not used the problem disappear.

> but this would mean keeping all the compatibility with the existing
ViscElPhys_Basic law (prescription possible of (young, poisson) or (kn,
ks) with (en, et) or (cn, cs))

"Merging" A and B does not mean to break B to make it fit in A by brute force. 
;)
Of course the point is to keep the good stuff, else no point merging.

I guess the straight way is to make a viscous version inheriting from 
CundallStrack. It may need to reorganize a little bit the code of CundallStrack 
to let sections of the go() function be used as standalone functions (mostly a 
way to avoid code duplication in the inheriting classes).
Not sure yet if the law functor is enough or if some Ig/Ip functors would need 
something as well.

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Raphaël Maurin]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Status: Answered => Open

Raphaël Maurin is still having a problem:
Hi Bruno,

Thanks for your answer, there are still some points that are unclear to me, I 
probably didn't understand everything:
- How do you get the normal component of the relative velocity only from the 
shear one ? I would need to evaluate the relative velocity for that and if I 
want to account for the ratcheting effect correction,  it seems to me that I 
need to call again getIncidentVel, maybe I missed something. 

>I'm not following you fully. 1/ I don't see a trick/approximation, to
me it is just ordinary code. So when you ask "why" I tend to just think
"why not?". 2/ why do you think it is irrelevant to include the 0.5*un
term in the branch between e.g. sphere and box or spheres of different
sizes? How would you like it? Did you realize that the same 0.5 is also
used between facets and spheres whatever the law functor?

>From what I understand the contact point is then always considered to be at 
>half the penetration depth. Two things. 1/Why not from a modeling point of 
>view as you say, but from a physical point of view I believe that the question 
>of the contact point is not that simple and is fundamental. If you consider an 
>interaction between two objects that have a different curvature, the 
>definition of the contact point is not straightforward and you should consider 
>the position of the contact point which gives you the right physical behavior 
>for the normal and tangential contact force. Defining it as half the 
>penetration depth is a choice (which I do not know if it is well reproducing 
>the "physical"/expected behavior). 2/ If you say that the contact point is 
>always considered as half the penetration depth, then I do not understand why 
>that:
const State& de1 = *static_cast(bodies[id1]->state.get());
const State& de2 = *static_cast(bodies[id2]->state.get());
const Vector3r shift2 = scene->isPeriodic ? 
scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero();
const Vector3r c1x = (geom.contactPoint - de1.pos);
const Vector3r c2x = (geom.contactPoint - de2.pos - shift2);
gives different results from that:
const Vector3r c1x = geom->radius1-0.5*geom->penetrationDepth*geom->normal;
const Vector3r c2x = -geom->radius2-0.5*geom->penetrationDepth*geom->normal;

?

> However, guess why ViscElPhys_Basic has a bug that CundallStrack does
not have? Because they are two different functors instead of one, i.e.
more lines of code with less eyes/fingers for each of them. Hence if
someone wants to do something very useful for the years to come he could
merge some of the existing functors. That was my point.

I am ok with that, but this would mean keeping all the compatibility
with the existing ViscElPhys_Basic law (prescription possible of (young,
poisson) or (kn, ks) with (en, et) or (cn, cs)) in the new one and
suppressing the ViscElPhys_Basic contact law. Is this possible ? In
which file(s) should it be included ?

Raphael

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Bruno Chareyre]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Status: Open => Answered

Bruno Chareyre proposed the following answer:
>  I also want the normal one

Very easy to find. Will you really waste cpu time in calculating the
incident shear twice when you actually only need the normal one (you
already have the shear one with my suggestion).

>contact point is always considered to be situated at the middle of the
penetration depth. While this is ok for particles of the >same size,
this to my opinion does not apply when considering the interaction of a
particle with a wall. My question was (and >still is): why is it
necessary to use this trick/approximation to avoid the bug ?

I'm not following you fully. 1/ I don't see a trick/approximation, to me
it is just ordinary code. So when you ask "why" I tend to just think
"why not?".  2/ why do you think it is irrelevant to include the 0.5*un
term in the branch between e.g. sphere and box or spheres of different
sizes? How would you like it? Did you realize that the same 0.5 is also
used between facets and spheres whatever the law functor?

> I believe it is not a good idea to leave this contact law with a
mistake

You are right there is nothing wrong in fixing existing code and you are very 
welcome to do it.
However, guess why ViscElPhys_Basic has a bug that CundallStrack does not have? 
Because they are two different functors instead of one, i.e. more lines of code 
with less eyes/fingers for each of them. Hence if someone wants to do something 
very useful for the years to come he could merge some of the existing functors. 
That was my point.
Cheers
Bruno

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Raphaël Maurin]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Status: Answered => Open

Raphaël Maurin is still having a problem:
Hi Bruno,

Not sure I completely agree with you on all what you are saying. 
- Regarding getIncidentVel, I thought at using shearInc but this gives me only 
the shear component and I also want the normal one so it does not work. I need 
to use getIncidentVel
- Regarding the periodicity,  it is implemented in 
Law2_ScGeom_ViscElPhys_Basic. Indeed, the shift which is necessary to evaluate 
the geometry of the contact between two particles 
(scene->cell->intrShiftPos(cellDist), usually called shift2 in the code) and 
that you are describing, is implemented in it. 
However, the problem does not seem to rely on the periodicity, at least for 
spheres, but on the periodicity for interaction between spheres and other 
objects. This can be very well identified by comparing the contact law with 
Law2_ScGeom_FrictPhys_CundallStrack. In the latter case, a "trick"/an 
approximation is used when considering periodic boundary conditions. It is the 
one I reported: 
if (!scene->isPeriodic && !sphericalBodies) {
  State* de1 = Body::byId(id1,scene)->state.get();
  State* de2 = Body::byId(id2,scene)->state.get();
  applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, 
id1, de1->se3.position, id2, de2->se3.position);}
 else {//we need to use correct branches in the periodic case, the following 
apply for spheres only
  Vector3r force = -phys->normalForce-shearForce;
  scene->forces.addForce(id1,force);
  scene->forces.addForce(id2,-force);
  scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* 
geom->normal.cross(force));
  scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* 
geom->normal.cross(force));
 }
Where we can see that in the case where the conditions are periodic the first 
condition is never fulfilled, and the part following the else is used. In that 
case, what is applied is that the contact point is always considered to be 
situated at the middle of the penetration depth. While this is ok for particles 
of the same size, this to my opinion does not apply when considering the 
interaction of a particle with a wall. My question was (and still is): why is 
it necessary to use this trick/approximation to avoid the bug ? 
- For the fact to implement a potential correction (which I will probably do) 
in CundallStrack, I am not sure it is relevant, especially considering that to 
my knowledge a non-negligible number of people use this contact law and I 
believe it is not a good idea to leave this contact law with a mistake. I could 
simply adjust the formulation to make it similar to CundallStrack, in order to 
have a coherence and correct the mistake (if there is any real one).

Cheers
Raphael

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Bruno Chareyre]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Status: Open => Answered

Bruno Chareyre proposed the following answer:
Hi,
In short you found that Law2_ScGeom_ViscElPhys_Basic does not implement 
periodicity and that it uses inexact expressions of the relative displacement.
I can confirm that. :)

The problem with the branch vector is because when two computational particles 
interact they have centers which can be arbitrarily far away - say, in 
different periodic windows of the periodic space. A consistent interaction 
geometry thus requires in principle to shift one of them to the same period 
where the other one is (so that (center - contactPoint) is a consistent branch 
vector). It is done for instance in [1]:
applyForceAtContactPoint(-phys->normalForce-shearForce, 
geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + 
(scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): 
Vector3r::Zero()));

However if the branch vector can be defined independently (like with
just radius and normal, the case of CundallStrack) there is no need for
a shift. The two approaches are mathematically equivalent.

You can access interaction->geom->getIncidentVel(...) where CundallStrack 
access interaction->geom->shearInc.
However shearInc is a cached value while getIncidentVel triggers complex 
calculations and needs input parameters (with possible mistakes); so I would 
suggest to use shearInc/scene->dt instead of getIncidentVel(...).

If improvements are in view, it may be much simpler and less energy-dissipative 
to introduce the viscous term directly into the CundallStrack family of 
functors (and leave ViscElPhys_Basic behind as it is). Since you only found 
"some" of the problems with this law functor.
Cheers
Bruno

[1]
https://github.com/yade/trunk/blob/master/pkg/dem/CohesiveFrictionalContactLaw.cpp#L163

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Raphaël Maurin]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Raphaël Maurin gave more information on the question:
Hi,

Even though I am not quite sure to understand everything, I found  the
origin of the problem by comparing the Law2_ScGeom_ViscElPhys_Basic to
the Law2_ScGeom_FrictPhys_CundallStrack. In the following of the
message, I focus only on the elastic part of
Law2_ScGeom_ViscElPhys_Basic, considering a case with zero viscous
damping.

When it is not periodic, there are not much difference in the evaluation of the 
normal and shear elastic contact force between the two contact laws, even if 
sign conventions are different. There is only a difference in the application 
of the contact force to the interacting elements. 
In the Law2_ScGeom_ViscElPhys_Basic, force and torques are evaluated and 
directly added to the particles following:
force = phys.normalForce + shearForce + shearForceVisc;
torque1 = -c1x.cross(force);
torque2 = c2x.cross(force);

with c1x = (geom.contactPoint - de1.pos) defined between the position of the 
particle (de1.pos) and the contact point.
and then

addForce (id1,-force,scene);
addForce (id2, force,scene);
addTorque(id1, torque1,scene);
addTorque(id2, torque2,scene);


In the  Law2_ScGeom_FrictPhys_CundallStrack it is done as: 
if (!scene->isPeriodic && !sphericalBodies) {
State* de1 = Body::byId(id1,scene)->state.get();
State* de2 = Body::byId(id2,scene)->state.get();
applyForceAtContactPoint(-phys->normalForce-shearForce, 
geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);}
else {//we need to use correct branches in the periodic case, the 
following apply for spheres only
Vector3r force = -phys->normalForce-shearForce;
scene->forces.addForce(id1,force);
scene->forces.addForce(id2,-force);

scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* 
geom->normal.cross(force));

scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* 
geom->normal.cross(force));
}

Different things regarding that:
- I didn't succeed to find the location of the function 
applyForceAtContactPoint, but I believe that it is doing exactly the same as in 
the viscous law where the force and torque are applied taking into account the 
branch vectors c1x and c2x defined with respect to the contact point. 
- We can see that when the scene is periodic (or when we consider spheres), the 
force is however not applied with respect to the contact point in the 
Law2_ScGeom_FrictPhys_CundallStrack, but with respect to the radius and the 
penetration depth. I think these two methods are completely equivalent for two 
spheres of the same size. 

If I am right with these two points, there are no differences between
the two contact laws regarding the elastic contribution in non-periodic
frameworks and when considering monodisperse particles.

However, in the periodic case, while there is no change in the
evaluation of the branch vector in Law2_ScGeom_ViscElPhys_Basic, it is
different in the Law2_ScGeom_FrictPhys_CundallStrack where it is assumed
that the second part of the if statement apply for all interactions.
While this is a priori not appropriate for all interactions (such as an
interaction between a sphere and a wall), this formulation seems to
avoid some troubles. Indeed, with the Law2_ScGeom_ViscElPhys_Basic not
modified formulation, the partices tend to spin a lot leading to
unphysical behavior.

Could someone explain me why this trick is used for periodic BC and how
it solves this spinning problem ?

Replacing the branch vector formulation based on the contact point by
the one based on the penetration depth such as done in
Law2_ScGeom_FrictPhys_CundallStrack, solves the whole problem, even when
re-activating the viscous damping.

Last remark: The way the contact force is evaluated in
Law2_ScGeom_ViscElPhys_Basic does not take into account the implemented
ratcheting effect correction, and I believe it would be better to base
it on the same formulation as the CundallStrack contact law to account
for that. I can do it, but in order to do that, I need to know how to
access the function getIncidentVel of ScGeom inside the contact law. Can
someone tell me how to do that ?

Thanks !
Cheers, 

Raphael

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Raphaël Maurin]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Status: Answered => Open

Raphaël Maurin is still having a problem:
Hi Bruno,

Thanks for the answer and for the trick to split the cylinder in 3. It
works very well when considering all element FrictMat and the
Law2_ScGeom_FrictPhys_CundallStrack. However, when I use ViscElMat and
Law2_ScGeom_ViscElPhys_Basic, I still get the same problem. Is it also
the case for you ?

A pragmatic solution is to use a cylinder that is FrictMat and interact with 
the particles with Law2_ScGeom_FrictPhys_CundallStrack, then it works. 
However, it would be nice to understand why it is not working with the 
Law2_ScGeom_ViscElPhys_Basic contact law, and try to fix this problem. 

I will check that further but if you have any ideas in the mean time, it
is more than welcome !

Raphael

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Bruno Chareyre]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Status: Open => Answered

Bruno Chareyre proposed the following answer:
My Christmas gift below. :)

The cylinders are ok in periodic boundary conditions overall. It has
been applied to sheared fiber suspensions for instance, where basically
every fiber can have different nodes in different periods (if we admit
that "being in different periods" makes any sense...). François also
generated periodic grids IIRC.

However having bodies larger than period is the door open to remaining bugs in 
interaction detection, that's most likely the source of this new problem.
If you create a cylinder with length equal to period and discretized in three 
segments then none of the segments are too large and everything is ok. And you 
can still use the "largeBody" trick to insert the floor.

Bruno 
__
from yade import pack
from yade.gridpfacet import *

dp= 6e-3

periodic = 1
cylinderBottom = 0
cylindersBottom = 1

ground = 0.2*dp
H = 60*dp
L = 20*dp

O.engines = [
 ForceResetter(),
 
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_GridConnection_Aabb()],allowBiggerThanPeriod
 = True),
 InteractionLoop(

[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Sphere_GridConnection_ScGridCoGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
# [Ip2_FrictMat_FrictMat_FrictPhys()],
# [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 PyRunner(command = 'check()', virtPeriod = 0.1),
 NewtonIntegrator(damping=0., gravity = (0,0,-9.81))
 ]

#Material creation
O.materials.append(ViscElMat(en=0.5, et=1., young=5e6, poisson=0.5, 
density=2500., frictionAngle=0.4, label='Mat'))
#O.materials.append(FrictMat(young = 5e6, poisson = 0.5, 
density=2500,frictionAngle=0.4, label='Mat'))

#Periodic Cell or containing box
if periodic:
 O.periodic = True
 O.cell.setBox(L,L,H)
 O.bodies.append(box(center= 
(L/2.,L/2.,ground),extents=(100,100,0),fixed=True,color = (0.,1.,0.),wire = 
True,material = 'Mat'))
else:
 O.bodies.append(box(center= (L/2.,L/2.,ground),extents=(100,100,0),wire = 
True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see 
below
 O.bodies.append(box(center= (0,L/2.,H/2.),extents=(0,L/2.,H/2.),wire = 
True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see
 O.bodies.append(box(center= 
(L,L/2.,H/2.),extents=(0,L/2.,H/2.),fixed=True,color = (0.,1.,0.),material = 
'Mat'))
 O.bodies.append(box(center= (L/2.,0,H/2.),extents=(L/2.,0,H/2.),wire = 
True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see
 O.bodies.append(box(center= 
(L/2.,L,H/2.),extents=(L/2.,0,H/2.),fixed=True,color = (0.,1.,0.),material = 
'Mat'))

#Bottom fixed cylinder or fixed particles
if cylinderBottom:
 n = len(O.bodies)
 cylinder(begin = (L/2.,100*L+0.001*dp,ground+dp/2.),end 
=(L/2.,-100*L+0.001*dp,ground+dp/2.),radius = dp/2.,fixed = True,color = 
(0,0,1),intMaterial = 'Mat',extMaterial = 'Mat')#Made invisible to see inside
else:
if cylindersBottom:
ep=0.001*dp #a small thing

nodes=[[L/2.,ep,ground+dp/2.],[L/2.,0.33*L+ep,ground+dp/2.],[L/2.,0.33*2*L+ep,ground+dp/2.],[L/2.,L+ep,ground+dp/2.]]

cylinderConnection([[L/2.,ep,ground+dp/2.],[L/2.,0.33*L+ep,ground+dp/2.],[L/2.,0.33*2*L+ep,ground+dp/2.],[L/2.,L+ep,ground+dp/2.]],dp/2.,fixed
 = True)
else:
for n in range(0,int(L/dp)+1):
O.bodies.append(sphere(center = 
(L/2.,n*dp,ground+dp/2.),radius = dp/2.,fixed = True,color = (0,0,1),wire = 
True,material = 'Mat'))#Made invisible to see inside

#Particle cloud for gravity deposition
partCloud = pack.SpherePack()
partCloud.makeCloud(minCorner=(4*L/10.,4*L/10.,ground+3.*dp),maxCorner=(6*L/10.,6*L/10.,0.9*H),rRelFuzz=0.,
 rMean=dp/2., num = 200)
partCloud.toSimulation(material='Mat')


O.dt = 1e-5
O.saveTmp()

#Function to check if the center of a particle is contained inside the cylinder 
or pseudo cylinder
def check():
 for b in O.bodies:
  if b.dynamic and b.state.pos[2]https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Raphaël Maurin]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Status: Answered => Open

Raphaël Maurin is still having a problem:
Hi Klaus,

Thank you for the answer. Can anyone confirm that ? 
If it is the case, is it complicated to implement ? 

Raphael

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Klaus Thoeni]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Status: Open => Answered

Klaus Thoeni proposed the following answer:
Hi Raphael,

I am unable to test at the moment but as far as I can remember from the
code periodic boundaries are not yet supported for cylinders.

Klaus

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

[Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Raphaël Maurin]

New question #429604 on Yade:
https://answers.launchpad.net/yade/+question/429604

Hi all, 

I am having a problem when I use a cylinder with periodic boundary condition. 

I made a simple script to reproduce the problem. The script is a column of 
grain falling inside a box under gravity, with a cylinder fixed on the lower 
(infinite) plane of the box. When I do that, everything works fine (option 
periodic = 0 and cylinderBottom = 1 in the script below). 
If now, I remove the sides of the box and replace it  with periodic boundary 
condition (periodic = 1 and cylinderBottom = 1 in script below), some particles 
literally goes inside the cylinder. And this whatever the stiffness of the 
interaction, the time step used, the contact law used 
(Law2_ScGeom_ViscElPhys_Basic or Law2_ScGeom_FrictPhys_CundallStrack) and the 
length of the cylinder (which I made much larger than the periodic cell).  
This does not happen if I consider instead of the cylinder a row of fixed 
spheres (periodic = 1 and cylinderBottom = 0 in the script below). 

There is also a strange behavior when defining the cylinder position and 
extents, which might be linked to the present problem:
- It is not possible to define a cylinder which extents corresponds to an 
integer of the periodic cell extent, e.g. cylinder(begin = 
(L/2.,-100*L,ground+dp/2.),end =(L/2.,100*L,ground+dp/2.), It considers 
that the length is zero at this moment, probably applying
However, if you interchange begin and end, it works, i.e. cylinder(begin = 
(L/2.,100*L,ground+dp/2.),end =(L/2.,-100*L,ground+dp/2.),
- The 3D view shows very strange position of the cylinder, which do not 
correspond to its actual position, when you consider the position of the 
cylinder nodes.

Do you have any idea what the problem could be linked to  ? 

Raphael


Script:
from yade import pack
from yade.gridpfacet import *

dp= 6e-3

periodic = 1
cylinderBottom = 1

ground = 0.2*dp
H = 300*dp
L = 20*dp

O.engines = [
 ForceResetter(),
 
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_GridConnection_Aabb()],allowBiggerThanPeriod
 = True),
 InteractionLoop(

[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Sphere_GridConnection_ScGridCoGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
#[Ip2_FrictMat_FrictMat_FrictPhys()],
#[Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 PyRunner(command = 'check()', virtPeriod = 0.1),
 NewtonIntegrator(damping=0., gravity = (0,0,-9.81))
 ]

#Material creation
O.materials.append(ViscElMat(en=0.5, et=1., young=5e6, poisson=0.5, 
density=2500., frictionAngle=0.4, label='Mat'))
#O.materials.append(FrictMat(young = 5e6, poisson = 0.5, 
density=2500,frictionAngle=0.4, label='Mat'))  

#Periodic Cell or containing box
if periodic:
O.periodic = True 
O.cell.setBox(L,L,H)
O.bodies.append(box(center= 
(L/2.,L/2.,ground),extents=(100,100,0),fixed=True,color = (0.,1.,0.),wire = 
True,material = 'Mat'))
else:
O.bodies.append(box(center= (L/2.,L/2.,ground),extents=(100,100,0),wire 
= True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see 
below
O.bodies.append(box(center= (0,L/2.,H/2.),extents=(0,L/2.,H/2.),wire = 
True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see
O.bodies.append(box(center= 
(L,L/2.,H/2.),extents=(0,L/2.,H/2.),fixed=True,color = (0.,1.,0.),material = 
'Mat'))
O.bodies.append(box(center= (L/2.,0,H/2.),extents=(L/2.,0,H/2.),wire = 
True,fixed=True,color = (0.,1.,0.),material = 'Mat'))#Made invisible to see
O.bodies.append(box(center= 
(L/2.,L,H/2.),extents=(L/2.,0,H/2.),fixed=True,color = (0.,1.,0.),material = 
'Mat'))

#Bottom fixed cylinder or fixed particles
if cylinderBottom:
n = len(O.bodies)
cylinder(begin = (L/2.,100*L,ground+dp/2.),end 
=(L/2.,-100*L,ground+dp/2.),radius = dp/2.,fixed = True,color = 
(0,0,1),intMaterial = 'Mat',extMaterial = 'Mat')#Made invisible to see inside
else:
for n in range(0,int(L/dp)+1):
O.bodies.append(sphere(center = (L/2.,n*dp,ground+dp/2.),radius 
= dp/2.,fixed = True,color = (0,0,1),wire = True,material = 'Mat'))#Made 
invisible to see inside

#Particle cloud for gravity deposition
partCloud = pack.SpherePack()
partCloud.makeCloud(minCorner=(4*L/10.,4*L/10.,ground+dp),maxCorner=(6*L/10.,6*L/10.,H),rRelFuzz=0.,
 rMean=dp/2., num = 2000)
partCloud.toSimulation(material='Mat')

O.saveTmp()
O.dt = 1e-5

#Function to check if the center of a particle is contained inside the cylinder 
or pseudo cylinder
def check():
for b in O.bodies:
if b.dynamic and b.state.pos[2]https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Re: [Yade-users] [Question #429604]: Cylinder and periodic boundary conditions

[Raphaël Maurin]

Question #429604 on Yade changed:
https://answers.launchpad.net/yade/+question/429604

Raphaël Maurin gave more information on the question:
I forgot to precise: this happens independently when using yadedaily
(2016.06a-24-0557faf~xenial)  and the latest yade version from source
code (git-54c46f3)

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

___
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp