Question #695558 on Yade changed: https://answers.launchpad.net/yade/+question/695558
Description changed to: 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 from yade import plot, utils from yade.gridpfacet import * # user define function # ---------------------------------------------------------------------------------------------- add motion engines def add_motion_engines(ids, linear_velocity, angular_velocity, center_of_rotation): ''' 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(): combined_motion_engine.dead = True translation_engine.dead = True rotation_engine.dead = True 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_radius = 1e-3 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 bristle_friction = radians(30) covar = 0 bristle_tip_spread_covariance = [[covar, 0], [0, covar]] target_dimension = [0.1,0.1,0.1] target_position = [0,0,0] grid_radius = 1e-3 target_ori = Quaternion((1,0,0), pi/4) target_young = 30e9 target_poisson = 0.3 target_density = 1000 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(), 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, 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 + bristle_length + bristle_radius + grid_radius], [-x_len, 0, 0 + bristle_length + bristle_radius + grid_radius], [0.0, y_len, z_len + bristle_length + bristle_radius + grid_radius], [0.0,-y_len, z_len + bristle_length + bristle_radius + grid_radius] ] 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[2], nodes[1]], [nodes[0], nodes[1], nodes[3]] ] 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 # # ---------------------------------------------------------------------- motion engine target_lin_vel = [0,0,-1] target_ang_vel = [0,0,0] add_motion_engines( 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 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() -- 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