# yade-users team mailing list archive

## [Question #695558]: Dynamics of pfacet object does not appear to be correct

```New question #695558 on Yade:

Hello,

I am simulating the interaction between a brush (bristles made of grid nodes and grid connections) and an object (made of pfacet). The bristles are arranged to form a square at its  base. The objects takes the form of a wedge placed symmetrically between the bristles. The object is made to move into the brush

I expect the object to bounce back up, but it is starting to rotate when it contacts the brush. Since the system is symmetric, I do not expect this rotation to happen. I could not find the bug in my code. I suppose there could be something in the dynamics I am missing.

Kind regards,
Rohit John
The minimum working code
# ---------------------------------------------------------------------------------------------------------------------------------------

# user define function
'''
Adds the engines which sets the initial velocities of a target. This is done because for some reason we cannot graph the velocities for non spherical objects by manually setting the velocities
'''
linear_velocity_magnitude = L2_norm(linear_velocity)
linear_velocity_vector    = normalise(linear_velocity)

angular_velocity_magnitude = L2_norm(angular_velocity)
angular_velocity_vector    = normalise(angular_velocity)

O.engines += [
CombinedKinematicEngine(
ids   = ids,
label ='combined_motion_engine') +

TranslationEngine(
translationAxis = linear_velocity_vector,
velocity        = linear_velocity_magnitude,
label           = 'translation_engine') +

RotationEngine(
rotationAxis     = angular_velocity_vector,
angularVelocity  = angular_velocity_magnitude,
rotateAroundZero = True,
zeroPoint        = center_of_rotation,
label            = 'rotation_engine'
),

PyRunner(command = "motion_engine_stopper()", iterPeriod = 1, label = "stopper")

]

# ---------------------------------------------------------------------------------------------- motion_engine_stopper
def motion_engine_stopper():
stopper.iterPeriod      = int(1e14)
print("Engines stopped")

# ---------------------------------------------------------------------------------------------- L2_norm
def L2_norm(vector):
mag = 0
for i in vector:
mag += i**2

return mag**0.5

# ---------------------------------------------------------------------------------------------- normalise
def normalise(vector):
mag = L2_norm(vector)
if mag == 0:
return [0,0,0]
res = [i/mag for i in vector]
return res

# ---------------------------------------------------------------------- defining material
bristle_length = 0.1
brush_x_dim    = 0.1
brush_y_dim    = 0.02
brush_position = [0,0.0025,0]
brush_ori      = Quaternion((1,0,0), 0)

bristle_x_no   = 10
bristle_y_no   = 2
bristle_seg_no = 5
bristle_x_density = bristle_x_no / brush_x_dim
bristle_y_density = bristle_y_no / brush_y_dim

bristle_young    = 3e9
bristle_poisson  = 0.3
bristle_density  = 1000

covar = 0
bristle_tip_spread_covariance = [[covar, 0], [0, covar]]

target_dimension = [0.1,0.1,0.1]
target_position  = [0,0,0]
target_ori       = Quaternion((1,0,0), pi/4)
target_young    = 30e9
target_poisson  = 0.3
target_density  = 1000

# ---------------------------------------------------------------------- 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(),
Bo1_Sphere_Aabb(),
]),

InteractionLoop(
[
Ig2_PFacet_PFacet_ScGeom(),
Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
Ig2_GridNode_GridNode_GridNodeGeom6D(),
Ig2_GridConnection_PFacet_ScGeom(),
Ig2_Sphere_PFacet_ScGridCoGeom(),
Ig2_Sphere_Sphere_ScGeom()
],
[
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(),
],
),
NewtonIntegrator(gravity = [0.0,0.0,0.0], damping = 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,
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,
]

x_len = 0.03
y_len = 0.05
z_len = 0.05
nodes = [
]

target_pos = [0,0,0] # centroid

for j in range(3):
for i in nodes:
target_pos[j] += i[j]
target_pos[j] / 4

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

pnode  = []
pcyl   = []
pfacet = []
for i in pfacet_nodes:
print(i)
pfacetCreator1(
i,
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
# # ---------------------------------------------------------------------- motion engine
target_lin_vel = [0,0,-1]
target_ang_vel = [0,0,0]

ids = target_ids,
linear_velocity  = target_lin_vel,
angular_velocity = target_ang_vel,
center_of_rotation = target_pos
)

# ----------------------------------------------------------------------------- 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

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()

--