testDeformableBodies.py 4.23 KB
Newer Older
Burak Er's avatar
Burak Er committed
1 2 3 4 5 6 7 8 9
#!/usr/bin/python
# -*- coding: utf-8 -*-
# Created by Burak ER
from yade.deformableelementsutils import *;

O=Omega()

## Create deformable elements internal materials and interaction element material

10 11 12 13 14 15
try:
	mat=LinIsoRayleighDampElastMat(label='aluminiummaterial');
except NameError:
	print("\nDid you compile yade with cmake option -DENABLE_FEMLIKE=1 ?\n")
	import sys
	sys.exit()
Burak Er's avatar
Burak Er committed
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

intermat=LinCohesiveStiffPropDampElastMat(label='intermat')

intermat.youngmodulus=70e11;

mat.density=2700;

mat.youngmodulus=70e9

O.materials.append(intermat)

O.materials.append(mat)

## Generate deformable element mesh from GMSH2 format

mesh=tetrahedral_mesh_generator('model.msh',Lin4NodeTetra,'aluminiummaterial',Lin4NodeTetra_Lin4NodeTetra_InteractionElement,'intermat')

## Define the body boundary force

# x position limits of bodies that are subject to the force
force_tail_body_lowerlimit=-1.01;

force_tail_body_upperlimit=-0.9999;

# Angular frequency of the force
amplitude=1000;

period=(1e-4)

45
omega=2*pi/period;
Burak Er's avatar
Burak Er committed
46

47
applicationperiod=period/2;
Burak Er's avatar
Burak Er committed
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105

## Define the fixed boundary

# x position limits of bodies that are subject to the fixed boundary

fixed_tail_body_lowerlimit=0.99;

fixed_tail_body_upperlimit=1.01;


def getinitialpos():
	positions=[];
	for body in O.bodies:
		positions.append(body.state.pos)
	return positions

def getboundarybodies():
	bdy=[];
	for body in O.bodies:
			if(body.shape.dispIndex==12):
				if(body.state.pos[0]<fixed_tail_body_upperlimit):
					if(body.state.pos[0]>fixed_tail_body_lowerlimit):
						bdy.append(body.id)
	return bdy
						
def getforcebodies():
	bdy=[];
	for body in O.bodies:
			if(body.shape.dispIndex==12):
				if(body.state.pos[0]<force_tail_body_upperlimit):
					if(body.state.pos[0]>force_tail_body_lowerlimit):
						if(O.time<applicationperiod):
							bdy.append(body.id)	
	return bdy						

initialpositions=getinitialpos();

forcebodies=getforcebodies();

boundarybodies=getboundarybodies();


## definition of functions
def fixboundaryelements():
	for i in boundarybodies:
		O.forces.addF(i,-O.forces.f(i))
		O.forces.addT(i,-O.forces.m(i))

def applyforcetoelements():
	for i in forcebodies:
		if(O.time<applicationperiod):
#			print 'apply force @'+str(O.time)
			O.forces.addF(i,-1*O.forces.f(i));
	
			O.forces.addF(i,Vector3(sin(omega*O.time)*amplitude,0,0));
#			plot.addData(force=O.forces.f(i));

def addplot():
106
	plot.addData(force=O.forces.f(forcebodies[1])[0],pos=(O.bodies[forcebodies[1]].state.pos[0]-initialpositions[forcebodies[1]][0]),vel=O.bodies[forcebodies[1]].state.vel[0],t=O.time,time=O.time,tm=O.time)
Burak Er's avatar
Burak Er committed
107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
		
for i in forcebodies:
	O.bodies[i].state.pos+=Vector3(0,0,0.05);

## Adaptive integration engine initialization: we have given the engines in ordered, they are runned sequentally for every substep of the integrator

#integratoreng=RungeKuttaCashKarp54Integrator([	## Resets forces and momenta the act on bodies
#						ForceResetter(),	
#						## Apply internal force to the deformable elements and internal force of the interaction element	 
#						FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
#						PyRunner(virtPeriod=1e-99,command='applyforcetoelements()'),
#
#])
# Integration tolerances of the RungeKuttaCashKarp54Integrator
#~ 
#~ integratoreng.rel_err=1e-6;
#~ 
#~ integratoreng.abs_err=1e-6;

# We use only the integrator engine

# Time step determines the exiting period of the integrator since the integrator performs one step from current_time to current_time+dt; using many substeps for any value of dt; then stops. 

O.dt=1e-8;

O.engines=[
133 134 135 136 137
	ForceResetter(),	
	## Apply internal force to the deformable elements and internal force of the interaction element	 
	FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
	PyRunner(iterPeriod=1,command='applyforcetoelements()'),
	NewtonIntegrator(damping=0,gravity=[0,0,0]),
Burak Er's avatar
Burak Er committed
138
#	    ## Plotting data: adds plots after one step of the integrator engine
139
	PyRunner(iterPeriod=1,command='addplot()')
Burak Er's avatar
Burak Er committed
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156


];
#Tolerances can be set for the optimum accuracy


from yade import plot

plot.plots={'t':'vel','time':'pos','tm':'force'}
plot.plot(subPlots=True)

try:
	from yade import qt
	qt.View()
	qt.Controller()
except ImportError: pass