yade-users team mailing list archive

[Loic Dugelas <question632960@xxxxxxxxxxxxxxxxxxxxx>]

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

Loic Dugelas posted a new comment:
Here you can find my script

#-*- coding: utf-8 -*-
from yade import qt, plot
from yade.gridpfacet import *
import gts, os, locale, sys
import numpy as np


################ Engine ############
	
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_GridConnection_Aabb(),]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_GridNode_GridNode_GridNodeGeom6D(),Ig2_Sphere_GridConnection_ScGridCoGeom(),Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),],
		[Ip2_WireMat_WireMat_WirePhys(),Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),],
		[ Law2_ScGeom_WirePhys_WirePM(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGridCoGeom_CohFrictPhys_CundallStrack(),Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),Law2_ScGridCoGeom_FrictPhys_CundallStrack(),Law2_ScGeom_FrictPhys_CundallStrack(),]
	),
	NewtonIntegrator(gravity=(0,0,-9.81),damping=0.9,label='newton'),
	PyRunner(initRun=True,iterPeriod=200,command='Force()'),
	PyRunner(virtPeriod=0.05, virtLast = 3, command='MEP_Pret()'),
]

################## Material & geometry ##########
youngc = 60e9
strenruptc=1400e9
densityc=7850.
poissonc = 0.3
shearCohc = 1e100
Rc = 0.008
cablemat = O.materials.append(CohFrictMat(young=youngc,poisson=poissonc,density=densityc,frictionAngle=radians(10),normalCohesion=pi*strenruptc,isCohesive=True,shearCohesion=shearCohc,momentRotationLaw=False, alphaKr = 1.0))
cablecomat = O.materials.append(FrictMat(young=youngc,poisson=poissonc,density=densityc,frictionAngle=radians(10)))

####################### Adding cable ##############
def Cable_droit_spheres(start, end, lseg, mtxnode, rayon): #just the meshing function

	lenght = end-start
	vecunit = lenght/(abs(lenght))
	lfinal = abs(lenght)
	cable = []
	nbseg = int(lfinal/lseg)
	for i in range(nbseg+1):
		x = start[0] + (vecunit*lfinal)[0]*i/nbseg
		y = start[1] + (vecunit*lfinal)[1]*i/nbseg
		z = start[2] + (vecunit*lfinal)[2]*i/nbseg
		cable.append(O.bodies.append(gridNode([x,y,z],rayon,wire=False,fixed=False,material=mtxnode)))
	
	return [cable]

#mapping the spheres of the cable
[c2a_1] = Cable_droit_spheres(start = Vector3(-1,0,0) , end = Vector3(-0.25, 0, 0), lseg = 0.1, mtxnode = cablemat,rayon = Rc)
[c2a_2] = Cable_droit_spheres(start = Vector3(-.20,0,-0.1), end = Vector3(0.20,0,-0.1), lseg = 0.1,mtxnode = cablemat,rayon = Rc)
[c2a_3] = Cable_droit_spheres(start = Vector3(0.25,0,0), end = Vector3(1,0,0), lseg = 0.1,mtxnode = cablemat,rayon = Rc)

cable = c2a_1 + c2a_2 + c2a_3


#and the link
concable = []
for i,j in zip(cable[:-1],cable[1:]):
	concable.append(O.bodies.append(gridConnection(O.bodies[i].id,O.bodies[j].id,Rc,wire=False, material = cablecomat) ))

O.bodies[cable[0]].state.blockedDOFs='xyzXYZ'
O.bodies[cable[-1]].state.blockedDOFs='xyzXYZ'

plot.plots={'t':('Fcable'),}
plot.plot(noShow=False, subPlots=False)
O.dt = 1e-6
#O.run()
################### Functions ####################

def MEP_Pret():
	
	Lcable = 0
	for i in range(len(cable)-1):
		Lcable += np.linalg.norm(O.bodies[cable[i]].state.pos-O.bodies[cable[i+1]].state.pos)

	if O.time > 0.6 :
		O.pause()
		pos1 = O.bodies[cable[-1]].state.pos

		for i in O.bodies[cable[-1]].intrs(): #deleting interaction with the extremity of the cable
			O.interactions.erase(i.id1,i.id2)
		for i in O.bodies[concable[-1]].intrs():
			O.interactions.erase(i.id1,i.id2)

                O.bodies.erase(concable[-1]) #deleting the cylinder (the
crash comes from here)

		O.bodies[cable[-1]].state.pos = Vector3(pos1[0]-0.10,pos1[1],pos1[2])
		concable.append(O.bodies.append( gridConnection(cable[-2],cable[-1], Rc,material = cablecomat) ))
		O.bodies[cable[-1]].state.pos = pos1
		O.run()
		
def Force():

	Fcable = np.linalg.norm(O.forces.f(cable[0]))
	t = O.time
	plot.addData(
	Fcable = Fcable/1000,
	t = O.time)
	if t> 3.0 :
		O.pause()

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