GS_program_DiM_(v0.4.1) -- Sage

%auto 

##########################
## Utility functions:
##########################


def bind_1st (f, a) :
    return lambda y : f (a, y)

def bind_2nd (f, b) :
    return lambda x : f (x, b)


def apply_to_ith (f, i) :
    def f_ (l) :
        ll = list (l)
        ll[i] = f(ll[i])
        return ll
    return f_
    
def apply_to_ith_t (f, i) :
    def f_ (t) :
        l = list (t)
        l[i] = f(l[i])
        return tuple (l)
    return f_


def map2 (f, l2) :
    return map (bind_1st (map, f), l2)
     

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

def print_indexed (l) :
    nl = l.nrows() if sage.structure.element.is_Matrix (l) else len (l)
    for i in xrange (nl) :
        print '{0:>3}: {1}'.format (i, l[i])



##########################
## Changing fields:
##########################


def RealFieldD (sig_digits_10) :
    return RealField (round (sig_digits_10 * log (10., 2.)))
    
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 tol
        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 sig_digits (ndigits, to_zero = True, tolerance = None) :
    def f_ (x) :
        return ND(ndigits) (chop(tolerance) (x)) if to_zero else ND(ndigits) (x)
    return f_
    

default_field = RDF

def list_change_field (l, field = default_field) :
    return map (field, l)

def tuple_change_field (t, field = default_field) :
    return tuple (map (field, t))
    
def list2_change_field (l2, field = default_field) :
    return map2 (field, l2)
    
def vector_change_field (v, field = default_field) :
    return v.change_ring (field)
    
def matrix_change_field (A, field = default_field) :  
    return A.change_ring (field)


def nodes_change_field (nds, field = default_field) :
    return list2_change_field (nds, field)

def characteristics_change_field (chrs, field = default_field) :
    return list2_change_field (chrs, field)
    
def nodal_loads_change_field (ndlds, field = default_field) :
    return map (apply_to_ith (bind_1st (map, field), 1), ndlds)
    
def element_loads_change_field (ellds, field = default_field) :
    return [[ld[0], ld[1], [field (val) for val in ld[2]]] for ld in ellds]
    


##########################
## Cross section:
##########################


def A_rect (b, h) :
    return b * h
    
def A_circ (d) :
    pi_n = parent(d) (pi)
    return pi_n * d^2 / 4
    
def I_rect (b, h) :
    return b * h^3 / 12
    
def I_circ (d) :
    pi_n = parent(d) (pi)
    return pi_n * d^4 / 64

    

##########################
## Nodes:
##########################


def ith_node (ni, nodes) :
    return nodes[ni]
    
def node_coords (nd) :
    return nd
    
def ith_node_coords (ni, nodes) :
    return node_coords (ith_node (ni, nodes))
    
def coords_left_to_right (nodes) :
    if parent (nodes[0]) == tuple :
        return map (apply_to_ith_t (lambda x : -x, -1), nodes)
    else :
        return map (apply_to_ith (lambda x : -x, -1), nodes)
    

    
##########################    
## Elements:
##########################

    
def ith_element (ei, elems) :
    return elems[ei]
    
def element_node_indices (el) :
    return el
    
def ith_element_node_indices (ei, elems) :
    return element_node_indices (ith_element (ei, elems))
    
def element_node_coords (el, nodes) :
    ei = element_node_indices (el)
    return [ith_node_coords (i, nodes) for i in ei]
    
def ith_element_node_coords (ei, elems, nodes) :
    return element_node_coords (ith_element (ei, elems), nodes)

 
def ith_element_chars_index (ei, echs) :
    return echs[ei]
    
def ith_characteristics (ci, chs) :
    return chs[ci]
    
def ith_element_characteristics (ei, echs, chs) :
    return ith_characteristics (ith_element_chars_index (ei, echs), chs)   


def node_element_incidences (nodes, elems) :
    nds_els = [[] for nd in nodes]
    for j in xrange (len (elems)) :
        i0, i1 = element_node_indices (elems[j])
        nds_els[i0].append (j)
        nds_els[i1].append (j)
    return nds_els



##########################    
## Supports:
##########################


def support_node_index (support) :
    return support[0]
    
def ith_support_node_index (si, supports) :
    return support_node_index (supports[si])
    
def support_conditions (support) :
    return support[1]
    
def ith_support_conditions (si, supports) :
    return support_conditions (supports[si])



##########################    
## Loads:
##########################


def loaded_node_index (nd_ld) :
    return nd_ld[0]
    
def ith_loaded_node_index (ni, nd_lds) :
    return loaded_node_index (nd_lds[ni])
    
def nodal_load_values (nd_ld) :
    return nd_ld[1]
    
def ith_nodal_load_values (ni, nd_lds) :
    return nodal_load_values (nd_lds[ni])


def loaded_element_index (el_ld) :
    return el_ld[0]
    
def ith_loaded_element_index (ei, el_lds) :
    return loaded_element_index (el_lds[ei])
    
def element_load_type (el_ld) :
    return el_ld[1]
    
def ith_element_load_type (ei, el_lds) :
    return element_load_type (el_lds[ei])
    
def element_load_values (el_ld) :
    return el_ld[2]

def ith_element_load_values (ei, el_lds) :
    return element_load_values (el_lds[ei])
    

    
#############################################
##  Simple generators:
#############################################


def make_nodes (nnds, nd1, dd) :
    return [map (lambda x, y : x + y, nd1, map (lambda x : i * x, dd)) for i in xrange (nnds)]

def make_elements (nels, el1nds, steps) :
    return [map (lambda x, y : x + y, el1nds, map (lambda x : i * x, steps))
            for i in xrange (nels)]       

    

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


def nodal_DsOF (n_nodes, supports, hinges = []) :

    n_dof_per_node = len (ith_support_conditions (0, supports))
    
    dof_table = ones_matrix (ZZ, n_nodes, n_dof_per_node)
    for s in supports :
        dof_table[support_node_index(s)] = support_conditions(s)
    for h in hinges :
        dof_table[h, 2] = 0
    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 (el, doft) :
    ei = element_node_indices (el)
    nn = len (ei)
    nd = doft.ncols()
    ix = matrix (ZZ, nn, nd)
    for i in xrange (nn) :
        ix[i] = doft[ ei[i] ]
    return ix.list()
    


############################    
## Hinges & end conditions:
############################    


def end_conditions (elems, releases, hinges, ne_inc) :
    rels = dict (releases)
    for h in hinges :
        els = ne_inc[h]
        for el in els :
            n0, n1 = ith_element_node_indices (el, elems)
            re = [0, 0]
            if el in rels :
                 re = list (rels[el])
            if h == n0 :
                re[0] = 1
                rels[el] = tuple (re)
            else :
                re[1] = 1
                rels[el] = tuple (re)

    ec = [(0, 0) for ee in range (len (elems))]
    for r in rels :
        ec[r] = rels[r]
    return ec
    
    
def all_hinges (hinges, releases, elems, ne_inc) :
    hh = {}
    hinges2 = copy (hinges)
    for r in releases :
        el = r[0]
        r0, r1 = r[1]
        n0, n1 = ith_element_node_indices (el, elems)
        if r0 == 1 :
            if n0 in hh :
                hh[n0].remove (el)
            else :
                ne = copy (ne_inc[n0])
                ne.remove (el)
                hh[n0] = ne
        if r1 == 1 :
            if n1 in hh :
                hh[n1].remove (el)
            else :
                ne = copy (ne_inc[n1])
                ne.remove (el)
                hh[n1] = ne
    for h in hh :
        if hh[h] == [] :
            if h not in hinges2 :
                hinges2.append (h)
    return hinges2                



############################    
## 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 ((x1 - x0)^2 + (z1 - z0)^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_cos1 (el_nd_coords) :
    [[x0, y0, z0], [x1, y1, z1]] = el_nd_coords
    return (x1 - x0)/element3D_length (el_nd_coords)
    
def element3D_cos2 (el_nd_coords) :
    [[x0, y0, z0], [x1, y1, z1]] = el_nd_coords
    return (y1 - y0)/element3D_length (el_nd_coords)

def element3D_cos3 (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 ((x1 - x0)^2 + (y1 - y0)^2 + (z1 - z0)^2)
    return (dx/l, dy/l, dz/l)
    


##############################    
## Element stiffness matrices:
##############################


def frame2D_element_stiffness_matrix_ (ee, aa, ii, ll) :
    k = zero_matrix (parent (ll), 6, 6)
    k[0,0] = k[3,3] = ee * aa / ll
    k[0,3] = k[3,0] = -k[0,0]
    k[1,1] = 12 * ee * ii / ll^3
    k[1,2] = k[2,1] = -6 * ee * ii /ll^2
    k[1,4] = k[4,1] = -k[1,1]
    k[1,5] = k[5,1] = k[1,2]
    k[2,2] = 4 * ee * ii / ll
    k[2,4] = k[4,2] = -k[1,2]
    k[2,5] = k[5,2] = 2 * ee * ii / ll
    k[4,4] = k[1,1]
    k[4,5] = k[5,4] = -k[1,2]
    k[5,5] = k[2,2]
    return k

def frame2D_element_stiffness_matrix (ch, nc) :
    (ee, aa, ii) = ch
    ll = element2D_length (nc)
    return frame2D_element_stiffness_matrix_ (ee, aa, ii, ll)
    

def truss2D_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
    
def truss2D_element_stiffness_matrix (ch, nc) :
    (ee, aa) = ch
    ll = element2D_length (nc)
    return truss2D_element_stiffness_matrix_ (ee, aa, ll)


def truss3D_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
    
def truss3D_element_stiffness_matrix (ch, nc) :
    (ee, aa) = ch
    ll = element3D_length (nc)
    return truss3D_element_stiffness_matrix_ (ee, aa, ll)
    

def local_stiffness_matrix (ei, fk, nodes, elems, echs, chars) :
    enc = ith_element_node_coords (ei, els, nds)
    ech = ith_element_characteristics (ei, echs, chars)
    ekl = fk (ech, enc)
    return ekl



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


def frame2D_element_global_to_local (el_nd_coords) :
    (c, s) = element2D_cos_sin (el_nd_coords)
    rr = zero_matrix (parent (c), 6, 6)
    rr[0,0] = rr[1,1] = rr[3,3] = rr[4,4] = c
    rr[2,2] = rr[5,5] = 1
    rr[0,1] = rr[3,4] = -s
    rr[1,0] = rr[4,3] = s
    return rr

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) :
    (c1, c2, c3) = element3D_cos_all (el_nd_coords)
    rr = zero_matrix (parent (c1), 2, 6)
    rr[0,0] = rr[1,3] = c1
    rr[0,1] = rr[1,4] = c2
    rr[0,2] = rr[1,5] = c3
    return rr


def transform_global_to_local (ei, fr, nodes, elems) :
    return fr (ith_element_node_coords (ei, elems, nodes))



##############################    
## Fixed end forces:
##############################


def distributed_force2D_ (qx, qz, ll) :
    f = vector (parent(ll), 6)
    f[0] = f[3] = -qx * ll / 2
    f[1] = f[4] = -qz * ll / 2
    f[2] = qz * ll^2 / 12
    f[5] = -f[2]
    return f
    
def distributed_force2D (qq, ll) :
    [qx, qz] = qq
    return distributed_force2D_ (qx, qz, ll)
    
    
def concentrated_force2D_ (fx, fz, m, al, ll) :
    f = vector (parent(ll), 6)
    a = al*ll
    b = ll - a
    f[0] = -fx * b / ll
    f[3] = -fx * a / ll
    f[1] = -fz * (3*a + b) * b^2 / ll^3 - 6*m*a*b / ll^3
    f[4] = -fz * (a + 3*b) * a^2 / ll^3 + 6*m*a*b / ll^3
    f[2] = fz * a * b^2 / ll^2 + m * b * (3*a - ll) / ll^2
    f[5] = -fz * a^2 * b / ll^2 + m * a * (3*b - ll) / ll^2
    return f
    
def concentrated_force2D (ff, ll) :
    [fx, fz, m, al] = ff
    return concentrated_force2D_ (fx, fz, m, al, ll)
    
    
def fixed_end_force_vector (f, ff, enc) :
    ll = element2D_length (enc)
    return f (ff, ll)

def local_fixed_end_forces (ei, elem_loads, nodes, elems) :
    elds = []
    for el in elem_loads :
        if ei == loaded_element_index (el) :
            elds.append (el)  
    if elds == [] :
        return zero_vector (6)
    enc = ith_element_node_coords (ei, elems, nodes)
    ff = vector (parent(enc[0][0]), 6)
    for el in elds :
        f = element_load_type (el)
        ldv = element_load_values (el)
        ff += fixed_end_force_vector (f, ldv, enc)
    return ff



################################    
## Static condensation:
################################


def condense_stiffness_matrix (kk, dof) :
    crw = kk.row (dof)
    ccl = kk.column (dof)
    ccl /= ccl[dof]
    for i in xrange (kk.nrows()) :
        kk[i] -= ccl[i] * crw
    return ccl

def condense_force_vector (ff, dof, kc) :
    fc = ff[dof]
    for i in xrange (ff.degree()) :
        ff[i] -= fc * kc[i]
    return

    
################################    
## Structure stiffness matrix:
################################


def structure_stiffness_matrix_ (fk, fr, ndof, tdof, nodes, elems, echs, chars) :
    field = parent (ith_node_coords (0, nodes)[0])
    sk = zero_matrix (field, ndof, ndof)
    for j in xrange (len (elems)) :
        enc = ith_element_node_coords (j, elems, nodes)
        ech = ith_element_characteristics (j, echs, chars)
        ekl = fk (ech, enc)
        er = fr (enc)
        ekg = er.T * ekl * er
        eix = element_DsOF (ith_element_node_indices (j, elems), 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


def structure_stiffness_matrix_c_ (fk, fr, ndof, tdof, nodes, elems, echs, chars, end_conds) :
    field = parent (ith_node_coords (0, nodes)[0])
    sk = zero_matrix (field, ndof, ndof)
    cond = []
    for j in xrange (len (elems)) :
        enc = ith_element_node_coords (j, elems, nodes)
        ech = ith_element_characteristics (j, echs, chars)
        ekl = fk (ech, enc)

        ec0, ec1 = end_conds[j]
        if ec0 == 1 :
            cr0 = condense_stiffness_matrix (ekl, 2)
        else :
            cr0 = []
        if ec1 == 1 :
            cr1 = condense_stiffness_matrix (ekl, 5)
        else :
            cr1 = []
        cond.append ((cr0, cr1))

        er = fr (enc)
        ekg = er.T * ekl * er
        eix = element_DsOF (ith_element_node_indices (j, elems), 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, cond


def structure_stiffness_matrix (fk, fr, ndof, tdof, nodes, elems, echs, chars, end_cs = []) :
    if end_cs == [] :
        return structure_stiffness_matrix_ (fk, fr, ndof, tdof, nodes, elems, echs, chars)
    else :
        return structure_stiffness_matrix_c_ (fk, fr, ndof, tdof, nodes, elems, echs, chars, end_cs)


    
##############################    
## Structure force vector:
##############################


def structure_force_vector_nodes (ndof, tdof, nodal_loads) :
    field = parent (ith_nodal_load_values (0, nodal_loads)[0])
    sf = zero_vector (field, ndof)
    for ndl in nodal_loads :
        for j in xrange (tdof.ncols()) :
            d = tdof[loaded_node_index (ndl), j]
            if d != -1 :
                sf[d] += nodal_load_values (ndl)[j]  
    return sf
     

def structure_force_vector_elems_ (ndof, tdof, elem_loads, elems, nodes, fr) :
    field = parent (ith_element_load_values (0, elem_loads)[0])
    sf = zero_vector (field, ndof)
    for eld in elem_loads :
        ei = loaded_element_index (eld)
        f = element_load_type (eld)
        ldv = element_load_values (eld)
        el = ith_element (ei, elems)
        enc = element_node_coords (el, nodes)
        efl = fixed_end_force_vector (f, ldv, enc)
        er = fr (enc)
        efg = er.T * efl
        ix = element_DsOF (element_node_indices (el), tdof)
        for j in xrange (len(ix)) :
            if ix[j] != -1 :
                sf[ ix[j] ] -= efg[j]
    return sf


def structure_force_vector_elems_c_ (ndof, tdof, elem_loads, elems, nodes, fr, cond) :
    field = parent (ith_element_load_values (0, elem_loads)[0])
    sf = zero_vector (field, ndof)
    for eld in elem_loads :
        ei = loaded_element_index (eld)
        f = element_load_type (eld)
        ldv = element_load_values (eld)
        el = ith_element (ei, elems)
        enc = element_node_coords (el, nodes)
        efl = fixed_end_force_vector (f, ldv, enc)

        ec0, ec1 = cond[ei]
        if ec0 != [] :
            condense_force_vector (efl, 2, ec0)
        if ec1 != [] :
            condense_force_vector (efl, 5, ec1)         

        er = fr (enc)
        efg = er.T * efl
        ix = element_DsOF (element_node_indices (el), tdof)
        for j in xrange (len(ix)) :
            if ix[j] != -1 :
                sf[ ix[j] ] -= efg[j]
    return sf
 

def structure_force_vector_elems (ndof, tdof, elem_loads, elems, nodes, fr, cond = []) :
    if cond == [] :
        return structure_force_vector_elems_ (ndof, tdof, elem_loads, elems, nodes, fr)
    else :
        return structure_force_vector_elems_c_ (ndof, tdof, elem_loads, elems, nodes, fr, cond)


def structure_force_vector (ndof, tdof, nodal_loads = [],
                            elem_loads = [], elems = [], nodes = [], fr = [], cond = []) :
    if nodal_loads != [] :
        sfn = structure_force_vector_nodes (ndof, tdof, nodal_loads)
    if elem_loads != [] :
        sfe = structure_force_vector_elems (ndof, tdof, elem_loads, elems, nodes, fr, cond)
    if nodal_loads != [] and elem_loads != [] :
        sf = sfn + sfe
    elif nodal_loads != [] :
        sf = sfn
    elif elem_loads != [] :
        sf = sfe
    return sf



##############################   
## Kinematic condensation:
##############################


def transl_DsOF (tdof) :
    ttdof = -ones_matrix (ZZ, tdof.nrows(), 2)
    ntdof = 0
    for j in xrange (tdof.nrows()) :
        if tdof[j, 0] != -1 :
            ttdof[j, 0] = ntdof
            ntdof += 1
        if tdof[j, 1] != -1 :
            ttdof[j, 1] = ntdof
            ntdof += 1
    return (ntdof, ttdof)
    

def kinematic_constraints_matrix (ntdof, ttdof, nodes, elems) :
    field = parent (ith_node_coords (0, nodes)[0])
    cnstr = zero_matrix (field, len (elems), ntdof)
    for j in xrange (len (elems)) :
        [e0, e1] = ith_element_node_indices (j, elems)
        enc = ith_element_node_coords (j, elems, nodes)
        (c, s) = element2D_cos_sin (enc)
        i = ttdof[e0, 0]
        if i != -1 :
            cnstr[j, i] = c
        i = ttdof[e0, 1]
        if i != -1 :
            cnstr[j, i] = -s
        i = ttdof[e1, 0]
        if i != -1 :
            cnstr[j, i] = -c
        i = ttdof[e1, 1]
        if i != -1 :
            cnstr[j, i] = s
    return cnstr


def masters (Kc) :
    return Kc.nonpivots()
    
def slaves (Kc) :
    return Kc.pivots()

def print_masters (ms, ttdof) :
    mss = set (ms)
    print 'masters:'
    for i in xrange (ttdof.nrows()) :
        for j in xrange (ttdof.ncols()) :
            if ttdof[i,j] in mss :
                print '   node', i, ': ',
                if j == 0 :
                    print 'horizontal displacement'
                else :
                    print 'vertical displacement'
    return


def kinematic_constraints_kernel (kc) :
    return column_matrix (kc . right_kernel (basis = 'pivot') . matrix())


def kinematic_condensation_matrix_ (ndof, ntdof, tdof, ttdof, ms, hh) :
    nrdof = ndof - ntdof
    nms = len (ms)
    cc = zero_matrix (hh.base_ring(), ndof, nrdof + nms)
    lt = []
    lm = []
    i = 0
    m = ms[0]
    cl = 0
    cc = zero_matrix (hh.base_ring(), ndof, nrdof + nms)
    for j in xrange (tdof.nrows()) :
        (d0, d1, d2) = tdof[j]
        (td0, td1) = ttdof[j]
        # if i < nms :
        #    m = ms[i]
        if d0 != -1 :
            lt.append (d0)
            if td0 == m :
                lm.append (cl)
                cl += 1
                i += 1
                if i < nms :
                    m = ms[i]    
        if d1 != -1 :
            lt.append (d1)
            if td1 == m :
                lm.append (cl)
                cl += 1
                i += 1
                if i < nms :
                    m = ms[i]
        if d2 != -1 :
            cc[d2, cl] = 1
            cl += 1
    for i in xrange (len (lt)) :
        for j in xrange (len (lm)) :
            cc[lt[i], lm[j]] = hh[i, j]
    return cc


def kinematic_condensation_matrix (ndof, tdof, nodes, elems) :
    ntdof, ttdof = transl_DsOF (tdof)
    kc = kinematic_constraints_matrix (ntdof, ttdof, nodes, elems)
    hh = kinematic_constraints_kernel (kc)
    ms = masters (kc)
    kk = kinematic_condensation_matrix_ (ndof, ntdof, tdof, ttdof, ms, hh)
    return ms, kk

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


def table_of_nodal_displacements (tdof, uu, uf = []) :
    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
    
def nodal_displacements (sk, sf, tdof, uf = []) :
    d = sk.solve_right (sf)
    return table_of_nodal_displacements (tdof, d, uf)
    
    
    
##################################    
## Displacements of element ends:
##################################


def element_displacements_global (elni, uu) :
    nn = len (elni)
    nd = uu.ncols()
    field = parent (uu[0,0])
    edg = zero_matrix (field, nn, nd)
    for i in xrange (nn) :
        edg[i] = uu[ elni[i] ]
    return vector (edg.list())
    
def element_displacements_local (elni, uu, rr) :
    edg = element_displacements_global (elni, uu)
    return rr * edg
    
    

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


def all_end_forces_ (fk, fr, uu, nodes, elems, echs, chars, elem_loads = []) :
    efs = []
    for j in xrange (len (elems)) :
        el = ith_element (j, elems)
        enc = element_node_coords (el, nodes)
        ec = ith_element_characteristics (j, echs, chars)
        ekl = fk (ec, enc)
        er = fr (enc)
        eni = element_node_indices (el)
        edl = element_displacements_local (eni, uu, er)
        efs.append (ekl * edl)

    if elem_loads != [] :
        for eld in elem_loads :
            ei = loaded_element_index (eld)
            f = element_load_type (eld)
            ldv = element_load_values (eld)
            enc = ith_element_node_coords (ei, elems, nodes)
            efl = fixed_end_force_vector (f, ldv, enc)
            efs[ei] += efl

    return matrix (efs)
    
    
def all_end_forces_c_ (fk, fr, uu, nodes, elems, echs, chars, end_conds, elem_loads = []) :
    efs = []
    cond = []
    for j in xrange (len (elems)) :
        el = ith_element (j, elems)
        enc = element_node_coords (el, nodes)
        ec = ith_element_characteristics (j, echs, chars)
        ekl = fk (ec, enc)
        ec0, ec1 = end_conds[j]
        if ec0 == 1 :
            cr0 = condense_stiffness_matrix (ekl, 2)
        else :
            cr0 = []
        if ec1 == 1 :
            cr1 = condense_stiffness_matrix (ekl, 5)
        else :
            cr1 = []
        cond.append ((cr0, cr1))
        er = fr (enc)
        eni = element_node_indices (el)
        edl = element_displacements_local (eni, uu, er)
        efs.append (ekl * edl)

    if elem_loads != [] :
        for eld in elem_loads :
            ei = loaded_element_index (eld)
            f = element_load_type (eld)
            ldv = element_load_values (eld)
            enc = ith_element_node_coords (ei, elems, nodes)
            efl = fixed_end_force_vector (f, ldv, enc)
            ec0, ec1 = cond[ei]
            if ec0 != [] :
                condense_force_vector (efl, 2, ec0)
            if ec1 != [] :
                condense_force_vector (efl, 5, ec1)         
            efs[ei] += efl

    return matrix (efs)
    
    
def all_end_forces (fk, fr, uu, nodes, elems, echs, chars, end_cs = [], elem_loads = []) :
    if end_cs == [] :
        return all_end_forces_ (fk, fr, uu, nodes, elems, echs, chars, elem_loads)
    else :
        return all_end_forces_c_ (fk, fr, uu, nodes, elems, echs, chars, end_cs, elem_loads)
    

    
##########################
##  Main function (DiM):
##########################
    
    
def DiM (kk, rr, nodes, elems, supports, chars, map_el_ch, hinges = [], releases = [],
         nodal_loads = [], elem_loads = [], field = default_field) :

    nds2 = coords_left_to_right (nodes_change_field (nodes, field))
    chrs2 = characteristics_change_field (chars, field)
    ndlds2 = nodal_loads_change_field (nodal_loads, field)
    ellds2 = element_loads_change_field (elem_loads, field)

    if hinges == [] and releases == [] :
        endc = []
        nn, tbl = nodal_DsOF (len (nds2), supports)
        sk = structure_stiffness_matrix_ (kk, rr, nn, tbl, nds2, elems, map_el_ch, chrs2)
        sf = structure_force_vector (nn, tbl, ndlds2, ellds2, elems, nds2, rr)
    else :
        ne_inc = node_element_incidences (nds2, elems)
        endc = end_conditions (elems, releases, hinges, ne_inc)
        hinges = all_hinges (hinges, releases, elems, ne_inc)
        nn, tbl = nodal_DsOF (len (nds2), supports, hinges)
        sk, cn = structure_stiffness_matrix_c_ (kk, rr, nn, tbl, nds2, elems, map_el_ch, chrs2, endc)
        sf = structure_force_vector (nn, tbl, ndlds2, ellds2, elems, nds2, rr, cn)

    uu = nodal_displacements (sk, sf, tbl)
    ef = all_end_forces (kk, rr, uu, nds2, elems, map_el_ch, chrs2, endc, ellds2)

    return (uu, ef)


def IDiM (kk, rr, nodes, elems, supports, chars, map_el_ch, hinges = [], releases = [],
         nodal_loads = [], elem_loads = [], field = default_field) :

    nds2 = coords_left_to_right (nodes_change_field (nodes, field))
    chrs2 = characteristics_change_field (chars, field)
    ndlds2 = nodal_loads_change_field (nodal_loads, field)
    ellds2 = element_loads_change_field (elem_loads, field)

    if hinges == [] and releases == [] :
        endc = []
        nn, tbl = nodal_DsOF (len (nds2), supports)
        sk = structure_stiffness_matrix_ (kk, rr, nn, tbl, nds2, elems, map_el_ch, chrs2)
        sf = structure_force_vector (nn, tbl, ndlds2, ellds2, elems, nds2, rr)
    else :
        ne_inc = node_element_incidences (nds2, elems)
        endc = end_conditions (elems, releases, hinges, ne_inc)
        hinges = all_hinges (hinges, releases, elems, ne_inc)
        nn, tbl = nodal_DsOF (len (nds2), supports, hinges)
        sk, cn = structure_stiffness_matrix_c_ (kk, rr, nn, tbl, nds2, elems, map_el_ch, chrs2, endc)
        sf = structure_force_vector (nn, tbl, ndlds2, ellds2, elems, nds2, rr, cn)
    
    ms, C = kinematic_condensation_matrix (nn, tbl, nds2, els)
    CtKC = C.T * sk * C
    Ctq = C.T * sf
    um = CtKC.solve_right (Ctq)
    ua = C * um
    uu = table_of_nodal_displacements (tbl, ua)
    ef = all_end_forces (kk, rr, uu, nds2, elems, map_el_ch, chrs2, endc, ellds2)

    return (uu, ef)

    
    
##########################
##  Graphics:
##########################


def plot_node (pt, rr, color = 'teal') :
    return polygon2d ([(pt[0]-rr, pt[1]-rr), (pt[0]+rr, pt[1]-rr),
                       (pt[0]+rr, pt[1]+rr), (pt[0]-rr, pt[1]+rr)],
                      color = color, fill = True, aspect_ratio = 1, zorder = 10)

def plot_nodes (pts, rr, color = 'teal') :
    return sum ([plot_node (pt, rr, color) for pt in pts])


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_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_end_releases (releases, elems, nodes, hinges, rr, edgecolor = 'teal') :
    prls = []
    rc = 5*rr/2
    for r in releases :
        ei = r[0]
        r0, r1 = r[1]
        n0, n1 = ith_element_node_indices (ei, elems)
        pt0, pt1 = ith_element_node_coords (ei, elems, nodes)
        dx = pt1[0] - pt0[0]
        dy = pt1[1] - pt0[1]
        ll = sqrt (dx^2 + dy^2)
        dx /= ll
        dy /= ll
        if r0 == 1 :
            if n0 not in hinges :
                h = (pt0[0] + rc*dx, pt0[1] + rc*dy)
                prls.append (plot_hinge (h, rr, edgecolor))
        if r1 == 1 :
            if n1 not in hinges :
                h = (pt1[0] - rc*dx, pt1[1] - rc*dy)
                prls.append (plot_hinge (h, rr, edgecolor))
    return sum (prls)          


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_support_r (pt, rr, color = 'green', thickness = 3) :
    return line2d ((pt, (pt[0]+rr, pt[1])), color = color, thickness = thickness, aspect_ratio = 1) \
           + arc (pt, rr, sector = (0, 7*pi/4), color = color, thickness = thickness,
                  aspect_ratio = 1)

def plot_supports (supports, nodes, length, color = 'green', thickness = 3) :
    psps = []
    for sp in supports :
        spn = ith_node_coords (support_node_index (sp), nodes)
        spc = support_conditions (sp)
        if spc[0] == 0:
            psps.append (plot_support_x (spn, length, color = color, thickness = thickness))
        if spc[1] == 0:
            psps.append (plot_support_y (spn, length, color = color, thickness = thickness))
        if len (spc) == 3 and spc[2] == 0:
            psps.append (plot_support_r (spn, 2*length/3, color = color, thickness = thickness))
    return sum (psps)


def plot_force (ff, head, scale, color = 'crimson', width = 5) :
    tail = head[0] - scale*ff[0], head[1] + scale*ff[1]
    return arrow2d (tail, head, width = width, color = color, aspect_ratio = 1, zorder = 12)
    
def plot_moment (mm, pt, scale, color = 'crimson', width = 4) :
    rr = scale*mm
    bc = 0.551915024494
    br = rr*bc
    if mm > 0 :
        bpp = [[(rr, 0), (rr, br), (br, rr), (0, rr)],
               [(-br, rr), (-rr, br), (-rr, 0)],
               [(-rr, -br), (-br, -rr), (0, -rr)]]
    else :
        bpp = [[(-rr, 0), (-rr, br), (-br, rr), (0, rr)],
               [(br, rr), (rr, br), (rr, 0)],
               [(rr, -br), (br, -rr), (0, -rr)]]
    bp = map2 (lambda x : (x[0]+pt[0], x[1]+pt[1]), bpp)
    return arrow2d (path = bp, width = width, color = color, aspect_ratio = 1, zorder = 12)

def plot_nodal_loads (loads, nodes, scale_f, scale_m = None, color = 'crimson', width = 5) :
    pls = []
    if scale_m == None :
        scale_m = 3*scale_f/4
    for li in loads :
        head = ith_node_coords (li[0], nodes)
        ff = li[1]
        if ff[0] != 0 or ff[1] != 0 :
            pls.append (plot_force ((ff[0], ff[1]), head, scale_f, color = color, width = width))
        if len (ff) == 3 and ff[2] != 0 :
            pls.append (plot_moment (ff[2], head, scale_m, color = color, width = 4*width/5))
    return sum (pls)

def plot_distributed (ld, elem, scale, color = 'crimson') :
    pt0, pt1 = elem
    dx = pt1[0] - pt0[0]
    dy = pt1[1] - pt0[1]
    ll = sqrt (dx^2 + dy^2)
    dx /= ll
    dy /= ll
    hx = 4*scale*ld[0]/3
    hy = 4*scale*ld[1]/3
    pt2 = (pt1[0] - hy*dy - hx*dx, pt1[1] + hy*dx - hx*dy)
    pt3 = (pt0[0] - hy*dy - hx*dx, pt0[1] + hy*dx - hx*dy)
    if ld[1] != 0 :
        cl = Color (color) . lighter (0.125)
        pld = polygon2d ([pt0, pt1, pt2, pt3], color = cl, alpha = 0.5, fill = True, aspect_ratio = 1)
    else :
        cl = Color (color) . lighter (0.325)
        pld = arrow2d (pt3, pt2, color = cl, width = 2, aspect_ratio = 1, zorder = 21)
    return pld
    
def plot_concentrated (ff, elem, scale, color = 'crimson', width = 5) :
    [fx, fz, m, al] = ff
    pt0, pt1 = elem
    dx = pt1[0] - pt0[0]
    dy = pt1[1] - pt0[1]
    ll = sqrt (dx^2 + dy^2)
    pt2 = (pt0[0] + al*dx, pt0[1] + al*dy)
    plcf = plcm = None
    if fx != 0 or fz != 0 :
        dx /= ll
        dy /= ll
        hx = scale*fx
        hy = scale*fz
        pt3 = (pt2[0] - hy*dy - hx*dx, pt2[1] + hy*dx - hx*dy)
        plcf = arrow2d (pt3, pt2, width = width, color = color, aspect_ratio = 1, zorder = 20)
    if m != 0 :
        scale_m = 3*scale/4
        plcm = plot_moment (m, pt2, scale_m, color = color, width = 4*width/5)
    if plcf != None and plcm != None :
        plc = plcf + plcm
    elif plcf != None :
        plc = plcf
    elif plcm != 0 :
        plc = plcm
    return plc
    
def plot_elem_loads (loads, nodes, elems, scale, color = 'crimson', width = 5) :
    pls = []
    for ld in loads :
        ei, f, v = ld
        ec = ith_element_node_coords (ei, elems, nodes)
        if f.func_name == distributed_force2D.func_name :
            pls.append (plot_distributed (v, ec, scale, color = color))
        elif f.func_name == concentrated_force2D.func_name :
            pls.append (plot_concentrated (v, ec, scale, color = color, width = width))
    return sum (pls)
    
    
def plot_frame (nodes, elems, supports = [], hinges = [], releases = [],
                nodal_loads = [], elem_loads = [], load_scale = 0.01) :
    prr = point2d (nodes, aspect_ratio = 1)
    rrx = (prr.xmax() - prr.xmin())
    rry = (prr.ymax() - prr.ymin())
    rr = rrx if rrx > rry else rry*3/2
    rn = rr/125
    rh = 4*rn/3
    re = 4*rn/5
    lsp = rr/29
    pt = plot_nodes (nodes, rn)
    pt += plot_elements (elems, nodes)
    if supports != [] :
        pt += plot_supports (supports, nodes, lsp)
    if hinges != [] :
        hngs = [ith_node_coords (hi, nodes) for hi in hinges]
        pt += plot_hinges (hngs, rh)
    if releases != [] :
        pt += plot_end_releases (releases, elems, nodes, hinges, re)
    if nodal_loads != [] :
        pt += plot_nodal_loads (nodal_loads, nodes, load_scale)
    if elem_loads != [] :
        pt += plot_elem_loads (elem_loads, nodes, elems, load_scale)
    return pt    

    
def plot_truss (nodes, elems, supports = [], nodal_loads = [], load_scale = 0.01) :
    prr = point2d (nodes, aspect_ratio = 1)
    rrx = (prr.xmax() - prr.xmin())
    rry = (prr.ymax() - prr.ymin())
    rr = rrx if rrx > rry else rry*3/2
    rn = rr/125
    lsp = rr/23
    pt = plot_hinges (nodes, rn)
    pt += plot_elements (elems, nodes)
    if supports != [] :
        pt += plot_supports (supports, nodes, lsp)
    if nodal_loads != [] :
        pt += plot_nodal_loads (nodal_loads, nodes, load_scale)
    return pt

%auto

# Kratice:

f2k = frame2D_element_stiffness_matrix
f2r = frame2D_element_global_to_local
t2k = truss2D_element_stiffness_matrix
t2r = truss2D_element_global_to_local
t3k = truss3D_element_stiffness_matrix
t3r = truss3D_element_global_to_local

df2 = distributed_force2D
cf2 = concentrated_force2D

nds_xy = [(0., 0.), (8., 0.), (3., 4.), (8., 4.)]
print 'nodes:'
print_indexed (nds_xy)
print 

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

supps = [(0, (0, 0, 1)), (1, (0, 0, 0))]
print 'supports:', supps
print

chars = [(3e7, 0.3*0.3, 0.3^4/12), (3e7, 0.3*0.6, 0.3*0.6^3/12)]
print 'characteristics:'
print_list (chars, '   ')
print

echs = [0, 1, 0]
print 'elements <-> chars:', echs
print 

ndlds = [(3, (0., 0., 50.))]
print 'nodal loads:'
print_list (ndlds, '   ')
print

ellds = [(1, df2, (10., 25.)), (0, cf2, (3/5*75., 4/5*75., 0., 2/3))]
print 'element loads:'
print_list (ellds, '   ')

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (8.00000000000000, 0.000000000000000)
  2: (3.00000000000000, 4.00000000000000)
  3: (8.00000000000000, 4.00000000000000)

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

supports: [(0, (0, 0, 1)), (1, (0, 0, 0))]

characteristics:
   (3.00000000000000e7, 0.0900000000000000, 0.000675000000000000)
   (3.00000000000000e7, 0.180000000000000, 0.00540000000000000)

elements <-> chars: [0, 1, 0]

nodal loads:
   (3, (0.000000000000000, 0.000000000000000, 50.0000000000000))

element loads:
   (1, <function distributed_force2D at 0x7efa4e4a6b18>,
(10.0000000000000, 25.0000000000000))
   (0, <function concentrated_force2D at 0x7efa4e4a6c08>,
(45.0000000000000, 60.0000000000000, 0.000000000000000, 2/3))

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (8.00000000000000, 0.000000000000000)
  2: (3.00000000000000, 4.00000000000000)
  3: (8.00000000000000, 4.00000000000000)

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

supports: [(0, (0, 0, 1)), (1, (0, 0, 0))]

characteristics:
   (3.00000000000000e7, 0.0900000000000000, 0.000675000000000000)
   (3.00000000000000e7, 0.180000000000000, 0.00540000000000000)

elements <-> chars: [0, 1, 0]

nodal loads:
   (3, (0.000000000000000, 0.000000000000000, 50.0000000000000))

element loads:
   (1, <function distributed_force2D at 0x7efa4e4a6b18>, (10.0000000000000, 25.0000000000000))
   (0, <function concentrated_force2D at 0x7efa4e4a6c08>, (45.0000000000000, 60.0000000000000, 0.000000000000000, 2/3))

plot_frame (nds_xy, els, supps, nodal_loads = ndlds, elem_loads = ellds)

uu, ff = DiM (f2k, f2r, nds_xy, els, supps, chars, echs, nodal_loads = ndlds, elem_loads = ellds)

print_indexed (uu.apply_map (ND(6)))

  0: (0.000000, 0.000000, -0.0110758)
  1: (0.000000, 0.000000, 0.000000)
  2: (0.0220305, 0.0165630, 0.00287090)
  3: (0.0219583, 0.000156406, 0.00258386)

  0: (0.000000, 0.000000, -0.0110758)
  1: (0.000000, 0.000000, 0.000000)
  2: (0.0220305, 0.0165630, 0.00287090)
  3: (0.0219583, 0.000156406, 0.00258386)

print_indexed (ff.apply_map (ND(6)))

  0: (2.33707, -29.2604, -2.48690e-14, -47.3371, -30.7396, 46.3019)
  1: (52.9939, -19.4259, -46.3019, -102.994, -105.574, -169.069)
  2: (105.574, -102.994, 192.907, -105.574, 102.994, 219.069)

  0: (2.33707, -29.2604, -2.48690e-14, -47.3371, -30.7396, 46.3019)
  1: (52.9939, -19.4259, -46.3019, -102.994, -105.574, -169.069)
  2: (105.574, -102.994, 192.907, -105.574, 102.994, 219.069)

print_indexed (ff.apply_map (NDC(6)))

  0: (2.33707, -29.2604, 0.000000, -47.3371, -30.7396, 46.3019)
  1: (52.9939, -19.4259, -46.3019, -102.994, -105.574, -169.069)
  2: (105.574, -102.994, 192.907, -105.574, 102.994, 219.069)

  0: (2.33707, -29.2604, 0.000000, -47.3371, -30.7396, 46.3019)
  1: (52.9939, -19.4259, -46.3019, -102.994, -105.574, -169.069)
  2: (105.574, -102.994, 192.907, -105.574, 102.994, 219.069)

nds = coords_left_to_right (nds_xy)
print_indexed (nds)

  0: (0.000000000000000, -0.000000000000000)
  1: (8.00000000000000, -0.000000000000000)
  2: (3.00000000000000, -4.00000000000000)
  3: (8.00000000000000, -4.00000000000000)

  0: (0.000000000000000, -0.000000000000000)
  1: (8.00000000000000, -0.000000000000000)
  2: (3.00000000000000, -4.00000000000000)
  3: (8.00000000000000, -4.00000000000000)

ndof, tbl_dof = nodal_DsOF (len (nds), supps)
print ndof
print_indexed (tbl_dof)

7
  0: (-1, -1, 0)
  1: (-1, -1, -1)
  2: (1, 2, 3)
  3: (4, 5, 6)

7
  0: (-1, -1, 0)
  1: (-1, -1, -1)
  2: (1, 2, 3)
  3: (4, 5, 6)

K = structure_stiffness_matrix (f2k, f2r, ndof, tbl_dof, nds, els, echs, chars)
show (K.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrr} 16200.0 & 3888.00 & 2916.00 & 8100.00 & 0.000000 & 0.000000 & 0.000000 \\ 3888.00 & 1.27564 \times 10^{6} & -258267. & 3888.00 & -1.08000 \times 10^{6} & 0.000000 & 0.000000 \\ 2916.00 & -258267. & 361852. & -35964.0 & 0.000000 & -15552.0 & -38880.0 \\ 8100.00 & 3888.00 & -35964.0 & 145800. & 0.000000 & 38880.0 & 64800.0 \\ 0.000000 & -1.08000 \times 10^{6} & 0.000000 & 0.000000 & 1.08380 \times 10^{6} & 0.000000 & 7593.75 \\ 0.000000 & 0.000000 & -15552.0 & 38880.0 & 0.000000 & 690552. & 38880.0 \\ 0.000000 & 0.000000 & -38880.0 & 64800.0 & 7593.75 & 38880.0 & 149850. \end{array}\right)$

q = structure_force_vector (ndof, tbl_dof, nodal_loads = ndlds)
show (q.apply_map (ND(4)))

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

xx = K.solve_right (q)
xx

(-0.000139709765309643, 0.000508162731905217, 0.000402196975359017,
-0.000109289491824465, 0.000503141253890175, -0.0000108341675260316,
0.000462595135194311)

(-0.000139709765309643, 0.000508162731905217, 0.000402196975359017, -0.000109289491824465, 0.000503141253890175, -0.0000108341675260316, 0.000462595135194311)

uu = table_of_nodal_displacements (tbl_dof, xx)
print_indexed (uu.apply_map (ND(8)))

  0: (0.00000000, 0.00000000, -0.00013970977)
  1: (0.00000000, 0.00000000, 0.00000000)
  2: (0.00050816273, 0.00040219698, -0.00010928949)
  3: (0.00050314125, -0.000010834168, 0.00046259514)

  0: (0.00000000, 0.00000000, -0.00013970977)
  1: (0.00000000, 0.00000000, 0.00000000)
  2: (0.00050816273, 0.00040219698, -0.00010928949)
  3: (0.00050314125, -0.000010834168, 0.00046259514)

ff = all_end_forces (f2k, f2r, uu, nds, els, echs, chars)
print_indexed (ff.apply_map (NDC(6)))

  0: (9.10437, -0.0492808, 0.000000, -9.10437, 0.0492808, 0.246404)
  1: (5.42320, -7.31306, -0.246404, -5.42320, 7.31306, 36.8117)
  2: (-7.31306, -5.42320, 8.50450, 7.31306, 5.42320, 13.1883)

  0: (9.10437, -0.0492808, 0.000000, -9.10437, 0.0492808, 0.246404)
  1: (5.42320, -7.31306, -0.246404, -5.42320, 7.31306, 36.8117)
  2: (-7.31306, -5.42320, 8.50450, 7.31306, 5.42320, 13.1883)

ffn = ff

q = structure_force_vector (ndof, tbl_dof, elem_loads = ellds, elems = els, nodes = nds, fr = f2r)
show (q.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(-22.2222,\,78.5556,\,65.1667,\,-7.63889,\,25.0000,\,62.5000,\,52.0833\right)$

xx = K \ q
xx

(-0.0109361315963667, 0.0215223076981514, 0.0161607879460245,
0.00298018761208612, 0.0214551125721685, 0.000167240265102427,
0.00212126492365365)

(-0.0109361315963667, 0.0215223076981514, 0.0161607879460245, 0.00298018761208612, 0.0214551125721685, 0.000167240265102427, 0.00212126492365365)

uu = table_of_nodal_displacements (tbl_dof, xx)
print_indexed (uu.apply_map (ND(8)))

  0: (0.00000000, 0.00000000, -0.010936132)
  1: (0.00000000, 0.00000000, 0.00000000)
  2: (0.021522308, 0.016160788, 0.0029801876)
  3: (0.021455113, 0.00016724027, 0.0021212649)

  0: (0.00000000, 0.00000000, -0.010936132)
  1: (0.00000000, 0.00000000, 0.00000000)
  2: (0.021522308, 0.016160788, 0.0029801876)
  3: (0.021455113, 0.00016724027, 0.0021212649)

ff = all_end_forces (f2k, f2r, uu, nds, els, echs, chars, elem_loads = ellds)
print_indexed (ff.apply_map (NDC(6)))

  0: (-6.76730, -29.2111, 0.000000, -38.2327, -30.7889, 46.0555)
  1: (47.5707, -12.1128, -46.0555, -97.5707, -112.887, -205.880)
  2: (112.887, -97.5707, 184.403, -112.887, 97.5707, 205.880)

  0: (-6.76730, -29.2111, 0.000000, -38.2327, -30.7889, 46.0555)
  1: (47.5707, -12.1128, -46.0555, -97.5707, -112.887, -205.880)
  2: (112.887, -97.5707, 184.403, -112.887, 97.5707, 205.880)

ffe = ff

q = structure_force_vector (ndof, tbl_dof, nodal_loads = ndlds, 
                            elem_loads = ellds, elems = els, nodes = nds, fr = f2r)
show (q.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(-22.2222,\,78.5556,\,65.1667,\,-7.63889,\,25.0000,\,62.5000,\,102.083\right)$

xx = K \ q
uu = table_of_nodal_displacements (tbl_dof, xx)
print_indexed (uu.apply_map (ND(8)))

  0: (0.00000000, 0.00000000, -0.011075841)
  1: (0.00000000, 0.00000000, 0.00000000)
  2: (0.022030470, 0.016562985, 0.0028708981)
  3: (0.021958254, 0.00015640610, 0.0025838601)

  0: (0.00000000, 0.00000000, -0.011075841)
  1: (0.00000000, 0.00000000, 0.00000000)
  2: (0.022030470, 0.016562985, 0.0028708981)
  3: (0.021958254, 0.00015640610, 0.0025838601)

ff = all_end_forces (f2k, f2r, uu, nds, els, echs, chars, elem_loads = ellds)
print_indexed (ff.apply_map (NDC(6)))

  0: (2.33707, -29.2604, 0.000000, -47.3371, -30.7396, 46.3019)
  1: (52.9939, -19.4259, -46.3019, -102.994, -105.574, -169.069)
  2: (105.574, -102.994, 192.907, -105.574, 102.994, 219.069)

  0: (2.33707, -29.2604, 0.000000, -47.3371, -30.7396, 46.3019)
  1: (52.9939, -19.4259, -46.3019, -102.994, -105.574, -169.069)
  2: (105.574, -102.994, 192.907, -105.574, 102.994, 219.069)

print_indexed ((ffn + ffe).apply_map (NDC(6)))

  0: (2.33707, -29.2604, 0.000000, -47.3371, -30.7396, 46.3019)
  1: (52.9939, -19.4259, -46.3019, -102.994, -105.574, -169.069)
  2: (105.574, -102.994, 192.907, -105.574, 102.994, 219.069)

  0: (2.33707, -29.2604, 0.000000, -47.3371, -30.7396, 46.3019)
  1: (52.9939, -19.4259, -46.3019, -102.994, -105.574, -169.069)
  2: (105.574, -102.994, 192.907, -105.574, 102.994, 219.069)

ek2 = local_stiffness_matrix (2, f2k, nds, els, echs, chars) 
show (ek2.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 675000. & 0.000000 & 0.000000 & -675000. & 0.000000 & 0.000000 \\ 0.000000 & 3796.88 & -7593.75 & 0.000000 & -3796.88 & -7593.75 \\ 0.000000 & -7593.75 & 20250.0 & 0.000000 & 7593.75 & 10125.0 \\ -675000. & 0.000000 & 0.000000 & 675000. & 0.000000 & 0.000000 \\ 0.000000 & -3796.88 & 7593.75 & 0.000000 & 3796.88 & 7593.75 \\ 0.000000 & -7593.75 & 10125.0 & 0.000000 & 7593.75 & 20250.0 \end{array}\right)$

er2 = transform_global_to_local (2, f2r, nds, els)
show (er2.apply_map (ND(6)))

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

ek2g = er2.T * ek2 * er2
show (ek2g.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 3796.88 & 0.000000 & -7593.75 & -3796.88 & 0.000000 & -7593.75 \\ 0.000000 & 675000. & 0.000000 & 0.000000 & -675000. & 0.000000 \\ -7593.75 & 0.000000 & 20250.0 & 7593.75 & 0.000000 & 10125.0 \\ -3796.88 & 0.000000 & 7593.75 & 3796.88 & 0.000000 & 7593.75 \\ 0.000000 & -675000. & 0.000000 & 0.000000 & 675000. & 0.000000 \\ -7593.75 & 0.000000 & 10125.0 & 7593.75 & 0.000000 & 20250.0 \end{array}\right)$

ef0 = local_fixed_end_forces (0, ellds, nds, els)
ef0.apply_map (ND(6))

(-15.0000, -15.5556, 22.2222, -30.0000, -44.4444, -44.4444)

(-15.0000, -15.5556, 22.2222, -30.0000, -44.4444, -44.4444)

er0 = transform_global_to_local (0, f2r, nds, els)
ef0g = er0.T * ef0
ef0g.apply_map (ND(6))

(-21.4444, 2.66667, 22.2222, -53.5556, -2.66667, -44.4444)

(-21.4444, 2.66667, 22.2222, -53.5556, -2.66667, -44.4444)

ef1 = local_fixed_end_forces (1, ellds, nds, els)
ef1.apply_map (ND(6))

(-25.0000, -62.5000, 52.0833, -25.0000, -62.5000, -52.0833)

(-25.0000, -62.5000, 52.0833, -25.0000, -62.5000, -52.0833)

er1 = transform_global_to_local (1, f2r, nds, els)
ef1g = er1.T * ef1
ef1g.apply_map (ND(6))

(-25.0000, -62.5000, 52.0833, -25.0000, -62.5000, -52.0833)

(-25.0000, -62.5000, 52.0833, -25.0000, -62.5000, -52.0833)

ef2 = local_fixed_end_forces (2, ellds, nds, els)
ef2

(0, 0, 0, 0, 0, 0)

(0, 0, 0, 0, 0, 0)

nds_xy = [(0.0, 0.0), (8.3, 2.1), (0.0, 4.0), (8.3, 4.0), (0.0, 7.0), (4.6, 7.0), (8.3, 7.0), (10.6, 7.0)]
print 'nodes:'
print_indexed (nds_xy)
print 

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

supps = [(0, (0, 0, 1)), (1, (1, 0, 1))]
print 'supports:', supps
print

hinges = [5]
print 'hinges:', hinges
print

rels = [(2, (1, 1))]
print 'end releases:'
print_list (rels, '   ')
print

chars = [(3e7, A_rect (0.3, 0.3), I_rect (0.3, 0.3)), (3e7, A_rect (0.3, 0.6), I_rect (0.3, 0.6)), 
         (2e8, A_circ (0.05), I_circ (0.05))]
print 'characteristics:'
print_list (chars, '   ')
print

echs = [0, 0, 2, 0, 0, 1, 1, 1]
print 'elements <-> chars:'
print_indexed (echs)
print 

ndlds = [(2, (0., 0., -80.)), (5, (0., 80., 0.)), (7, (0., 0., 60.))]
print 'nodal loads:'
print_list (ndlds, '   ')
print 

ellds = [(5, df2, (0., 18.)), (6, df2, (0., 18.)), (6, cf2, (0., 75., 0., 1/2))]
print 'element loads:'
print_list (ellds, '   ')

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (8.30000000000000, 2.10000000000000)
  2: (0.000000000000000, 4.00000000000000)
  3: (8.30000000000000, 4.00000000000000)
  4: (0.000000000000000, 7.00000000000000)
  5: (4.60000000000000, 7.00000000000000)
  6: (8.30000000000000, 7.00000000000000)
  7: (10.6000000000000, 7.00000000000000)

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

supports: [(0, (0, 0, 1)), (1, (1, 0, 1))]

hinges: [5]

end releases:
   (2, (1, 1))

characteristics:
   (3.00000000000000e7, 0.0900000000000000, 0.000675000000000000)
   (3.00000000000000e7, 0.180000000000000, 0.00540000000000000)
   (2.00000000000000e8, 0.00196349540849362, 3.06796157577128e-7)

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

nodal loads:
   (2, (0.000000000000000, 0.000000000000000, -80.0000000000000))
   (5, (0.000000000000000, 80.0000000000000, 0.000000000000000))
   (7, (0.000000000000000, 0.000000000000000, 60.0000000000000))

element loads:
   (5, <function distributed_force2D at 0x663b758>, (0.000000000000000,
18.0000000000000))
   (6, <function distributed_force2D at 0x663b758>, (0.000000000000000,
18.0000000000000))
   (6, <function concentrated_force2D at 0x663b848>, (0.000000000000000,
75.0000000000000, 0.000000000000000, 1/2))

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (8.30000000000000, 2.10000000000000)
  2: (0.000000000000000, 4.00000000000000)
  3: (8.30000000000000, 4.00000000000000)
  4: (0.000000000000000, 7.00000000000000)
  5: (4.60000000000000, 7.00000000000000)
  6: (8.30000000000000, 7.00000000000000)
  7: (10.6000000000000, 7.00000000000000)

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

supports: [(0, (0, 0, 1)), (1, (1, 0, 1))]

hinges: [5]

end releases:
   (2, (1, 1))

characteristics:
   (3.00000000000000e7, 0.0900000000000000, 0.000675000000000000)
   (3.00000000000000e7, 0.180000000000000, 0.00540000000000000)
   (2.00000000000000e8, 0.00196349540849362, 3.06796157577128e-7)

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

nodal loads:
   (2, (0.000000000000000, 0.000000000000000, -80.0000000000000))
   (5, (0.000000000000000, 80.0000000000000, 0.000000000000000))
   (7, (0.000000000000000, 0.000000000000000, 60.0000000000000))

element loads:
   (5, <function distributed_force2D at 0x663b758>, (0.000000000000000, 18.0000000000000))
   (6, <function distributed_force2D at 0x663b758>, (0.000000000000000, 18.0000000000000))
   (6, <function concentrated_force2D at 0x663b848>, (0.000000000000000, 75.0000000000000, 0.000000000000000, 1/2))

plot_frame (nds_xy, els, supps, hinges, rels, ndlds, ellds)

uu, ff = DiM (f2k, f2r, nds_xy, els, supps, chars, echs, hinges, rels, ndlds, ellds)

print_indexed (ff.apply_map (NDC(6)))

  0: (124.670, 0.000000, 0.000000, -124.670, 0.000000, 0.000000)
  1: (124.670, 154.347, -80.0000, -124.670, -154.347, -383.041)
  2: (-154.347, 0.000000, 0.000000, 154.347, 0.000000, 0.000000)
  3: (179.730, 0.000000, 0.000000, -179.730, 0.000000, 0.000000)
  4: (179.730, -154.347, 0.000000, -179.730, 154.347, 463.041)
  5: (154.347, -124.670, 383.041, -154.347, 41.8699, 0.000000)
  6: (154.347, 38.1301, 0.000000, -154.347, -179.730, -403.041)
  7: (0.000000, 0.000000, -60.0000, 0.000000, 0.000000, 60.0000)

  0: (124.670, 0.000000, 0.000000, -124.670, 0.000000, 0.000000)
  1: (124.670, 154.347, -80.0000, -124.670, -154.347, -383.041)
  2: (-154.347, 0.000000, 0.000000, 154.347, 0.000000, 0.000000)
  3: (179.730, 0.000000, 0.000000, -179.730, 0.000000, 0.000000)
  4: (179.730, -154.347, 0.000000, -179.730, 154.347, 463.041)
  5: (154.347, -124.670, 383.041, -154.347, 41.8699, 0.000000)
  6: (154.347, 38.1301, 0.000000, -154.347, -179.730, -403.041)
  7: (0.000000, 0.000000, -60.0000, 0.000000, 0.000000, 60.0000)

nds_xy = [(0.0, 0.0), (8.3, 2.1), (0.0, 4.0), (8.3, 4.0), (0.0, 7.0), (4.6, 7.0), (8.3, 7.0), (10.6, 7.0), (6.45, 7.0)]
print 'nodes:'
print_indexed (nds_xy)
print 

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

echs = [0, 0, 2, 0, 0, 1, 1, 1, 1]

ellds = [(5, df2, (0., 18.)), (6, df2, (0., 18.)), (7, df2, (0., 18.)), (6, cf2, (0., 75., 0., 1))]
print 'element loads:'
print_list (ellds, '   ')

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (8.30000000000000, 2.10000000000000)
  2: (0.000000000000000, 4.00000000000000)
  3: (8.30000000000000, 4.00000000000000)
  4: (0.000000000000000, 7.00000000000000)
  5: (4.60000000000000, 7.00000000000000)
  6: (8.30000000000000, 7.00000000000000)
  7: (10.6000000000000, 7.00000000000000)
  8: (6.45000000000000, 7.00000000000000)

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

element loads:
   (5, <function distributed_force2D at 0x663b758>, (0.000000000000000,
18.0000000000000))
   (6, <function distributed_force2D at 0x663b758>, (0.000000000000000,
18.0000000000000))
   (7, <function distributed_force2D at 0x663b758>, (0.000000000000000,
18.0000000000000))
   (6, <function concentrated_force2D at 0x663b848>, (0.000000000000000,
75.0000000000000, 0.000000000000000, 1))

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (8.30000000000000, 2.10000000000000)
  2: (0.000000000000000, 4.00000000000000)
  3: (8.30000000000000, 4.00000000000000)
  4: (0.000000000000000, 7.00000000000000)
  5: (4.60000000000000, 7.00000000000000)
  6: (8.30000000000000, 7.00000000000000)
  7: (10.6000000000000, 7.00000000000000)
  8: (6.45000000000000, 7.00000000000000)

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

element loads:
   (5, <function distributed_force2D at 0x663b758>, (0.000000000000000, 18.0000000000000))
   (6, <function distributed_force2D at 0x663b758>, (0.000000000000000, 18.0000000000000))
   (7, <function distributed_force2D at 0x663b758>, (0.000000000000000, 18.0000000000000))
   (6, <function concentrated_force2D at 0x663b848>, (0.000000000000000, 75.0000000000000, 0.000000000000000, 1))

plot_frame (nds_xy, els, supps, hinges, rels, ndlds, ellds)

uu, ff = DiM (f2k, f2r, nds_xy, els, supps, chars, echs, hinges, rels, ndlds, ellds)

print_indexed (ff.apply_map (NDC(6)))

  0: (124.670, 0.000000, 0.000000, -124.670, 0.000000, 0.000000)
  1: (124.670, 154.347, -80.0000, -124.670, -154.347, -383.041)
  2: (-154.347, 0.000000, 0.000000, 154.347, 0.000000, 0.000000)
  3: (179.730, 0.000000, 0.000000, -179.730, 0.000000, 0.000000)
  4: (179.730, -154.347, 0.000000, -179.730, 154.347, 463.041)
  5: (154.347, -124.670, 383.041, -154.347, 41.8699, 0.000000)
  6: (154.347, 38.1301, 0.000000, -154.347, -146.430, -101.343)
  7: (154.347, 146.430, 101.343, -154.347, -179.730, -403.041)
  8: (0.000000, 0.000000, -60.0000, 0.000000, 0.000000, 60.0000)

  0: (124.670, 0.000000, 0.000000, -124.670, 0.000000, 0.000000)
  1: (124.670, 154.347, -80.0000, -124.670, -154.347, -383.041)
  2: (-154.347, 0.000000, 0.000000, 154.347, 0.000000, 0.000000)
  3: (179.730, 0.000000, 0.000000, -179.730, 0.000000, 0.000000)
  4: (179.730, -154.347, 0.000000, -179.730, 154.347, 463.041)
  5: (154.347, -124.670, 383.041, -154.347, 41.8699, 0.000000)
  6: (154.347, 38.1301, 0.000000, -154.347, -146.430, -101.343)
  7: (154.347, 146.430, 101.343, -154.347, -179.730, -403.041)
  8: (0.000000, 0.000000, -60.0000, 0.000000, 0.000000, 60.0000)

nds = coords_left_to_right (nds_xy)
print_list (nds)

(0.000000000000000, -0.000000000000000)
(8.30000000000000, -2.10000000000000)
(0.000000000000000, -4.00000000000000)
(8.30000000000000, -4.00000000000000)
(0.000000000000000, -7.00000000000000)
(4.60000000000000, -7.00000000000000)
(8.30000000000000, -7.00000000000000)
(10.6000000000000, -7.00000000000000)
(6.45000000000000, -7.00000000000000)

(0.000000000000000, -0.000000000000000)
(8.30000000000000, -2.10000000000000)
(0.000000000000000, -4.00000000000000)
(8.30000000000000, -4.00000000000000)
(0.000000000000000, -7.00000000000000)
(4.60000000000000, -7.00000000000000)
(8.30000000000000, -7.00000000000000)
(10.6000000000000, -7.00000000000000)
(6.45000000000000, -7.00000000000000)

ne_inc = node_element_incidences (nds, els)
print_indexed (ne_inc)

  0: [0]
  1: [3]
  2: [0, 1, 2]
  3: [2, 3, 4]
  4: [1, 5]
  5: [5, 6]
  6: [4, 7, 8]
  7: [8]
  8: [6, 7]

  0: [0]
  1: [3]
  2: [0, 1, 2]
  3: [2, 3, 4]
  4: [1, 5]
  5: [5, 6]
  6: [4, 7, 8]
  7: [8]
  8: [6, 7]

endc = end_conditions (els, rels, hinges, ne_inc)
print_indexed (endc)

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

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

ndof, tbl_dof = nodal_DsOF (len (nds), supps, hinges)
print ndof
print_indexed (tbl_dof)

23
  0: (-1, -1, 0)
  1: (1, -1, 2)
  2: (3, 4, 5)
  3: (6, 7, 8)
  4: (9, 10, 11)
  5: (12, 13, -1)
  6: (14, 15, 16)
  7: (17, 18, 19)
  8: (20, 21, 22)

23
  0: (-1, -1, 0)
  1: (1, -1, 2)
  2: (3, 4, 5)
  3: (6, 7, 8)
  4: (9, 10, 11)
  5: (12, 13, -1)
  6: (14, 15, 16)
  7: (17, 18, 19)
  8: (20, 21, 22)

K, cc = structure_stiffness_matrix (f2k, f2r, ndof, tbl_dof, nds, els, echs, chars, endc)
show (K.apply_map (ND(4)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrrrrrrrrrrrrrrrrr} 20250. & 0.0000 & 0.0000 & 7594. & 0.0000 & 10120. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 35430. & -33660. & 0.0000 & 0.0000 & 0.0000 & -35430. & 0.0000 & -33660. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & -33660. & 42630. & 0.0000 & 0.0000 & 0.0000 & 33660. & 0.0000 & 21320. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 7594. & 0.0000 & 0.0000 & 60110. & 0.0000 & -5906. & -47310. & 0.0000 & 0.0000 & -9000. & 0.0000 & -13500. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 1.575 \times 10^{6} & -8.882 \times 10^{-16} & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -900000. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 10120. & 0.0000 & 0.0000 & -5906. & 0.0000 & 47250. & 0.0000 & 0.0000 & 0.0000 & 13500. & 0.0000 & 13500. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & -35430. & 33660. & -47310. & 0.0000 & 0.0000 & 91740. & 0.0000 & 20160. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -9000. & 0.0000 & -13500. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 8.882 \times 10^{-16} & 0.0000 & 2.321 \times 10^{6} & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -900000. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & -33660. & 21320. & 0.0000 & 0.0000 & 0.0000 & 20160. & 0.0000 & 69630. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 13500. & 0.0000 & 13500. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & -9000. & 0.0000 & 13500. & 0.0000 & 0.0000 & 0.0000 & 1.183 \times 10^{6} & 0.0000 & 13500. & -1.174 \times 10^{6} & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & -900000. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 905000. & -22970. & 0.0000 & -4993. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & -13500. & 0.0000 & 13500. & 0.0000 & 0.0000 & 0.0000 & 13500. & -22970. & 132700. & 0.0000 & 22970. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -1.174 \times 10^{6} & 0.0000 & 0.0000 & 4.093 \times 10^{6} & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -2.919 \times 10^{6} & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -4993. & 22970. & 0.0000 & 81750. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -76760. & -142000. \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -9000. & 0.0000 & 13500. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 5.276 \times 10^{6} & 0.0000 & 13500. & -2.348 \times 10^{6} & 0.0000 & 0.0000 & -2.919 \times 10^{6} & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -900000. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 1.367 \times 10^{6} & 100300. & 0.0000 & -159800. & -183700. & 0.0000 & -307000. & 284000. \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -13500. & 0.0000 & 13500. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 13500. & 100300. & 659000. & 0.0000 & 183700. & 140900. & 0.0000 & -284000. & 175100. \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -2.348 \times 10^{6} & 0.0000 & 0.0000 & 2.348 \times 10^{6} & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -159800. & 183700. & 0.0000 & 159800. & 183700. & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -183700. & 140900. & 0.0000 & 183700. & 281700. & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -2.919 \times 10^{6} & 0.0000 & -2.919 \times 10^{6} & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 5.838 \times 10^{6} & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -76760. & 0.0000 & -307000. & -284000. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 383800. & -142000. \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -142000. & 0.0000 & 284000. & 175100. & 0.0000 & 0.0000 & 0.0000 & 0.0000 & -142000. & 613000. \end{array}\right)$

q = structure_force_vector (ndof, tbl_dof, ndlds, ellds, els, nds, f2r, cc)
show (q.apply_map (ND(4)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.0000,\,0.0000,\,0.0000,\,0.0000,\,0.0000,\,-80.00,\,0.0000,\,0.0000,\,0.0000,\,0.0000,\,51.75,\,-47.61,\,0.0000,\,123.5,\,0.0000,\,16.65,\,5.134,\,0.0000,\,0.0000,\,60.00,\,0.0000,\,112.5,\,2.567\right)$

uu = nodal_displacements (K, q, tbl_dof)
print_indexed (uu.apply_map (ND(8)))

  0: (0.00000000, 0.00000000, 0.0048548963)
  1: (-0.036903300, 0.00000000, -0.010918927)
  2: (-0.019419585, 0.00018469612, 0.0048548963)
  3: (-0.016157338, 0.00012647675, -0.010918927)
  4: (-0.017462685, 0.00032321821, -0.017592618)
  5: (-0.017594166, 0.095853643, 0.00000000)
  6: (-0.017699923, 0.00032617689, 0.023380439)
  7: (-0.017699923, -0.054428462, 0.024232291)
  8: (-0.017647044, 0.046720889, 0.026201787)

  0: (0.00000000, 0.00000000, 0.0048548963)
  1: (-0.036903300, 0.00000000, -0.010918927)
  2: (-0.019419585, 0.00018469612, 0.0048548963)
  3: (-0.016157338, 0.00012647675, -0.010918927)
  4: (-0.017462685, 0.00032321821, -0.017592618)
  5: (-0.017594166, 0.095853643, 0.00000000)
  6: (-0.017699923, 0.00032617689, 0.023380439)
  7: (-0.017699923, -0.054428462, 0.024232291)
  8: (-0.017647044, 0.046720889, 0.026201787)

ff = all_end_forces (f2k, f2r, uu, nds, els, echs, chars, endc)
print_indexed (ff.apply_map (NDC(6)))

  0: (124.670, 0.000000, 0.000000, -124.670, 0.000000, 0.000000)
  1: (124.670, 154.347, -80.0000, -124.670, -154.347, -383.041)
  2: (-154.347, 0.000000, 0.000000, 154.347, 0.000000, 0.000000)
  3: (179.730, 0.000000, 0.000000, -179.730, 0.000000, 0.000000)
  4: (179.730, -154.347, 0.000000, -179.730, 154.347, 463.041)
  5: (154.347, -72.9199, 335.431, -154.347, 72.9199, 0.000000)
  6: (154.347, 50.6176, 0.000000, -154.347, -50.6176, -93.6426)
  7: (154.347, 163.080, 96.2095, -154.347, -163.080, -397.908)
  8: (0.000000, 0.000000, -60.0000, 0.000000, 0.000000, 60.0000)

  0: (124.670, 0.000000, 0.000000, -124.670, 0.000000, 0.000000)
  1: (124.670, 154.347, -80.0000, -124.670, -154.347, -383.041)
  2: (-154.347, 0.000000, 0.000000, 154.347, 0.000000, 0.000000)
  3: (179.730, 0.000000, 0.000000, -179.730, 0.000000, 0.000000)
  4: (179.730, -154.347, 0.000000, -179.730, 154.347, 463.041)
  5: (154.347, -72.9199, 335.431, -154.347, 72.9199, 0.000000)
  6: (154.347, 50.6176, 0.000000, -154.347, -50.6176, -93.6426)
  7: (154.347, 163.080, 96.2095, -154.347, -163.080, -397.908)
  8: (0.000000, 0.000000, -60.0000, 0.000000, 0.000000, 60.0000)

ek6 = local_stiffness_matrix (6, f2k, nds, els, echs, chars) 
show (ek6.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 2.91892 \times 10^{6} & 0.000000 & 0.000000 & -2.91892 \times 10^{6} & 0.000000 & 0.000000 \\ 0.000000 & 307030. & -284003. & 0.000000 & -307030. & -284003. \\ 0.000000 & -284003. & 350270. & 0.000000 & 284003. & 175135. \\ -2.91892 \times 10^{6} & 0.000000 & 0.000000 & 2.91892 \times 10^{6} & 0.000000 & 0.000000 \\ 0.000000 & -307030. & 284003. & 0.000000 & 307030. & 284003. \\ 0.000000 & -284003. & 175135. & 0.000000 & 284003. & 350270. \end{array}\right)$

cr6 = condense_stiffness_matrix (ek6, 2)
show (ek6.apply_map (NDC(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 2.91892 \times 10^{6} & 0.000000 & 0.000000 & -2.91892 \times 10^{6} & 0.000000 & 0.000000 \\ 0.000000 & 76757.5 & 0.000000 & 0.000000 & -76757.5 & -142001. \\ 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ -2.91892 \times 10^{6} & 0.000000 & 0.000000 & 2.91892 \times 10^{6} & 0.000000 & 0.000000 \\ 0.000000 & -76757.5 & 0.000000 & 0.000000 & 76757.5 & 142001. \\ 0.000000 & -142001. & 0.000000 & 0.000000 & 142001. & 262703. \end{array}\right)$

ef6 = local_fixed_end_forces (6, ellds, nds, els)
ef6.apply_map (ND(6))

(0.000000, -16.6500, 5.13375, 0.000000, -91.6500, -5.13375)

(0.000000, -16.6500, 5.13375, 0.000000, -91.6500, -5.13375)

condense_force_vector (ef6, 2, cr6) 
ef6.apply_map (ND(6))

(0.000000, -12.4875, 0.000000, 0.000000, -95.8125, -7.70062)

(0.000000, -12.4875, 0.000000, 0.000000, -95.8125, -7.70062)

ek2 = local_stiffness_matrix (2, f2k, nds, els, echs, chars) 
show (ek2.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 47313.1 & 0.000000 & 0.000000 & -47313.1 & 0.000000 & 0.000000 \\ 0.000000 & 1.28774 & -5.34410 & 0.000000 & -1.28774 & -5.34410 \\ 0.000000 & -5.34410 & 29.5707 & 0.000000 & 5.34410 & 14.7854 \\ -47313.1 & 0.000000 & 0.000000 & 47313.1 & 0.000000 & 0.000000 \\ 0.000000 & -1.28774 & 5.34410 & 0.000000 & 1.28774 & 5.34410 \\ 0.000000 & -5.34410 & 14.7854 & 0.000000 & 5.34410 & 29.5707 \end{array}\right)$

condense_stiffness_matrix (ek2, 2)
show (ek2.apply_map (NDC(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 47313.1 & 0.000000 & 0.000000 & -47313.1 & 0.000000 & 0.000000 \\ 0.000000 & 0.321934 & 0.000000 & 0.000000 & -0.321934 & -2.67205 \\ 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 & 0.000000 \\ -47313.1 & 0.000000 & 0.000000 & 47313.1 & 0.000000 & 0.000000 \\ 0.000000 & -0.321934 & 0.000000 & 0.000000 & 0.321934 & 2.67205 \\ 0.000000 & -2.67205 & 0.000000 & 0.000000 & 2.67205 & 22.1780 \end{array}\right)$

condense_stiffness_matrix (ek2, 5)
show (ek2.apply_map (NDC(6)))

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

nds_xy = [(0., 0.), (3., 0.), (6., 0.), (1.5, 3.), (4.5, 3.)]
print "joints:"
print_indexed (nds_xy)
print

els = [(0,1), (1,2), (0,3), (3,1), (1,4), (4,2), (3,4)]
print "bars:"
print_indexed (els)
print

supps = [(0, (0, 0)), (2, (1, 0))]
print "supports:", supps
print

chars = [(2*1e8, 0.001)]
print 'characteristics:', chars
print

echs = [0, 0, 0, 0, 0, 0, 0]

ndlds = [(1, (0., 100.)), (4, (75., 0.))]
print "loads:" 
print_list (ndlds, '   ')
print

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

bars:
  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: [(2.00000000000000e8, 0.00100000000000000)]

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

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

bars:
  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: [(2.00000000000000e8, 0.00100000000000000)]

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

plot_truss (nds_xy, els, supps, ndlds)

uu, ft = DiM (t2k, t2r, nds_xy, els, supps, chars, echs, nodal_loads = ndlds)

print_indexed (ft.apply_map (NDC(6)))

  0: (-81.2500, 81.2500)
  1: (-43.7500, 43.7500)
  2: (13.9754, -13.9754)
  3: (-13.9754, 13.9754)
  4: (-97.8280, 97.8280)
  5: (97.8280, -97.8280)
  6: (12.5000, -12.5000)

  0: (-81.2500, 81.2500)
  1: (-43.7500, 43.7500)
  2: (13.9754, -13.9754)
  3: (-13.9754, 13.9754)
  4: (-97.8280, 97.8280)
  5: (97.8280, -97.8280)
  6: (12.5000, -12.5000)

nds = coords_left_to_right (nds_xy)

ndof, tbl_dof = nodal_DsOF (len (nds), supps)
print ndof
print_indexed (tbl_dof)

7
  0: (-1, -1)
  1: (0, 1)
  2: (2, -1)
  3: (3, 4)
  4: (5, 6)

7
  0: (-1, -1)
  1: (0, 1)
  2: (2, -1)
  3: (3, 4)
  4: (5, 6)

ek4 = local_stiffness_matrix (4, t2k, nds, els, echs, chars) 
show (ek4.apply_map (ND(6)))

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

er4 = transform_global_to_local (4, t2r, nds, els)
show (er4.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 0.447214 & -0.894427 & 0.000000 & 0.000000 \\ 0.000000 & 0.000000 & 0.447214 & -0.894427 \end{array}\right)$

ek4g = er4.T * ek4 * er4
show (ek4g.apply_map (ND(6)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 11925.7 & -23851.4 & -11925.7 & 23851.4 \\ -23851.4 & 47702.8 & 23851.4 & -47702.8 \\ -11925.7 & 23851.4 & 11925.7 & -23851.4 \\ 23851.4 & -47702.8 & -23851.4 & 47702.8 \end{array}\right)$

print "joints:"
print_indexed (nds_xy)
print

print "bars:"
print_indexed (els)
print

supps = [(0, (0, 0, 1)), (2, (1, 0, 1))]
print "supports:", supps
print

chars = [(2*1e8, 0.001, 1.2e-6)]
print 'characteristics:', chars
print

ndlds = [(1, (0., 100., 0.)), (4, (75., 0., 0.))]
print "loads:" 
print_list (ndlds, '   ')
print

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

bars:
  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, 1)), (2, (1, 0, 1))]

characteristics: [(2.00000000000000e8, 0.00100000000000000,
1.20000000000000e-6)]

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

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

bars:
  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, 1)), (2, (1, 0, 1))]

characteristics: [(2.00000000000000e8, 0.00100000000000000, 1.20000000000000e-6)]

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

plot_frame (nds_xy, els, supps, nodal_loads = ndlds)

uu, ff = DiM (f2k, f2r, nds_xy, els, supps, chars, echs, nodal_loads = ndlds)

print_indexed (ff.apply_map (ND(6)))

  0: (-81.1876, -0.0644775, 0.00972327, 81.1876, 0.0644775, 0.183709)
  1: (-43.7521, 0.182171, -0.345888, 43.7521, -0.182171, -0.200625)
  2: (13.8898, -0.0269820, -0.00972327, -13.8898, 0.0269820, 0.100224)
  3: (-13.9853, -0.0623117, 0.0816645, 13.9853, 0.0623117, 0.127335)
  4: (-97.5227, -0.0231945, 0.0348437, 97.5227, 0.0231945, 0.0429529)
  5: (97.6660, -0.0833106, 0.0788072, -97.6660, 0.0833106, 0.200625)
  6: (12.3863, 0.101216, -0.181888, -12.3863, -0.101216, -0.121760)

  0: (-81.1876, -0.0644775, 0.00972327, 81.1876, 0.0644775, 0.183709)
  1: (-43.7521, 0.182171, -0.345888, 43.7521, -0.182171, -0.200625)
  2: (13.8898, -0.0269820, -0.00972327, -13.8898, 0.0269820, 0.100224)
  3: (-13.9853, -0.0623117, 0.0816645, 13.9853, 0.0623117, 0.127335)
  4: (-97.5227, -0.0231945, 0.0348437, 97.5227, 0.0231945, 0.0429529)
  5: (97.6660, -0.0833106, 0.0788072, -97.6660, 0.0833106, 0.200625)
  6: (12.3863, 0.101216, -0.181888, -12.3863, -0.101216, -0.121760)

print_indexed (ft.apply_map (ND(6)))

  0: (-81.2500, 81.2500)
  1: (-43.7500, 43.7500)
  2: (13.9754, -13.9754)
  3: (-13.9754, 13.9754)
  4: (-97.8280, 97.8280)
  5: (97.8280, -97.8280)
  6: (12.5000, -12.5000)

  0: (-81.2500, 81.2500)
  1: (-43.7500, 43.7500)
  2: (13.9754, -13.9754)
  3: (-13.9754, 13.9754)
  4: (-97.8280, 97.8280)
  5: (97.8280, -97.8280)
  6: (12.5000, -12.5000)

ndof, tbl_dof = nodal_DsOF (len (nds), supps)
print ndof
print_indexed (tbl_dof)

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

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

nds_xy = make_nodes (17, (0., 0.), (0., 3.25)) 
nds_xy += make_nodes (17, (5.5, 0.), (0., 3.25)) 
print_indexed (nds_xy)

  0: [0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 3.25000000000000]
  2: [0.000000000000000, 6.50000000000000]
  3: [0.000000000000000, 9.75000000000000]
  4: [0.000000000000000, 13.0000000000000]
  5: [0.000000000000000, 16.2500000000000]
  6: [0.000000000000000, 19.5000000000000]
  7: [0.000000000000000, 22.7500000000000]
  8: [0.000000000000000, 26.0000000000000]
  9: [0.000000000000000, 29.2500000000000]
 10: [0.000000000000000, 32.5000000000000]
 11: [0.000000000000000, 35.7500000000000]
 12: [0.000000000000000, 39.0000000000000]
 13: [0.000000000000000, 42.2500000000000]
 14: [0.000000000000000, 45.5000000000000]
 15: [0.000000000000000, 48.7500000000000]
 16: [0.000000000000000, 52.0000000000000]
 17: [5.50000000000000, 0.000000000000000]
 18: [5.50000000000000, 3.25000000000000]
 19: [5.50000000000000, 6.50000000000000]
 20: [5.50000000000000, 9.75000000000000]
 21: [5.50000000000000, 13.0000000000000]
 22: [5.50000000000000, 16.2500000000000]
 23: [5.50000000000000, 19.5000000000000]
 24: [5.50000000000000, 22.7500000000000]
 25: [5.50000000000000, 26.0000000000000]
 26: [5.50000000000000, 29.2500000000000]
 27: [5.50000000000000, 32.5000000000000]
 28: [5.50000000000000, 35.7500000000000]
 29: [5.50000000000000, 39.0000000000000]
 30: [5.50000000000000, 42.2500000000000]
 31: [5.50000000000000, 45.5000000000000]
 32: [5.50000000000000, 48.7500000000000]
 33: [5.50000000000000, 52.0000000000000]

  0: [0.000000000000000, 0.000000000000000]
  1: [0.000000000000000, 3.25000000000000]
  2: [0.000000000000000, 6.50000000000000]
  3: [0.000000000000000, 9.75000000000000]
  4: [0.000000000000000, 13.0000000000000]
  5: [0.000000000000000, 16.2500000000000]
  6: [0.000000000000000, 19.5000000000000]
  7: [0.000000000000000, 22.7500000000000]
  8: [0.000000000000000, 26.0000000000000]
  9: [0.000000000000000, 29.2500000000000]
 10: [0.000000000000000, 32.5000000000000]
 11: [0.000000000000000, 35.7500000000000]
 12: [0.000000000000000, 39.0000000000000]
 13: [0.000000000000000, 42.2500000000000]
 14: [0.000000000000000, 45.5000000000000]
 15: [0.000000000000000, 48.7500000000000]
 16: [0.000000000000000, 52.0000000000000]
 17: [5.50000000000000, 0.000000000000000]
 18: [5.50000000000000, 3.25000000000000]
 19: [5.50000000000000, 6.50000000000000]
 20: [5.50000000000000, 9.75000000000000]
 21: [5.50000000000000, 13.0000000000000]
 22: [5.50000000000000, 16.2500000000000]
 23: [5.50000000000000, 19.5000000000000]
 24: [5.50000000000000, 22.7500000000000]
 25: [5.50000000000000, 26.0000000000000]
 26: [5.50000000000000, 29.2500000000000]
 27: [5.50000000000000, 32.5000000000000]
 28: [5.50000000000000, 35.7500000000000]
 29: [5.50000000000000, 39.0000000000000]
 30: [5.50000000000000, 42.2500000000000]
 31: [5.50000000000000, 45.5000000000000]
 32: [5.50000000000000, 48.7500000000000]
 33: [5.50000000000000, 52.0000000000000]

els = make_elements (16, (0, 1), (1, 1)) 
els += make_elements (16, (17, 18), (1, 1)) 
els += make_elements (16, (1, 18), (1, 1)) 
print_indexed (els)

  0: [0, 1]
  1: [1, 2]
  2: [2, 3]
  3: [3, 4]
  4: [4, 5]
  5: [5, 6]
  6: [6, 7]
  7: [7, 8]
  8: [8, 9]
  9: [9, 10]
 10: [10, 11]
 11: [11, 12]
 12: [12, 13]
 13: [13, 14]
 14: [14, 15]
 15: [15, 16]
 16: [17, 18]
 17: [18, 19]
 18: [19, 20]
 19: [20, 21]
 20: [21, 22]
 21: [22, 23]
 22: [23, 24]
 23: [24, 25]
 24: [25, 26]
 25: [26, 27]
 26: [27, 28]
 27: [28, 29]
 28: [29, 30]
 29: [30, 31]
 30: [31, 32]
 31: [32, 33]
 32: [1, 18]
 33: [2, 19]
 34: [3, 20]
 35: [4, 21]
 36: [5, 22]
 37: [6, 23]
 38: [7, 24]
 39: [8, 25]
 40: [9, 26]
 41: [10, 27]
 42: [11, 28]
 43: [12, 29]
 44: [13, 30]
 45: [14, 31]
 46: [15, 32]
 47: [16, 33]

  0: [0, 1]
  1: [1, 2]
  2: [2, 3]
  3: [3, 4]
  4: [4, 5]
  5: [5, 6]
  6: [6, 7]
  7: [7, 8]
  8: [8, 9]
  9: [9, 10]
 10: [10, 11]
 11: [11, 12]
 12: [12, 13]
 13: [13, 14]
 14: [14, 15]
 15: [15, 16]
 16: [17, 18]
 17: [18, 19]
 18: [19, 20]
 19: [20, 21]
 20: [21, 22]
 21: [22, 23]
 22: [23, 24]
 23: [24, 25]
 24: [25, 26]
 25: [26, 27]
 26: [27, 28]
 27: [28, 29]
 28: [29, 30]
 29: [30, 31]
 30: [31, 32]
 31: [32, 33]
 32: [1, 18]
 33: [2, 19]
 34: [3, 20]
 35: [4, 21]
 36: [5, 22]
 37: [6, 23]
 38: [7, 24]
 39: [8, 25]
 40: [9, 26]
 41: [10, 27]
 42: [11, 28]
 43: [12, 29]
 44: [13, 30]
 45: [14, 31]
 46: [15, 32]
 47: [16, 33]

stup1_e = range (0, 16) 
stup1d_e = range (0, 8) 
stup1g_e = range (8, 16) 
stup2_e = range (16, 32) 
stup2d_e = range (16, 24) 
stup2g_e = range (24, 32) 
grede_e = range (32, 48) 
stup1_n = range (0, 17)
stup2_n = range (17, 34)

supps = [(0, (0, 0, 0)), (17, (0, 0, 0))]

ch_stupd = (3e7, 0.6*0.8, 0.6*0.8^3/12) 
ch_stupg = (3e7, 0.4*0.6, 0.4*0.6^3/12) 
ch_greda = (3e7, 0.4*0.8, 0.4*0.8^3/12) 
chars = [ch_stupd, ch_stupg, ch_greda]

echs = [0 for i in stup1d_e] + [1 for i in stup1g_e] 
echs += [0 for i in stup2d_e] + [1 for i in stup2g_e] 
echs += [2 for i in grede_e] 
print_indexed (echs)

ndlds = [(i, (45./16*i, 0., 0.)) for i in stup1_n[1:]] 
print_list (ndlds)

(1, (2.81250000000000, 0.000000000000000, 0.000000000000000))
(2, (5.62500000000000, 0.000000000000000, 0.000000000000000))
(3, (8.43750000000000, 0.000000000000000, 0.000000000000000))
(4, (11.2500000000000, 0.000000000000000, 0.000000000000000))
(5, (14.0625000000000, 0.000000000000000, 0.000000000000000))
(6, (16.8750000000000, 0.000000000000000, 0.000000000000000))
(7, (19.6875000000000, 0.000000000000000, 0.000000000000000))
(8, (22.5000000000000, 0.000000000000000, 0.000000000000000))
(9, (25.3125000000000, 0.000000000000000, 0.000000000000000))
(10, (28.1250000000000, 0.000000000000000, 0.000000000000000))
(11, (30.9375000000000, 0.000000000000000, 0.000000000000000))
(12, (33.7500000000000, 0.000000000000000, 0.000000000000000))
(13, (36.5625000000000, 0.000000000000000, 0.000000000000000))
(14, (39.3750000000000, 0.000000000000000, 0.000000000000000))
(15, (42.1875000000000, 0.000000000000000, 0.000000000000000))
(16, (45.0000000000000, 0.000000000000000, 0.000000000000000))

(1, (2.81250000000000, 0.000000000000000, 0.000000000000000))
(2, (5.62500000000000, 0.000000000000000, 0.000000000000000))
(3, (8.43750000000000, 0.000000000000000, 0.000000000000000))
(4, (11.2500000000000, 0.000000000000000, 0.000000000000000))
(5, (14.0625000000000, 0.000000000000000, 0.000000000000000))
(6, (16.8750000000000, 0.000000000000000, 0.000000000000000))
(7, (19.6875000000000, 0.000000000000000, 0.000000000000000))
(8, (22.5000000000000, 0.000000000000000, 0.000000000000000))
(9, (25.3125000000000, 0.000000000000000, 0.000000000000000))
(10, (28.1250000000000, 0.000000000000000, 0.000000000000000))
(11, (30.9375000000000, 0.000000000000000, 0.000000000000000))
(12, (33.7500000000000, 0.000000000000000, 0.000000000000000))
(13, (36.5625000000000, 0.000000000000000, 0.000000000000000))
(14, (39.3750000000000, 0.000000000000000, 0.000000000000000))
(15, (42.1875000000000, 0.000000000000000, 0.000000000000000))
(16, (45.0000000000000, 0.000000000000000, 0.000000000000000))

ellds = [(i, df2, (0., 50.)) for i in grede_e] 
print_list (ellds)

(32, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(33, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(34, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(35, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(36, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(37, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(38, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(39, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(40, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(41, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(42, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(43, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(44, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(45, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(46, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))
(47, <function distributed_force2D at 0x663b758>, (0.000000000000000,
50.0000000000000))

(32, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(33, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(34, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(35, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(36, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(37, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(38, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(39, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(40, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(41, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(42, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(43, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(44, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(45, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(46, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))
(47, <function distributed_force2D at 0x663b758>, (0.000000000000000, 50.0000000000000))

plot_frame (nds_xy, els, supps, nodal_loads = ndlds, elem_loads = ellds)

plot_frame (nds_xy, els, supps, nodal_loads = ndlds, elem_loads = ellds, load_scale = 0.5)

uu, ff = DiM (f2k, f2r, nds_xy, els, supps, chars, echs, nodal_loads = ndlds, elem_loads = ellds)

print_indexed (uu.apply_map (NDC(6)))

  0: (0.000000, 0.000000, 0.000000)
  1: (0.00236427, -0.0000193372, -0.00107080)
  2: (0.00697908, -0.0000315969, -0.00142277)
  3: (0.0124689, -0.0000280841, -0.00160954)
  4: (0.0184441, -7.49419e-6, -0.00173194)
  5: (0.0247409, 0.0000294104, -0.00182024)
  6: (0.0312558, 0.0000810878, -0.00188260)
  7: (0.0378983, 0.000145535, -0.00191533)
  8: (0.0446066, 0.000220401, -0.00196336)
  9: (0.0527081, 0.000385886, -0.00208554)
 10: (0.0607584, 0.000560595, -0.00208959)
 11: (0.0686236, 0.000737412, -0.00207715)
 12: (0.0761909, 0.000908489, -0.00203384)
 13: (0.0833531, 0.00106523, -0.00196466)
 14: (0.0900153, 0.00119831, -0.00187808)
 15: (0.0960890, 0.00129763, -0.00174516)
 16: (0.101549, 0.00135249, -0.00177972)
 17: (0.000000, 0.000000, 0.000000)
 18: (0.00237212, 0.00101239, -0.000967134)
 19: (0.00697593, 0.00195564, -0.00134197)
 20: (0.0124667, 0.00282105, -0.00152555)
 21: (0.0184408, 0.00360732, -0.00164831)
 22: (0.0247369, 0.00431521, -0.00173673)
 23: (0.0312502, 0.00494626, -0.00179768)
 24: (0.0378971, 0.00550247, -0.00184113)
 25: (0.0445902, 0.00598620, -0.00183094)
 26: (0.0527047, 0.00681377, -0.00180705)
 27: (0.0607497, 0.00750798, -0.00183667)
 28: (0.0686149, 0.00807596, -0.00182016)
 29: (0.0761812, 0.00852554, -0.00177729)
 30: (0.0833428, 0.00886532, -0.00170930)
 31: (0.0900027, 0.00910465, -0.00161470)
 32: (0.0960845, 0.00925359, -0.00153183)
 33: (0.101511, 0.00932286, -0.00127586)

  0: (0.000000, 0.000000, 0.000000)
  1: (0.00236427, -0.0000193372, -0.00107080)
  2: (0.00697908, -0.0000315969, -0.00142277)
  3: (0.0124689, -0.0000280841, -0.00160954)
  4: (0.0184441, -7.49419e-6, -0.00173194)
  5: (0.0247409, 0.0000294104, -0.00182024)
  6: (0.0312558, 0.0000810878, -0.00188260)
  7: (0.0378983, 0.000145535, -0.00191533)
  8: (0.0446066, 0.000220401, -0.00196336)
  9: (0.0527081, 0.000385886, -0.00208554)
 10: (0.0607584, 0.000560595, -0.00208959)
 11: (0.0686236, 0.000737412, -0.00207715)
 12: (0.0761909, 0.000908489, -0.00203384)
 13: (0.0833531, 0.00106523, -0.00196466)
 14: (0.0900153, 0.00119831, -0.00187808)
 15: (0.0960890, 0.00129763, -0.00174516)
 16: (0.101549, 0.00135249, -0.00177972)
 17: (0.000000, 0.000000, 0.000000)
 18: (0.00237212, 0.00101239, -0.000967134)
 19: (0.00697593, 0.00195564, -0.00134197)
 20: (0.0124667, 0.00282105, -0.00152555)
 21: (0.0184408, 0.00360732, -0.00164831)
 22: (0.0247369, 0.00431521, -0.00173673)
 23: (0.0312502, 0.00494626, -0.00179768)
 24: (0.0378971, 0.00550247, -0.00184113)
 25: (0.0445902, 0.00598620, -0.00183094)
 26: (0.0527047, 0.00681377, -0.00180705)
 27: (0.0607497, 0.00750798, -0.00183667)
 28: (0.0686149, 0.00807596, -0.00182016)
 29: (0.0761812, 0.00852554, -0.00177729)
 30: (0.0833428, 0.00886532, -0.00170930)
 31: (0.0900027, 0.00910465, -0.00161470)
 32: (0.0960845, 0.00925359, -0.00153183)
 33: (0.101511, 0.00932286, -0.00127586)

print_indexed (ff.apply_map (NDC(6)))

  0: (-85.6788, -167.585, 525.363, 85.6788, 167.585, 19.2877)
  1: (-54.3196, -151.084, 328.684, 54.3196, 151.084, 162.338)
  2: (15.5641, -150.961, 289.447, -15.5641, 150.961, 201.175)
  3: (91.2293, -146.389, 266.806, -91.2293, 146.389, 208.958)
  4: (163.516, -140.821, 249.699, -163.516, 140.821, 207.968)
  5: (228.970, -133.631, 231.888, -228.970, 133.631, 202.413)
  6: (285.553, -126.427, 213.177, -285.553, 126.427, 197.710)
  7: (331.711, -108.841, 188.217, -331.711, 108.841, 165.517)
  8: (366.614, -114.926, 194.876, -366.614, 114.926, 178.635)
  9: (387.048, -95.5672, 155.566, -387.048, 95.5672, 155.027)
 10: (391.717, -82.6259, 133.441, -391.717, 82.6259, 135.094)
 11: (379.001, -66.9688, 105.945, -379.001, 66.9688, 111.703)
 12: (347.249, -50.1784, 76.9425, -347.249, 50.1784, 86.1375)
 13: (294.806, -31.5420, 45.5017, -294.806, 31.5420, 57.0099)
 14: (220.041, -14.0447, 13.9884, -220.041, 14.0447, 31.6570)
 15: (121.529, 20.2584, -30.6230, -121.529, -20.2584, -35.2169)
 16: (4485.68, -214.915, 577.778, -4485.68, 214.915, 120.696)
 17: (4179.32, -228.604, 460.059, -4179.32, 228.604, 282.904)
 18: (3834.44, -223.102, 405.921, -3834.44, 223.102, 319.160)
 19: (3483.77, -219.236, 385.268, -3483.77, 219.236, 327.250)
 20: (3136.48, -213.554, 367.920, -3136.48, 213.554, 326.132)
 21: (2796.03, -206.681, 350.261, -2796.03, 206.681, 321.453)
 22: (2464.45, -197.011, 330.410, -2464.45, 197.011, 309.875)
 23: (2143.29, -194.909, 314.319, -2143.29, 194.909, 319.134)
 24: (1833.39, -166.324, 268.688, -1833.39, 166.324, 271.864)
 25: (1537.95, -160.370, 262.571, -1537.95, 160.370, 258.633)
 26: (1258.28, -145.187, 234.831, -1258.28, 145.187, 237.026)
 27: (995.999, -129.906, 208.249, -995.999, 129.906, 213.946)
 28: (752.751, -112.947, 179.019, -752.751, 112.947, 188.057)
 29: (530.194, -95.0205, 148.121, -530.194, 95.0205, 160.695)
 30: (329.959, -73.1428, 113.349, -329.959, 73.1428, 124.365)
 31: (153.471, -65.2584, 89.0331, -153.471, 65.2584, 123.057)
 32: (-13.6890, 31.3592, -347.971, 13.6890, -306.359, -580.754)
 33: (5.50222, 69.8836, -451.785, -5.50222, -344.884, -688.825)
 34: (3.86565, 75.6652, -467.981, -3.86565, -350.665, -704.427)
 35: (5.68184, 72.2866, -458.657, -5.68184, -347.287, -695.170)
 36: (6.87308, 65.4544, -439.857, -6.87308, -340.454, -676.392)
 37: (9.67047, 56.5825, -415.591, -9.67047, -331.582, -651.863)
 38: (2.10194, 46.1584, -385.926, -2.10194, -321.158, -624.195)
 39: (28.5853, 34.9027, -360.393, -28.5853, -309.903, -587.822)
 40: (5.95319, 20.4338, -334.201, -5.95319, -295.434, -534.435)
 41: (15.1837, 4.66933, -288.468, -15.1837, -279.669, -493.463)
 42: (15.2804, -12.7157, -241.039, -15.2804, -262.284, -445.275)
 43: (16.9596, -31.7525, -188.646, -16.9596, -243.247, -392.965)
 44: (17.9261, -52.4422, -131.639, -17.9261, -222.558, -336.178)
 45: (21.8777, -74.7650, -70.9983, -21.8777, -200.235, -274.044)
 46: (7.88434, -98.5123, -1.03404, -7.88434, -176.488, -213.398)
 47: (65.2584, -121.529, 35.2169, -65.2584, -153.471, -123.057)

  0: (-85.6788, -167.585, 525.363, 85.6788, 167.585, 19.2877)
  1: (-54.3196, -151.084, 328.684, 54.3196, 151.084, 162.338)
  2: (15.5641, -150.961, 289.447, -15.5641, 150.961, 201.175)
  3: (91.2293, -146.389, 266.806, -91.2293, 146.389, 208.958)
  4: (163.516, -140.821, 249.699, -163.516, 140.821, 207.968)
  5: (228.970, -133.631, 231.888, -228.970, 133.631, 202.413)
  6: (285.553, -126.427, 213.177, -285.553, 126.427, 197.710)
  7: (331.711, -108.841, 188.217, -331.711, 108.841, 165.517)
  8: (366.614, -114.926, 194.876, -366.614, 114.926, 178.635)
  9: (387.048, -95.5672, 155.566, -387.048, 95.5672, 155.027)
 10: (391.717, -82.6259, 133.441, -391.717, 82.6259, 135.094)
 11: (379.001, -66.9688, 105.945, -379.001, 66.9688, 111.703)
 12: (347.249, -50.1784, 76.9425, -347.249, 50.1784, 86.1375)
 13: (294.806, -31.5420, 45.5017, -294.806, 31.5420, 57.0099)
 14: (220.041, -14.0447, 13.9884, -220.041, 14.0447, 31.6570)
 15: (121.529, 20.2584, -30.6230, -121.529, -20.2584, -35.2169)
 16: (4485.68, -214.915, 577.778, -4485.68, 214.915, 120.696)
 17: (4179.32, -228.604, 460.059, -4179.32, 228.604, 282.904)
 18: (3834.44, -223.102, 405.921, -3834.44, 223.102, 319.160)
 19: (3483.77, -219.236, 385.268, -3483.77, 219.236, 327.250)
 20: (3136.48, -213.554, 367.920, -3136.48, 213.554, 326.132)
 21: (2796.03, -206.681, 350.261, -2796.03, 206.681, 321.453)
 22: (2464.45, -197.011, 330.410, -2464.45, 197.011, 309.875)
 23: (2143.29, -194.909, 314.319, -2143.29, 194.909, 319.134)
 24: (1833.39, -166.324, 268.688, -1833.39, 166.324, 271.864)
 25: (1537.95, -160.370, 262.571, -1537.95, 160.370, 258.633)
 26: (1258.28, -145.187, 234.831, -1258.28, 145.187, 237.026)
 27: (995.999, -129.906, 208.249, -995.999, 129.906, 213.946)
 28: (752.751, -112.947, 179.019, -752.751, 112.947, 188.057)
 29: (530.194, -95.0205, 148.121, -530.194, 95.0205, 160.695)
 30: (329.959, -73.1428, 113.349, -329.959, 73.1428, 124.365)
 31: (153.471, -65.2584, 89.0331, -153.471, 65.2584, 123.057)
 32: (-13.6890, 31.3592, -347.971, 13.6890, -306.359, -580.754)
 33: (5.50222, 69.8836, -451.785, -5.50222, -344.884, -688.825)
 34: (3.86565, 75.6652, -467.981, -3.86565, -350.665, -704.427)
 35: (5.68184, 72.2866, -458.657, -5.68184, -347.287, -695.170)
 36: (6.87308, 65.4544, -439.857, -6.87308, -340.454, -676.392)
 37: (9.67047, 56.5825, -415.591, -9.67047, -331.582, -651.863)
 38: (2.10194, 46.1584, -385.926, -2.10194, -321.158, -624.195)
 39: (28.5853, 34.9027, -360.393, -28.5853, -309.903, -587.822)
 40: (5.95319, 20.4338, -334.201, -5.95319, -295.434, -534.435)
 41: (15.1837, 4.66933, -288.468, -15.1837, -279.669, -493.463)
 42: (15.2804, -12.7157, -241.039, -15.2804, -262.284, -445.275)
 43: (16.9596, -31.7525, -188.646, -16.9596, -243.247, -392.965)
 44: (17.9261, -52.4422, -131.639, -17.9261, -222.558, -336.178)
 45: (21.8777, -74.7650, -70.9983, -21.8777, -200.235, -274.044)
 46: (7.88434, -98.5123, -1.03404, -7.88434, -176.488, -213.398)
 47: (65.2584, -121.529, 35.2169, -65.2584, -153.471, -123.057)

def nodes3D_as_supports (nd_ids) : 
    return [(ni, (0, 0, 0)) for ni in nd_ids]

nds_xy = [(0., 0., 0.), (2., 0., 0.), (0., 2., 0.), (0., 0., 2.)]
print "joints:"
print_indexed (nds_xy)
print

els = [(0, 3), (1, 3), (2, 3)]
print "bars:"
print_indexed (els)
print

supps = nodes3D_as_supports ([0, 1, 2])
print "supports:", supps
print

chars = [(2*1e8, 0.001)]
print 'characteristics:', chars
print

echs = [0, 0, 0]

ndlds = [(3, (100., 100., 0.))]
print "loads:" 
print_list (ndlds, '   ')
print

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)

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

characteristics: [(2.00000000000000e8, 0.00100000000000000)]

loads:
   (3, (100.000000000000, 100.000000000000, 0.000000000000000))

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)

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

characteristics: [(2.00000000000000e8, 0.00100000000000000)]

loads:
   (3, (100.000000000000, 100.000000000000, 0.000000000000000))

uu, ft = DiM (t3k, t3r, nds_xy, els, supps, chars, echs, nodal_loads = ndlds)

print_indexed (ft.apply_map (NDC(6)))

  0: (-200.000, 200.000)
  1: (141.421, -141.421)
  2: (141.421, -141.421)

  0: (-200.000, 200.000)
  1: (141.421, -141.421)
  2: (141.421, -141.421)

nds_xy = [(0., 0.), (3., 4.), (8., 4.), (8., 0.)]
print 'nodes:'
print_indexed (nds_xy)
print 

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

supps = [(0, (0, 0, 0)), (3, (0, 0, 0))]
print 'supports:', supps
print

chars = [(3e7, 0.5*0.5, 0.5^4/12)]
print 'characteristics:'
print_list (chars, '   ')
print

echs = [0, 0, 0]
print 'elements <-> chars:', echs
print 

ndlds = [(1, (100., 0., 0.))]
print 'nodal loads:'
print_list (ndlds, '   ')
print

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (3.00000000000000, 4.00000000000000)
  2: (8.00000000000000, 4.00000000000000)
  3: (8.00000000000000, 0.000000000000000)

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

supports: [(0, (0, 0, 0)), (3, (0, 0, 0))]

characteristics:
   (3.00000000000000e7, 0.250000000000000, 0.00520833333333333)

elements <-> chars: [0, 0, 0]

nodal loads:
   (1, (100.000000000000, 0.000000000000000, 0.000000000000000))

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (3.00000000000000, 4.00000000000000)
  2: (8.00000000000000, 4.00000000000000)
  3: (8.00000000000000, 0.000000000000000)

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

supports: [(0, (0, 0, 0)), (3, (0, 0, 0))]

characteristics:
   (3.00000000000000e7, 0.250000000000000, 0.00520833333333333)

elements <-> chars: [0, 0, 0]

nodal loads:
   (1, (100.000000000000, 0.000000000000000, 0.000000000000000))

plot_frame (nds_xy, els, supps, nodal_loads = ndlds)

uu, ff = DiM (f2k, f2r, nds_xy, els, supps, chars, echs, nodal_loads = ndlds)

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

  0: (0.00000000, 0.00000000, 0.00000000)
  1: (0.0017851950, 0.0012900092, -0.000095396427)
  2: (0.0017577034, 0.000015604534, -0.00017506838)
  3: (0.00000000, 0.00000000, 0.00000000)

  0: (0.00000000, 0.00000000, 0.00000000)
  1: (0.0017851950, 0.0012900092, -0.000095396427)
  2: (0.0017577034, 0.000015604534, -0.00017506838)
  3: (0.00000000, 0.00000000, 0.00000000)

print_indexed (ff.apply_map (NDC(6)))

  0: (-58.6644, -29.4551, 76.6188, 58.6644, 29.4551, 70.6565)
  1: (41.2373, 29.2585, -70.6565, -41.2373, -29.2585, -75.6360)
  2: (29.2585, -41.2373, 75.6360, -29.2585, 41.2373, 89.3132)

  0: (-58.6644, -29.4551, 76.6188, 58.6644, 29.4551, 70.6565)
  1: (41.2373, 29.2585, -70.6565, -41.2373, -29.2585, -75.6360)
  2: (29.2585, -41.2373, 75.6360, -29.2585, 41.2373, 89.3132)

uu, ff = IDiM (f2k, f2r, nds_xy, els, supps, chars, echs, nodal_loads = ndlds)

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

  0: (0.00000000, 0.00000000, 0.00000000)
  1: (0.0017463184, 0.0013097388, -0.000088599976)
  2: (0.0017463184, -0.00000000, -0.00016949561)
  3: (0.00000000, 0.00000000, 0.00000000)

  0: (0.00000000, 0.00000000, 0.00000000)
  1: (0.0017463184, 0.0013097388, -0.000088599976)
  2: (0.0017463184, -0.00000000, -0.00016949561)
  3: (0.00000000, 0.00000000, 0.00000000)

print_indexed (ff.apply_map (NDC(6)))

  0: (0.000000, -29.4210, 76.3212, 0.000000, 29.4210, 70.7837)
  1: (0.000000, 29.3247, -70.7837, 0.000000, -29.3247, -75.8397)
  2: (0.000000, -41.2303, 75.8397, 0.000000, 41.2303, 89.0815)

  0: (0.000000, -29.4210, 76.3212, 0.000000, 29.4210, 70.7837)
  1: (0.000000, 29.3247, -70.7837, 0.000000, -29.3247, -75.8397)
  2: (0.000000, -41.2303, 75.8397, 0.000000, 41.2303, 89.0815)

nds = coords_left_to_right (nds_xy)

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

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

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

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

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 2.04960 \times 10^{6} & -712800. & 30000.0 & -1.50000 \times 10^{6} & 0.000000 & 0.000000 \\ -712800. & 980400. & -15000.0 & 0.000000 & -15000.0 & -37500.0 \\ 30000.0 & -15000.0 & 250000. & 0.000000 & 37500.0 & 62500.0 \\ -1.50000 \times 10^{6} & 0.000000 & 0.000000 & 1.52930 \times 10^{6} & 0.000000 & 58593.8 \\ 0.000000 & -15000.0 & 37500.0 & 0.000000 & 1.89000 \times 10^{6} & 37500.0 \\ 0.000000 & -37500.0 & 62500.0 & 58593.8 & 37500.0 & 281250. \end{array}\right)$

q = structure_force_vector (ndof, tdof, nodal_loads = ndlds)
show (q.apply_map (ND(4)))

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

ntdof, ttdof = transl_DsOF (tdof)
ntdof; ttdof

4
[-1 -1]
[ 0  1]
[ 2  3]
[-1 -1]

4
[-1 -1]
[ 0  1]
[ 2  3]
[-1 -1]

ms, C = kinematic_condensation_matrix (ndof, tdof, nds, els)
show (C.apply_map (ND(4)))

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

print_masters (ms, ttdof)

masters:
   node 2 :  horizontal displacement

masters:
   node 2 :  horizontal displacement

CtKC = C.T * K * C
show (CtKC.apply_map (ND(8)))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 250000.00 & 18750.000 & 62500.000 \\ 18750.000 & 61171.875 & 30468.750 \\ 62500.000 & 30468.750 & 281250.00 \end{array}\right)$

Ctq = C.T * q
show (Ctq.apply_map (ND(4)))

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

um = CtKC.solve_right (Ctq)
um

(-0.0000885999759239197, 0.00174631836603668, -0.000169495606115325)

(-0.0000885999759239197, 0.00174631836603668, -0.000169495606115325)

ua = C * um
ua

(0.00174631836603668, 0.00130973877452751, -0.0000885999759239197,
0.00174631836603668, 0.000000000000000, -0.000169495606115325)

(0.00174631836603668, 0.00130973877452751, -0.0000885999759239197, 0.00174631836603668, 0.000000000000000, -0.000169495606115325)

uu = table_of_nodal_displacements (tdof, ua)
print_indexed (uu.apply_map (ND(8)))

  0: (0.00000000, 0.00000000, 0.00000000)
  1: (0.0017463184, 0.0013097388, -0.000088599976)
  2: (0.0017463184, 0.00000000, -0.00016949561)
  3: (0.00000000, 0.00000000, 0.00000000)

  0: (0.00000000, 0.00000000, 0.00000000)
  1: (0.0017463184, 0.0013097388, -0.000088599976)
  2: (0.0017463184, 0.00000000, -0.00016949561)
  3: (0.00000000, 0.00000000, 0.00000000)

ff = all_end_forces (f2k, f2r, uu, nds, els, echs, chars)
print_indexed (ff.apply_map (NDC(6)))

  0: (0.000000, -29.4210, 76.3212, 0.000000, 29.4210, 70.7837)
  1: (0.000000, 29.3247, -70.7837, 0.000000, -29.3247, -75.8397)
  2: (0.000000, -41.2303, 75.8397, 0.000000, 41.2303, 89.0815)

  0: (0.000000, -29.4210, 76.3212, 0.000000, 29.4210, 70.7837)
  1: (0.000000, 29.3247, -70.7837, 0.000000, -29.3247, -75.8397)
  2: (0.000000, -41.2303, 75.8397, 0.000000, 41.2303, 89.0815)

nds_xy = [(0., 0.), (10., 1.), (2., 4.), (7., 5.), (2., 8.), (7., 8.)]
print 'nodes:'
print_indexed (nds_xy)
print 

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

supps = [(0, (0, 0, 0)), (1, (0, 0, 0))]
print 'supports:', supps
print

chars = [(3e7, 0.5*0.5, 0.5^4/12)]
print 'characteristics:'
print_list (chars, '   ')
print

echs = [0, 0, 0, 0, 0, 0]
print 'elements <-> chars:', echs
print 

ndlds = [(2, (100., 0., 0.)), (4, (100., 0., 0.))]
print 'nodal loads:'
print_list (ndlds, '   ')
print

ellds = [(1, cf2, (0., 50., 0., 1/3)), (1, cf2, (0., 50., 0., 2/3)), (4, cf2, (0., 50., 0., 1/2))]
print 'element loads:'
print_list (ellds, '   ')

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (10.0000000000000, 1.00000000000000)
  2: (2.00000000000000, 4.00000000000000)
  3: (7.00000000000000, 5.00000000000000)
  4: (2.00000000000000, 8.00000000000000)
  5: (7.00000000000000, 8.00000000000000)

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

supports: [(0, (0, 0, 0)), (1, (0, 0, 0))]

characteristics:
   (3.00000000000000e7, 0.250000000000000, 0.00520833333333333)

elements <-> chars: [0, 0, 0, 0, 0, 0]

nodal loads:
   (2, (100.000000000000, 0.000000000000000, 0.000000000000000))
   (4, (100.000000000000, 0.000000000000000, 0.000000000000000))

element loads:
   (1, <function concentrated_force2D at 0xadc6725c>, (0.000000000000000,
50.0000000000000, 0.000000000000000, 1/3))
   (1, <function concentrated_force2D at 0xadc6725c>, (0.000000000000000,
50.0000000000000, 0.000000000000000, 2/3))
   (4, <function concentrated_force2D at 0xadc6725c>, (0.000000000000000,
50.0000000000000, 0.000000000000000, 1/2))

nodes:
  0: (0.000000000000000, 0.000000000000000)
  1: (10.0000000000000, 1.00000000000000)
  2: (2.00000000000000, 4.00000000000000)
  3: (7.00000000000000, 5.00000000000000)
  4: (2.00000000000000, 8.00000000000000)
  5: (7.00000000000000, 8.00000000000000)

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

supports: [(0, (0, 0, 0)), (1, (0, 0, 0))]

characteristics:
   (3.00000000000000e7, 0.250000000000000, 0.00520833333333333)

elements <-> chars: [0, 0, 0, 0, 0, 0]

nodal loads:
   (2, (100.000000000000, 0.000000000000000, 0.000000000000000))
   (4, (100.000000000000, 0.000000000000000, 0.000000000000000))

element loads:
   (1, <function concentrated_force2D at 0xadc6725c>, (0.000000000000000, 50.0000000000000, 0.000000000000000, 1/3))
   (1, <function concentrated_force2D at 0xadc6725c>, (0.000000000000000, 50.0000000000000, 0.000000000000000, 2/3))
   (4, <function concentrated_force2D at 0xadc6725c>, (0.000000000000000, 50.0000000000000, 0.000000000000000, 1/2))

plot_frame (nds_xy, els, supps, nodal_loads = ndlds, elem_loads = ellds)

uu, ff = DiM (f2k, f2r, nds_xy, els, supps, chars, echs, nodal_loads = ndlds, elem_loads = ellds)

print_indexed (uu.apply_map (NDC(8)))

  0: (0.00000000, 0.00000000, 0.00000000)
  1: (0.00000000, 0.00000000, 0.00000000)
  2: (0.0027947712, 0.0013738427, -0.00038016241)
  3: (0.0021893491, -0.0014557156, -0.00014215334)
  4: (0.0037396186, 0.0013674577, 0.00018142955)
  5: (0.0036836430, -0.0014309268, -0.000047993723)

  0: (0.00000000, 0.00000000, 0.00000000)
  1: (0.00000000, 0.00000000, 0.00000000)
  2: (0.0027947712, 0.0013738427, -0.00038016241)
  3: (0.0021893491, -0.0014557156, -0.00014215334)
  4: (0.0037396186, 0.0013674577, 0.00018142955)
  5: (0.0036836430, -0.0014309268, -0.000047993723)

print_indexed (ff.apply_map (NDC(6)))

  0: (-35.3143, -47.4616, 119.410, 35.3143, 47.4616, 92.8451)
  1: (56.9863, 9.75403, -102.981, -56.9863, -109.754, -201.706)
  2: (223.555, -34.0429, 80.6649, -223.555, 34.0429, 89.5495)
  3: (-11.9719, -16.0366, 10.1360, 11.9719, 16.0366, 54.0103)
  4: (83.9634, 11.9719, -54.0103, -83.9634, -61.9719, -130.849)
  5: (61.9719, -83.9634, 130.849, -61.9719, 83.9634, 121.041)

  0: (-35.3143, -47.4616, 119.410, 35.3143, 47.4616, 92.8451)
  1: (56.9863, 9.75403, -102.981, -56.9863, -109.754, -201.706)
  2: (223.555, -34.0429, 80.6649, -223.555, 34.0429, 89.5495)
  3: (-11.9719, -16.0366, 10.1360, 11.9719, 16.0366, 54.0103)
  4: (83.9634, 11.9719, -54.0103, -83.9634, -61.9719, -130.849)
  5: (61.9719, -83.9634, 130.849, -61.9719, 83.9634, 121.041)

uu, ff = IDiM (f2k, f2r, nds_xy, els, supps, chars, echs, nodal_loads = ndlds, elem_loads = ellds)

print_indexed (uu.apply_map (NDC(8)))

  0: (0.00000000, 0.00000000, 0.00000000)
  1: (0.00000000, 0.00000000, 0.00000000)
  2: (0.0027001656, 0.0013500828, -0.00035387188)
  3: (0.0021131731, -0.0015848798, -0.00013200118)
  4: (0.0035509470, 0.0013500828, 0.00020633752)
  5: (0.0035509470, -0.0015848798, -0.000022808979)

  0: (0.00000000, 0.00000000, 0.00000000)
  1: (0.00000000, 0.00000000, 0.00000000)
  2: (0.0027001656, 0.0013500828, -0.00035387188)
  3: (0.0021131731, -0.0015848798, -0.00013200118)
  4: (0.0035509470, 0.0013500828, 0.00020633752)
  5: (0.0035509470, -0.0015848798, -0.000022808979)

print_indexed (ff.apply_map (NDC(6)))

  0: (0.000000, -46.6974, 116.782, 0.000000, 46.6974, 92.0548)
  1: (0.000000, 9.85065, -102.733, 0.000000, -109.851, -202.447)
  2: (0.000000, -34.6720, 82.5548, 0.000000, 34.6720, 90.8049)
  3: (0.000000, -16.2806, 10.6781, 0.000000, 16.2806, 54.4445)
  4: (0.000000, 12.1421, -54.4445, 0.000000, -62.1421, -131.266)
  5: (0.000000, -83.7194, 131.266, 0.000000, 83.7194, 119.892)

  0: (0.000000, -46.6974, 116.782, 0.000000, 46.6974, 92.0548)
  1: (0.000000, 9.85065, -102.733, 0.000000, -109.851, -202.447)
  2: (0.000000, -34.6720, 82.5548, 0.000000, 34.6720, 90.8049)
  3: (0.000000, -16.2806, 10.6781, 0.000000, 16.2806, 54.4445)
  4: (0.000000, 12.1421, -54.4445, 0.000000, -62.1421, -131.266)
  5: (0.000000, -83.7194, 131.266, 0.000000, 83.7194, 119.892)

GS_program_DiM_(v0.4.1)

53 days ago by fresl