Rotate cylinder

Asked by Mithushan Soundaranathan on 2021-03-25

Hi,

How can I rotate this half cylinder 90 degrees in y-axis:

#Cylinderpunch
cyl1 =geom.facetCylinder((0,-0.15,0),8*Diameter,0.07,Quaternion((0,0,1),.5*pi),24,wallMask=4)
cyl1 = [f for f in cyl1 if f.state.pos[0] < 0] # "half-cylinder"
cyl1=O.bodies.append(cyl1)

for j in cyl1:
     body= O.bodies[j]
     body.state.rot=(0,90,0)

I tried this, but it did not move my cylinder.

Here is my code:

#!/usr/bin/env python
#encoding: ascii

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test

from yade import utils, plot, timing
from yade import pack
import pandas as pd

o = Omega()

# Physical parameters
fr = 0.52
rho = 1050 #kg/m3
Diameter = 0.0079 #79 micrometer
r1 = Diameter
r2 = Diameter
k1 = 126000.0
kp = 12.0*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

o.dt = 1.0e-5

particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

#*************************************

# Add material
mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression
# Spheres for compression

sp=pack.SpherePack()
sp.makeCloud((-9.5*Diameter,-9.5*Diameter,-35.0*Diameter),(9.5*Diameter,9.5*Diameter,15.0*Diameter), rMean=Diameter/2.0)

#Cylinder cloud sphere pack
cyl = pack.inCylinder((0,0,0),(0,0,20*Diameter),8*Diameter-0.006)
sp = pack.filterSpherePack(cyl,sp,True,material=mat1)
sp.toSimulation(material=mat1)

##No of particles
count = 0
for b in O.bodies:
    if isinstance(b.shape, Sphere):
        count +=1
##Mass of particles
n_p=count

######################################################################
walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=8*Diameter,height=0.066,segmentsNumber=20,wallMask=6,material=mat1))
#height=0.067
##Top punch##
#find max Z coordinate of your cloud
max_z = aabbExtrema()[1][2]
min_z = aabbExtrema()[0][2]

#Top_punch=O.bodies.append(wall((0,0,max_z), axis=2, sense=-1, material=mat1))
#O.bodies[Top_punch].state.vel = (0, 0, -0.01)

#Cylinderpunch
cyl1 =geom.facetCylinder((0,-0.15,0),8*Diameter,0.07,Quaternion((0,0,1),.5*pi),24,wallMask=4)
cyl1 = [f for f in cyl1 if f.state.pos[0] < 0] # "half-cylinder"
cyl1=O.bodies.append(cyl1)

for j in cyl1:
 body= O.bodies[j]
 body.state.rot=(0,180,0)

for i in cyl1:
 body= O.bodies[i]
 body.state.vel = (0,1,0)

df = pd.DataFrame(columns=['time','Porosity','Maxoverlap','Stress'])
df.to_csv('PH101_Gao_speed_5.csv')
from csv import writer
def append_list_as_row(file_name, list_of_elem):
    # Open file in append mode
    with open(file_name, 'a+', newline='') as write_obj:
        # Create a writer object from csv module
        csv_writer = writer(write_obj)
        # Add contents of list as last row in the csv file
        csv_writer.writerow(list_of_elem)

# Add engines
o.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
                         Bo1_Wall_Aabb(),
                         Bo1_Facet_Aabb()
                         ]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
     Ig2_Facet_Sphere_ScGeom(),
     Ig2_Wall_Sphere_ScGeom()],
    [Ip2_LudingMat_LudingMat_LudingPhys()],
    [Law2_ScGeom_LudingPhys_Basic()]
  ),
  NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
  #VTKRecorder(fileName='vtk-',recorders=['all'],iterPeriod=10000),
  #PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
  #PyRunner(command='print(voxelPorosity(resolution=600,start=aabbExtrema()[0],end=aabbExtrema()[1]))', realPeriod=15),
  PyRunner(command='checkdata()', realPeriod=3),
  PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
  DeformControl(label="DefControl")
]

def checkdata():
 if O.iter < 20000:
  return
 else:
  time=o.realtime
  porosity=voxelPorosity(resolution=600,start=aabbExtrema()[0],end=aabbExtrema()[1])
  Maxoverlap=utils.maxOverlapRatio()
  stress=getStress().trace()/3.
  row_contents = [time,time,porosity,Maxoverlap,stress]
  append_list_as_row('PH101_Gao_speed_5.csv', row_contents)

def checkForce():
    # at the very start, unbalanced force can be low as there is only few
    # contacts, but it does not mean the packing is stable
    if O.iter < 20000:
        return
    # the rest will be run only if unbalanced is < .1 (stabilized packing)
    timing.reset()
    if unbalancedForce() > 0.2:
        return
    # add plate at upper box side

    highSphere = 0.0
    for b in O.bodies:
        if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
            highSphere = b.state.pos[2]
        else:
            pass

    O.bodies.append(wall(highSphere+0.5*Diameter, axis=2, sense=-1, material=mat1))
    # without this line, the plate variable would only exist inside this
    # function
    global plate
    plate = O.bodies[-1] # the last particles is the plate
    # Wall objects are "fixed" by default, i.e. not subject to forces
    # prescribing a velocity will therefore make it move at constant velocity
    # (downwards)
    plate.state.vel = (0, 0, -0.01)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2021-03-30
Last query:
2021-03-30
Last reply:
2021-03-26
Best Jan Stránský (honzik) said : #1

Hello,

> Here is my code:

most of the code is irrelevant to the problem.
Please try to create a MWE [1], M = minimal (unnecessary code removed, focusing on the problem)
Your "initial" code would be just fine (with "Diamater" variable defined), showing the entire problem ("I intended to rotate these bodies, but they are not rotated")
The rest just scares/discourages potential helpers :-)

state.rot is a method to get relative rotation from initial rotation [1], something like state.displ().
With your code, you just overwritten the method to the (0,90,0) tuple, try this:
###
Diameter = 0.0079 #79 micrometer
cyl1 =geom.facetCylinder((0,-0.15,0),8*Diameter,0.07,Quaternion((0,0,1),.5*pi),24,wallMask=4)
cyl1 = [f for f in cyl1 if f.state.pos[0] < 0] # "half-cylinder"
cyl1=O.bodies.append(cyl1)

for j in cyl1:
    body= O.bodies[j]
    print(body.state.rot) # <bound method State.rot ...
    body.state.rot=(0,90,0)
    print(body.state.rot) # (0,90,0)
###

Just a note, in Yade the angles are rarely (if at all) in degrees..

To rotate a body, you have to modify state.ori. If the center of rotation is not state.pos, you have to modify state.pos as well:
###
Diameter = 0.0079 #79 micrometer
cyl1 =geom.facetCylinder((0,-0.15,0),8*Diameter,0.07,Quaternion((0,0,1),.5*pi),24,wallMask=4)
cyl1 = [f for f in cyl1 if f.state.pos[0] < 0] # "half-cylinder"
cyl1=O.bodies.append(cyl1)

rot = Quaternion((0,0,1),.5*pi) # 90 degrees around z axis
center = Vector3(0,0,0) # center of rotation

for j in cyl1:
    body= O.bodies[j]
    # change rotation
    body.state.ori *= rot
    # shift the body
    p = body.state.pos - center # relative vector from center to pos
    p = rot * p # rotated relative vector
    body.state.pos = center + p # new pos
###

cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.State.rot

Jérôme Duriez (jduriez) said : #2

Depending on your needs (change the default initial orientation by a 90 degrees rotation, or impose a gradual rotation up to 90 degrees), it is a good YADE habit (for the 2nd case) to go through angular velocities, resp. velocities, and not update directly orientation, resp. position.
See https://yade-dem.org/doc/user.html#imposed-velocity (the same may apply for angular velocities and orientation)

Thanks Jan Stránský, that solved my question.