#!/usr/bin/env python ''' fem.py - a program that handles truss members (t) with frame members, i.e. beams (b). ''' # # (c) Copyright 2006, 2008 Algis Kabaila. All rights reserved. # # This document is free; you can redistribute it and/or modify it under the terms # of the GNU General Public License as published by the Free Software Foundation; # either version 2 of the License, or (at your option) any later version. # # This document is distributed in the hope that it will be useful, but WITHOUT # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public # License for more details. # # You should have received a copy of the GNU General Public License along with # this document; if not, write to the Free Software Foundation, Inc., # 675 Mass Ave, Cambridge, MA 02139, USA. # # Version 0.4.4-x x = revision number = $Rev$ # # $Id$ # Trying vcs of bazaar with small changes. # import sys from math import sqrt from core_mods import matrixM from core_mods import core_mod as cm_ from core_mods import joint_mod as jm_ # select any available matrix-linear algebra module m_ = matrixM.Matrix() # Global functions ------------------------------------------------------- def d_print(line): '''Print string line on the console if debug is True.''' if debug: print '*d* ', line # Definition of object classes ------------------------------------------- class Prescribes: "Prescribed non-zero displacements for frames." def __init__(self, glob, nodeList): self.glob = glob self.nodeList = nodeList self.prescribesDict = {} # Format prescribesDict = {freedomNo(world) : [vals LC1, vals LC2..]) # vals are the magnitude of the prescribe - # all in world numbering 1, 2 , 3 etc, def modify(self, Ks, Ps, nodeList): ''' Modify Ks and Ps for effects of prescribed displacements. Called by getKsPsUs. ''' glob = self.glob if glob.noPrescribes == 0: return Ks, Ps # Quit without modifying Ks and Ps f = glob.f fo = glob.fo nn = glob.noFreedoms bigNo = 0 # pick the larges diagonal term in the stiffness matrix for i in range(nn): if Ks[i,i] > bigNo: bigNo = Ks[i,i] bigNo = bigNo * 1e20 prescrdict = self.prescribesDict # alias freedomvec = prescrdict.keys() for i in range(len(freedomvec)): freedom = self.prescribedFreedoms[i] - 1 Ks[freedom, freedom] = bigNo for ii in range(glob.noLoadCases): try: value = prescrdict[freedom + 1][ii] except KeyError: cm_.printer(fo, 'no such key') raise KeyError Ps[freedom, ii] = value * bigNo return Ks, Ps def readPrescribes(self): ''' Create dictionary "prescribeDict" with key = freedom and value = list of values, noLoads (load cases) long. Called by main(). ''' glob = self.glob f = glob.f fo = glob.fo if self.glob.noPrescribes == 0: cm_.getTitle(f, fo, 'end') return else: d = self.prescribesDict prescribedFreedoms = d.keys() # self.prescribedFreedoms title = cm_.getTitle(f, fo, 'prescribes') for ii in range(glob.noLoadCases): for i in range(glob.noPrescribes): lst = cm_.listRead(f,fo) if lst[0] == 'b': line = cm_.listProcess(fo, lst, 'sisssfff')[1:] # b 3 n y n 0 -10e-3 0 node = self.nodeList.ds[line[0]] # directory format: # {nodeNo : [bt, ms, x, y, # flagX='f', flagY='f', flagXY='f']} freedoms = [node[4], node[5], node[6]] # world units useFlags = [line[1], line[2], line[3]] values = [line[4], line[5], line[6]] elif lst[0] == 't': line = cm_.listProcess(fo, lst, 'sissff')[1:] node = self.nodeList.ds[line[0]] # directory format: # {nodeNo : [bt, ms, x, y, flagX='f', flagY='f']} freedoms = [node[4], node[5]] # world units useFlags = [line[1], line[2]] values = [line[3], line[4]] else: cm_.printer(fo, lst.__str__()) cm_.printer(fo, 'Found type = ' + lst[0]) cm_.abandon('Type must be either (b) or (t)') if ii == 0: for j in range(len(useFlags)): if useFlags[j] == 'y': d[freedoms[j]] = [values[j]] prescribedFreedoms.append(freedoms[j]) elif ii > 0: for j in range(len(useFlags)): if useFlags[j] == 'y': d[freedoms[j]].append(values[j]) self.prescribedFreedoms = prescribedFreedoms # world values self.prescribesDict = d # world values title = cm_.getTitle(f, fo, 'end') # End of Definition of object classes ------------------------------------- # start of frame functions--------------------------------------------- def get_beam_Bkbar(E, dx, dy, area, Ix): ''' Get frame member kbar stiffness and B matrix and member length L. doLocationVectors must be called **before** this method ''' B = m_.zerosM(6,3) # B for frames is (6x3) L = sqrt(dx * dx + dy * dy) a = dx/L # cos alpha b = dy/L # cos beta B[0,0] = -a B[0,1] = b B[1,0] = -b B[1,1] = -a B[2,1] = -L B[2,2] = -1 # ---------- B[3,0] = a B[3,1] = -b B[4,0] = b B[4,1] = a B[5,2] = 1 # Aliases: A = area # area kbar = m_.zerosM(3, 3) kbar[0,0] = A / Ix kbar[1,1] = 12./(L * L) kbar[1,2] = -6. / L kbar[2,1] = -6. / L kbar[2,2] = 4 kbar = (E * Ix / L) * kbar return L, B, kbar def getk(L, B, kbar): Bt = m_.transposeM(B) temp= m_.multM(B,kbar) return m_.multM(temp,Bt) # k for a frame member def frameStressMatrix(L, B, kbar): 'Find Stress Matrix for a beam (frame member).' A = m_.zerosM(4, 3) A[0,0] = A[1,1] = A[2,2] = A[3,2] = 1 A[2,1] = L Bt = m_.transposeM(B) t = m_.multM(kbar, Bt) stressMatrix = m_.multM(A, t) return stressMatrix # end of frame functions--------------------------------------------- # start of truss functions ------------------------------------------ def B_kbar(dx, dy, E, area): "Compute matrix B and member stiffness kbar for truss bar." B = m_.zerosM(4, 1) l = sqrt(dx * dx + dy * dy) cc = dx/l # cos ss = dy/l # sin B[0, 0] = -cc B[1, 0] = -ss B[2, 0] = cc B[3, 0] = ss kbar = (E * area) / l return B, kbar def k_truss_element(dx, dy, E, area): 'Compute k for truss bar.' B, kbar = B_kbar(dx, dy, E, area) Bt = m_.transposeM(B) k = kbar * m_.multM(B, Bt) return k def trussStressMatrix(B, kbar): '''Returns truss truss bar's stress matrix.''' Bt = m_.transposeM(B) stressMatrix = kbar * Bt return stressMatrix # end of truss functions -------------------------------------------- # start of 3d truss functions --------------------------------------- def B_kbar3d(dx, dy, dz, E, area): '''determine l, B, kbar for a 3d truss member.''' B = m_.zerosM(6, 1) l = sqrt(dx * dx + dy * dy + dz * dz) cc = dx/l # cos(x) = dir cos with x axis ss = dy/l # cos with y axis dd = dz/l # cos with z axis B[0, 0] = -cc B[1, 0] = -ss B[2, 0] = -dd B[3, 0] = cc B[4, 0] = ss B[5, 0] = dd kbar = (E * area) / l return B, kbar def k_truss3d(dx, dy, dz, E, area): '''Find k for a 3d truss member.''' B, kbar = B_kbar3d(dx, dy, dz, E, area) # print '*d* E, area = ' , E, area # print '*d* kbar = ', kbar Bt = m_.transposeM(B) # print '*d* B = ', B k = kbar * m_.multM(B, Bt) # print '*d* k = ', k return k def truss3dStressMatrix(B, kbar): '''Returns truss truss bar's stress matrix.''' Bt = m_.transposeM(B) stressMatrix = kbar * Bt return stressMatrix # end of 3d truss functions ----------------------------------------- class StructureStiffness: def __init__(self, glob, memberList, materialList, sectionList, loadMatrix): self.glob = glob self.memberList = memberList self.materialList = materialList self.sectionList = sectionList self.loadMatrix = loadMatrix def assembleStiffness(self): glob = self.glob noFreedoms = glob.noFreedoms noLoadCases = glob.noLoadCases memberList = self.memberList materialList = self.materialList sectionList = self.sectionList fo = glob.fo Ks = m_.zerosM(noFreedoms, noFreedoms) # Calculate member stiffnesses for all membeers and assemble to Ks. i_list = memberList.ds.keys() for i in i_list: member = memberList.ds[i] d_print(('member')) d_print(member) tag = member[0] if not(tag in ['t', 'b', 'u']): print 'unrecognized member tag in assembleStiffness method of structureStiffness class' sys.exit(['tag not recognized']) d_print(tag) E = materialList.dictE[member[-2]] d_print('E = %12.5E' % (E)) section = sectionList.ds[member[-1]] d_print(('section = ', section)) area = section[0] Ix = section[1] d_print('area, Ix = %12.5E %12.5E' % (area, Ix)) properties = memberList.properties[i] d_print(('properties = ', properties)) dx = properties[0] dy = properties[1] if tag == 'u': dz = properties[2] d_print(('tag,dx,dy,dz', tag, dx, dy, dz)) loc_vector = memberList.properties[i][3] else: d_print((' tag, dx, dy =', tag, dx, dy)) loc_vector = memberList.properties[i][2] d_print(('Location vector = ', loc_vector)) if tag == 't': k = k_truss_element(dx, dy, E, area) elif tag == 'b': l, B, kbar = get_beam_Bkbar(E, dx, dy, area, Ix) k = getk(l, B, kbar) elif tag == 'u': k = k_truss3d(dx, dy, dz, E, area) Ks = cm_.assemble(Ks, k, loc_vector) cm_.printer(fo, '\nStructure Stiffness Matrix, Ks.') m_.printer_m(fo, Ks) self.Ks = Ks return Ks def solve(self): glob = self.glob fo = glob.fo loadMatrix = self.loadMatrix memberList = self.memberList materialList = self.materialList sectionList = self.sectionList Ks = self.Ks us = m_.solveM(Ks, loadMatrix) cm_.printer(fo, 'Structure displacements (transposed matrix).') uss = m_.transposeM(us) m_.printer_m(fo, uss) i_list = memberList.ds.keys() cm_.printer(fo, " Member stress resultants.") cm_.printer(fo, " Member LC Axial Shear B.M. 1 B.M. 2") i_list.sort() for i in i_list: member = memberList.ds[i] tag = member[0] assert tag in ['t', 'b', 'u'], 'Unknown element tag in solve subroutine' E = materialList.dictE[member[-2]] section = sectionList.ds[member[-1]] area = section[0] Ix = section[1] properties = memberList.properties[i] dx = properties[0] dy = properties[1] if tag in ['t', 'b']: # loc_vector = memberList.properties[i][2] loc_vector = properties[2] elif tag == 'u': # loc_vector = memberList.properties[i][3] loc_vector = properties[3] dz = properties[2] # get the element displacements ue = findElDisps(glob.noLoadCases, loc_vector, us) if tag == 't': B, kbar = B_kbar(dx, dy, E, area) stressMatrix = trussStressMatrix(B, kbar) stress = m_.multM(stressMatrix, ue) ss = stress for ii in range(len(stress[0])): cm_.printer(fo, "%5d %5d %12.5E " %\ (i, ii+1, ss[0,ii])) cm_.printer(fo, ' ') elif tag == 'b': l, B, kbar = get_beam_Bkbar(E, dx, dy, area, Ix) stressMatrix = frameStressMatrix(l, B, kbar) stress = m_.multM(stressMatrix, ue) ss = stress for ii in range(len(stress[0])): cm_.printer(fo, "%5d %5d %12.5E %12.5E %12.5E %12.5E" %\ (i, ii+1, ss[0,ii], ss[1,ii], ss[2,ii], ss[3,ii])) cm_.printer(fo, ' ') elif tag == 'u': B, kbar = B_kbar3d(dx, dy, dz, E, area) stressMatrix = truss3dStressMatrix(B, kbar) stress = m_.multM(stressMatrix, ue) ss = stress for ii in range(len(stress[0])): cm_.printer(fo, "%5d %5d %12.5E " %\ (i, ii+1, ss[0,ii])) cm_.printer(fo, ' ') return us # ------------------------------------------------------------------- def findElDisps(noLoadCases, loc_vector, us): '''Extract element displacements, ue, from structure displacements us. ''' # print '*d* loc_vector = ', loc_vector size = len(loc_vector) d_print('\n*d* loc_vector = ' + loc_vector.__str__() ) ue = m_.zerosM(size, noLoadCases) row = -1 for v in loc_vector: i = abs(v) - 1 # Convert from world to program units. row += 1 for j in range(noLoadCases): if i < 0: ue[row, j] = 0 else: ue[row,j] = cm_.sgn(v) * us[i , j] return ue # ------------------------------------------------------------------- def main(): '''This is the mainline for combination of trusses with frames.''' f, fo = cm_.prolog('beam1a') glob = jm_.GlobalVariables(f, fo) print 'glob, global variable.' glob.readGlobals() print glob materialList = jm_.MaterialList(glob) materialList.readMaterialList() sectionList = jm_.SectionList(glob) sectionList.readSectionList() nodeList = jm_.NodeList(glob) nodeList.readNodeList() nodeList.findNoFreedoms() print glob memberList = jm_.MemberList(glob, nodeList) memberList.readMemberList() memberList.doLocationVectors() loadCasesList = jm_.LoadCasesList(glob, nodeList) loadCasesList.readLoadCase() loadMatrix = loadCasesList.doLoadMatrix() prescribes = Prescribes(glob, nodeList) prescribes.readPrescribes() cm_.printer(fo, '\nAll data read in. Solution follows. ------------------->\n') cm_.printer(fo, 'Load Matrix (transposed):') m_.printer_m(fo, m_.transposeM(loadMatrix)) print materialList print sectionList nodeList.printNodeList() memberList.printMemberList() memberList.print_member_properties() d_print('\nLoad Matrix') d_print(loadMatrix) structureStiffness = StructureStiffness(glob, memberList, materialList,\ sectionList, loadMatrix) Ks = structureStiffness.assembleStiffness() prescribes.modify(Ks, loadMatrix, nodeList) us = structureStiffness.solve() f.close() fo.close() # ------------------------------------------------------------------- if __name__ == '__main__': debug = False # True main()