Dynamics of pfacet object does not appear to be correct

I found that the order of defining the pfacet faces is affecting the rotation. For
pfacet_nodes = [
    [nodes[0], nodes[2], nodes[1]],
    [nodes[0], nodes[1], nodes[3]]
]

The object is rotating about the z-axis. For
pfacet_nodes = [
    [nodes[0], nodes[1], nodes[2]],
    [nodes[0], nodes[1], nodes[3]]
]

The object is rotating about the y-axis. This behaviour seems very unnatural.

Hello,

I found that the engine Law2_ScGridCoGeom_FrictPhys_CundallStrack(), was part of the problem. This got rid of the rotation that was caused at the start of the contact. But after some large deformation, there are still some problems. The dynamics depending on the order of defining the pfacet. I think the error may lay with the engines that I am using. Has anywork worked with pfacet interacting with gridnode/ grid connection? I did not find an example in the documentation.

Kind regards,
Rohit John

Hello,

Please disregards, the comment I made above (#2), now I realize that I need Law2_ScGridCoGeom_FrictPhys_CundallStrack() to handle the interaction between the end of the gridConnection, gridnode flexible rod and the pfacet. Please help.

Kind regards,
Rohit John

He John,

you can find an example on how to use gridConnections and pFacets here [1]. Just as a reminder, the end nodes of a grdiConnection can rotate. Regarding symmetry, one little instability is enough to change the system behaviour. Consider using some damping in Newtonintegrator and you will see it's behaving as expected.

HTH
Klaus

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/cylinders/mikado.py

Dear Klaus Thoeni,

Thanks for your insight into the matter. I tried your solution and it worked. However, project involves studying the dynamics of this interaction, specifically the change in velocity and angular velocity of an object during an impact with an array of flexible rods (brush). Since the damping artificial [1], I believe it will affect the dynamics. Is there any other way to reduce the instability.

Do you know if this is caused by some round off error in the solver? If yes is there a way to verify if this is caused by round off?

Kind regards,
Rohit K. John

[1] https://yade-dem.org/doc/formulation.html#numericaldamping

Yes, most likely the instability comes from numerical rounding errors and the error made by the integration of the equation of motion.

Now, you are right. Newton damping is artificial but from my experience a very low value already helps and the dynamic behaviour might not be so wrong after all , unless you are interested of an accurate time scale. Alternatively, you could use viscous or hysteresis damping (or coefficients of restitution). We have the relevant contact models but not sure if they will work with grids and pFacets.

Klaus

Dear Klaus Thoeni,

Thank you for your input on my suspicions regarding the rounding error.

Is it possible for you to give a value (or a range of values) for the newton damping which worked for you? I tried different values but, I do not have a way to verify which one is close to an actual response.

> Alternatively, you could use viscous or hysteresis damping (or coefficients of restitution). We have the relevant contact models but not sure if they will work with grids and pFacets.
- Do you mean FrictViscoMat? or do you have a different model in mind?

Kind regards,
Rohit John

Hi John,

apparently there seems to be a bug in the distribution of the forces onto the gridNodes [1]. I have to look into it but maybe once fixed it might also fix part of your problem. What it won't fix is the damping/restitution coefficient issue. Several materials/contact laws in Yade support viscous damping/restitution coefficients. FrictViscoMat is one of them but it only has damping in normal direction. ViscElMat is more general or you could also use Ip2_FrictMat_FrictMat_MindlinPhys

I aim to fix the bug next week.

Klaus

[1] https://gitlab.com/yade-dev/trunk/-/issues/184

Dear Klaus,

> Apparently there seems to be a bug in the distribution of the forces onto the gridNodes [1].
- Thanks for this information. Please let me know if I can help you fix this issue in any way.

> ViscElMat is more general or you could also use Ip2_FrictMat_FrictMat_MindlinPhys.
- I shall look into this. Thanks for this information.

Kind regards,
Rohit John

Dear Klaus,

I was testing the mikado example [1], and I found an bug. This example has many vertical cylinders on a pfacet square. In [1], the lengths are random and the cylinders start in contact with the pfacet. I modified this slightly so that the cylinders have the same length and start from a height above the pfacet. The cylinders contacting the diagonal of the square bounce differently compared to the rest. Since all cylinders have the same height and dimension, I expect all of them to bounce uniformly.

Kind regards,
Rohit K. John

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/cylinders/mikado.py

# ---------------------------------------- code
# encoding: utf-8
"An example showing how two create cylinders with random length."

# from builtins import range
from yade.gridpfacet import *

#### Parameters ####
Lmin=3. # minimum length
Lmax=6. # maximum length
n=10 # number of cylinders in one direction
r=0.5 # radius of the cylinder element
phi=30. # friction angle
E=1e6 # Young's modulus

#### Engines ####
O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_GridConnection_Aabb(),
  Bo1_PFacet_Aabb(),
 ]),
 InteractionLoop([
  Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), # cylinder-cylinder interaction
  Ig2_Sphere_PFacet_ScGridCoGeom(), # needed for GridNode-pFacet interaction (why is this not included in Ig2_GridConnection_PFacet_ScGeom???)
  Ig2_GridConnection_PFacet_ScGeom(), # Cylinder-pFcet interaction
 ],
 [
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),
  Ip2_FrictMat_FrictMat_FrictPhys()
 ],
 [
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), # contact law for "internal" cylider forces
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(), # contact law for Cylinder-pFacet interaction
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack() # contact law for cylinder-cylinder interaction
 ]
 ),
 GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.1,label='ts'),
 NewtonIntegrator(gravity=(0.,0,-10),damping=0.,label='newton'),
 #NewtonIntegrator(gravity=(-1.,0,-10),damping=0.5,label='newton'),
]

#### Creat materials ####
O.materials.append( CohFrictMat( young=E,poisson=0.3,density=1000,frictionAngle=radians(phi),normalCohesion=1e10,shearCohesion=1e10,momentRotationLaw=True,label='cMat' ) ) # material to create the gridConnections
O.materials.append( FrictMat( young=E,poisson=0.3,density=1000,frictionAngle=radians(phi),label='fMat' ) ) # material for general interactions

#### Create cylinders ####
nodesIds=[]
cylIds=[]
ext=12.
dxy=ext/(n-1)
dL=Lmax-Lmin
random.seed( 10 )
for i in range(0,n):
 y=-ext/2+i*dxy
 for i in range(0,n):
  x=-ext/2+i*dxy
  L=Lmin+dL
  color=[random.random(),random.random(),random.random()]
  cylinder((x,y,3*r),(x,y,L+3*r),radius=r,nodesIds=nodesIds,cylIds=cylIds,
     fixed=False,color=color,intMaterial='cMat',extMaterial='fMat')

#### Creat ground with pFacets ####
color=[255./255.,102./255.,0./255.]
n0=O.bodies.append( gridNode([-10,-10,0],r,wire=False,fixed=True,material='cMat',color=color) )
n1=O.bodies.append( gridNode([10,-10,0],r,wire=False,fixed=True,material='cMat',color=color) )
n2=O.bodies.append( gridNode([10,10,0],r,wire=False,fixed=True,material='cMat',color=color) )
n3=O.bodies.append( gridNode([-10,10,0],r,wire=False,fixed=True,material='cMat',color=color) )
O.bodies.append( gridConnection(n0,n1,r,color=color,material='fMat') )
O.bodies.append( gridConnection(n1,n2,r,color=color,material='fMat') )
O.bodies.append( gridConnection(n2,n0,r,color=color,material='fMat') )
O.bodies.append( gridConnection(n2,n3,r,color=color,material='fMat') )
O.bodies.append( gridConnection(n3,n0,r,color=color,material='fMat') )
O.bodies.append( pfacet(n0,n1,n2,wire=False,material='fMat',color=color) )
O.bodies.append( pfacet(n0,n2,n3,wire=False,material='fMat',color=color) )

#### For viewing ####
from yade import qt
qt.View()
Gl1_Sphere.stripes=True

#### Allows to reload the simulation ####
O.saveTmp()

In the comment above, I have given the code which highlights the bug. I've modified the mikado example [1]. The newton damping is set to 0, but you can still see this bug when it is not 0. I removed the random lenghts and raised the position of the cylinders

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/cylinders/mikado.py

Thanks for that. You are right. It is even more obvious is you use spheres instead of the cylinders and and extent them over the whole area of the pfactes. The spheres in contact with the cylinder in the middle bounce differently. It might be related to this issue [1]. I will have to look into it. Thanks for reporting.

Klaus

[1] https://gitlab.com/yade-dev/trunk/-/issues/184

Dear Klaus Thoeni,

I was playing around with the pfacet and I found a discrepancy in mass of the nodes in a pfacet. I calculated the mass of the nodes using the pfacet radius, the density of the material. The mass observed (O.bodies[id].state.mass) and mass I calculated is not equal. Moreover, changing the size of the pfacet changed the mass and the mass is not equally distributed to the nodes. [1] says that "The mass of the PFacet is equally lumped into the 3 nodes." so I thought the mass should be equal for all nodes.

[1] Effeindzourou, A., Thoeni, K., Chareyre, B., & Giacomini, A. (2015). A general method for modelling deformable structures in DEM. Proceedings of the 4th International Conference on Particle-Based Methods - Fundamentals and Applications, PARTICLES 2015, October 2018, 744–754.

Kind regards,
Rohit John

Please find the code for this
# -----------------------------------------------------------------------------------------------------------
# encoding: utf-8
from yade import qt
from yade.gridpfacet import *

###########################
##### ENGINES #####
###########################

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(),
  Bo1_Wall_Aabb(),
  Bo1_PFacet_Aabb(),
 ],sortThenCollide=True),
 InteractionLoop(
 [
        Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_Sphere_PFacet_ScGridCoGeom(),
  Ig2_Wall_PFacet_ScGeom(),
        Ig2_Wall_Sphere_ScGeom()
 ],
 [
        Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),
  Ip2_FrictMat_FrictMat_FrictPhys()],
 [
        Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
  Law2_ScGeom_FrictPhys_CundallStrack(),
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack()
 ]),
    GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.8,label='ts'),
 NewtonIntegrator(gravity=(0,-9.81,0),damping=0.1,label='newton')
]

young = 1e9
density = 1e3
poisson = 0.3
friction= 30

O.materials.append(
    CohFrictMat(
        young=young,
        poisson=poisson,
        density=density,
        frictionAngle=radians(friction),
        normalCohesion=3e7,
        shearCohesion=3e7,
        momentRotationLaw=True,
        label='gridNodeMat'
        ))
O.materials.append(
    FrictMat(
        young=young,
        poisson=poisson,
        density=density,
        frictionAngle=radians(friction),
        label='gridConnectionMat'
        ))

###################################
##### PFacet creators #####
###################################

fixed = False
color=[255./255.,102./255.,0./255.]

nodesIds = []
cylIds = []
pfIds = []

#position of the node in the middle
r = 0.05

## Option 2: pfacetCreator1(vertices) -> based on 3 vertices
scale = 2

v1=Vector3(2, 0, 0) * scale
v2=Vector3(3, 0, 0) * scale
v3=Vector3(2.5, 2, 0) * scale
vertices=[v1,v2,v3]

pfacetCreator1(
    vertices,
    r,
    nodesIds=nodesIds,
    cylIds=cylIds,
    pfIds=pfIds,
    wire=False,
    color=color,
    fixed=fixed,
    materialNodes='gridNodeMat',
    material='gridConnectionMat')

print('Mass of node 0: ', O.bodies[nodesIds[0]].state.mass)
print('Mass of node 1: ', O.bodies[nodesIds[1]].state.mass)
print('Mass of node 2: ', O.bodies[nodesIds[2]].state.mass)
print('Calculated Mass: ', 4.0/3.0 * pi * r**3 * density)
#####################
##### Wall ###
#####################

O.bodies.append(utils.wall(position=-1,sense=0, axis=1,color=Vector3(1,0,0),material='gridConnectionMat'))

qt.View()
O.saveTmp()

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Dear Klaus Thoeni,

I used the fix mentioned in [1] to fix the problem I was facing. And it does fix the dynamics my simulation (see code in the first description of this problem). However, I ran into a problem when I clumped the pfacet together (I added the nodes, gridconnection and pfacets into one clump). The clump is having an unexpected rotation. I have put the code below.

I did not download a new version of yade. But I pasted the fix given in [1] into the YADE source code I am using and recompiled it. Is this a correct way to use this fix? I am using Ubuntu 18.04.5 LTS.

Kind regards,
Rohit K. John

[1] https://gitlab.com/yade-dev/trunk/-/merge_requests/628/diffs?commit_id=fc8e01415108ffd30503bbc37234d51d8d0b0499

# ---------------------------------------------------------------------------------------------------------------- code
from yade import plot, utils
from yade.gridpfacet import *

# ---------------------------------------------------------------------- defining material
bristle_radius = 1e-3
bristle_length = 0.1
bristle_seg_no = 5

bristle_young = 3e9
bristle_poisson = 0.3
bristle_density = 1000
bristle_friction = radians(30)

grid_radius = 1e-3
target_young = 30e9
target_poisson = 0.3
target_density = 10000
target_friction = radians(30)

# ---------------------------------------------------------------------- material
O.materials.append(
    FrictMat(
        young = bristle_young,
        poisson = bristle_poisson,
        density = bristle_density,
        label = 'cyl_ext_mat',
        frictionAngle = bristle_friction,
    )
)

O.materials.append(
    CohFrictMat(
        young = bristle_young,
        poisson = bristle_poisson,
        density = bristle_density,
        label = 'cyl_int_mat',

        frictionAngle = bristle_friction,
        momentRotationLaw = True,
        normalCohesion = 1e40,
        shearCohesion = 1e40,
    )
)

O.materials.append(
    FrictMat(
        young = target_young,
        poisson = target_poisson,
        density = target_density,
        label = 'pfacet_ext_mat',
        frictionAngle = target_friction,
    )
)

O.materials.append(
    CohFrictMat(
        young = 100*target_young,
        poisson = target_poisson,
        density = target_density,
        label = 'pfacet_int_mat',

        frictionAngle = target_friction,
        momentRotationLaw = True,
        normalCohesion = 1e40,
        shearCohesion = 1e40,
    )
)

# ---------------------------------------------------------------------- Engines
O.engines = [
    ForceResetter(),
    InsertionSortCollider([
        Bo1_GridConnection_Aabb(),
        Bo1_PFacet_Aabb(),
    ]),
    InteractionLoop(
        [
            Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
            Ig2_GridNode_GridNode_GridNodeGeom6D(),
            Ig2_PFacet_PFacet_ScGeom(),
            Ig2_GridConnection_PFacet_ScGeom(),
            Ig2_Sphere_PFacet_ScGridCoGeom(),
        ],
        [
            Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow = True, setCohesionOnNewContacts = False),
            Ip2_FrictMat_FrictMat_FrictPhys()
        ],
        [
            Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
            Law2_ScGeom_FrictPhys_CundallStrack(),
            Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
            Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
        ]
    ),
    NewtonIntegrator(gravity = (0,0,-10),damping = 0.0)
]
# ---------------------------------------------------------------------- bodies
# ------------------------------------- custome brush
nodesIds=[]
cylIds=[]

# bristle root location
base_pos = [
    [0.003, 0.003, 0],
    [0.003, -0.003, 0],
    [-0.003, 0.003, 0],
    [-0.003, -0.003, 0],
    ]
dz = bristle_length / bristle_seg_no # grid connection length

for i in base_pos:
    bristle = []
    for j in range(bristle_seg_no + 1):
        bristle.append([
            i[0],
            i[1],
            bristle_length - (i[2] + dz * j)
        ])
    cylinderConnection(
        bristle,
        bristle_radius,
        nodesIds,
        cylIds,
        color=[1,0,0],
        intMaterial='cyl_int_mat',
        extMaterial='cyl_ext_mat'
    )
    O.bodies[nodesIds[-1]].state.blockedDOFs='xyzXYZ'

# ------------------------------------- manaul target
origin = [
    0,
    0,
    bristle_length + bristle_radius + grid_radius
]

x_len = 0.03
y_len = 0.05
z_len = 0.05
nodes = [
    [x_len, 0, 0 + origin[2]],
    [-x_len, 0, 0 + origin[2]],
    [0.0, y_len, z_len + origin[2]],
    [0.0,-y_len, z_len + origin[2]],
]

pfacet_nodes = [
    [nodes[0], nodes[1], nodes[2]],
    [nodes[0], nodes[3], nodes[1]],
]

pnode = []
pcyl = []
pfacet = []
for i in pfacet_nodes:
    print(i)
    pfacetCreator1(
        i,
        grid_radius,
        nodesIds = pnode,
        cylIds = pcyl,
        pfIds = pfacet,
        wire = False,
        fixed = False,
        color = [0.5,0.5,0.5],
        materialNodes = 'pfacet_int_mat',
        material = 'pfacet_ext_mat',
        )

target_ids = pnode + pcyl + pfacet
target_clump_id = O.bodies.clump(target_ids)

# ---------------------------------------------------------------------- graphing and additional engines
graph_iter = 1000
plot.plots = {
    't':(
        'wx_avg', 'wy_avg', 'wz_avg',
        )}

# ----------------- initialising the variables
# These variables will be added each iteration and its average will be taken (dividing my the number of iteration) when the graph is plotted
wx_avg = 0
wy_avg = 0
wz_avg = 0
# bodyList = target_body.node_id_list
bodyList = target_ids

def get_avg_state():
    global wx_avg, wy_avg, wz_avg

    wx_temp = 0
    wy_temp = 0
    wz_temp = 0

    for i in bodyList:
        wx_temp += O.bodies[i].state.angVel[0]
        wy_temp += O.bodies[i].state.angVel[1]
        wz_temp += O.bodies[i].state.angVel[2]

    # averaging for all bodies
    no_of_bodies = len(bodyList)

    wx_temp /= no_of_bodies
    wy_temp /= no_of_bodies
    wz_temp /= no_of_bodies

    # addind to the cumulative states. This will later be divided by graph_iter to get the average

    wx_avg += wx_temp
    wy_avg += wy_temp
    wz_avg += wz_temp

def graph():
    global wx_avg, wy_avg, wz_avg

    wx_avg /= graph_iter
    wy_avg /= graph_iter
    wz_avg /= graph_iter

    plot.addData(t = O.time,
        wx_avg = wx_avg, wy_avg = wy_avg, wz_avg = wz_avg,
        )

plot.plot()
# ---------------------------------------------------------------------- plotting and simulation stop engines
end_time = 4e-1
O.engines += [
    PyRunner(command = 'graph()', iterPeriod = graph_iter),
    PyRunner(command = 'get_avg_state()', iterPeriod = 1),
]

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

Thanks for more feedback.

> I did not download a new version of yade. But I pasted the fix given in [1]

If you did it carefully it should be the ok but it is really not a recommended method (because there might be other bugfixes in the meantime, for instance).

Best method is to download the sources with git [1], then you update with "cd trunk; git pull".
Please confirm the issue with current trunk version first, else we don't know really what we are speaking about.
Bruno

[1] https://yade-dem.org/doc/installation.html#source-code

I tried you example (nice one!) and I did not see something obviously wrong.

It looks a bit like a ball-on-a-tip problem. This triangular object has to fall one way or the other since it is in very unstable equilibrium, so I'm not surprised that it rotates (which doesn't mean there is no problem, maybe you have more insight?).

Note: there is a display problem which makes it look as if the clumped objects were splitting appart. It goes back to normal when you pause the simulation. Nothing to worry about. (reason is: the cylinders are not dynamic objects and there orientation is not used anywhere in the code, so it is not updated at runtime - only when 3D view is refreshed orientation is recomputed for display).

Bruno

with newton damping 0.2 and
base_pos = [
    [0.011, 0.011, 0],
    [0.011, -0.011, 0],
    [-0.011, 0.011, 0],
    [-0.011, -0.011, 0],
    ]

The solution is stable and symmetric. My impression is you defined a very unstable situation and simulation reflects that.

Dear Bruno,

Thanks for you swift reply.
> Best method is to download the sources with git [1], then you update with "cd trunk; git pull".
- I shall try this and get back to you. In [2] it is mentioned that Ubuntu18.04 has YADE 2018. Will there be any errors when I use the latest version compiled in my computer?

> It looks a bit like a ball-on-a-tip problem. This triangular object has to fall one way or the other since it is in very unstable equilibrium, so I'm not surprised that it rotates (which doesn't mean there is no problem, maybe you have more insight?).
- If you do not clump the pfacets, triangular object remains on the bristles. I expected the simulation with the clump and one without the clump to have similar dynamics.

> Note: there is a display problem which makes it look as if the clumped objects were splitting appart.
- Thanks for your advice. It is not just this.

In the following I have put a code simulating two collisions (sphere-sphere and sphere-pfacet). I have made it so that the sphere contacts the pfacet at its centroid (centre of mass) and the pfacet body and spheres have the same mass. I expect the pfacet to move without any rotation, but it does. I have two observations.

1) The mass of the pfacet is not equally distributed to the node - this is causing an unexpected rotation. This can be solved if the pfacet is an equilateral triangle (remove the '2*' in the line Vector3([2*pfacet_circum_rad * cos(node_angle[0]), pfacet_circum_rad * sin(node_angle[0]), 0]),).

2) If you clump the pfacet together, then the collision dynamics appears to be different. But if they are not clumped together, the motion of the pfacet and sphere after collision are the same (as expected)

Kind regards,
Rohit K. John

[2] https://yade-dem.org/doc/installation.html

from yade import utils
from yade.gridpfacet import *

# ------------------------------------------------------------------- Data
sphere_young = 1e9
sphere_poisson = 0.3
sphere_friction = 30

sphere_radius = 0.1
sphere_1_position = [0, 0, 3 * sphere_radius]
sphere_2_position = [0, 5 * sphere_radius, 0]
sphere_3_position = [0, 5 * sphere_radius, 3 * sphere_radius]

pfacet_side = 0.5
pfacet_circum_rad = pfacet_side/sqrt(3)
pfacet_grid_rad = sphere_radius
node_angle = [radians(0), radians(120), radians(240)]
pfacet_nodes = [
    Vector3([2*pfacet_circum_rad * cos(node_angle[0]), pfacet_circum_rad * sin(node_angle[0]), 0]),
    Vector3([pfacet_circum_rad * cos(node_angle[1]), pfacet_circum_rad * sin(node_angle[1]), 0]),
    Vector3([pfacet_circum_rad * cos(node_angle[2]), pfacet_circum_rad * sin(node_angle[2]), 0]),
]

pfacet_centroid = (pfacet_nodes[0] + pfacet_nodes[1] + pfacet_nodes[2])/3
pfacet_nodes = [
    pfacet_nodes[0] - pfacet_centroid,
    pfacet_nodes[1] - pfacet_centroid,
    pfacet_nodes[2] - pfacet_centroid
    ]

pfacet_young = sphere_young
pfacet_poisson = sphere_poisson
pfacet_density = 1000
pfacet_friction = sphere_friction
# ------------------------------------------------------------------- Material
pfacet_int_mat = "pfacet_int_mat"
pfacet_ext_mat = "pfacet_ext_mat"

O.materials.append(
    FrictMat(
        young = pfacet_young,
        poisson = pfacet_poisson,
        density = pfacet_density,
        label = 'pfacet_ext_mat',
        frictionAngle = radians(pfacet_friction),
    )
)

O.materials.append(
    CohFrictMat(
        young = pfacet_young,
        poisson = pfacet_poisson,
        density = pfacet_density,
        label = 'pfacet_int_mat',

        frictionAngle = radians(pfacet_friction),
        momentRotationLaw = True,
        normalCohesion = 1e40,
        shearCohesion = 1e40,
    )
)
# ------------------------------------------------------------------- Engines
O.engines = [
 ForceResetter(),
 InsertionSortCollider([
        Bo1_Sphere_Aabb(),
  Bo1_GridConnection_Aabb(),
  Bo1_PFacet_Aabb(),
 ]),
 InteractionLoop([
  Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), # cylinder-cylinder interaction
  Ig2_Sphere_PFacet_ScGridCoGeom(), # needed for GridNode-pFacet interaction (why is this not included in Ig2_GridConnection_PFacet_ScGeom???)
  Ig2_GridConnection_PFacet_ScGeom(), # Cylinder-pFcet interaction
        Ig2_Sphere_Sphere_ScGeom(),
 ],
 [
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),
  Ip2_FrictMat_FrictMat_FrictPhys()
 ],
 [
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), # contact law for "internal" cylider forces
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(), # contact law for Cylinder-pFacet interaction
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), # contact law for cylinder-cylinder interaction
  Law2_ScGeom_FrictPhys_CundallStrack(),
 ]
 ),
 NewtonIntegrator(gravity=(0.,0,-0),damping=0.0,label='newton'),
    PyRunner(command = 'graph()', iterPeriod = 100)
]

# ------------------------------------------------------------------- Body
# ------------------ pfacet
pnode = []
pcyl = []
pfacet = []
pfacetCreator1(
    pfacet_nodes,
    pfacet_grid_rad,
    nodesIds = pnode,
    cylIds = pcyl,
    pfIds = pfacet,
    wire = False,
    fixed = False,
    color = [0.5,0.5,0.5],
    materialNodes = pfacet_int_mat,
    material = pfacet_ext_mat,
    )

target_ids = pnode + pcyl + pfacet
mass_node = O.bodies[pnode[0]].state.mass
volume_node = 4.0/3.0 * pi * pfacet_grid_rad**3
print('mass of node 0', O.bodies[pnode[0]].state.mass)
print('mass of node 1', O.bodies[pnode[1]].state.mass)
print('mass of node 2', O.bodies[pnode[2]].state.mass)
# target_clump_id = O.bodies.clump(target_ids)

# ------------------ sphere
sphere_density = 3.0 * mass_node / volume_node
sphere_mat = "sph_mat"

O.materials.append(
    FrictMat(
        young = sphere_young,
        poisson = sphere_poisson,
        density = sphere_density,
        label = sphere_mat,
        frictionAngle = radians(sphere_friction),
    )
)

sph_1 = sphere(sphere_1_position, sphere_radius, material = sphere_mat)
sph_2 = sphere(sphere_2_position, sphere_radius, material = sphere_mat)
sph_3 = sphere(sphere_3_position, sphere_radius, material = sphere_mat)

sph_1_ID = O.bodies.append(sph_1)
sph_2_ID = O.bodies.append(sph_2)
sph_3_ID = O.bodies.append(sph_3)
O.bodies[sph_1_ID].state.vel = [0,0,-15]
O.bodies[sph_3_ID].state.vel = [0,0,-15]
# ------------------------------------------------------------------- Other stuff
print('pfacet mass: ', O.bodies[pnode[0]].state.mass)
print('sphere mass: ', O.bodies[sph_1_ID].state.mass)
# ---------------------------------------------- plotting
from yade import plot
plot.plots = {
    't':('fz_sph', 'fz_pfacet')
}

def graph():
    t = O.time
    fz_pfacet = Vector3([0,0,0])
    for id in pnode:
        fz_pfacet = fz_pfacet + O.forces.f(id)

    fz_sph = O.forces.f(sph_2_ID)
    plot.addData(
        t = t,
        fz_sph = fz_sph[2], fz_pfacet = fz_pfacet[2]
    )

plot.plot()

# ---------------------------------------------- simulations
O.dt = 2.5e-08 # utils.PWaveTimeStep()
O.saveTmp()

Dear Bruno,

I think I found a solution. I found that setting the mass of the grdiConnections to 0 seemed to fix the problem of the triangular object rotation. It also seem to fix the example I had given above. Please let me know if I should add the whole code. I just added the lines given below after creating the pfacet.

# add this before clumping
for i in pcyl:
    O.bodies[i].state.mass = 0

Moreover, setting the mass of all the nodes to be equal also fixed the other issue (observation 1 in the above comment). Please let me know what you think.

Kind regards,
Rohit K. John

Very interesting. I'm unsure where to improve that in source code though.
Mass-inertia assignements are done mainly in some python helper function. As you found out you can always overwrite what they give (fortunately).

I would have to check again in details to explain precisely what you found but for sure it happens in one/some of these functions. Having assymetric mass with symmetric geometry is not a nice behaviour obviously. And we should improve.

Clumping is also something which uses python wrappers to update mass-inertia of the produced clump. I realize it will count the mass of cylinders whereas when it's not clumped mass of cylinders doesn't matter. So the total mass probably isn't the same before and after clumping (except when setting their mass to zero first as you suggest).

But still, if the shape is symmetric the mass of the cylinders should be symmetric too, so there is something to improve.

Bruno

Ah... the clumping step is using positions+mass to deduce total mass and center of mass while "position" of the cylinder is one of its end nodes IIRC (not the middle point as one would expect) . Hence some randomness depending on the ordering of nodes and cylinders.

It would explain why setting those masses to zero helped. So, setting cylinders mass to zero sounds like a good idea. It should only happen after assigning half of its actual mass to each of the nodes though.
The facets should do the same to avoid the "change mass by clumping" problem: compute self mass, then split in three and assign to the nodes, and finally set facet mass to zero.

If you check what's behind functions like cylinderConnection() or pFacetCreator1() you should quickly reach where such things happen and change behavior. If you are able to do that feel free to contribute it to trunk. Else maybe we should open a bug to not forget.

Bruno

Dear Bruno,

Thanks for your explanation. I will look into it. How do I open a bug and what information should I put into it?

Kind regards,
Rohit

Dear Bruno,

I have found one more bug. The angular momentum value of the clump returns (0,0,0) even though inertia is not 0 and angular velocity is not 0. This happens when I set the mass of the gridConnections to 0. However, if I do not do this, I get a non zero value for the angular momentum.

In the simulation, two sphere moving the opposite direction. They are not collinear, they are offset. One is offset in the +y while the other in the -y direction. The are set to contact a pfacet cube at the same moment of time on opposite faces. I wanted to check if the angular momentum was conserved. It turns out that it is not conserved.

I am using ubuntu 18.04. I found this error in yade-2021-03-22.git-096f9a3, yade-2020-11-05.git-7d0ae90
 and Yade 2018.02b. All of them showed this bug

Kind regards,
Rohit K. John

#---------------------------------------------------- GTS FILE. File name = 'cube.gts'
14 36 24 GtsSurface GtsFace GtsEdge GtsVertex
0.5 0.5 0.5
0.5 0.5 -0.5
0.5 -0.5 0.5
0.5 -0.5 -0.5
-0.5 0.5 0.5
-0.5 0.5 -0.5
-0.5 -0.5 0.5
-0.5 -0.5 -0.5
0.5 0.0 0.0
0.0 -0.5 0.0
0.0 0.0 0.5
-0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 -0.5
6 8
2 6
1 2
8 7
3 4
5 6
3 7
1 3
8 4
7 5
5 1
4 2
9 3
2 9
4 9
1 9
10 7
4 10
8 10
3 10
11 3
5 11
7 11
1 11
12 6
7 12
8 12
5 12
13 2
5 13
6 13
1 13
14 8
2 14
6 14
4 14
7 21 23
7 17 20
6 25 28
9 33 36
8 13 16
3 29 32
12 14 15
5 15 13
3 16 14
9 18 19
4 19 17
5 20 18
11 22 24
10 23 22
8 24 21
4 26 27
1 27 25
10 28 26
6 30 31
2 31 29
11 32 30
2 34 35
1 35 33
12 36 34

#------------------------------------------------------------------------------------------- YADE Script
from yade.gridpfacet import *
from yade import geom, utils
from yade import plot
import sys, os
sys.path.append(".")

# ---------------------------------------------------------------------------- input parameter
# ----------------------------------------------------- target
target_young = 50e9
target_density = 1000
target_poisson = 0.3
target_friction = radians(30)

p_radius = 5e-2

# ---------------------------------------------------------------------------------------------------------- Materials
target_int_mat = 'pfacet_int_mat'
target_ext_mat = 'pfacet_ext_mat'

O.materials.append(
    FrictMat(
        young = target_young,
        poisson = target_poisson,
        density = target_density,
        label = target_ext_mat,
        frictionAngle = target_friction,
    )
)

O.materials.append(
    CohFrictMat(
        young = target_young,
        poisson = target_poisson,
        density = target_density,
        label = target_int_mat,

        frictionAngle = target_friction,
        normalCohesion = 3e100,
        shearCohesion = 3e100,
        momentRotationLaw = True,
    )
)
# ---------------------------------------------------------------------------------------------------------- Engines
O.engines = [
                ForceResetter(),

                InsertionSortCollider([
                    Bo1_GridConnection_Aabb(),
                    Bo1_PFacet_Aabb(),
                    Bo1_Sphere_Aabb(),
                ]),

                InteractionLoop(
                    [
                        Ig2_PFacet_PFacet_ScGeom(),
                        Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
                        Ig2_GridNode_GridNode_GridNodeGeom6D(),
                        Ig2_GridConnection_PFacet_ScGeom(),
                        Ig2_Sphere_PFacet_ScGridCoGeom(),
                    ],
                    [
                        Ip2_FrictMat_FrictMat_FrictPhys(),
                        Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(
                            setCohesionNow = True,
                            setCohesionOnNewContacts = False
                            ),
                    ],
                    [
                        Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
                        Law2_ScGeom_FrictPhys_CundallStrack(),
                        Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
                        Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
                    ],
                )
            ]
# ---------------------------------------------------------------------------------------------------------- objects
# ---------------------------------------------------------------------- target
(
pnode,
pcyl,
pfacet
) = gtsPFacet(
    'cube.gts',
    radius = p_radius,
    shift = (0,0,0),
    scale = 1,
    wire = False,
    fixed = False,
    color = [0.1,0.5,0.1],
    materialNodes = 'pfacet_int_mat',
    material = 'pfacet_ext_mat',
)

target_ids = pnode + pcyl + pfacet
for i in pcyl:
    O.bodies[i].state.mass = 0

target_clump_ID = O.bodies.clump(target_ids)
O.bodies[target_clump_ID].state.inertia = Vector3([0.01,0.01,0.01])
O.bodies[target_clump_ID].state.mass = 1
print(target_clump_ID)

# ---------------------------------------------------------------------- spheres
sp1 = sphere([ 0.7, 0.25, 0], 5e-2, material = target_ext_mat)
sp2 = sphere([-0.7, -0.25, 0], 5e-2, material = target_ext_mat)

sp1_ID = O.bodies.append(sp1)
sp2_ID = O.bodies.append(sp2)

O.bodies[sp1_ID].state.vel = [-1,0,0]
O.bodies[sp2_ID].state.vel = [ 1,0,0]
# ----------------------------------------------------------------------------- Additional engines
ids = target_clump_ID
O.engines += [
    NewtonIntegrator(gravity = [0,0,0], damping = 0),
    PyRunner(command = 'graph()', iterPeriod = 1000)
]

# ---------------------------------------------------------------------- plotting
plot.plots = {'t':('wx', 'wy', 'wz'), 't1':('Lx', 'Ly', 'Lz')}

def graph():
    w = O.bodies[ids].state.angVel
    pos_sp1 = O.bodies[sp1_ID].state.pos
    vel_sp1 = O.bodies[sp1_ID].state.vel
    angMom_sp1 = O.bodies[sp1_ID].state.angMom
    mass_sp1 = O.bodies[sp1_ID].state.mass
    L_sp1 = angMom_sp1 + mass_sp1 * pos_sp1.cross(vel_sp1)

    pos_sp2 = O.bodies[sp2_ID].state.pos
    vel_sp2 = O.bodies[sp2_ID].state.vel
    angMom_sp2 = O.bodies[sp2_ID].state.angMom
    mass_sp2 = O.bodies[sp2_ID].state.mass
    L_sp2 = angMom_sp2 + mass_sp2 * pos_sp2.cross(vel_sp2)

    angVel_clump = O.bodies[target_clump_ID].state.angVel
    mmoi_clump = O.bodies[target_clump_ID].state.inertia
    mmoi_mat_clump = Matrix3([
        mmoi_clump[0], 0, 0,
        0, mmoi_clump[1], 0,
        0, 0, mmoi_clump[2]
    ])
    L_clump = mmoi_mat_clump * angVel_clump

    L_tot = L_clump + L_sp1 + L_sp2
    angMom_clump = O.bodies[target_clump_ID].state.angMom
    print(L_tot)

    plot.addData(
        t = O.time, t1 = O.time,
        wx = w[0], wy = w[1], wz = w[2],
        Lx = L_tot[0], Ly = L_tot[1], Lz = L_tot[2],
        )

plot.plot()
# ----------------------------------------------------------------------------- Sim controller
O.dt = 2e-7#utils.PWaveTimeStep()
O.saveTmp()

Thanks.
Re-posted here: https://gitlab.com/yade-dev/trunk/-/issues/214

Bruno

In fact when it is a new bug better start a new thread (if on launchpad), or even better report an "issue" directly on gitlab:
https://gitlab.com/yade-dev/trunk/-/issues/new

B