| 33 | import rod |
| 34 | |
| 35 | class FE_model: |
| 36 | def ascending(self, i,j): # {{{ |
| 37 | return (i,j) if i < j else (j,i) |
| 38 | # }}} |
| 39 | def connectivity(self, triangle_mesh): # {{{ |
| 40 | """ |
| 41 | Given a mesh from triangle, returns nodes, (mesh's ['vertices']) |
| 42 | and conn, a set of elements defined by their node ID's, eg, |
| 43 | nodes = [[ 0.0, 0.0], |
| 44 | [10.0, 0.0], |
| 45 | [10.0, 1.0], ... ] |
| 46 | conn = {(0, 3), |
| 47 | (0, 18), |
| 48 | (1, 2), |
| 49 | (1, 30), ... } |
| 50 | Node ID's begin with 0. |
| 51 | """ |
| 52 | conn = [] |
| 53 | for i,j,k in triangle_mesh['triangles']: |
| 54 | conn.append( self.ascending(i, j) ) # add node ID tuples in |
| 55 | conn.append( self.ascending(j, k) ) # ascending order so that |
| 56 | conn.append( self.ascending(i, k) ) # set() removes duplicates |
| 57 | conn = set(conn) # remove duplicates |
| 58 | return triangle_mesh['vertices'], conn |
| 59 | # }}} |
| 60 | def __init__(self, triangle_mesh): # {{{ |
| 61 | self.E = 10.0e6 # psi |
| 62 | self.R = 0.1 # radius, inches |
| 63 | self.rho = 0.1 # density, lb/inches**3 |
| 64 | nodes, conn = self.connectivity(triangle_mesh) |
| 65 | |
| 66 | self.node = [] |
| 67 | for nid,(x,y) in enumerate(nodes): |
| 68 | self.node.append( rod.Node(id=nid, x=x, y=y) ) |
| 69 | # store coordinate min/max to set plot frame, |
| 70 | # help with autoscaling |
| 71 | self.xmin = min([_.x for _ in self.node]) |
| 72 | self.xmax = max([_.x for _ in self.node]) |
| 73 | self.ymin = min([_.y for _ in self.node]) |
| 74 | self.ymax = max([_.y for _ in self.node]) |
| 75 | |
| 76 | self.constrained_dof = [] |
| 77 | for nid in triangle_mesh['constrained_x']: |
| 78 | self.constrained_dof.append( 2*nid ) |
| 79 | for nid in triangle_mesh['constrained_y']: |
| 80 | self.constrained_dof.append( 2*nid + 1 ) |
| 81 | self.constrained_dof = set( self.constrained_dof ) |
| 82 | #print('INIT: constrained_dof=',self.constrained_dof) |
| 83 | |
| 84 | _elem = {} |
| 85 | for eid,(i,j) in enumerate(conn): |
| 86 | nA = rod.Node(id=i, x=nodes[i,0], y=nodes[i,1]) |
| 87 | nB = rod.Node(id=j, x=nodes[j,0], y=nodes[j,1]) |
| 88 | _elem[eid] = rod.Rod_Elem(eid, self.node[i], self.node[j], |
| 89 | self.R, self.E, self.rho) |
| 90 | self.elem = _elem |
| 91 | self.nDof = 2*len(self.node) |
| 92 | # }}} |