PPPK_program_DiM_2021_(v0.2.1)

###########################
## Nodes, elements etc.:
###########################
      
        
def element_node_coords (el, nds) : 
    return [nds[el[0]], nds[el[1]]]


def coords_left_to_right (nodes) :
    return map (apply_to_ith (lambda x : -x, -1), nodes)
        
            
def ith_element_characteristics (ei, echs, chars) : 
    return chars[echs[ei]]
    
      
def A_rect (b, h) :
    return b * h
    
def A_circ (d) :
    pi_n = parent(d) (pi)
    return pi_n * d^2 / 4
    
    
def support_list (d) : 
    sl = []
    for s in d.keys() :
        sl.append (s)
    return sl
    
def support2d_dict (l) :
    sd = {} 
    for s in l : 
        sd[s] = (0, 0)
    return sd
    
def support3d_dict (l) :
    sd = {} 
    for s in l : 
        sd[s] = (0, 0, 0)
    return sd
    

 
###############################    
## Element stiffness matrix:
###############################
  
  
def truss_element_stiffness_matrix (ee, aa, ll) :
    k = zero_matrix (parent (ll), 2, 2)
    k[0,0] = k[1,1] = ee * aa / ll
    k[0,1] = k[1,0] = -k[0,0]
    return k
    


##############################    
## Geometry of 2D elements:
##############################


def element2D_length (el_nd_coords) : 
    [[x0, z0], [x1, z1]] = el_nd_coords 
    return sqrt ((x1 - x0)^2 + (z1 - z0)^2)


def element2D_cos (el_nd_coords) :
    [[x0, z0], [x1, z1]] = el_nd_coords
    return (x1 - x0)/element2D_length (el_nd_coords)
    
def element2D_sin (el_nd_coords) :
    [[x0, z0], [x1, z1]] = el_nd_coords
    return -(z1 - z0)/element2D_length (el_nd_coords)

def element2D_cos_sin (el_nd_coords) :
    [[x0, z0], [x1, z1]] = el_nd_coords
    dx = x1 - x0
    dz = z1 - z0
    l = sqrt (dx^2 + dz^2)
    return (dx/l, -dz/l)



##############################    
## Geometry of 3D elements:
##############################


def element3D_length (el_nd_coords) :
    [[x0, y0, z0], [x1, y1, z1]] = el_nd_coords
    return sqrt ((x1 - x0)^2 + (y1 - y0)^2 + (z1 - z0)^2)
    
def element3D_cos_x (el_nd_coords) :
    [[x0, y0, z0], [x1, y1, z1]] = el_nd_coords
    return (x1 - x0)/element3D_length (el_nd_coords)
    
def element3D_cos_y (el_nd_coords) :
    [[x0, y0, z0], [x1, y1, z1]] = el_nd_coords
    return (y1 - y0)/element3D_length (el_nd_coords)

def element3D_cos_z (el_nd_coords) :
    [[x0, y0, z0], [x1, y1, z1]] = el_nd_coords
    return (z1 - z0)/element3D_length (el_nd_coords)

def element3D_cos_all (el_nd_coords) :
    [[x0, y0, z0], [x1, y1, z1]] = el_nd_coords
    dx = x1 - x0
    dy = y1 - y0
    dz = z1 - z0
    l = sqrt (dx^2 + dy^2 + dz^2)
    return (dx/l, dy/l, dz/l)
    


##############################    
## Transformation matrices:
##############################


def truss2D_element_global_to_local (el_nd_coords) : 
    (c, s) = element2D_cos_sin (el_nd_coords)
    rr = zero_matrix (parent (c), 2, 4)
    rr[0,0] = rr[1,2] = c
    rr[0,1] = rr[1,3] = -s
    return rr 
    
    
def truss3D_element_global_to_local (el_nd_coords) :
    (cx, cy, cz) = element3D_cos_all (el_nd_coords)
    rr = zero_matrix (parent (cx), 2, 6)
    rr[0,0] = rr[1,3] = cx
    rr[0,1] = rr[1,4] = cy
    rr[0,2] = rr[1,5] = cz
    return rr



##########################    
## Degrees of freedom:
##########################    


def nodal_DsOF (n_nodes, n_dof_per_node, supports) : 
    dof_table = ones_matrix (ZZ, n_nodes, n_dof_per_node)
    for s in supports.keys() :
        dof_table[s] = supports [s] 
    dof_table -= ones_matrix (ZZ, n_nodes, n_dof_per_node)
    ndof = 0
    for i in xrange (n_nodes) :
        for j in xrange (n_dof_per_node) :
            if dof_table[i, j] == 0 :
                dof_table[i, j] = ndof
                ndof += 1
    return (ndof, dof_table) 
 
    
def element_DsOF (ei, doft) : 
    eix = list (doft [ei[0]])
    eix.extend (doft [ei[1]])
    return eix
    
    
    
################################    
## Structure stiffness matrix:
################################


def structure_stiffness_matrix (ndof, tdof, nds, els, echs, chars) : 
    field = parent (nds[0][0])
    if tdof.ncols() == 2 :
        fl = element2D_length
        fr = truss2D_element_global_to_local
    else :
        fl = element3D_length
        fr = truss3D_element_global_to_local
    sk = matrix (field, ndof, ndof)
    for j in xrange (len (els)) :
        ej = els[j]
        enc = element_node_coords (ej, nds) 
        ll = fl (enc)
        (ee, aa) = ith_element_characteristics (j, echs, chars)
        ekl = truss_element_stiffness_matrix (ee, aa, ll) 
        er = fr (enc)
        ekg = er.T * ekl * er
        eix = element_DsOF (ej, tdof)
        for i in xrange (len (eix)) :
            if eix[i] != -1 : 
                for ii in xrange (len (eix)) : 
                    if eix[ii] != -1 :
                        sk[eix[i], eix[ii]] += ekg[i, ii] 
    return sk
    


##############################    
## Structure force vectors:
##############################


def structure_force_vector (ndof, tdof, nodal_loads) :
    field = parent (nodal_loads.values()[0][0])
    sf = zero_vector (field, ndof)
    for ni in nodal_loads.keys() :
        nlv = nodal_loads[ni]
        for j in xrange (tdof.ncols()) : 
            d = tdof[ni, j]
            if d != -1 :
                sf[d] += nlv[j]  
    return sf



##############################    
## Nodal displacements:
##############################    


def nodal_displacements (sk, sf, tdof, tolerance = None) :
    skf = block_matrix (1, 2, [sk, sf.column()])
    field = skf.base_ring()
    try :
        (field (1)).ulp()
    except : 
        is_zero_f = is_zero
    else :
        is_zero_f = is_zero_wtol (tolerance)
    skf, pl = ref_wpp (skf, is_zero_f, in_place = True)
    if pl[-1] == skf.ncols() - 1 :
        raise Exception ("system has NO solution")
    if len (pl) < skf.ncols() - 1 : 
        print "system has MULTIPLE solutions"
    u = ref_solution (skf, pl)
    return table_of_nodal_displacements (tdof, u)
    
        
def table_of_nodal_displacements (tdof, uu) :
    field = parent (uu[0])
    ndd = zero_matrix (field, tdof.nrows(), tdof.ncols())
    for i in xrange (tdof.nrows()) :
        for j in xrange (tdof.ncols()) : 
            if tdof[i,j] != -1 :
                ndd[i,j] = uu[tdof[i,j]]
    return ndd



###################################    
## Displacements of element ends:
###################################

    
def elem_ends_displacements_global (ei, uu) : 
    edg = list (uu [ei[0]])
    edg.extend (uu [ei[1]])
    return vector (edg)
    
def elem_ends_displacements_local (ei, uu, rr) :
    edg = elem_ends_displacements_global (ei, uu)
    return rr * edg
     

    
############################ 
## Element end forces:
############################


def end_forces (uu, nds, els, elchs, chars) : 
    field = parent (nds[0][0])
    if uu.ncols() == 2 :
        fl = element2D_length
        fr = truss2D_element_global_to_local
    else :
        fl = element3D_length
        fr = truss3D_element_global_to_local
    ne = len (els)
    efs = matrix (field, ne, 2)
    for j in xrange (ne) :
        ej = els[j]
        enc = element_node_coords (ej, nds) 
        ll = fl (enc)
        (ee, aa) = ith_element_characteristics (j, elchs, chars)
        ekl = truss_element_stiffness_matrix (ee, aa, ll) 
        er = fr (enc)
        edl = elem_ends_displacements_local (ej, uu, er)
        efs[j] = ekl * edl
    return efs
     
      
 
##########################
##  Main function (DiM): 
##########################
    
    
def DiM (nodes, elems, supports, chars, elchs, loads, field = RR, tolerance = None) :

    nds2 = coords_left_to_right (nodes_change_field (nodes, field))
    chrs2 = characteristics_change_field (chars, field)
    lds2 = nodal_loads_change_field (loads, field)
    
    ndim = len (nodes[0]) 
    
    if type (supports) == dict : 
        supps_d = supports
    else : 
        supps_d = support2d_dict (supports) if ndim == 2 else support3d_dict (supports)
    
    n, t = nodal_DsOF (len (nds2), ndim, supps_d) 
    sk = structure_stiffness_matrix (n, t, nds2, elems, elchs, chrs2)  
    sf = structure_force_vector (n, t, lds2)
    u = nodal_displacements (sk, sf, t, tolerance)
    ef = end_forces (u, nds2, elems, elchs, chrs2)

    return (u, ef)  
    
    
       
############################
## Row echelon form: 
############################


def is_zero (x) :
    return x == 0

def is_zero_wtol (tolerance = None) :
    def f_ (x) :
        t = sqrt ((parent(x)(1)).ulp()) if tolerance == None else tolerance
        return abs (x) <= t
    return f_


def ref_wpp (A, is_zero_f = is_zero_wtol(), in_place = False) :
    
    Ar = A if in_place == True else copy (A)
        
    r = Ar.nrows()
    c = Ar.ncols() 
    pivot = 0
    pivotlist = []
    
    for j in xrange (c) :
        if pivot == r :
            break
        maxrow = pivot
        for k in xrange (pivot + 1, r):
            if abs (Ar[k, j]) > abs (Ar[maxrow, j]) : 
                maxrow = k
        if pivot != maxrow :
            Ar.swap_rows (pivot, maxrow)
        if not is_zero_f (Ar[pivot, j]) :
            scale = Ar[pivot, j]
            Ar.rescale_row (pivot, 1/scale, j)
            for i in xrange (pivot + 1, r) :
                scale = Ar[i, j]
                Ar.add_multiple_of_row (i, pivot, -scale, j)
            pivotlist.append (j)
            pivot += 1
            
    return Ar, pivotlist


def ref_solution (Ar, pivl) :
     
    b = Ar.column (-1)
    nx = Ar.ncols() - 1
    np = len (pivl)
    x = zero_vector (Ar.base_ring(), nx)
    
    for i in xrange (np - 1, -1, -1) :
        x[pivl[i]] = b[i]
        for j in xrange (np - 1, i, -1) :
            jj = pivl[j]
            x[pivl[i]] -= Ar[i, jj] * x[jj]
            
    return x
   
    
           
##########################
## Utility functions:
##########################

        
def print_list (l, indent = '') : 
    for x in l :
        print '{0}{1}'.format (indent, x)

def print_indexed (l, width = 3) :
    for i in xrange (len(l)) :
        print '{0:>{w}}: {1}'.format (i, l[i], w = width)

def print_dict (d, width = 3) :
    for i in d.keys() :
        print '{0:>{w}}: {1}'.format (i, d[i], w = width)

   
def ND (d) :
    def f_ (x) :
        return N (x, digits = d)
    return f_  
    
    
def chop (tol = None) :
    def f_ (x) :
        t = sqrt ((parent(x) (1)).ulp()) if tol == None else tol
        return x if (abs (x) >= t) else parent(x) (0)
    return f_
    
    
def bind_1st (f, a) : 
    def f_ (y) : 
        return f (a, y)
    return f_
    
def bind_2nd (f, b) :
    def f_ (x) : 
        return f (x, b)
    return f_
   
   
def apply_to_ith (f, i) : 
    def f_ (l) : 
        ll = list (l)    # list() <= copy()
        ll[i] = f(ll[i])
        return ll
    return f_ 
    

def map2 (f, ll) :
    return map (bind_1st (map, f), ll)
 
def map_dict (f, d) :
    return { k : f (d[k]) for k in d.keys() }
    
def map_dict2 (f, dl) :
    return map_dict (bind_1st (map, f), dl)


def all_indices (l) :
    return range (len (l))
 
def others (l1, l2) :
    return [c for c in l2 if c not in l1]    

def list_indexed (inds, l) :
    return [l[i] for i in inds]

 
def nodes_change_field (nds, field = RR) :
    return map2 (field, nds)

def characteristics_change_field (chrs, field = RR) : 
    return map_dict2 (field, chrs)
    
def nodal_loads_change_field (ndlds, field = RR) : 
    return map_dict2 (field, ndlds)
    
    
    
##########################
##  Graphics (2D):
##########################


def plot_hinge (pt, rr, edgecolor = 'blue') :
    return circle (pt, rr, edgecolor = edgecolor, facecolor = 'white',
                   fill = True, aspect_ratio = 1, zorder = 15)

def plot_hinges (pts, rr, edgecolor = 'blue') :
    return sum ([plot_hinge (pt, rr, edgecolor) for pt in pts])


def plot_hinge_wi (i, pt, rr, edgecolor = 'blue', tcolor = 'blue') :
    c = circle (pt, rr, edgecolor = edgecolor, facecolor = 'white',
                fill = True, aspect_ratio = 1, zorder = 15)
    t = text (str (i), (pt[0] + 1.125*rr, pt[1] + 1.125*rr), 
              vertical_alignment = 'bottom', horizontal_alignment = 'left', color = tcolor, zorder = 23)
    return c + t
    
def plot_hinges_wi (pts, rr, edgecolor = 'blue', tcolor = 'blue') :
    return sum ([plot_hinge_wi (i, pts[i], rr, edgecolor, tcolor) for i in xrange (len (pts))]) 
   
    
def plot_element (elnds, color = 'teal', thickness = 4) :
    return line2d (elnds, color = color, thickness = thickness, aspect_ratio = 1)

def plot_elements (elems, nodes, color = 'teal', thickness = 4) :
    return sum ([plot_element (element_node_coords (el, nodes), thickness = thickness, color = color)
                 for el in elems])

def plot_element_wi (i, elnds, rr, color = 'teal', tcolor = 'orange', thickness = 4) :
    l = line2d (elnds, color = color, thickness = thickness, aspect_ratio = 1)
    n0, n1 = elnds
    tp = ((n0[0] + n1[0])/2 + rr, (n0[1] + n1[1])/2 - 1.125*rr)
    t = text (str (i), tp, fontstyle = 'oblique', fontweight = 'semibold', 
              vertical_alignment = 'top', horizontal_alignment = 'left', color = tcolor, zorder = 23)
    return l + t

def plot_elements_wi (elems, nodes, rr, color = 'teal', tcolor = 'orange', thickness = 4) :
    return sum ([plot_element_wi (i, element_node_coords (elems[i], nodes), rr, color, tcolor, thickness)
                 for i in xrange (len (elems))])


def plot_support_x (pt, length, color = 'green', thickness = 3) :
    return line2d ((pt, (pt[0]-length, pt[1])), color = color, thickness = thickness, aspect_ratio = 1)
    
def plot_support_y (pt, length, color = 'green', thickness = 3) :
    return line2d ((pt, (pt[0], pt[1]-length)), color = color, thickness = thickness, aspect_ratio = 1)

def plot_supports (supports, nodes, length, color = 'green', thickness = 2) :
    psps = []
    for ni in supports.keys() :
        sn = nodes[ni]
        sc = supports[ni]
        if sc[0] == 0:
            psps.append (plot_support_x (sn, length, color = color, thickness = thickness))
        if sc[1] == 0:
            psps.append (plot_support_y (sn, length, color = color, thickness = thickness))
    return sum (psps)


def plot_support_hinge (pt, rr) :
    return circle (pt, rr, edgecolor = 'black', facecolor = 'black',
                   fill = True, aspect_ratio = 1, zorder = 15)
                   
def plot_support_hinges (supports, nodes, rr) :
    psps = []
    for ni in supports :
        sn = nodes[ni]
        psps.append (plot_support_hinge (sn, rr))
    return sum (psps)
                   

def plot_force (ff, head, scale, color = 'crimson', width = 4) :
    tail = head[0] - scale*ff[0], head[1] + scale*ff[1]
    return arrow2d (tail, head, width = width, color = color, aspect_ratio = 1, zorder = 13)
    
def plot_forces (loads, nodes, scale, color = 'crimson', width = 4) :
    pls = []
    for ni in loads.keys() :
        head = nodes[ni]
        ff = loads[ni]
        pls.append (plot_force (ff, head, scale, color = color, width = width))
    return sum (pls)

    
def plot_truss2d (nodes, elems, supports, loads, load_scale = 0.01, indices = True) :
    prr = point2d (nodes, aspect_ratio = 1)
    rrx = (prr.xmax() - prr.xmin())
    rry = (prr.ymax() - prr.ymin())
    rr = rrx if rrx > rry else rry
    rn = rr/125
    lsp = rr/23
    if indices == True :
        pt = plot_hinges_wi (nodes, rn)
        pt += plot_elements_wi (elems, nodes, rn)
    else :
        pt = plot_hinges (nodes, rn)
        pt += plot_elements (elems, nodes)
    if type (supports) == list :
        if supports != [] :
            pt += plot_support_hinges (supports, nodes, rn)
    else : 
        if supports != {} :
            pt += plot_supports (supports, nodes, lsp)
    if loads != {} :
        pt += plot_forces (loads, nodes, load_scale)
    return pt
    
    
    
##########################
##  Graphics (3D):
##########################


def plot_hinge3d (pt, color = 'white', size = 1/11) : 
    return sphere (pt, size = size, color = color, aspect_ratio = 1, frame = False)

def plot_hinges3d (pts, color = 'white', size = 1/11) :
    return sum ([plot_hinge3d (pt, color, size) for pt in pts])  


def plot_support3d (pt, color = 'black', size = 1/13) : 
    return sphere (pt, size = size, color = color, aspect_ratio = 1, frame = False)
  
def plot_supports3d (pts, color = 'black', size = 1/13) : 
    return sum ([plot_support3d (pt, color, size) for pt in pts])
    
    
def plot_bar3d (joints, color = 'slategray', radius_s = 23) : 
    ll = element3D_length (joints)
    rr = ll/radius_s
    return line3d (joints, radius = rr, color = color, aspect_ratio = 1, frame = False)
    
def plot_bars3d (bars, joints, color = 'slategray', radius_s = 25) : 
    return sum ([plot_bar3d ((joints[b[0]], joints[b[1]]), color = color, radius_s = radius_s) for b in bars])
    
    
def plot_force3d (ff, joint = (0, 0, 0), scale = 0.01, pull = True, radius_s = 27, color = 'purple', displ = 0.1325) :
    fs = scale * vector (ff)
    nf = norm (fs)
    rd = nf / radius_s
    hl = 7.625 * rd
    hr = 3.125 * rd
    if pull == True :
        tail = vector (joint) + 1.125 * displ * fs
        head = tail + fs
    else :
        head = vector (joint) - displ * fs
        tail = head - fs
    ar = arrow3d (tail, head, radius = rd, head_len = hl, head_radius = hr, color = color, aspect_ratio = 1, frame = False)
    return ar

def plot_forces3d_b__ (forces, joints, scale, pull, radius_s, color, displ) : 
    pfs = []
    for li in forces.keys() :
       pfs.append (plot_force3d (forces[li], joints[li], scale, pull, radius_s, color, displ))
    return sum (pfs)
    
def plot_forces3d_d__ (forces, joints, scale, pull, radius_s, color, displ) : 
    pfs = []
    for li in forces.keys() :
        pfs.append (plot_force3d (forces[li], joints[li], scale, pull[li], radius_s, color, displ))
    return sum (pfs)

def plot_forces3d (forces, joints, scale = 0.01, pull = True, radius_s = 27, color = 'purple', displ = 0.1325) : 
    if type (pull) == bool :
        pfs = plot_forces3d_b__ (forces, joints, scale, pull, radius_s, color, displ)
    else :
        pfs = plot_forces3d_d__ (forces, joints, scale, pull, radius_s, color, displ)
    return pfs 
    
    
def plot_truss3d (joints, bars, supports, forces = {}, force_scale = 0.01, pull = True) : 
    pfjs = plot_hinges3d (list_indexed (others (supports, all_indices (joints)), joints))
    psjs = plot_supports3d (list_indexed (supports, joints))
    pbs = plot_bars3d (bars, joints)
    pt = pfjs + psjs + pbs
    if forces != {} :
        pls = plot_forces3d (forces, joints, force_scale, pull)
        pt += pls
    return pt

nds_xy = [(0., 0.), (3., 0.), (6., 0.), (1.5, 2.5), (4.5, 2.5)]
print 'nodes:'
print_indexed (nds_xy)
print

els = [(0, 1), (1, 2), (0, 3), (3, 1), (1, 4), (4, 2), (3, 4)]
print 'elements:'
print_indexed (els)
print

supps = { 0 : (0, 0), 2 : (1, 0) }
print 'supports:'
print_dict (supps)
print

chars = { 0 : (2*1e8, 0.05*0.05) }
print 'characteristics:'
print_dict (chars)
print

echs = [0] * 7
print 'elements <-> chars:'
print_indexed (echs)
print

lds = { 1 : (0., 100.), 4 : (75., 0.) }
print 'loads:'
print_dict (lds)
print

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (3.00000000000000, 0.000000000000000)
  2: (6.00000000000000, 0.000000000000000)
  3: (1.50000000000000, 2.50000000000000)
  4: (4.50000000000000, 2.50000000000000)

elements:
  0: (0, 1)
  1: (1, 2)
  2: (0, 3)
  3: (3, 1)
  4: (1, 4)
  5: (4, 2)
  6: (3, 4)

supports:
  0: (0, 0)
  2: (1, 0)

characteristics:
  0: (2.00000000000000e8, 0.00250000000000000)

elements <-> chars:
  0: 0
  1: 0
  2: 0
  3: 0
  4: 0
  5: 0
  6: 0

loads:
  1: (0.000000000000000, 100.000000000000)
  4: (75.0000000000000, 0.000000000000000)

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (3.00000000000000, 0.000000000000000)
  2: (6.00000000000000, 0.000000000000000)
  3: (1.50000000000000, 2.50000000000000)
  4: (4.50000000000000, 2.50000000000000)

elements:
  0: (0, 1)
  1: (1, 2)
  2: (0, 3)
  3: (3, 1)
  4: (1, 4)
  5: (4, 2)
  6: (3, 4)

supports:
  0: (0, 0)
  2: (1, 0)

characteristics:
  0: (2.00000000000000e8, 0.00250000000000000)

elements <-> chars:
  0: 0
  1: 0
  2: 0
  3: 0
  4: 0
  5: 0
  6: 0

loads:
  1: (0.000000000000000, 100.000000000000)
  4: (75.0000000000000, 0.000000000000000)

plot_truss2d (nds_xy, els, supps, lds, indices = False)

plot_truss2d (nds_xy, els, supps, lds, 0.0075)

nds = coords_left_to_right (nds_xy)
print_list (nds)

[0.000000000000000, -0.000000000000000]
[3.00000000000000, -0.000000000000000]
[6.00000000000000, -0.000000000000000]
[1.50000000000000, -2.50000000000000]
[4.50000000000000, -2.50000000000000]

[0.000000000000000, -0.000000000000000]
[3.00000000000000, -0.000000000000000]
[6.00000000000000, -0.000000000000000]
[1.50000000000000, -2.50000000000000]
[4.50000000000000, -2.50000000000000]

ndof, tdof = nodal_DsOF (len (nds), 2, supps)
ndof; tdof

7
[-1 -1]
[ 0  1]
[ 2 -1]
[ 3  4]
[ 5  6]

7
[-1 -1]
[ 0  1]
[ 2 -1]
[ 3  4]
[ 5  6]

ec2 = element_node_coords (els[2], nds)
ec2

[[0.000000000000000, -0.000000000000000],
 [1.50000000000000, -2.50000000000000]]

[[0.000000000000000, -0.000000000000000],
 [1.50000000000000, -2.50000000000000]]

l2 = element2D_length (ec2)
l2

2.91547594742265

2.91547594742265

ee, aa = chars[0]
ee; aa

2.00000000000000e8
0.00250000000000000

2.00000000000000e8
0.00250000000000000

k2 = truss_element_stiffness_matrix (ee, aa, l2)
show (k2)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 171498.585142509 & -171498.585142509 \\ -171498.585142509 & 171498.585142509 \end{array}\right)$

r2 = truss2D_element_global_to_local (ec2)
show (r2)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 0.514495755427526 & -0.857492925712544 & 0.000000000000000 & 0.000000000000000 \\ 0.000000000000000 & 0.000000000000000 & 0.514495755427526 & -0.857492925712544 \end{array}\right)$

k2g = r2.T * k2 * r2
show (k2g)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 45396.6843024288 & -75661.1405040480 & -45396.6843024288 & 75661.1405040480 \\ -75661.1405040480 & 126101.900840080 & 75661.1405040480 & -126101.900840080 \\ -45396.6843024288 & 75661.1405040480 & 45396.6843024288 & -75661.1405040480 \\ 75661.1405040480 & -126101.900840080 & -75661.1405040480 & 126101.900840080 \end{array}\right)$

K = structure_stiffness_matrix (ndof, tdof, nds, els, echs, chars)
show (K.apply_map (ND(8)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 424126.70 & 0.00000000 & -166666.67 & -45396.684 & -75661.141 & -45396.684 & 75661.141 \\ 0.00000000 & 252203.80 & 0.00000000 & -75661.141 & -126101.90 & 75661.141 & -126101.90 \\ -166666.67 & 0.00000000 & 212063.35 & 0.00000000 & 0.00000000 & -45396.684 & -75661.141 \\ -45396.684 & -75661.141 & 0.00000000 & 257460.04 & 0.00000000 & -166666.67 & 0.00000000 \\ -75661.141 & -126101.90 & 0.00000000 & 0.00000000 & 252203.80 & 0.00000000 & 0.00000000 \\ -45396.684 & 75661.141 & -45396.684 & -166666.67 & 0.00000000 & 257460.04 & 0.00000000 \\ 75661.141 & -126101.90 & -75661.141 & 0.00000000 & 0.00000000 & 0.00000000 & 252203.80 \end{array}\right)$

q = structure_force_vector (ndof, tdof, lds)
show (q.apply_map (ND(8)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.00000000,\,100.00000,\,0.00000000,\,0.00000000,\,0.00000000,\,75.000000,\,0.00000000\right)$

u = nodal_displacements (K, q, tdof)
print_indexed (u.rows())

  0: (0.000000000000000, 0.000000000000000)
  1: (0.000517500000000000, 0.00111700945769896)
  2: (0.000810000000000000, 0.000000000000000)
  3: (0.000941775759218208, 0.000713754728849480)
  4: (0.000806775759218209, 0.000646254728849480)

  0: (0.000000000000000, 0.000000000000000)
  1: (0.000517500000000000, 0.00111700945769896)
  2: (0.000810000000000000, 0.000000000000000)
  3: (0.000941775759218208, 0.000713754728849480)
  4: (0.000806775759218209, 0.000646254728849480)

f = end_forces (u, nds, els, echs, chars)
print_indexed (f.rows())

  0: (-86.2500000000000, 86.2500000000000)
  1: (-48.7500000000000, 48.7500000000000)
  2: (21.8660696056698, -21.8660696056698)
  3: (-21.8660696056699, 21.8660696056699)
  4: (-94.7529682912361, 94.7529682912361)
  5: (94.7529682912361, -94.7529682912361)
  6: (22.4999999999999, -22.4999999999999)

  0: (-86.2500000000000, 86.2500000000000)
  1: (-48.7500000000000, 48.7500000000000)
  2: (21.8660696056698, -21.8660696056698)
  3: (-21.8660696056699, 21.8660696056699)
  4: (-94.7529682912361, 94.7529682912361)
  5: (94.7529682912361, -94.7529682912361)
  6: (22.4999999999999, -22.4999999999999)

uu, ff = DiM (nds_xy, els, supps, chars, echs, lds)

print_indexed (uu.apply_map (ND(8)).rows())

  0: (0.00000000, 0.00000000)
  1: (0.00051750000, 0.0011170095)
  2: (0.00081000000, 0.00000000)
  3: (0.00094177576, 0.00071375473)
  4: (0.00080677576, 0.00064625473)

  0: (0.00000000, 0.00000000)
  1: (0.00051750000, 0.0011170095)
  2: (0.00081000000, 0.00000000)
  3: (0.00094177576, 0.00071375473)
  4: (0.00080677576, 0.00064625473)

print_indexed (ff.apply_map (ND(8)).rows())

  0: (-86.250000, 86.250000)
  1: (-48.750000, 48.750000)
  2: (21.866070, -21.866070)
  3: (-21.866070, 21.866070)
  4: (-94.752968, 94.752968)
  5: (94.752968, -94.752968)
  6: (22.500000, -22.500000)

  0: (-86.250000, 86.250000)
  1: (-48.750000, 48.750000)
  2: (21.866070, -21.866070)
  3: (-21.866070, 21.866070)
  4: (-94.752968, 94.752968)
  5: (94.752968, -94.752968)
  6: (22.500000, -22.500000)

uQQ, fQQ = DiM (nds_xy, els, supps, chars, echs, lds, field = QQ)

Traceback (click to the left of this block for traceback)
...
TypeError: unable to convert
16012798675095600000000/1028377514911694911*sqrt(17/2) to a rational

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_22.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("dVFRLCBmUVEgPSBEaU0gKG5kc194eSwgZWxzLCBzdXBwcywgY2hhcnMsIGVjaHMsIGxkcywgZmllbGQgPSBRUSk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpein5Pm/___code___.py", line 2, in <module>
    exec compile(u'uQQ, fQQ = DiM (nds_xy, els, supps, chars, echs, lds, field = QQ)
  File "", line 1, in <module>
    
  File "/tmp/tmpgt8if5/___code___.py", line 314, in DiM
    sk = structure_stiffness_matrix (n, t, nds2, elems, elchs, chrs2)  
  File "/tmp/tmpgt8if5/___code___.py", line 194, in structure_stiffness_matrix
    sk[eix[i], eix[ii]] += ekg[i, ii] 
  File "sage/matrix/matrix0.pyx", line 1424, in sage.matrix.matrix0.Matrix.__setitem__ (build/cythonized/sage/matrix/matrix0.c:8803)
  File "sage/matrix/matrix0.pyx", line 1529, in sage.matrix.matrix0.Matrix._coerce_element (build/cythonized/sage/matrix/matrix0.c:10280)
  File "sage/structure/parent.pyx", line 900, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9197)
  File "sage/structure/coerce_maps.pyx", line 162, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4556)
  File "sage/structure/coerce_maps.pyx", line 157, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4448)
  File "sage/rings/rational.pyx", line 531, in sage.rings.rational.Rational.__init__ (build/cythonized/sage/rings/rational.c:6395)
  File "sage/rings/rational.pyx", line 619, in sage.rings.rational.Rational.__set_value (build/cythonized/sage/rings/rational.c:7516)
  File "sage/symbolic/expression.pyx", line 1163, in sage.symbolic.expression.Expression._rational_ (build/cythonized/sage/symbolic/expression.cpp:9379)
TypeError: unable to convert 16012798675095600000000/1028377514911694911*sqrt(17/2) to a rational

gr = (1 + sqrt (5))/2
gr

1/2*sqrt(5) + 1/2

1/2*sqrt(5) + 1/2

grAA = AA (gr)
grAA

1.618033988749895?

1.618033988749895?

grAA.minpoly()

x^2 - x - 1

x^2 - x - 1

gr2 = (1 - sqrt (5))/2
AA (gr2); (AA (gr2)).minpoly()

-0.618033988749895?
x^2 - x - 1

-0.618033988749895?
x^2 - x - 1

RIF (grAA)

1.618033988749895?

1.618033988749895?

(RIF (grAA)).str (style='brackets')

'[1.6180339887498946 .. 1.6180339887498950]'

'[1.6180339887498946 .. 1.6180339887498950]'

uAA, fAA = DiM (nds_xy, els, supps, chars, echs, lds, field = AA)

Traceback (click to the left of this block for traceback)
...
TypeError: Illegal initializer for algebraic number

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_29.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("dUFBLCBmQUEgPSBEaU0gKG5kc194eSwgZWxzLCBzdXBwcywgY2hhcnMsIGVjaHMsIGxkcywgZmllbGQgPSBBQSk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpTuANht/___code___.py", line 2, in <module>
    exec compile(u'uAA, fAA = DiM (nds_xy, els, supps, chars, echs, lds, field = AA)
  File "", line 1, in <module>
    
  File "/tmp/tmpgt8if5/___code___.py", line 302, in DiM
    nds2 = coords_left_to_right (nodes_change_field (nodes, field))
  File "/tmp/tmpgt8if5/___code___.py", line 456, in nodes_change_field
    return map2 (field, nds)
  File "/tmp/tmpgt8if5/___code___.py", line 436, in map2
    return map (bind_1st (map, f), ll)
  File "/tmp/tmpgt8if5/___code___.py", line 418, in f_
    return f (a, y)
  File "sage/structure/parent.pyx", line 900, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9197)
  File "sage/structure/coerce_maps.pyx", line 162, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4556)
  File "sage/structure/coerce_maps.pyx", line 157, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4448)
  File "/opt/SageMath/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 760, in _element_constructor_
    return AlgebraicReal(x)
  File "/opt/SageMath/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 4739, in __init__
    AlgebraicNumber_base.__init__(self, AA, x)
  File "/opt/SageMath/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3106, in __init__
    raise TypeError("Illegal initializer for algebraic number")
TypeError: Illegal initializer for algebraic number

AA (1/2)

1/2

1/2

AA (0.5)

Traceback (click to the left of this block for traceback)
...
TypeError: Illegal initializer for algebraic number

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_31.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("QUEgKDAuNSk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpZjg7r7/___code___.py", line 3, in <module>
    exec compile(u'AA (_sage_const_0p5 )
  File "", line 1, in <module>
    
  File "sage/structure/parent.pyx", line 900, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9197)
  File "sage/structure/coerce_maps.pyx", line 162, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4556)
  File "sage/structure/coerce_maps.pyx", line 157, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4448)
  File "/opt/SageMath/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 760, in _element_constructor_
    return AlgebraicReal(x)
  File "/opt/SageMath/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 4739, in __init__
    AlgebraicNumber_base.__init__(self, AA, x)
  File "/opt/SageMath/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3106, in __init__
    raise TypeError("Illegal initializer for algebraic number")
TypeError: Illegal initializer for algebraic number

QQ (0.5)

1/2

1/2

AA (QQ (0.5))

1/2

1/2

nds_xyQQ = nodes_change_field (nds_xy, QQ)
charsQQ = characteristics_change_field (chars, QQ)
ldsQQ = nodal_loads_change_field (lds, QQ)
print_indexed (nds_xyQQ)
print
print_dict (charsQQ)
print
print_dict (ldsQQ)

  0: [0, 0]
  1: [3, 0]
  2: [6, 0]
  3: [3/2, 5/2]
  4: [9/2, 5/2]

  0: [200000000, 8895999263942/3558399705576799]

  1: [0, 100]
  4: [75, 0]

  0: [0, 0]
  1: [3, 0]
  2: [6, 0]
  3: [3/2, 5/2]
  4: [9/2, 5/2]

  0: [200000000, 8895999263942/3558399705576799]

  1: [0, 100]
  4: [75, 0]

(AA (1)).ulp()

Traceback (click to the left of this block for traceback)
...
AttributeError: 'AlgebraicReal' object has no attribute 'ulp'

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_35.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("KEFBICgxKSkudWxwKCk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpzQgPTQ/___code___.py", line 3, in <module>
    exec compile(u'(AA (_sage_const_1 )).ulp()
  File "", line 1, in <module>
    
  File "sage/structure/element.pyx", line 489, in sage.structure.element.Element.__getattr__ (build/cythonized/sage/structure/element.c:4607)
  File "sage/structure/element.pyx", line 502, in sage.structure.element.Element.getattr_from_category (build/cythonized/sage/structure/element.c:4716)
  File "sage/cpython/getattr.pyx", line 394, in sage.cpython.getattr.getattr_from_other_class (build/cythonized/sage/cpython/getattr.c:2614)
AttributeError: 'AlgebraicReal' object has no attribute 'ulp'

uAA, fAA = DiM (nds_xyQQ, els, supps, charsQQ, echs, ldsQQ, field = AA)

print_indexed (uAA.rows())

  0: (0, 0)
  1: (0.0005174999999999999?, 0.001117009457698961?)
  2: (0.0008099999999999998?, 0)
  3: (0.0009417757592182085?, 0.0007137547288494802?)
  4: (0.0008067757592182086?, 0.000646254728849481?)

  0: (0, 0)
  1: (0.0005174999999999999?, 0.001117009457698961?)
  2: (0.0008099999999999998?, 0)
  3: (0.0009417757592182085?, 0.0007137547288494802?)
  4: (0.0008067757592182086?, 0.000646254728849481?)

for i in xrange (uAA.nrows()) :
    print i, ': ', 
    for j in xrange (uAA.ncols()) :
        print (RIF (uAA[i,j])).str (style='brackets'), 
    print

0 :  [0.0000000000000000 .. -0.0000000000000000] [0.0000000000000000 ..
-0.0000000000000000]
1 :  [0.00051749999999999983 .. 0.00051749999999999995] [0.0011170094576989605
.. 0.0011170094576989608]
2 :  [0.00080999999999999973 .. 0.00080999999999999985] [0.0000000000000000 ..
-0.0000000000000000]
3 :  [0.00094177575921820843 .. 0.00094177575921820854] [0.00071375472884948015
.. 0.00071375472884948027]
4 :  [0.00080677575921820851 .. 0.00080677575921820862] [0.00064625472884948019
.. 0.00064625472884948031]

0 :  [0.0000000000000000 .. -0.0000000000000000] [0.0000000000000000 .. -0.0000000000000000]
1 :  [0.00051749999999999983 .. 0.00051749999999999995] [0.0011170094576989605 .. 0.0011170094576989608]
2 :  [0.00080999999999999973 .. 0.00080999999999999985] [0.0000000000000000 .. -0.0000000000000000]
3 :  [0.00094177575921820843 .. 0.00094177575921820854] [0.00071375472884948015 .. 0.00071375472884948027]
4 :  [0.00080677575921820851 .. 0.00080677575921820862] [0.00064625472884948019 .. 0.00064625472884948031]

print_indexed (uAA.apply_map (RR).rows())

  0: (0.000000000000000, 0.000000000000000)
  1: (0.000517500000000000, 0.00111700945769896)
  2: (0.000810000000000000, 0.000000000000000)
  3: (0.000941775759218208, 0.000713754728849480)
  4: (0.000806775759218209, 0.000646254728849480)

  0: (0.000000000000000, 0.000000000000000)
  1: (0.000517500000000000, 0.00111700945769896)
  2: (0.000810000000000000, 0.000000000000000)
  3: (0.000941775759218208, 0.000713754728849480)
  4: (0.000806775759218209, 0.000646254728849480)

print_indexed (uAA.apply_map (lambda x : x.minpoly()).rows())

  0: (x, x)
  1: (x - 736588739054397393/1423359882230720000000, x^2 -
288230376151720719/444799963197100000000*x -
414598691758003837870589603070027943/791388029040566057813793640000000000000000)
  2: (x - 288230376151720719/355839970557680000000, x)
  3: (x^2 - 1505203075458985977/1423359882230720000000*x +
7948840647976552156658803402646178161/729343207563785678881192218624000000000000\
00, x^2 - 4515609226376957931/7116799411153600000000*x -
2292100314062306347234409167093222619/405190670868769821600662343680000000000000\
00)
  4: (x^2 - 32025597350191191/40667425206592000000*x -
45368692928837545462856871630574783/2917372830255142715524768874496000000000000,
x^2 - 3554841305871222201/7116799411153600000000*x -
3842866309144804775346886639138925891/405190670868769821600662343680000000000000\
00)

  0: (x, x)
  1: (x - 736588739054397393/1423359882230720000000, x^2 - 288230376151720719/444799963197100000000*x - 414598691758003837870589603070027943/791388029040566057813793640000000000000000)
  2: (x - 288230376151720719/355839970557680000000, x)
  3: (x^2 - 1505203075458985977/1423359882230720000000*x + 7948840647976552156658803402646178161/72934320756378567888119221862400000000000000, x^2 - 4515609226376957931/7116799411153600000000*x - 2292100314062306347234409167093222619/40519067086876982160066234368000000000000000)
  4: (x^2 - 32025597350191191/40667425206592000000*x - 45368692928837545462856871630574783/2917372830255142715524768874496000000000000, x^2 - 3554841305871222201/7116799411153600000000*x - 3842866309144804775346886639138925891/40519067086876982160066234368000000000000000)

print_indexed (fAA.rows())

  0: (-86.2500000000000?, 86.2500000000000?)
  1: (-48.75000000000000?, 48.75000000000000?)
  2: (21.86606960566988?, -21.86606960566988?)
  3: (-21.86606960566988?, 21.86606960566988?)
  4: (-94.7529682912362?, 94.7529682912362?)
  5: (94.7529682912362?, -94.7529682912362?)
  6: (22.50000000000000?, -22.50000000000000?)

  0: (-86.2500000000000?, 86.2500000000000?)
  1: (-48.75000000000000?, 48.75000000000000?)
  2: (21.86606960566988?, -21.86606960566988?)
  3: (-21.86606960566988?, 21.86606960566988?)
  4: (-94.7529682912362?, 94.7529682912362?)
  5: (94.7529682912362?, -94.7529682912362?)
  6: (22.50000000000000?, -22.50000000000000?)

for i in xrange (fAA.nrows()) :
    print i, ': ', 
    for j in xrange (fAA.ncols()) :
        print (RIF (fAA[i,j])).str (style='brackets'), 
    print

0 :  [-86.250000000000015 .. -86.249999999999985] [86.249999999999985 ..
86.250000000000015]
1 :  [-48.750000000000008 .. -48.749999999999992] [48.749999999999992 ..
48.750000000000008]
2 :  [21.866069605669874 .. 21.866069605669879] [-21.866069605669879 ..
-21.866069605669874]
3 :  [-21.866069605669879 .. -21.866069605669874] [21.866069605669874 ..
21.866069605669879]
4 :  [-94.752968291236144 .. -94.752968291236115] [94.752968291236115 ..
94.752968291236144]
5 :  [94.752968291236129 .. 94.752968291236144] [-94.752968291236144 ..
-94.752968291236129]
6 :  [22.499999999999996 .. 22.500000000000004] [-22.500000000000004 ..
-22.499999999999996]

0 :  [-86.250000000000015 .. -86.249999999999985] [86.249999999999985 .. 86.250000000000015]
1 :  [-48.750000000000008 .. -48.749999999999992] [48.749999999999992 .. 48.750000000000008]
2 :  [21.866069605669874 .. 21.866069605669879] [-21.866069605669879 .. -21.866069605669874]
3 :  [-21.866069605669879 .. -21.866069605669874] [21.866069605669874 .. 21.866069605669879]
4 :  [-94.752968291236144 .. -94.752968291236115] [94.752968291236115 .. 94.752968291236144]
5 :  [94.752968291236129 .. 94.752968291236144] [-94.752968291236144 .. -94.752968291236129]
6 :  [22.499999999999996 .. 22.500000000000004] [-22.500000000000004 .. -22.499999999999996]

print_indexed (fAA.apply_map (RR).rows())

  0: (-86.2500000000000, 86.2500000000000)
  1: (-48.7500000000000, 48.7500000000000)
  2: (21.8660696056699, -21.8660696056699)
  3: (-21.8660696056699, 21.8660696056699)
  4: (-94.7529682912361, 94.7529682912361)
  5: (94.7529682912361, -94.7529682912361)
  6: (22.5000000000000, -22.5000000000000)

  0: (-86.2500000000000, 86.2500000000000)
  1: (-48.7500000000000, 48.7500000000000)
  2: (21.8660696056699, -21.8660696056699)
  3: (-21.8660696056699, 21.8660696056699)
  4: (-94.7529682912361, 94.7529682912361)
  5: (94.7529682912361, -94.7529682912361)
  6: (22.5000000000000, -22.5000000000000)

joints = [ (1.6, 3.0), (3*1.6, 3.0), (5*1.6, 3.0), (0.0, 0.0), (3.2, 0.0), (2*3.2, 0.0), (3*3.2, 0.0) ]
print 'joints:'
print_indexed (joints)
print

bars = [ (0, 1), (1, 2), (3, 0), (0, 4), (4, 1), (1, 5), (5, 2), (2, 6), (3, 4), (4, 5), (5, 6) ] 
print 'bars:' 
print_indexed (bars)
print

supports = { 3 : (0, 0), 5 : (1, 0) }
print 'supports:'
print_dict (supports)
print

chars = { 0 : (2*1e8, 0.05*0.05) } 
print 'characteristics:'
print_dict (chars)
print

echs = [0] * len (bars)
print 'elements <-> chars:'
print_indexed (echs)
print

loads = { 1 : (0.0, -22.0), 6 : (22.4, 16.8) } 
print 'loads:'
print_dict (loads)
print

joints:
  0: (1.60000000000000, 3.00000000000000)
  1: (4.80000000000000, 3.00000000000000)
  2: (8.00000000000000, 3.00000000000000)
  3: (0.000000000000000, 0.000000000000000)
  4: (3.20000000000000, 0.000000000000000)
  5: (6.40000000000000, 0.000000000000000)
  6: (9.60000000000000, 0.000000000000000)

bars:
  0: (0, 1)
  1: (1, 2)
  2: (3, 0)
  3: (0, 4)
  4: (4, 1)
  5: (1, 5)
  6: (5, 2)
  7: (2, 6)
  8: (3, 4)
  9: (4, 5)
 10: (5, 6)

supports:
  3: (0, 0)
  5: (1, 0)

characteristics:
  0: (2.00000000000000e8, 0.00250000000000000)

elements <-> chars:
  0: 0
  1: 0
  2: 0
  3: 0
  4: 0
  5: 0
  6: 0
  7: 0
  8: 0
  9: 0
 10: 0

loads:
  1: (0.000000000000000, -22.0000000000000)
  6: (22.4000000000000, 16.8000000000000)

joints:
  0: (1.60000000000000, 3.00000000000000)
  1: (4.80000000000000, 3.00000000000000)
  2: (8.00000000000000, 3.00000000000000)
  3: (0.000000000000000, 0.000000000000000)
  4: (3.20000000000000, 0.000000000000000)
  5: (6.40000000000000, 0.000000000000000)
  6: (9.60000000000000, 0.000000000000000)

bars:
  0: (0, 1)
  1: (1, 2)
  2: (3, 0)
  3: (0, 4)
  4: (4, 1)
  5: (1, 5)
  6: (5, 2)
  7: (2, 6)
  8: (3, 4)
  9: (4, 5)
 10: (5, 6)

supports:
  3: (0, 0)
  5: (1, 0)

characteristics:
  0: (2.00000000000000e8, 0.00250000000000000)

elements <-> chars:
  0: 0
  1: 0
  2: 0
  3: 0
  4: 0
  5: 0
  6: 0
  7: 0
  8: 0
  9: 0
 10: 0

loads:
  1: (0.000000000000000, -22.0000000000000)
  6: (22.4000000000000, 16.8000000000000)

plot_truss2d (joints, bars, supports, loads, 0.0425)

uu, ff = DiM (joints, bars, supports, chars, echs, loads)

print_indexed (ff.rows())

  0: (-14.8266666666667, 14.8266666666667)
  1: (-17.9200000000000, 17.9200000000000)
  2: (-15.7533333333333, 15.7533333333333)
  3: (15.7533333333333, -15.7533333333333)
  4: (-15.7533333333333, 15.7533333333333)
  5: (-9.18000000000000, 9.18000000000000)
  6: (19.0400000000000, -19.0400000000000)
  7: (-19.0400000000000, 19.0400000000000)
  8: (-14.9866666666667, 14.9866666666667)
  9: (-0.159999999999982, 0.159999999999982)
 10: (-13.4400000000000, 13.4400000000000)

  0: (-14.8266666666667, 14.8266666666667)
  1: (-17.9200000000000, 17.9200000000000)
  2: (-15.7533333333333, 15.7533333333333)
  3: (15.7533333333333, -15.7533333333333)
  4: (-15.7533333333333, 15.7533333333333)
  5: (-9.18000000000000, 9.18000000000000)
  6: (19.0400000000000, -19.0400000000000)
  7: (-19.0400000000000, 19.0400000000000)
  8: (-14.9866666666667, 14.9866666666667)
  9: (-0.159999999999982, 0.159999999999982)
 10: (-13.4400000000000, 13.4400000000000)

joints = [ (1.6, 3.0), (3*1.6, 3.0), (5*1.6, 3.0), (0.0, 0.0), (3.2, 0.0), (2*3.2, 0.0), (3*3.2, 0.0), (-0.725, 0.0), (0.0, -0.725), (2*3.2, -0.725) ]
print 'joints:'
print_indexed (joints)
print

bars = [ (0, 1), (1, 2), (3, 0), (0, 4), (4, 1), (1, 5), (5, 2), (2, 6), (3, 4), (4, 5), (5, 6), (7, 3), (3, 8), (9, 5) ]
print 'bars:' 
print_indexed (bars)
print

supports = { 7 : (0, 0), 8 : (0, 0), 9 : (0, 0) }
print 'supports:' 
print_dict (supports)
print

chars = { 0 : (2*1e8, 0.05*0.05) } 
print 'characteristics:'
print_dict (chars)
print

echs = [0] * len (bars)
print 'elements <-> chars:'
print_indexed (echs)
print

loads = { 1 : (0.0, -22.0), 6 : (22.4, 16.8) } 
print 'loads:'
print_dict (loads)
print

joints:
  0: (1.60000000000000, 3.00000000000000)
  1: (4.80000000000000, 3.00000000000000)
  2: (8.00000000000000, 3.00000000000000)
  3: (0.000000000000000, 0.000000000000000)
  4: (3.20000000000000, 0.000000000000000)
  5: (6.40000000000000, 0.000000000000000)
  6: (9.60000000000000, 0.000000000000000)
  7: (-0.725000000000000, 0.000000000000000)
  8: (0.000000000000000, -0.725000000000000)
  9: (6.40000000000000, -0.725000000000000)

bars:
  0: (0, 1)
  1: (1, 2)
  2: (3, 0)
  3: (0, 4)
  4: (4, 1)
  5: (1, 5)
  6: (5, 2)
  7: (2, 6)
  8: (3, 4)
  9: (4, 5)
 10: (5, 6)
 11: (7, 3)
 12: (3, 8)
 13: (9, 5)

supports:
  8: (0, 0)
  9: (0, 0)
  7: (0, 0)

characteristics:
  0: (2.00000000000000e8, 0.00250000000000000)

elements <-> chars:
  0: 0
  1: 0
  2: 0
  3: 0
  4: 0
  5: 0
  6: 0
  7: 0
  8: 0
  9: 0
 10: 0
 11: 0
 12: 0
 13: 0

loads:
  1: (0.000000000000000, -22.0000000000000)
  6: (22.4000000000000, 16.8000000000000)

joints:
  0: (1.60000000000000, 3.00000000000000)
  1: (4.80000000000000, 3.00000000000000)
  2: (8.00000000000000, 3.00000000000000)
  3: (0.000000000000000, 0.000000000000000)
  4: (3.20000000000000, 0.000000000000000)
  5: (6.40000000000000, 0.000000000000000)
  6: (9.60000000000000, 0.000000000000000)
  7: (-0.725000000000000, 0.000000000000000)
  8: (0.000000000000000, -0.725000000000000)
  9: (6.40000000000000, -0.725000000000000)

bars:
  0: (0, 1)
  1: (1, 2)
  2: (3, 0)
  3: (0, 4)
  4: (4, 1)
  5: (1, 5)
  6: (5, 2)
  7: (2, 6)
  8: (3, 4)
  9: (4, 5)
 10: (5, 6)
 11: (7, 3)
 12: (3, 8)
 13: (9, 5)

supports:
  8: (0, 0)
  9: (0, 0)
  7: (0, 0)

characteristics:
  0: (2.00000000000000e8, 0.00250000000000000)

elements <-> chars:
  0: 0
  1: 0
  2: 0
  3: 0
  4: 0
  5: 0
  6: 0
  7: 0
  8: 0
  9: 0
 10: 0
 11: 0
 12: 0
 13: 0

loads:
  1: (0.000000000000000, -22.0000000000000)
  6: (22.4000000000000, 16.8000000000000)

plot_truss2d (joints, bars, supports, loads, 0.0425)

uu, ff = DiM (joints, bars, supports, chars, echs, loads)

print_indexed (ff.rows())

  0: (-14.8266666666667, 14.8266666666667)
  1: (-17.9200000000000, 17.9200000000000)
  2: (-15.7533333333333, 15.7533333333333)
  3: (15.7533333333333, -15.7533333333333)
  4: (-15.7533333333333, 15.7533333333333)
  5: (-9.17999999999999, 9.17999999999999)
  6: (19.0400000000000, -19.0400000000000)
  7: (-19.0400000000000, 19.0400000000000)
  8: (-14.9866666666667, 14.9866666666667)
  9: (-0.159999999999989, 0.159999999999989)
 10: (-13.4400000000000, 13.4400000000000)
 11: (-22.4000000000000, 22.4000000000000)
 12: (-13.9000000000000, 13.9000000000000)
 13: (8.70000000000000, -8.70000000000000)

  0: (-14.8266666666667, 14.8266666666667)
  1: (-17.9200000000000, 17.9200000000000)
  2: (-15.7533333333333, 15.7533333333333)
  3: (15.7533333333333, -15.7533333333333)
  4: (-15.7533333333333, 15.7533333333333)
  5: (-9.17999999999999, 9.17999999999999)
  6: (19.0400000000000, -19.0400000000000)
  7: (-19.0400000000000, 19.0400000000000)
  8: (-14.9866666666667, 14.9866666666667)
  9: (-0.159999999999989, 0.159999999999989)
 10: (-13.4400000000000, 13.4400000000000)
 11: (-22.4000000000000, 22.4000000000000)
 12: (-13.9000000000000, 13.9000000000000)
 13: (8.70000000000000, -8.70000000000000)

supps_l = support_list (supports)
supps_l

[8, 9, 7]

[8, 9, 7]

print_dict (support2d_dict (supps_l))

  8: (0, 0)
  9: (0, 0)
  7: (0, 0)

  8: (0, 0)
  9: (0, 0)
  7: (0, 0)

plot_truss2d (joints, bars, supps_l, loads, 0.0425)

uu, ff = DiM (joints, bars, supps_l, chars, echs, loads)

print_indexed (ff.rows())

  0: (-14.8266666666667, 14.8266666666667)
  1: (-17.9200000000000, 17.9200000000000)
  2: (-15.7533333333333, 15.7533333333333)
  3: (15.7533333333333, -15.7533333333333)
  4: (-15.7533333333333, 15.7533333333333)
  5: (-9.17999999999999, 9.17999999999999)
  6: (19.0400000000000, -19.0400000000000)
  7: (-19.0400000000000, 19.0400000000000)
  8: (-14.9866666666667, 14.9866666666667)
  9: (-0.159999999999989, 0.159999999999989)
 10: (-13.4400000000000, 13.4400000000000)
 11: (-22.4000000000000, 22.4000000000000)
 12: (-13.9000000000000, 13.9000000000000)
 13: (8.70000000000000, -8.70000000000000)

  0: (-14.8266666666667, 14.8266666666667)
  1: (-17.9200000000000, 17.9200000000000)
  2: (-15.7533333333333, 15.7533333333333)
  3: (15.7533333333333, -15.7533333333333)
  4: (-15.7533333333333, 15.7533333333333)
  5: (-9.17999999999999, 9.17999999999999)
  6: (19.0400000000000, -19.0400000000000)
  7: (-19.0400000000000, 19.0400000000000)
  8: (-14.9866666666667, 14.9866666666667)
  9: (-0.159999999999989, 0.159999999999989)
 10: (-13.4400000000000, 13.4400000000000)
 11: (-22.4000000000000, 22.4000000000000)
 12: (-13.9000000000000, 13.9000000000000)
 13: (8.70000000000000, -8.70000000000000)

aa = tan (RR (pi)/3.)

joints = [(0., -1., 0.), (-1., 0., 0.), (0., 1., 0.), (2., -1., 0.), (2., 1., 0.), (0., 0., aa), (2., 0., aa)]
print 'joints:'
print_indexed (joints)
print

bars = [(0, 5), (1, 5), (2, 5), (3, 6), (4, 6), (5, 6)]
print 'bars:'
print_indexed (bars)
print

supports = [0, 1, 2, 3, 4]
print 'supports:'
print_list (supports, '  ')
print

chars = { 0 : (2*1e8, 0.05*0.05) }
print 'characteristics:'
print_dict (chars)
print

echs = [0] * len (bars)
print 'elements <-> chars:'
print_indexed (echs)
print

loads = { 6 : (-100., 0., 0.) }
print 'loads:'
print_dict (loads)
print

joints:
  0: (0.000000000000000, -1.00000000000000, 0.000000000000000)
  1: (-1.00000000000000, 0.000000000000000, 0.000000000000000)
  2: (0.000000000000000, 1.00000000000000, 0.000000000000000)
  3: (2.00000000000000, -1.00000000000000, 0.000000000000000)
  4: (2.00000000000000, 1.00000000000000, 0.000000000000000)
  5: (0.000000000000000, 0.000000000000000, 1.73205080756888)
  6: (2.00000000000000, 0.000000000000000, 1.73205080756888)

bars:
  0: (0, 5)
  1: (1, 5)
  2: (2, 5)
  3: (3, 6)
  4: (4, 6)
  5: (5, 6)

supports:
  0
  1
  2
  3
  4

characteristics:
  0: (2.00000000000000e8, 0.00250000000000000)

elements <-> chars:
  0: 0
  1: 0
  2: 0
  3: 0
  4: 0
  5: 0

loads:
  6: (-100.000000000000, 0.000000000000000, 0.000000000000000)

joints:
  0: (0.000000000000000, -1.00000000000000, 0.000000000000000)
  1: (-1.00000000000000, 0.000000000000000, 0.000000000000000)
  2: (0.000000000000000, 1.00000000000000, 0.000000000000000)
  3: (2.00000000000000, -1.00000000000000, 0.000000000000000)
  4: (2.00000000000000, 1.00000000000000, 0.000000000000000)
  5: (0.000000000000000, 0.000000000000000, 1.73205080756888)
  6: (2.00000000000000, 0.000000000000000, 1.73205080756888)

bars:
  0: (0, 5)
  1: (1, 5)
  2: (2, 5)
  3: (3, 6)
  4: (4, 6)
  5: (5, 6)

supports:
  0
  1
  2
  3
  4

characteristics:
  0: (2.00000000000000e8, 0.00250000000000000)

elements <-> chars:
  0: 0
  1: 0
  2: 0
  3: 0
  4: 0
  5: 0

loads:
  6: (-100.000000000000, 0.000000000000000, 0.000000000000000)

supps_d = support3d_dict (supports)
print_dict (supps_d)

  0: (0, 0, 0)
  1: (0, 0, 0)
  2: (0, 0, 0)
  3: (0, 0, 0)
  4: (0, 0, 0)

  0: (0, 0, 0)
  1: (0, 0, 0)
  2: (0, 0, 0)
  3: (0, 0, 0)
  4: (0, 0, 0)

plot_truss3d (joints, bars, supports, loads, force_scale = 0.00825, pull = False)

Click on 3-D image to make it live. Right-click on live image for a control menu.

u, f = DiM (joints, bars, supports, chars, echs, loads)

print_indexed (f.rows())

  0: (-100.000000000000, 100.000000000000)
  1: (200.000000000000, -200.000000000000)
  2: (-100.000000000000, 100.000000000000)
  3: (0.000000000000000, 0.000000000000000)
  4: (0.000000000000000, 0.000000000000000)
  5: (100.000000000000, -100.000000000000)

  0: (-100.000000000000, 100.000000000000)
  1: (200.000000000000, -200.000000000000)
  2: (-100.000000000000, 100.000000000000)
  3: (0.000000000000000, 0.000000000000000)
  4: (0.000000000000000, 0.000000000000000)
  5: (100.000000000000, -100.000000000000)

loads2 = { 5 : (0., -100., 0.) }
print_dict (loads2)

  5: (0.000000000000000, -100.000000000000, 0.000000000000000)

  5: (0.000000000000000, -100.000000000000, 0.000000000000000)

plot_truss3d (joints, bars, supports, loads2, force_scale = 0.00825)

u2, f2 = DiM (joints, bars, supports, chars, echs, loads2)

print_indexed (f2.apply_map (chop()).rows())

  0: (100.000000000000, -100.000000000000)
  1: (0.000000000000000, 0.000000000000000)
  2: (-100.000000000000, 100.000000000000)
  3: (0.000000000000000, 0.000000000000000)
  4: (0.000000000000000, 0.000000000000000)
  5: (0.000000000000000, 0.000000000000000)

  0: (100.000000000000, -100.000000000000)
  1: (0.000000000000000, 0.000000000000000)
  2: (-100.000000000000, 100.000000000000)
  3: (0.000000000000000, 0.000000000000000)
  4: (0.000000000000000, 0.000000000000000)
  5: (0.000000000000000, 0.000000000000000)

loads3 = { 6 : (100., 0., 0.), 5 : (0., -100., 0.) }
print_dict (loads3)

  5: (0.000000000000000, -100.000000000000, 0.000000000000000)
  6: (100.000000000000, 0.000000000000000, 0.000000000000000)

  5: (0.000000000000000, -100.000000000000, 0.000000000000000)
  6: (100.000000000000, 0.000000000000000, 0.000000000000000)

supps_d = support3d_dict (supports)
print_dict (supps_d)

  0: (0, 0, 0)
  1: (0, 0, 0)
  2: (0, 0, 0)
  3: (0, 0, 0)
  4: (0, 0, 0)

  0: (0, 0, 0)
  1: (0, 0, 0)
  2: (0, 0, 0)
  3: (0, 0, 0)
  4: (0, 0, 0)

lpd3 = { 5 : False, 6 : True }

plot_truss3d (joints, bars, supports, loads3, force_scale = 0.00825, pull = lpd3)

u3, f3 = DiM (joints, bars, supps_d, chars, echs, loads3)

print_indexed (f3.rows())

  0: (200.000000000000, -200.000000000000)
  1: (-200.000000000000, 200.000000000000)
  2: (-2.71050543121376e-14, 2.71050543121376e-14)
  3: (0.000000000000000, 0.000000000000000)
  4: (0.000000000000000, 0.000000000000000)
  5: (-100.000000000000, 100.000000000000)

  0: (200.000000000000, -200.000000000000)
  1: (-200.000000000000, 200.000000000000)
  2: (-2.71050543121376e-14, 2.71050543121376e-14)
  3: (0.000000000000000, 0.000000000000000)
  4: (0.000000000000000, 0.000000000000000)
  5: (-100.000000000000, 100.000000000000)

loads4 = { 5 : (-100., 100., 100. * aa), 6 : (100., 0., -100. * aa) }
print_dict (loads4)

  5: (-100.000000000000, 100.000000000000, 173.205080756888)
  6: (100.000000000000, 0.000000000000000, -173.205080756888)

  5: (-100.000000000000, 100.000000000000, 173.205080756888)
  6: (100.000000000000, 0.000000000000000, -173.205080756888)

plot_truss3d (joints, bars, supports, loads4, force_scale = 0.00425)

plot_truss3d (joints, bars, supports, loads4, force_scale = 0.00425, pull = False)

lpd4 = { 5 : True, 6 : False }

plot_truss3d (joints, bars, supports, loads4, force_scale = 0.00425, pull = lpd4)

u4, f4 = DiM (joints, bars, supps_d, chars, echs, loads4)

print_indexed (f4.apply_map (chop()).rows())

  0: (0.000000000000000, 0.000000000000000)
  1: (0.000000000000000, 0.000000000000000)
  2: (200.000000000000, -200.000000000000)
  3: (-100.000000000000, 100.000000000000)
  4: (-100.000000000000, 100.000000000000)
  5: (-100.000000000000, 100.000000000000)

  0: (0.000000000000000, 0.000000000000000)
  1: (0.000000000000000, 0.000000000000000)
  2: (200.000000000000, -200.000000000000)
  3: (-100.000000000000, 100.000000000000)
  4: (-100.000000000000, 100.000000000000)
  5: (-100.000000000000, 100.000000000000)

PPPK_program_DiM_2021_(v0.2.1)_(Truss)

255 days ago by fresl