NS_Truss_1_(v0.6) -- Sage

%auto

def map2 (f, l2) :
    return map (lambda y : map (f, y), l2)

def map_dict (f, d) :
    return dict (zip (d.keys(), map (f, d.values())))
    
def map2_dict (f, d) :
    return dict (zip (d.keys(), map2 (f, d.values())))

def filter_dict (f, d) : 
    fd = {}
    for k in d.keys():
        if f (d[k]) : 
            fd[k] = d[k]
    return fd
    

def ND (sig_digits_10) : 
    return lambda x : N (x, digits = sig_digits_10)
    
def chop (tolerance = None) :
    def f_ (x) :
        t = sqrt ((parent(x)(1)).ulp()) if tolerance == None else tolerance
        return x if (abs(x) >= t) else parent(x)(0)
    return f_
 
def NDC (sig_digits_10, tolerance = None) : 
    return lambda x : N (chop(tolerance) (x), digits = sig_digits_10)
   
def sqrt_ulp (field = RDF) :
    return sqrt (field(1).ulp())
    

def dict_joint_bars (free_nodes, elems) : 
    nds_els = {}
    for fni in free_nodes :
        nds_els[fni] = []
    for (ei, ec) in elems.iteritems() : 
        (i0, i1) = ec
        if i0 in free_nodes :
            nds_els[i0].append (ei)
        if i1 in free_nodes :
            nds_els[i1].append (ei)
    return nds_els


def dict_equation_index (free_nodes, dimen = 3) : 
    return dict (zip (sorted (free_nodes), xrange (0, dimen * len (free_nodes), dimen)))  
      
def dict_index_equation (free_nodes, dimen = 3) : 
    return dict (zip (xrange (0, dimen * len (free_nodes), dimen), sorted (free_nodes)))  
      
def dict_unknown_index (elems) : 
    return dict (zip (sorted (elems.keys()), xrange (len (elems))))   

def dict_index_unknown (elems) : 
    return dict (zip (xrange (len (elems)), sorted (elems.keys())))   

def dict_label_index (lbls) : 
    return dict (zip (sorted (lbls), xrange (len (lbls))))   

def dict_index_label (lbls) : 
    return dict (zip (xrange (len (lbls)), sorted (lbls)))   

        
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 all_indices (l) :
    return xrange (len (l))
    

def maxwells_rule (joints, supports, bars) : 
    dim = len (joints.values()[0]) 
    free_joints = others (supports, joints)
    mxw = dim*len (free_joints) - len (bars)
    print dim, '*', len (free_joints), '-', len (bars), ' == ', mxw, 
    if mxw != 0 :
        print ' !=  0'
    else :
        print 
    return 


def base_vectors (dim, field = RDF) :
    bv = zero_vector (field, dim)
    bvs = [copy (bv) for i in xrange (dim)]
    for i in xrange (dim) :
        bvs[i][i] = field (1)
    return bvs


def bar_length (xxi, xxj) : 
    return sqrt (sum (map (lambda xi, xj : (xj - xi)^2, xxi, xxj)))


def equilibrium_matrix (nodes, elems, free_nodes, field = RDF) : 
    nfj = len (free_nodes)
    nb = len (elems)
    dim = len (nodes.values()[0]) 
    ne_inc = dict_joint_bars (free_nodes, elems)
    ukws = dict_unknown_index (elems)
    eqs = dict_equation_index (free_nodes, dim)
    A = zero_matrix (field, dim*nfj, nb)
    
    for ni in free_nodes : 
        ei = eqs[ni]
        bni = ne_inc[ni]
        for bi in bni :
            if ni == elems[bi][0] :
                nj = elems[bi][1]
            else :
                nj = elems[bi][0]
            xxi = map (field, nodes[ni]) 
            xxj = map (field, nodes[nj]) 
            # if field == SR :
            #     if ni < nj :
            #         l = var ('l_%d%d' % (ni, nj))
            #     else :
            #         l = var ('l_%d%d' % (nj, ni))
            # else :
            l = sqrt (sum (map (lambda xi, xj : (xj - xi)^2, xxi, xxj)))
            cc = map (lambda xi, xj : (xj - xi) / l, xxi, xxj)
            ui = ukws[bi]
            A[ei : (ei+dim), ui] += matrix (dim, 1, cc)
            
    return A


def load_vector (loads, free_nodes, field = RDF) :
    dim = len (loads.values()[0]) 
    nfj = len (free_nodes)
    eqs = dict_equation_index (free_nodes, dim)
    b = zero_vector (field, dim*nfj)
    for li, ff in loads.iteritems() :
        ei = eqs[li]
        b[ei : (ei+dim)] = map (field, ff)
    return b


def augmented_equil_matrix (nodes, elems, loads, free_nodes, field = RDF) : 
    nfj = len (free_nodes)
    nb = len (elems)
    dim = len (nodes.values()[0]) 
    ne_inc = dict_joint_bars (free_nodes, elems)
    ukws = dict_unknown_index (elems)
    eqs = dict_equation_index (free_nodes, dim)
    A = zero_matrix (field, dim*nfj, nb + 1)
    
    for ni in free_nodes : 
        ei = eqs[ni]
        bni = ne_inc[ni]
        for bi in bni :
            if ni == elems[bi][0] :
                nj = elems[bi][1]
            else :
                nj = elems[bi][0]
            xxi = map (field, nodes[ni]) 
            xxj = map (field, nodes[nj]) 
            l = sqrt (sum (map (lambda xi, xj : (xj - xi)^2, xxi, xxj)))
            cc = map (lambda xi, xj : (xj - xi) / l, xxi, xxj)
            ui = ukws[bi]
            A[ei : (ei+dim), ui] += matrix (dim, 1, cc)
    
    for li, ff in loads.iteritems() :
        ei = eqs[li]
        ej = nb
        A[ei : (ei+dim), ej] = matrix (dim, 1, map (lambda x : -field(x), ff))
        
    A.subdivide ([], [nb])
    return A


def dict_bar_force (elems, forces) :
    ukws = dict_label_index (elems.keys()) 
    return dict ([(b, forces[i]) for b, i in ukws.iteritems()]) 


def truss_skd (nodes, elems, supports, loads, field = RDF) : 
    free_nodes = others (supports, nodes)
    A = equilibrium_matrix (nodes, elems, free_nodes, field)
    b = load_vector (loads, free_nodes, field)
    ff = A.solve_right (-b)
    return dict_bar_force (elems, ff)
    

def is_zz (x) :
    return x == 0

def is_zz_tol (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 (m, is_zero_f = is_zz_tol(), in_place = False) :
    pivot = 0
    if in_place == True : 
        r = m
    else :
        r = copy (m)
    pivotlist = []
    rnr = r.nrows()
    for i in xrange (r.ncols()) :
        if pivot == rnr :
            break
        maxrow = pivot
        for k in xrange (pivot+1, r.nrows()):
            if abs (r[k][i]) > abs (r[maxrow][i]) : 
                maxrow = k
        if pivot != maxrow :
            r.swap_rows (pivot, maxrow)
        if not is_zero_f (r[pivot, i]) :
            scale = r[pivot, i]
            r.rescale_row (pivot, 1/scale, i)
            for j in xrange (pivot+1, r.nrows()) :
                scale = r[j, i]
                r.add_multiple_of_row (j, pivot, -scale, i)
            pivotlist.append (i)
            pivot += 1
    return r, pivotlist

def ref_solution (Ar, piv) : 
    v = vector (Ar.base_ring(), Ar.ncols() - 1)
    np = len (piv)
    for j in xrange (np) :
        v[piv[j]] = Ar[j,-1]
    for j in xrange (np-1, -1, -1) :
        smj = 0 
        for jj in xrange (np-1, j, -1) :
            smj += Ar[j, piv[jj]] * v[piv[jj]]
        v[piv[j]] -= smj
    return v


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

def rref_solution (Arr, piv) :
    v = vector (Arr.base_ring(), Arr.ncols() - 1)
    np = len (piv)
    for j in xrange (np) :
        v[piv[j]] = Arr[j,-1]
    return v
    
 
def truss (nodes, elems, supports, loads, is_zero_f = is_zz_tol(), field = RDF) : 
    free_nodes = others (supports, nodes)
    Ab = augmented_equil_matrix (nodes, elems, loads, free_nodes, field)
    Ab, piv = rref_wpp (Ab, is_zero_f = is_zero_f, in_place = True)
    return dict_bar_force (elems, rref_solution (Ab, piv))

    
def pivot_columns (mr) : 
    nr = mr.nrows()
    pc = []
    for i in xrange (mr.nrows()) : 
        for j in xrange (i, mr.ncols()) :
            if mr[i, j] != 0 :
                pc.append (j)
                break
    return pc


def ref_ext_wpp (m, b = None, is_zero_f = is_zz_tol()) :
    pivot = 0
    if b == None :
        r = block_matrix (1, 2, [m, identity_matrix (m.base_ring(), m.nrows())])
    else :
        r = block_matrix (1, 3, [m, identity_matrix (m.base_ring(), m.nrows()), b])
    mncols = m.ncols()
    pivotlist = []
    for i in xrange (mncols) :
        if pivot == r.nrows() :
            break
        maxrow = pivot
        for k in xrange (pivot+1, r.nrows()):
            if abs (r[k][i]) > abs (r[maxrow][i]) : 
                maxrow = k
        if pivot != maxrow :
            r.swap_rows (pivot, maxrow)
        if not is_zero_f (r[pivot, i]) :
            scale = r[pivot, i]
            r.rescale_row (pivot, 1/scale, i)
            for j in xrange (pivot+1, r.nrows()) :
                scale = r[j, i]
                r.add_multiple_of_row (j, pivot, -scale, i)
            pivotlist.append (i)
            pivot += 1
    return r, pivotlist


def rref_ext_wpp (m, b = None, is_zero_f = is_zz_tol()) :
    pivot = 0
    if b == None :
        r = block_matrix (1, 2, [m, identity_matrix (m.base_ring(), m.nrows())])
    else :
        r = block_matrix (1, 3, [m, identity_matrix (m.base_ring(), m.nrows()), b])
    mncols = m.ncols()
    pivotlist = []
    for i in xrange (mncols) :
        if pivot == r.nrows() :
            break
        maxrow = pivot
        for k in xrange (pivot+1, r.nrows()):
            if abs (r[k][i]) > abs (r[maxrow][i]) : 
                maxrow = k
        if pivot != maxrow :
            r.swap_rows (pivot, maxrow)
        if not is_zero_f (r[pivot, i]) :
            scale = r[pivot, i]
            r.rescale_row (pivot, 1/scale, i)
            for j in xrange (r.nrows()) :
                if j == pivot :
                    continue
                scale = r[j, i]
                r.add_multiple_of_row (j, pivot, -scale, i)
            pivotlist.append (i)
            pivot += 1
    return r, pivotlist


def column_space (m, piv, is_zero_f = is_zz_tol(), rref = False) :
    cs = []
    for p in piv :
        cs.append (m.column (p))
    if rref == False : 
        return cs
    mcs = matrix (cs)
    mcsr, pp = rref_wpp (mcs, is_zero_f = is_zero_f, in_place = True)
    return mcsr.rows()
    
def row_space (mr, piv) :
    rs = []
    m = mr 
    for i in xrange (len (piv)) :
        rs.append (m.row (i))
    return rs

def left_kernel (mre, piv) :
    lk = []
    ii = mre 
    nr = ii.nrows()
    r0 = len (piv)
    for i in xrange (r0, nr) :
        lk.append (ii.row (i))
    return lk

def subs_r_ (m, nc, piv, np, i) : 
    v = vector (m.base_ring(), nc)
    v[i] = 1
    for j in xrange (np) :
        v[piv[j]] = -m[j,i]
    return v
    
def backsubs_ (m, nc, piv, np, i) : 
    v = vector (m.base_ring(), nc)
    v[i] = 1
    for j in xrange (np-1, -1, -1) :
        smj = 0 
        for jj in xrange (np-1, j, -1) :
            smj += m[j, piv[jj]] * v[piv[jj]]
        v[piv[j]] = -m[j,i] - smj
    return v

def kernel (mr, piv, backsubs = backsubs_) :
    k = []
    m = mr 
    nc = m.ncols()
    np = len (piv)
    piv_s = set (piv)
    for i in xrange (nc) :
        if not i in piv_s : 
            k.append (backsubs (m, nc, piv, np, i))
    return k

def subspace_bases (m, is_zero_f = is_zz_tol()) :
    mre, piv = ref_ext_wpp (m, is_zero_f = is_zero_f)
    rs = row_space (mre.subdivision (0, 0), piv)
    cs = column_space (m, piv, is_zero_f = is_zero_f)
    rk = kernel (mre.subdivision (0, 0), piv)
    lk = left_kernel (mre.subdivision (0, 1), piv)
    return rs, cs, rk, lk

def subspace_bases_r (m, is_zero_f = is_zz_tol()) :
    mre, piv = rref_ext_wpp (m, is_zero_f = is_zero_f)
    rs = row_space (mre.subdivision (0, 0), piv)
    cs = column_space (m, piv, is_zero_f = is_zero_f)
    rk = kernel (mre.subdivision (0, 0), piv)
    lk = left_kernel (mre.subdivision (0, 1), piv)
    return rs, cs, rk, lk

def subspaces (m, is_zero_f = is_zz_tol(), return_pivots = False) :
    mre, piv = ref_ext_wpp (m, is_zero_f = is_zero_f)
    nj3 = m.nrows()
    nb = m.ncols()
    rm = len (piv)
    rsb = row_space (mre.subdivision (0, 0), piv)
    rsn = nb
    rsd = rm
    rkb = kernel (mre.subdivision (0, 0), piv)
    rkn = nb
    rkd = nb - rm
    csb = column_space (m, piv)
    csn = nj3
    csd = rm
    lkb = left_kernel (mre.subdivision (0, 1), piv)
    lkn = nj3
    lkd = nj3 - rm
    rs = {'name' : 'Row space', 'degree' : rsn, 'dimension' : rsd, 'basis' : rsb} 
    cs = {'name' : 'Column space', 'degree' : csn, 'dimension' : csd, 'basis' : csb} 
    rk = {'name' : 'Kernel', 'degree' : rkn, 'dimension' : rkd, 'basis' : rkb} 
    lk = {'name' : 'Left kernel', 'degree' : lkn, 'dimension' : lkd, 'basis' : lkb} 
    if return_pivots == True :
        r = rs, cs, rk, lk, piv
    else :
        r = rs, cs, rk, lk
    return r

def subspaces_r (m, is_zero_f = is_zz_tol(), return_pivots = False) :
    mre, piv = rref_ext_wpp (m, is_zero_f = is_zero_f)
    nj3 = m.nrows()
    nb = m.ncols()
    rm = len (piv)
    rsb = row_space (mre.subdivision (0, 0), piv)
    rsn = nb
    rsd = rm
    rkb = kernel (mre.subdivision (0, 0), piv)
    rkn = nb
    rkd = nb - rm
    csb = column_space (m, piv, rref = True, is_zero_f = is_zero_f)
    csn = nj3
    csd = rm
    lkb = left_kernel (mre.subdivision (0, 1), piv)
    lkn = nj3
    lkd = nj3 - rm
    rs = {'name' : 'Row space', 'degree' : rsn, 'dimension' : rsd, 'basis' : rsb} 
    cs = {'name' : 'Column space', 'degree' : csn, 'dimension' : csd, 'basis' : csb} 
    rk = {'name' : 'Kernel', 'degree' : rkn, 'dimension' : rkd, 'basis' : rkb} 
    lk = {'name' : 'Left kernel', 'degree' : lkn, 'dimension' : lkd, 'basis' : lkb} 
    if return_pivots == True :
        r = rs, cs, rk, lk, piv
    else :
        r = rs, cs, rk, lk
    return r

def print_subspace (ss, sb = True) :
    print '{0}{1}'.format (ss['name'], ':')
    print '   Vector space of degree', ss['degree'], 'and dimension', ss['dimension']
    if sb == True : 
        print '   Basis:'
        print_list (ss['basis'], '     ')
        
def are_linearly_dependent (vl, is_zero_f = is_zz_tol()) :
    m = matrix (vl) 
    mre, piv = rref_wpp (m, is_zero_f = is_zero_f)
    return len (vl) != len (piv)

def is_in_subspace (v, sp, is_zero_f = is_zz_tol()) : 
    vl = [v]
    vl.extend (sp)
    return are_linearly_dependent (vl, is_zero_f = is_zero_f)

def has_component_in_subspace (v, sp, is_zero_f = is_zz_tol()) :
    r = False
    for vv in sp :
        s = v * vv
        if not is_zero_f (s) :
            r = True
            break
    return r


def all_bar_forces (primary_bars, redundant_bars, primary_forces) : 
    bfs = vector (primary_forces.base_ring(), len (primary_bars) + len (redundant_bars))
    for i in xrange (len (primary_bars)) :
        bfs[primary_bars[i]] = primary_forces[i]
    return bfs    


def diag_flex (elems, nodes, E, F, field = RDF) :
    ldelta = []
    for el in elems :
        ni, nj = elems[el]
        xi = nodes[ni] 
        xj = nodes[nj] 
        ldelta.append (bar_length (xi, xj) / (E * F))
    return diagonal_matrix (RDF, ldelta)


def diag_stiff (elems, nodes, E, F, field = RDF) :
    lk = []
    for el in elems :
        ni, nj = elems[el]
        xi = nodes[ni] 
        xj = nodes[nj] 
        lk.append (E * F / bar_length (xi, xj))
    return diagonal_matrix (RDF, lk)


def make_bars__ (nb, b0j, steps) :
    return [map (lambda x, y : x + y, b0j, map (lambda x : i * x, steps)) for i in xrange (nb)]


def schwedler_1 (r0, hh, n) :
    if not is_even (n) :
        print 'n must be even'
        return

    jsl = []
    bsl = []
    ssl = []
    rsl = [r0]
    
    h0 = hh[-1]
    rr = r0^2/(2*h0) + h0/2
    aa = rr - h0
    for h in hh[:-1] :
        rsl.append (sqrt (rr^2 - (aa+h)^2))
    hh.insert (0, 0.0)
    alpha = RDF (2*pi/n)
    for i in xrange (len (rsl)) :
        for j in xrange (n) :
            x = rsl[i] * cos (j*alpha)
            y = rsl[i] * sin (j*alpha)
            z = hh[i]
            jsl.append ((x, y, z))
    jsd = dict (zip (xrange (len (jsl)), jsl))
    
    for i in xrange (0, len (hh) - 2) :
        n0 = i*n
        for j in xrange (n) :
            bsl.append ((n0+j, n0+j+n))
        for j in xrange (n) :
            bsl.append ((n0+n+j, n*(i+1)+(n0+j+1)%n))
        for j in xrange (0, n, 2) :
            bsl.append ((n0+j, (n*(i+1)+(n0+j+1)%n)))
        for j in xrange (0, n, 2) :
            bsl.append ((n0+j, (n*(i+1)+(n0+j-1)%n)))
    bsd = dict (zip (xrange (len (bsl)), bsl))
           
    ssl = range (n)
    return jsd, bsd, ssl


def schwedler_2 (r0, hh, n, ccw = True) :        
    jsl = []
    bsl = []
    ssl = []
    rsl = [r0]
    
    h0 = hh[-1]
    rr = r0^2/(2*h0) + h0/2
    aa = rr - h0
    for h in hh[:-1] :
        rsl.append (sqrt (rr^2 - (aa+h)^2))
    hh.insert (0, 0.0)
    alpha = RDF (2*pi/n) 
    for i in xrange (len (rsl)) :
        for j in xrange (n) : 
            x = rsl[i] * cos (j*alpha)
            y = rsl[i] * sin (j*alpha)
            z = hh[i]
            jsl.append ((x, y, z))
    jsd = dict (zip (xrange (len (jsl)), jsl))
    
    for i in xrange (0, len (hh) - 2) : 
        n0 = i*n 
        for j in xrange (n) : 
            bsl.append ((n0+j, n0+j+n)) 
        for j in xrange (n) : 
            bsl.append ((n0+n+j, n*(i+1)+(n0+j+1)%n))
        if ccw == True :
            for j in xrange (0, n) : 
                bsl.append ((n0+j, (n*(i+1)+(n0+j+1)%n)))
        else :
            for j in xrange (0, n) : 
                bsl.append ((n0+j, (n*(i+1)+(n0+j-1)%n)))
    bsd = dict (zip (xrange (len (bsl)), bsl))
           
    ssl = range (n)
    return jsd, bsd, ssl


def schwedler_kind (r0, hh, n) :        
    jsl = []
    bsl = []
    ssl = []
    rsl = [r0]
    
    h0 = hh[-1]
    rr = r0^2/(2*h0) + h0/2
    aa = rr - h0
    for h in hh[:-1] :
        rsl.append (sqrt (rr^2 - (aa+h)^2))
    hh.insert (0, 0.0)
    alpha = RDF (2*pi/n) 
    for i in xrange (len (rsl)) :
        for j in xrange (n) : 
            x = rsl[i] * cos (j*alpha)
            y = rsl[i] * sin (j*alpha)
            z = hh[i]
            jsl.append ((x, y, z))
    jsd = dict (zip (xrange (len (jsl)), jsl))
    
    for i in xrange (0, len (hh) - 2) : 
        n0 = i*n 
        for j in xrange (n) : 
            bsl.append ((n0+j, n0+j+n)) 
        for j in xrange (n) : 
            bsl.append ((n0+n+j, n*(i+1)+(n0+j+1)%n))
    bsd = dict (zip (xrange (len (bsl)), bsl))
           
    ssl = range (n)
    return jsd, bsd, ssl


def schwedler_sind (r0, hh, n) :        
    jsl = []
    bsl = []
    ssl = []
    rsl = [r0]
    
    h0 = hh[-1]
    rr = r0^2/(2*h0) + h0/2
    aa = rr - h0
    for h in hh[:-1] :
        rsl.append (sqrt (rr^2 - (aa+h)^2))
    hh.insert (0, 0.0)
    alpha = RDF (2*pi/n) 
    for i in xrange (len (rsl)) :
        for j in xrange (n) : 
            x = rsl[i] * cos (j*alpha)
            y = rsl[i] * sin (j*alpha)
            z = hh[i]
            jsl.append ((x, y, z))
    jsd = dict (zip (xrange (len (jsl)), jsl))
    
    for i in xrange (0, len (hh) - 2) : 
        n0 = i*n 
        for j in xrange (n) : 
            bsl.append ((n0+j, n0+j+n)) 
        for j in xrange (n) : 
            bsl.append ((n0+n+j, n*(i+1)+(n0+j+1)%n))
        for j in xrange (0, n) : 
            bsl.append ((n0+j, (n*(i+1)+(n0+j+1)%n)))
        for j in xrange (0, n) : 
            bsl.append ((n0+j, (n*(i+1)+(n0+j-1)%n)))
    bsd = dict (zip (xrange (len (bsl)), bsl))
           
    ssl = range (n)
    return jsd, bsd, ssl


def hypar_xy (x, y, z_scale) :
    return x * y * z_scale 
    
def hypar_x2y2 (x, y, z_scale) :
    return (x^2 - y^2) * z_scale 
    
def surface_net (f, lx, ly, m, n) : 
    jsl = []
    bsl = []
    ssl = []

    m2 = m + 2
    n2 = n + 2
    xmin = -lx/2
    ymin = -ly/2
    hx = lx/(m2-1)
    hy = ly/(n2-1) 
    
    for j in xrange (n2) : 
        for i in xrange (m2) : 
            xx = xmin + i*hx
            yy = ymin + j*hy
            jsl.append ((xx, yy, f (xx, yy)))
    jsd = dict (zip (xrange (len (jsl)), jsl))

    for j in xrange (1, n2-1) : 
        n0 = j*m2 
        for i in xrange (m2-1) : 
            bsl.append ((n0+i, n0+i+1)) 
    for i in xrange (1, m2-1) : 
        for j in xrange (0, (n2-1)*m2, m2) : 
            bsl.append ((i+j, i+j+m2)) 
    bsd = dict (zip (xrange (len (bsl)), bsl))
        
    ssl = range (m2)
    for i in xrange (1, n2-1) :
        ssl.append (i*m2)
        ssl.append (i*m2 + m2 - 1)
    ssl.extend (range ((n2-1)*m2, n2*m2))

    return jsd, bsd, ssl


def make_loads__ (nl, n0, step, ld) :
    return zip (range (n0, n0 + nl*step, step), [ld for i in xrange (nl)])
    
def make_loads (generator) :
    ll = []
    for gg in generator :
        if len (gg) == 3 : 
            nl, n0, ld = gg
            step = 1
        else :
            nl, n0, step, ld = gg
        ll.extend (make_loads__ (nl, n0, step, ld))
    return dict (ll)
    
    
def tension_compression (bar_forces, tolerance = 1.5e-8) :
    tt = [] 
    cc = []
    zz = []
    for (bi, f) in bar_forces.iteritems() : 
        if abs (f) > tolerance : 
            if f > 0 :
                tt.append ((bi, f))
            else :
                cc.append ((bi, f))
        else :
            zz.append (bi)
    tt = sorted (tt, key = lambda x : x[1])
    cc = sorted (cc, key = lambda x : x[1])
    return (tt, cc, zz)

def unique_lp (l, tolerance = 1e-8) : 
    li = []
    if l != []: 
        (ai, av) = l[0]
        ll = [ai]
        for (bi, bv) in l[1:] :
            if abs (av - bv) < tolerance :
                ll.append (bi)
            else : 
                li.append (copy (ll))
                ll = [bi] 
                av = bv
        li.append (copy (ll))
    return li   


def plot_free_joints2d (pts, rr, color = 'black') :
    return sum ([circle (pt, rr, edgecolor = color, facecolor = 'white',
                 fill = True, aspect_ratio = 1) for pt in pts])

def plot_supports2d (pts, rr) :
    return sum ([circle (pt, rr, edgecolor = 'black', facecolor = 'black',
                 fill = True, aspect_ratio = 1) for pt in pts])

def plot_bars2d (bars, joints, color = 'black', thickness = 2) :
    return sum ([line2d ((joints[b[0]], joints[b[1]]), thickness = thickness,
                 color = color, aspect_ratio = 1) for b in bars])

def plot_loads2d (loads, joints, load_scale) :
    pls = []
    for li in loads :
        tail = joints[li]
        ff = loads[li]
        head = [tail[0] + load_scale*ff[0], tail[1] + load_scale*ff[1]]
        pls.append (arrow2d (tail, head, width = 3, color = 'crimson', aspect_ratio = 1))
    return sum (pls)
                                                 
def plot_truss2d (joints, bars, supports, loads = [], load_scale = 0.01) :
    prr = point2d (joints.values(), aspect_ratio = 1)
    rrx = (prr.xmax() - prr.xmin())
    rry = (prr.ymax() - prr.ymin())
    rr = rrx if rrx > rry else rry*3/2
    pfjs = plot_free_joints2d (list_indexed (others (supports, joints), joints), rr/120)
    psjs = plot_supports2d (list_indexed (supports, joints), rr/150)
    pbs = plot_bars2d (bars.values(), joints)
    pt = pfjs + psjs + pbs
    if loads != [] :
        pls = plot_loads2d (loads, joints, load_scale)
        pt += pls
    return pt
                 
def plot_truss_with_forces2d (joints, bars, supports, forces) :
    prr = point2d (joints.values(), aspect_ratio = 1)
    rrx = (prr.xmax() - prr.xmin())
    rry = (prr.ymax() - prr.ymin())
    rr = rrx if rrx > rry else rry*3/2
    pfjs = plot_free_joints2d (list_indexed (others (supports, joints), joints), rr/120)
    psjs = plot_supports2d (list_indexed (supports, joints), rr/150)
    pt = pfjs + psjs
    tt, cc, zi = tension_compression (forces)
    ti = unique_lp (tt)
    ti.reverse()
    ci = unique_lp (cc)
    if zi != [] :
        cl = Color ('orange')
        cl = cl.lighter (0.425)
        pbz = plot_bars2d (list_indexed (zi, bars), joints, color = cl, thickness = 3)
        pt += pbz
    if ti != [] :
        step = 1./len (ti)
        cl = Color ('darkblue')
        pbt = plot_bars2d (list_indexed (ti[0], bars), joints, color = cl, thickness = 3)
        for i in xrange (1, len (ti)) :
            cl = cl.lighter (step)
            pbt += plot_bars2d (list_indexed (ti[i], bars), joints, color = cl, thickness = 3)
        pt += pbt
    if ci != [] :
        step = 1./len (ci)
        cl = Color ('red')
        pbt = plot_bars2d (list_indexed (ci[0], bars), joints, color = cl, thickness = 3)
        for i in xrange (1, len (ci)) :
            cl = cl.lighter (step)
            pbt += plot_bars2d (list_indexed (ci[i], bars), joints, color = cl, thickness = 3)
        pt += pbt
    return pt

def plot_truss_with_displacements2d (joints, bars, supports, displs, displ_scale = 1) : 
    prr = point2d (joints.values(), aspect_ratio = 1)
    rrx = (prr.xmax() - prr.xmin())
    rry = (prr.ymax() - prr.ymin())
    rr = rrx if rrx > rry else rry*3/2
    fjis = others (supports, joints)
    fjs = list_indexed (fjis, joints)
    sjs = list_indexed (supports, joints)
    pfjs = plot_free_joints2d (fjs, rr/180, 'grey')
    psjs = plot_supports2d (sjs, rr/150)
    pbs = plot_bars2d (bars.values(), joints, color = 'grey', thickness = 1)
    pt = pfjs + psjs + pbs
    ds = displ_scale
    dfjs = copy (fjs)
    dispc = map (chop(), displs)
    for i in xrange (len (dfjs)) :
        jnt = dfjs[i]
        dfjs[i] = (jnt[0] + ds*dispc[i*2], jnt[1] + ds*dispc[i*2 + 1])
    djs = copy (joints)
    for i in xrange (len (fjis)) :
        djs[fjis[i]] = dfjs[i] 
    pdfjs = plot_free_joints2d (dfjs, rr/120)
    pdbs = plot_bars2d (bars.values(), djs)
    pt = pbs + pfjs + pdbs + pdfjs + psjs
    return pt


def plot_free_joints3d (pts, color = 'white', size = 27) : 
    return point3d (pts, color = color, size = size, aspect_ratio = 1, frame = False)
  
def plot_supports3d (pts, color = 'black', size = 23) : 
    return point3d (pts, color = color, size = size, aspect_ratio = 1, frame = False)
  
def plot_bars3d (bars, joints, color = 'black', thickness = 7) : 
    return sum ([line3d ((joints[b[0]], joints[b[1]]), thickness = thickness,
                 color = color, aspect_ratio = 1, frame = False) for b in bars]) 


def plot_arrow3d (ff, tail = (0, 0, 0), scale = 1, color = 'blue', thickness = 1) :
    fs = scale * ff
    nf = norm (fs)
    if nf != 0 :
        fw = (3.425 * thickness * nf)^(2/5)
        # fw = sqrt (0.825 * radius * nf)
        # fw = radius * 3.25
        rd = fw / 50
        hl = 6.425 * rd
        hr = 1.625 * rd
        vtl = vector (tail)
        ar = arrow3d (vtl, vtl + fs, radius = rd, head_len = hl, head_radius = hr, color = color, aspect_ratio = 1, frame = False)
        # ar = arrow3d (vtl, vtl + fs, width = fw, head_len = hl, head_radius = hr, color = color)
    else :
        ar = Sphere (0.001, color = 'white').translate (tail)
    return ar

def plot_loads3d (loads, joints, load_scale, pull = True, thickness = 1, color = 'crimson') : 
    pls = []
    for li in loads :
        ff = vector (loads[li]) 
        if pull == True :
            tail = joints[li]
        else : 
            head = joints[li]
            tail = [head[0] - load_scale*ff[0], head[1] - load_scale*ff[1], head[2] - load_scale*ff[2]]
        pls.append (plot_arrow3d (ff, tail, color = color, thickness = thickness, scale = load_scale))
    return sum (pls)

    
def plot_truss3d (joints, bars, supports, loads = [], load_scale = 0.01, pull = True) : 
    pfjs = plot_free_joints3d (list_indexed (others (supports, joints), joints))
    psjs = plot_supports3d (list_indexed (supports, joints))
    pbs = plot_bars3d (bars.values(), joints)
    pt = pfjs + psjs + pbs
    if loads != [] :
        pls = plot_loads3d (loads, joints, load_scale, pull)
        pt += pls
    return pt

   
def plot_truss_with_forces3d (joints, bars, supports, forces) : 
    pfjs = plot_free_joints3d (list_indexed (others (supports, joints), joints))
    psjs = plot_supports3d (list_indexed (supports, joints))
    pt = pfjs + psjs 
    tt, cc, zi = tension_compression (forces)
    ti = unique_lp (tt)
    ti.reverse() 
    ci = unique_lp (cc)
    if zi != [] :
        cl = Color ('orange')
        cl = cl.lighter (0.625)
        pbz = plot_bars3d (list_indexed (zi, bars), joints, color = cl, thickness = 7) 
        pt += pbz
    if ti != [] :
        step = 1./len (ti) 
        cl = Color ('indigo')
        pbt = plot_bars3d (list_indexed (ti[0], bars), joints, color = cl, thickness = 7) 
        for i in xrange (1, len (ti)) :
            cl = cl.lighter (step)
            pbt += plot_bars3d (list_indexed (ti[i], bars), joints, color = cl, thickness = 7) 
        pt += pbt
    if ci != [] :
        step = 1./len (ci) 
        cl = Color ('red')
        pbt = plot_bars3d (list_indexed (ci[0], bars), joints, color = cl, thickness = 7) 
        for i in xrange (1, len (ci)) :
            cl = cl.lighter (step)
            pbt += plot_bars3d (list_indexed (ci[i], bars), joints, color = cl, thickness = 7) 
        pt += pbt
    return pt


def plot_disp_lines3d (joints, djoints, color = 'orange', thickness = 3) : 
    return sum ([line3d ((joints[i], djoints[i]), arrow_head = True, thickness = thickness,
                 color = color, aspect_ratio = 1, frame = False) for i in xrange (len (joints))]) 
                         
def plot_truss_with_displacements3d (joints, bars, supports, displs, 
                                     displ_scale = 1, plot_dlines = True) : 
    fjis = others (supports, joints)
    fjs = list_indexed (fjis, joints)
    pfjs = plot_free_joints3d (fjs, color = 'whitesmoke', size = 18)
    psjs = plot_supports3d (list_indexed (supports, joints))
    pbs = plot_bars3d (bars.values(), joints, color = 'whitesmoke', thickness = 4)
    ds = displ_scale
    dfjs = copy (fjs)
    dispc = map (chop(), displs)
    for i in xrange (len (dfjs)) :
        jnt = dfjs[i]
        dfjs[i] = (jnt[0] + ds*dispc[i*3], jnt[1] + ds*dispc[i*3 + 1], jnt[2] + ds*dispc[i*3 + 2])
    djs = copy (joints)
    for i in xrange (len (fjis)) :
        djs[fjis[i]] = dfjs[i] 
    pdfjs = plot_free_joints3d (dfjs)
    pdbs = plot_bars3d (bars.values(), djs)
    pt = pfjs + psjs + pbs + pdfjs + pdbs
    if plot_dlines == True :
        pdl = plot_disp_lines3d (fjs, dfjs)
        pt += pdl
    return pt


def plot_vectors3d (vs, color = 'green', thickness = 0.1) : 
    pvs = []
    coord_axes = []
    coord_axes.append (plot_arrow3d (vector ((1, 0, 0)), (0, 0, 0), thickness = 0.01))
    coord_axes.append (plot_arrow3d (vector ((0, 1, 0)), (0, 0, 0), thickness = 0.01))
    coord_axes.append (plot_arrow3d (vector ((0, 0, 1)), (0, 0, 0), thickness = 0.01))
    for v in vs :
        pvs.append (plot_arrow3d (v, (0, 0, 0), color = color, thickness = thickness))
    return sum (coord_axes) + sum (pvs)
    
    
def plot_vectors2d (vs, color = 'green', thickness = 1) :
    pvs = [] 
    for v in vs :
        pvs.append (v.plot (color = color, thickness = thickness, aspect_ratio = 1))
    return sum (pvs)


def print_indexed (l, i0 = 0, idx_len = 4, indent = '') :
    fstr = '{0}{1:<' + str (idx_len) + '} {2}'
    for i, e in enumerate (l) :
        print fstr.format (indent, str(i+i0) + ':', e)

def print_list (l, indent = '', braces = False) : 
    if l == [] :
        print '{0}{1}'.format (indent, '[]')
    else :
        if braces == True :
            print '{0}{1}'.format (indent, '[')
        for x in l :
            print '{0}{1}'.format (indent, x)
        if braces == True :
            print '{0}{1}'.format (indent, ']')
        
def print_dict (d, key_len = 4, sort_keys = True, indent = '') : 
    fstr = '{0}{1:<' + str (key_len) + '} {2}'
    for k in sorted (d.keys()) if sort_keys == True else d.keys():
        print fstr.format (indent, str(k) + ':', d[k])

joints = { 0: (0.0, 0.0, 0.0), 1: (4.0, -1.5, 0.0), 2: (2.0, 3.5, 0.0), 3: (1.5, 2.0, 3.0) }
print 'joints:'
print_dict (joints, indent = '   ')
print

bars = { 0: (0, 3), 1: (1, 3), 2: (2, 3) }
print 'bars:' 
print_dict (bars, indent = '   ')
print

supports = [0, 1, 2]
print 'supports:', supports
print

loads = { 3: (75.0, 100.0, -25.0) }
print 'loads:' 
print_dict (loads, indent = '   ')

joints:
   0:   (0.000000000000000, 0.000000000000000, 0.000000000000000)
   1:   (4.00000000000000, -1.50000000000000, 0.000000000000000)
   2:   (2.00000000000000, 3.50000000000000, 0.000000000000000)
   3:   (1.50000000000000, 2.00000000000000, 3.00000000000000)

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

supports: [0, 1, 2]

loads:
   3:   (75.0000000000000, 100.000000000000, -25.0000000000000)

joints:
   0:   (0.000000000000000, 0.000000000000000, 0.000000000000000)
   1:   (4.00000000000000, -1.50000000000000, 0.000000000000000)
   2:   (2.00000000000000, 3.50000000000000, 0.000000000000000)
   3:   (1.50000000000000, 2.00000000000000, 3.00000000000000)

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

supports: [0, 1, 2]

loads:
   3:   (75.0000000000000, 100.000000000000, -25.0000000000000)

plot_truss3d (joints, bars, supports, loads)

bar_forces = truss (joints, bars, supports, loads)
bar_forces

{0: 121.55658196570405, 1: -22.492836817374567, 2: -119.27259222775442}

{0: 121.55658196570405, 1: -22.492836817374567, 2: -119.27259222775442}

print_dict (bar_forces)

0:   121.55658196570405
1:   -22.492836817374567
2:   -119.27259222775442

0:   121.55658196570405
1:   -22.492836817374567
2:   -119.27259222775442

plot_truss_with_forces3d (joints, bars, supports, bar_forces)

free_joints = others (supports, joints)
free_joints

[3]

[3]

maxwells_rule (joints, supports, bars)

3 * 1 - 3  ==  0

3 * 1 - 3  ==  0

A = equilibrium_matrix (joints, bars, free_joints)
show (A)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -0.3841106397986879 & 0.4767312946227962 & 0.14744195615489714 \\ -0.5121475197315839 & -0.6674238124719146 & 0.4423258684646914 \\ -0.7682212795973759 & -0.5720775535473553 & -0.8846517369293828 \end{array}\right)$

b = load_vector (loads, free_joints)
show (b)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(75.0,\,100.0,\,-25.0\right)$

x = -A \ b
show (x)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(121.55658196570404,\,-22.492836817374574,\,-119.2725922277544\right)$

bfs = dict_bar_force (bars, x)
print_dict (bfs)

0:   121.55658196570404
1:   -22.492836817374574
2:   -119.2725922277544

0:   121.55658196570404
1:   -22.492836817374574
2:   -119.2725922277544

Ab = augmented_equil_matrix (joints, bars, loads, free_joints)
show (Ab)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} -0.3841106397986879 & 0.4767312946227962 & 0.14744195615489714 & -75.0 \\ -0.5121475197315839 & -0.6674238124719146 & 0.4423258684646914 & -100.0 \\ -0.7682212795973759 & -0.5720775535473553 & -0.8846517369293828 & 25.0 \end{array}\right)$

Abr, pivots = ref_wpp (Ab)
show (Abr)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} 1.0 & 0.7446780878644506 & 1.1515584902738285 & -32.542706982944395 \\ 0.0 & 0.9999999999999999 & 0.7731921410338582 & -114.71346776861033 \\ 0.0 & -5.551115123125783 \times 10^{-17} & 1.0 & -119.27259222775442 \end{array}\right)$

show (Abr.apply_map (NDC (8)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} 1.0000000 & 0.74467809 & 1.1515585 & -32.542707 \\ 0.00000000 & 1.0000000 & 0.77319214 & -114.71347 \\ 0.00000000 & 0.00000000 & 1.0000000 & -119.27259 \end{array}\right)$

pivots

[0, 1, 2]

[0, 1, 2]

xx = ref_solution (Abr, pivots)
xx

(121.55658196570405, -22.492836817374567, -119.27259222775442)

(121.55658196570405, -22.492836817374567, -119.27259222775442)

print_dict (dict_bar_force (bars, xx))

0:   121.55658196570405
1:   -22.492836817374567
2:   -119.27259222775442

0:   121.55658196570405
1:   -22.492836817374567
2:   -119.27259222775442

rs, cs, rk, lk = subspaces (A)

print_subspace (cs)

Column space:
   Vector space of degree 3 and dimension 3
   Basis:
     (-0.3841106397986879, -0.5121475197315839, -0.7682212795973759)
     (0.4767312946227962, -0.6674238124719146, -0.5720775535473553)
     (0.14744195615489714, 0.4423258684646914, -0.8846517369293828)

Column space:
   Vector space of degree 3 and dimension 3
   Basis:
     (-0.3841106397986879, -0.5121475197315839, -0.7682212795973759)
     (0.4767312946227962, -0.6674238124719146, -0.5720775535473553)
     (0.14744195615489714, 0.4423258684646914, -0.8846517369293828)

show (A)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -0.3841106397986879 & 0.4767312946227962 & 0.14744195615489714 \\ -0.5121475197315839 & -0.6674238124719146 & 0.4423258684646914 \\ -0.7682212795973759 & -0.5720775535473553 & -0.8846517369293828 \end{array}\right)$

print_subspace (rs)

Row space:
   Vector space of degree 3 and dimension 3
   Basis:
     (1.0, 0.7446780878644506, 1.1515584902738285)
     (0.0, 0.9999999999999999, 0.7731921410338582)
     (0.0, -5.551115123125783e-17, 1.0)

Row space:
   Vector space of degree 3 and dimension 3
   Basis:
     (1.0, 0.7446780878644506, 1.1515584902738285)
     (0.0, 0.9999999999999999, 0.7731921410338582)
     (0.0, -5.551115123125783e-17, 1.0)

show (Abr.subdivision (0, 0))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1.0 & 0.7446780878644506 & 1.1515584902738285 \\ 0.0 & 0.9999999999999999 & 0.7731921410338582 \\ 0.0 & -5.551115123125783 \times 10^{-17} & 1.0 \end{array}\right)$

print_subspace (rk)

Kernel:
   Vector space of degree 3 and dimension 0

Kernel:
   Vector space of degree 3 and dimension 0

print_subspace (lk)

Left kernel:
   Vector space of degree 3 and dimension 0

Left kernel:
   Vector space of degree 3 and dimension 0

Abrr, pivots = rref_wpp (Ab)
show (Abrr)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} 1.0 & 1.1102230246251565 \times 10^{-16} & 0.0 & 121.55658196570405 \\ 0.0 & 0.9999999999999999 & 0.0 & -22.492836817374567 \\ 0.0 & -5.551115123125783 \times 10^{-17} & 1.0 & -119.27259222775442 \end{array}\right)$

show (Abrr.apply_map (NDC (8)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} 1.0000000 & 0.00000000 & 0.00000000 & 121.55658 \\ 0.00000000 & 1.0000000 & 0.00000000 & -22.492837 \\ 0.00000000 & 0.00000000 & 1.0000000 & -119.27259 \end{array}\right)$

print_dict (dict_bar_force (bars, rref_solution (Abrr, pivots)))

0:   121.55658196570405
1:   -22.492836817374567
2:   -119.27259222775442

0:   121.55658196570405
1:   -22.492836817374567
2:   -119.27259222775442

joints = { 0: (-2.0, -1.0, 0.0), 1: (0.0, 0.0, 0.0), 2: (2.0, 1.0, 0.0), 3: (0.0, 0.0, 2.0) }
print 'joints:'
print_dict (joints, indent = '   ')
print

bars = { 0: (0, 3), 1: (1, 3), 2: (2, 3) }
print 'bars:' 
print_dict (bars, indent = '   ')
print

supports = [0, 1, 2]
print 'supports:', supports
print

loads = { 3: (100.0, 50.0, 0.0) }
print 'loads:' 
print_dict (loads, indent = '   ')

joints:
   0:   (-2.00000000000000, -1.00000000000000, 0.000000000000000)
   1:   (0.000000000000000, 0.000000000000000, 0.000000000000000)
   2:   (2.00000000000000, 1.00000000000000, 0.000000000000000)
   3:   (0.000000000000000, 0.000000000000000, 2.00000000000000)

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

supports: [0, 1, 2]

loads:
   3:   (100.000000000000, 50.0000000000000, 0.000000000000000)

joints:
   0:   (-2.00000000000000, -1.00000000000000, 0.000000000000000)
   1:   (0.000000000000000, 0.000000000000000, 0.000000000000000)
   2:   (2.00000000000000, 1.00000000000000, 0.000000000000000)
   3:   (0.000000000000000, 0.000000000000000, 2.00000000000000)

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

supports: [0, 1, 2]

loads:
   3:   (100.000000000000, 50.0000000000000, 0.000000000000000)

plot_truss3d (joints, bars, supports, loads)

free_joints = others (supports, joints)
free_joints

[3]

[3]

maxwells_rule (joints, supports, bars)

3 * 1 - 3  ==  0

3 * 1 - 3  ==  0

Ab = augmented_equil_matrix (joints, bars, loads, free_joints)
show (Ab)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} -0.6666666666666666 & 0.0 & 0.6666666666666666 & -100.0 \\ -0.3333333333333333 & 0.0 & 0.3333333333333333 & -50.0 \\ -0.6666666666666666 & -1.0 & -0.6666666666666666 & -0.0 \end{array}\right)$

Abr, pivots = ref_wpp (Ab)
show (Abr)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} 1.0 & -0.0 & -1.0 & 150.0 \\ 0.0 & 1.0 & 1.3333333333333333 & -100.0 \\ 0.0 & 0.0 & 0.0 & 0.0 \end{array}\right)$

pivots

[0, 1]

[0, 1]

primary_system_bars = pivots
primary_system_bars

[0, 1]

[0, 1]

redundant_bars = others (pivots, bars)
redundant_bars

[2]

[2]

bar_forces = dict_bar_force (bars, ref_solution (Abr, pivots))
print_dict (bar_forces)

0:   150.0
1:   -100.0
2:   0.0

0:   150.0
1:   -100.0
2:   0.0

plot_truss_with_forces3d (joints, bars, supports, bar_forces)

rs, cs, rk, lk = subspaces (Ab.subdivision (0, 0))

print_subspace (cs)

Column space:
   Vector space of degree 3 and dimension 2
   Basis:
     (-0.6666666666666666, -0.3333333333333333, -0.6666666666666666)
     (0.0, 0.0, -1.0)

Column space:
   Vector space of degree 3 and dimension 2
   Basis:
     (-0.6666666666666666, -0.3333333333333333, -0.6666666666666666)
     (0.0, 0.0, -1.0)

show (Ab.subdivision (0, 0))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -0.6666666666666666 & 0.0 & 0.6666666666666666 \\ -0.3333333333333333 & 0.0 & 0.3333333333333333 \\ -0.6666666666666666 & -1.0 & -0.6666666666666666 \end{array}\right)$

pivots

[0, 1]

[0, 1]

print_subspace (rs)

Row space:
   Vector space of degree 3 and dimension 2
   Basis:
     (1.0, -0.0, -1.0)
     (0.0, 1.0, 1.3333333333333333)

Row space:
   Vector space of degree 3 and dimension 2
   Basis:
     (1.0, -0.0, -1.0)
     (0.0, 1.0, 1.3333333333333333)

show (Abr.subdivision (0, 0))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1.0 & -0.0 & -1.0 \\ 0.0 & 1.0 & 1.3333333333333333 \\ 0.0 & 0.0 & 0.0 \end{array}\right)$

print_subspace (rk)

Kernel:
   Vector space of degree 3 and dimension 1
   Basis:
     (1.0, -1.3333333333333333, 1.0)

Kernel:
   Vector space of degree 3 and dimension 1
   Basis:
     (1.0, -1.3333333333333333, 1.0)

bar_forces_ss = dict_bar_force (bars, rk['basis'][0])
print_dict (bar_forces_ss)

0:   1.0
1:   -1.3333333333333333
2:   1.0

0:   1.0
1:   -1.3333333333333333
2:   1.0

plot_truss_with_forces3d (joints, bars, supports, bar_forces_ss)

print_subspace (lk)

Left kernel:
   Vector space of degree 3 and dimension 1
   Basis:
     (-0.5, 1.0, 0.0)

Left kernel:
   Vector space of degree 3 and dimension 1
   Basis:
     (-0.5, 1.0, 0.0)

plot_truss_with_displacements3d (joints, bars, supports, lk['basis'][0])

b = load_vector (loads, free_joints)
b

(100.0, 50.0, 0.0)

(100.0, 50.0, 0.0)

is_in_subspace (b, cs['basis'])

True

True

c = vector ([-50, 100, 0])
c

(-50, 100, 0)

(-50, 100, 0)

is_in_subspace (c, cs['basis'])

False

False

is_in_subspace (c, lk['basis'])

True

True

c2 = vector ([100, 0, 0])
c2

(100, 0, 0)

(100, 0, 0)

is_in_subspace (c2, cs['basis'])

False

False

is_in_subspace (c2, lk['basis'])

False

False

has_component_in_subspace (c2, lk['basis'])

True

True

loads2 = { 3: (-50.0, 100.0, 0.0) }

plot_truss3d (joints, bars, supports, loads2)

Ab12 = augmented_equil_matrix (joints, bars, loads2, free_joints)
show (Ab12)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} -0.6666666666666666 & 0.0 & 0.6666666666666666 & 50.0 \\ -0.3333333333333333 & 0.0 & 0.3333333333333333 & -100.0 \\ -0.6666666666666666 & -1.0 & -0.6666666666666666 & -0.0 \end{array}\right)$

Abrr12, pivots12 = rref_wpp (Ab12)
show (Abrr12)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} 1.0 & 0.0 & -1.0 & 0.0 \\ 0.0 & 1.0 & 1.3333333333333333 & 0.0 \\ 0.0 & 0.0 & 0.0 & 1.0 \end{array}\right)$

pivots12

[0, 1, 3]

[0, 1, 3]

loads3 = { 3: (100.0, 0.0, 0.0) }

plot_truss3d (joints, bars, supports, loads3)

Ab13 = augmented_equil_matrix (joints, bars, loads3, free_joints)
show (Ab13)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} -0.6666666666666666 & 0.0 & 0.6666666666666666 & -100.0 \\ -0.3333333333333333 & 0.0 & 0.3333333333333333 & -0.0 \\ -0.6666666666666666 & -1.0 & -0.6666666666666666 & -0.0 \end{array}\right)$

Abrr13, pivots13 = rref_wpp (Ab13)
show (Abrr13)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} 1.0 & 0.0 & -1.0 & 0.0 \\ 0.0 & 1.0 & 1.3333333333333333 & 0.0 \\ 0.0 & 0.0 & 0.0 & 1.0 \end{array}\right)$

pivots13

[0, 1, 3]

[0, 1, 3]

bars2 = { 0: (0, 3), 1: (2, 3), 2: (1, 3) }
print 'bars:' 
print_dict (bars2, indent = '   ')
print

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

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

Ab2 = augmented_equil_matrix (joints, bars2, loads, free_joints)
show (Ab2)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} -0.6666666666666666 & 0.6666666666666666 & 0.0 & -100.0 \\ -0.3333333333333333 & 0.3333333333333333 & 0.0 & -50.0 \\ -0.6666666666666666 & -0.6666666666666666 & -1.0 & -0.0 \end{array}\right)$

show (Ab)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} -0.6666666666666666 & 0.0 & 0.6666666666666666 & -100.0 \\ -0.3333333333333333 & 0.0 & 0.3333333333333333 & -50.0 \\ -0.6666666666666666 & -1.0 & -0.6666666666666666 & -0.0 \end{array}\right)$

Abr2, p2 = ref_wpp (Ab2)
show (Abr2)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} 1.0 & -1.0 & -0.0 & 150.0 \\ 0.0 & 1.0 & 0.75 & -75.0 \\ 0.0 & 0.0 & 0.0 & 0.0 \end{array}\right)$

bar_forces2 = dict_bar_force (bars2, ref_solution (Abr2, p2))
print_dict (bar_forces2)

0:   75.0
1:   -75.0
2:   0.0

0:   75.0
1:   -75.0
2:   0.0

plot_truss_with_forces3d (joints, bars2, supports, bar_forces2)

rs2, cs2, rk2, lk2 = subspaces (Ab2.subdivision (0, 0))

print_subspace (cs2)

Column space:
   Vector space of degree 3 and dimension 2
   Basis:
     (-0.6666666666666666, -0.3333333333333333, -0.6666666666666666)
     (0.6666666666666666, 0.3333333333333333, -0.6666666666666666)

Column space:
   Vector space of degree 3 and dimension 2
   Basis:
     (-0.6666666666666666, -0.3333333333333333, -0.6666666666666666)
     (0.6666666666666666, 0.3333333333333333, -0.6666666666666666)

show (Ab2.subdivision (0, 0))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -0.6666666666666666 & 0.6666666666666666 & 0.0 \\ -0.3333333333333333 & 0.3333333333333333 & 0.0 \\ -0.6666666666666666 & -0.6666666666666666 & -1.0 \end{array}\right)$

print_subspace (rs2)

Row space:
   Vector space of degree 3 and dimension 2
   Basis:
     (1.0, -1.0, -0.0)
     (0.0, 1.0, 0.75)

Row space:
   Vector space of degree 3 and dimension 2
   Basis:
     (1.0, -1.0, -0.0)
     (0.0, 1.0, 0.75)

show (Abr2.subdivision (0, 0))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1.0 & -1.0 & -0.0 \\ 0.0 & 1.0 & 0.75 \\ 0.0 & 0.0 & 0.0 \end{array}\right)$

bar_forces_ss2 = dict_bar_force (bars2, rk2['basis'][0])
print_dict (bar_forces_ss2)

0:   -0.75
1:   -0.75
2:   1.0

0:   -0.75
1:   -0.75
2:   1.0

plot_truss_with_forces3d (joints, bars2, supports, bar_forces_ss2)

joints = { 0: (0.0, 0.0, 0.0), 1: (2.0, 0.0, 0.0), 2: (0.0, 0.0, 2.0) }
print 'joints:'
print_dict (joints, indent = '   ')
print

bars = { 0: (0, 2), 1: (1, 2) }
print 'bars:' 
print_dict (bars, indent = '   ')
print

supports = [0, 1]
print 'supports:', supports
print

loads = { 2: (100.0, 0.0, 0.0) }
print 'loads:' 
print_dict (loads, indent = '   ')

joints:
   0:   (0.000000000000000, 0.000000000000000, 0.000000000000000)
   1:   (2.00000000000000, 0.000000000000000, 0.000000000000000)
   2:   (0.000000000000000, 0.000000000000000, 2.00000000000000)

bars:
   0:   (0, 2)
   1:   (1, 2)

supports: [0, 1]

loads:
   2:   (100.000000000000, 0.000000000000000, 0.000000000000000)

joints:
   0:   (0.000000000000000, 0.000000000000000, 0.000000000000000)
   1:   (2.00000000000000, 0.000000000000000, 0.000000000000000)
   2:   (0.000000000000000, 0.000000000000000, 2.00000000000000)

bars:
   0:   (0, 2)
   1:   (1, 2)

supports: [0, 1]

loads:
   2:   (100.000000000000, 0.000000000000000, 0.000000000000000)

plot_truss3d (joints, bars, supports, loads)

free_joints = others (supports, joints)
free_joints

[2]

[2]

maxwells_rule (joints, supports, bars)

3 * 1 - 2  ==  1  !=  0

3 * 1 - 2  ==  1  !=  0

Ab = augmented_equil_matrix (joints, bars, loads, free_joints)
show (Ab)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr|r} 0.0 & 0.7071067811865475 & -100.0 \\ 0.0 & 0.0 & -0.0 \\ -1.0 & -0.7071067811865475 & -0.0 \end{array}\right)$

Abr, pivots = ref_wpp (Ab)
show (Abr)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr|r} 1.0 & 0.7071067811865475 & 0.0 \\ 0.0 & 1.0 & -141.4213562373095 \\ 0.0 & 0.0 & 0.0 \end{array}\right)$

bar_forces = dict_bar_force (bars, ref_solution (Abr, pivots))
print_dict (bar_forces)

0:   100.0
1:   -141.4213562373095

0:   100.0
1:   -141.4213562373095

pivots

[0, 1]

[0, 1]

plot_truss_with_forces3d (joints, bars, supports, bar_forces)

rs, cs, rk, lk = subspaces (Ab.subdivision (0, 0))

print_subspace (cs)

Column space:
   Vector space of degree 3 and dimension 2
   Basis:
     (0.0, 0.0, -1.0)
     (0.7071067811865475, 0.0, -0.7071067811865475)

Column space:
   Vector space of degree 3 and dimension 2
   Basis:
     (0.0, 0.0, -1.0)
     (0.7071067811865475, 0.0, -0.7071067811865475)

show (Ab.subdivision (0, 0))

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

b = load_vector (loads, free_joints)
show (b)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(100.0,\,0.0,\,0.0\right)$

is_in_subspace (b, cs['basis'])

True

True

print_subspace (rs)

Row space:
   Vector space of degree 2 and dimension 2
   Basis:
     (1.0, 0.7071067811865475)
     (0.0, 1.0)

Row space:
   Vector space of degree 2 and dimension 2
   Basis:
     (1.0, 0.7071067811865475)
     (0.0, 1.0)

show (Abr.subdivision (0, 0))

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

print_subspace (rk)

Kernel:
   Vector space of degree 2 and dimension 0

Kernel:
   Vector space of degree 2 and dimension 0

print_subspace (lk)

Left kernel:
   Vector space of degree 3 and dimension 1
   Basis:
     (0.0, 1.0, 0.0)

Left kernel:
   Vector space of degree 3 and dimension 1
   Basis:
     (0.0, 1.0, 0.0)

plot_truss_with_displacements3d (joints, bars, supports, lk['basis'][0], displ_scale = 1.25)

joints = { 0: (-1.0, 0.0), 1: (0.0, 0.0), 2: (4.0, 0.0), 3: (0.0, 2.0) }
print 'joints:'
print_dict (joints, indent = '   ')
print

bars = { 0: (0, 3), 1: (1, 3), 2: (2, 3) }
print 'bars:' 
print_dict (bars, indent = '   ')
print

supports = [0, 1, 2]
print 'supports:', supports
print

loads = { 3: (100.0, 0.0) }
print 'loads:' 
print_dict (loads, indent = '   ')

joints:
   0:   (-1.00000000000000, 0.000000000000000)
   1:   (0.000000000000000, 0.000000000000000)
   2:   (4.00000000000000, 0.000000000000000)
   3:   (0.000000000000000, 2.00000000000000)

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

supports: [0, 1, 2]

loads:
   3:   (100.000000000000, 0.000000000000000)

joints:
   0:   (-1.00000000000000, 0.000000000000000)
   1:   (0.000000000000000, 0.000000000000000)
   2:   (4.00000000000000, 0.000000000000000)
   3:   (0.000000000000000, 2.00000000000000)

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

supports: [0, 1, 2]

loads:
   3:   (100.000000000000, 0.000000000000000)

plot_truss2d (joints, bars, supports, loads)

free_joints = others (supports, joints)
free_joints

[3]

[3]

maxwells_rule (joints, supports, bars)

2 * 1 - 3  ==  -1  !=  0

2 * 1 - 3  ==  -1  !=  0

Ab = augmented_equil_matrix (joints, bars, loads, free_joints)
show (Ab)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} -0.4472135954999579 & 0.0 & 0.8944271909999159 & -100.0 \\ -0.8944271909999159 & -1.0 & -0.4472135954999579 & -0.0 \end{array}\right)$

Abr, pivots = ref_wpp (Ab)
show (Abr)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr|r} 1.0 & 1.118033988749895 & 0.5 & 0.0 \\ 0.0 & 1.0 & 2.23606797749979 & -200.0 \end{array}\right)$

primary_system_bars = pivots
primary_system_bars

[0, 1]

[0, 1]

redundant_bars = others (pivots, bars)
redundant_bars

[2]

[2]

x = ref_solution (Abr, pivots)
x

(223.60679774997897, -200.0, 0.0)

(223.60679774997897, -200.0, 0.0)

bar_forces = dict_bar_force (bars, ref_solution (Abr, pivots))
print_dict (bar_forces)

0:   223.60679774997897
1:   -200.0
2:   0.0

0:   223.60679774997897
1:   -200.0
2:   0.0

plot_truss_with_forces2d (joints, bars, supports, bar_forces)

rs, cs, rk, lk = subspaces (Ab.subdivision (0, 0))

print_subspace (cs)

Column space:
   Vector space of degree 2 and dimension 2
   Basis:
     (-0.4472135954999579, -0.8944271909999159)
     (0.0, -1.0)

Column space:
   Vector space of degree 2 and dimension 2
   Basis:
     (-0.4472135954999579, -0.8944271909999159)
     (0.0, -1.0)

show (Ab.subdivision (0, 0))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -0.4472135954999579 & 0.0 & 0.8944271909999159 \\ -0.8944271909999159 & -1.0 & -0.4472135954999579 \end{array}\right)$

print_subspace (rs)

Row space:
   Vector space of degree 3 and dimension 2
   Basis:
     (1.0, 1.118033988749895, 0.5)
     (0.0, 1.0, 2.23606797749979)

Row space:
   Vector space of degree 3 and dimension 2
   Basis:
     (1.0, 1.118033988749895, 0.5)
     (0.0, 1.0, 2.23606797749979)

show (Abr.subdivision (0, 0))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1.0 & 1.118033988749895 & 0.5 \\ 0.0 & 1.0 & 2.23606797749979 \end{array}\right)$

print_subspace (rk)

Kernel:
   Vector space of degree 3 and dimension 1
   Basis:
     (2.0000000000000004, -2.23606797749979, 1.0)

Kernel:
   Vector space of degree 3 and dimension 1
   Basis:
     (2.0000000000000004, -2.23606797749979, 1.0)

bar_forces_ss = dict_bar_force (bars, rk['basis'][0])
print_dict (bar_forces_ss)

0:   2.0000000000000004
1:   -2.23606797749979
2:   1.0

0:   2.0000000000000004
1:   -2.23606797749979
2:   1.0

plot_truss_with_forces2d (joints, bars, supports, bar_forces_ss)

print_subspace (lk)

Left kernel:
   Vector space of degree 2 and dimension 0
   Basis:
     []

Left kernel:
   Vector space of degree 2 and dimension 0
   Basis:
     []

joints = { 0: (-2.0, 0.0, 0.0), 1: (0.0, 0.0, 0.0), 2: (2.0, 0.0, 0.0), 3: (0.0, -2.0, 0.0), 4: (0.0, 2.0, 0.0), 5: (0.0, 0.0, 2.0) }
print 'joints:'
print_dict (joints, indent = '   ')
print

bars = { 0: (0, 5), 1: (1, 5), 2: (2, 5), 3: (3, 5), 4: (4, 5) }
print 'bars:' 
print_dict (bars, indent = '   ')
print

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

loads = { 5: (100.0, 100.0, 0.0) }
print 'loads:' 
print_dict (loads, indent = '   ')

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

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

supports: [0, 1, 2, 3, 4]

loads:
   5:   (100.000000000000, 100.000000000000, 0.000000000000000)

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

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

supports: [0, 1, 2, 3, 4]

loads:
   5:   (100.000000000000, 100.000000000000, 0.000000000000000)

plot_truss3d (joints, bars, supports, loads)

free_joints = others (supports, joints)
free_joints

[5]

[5]

maxwells_rule (joints, supports, bars)

3 * 1 - 5  ==  -2  !=  0

3 * 1 - 5  ==  -2  !=  0

Ab = augmented_equil_matrix (joints, bars, loads, free_joints)
show (Ab)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr|r} -0.7071067811865475 & 0.0 & 0.7071067811865475 & 0.0 & 0.0 & -100.0 \\ 0.0 & 0.0 & 0.0 & -0.7071067811865475 & 0.7071067811865475 & -100.0 \\ -0.7071067811865475 & -1.0 & -0.7071067811865475 & -0.7071067811865475 & -0.7071067811865475 & -0.0 \end{array}\right)$

Abr, pivots = ref_wpp (Ab)
show (Abr)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr|r} 1.0 & -0.0 & -1.0 & -0.0 & -0.0 & 141.4213562373095 \\ 0.0 & 1.0 & 1.414213562373095 & 0.7071067811865475 & 0.7071067811865475 & -100.0 \\ 0.0 & 0.0 & 0.0 & 1.0 & -1.0 & 141.4213562373095 \end{array}\right)$

primary_system_bars = pivots
primary_system_bars

[0, 1, 3]

[0, 1, 3]

redundant_bars = others (pivots, bars)
redundant_bars

[2, 4]

[2, 4]

bar_forces = dict_bar_force (bars, ref_solution (Abr, pivots))
print_dict (bar_forces)

0:   141.4213562373095
1:   -200.0
2:   0.0
3:   141.4213562373095
4:   0.0

0:   141.4213562373095
1:   -200.0
2:   0.0
3:   141.4213562373095
4:   0.0

plot_truss_with_forces3d (joints, bars, supports, bar_forces)

A = equilibrium_matrix (joints, bars, free_joints)
show (A)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} -0.7071067811865475 & 0.0 & 0.7071067811865475 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & -0.7071067811865475 & 0.7071067811865475 \\ -0.7071067811865475 & -1.0 & -0.7071067811865475 & -0.7071067811865475 & -0.7071067811865475 \end{array}\right)$

rs, cs, rk, lk = subspaces (A)

print_subspace (rs)

Row space:
   Vector space of degree 5 and dimension 3
   Basis:
     (1.0, -0.0, -1.0, -0.0, -0.0)
     (0.0, 1.0, 1.414213562373095, 0.7071067811865475, 0.7071067811865475)
     (0.0, 0.0, 0.0, 1.0, -1.0)

Row space:
   Vector space of degree 5 and dimension 3
   Basis:
     (1.0, -0.0, -1.0, -0.0, -0.0)
     (0.0, 1.0, 1.414213562373095, 0.7071067811865475, 0.7071067811865475)
     (0.0, 0.0, 0.0, 1.0, -1.0)

print_subspace (cs)

Column space:
   Vector space of degree 3 and dimension 3
   Basis:
     (-0.7071067811865475, 0.0, -0.7071067811865475)
     (0.0, 0.0, -1.0)
     (0.0, -0.7071067811865475, -0.7071067811865475)

Column space:
   Vector space of degree 3 and dimension 3
   Basis:
     (-0.7071067811865475, 0.0, -0.7071067811865475)
     (0.0, 0.0, -1.0)
     (0.0, -0.7071067811865475, -0.7071067811865475)

print_subspace (rk)

Kernel:
   Vector space of degree 5 and dimension 2
   Basis:
     (1.0, -1.414213562373095, 1.0, -0.0, 0.0)
     (0.0, -1.414213562373095, 0.0, 1.0, 1.0)

Kernel:
   Vector space of degree 5 and dimension 2
   Basis:
     (1.0, -1.414213562373095, 1.0, -0.0, 0.0)
     (0.0, -1.414213562373095, 0.0, 1.0, 1.0)

bar_forces_ss0 = dict_bar_force (bars, rk['basis'][0])
print_dict (bar_forces_ss0)

0:   1.0
1:   -1.414213562373095
2:   1.0
3:   -0.0
4:   0.0

0:   1.0
1:   -1.414213562373095
2:   1.0
3:   -0.0
4:   0.0

plot_truss_with_forces3d (joints, bars, supports, bar_forces_ss0)

bar_forces_ss1 = dict_bar_force (bars, rk['basis'][1])
print_dict (bar_forces_ss1)

0:   0.0
1:   -1.414213562373095
2:   0.0
3:   1.0
4:   1.0

0:   0.0
1:   -1.414213562373095
2:   0.0
3:   1.0
4:   1.0

plot_truss_with_forces3d (joints, bars, supports, bar_forces_ss1)

print_subspace (lk)

Left kernel:
   Vector space of degree 3 and dimension 0
   Basis:
     []

Left kernel:
   Vector space of degree 3 and dimension 0
   Basis:
     []

joints = { 0: (0.0, 0.0, 0.0), 1: (2.0, 0.0, 0.0), 2: (0.0, 2.0, 0.0), 3: (0.0, 0.0, 2.0) }
print 'joints:'
print_dict (joints, indent = '   ')
print

bars = { 0: (0, 3), 1: (1, 3), 2: (2, 3), 3: (0, 1), 4: (1, 2), 5: (2, 0) }
print 'bars:' 
print_dict (bars, indent = '   ')
print

supports = []

joints:
   0:   (0.000000000000000, 0.000000000000000, 0.000000000000000)
   1:   (2.00000000000000, 0.000000000000000, 0.000000000000000)
   2:   (0.000000000000000, 2.00000000000000, 0.000000000000000)
   3:   (0.000000000000000, 0.000000000000000, 2.00000000000000)

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

joints:
   0:   (0.000000000000000, 0.000000000000000, 0.000000000000000)
   1:   (2.00000000000000, 0.000000000000000, 0.000000000000000)
   2:   (0.000000000000000, 2.00000000000000, 0.000000000000000)
   3:   (0.000000000000000, 0.000000000000000, 2.00000000000000)

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

maxwells_rule (joints, supports, bars)

3 * 4 - 6  ==  6  !=  0

3 * 4 - 6  ==  6  !=  0

free_joints = others (supports, joints)
free_joints

[0, 1, 2, 3]

[0, 1, 2, 3]

plot_truss3d (joints, bars, supports)

A = equilibrium_matrix (joints, bars, free_joints)
show (A.apply_map (ND(4)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 0.0000 & 0.0000 & 0.0000 & 1.000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 1.000 \\ 1.000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & -0.7071 & 0.0000 & -1.000 & -0.7071 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.7071 & 0.0000 \\ 0.0000 & 0.7071 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.7071 & 0.0000 \\ 0.0000 & 0.0000 & -0.7071 & 0.0000 & -0.7071 & -1.000 \\ 0.0000 & 0.0000 & 0.7071 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.7071 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.7071 & 0.0000 & 0.0000 & 0.0000 \\ -1.000 & -0.7071 & -0.7071 & 0.0000 & 0.0000 & 0.0000 \end{array}\right)$

rs, cs, rk, lk = subspaces (A)

print_subspace (rs)

Row space:
   Vector space of degree 6 and dimension 6
   Basis:
     (1.0, 0.0, 0.0, 0.0, 0.0, 0.0)
     (0.0, 1.0, -0.0, 1.4142135623730951, 1.0, -0.0)
     (0.0, 0.0, 1.0, -0.0, 1.0, 1.4142135623730951)
     (0.0, 0.0, 0.0, 1.0, 0.7071067811865475, -0.0)
     (0.0, 0.0, 0.0, 0.0, 1.0, 0.0)
     (0.0, 0.0, 0.0, 0.0, 0.0, 1.0)

Row space:
   Vector space of degree 6 and dimension 6
   Basis:
     (1.0, 0.0, 0.0, 0.0, 0.0, 0.0)
     (0.0, 1.0, -0.0, 1.4142135623730951, 1.0, -0.0)
     (0.0, 0.0, 1.0, -0.0, 1.0, 1.4142135623730951)
     (0.0, 0.0, 0.0, 1.0, 0.7071067811865475, -0.0)
     (0.0, 0.0, 0.0, 0.0, 1.0, 0.0)
     (0.0, 0.0, 0.0, 0.0, 0.0, 1.0)

print_subspace (cs)

Column space:
   Vector space of degree 12 and dimension 6
   Basis:
     (0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0)
     (0.0, 0.0, 0.0, -0.7071067811865475, 0.0, 0.7071067811865475, 0.0, 0.0,
0.0, 0.7071067811865475, 0.0, -0.7071067811865475)
     (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.7071067811865475,
0.7071067811865475, 0.0, 0.7071067811865475, -0.7071067811865475)
     (1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
     (0.0, 0.0, 0.0, -0.7071067811865475, 0.7071067811865475, 0.0,
0.7071067811865475, -0.7071067811865475, 0.0, 0.0, 0.0, 0.0)
     (0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0)

Column space:
   Vector space of degree 12 and dimension 6
   Basis:
     (0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0)
     (0.0, 0.0, 0.0, -0.7071067811865475, 0.0, 0.7071067811865475, 0.0, 0.0, 0.0, 0.7071067811865475, 0.0, -0.7071067811865475)
     (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.7071067811865475, 0.7071067811865475, 0.0, 0.7071067811865475, -0.7071067811865475)
     (1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
     (0.0, 0.0, 0.0, -0.7071067811865475, 0.7071067811865475, 0.0, 0.7071067811865475, -0.7071067811865475, 0.0, 0.0, 0.0, 0.0)
     (0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0)

print_subspace (rk)

Kernel:
   Vector space of degree 6 and dimension 0
   Basis:
     []

Kernel:
   Vector space of degree 6 and dimension 0
   Basis:
     []

print_subspace (lk)

Left kernel:
   Vector space of degree 12 and dimension 6
   Basis:
     (0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0)
     (1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
     (0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0)
     (0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0)
     (0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0)
     (0.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0)

Left kernel:
   Vector space of degree 12 and dimension 6
   Basis:
     (0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0)
     (1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
     (0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0)
     (0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0)
     (0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0)
     (0.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0)

plot_truss_with_displacements3d (joints, bars, supports, lk['basis'][0], displ_scale = 0.325)

plot_truss_with_displacements3d (joints, bars, supports, lk['basis'][1], displ_scale = 0.325)

plot_truss_with_displacements3d (joints, bars, supports, lk['basis'][2], displ_scale = 0.325)

plot_truss_with_displacements3d (joints, bars, supports, lk['basis'][3], displ_scale = 0.325)

plot_truss_with_displacements3d (joints, bars, supports, lk['basis'][4], displ_scale = 0.325)

plot_truss_with_displacements3d (joints, bars, supports, lk['basis'][5], displ_scale = 0.325)

NS_Truss_1_(v0.6)

986 days ago by fresl