#(c) Rory Clune 05/16/2011 from numpy import * from numpy.linalg import inv import pygame class structure: def __init__(self, geom = None, connect = None, fix = array([0,1,2]), load = None, A = None, I = None): global connectivity, ndof, nbeams, fixity self.geometry = geom self.connectivity = connect self.fixity = fix self.F = load if A == None: self.A = 10*ones((connect.shape[0])) else: self.A = A if I == None: self.I = 10*ones((connect.shape[0])) else: self.I = I self.displacement_limit = 0.00 def analyze(self): geometry = self.geometry; connectivity = self.connectivity; fixity = self.fixity; F = self.F #-------- # 'MESHING' - TO BE ADDED - need to generate new connectivity and geometry for higher resolution #-------- # nbeams = connectivity.shape[0] # div = 3 # conn = zeros((nbeams * div, 2)) # geom = zeros((geometry.shape[0] + nbeams*(div-1), 2)) # for i, e in enumerate(connectivity): # #two entries are known # conn[e[0]*div,0] = e[0]*div # conn[e[0]*div+div-1,1] = e[1]*div # geom[e[0]*div,:] = geometry[e[0],:] # geom[e[1]*div,:] = geometry[e[1],:] # #the rest (new nodes) need to be filled in, using various types of interpolation # conn[e[0]*div,1] = conn[e[0]*div,0]+1 # conn[e[0]*div+div-1,0] = conn[e[0]*div+div-1,1]-1 # for j in range(1,div-1,1): # conn[e[0]*div+j,0] = conn[e[0]*div+j-1,1] # conn[e[0]*div+j,1] = conn[e[0]*div+j,0] + 1 # for j in range(1,div,1): # geom[e[0]*div+j,:] = geometry[e[0],:] + (float(j)/div)*(geometry[e[1],:] - geometry[e[0],:]) # ndof = 3*geometry.shape[0] # newF = zeros((ndof)) # for i in range(0, self.F.shape[0], 1): # newF[div*(i/div)+i] = self.F[i] # self.F = newF # newFixity = zeros((ndof)) # for i in range(0,self.fixity.shape[0],1): # newFixity[div(i/div)+i] = self.fixity[i] # self.fixity = newFixity # raw_input('paused after subdivision') #-------- # BUILD EFT from connectivity matrix #-------- EFT = zeros((connectivity.shape[0],6), dtype = int) for i in range(0,connectivity.shape[0],1): EFT[i,0:3] = 3*(connectivity[i,0]+1)-arange(3,0,-1) EFT[i,3:6] = 3*(connectivity[i,1]+1)-arange(3,0,-1) #-------- # BUILD STIFFNESS MATRIX from EFT & element info #-------- K = zeros((ndof, ndof)) for i, e in enumerate(connectivity): # build local (element) stiffness matrix: l = sum((geometry[e][0]-geometry[e][1])**2)**0.5 c = (geometry[e][1][0] - geometry[e][0][0])/l s = (geometry[e][1][1] - geometry[e][0][1])/l e = 500; a = self.A[i]; I = self.I[i] k_local = array( [ [e*a*c**2/l + 12*e*I*s**2/l**3, e*a*c*s/l - 12*e*I*c*s/l**3, -6*e*I*s/l**2, -e*a*c**2/l - 12*e*I*s**2/l**3, -e*a*c*s/l + 12*e*I*c*s/l**3, -6*e*I*s/ l**2], [e*a*c*s/l - 12*e*I*c*s/l**3, e*a*s**2/l + 12*e*I*c**2/l**3, 6*e*I*c/l**2, -e*a*c*s/l + 12*e*I*c*s/l**3, -e*a*s**2/l - 12*e*I*c**2/l**3, 6*e*I*c/l**2], [6*e*I*s/l**2, 6*e*I*c/l**2, 4*e*I/l, 6*e*I*s/l**2, -6*e*I*c/l**2, 2*e*I/l], [-e*a*c**2/l - 12*e*I*s**2/l**3, -e*a*c*s/l + 12*e*I*c*s/l**3, 6*e*I*s/l**2, e*a*c**2/l + 12*e*I*s**2/l**3, e*a*c*s/l - 12*e*I*c*s/l**3, 6*e*I*s/l**2], [-e*a*c*s/l + 12*e*I*c*s/l**3, -e*a*s*s/l - 12*e*I*c**2/l**3, -6*e*I*c/l**2, e*a*c*s/l - 12*e*I*c*s/l**3, e*a*s**2/l + 12*e*I*c**2/l**3, -6*e*I*c/l**2], [-6*e*I*s/l**2, 6*e*I*c/l**2, 2*e*I/l, 6*e*I*s/l**2, -6*e*I*c/l**2, 4*e*I/l] ]) # assemble global stiffness matrix: K[EFT[i,0]:EFT[i,2]+1, EFT[i,0]:EFT[i,2]+1 ] += k_local[0:3,0:3] #DO THIS IN FEWER STEPS? using enumeration? K[EFT[i,3]:EFT[i,5]+1, EFT[i,3]:EFT[i,5]+1 ] += k_local[3:6,3:6] K[EFT[i,0]:EFT[i,2]+1, EFT[i,3]:EFT[i,5]+1 ] += k_local[0:3,3:6] K[EFT[i,3]:EFT[i,5]+1, EFT[i,0]:EFT[i,2]+1 ] += k_local[3:6,0:3] #-------- # APPLY BOUNDARY CONDITIONS #-------- K[fixity[:],:] = 0; K[:,fixity[:]] = 0 K[fixity[:],fixity[:]] = 1 F[fixity[:]] = 0 self.K = K u = dot(linalg.inv(K), F) return u def EnergyObjective(self, x): self.input_x(self, x) u = self.analyze() E = dot(dot(transpose(u), self.K), u) # E = dot(dot(transpose(self.F), transpose(linalg.inv(self.K))), self.F) return E def VolumeObjective(self, x): self.input_x(self, x) C = 0 for i, e in enumerate(self.connectivity): l = sum((self.geometry[e][0]-self.geometry[e][1])**2)**0.5 C += self.A[i]*l return C def Constraint(self, x): self.input_x(self, x) ans = maximum(0, fabs(self.analyze())-self.displacement_limit) # vector of displacements return ans #----------- # DISPLAY A STRUCTURE USING PYGAME #----------- def show(self, xscale = 40, yscale = 40, xshift = 20, yshift = 100, load_pos = 1): '''Animates the loading of a structure using pygame''' screen = pygame.display.set_mode((640, 480)) k = 0; maxk = fabs(self.F[load_pos[0]]); dk = int(maxk/10); running = 1; cycles= -1; count = 0 while running: count+=1 # Pygame stuff: pygame.time.wait(10) event = pygame.event.poll() if event.type == pygame.QUIT: running = 0 screen.fill((255, 255, 255)) # Set the load: load = k self.F[load_pos] = k # Find the displacements u = self.analyze() # Draw the structure: for i, pt in enumerate(self.connectivity): x1 = xscale * (self.geometry[pt[0],0] + u[3*(pt[0]+1)-3]) + xshift y1 = yscale*(self.geometry[pt[0],1] + u[3*(pt[0]+1)-2])+yshift x2 = xscale*(self.geometry[pt[1],0] + u[3*(pt[1]+1)-3])+xshift y2 = yscale*(self.geometry[pt[1],1] + u[3*(pt[1]+1)-2])+yshift pygame.draw.line(screen, (0,0,255), (x1, y1), (x2, y2), 2) pygame.draw.circle(screen, (255,140,0), (int(x1),int(y1)), 4, 0) pygame.draw.circle(screen, (255,140,0), (int(x2),int(y2)), 4, 0) if k >= maxk or k <= -maxk: dk *= -1 cycles += 1 k += dk #Output image at each frame for movie making: # pygame.image.save(screen, '%d.png' % count) pygame.display.flip() if(cycles == -1): pygame.time.wait(1500); cycles += 1 pygame.image.save(screen, 'Structure.png') if cycles>15: cycles = -1 running = 0